## Copyright (C) 2018 Avinoam Kalma ## Copyright (C) 2018 Carnë Draug ## ## 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 . ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{level}, @var{sep}] =} otsuthresh (@var{hist}) ## Compute global image threshold for histogram, using Otsu's method. ## ## Given an image histogram @var{hist} finds the optimal threshold ## value @var{level} for conversion to a binary image with ## @code{im2bw}. ## ## This function uses ``graythresh'' function which 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}. ## The 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}. It is a value within the range [0 1], the ## lower bound (zero) being attainable by, and only by, histograms having a ## single constant gray level, and the upper bound being attainable by, and only ## by, two-valued pictures. ## ## @seealso{graythresh, im2bw} ## @end deftypefn function [varargout] = otsuthresh (hist) if (nargin != 1) print_usage(); endif if (! isvector (hist) || ! isnumeric(hist) || ! isreal (hist) || issparse (hist) || any (hist == Inf) || any (hist == NaN) || any (hist < 0) || any (hist != fix (hist))) error ("otsuthresh: HIST must be a vector of non-negative integers"); endif if (! strcmp (class (hist), "double")) hist = double (hist); end [varargout{1:nargout}] = graythresh (hist(:).', "otsu"); endfunction %!test %!shared h %! histo = zeros (1, 256); %! histo([ 29 33 37 41 46 50 54 58 62 66 70 74 78 82 ... %! 86 90 94 98 102 106 110 114 118 122 126 131 135 139 ... %! 143 147 151 155 159 163 167 171 175 179 183 187 191 195 ... %! 199 203 207 211 216 220 224 228 232 236 240 244 248 252]) = ... %! [2 27 51 144 132 108 43 29 22 21 22 20 10 16 17 12 13 14 12 13 ... %! 15 25 19 20 23 37 23 65 92 84 87 54 50 54 33 73 76 64 57 58 47 ... %! 48 30 27 22 20 20 11 12 12 11 7 17 31 37 31]; %! assert (otsuthresh (histo), 114.5/255); %! %! I = max (phantom (), 0); %! h = imhist (I); %! assert (otsuthresh (h), 178/255); %! assert (otsuthresh (h'), 178/255); %! h = imhist (I, 10); %! assert (otsuthresh (h), 170/255); %! h(1) = NaN; %!error otsuthresh (h); %! h(1) = Inf; %!error otsuthresh (h); %!test %! H(1) = 100; %! assert (otsuthresh (H), 0); %!test %! H = zeros (1, 256); %! assert (otsuthresh (H), 0); %! assert (otsuthresh (H'), 0); %!test %! H = zeros (1, 5); %! assert (otsuthresh (H), 0); %! assert (otsuthresh (H'), 0); %!test %! H = uint8([10 20 30]); %! assert (otsuthresh (H), 0.5); %! assert (otsuthresh (H'), 0.5); %!test %! H = int32([100 200 300]); %! assert (otsuthresh (H), 0.5); %! assert (otsuthresh (H'), 0.5); %!test %! H = int32([100 200]); %! assert (otsuthresh (H), 0); %! assert (otsuthresh (H'), 0); %!test %! H = single([10 20 30 40]); %! assert (otsuthresh (H), 1/3); %! assert (otsuthresh (H'), 1/3); %!test %! H = uint16([10 20 30 40 50 60 70 80 90 100]); %! assert (otsuthresh (H), 5/9); %! assert (otsuthresh (H'), 5/9); %!test %! H = int16([10 20 30 40 50 60 70 80 90 100]); %! assert (otsuthresh (H), 5/9); %! assert (otsuthresh (H'), 5/9); %!test %! H = int16(1:255); %! assert (otsuthresh (H), 0.614173228346457, 1e-15); %! assert (otsuthresh (H'), 0.614173228346457, 1e-15); %!test %! H = int16(1:1023); %! assert (otsuthresh (H), 0.617416829745597, 1e-15); %! assert (otsuthresh (H'), 0.617416829745597, 1e-15); %!test %! H = int8(1:1023); %! assert (otsuthresh (H), 0.529354207436399, 1e-15); %! assert (otsuthresh (H'), 0.529354207436399, 1e-15); %!test %! q = warning ("query", "Octave:data-file-in-path"); %! warning ("off", "Octave:data-file-in-path"); %! S = load ("penny.mat"); %! warning (q.state, "Octave:data-file-in-path"); %! h = imhist (uint8 (S.P)); %! assert (otsuthresh (h), 94/255); %!test %! I = max (phantom (), 0); %! h = imhist (I, 5); %! assert (otsuthresh (h), 0.625); %!error id=Octave:invalid-fun-call otsuthresh () %!error id=Octave:invalid-fun-call otsuthresh (ones(10), 5) %!error id=Octave:undefined-function otsuthresh (xxx) %!error id=Octave:undefined-function otsuthresh (xxx, 5) %!error id=Octave:undefined-function otsuthresh ([xxx], 5) %!error otsuthresh ([]) %!error otsuthresh ([]) %!error otsuthresh ([Inf]) %!error otsuthresh ([Inf 10]) %!error otsuthresh ([Inf 10]) %!error otsuthresh ([NaN]) %!error otsuthresh ([10 NaN]) %!error otsuthresh (zeros (5)) %!error otsuthresh ([10 -10]) %!error otsuthresh ([10 20; 30 40]) %!error otsuthresh (sparse ([10 10 0 0 0])) %!demo %! I = max (phantom (), 0); %! figure; imshow (I); %! title ("Original image"); %! h = imhist (I); %! t = otsuthresh (h); %! J = im2bw (I); %! figure; imshow (J); %! title_line = sprintf ("Black and white image after thresholding, t=%g", %! t*255); %! title (title_line); %!demo %! q = warning ("query", "Octave:data-file-in-path"); %! warning ("off", "Octave:data-file-in-path"); %! S = load ("penny.mat"); %! warning (q.state, "Octave:data-file-in-path"); %! I = uint8 (S.P); %! figure; imshow (I); %! title ("Original penny image"); %! h = imhist (I); %! t = otsuthresh (h); %! J = im2bw (I); %! figure; imshow (J); %! title_line = sprintf ("Black and white penny image after thresholding, t=%g", %! t*255); %! title (title_line); %! I = 255 - I; %! figure; imshow(I); %! title ("Negative penny image"); %! h = imhist (I); %! t = otsuthresh (h); %! J = im2bw (I); %! figure; imshow (J); %! title_line = sprintf ("Black and white negative penny image after thresholding, t=%g", %! t*255); %! title (title_line);