octave_packages/image-1.0.15/graythresh.m

   1 ## Copyright (C) 2005  Barre-Piquot
   2 ##
   3 ## This program is free software; you can redistribute it and/or
   4 ## modify it under the terms of the GNU General Public License
   5 ## as published by the Free Software Foundation; either version 2
   6 ## of the License, or (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12
  13 ## -*- texinfo -*-
  14 ## @deftypefn {Function File} {@var{level}=} graythresh (@var{I})
  15 ## Compute global image threshold using Otsu's method.
  16 ##
  17 ## The output @var{level} is a global threshold (level) that can be used to convert
  18 ## an intensity image to a binary image with @code{im2bw}.
  19 ## @var{level} is a normalized intensity value that lies in the range [0, 1].
  20 ##
  21 ## The function uses Otsu's method, which chooses the threshold to
  22 ## minimize the intraclass variance of the black and white pixels.
  23 ##
  24 ## Color images are converted grayscale before @var{level} is computed.
  25 ## @seealso{im2bw}
  26 ## @end deftypefn
  27
  28 ## Note:
  29 ## This function is taken from
  30 ## http://www.irit.fr/recherches/SAMOVA/MEMBERS/JOLY/Homepage_files/IRR05/Barre-Piquot/graythresh.m
  31 ## I added texinfo documentation, error checking and sanitised the code.
  32 ##    -- Søren Hauberg
  33
  34 function level = graythresh (I)
  35     ## Input checking
  36     if (nargin != 1)
  37       print_usage();
  38     endif
  39     if (!isgray(I) && !isrgb(I))
  40       error("graythresh: input must be an image");
  41     endif
  42
  43     ## If the image is RGB convert it to grayscale
  44     if (isrgb(I))
  45       I = rgb2gray(I);
  46     endif
  47
  48     ## Calculation of the normalized histogram
  49     n = 256;
  50     h = hist(I(:), 1:n);
  51     h = h/(length(I(:))+1);
  52
  53     ## Calculation of the cumulated histogram and the mean values
  54     w = cumsum(h);
  55     mu = zeros(n, 1); mu(1) = h(1);
  56     for i=2:n
  57         mu(i) = mu(i-1) + i*h(i);
  58     end
  59
  60     ## Initialisation of the values used for the threshold calculation
  61     level = find (h > 0, 1);
  62     w0 = w(level);
  63     w1 = 1-w0;
  64     mu0 = mu(level)/w0;
  65     mu1 = (mu(end)-mu(level))/w1;
  66     max = w0*w1*(mu1-mu0)*(mu1-mu0);
  67
  68     ## For each step of the histogram, calculation of the threshold and storing of the maximum
  69     for i = find (h > 0)
  70         w0 = w(i);
  71         w1 = 1-w0;
  72         mu0 = mu(i)/w0;
  73         mu1 = (mu(end)-mu(i))/w1;
  74         s = w0*w1*(mu1-mu0)*(mu1-mu0);
  75         if (s > max)
  76             max = s;
  77             level = i;
  78         endif
  79     endfor
  80
  81     ## Normalisation of the threshold
  82     level /= n;
  83 endfunction
  84