因昨日写的匆忙，有不少错误，特重写并增加了适当的解释，并且觉得对于搞图像处理的朋友来说，这个可能在很多场合还是用得着的，因此删除了原随笔，又重新发布了下，非有意刷屏。
原来一直认为彩色图像的去色算法没啥研究价值，网络上已经有很多类似的算法了，比如著名的Gray = R*0.299 + G*0.587 + B*0.114公式，或者LAB颜色通道的L值，HSL的L通道等等，直到最近看一些论文，发现原来也有很多人对这个过程进行过详细的研究和创新。在学习这些算法的同时，使我也认识到，任何你认为简单的东西在背后都有可能有着复杂的机理，只是你没有发现而已。
我先拿两个图片的例子来说明为什么要对彩图去色的过程进行研究，同时也说明了常用的方法存在的一些缺陷。
下面水平4副图从左至右依次为原图、PS的彩色转灰度命令、PS的去色命令以及matlab的rgb2gray命令的效果。
上面的效果有什么问题，不错，很多人包括我开始也觉得似乎没啥大毛病，可仔细看看，你就会发现第一幅图中的晚霞以及晚霞在水中的倒影和太阳在去色后已经很难找到其踪影了，而第二副图的红花和绿叶去色后基本变得一致了，这其实即是所谓的对比度丢失，这种丢失对于普通的图像处理用户也许问题不大，不过对于图像分析方面是很不利的。
针对这个问题，不少作者都提出了解决方案，我这里列出三篇论文供有兴趣的朋友参考：
（1）Color2Gray: Salience-Preserving Color Removal ，作者Amy A. Gooch Sven C. Olsen Jack Tumblin Bruce Gooch
论文中提及的源码下载地址已经无效了，我也不记得我这是从哪里下载到的，但确实是和原文匹配的代码：http://files.cnblogs.com/Imageshop/Color2GrayMatlabCode.rar
这个代码仅仅具有学习价值，因为作者在论文中说100*100大小的图像算法用时12.7秒，这么长的时间那还搞个屁啊。
（2）Contrast Preserving Decolorization ，作者 Cewu Lu Li Xu Jiaya Jia
以及：
（3）Real-time Contrast Preserving Decolorization，作者 Cewu Lu Li Xu Jiaya Jia
后两篇文章配套的下载地址为：http://appsrv.cse.cuhk.edu.hk/~xuli/mypapers/cpcolor2gray_v2.zip
后两篇文章都是我们国人的杰作，不过不是大陆的，是香港大学的几位高人的作品，我特别关注的是Jiaya Jia，。
特别是最后一篇文章的算法，效果好且执行速度快，能满足实时的要求。
这3篇论文都是以最小化一个能量函数为目标的：
其中gx,gy为灰度化后的像素值。而δx,y则表示颜色对比度，三篇论文中开始的时候都是用的LAB颜色空间的一些相关计算公式。
在第三篇论文中，作者进一步将目标函数改写为：
式中：Nб是高斯分布函数，△gx,y=gx-gy ;
然后优化的目标就是求上式的最大值。
在第三篇论文中提出算法最为实用，他没有像第二篇那样采用了二维的模型，而是简化为一维模型，类似于公式Gray = R*0.299 + G*0.587 + B*0.114，论文中也是使用W1，W2，W3三个系数来控制结果值，但是这三个系数不是固定的，而是同用户输入的图像自适应的。同时有约束条件W1>0;W2>0;W3>0，以及W1+W2+W3=1;满足这几个条件的W1，W2，W3的值的组合还是有无限个的，但是作者注意到系数的微小变化对于输出的结果的影响不是特别大，因此，论文中提出了将每个系数在[0,1]之间量化为10等份，即只取0、0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8，0.9、1.0这11个值，在满足约束条件的情况下，W1、W2、W2的组合总共只会有11*（11+1）/2 =66 种，如下所示：
W = [ 0 0 1.0000 0 0.1000 0.9000 0 0.2000 0.8000 0 0.3000 0.7000 0 0.4000 0.6000 0 0.5000 0.5000 0 0.6000 0.4000 0 0.7000 0.3000 0 0.8000 0.2000 0 0.9000 0.1000 0 1.0000 0 0.1000 0 0.9000 0.1000 0.1000 0.8000 0.1000 0.2000 0.7000 0.1000 0.3000 0.6000 0.1000 0.4000 0.5000 0.1000 0.5000 0.4000 0.1000 0.6000 0.3000 0.1000 0.7000 0.2000 0.1000 0.8000 0.1000 0.1000 0.9000 0 0.2000 0 0.8000 0.2000 0.1000 0.7000 0.2000 0.2000 0.6000 0.2000 0.3000 0.5000 0.2000 0.4000 0.4000 0.2000 0.5000 0.3000 0.2000 0.6000 0.2000 0.2000 0.7000 0.1000 0.2000 0.8000 0 0.3000 0 0.7000 0.3000 0.1000 0.6000 0.3000 0.2000 0.5000 0.3000 0.3000 0.4000 0.3000 0.4000 0.3000 0.3000 0.5000 0.2000 0.3000 0.6000 0.1000 0.3000 0.7000 0.0000 0.4000 0 0.6000 0.4000 0.1000 0.5000 0.4000 0.2000 0.4000 0.4000 0.3000 0.3000 0.4000 0.4000 0.2000 0.4000 0.5000 0.1000 0.4000 0.6000 0.0000 0.5000 0 0.5000 0.5000 0.1000 0.4000 0.5000 0.2000 0.3000 0.5000 0.3000 0.2000 0.5000 0.4000 0.1000 0.5000 0.5000 0 0.6000 0 0.4000 0.6000 0.1000 0.3000 0.6000 0.2000 0.2000 0.6000 0.3000 0.1000 0.6000 0.4000 0.0000 0.7000 0 0.3000 0.7000 0.1000 0.2000 0.7000 0.2000 0.1000 0.7000 0.3000 0.0000 0.8000 0 0.2000 0.8000 0.1000 0.1000 0.8000 0.2000 0.0000 0.9000 0 0.1000 0.9000 0.1000 0.0000 1.0000 0 0];
优化的目的就是从输入图像的数据确定最优的权值使得式（2）的值最大。
考虑到图片内在的颜色的冗余性，为进一步提高速度，论文提出不在原始的图像中搜索最优系数，而是现将图像缩小到一定范围，在缩放后的图中找寻，作者建议缩小到64*64（原图为等宽高的情况）大小，此时即能提高速度，又不会降低精度，但是注意一点就是最好是采用最近邻插值，因为这不会产生新的像素。
上面两种优化方式，特别是搜索空间的这种离散化，我想在很多算法中都可以去模仿的。
贴些和第三篇文章基本对应的matlab代码。
1 function img = rtcprgb2gray(im) 2 3 %% Proprocessing 4 [n,m,ch] = size(im); 5 sigma = 0.05; 6 W = wei(); 7 8 9 %% Global and Local Contrast Computing 10 ims = imresize(im, round(64/sqrt(n*m)*[n,m]),'nearest'); 11 R = ims(:,:,1);G = ims(:,:,2);B = ims(:,:,3); 12 imV = [R(:),G(:),B(:)]; 13 defaultStream = RandStream.getDefaultStream; savedState = defaultStream.State; 14 t1 = randperm(size(imV,1)); 15 Pg = [imV - imV(t1,:)]; 16 17 ims = imresize(ims, 0.5 ,'nearest'); 18 Rx = ims(:,1:end-1,1) - ims(:,2:end,1); 19 Gx = ims(:,1:end-1,2) - ims(:,2:end,2); 20 Bx = ims(:,1:end-1,3) - ims(:,2:end,3); 21 22 Ry = ims(1:end-1,:,1) - ims(2:end,:,1); 23 Gy = ims(1:end-1,:,2) - ims(2:end,:,2); 24 By = ims(1:end-1,:,3) - ims(2:end,:,3); 25 Pl = [[Rx(:),Gx(:),Bx(:)];[Ry(:),Gy(:),By(:)]]; 26 27 P = [Pg;Pl ]; 28 29 det = sqrt(sum(P.^2,2))/1.41 ; 30 31 P( (det < 0.05),:) = []; det( (det < 0.05)) = []; 32 detM = repmat(det,[1,size(W,1)]); L = P*W'; 33 34 %% Energy optimization 35 36 U = log(exp(- (L + detM ).^2/sigma.^2) + exp(- (L- detM).^2/sigma.^2)); 37 Es = mean(U); 38 39 40 41 %% Output 42 [NULLval,bw] = max(Es); 43 img = imlincomb(W(bw,1),im(:,:,1) , W(bw,2),im(:,:,2) , W(bw,3),im(:,:,3)); 44 45 end
第10行即为减少搜索样本。第14行产生不重复的随即对，第27行以及下面计算det的代码未和原文对应，不过似乎不影响结果，第36行计算能量，然后取能量最大的哪个作为权值。
上述代码中PL值其实即计算水平和垂直方向的梯度值，而再次使用ims = imresize(ims, 0.5 ,'nearest');这样的语句也是为了加快计算速度，对为什么用这样的计算方法，原始论文也没有提及，编码的时候似乎去掉这个pl效果也没啥区别。
上述代码如果要在项目中使用，还必须转换为其他的语种，转换应该不存在困难的地方，基本就是一些矩阵或者说数组的操作，唯一比较困难的就是randperm函数，他产生一定范围内的不重复的随机数，这个可以参考http://blog.csdn.net/devfun/article/details/6534476一文，抑或是复制下面的代码：
private static int[] RandPerm(int N) { //http://blog.csdn.net/devfun/article/details/6534476 int X; Random Rnd = new Random(Environment.TickCount); int[] TempArray = new int[N]; int[] Value = new int[N]; for (X = 0; X < N; X++) TempArray[X] = X; for (X = 0; X < N; X++) { int Seed = Rnd.Next(0, N - X); //从剩下的随机数里生成 Value[X] = TempArray[Seed]; //赋值给结果数组 TempArray[Seed] = TempArray[N - X - 1]; //把随机数产生过的位置替换为未被选中的值。 } return Value; }
至于算法的更加原理性的东西，只能看论文本身了吧，实用为王，对于那些提供了参考代码的论文，基本知道原理的意思，然后关键就是掌握代码的改写了。
至于速度，在我I3的笔记本CPU上，参考上述M代码，用C#编制的程序，对于1024*768的图片耗时约为40ms（未优化代码，优化应该能达到20ms)。
而算法效果，我们可以通过以下十几副图像的比较得出结论。
下面水平4副图从左至右依次为原图、PS的彩色转灰度命令、PS的去色命令以论文3的结果。
原图 PS的彩色转灰度 PS的去色命令 论文3的结果
由以上举的一些例子，可以明显的看出论文中的效果要很多，有些图PS转换为灰度后丢失了太多的信息了就变得毫无意义了。
在计算机识别领域，有很多算法需要将彩色先转换为灰度，我想这个时候上述论文中的这种转换效果就非常有意义了。
提供一个C#的测试程序：http://files.cnblogs.com/Imageshop/Decolorization.rar
关于去色，仔细找找，其实还有很多其他的比较成功和优秀的算法，因此，不要小看一些基础的算法，只要你仔细研究，总会有新发现。
*********************************作者： laviewpbt 时间： 2013.11.19 联系QQ: 33184777 转载请保留本行信息************************