RGB图像转换到CIELab空间的研究及优化

第一次写博客，没有太多经验，初入图像处理方向的炒鸡菜鸟，看着大神们都有自己的博客，而且总是学习大神的文章，心里除了满满的敬佩，还有一丝丝的失落，后悔自己当初读本科时没学计算机方向，现在研究僧半路出家学图像，没有太多计算机基础，很多东西都是自己重新学习摸索，非常吃力，编程神马的真是效率超低，本文的这个小程序整整花费一天才写好，，，无语了。。。

废话不多说，作为处女作，难免会有很多写的不周到的地方，还请大家批评指正。

[b]************************[/b] 忧桑的分割线 [b]***************************[/b]

此文主要参照了laviewpbt大神的两篇文章，目前仅仅做了RGB到Lab的转换，反变换暂时未做

CIELAB颜色空间的更多原理说明见：

http://en.wikipedia.org/wiki/Lab_color_space

laviewpbt大神的两篇文章:

http://www.cnblogs.com/Imageshop/archive/2013/01/31/2888097.html

http://www.cnblogs.com/Imageshop/archive/2013/02/02/2889897.html

下面仅仅是简单说明一下，因为细节问题在以上两篇博文里已经说的很详细了

（1）主要流程: RGB->XYZ->Lab

（2）公式：

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪R=γ(r255.0)G=γ(g255.0)B=γ(b255.0)r,g,b∈[0,255](1)

γ(x)=⎧⎩⎨⎪⎪⎪⎪⎪⎪(x+0.0551.055)2.4x12.92x>0.04045otherwise(2)

⎡⎣⎢XYZ⎤⎦⎥=M∗⎡⎣⎢RGB⎤⎦⎥(3)

M=⎡⎣⎢0.4339530.2126710.0177580.3762190.7151600.1094770.1898280.0721690.872765⎤⎦⎥(4)

M的取值是直接取了laviewpbt大神文章中归一化好的值。

L=116f(YYN)−16a=500[f(XXN)−f(YYN)]b=200[f(YYN)−f(ZZN)](5)

f(t)={t√37.787tt>0.008856otherwise(6)

XYZ->Lab的公式参考 冈萨雷斯的 数字图像处理 第3版 阮秋琦译 p273

(3) 优化

通过以上公式可以发现，整个计算过程中涉及到大量的浮点运算，对于PC应该还可以顶住，如果你的程序要在嵌入式系统中跑，这么大量的浮点运算对CPU是很大的开销，所以需要对其进行优化。优化的方法首先想到就是查表，用空间换时间，其次就是用移位代替除法，用整形替换浮点型。主要的需要用到优化的有以下几点：

公式(1)和公式(2)的自变量有256个固定的值，即[0,255]，故公式中的除法与乘方运算可以将浮点型变成整型，做成查找表；

M矩阵中全部为浮点型数据，需要转换成整形计算；

公式(6)又涉及到浮点型数据运算和有限个自变量的问题，依然需要转换类型并做成查找表；

RGB三通道的值在[0,255]之间变化，而且公式(2)的值域在[0,1]之间，为将浮点型转换成整型，对其放大1024倍，即左移10位，实现如下：

int gammatable[256];
memset(gammatable, 0, 256 * sizeof(int));
for(i = 0; i < 11; i++)
gammatable[i] = (i / 3294.6) * 1024;
for(; i < 256; i++)
gammatable[i] = (pow((i + 14.025) / 269.025, 2.4)) * 1024;

M矩阵放大220倍，即左移20位，用Matlab计算得到数据如下（我承认我偷懒了）：

M3 =1.0e+05 *
4.5503    3.9449    1.9905
2.2300    7.4990    0.7567
0.1862    1.1479    9.1516

f(t)的计算可以用查表，以为前面的计算总共移位30位，在计算出[X,Y,Z]T后需要先将处理M时所乘的220移回来的，f(t)的值域为[0,1]，充分利用还剩下的10位未移回来的位，将其放大1024倍，即：XYZ的三个分量仅有1024种取值，可以作为Lab的输入。

int ftable[1024];
memset(ftable, 0, 1024 * sizeof(int));
for(i = 0; i < 1024; i++)
ftable[i] = (i > 9) ? pow((float)i / 1024, 1.0 / 3) * 1024 + 0.5 : (7.787 * i + 141.2413);

L的值域为[0,100]，a和b的值域分别为[-169,169]，[-160,160]，至于怎么知道的，我也忘了从哪里看的，总之它是对的，有了这个值域，就可以做查找表了，而且在之后的归一化到[0,255]范围内以方便显示查看是有帮助的。

int Ltable[100];
int atable[338];
int btable[320];
memset(Ltable, 0, 100 * sizeof(int));
memset(atable, 0, 338 * sizeof(int));
memset(btable, 0, 320 * sizeof(int));
for(i = 0; i < 100; i++)
Ltable[i] = (float)i / 100 * 255;
for(i = 0; i < 338; i++)
atable[i] = (float)i / 338 * 255;
for(i = 0; i < 320; i++)
btable[i] = (float)i / 320 * 255;

如此全部的优化工作做完，可以开森的查表了（可能有人会说：L，a，b的计算中只有整型除法，没必要查表，，，对于一个强迫症晚期的人来说，这种缺陷是无法容忍的）

(4) 实现

全部工作做完，给出全部代码，纯C版（其实我对C++很凌乱，不要鄙视我，我是半路出家学这玩意儿的，刚入门）

/**
* @brief RGB空间转换到Lab空间，结果已经归一化到[0,255]以方便显示，效果与OpenCV自带函数效果一致
* @param rgb
* @param width
* @param height
* @return lab空间图像数据首地址
*/
unsigned char* rgb2lab(unsigned char* rgb, int width, int height)
{
unsigned char *lab;
int gammatable[256];
int ftable[1024];
int Ltable[100];
int atable[338];
int btable[320];
int L, a1, b1;
long long *X, *Y, *Z;
int i, r, g, b;
unsigned char *pCur, *pEnd;
if((lab = (unsigned char*)malloc(width * height * 3 * sizeof(unsigned char))) == NULL)
{
printf("No space distributed for CIE_Lab!\n");
return NULL;
}
if((X = (long long*)malloc(width * height * sizeof(long long))) == NULL)
{
printf("No space distributed for X in CIE_XYZ!\n");
return NULL;
}
if((Y = (long long*)malloc(width * height * sizeof(long long))) == NULL)
{
printf("No space distributed for Y in CIE_XYZ!\n");
return NULL;
}
if((Z = (long long*)malloc(width * height * sizeof(long long))) == NULL)
{
printf("No space distributed for Z in CIE_XYZ!\n");
return NULL;
}
//创建伽马校正的查找表,将浮点数转变成整数，放大1024倍
memset(gammatable, 0, 256 * sizeof(int));
for(i = 0; i < 11; i++)
gammatable[i] = (i / 3294.6) * 1024;
for(; i < 256; i++)
gammatable[i] = (pow((i + 14.025) / 269.025, 2.4)) * 1024;
// XYZ2LAB查找表
memset(ftable, 0, 1024 * sizeof(int));
for(i = 0; i < 1024; i++)
ftable[i] = (i > 9) ? pow((float)i / 1024, 1.0 / 3) * 1024 + 0.5 : (7.787 * i + 141.2413);
//归一化到[0,255]的查找表（由于涉及到浮点除法）
memset(Ltable, 0, 100 * sizeof(int));
memset(atable, 0, 338 * sizeof(int));
memset(btable, 0, 320 * sizeof(int));
for(i = 0; i < 100; i++)
Ltable[i] = (float)i / 100 * 255;
for(i = 0; i < 338; i++)
atable[i] = (float)i / 338 * 255;
for(i = 0; i < 320; i++)
btable[i] = (float)i / 320 * 255;
for(i = 0, pCur = rgb, pEnd = rgb + width * height * 3; pCur < pEnd; i++)
{
b = *pCur++;
g = *pCur++;
r = *pCur++;
//M矩阵系数被放大了2^20倍，结果需要移位20位
X[i] = (gammatable * 199049 + gammatable[g] * 394494 + gammatable[r] * 455033 + 524288) >> 20;
Y[i] = (gammatable[b] * 75675 + gammatable[g] * 749900 + gammatable[r] * 223002 + 524288) >> 20;
Z[i] = (gammatable[b] * 915161 + gammatable[g] * 114795 + gammatable[r] * 18621 + 524288) >> 20;
}
for(i = 0, pCur = lab, pEnd = lab + width * height * 3; pCur < pEnd; i++)
{
L = (116 * ftable[Y[i]] - 16 * 1024) >> 10;
a1 = (500 * (ftable[X[i]] - ftable[Y[i]])) >> 10;
b1 = (200 * (ftable[Y[i]] - ftable[Z[i]])) >> 10;
*pCur++ = Ltable[L];
*pCur++ = atable[a1 + 169];
*pCur++ = btable[b1 + 160];
}
free(X);
free(Y);
free(Z);
return lab;
}

原图如下，来源于MSRA10K



与OpenCV中的rgb2lab比较，可以得出如下结果：



OpenCV_lab是OpenCV版的结果，mxnLib_lab是本文中代码实现的结果，可以看出二者还是很一致的！

[b][b]**********************************[/b] 我是分割线 [b]************************************[/b]

今天就写到这里把，总之都是第一次，第一次在CSDN写博客，第一次使用CSDN的markdown编辑器，在线学习了LaTex公式格式，很不熟练，花费了整整一晚上的时间才码完。欢迎各位提出意见，如果错误，请提出批评指正。