/usr/share/octave/packages/image-2.2.2/graythresh.m is in octave-image 2.2.2-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 | ## Copyright (C) 2004 Antti Niemistö <antti.niemisto@tut.fi>
## Copyright (C) 2007 Søren Hauberg
## Copyright (C) 2012-2013 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{method}, @var{options})
## @deftypefnx {Function File} {[@var{level}, @var{sep}] =} graythresh (@var{hist}, @dots{})
## 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. An image histogram @var{hist}
## can also be used to allow for preprocessing of the histogram.
##
## The optional argument @var{method} is the algorithm to be used (default's to
## Otsu). Some methods may have other @var{options} and/or return an extra
## value @var{sep} (see each entry for details). The available @var{method}s 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}. 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.
##
## @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 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 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 fraction of pixels (set at @var{options}) to be background.
## If no value is given, assumes 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}
## @end table
##
## @seealso{im2bw}
## @end deftypefn
## Notes:
## * The following methods were adapted from http://www.cs.tut.fi/~ant/histthresh/
## intermodes percentile minimum
## MaxEntropy MaxLikelihood intermeans
## moments minerror concavity
## * Carnë Draug implemented and vectorized the Otsu's method
## * Carnë Draug vectorized percentile and moments.
## * missing methods from ImageJ
## Yen triangle RenyiEntropy
## Shanbhag Li Huang
## ImageJ
function [varargout] = graythresh (img, algo = "otsu", varargin)
## Input checking
if (nargin < 1 || nargin > 3)
print_usage();
elseif (nargin > 2 && !any (strcmpi (algo, {"percentile"})))
error ("graythresh: algorithm `%s' does not accept any options.", algo);
else
hist_in = false;
## If the image is RGB convert it to grayscale
if (isrgb (img))
img = rgb2gray (img);
elseif (isgray (img))
## do nothing
elseif (isvector (img) && !issparse (img) && isreal (img) && all (img >= 0))
hist_in = true;
ihist = img;
else
error ("graythresh: input must be an image or an histogram.");
endif
endif
## the "mean" is the simplest of all, we can get rid of it right here
if (strcmpi (algo, "mean"))
varargout{1} = mean (im2double(img)(:));
return
endif
## we only need to do this if the input is an image. If an histogram, is
## supplieed, no need to do any of this
if (!hist_in)
## if image is uint we do nothing. If it's int, then we convert to uint since
## it may mess up calculations later. If it's double we need some bins so we
## choose uint8 by default... but should we adjust the histogram before?
if (isa (img, "uint16") || isa (img, "uint8"))
## do nothing
elseif (isa (img, "int16"))
img = im2uint16 (img);
img_max = 65535;
else
img = im2uint8 (img);
endif
ihist = hist (img(:), 0:intmax (class (img)));
endif
switch tolower (algo)
case {"concavity"}, thresh = concavity (ihist);
case {"intermeans"}, thresh = intermeans (ihist, floor (mean (img (:))));
case {"intermodes"}, thresh = intermodes (ihist);
case {"maxentropy"}, thresh = maxentropy (ihist);
case {"maxlikelihood"}, thresh = maxlikelihood (ihist);
case {"minerror"}, thresh = minerror_iter (ihist, floor (mean (img (:))));
case {"minimum"}, thresh = minimum (ihist);
case {"moments"}, thresh = moments (ihist);
case {"otsu"}, thresh = otsu (ihist, nargout > 1);
case {"percentile"}, thresh = percentile (ihist, varargin{:});
otherwise, error ("graythresh: unknown method '%s'", algo);
endswitch
## normalize the threshold value to the [0 1] range
thresh{1} = double (thresh{1}) / (numel (ihist) - 1);
## some algorithms may return more than one value...
for i = 1:numel (thresh)
varargout{i} = thresh{i};
endfor
endfunction
function [thresh] = otsu (ihist, compute_good)
## this method is quite well explained at
## http://www.labbookpages.co.uk/software/imgProc/otsuThreshold.html
##
## It does not, however, explain how to compute the goodness of threshold and
## there's not many pages explaining it either. For that, one really needs to
## check the paper.
##
## The implemententation on the link above assumes that threshold is to be
## made for values "greater or equal than" but that is not the case (in im2bw
## and also not ImageJ) so we subtract 1 at the end.
bins = 0:(numel (ihist) - 1);
total = sum (ihist);
## b = black, w = white
b_totals = cumsum ([0 ihist(1:end-1)]);
b_weights = b_totals / total;
b_means = [0 cumsum(bins(1:end-1) .* ihist(1:end-1))] ./ b_totals;
w_totals = total - b_totals;
w_weights = w_totals / total;
w_means = (cumsum (bins(end:-1:1) .* ihist(end:-1:1)) ./ w_totals(end:-1:1))(end:-1:1);
## between class variance (its maximum is the best threshold)
bcv = b_weights .* w_weights .* (b_means - w_means).^2;
## in case there's more than one place with best maximum (for example, a group
## of empty bins, we select the one in the center (this is compatible with ImageJ)
thresh{1} = ceil (mean (find (bcv == max (bcv)))) - 2;
## we subtract 2, once for the 1 based indexes and another for the greater
## than or equal problem
if (compute_good)
## basically we need to divide the between class variance by the total
## variance which is a single value, independent of the threshold. From the
## paper, last of the equation 12, eta = sigma²b / sigma²t
##(where b = between and t = total)
norm_hist = ihist / total;
total_mean = sum (bins .* norm_hist);
total_variance = sum (((bins - total_mean).^2) .* norm_hist);
thresh{2} = max (bcv) / total_variance;
endif
endfunction
function level = moments (y)
n = numel (y) - 1;
## The threshold is chosen such that partial_sumA(y,t)/partial_sumA(y,n)
## is closest to x0.
sumY = sum (y);
Avec = cumsum (y) / sumY;
sumB = partial_sumB (y,n);
sumC = partial_sumC (y,n);
sumD = partial_sumD (y,n);
## The following finds x0.
x2 = (sumB*sumC - sumY*sumD) / (sumY*sumC - sumB^2);
x1 = (sumB*sumD - sumC^2) / (sumY*sumC - sumB^2);
x0 = .5 - (sumB/sumY + x2/2) / sqrt (x2^2 - 4*x1);
## And finally the threshold
[~, ind] = min (abs (Avec-x0));
level{1} = ind-1;
endfunction
function T = maxentropy(y)
n = numel (y) - 1;
dz_warn = warning ("query", "divide-by-zero");
unwind_protect
warning ("off", "Octave:divide-by-zero");
## The threshold is chosen such that the following expression is minimized.
sumY = sum (y);
negY = negativeE (y, n);
for j = 0:n
sumA = partial_sumA (y, j);
negE = negativeE (y, j);
sum_diff = sumY - sumA;
vec(j+1) = negE/sumA - log10 (sumA) + (negY-negE)/(sum_diff) - log10 (sum_diff);
end
unwind_protect_cleanup
## restore broadcats warning status
warning (dz_warn.state, "Octave:divide-by-zero");
end_unwind_protect
[~,ind] = min (vec);
T{1} = ind-1;
endfunction
function [T] = intermodes (y)
## checked with ImageJ and is slightly different but not by much
n = numel (y) - 1;
% 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{1} = 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{1} = floor(mean(TT));
endfunction
## The threshold is chosen such that 50% (in case of p = 0.5) of
## pixels lie in each category.
function [T] = percentile (y, p = 0.5)
Avec = cumsum (y) / sum (y);
[~, ind] = min (abs (Avec - p));
T{1} = ind -1;
endfunction
function T = minimum(y)
n = numel (y) - 1;
% 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{1} = 0;
return
end
end
peakfound = false;
for k = 2:n
if y(k-1) < y(k) && y(k+1) < y(k)
peakfound = true;
end
if peakfound && y(k-1) >= y(k) && y(k+1) >= y(k)
T{1} = k-1;
return
end
end
endfunction
function [Tout] = minerror_iter (y, T)
n = numel (y) - 1;
Tprev = NaN;
dz_warn = warning ("query", "divide-by-zero");
unwind_protect
warning ("off", "Octave:divide-by-zero");
sumA = partial_sumA (y, n);
sumB = partial_sumB (y, n);
sumC = partial_sumC (y, n);
while T ~= Tprev
% Calculate some statistics.
sumAT = partial_sumA (y, T);
sumBT = partial_sumB (y, T);
sumCT = partial_sumC (y, T);
sumAdiff = sumA - sumAT;
mu = sumBT/sumAT;
nu = (sumB-sumBT)/(sumAdiff);
p = sumAT/sumA;
q = (sumAdiff) / sumA;
sigma2 = sumCT/sumAT-mu^2;
tau2 = (sumC-sumCT) / (sumAdiff) - 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.')
break
endif
% 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
endwhile
unwind_protect_cleanup
## restore broadcats warning status
warning (dz_warn.state, "Octave:divide-by-zero");
end_unwind_protect
Tout{1} = T;
endfunction
#{
## this is an implementation of the original minerror algorithm but seems
## to converge less often than the iterative version. This one 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);
% 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
vec(vec==-inf) = NaN;
[minimum,ind] = min(vec);
T = ind-1;
endfunction
#}
function Tout = maxlikelihood (y)
n = numel (y) - 1;
## initial estimate for the threshold is found with the minimum algorithm
T = minimum (y){1};
sumY = sum (y);
sumB = partial_sumB (y, n);
sumC = partial_sumC (y, n);
sumAT = partial_sumA (y, T);
sumBT = partial_sumB (y, T);
sumCT = partial_sumC (y, T);
## initial values for the statistics
mu = sumBT / sumAT;
nu = (sumB - sumBT) / (sumY - sumAT);
p = sumAT / sumY;
q = (sumY - sumAT) / sumY;
sigma2 = sumCT / sumAT - mu^2;
tau2 = (sumC - sumCT) / (sumY - sumAT) - nu^2;
## Return if sigma2 or tau2 are zero, to avoid division by zero
if (sigma2 == 0 || tau2 == 0)
Tout{1} = T;
return
endif
do
## we store the previous values for comparison at the end (we will stop when
## they stop changing)
mu_prev = mu;
nu_prev = nu;
p_prev = p;
q_prev = q;
sigma2_prev = nu;
tau2_prev = nu;
for i = 0:n
phi(i+1) = p/sqrt((sigma2)) * exp(-((i-mu)^2) / (2*sigma2)) / ...
(p/sqrt(sigma2) * exp(-((i-mu)^2) / (2*sigma2)) + ...
(q/sqrt(tau2)) * exp(-((i-nu)^2) / (2*tau2)));
endfor
ind = 0:n;
gamma = 1-phi;
F = phi*y';
G = gamma*y';
mu = ind.*phi*y'/F;
nu = ind.*gamma*y'/G;
p = F / sumY;
q = G / sumY;
sigma2 = ind.^2.*phi*y'/F - mu^2;
tau2 = ind.^2.*gamma*y'/G - nu^2;
until (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)
## 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)
Tout{1} = 0;
return
endif
## The threshold is the integer part of the solution of the quadratic equation
Tout{1} = floor((w1+sqrt(sqterm))/w0);
endfunction
function Tout = intermeans (y, T)
n = numel (y) - 1;
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.
sumY = sum (y);
sumB = partial_sumB (y, n);
while T ~= Tprev
sumAT = partial_sumA (y, T);
sumBT = partial_sumB (y, T);
mu = sumBT/sumAT;
nu = (sumB-sumBT)/(sumY-sumAT);
Tprev = T;
T = floor((mu+nu)/2);
end
Tout{1} = T;
endfunction
function T = concavity (h)
n = numel (h) - 1;
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{1} = 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
## these were tested with ImageJ
%!shared img, histo
%! ## this is the old default.img that came with GNU Octave. While the current
%! ## is very very similar, is off just enough for us to get precision errors
%! img = uint8 (reshape ([138 138 138 142 142 138 142 138 138 117 105 81 69 61 53 40 49 45 40 36 40 45 53 49 65 73 121 166 210 243 247 247 247 239 235 178 154 170 150 150 162 174 190 190 194 186 178 170 154 182 198 174 117 138 138 142 138 142 142 146 142 138 138 130 109 97 81 73 69 57 53 53 57 61 61 69 73 77 105 121 158 219 243 243 247 243 243 243 206 150 158 158 158 150 158 182 186 190 194 186 174 190 206 198 162 138 142 138 142 146 138 142 142 138 146 142 134 142 130 121 101 97 85 85 81 81 81 85 93 85 73 57 61 93 150 194 215 239 243 243 243 223 166 138 158 158 154 142 162 178 190 190 198 186 182 186 174 162 182 146 142 138 142 142 146 142 146 146 146 146 142 142 142 134 125 101 85 73 65 69 73 73 57 40 53 49 57 69 85 125 166 182 178 178 174 150 130 121 146 146 150 142 166 182 190 182 174 166 162 170 194 198 138 138 146 146 138 146 146 146 146 142 150 146 146 142 130 93 65 45 45 49 45 40 49 40 49 49 49 49 61 81 113 142 150 154 154 146 142 134 125 125 138 134 125 146 162 178 178 178 166 186 202 206 186 142 142 142 134 142 146 142 150 142 146 142 146 146 130 81 53 49 49 45 49 40 36 36 32 36 36 36 53 73 89 125 150 146 134 138 146 138 146 138 142 117 117 113 117 146 166 174 178 182 178 178 170 146 142 142 138 142 146 142 142 146 150 138 146 142 130 73 49 40 49 57 65 69 73 61 61 53 57 53 61 77 77 97 113 138 134 130 138 142 150 146 150 134 138 121 121 101 121 150 158 154 142 150 162 166 178 138 138 146 142 142 142 142 146 146 142 142 130 73 57 49 36 49 65 77 85 89 85 81 81 81 85 93 93 97 105 117 125 150 158 154 162 162 166 154 134 150 130 125 113 138 182 174 154 130 178 227 239 239 134 138 142 138 142 142 146 146 138 150 125 61 49 32 32 45 49 57 65 85 101 105 101 101 109 125 117 113 109 138 134 125 166 178 170 162 150 170 162 170 150 146 150 138 125 162 186 182 142 206 247 247 243 138 138 138 138 142 142 146 146 146 130 85 45 45 36 40 53 45 57 69 97 125 130 130 134 138 146 142 134 142 158 138 117 146 174 170 174 178 170 174 170 166 154 162 158 130 134 170 178 158 190 243 247 247 142 142 142 142 142 146 146 142 138 89 53 45 40 45 45 49 57 77 93 125 138 150 154 158 158 162 154 150 166 174 142 73 125 174 178 174 182 182 178 178 174 166 174 174 162 125 154 170 174 170 227 247 251 142 138 142 142 142 142 142 138 105 61 40 40 32 40 40 49 61 89 117 146 154 158 162 170 170 174 162 166 174 182 150 65 146 166 174 186 198 198 198 190 178 178 174 174 158 134 154 198 194 174 202 251 251 146 142 142 142 146 150 138 134 69 40 40 36 32 40 45 45 65 101 134 150 158 166 174 178 174 174 174 170 170 174 142 73 150 162 178 194 202 202 194 194 178 178 154 134 125 138 154 198 194 186 190 243 251 150 146 146 146 146 150 130 109 53 45 28 40 40 36 32 49 73 101 130 154 162 170 170 170 178 182 178 178 174 158 142 121 146 158 178 174 186 190 186 186 174 146 105 109 113 130 150 178 202 190 186 243 251 146 146 146 146 150 142 109 73 49 40 32 40 40 45 40 53 69 93 130 154 162 170 174 178 182 182 186 182 178 154 146 130 138 142 150 170 182 178 174 166 150 117 97 105 113 130 150 150 174 182 190 243 251 146 146 154 146 150 134 105 53 40 45 45 40 40 36 36 40 69 105 134 162 170 174 178 182 182 182 186 190 186 178 170 158 154 150 162 182 182 174 174 174 150 113 109 113 113 130 150 162 186 186 190 239 251 154 150 146 150 146 125 77 49 36 40 36 40 36 28 40 36 77 113 138 150 170 170 174 186 190 190 190 194 190 186 194 190 170 162 174 194 174 182 170 170 158 121 113 113 113 146 158 170 210 215 215 206 243 150 146 150 150 150 113 57 49 40 45 45 49 49 40 32 45 85 113 142 170 178 174 182 194 190 194 194 198 198 198 210 210 182 162 170 190 182 186 170 170 162 130 121 113 121 146 154 150 198 215 206 210 215 150 150 150 150 150 105 49 45 40 49 49 57 40 49 49 53 85 121 158 182 178 174 182 198 194 194 194 194 202 202 194 186 174 154 162 166 178 174 170 170 170 158 117 113 130 150 154 121 182 194 206 215 206 158 150 150 150 146 97 45 36 49 49 49 40 40 49 49 65 97 130 154 174 174 174 186 194 194 194 194 198 198 186 170 158 154 158 138 158 162 170 190 182 174 170 138 138 142 154 134 142 146 170 206 219 215 150 150 158 158 150 85 36 40 40 40 40 45 45 49 49 65 97 130 146 166 166 174 182 190 194 194 194 194 190 182 162 158 150 158 182 186 178 198 206 198 190 174 154 174 174 142 142 170 170 166 202 223 219 158 150 150 150 146 85 40 45 40 40 36 45 53 45 49 53 93 117 130 154 162 174 190 186 194 194 194 190 186 178 162 162 170 174 182 198 210 206 210 198 198 182 170 178 174 158 154 194 194 174 198 210 215 150 154 158 150 150 85 49 45 40 40 32 36 53 40 45 53 81 109 142 158 158 174 178 182 190 190 194 190 190 178 170 174 178 186 190 190 206 215 202 206 194 186 178 182 174 154 170 198 210 186 186 202 215 150 154 150 154 150 97 45 40 40 40 36 36 45 40 45 73 89 113 142 158 158 174 174 182 186 186 194 186 182 178 174 170 105 166 206 186 190 202 198 194 190 182 182 174 166 154 162 198 215 202 182 202 219 154 150 154 150 146 117 61 45 45 45 36 53 53 49 53 77 93 101 125 158 162 174 174 178 174 186 190 182 182 186 182 182 77 125 198 194 186 190 190 178 178 178 162 162 162 154 186 210 227 210 190 206 223 154 150 154 150 154 138 65 45 45 45 40 49 49 40 53 65 77 89 113 150 158 166 166 170 178 182 186 182 170 170 170 162 81 117 186 190 186 182 178 186 174 166 162 150 130 154 194 227 227 219 202 202 219 154 154 150 154 146 146 89 45 40 45 40 49 49 36 40 57 65 89 109 138 146 158 158 170 170 178 182 178 162 150 158 154 113 146 186 182 178 182 178 170 170 162 146 138 138 146 202 223 231 219 210 190 215 130 130 130 130 130 130 109 45 53 40 32 36 40 45 53 61 65 81 97 117 130 138 150 158 158 178 170 162 158 138 142 150 146 166 178 174 174 170 170 170 162 158 138 117 117 142 202 223 239 223 215 186 206 61 61 65 69 69 65 57 36 40 36 32 40 40 53 57 53 57 69 93 105 109 130 138 142 154 162 150 138 142 125 121 150 162 170 170 166 170 170 170 166 162 138 121 113 130 170 202 223 227 231 202 178 182 45 49 45 40 40 40 45 45 45 45 36 40 32 49 61 61 57 65 73 81 101 109 121 130 142 146 121 89 93 117 113 134 154 174 166 162 166 170 170 162 154 150 142 150 223 186 194 215 231 227 206 182 174 49 40 45 45 49 49 45 49 49 49 49 40 36 45 57 69 65 61 65 69 85 93 109 109 117 109 89 57 57 81 97 113 154 162 166 162 170 158 158 162 154 162 174 231 239 178 186 210 231 239 210 194 178 49 36 49 45 49 49 49 45 45 49 49 36 40 40 45 36 53 53 53 57 57 69 69 73 69 61 57 45 45 65 89 105 125 142 146 150 150 154 162 170 174 223 235 247 231 178 178 206 227 227 223 198 190 40 53 36 45 40 40 40 40 45 40 40 45 45 45 45 40 53 49 49 45 53 45 32 36 36 36 36 40 49 45 61 73 89 93 97 113 125 142 186 202 239 239 243 251 239 198 166 194 215 235 227 215 202 40 45 36 32 36 40 40 45 40 40 45 49 45 49 45 49 40 40 45 49 40 45 45 45 49 49 32 40 49 40 49 57 69 81 101 134 170 206 235 243 243 239 247 251 247 210 170 186 202 231 231 227 210 49 45 49 40 40 40 49 45 40 40 45 45 45 40 45 45 45 49 40 49 40 49 45 45 36 40 40 45 45 45 45 65 121 150 210 239 243 243 247 243 243 247 251 251 239 223 178 174 194 219 239 231 219 36 45 45 40 40 49 40 45 49 49 40 40 45 49 40 40 45 49 45 40 49 45 40 40 40 49 40 45 40 49 49 121 162 215 247 247 247 247 247 243 247 251 251 251 247 239 223 194 186 202 215 210 210 36 45 45 40 40 49 40 45 32 36 49 36 45 49 40 40 45 40 36 40 45 45 40 40 40 36 45 32 40 49 57 121 142 215 243 247 243 247 243 247 251 251 251 251 247 247 247 227 186 194 190 190 182 40 32 45 32 45 40 45 45 49 45 40 45 49 36 40 45 32 40 45 45 49 45 45 45 45 53 49 53 45 45 40 69 97 186 239 243 247 247 247 251 251 251 251 251 243 243 231 202 202 206 206 186 170 53 40 40 40 40 40 36 32 32 36 45 53 49 32 36 32 36 32 40 49 40 40 45 40 40 53 45 49 49 40 32 40 49 138 219 235 247 247 251 251 251 251 251 247 243 235 198 206 210 198 190 186 186 73 69 61 57 61 49 53 40 49 45 40 49 49 49 57 57 53 49 53 53 45 40 45 40 45 49 45 49 45 40 32 53 69 101 215 231 247 247 247 247 251 251 251 243 235 219 194 202 202 186 186 190 194], [53 40]));
%!assert (graythresh (img, "percentile"), 142/255);
%!assert (graythresh (img, "percentile", 0.5), 142/255);
%!assert (graythresh (img, "moments"), 142/255);
%!assert (graythresh (img, "minimum"), 93/255);
%!assert (graythresh (img, "maxentropy"), 150/255);
%!assert (graythresh (img, "intermodes"), 99/255);
%!assert (graythresh (img, "otsu"), 115/255);
%! histo = hist (img(:), 0:255);
%!assert (graythresh (histo, "otsu"), 115/255);
## for the mean our results differ from matlab because we do not calculate it
## from the histogram. Our results should be more accurate.
%!assert (graythresh (img, "mean"), 0.51445615982, 0.000000001); # here our results differ from ImageJ
|