isp算法学习CCM（Color Correction Matrix）
1700 Words|Read in about 8 Min|本文总阅读量次
本文是isp算法学习的第十二篇章，主要是讲述CCM。
0开篇介绍
0.1人眼波长反应曲线
人眼之所以能感受到自然界的颜色是因为人眼的视锥细胞在起作用。人眼主要通过三种视锥细胞感受三种不同波长的光从而感受颜色。如图所示是人眼感受不同波长的反应曲线，分别对应三种不同的波长，所以通常用RGB三原色来表示颜色。
0.2sensor波长响应曲线
如图是IMX415C和ICX262AQ两种sensor的感光特性，可以看出来和人眼的感光曲线有很大的不同，而且同样SONY生产的sensor各自的感光曲线的特性也有很大差异，所以如果直接用sensor感光的特性来表示颜色回合人眼有很大的差异，且同一种颜色通过不同sensor感光后得到的数据是不一样的，而对于显示器而言就会表现出不同的颜色，这是我们不希望看到的，我们希望同一种颜色即使使用不同的sensor也能得到相同的RGB数据且和人眼感受的颜色一致或者接近。
IMX415C的sensor的感光特性图
ICX262AQsensor的感光特性图
为了达到上述的要求，我们可以理解为人眼对物体感受的颜色是我们的目标，那么就需要将sensor感光数据经过某种变换达到我们的目标。假设人眼能感受到的颜色种类有m种，那么自然界的颜色就可以表示为一个3Xm的矩阵，同理sensor对自然界的感光也可以得到一个3Xm的矩阵。
那么我们需要做的就是将右侧sensor感光的数据转换到左侧人眼感光的数据上来。（左侧为人眼感光数据矩阵，右侧为sensor感光数据矩阵），其中M的9个元素为3X3的CCM矩阵。
0.3CCM
CCM(Color Correction Matrix)的作用就是通过一个3X3的矩阵使得颜色更接近人眼所感受的颜色。
1颜色学基础
1.1CIE RGB
这里就不得不提到颜色匹配实验了。
在20世纪初，大卫·赖特和约翰·吉尔德独立地与人类观察者进行了实验，在实验中，测试对象试图将红、绿、蓝光源的光谱颜色以10纳米的增量相匹配。这个实验非常类似于我们用滑块来匹配目标颜色，作为其他三种颜色的组合。
该实验原理如下图所示，右侧有两个屏幕，一侧有一个可以改变波长的光源，一侧是固定三种波长（红绿蓝）的光源，然后人眼通过一个角度去看这两个屏幕。f光源改变不同的波长呈现出不同的颜色，然后通过P1，P2，P3三个波长的光源不同的强度的混合，使得人眼感受到两个屏幕的颜色一致，然后记录下此时的三色值，那么该三色值就表示f测波长对应的颜色值。
国际照明委员会（CIE）于1931年提出了一个RGB颜色模型
从下图中可以看到CIE RGB有数据是负数，这个是因为某个波长r(λ)颜色的光通过P1，P2，P3在一侧混合无法得到，那么就需要将某个光源放到f测去才能到达这种效果，那么此时就相当于对改颜色做了减法，那么就出现了负数。
光谱三刺激值
如果待测色光是辐射能量都相同的单色光（光谱色，光波段为380-760nm），我们把得到的三原色光强度值称之为：光谱三刺激值，用 r(λ),g(λ),r(λ)表示。颜色匹配过程可用公式表示如下： $$ C_\lambda=\overline{r}(R)+\overline{g}(G)+\overline{b}(B) $$
光谱三刺激值曲线图
CIE将实验的参考光源CIE R标准化为波长为700.0 nm的单色光，参考光源CIE G标准化为546.1 nm，参考光源CIE B标准化为435.8 nm，并给出了它们之间的相对强度比。这是我们讨论颜色的关键一步——我们最终基于一些物理上可测量的属性。
让我们看看黄色波长570nm。从颜色匹配图中我们可以看出，颜色匹配需要0.16768个单位的CIE R, 0.17087个单位的CIE G和- 0.00135个单位的CIE B。一组三个RGB坐标只是为了三维的呈现。如果我们读取每个波长的颜色匹配值，我们得到下面的图:
Q:为什么光谱曲线没有经过CIE RGB立方体的纯红、绿或蓝端点，
A:它只是委员会对颜色匹配函数进行规范化和缩放的结果。在合理的范围内，光源可以根据需要被构造成强大的，所以什么构成“1.0”总是被任意地定义。
1.2CIE XYZ
因为发现这个上述RGB模型存在缺陷，其中红色响应在435.1~546.1波段出现了负值，给配色造成了极大不便。CIE RGB色彩空间植根于单色光的具体物理特性，我们可以用它来定义任何颜色感觉。然而，委员会努力创造一个衍生空间，将有一些有用的属性，其中两个值得注意:
- 没有光谱颜色的负坐标
- 色度(色调和色彩)与亮度(亮度的视觉感知)的分离
经过一番考虑，我们得到了下面的一组方程:
$$
X = 2.769×R_{CIE} + 1.752×G_{CIE} + 1.130×B_{CIE}\
Y = 1.000×R_{CIE} + 4.591×G_{CIE} + 0.061×B_{CIE}\
Z = 0.000×R_{CIE} + 0.057×G_{CIE} + 5.594×B_{CIE}
$$
这些因素可能看起来是任意的，在某种意义上它们是，但这只是另一个从一个空间到另一个空间的值映射，类似于我们用RGB到R’G’B’变换所做的。选择Y分量作为亮度分量，而X和Z定义色度。
国际照明委员会（CIE）于1931年首次用纯数学的方法定义了一种颜色空间模型，即CIE 1931 XYZ 颜色空间，它覆盖了人眼所能感知的全部颜色，而且不依赖于任何特定的物理实现。当一个颜色以XYZ 形式表示时，在任何支持XYZ标准的显示器上都会渲染出一致的颜色，与显示器的具体物理特性无关。
你可以看到XYZ空间是如何围绕光谱曲线设计的，以确保它符合正八分区。它的小尺寸并不是一个真正的问题，通过确保感知到的最亮波长555nm的Y值为1，缩放使事情更方便使用。作为最后一步，我们可以单独显示XYZ空间和光谱颜色:
由于CIE XYZ空间源自CIE RGB，它也基于可测量的物理量。XYZ色彩空间是用于所有基于矩阵的RGB显示色彩空间转换的基本色彩空间。任何两个色彩空间之间的映射都是通过一个通用的CIE XYZ中间体完成的。这简化了颜色变换，因为我们不需要知道如何在任意两个颜色空间之间转换颜色，我们只需要能够将它们都转换到XYZ空间。
CIE 1931 XYZ系统光谱三刺激值曲线图
1.3CIE xy色度图
三维图表玩起来很有趣，但它们通常不是特别实用——2D图表通常更容易处理和推理。由于XYZ空间的Y分量没有任何色彩和色调，所以我们只剩下两个影响颜色色度的分量。我们可以执行三种旨在降低空间维数的操作（这里归一化处理）:
$$
x = X / (X + Y + Z)\
y = Y / (X + Y + Z)\
z = Z / (X + Y + Z)
$$
小写和大写的区别在这里很重要——x和X不一样。这些看似任意的方程有一个简单的视觉解释——它是在XYZ坐标[1.0,0.0,0.0]和[0.0,1.0,0.0]和[0.0,0.0,1.0]之间的三角形上的投影。
注意到x, y, z加起来等于1，所以我们可以去掉最后一个分量，因为它是冗余的——我们总是可以通过从1中减去x和y来重新生成它。拒绝z等于在xy平面上的平面投影:
如果我们对光谱颜色的每个组合重复这一步骤，我们最终可以呈现出被称为CIE xy色度图的2D图。你以前可能见过这个马蹄形:
补充说明
1.色觉异常的色度学表现
大部分色觉异常者仅有一种锥细胞无法正常工作，而剩余的锥细胞依然能提供“彩色”的视觉信息。因此，色觉异常者眼中的世界并非只有黑白两色，只不过他们的辨色能力比普通人更弱，他们所“看见”的世界也和普通人的不同。
由于辨色能力的缺陷，很多普通人看来明显不同的颜色，色觉异常者却认为几乎一致，这种现象称为颜色混同（color confusion）**。**研究者在CIE1931色品图上绘制出了混同线，来标识色觉异常者难以辨别的颜色。从图中可以看出，红、绿色盲都难以分辨红色和绿色，红色盲还无法分辨蓝色和紫色，可能将青色视为灰色，而绿色盲容易混淆青色和紫色。
2.等色温线
色温是由英国物理学家洛德·开尔文于19世纪末提出的，其通过对比光源的颜色和理想黑体来确定。当某光源所呈现的颜色与黑体在某一温度下的颜色相同时，则此时黑体的温度就称为该光源的颜色温度，简称色温。由于某些光源的光谱分布与黑体相差较远，它们在某一温度时的相对光谱功率分布所决定的色品坐标不一定准确地落在色品图的黑体温度轨迹上，这时不能用色温来描述其颜色，而采用相关色温来表征这类光源的光色特性。相关色温定义为当光源的颜色与黑体在某一温度下的颜色最接近时的黑体温度。一般，日常生活中所提到的色温是指相关色温。（不同色温下的黑体轨迹和等色温线）
1.4sRGB色域
sRGB代表了标准的红、绿、蓝，即CRT显示器、LCD面板、投影机、打印机以及其他设备中色彩再现所使用的三个基本色素。sRGB的色彩空间基于独立的色彩坐标，可以使色彩在不同的设备使用传输中对应于同一的色彩坐标体系，而不受这些设备各自具有的不同色彩坐标的影响。 投影显示系统的sRGB是微软公司与精工爱普生公司、三菱公司合作开发的，目的是建立一个可以满足计算机和投影显示需求的色彩管理标准，使得显示设备无须经过特别的色彩信息分析，就可以正确地表现出图象文件。sRGB消除了不同显示系统在色彩还原上的原有差异。不同显示设备间的RGB色彩，自然会发生一些变化，因而经过不同的显示设备后就无法正确地再现色彩。
因为现实器也不可能现实自然界所有的颜色，所以就在xy坐标系中选取一个三角形范围作为现实器的显示色域。同时通过定义不同的区域还有Adobe RGB色域空间和其他。
常用的设备无关颜色空间是CIE XYZ颜色空间。而生活中常用的设备无关颜色空间还有sRGB、AdobeRGB和ProPhoto RGB颜色空间。通常浏览网页时，我们看到的图像使用的就是sRGB颜色空间。目前ISP中CCM的目标通常是以sRGB为常用的设备无关颜色空间。
颜色管理流程示例
之所以图像在跨设备进行传输时产生了变化，主要的原因是两台设备采用设备相关的RGB传输，而输入相同的RGB值在不同设备中却对应不同的实际颜色所致。
在设备相关颜色空间与设备无关颜色空间之间尽可能准确地建立数学映射关系就是颜色特征化。简单来说，颜色特征化能通过设备的数字驱动值RGB精确预测出设备显现的实际颜色。因此，设备的颜色特征化是颜色管理的基础步骤和关键环节，它可以精确控制颜色信息在图像设备之间的复制与传递。
1.5L*A*B和LUV
生活中我们经常使用设备相关的RGB颜色空间进行颜色的处理和表征。但是，发现RGB间距相同的颜色组合所产生的颜色感知差异是不同的，这是因为在该颜色空间中，各色区的颜色感知差异容限大小不等。虽然在驱动值RGB间距相等的颜色对，但是人眼感知还是有差异，左边会觉得差异更大。
1.5.1颜色空间的非均匀性
上述问题同样存在于国际照明委员会（CIE）推荐的CIE1931和CIE1964标准色度系统中，这也使得在工业生产中难以使用标准色度系统对颜色进行精确的比较和评价。
1942年麦克亚当对CIE-xy色品图的非均匀性进行了详细的实验研究，他选定了25个颜色中心，在5~9个方向上用加混色进行多次颜色匹配测试。在实验中，他采用2°观察视场并保持颜色亮度不变。实验表明，即使对某一确定的方向，每次实验的结果也是不同的，并在视觉可识别的界限内变动。麦克亚当将不同方向上该变动的标准偏差在色品图上描点，得到下图的麦克亚当椭圆 (图中椭圆放大了10 倍)。
通过麦克亚当椭圆在CIE-xy色品图上的分布，我们可以看出CIE-xy色品图是非常不均匀的。在色品图左下方的蓝紫色区，色品坐标稍有变化就能识别，而对于上方的绿色区，色品坐标约变化前者的10倍才能分辨。理想的均匀色品图应满足图上任何位置的麦克亚当椭圆都是半径相等的圆。
1.5.2均匀颜色空间
国际照明委员会根据麦克亚当的工作对CIE-xy色品图进行适当的变换，对色品图的非均匀性进行了一定程度的改善。同时，在结合明度与亮度的非线性关系，对亮度因数进行修正后，于1976年推荐了CIE1976 L*u*v颜色空间（也称CIELUV颜色空间）和CIE1976 L*a*b颜色空间（也称CIELAB颜色空间），分别用于电视工业等加混色和表面色料工业等减混色的表示和评价。
1.5.2.1CIE1976 L*u*v*颜色空间
在CIE1976 L*u*v颜色空间中，L*为明度，u*、v*表示色品坐标，其计算公式为
$$
\begin{cases}
L^=116(Y/Y_n)^{1/3}-16,&&Y/Y_n>(24/116)^3\
u^=13L(u'-u'_n)\
v^=13L(v'-v'_n)
\end{cases}
$$
式中，u'、v’和u'n、v'n分别是颜色样品和CIE标准照明体的CIE1976 UCS色品坐标；X、Y、Z为颜色样品的三刺激值，Xn、Yn、Zn为CIE标准照明体照射在完全漫反射体上，然后反射到观察者眼中的三刺激值，其中Yn=100。
如果 u* 和 v* 都是正值，则 huv 在 0° 和 90° 之间；
如果 v* 是正值，而 u* 是负值，则 huv 在 90° 和 180° 之间；
如果 v和 u 都是负值，则 huv在 180° 和 270° 之间；
如果 v* 是负值，而 u* 是正值，则 huv 在 270° 和 360° 之间。
1.5.2.2CIE1976 L*a*b*颜色空间
CIE L*a*b*又称CIELAB颜色空间，是目前应用最广的颜色空间之一，它的主要优点是颜色之间的距离与人的知觉更好地符合线性关系，尤其是描述较暗的颜色准确度更高。轻微的不足是描述黄色区域时线性关系会有一定变化，即颜色宽容量的圆圈在黄色附近直径会有变化。
CIELAB 颜色空间在知觉均匀性方面与CIE 1976 UCS类似，两种标准都不是十全十美的，但都取得了广泛的应用。关于这两个颜色空间哪一个会更均匀一些，学术界还没有形成明确的结论。目前似乎是CIELAB的应用会更多一些。
在CIE1976 L*a*b颜色空间中，L*为明度，a*、b*表示色品坐标，其计算公式为
$$
\begin{cases}
L^=116(Y/Y_n)^{1/3}-16,&&Y/Y_n>(24/116)^3\
a^=500[f(X/X_n)-f(Y/Y_n)]\
b^=200[f(Y/Y_n)-f(Z/Z_n)]
\end{cases}
$$
式中，X、Y、Z和Xn、Yn、Zn的含义与CIELUV颜色空间的对应参数相同。
如果 a* 和 b* 均为正值，hab 位于 0° 和 90° 之间；
如果 b* 为正值，a* 为负值，hab 位于 90° 和 180° 之间；
如果 b* 和 a* 均为负值，hab 位于 180° 和 270° 之间；
如果 b* 为负值，a* 为正值，hab 位于 270° 和 360° 之间。
此外，在CIELUV和CIELAB颜色空间中还各有一组与心理量近似对应的感知属性，即明度（L*，如上）以及彩度和色调角，以下是具体计算公式，并以下标“uv”和“ab”来区分。
1.5.2.3彩度
$$
\begin{cases}
C_{uv}^*=(u^{*2}+v^{*2})^{1/2}\
C_{ab}^*=(a^{*2}+b^{*2})^{1/2}
\end{cases}
$$
1.5.2.4色调角
$$
\begin{cases}
h_{uv}^*=arctan(v^*/u^*)\
h_{ab}^*=arctan(b^*/a^*)
\end{cases}
$$
1.5.2.5饱和度
另外，在CIELUV颜色空间中还有一个心理相关量即饱和度，其计算公式为 $$ S_{uv}=13[(u'-u'_n)^2+(v'-v'_n)^2]^{1/2}=C_{uv}^*/L^* $$ CIELUV和CIELAB颜色空间由CIE XYZ系统的非线性变换推导而来，以其转换公式简单、颜色空间相对均匀等优点，在色彩行业得到了广泛的应用
1.5.3计算色差
色差公式旨在提供成对颜色样品之间感知色差的定量表示，理想的色差公式应基于真正视觉感知均匀的颜色空间，其预测的色差应与目视判断具有良好的一致性。颜色科学工作者对色差评价进行了大量研究，并提出了多种色差公式，其中应用广泛的典型色差评价模型包括CIELAB、CMC(l:c)、CIE94以及CIEDE2000色差公式。
1.5.3.1CIELUV和CIELAB色差公式
CIELUV和CIELAB色差公式以被比较颜色点在对应颜色空间中的欧氏距离表示色差，即
$$
\begin{cases}
\Delta{E_{uv}^*}=(\Delta{L^*}^2+\Delta{u^*}^2+\Delta{v^*}^2)^{1/2}\
\Delta{E_{ab}^*}=(\Delta{L^*}^2+\Delta{a^*}^2+\Delta{b^*}^2)^{1/2}
\end{cases}
$$
然而，随着更多心理物理学实验的开展，人们逐渐发现上述两个色差方程并不能准确量化中小幅度的色差，并且在u*v*、a*b*平面中麦克亚当椭圆的形状和分布也不均匀规整。(色差相同的颜色对呈现出不一致的感知差异，很明显左边的色差更大)
1.5.3.2CMC(l:c)色差公式
CMC(l:c)色差公式克服了当时其它色差公式在深色及中性色范围计算结果与视觉感知偏差较大的现象，具有比CIELAB公式更好的视觉一致性，并于1995年成为纺织行业的国际标准。CMC(l:c)色差模型的具体计算公式如下：
$$
\Delta{E_{CMC(l:c)}}={[\Delta{L_{ab}^*}/(lS_L)]^2+[\Delta{C_{ab}^*}/(cS_C)]^2+[\Delta{H_{ab}^*}/S_H]^2}^{1/2}
$$
式中，△L*ab、△C*ab和△H*ab分别对应于CIELAB颜色空间的明度差、彩度差和色调差；、和定义容差椭球体的半轴长度，分别对应明度、彩度和色调方向；明度权重因子*l*和彩度权重因子*c*分别用来调整明度和彩度对总色差的影响程度，所以在不同的应用场合应取其不同的比值。大量实验表明，对色差的可接受性评价时，推荐采用l:c=2:1，如在纺织业界对产品的质量控制大多采用CMC(2:1)公式；而对色差的可察觉性评价时，推荐采用l:c=1:1，如对数字系统的色度校正以及涂料或塑料等行业一般采用CMC(1:1)公式。
1.5.3.3CIE94色差公式
国际照明委员会于1994年推荐了CIE94色差公式。它具有与 CMC(l:c)色差公式相似的结构，但具有更简单的加权函数。该公式探讨了色差距离的均匀性和参量对知觉色差的影响，并对参照条件进行了严格的说明。然而，在实际观察条件与其不符时，需要调整明度、彩度和色调分色差在总色差中的相对权重。CIE94色差模型的具体计算公式如下：
$$
\Delta{E_{94}^*}={[\Delta{L^*}/(K_LS_L)]^2+[\Delta{C_{ab}^*/(K_CS_C)}]^2+[\Delta{H_{ab}^*/(K_HS_H)}]^2}^{1/2}
$$
式中，参数因子KL、KC和KH用以调整明度、彩度和色调分色差的相对权重，SL、SC和SH的意义与CMC(*l*:*c*)色差公式相似。对于除纺织业以外的一般应用，建议所有参数因子的比值均为1；对于纺织行业，建议采用2:1:1的比值。
1.5.3.4CIEDE2000色差公式
为进一步改善工业色差评价的视觉一致性，国际照明委员会成立了TC1-47技术委员会，对已有色差模型和视觉评价数据进行大量分析，并于2000年推荐了新的色差评价公式，即
$$
\Delta{E_{00}}={[\Delta{L'}/(k_LS_L)]^2+[\Delta{C'_{ab}}/(k_CS_C)]^2+[\Delta{H'_{ab}}/(k_HS_H)]^2\
+R_T[\Delta{C'_{ab}}/(k_CS_C)]}^{1/2}
$$
式中，为RT为旋转项，其余因子意义与上述公式相似。CIEDE2000色差公式不仅包括明度、彩度和色调加权函数，还包括彩度差和色调差的交互项，提高了蓝色区域的色差预测性能，也是**CIE目前推荐的最新色差评价公式**。
1.5.3.4不同色差公式的计算结果
1.6 24色卡
ColorChecker Classic色卡由4行6列共计24种不同的色块组成，这些色块涵盖了经科学配制的自然色、原色以及中性色，可在不同光源下与其代表的物体色相匹配，以实现颜色的准确再现。其中，第一行代表了常见的物体色（natural colors），例如深肤色、浅肤色、树叶、天空等；第二行则包含了自然界除上述典型物体色之外的其它颜色（miscellaneous colors）；第三行的色块为饱和度较高的三基色或其组合色（primary and secondary colors），通常用于测试、校准设备的色域边界；最后一行均为中性色块（grayscale colors），是白平衡校正、灰阶校正的重要参考依据，如下图所示。
下表列出了ColorChecker Classic各色块的色名以及所对应的sRGB颜色外貌与参考数值，其中序号20至23中性色块色名中的数字表示其孟塞尔明度值。以第22号色块（neutral 5）为例，其孟塞尔标号为N5，表面反射率约为18.4%。
24个色块在CIE1931 xy色品图上的分布情况，可见其色域基本覆盖了sRGB，适用于绝大多数颜色管理场景。
另外，航空航天中视频画面传输也会被用到，红色区域所示即为24色卡
2算法与案例精选
2.1关于CCM
2.1.1卢瑟条件
卢瑟条件（Luther condition）大意是说三个线性无关的原色，经过混合能够表示任意一种颜色。
2.1.2What is CCM & Why CCM
由于相机的颜色匹配特性通常不满足卢瑟条件（也就是说sensor的RGB响应通常不是线性无关的），即相机的颜色匹配特性与CIE标准观察者之间并不存在线性关系。因此，我们需要某种方法来校正相机的特性，使其接近标准观察者。在实际的ISP处理中，这种方法通常以3x3矩阵的形式出现，称为色彩校正矩阵（colour correction matrix）
CCM的作用主要有两个：
- 完成了sensor_RGB色彩空间到设备无关的颜色空间sRGB色彩空间的转换
- 使得相机的颜色匹配特性满足卢瑟条件，具体效果如下图所示
2.1.3CCM公式的物理意义
CCM公式的物理意义是从一种颜色种减除另外两种颜色的成分，以增加该颜色的饱和度，使变换的结果接近人的视觉感受，或者更符合人的主观审美。由于输入颜色可以有上千万种组合，而CCM参数却只有9个，所以CCM实际上只能优先保证几个最重要的颜色在人看来是“正确”的，而不可能面面俱到地保证所有颜色在所有条件下都是最优的。
2.1.4CCM基本约束
CCM公式的一个基本约束就是不能破坏白平衡，即对于任何R=G=B的输入，必须保证输出满足R'=G'=B'。正式由于这个原因，颜色校正操作只能放在白平衡调整之后执行。通常归一化之后，每一行都会是固定值，比如每一行和为1。
实践中通常使用X-Rite 24色卡上的18个彩色色块为标准计算校正系数，基本原理是用摄像机拍摄色卡，提取18个色块的平均颜色（Rn, Gn, Bn），n=1..18 构成输入矩阵
$$
S=
\begin{bmatrix}
R_1 & R_2 & R_3 &…&R_n\
G_1 & G_2 & G_3 &…&G_n\
B_1 & B_2 & B_3 &…&B_n
\end{bmatrix}
$$
用24色卡上的18个彩色色块的标准sRGB值构成目标矩阵
$$
S'=
\begin{bmatrix}
R_1' & R_2' & R_3' &…&R_n'\
G_1' & G_2' & G_3' &…&G_n'\
B_1' & B_2' & B_3' &…&B_n'
\end{bmatrix}
$$
则有关于CCM的矩阵方程
$$
S'=M·S\
\begin{bmatrix}
R_1' & R_2' & R_3' &…&R_n'\
G_1' & G_2' & G_3' &…&G_n'\
B_1' & B_2' & B_3' &…&B_n'
\end{bmatrix}=
\begin{bmatrix}
a_{11} & a_{12} & a_{13}\
a_{21} & a_{22} & a_{23}\
a_{31} & a_{32} & a_{33}
\end{bmatrix}
\begin{bmatrix}
R_1 & R_2 & R_3 &…&R_n\
G_1 & G_2 & G_3 &…&G_n\
B_1 & B_2 & B_3 &…&B_n
\end{bmatrix}\
其中\begin{bmatrix}
a_{11} & a_{12} & a_{13}\
a_{21} & a_{22} & a_{23}\
a_{31} & a_{32} & a_{33}
\end{bmatrix}=
\begin{bmatrix}
r_1 & r_2 & r_3\
g_1 & g_2 & g_3\
b_1 & b_2 & b_3
\end{bmatrix}=
\begin{bmatrix}
RR & GR & BR\
RG & GG & BG\
RB & GB & BB
\end{bmatrix}
$$
求解M的过程为
$$
S'·S^T=M·S·S^T\
M=S'·S^T(S·S^T)^{-1}
$$
上述过程需要使用某种色卡，因此称为patch-based方法。
2.1.5CCM常见的调试方法
1）将参数在L*a*b* 色域中标注出来:
2）重点关注RGB三原色的色相，往哪边偏就把对应的数值减小，把相邻的数值加大
例如，R向M方向偏，就需要把①对应的数值减小，把②对应的数值加大。
3）用主色分量上的值补偿之前的修改，使每行之和恒等于1
4）饱和度
单一颜色的饱和度可以将主色分量所在列的另外两个分量减小或加大相同的比例。
例如，要加大R的饱和度，需要将g1，b1减小相同的比例；
要减小R的饱和度，需要将g1，b1加大相同的比例。
5）调试小技巧
| 问题现象 | 调试方向 | 备注 |
|---|---|---|
| 红色偏深 | 减rg和减rb | 红色比例增多，相应的要增加gg，bb，则绿色和蓝色饱和度会提高 |
| 绿色更嫩 | 减gb和减gr | gb越小，绿色越黄；gr越小，绿色越青 |
| 肤色偏黄 | 加gb | 蓝和黄相对时，蓝色会变淡 |
| 肤色变红 | 减br或加bg | br越小，蓝色越青 |
| 肤色偏红 | 加br | br越大，蓝色越紫 |
| 红色偏桔 | 减gr或加rb | rg越小，红色越粉 |
| 红色偏粉 | 加rg或减rb | rg越大，红色越桔，红色系饱和度回收 |
| 黄色偏红 | 加br或减bg | br越大，黄色越绿 |
| 黄色偏绿 | 减br或加bg | br越小，黄色越红 |
| 绿色偏青 | 加gr | |
| 绿色偏黄 | 加gb | 蓝黄相对，蓝色变淡 |
| 蓝色变淡 | 减gb | gb越小，绿色越黄，且绿色饱和度下降，反之增加 |
| 蓝色偏青 | 加br或减bg | br越大，蓝色越紫 |
| 蓝色偏紫 | 减br或加bg | br越小，蓝色越青；加bg，蓝色和绿色饱和度都下降，但整体饱和度增加 |
2.2关于CCM算法
这里叫颜色校正算法其实主要是针对转换过程而言，因为CCM只是转换算法的一种。上面我们已经知道我们整个过程其实就是将一个矩阵转换为另外一个矩阵，那么我们首先能想到的就是LUT，就是将两个矩阵表示的颜色都通过一张表来表示对应的关系，那么进来一个颜色就可以通过查表快速得到想要的颜色。
2.2.1 3D-LUT
之前使用过2D的LUT（Look up table），现在用到的是三维空间的查找法。左边与设备无关，RGB通常是三维的空间，三个值都是0-255；和sensor相关的就是右边的颜色空间，通过3D-LUT来获取对应的颜值值。如果新空间到老的空间有对应的点就直接复制，没有的话通过插值的方式来获取。这里会使用3D插值，四个点组成一个四棱锥，获取四棱锥里面四个点距离通过比例进行三维插值。
这种方式用的比较少。
2.2多项式拟合
实际中大多数厂商用得都是这种，求的是CCM矩阵。多项式拟合有分成好多种，这里只介绍典型的。
2.2.1算法
imatest官网有对应的算法，具体网址可以点击这里
1function ccm_test
2% Test program for applying Color Correction matrix.
3% Output for display gamma = 2.2 -- sRGB and Adobe RGB color spaces.
4
5% Read image file
6disp('Enter the file for the Color Correction Matrix (CCM)');
7imagefile = input('','s');
8
9gamma = input('Enter encoding gamma (typically 0.5 to 1) ');
10
11imageRGB = imread(imagefile);
12imtype = class(imageRGB); % Image type: typically uint8 or uint16.
13imageRGB = double(imageRGB)/double(intmax(imtype)); % Normalize to maximum for image type.
14figure; image(imageRGB);
15title(['Original image: gamma = ' num2str(gamma)]); % Works without Image Processing Toolbox.
16
17linearRGB = imageRGB.^(1/gamma); % Linearize the image (apply inverse of encoding gamma).
18
19% The CCM has 3 lines, which makes it tricky to enter.
20disp('(Calculate a 3x3 CCM in Color/Tone Interactive, then press File, ''Copy Color Matrix to Clipboard'').');
21disp('Press ''enter'' when the the CCM is in the clipboard');
22inp = input('','s'); % pause fails after copy to clipboard. inp is not used.
23
24ccm = clipboard('paste');
25ccm = strrep(ccm,10,';'); % Bring back semicolons.
26ccmnum = str2num(ccm), % str2double doesn't work.
27% Change to 2D to apply matrix; then change back.
28[my, mx, mc] = size(linearRGB); % rows, columns, colors (3)
29linearRGB = reshape(linearRGB,my*mx,mc);
30correctedRGB = linearRGB*ccmnum;
31correctedRGB = min(correctedRGB,1); correctedRGB = max(correctedRGB,0); % Place limits on output.
32correctedRGB = correctedRGB.^(1/2.2); % Apply gamma for sRGB, Adobe RGB color space.
33% Deal with saturated pixels. Not perfect, but this is what cameras do. Related to "purple fringing".
34correctedRGB(linearRGB==1) = 1; % Don't change saturated pixels. (We don't know HOW saturated.)
35correctedRGB = reshape(correctedRGB, my, mx, mc);
36
37figure; image(correctedRGB); title('Corrected image');
2.2.2其他形式
多项式拟合还存在一种根多项式和多次项式
2.2.3多项式拟合原理
其中我们sensorRGB与sRGB之间有下面的关系
- sensorRGB到sensor XYZ颜色空间，需要矩阵M1
- sensor XYZ颜色空间到sensor
L*a*b*颜色空间，需要矩阵M2 - sensor
L*a*b*颜色空间到目标的L*a*b*颜色空间，需要矩阵M3 - 目标的
L*a*b*颜色空间到目标的XYZ颜色空间，需要矩阵M4 - 目标的XYZ颜色空间到sRGB颜色空间，需要矩阵M5
上述M2，M4，M5都是已知的，通过sensor L*a*b*颜色空间到目标的 L*a*b*颜色空间，可以求出M3，再通过sensorRGB到sensor XYZ颜色空间，需要矩阵M1，最终M1*M2*M3*M4*M5的综合为所需的CCM矩阵。
3.3神经网络拟合
当然实际中没法把所有的颜色的数据都得到，那么就通过采样的到部分坐标的数据形成一个三维表，其他不在采样点位置的数据就可以通过插值的方式求出来。
除了上述的LUT的方式还有神经网络的方式，因为其实输入就是一个矩阵，输出也是一个矩阵，那么中间通过网络连接，然后通过数据训练也能得到一个转换的网络参数。
3算法代码实现
这里参考QiuJueqin的CCM实现，添加CCM算法，可框选24色卡,，具体点击这里
3.1demo1
这个demo1主要是140个色块的CCM校正
1clear; close all; clc;
2
3%% data preparation
4
5% load DSG color checker's spectral reflectance data
6% 加载光谱反射数据
7load('spectral_reflectance_data.mat');
8spectra = spectral_reflectance_data.XRite_DSG;
9
10% calculate the XYZ values for the color checker under D65
11%这里计算出来的是目标的XYZ颜色空间
12XYZ = spectra2colors(spectra, 400:5:700,...
13 'spd', 'D65');
14
15% split the left half as training samples and the right half as validation
16%将左半部分作为训练样本，右半部分作为验证样本
17XYZ = reshape(XYZ, 14, 10, 3);
18XYZ_train = reshape(XYZ(1:7, :, :), [], 3);
19XYZ_val = reshape(XYZ(8:14, :, :), [], 3);
20
21% visualize the target colors
22figureFullScreen('color', 'w');
23ax1 = subplot(1,2,1);
24%把点变成像色卡一样的块状填充到画面上，形成拟合出来的色卡
25colors2checker(xyz2rgb(XYZ_train),...
26 'layout', [10, 7],...
27 'squaresize', 100,...
28 'parent', ax1);
29title('Training samples (target sRGB)');
30
31ax2 = subplot(1,2,2);
32colors2checker(xyz2rgb(XYZ_val),...
33 'layout', [10, 7],...
34 'squaresize', 100,...
35 'parent', ax2);
36title('Validation samples (target sRGB)');
37
38% read color checker image and extract color responses
39% the darkness level and the spatial non-uniformity has been corrected for
40% this image
41%真实sensor显示的RGB
42dsg_img = imread('dsg_Nikon_D3x.tiff');
43
44% you can also manually select roi for the color checker
45%选择的区域，也可以手动选择ROI区域
46roi = [519,310,65,65;669,309,65,65;819,308,65,65;969,307,65,65;1120,306,65,65;1270,306,65,65;1421,305,65,65;1572,304,65,65;1722,303,65,65;1873,302,65,65;2025,301,65,65;2176,300,65,65;2327,299,65,65;2479,299,65,65;519,460,65,65;669,459,65,65;819,458,65,65;970,457,65,65;1120,456,65,65;1271,456,65,65;1421,455,65,65;1572,454,65,65;1723,453,65,65;1874,452,65,65;2025,452,65,65;2176,451,65,65;2328,450,65,65;2479,449,65,65;519,650,65,65;669,609,65,65;820,608,65,65;970,607,65,65;1120,607,65,65;1271,606,65,65;1422,605,65,65;1572,604,65,65;1723,604,65,65;1874,603,65,65;2026,602,65,65;2177,601,65,65;2328,601,65,65;2480,600,65,65;520,760,65,65;670,759,65,65;820,758,65,65;970,758,65,65;1121,757,65,65;1271,756,65,65;1422,755,65,65;1573,755,65,65;1724,754,65,65;1875,753,65,65;2026,753,65,65;2177,752,65,65;2329,751,65,65;2480,751,65,65;520,910,65,65;670,909,65,65;820,908,65,65;971,908,65,65;1121,907,65,65;1272,907,65,65;1422,906,65,65;1573,905,65,65;1724,905,65,65;1875,904,65,65;2027,903,65,65;2178,903,65,65;2329,902,65,65;2481,901,65,65;520,1060,65,65;670,1059,65,65;820,1059,65,65;971,1058,65,65;1121,1057,65,65;1272,1057,65,65;1423,1056,65,65;1574,1056,65,65;1725,1055,65,65;1876,1055,65,65;2027,1054,65,65;2179,1053,65,65;2330,1053,65,65;2482,1052,65,65;520,1210,65,65;670,1209,65,65;821,1209,65,65;971,1208,65,65;1122,1208,65,65;1272,1207,65,65;1423,1207,65,65;1574,1206,65,65;1725,1206,65,65;1876,1205,65,65;2028,1205,65,65;2179,1204,65,65;2331,1204,65,65;2482,1203,65,65;520,1360,65,65;671,1360,65,65;821,1359,65,65;972,1359,65,65;1122,1358,65,65;1273,1358,65,65;1424,1357,65,65;1575,1357,65,65;1726,1357,65,65;1877,1356,65,65;2028,1356,65,65;2180,1355,65,65;2331,1355,65,65;2483,1354,65,65;521,1510,65,65;671,1510,65,65;821,1510,65,65;972,1509,65,65;1122,1509,65,65;1273,1508,65,65;1424,1508,65,65;1575,1508,65,65;1726,1507,65,65;1877,1507,65,65;2029,1507,65,65;2180,1506,65,65;2332,1506,65,65;2484,1505,65,65;521,1661,65,65;671,1660,65,65;822,1660,65,65;972,1660,65,65;1123,1659,65,65;1274,1659,65,65;1424,1659,65,65;1576,1658,65,65;1727,1658,65,65;1878,1658,65,65;2029,1657,65,65;2181,1657,65,65;2332,1657,65,65;2484,1656,65,65];
47%计算选择的140个色块区域平均的RGB值
48RGB = checker2colors(dsg_img, [10, 14],...
49 'roi', roi,...
50 'scale', 4); % scaling only for better visualization
51
52% split the left half as training samples and the right half as validation
53%跟上述类似，一半训练，一半验证
54RGB = reshape(RGB, 14, 10, 3);
55RGB_train = reshape(RGB(1:7, :, :), [], 3);
56RGB_val = reshape(RGB(8:14, :, :), [], 3);
57
58% visualize the camera colors
59figureFullScreen('color', 'w');
60ax1 = subplot(1,2,1);
61
62colors2checker(RGB_train,...
63 'layout', [10, 7],...
64 'squaresize', 100,...
65 'parent', ax1);
66title('Training samples (camera RGB, w/o gamma)');
67
68ax2 = subplot(1,2,2);
69colors2checker(RGB_val,...
70 'layout', [10, 7],...
71 'squaresize', 100,...
72 'parent', ax2);
73title('Validation samples (camera RGB, w/o gamma)');
74
75clearvars -except XYZ_train XYZ_val RGB_train RGB_val
76
77%% color correction
78
79% training
80%可以定义3x3模型，也可以定义其他模型，这里是6x3模型
81model = 'root6x3'; % color correction model
82%根据光谱特性，sRGB转换到了XYZ颜色空间
83targetcolorspace = 'XYZ';
84%这里用到的计算方式是CIE △E2000的方式求解
85% OPTIONAL PARAMETERS:
86% loss: loss function to be minimized.
87% 'mse' | 'ciede00' (default) | 'ciede94' |
88% 'ciedelab' | 'cmcde'
89% model: color correction model, based on which the camera
90% responses will be expanded.
91% 'linear3x3' (default) | 'root6x3' | 'root13x3' |
92% 'poly4x3' | 'poly6x3' | 'poly7x3' | 'poly9x3'
93metrics = {'mse', 'ciede00', 'ciedelab'}; % only for evaluation
94
95% default loss function: ciede00
96% if RGB_train has been white-balanced, consider using 'preservewhite' and
97% 'whitepoint' parameters.
98[matrix, scale, XYZ_train_pred, errs_train] = ccmtrain(RGB_train,...
99 XYZ_train,...
100 'model', model,...
101 'bias', true,...
102 'targetcolorspace', targetcolorspace,...
103 'metric', metrics);
104% validation
105[XYZ_val_pred, errs_val] = ccmvalidate(RGB_val,...
106 XYZ_val,...
107 model,...
108 matrix,...
109 scale,...
110 'targetcolorspace', targetcolorspace,...
111 'metric', metrics);
112% visualize final results
113figureFullScreen('color', 'w');
114ax1 = subplot(1,2,1);
115colors2checker({xyz2rgb(XYZ_train), xyz2rgb(XYZ_train_pred)},...
116 'layout', [10, 7],...
117 'squaresize', 100,...
118 'legend', {'Ground-truth', 'Predicted'},...
119 'parent', ax1);
120title('Training samples comparison');
121
122ax2 = subplot(1,2,2);
123colors2checker({xyz2rgb(XYZ_val), xyz2rgb(XYZ_val_pred)},...
124 'layout', [10, 7],...
125 'squaresize', 100,...
126 'legend', {'Ground-truth', 'Predicted'},...
127 'parent', ax2);
128title('Validation samples comparison');
框选的区域是DIGITAL SG的卡，有140个色块
140分成两部分，其中70为训练集，另外70个为实际结果的验证
3.2demo2
跟demo1区别，加了白点的保持
1clear; close all; clc;
2
3%% data preparation
4
5% load DSG color checker's spectral reflectance data
6load('spectral_reflectance_data.mat');
7spectra = spectral_reflectance_data.XRite_DSG;
8
9% calculate the XYZ values for the color checker under D65
10XYZ = spectra2colors(spectra, 400:5:700,...
11 'spd', 'D65');
12
13% read color checker image and extract color responses
14% the darkness level and the spatial non-uniformity has been corrected for
15% this image
16dsg_img = imread('dsg_Nikon_D3x.tiff');
17
18% you can also manually select roi for the color checker
19roi = [519,310,65,65;669,309,65,65;819,308,65,65;969,307,65,65;1120,306,65,65;1270,306,65,65;1421,305,65,65;1572,304,65,65;1722,303,65,65;1873,302,65,65;2025,301,65,65;2176,300,65,65;2327,299,65,65;2479,299,65,65;519,460,65,65;669,459,65,65;819,458,65,65;970,457,65,65;1120,456,65,65;1271,456,65,65;1421,455,65,65;1572,454,65,65;1723,453,65,65;1874,452,65,65;2025,452,65,65;2176,451,65,65;2328,450,65,65;2479,449,65,65;519,650,65,65;669,609,65,65;820,608,65,65;970,607,65,65;1120,607,65,65;1271,606,65,65;1422,605,65,65;1572,604,65,65;1723,604,65,65;1874,603,65,65;2026,602,65,65;2177,601,65,65;2328,601,65,65;2480,600,65,65;520,760,65,65;670,759,65,65;820,758,65,65;970,758,65,65;1121,757,65,65;1271,756,65,65;1422,755,65,65;1573,755,65,65;1724,754,65,65;1875,753,65,65;2026,753,65,65;2177,752,65,65;2329,751,65,65;2480,751,65,65;520,910,65,65;670,909,65,65;820,908,65,65;971,908,65,65;1121,907,65,65;1272,907,65,65;1422,906,65,65;1573,905,65,65;1724,905,65,65;1875,904,65,65;2027,903,65,65;2178,903,65,65;2329,902,65,65;2481,901,65,65;520,1060,65,65;670,1059,65,65;820,1059,65,65;971,1058,65,65;1121,1057,65,65;1272,1057,65,65;1423,1056,65,65;1574,1056,65,65;1725,1055,65,65;1876,1055,65,65;2027,1054,65,65;2179,1053,65,65;2330,1053,65,65;2482,1052,65,65;520,1210,65,65;670,1209,65,65;821,1209,65,65;971,1208,65,65;1122,1208,65,65;1272,1207,65,65;1423,1207,65,65;1574,1206,65,65;1725,1206,65,65;1876,1205,65,65;2028,1205,65,65;2179,1204,65,65;2331,1204,65,65;2482,1203,65,65;520,1360,65,65;671,1360,65,65;821,1359,65,65;972,1359,65,65;1122,1358,65,65;1273,1358,65,65;1424,1357,65,65;1575,1357,65,65;1726,1357,65,65;1877,1356,65,65;2028,1356,65,65;2180,1355,65,65;2331,1355,65,65;2483,1354,65,65;521,1510,65,65;671,1510,65,65;821,1510,65,65;972,1509,65,65;1122,1509,65,65;1273,1508,65,65;1424,1508,65,65;1575,1508,65,65;1726,1507,65,65;1877,1507,65,65;2029,1507,65,65;2180,1506,65,65;2332,1506,65,65;2484,1505,65,65;521,1661,65,65;671,1660,65,65;822,1660,65,65;972,1660,65,65;1123,1659,65,65;1274,1659,65,65;1424,1659,65,65;1576,1658,65,65;1727,1658,65,65;1878,1658,65,65;2029,1657,65,65;2181,1657,65,65;2332,1657,65,65;2484,1656,65,65];
20RGB = checker2colors(dsg_img, [10, 14],...
21 'roi', roi,...
22 'show', false,...
23 'scale', 4); % scaling only for better visualization
24
25clearvars -except XYZ RGB
26
27
28%% color correction with white point preserved
29
30% white balancing RGB values
31%增加了白平衡RGB值
32neutral_patches_idx = [61, 62, 63, 64, 65];
33gains = [RGB(neutral_patches_idx, 1) \ RGB(neutral_patches_idx, 2),...
34 1,...
35 RGB(neutral_patches_idx, 3) \ RGB(neutral_patches_idx, 2)];
36
37RGB_wb = RGB .* gains;
38
39% use D65's XYZ value as white point to be preserved after color correction
40white_point = whitepoint('d65');
41
42% training
43model = 'root6x3';
44[matrix, scale, XYZ_pred, errs_train] = ccmtrain(RGB_wb,...
45 XYZ,...
46 'model', model,...
47 'targetcolorspace', 'xyz',...
48 'whitepoint', white_point);
49
50% check if [1, 1, 1] has been preserved as [0.9505, 1.0000, 1.0888]
51predicted_white_point = ccmapply([1, 1, 1],...
52 model,...
53 matrix);
54
55white_point_err = sqrt(sum((predicted_white_point - white_point).^2));
56fprintf('residual error between user-specified and predicted white points: %.3e\n', white_point_err);
57
58% visualization
59figureFullScreen('color', 'w');
60ax1 = subplot(1,2,1);
61colors2checker(RGB_wb .^ (1/2.2),...
62 'layout', [10, 14],...
63 'squaresize', 100,...
64 'parent', ax1);
65title('White balaned camera responses before color correction (gamma = 1/2.2)');
66
67ax2 = subplot(1,2,2);
68colors2checker(xyz2rgb(XYZ_pred),...
69 'layout', [10, 14],...
70 'squaresize', 100,...
71 'parent', ax2);
72title('Camera responses (sRGB) after color correction');
参考
[1] Vinegar. CCM调试的理论依据, 2023.
[2] piliweant. CCM 调节笔记, 2021.
[4] 烫手的洋芋. Color Correction Tuning Guide, 2018.