我想将 RAW 图像数据 (RGGB) 转换为 sRGB 图像。有许多专门的方法可以做到这一点，但首先要了解基础知识，我已经实现了一些简单的算法，例如通过降低分辨率进行去拜耳法。我目前的管道是：
- 通过blacklevel和whitelevel 重新缩放u16 输入数据
- 应用白平衡系数
- 尺寸减小的 Debayer，G 的平均值：g=((g0+g1)/2)
- 计算 D65 光源 XYZ_TO_CAM 的伪逆（来自 Adobe DNG）
- 通过 CAM_TO_XYZ 将 debayered RGB 数据转换为 XYZ
- 将 XYZ 转换为 D65 sRGB（取自 Bruce Lindbloom 的矩阵）
- 应用伽玛校正（现在简单的例程，应该用 sRGB 伽玛代替）
- 从 [minval..maxval] 重新缩放到 [0..1] 并将 f32 转换为 u16
- 另存为 tiff
问题是，如果我跳过白平衡系数乘法（或者只是将它们替换为 1.0），输出图像看起来已经可以接受了。如果我应用系数（取自 DNG 中的 AsShot），输出会有很大的色偏。而且我不确定是否必须乘以coef或1/coef。
第一个图像是 wb_coefs 设置为 1.0 的管道的结果。
第二个图像是具有“正确” wb_coefs 的结果。
我的管道有什么问题？
附加问题：
- 我不确定重新缩放过程。我是否必须在每一步之后重新缩放到 [0..1] 还是在 u16 转换作为最后阶段时重新缩放是否足够？
完整代码：
macro_rules! max {
($x: expr) => ($x);
($x: expr, $($z: expr),+) => {{
let y = max!($($z),*);
if $x > y {
$x
} else {
y
}
}}
}
macro_rules! min {
($x: expr) => ($x);
($x: expr, $($z: expr),+) => {{
let y = min!($($z),*);
if $x < y {
$x
} else {
y
}
}}
}
/// sRGB D65
const XYZD65_TO_SRGB: [[f32; 3]; 4] = [
[3.2404542, -1.5371385, -0.4985314],
[-0.9692660, 1.8760108, 0.0415560],
[0.0556434, -0.2040259, 1.0572252],
[0.0, 0.0, 0.0],
];
// buf: RAW image data
fn to_srgb(buf: &Vec<u16>, width: usize, height: usize) {
let w = width / 2;
let h = height / 2;
let blacklevel: [u16; 4] = [511, 511, 511, 511];
let whitelevel: [u16; 4] = [12735, 12735, 12735, 12735];
let xyz2cam_d65: [[i32; 3]; 4] = [[6722, -635, -963], [-4287, 12460, 2028], [-908, 2162, 5668], [0, 0, 0]];
let cam2xyz = convert_matrix::<4>(xyz2cam_d65);
eprintln!("CAM_TO_XYZ: {:?}", cam2xyz);
// from DNG
// As Shot Neutral: 0.518481 1 0.545842
//let wb_coef = [1.0/0.518481, 1.0, 1.0, 1.0/0.545842];
//let wb_coef = [0.518481, 1.0, 1.0, 0.545842];
let wb_coef = [1.0, 1.0, 1.0, 1.0];
// b/w level correction, rescale, debayer
let mut rgb = vec![0.0_f32; width / 2 * height / 2 * 3];
for row in 0..h {
for col in 0..w {
let r0 = buf[(row * 2 + 0) * width + (col * 2) + 0];
let g0 = buf[(row * 2 + 0) * width + (col * 2) + 1];
let g1 = buf[(row * 2 + 1) * width + (col * 2) + 0];
let b0 = buf[(row * 2 + 1) * width + (col * 2) + 1];
let r0 = ((r0.saturating_sub(blacklevel[0])) as f32 / (whitelevel[0] - blacklevel[0]) as f32) * wb_coef[0];
let g0 = ((g0.saturating_sub(blacklevel[1])) as f32 / (whitelevel[1] - blacklevel[1]) as f32) * wb_coef[1];
let g1 = ((g1.saturating_sub(blacklevel[2])) as f32 / (whitelevel[2] - blacklevel[2]) as f32) * wb_coef[2];
let b0 = ((b0.saturating_sub(blacklevel[3])) as f32 / (whitelevel[3] - blacklevel[3]) as f32) * wb_coef[3];
rgb[row * w * 3 + (col * 3) + 0] = r0;
rgb[row * w * 3 + (col * 3) + 1] = (g0 + g1) / 2.0;
rgb[row * w * 3 + (col * 3) + 2] = b0;
}
}
// Convert to XYZ by CAM_TO_XYZ from D65 illuminant
let mut xyz = vec![0.0_f32; w * h * 3];
for row in 0..h {
for col in 0..w {
let r = rgb[row * w * 3 + (col * 3) + 0];
let g = rgb[row * w * 3 + (col * 3) + 1];
let b = rgb[row * w * 3 + (col * 3) + 2];
xyz[row * w * 3 + (col * 3) + 0] = cam2xyz[0][0] * r + cam2xyz[0][1] * g + cam2xyz[0][2] * b;
xyz[row * w * 3 + (col * 3) + 1] = cam2xyz[1][0] * r + cam2xyz[1][1] * g + cam2xyz[1][2] * b;
xyz[row * w * 3 + (col * 3) + 2] = cam2xyz[2][0] * r + cam2xyz[2][1] * g + cam2xyz[2][2] * b;
}
}
// Track min/max value for rescaling/clipping
let mut maxval = 1.0;
let mut minval = 0.0;
// Convert to sRGB from XYZ
let mut srgb = vec![0.0; w * h * 3];
for row in 0..h {
for col in 0..w {
let r = xyz[row * w * 3 + (col * 3) + 0] as f32;
let g = xyz[row * w * 3 + (col * 3) + 1] as f32;
let b = xyz[row * w * 3 + (col * 3) + 2] as f32;
srgb[row * w * 3 + (col * 3) + 0] = XYZD65_TO_SRGB[0][0] * r + XYZD65_TO_SRGB[0][1] * g + XYZD65_TO_SRGB[0][2] * b;
srgb[row * w * 3 + (col * 3) + 1] = XYZD65_TO_SRGB[1][0] * r + XYZD65_TO_SRGB[1][1] * g + XYZD65_TO_SRGB[1][2] * b;
srgb[row * w * 3 + (col * 3) + 2] = XYZD65_TO_SRGB[2][0] * r + XYZD65_TO_SRGB[2][1] * g + XYZD65_TO_SRGB[2][2] * b;
let r = srgb[row * w * 3 + (col * 3) + 0];
let g = srgb[row * w * 3 + (col * 3) + 1];
let b = srgb[row * w * 3 + (col * 3) + 2];
maxval = max!(maxval, r, g, b);
minval = min!(minval, r, g, b);
}
}
gamma_corr(&mut srgb, w, h, 2.2);
let mut output = vec![0_u16; w * h * 3];
for row in 0..h {
for col in 0..w {
let r = srgb[row * w * 3 + (col * 3) + 0];
let g = srgb[row * w * 3 + (col * 3) + 1];
let b = srgb[row * w * 3 + (col * 3) + 2];
output[row * w * 3 + (col * 3) + 0] = (clip(r, minval, maxval) * (u16::MAX as f32)) as u16;
output[row * w * 3 + (col * 3) + 1] = (clip(g, minval, maxval) * (u16::MAX as f32)) as u16;
output[row * w * 3 + (col * 3) + 2] = (clip(b, minval, maxval) * (u16::MAX as f32)) as u16;
}
}
let img = DynamicImage::ImageRgb16(ImageBuffer::from_raw(w as u32, h as u32, output).unwrap());
img.save_with_format("/tmp/test.tif", image::ImageFormat::Tiff).unwrap();
}
fn pseudoinverse<const N: usize>(matrix: [[f32; 3]; N]) -> [[f32; 3]; N] {
let mut result: [[f32; 3]; N] = [Default::default(); N];
let mut work: [[f32; 6]; 3] = [Default::default(); 3];
let mut num: f32 = 0.0;
for i in 0..3 {
for j in 0..6 {
work[i][j] = if j == i + 3 { 1.0 } else { 0.0 };
}
for j in 0..3 {
for k in 0..N {
work[i][j] += matrix[k][i] * matrix[k][j];
}
}
}
for i in 0..3 {
num = work[i][i];
for j in 0..6 {
work[i][j] /= num;
}
for k in 0..3 {
if k == i {
continue;
}
num = work[k][i];
for j in 0..6 {
work[k][j] -= work[i][j] * num;
}
}
}
for i in 0..N {
for j in 0..3 {
result[i][j] = 0.0;
for k in 0..3 {
result[i][j] += work[j][k + 3] * matrix[i][k];
}
}
}
result
}
fn convert_matrix<const N: usize>(adobe_xyz_to_cam: [[i32; 3]; N]) -> [[f32; N]; 3] {
let mut xyz_to_cam: [[f32; 3]; N] = [[0.0; 3]; N];
let mut cam_to_xyz: [[f32; N]; 3] = [[0.0; N]; 3];
for i in 0..N {
for j in 0..3 {
xyz_to_cam[i][j] = adobe_xyz_to_cam[i][j] as f32 / 10000.0;
}
}
eprintln!("XYZ_TO_CAM: {:?}", xyz_to_cam);
let inverse = pseudoinverse::<N>(xyz_to_cam);
for i in 0..3 {
for j in 0..N {
cam_to_xyz[i][j] = inverse[j][i];
}
}
cam_to_xyz
}
fn clip(v: f32, minval: f32, maxval: f32) -> f32 {
(v + minval.abs()) / (maxval + minval.abs())
}
// https://kosinix.github.io/raster/docs/src/raster/filter.rs.html#339-359
fn gamma_corr(rgb: &mut Vec<f32>, w: usize, h: usize, gamma: f32) {
for row in 0..h {
for col in 0..w {
let r = rgb[row * w * 3 + (col * 3) + 0];
let g = rgb[row * w * 3 + (col * 3) + 1];
let b = rgb[row * w * 3 + (col * 3) + 2];
rgb[row * w * 3 + (col * 3) + 0] = r.powf(1.0 / gamma);
rgb[row * w * 3 + (col * 3) + 1] = g.powf(1.0 / gamma);
rgb[row * w * 3 + (col * 3) + 2] = b.powf(1.0 / gamma);
}
}
}
此示例的 DNG 可在以下位置找到：https ://chaospixel.com/pub/misc/dng/sample.dng (~40 MiB)。
