Skip to content

Commit 7da853e

Browse files
dleglanddlegland
committed
add geodesic distance map
1 parent 9864671 commit 7da853e

File tree

4 files changed

+431
-3
lines changed

4 files changed

+431
-3
lines changed

‎src/@Image/distanceMap.m‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -12,11 +12,11 @@
1212
% show(map)
1313
%
1414
% See also
15-
%
15+
% skeleton, geodesicDistanceMap
1616

1717
% ------
1818
% Author: David Legland
19-
% e-mail: david.legland@inra.fr
19+
% e-mail: david.legland@inrae.fr
2020
% Created: 2011-03-27, using Matlab 7.9.0.529 (R2009b)
2121
% Copyright 2011 INRA - Cepia Software Platform.
2222

‎src/@Image/geodesicDistanceMap.m‎

Lines changed: 322 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,322 @@
1+
function dist = geodesicDistanceMap(marker, mask, varargin)
2+
% Geodesic distance transform for binary or label images.
3+
%
4+
% RES = geodesicDistanceMap(MARKER, MASK);
5+
% where MASK and MARKER are binary images the same size, computes for
6+
% each foreground pixel the minimum distance to the marker, using a path
7+
% that is contained in the foreground. If the markers can not be reached
8+
% from a foreground pixel, the corresponding result is Inf.
9+
% The function propagates distances using chamfer weights, that
10+
% approximate the Euclidean distance to neighbor pixels.
11+
%
12+
% RES = geodesicDistanceMap(..., WEIGHTS);
13+
% Specifies different weights for computing distance between 2 pixels.
14+
% WEIGHTS is a 2 elements array, with WEIGHTS(1) corresponding to the
15+
% distance between two orthonal pixels, and WEIGHTS(2) corresponding to
16+
% the distance between two diagonal pixels.
17+
% Possible choices
18+
% WEIGHTS = [5 7 11] -> best choice for 5x5 chamfer masks (default)
19+
% WEIGHTS = [1 sqrt(2)] -> quasi-euclidean distance
20+
% WEIGHTS = [1 Inf] -> "Manhattan" or "cityblock" distance
21+
% WEIGHTS = [1 1] -> "Chessboard" distance
22+
% WEIGHTS = [3 4] -> Borgerfors' weights
23+
% WEIGHTS = [5 7] -> close approximation of sqrt(2)
24+
%
25+
% Note: when specifying weights, the result has the same class/data type
26+
% than the array of weights. It is possible to balance between speed and
27+
% memory usage:
28+
% - if weights are double (the default), the memory usage is high, but
29+
% the result can be given in pixel units
30+
% - if weights are integer (for Borgefors weights, for example), the
31+
% memory usage is reduced, but representation limit of datatype can
32+
% be easily reached. One needs to divide by the first weight to get
33+
% result comparabale with natural distances.
34+
% For uint8, using [3 4] weigths, the maximal computable distance is
35+
% around 255/3 = 85 pixels. Using 'int16' seems to be a good
36+
% tradeoff, the maximal distance with [3 4] weights is around 11000
37+
% pixels.
38+
%
39+
% RES = geodesicDistanceMap(..., 'verbose', true);
40+
% Displays info on iterations.
41+
%
42+
% Example
43+
% % computes distance function inside a complex particle
44+
% mask = Image.read('circles.png');
45+
% marker = Image.false(size(mask));
46+
% marker(80, 80) = 1;
47+
% % compute using quasi-enclidean weights
48+
% dist = geodesicDistanceMap(marker, mask);
49+
% figure; imshow(dist, []);
50+
% colormap(jet); title('Quasi-euclidean distance');
51+
% % compute using integer weights, giving integer results
52+
% dist34 = geodesicDistanceMap(marker, mask, int16([3 4]));
53+
% figure; imshow(double(dist34)/3, [0 max(dist34(mask))/3]);
54+
% colormap(jet); title('Borgefors 3-4 weights');
55+
%
56+
%
57+
% % uses the examples from bwdist with different distances
58+
% img = ones(255, 255);
59+
% img(126, 126) = 0;
60+
% res1 = geodesicDistanceMap(img);
61+
% res2 = geodesicDistanceMap(img, [1 inf]);
62+
% res3 = geodesicDistanceMap(img, [1 1]);
63+
% res4 = geodesicDistanceMap(img, [1 1.5]);
64+
% figure;
65+
% subplot(221); subimage(mat2gray(res1));
66+
% hold on; imcontour(res1); title('quasi-euclidean');
67+
% subplot(222); subimage(mat2gray(res2));
68+
% hold on; imcontour(res2); title('city-block');
69+
% subplot(223); subimage(mat2gray(res3));
70+
% hold on; imcontour(res3); title('chessboard');
71+
% subplot(224); subimage(mat2gray(res4));
72+
% hold on; imcontour(res4); title('approx euclidean');
73+
%
74+
% The function uses scanning algorithm. Each iteration consists in a
75+
% sequence of a forward and a backward scan. Iterations stop when
76+
% stability is reached.
77+
%
78+
% See also
79+
% distanceMap, reconstruction
80+
%
81+
82+
% ------
83+
% Author: David Legland
84+
% e-mail: david.legland@inrae.fr
85+
% INRAE - BIA Research Unit - BIBS Platform (Nantes)
86+
% Created: 2020-10-19, using Matlab 9.8.0.1323502 (R2020a)
87+
% Copyright 2020 INRAE.
88+
89+
90+
%% Default options
91+
92+
% default weights for orthogonal, diagonal, and chess-knight neighbors
93+
weights = [5 7 11];
94+
95+
normalize = true;
96+
97+
% silent processing by default
98+
verbose = false;
99+
100+
101+
%% Process input arguments
102+
103+
% extract user-specified weights
104+
if ~isempty(varargin) && isnumeric(varargin{1})
105+
weights = varargin{1};
106+
varargin(1) = [];
107+
end
108+
109+
% extract verbosity option
110+
if length(varargin) > 1
111+
varName = varargin{1};
112+
if ~ischar(varName)
113+
error('Require options as name-value pairs');
114+
end
115+
116+
if strcmpi(varName, 'normalize')
117+
normalize = varargin{2};
118+
elseif strcmpi(varName, 'verbose')
119+
verbose = varargin{2};
120+
else
121+
error(['unknown option: ' varName]);
122+
end
123+
end
124+
125+
126+
%% Pre-Processing
127+
128+
if size(marker) ~= size(mask)
129+
error('Marker and mask image must have same size');
130+
end
131+
132+
% extract weights in specific directions
133+
w1 = weights(1);
134+
w2 = weights(2);
135+
136+
% small check up to avoid degenerate cases
137+
if w2 == 0
138+
w2 = 2 * w1;
139+
end
140+
141+
% initialize weight associated to chess-knight move
142+
if length(weights) == 2
143+
nShifts = 4;
144+
% the list of shifts for forward scans (shift in y, shift in x, weight)
145+
fwdShifts = [...
146+
-1 -1; ...
147+
-1 0; ...
148+
-1 +1; ...
149+
0 -1];
150+
151+
% the list of shifts for backward scans (shift in y, shift in x, weight)
152+
bckShifts = [...
153+
+1 +1; ...
154+
+1 0; ...
155+
+1 -1; ...
156+
0 +1];
157+
ws = [w2 w1 w2 w1];
158+
else
159+
w3 = weights(3);
160+
nShifts = 8;
161+
162+
% the list of shifts for forward scans (shift in y, shift in x, weight)
163+
fwdShifts = [...
164+
-2 -1; ...
165+
-2 +1; ...
166+
-1 -2; ...
167+
-1 -1; ...
168+
-1 0; ...
169+
-1 +1; ...
170+
-1 +2; ...
171+
0 -1 ];
172+
173+
% the list of shifts for backward scans (shift in y, shift in x, weight)
174+
bckShifts = [...
175+
+2 +1; ...
176+
+2 -1; ...
177+
+1 +2; ...
178+
+1 +1; ...
179+
+1 0; ...
180+
+1 -1; ...
181+
+1 -2; ...
182+
0 +1];
183+
ws = [w3 w3 w3 w2 w1 w2 w3 w1];
184+
end
185+
186+
% determines type of output from type of weights
187+
outputType = class(w1);
188+
189+
% allocate memory for result
190+
dist = ones(size(mask), outputType);
191+
192+
% init result: either max value, or 0 for marker pixels
193+
if isinteger(w1)
194+
dist(:) = intmax(outputType);
195+
else
196+
dist(:) = inf;
197+
end
198+
dist(marker) = 0;
199+
200+
% size of image
201+
[D1, D2] = size(mask);
202+
203+
204+
%% Iterations until no more changes
205+
206+
% apply forward and backward iteration until no more changes are made
207+
modif = true;
208+
nIter = 1;
209+
while modif
210+
modif = false;
211+
212+
%% Forward iteration
213+
214+
if verbose
215+
disp(sprintf('Forward iteration %d', nIter)); %#ok<DSPS>
216+
end
217+
218+
% forward iteration on lines
219+
for i = 1:D1
220+
% process all pixels inside the current line
221+
for j = 1:D2
222+
% computes only for pixel inside mask
223+
if mask.Data(i, j) == 0
224+
continue;
225+
end
226+
227+
% compute minimal propagated distance around neighbors in
228+
% forward mask
229+
newVal = dist(i, j);
230+
for k = 1:nShifts
231+
% coordinate of neighbor
232+
i2 = i + fwdShifts(k, 1);
233+
j2 = j + fwdShifts(k, 2);
234+
235+
% check bounds of neighbor
236+
if i2 < 1 || i2 > D1
237+
continue;
238+
end
239+
if j2 < 1 || j2 > D2
240+
continue;
241+
end
242+
243+
% compute new value of local distance map
244+
if mask.Data(i2, j2) == mask.Data(i, j)
245+
% neighbor in same region
246+
% -> add offset weight to neighbor distance
247+
newVal = min(newVal, dist(i2, j2) + ws(k));
248+
end
249+
end
250+
251+
% if distance was changed, update result, and toggle flag
252+
if newVal < dist(i,j)
253+
modif = true;
254+
dist(i,j) = newVal;
255+
end
256+
end
257+
258+
end % iteration on lines
259+
260+
% check end of iteration
261+
if ~modif && nIter ~= 1
262+
break;
263+
end
264+
265+
%% Backward iteration
266+
modif = false;
267+
268+
if verbose
269+
disp(sprintf('Backward iteration %d', nIter)); %#ok<DSPS>
270+
end
271+
272+
% backward iteration on lines
273+
for i = D1:-1:1
274+
% process all pixels inside the current line
275+
for j = D2:-1:1
276+
% computes only for pixel inside mask
277+
if mask.Data(i, j) == 0
278+
continue;
279+
end
280+
281+
% compute minimal propagated distance around neighbors in
282+
% backward mask
283+
newVal = dist(i, j);
284+
for k = 1:nShifts
285+
% coordinate of neighbor
286+
i2 = i + bckShifts(k, 1);
287+
j2 = j + bckShifts(k, 2);
288+
289+
% check bounds of neighbor
290+
if i2 < 1 || i2 > D1
291+
continue;
292+
end
293+
if j2 < 1 || j2 > D2
294+
continue;
295+
end
296+
297+
% compute new value of local distance map
298+
if mask.Data(i2, j2) == mask.Data(i, j)
299+
% neighbor in same region
300+
% -> add offset weight to neighbor distance
301+
newVal = min(newVal, dist(i2, j2) + ws(k));
302+
end
303+
end
304+
305+
% if distance was changed, update result, and toggle flag
306+
if newVal < dist(i,j)
307+
modif = true;
308+
dist(i,j) = newVal;
309+
end
310+
end
311+
312+
end % line iteration
313+
314+
nIter = nIter+1;
315+
end % until no more modif
316+
317+
% normalize map
318+
if normalize
319+
dist(mask > 0) = dist(mask > 0) / w1;
320+
end
321+
322+
dist = Image('Data', dist, 'Parent', mask, 'Type', 'Intensity');

‎src/@Image/reconstruction.m‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,8 @@
2020
% show(rec);
2121
%
2222
% See also
23-
% killBorders, fillHoles, extendedMinima, floodFill, imreconstruct
23+
% killBorders, fillHoles, extendedMinima, floodFill, imreconstruct,
24+
% geodesicDistanceMap
2425

2526
% ------
2627
% Author: David Legland

0 commit comments

Comments
 (0)