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

distancetransform.hxx VIGRA

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 */
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 
BasicImage< float > FImage
Definition: stdimage.hxx:143
Definition: accessor.hxx:43
doxygen_overloaded_function(template<... > void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition: fftw3.hxx:1037
BasicImageIterator< float, float ** > Iterator
Definition: basicimage.hxx:530
void distanceTransform(...)
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616
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.10.0