因昨日寫的匆忙，有不少錯誤，特重寫並增加了適當的解釋，並且覺得對於搞圖像處理的朋友來說，這個可能在很多場合還是用得着的，因此刪除了原隨筆，又重新發布了下，非有意刷屏。
原來一直認為彩色圖像的去色算法沒啥研究價值，網絡上已經有很多類似的算法了，比如著名的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篇論文都是以最小化一個能量函數為目標的：
（1）
其中gx,gy為灰度化后的像素值。而δx,y則表示顏色對比度，三篇論文中開始的時候都是用的LAB顏色空間的一些相關計算公式。
在第三篇論文中，作者進一步將目標函數改寫為：
（2）
式中：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 轉載請保留本行信息************************