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

第一次写博客，没有太多经验，初入图像处理方向的炒鸡菜鸟，看着大神们都有自己的博客，而且总是学习大神的文章，心里除了满满的敬佩，还有一丝丝的失落，后悔自己当初读本科时没学计算机方向，现在研究僧半路出家学图像，没有太多计算机基础，很多东西都是自己重新学习摸索，非常吃力，编程神马的真是效率超低，本文的这个小程序整整花费一天才写好，，，无语了。。。
废话不多说，作为处女作，难免会有很多写的不周到的地方，还请大家批评指正。
************************ 忧桑的分割线 ***************************
此文主要参照了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)

\left\{ {\begin{array}{{c}}{R = \gamma (\frac{r}{{255.0}})}\\{G = \gamma (\frac{g}{{255.0}})}\\{B = \gamma (\frac{b}{{255.0}})}\end{array}}\right. \tag{1}

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

\gamma (x) = \left\{ {\begin{array}{*{20}{c}} {{{\left( {\frac{{x + 0.055}}{{1.055}}} \right)}^{2.4}}}\\ {\frac{x}{{12.92}}}\end{array}} \right.\tag{2}

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

\left[ {\begin{array}{*{20}{c}}X\\Y\\Z\end{array}} \right] = M * \left[ {\begin{array}{*{20}{c}}R\\G\\B\end{array}} \right]\tag{3}

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

M = \left[ {\begin{array}{*{20}{c}}{{\text{0}}{\text{.433953}}}}{\text{.376219}}}}{\text{.189828}}} \\ {{\text{0}}{\text{.212671}}}}{\text{.715160}}}}{\text{.072169}}} \\ {{\text{0}}{\text{.017758}}}}{\text{.109477}}}}{\text{.872765}}} \end{array}} \right]\tag{4}
M的取值是直接取了laviewpbt大神文章中归一化好的值。

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

\begin{gathered}L = 116f\left( {\frac{Y}{{{Y_N}}}} \right) - 16 \\a = 500\left[ {f\left( {\frac{X}{{{X_N}}}} \right) - f\left( {\frac{Y}{{{Y_N}}}} \right)} \right] \\b = 200\left[ {f\left( {\frac{Y}{{{Y_N}}}} \right) - f\left( {\frac{Z}{{{Z_N}}}} \right)} \right] \\ \end{gathered}\tag{5}

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

f(t) = \left\{ {\begin{array}{*{20}{c}}{\begin{array}{*{20}{c}}{\sqrt[3]{t}} \\ {7.787t} \end{array}}{*{20}{c}}{t > 0.008856} \\ {otherwise} \end{array}} \end{array}} \right.\tag{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矩阵放大2202^{20}倍，即左移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)f(t)的计算可以用查表，以为前面的计算总共移位30位，在计算出[X,Y,Z]T[X,Y,Z]^T后需要先将处理M时所乘的2202^{20}移回来的，f(t)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);

LL的值域为[0,100]，aa和b<script type="math/tex" id="MathJax-Element-14">b</script>的值域分别为[-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[b] * 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是本文中代码实现的结果，可以看出二者还是很一致的！
********************************** 我是分割线 ************************************
今天就写到这里把，总之都是第一次，第一次在CSDN写博客，第一次使用CSDN的markdown编辑器，在线学习了LaTex公式格式，很不熟练，花费了整整一晚上的时间才码完。欢迎各位提出意见，如果错误，请提出批评指正。