[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

distancetransform.hxx
1/************************************************************************/
2/* */
3/* Copyright 1998-2002 by Ullrich Koethe */
4/* */
5/* This file is part of the VIGRA computer vision library. */
6/* The VIGRA Website is */
7/* http://hci.iwr.uni-heidelberg.de/vigra/ */
8/* Please direct questions, bug reports, and contributions to */
9/* ullrich.koethe@iwr.uni-heidelberg.de or */
10/* vigra@informatik.uni-hamburg.de */
11/* */
12/* Permission is hereby granted, free of charge, to any person */
13/* obtaining a copy of this software and associated documentation */
14/* files (the "Software"), to deal in the Software without */
15/* restriction, including without limitation the rights to use, */
16/* copy, modify, merge, publish, distribute, sublicense, and/or */
17/* sell copies of the Software, and to permit persons to whom the */
18/* Software is furnished to do so, subject to the following */
19/* conditions: */
20/* */
21/* The above copyright notice and this permission notice shall be */
22/* included in all copies or substantial portions of the */
23/* Software. */
24/* */
25/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26/* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27/* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28/* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29/* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30/* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31/* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32/* OTHER DEALINGS IN THE SOFTWARE. */
33/* */
34/************************************************************************/
35
36
37#ifndef VIGRA_DISTANCETRANSFORM_HXX
38#define VIGRA_DISTANCETRANSFORM_HXX
39
40#include <cmath>
41#include "stdimage.hxx"
42#include "multi_shape.hxx"
43
44 namespace vigra {
45
46/*
47 * functors to determine the distance norm
48 * these functors assume that dx and dy are positive
49 * (this is OK for use in internalDistanceTransform())
50 */
51
52// chessboard metric
53 struct InternalDistanceTransformLInifinityNormFunctor
54{
55 float operator()(float dx, float dy) const
56 {
57 return (dx < dy) ? dy : dx;
58 }
59};
60
61// Manhattan metric
62 struct InternalDistanceTransformL1NormFunctor
63{
64 float operator()(float dx, float dy) const
65 {
66 return dx + dy;
67 }
68};
69
70// Euclidean metric
71 struct InternalDistanceTransformL2NormFunctor
72{
73 float operator()(float dx, float dy) const
74 {
75 return VIGRA_CSTD::sqrt(dx*dx + dy*dy);
76 }
77};
78
79
80 template <class SrcImageIterator, class SrcAccessor,
81 class DestImageIterator, class DestAccessor,
82 class ValueType, class NormFunctor>
83 void
84internalDistanceTransform(SrcImageIterator src_upperleft,
85 SrcImageIterator src_lowerright, SrcAccessor sa,
86 DestImageIterator dest_upperleft, DestAccessor da,
87 ValueType background, NormFunctor norm)
88{
89 int w = src_lowerright.x - src_upperleft.x;
90 int h = src_lowerright.y - src_upperleft.y;
91
92 FImage xdist(w,h), ydist(w,h);
93
94 xdist = (FImage::value_type)w; // init x and
95 ydist = (FImage::value_type)h; // y distances with 'large' values
96
97 SrcImageIterator sy = src_upperleft;
98 DestImageIterator ry = dest_upperleft;
99 FImage::Iterator xdy = xdist.upperLeft();
100 FImage::Iterator ydy = ydist.upperLeft();
101 SrcImageIterator sx = sy;
102 DestImageIterator rx = ry;
103 FImage::Iterator xdx = xdy;
104 FImage::Iterator ydx = ydy;
105
106 const Diff2D left(-1, 0);
107 const Diff2D right(1, 0);
108 const Diff2D top(0, -1);
109 const Diff2D bottom(0, 1);
110
111 int x,y;
112 if(sa(sx) != background) // first pixel
113 {
114 *xdx = 0.0;
115 *ydx = 0.0;
116 da.set(0.0, rx);
117 }
118 else
119 {
120 da.set(norm(*xdx, *ydx), rx);
121 }
122
123
124 for(x=1, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x;
125 x<w;
126 ++x, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x) // first row left to right
127 {
128 if(sa(sx) != background)
129 {
130 *xdx = 0.0;
131 *ydx = 0.0;
132 da.set(0.0, rx);
133 }
134 else
135 {
136 *xdx = xdx[left] + 1.0f; // propagate x and
137 *ydx = ydx[left]; // y components of distance from left pixel
138 da.set(norm(*xdx, *ydx), rx); // calculate distance from x and y components
139 }
140 }
141 for(x=w-2, xdx.x -= 2, ydx.x -= 2, sx.x -= 2, rx.x -= 2;
142 x>=0;
143 --x, --xdx.x, --ydx.x, --sx.x, --rx.x) // first row right to left
144 {
145 float d = norm(xdx[right] + 1.0f, ydx[right]);
146
147 if(da(rx) < d) continue;
148
149 *xdx = xdx[right] + 1.0f;
150 *ydx = ydx[right];
151 da.set(d, rx);
152 }
153 for(y=1, ++xdy.y, ++ydy.y, ++sy.y, ++ry.y;
154 y<h;
155 ++y, ++xdy.y, ++ydy.y, ++sy.y, ++ry.y) // top to bottom
156 {
157 sx = sy;
158 rx = ry;
159 xdx = xdy;
160 ydx = ydy;
161
162 if(sa(sx) != background) // first pixel of current row
163 {
164 *xdx = 0.0f;
165 *ydx = 0.0f;
166 da.set(0.0, rx);
167 }
168 else
169 {
170 *xdx = xdx[top];
171 *ydx = ydx[top] + 1.0f;
172 da.set(norm(*xdx, *ydx), rx);
173 }
174
175 for(x=1, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x;
176 x<w;
177 ++x, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x) // current row left to right
178 {
179 if(sa(sx) != background)
180 {
181 *xdx = 0.0f;
182 *ydx = 0.0f;
183 da.set(0.0, rx);
184 }
185 else
186 {
187 float d1 = norm(xdx[left] + 1.0f, ydx[left]);
188 float d2 = norm(xdx[top], ydx[top] + 1.0f);
189
190 if(d1 < d2)
191 {
192 *xdx = xdx[left] + 1.0f;
193 *ydx = ydx[left];
194 da.set(d1, rx);
195 }
196 else
197 {
198 *xdx = xdx[top];
199 *ydx = ydx[top] + 1.0f;
200 da.set(d2, rx);
201 }
202 }
203 }
204 for(x=w-2, xdx.x -= 2, ydx.x -= 2, sx.x -= 2, rx.x -= 2;
205 x>=0;
206 --x, --xdx.x, --ydx.x, --sx.x, --rx.x) // current row right to left
207 {
208 float d1 = norm(xdx[right] + 1.0f, ydx[right]);
209
210 if(da(rx) < d1) continue;
211
212 *xdx = xdx[right] + 1.0f;
213 *ydx = ydx[right];
214 da.set(d1, rx);
215 }
216 }
217 for(y=h-2, xdy.y -= 2, ydy.y -= 2, sy.y -= 2, ry.y -= 2;
218 y>=0;
219 --y, --xdy.y, --ydy.y, --sy.y, --ry.y) // bottom to top
220 {
221 sx = sy;
222 rx = ry;
223 xdx = xdy;
224 ydx = ydy;
225
226 float d = norm(xdx[bottom], ydx[bottom] + 1.0f);
227 if(d < da(rx)) // first pixel of current row
228 {
229 *xdx = xdx[bottom];
230 *ydx = ydx[bottom] + 1.0f;
231 da.set(d, rx);
232 }
233
234 for(x=1, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x;
235 x<w;
236 ++x, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x) // current row left to right
237 {
238 float d1 = norm(xdx[left] + 1.0f, ydx[left]);
239 float d2 = norm(xdx[bottom], ydx[bottom] + 1.0f);
240
241 if(d1 < d2)
242 {
243 if(da(rx) < d1) continue;
244 *xdx = xdx[left] + 1.0f;
245 *ydx = ydx[left];
246 da.set(d1, rx);
247 }
248 else
249 {
250 if(da(rx) < d2) continue;
251 *xdx = xdx[bottom];
252 *ydx = ydx[bottom] + 1.0f;
253 da.set(d2, rx);
254 }
255 }
256 for(x=w-2, xdx.x -= 2, ydx.x -= 2, sx.x -= 2, rx.x -= 2;
257 x>=0;
258 --x, --xdx.x, --ydx.x, --sx.x, --rx.x) // current row right to left
259 {
260 float d1 = norm(xdx[right] + 1.0f, ydx[right]);
261
262 if(da(rx) < d1) continue;
263 *xdx = xdx[right] + 1.0f;
264 *ydx = ydx[right];
265 da.set(d1, rx);
266 }
267 }
268}
269
270/********************************************************/
271/* */
272/* distanceTransform */
273/* */
274/********************************************************/
275
276/** \addtogroup DistanceTransform Distance Transform
277*/
278//@{
279
280/** For all background pixels, calculate the distance to
281 the nearest object or contour. The label of the pixels to be considered
282 background in the source image is passed in the parameter 'background'.
283 Source pixels with other labels will be considered objects. In the
284 destination image, all pixels corresponding to background will be assigned
285 the their distance value, all pixels corresponding to objects will be
286 assigned 0.
287
288 The parameter 'norm' gives the distance norm to be used:
289
290 <ul>
291
292 <li> norm == 0: use chessboard distance (L-infinity norm)
293 <li> norm == 1: use Manhattan distance (L1 norm)
294 <li> norm == 2: use Euclidean distance (L2 norm)
295
296 </ul>
297
298 If you use the L2 norm, the destination pixels must be real valued to give
299 correct results.
300
301 <b> Declarations:</b>
302
303 pass 2D array views:
304 \code
305 namespace vigra {
306 template <class T1, class S1,
307 class T2, class S2,
308 class ValueType>
309 void
310 distanceTransform(MultiArrayView<2, T1, S1> const & src,
311 MultiArrayView<2, T2, S2> dest,
312 ValueType background, int norm);
313 }
314 \endcode
315
316 \deprecatedAPI{distanceTransform}
317 pass \ref ImageIterators and \ref DataAccessors :
318 \code
319 namespace vigra {
320 template <class SrcImageIterator, class SrcAccessor,
321 class DestImageIterator, class DestAccessor,
322 class ValueType>
323 void distanceTransform(SrcImageIterator src_upperleft,
324 SrcImageIterator src_lowerright, SrcAccessor sa,
325 DestImageIterator dest_upperleft, DestAccessor da,
326 ValueType background, int norm);
327 }
328 \endcode
329 use argument objects in conjunction with \ref ArgumentObjectFactories :
330 \code
331 namespace vigra {
332 template <class SrcImageIterator, class SrcAccessor,
333 class DestImageIterator, class DestAccessor,
334 class ValueType>
335 void distanceTransform(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
336 pair<DestImageIterator, DestAccessor> dest,
337 ValueType background, int norm);
338 }
339 \endcode
340 \deprecatedEnd
341
342 <b> Usage:</b>
343
344 <b>\#include</b> <vigra/distancetransform.hxx><br>
345 Namespace: vigra
346
347 \code
348 MultiArray<2, unsigned char> src(w,h), edges(w,h);
349 MultiArray<2, float> distance(w, h);
350 ...
351
352 // detect edges in src image (edges will be marked 1, background 0)
353 differenceOfExponentialEdgeImage(src, edges, 0.8, 4.0);
354
355 // find distance of all pixels from nearest edge
356 distanceTransform(edges, distance, 0, 2);
357 // ^ background label ^ norm (Euclidean)
358 \endcode
359
360 \deprecatedUsage{distanceTransform}
361 \code
362
363 vigra::BImage src(w,h), edges(w,h);
364 vigra::FImage distance(w, h);
365
366 // empty edge image
367 edges = 0;
368 ...
369
370 // detect edges in src image (edges will be marked 1, background 0)
371 vigra::differenceOfExponentialEdgeImage(srcImageRange(src), destImage(edges),
372 0.8, 4.0);
373
374 // find distance of all pixels from nearest edge
375 vigra::distanceTransform(srcImageRange(edges), destImage(distance),
376 0, 2);
377 // ^ background label ^ norm (Euclidean)
378 \endcode
379 <b> Required Interface:</b>
380 \code
381
382 SrcImageIterator src_upperleft, src_lowerright;
383 DestImageIterator dest_upperleft;
384
385 SrcAccessor sa;
386 DestAccessor da;
387
388 ValueType background;
389 float distance;
390
391 sa(src_upperleft) != background;
392 da(dest_upperleft) < distance;
393 da.set(distance, dest_upperleft);
394
395 \endcode
396 \deprecatedEnd
397*/
398doxygen_overloaded_function(template <...> void distanceTransform)
399
400 template <class SrcImageIterator, class SrcAccessor,
401 class DestImageIterator, class DestAccessor,
402 class ValueType>
403 inline void
404 distanceTransform(SrcImageIterator src_upperleft,
405 SrcImageIterator src_lowerright, SrcAccessor sa,
406 DestImageIterator dest_upperleft, DestAccessor da,
407 ValueType background, int norm)
408{
409 if(norm == 1)
410 {
411 internalDistanceTransform(src_upperleft, src_lowerright, sa,
412 dest_upperleft, da, background,
413 InternalDistanceTransformL1NormFunctor());
414 }
415 else if(norm == 2)
416 {
417 internalDistanceTransform(src_upperleft, src_lowerright, sa,
418 dest_upperleft, da, background,
419 InternalDistanceTransformL2NormFunctor());
420 }
421 else
422 {
423 internalDistanceTransform(src_upperleft, src_lowerright, sa,
424 dest_upperleft, da, background,
425 InternalDistanceTransformLInifinityNormFunctor());
426 }
427}
428
429 template <class SrcImageIterator, class SrcAccessor,
430 class DestImageIterator, class DestAccessor,
431 class ValueType>
432 inline void
433 distanceTransform(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
434 pair<DestImageIterator, DestAccessor> dest,
435 ValueType background, int norm)
436{
437 distanceTransform(src.first, src.second, src.third,
438 dest.first, dest.second, background, norm);
439}
440
441 template <class T1, class S1,
442 class T2, class S2,
443 class ValueType>
444 inline void
445 distanceTransform(MultiArrayView<2, T1, S1> const & src,
446 MultiArrayView<2, T2, S2> dest,
447 ValueType background, int norm)
448{
449 vigra_precondition(src.shape() == dest.shape(),
450 "distanceTransform(): shape mismatch between input and output.");
451 distanceTransform(srcImageRange(src),
452 destImage(dest), background, norm);
453}
454
455//@}
456
457} // namespace vigra
458
459#endif // VIGRA_DISTANCETRANSFORM_HXX
460
vigra::BasicImage< float >::value_type
float value_type
Definition basicimage.hxx:481
vigra::BasicImage< float >::Iterator
BasicImageIterator< float, float ** > Iterator
Definition basicimage.hxx:532
vigra::Diff2D
Two dimensional difference vector.
Definition diff2d.hxx:186
vigra::MultiArrayView
Base class for, and view to, MultiArray.
Definition multi_array.hxx:705
vigra
vigra::FImage
BasicImage< float > FImage
Definition stdimage.hxx:143
vigra::norm
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition fftw3.hxx:1037
vigra::distanceTransform
void distanceTransform(...)

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python vigra 1.12.4