如何将白平衡系数应用于sRGB输出的RAW图像

html5 • • 问答

我想将 RAW 图像数据 (RGGB) 转换为 sRGB 图像。有许多专门的方法可以做到这一点，但首先要了解基础知识，我已经实现了一些简单的算法，例如通过降低分辨率进行去拜耳化。我目前的管道是：

  • 通过重新调节U16输入数据的黑和whitelevel
  • 应用白平衡系数
  • 尺寸减小的 Debayer，G 的平均值：g=((g0+g1)/2)
  • 计算 D65 光源 XYZ_TO_CAM 的伪逆（来自 Adob​​e 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)。

回答

得到错误颜色的主要原因是我们必须将rgb2cam矩阵的行标准化为1，如以下指南中所述。

根据DNG 规范：

ColorMatrix1 定义了一个转换矩阵，该矩阵在第一个校准光源下将 XYZ 值转换为参考相机原始色彩空间值。

这意味着如果校准光源是 D65，则 ColorMatrix 将 XYZ 转换为“相机 RGB”。
（按原样转换，不使用任何白平衡缩放系数）。

  • 逆 ColorMatrix，从“相机 RGB”转换为 XYZ。
    将 XYZ 转换为 sRGB 后，结果是色彩平衡的sRGB。
    结论是，ColorMatrix 中包含了while 平衡系数（白平衡系数适用于D65 光源）。
  • 对 的行进行归一化rgb2cam会1抵消 while 平衡系数，并且只保留“颜色校正矩阵”（数学有点复杂）。
  • 没有对行进行标准化，我们将同时平衡乘数两次：
  1. 来自 ColorMatrix 的比例系数，用于平衡 D65 的输入。
  2. 取自 AsShotNatural 的比例系数，用于平衡场景光源的输入（场景光源接近 D65）。
    缩放两次的结果是极端的色偏。

跟踪最大值以避免“高光中的洋红色”：

我们假设跟踪“理论最大颜色值”，而不是跟踪输入图像中的实际最大颜色值。

  • whitelevel - blacklevel通过白平衡乘数进行拍摄和缩放。
    跟踪结果...

指导原则是两种情况下的颜色应该相同：

  • 将处理应用于图像的小块，并将这些块放在一起（我们无法跟踪全局最小值和最大值）。
  • 将处理应用于整个图像。

我想您必须跟踪缩放的最大值whitelevel - blacklevel，只有当白平衡乘数小于1.
当所有乘数都为1或以上时，我们可以将结果裁剪为1.0，而不跟踪最大值。
笔记：

  • 缩小规模并跟踪最大值可能有优势，但我不知道这个主题。
    在我的解决方案中，我们只是乘以 upper（高于 1.0），然后裁剪结果。

该解决方案基于在 MATLAB 中处理 RAW 图像指南中。

我发布了 MATLAB 实现和 Python 实现（但没有 Rust 实现）。

第一步是sample.dng使用dcraw命令行提取原始拜耳图像：

dcraw -4 -D -T sample.dng  

dcraw -4 -D -T sample.dng

将 tiff 输出重命名为sample__lin_bayer.tif.

转换过程：

  • 重新调整uint16由输入数据blacklevel和whitelevel（减去blacklevel从所有像素和规模通过whitelevel - blacklevel）。
  • 应用白平衡缩放系数。
    缩放系数等于1./AsShotNatural。
    按红色缩放系数缩放拜耳对齐中的红色像素，按绿色缩放比例缩放绿色，按蓝色缩放比例缩放蓝色。
    假设：最小比例是1.0，其他人在上面1.0（我们除以最小比例以确保）。
  • 将缩放结果裁剪为 [0, 1]（由于demosaic实现限制，需要裁剪）。
  • 使用 MATLABdemosaic函数或cv2.cvtColorPython 进行去马赛克 (Debayer) 。
  • 计算rgb2cam矩阵：rgb2cam = ColorMatrix * rgb2xyz。
    rgb2xyz矩阵取自Bruce Lindbloom 站点。
  • 规范化rgb2cam矩阵的行，使每行的总和等于1（每行除以行的总和）。
  • cam2rgb通过求逆计算矩阵rgb2cam: cam2rgb = inv(rgb2cam)。
    cam2rgb是“CCM矩阵”（色彩校正矩阵）。
  • 左乘矩阵乘以cam2rgb每个 RGB 元组（应用颜色校正）。
  • 应用伽马校正（使用 sRGB 标准伽马）。
  • 转换uint8并保存为 PNG（PNG 格式用于在 SO 网站上发布）。

MATLAB 实现：

Python实现：

filename = 'sample__lin_bayer.tif'; % Output of: dcraw -4 -D -T sample.dng

% Exif information:
blacklevel = 511; % blacklevel = meta_info.SubIFDs{1}.BlackLevel(1);
whitelevel = 12735; % whitelevel = meta_info.SubIFDs{1}.WhiteLevel;
AsShotNeutral = [0.5185 1 0.5458];
ColorMatrix = [ 0.6722   -0.0635   -0.0963
               -0.4287    1.2460    0.2028
               -0.0908    0.2162    0.5668];
bayer_type = 'rggb';


% Constant matrix for converting sRGB to XYZ(D65):
% http://www.brucelindbloom.com/Eqn_RGB_XYZ_Matrix.html
rgb2xyz = [0.4124564  0.3575761  0.1804375
           0.2126729  0.7151522  0.0721750
           0.0193339  0.1191920  0.9503041];


% Read input image (Bayer mosaic alignment, pixels data type is uint16):
raw = imread(filename);

% "Linearizing":
% There is no LinearizationTable so we only have to subtract the black level.
% Convert from range [blacklevel, whitelevel] to range [0, 1].
lin_bayer = (double(raw) - blacklevel) / (whitelevel - blacklevel);
lin_bayer = max(0, min(lin_bayer, 1));

% The White Balance multipliers are 1./AsShotNeutral
wb_multipliers = 1./AsShotNeutral;
r_scale = wb_multipliers(1); % Assume value is above 1
g_scale = wb_multipliers(2); % Assume value = 1
b_scale = wb_multipliers(3); % Assume value is above 1

% Bayer alignment is RGGB:
% R G
% G B
%
% Apply white balancing to linear Bayer image.
balanced_bayer = lin_bayer;
balanced_bayer(1:2:end, 1:2:end) = balanced_bayer(1:2:end, 1:2:end)*r_scale; % Red   (indices [1, 3, 5,... ; 1, 3, 5,... ])
balanced_bayer(1:2:end, 2:2:end) = balanced_bayer(1:2:end, 2:2:end)*g_scale; % Green (indices [1, 3, 5,... ; 2, 4, 6,... ])
balanced_bayer(2:2:end, 1:2:end) = balanced_bayer(2:2:end, 1:2:end)*g_scale; % Green (indices [2, 4, 6,... ; 1, 3, 5,... ])
balanced_bayer(2:2:end, 2:2:end) = balanced_bayer(2:2:end, 2:2:end)*b_scale; % Blue  (indices [2, 4, 6,... ; 2, 4, 6,... ])

% Clip to range [0, 1] for avoiding "pinkish highlights" (avoiding "magenta cast" in the highlights).
balanced_bayer = min(balanced_bayer, 1);

% Demosaicing
temp = uint16(balanced_bayer*(2^16-1)); % Convert from double to uint16, because MATLAB demosaic() function requires a uint8 or uint16 input.
lin_rgb = double(demosaic(temp, bayer_type))/(2^16-1); % Apply Demosaicing and convert range back type double and range [0, 1].

% Color Space Conversion
xyz2cam = ColorMatrix; % ColorMatrix applies XYZ(D65) to CAM_rgb
rgb2cam = xyz2cam * rgb2xyz;

% Result:
% rgb2cam = [0.2619    0.1835    0.0252
%            0.0921    0.7620    0.2053
%            0.0195    0.1897    0.5379]

% Normalize rows to 1. MATLAB shortcut: rgb2cam = rgb2cam ./ repmat(sum(rgb2cam,2),1,3);

rows_sum = sum(rgb2cam, 2);
% Result:
% rows_sum = [0.4706
%             1.0593
%             0.7470]

% Divide element of every row by the sum of the row:
rgb2cam(1, :) = rgb2cam(1, :) / rows_sum(1); % Divide top row
rgb2cam(2, :) = rgb2cam(2, :) / rows_sum(2); % Divide center row
rgb2cam(3, :) = rgb2cam(3, :) / rows_sum(3); % Divide bottom row
% Result (sum of every row is 1):
% rgb2cam = [0.5566    0.3899    0.0535
%            0.0869    0.7193    0.1938
%            0.0261    0.2539    0.7200]

cam2rgb = inv(rgb2cam); % Invert matrix
% Result:
% cam2rgb = [ 1.9644   -1.1197    0.1553
%            -0.2412    1.6738   -0.4326
%             0.0139   -0.5498    1.5359]

R = lin_rgb(:, :, 1);
G = lin_rgb(:, :, 2);
B = lin_rgb(:, :, 3);

% Left multiply matrix cam2rgb by each RGB tuple (convert from "camera RGB" to "linear sRGB").
sR = cam2rgb(1,1)*R + cam2rgb(1,2)*G + cam2rgb(1,3)*B;
sG = cam2rgb(2,1)*R + cam2rgb(2,2)*G + cam2rgb(2,3)*B;
sB = cam2rgb(3,1)*R + cam2rgb(3,2)*G + cam2rgb(3,3)*B;

lin_srgb = cat(3, sR, sG, sB);
lin_srgb = max(min(lin_srgb, 1), 0); % Clip to range [0, 1]

% Convet from "Linear sRGB" to sRGB (apply gamma)
sRGB = lin2rgb(lin_srgb); % lin2rgb MATLAB functions uses the exact formula [you may approximate it to power of (1/gamma)].

% Show the result, and save to sRGB.png
figure;imshow(sRGB);impixelinfo;title('sRGB'); 
imwrite(im2uint8(sRGB), 'sRGB.png');



% Inverting 3x3 matrix (some help of MATLAB Symbolic Toolbox):
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Assume:
% A = [ a11, a12, a13]
%     [ a21, a22, a23]
%     [ a31, a32, a33]
%
% 1. Compute determinant of A:
% detA = a11*a22*a33 - a11*a23*a32 - a12*a21*a33 + a12*a23*a31 + a13*a21*a32 - a13*a22*a31
%
% 2. Compute the inverse of the matrix A:
% invA = [  (a22*a33 - a23*a32)/detA, -(a12*a33 - a13*a32)/detA,  (a12*a23 - a13*a22)/detA
%          -(a21*a33 - a23*a31)/detA,  (a11*a33 - a13*a31)/detA, -(a11*a23 - a13*a21)/detA
%           (a21*a32 - a22*a31)/detA, -(a11*a32 - a12*a31)/detA,  (a11*a22 - a12*a21)/detA]

结果（缩小）：

RawTherapee 的结果（禁用所有增强功能）：

MATLAB 结果：

蟒蛇结果：

笔记：

  • 由于曝光不足（并且因为我们没有应用任何亮度校正），结果看起来很暗。
以上是如何将白平衡系数应用于sRGB输出的RAW图像的全部内容。
THE END
分享
二维码
< <上一篇
下一篇>>