使用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。

- diegoaguilar

Answer 1

在计算直方图时，每个强度值的频率计算是正确的，尽管存在一些小错误...稍后会详细说明。此外，我个人建议不要在这里使用循环。请参见我的帖子末尾的小注释。

然而，您的代码存在三个问题：

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

histogram_values 应该包含您的直方图，但您正在通过 0:255 的向量来初始化直方图。每个强度值应该从 计数为 0 开始，所以您实际上需要这样做：

histogram_values = zeros(256,1);

问题＃2-for循环中的轻微错误

你的强度范围是0到255，但MATLAB从1开始索引。如果您得到强度为0的强度值，那么将会出现越界错误。因此，正确的做法是将 p 加1，以便从1开始索引。

然而，我需要指出的一个复杂点是，如果您有一个uint8 精度图像，将255的强度值加1将简单地使该值饱和为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向量，然后使用eq（等于）函数的bsxfun进行广播，以便获得一个3D矩阵，其中每个切片的大小与灰度图像相同，并提供与感兴趣的强度相等的位置。具体而言，第一个切片告诉您元素等于0的位置，第二个切片告诉您元素等于1，直到最后一个切片，告诉您元素等于255的位置。

对于第二行代码，一旦我们计算出这个3D矩阵，我们分别对每一行进行求和，然后对这个中间结果的每一列进行求和。然后，我们得到每个强度的总和，因此这是一个3D向量，因此我们将其reshape回一个单独的1D向量以完成计算。

---

为了显示直方图，我将使用histc标志对bar进行设置。 如果我们使用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');

我们得到这个：

- rayryeng

---

抱歉我有点笨，但我从来不理解“permute”。为什么要排列[1 3 2]？为什么不是例如[3 2 1]？由于使用了三个值（即[1 3 2]），您正在创建一个三维矩阵。为什么需要三个维度？im是二维的，为什么不是两个维度？ - kkuilla

谢谢@rayryeng，你能检查一下我的编辑吗？当然，我并不期望得到与imhist生成的图形完全相同的克隆版本，但它甚至看起来都不像是图像直方图。我该如何改进它呢？此外，很明显它需要在y轴上进行调整。 - diegoaguilar

@diegoaguilar - 我不会使用 histogram。我会使用带有 histc 标志的 bar。我会更新我的帖子。 - rayryeng

谢谢@rayryeng。我已经实现了，有一个问题，我想展示实现、内置函数和最终图像在中心的效果，我该怎么做？这是我正在尝试的：http://kopy.io/1oiuZ - diegoaguilar

是的，我想知道哪一个更快。我打赌是 a(1,1,:) = 0:255;。 - kkuilla

Answer 2

你的代码看起来是正确的。问题在于调用直方图时需要提供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

- bendervader

---

小注：histogram_values 应该最初全部为零，但在您的次要函数中，这是正确的。 - rayryeng