【发布时间】：2021-10-16 00:19:55
【问题描述】：
我想将 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）
- 应用 gamma 校正（现在是简单的例程，应该替换为 sRGB gamma）
- 从 [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)。
【问题讨论】：
-
根据下面的guide，你有tp使用
1/AsShotNeutral（但指南不准确）。您可以关注Developing A RAW Photo By Hand 第 1 部分和第 2 部分。注意：该站点看起来已损坏 - 为了查看图像，我必须在浏览器中放大和缩小。如果实现是在 Python 或 MATLAB 中，我本可以为您提供更好的帮助。 -
我试图按照你的代码。我认为问题出在“从 [minval..maxval] 重新调整为 [0..1]”阶段。您不必跟踪最小值/最大值：
maxval = max!(maxval, r, g, b);。代替clip(r, minval, maxval)，使用clip(r, 0.0, 1.0)（g和b相同）。 -
@Rotem 跟踪是因为我已经读到仅剪切 0..1 可能会导致高光中出现洋红色，因此您必须跟踪最大值并将所有通道减少相同的数量只剪掉一两个溢出的通道。不确定我的方法是否正确。在您链接的文章中，有一个很好的提示：由于 WB 为 ~1.9、1.0、~1.8，因此必须对这些进行缩放，因此没有一个 >1.0。我已经添加了这个并且还添加了 clip(..., 0.0, 1.0) 但结果仍然有一个极端的偏色。现在它是红色/洋红色而不是绿色，因为 1,9、1,0、1,8 的绿色部分要低得多。
-
我猜问题出在“计算 D65 光源 XYZ_TO_CAM 的伪逆”阶段。根据odelama，您应该使用
ForwardMatrix1和ForwardMatrix2（尝试仅使用ForwardMatrix2）。首先将其从 XYZD50 转换为 XYZD65。 -
“从 [minval..maxval] 重新调整为 [0..1]”意味着您撤消白平衡和黑电平校正的任何影响。它基本上使您早期的大部分计算步骤变得无关紧要。如果您确实想缩放以避免剪裁，请平等缩放所有通道，并在原处保留 0（将负值剪裁为 0）。
标签： image image-processing rust rgb raw