MATLAB: Two-dimensional Otsu’s method

otsu

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

Best Answer

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).
Related Question