用Matlab实现图像直方图
Image histogram implementation with Matlab
我正在尝试在 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。
您的代码看起来是正确的。问题在于对直方图的调用。您需要在调用 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
在计算直方图方面，每个强度的频率计算是正确的，尽管存在轻微错误……稍后会详细介绍。另外，我个人会避免在这里使用循环。 post.
见文末我的小贴士
然而，您的代码存在三个问题：
问题 #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 and two sum 个调用来尝试：
mat = bsxfun(@eq, permute(0:255, [1 3 2]), im);
h = reshape(sum(sum(mat, 2), 1), 256, 1);
如果你想更详细地解释这段代码是如何在幕后工作的，请看我和 kkuilla 之间的对话：https://chat.whosebug.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');
我们得到这个：
我正在尝试在 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。
您的代码看起来是正确的。问题在于对直方图的调用。您需要在调用 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
在计算直方图方面，每个强度的频率计算是正确的，尽管存在轻微错误……稍后会详细介绍。另外，我个人会避免在这里使用循环。 post.
见文末我的小贴士然而，您的代码存在三个问题：
问题 #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 and two sum 个调用来尝试：
mat = bsxfun(@eq, permute(0:255, [1 3 2]), im);
h = reshape(sum(sum(mat, 2), 1), 256, 1);
如果你想更详细地解释这段代码是如何在幕后工作的，请看我和 kkuilla 之间的对话：https://chat.whosebug.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');
我们得到这个：