Skip to content

Commit e8f5207

Browse files
dleglanddlegland
committed
add multilevelOtsuThreshold.m
1 parent 2f471fe commit e8f5207

File tree

3 files changed

+226
-0
lines changed

3 files changed

+226
-0
lines changed

‎src/@Image/multilevelOtsuThreshold.m‎

Lines changed: 187 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,187 @@
1+
function [regions, thresholds, maxSig] = multilevelOtsuThreshold(obj, nClasses, varargin)
2+
% Multilevel Thresholding using Otsu Method.
3+
%
4+
% CLASSES = multilevelOtsuThreshold(IMG, NC)
5+
% Computes a segmentation of the input grayscale image IMG into NC
6+
% classes. The number of classes must be comprised between 2 and 5.
7+
%
8+
% [CLASSES, THRESHOLDS] = multilevelOtsuThreshold(IMG, NC)
9+
% Also returns the list of threshold. The number of threshold values is
10+
% the number of classes minus one.
11+
%
12+
% Example
13+
% % segment cameraman image into three classes
14+
% img = Image.read('cameraman.tif');
15+
% [classes, threshs] = multilevelOtsuThreshold(img, 3);
16+
% rgb = label2rgb(classes, 'jet', [1 0 0], 'shuffle');
17+
% figure; show(rgb);
18+
% % also displays histogram with threshold levels
19+
% figure; histogram(img);
20+
% hold on;
21+
% for i = 1:2
22+
% plot(threshs([i i]), [0 1800], 'color', 'r', 'linewidth', 2);
23+
% end
24+
%
25+
% References
26+
% * Ping-Sung Liao, Tse-Sheng Chen, Pau-Choo Chung (2001). "A Fast
27+
% Algorithm for Multilevel Thresholding". Journal of Information Science
28+
% and Engineering, Vol. 17 No. 5, pp. 713-727.
29+
% https://jise.iis.sinica.edu.tw/JISESearch/pages/View/PaperView.jsf?keyId=86_1302#
30+
% * https://imagej.net/Multi_Otsu_Threshold
31+
%
32+
% See also
33+
% otsuThreshold, histogram, maxEntropyThreshold
34+
%
35+
36+
% ------
37+
% Author: David Legland
38+
% e-mail: david.legland@inrae.fr
39+
% INRAE - BIA Research Unit - BIBS Platform (Nantes)
40+
% Created: 2021-02-02, using Matlab 9.8.0.1323502 (R2020a)
41+
% Copyright 2021 INRAE.
42+
43+
44+
%% Input arguments
45+
46+
if ~isGrayscaleImage(obj)
47+
error('Requires a grayscale uint8 image as first input');
48+
end
49+
50+
if nClasses < 2 || nClasses > 5
51+
error('The number of classes must be comprised between 2 and 5.');
52+
end
53+
54+
55+
%% Initialisations
56+
57+
% compute histogram, and convert into probability density
58+
[h, levels] = histogram(obj, varargin{:});
59+
h = h / sum(h);
60+
61+
% number of gray levels
62+
nLevels = length(levels);
63+
64+
Puv = zeros(nLevels, nLevels);
65+
Suv = zeros(nLevels, nLevels);
66+
Huv = zeros(nLevels, nLevels);
67+
68+
% initialize diagonal terms
69+
for i = 1:nLevels
70+
Puv(i, i) = h(i);
71+
Suv(i, i) = h(i) * i;
72+
end
73+
74+
% initialize first rows
75+
for i = 2:nLevels
76+
Puv(1, i) = Puv(1, i-1) + h(i);
77+
Suv(1, i) = Suv(1, i-1) + i*h(i);
78+
end
79+
% initialize the rest of the matrices
80+
for u = 2:nLevels
81+
for v = u+1:nLevels
82+
Puv(u, v) = Puv(1, v) - Puv(1, u-1);
83+
Suv(u, v) = Suv(1, v) - Suv(1, u-1);
84+
end
85+
end
86+
87+
% now calculate Huv
88+
for u = 1:nLevels
89+
for v = u+1:nLevels
90+
if Puv(u,v) > 0
91+
Huv(u,v) = Suv(u,v) * Suv(u,v) / Puv(u,v);
92+
end
93+
end
94+
end
95+
96+
97+
%% Main processing
98+
99+
% create array for threshold values (first one is zeros, last one is inf)
100+
thresholds = zeros(nClasses+1, 1);
101+
thresholds(end) = inf;
102+
103+
% iterate over all possible combinations
104+
maxSig = 0.0;
105+
switch nClasses
106+
case 2
107+
% two classes -> only need to find second threshold value
108+
for i = 1:nLevels-1
109+
Sq = Huv(1,i) + Huv(i+1, end);
110+
if Sq >= maxSig
111+
thresholds(2) = i;
112+
maxSig = Sq;
113+
end
114+
end
115+
116+
case 3
117+
% three classes
118+
for i = 1:nLevels-2
119+
for j = i+1:nLevels-1
120+
Sq = Huv(1,i) + Huv(i+1, j) + Huv(j+1, end);
121+
if Sq >= maxSig
122+
thresholds(2) = i;
123+
thresholds(3) = j;
124+
maxSig = Sq;
125+
end
126+
end
127+
end
128+
129+
case 4
130+
% four classes
131+
for i = 1:nLevels-3
132+
for j = i+1:nLevels-2
133+
for k = j+1:nLevels-1
134+
Sq = Huv(1,i) + Huv(i+1, j) + Huv(j+1, k) + Huv(k+1, end);
135+
if Sq >= maxSig
136+
thresholds(2) = i;
137+
thresholds(3) = j;
138+
thresholds(4) = k;
139+
maxSig = Sq;
140+
end
141+
end
142+
end
143+
end
144+
145+
case 5
146+
% five classes
147+
for i = 1:nLevels-3
148+
for j = i+1:nLevels-2
149+
for k = j+1:nLevels-1
150+
for m = k+1:nLevels-1
151+
Sq = Huv(1,i) + Huv(i+1, j) + Huv(j+1, k) + Huv(k+1, m) + Huv(m+1, end);
152+
if Sq >= maxSig
153+
thresholds(2) = i;
154+
thresholds(3) = j;
155+
thresholds(4) = k;
156+
thresholds(5) = m;
157+
maxSig = Sq;
158+
end
159+
end
160+
end
161+
end
162+
end
163+
164+
otherwise
165+
error('Can not manage %d number of classes', nClasses);
166+
end
167+
168+
169+
%% Create region image
170+
171+
% compute array of region labels
172+
regData = zeros(size(obj), 'uint8');
173+
for i = 1:nClasses
174+
inds = obj >= thresholds(i) & obj < thresholds(i+1);
175+
regData(inds) = i;
176+
end
177+
178+
% create result image
179+
newName = createNewName(obj, '%s-OtsuClasses');
180+
regions = Image(...
181+
'Data', regData, ...
182+
'Parent', obj, ...
183+
'Type', 'Label', ...
184+
'Name', newName);
185+
186+
% format threshold result
187+
thresholds([1 end]) = [];

‎src/@Image/otsuThreshold.m‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,7 @@
2626
%
2727
% See also
2828
% binarize, otsuThresholdValue, maxEntropyThreshold, lt, le, histogram
29+
% multilevelOtsuThreshold
2930
%
3031

3132
% ------

‎tests/image/test_multilevelOtsuThreshold.m‎

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,38 @@
1+
function tests = test_multilevelOtsuThreshold
2+
% Test suite for the file multilevelOtsuThreshold.
3+
%
4+
% Test suite for the file multilevelOtsuThreshold
5+
%
6+
% Example
7+
% test_multilevelOtsuThreshold
8+
%
9+
% See also
10+
% multilevelOtsuThreshold
11+
12+
% ------
13+
% Author: David Legland
14+
% e-mail: david.legland@inrae.fr
15+
% Created: 2021-02-02, using Matlab 9.8.0.1323502 (R2020a)
16+
% Copyright 2021 INRAE - BIA-BIBS.
17+
18+
tests = functiontests(localfunctions);
19+
20+
function test_Simple(testCase) %#ok<*DEFNU>
21+
% Test call of function without argument.
22+
23+
img = Image.read('cameraman.tif');
24+
25+
classes = multilevelOtsuThreshold(img, 3);
26+
27+
assertEqual(testCase, size(classes), size(img));
28+
29+
30+
function test_getLevels(testCase) %#ok<*DEFNU>
31+
% Test call of function without argument.
32+
33+
img = Image.read('cameraman.tif');
34+
35+
[classes, levels] = multilevelOtsuThreshold(img, 3);
36+
37+
assertEqual(testCase, size(classes), size(img));
38+
assertEqual(testCase, length(levels), 2);

0 commit comments

Comments
 (0)