two-dimensional Otsu's method

I am trying to implement 2D Otsu segmentation method, but i am facing problem in my code. in 2D otsu the gray-level value of each pixel as well as the average value of its immediate neighborhood is studied so that the binarization results are greatly improved. I am attaching code,but it is not working and also hyperlink http://en.wikipedia.org/wiki/Otsu%27s_method

function inputs and output:

%hists is a 256\times 256 2D-histogram of grayscale value and neighborhood average grayscale value pair.
%total is the number of pairs in the given image.
%threshold is the threshold obtained.

function threshold = 2D_otsu(hists, total)
maximum = 0.0;
threshold = 0;
helperVec = 0:255;
mu_t0 = sum(sum(repmat(helperVec',1,256).*hists));
mu_t1 = sum(sum(repmat(helperVec,256,1).*hists));
p_0 = zeros(256);
mu_i = p_0;
mu_j = p_0;
for ii = 1:256
    for jj = 1:256
        if jj == 1
            if ii == 1
                p_0(1,1) = hists(1,1);
            else
                p_0(ii,1) = p_0(ii-1,1) + hists(ii,1);
                mu_i(ii,1) = mu_i(ii-1,1)+(ii-1)*hists(ii,1);
                mu_j(ii,1) = mu_j(ii-1,1);
            end
        else
            p_0(ii,jj) = p_0(ii,jj-1)+p_0(ii-1,jj)-p_0(ii-1,jj-1)+hists(ii,jj);
            mu_i(ii,jj) = mu_i(ii,jj-1)+mu_i(ii-1,jj)-mu_i(ii-1,jj-1)+(ii-1)*hists(ii,jj);
            mu_j(ii,jj) = mu_j(ii,jj-1)+mu_j(ii-1,jj)-mu_j(ii-1,jj-1)+(jj-1)*hists(ii,jj);
        end

          if (p_0(ii,jj) == 0)
              continue;
          end
          if (p_0(ii,jj) == total)
              break;
          end
          tr = ((mu_i(ii,jj)-p_0(ii,jj)*mu_t0)^2 + (mu_j(ii,jj)-p_0(ii,jj)*mu_t0)^2)/(p_0(ii,jj)*(1-p_0(ii,jj)));

          if ( tr >= maximum )
              threshold = ii;
              maximum = tr;
          end
      end
    end

[Geoff Hayes]

pramod - please clarify what you mean by it is not working. Are you observing an error, and if so, what is it (please copy and paste the full error message)? Or, are the results not what you expect. If this is the case, please describe your inputs (which you can attach as a mat file), what the outputs from the above are, and what they should be. Also, in the future, rather than copying and pasting your code into the question body, just attach the code (m file).

[Harsha]

Hi Geoff Hayes thank you for your reply, I want to segment lung using 2d Otsu method. Because 2d Otsu is improved one compare to 1D Otsu (graythresh() and multithresh()),2D Otsu threshold by averaging image. I am referring : https://en.wikipedia.org/wiki/Otsu's_method

Here I am attaching ZIP file contain test.m (main file) Otsu_2D(2D otsu logic) 1.tif(image)

Following Error:-

Attempted to access p_0(0,2); index must be a positive integer or logical.

Error in otsu_2D (line 26)
             p_0(ii,jj) = p_0(ii,jj-1)+p_0(ii-1,jj)-p_0(ii-1,jj-1)+hists(ii,jj);

Error in t2 (line 56)
 thresh = otsu_2D(hist2d, total) %threshold obtained 2D Otsu

[Prabhu Bevinamarad]

Hi,

I am new to Matlab. I wanted to apply two-dimensional Otsu's method for thresholding. I am not getting how to find hists i.e. 2D-histogram of grayscale value and neighborhood average grayscale value pair and total is the number of pairs in the given image which are passed as a parameters for otsu_2D function.Can anybody suggest which functions are used.

Thank you

Answer 1

[Geoff Hayes]

pramod - the error message is telling you that the code is trying to access an element within your matrix using an index that is zero. I realize that you are considering the neighbourhood surrounding each cell within your 256x256 matrix, but your code will have to account for the edge cases when subtracting one from the index leads to a zero which "falls" outside of your matrix. You will have to guard against this code as follows

 for ii = 1:256
    for jj = 1:256
        if jj == 1
            if ii == 1
                p_0(1,1) = hists(1,1);
            else
                p_0(ii,1) = p_0(ii-1,1) + hists(ii,1);
                mu_i(ii,1) = mu_i(ii-1,1)+(ii-1)*hists(ii,1);
                mu_j(ii,1) = mu_j(ii-1,1);
            end
        else
            % here is the new code  ***** guard against ii==1
            if ii==1
                p_0(ii,jj) = p_0(ii,jj-1)+hists(ii,jj);
                mu_i(ii,jj) = mu_i(ii,jj-1)+(ii-1)*hists(ii,jj);
                mu_j(ii,jj) = mu_j(ii,jj-1)+(jj-1)*hists(ii,jj);     
            else
                p_0(ii,jj) = p_0(ii,jj-1)+p_0(ii-1,jj)-p_0(ii-1,jj-1)+hists(ii,jj);
                mu_i(ii,jj) = mu_i(ii,jj-1)+mu_i(ii-1,jj)-mu_i(ii-1,jj-1)+(ii-1)*hists(ii,jj);
                mu_j(ii,jj) = mu_j(ii,jj-1)+mu_j(ii-1,jj)-mu_j(ii-1,jj-1)+(jj-1)*hists(ii,jj);
            end
        end

I don't know for sure if the above is correct but it will guard against the case where jj is greater than one and ii is equal to one (the pattern is similar to what is already there for the case where jj is one and ii is one). A note should be added to the Wikipedia article indicating that there is a bug with the code.

ImageAnalyst would have a better idea as to whether the above makes sense or not.

And, rather than posting a binary zip file which some of us may be reluctant to open, just attach each file individually (since there are just three).

[Harsha]

Thank you Geoff Hayes, but I am getting threshold as 0 for above change. My aim is to get threshold from 2D Otsu method.Here I am attaching 2D Otsu method computation,can you pleas help..

[Geoff Hayes]

pramod - you will need to attach your inputs, the hists and total so that we can try to figure out what is incorrect with the above implementation.

[Harsha]

hists is a 256x256 because (L x L) 2D-histogram of grayscale value and neighborhood average grayscale value pair. total is the number of pairs in the given image(256x256 because image size).Here I am attaching test.m having hists and total variable calculated but as 2D Otsu computation I unable to get threshold,can you please help.I calculated hists by hist2d as 2D histogram of two images (image and average of that).

[Geoff Hayes]

pramod - I was hoping for the inputs to the function (the hists and the total). I don't have the Image Processing Toolbox so cannot run the code that you provided in test.m.

[Geoff Hayes]

pramod - since the only answer you are getting is a threshold value of zero, then we should try to figure out why that may be. This variable is set as

        if ( tr >= maximum )
            threshold = ii;
            maximum = tr;
        end

So threshold will be set to ii (which is an integer in the interval [1,256]) only if the trace, tr, is greater than maximum which is initialized to zero. This implies that tr is never greater than zero and so is either zero or negative. The trace is initialized as

 tr = ((mu_i(ii,jj)-p_0(ii,jj)*mu_t0)^2 + ...
       (mu_j(ii,jj)-p_0(ii,jj)*mu_t0)^2)/(p_0(ii,jj)*(1-p_0(ii,jj)));

Since we are squaring the first two terms and adding them together, then the numerator will always be positive. The problem then is with the denominator

 (p_0(ii,jj)*(1-p_0(ii,jj)))

The name p_0 suggests that this variable should be a probability matrix with values between 0 and 1. This should guarantee that 1-p_0(ii,jj) should be positive but... p_0 is initialized using hists whose name suggests that its elements are integer counts representing frequency whose values are 0,1,2,... etc. which is exactly how you have constructed this input in your test.m code. Try doing the following in test.m

thresh = otsu_2D(hist2d./total, total);

so that the first input to this function is a probability matrix (this should (?) be the the joint probability mass function in 2-dimensional histogram from the Wikipedia entry). Then in the otsu_2d code, replace

 if (p_0(ii,jj) == total)

with

 if fabs(p_0(ii,jj) - 1.0) < eps

or comment it out altogether. This should ensure that the trace is positive and that you should obtain a non-zero threshold.

I think that there could be other flaws with this implementation of the algorithm. It is unfortunate that the author did not comment the code better or even verify that it works before posting it to Wikipedia. (Or do a better job of mapping the discussed algorithm to the code to make it easier to follow.)

[Harsha]

Geoff Hayes- for above changes I am getting threshold = 256,but I am trying to map the discussed algorithm to the code. For a reference I am attaching results for above changes.

hists is an 2d-histogram of image with size M×N can be represented by a 2D gray level intensity function f(x, y). The value of f(x, y) is the gray level, ranging from 0 to L-1, where L is the number of distinct gray levels. In a 2D thresholding method, the gray level of a pixel and its local average gray level are both used. The local average gray level is also divided into the same L values, let g(x, y) be the function of the local average gray level.f(x,y) is an input image and g(x,y) average gray level of input image. total is the number of pixels in the input image(i.e M x N).

[Harsha]

Hi Geoff Hayes- finally I got something, using 1D Otsu(i.e multithresh in MATLAB) I got threshold as 88.using 2D Otsu method I got threshold as 85.From this, one of the advantage I got is When two-dimensional Otsu method is adopted, gray mean is considered, thus it can better segment the images with Gauss noise. Here I attaching 2D Otsu with results,test.m is same as above attached one. Can you please check that for me and reply.

[Klaus Grouleff]

Hi Harsha, how did you generate the hist input for otsu_2D.m?

[Harsha]

hists is an 2d-histogram of image with size M×N can be represented by a 2D gray level intensity function f(x, y). The value of f(x, y) is the gray level, ranging from 0 to L-1, where L is the number of distinct gray levels. In a 2D thresholding method, the gray level of a pixel and its local average gray level are both used. The local average gray level is also divided into the same L values, let g(x, y) be the function of the local average gray level.f(x,y) is an input image and g(x,y) average gray level of input image. total is the number of pixels in the input image(i.e M x N).

[Image Analyst]

No idea what this means exactly but as far as I can tell, it's this:

% Read in input image f(x,y) with L = 256 gray levels:
fxy = imread('cameraman.tif');
% No idea what "The local average gray level is also divided into the same L values" means,
% but to get an image where each pixel is the average in a square window around the pixel, do thi
windowWidth = 3; % Whatever you want.
kernel = ones(windowWidth, windowWidth) / windowWidth^2;
gxy = imfilter(fxy, kernel);
% Now get the total number of pixels in fxy and gxy:
total = numel(fxy)

Note that images are indexed fxy(y, x), which is fxy(row, column), and NOT as fxy(x, y).

But gxy is simply a blurred input image -- the local average. It is not a 2-D histogram. I'm wondering if he wants 256 histograms, one for the neighbors of each pixel (not including the pixel itself), like

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 18;
fprintf('Beginning to run %s.m ...\n', mfilename);
echo off;

%===============================================================================
% Get the name of the demo image the user wants to use.
% Let's let the user select from a list of all the demo images that ship with the Image Processing Toolbox.
folder = fileparts(which('cameraman.tif')); % Determine where demo folder is (works with all versions).
% Demo images have extensions of TIF, PNG, and JPG.  Get a list of all of them.
imageFiles = [dir(fullfile(folder,'*.TIF')); dir(fullfile(folder,'*.PNG')); dir(fullfile(folder,'*.jpg'))];
for k = 1 : length(imageFiles)
	% 	fprintf('%d: %s\n', k, files(k).name);
	[~, baseFileName, extension] = fileparts(imageFiles(k).name);
	ca{k} = [baseFileName, extension];
end
% Sort the base file names alphabetically.
[ca, sortOrder] = sort(ca);
imageFiles = imageFiles(sortOrder);
button = menu('Use which gray scale demo image?', ca); % Display all image file names in a popup menu.
% Get the base filename.
baseFileName = imageFiles(button).name; % Assign the one on the button that they clicked on.
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);

grayImage = imread(fullFileName);
[rows, columns, numberOfColorChannels] = size(grayImage);
if numberOfColorChannels == 3
	grayImage = rgb2gray(grayImage);
end
subplot(2, 1, 1);
imshow(grayImage);
caption = sprintf('Original gray scale image : %s', baseFileName);
title('Original gray scale image');
drawnow;
% Intantiate 2-D histogram
hist2d = zeros(256, 256);
for col = 2 : columns-1
	fprintf('Processing column %d of original image...\n', col);
	for row = 2 : rows-1
		% Get 3x3 subimage, except for center pixel, so that's 8 pixels.
		subImage = [grayImage(row-1, col-1);...
			grayImage(row-1, col);...
			grayImage(row-1, col+1);...
			grayImage(row, col-1);...
			grayImage(row, col+1);...
			grayImage(row+1, col+1);...
			grayImage(row+1, col+1);...
			grayImage(row+1, col+1)];
		% Get the gray levle of the center pixel
		centerGL = grayImage(row, col) + 1; % Add 1 so we don't get 0 because matrices can't have a zeroeth row or column.
		for k = 1 : length(subImage)
			% Get this gray level
			thisGL = subImage(k) + 1; % Add 1 so we don't get 0 because matrices can't have a zeroeth row or column.
			% Increment the histogram at the column for the center pixel's gray level
			hist2d(centerGL, thisGL) = hist2d(centerGL, thisGL) + 1;
		end
	end
end
% Display the 2-D histogram
subplot(2, 1, 2);
imshow(hist2d, [], 'Colormap', jet(256));
title('2D histogram');
colorbar

Answer 2

[Harsha]

hists is an 2d-histogram of image with size M×N can be represented by a 2D gray level intensity function f(x, y). The value of f(x, y) is the gray level, ranging from 0 to L-1, where L is the number of distinct gray levels. In a 2D thresholding method, the gray level of a pixel and its local average gray level are both used. The local average gray level is also divided into the same L values, let g(x, y) be the function of the local average gray level.f(x,y) is an input image and g(x,y) average gray level of input image. total is the number of pixels in the input image(i.e M x N).

[Madhava Teja Munagala]

Hii harsha bro, im doing project on otsu 2d,,, im not understanding code,how to implement,plz help me

[Harsha]

Here I am attaching 2D otsu MATLAB code, I used this code in my research

https://in.mathworks.com/matlabcentral/fileexchange/66054-fully-automated-segmentation-of-lung-parenchyma

Complete Code will be available on October 2018. For more information please refer following publications: [1] S. Pramod Kumar and M. V. Latte, “Fully automated segmentation of lung parenchyma using break and repair strategy”, Journal of Intelligent Systems. De Gruyter (2017). DOI:10.1515/jisys-2017-0020

https://www.degruyter.com/view/j/jisys.ahead-of-print/jisys-2017-0020/jisys-2017-0020.xml?format=INT

[2] S. Pramod Kumar and M. V. Latte,“Modified and Optimized Method for Segmenting Pulmonary Parenchyma in CT Lung Images, Based on Fractional Calculus and Natural Selection”, Journal of Intelligent Systems. De Gruyter (2017). DOI:10.1515/jisys-2017-0028 https://www.degruyter.com/view/j/jisys.ahead-of-print/jisys-2017-0028/jisys-2017-0028.xml?format=INT

[Madhava Teja Munagala]

thank u harsha bro, i get output, how to analyze loop

[Image Analyst]

Madhava, I don't understand what you said.

Please elaborate after reading this link.

[MAS]

It works absolutely fine for two clusters. Any suggestions to update it for handling multi-level thresholding!?