[ 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
84 internalDistanceTransform(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 Perform a distance transform using either the Euclidean, Manhattan,
278 or chessboard metrics.
279
280 See also: \ref MultiArrayDistanceTransform "multi-dimensional distance transforms"
281 */
282 //@{
283
284 /** For all background pixels, calculate the distance to
285 the nearest object or contour. The label of the pixels to be considered
286 background in the source image is passed in the parameter 'background'.
287 Source pixels with other labels will be considered objects. In the
288 destination image, all pixels corresponding to background will be assigned
289 the their distance value, all pixels corresponding to objects will be
290 assigned 0.
291
292 The parameter 'norm' gives the distance norm to be used:
293
294 <ul>
295
296 <li> norm == 0: use chessboard distance (L-infinity norm)
297 <li> norm == 1: use Manhattan distance (L1 norm)
298 <li> norm == 2: use Euclidean distance (L2 norm)
299
300 </ul>
301
302 If you use the L2 norm, the destination pixels must be real valued to give
303 correct results.
304
305 <b> Declarations:</b>
306
307 pass 2D array views:
308 \code
309 namespace vigra {
310 template <class T1, class S1,
311 class T2, class S2,
312 class ValueType>
313 void
314 distanceTransform(MultiArrayView<2, T1, S1> const & src,
315 MultiArrayView<2, T2, S2> dest,
316 ValueType background, int norm);
317 }
318 \endcode
319
320 \deprecatedAPI{distanceTransform}
321 pass \ref ImageIterators and \ref DataAccessors :
322 \code
323 namespace vigra {
324 template <class SrcImageIterator, class SrcAccessor,
325 class DestImageIterator, class DestAccessor,
326 class ValueType>
327 void distanceTransform(SrcImageIterator src_upperleft,
328 SrcImageIterator src_lowerright, SrcAccessor sa,
329 DestImageIterator dest_upperleft, DestAccessor da,
330 ValueType background, int norm);
331 }
332 \endcode
333 use argument objects in conjunction with \ref ArgumentObjectFactories :
334 \code
335 namespace vigra {
336 template <class SrcImageIterator, class SrcAccessor,
337 class DestImageIterator, class DestAccessor,
338 class ValueType>
339 void distanceTransform(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
340 pair<DestImageIterator, DestAccessor> dest,
341 ValueType background, int norm);
342 }
343 \endcode
344 \deprecatedEnd
345
346 <b> Usage:</b>
347
348 <b>\#include</b> <vigra/distancetransform.hxx><br>
349 Namespace: vigra
350
351 \code
352 MultiArray<2, unsigned char> src(w,h), edges(w,h);
353 MultiArray<2, float> distance(w, h);
354 ...
355
356 // detect edges in src image (edges will be marked 1, background 0)
357 differenceOfExponentialEdgeImage(src, edges, 0.8, 4.0);
358
359 // find distance of all pixels from nearest edge
360 distanceTransform(edges, distance, 0, 2);
361 // ^ background label ^ norm (Euclidean)
362 \endcode
363
364 \deprecatedUsage{distanceTransform}
365 \code
366
367 vigra::BImage src(w,h), edges(w,h);
368 vigra::FImage distance(w, h);
369
370 // empty edge image
371 edges = 0;
372 ...
373
374 // detect edges in src image (edges will be marked 1, background 0)
375 vigra::differenceOfExponentialEdgeImage(srcImageRange(src), destImage(edges),
376 0.8, 4.0);
377
378 // find distance of all pixels from nearest edge
379 vigra::distanceTransform(srcImageRange(edges), destImage(distance),
380 0, 2);
381 // ^ background label ^ norm (Euclidean)
382 \endcode
383 <b> Required Interface:</b>
384 \code
385
386 SrcImageIterator src_upperleft, src_lowerright;
387 DestImageIterator dest_upperleft;
388
389 SrcAccessor sa;
390 DestAccessor da;
391
392 ValueType background;
393 float distance;
394
395 sa(src_upperleft) != background;
396 da(dest_upperleft) < distance;
397 da.set(distance, dest_upperleft);
398
399 \endcode
400 \deprecatedEnd
401 */
402 doxygen_overloaded_function(template <...> void distanceTransform)
403
404 template <class SrcImageIterator, class SrcAccessor,
405 class DestImageIterator, class DestAccessor,
406 class ValueType>
407 inline void
408 distanceTransform(SrcImageIterator src_upperleft,
409 SrcImageIterator src_lowerright, SrcAccessor sa,
410 DestImageIterator dest_upperleft, DestAccessor da,
411 ValueType background, int norm)
412 {
413 if(norm == 1)
414 {
415 internalDistanceTransform(src_upperleft, src_lowerright, sa,
416 dest_upperleft, da, background,
417 InternalDistanceTransformL1NormFunctor());
418 }
419 else if(norm == 2)
420 {
421 internalDistanceTransform(src_upperleft, src_lowerright, sa,
422 dest_upperleft, da, background,
423 InternalDistanceTransformL2NormFunctor());
424 }
425 else
426 {
427 internalDistanceTransform(src_upperleft, src_lowerright, sa,
428 dest_upperleft, da, background,
429 InternalDistanceTransformLInifinityNormFunctor());
430 }
431 }
432
433 template <class SrcImageIterator, class SrcAccessor,
434 class DestImageIterator, class DestAccessor,
435 class ValueType>
436 inline void
437 distanceTransform(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
438 pair<DestImageIterator, DestAccessor> dest,
439 ValueType background, int norm)
440 {
441 distanceTransform(src.first, src.second, src.third,
442 dest.first, dest.second, background, norm);
443 }
444
445 template <class T1, class S1,
446 class T2, class S2,
447 class ValueType>
448 inline void
449 distanceTransform(MultiArrayView<2, T1, S1> const & src,
450 MultiArrayView<2, T2, S2> dest,
451 ValueType background, int norm)
452 {
453 vigra_precondition(src.shape() == dest.shape(),
454 "distanceTransform(): shape mismatch between input and output.");
455 distanceTransform(srcImageRange(src),
456 destImage(dest), background, norm);
457 }
458
459 //@}
460
461 } // namespace vigra
462
463 #endif // VIGRA_DISTANCETRANSFORM_HXX
464
vigra::FImage
BasicImage< float > FImage
Definition: stdimage.hxx:143
vigra
Definition: accessor.hxx:43
vigra::doxygen_overloaded_function
doxygen_overloaded_function(template<... > void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
vigra::norm
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition: fftw3.hxx:1037
vigra::BasicImage< float >::Iterator
BasicImageIterator< float, float **> Iterator
Definition: basicimage.hxx:530
vigra::distanceTransform
void distanceTransform(...)
vigra::sqrt
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616
vigra::BasicImage< float >::value_type
float value_type
Definition: basicimage.hxx:479

© 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.11.0