--- a/main/image/NEWS	Thu Sep 27 12:51:44 2012 +0000
+++ b/main/image/NEWS	Thu Sep 27 13:18:13 2012 +0000
@@ -50,7 +50,7 @@
       minimum
       MinError
       moments
-      otsu
+      Otsu
       percentile
 
  ** The following functions have been deprecated (see their help text

--- a/main/image/inst/bwhitmiss.m	Thu Sep 27 12:51:44 2012 +0000
+++ b/main/image/inst/bwhitmiss.m	Thu Sep 27 13:18:13 2012 +0000
@@ -1,67 +1,67 @@
-## Copyright (C) 2008 Soren Hauberg <soren@hauberg.org>
-##
-## 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 3 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.
-##
-## You should have received a copy of the GNU General Public License along with
-## this program; if not, see <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn {Function File} @var{bw2} = bwhitmiss (@var{bw1}, @var{se1}, @var{se1})
-## @deftypefnx{Function File} @var{bw2} = bwhitmiss (@var{bw1}, @var{interval})
-## Perform the binary hit-miss operation.
-##
-## If two structuring elements @var{se1} and @var{se1} are given, the hit-miss
-## operation is defined as
-## @example
-## bw2 = erode(bw1, se1) & erode(!bw1, se2);
-## @end example
-## If instead an 'interval' array is given, two structuring elements are computed
-## as
-## @example
-## se1 = (interval ==  1)
-## se2 = (interval == -1)
-## @end example
-## and then the operation is defined as previously.
-## @seealso{bwmorph}
-## @end deftypefn
-
-function bw = bwhitmiss(im, varargin)
-  ## Checkinput
-  if (nargin != 2 && nargin != 3)
-    print_usage();
-  endif
-  if (!ismatrix(im) || !isreal(im))
-    error("bwhitmiss: first input argument must be a real matrix");
-  endif
-
-  ## Get structuring elements
-  if (nargin == 2) # bwhitmiss (im, interval)
-    interval = varargin{1};
-    if (!isreal(interval))
-      error("bwhitmiss: second input argument must be a real matrix");
-    endif
-    if (!all( (interval(:) == 1) | (interval(:) == 0) | (interval(:) == -1) ))
-      error("bwhitmiss: second input argument can only contain the values -1, 0, and 1");
-    endif
-    se1 = (interval ==  1);
-    se2 = (interval == -1);
-  else # bwhitmiss (im, se1, se2)
-    se1 = varargin{1};
-    se2 = varargin{2};
-    if (!all((se1(:) == 1) | (se1(:) == 0)) || !all((se2(:) == 1) | (se2(:) == 0)))
-      error("bwhitmiss: structuring elements can only contain zeros and ones.");
-    endif
-  endif
-  
-  ## Perform filtering
-  bw = erode(im, se1) & erode(!im, se2);
-
-endfunction
+## Copyright (C) 2008 Soren Hauberg <soren@hauberg.org>
+##
+## 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 3 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.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} @var{bw2} = bwhitmiss (@var{bw1}, @var{se1}, @var{se1})
+## @deftypefnx{Function File} @var{bw2} = bwhitmiss (@var{bw1}, @var{interval})
+## Perform the binary hit-miss operation.
+##
+## If two structuring elements @var{se1} and @var{se1} are given, the hit-miss
+## operation is defined as
+## @example
+## bw2 = erode(bw1, se1) & erode(!bw1, se2);
+## @end example
+## If instead an 'interval' array is given, two structuring elements are computed
+## as
+## @example
+## se1 = (interval ==  1)
+## se2 = (interval == -1)
+## @end example
+## and then the operation is defined as previously.
+## @seealso{bwmorph}
+## @end deftypefn
+
+function bw = bwhitmiss(im, varargin)
+  ## Checkinput
+  if (nargin != 2 && nargin != 3)
+    print_usage();
+  endif
+  if (!ismatrix(im) || !isreal(im))
+    error("bwhitmiss: first input argument must be a real matrix");
+  endif
+
+  ## Get structuring elements
+  if (nargin == 2) # bwhitmiss (im, interval)
+    interval = varargin{1};
+    if (!isreal(interval))
+      error("bwhitmiss: second input argument must be a real matrix");
+    endif
+    if (!all( (interval(:) == 1) | (interval(:) == 0) | (interval(:) == -1) ))
+      error("bwhitmiss: second input argument can only contain the values -1, 0, and 1");
+    endif
+    se1 = (interval ==  1);
+    se2 = (interval == -1);
+  else # bwhitmiss (im, se1, se2)
+    se1 = varargin{1};
+    se2 = varargin{2};
+    if (!all((se1(:) == 1) | (se1(:) == 0)) || !all((se2(:) == 1) | (se2(:) == 0)))
+      error("bwhitmiss: structuring elements can only contain zeros and ones.");
+    endif
+  endif
+  
+  ## Perform filtering
+  bw = erode(im, se1) & erode(!im, se2);
+
+endfunction

--- a/main/image/inst/graythresh.m	Thu Sep 27 12:51:44 2012 +0000
+++ b/main/image/inst/graythresh.m	Thu Sep 27 13:18:13 2012 +0000
@@ -1,726 +1,728 @@
-## Copyright (C) 2004 Antti Niemistö <antti.niemisto@tut.fi>
-## Copyright (C) 2005 Barre-Piquot
-## Copyright (C) 2007 Søren Hauberg
-## Copyright (C) 2012 Carnë Draug <carandraug+dev@gmail.com>
-##
-## 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 3 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.
-##
-## You should have received a copy of the GNU General Public License along with
-## this program; if not, see <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn {Function File} {[@var{level}, @var{sep}] =} graythresh (@var{img})
-## @deftypefnx {Function File} {[@var{level}, @var{sep}] =} graythresh (@var{img}, @var{algo})
-## @deftypefnx {Function File} {[@var{level}] =} graythresh (@var{img}, "percentile", @var{fraction})
-## Compute global image threshold.
-##
-## Given an image @var{img} finds the optimal threshold value @var{level} for
-## conversion to a binary image with @code{im2bw}. Color images are converted
-## to grayscale before @var{level} is computed.
-##
-## The optional argument @var{method} is the algorithm to be used (default's to
-## Otsu). The available algorithms are:
-##
-## @table @asis
-## @item Otsu (default)
-## Implements Otsu's method as described in @cite{Nobuyuki Otsu (1979). "A
-## threshold selection method from gray-level histograms", IEEE Trans. Sys.,
-## Man., Cyber. 9 (1): 62-66}. This algorithm chooses the threshold to minimize
-## the intraclass variance of the black and white pixels. The second output,
-## @var{sep} represents the ``goodness'' (or separability) of the threshold at
-## @var{level}.
-##
-## @item concavity
-## Find a global threshold for a grayscale image by choosing the threshold to
-## be in the shoulder of the histogram @cite{A. Rosenfeld, and P. De La Torre
-## (1983). "Histogram concavity analysis as an aid in threshold selection", IEEE
-## Transactions on Systems, Man, and Cybernetics, 13: 231-235}.
-##
-## @item Huang
-## Not yet implemented.
-##
-## @item ImageJ
-## A variation of the intermeans algorithm, this is the default for ImageJ.
-## Not yet implemented.
-##
-## @item intermodes
-## This assumes a bimodal histogram and chooses the threshold to be the mean of
-## the two peaks of the bimodal histogram @cite{J. M. S. Prewitt, and M. L.
-## Mendelsohn (1966). "The analysis of cell images", Annals of the New York
-## Academy of Sciences, 128: 1035-1053}.
-##
-## Images with histograms having extremely unequal peaks or a broad and flat
-## valley are unsuitable for this method.
-##
-## @item intermeans
-## Iterative procedure based on the iterative intermeans algorithm of @cite{T.
-## Ridler, and S. Calvard (1978). "Picture thresholding using an iterative
-## selection method", IEEE Transactions on Systems, Man, and Cybernetics, 8: 630-632}
-## and @cite{H. J. Trussell (1979). "Comments on 'Picture thresholding using an
-## iterative selection method'", IEEE Transactions on Systems, Man, and Cybernetics,
-## 9: 311}.
-##
-## Note that several implementations of this method exist. See the source code
-## for details.
-##
-## @item Li
-## Not yet implemented.
-##
-## @item MaxEntropy
-## Implements Kapur-Sahoo-Wong (Maximum Entropy) thresholding method based on the
-## entropy of the image histogram @cite{J. N. Kapur, P. K. Sahoo, and A. C. K. Wong
-## (1985). "A new method for gray-level picture thresholding using the entropy
-## of the histogram", Graphical Models and Image Processing, 29(3): 273-285}.
-##
-## @item MaxLikelihood
-## Find a global threshold for a grayscale image using the maximum likelihood
-## via expectation maximization method @cite{A. P. Dempster, N. M. Laird, and D. B.
-## Rubin (1977). "Maximum likelihood from incomplete data via the EM algorithm",
-## Journal of the Royal Statistical Society, Series B, 39:1-38}.
-##
-## @item mean
-## The mean intensity value. It is mostly used by other methods as a first guess
-## threshold.
-##
-## @item MinError
-## An iterative implementation of Kittler and Illingworth's Minimum Error
-## thresholding @cite{J. Kittler, and J. Illingworth (1986). "Minimum error
-## thresholding", Pattern recognition, 19: 41-47}.
-##
-## This implementation seems to converge more often than the original.
-## Nevertheless, sometimes the algorithm does not converge to a solution. In
-## that case a warning is displayed and defaults to the initial estimate of the
-## mean method.
-##
-## @item minimum
-## This assumes a bimodal histogram and chooses the threshold to be in the
-## valley of the bimodal histogram.  This method is also known as the mode
-## method @cite{J. M. S. Prewitt, and M. L. Mendelsohn (1966). "The analysis of
-## cell images", Annals of the New York Academy of Sciences, 128: 1035-1053}.
-##
-## Images with histograms having extremely unequal peaks or a broad and flat
-## valley are unsuitable for this method.
-##
-## @item moments
-## Find a global threshold for a grayscale image using moment preserving
-## thresholding method @cite{W. Tsai (1985). "Moment-preserving thresholding:
-## a new approach", Computer Vision, Graphics, and Image Processing, 29: 377-393}
-##
-## @item percentile
-## Assumes a specific @var{fraction} of pixels to be background.  If no value is
-## given, assumes a value of 0.5 (equal distribution of background and foreground)
-## @cite{W Doyle (1962). "Operation useful for similarity-invariant pattern
-## recognition", Journal of the Association for Computing Machinery 9: 259-267}
-##
-## @item RenyiEntropy
-## Not yet implemented.
-##
-## @item Shanbhag
-## Not yet implemented.
-##
-## @item triangle
-## Not yet implemented.
-##
-## @item Yen
-## Not yet implemented.
-## @end table
-##
-## @seealso{im2bw}
-## @end deftypefn
-
-## Notes:
-##  * Otsu's method is a function taken from
-##    http://www.irit.fr/~Philippe.Joly/Teaching/M2IRR/IRR05/ Søren Hauberg
-##    added texinfo documentation, error checking and sanitised the code.
-##  * The following methods were adapted from http://www.cs.tut.fi/~ant/histthresh/
-##      intermodes    percentile      minimum
-##      MaxEntropy    MaxLikelihood   intermeans
-##      moments       minerror        concavity
-
-
-function [level, good] = graythresh (img, algo = "otsu", varargin)
-  ## Input checking
-  if (nargin < 1 || nargin > 3)
-    print_usage();
-  elseif (!isgray (img) && !isrgb (img))
-    error ("graythresh: input must be an image");
-  endif
-
-  ## If the image is RGB convert it to grayscale
-  if (isrgb (img))
-    img = rgb2gray (img);
-  endif
-
-  switch tolower (algo)
-    case {"concavity"},     [level]       = concavity (img);
-    case {"intermeans"},    [level]       = intermeans (img);
-    case {"intermodes"},    [level]       = intermodes (img);
-    case {"maxlikelihood"}, [level]       = maxlikelihood (img);
-    case {"maxentropy"},    [level]       = maxentropy (img);
-    case {"mean"},          [level]       = mean (img(:));
-    case {"minimum"},       [level]       = minimum (img);
-    case {"minerror"},      [level]       = minerror_iter (img);
-    case {"moments"},       [level]       = moments (img);
-    case {"otsu"},          [level, good] = otsu (img);
-    case {"percentile"},    [level]       = percentile (img, varargin{1});
-    otherwise, error ("graythresh: unknown method '%s'", algo);
-  endswitch
-endfunction
-
-function [level, good] = otsu (img)
-  ## Calculation of the normalized histogram
-  n = double (intmax (class (img))) + 1;
-  h = hist (img(:), 1:n);
-  h = h / (length (img(:)) + 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);
-  endfor
-
-  ## Initialisation of the values used for the threshold calculation
-  level = find (h > 0, 1);
-  w0  = w(level);
-  w1  = 1 - w0;
-  mu0 = mu(level) / w0;
-  mu1 = (mu(end) - mu(level)) / w1;
-  good = w0 * w1 * (mu1 - mu0) * (mu1 - mu0);
-
-  ## For each step of the histogram, calculation of the threshold
-  ## and storing of the maximum
-  for i = find (h > 0)
-    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 > good)
-      good  = s;
-      level = i;
-    endif
-  endfor
-
-  ## Normalisation of the threshold
-  level /= n;
-endfunction
-
-#{
-  ##The following is another implementation of Otsu's method from the HistThresh
-  ##toolbox
-  if nargin == 1
-    n = 255;
-  end
-
-  % This algorithm is implemented in the Image Processing Toolbox.
-  %I = uint8(I);
-  %T = n*graythresh(I);
-
-  % The implementation below uses the notations from the paper, but is
-  % significantly slower.
-
-  I = double(I);
-
-  % Calculate the histogram.
-  y = hist(I(:),0:n);
-
-  warning off MATLAB:divideByZero
-  for j = 0:n
-    mu = partial_sumB(y,j)/partial_sumA(y,j);
-    nu = (partial_sumB(y,n)-partial_sumB(y,j))/(partial_sumA(y,n)-partial_sumA(y,j));
-    vec(j+1) = partial_sumA(y,j)*(partial_sumA(y,n)-partial_sumA(y,j))*(mu-nu)^2;
-  end
-  warning on MATLAB:divideByZero
-
-  [maximum,ind] = max(vec);
-  T = ind-1;
-#}
-
-function level = moments (I, n)
-
-  if nargin == 1
-    n = 255;
-  end
-
-  I = double(I);
-
-  % Calculate the histogram.
-  y = hist (I(:), 0:n);
-
-  % The threshold is chosen such that partial_sumA(y,t)/partial_sumA(y,n) is closest to x0.
-  Avec = zeros (1, n+1);
-
-  for t = 0:n
-    Avec(t+1) = partial_sumA (y, t) / partial_sumA (y, n);
-  end
-
-  % The following finds x0.
-  x2 = (partial_sumB(y,n)*partial_sumC(y,n) - partial_sumA(y,n)*partial_sumD(y,n)) / (partial_sumA(y,n)*partial_sumC(y,n) - partial_sumB(y,n)^2);
-  x1 = (partial_sumB(y,n)*partial_sumD(y,n) - partial_sumC(y,n)^2) / (partial_sumA (y,n)*partial_sumC(y,n) - partial_sumB(y,n)^2);
-  x0 = .5 - (partial_sumB(y,n)/partial_sumA(y,n) + x2/2) / sqrt (x2^2 - 4*x1);
-
-  % And finally the threshold.
-  [minimum, ind] = min (abs (Avec-x0));
-  level = ind-1;
-
-  ## Normalisation of the threshold
-  level /= n;
-endfunction
-
-function T = maxentropy(I,n)
-  if nargin == 1
-    n = 255;
-  end
-
-  I = double(I);
-
-  % Calculate the histogram.
-  y = hist(I(:),0:n);
-
-  warning off
-  % The threshold is chosen such that the following expression is minimized.
-  for j = 0:n
-    vec(j+1) = negativeE(y,j)/partial_sumA(y,j) - log10(partial_sumA(y,j)) + ...
-        (negativeE(y,n)-negativeE(y,j))/(partial_sumA(y,n)-partial_sumA(y,j)) - log10(partial_sumA(y,n)-partial_sumA(y,j));
-  end
-  warning on
-
-  [minimum,ind] = min(vec);
-  T = ind-1;
-endfunction
-
-function [T] = intermodes (I,n)
-  if nargin == 1
-    n = 255;
-  end
-
-  I = double(I);
-
-  % Calculate the histogram.
-  y = hist(I(:),0:n);
-
-  % Smooth the histogram by iterative three point mean filtering.
-  iter = 0;
-  while ~bimodtest(y)
-    h = ones(1,3)/3;
-    y = conv2(y,h,'same');
-    iter = iter+1;
-    % If the histogram turns out not to be bimodal, set T to zero.
-    if iter > 10000;
-      T = 0;
-      return
-    end
-  end
-
-  % The threshold is the mean of the two peaks of the histogram.
-  ind = 0;
-  for k = 2:n
-    if y(k-1) < y(k) && y(k+1) < y(k)
-      ind = ind+1;
-      TT(ind) = k-1;
-    end
-  end
-  T = floor(mean(TT));
-endfunction
-
-function [T] = percentile (I, p, n)
-  if nargin == 1
-    p = 0.5;
-    n = 255;
-  elseif nargin == 2
-    n = 255;
-  end
-
-  I = double(I);
-
-  % Calculate the histogram.
-  y = hist(I(:),0:n);
-
-  % The threshold is chosen such that 50% of pixels lie in each category.
-  Avec = zeros(1,n+1);
-  for t = 0:n
-    Avec(t+1) = partial_sumA(y,t)/partial_sumA(y,n);
-  end
-
-  [minimum,ind] = min(abs(Avec-p));
-  T = ind-1;
-endfunction
-
-
-function T = minimum(I,n);
-  if nargin == 1
-    n = 255;
-  end
-
-  I = double(I);
-
-  % Calculate the histogram.
-  y = hist(I(:),0:n);
-
-  % Smooth the histogram by iterative three point mean filtering.
-  iter = 0;
-  while ~bimodtest(y)
-    h = ones(1,3)/3;
-    y = conv2(y,h,'same');
-    iter = iter+1;
-    % If the histogram turns out not to be bimodal, set T to zero.
-    if iter > 10000;
-      T = 0;
-      return
-    end
-  end
-
-  % The threshold is the minimum between the two peaks.
-  for k = 2:n
-    if y(k-1) > y(k) && y(k+1) > y(k)
-      T = k-1;
-    end
-  end
-endfunction
-
-function [T] = minerror_iter (I,n)
-  if nargin == 1
-    n = 255;
-  end
-
-  I = double(I);
-
-  % Calculate the histogram.
-  y = hist(I(:),0:n);
-
-  % The initial estimate for the threshold is found with the MEAN algorithm.
-  T = floor (mean (I, img(:)));
-  Tprev = NaN;
-
-  while T ~= Tprev
-
-    % Calculate some statistics.
-    mu = partial_sumB(y,T)/partial_sumA(y,T);
-    nu = (partial_sumB(y,n)-partial_sumB(y,T))/(partial_sumA(y,n)-partial_sumA(y,T));
-    p = partial_sumA(y,T)/partial_sumA(y,n);
-    q = (partial_sumA(y,n)-partial_sumA(y,T)) / partial_sumA(y,n);
-    sigma2 = partial_sumC(y,T)/partial_sumA(y,T)-mu^2;
-    tau2 = (partial_sumC(y,n)-partial_sumC(y,T)) / (partial_sumA(y,n)-partial_sumA(y,T)) - nu^2;
-
-    % The terms of the quadratic equation to be solved.
-    w0 = 1/sigma2-1/tau2;
-    w1 = mu/sigma2-nu/tau2;
-    w2 = mu^2/sigma2 - nu^2/tau2 + log10((sigma2*q^2)/(tau2*p^2));
-
-    % If the next threshold would be imaginary, return with the current one.
-    sqterm = w1^2-w0*w2;
-    if sqterm < 0
-      warning('MINERROR:NaN','Warning: th_minerror_iter did not converge.')
-      return
-    end
-
-    % The updated threshold is the integer part of the solution of the
-    % quadratic equation.
-    Tprev = T;
-    T = floor((w1+sqrt(sqterm))/w0);
-
-    % If the threshold turns out to be NaN, return with the previous threshold.
-    if isnan(T)
-      warning('MINERROR:NaN','Warning: th_minerror_iter did not converge.')
-      T = Tprev;
-    end
-
-  end
-endfunction
-#{
-  ## this was not an implementatino of the original minerror algorithm but seems
-  ## to converg more often than the original. The original is (also from the
-  ## HistThresh toolbox
-  function T = th_minerror(I,n)
-    if nargin == 1
-      n = 255;
-    end
-
-    I = double(I);
-
-    % Calculate the histogram.
-    y = hist(I(:),0:n);
-
-    warning off
-    % The threshold is chosen such that the following expression is minimized.
-    for j = 0:n
-      mu = partial_sumB(y,j)/partial_sumA(y,j);
-      nu = (partial_sumB(y,n)-partial_sumB(y,j))/(partial_sumA(y,n)-partial_sumA(y,j));
-      p = partial_sumA(y,j)/partial_sumA(y,n);
-      q = (partial_sumA(y,n)-partial_sumA(y,j)) / partial_sumA(y,n);
-      sigma2 = partial_sumC(y,j)/partial_sumA(y,j)-mu^2;
-      tau2 = (partial_sumC(y,n)-partial_sumC(y,j)) / (partial_sumA(y,n)-partial_sumA(y,j)) - nu^2;
-      vec(j+1) = p*log10(sqrt(sigma2)/p) + q*log10(sqrt(tau2)/q);
-    end
-    warning on
-
-    vec(vec==-inf) = NaN;
-    [minimum,ind] = min(vec);
-    T = ind-1;
-  endfunction
-#}
-
-function T = maxlikelihood (I,n)
-  if nargin == 1
-    n = 255;
-  end
-
-  I = double(I);
-
-  % Calculate the histogram.
-  y = hist(I(:),0:n);
-
-  % The initial estimate for the threshold is found with the MINIMUM
-  % algorithm.
-  T = th_minimum(I,n);
-
-  % Calculate initial values for the statistics.
-  mu = partial_sumB(y,T)/partial_sumA(y,T);
-  nu = (partial_sumB(y,n)-partial_sumB(y,T))/(partial_sumA(y,n)-partial_sumA(y,T));
-  p = partial_sumA(y,T)/partial_sumA(y,n);
-  q = (partial_sumA(y,n)-partial_sumA(y,T)) / partial_sumA(y,n);
-  sigma2 = partial_sumC(y,T)/partial_sumA(y,T)-mu^2;
-  tau2 = (partial_sumC(y,n)-partial_sumC(y,T)) / (partial_sumA(y,n)-partial_sumA(y,T)) - nu^2;
-
-  mu_prev = NaN;
-  nu_prev = NaN;
-  p_prev = NaN;
-  q_prev = NaN;
-  sigma2_prev = NaN;
-  tau2_prev = NaN;
-
-  while abs(mu-mu_prev) > eps || abs(nu-nu_prev) > eps || ...
-        abs(p-p_prev) > eps || abs(q-q_prev) > eps || ...
-        abs(sigma2-sigma2_prev) > eps || abs(tau2-tau2_prev) > eps
-    for i = 0:n
-      phi(i+1) = p/q * exp(-((i-mu)^2) / (2*sigma2)) / ...
-          (p/sqrt(sigma2) * exp(-((i-mu)^2) / (2*sigma2)) + ...
-           (q/sqrt(tau2)) * exp(-((i-nu)^2) / (2*tau2)));
-    end
-    ind = 0:n;
-    gamma = 1-phi;
-    F = phi*y';
-    G = gamma*y';
-    p_prev = p;
-    q_prev = q;
-    mu_prev = mu;
-    nu_prev = nu;
-    sigma2_prev = nu;
-    tau2_prev = nu;
-    p = F/partial_sumA(y,n);
-    q = G/partial_sumA(y,n);
-    mu = ind.*phi*y'/F;
-    nu = ind.*gamma*y'/G;
-    sigma2 = ind.^2.*phi*y'/F - mu^2;
-    tau2 = ind.^2.*gamma*y'/G - nu^2;
-  end
-
-  % The terms of the quadratic equation to be solved.
-  w0 = 1/sigma2-1/tau2;
-  w1 = mu/sigma2-nu/tau2;
-  w2 = mu^2/sigma2 - nu^2/tau2 + log10((sigma2*q^2)/(tau2*p^2));
-
-  % If the threshold would be imaginary, return with threshold set to zero.
-  sqterm = w1^2-w0*w2;
-  if sqterm < 0;
-    T = 0;
-    return
-  end
-
-  % The threshold is the integer part of the solution of the quadratic
-  % equation.
-  T = floor((w1+sqrt(sqterm))/w0);
-endfunction
-
-function T = intermeans (I,n)
-  if nargin == 1
-    n = 255;
-  end
-
-  I = double(I);
-
-  % Calculate the histogram.
-  y = hist(I(:),0:n);
-
-  % The initial estimate for the threshold is found with the MEAN algorithm.
-  T = floor (mean (I, img(:)));
-  Tprev = NaN;
-
-  % The threshold is found iteratively. In each iteration, the means of the
-  % pixels below (mu) the threshold and above (nu) it are found. The
-  % updated threshold is the mean of mu and nu.
-  while T ~= Tprev
-    mu = partial_sumB(y,T)/partial_sumA(y,T);
-    nu = (partial_sumB(y,n)-partial_sumB(y,T))/(partial_sumA(y,n)-partial_sumA(y,T));
-    Tprev = T;
-    T = floor((mu+nu)/2);
-  end
-endfunction
-
-function T = concavity (I,n)
-  if nargin == 1
-    n = 255;
-  end
-
-  I = double(I);
-
-  % Calculate the histogram and its convex hull.
-  h = hist(I(:),0:n);
-  H = hconvhull(h);
-
-  % Find the local maxima of the difference H-h.
-  lmax = flocmax(H-h);
-
-  % Find the histogram balance around each index.
-  for k = 0:n
-    E(k+1) = hbalance(h,k);
-  end
-
-  % The threshold is the local maximum with highest balance.
-  E = E.*lmax;
-  [dummy ind] = max(E);
-  T = ind-1;
-endfunction
-
-################################################################################
-## Auxiliary functions from HistThresh toolbox http://www.cs.tut.fi/~ant/histthresh/
-################################################################################
-
-## partial sums from C. A. Glasbey, "An analysis of histogram-based thresholding
-## algorithms," CVGIP: Graphical Models and Image Processing, vol. 55, pp. 532-537, 1993.
-function x = partial_sumA (y, j)
-  x = sum (y(1:j+1));
-endfunction
-function x = partial_sumB (y, j)
-  ind = 0:j;
-  x   = ind*y(1:j+1)';
-endfunction
-function x = partial_sumC (y, j)
-  ind = 0:j;
-  x = ind.^2*y(1:j+1)';
-endfunction
-function x = partial_sumD (y, j)
-  ind = 0:j;
-  x = ind.^3*y(1:j+1)';
-endfunction
-
-## Test if a histogram is bimodal.
-function b = bimodtest(y)
-  len = length(y);
-  b = false;
-  modes = 0;
-
-  % Count the number of modes of the histogram in a loop. If the number
-  % exceeds 2, return with boolean return value false.
-  for k = 2:len-1
-    if y(k-1) < y(k) && y(k+1) < y(k)
-      modes = modes+1;
-      if modes > 2
-        return
-      end
-    end
-  end
-
-  % The number of modes could be less than two here
-  if modes == 2
-    b = true;
-  end
-endfunction
-
-## Find the local maxima of a vector using a three point neighborhood.
-function y = flocmax(x)
-%  y    binary vector with maxima of x marked as ones
-
-  len = length(x);
-  y = zeros(1,len);
-
-  for k = 2:len-1
-    [dummy,ind] = max(x(k-1:k+1));
-    if ind == 2
-      y(k) = 1;
-    end
-  end
-endfunction
-
-## Calculate the balance measure of the histogram around a histogram index.
-function E = hbalance(y,ind)
-%  y    histogram
-%  ind  index about which balance is calculated
-%
-% Out:
-%  E    balance measure
-%
-% References:
-%
-% A. Rosenfeld and P. De La Torre, "Histogram concavity analysis as an aid
-% in threhold selection," IEEE Transactions on Systems, Man, and
-% Cybernetics, vol. 13, pp. 231-235, 1983.
-%
-% P. K. Sahoo, S. Soltani, and A. K. C. Wong, "A survey of thresholding
-% techniques," Computer Vision, Graphics, and Image Processing, vol. 41,
-% pp. 233-260, 1988.
-
-  n = length(y)-1;
-  E = partial_sumA(y,ind)*(partial_sumA(y,n)-partial_sumA(y,ind));
-endfunction
-
-## Find the convex hull of a histogram.
-function H = hconvhull(h)
-  % In:
-  %  h    histogram
-  %
-  % Out:
-  %  H    convex hull of histogram
-  %
-  % References:
-  %
-  % A. Rosenfeld and P. De La Torre, "Histogram concavity analysis as an aid
-  % in threhold selection," IEEE Transactions on Systems, Man, and
-  % Cybernetics, vol. 13, pp. 231-235, 1983.
-
-  len = length(h);
-  K(1) = 1;
-  k = 1;
-
-  % The vector K gives the locations of the vertices of the convex hull.
-  while K(k)~=len
-
-    theta = zeros(1,len-K(k));
-    for i = K(k)+1:len
-      x = i-K(k);
-      y = h(i)-h(K(k));
-      theta(i-K(k)) = atan2(y,x);
-    end
-
-    maximum = max(theta);
-    maxloc = find(theta==maximum);
-    k = k+1;
-    K(k) = maxloc(end)+K(k-1);
-
-  end
-
-  % Form the convex hull.
-  H = zeros(1,len);
-  for i = 2:length(K)
-    H(K(i-1):K(i)) = h(K(i-1))+(h(K(i))-h(K(i-1)))/(K(i)-K(i-1))*(0:K(i)-K(i-1));
-  end
-endfunction
-
-## Entroy function. Note that the function returns the negative of entropy.
-function x = negativeE(y,j)
-  ## used by the maxentropy method only
-  y = y(1:j+1);
-  y = y(y~=0);
-  x = sum(y.*log10(y));
-endfunction
+## Copyright (C) 2004 Antti Niemistö <antti.niemisto@tut.fi>
+## Copyright (C) 2005 Barre-Piquot
+## Copyright (C) 2007 Søren Hauberg
+## Copyright (C) 2012 Carnë Draug <carandraug+dev@gmail.com>
+##
+## 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 3 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.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{level}, @var{sep}] =} graythresh (@var{img})
+## @deftypefnx {Function File} {[@var{level}, @var{sep}] =} graythresh (@var{img}, @var{algo})
+## @deftypefnx {Function File} {[@var{level}] =} graythresh (@var{img}, "percentile", @var{fraction})
+## Compute global image threshold.
+##
+## Given an image @var{img} finds the optimal threshold value @var{level} for
+## conversion to a binary image with @code{im2bw}. Color images are converted
+## to grayscale before @var{level} is computed.
+##
+## The optional argument @var{method} is the algorithm to be used (default's to
+## Otsu). The available algorithms are:
+##
+## @table @asis
+## @item Otsu (default)
+## Implements Otsu's method as described in @cite{Nobuyuki Otsu (1979). "A
+## threshold selection method from gray-level histograms", IEEE Trans. Sys.,
+## Man., Cyber. 9 (1): 62-66}. This algorithm chooses the threshold to minimize
+## the intraclass variance of the black and white pixels. The second output,
+## @var{sep} represents the ``goodness'' (or separability) of the threshold at
+## @var{level}.
+##
+## @item concavity
+## Find a global threshold for a grayscale image by choosing the threshold to
+## be in the shoulder of the histogram @cite{A. Rosenfeld, and P. De La Torre
+## (1983). "Histogram concavity analysis as an aid in threshold selection", IEEE
+## Transactions on Systems, Man, and Cybernetics, 13: 231-235}.
+##
+## @item Huang
+## Not yet implemented.
+##
+## @item ImageJ
+## A variation of the intermeans algorithm, this is the default for ImageJ.
+## Not yet implemented.
+##
+## @item intermodes
+## This assumes a bimodal histogram and chooses the threshold to be the mean of
+## the two peaks of the bimodal histogram @cite{J. M. S. Prewitt, and M. L.
+## Mendelsohn (1966). "The analysis of cell images", Annals of the New York
+## Academy of Sciences, 128: 1035-1053}.
+##
+## Images with histograms having extremely unequal peaks or a broad and flat
+## valley are unsuitable for this method.
+##
+## @item intermeans
+## Iterative procedure based on the iterative intermeans algorithm of @cite{T.
+## Ridler, and S. Calvard (1978). "Picture thresholding using an iterative
+## selection method", IEEE Transactions on Systems, Man, and Cybernetics, 8: 630-632}
+## and @cite{H. J. Trussell (1979). "Comments on 'Picture thresholding using an
+## iterative selection method'", IEEE Transactions on Systems, Man, and Cybernetics,
+## 9: 311}.
+##
+## Note that several implementations of this method exist. See the source code
+## for details.
+##
+## @item Li
+## Not yet implemented.
+##
+## @item MaxEntropy
+## Implements Kapur-Sahoo-Wong (Maximum Entropy) thresholding method based on the
+## entropy of the image histogram @cite{J. N. Kapur, P. K. Sahoo, and A. C. K. Wong
+## (1985). "A new method for gray-level picture thresholding using the entropy
+## of the histogram", Graphical Models and Image Processing, 29(3): 273-285}.
+##
+## @item MaxLikelihood
+## Find a global threshold for a grayscale image using the maximum likelihood
+## via expectation maximization method @cite{A. P. Dempster, N. M. Laird, and D. B.
+## Rubin (1977). "Maximum likelihood from incomplete data via the EM algorithm",
+## Journal of the Royal Statistical Society, Series B, 39:1-38}.
+##
+## @item mean
+## The mean intensity value. It is mostly used by other methods as a first guess
+## threshold.
+##
+## @item MinError
+## An iterative implementation of Kittler and Illingworth's Minimum Error
+## thresholding @cite{J. Kittler, and J. Illingworth (1986). "Minimum error
+## thresholding", Pattern recognition, 19: 41-47}.
+##
+## This implementation seems to converge more often than the original.
+## Nevertheless, sometimes the algorithm does not converge to a solution. In
+## that case a warning is displayed and defaults to the initial estimate of the
+## mean method.
+##
+## @item minimum
+## This assumes a bimodal histogram and chooses the threshold to be in the
+## valley of the bimodal histogram.  This method is also known as the mode
+## method @cite{J. M. S. Prewitt, and M. L. Mendelsohn (1966). "The analysis of
+## cell images", Annals of the New York Academy of Sciences, 128: 1035-1053}.
+##
+## Images with histograms having extremely unequal peaks or a broad and flat
+## valley are unsuitable for this method.
+##
+## @item moments
+## Find a global threshold for a grayscale image using moment preserving
+## thresholding method @cite{W. Tsai (1985). "Moment-preserving thresholding:
+## a new approach", Computer Vision, Graphics, and Image Processing, 29: 377-393}
+##
+## @item percentile
+## Assumes a specific @var{fraction} of pixels to be background.  If no value is
+## given, assumes a value of 0.5 (equal distribution of background and foreground)
+## @cite{W Doyle (1962). "Operation useful for similarity-invariant pattern
+## recognition", Journal of the Association for Computing Machinery 9: 259-267}
+##
+## @item RenyiEntropy
+## Not yet implemented.
+##
+## @item Shanbhag
+## Not yet implemented.
+##
+## @item triangle
+## Not yet implemented.
+##
+## @item Yen
+## Not yet implemented.
+## @end table
+##
+## @seealso{im2bw}
+## @end deftypefn
+
+## Notes:
+##  * Otsu's method is a function taken from
+##    http://www.irit.fr/~Philippe.Joly/Teaching/M2IRR/IRR05/ Søren Hauberg
+##    added texinfo documentation, error checking and sanitised the code.
+##  * The following methods were adapted from http://www.cs.tut.fi/~ant/histthresh/
+##      intermodes    percentile      minimum
+##      MaxEntropy    MaxLikelihood   intermeans
+##      moments       minerror        concavity
+
+
+function [level, good] = graythresh (img, algo = "otsu", varargin)
+  ## Input checking
+  if (nargin < 1 || nargin > 3)
+    print_usage();
+  elseif (nargin > 2 && !any (strcmpi (varargin{1}, {"percentile"})))
+    error ("graythresh: algorithm `%s' does not accept any options.", algo);
+  elseif (!isgray (img) && !isrgb (img))
+    error ("graythresh: input must be an image");
+  endif
+
+  ## If the image is RGB convert it to grayscale
+  if (isrgb (img))
+    img = rgb2gray (img);
+  endif
+
+  switch tolower (algo)
+    case {"concavity"},     [level]       = concavity (img);
+    case {"intermeans"},    [level]       = intermeans (img);
+    case {"intermodes"},    [level]       = intermodes (img);
+    case {"maxlikelihood"}, [level]       = maxlikelihood (img);
+    case {"maxentropy"},    [level]       = maxentropy (img);
+    case {"mean"},          [level]       = mean (img(:));
+    case {"minimum"},       [level]       = minimum (img);
+    case {"minerror"},      [level]       = minerror_iter (img);
+    case {"moments"},       [level]       = moments (img);
+    case {"otsu"},          [level, good] = otsu (img);
+    case {"percentile"},    [level]       = percentile (img, varargin{1});
+    otherwise, error ("graythresh: unknown method '%s'", algo);
+  endswitch
+endfunction
+
+function [level, good] = otsu (img)
+  ## Calculation of the normalized histogram
+  n = double (intmax (class (img))) + 1;
+  h = hist (img(:), 1:n);
+  h = h / (length (img(:)) + 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);
+  endfor
+
+  ## Initialisation of the values used for the threshold calculation
+  level = find (h > 0, 1);
+  w0  = w(level);
+  w1  = 1 - w0;
+  mu0 = mu(level) / w0;
+  mu1 = (mu(end) - mu(level)) / w1;
+  good = w0 * w1 * (mu1 - mu0) * (mu1 - mu0);
+
+  ## For each step of the histogram, calculation of the threshold
+  ## and storing of the maximum
+  for i = find (h > 0)
+    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 > good)
+      good  = s;
+      level = i;
+    endif
+  endfor
+
+  ## Normalisation of the threshold
+  level /= n;
+endfunction
+
+#{
+  ##The following is another implementation of Otsu's method from the HistThresh
+  ##toolbox
+  if nargin == 1
+    n = 255;
+  end
+
+  % This algorithm is implemented in the Image Processing Toolbox.
+  %I = uint8(I);
+  %T = n*graythresh(I);
+
+  % The implementation below uses the notations from the paper, but is
+  % significantly slower.
+
+  I = double(I);
+
+  % Calculate the histogram.
+  y = hist(I(:),0:n);
+
+  warning off MATLAB:divideByZero
+  for j = 0:n
+    mu = partial_sumB(y,j)/partial_sumA(y,j);
+    nu = (partial_sumB(y,n)-partial_sumB(y,j))/(partial_sumA(y,n)-partial_sumA(y,j));
+    vec(j+1) = partial_sumA(y,j)*(partial_sumA(y,n)-partial_sumA(y,j))*(mu-nu)^2;
+  end
+  warning on MATLAB:divideByZero
+
+  [maximum,ind] = max(vec);
+  T = ind-1;
+#}
+
+function level = moments (I, n)
+
+  if nargin == 1
+    n = 255;
+  end
+
+  I = double(I);
+
+  % Calculate the histogram.
+  y = hist (I(:), 0:n);
+
+  % The threshold is chosen such that partial_sumA(y,t)/partial_sumA(y,n) is closest to x0.
+  Avec = zeros (1, n+1);
+
+  for t = 0:n
+    Avec(t+1) = partial_sumA (y, t) / partial_sumA (y, n);
+  end
+
+  % The following finds x0.
+  x2 = (partial_sumB(y,n)*partial_sumC(y,n) - partial_sumA(y,n)*partial_sumD(y,n)) / (partial_sumA(y,n)*partial_sumC(y,n) - partial_sumB(y,n)^2);
+  x1 = (partial_sumB(y,n)*partial_sumD(y,n) - partial_sumC(y,n)^2) / (partial_sumA (y,n)*partial_sumC(y,n) - partial_sumB(y,n)^2);
+  x0 = .5 - (partial_sumB(y,n)/partial_sumA(y,n) + x2/2) / sqrt (x2^2 - 4*x1);
+
+  % And finally the threshold.
+  [minimum, ind] = min (abs (Avec-x0));
+  level = ind-1;
+
+  ## Normalisation of the threshold
+  level /= n;
+endfunction
+
+function T = maxentropy(I,n)
+  if nargin == 1
+    n = 255;
+  end
+
+  I = double(I);
+
+  % Calculate the histogram.
+  y = hist(I(:),0:n);
+
+  warning off
+  % The threshold is chosen such that the following expression is minimized.
+  for j = 0:n
+    vec(j+1) = negativeE(y,j)/partial_sumA(y,j) - log10(partial_sumA(y,j)) + ...
+        (negativeE(y,n)-negativeE(y,j))/(partial_sumA(y,n)-partial_sumA(y,j)) - log10(partial_sumA(y,n)-partial_sumA(y,j));
+  end
+  warning on
+
+  [minimum,ind] = min(vec);
+  T = ind-1;
+endfunction
+
+function [T] = intermodes (I,n)
+  if nargin == 1
+    n = 255;
+  end
+
+  I = double(I);
+
+  % Calculate the histogram.
+  y = hist(I(:),0:n);
+
+  % Smooth the histogram by iterative three point mean filtering.
+  iter = 0;
+  while ~bimodtest(y)
+    h = ones(1,3)/3;
+    y = conv2(y,h,'same');
+    iter = iter+1;
+    % If the histogram turns out not to be bimodal, set T to zero.
+    if iter > 10000;
+      T = 0;
+      return
+    end
+  end
+
+  % The threshold is the mean of the two peaks of the histogram.
+  ind = 0;
+  for k = 2:n
+    if y(k-1) < y(k) && y(k+1) < y(k)
+      ind = ind+1;
+      TT(ind) = k-1;
+    end
+  end
+  T = floor(mean(TT));
+endfunction
+
+function [T] = percentile (I, p, n)
+  if nargin == 1
+    p = 0.5;
+    n = 255;
+  elseif nargin == 2
+    n = 255;
+  end
+
+  I = double(I);
+
+  % Calculate the histogram.
+  y = hist(I(:),0:n);
+
+  % The threshold is chosen such that 50% of pixels lie in each category.
+  Avec = zeros(1,n+1);
+  for t = 0:n
+    Avec(t+1) = partial_sumA(y,t)/partial_sumA(y,n);
+  end
+
+  [minimum,ind] = min(abs(Avec-p));
+  T = ind-1;
+endfunction
+
+
+function T = minimum(I,n);
+  if nargin == 1
+    n = 255;
+  end
+
+  I = double(I);
+
+  % Calculate the histogram.
+  y = hist(I(:),0:n);
+
+  % Smooth the histogram by iterative three point mean filtering.
+  iter = 0;
+  while ~bimodtest(y)
+    h = ones(1,3)/3;
+    y = conv2(y,h,'same');
+    iter = iter+1;
+    % If the histogram turns out not to be bimodal, set T to zero.
+    if iter > 10000;
+      T = 0;
+      return
+    end
+  end
+
+  % The threshold is the minimum between the two peaks.
+  for k = 2:n
+    if y(k-1) > y(k) && y(k+1) > y(k)
+      T = k-1;
+    end
+  end
+endfunction
+
+function [T] = minerror_iter (I,n)
+  if nargin == 1
+    n = 255;
+  end
+
+  I = double(I);
+
+  % Calculate the histogram.
+  y = hist(I(:),0:n);
+
+  % The initial estimate for the threshold is found with the MEAN algorithm.
+  T = floor (mean (I, img(:)));
+  Tprev = NaN;
+
+  while T ~= Tprev
+
+    % Calculate some statistics.
+    mu = partial_sumB(y,T)/partial_sumA(y,T);
+    nu = (partial_sumB(y,n)-partial_sumB(y,T))/(partial_sumA(y,n)-partial_sumA(y,T));
+    p = partial_sumA(y,T)/partial_sumA(y,n);
+    q = (partial_sumA(y,n)-partial_sumA(y,T)) / partial_sumA(y,n);
+    sigma2 = partial_sumC(y,T)/partial_sumA(y,T)-mu^2;
+    tau2 = (partial_sumC(y,n)-partial_sumC(y,T)) / (partial_sumA(y,n)-partial_sumA(y,T)) - nu^2;
+
+    % The terms of the quadratic equation to be solved.
+    w0 = 1/sigma2-1/tau2;
+    w1 = mu/sigma2-nu/tau2;
+    w2 = mu^2/sigma2 - nu^2/tau2 + log10((sigma2*q^2)/(tau2*p^2));
+
+    % If the next threshold would be imaginary, return with the current one.
+    sqterm = w1^2-w0*w2;
+    if sqterm < 0
+      warning('MINERROR:NaN','Warning: th_minerror_iter did not converge.')
+      return
+    end
+
+    % The updated threshold is the integer part of the solution of the
+    % quadratic equation.
+    Tprev = T;
+    T = floor((w1+sqrt(sqterm))/w0);
+
+    % If the threshold turns out to be NaN, return with the previous threshold.
+    if isnan(T)
+      warning('MINERROR:NaN','Warning: th_minerror_iter did not converge.')
+      T = Tprev;
+    end
+
+  end
+endfunction
+#{
+  ## this was not an implementatino of the original minerror algorithm but seems
+  ## to converg more often than the original. The original is (also from the
+  ## HistThresh toolbox
+  function T = th_minerror(I,n)
+    if nargin == 1
+      n = 255;
+    end
+
+    I = double(I);
+
+    % Calculate the histogram.
+    y = hist(I(:),0:n);
+
+    warning off
+    % The threshold is chosen such that the following expression is minimized.
+    for j = 0:n
+      mu = partial_sumB(y,j)/partial_sumA(y,j);
+      nu = (partial_sumB(y,n)-partial_sumB(y,j))/(partial_sumA(y,n)-partial_sumA(y,j));
+      p = partial_sumA(y,j)/partial_sumA(y,n);
+      q = (partial_sumA(y,n)-partial_sumA(y,j)) / partial_sumA(y,n);
+      sigma2 = partial_sumC(y,j)/partial_sumA(y,j)-mu^2;
+      tau2 = (partial_sumC(y,n)-partial_sumC(y,j)) / (partial_sumA(y,n)-partial_sumA(y,j)) - nu^2;
+      vec(j+1) = p*log10(sqrt(sigma2)/p) + q*log10(sqrt(tau2)/q);
+    end
+    warning on
+
+    vec(vec==-inf) = NaN;
+    [minimum,ind] = min(vec);
+    T = ind-1;
+  endfunction
+#}
+
+function T = maxlikelihood (I,n)
+  if nargin == 1
+    n = 255;
+  end
+
+  I = double(I);
+
+  % Calculate the histogram.
+  y = hist(I(:),0:n);
+
+  % The initial estimate for the threshold is found with the MINIMUM
+  % algorithm.
+  T = th_minimum(I,n);
+
+  % Calculate initial values for the statistics.
+  mu = partial_sumB(y,T)/partial_sumA(y,T);
+  nu = (partial_sumB(y,n)-partial_sumB(y,T))/(partial_sumA(y,n)-partial_sumA(y,T));
+  p = partial_sumA(y,T)/partial_sumA(y,n);
+  q = (partial_sumA(y,n)-partial_sumA(y,T)) / partial_sumA(y,n);
+  sigma2 = partial_sumC(y,T)/partial_sumA(y,T)-mu^2;
+  tau2 = (partial_sumC(y,n)-partial_sumC(y,T)) / (partial_sumA(y,n)-partial_sumA(y,T)) - nu^2;
+
+  mu_prev = NaN;
+  nu_prev = NaN;
+  p_prev = NaN;
+  q_prev = NaN;
+  sigma2_prev = NaN;
+  tau2_prev = NaN;
+
+  while abs(mu-mu_prev) > eps || abs(nu-nu_prev) > eps || ...
+        abs(p-p_prev) > eps || abs(q-q_prev) > eps || ...
+        abs(sigma2-sigma2_prev) > eps || abs(tau2-tau2_prev) > eps
+    for i = 0:n
+      phi(i+1) = p/q * exp(-((i-mu)^2) / (2*sigma2)) / ...
+          (p/sqrt(sigma2) * exp(-((i-mu)^2) / (2*sigma2)) + ...
+           (q/sqrt(tau2)) * exp(-((i-nu)^2) / (2*tau2)));
+    end
+    ind = 0:n;
+    gamma = 1-phi;
+    F = phi*y';
+    G = gamma*y';
+    p_prev = p;
+    q_prev = q;
+    mu_prev = mu;
+    nu_prev = nu;
+    sigma2_prev = nu;
+    tau2_prev = nu;
+    p = F/partial_sumA(y,n);
+    q = G/partial_sumA(y,n);
+    mu = ind.*phi*y'/F;
+    nu = ind.*gamma*y'/G;
+    sigma2 = ind.^2.*phi*y'/F - mu^2;
+    tau2 = ind.^2.*gamma*y'/G - nu^2;
+  end
+
+  % The terms of the quadratic equation to be solved.
+  w0 = 1/sigma2-1/tau2;
+  w1 = mu/sigma2-nu/tau2;
+  w2 = mu^2/sigma2 - nu^2/tau2 + log10((sigma2*q^2)/(tau2*p^2));
+
+  % If the threshold would be imaginary, return with threshold set to zero.
+  sqterm = w1^2-w0*w2;
+  if sqterm < 0;
+    T = 0;
+    return
+  end
+
+  % The threshold is the integer part of the solution of the quadratic
+  % equation.
+  T = floor((w1+sqrt(sqterm))/w0);
+endfunction
+
+function T = intermeans (I,n)
+  if nargin == 1
+    n = 255;
+  end
+
+  I = double(I);
+
+  % Calculate the histogram.
+  y = hist(I(:),0:n);
+
+  % The initial estimate for the threshold is found with the MEAN algorithm.
+  T = floor (mean (I, img(:)));
+  Tprev = NaN;
+
+  % The threshold is found iteratively. In each iteration, the means of the
+  % pixels below (mu) the threshold and above (nu) it are found. The
+  % updated threshold is the mean of mu and nu.
+  while T ~= Tprev
+    mu = partial_sumB(y,T)/partial_sumA(y,T);
+    nu = (partial_sumB(y,n)-partial_sumB(y,T))/(partial_sumA(y,n)-partial_sumA(y,T));
+    Tprev = T;
+    T = floor((mu+nu)/2);
+  end
+endfunction
+
+function T = concavity (I,n)
+  if nargin == 1
+    n = 255;
+  end
+
+  I = double(I);
+
+  % Calculate the histogram and its convex hull.
+  h = hist(I(:),0:n);
+  H = hconvhull(h);
+
+  % Find the local maxima of the difference H-h.
+  lmax = flocmax(H-h);
+
+  % Find the histogram balance around each index.
+  for k = 0:n
+    E(k+1) = hbalance(h,k);
+  end
+
+  % The threshold is the local maximum with highest balance.
+  E = E.*lmax;
+  [dummy ind] = max(E);
+  T = ind-1;
+endfunction
+
+################################################################################
+## Auxiliary functions from HistThresh toolbox http://www.cs.tut.fi/~ant/histthresh/
+################################################################################
+
+## partial sums from C. A. Glasbey, "An analysis of histogram-based thresholding
+## algorithms," CVGIP: Graphical Models and Image Processing, vol. 55, pp. 532-537, 1993.
+function x = partial_sumA (y, j)
+  x = sum (y(1:j+1));
+endfunction
+function x = partial_sumB (y, j)
+  ind = 0:j;
+  x   = ind*y(1:j+1)';
+endfunction
+function x = partial_sumC (y, j)
+  ind = 0:j;
+  x = ind.^2*y(1:j+1)';
+endfunction
+function x = partial_sumD (y, j)
+  ind = 0:j;
+  x = ind.^3*y(1:j+1)';
+endfunction
+
+## Test if a histogram is bimodal.
+function b = bimodtest(y)
+  len = length(y);
+  b = false;
+  modes = 0;
+
+  % Count the number of modes of the histogram in a loop. If the number
+  % exceeds 2, return with boolean return value false.
+  for k = 2:len-1
+    if y(k-1) < y(k) && y(k+1) < y(k)
+      modes = modes+1;
+      if modes > 2
+        return
+      end
+    end
+  end
+
+  % The number of modes could be less than two here
+  if modes == 2
+    b = true;
+  end
+endfunction
+
+## Find the local maxima of a vector using a three point neighborhood.
+function y = flocmax(x)
+%  y    binary vector with maxima of x marked as ones
+
+  len = length(x);
+  y = zeros(1,len);
+
+  for k = 2:len-1
+    [dummy,ind] = max(x(k-1:k+1));
+    if ind == 2
+      y(k) = 1;
+    end
+  end
+endfunction
+
+## Calculate the balance measure of the histogram around a histogram index.
+function E = hbalance(y,ind)
+%  y    histogram
+%  ind  index about which balance is calculated
+%
+% Out:
+%  E    balance measure
+%
+% References:
+%
+% A. Rosenfeld and P. De La Torre, "Histogram concavity analysis as an aid
+% in threhold selection," IEEE Transactions on Systems, Man, and
+% Cybernetics, vol. 13, pp. 231-235, 1983.
+%
+% P. K. Sahoo, S. Soltani, and A. K. C. Wong, "A survey of thresholding
+% techniques," Computer Vision, Graphics, and Image Processing, vol. 41,
+% pp. 233-260, 1988.
+
+  n = length(y)-1;
+  E = partial_sumA(y,ind)*(partial_sumA(y,n)-partial_sumA(y,ind));
+endfunction
+
+## Find the convex hull of a histogram.
+function H = hconvhull(h)
+  % In:
+  %  h    histogram
+  %
+  % Out:
+  %  H    convex hull of histogram
+  %
+  % References:
+  %
+  % A. Rosenfeld and P. De La Torre, "Histogram concavity analysis as an aid
+  % in threhold selection," IEEE Transactions on Systems, Man, and
+  % Cybernetics, vol. 13, pp. 231-235, 1983.
+
+  len = length(h);
+  K(1) = 1;
+  k = 1;
+
+  % The vector K gives the locations of the vertices of the convex hull.
+  while K(k)~=len
+
+    theta = zeros(1,len-K(k));
+    for i = K(k)+1:len
+      x = i-K(k);
+      y = h(i)-h(K(k));
+      theta(i-K(k)) = atan2(y,x);
+    end
+
+    maximum = max(theta);
+    maxloc = find(theta==maximum);
+    k = k+1;
+    K(k) = maxloc(end)+K(k-1);
+
+  end
+
+  % Form the convex hull.
+  H = zeros(1,len);
+  for i = 2:length(K)
+    H(K(i-1):K(i)) = h(K(i-1))+(h(K(i))-h(K(i-1)))/(K(i)-K(i-1))*(0:K(i)-K(i-1));
+  end
+endfunction
+
+## Entroy function. Note that the function returns the negative of entropy.
+function x = negativeE(y,j)
+  ## used by the maxentropy method only
+  y = y(1:j+1);
+  y = y(y~=0);
+  x = sum(y.*log10(y));
+endfunction

--- a/main/image/inst/imbothat.m	Thu Sep 27 12:51:44 2012 +0000
+++ b/main/image/inst/imbothat.m	Thu Sep 27 13:18:13 2012 +0000
@@ -1,55 +1,55 @@
-## Copyright (C) 2005  Carvalho-Mariel
-## Copyright (C) 2010-2011 Carnë Draug <carandraug+dev@gmail.com>
-##
-## 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 3 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.
-##
-## You should have received a copy of the GNU General Public License along with
-## this program; if not, see <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn {Function File} @var{B} = imbothat (@var{A}, @var{se})
-## Perform morphological bottom hat filtering.
-##
-## The image @var{A} must be a grayscale or binary image, and @var{se} must be a
-## structuring element. Both must have the same class, e.g., if @var{A} is a
-## logical matrix, @var{se} must also be logical.
-##
-## A bottom-hat transform is also known as 'black', or 'closing', top-hat
-## transform.
-##
-## @seealso{imerode, imdilate, imopen, imclose, imtophat, mmgradm}
-## @end deftypefn
-
-function retval = imbothat (im, se)
-
-  ## Checkinput
-  if (nargin != 2)
-    print_usage();
-  endif
-  if (!ismatrix(im) || !isreal(im))
-    error("imbothat: first input argument must be a real matrix");
-  elseif (!ismatrix(se) || !isreal(se))
-    error("imbothat: second input argument must be a real matrix");
-  elseif ( !strcmp(class(im), class(se)) )
-    error("imbothat: image and structuring element must have the same class");
-  endif
-
-  ## Perform filtering
-  ## Note that in case that the transform is to applied to a logical image,
-  ## subtraction must be handled in a different way (x & !y) instead of (x - y)
-  ## or it will return a double precision matrix
-  if (islogical(im))
-    retval = imclose(im, se) & !im;
-  else
-    retval = imclose(im, se) - im;
-  endif
-
-endfunction
+## Copyright (C) 2005  Carvalho-Mariel
+## Copyright (C) 2010-2011 Carnë Draug <carandraug+dev@gmail.com>
+##
+## 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 3 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.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} @var{B} = imbothat (@var{A}, @var{se})
+## Perform morphological bottom hat filtering.
+##
+## The image @var{A} must be a grayscale or binary image, and @var{se} must be a
+## structuring element. Both must have the same class, e.g., if @var{A} is a
+## logical matrix, @var{se} must also be logical.
+##
+## A bottom-hat transform is also known as 'black', or 'closing', top-hat
+## transform.
+##
+## @seealso{imerode, imdilate, imopen, imclose, imtophat, mmgradm}
+## @end deftypefn
+
+function retval = imbothat (im, se)
+
+  ## Checkinput
+  if (nargin != 2)
+    print_usage();
+  endif
+  if (!ismatrix(im) || !isreal(im))
+    error("imbothat: first input argument must be a real matrix");
+  elseif (!ismatrix(se) || !isreal(se))
+    error("imbothat: second input argument must be a real matrix");
+  elseif ( !strcmp(class(im), class(se)) )
+    error("imbothat: image and structuring element must have the same class");
+  endif
+
+  ## Perform filtering
+  ## Note that in case that the transform is to applied to a logical image,
+  ## subtraction must be handled in a different way (x & !y) instead of (x - y)
+  ## or it will return a double precision matrix
+  if (islogical(im))
+    retval = imclose(im, se) & !im;
+  else
+    retval = imclose(im, se) - im;
+  endif
+
+endfunction

--- a/main/image/inst/imclose.m	Thu Sep 27 12:51:44 2012 +0000
+++ b/main/image/inst/imclose.m	Thu Sep 27 13:18:13 2012 +0000
@@ -1,44 +1,44 @@
-## Copyright (C) 2008 Soren Hauberg <soren@hauberg.org>
-##
-## 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 3 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.
-##
-## You should have received a copy of the GNU General Public License along with
-## this program; if not, see <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn {Function File} @var{B} = imclose (@var{A}, @var{se})
-## Perform morphological closing on a given image.
-## The image @var{A} must be a grayscale or binary image, and @var{se} must be a
-## structuring element.
-##
-## The closing corresponds to a dilation followed by an erosion of the image, i.e.
-## @example
-## B = imerode(imdilate(A, se), se);
-## @end example
-## @seealso{imdilate, imerode, imclose}
-## @end deftypefn
-
-function retval = imclose (im, se)
-  ## Checkinput
-  if (nargin != 2)
-    print_usage();
-  endif
-  if (!ismatrix(im) || !isreal(im))
-    error("imclose: first input argument must be a real matrix");
-  endif
-  if (!ismatrix(se) || !isreal(se))
-    error("imclose: second input argument must be a real matrix");
-  endif
-  
-  ## Perform filtering
-  retval = imerode(imdilate(im, se), se);
-
-endfunction
+## Copyright (C) 2008 Soren Hauberg <soren@hauberg.org>
+##
+## 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 3 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.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} @var{B} = imclose (@var{A}, @var{se})
+## Perform morphological closing on a given image.
+## The image @var{A} must be a grayscale or binary image, and @var{se} must be a
+## structuring element.
+##
+## The closing corresponds to a dilation followed by an erosion of the image, i.e.
+## @example
+## B = imerode(imdilate(A, se), se);
+## @end example
+## @seealso{imdilate, imerode, imclose}
+## @end deftypefn
+
+function retval = imclose (im, se)
+  ## Checkinput
+  if (nargin != 2)
+    print_usage();
+  endif
+  if (!ismatrix(im) || !isreal(im))
+    error("imclose: first input argument must be a real matrix");
+  endif
+  if (!ismatrix(se) || !isreal(se))
+    error("imclose: second input argument must be a real matrix");
+  endif
+  
+  ## Perform filtering
+  retval = imerode(imdilate(im, se), se);
+
+endfunction

--- a/main/image/inst/imdilate.m	Thu Sep 27 12:51:44 2012 +0000
+++ b/main/image/inst/imdilate.m	Thu Sep 27 13:18:13 2012 +0000
@@ -1,69 +1,69 @@
-## Copyright (C) 2004 Josep Mones i Teixidor <jmones@puntbarra.com>
-## Copyright (C) 2008 Soren Hauberg <soren@hauberg.org>
-## Copyright (C) 2010 Carnë Draug <carandraug+dev@gmail.com>
-##
-## 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 3 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.
-##
-## You should have received a copy of the GNU General Public License along with
-## this program; if not, see <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn {Function File} {@var{B} =} imdilate (@var{A}, @var{se})
-## Perform morphological dilation on a given image.
-##
-## The image @var{A} must be a grayscale or binary image, and @var{se} a
-## structuring element. Both must have the same class, e.g., if @var{A} is a
-## logical matrix, @var{se} must also be logical. Note that the erosion
-## algorithm is different for each class, being much faster for logical
-## matrices. As such, if you have a binary matrix, you should use @code{logical}
-## first. This will also reduce the memory usage of your code.
-##
-## The center of @var{SE} is calculated using floor((size(@var{SE})+1)/2).
-##
-## Pixels outside the image are considered to be 0.
-##
-## @seealso{imerode, imopen, imclose}
-## @end deftypefn
-
-function retval = imdilate(im, se)
-  ## Checkinput
-  if (nargin != 2)
-    print_usage();
-  endif
-  if (!ismatrix(im) || !isreal(im))
-    error("imdilate: first input argument must be a real matrix");
-  elseif (!ismatrix(se) || !isreal(se))
-    error("imdilate: second input argument must be a real matrix");
-  elseif ( !strcmp(class(im), class(se)) )
-    error("imdilate: image and structuring element must have the same class");
-  endif
-
-  ## Perform filtering
-  ## Filtering must be done with the reflection of the structuring element (they
-  ## are not always symmetrical)
-  se      = imrotate(se, 180);
-
-  ## If image is binary/logical, try to use filter2 (much faster)
-  if (islogical(im))
-    retval  = filter2(se,im)>0;
-  else
-    retval  = ordfiltn(im, sum(se(:)!=0), se, 0);
-  endif
-
-endfunction
-
-%!demo
-%! imdilate(eye(5),ones(2,2))
-%! % returns a thick diagonal.
-
-%!assert(imdilate(eye(3),[1])==eye(3));                     # using [1] as a mask returns the same value
-%!assert(imdilate(eye(3),[1,0,0])==[0,0,0;1,0,0;0,1,0]);    # check if it works with non-symmetric SE
-%!assert(imdilate(eye(5),[1,0,0,0])==[0,0,0,0,0;1,0,0,0,0;0,1,0,0,0;0,0,1,0,0;0,0,0,1,0]); # test if center is correctly calculated on even masks
+## Copyright (C) 2004 Josep Mones i Teixidor <jmones@puntbarra.com>
+## Copyright (C) 2008 Soren Hauberg <soren@hauberg.org>
+## Copyright (C) 2010 Carnë Draug <carandraug+dev@gmail.com>
+##
+## 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 3 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.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{B} =} imdilate (@var{A}, @var{se})
+## Perform morphological dilation on a given image.
+##
+## The image @var{A} must be a grayscale or binary image, and @var{se} a
+## structuring element. Both must have the same class, e.g., if @var{A} is a
+## logical matrix, @var{se} must also be logical. Note that the erosion
+## algorithm is different for each class, being much faster for logical
+## matrices. As such, if you have a binary matrix, you should use @code{logical}
+## first. This will also reduce the memory usage of your code.
+##
+## The center of @var{SE} is calculated using floor((size(@var{SE})+1)/2).
+##
+## Pixels outside the image are considered to be 0.
+##
+## @seealso{imerode, imopen, imclose}
+## @end deftypefn
+
+function retval = imdilate(im, se)
+  ## Checkinput
+  if (nargin != 2)
+    print_usage();
+  endif
+  if (!ismatrix(im) || !isreal(im))
+    error("imdilate: first input argument must be a real matrix");
+  elseif (!ismatrix(se) || !isreal(se))
+    error("imdilate: second input argument must be a real matrix");
+  elseif ( !strcmp(class(im), class(se)) )
+    error("imdilate: image and structuring element must have the same class");
+  endif
+
+  ## Perform filtering
+  ## Filtering must be done with the reflection of the structuring element (they
+  ## are not always symmetrical)
+  se      = imrotate(se, 180);
+
+  ## If image is binary/logical, try to use filter2 (much faster)
+  if (islogical(im))
+    retval  = filter2(se,im)>0;
+  else
+    retval  = ordfiltn(im, sum(se(:)!=0), se, 0);
+  endif
+
+endfunction
+
+%!demo
+%! imdilate(eye(5),ones(2,2))
+%! % returns a thick diagonal.
+
+%!assert(imdilate(eye(3),[1])==eye(3));                     # using [1] as a mask returns the same value
+%!assert(imdilate(eye(3),[1,0,0])==[0,0,0;1,0,0;0,1,0]);    # check if it works with non-symmetric SE
+%!assert(imdilate(eye(5),[1,0,0,0])==[0,0,0,0,0;1,0,0,0,0;0,1,0,0,0;0,0,1,0,0;0,0,0,1,0]); # test if center is correctly calculated on even masks

--- a/main/image/inst/imerode.m	Thu Sep 27 12:51:44 2012 +0000
+++ b/main/image/inst/imerode.m	Thu Sep 27 13:18:13 2012 +0000
@@ -1,65 +1,65 @@
-## Copyright (C) 2004 Josep Mones i Teixidor <jmones@puntbarra.com>
-## Copyright (C) 2008 Soren Hauberg <soren@hauberg.org>
-## Copyright (C) 2011 Carnë Draug <carandraug+dev@gmail.com>
-##
-## 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 3 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.
-##
-## You should have received a copy of the GNU General Public License along with
-## this program; if not, see <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn {Function File} {@var{B} =} imerode (@var{A}, @var{se})
-## Perform morphological erosion on a given image.
-##
-## The image @var{A} must be a grayscale or binary image, and @var{se} a
-## structuring element. Both must have the same class, e.g., if @var{A} is a
-## logical matrix, @var{se} must also be logical. Note that the erosion
-## algorithm is different for each class, being much faster for logical
-## matrices. As such, if you have a binary matrix, you should use @code{logical}
-## first. This will also reduce the memory usage of your code.
-##
-## The center of @var{SE} is calculated using floor((size(@var{SE})+1)/2).
-##
-## Pixels outside the image are considered to be 0.
-##
-## @seealso{imdilate, imopen, imclose}
-## @end deftypefn
-
-function retval = imerode(im, se)
-  ## Checkinput
-  if (nargin != 2)
-    print_usage();
-  endif
-  if (!ismatrix(im) || !isreal(im))
-    error("imerode: first input argument must be a real matrix");
-  elseif (!ismatrix(se) || !isreal(se))
-    error("imerode: second input argument must be a real matrix");
-  elseif ( !strcmp(class(im), class(se)) )
-    error("imerode: image and structuring element must have the same class");
-  endif
-
-  ## Perform filtering
-  ## If image is binary/logical, try to use filter2 (much faster)
-  if (islogical(im))
-    thr     = sum(se(:));
-    retval  = filter2(se,im) == thr;
-  else
-    retval  = ordfiltn(im, 1, se, 0);
-  endif
-
-endfunction
-
-%!demo
-%! imerode(ones(5,5),ones(3,3))
-%! % creates a zeros border around ones.
-
-%!assert(imerode(eye(3),[1])==eye(3));                                           # using [1] as a mask returns the same value
-%!assert(imerode([0,1,0;1,1,1;0,1,0],[0,0,0;0,0,1;0,1,1])==[1,0,0;0,0,0;0,0,0]); # check if it works with non-symmetric SE
+## Copyright (C) 2004 Josep Mones i Teixidor <jmones@puntbarra.com>
+## Copyright (C) 2008 Soren Hauberg <soren@hauberg.org>
+## Copyright (C) 2011 Carnë Draug <carandraug+dev@gmail.com>
+##
+## 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 3 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.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{B} =} imerode (@var{A}, @var{se})
+## Perform morphological erosion on a given image.
+##
+## The image @var{A} must be a grayscale or binary image, and @var{se} a
+## structuring element. Both must have the same class, e.g., if @var{A} is a
+## logical matrix, @var{se} must also be logical. Note that the erosion
+## algorithm is different for each class, being much faster for logical
+## matrices. As such, if you have a binary matrix, you should use @code{logical}
+## first. This will also reduce the memory usage of your code.
+##
+## The center of @var{SE} is calculated using floor((size(@var{SE})+1)/2).
+##
+## Pixels outside the image are considered to be 0.
+##
+## @seealso{imdilate, imopen, imclose}
+## @end deftypefn
+
+function retval = imerode(im, se)
+  ## Checkinput
+  if (nargin != 2)
+    print_usage();
+  endif
+  if (!ismatrix(im) || !isreal(im))
+    error("imerode: first input argument must be a real matrix");
+  elseif (!ismatrix(se) || !isreal(se))
+    error("imerode: second input argument must be a real matrix");
+  elseif ( !strcmp(class(im), class(se)) )
+    error("imerode: image and structuring element must have the same class");
+  endif
+
+  ## Perform filtering
+  ## If image is binary/logical, try to use filter2 (much faster)
+  if (islogical(im))
+    thr     = sum(se(:));
+    retval  = filter2(se,im) == thr;
+  else
+    retval  = ordfiltn(im, 1, se, 0);
+  endif
+
+endfunction
+
+%!demo
+%! imerode(ones(5,5),ones(3,3))
+%! % creates a zeros border around ones.
+
+%!assert(imerode(eye(3),[1])==eye(3));                                           # using [1] as a mask returns the same value
+%!assert(imerode([0,1,0;1,1,1;0,1,0],[0,0,0;0,0,1;0,1,1])==[1,0,0;0,0,0;0,0,0]); # check if it works with non-symmetric SE

--- a/main/image/inst/imopen.m	Thu Sep 27 12:51:44 2012 +0000
+++ b/main/image/inst/imopen.m	Thu Sep 27 13:18:13 2012 +0000
@@ -1,44 +1,44 @@
-## Copyright (C) 2008 Soren Hauberg <soren@hauberg.org>
-##
-## 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 3 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.
-##
-## You should have received a copy of the GNU General Public License along with
-## this program; if not, see <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn {Function File} @var{B} = imopen (@var{A}, @var{se})
-## Perform morphological opening on a given image.
-## The image @var{A} must be a grayscale or binary image, and @var{se} must be a
-## structuring element.
-##
-## The opening corresponds to an erosion followed by a dilation of the image, i.e.
-## @example
-## B = imdilate(imerode(A, se), se);
-## @end example
-## @seealso{imdilate, imerode, imclose}
-## @end deftypefn
-
-function retval = imopen(im, se)
-  ## Checkinput
-  if (nargin != 2)
-    print_usage();
-  endif
-  if (!ismatrix(im) || !isreal(im))
-    error("imopen: first input argument must be a real matrix");
-  endif
-  if (!ismatrix(se) || !isreal(se))
-    error("imopen: second input argument must be a real matrix");
-  endif
-  
-  ## Perform filtering
-  retval = imdilate(imerode(im, se), se);
-
-endfunction
+## Copyright (C) 2008 Soren Hauberg <soren@hauberg.org>
+##
+## 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 3 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.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} @var{B} = imopen (@var{A}, @var{se})
+## Perform morphological opening on a given image.
+## The image @var{A} must be a grayscale or binary image, and @var{se} must be a
+## structuring element.
+##
+## The opening corresponds to an erosion followed by a dilation of the image, i.e.
+## @example
+## B = imdilate(imerode(A, se), se);
+## @end example
+## @seealso{imdilate, imerode, imclose}
+## @end deftypefn
+
+function retval = imopen(im, se)
+  ## Checkinput
+  if (nargin != 2)
+    print_usage();
+  endif
+  if (!ismatrix(im) || !isreal(im))
+    error("imopen: first input argument must be a real matrix");
+  endif
+  if (!ismatrix(se) || !isreal(se))
+    error("imopen: second input argument must be a real matrix");
+  endif
+  
+  ## Perform filtering
+  retval = imdilate(imerode(im, se), se);
+
+endfunction

--- a/main/image/inst/imtophat.m	Thu Sep 27 12:51:44 2012 +0000
+++ b/main/image/inst/imtophat.m	Thu Sep 27 13:18:13 2012 +0000
@@ -1,78 +1,78 @@
-## Copyright (C) 2005 Carvalho-Mariel
-## Copyright (C) 2010, 2011 Carnë Draug <carandraug+dev@gmail.com>
-##
-## 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 3 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.
-##
-## You should have received a copy of the GNU General Public License along with
-## this program; if not, see <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn {Function File} @var{B} = imtophat (@var{A}, @var{se})
-## Perform morphological top hat filtering.
-##
-## The image @var{A} must be a grayscale or binary image, and @var{se} must be a
-## structuring element. Both must have the same class, e.g., if @var{A} is a
-## logical matrix, @var{se} must also be logical.
-##
-## A 'black', or 'closing', top-hat transform is also known as bottom-hat
-## transform and so that is the same @code{imbothat}.
-##
-## @seealso{imerode, imdilate, imopen, imclose, imbothat, mmgradm}
-## @end deftypefn
-
-function retval = imtophat (im, se, trans)
-
-  ## Checkinput
-  if (nargin < 2 || nargin > 3)
-    print_usage();
-  elseif (nargin == 2)
-    trans = "white";
-  elseif (nargin == 3 && ischar (trans))
-    warning ("Use of the 'trans' argument in imtophat to specify a closing or opening top-hat transform has been deprecated in favor of using the functions 'imtophat' or 'imbothat'. This option will be removed from future versions of the 'imtophat' function");
-  endif
-  if (!ismatrix(im) || !isreal(im))
-    error("imtophat: first input argument must be a real matrix");
-  elseif (!ismatrix(se) || !isreal(se))
-    error("imtophat: second input argument must be a real matrix");
-  elseif ( !strcmp(class(im), class(se)) )
-    error("imtophat: image and structuring element must have the same class");
-  elseif (!ischar (trans))
-    error("imtophat: third input argument must be a string with 'white', 'open, 'black', or 'close'");
-  endif
-
-  ## Perform filtering
-  ## Note that in case that the transform is to applied to a logical image,
-  ## subtraction must be handled in a different way (x & !y) instead of (x - y)
-  ## or it will return a double precision matrix
-  if ( strcmpi(trans, "white") || strcmpi(trans, "open") )
-    if (islogical(im))
-      retval = im & !imopen(im,se);
-    else
-      retval = im - imopen(im, se);
-    endif
-  elseif ( strcmpi(trans, "black") || strcmpi(trans, "close") )
-    retval = imbothat (im, se);
-  else
-    error ("Unexpected type of top-hat transform");
-  endif
-
-endfunction
-
-%!test
-%! I = [1 1 1; 1 1 1; 1 1 1;];
-%! se = [1 1; 0 1;];
-%! ## class of input should be the same as the output
-%! result = imtophat(logical(I), logical(se));
-%! expected = 0.5 < [0 0 1; 0 0 1; 1 1 1];
-%! assert(expected, result);
-%! result = imtophat((I), (se));
-%! expected = [0 0 1; 0 0 1; 1 1 1];
-%! assert(expected, result);
+## Copyright (C) 2005 Carvalho-Mariel
+## Copyright (C) 2010, 2011 Carnë Draug <carandraug+dev@gmail.com>
+##
+## 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 3 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.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} @var{B} = imtophat (@var{A}, @var{se})
+## Perform morphological top hat filtering.
+##
+## The image @var{A} must be a grayscale or binary image, and @var{se} must be a
+## structuring element. Both must have the same class, e.g., if @var{A} is a
+## logical matrix, @var{se} must also be logical.
+##
+## A 'black', or 'closing', top-hat transform is also known as bottom-hat
+## transform and so that is the same @code{imbothat}.
+##
+## @seealso{imerode, imdilate, imopen, imclose, imbothat, mmgradm}
+## @end deftypefn
+
+function retval = imtophat (im, se, trans)
+
+  ## Checkinput
+  if (nargin < 2 || nargin > 3)
+    print_usage();
+  elseif (nargin == 2)
+    trans = "white";
+  elseif (nargin == 3 && ischar (trans))
+    warning ("Use of the 'trans' argument in imtophat to specify a closing or opening top-hat transform has been deprecated in favor of using the functions 'imtophat' or 'imbothat'. This option will be removed from future versions of the 'imtophat' function");
+  endif
+  if (!ismatrix(im) || !isreal(im))
+    error("imtophat: first input argument must be a real matrix");
+  elseif (!ismatrix(se) || !isreal(se))
+    error("imtophat: second input argument must be a real matrix");
+  elseif ( !strcmp(class(im), class(se)) )
+    error("imtophat: image and structuring element must have the same class");
+  elseif (!ischar (trans))
+    error("imtophat: third input argument must be a string with 'white', 'open, 'black', or 'close'");
+  endif
+
+  ## Perform filtering
+  ## Note that in case that the transform is to applied to a logical image,
+  ## subtraction must be handled in a different way (x & !y) instead of (x - y)
+  ## or it will return a double precision matrix
+  if ( strcmpi(trans, "white") || strcmpi(trans, "open") )
+    if (islogical(im))
+      retval = im & !imopen(im,se);
+    else
+      retval = im - imopen(im, se);
+    endif
+  elseif ( strcmpi(trans, "black") || strcmpi(trans, "close") )
+    retval = imbothat (im, se);
+  else
+    error ("Unexpected type of top-hat transform");
+  endif
+
+endfunction
+
+%!test
+%! I = [1 1 1; 1 1 1; 1 1 1;];
+%! se = [1 1; 0 1;];
+%! ## class of input should be the same as the output
+%! result = imtophat(logical(I), logical(se));
+%! expected = 0.5 < [0 0 1; 0 0 1; 1 1 1];
+%! assert(expected, result);
+%! result = imtophat((I), (se));
+%! expected = [0 0 1; 0 0 1; 1 1 1];
+%! assert(expected, result);

--- a/main/image/inst/rgbplot.m	Thu Sep 27 12:51:44 2012 +0000
+++ b/main/image/inst/rgbplot.m	Thu Sep 27 13:18:13 2012 +0000
@@ -1,48 +1,48 @@
-## Copyright (C) 2005 Berge-Gladel
-##
-## 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 3 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.
-##
-## You should have received a copy of the GNU General Public License along with
-## this program; if not, see <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn {Function File} {@var{h} =} rgbplot (@var{cmap})
-## Plot a given color map.
-##
-## The plot will display 3 lines, red, green and blue, showing the variation of
-## RGB values on @var{cmap}.  While input matrix must be a color map (see
-## @code{iscolormap}), it can also be used to see a colored line profile
-## (intensity values will have to be converted to class double with
-## @code{im2double}).
-##
-## If an output is requested, the graphics handle @var{h} to the plot is returned.
-##
-## @seealso{colormap, iscolormap}
-## @end deftypefn
-
-function h_out = rgbplot(map)
-
-  if (nargin != 1)
-    print_usage;
-  elseif (!iscolormap (map))
-    error("rgbplot: input must be a colormap");
-  endif
-
-  h = plot (map(:,1), "-r", map(:,2), "g-", map(:,3), "b-");
-  if (nargout > 0)
-    h_out = h;
-  endif
-endfunction
-
-%!demo
-%! ## look at the distribution of RGB values for the jet colormap
-%! cmap = jet (64);
-%! rgbplot (cmap);
+## Copyright (C) 2005 Berge-Gladel
+##
+## 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 3 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.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{h} =} rgbplot (@var{cmap})
+## Plot a given color map.
+##
+## The plot will display 3 lines, red, green and blue, showing the variation of
+## RGB values on @var{cmap}.  While input matrix must be a color map (see
+## @code{iscolormap}), it can also be used to see a colored line profile
+## (intensity values will have to be converted to class double with
+## @code{im2double}).
+##
+## If an output is requested, the graphics handle @var{h} to the plot is returned.
+##
+## @seealso{colormap, iscolormap}
+## @end deftypefn
+
+function h_out = rgbplot(map)
+
+  if (nargin != 1)
+    print_usage;
+  elseif (!iscolormap (map))
+    error("rgbplot: input must be a colormap");
+  endif
+
+  h = plot (map(:,1), "-r", map(:,2), "g-", map(:,3), "b-");
+  if (nargout > 0)
+    h_out = h;
+  endif
+endfunction
+
+%!demo
+%! ## look at the distribution of RGB values for the jet colormap
+%! cmap = jet (64);
+%! rgbplot (cmap);

--- a/main/image/inst/wavelength2rgb.m	Thu Sep 27 12:51:44 2012 +0000
+++ b/main/image/inst/wavelength2rgb.m	Thu Sep 27 13:18:13 2012 +0000
@@ -1,118 +1,118 @@
-## Copyright (C) 2011 William Krekeler <WKrekeler@cleanearthtech.com>
-## Copyright (C) 2012 Carnë Draug <carandraug+dev@gmail.com>
-##
-## 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 3 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.
-##
-## You should have received a copy of the GNU General Public License along with
-## this program; if not, see <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn{Function File} {@var{rgb} =} wavelength2rgb (@var{wavelength})
-## @deftypefnx{Function File} {@var{rgb} =} wavelength2rgb (@var{wavelength}, @var{intensity_max})
-## @deftypefnx{Function File} {@var{rgb} =} wavelength2rgb (@var{wavelength}, @var{intensity_max}, @var{gamma})
-## Convert wavelength in nm into an RGB value set.
-##
-## Output:
-## @itemize @bullet
-## @item @var{rgb}: value set on scale of 1:@var{intensity_max}
-## @end itemize
-##
-## Input:
-## Output:
-## @itemize @bullet
-## @item @var{wavelength}: wavelength of input light value in nm. Must be a positive
-## numeric scalar.
-## @item @var{intensity_max}: max of integer scale to output values. Defaults to 255.
-## @item @var{gamma}: Controls luminance. Must be a value between 0 and 1. Defaults to 0.8.
-## @end itemize
-##
-## @example
-## X = 350:0.5:800;
-## Y = exp (log (X.^3));   # arbitrary
-## figure, plot(X, Y)
-## stepSize = 0.5 * max (diff (X));
-## for n = 1:numel(Y)
-##   RGB = wavelength2rgb (X(n), 1);
-##   if (sum (RGB) > 0 )   # ie don't fill black
-##     hold on, area( (X(n)-stepSize):(stepSize/1):(X(n)+stepSize), ...
-##                   repmat( Y(n), 1, numel((X(n)-stepSize):(stepSize/1):(X(n)+stepSize)) ), ...
-##                   'EdgeColor', RGB, 'FaceColor', RGB ); hold off
-##   endif
-## endfor
-## @end example
-##
-## Reference:
-## @itemize @bullet
-## @item @uref{http://stackoverflow.com/questions/2374959/algorithm-to-convert-any-positive-integer-to-an-rgb-value}
-## @item @uref{http://www.midnightkite.com/color.html} per Dan Bruton
-## @end itemize
-## @end deftypefn
-
-function rgb = wavelength2rgb (wavelength, intensity_max = 255, gamma = 0.8)
-
-  if (nargin < 1 || nargin > 3)
-    print_usage;
-  elseif (!isnumeric (wavelength) || !isscalar (wavelength) || wavelength <= 0)
-    error ("wavelength must a positive numeric scalar");
-  elseif (!isnumeric (intensity_max) || !isscalar (intensity_max) || intensity_max <= 0)
-    error ("intensity_max must a positive numeric scalar");
-  elseif (!isnumeric (gamma) || !isscalar (gamma) || gamma > 1 || gamma < 0)
-    error ("gamma must a numeric scalar between 1 and 0");
-  endif
-
-  rgb = zeros (3, numel (wavelength));  # initialize rgb, each rgb set stored by column
-
-  ## set the factor
-  if ( wavelength >= 380 && wavelength < 420 )
-     factor = 0.3 + 0.7*(wavelength - 380) / (420 - 380);
-  elseif ( wavelength >= 420 && wavelength < 701 )
-     factor = 1;
-  elseif ( wavelength >= 420 && wavelength < 701 )
-     factor = 0.3 + 0.7*(780 - wavelength) / (780 - 700);
-  else
-     factor = 0;
-  endif
-
-  ## initialize rgb
-  if ( wavelength >= 380 && wavelength < 440 )
-     rgb(1) = -(wavelength - 440) / (440 - 380);
-     rgb(2) = 0.0;
-     rgb(3) = 1.0;
-  elseif ( wavelength >= 440 && wavelength < 490 )
-     rgb(1) = 0.0;
-     rgb(2) = (wavelength - 440) / (490 - 440);
-     rgb(3) = 1.0;
-  elseif ( wavelength >= 490 && wavelength < 510 )
-     rgb(1) = 0.0;
-     rgb(2) = 1.0;
-     rgb(3) = -(wavelength - 510) / (510 - 490);
-  elseif ( wavelength >= 510 && wavelength < 580 )
-     rgb(1) = (wavelength - 510) / (580 - 510);
-     rgb(2) = 1.0;
-     rgb(3) = 0.0;
-  elseif ( wavelength >= 580 && wavelength < 645 )
-     rgb(1) = 1.0;
-     rgb(2) = -(wavelength - 645) / (645 - 580);
-     rgb(3) = 0.0;
-  elseif ( wavelength >= 645 && wavelength < 781 )
-     rgb(1) = 1.0;
-     rgb(2) = 0.0;
-     rgb(3) = 0.0;
-  else
-     rgb(1) = 0.0;
-     rgb(2) = 0.0;
-     rgb(3) = 0.0;
-  endif
-
-  ## correct rgb
-  rgb = intensity_max .* (rgb .* factor) .^gamma .* ( rgb > 0 );
-
-endfunction
+## Copyright (C) 2011 William Krekeler <WKrekeler@cleanearthtech.com>
+## Copyright (C) 2012 Carnë Draug <carandraug+dev@gmail.com>
+##
+## 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 3 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.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn{Function File} {@var{rgb} =} wavelength2rgb (@var{wavelength})
+## @deftypefnx{Function File} {@var{rgb} =} wavelength2rgb (@var{wavelength}, @var{intensity_max})
+## @deftypefnx{Function File} {@var{rgb} =} wavelength2rgb (@var{wavelength}, @var{intensity_max}, @var{gamma})
+## Convert wavelength in nm into an RGB value set.
+##
+## Output:
+## @itemize @bullet
+## @item @var{rgb}: value set on scale of 1:@var{intensity_max}
+## @end itemize
+##
+## Input:
+## Output:
+## @itemize @bullet
+## @item @var{wavelength}: wavelength of input light value in nm. Must be a positive
+## numeric scalar.
+## @item @var{intensity_max}: max of integer scale to output values. Defaults to 255.
+## @item @var{gamma}: Controls luminance. Must be a value between 0 and 1. Defaults to 0.8.
+## @end itemize
+##
+## @example
+## X = 350:0.5:800;
+## Y = exp (log (X.^3));   # arbitrary
+## figure, plot(X, Y)
+## stepSize = 0.5 * max (diff (X));
+## for n = 1:numel(Y)
+##   RGB = wavelength2rgb (X(n), 1);
+##   if (sum (RGB) > 0 )   # ie don't fill black
+##     hold on, area( (X(n)-stepSize):(stepSize/1):(X(n)+stepSize), ...
+##                   repmat( Y(n), 1, numel((X(n)-stepSize):(stepSize/1):(X(n)+stepSize)) ), ...
+##                   'EdgeColor', RGB, 'FaceColor', RGB ); hold off
+##   endif
+## endfor
+## @end example
+##
+## Reference:
+## @itemize @bullet
+## @item @uref{http://stackoverflow.com/questions/2374959/algorithm-to-convert-any-positive-integer-to-an-rgb-value}
+## @item @uref{http://www.midnightkite.com/color.html} per Dan Bruton
+## @end itemize
+## @end deftypefn
+
+function rgb = wavelength2rgb (wavelength, intensity_max = 255, gamma = 0.8)
+
+  if (nargin < 1 || nargin > 3)
+    print_usage;
+  elseif (!isnumeric (wavelength) || !isscalar (wavelength) || wavelength <= 0)
+    error ("wavelength must a positive numeric scalar");
+  elseif (!isnumeric (intensity_max) || !isscalar (intensity_max) || intensity_max <= 0)
+    error ("intensity_max must a positive numeric scalar");
+  elseif (!isnumeric (gamma) || !isscalar (gamma) || gamma > 1 || gamma < 0)
+    error ("gamma must a numeric scalar between 1 and 0");
+  endif
+
+  rgb = zeros (3, numel (wavelength));  # initialize rgb, each rgb set stored by column
+
+  ## set the factor
+  if ( wavelength >= 380 && wavelength < 420 )
+     factor = 0.3 + 0.7*(wavelength - 380) / (420 - 380);
+  elseif ( wavelength >= 420 && wavelength < 701 )
+     factor = 1;
+  elseif ( wavelength >= 420 && wavelength < 701 )
+     factor = 0.3 + 0.7*(780 - wavelength) / (780 - 700);
+  else
+     factor = 0;
+  endif
+
+  ## initialize rgb
+  if ( wavelength >= 380 && wavelength < 440 )
+     rgb(1) = -(wavelength - 440) / (440 - 380);
+     rgb(2) = 0.0;
+     rgb(3) = 1.0;
+  elseif ( wavelength >= 440 && wavelength < 490 )
+     rgb(1) = 0.0;
+     rgb(2) = (wavelength - 440) / (490 - 440);
+     rgb(3) = 1.0;
+  elseif ( wavelength >= 490 && wavelength < 510 )
+     rgb(1) = 0.0;
+     rgb(2) = 1.0;
+     rgb(3) = -(wavelength - 510) / (510 - 490);
+  elseif ( wavelength >= 510 && wavelength < 580 )
+     rgb(1) = (wavelength - 510) / (580 - 510);
+     rgb(2) = 1.0;
+     rgb(3) = 0.0;
+  elseif ( wavelength >= 580 && wavelength < 645 )
+     rgb(1) = 1.0;
+     rgb(2) = -(wavelength - 645) / (645 - 580);
+     rgb(3) = 0.0;
+  elseif ( wavelength >= 645 && wavelength < 781 )
+     rgb(1) = 1.0;
+     rgb(2) = 0.0;
+     rgb(3) = 0.0;
+  else
+     rgb(1) = 0.0;
+     rgb(2) = 0.0;
+     rgb(3) = 0.0;
+  endif
+
+  ## correct rgb
+  rgb = intensity_max .* (rgb .* factor) .^gamma .* ( rgb > 0 );
+
+endfunction