Image histogram implementation with Matlab

Question

I'm tyring to implement (I know there's a custom function for achieving it) the grayscale image histogram in Matlab, so far I've tried:

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)

However when I call the function, for example: histogram_matlab('Harris.png')

I obtain some graph like:

which is obviously not what I expect, x axis should go from 0 to 255 and y axis from 0 to whatever max value is stored in histogram_values.

I need to obtain something like what imhist offers:

How should I set it up? Am I doing a bad implementation?

Edit

I've changed my code to improvements and corrections suggested by @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)

However the histogram plot is not that expected:

Here it's noticeable that there's some issue or error on y axis as it definitely would reach MORE than 2.

Answer 1

In terms of calculating the histogram, the computation of the frequency per intensity is correct though there is a slight error... more on that later. Also, I would personally avoid using loops here. See my small note at the end of this post.

Nevertheless, there are three problems with your code:

Problem #1 - Histogram is not initialized properly

histogram_values should contain your histogram, yet you are initializing the histogram by a vector of 0:255. Each intensity value should start with a count of 0, and so you actually need to do this:

histogram_values = zeros(256,1);

Problem #2 - Slight error in for loop

Your intensities range from 0 to 255, yet MATLAB starts indexing at 1. If you ever get intensities that are 0, you will get an out-of-bounds error. As such, the proper thing to do is to take p and add it with 1 so that you start indexing at 1. However, one intricacy I need to point out is that if you have a uint8 precision image, adding 1 to an intensity of 255 will simply saturate the value to 255. It won't go to 256.... so it's also prudent that you cast to something like double to ensure that 256 will be reached.

Therefore:

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

Problem #3 - Not calling histogram right

You should override the behaviour of histogram and include the edges. Basically, do this:

histogram(histogram_values, 0:255);

The second vector specifies where we should place bars on the x-axis.

Small note

You can totally implement the histogram computation yourself without any for loops. You can try this with a combination of bsxfun, permute, reshape and two sum calls:

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

---

If you'd like a more detailed explanation of how this code works under the hood, see this conversation between kkuilla and myself: https://chat.stackoverflow.com/rooms/81987/conversation/explanation-of-computing-an-images-histogram-vectorized

However, the gist of it as below.

---

The first line of code creates a 3D vector of 1 column that ranges from 0 to 255 by permute, and then using bsxfun with the eq (equals) function, we use broadcasting so that we get a 3D matrix where each slice is the same size as the grayscale image and gives us locations that are equal to an intensity of interest. Specifically, the first slice tells you where elements are equal to 0, the second slice tells you where elements are equal to 1 up until the last slice where it tells you where elements are equal to 255.

For the second line of code, once we compute this 3D matrix, we compute two sums - first summing each row independently, then summing each column of this intermediate result. We then get the total sum per slice which tells us how many values there were for each intensity. This is consequently a 3D vector, and so we reshape this back into a single 1D vector to finish the computation.

---

In order to display a histogram, I would use bar with the histc flag. Here's a reproducible example if we use the cameraman.tif image:

%// 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');

We get this: