如何将白平衡系数应用于 RAW 图像以进行 sRGB 输出
How to apply white balance coefficents to RAW image for sRGB output
我想将 RAW 图像数据 (RGGB) 转换为 sRGB 图像。有很多专门的方法可以做到这一点，但为了首先了解基础知识，我已经实现了一些简单的算法，比如通过降低分辨率进行去拜耳化。
我当前的管道是：
- 通过 blacklevel 和 whitelevel
重新缩放 u16 输入数据
- 应用白平衡系数
- 尺寸减小的 Debayer，G 的平均值：g=((g0+g1)/2)
- 计算 D65 光源的伪逆 XYZ_TO_CAM（来自 Adobe DNG）
- 通过 CAM_TO_XYZ
将去拜耳化的 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，如下面 guide.
所述
根据DNG spec：
ColorMatrix1 defines a transformation matrix that converts XYZ values to reference camera native color space values, under the first calibration illuminant.
表示如果校准光源为D65，ColorMatrix将XYZ转换为“相机RGB”。
（按原样转换，不使用任何白平衡缩放系数）。
- 反向 ColorMatrix，从“相机 RGB”转换为 XYZ。
将 XYZ 转换为 sRGB 后，结果为 颜色平衡 sRGB。
结论是ColorMatrix里面包含了白平衡系数（白平衡系数适用D65光源）。
- 将
rgb2cam的行归一化为1，中和了while平衡系数，只保留“颜色校正矩阵”（数学有点复杂）。
- 在不对行进行标准化的情况下，我们将 while 平衡乘数缩放两次：
- 平衡 D65 输入的 ColorMatrix 缩放系数。
- 取自 AsShotNatural 的比例系数平衡场景光源的输入（场景光源接近 D65）。
缩放两次的结果是极度偏色。
跟踪最大值以避免“高光中的洋红色色偏”：
我们假设跟踪“理论最大颜色值”，而不是跟踪输入图像中的实际最大颜色值。
- 取
whitelevel - blacklevel 并按白平衡倍数缩放。
跟踪结果...
指导原则是两种情况下的颜色应该相同：
- 将处理应用于图像的小块，并将这些块放在一起（我们无法跟踪全局最小值和最大值）。
- 将处理应用到整个图像。
我想你必须跟踪缩放的最大值 whitelevel - blacklevel，只有当白平衡乘数小于 1。
当所有乘数都为 1 或以上时，我们可以将结果裁剪为 1.0，而不跟踪最大值。
注：
- 缩小和跟踪最大值可能有优势，但我不知道这个主题。
在我的解决方案中，我们只是乘以上限（1.0 以上），然后裁剪结果。
解决方案基于 Processing RAW Images in MATLAB 指南。
我发布了 MATLAB 实现和 Python 实现（但没有 Rust 实现）。
第一步是使用 dcraw 命令行从 sample.dng 中提取原始拜耳图像：
dcraw -4 -D -T sample.dng
将 tiff 输出重命名为 sample__lin_bayer.tif。
转换过程：
- 通过
blacklevel 和 whitelevel 重新缩放 uint16 输入数据（从所有像素中减去 blacklevel 并缩放 whitelevel - blacklevel）。
- 应用白平衡缩放系数。
缩放系数等于 1./AsShotNatural.
Bayer对齐中的红色像素按红色缩放系数缩放，绿色按绿色缩放，蓝色按蓝色缩放。
假设：最小缩放比例为 1.0，其他比例高于 1.0（我们除以最小缩放比例以确保）。
- 将缩放后的结果裁剪为 [0, 1]（由于
demosaic 实施限制，需要裁剪）。
- 使用 MATLAB
demosaic 函数或 Python 中的 cv2.cvtColor 去马赛克 (Debayer)。
- 计算
rgb2cam矩阵：rgb2cam = ColorMatrix * rgb2xyz.
rgb2xyz 矩阵取自 Bruce Lindbloom site.
- 标准化
rgb2cam 矩阵的行，使每行的总和等于 1（每行除以该行的总和）。
- 通过反转
rgb2cam 计算 cam2rgb 矩阵：cam2rgb = inv(rgb2cam).
cam2rgb是“CCM矩阵”（Color Correction Matrix）。
- 每个 RGB 元组左乘矩阵
cam2rgb（应用颜色校正）。
- 应用伽玛校正（使用 sRGB 标准伽玛）。
- 转换为
uint8并保存为PNG格式（PNG格式用于SO网站发布）。
MATLAB 实现：
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]
Python 实施：
import numpy as np
import cv2
def lin2rgb(im):
""" Convert im from "Linear sRGB" to sRGB - apply Gamma. """
# sRGB standard applies gamma = 2.4, Break Point = 0.00304 (and computed Slope = 12.92)
g = 2.4
bp = 0.00304
inv_g = 1/g
sls = 1 / (g/(bp**(inv_g - 1)) - g*bp + bp)
fs = g*sls / (bp**(inv_g - 1))
co = fs*bp**(inv_g) - sls*bp
srgb = im.copy()
srgb[im <= bp] = sls * im[im <= bp]
srgb[im > bp] = np.power(fs*im[im > bp], inv_g) - co
return srgb
filename = 'sample__lin_bayer.tif' # Output of: dcraw -4 -D -T sample.dng
# Exif information:
blacklevel = 511
whitelevel = 12735
AsShotNeutral = np.array([0.5185, 1, 0.5458])
ColorMatrix = np.array([[ 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 = np.array([[0.4124564, 0.3575761, 0.1804375],
[0.2126729, 0.7151522, 0.0721750],
[0.0193339, 0.1191920, 0.9503041]])
# Read input image (Bayer mosaic alignement, pixeles data type is np.uint16):
raw = cv2.imread(filename, cv2.IMREAD_UNCHANGED)
# "Linearizing":
# There is no LinearizationTable so we only have to subtract the black level.
# Convert from range [blacklevel, whitelevel] to range [0, 1] (convert type to np.float64).
lin_bayer = (raw.astype(np.float64) - blacklevel) / (whitelevel - blacklevel)
lin_bayer = lin_bayer.clip(0, 1)
# The White Balance multipliers are 1./AsShotNeutral
wb_multipliers = 1 / AsShotNeutral
r_scale = wb_multipliers[0] # Assume value is above 1
g_scale = wb_multipliers[1] # Assume value = 1
b_scale = wb_multipliers[2] # Assume value is above 1
# Bayer alignment is RGGB:
# R G
# G B
#
# Apply white balancing to linear Bayer image.
balanced_bayer = lin_bayer.copy()
balanced_bayer[0::2, 0::2] = balanced_bayer[0::2, 0::2]*r_scale # Red (indices [0, 2, 4,... ; 0, 2, 4,... ])
balanced_bayer[0::2, 1::2] = balanced_bayer[0::2, 1::2]*g_scale # Green (indices [0, 2, 4,... ; 1, 3, 5,... ])
balanced_bayer[1::2, 0::2] = balanced_bayer[1::2, 0::2]*g_scale # Green (indices [1, 3, 5,... ; 0, 2, 4,... ])
balanced_bayer[1::2, 1::2] = balanced_bayer[1::2, 1::2]*b_scale # Blue (indices [1, 3, 5,... ; 0, 2, 4,... ])
# Clip to range [0, 1] for avoiding "pinkish highlights" (avoiding "magenta cast" in the highlights).
balanced_bayer = np.minimum(balanced_bayer, 1)
# Demosaicing:
temp = np.round((balanced_bayer*(2**16-1))).astype(np.uint16) # Convert from double to np.uint16, because OpenCV demosaic() function requires a uint8 or uint16 input.
lin_rgb = cv2.cvtColor(temp, cv2.COLOR_BayerBG2RGB).astype(np.float64)/(2**16-1) # Apply Demosaicing and convert back to np.float64 in range [0, 1] (is there a bug in OpenCV Bayer naming?).
# 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 = np.sum(rgb2cam, 1)
# Result:
# rows_sum = [0.4706
# 1.0593
# 0.7470]
# Divide element of every row by the sum of the row:
rgb2cam[0, :] = rgb2cam[0, :] / rows_sum[0] # Divide top row
rgb2cam[1, :] = rgb2cam[1, :] / rows_sum[1] # Divide center row
rgb2cam[2, :] = rgb2cam[2, :] / rows_sum[2] # 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 = np.linalg.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[:, :, 0]
g = lin_rgb[:, :, 1]
b = lin_rgb[:, :, 2]
# Left multiply matrix cam2rgb by each RGB tuple (convert from "camera RGB" to "linear sRGB").
sr = cam2rgb[0, 0]*r + cam2rgb[0, 1]*g + cam2rgb[0, 2]*b
sg = cam2rgb[1, 0]*r + cam2rgb[1, 1]*g + cam2rgb[1, 2]*b
sb = cam2rgb[2, 0]*r + cam2rgb[2, 1]*g + cam2rgb[2, 2]*b
lin_srgb = np.dstack([sr, sg, sb])
lin_srgb = lin_srgb.clip(0, 1) # Clip to range [0, 1]
# Convert 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)].
# Save to sRGB.png
cv2.imwrite('sRGB.png', cv2.cvtColor((sRGB*255).astype(np.uint8), cv2.COLOR_RGB2BGR))
结果（缩小比例）：
RawTherapee的结果（禁用所有增强功能）：
MATLAB 结果：
Python 结果：
注：
- 由于曝光度低（并且因为我们没有应用任何亮度校正），结果看起来很暗。
我想将 RAW 图像数据 (RGGB) 转换为 sRGB 图像。有很多专门的方法可以做到这一点，但为了首先了解基础知识，我已经实现了一些简单的算法，比如通过降低分辨率进行去拜耳化。 我当前的管道是：
- 通过 blacklevel 和 whitelevel 重新缩放 u16 输入数据
- 应用白平衡系数
- 尺寸减小的 Debayer，G 的平均值：g=((g0+g1)/2)
- 计算 D65 光源的伪逆 XYZ_TO_CAM（来自 Adobe DNG）
- 通过 CAM_TO_XYZ 将去拜耳化的 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，如下面 guide.
根据DNG spec：
ColorMatrix1 defines a transformation matrix that converts XYZ values to reference camera native color space values, under the first calibration illuminant.
表示如果校准光源为D65，ColorMatrix将XYZ转换为“相机RGB”。
（按原样转换，不使用任何白平衡缩放系数）。
- 反向 ColorMatrix，从“相机 RGB”转换为 XYZ。
将 XYZ 转换为 sRGB 后，结果为 颜色平衡 sRGB。
结论是ColorMatrix里面包含了白平衡系数（白平衡系数适用D65光源）。 - 将
rgb2cam的行归一化为1，中和了while平衡系数，只保留“颜色校正矩阵”（数学有点复杂）。 - 在不对行进行标准化的情况下，我们将 while 平衡乘数缩放两次：
- 平衡 D65 输入的 ColorMatrix 缩放系数。
- 取自 AsShotNatural 的比例系数平衡场景光源的输入（场景光源接近 D65）。
缩放两次的结果是极度偏色。
跟踪最大值以避免“高光中的洋红色色偏”：
我们假设跟踪“理论最大颜色值”，而不是跟踪输入图像中的实际最大颜色值。
- 取
whitelevel - blacklevel并按白平衡倍数缩放。
跟踪结果...
指导原则是两种情况下的颜色应该相同：
- 将处理应用于图像的小块，并将这些块放在一起（我们无法跟踪全局最小值和最大值）。
- 将处理应用到整个图像。
我想你必须跟踪缩放的最大值 whitelevel - blacklevel，只有当白平衡乘数小于 1。
当所有乘数都为 1 或以上时，我们可以将结果裁剪为 1.0，而不跟踪最大值。
注：
- 缩小和跟踪最大值可能有优势，但我不知道这个主题。
在我的解决方案中，我们只是乘以上限（1.0 以上），然后裁剪结果。
解决方案基于 Processing RAW Images in MATLAB 指南。
我发布了 MATLAB 实现和 Python 实现（但没有 Rust 实现）。
第一步是使用 dcraw 命令行从 sample.dng 中提取原始拜耳图像：
dcraw -4 -D -T sample.dng
将 tiff 输出重命名为 sample__lin_bayer.tif。
转换过程：
- 通过
blacklevel和whitelevel重新缩放uint16输入数据（从所有像素中减去blacklevel并缩放whitelevel - blacklevel）。 - 应用白平衡缩放系数。
缩放系数等于1./AsShotNatural.
Bayer对齐中的红色像素按红色缩放系数缩放，绿色按绿色缩放，蓝色按蓝色缩放。
假设：最小缩放比例为1.0，其他比例高于1.0（我们除以最小缩放比例以确保）。 - 将缩放后的结果裁剪为 [0, 1]（由于
demosaic实施限制，需要裁剪）。 - 使用 MATLAB
demosaic函数或 Python 中的cv2.cvtColor去马赛克 (Debayer)。 - 计算
rgb2cam矩阵：rgb2cam = ColorMatrix * rgb2xyz.
rgb2xyz矩阵取自 Bruce Lindbloom site. - 标准化
rgb2cam矩阵的行，使每行的总和等于1（每行除以该行的总和）。 - 通过反转
rgb2cam计算cam2rgb矩阵：cam2rgb = inv(rgb2cam).
cam2rgb是“CCM矩阵”（Color Correction Matrix）。 - 每个 RGB 元组左乘矩阵
cam2rgb（应用颜色校正）。 - 应用伽玛校正（使用 sRGB 标准伽玛）。
- 转换为
uint8并保存为PNG格式（PNG格式用于SO网站发布）。
MATLAB 实现：
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]
Python 实施：
import numpy as np
import cv2
def lin2rgb(im):
""" Convert im from "Linear sRGB" to sRGB - apply Gamma. """
# sRGB standard applies gamma = 2.4, Break Point = 0.00304 (and computed Slope = 12.92)
g = 2.4
bp = 0.00304
inv_g = 1/g
sls = 1 / (g/(bp**(inv_g - 1)) - g*bp + bp)
fs = g*sls / (bp**(inv_g - 1))
co = fs*bp**(inv_g) - sls*bp
srgb = im.copy()
srgb[im <= bp] = sls * im[im <= bp]
srgb[im > bp] = np.power(fs*im[im > bp], inv_g) - co
return srgb
filename = 'sample__lin_bayer.tif' # Output of: dcraw -4 -D -T sample.dng
# Exif information:
blacklevel = 511
whitelevel = 12735
AsShotNeutral = np.array([0.5185, 1, 0.5458])
ColorMatrix = np.array([[ 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 = np.array([[0.4124564, 0.3575761, 0.1804375],
[0.2126729, 0.7151522, 0.0721750],
[0.0193339, 0.1191920, 0.9503041]])
# Read input image (Bayer mosaic alignement, pixeles data type is np.uint16):
raw = cv2.imread(filename, cv2.IMREAD_UNCHANGED)
# "Linearizing":
# There is no LinearizationTable so we only have to subtract the black level.
# Convert from range [blacklevel, whitelevel] to range [0, 1] (convert type to np.float64).
lin_bayer = (raw.astype(np.float64) - blacklevel) / (whitelevel - blacklevel)
lin_bayer = lin_bayer.clip(0, 1)
# The White Balance multipliers are 1./AsShotNeutral
wb_multipliers = 1 / AsShotNeutral
r_scale = wb_multipliers[0] # Assume value is above 1
g_scale = wb_multipliers[1] # Assume value = 1
b_scale = wb_multipliers[2] # Assume value is above 1
# Bayer alignment is RGGB:
# R G
# G B
#
# Apply white balancing to linear Bayer image.
balanced_bayer = lin_bayer.copy()
balanced_bayer[0::2, 0::2] = balanced_bayer[0::2, 0::2]*r_scale # Red (indices [0, 2, 4,... ; 0, 2, 4,... ])
balanced_bayer[0::2, 1::2] = balanced_bayer[0::2, 1::2]*g_scale # Green (indices [0, 2, 4,... ; 1, 3, 5,... ])
balanced_bayer[1::2, 0::2] = balanced_bayer[1::2, 0::2]*g_scale # Green (indices [1, 3, 5,... ; 0, 2, 4,... ])
balanced_bayer[1::2, 1::2] = balanced_bayer[1::2, 1::2]*b_scale # Blue (indices [1, 3, 5,... ; 0, 2, 4,... ])
# Clip to range [0, 1] for avoiding "pinkish highlights" (avoiding "magenta cast" in the highlights).
balanced_bayer = np.minimum(balanced_bayer, 1)
# Demosaicing:
temp = np.round((balanced_bayer*(2**16-1))).astype(np.uint16) # Convert from double to np.uint16, because OpenCV demosaic() function requires a uint8 or uint16 input.
lin_rgb = cv2.cvtColor(temp, cv2.COLOR_BayerBG2RGB).astype(np.float64)/(2**16-1) # Apply Demosaicing and convert back to np.float64 in range [0, 1] (is there a bug in OpenCV Bayer naming?).
# 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 = np.sum(rgb2cam, 1)
# Result:
# rows_sum = [0.4706
# 1.0593
# 0.7470]
# Divide element of every row by the sum of the row:
rgb2cam[0, :] = rgb2cam[0, :] / rows_sum[0] # Divide top row
rgb2cam[1, :] = rgb2cam[1, :] / rows_sum[1] # Divide center row
rgb2cam[2, :] = rgb2cam[2, :] / rows_sum[2] # 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 = np.linalg.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[:, :, 0]
g = lin_rgb[:, :, 1]
b = lin_rgb[:, :, 2]
# Left multiply matrix cam2rgb by each RGB tuple (convert from "camera RGB" to "linear sRGB").
sr = cam2rgb[0, 0]*r + cam2rgb[0, 1]*g + cam2rgb[0, 2]*b
sg = cam2rgb[1, 0]*r + cam2rgb[1, 1]*g + cam2rgb[1, 2]*b
sb = cam2rgb[2, 0]*r + cam2rgb[2, 1]*g + cam2rgb[2, 2]*b
lin_srgb = np.dstack([sr, sg, sb])
lin_srgb = lin_srgb.clip(0, 1) # Clip to range [0, 1]
# Convert 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)].
# Save to sRGB.png
cv2.imwrite('sRGB.png', cv2.cvtColor((sRGB*255).astype(np.uint8), cv2.COLOR_RGB2BGR))
结果（缩小比例）：
RawTherapee的结果（禁用所有增强功能）：
MATLAB 结果：
Python 结果：
注：
- 由于曝光度低（并且因为我们没有应用任何亮度校正），结果看起来很暗。