爱华网Matlab的图像工具箱中有计算直方图的函数imhist。然而，课程设计总是有很多限制，比如这次的图像处理课程设计，对于图像处理工具箱的使用是有限制的。
所以得自己写计算直方图的算法。我看了一下imhist的代码，发现它调用了MEX，所以速度很快。可是我对于如何编写MEX文件没有研究，手头资料又比较有限，而且时间也很仓促，这两周对付了六门考试和两个课程设计……
我以前只知道m语言的循环慢，但有多慢，不大了解。
一直在考虑如何避免使用循环，想到的一种办法是将图像的2D矩阵转换成向量，并和向量0:255利用meshgrid构成两个矩阵，然后对这两个矩阵进行‘==’运算，对每列求和即可得到直方图，但是我算了一下，发现这种方法内存消耗是非常惊人的，以至于在Matlab里试验时Matlab直接提示内存不足，而无法执行。
后来决定使用一种折中的方法，使用一个规模不是很大的循环来进行计算，在这种思想的指导下，产生了如下的代码：
function bars=hi stogram(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，所以不要使用下面的代码）functionbars=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，是我写的几个函数中速度最快的了！居然会有这种事情！我费尽心机不用循环，然而最简单的使用循环的代码却成了最快的代码！糊涂了。到底怎样才能写出最快的代码呢？
