用 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-轴上放置条形的位置。

小记

您完全可以自己实现直方图计算，无需任何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 和我之间的对话：https://chat.stackoverflow.com/rooms/81987/conversation/explanation-of-computing-an-images-histogram-vectorized

但是，它的要点如下。

---

第一行代码通过permute创建一个范围从0到255的1列的3D向量，然后使用bsxfun和eq（等于）函数，我们使用广播，所以我们得到一个3D 矩阵，其中每个切片与灰度图像大小相同，并为我们提供与感兴趣强度相等的位置。具体来说，第一个切片告诉你哪里元素等于 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');

我们得到这个：

Answer 2

您的代码看起来正确。问题在于对直方图的调用。您需要在调用 histogram 时提供 bin 的数量，否则它们将被自动计算。

试试这个简单的修改，它调用 stem 来获得正确的图，而不是依赖于直方图

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
stem(histogram_values); axis tight;

编辑：在对代码进行一些检查后，您会遇到 0/1 错误。如果你有一个像素值为零，那么histogram_value(p) 会给你一个索引错误

试试这个。这种简单的情况不需要向量化：

function hv = histogram_matlab_vec(I)

  assert(isa(I,'uint8')); % for now we assume uint8 with range [0, 255]

  hv = zeros(1,256);
  for i = 1 : numel(I)
    p = I(i);
    hv(p + 1) = hv(p + 1) + 1;
  end

  stem(hv); axis tight;
end