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 = γ (r 255.0) G = γ (g 255.0) B = γ (b 255.0) r, g, b \in [0, 255] (1)

$\left\{ {\begin{array}{{c}}{R = \gamma (\frac{r}{{255.0}})}&{}\\{G = \gamma (\frac{g}{{255.0}})}&{r,g,b \in [0,255]}\\{B = \gamma (\frac{b}{{255.0}})}&{}\end{array}}\right. \tag{1}$

γ (x) = ⎧ ⎩ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ (x + 0.055 1.055) 2.4 x 12.92 x > 0.04045 o t h e r w i s e (2)

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

⎡ ⎣ ⎢ X Y Z ⎤ ⎦ ⎥ = M * ⎡ ⎣ ⎢ R G B ⎤ ⎦ ⎥ (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 .433953 0 .212671 0 .017758 0 .376219 0 .715160 0 .109477 0 .189828 0 .072169 0 .872765 ⎤ ⎦ ⎥ (4)

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

L = 116 f (Y Y N) - 16 a = 500 [f (X X N) - f (Y Y N)] b = 200 [f (Y Y N) - f (Z Z N)] (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 \sqrt 3 7.787 t t > 0.008856 o t h e r w i s e (6)

$f(t) = \left\{ {\begin{array}{*{20}{c}}{\begin{array}{*{20}{c}}{\sqrt[3]{t}} \\ {7.787t} \end{array}}&{\begin{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矩阵放大 $2^{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)$ 的计算可以用查表，以为前面的计算总共移位30位，在计算出 $[X,Y,Z]^T$ 后需要先将处理M时所乘的 $2^{20}$ 移回来的， $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 <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
MSRA_0_0_77
与OpenCV中的rgb2lab比较，可以得出如下结果：
RGB2Lab
OpenCV_lab是OpenCV版的结果，mxnLib_lab是本文中代码实现的结果，可以看出二者还是很一致的！
********************************** 我是分割线 ************************************
今天就写到这里把，总之都是第一次，第一次在CSDN写博客，第一次使用CSDN的markdown编辑器，在线学习了LaTex公式格式，很不熟练，花费了整整一晚上的时间才码完。欢迎各位提出意见，如果错误，请提出批评指正。