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）公式：

$$\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}$$
$$\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}$$
$$\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 = \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大神文章中归一化好的值。
$$\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) = \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. 公式(1)和公式(2)的自变量有256个固定的值，即[0,255]，故公式中的除法与乘方运算可以将浮点型变成整型，做成查找表；
2. M矩阵中全部为浮点型数据，需要转换成整形计算；
3. 公式(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

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

OpenCV_lab是OpenCV版的结果，mxnLib_lab是本文中代码实现的结果，可以看出二者还是很一致的！
********************************** 我是分割线 ************************************
今天就写到这里把，总之都是第一次，第一次在CSDN写博客，第一次使用CSDN的markdown编辑器，在线学习了LaTex公式格式，很不熟练，花费了整整一晚上的时间才码完。欢迎各位提出意见，如果错误，请提出批评指正。