使用Matlab实现图像直方图

Question

我正在试图实现(我知道有一个实现它的自定义函数)Matlab中的灰度图像直方图,到目前为止我已经尝试过:

function h = histogram_matlab(imageSource)
openImage = rgb2gray(imread(imageSource));
[rows,cols] =  size(openImage);

histogram_values = [0:255];

for i = 1:rows
  for j = 1:cols
    p = openImage(i,j);
    histogram_values(p) = histogram_values(p) + 1;

  end
end
histogram(histogram_values)

但是,当我调用该函数时,例如: histogram_matlab('Harris.png')

我得到一些图形:

这显然不是我所期望的,x轴应该从0到255,y轴从0到存储的最大值histogram_values.

我需要获得类似的东西imhist:

我该如何设置？我做得不好吗？

编辑

我已将我的代码更改为@rayryeng建议的改进和更正:

function h = histogram_matlab(imageSource)
openImage = rgb2gray(imread(imageSource));
[rows,cols] =  size(openImage);

histogram_values = zeros(256,1)

for i = 1:rows
  for j = 1:cols
    p = double(openImage(i,j)) +1;
    histogram_values(p) = histogram_values(p) + 1;

  end
end
histogram(histogram_values, 0:255)

但直方图不是预期的:

这里很明显,y轴上存在一些问题或错误,因为它肯定会达到2以上.

Answer 1

在计算直方图方面,虽然存在轻微误差,但每个强度的频率计算是正确的......稍后会有更多.另外,我个人会避免在这里使用循环.请参阅本文末尾的小记.

然而,您的代码有三个问题:

问题#1 - 直方图未正确初始化

histogram_values应该包含你的直方图,但是你用矢量来初始化直方图0:255.每个强度值应以0的计数开始,因此您实际需要这样做:

histogram_values = zeros(256,1);

问题#2 - for循环中出现轻微错误

您的强度范围从0到255,但MATLAB开始索引为1.如果您的强度为0,则会出现越界误差.因此,正确的做法是采取p并添加1,以便您开始索引为1.然而,我需要指出的一个复杂性是,如果你有一个uint8精确的图像,添加1强度为255将简单地将值饱和到255.它不会达到256 ....所以你也要谨慎,double以确保达到256.

因此:

histogram_values = zeros(256,1);
for i = 1:rows
  for j = 1:cols
    p = double(openImage(i,j)) + 1;
    histogram_values(p) = histogram_values(p) + 1;

  end
end

问题#3 - 没有histogram正确的呼叫

您应该覆盖histogram边缘的行为并包括边缘.基本上,这样做:

histogram(histogram_values, 0:255);

第二个向量指定我们应该在x-axis 上放置条形图的位置.

小记

你可以自己完全实现直方图计算而不需要任何for循环.你可以用组合试试这个bsxfun,permute,reshape和两个sum电话:

mat = bsxfun(@eq, permute(0:255, [1 3 2]), im);
h = reshape(sum(sum(mat, 2), 1), 256, 1);

如果您想更详细地解释这段代码的工作原理,请参阅kkuilla和我之间的对话:http://chat.stackoverflow.com/rooms/81987/conversation/explanation-of-computing-an- 图像直方图,矢量

但是,它的要点如下.

代码的第一行创建了由范围从0到255 1列的3D向量permute,然后使用bsxfun与eq(等于)功能,我们使用广播使我们得到一个三维矩阵,其中每个切片的大小与灰度相同图像并为我们提供等于感兴趣强度的位置.具体来说,第一个切片告诉你元素在哪里等于0,第二个切片告诉你哪个元素等于1,直到最后一个切片告诉你元素等于255的位置.

对于第二行代码,一旦我们计算了这个3D矩阵,我们计算两个总和 - 首先独立地对每一行求和,然后对该中间结果的每一列求和.然后我们得到每个切片的总和,它告诉我们每个强度有多少个值.这是一个3D矢量,所以我们reshape回到单个1D矢量来完成计算.

为了显示直方图,我会用bar与histc标志.如果我们使用cameraman.tif图像,这是一个可重现的例子:

%// Read in grayscale image
openImage = imread('cameraman.tif');
[rows,cols] = size(openImage); 

%// Your code corrected
histogram_values = zeros(256,1);
for i = 1:rows
  for j = 1:cols
    p = double(openImage(i,j)) + 1;
    histogram_values(p) = histogram_values(p) + 1;    
  end
end

%// Show histogram
bar(0:255, histogram_values, 'histc');

我们得到这个: