## Copyright (C) 2005  Barre-Piquot
##
## This program is free software; you can redistribute it and/or
## modify it under the terms of the GNU General Public License
## as published by the Free Software Foundation; either version 2
## of the License, or (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.

## -*- texinfo -*-
## @deftypefn {Function File} {@var{level}=} graythresh (@var{I})
## Compute global image threshold using Otsu's method.
##
## The output @var{level} is a global threshold (level) that can be used to convert
## an intensity image to a binary image with @code{im2bw}.
## @var{level} is a normalized intensity value that lies in the range [0, 1]. 
## 
## The function uses Otsu's method, which chooses the threshold to
## minimize the intraclass variance of the black and white pixels. 
## 
## Color images are converted grayscale before @var{level} is computed.
## @seealso{im2bw}
## @end deftypefn

## Note:
## This function is taken from
## http://www.irit.fr/recherches/SAMOVA/MEMBERS/JOLY/Homepage_files/IRR05/Barre-Piquot/graythresh.m
## I added texinfo documentation, error checking and sanitised the code.
##    -- Søren Hauberg

function level = graythresh (I)
    ## Input checking
    if (nargin != 1)
      print_usage();
    endif
    if (!isgray(I) && !isrgb(I))
      error("graythresh: input must be an image");
    endif
    
    ## If the image is RGB convert it to grayscale
    if (isrgb(I))
      I = rgb2gray(I);
    endif

    ## Calculation of the normalized histogram
    n = 256;
    h = hist(I(:), n);        
    h = h/(length(I(:))+1);
    
    ## Calculation of the cumulated histogram and the mean values
    w = cumsum(h);
    mu = zeros(n, 1); mu(1) = h(1);
    for i=2:n
        mu(i) = mu(i-1) + i*h(i);
    end    
         
    ## Initialisation of the values used for the threshold calculation
    w0 = w(1);
    w1 = 1-w0;
    mu0 = mu(1)/w0;
    mu1 = (mu(end)-mu(1))/w1;
    max = w0*w1*(mu1-mu0)*(mu1-mu0);
    level = 1;
    
    ## For each step of the histogram, calculation of the threshold and storing of the maximum
    for i = 2:n
        w0 = w(i);
        w1 = 1-w0;
        mu0 = mu(i)/w0;
        mu1 = (mu(end)-mu(i))/w1;
        s = w0*w1*(mu1-mu0)*(mu1-mu0);
        if (s > max)
            max = s;
            level = i;
        endif
    endfor
    
    ## Normalisation of the threshold        
    level /= n;
endfunction