用MATLAB求影象直方圖的演算法
Matlab的影象工具箱中有計算直方圖的函式imhist。然而，課程設計總是有很多限制，比如這次的影象處理課程設計，對於影象處理工具箱的使用是有限制的。
所以得自己寫計算直方圖的演算法。我看了一下imhist的程式碼，發現它呼叫了MEX，所以速度很快。可是我對於如何編寫MEX檔案沒有研究，手頭資料又比較有限，而且時間也很倉促，這兩週對付了六門考試和兩個課程設計……
我以前只知道m語言的迴圈慢，但有多慢，不大瞭解。
一直在考慮如何避免使用迴圈，想到的一種辦法是將影象的2D矩陣轉換成向量，並和向量0:255利用meshgrid構成兩個矩陣，然後對這兩個矩陣進行‘==’運算，對每列求和即可得到直方圖，但是我算了一下，發現這種方法記憶體消耗是非常驚人的，以至於在Matlab裡試驗時Matlab直接提示記憶體不足，而無法執行。
後來決定使用一種折中的方法，使用一個規模不是很大的迴圈來進行計算，在這種思想的指導下，產生了如下的程式碼：
function bars=histogram(I)
%用==來提取某個灰度的畫素
%並用sum來計算個數
tic
bars=zeros(1,256);
for value=0:255
bars(value+1)=sum(value==I(:));
end
bars=bars./numel(I);
toc
tic和toc是用來設定計時器，以測試函式的效能。
以如下方式使用這個函式：
首先讀取一幅影象，例如：
>>RGB=imread('1.jpg');
轉換為灰度圖：
>>I=rgb2gray(RGB);
獲取直方圖：
>>bars=histogram(I);
顯示直方圖：
>>bar(0:255,bars);
有一點需要注意的是得到的直方圖向量是進行了歸一化的，即bars向量的每個分量只是對應的灰度在原圖中佔有的比例，而不是實際數值，這樣做是為了之後的直方圖均衡化處理更方便。
這段程式碼的效能如何？對一幅540*720的灰度圖象，耗時3.304s。
所以這段程式碼實在是太慢了，看來迴圈實在是m語言的一大忌諱。
那麼能不能找到一種不用迴圈的方法呢？我又想了一種方法：
首先將bars構造成一個1x256的0向量，然後使用如下語句計算直方圖：
bars(I(:))=bars(I(:))+1
我考慮也許這樣就可以了，但是事與願違，這樣只會使影象出現過的畫素值在bars的對應位置中置1，而不會計算總數。
今天中午吃飯的時候，我覺得可以利用diff（差分）的方法來避免迴圈，以下用例子來說明想法：
有一個影象矩陣I：
>>I=[1 2 0;3 2 3;4 3 1]
I =
1 2 0
3 2 3
4 3 1
將I變成一個向量並排序：
>> I=sort(I(:)).'
I =
0 1 1 2 2 3 3 3 4
對I求diff：
>> A=diff(I)
A =
1 0 1 0 1 0 0 1
在A的末尾補一個1：
>> A=[A 1]
A =
1 0 1 0 1 0 0 1 1
將A與向量1:numel(I)對應項相乘：
>> A=A.*(1:numel(I))
A =
1 0 3 0 5 0 0 8 9
將A中的非0元素取出：
>> A=nonzeros(A).'
A =
1 3 5 8 9
在A的首部補一個0：
>> A=[0 A]
A =
0 1 3 5 8 9
再求一次diff：
>> A=diff(A)
A =
1 2 2 3 1
通過A的結果我們可以看到，原影象中有一個0，兩個1，兩個2，三個3，一個4。即最後的A就是影象的直方圖。
利用這種思路，我重寫了histogram函式（注意現在的histogram還有一個嚴重的bug，所以不要使用下面的程式碼）function bars=histogram(I)
tic
bars=double(diff([0;nonzeros([diff(uint32(sort(I(:))));1].*uint32(1:numel(I)).')]).')./numel(I);
toc
和先前那個函式一樣，對結果進行了歸一化。
這段程式碼的效率怎樣呢？對付那幅540x720的影象，耗時0.12s，比開始的那個histogram快了27.5倍。
但是這段程式碼有一個嚴重的bug，吃晚飯的時候我想，如果影象中沒有某個灰度的畫素，那麼這段程式碼肯定會得到極其錯誤的結果，例如：
>> I=[0 1 1;3 1 4;1 0 1]
I =
0 1 1
3 1 4
1 0 1
這個矩陣裡沒有2。重複上面的過程，看我們得到什麼樣的結果：
>> I=sort(I(:)).'
I =
0 0 1 1 1 1 1 3 4
>> A=diff(I)
A =
0 1 0 0 0 0 2 1
>> A=[A 1]
A =
0 1 0 0 0 0 2 1 1
>> A=A.*(1:numel(I))
A =
0 2 0 0 0 0 14 8 9
>> A=nonzeros(A).'
A =
2 14 8 9
>> A=[0 A]
A =
0 2 14 8 9
>> A=diff(A)
A =
2 12 -6 1
可以看出，完全是牛頭不對馬嘴。為什麼會這樣，原因就在於原矩陣中缺乏灰度為2的畫素導致排序後的向量不是連續地變化。怎樣解決這個看似棘手的問題呢？吃完飯的時候我就想出來了：將原來的影象補充每種灰度的畫素一個，然後對求得的直方圖向量每個元素減一，即可，如：
>> I=[0 1 1;3 1 4;1 0 1]
I =
0 1 1
3 1 4
1 0 1
>> I=sort([I(:).' 0:4])
I =
0 0 0 1 1 1 1 1 1 2 3 3 4 4
>> A=diff(I)
A =
0 0 1 0 0 0 0 0 1 1 0 1 0
>> A=[A 1]
A =
0 0 1 0 0 0 0 0 1 1 0 1 0 1
>> A=A.*(1:numel(I))
A =
0 0 3 0 0 0 0 0 9 10 0 12 0 14
>> A=nonzeros(A).'
A =
3 9 10 12 14
>> A=[0 A]
A =
0 3 9 10 12 14
>> A=diff(A)
A =
3 6 1 2 2
>> A=A-1
A =
2 5 0 1 1
從A可看出，0有1個，1有5個，2有0個，3有一個，4有一個。也就是通過補充每種灰度的畫素各一個的方式使排序後的向量連續，最後減1即可去除加入的畫素的影響，所以最終實用的histogram函式如下：
function bars=histogram(I)
tic
I=[I(:);(0:255).'];
bars=(double(diff([0;nonzeros([diff(uint32(sort(I)));1].*uint32(1:numel(I)).')]).')-1)./(numel(I)-256);
toc
所以m語言的迴圈實在是……你看上面的程式碼，又是sort又是diff，還比使用了迴圈的程式碼快那麼多。
晚上的時候我想，如果用傳統的思路寫這個函式，會有多慢呢？試了一下：
function bars=histogram(I)
tic
bars=zeros(1,256);
I=uint16(I);
for v=1:numel(I)
bars(I(v)+1)=bars(I(v)+1)+1;
end
bars=bars./numel(I);
toc
我想它一定很慢，因為對付那幅540x720的影象，它要迴圈388800次。
可是結果卻出乎我的意料：耗時0.04s，是我寫的幾個函式中速度最快的了！居然會有這種事情！我費盡心機不用迴圈，然而最簡單的使用迴圈的程式碼卻成了最快的程式碼！糊塗了。到底怎樣才能寫出最快的程式碼呢？