Skip to content

Commit edf3ddc

Browse files
dleglanddlegland
committed
add chamferDistanceMap.m
1 parent 2aebf69 commit edf3ddc

File tree

3 files changed

+288
-3
lines changed

3 files changed

+288
-3
lines changed

‎src/@Image/chamferDistanceMap.m‎

Lines changed: 203 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,203 @@
1+
function map = chamferDistanceMap(obj, varargin)
2+
% Distance map of a binary image computed using chamfer mask.
3+
%
4+
% DISTMAP = chamferDistanceMap(IMG)
5+
% DISTMAP = chamferDistanceMap(IMG, WEIGHTS)
6+
% Computes the distance map of the input image using chamfer weights. The
7+
% aim of this function is similar to that of the "distanceMap" one, with
8+
% the following specificities:
9+
% * possibility to use 5-by-5 chamfer masks
10+
% * possibility to compute distance maps for label images with touching
11+
% regions.
12+
%
13+
% Example
14+
% chamferDistanceMap
15+
%
16+
% See also
17+
% distanceMap, geodesicDistanceMap
18+
%
19+
20+
% ------
21+
% Author: David Legland
22+
% e-mail: david.legland@inrae.fr
23+
% INRAE - BIA Research Unit - BIBS Platform (Nantes)
24+
% Created: 2021-11-18, using Matlab 9.10.0.1684407 (R2021a) Update 3
25+
% Copyright 2021 INRAE.
26+
27+
%% Process input arguments
28+
29+
% default weights for orthogonal or diagonal
30+
weights = [5 7 11];
31+
32+
normalize = true;
33+
34+
% extract user-specified weights
35+
if ~isempty(varargin)
36+
weights = varargin{1};
37+
varargin(1) = [];
38+
end
39+
40+
% extract verbosity option
41+
verbose = false;
42+
if length(varargin) > 1
43+
varName = varargin{1};
44+
if ~ischar(varName)
45+
error('Require options as name-value pairs');
46+
end
47+
48+
if strcmpi(varName, 'normalize')
49+
normalize = varargin{2};
50+
elseif strcmpi(varName, 'verbose')
51+
verbose = varargin{2};
52+
else
53+
error(['unknown option: ' varName]);
54+
end
55+
end
56+
57+
58+
%% Initialisations
59+
60+
% determines type of output from type of weights
61+
outputType = class(weights);
62+
63+
% small check up to avoid degenerate cases
64+
w1 = weights(1);
65+
w2 = weights(2);
66+
if w2 < w1
67+
w2 = 2 * w1;
68+
end
69+
70+
% shifts in directions i and j for (1) forward and (2) backward iterations
71+
if length(weights) == 2
72+
nShifts = 4;
73+
di1 = [-1 -1 -1 0];
74+
dj1 = [-1 0 1 -1];
75+
di2 = [+1 +1 +1 0];
76+
dj2 = [-1 0 1 +1];
77+
ws = [w2 w1 w2 w1];
78+
79+
elseif length(weights) == 3
80+
nShifts = 8;
81+
w3 = weights(3);
82+
di1 = [-2 -2 -1 -1 -1 -1 -1 0];
83+
dj1 = [-1 +1 -2 -1 0 1 +2 -1];
84+
di2 = [+2 +2 +1 +1 +1 +1 +1 0];
85+
dj2 = [-1 +1 +2 +1 0 -1 -2 +1];
86+
ws = [w3 w3 w3 w2 w1 w2 w3 w1];
87+
end
88+
89+
% allocate memory for result
90+
dist = ones(size(obj.Data), outputType);
91+
92+
% init result: either max value, or 0 for marker pixels
93+
if isinteger(w1)
94+
dist(:) = intmax(outputType);
95+
else
96+
dist(:) = inf;
97+
end
98+
dist(obj.Data == 0) = 0;
99+
100+
% size of image
101+
[D1, D2] = size(obj.Data);
102+
103+
104+
%% Forward iteration
105+
106+
if verbose
107+
disp('Forward iteration %d');
108+
end
109+
110+
for i = 1:D1
111+
for j = 1:D2
112+
% computes only for pixels within a region
113+
if obj.Data(i, j) == 0
114+
continue;
115+
end
116+
117+
% compute minimal propagated distance
118+
newVal = dist(i, j);
119+
for k = 1:nShifts
120+
% coordinate of neighbor
121+
i2 = i + di1(k);
122+
j2 = j + dj1(k);
123+
124+
% check bounds
125+
if i2 < 1 || i2 > D1 || j2 < 1 || j2 > D2
126+
continue;
127+
end
128+
129+
% compute new value
130+
if obj.Data(i2, j2) == obj.Data(i, j)
131+
% neighbor in same region
132+
% -> add offset weight to neighbor distance
133+
newVal = min(newVal, dist(i2, j2) + ws(k));
134+
else
135+
% neighbor in another region
136+
% -> initialize with the offset weight
137+
newVal = min(newVal, ws(k));
138+
end
139+
end
140+
141+
% if distance was changed, update result
142+
dist(i,j) = newVal;
143+
end
144+
145+
end % iteration on lines
146+
147+
148+
149+
%% Backward iteration
150+
151+
if verbose
152+
disp('Backward iteration');
153+
end
154+
155+
for i = D1:-1:1
156+
for j = D2:-1:1
157+
% computes only for foreground pixels
158+
if obj.Data(i, j) == 0
159+
continue;
160+
end
161+
162+
% compute minimal propagated distance
163+
newVal = dist(i, j);
164+
for k = 1:nShifts
165+
% coordinate of neighbor
166+
i2 = i + di2(k);
167+
j2 = j + dj2(k);
168+
169+
% check bounds
170+
if i2 < 1 || i2 > D1 || j2 < 1 || j2 > D2
171+
continue;
172+
end
173+
174+
% compute new value
175+
if obj.Data(i2, j2) == obj.Data(i, j)
176+
% neighbor in same region
177+
% -> add offset weight to neighbor distance
178+
newVal = min(newVal, dist(i2, j2) + ws(k));
179+
else
180+
% neighbor in another region
181+
% -> initialize with the offset weight
182+
newVal = min(newVal, ws(k));
183+
end
184+
end
185+
186+
% if distance was changed, update result
187+
dist(i,j) = newVal;
188+
end
189+
190+
end % line iteration
191+
192+
if normalize
193+
dist(obj.Data>0) = dist(obj.Data>0) / w1;
194+
end
195+
196+
newName = createNewName(obj, '%s-distMap');
197+
198+
% create new image
199+
map = Image('Data', dist, ...
200+
'Parent', obj, ...
201+
'Name', newName, ...
202+
'Type', 'intensity', ...
203+
'ChannelNames', {'distance'});

‎src/@Image/distanceMap.m‎

Lines changed: 8 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -3,22 +3,27 @@
33
%
44
% MAP = distanceMap(BIN)
55
% The distance transform is an operator applied to binary images, that
6-
% results in a graylevel image that contains, for each foregournd pixel,
6+
% results in an intensity image that contains, for each foreground pixel,
77
% the distance to the closest background pixel.
88
%
9+
% This function requires binary image as input. For label images that
10+
% contain adjacent regions, the function 'chamferDistanceMap' could be
11+
% more adapted.
12+
%
913
% Example
1014
% img = Image.read('circles.png');
1115
% map = distanceMap(img);
1216
% show(map)
1317
%
1418
% See also
15-
% skeleton, geodesicDistanceMap
19+
% skeleton, geodesicDistanceMap, chamferDistanceMap
1620

1721
% ------
1822
% Author: David Legland
1923
% e-mail: david.legland@inrae.fr
24+
% INRAE - BIA Research Unit - BIBS Platform (Nantes)
2025
% Created: 2011-03-27, using Matlab 7.9.0.529 (R2009b)
21-
% Copyright 2011 INRA - Cepia Software Platform.
26+
% Copyright 2021 INRAE.
2227

2328
% check type
2429
if ~strcmp(obj.Type, 'binary')

‎tests/image/test_chamferDistanceMap.m‎

Lines changed: 77 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,77 @@
1+
function tests = test_chamferDistanceMap
2+
% Test suite for the file chamferDistanceMap.
3+
%
4+
% Test suite for the file chamferDistanceMap
5+
%
6+
% Example
7+
% test_chamferDistanceMap
8+
%
9+
% See also
10+
% chamferDistanceMap
11+
12+
% ------
13+
% Author: David Legland
14+
% e-mail: david.legland@inrae.fr
15+
% Created: 2021-11-18, using Matlab 9.10.0.1684407 (R2021a) Update 3
16+
% Copyright 2021 INRAE - BIA-BIBS.
17+
18+
tests = functiontests(localfunctions);
19+
20+
function test_SimpleBinary(testCase) %#ok<*DEFNU>
21+
% Test call of function without argument.
22+
23+
% generate an image of a 10x10 square, with one pixel border
24+
img = Image.false([12, 12]);
25+
img(2:11, 2:11) = 1;
26+
27+
% compute distance map
28+
distMap = chamferDistanceMap(img);
29+
30+
assertEqual(testCase, size(img), size(distMap));
31+
% distance equal to 1 on the borders of the square
32+
assertEqual(testCase, 1, distMap(2, 5));
33+
assertEqual(testCase, 1, distMap(11, 5));
34+
assertEqual(testCase, 1, distMap(5, 2));
35+
assertEqual(testCase, 1, distMap(5, 11));
36+
% distance equal to 5 in the middle of the square
37+
assertEqual(testCase, 5, distMap(6, 6));
38+
39+
40+
function test_TouchingLabels(testCase)
41+
% Aim is to compute distance map within each label, even if some of them
42+
% touch each other.
43+
% Uses an image with a completely landlocked label region.
44+
45+
data = [...
46+
0 0 0 0 0 0 0 0 0 0 0; ...
47+
0 1 1 1 2 2 2 3 3 3 0; ...
48+
0 1 1 1 2 2 2 3 3 3 0; ...
49+
0 1 1 1 2 2 2 3 3 3 0; ...
50+
0 1 1 1 4 4 4 3 3 3 0; ...
51+
0 1 1 1 4 4 4 3 3 3 0; ...
52+
0 1 1 1 4 4 4 3 3 3 0; ...
53+
0 1 1 1 5 5 5 3 3 3 0; ...
54+
0 1 1 1 5 5 5 3 3 3 0; ...
55+
0 1 1 1 5 5 5 3 3 3 0; ...
56+
0 0 0 0 0 0 0 0 0 0 0; ...
57+
];
58+
img = Image(data, 'type', 'label');
59+
60+
61+
distMap = chamferDistanceMap(img);
62+
63+
exp = Image([...
64+
0 0 0 0 0 0 0 0 0 0 0; ...
65+
0 1 1 1 1 1 1 1 1 1 0; ...
66+
0 1 2 1 1 2 1 1 2 1 0; ...
67+
0 1 2 1 1 1 1 1 2 1 0; ...
68+
0 1 2 1 1 1 1 1 2 1 0; ...
69+
0 1 2 1 1 2 1 1 2 1 0; ...
70+
0 1 2 1 1 1 1 1 2 1 0; ...
71+
0 1 2 1 1 1 1 1 2 1 0; ...
72+
0 1 2 1 1 2 1 1 2 1 0; ...
73+
0 1 1 1 1 1 1 1 1 1 0; ...
74+
0 0 0 0 0 0 0 0 0 0 0; ...
75+
]);
76+
assertEqual(testCase, size(img), size(distMap));
77+
assertEqual(testCase, distMap.Data, exp.Data);

0 commit comments

Comments
 (0)