用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，是我写的几个函数中速度最快的了！居然会有这种事情！我费尽心机不用循环，然而最简单的使用循环的代码却成了最快的代码！糊涂了。到底怎样才能写出最快的代码呢？