[ 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 
43 namespace vigra {
44 
45 /*
46  * functors to determine the distance norm
47  * these functors assume that dx and dy are positive
48  * (this is OK for use in internalDistanceTransform())
49  */
50 
51 // chessboard metric
52 struct InternalDistanceTransformLInifinityNormFunctor
53 {
54  float operator()(float dx, float dy) const
55  {
56  return (dx < dy) ? dy : dx;
57  }
58 };
59 
60 // Manhattan metric
61 struct InternalDistanceTransformL1NormFunctor
62 {
63  float operator()(float dx, float dy) const
64  {
65  return dx + dy;
66  }
67 };
68 
69 // Euclidean metric
70 struct InternalDistanceTransformL2NormFunctor
71 {
72  float operator()(float dx, float dy) const
73  {
74  return VIGRA_CSTD::sqrt(dx*dx + dy*dy);
75  }
76 };
77 
78 
79 template <class SrcImageIterator, class SrcAccessor,
80  class DestImageIterator, class DestAccessor,
81  class ValueType, class NormFunctor>
82 void
83 internalDistanceTransform(SrcImageIterator src_upperleft,
84  SrcImageIterator src_lowerright, SrcAccessor sa,
85  DestImageIterator dest_upperleft, DestAccessor da,
86  ValueType background, NormFunctor norm)
87 {
88  int w = src_lowerright.x - src_upperleft.x;
89  int h = src_lowerright.y - src_upperleft.y;
90 
91  FImage xdist(w,h), ydist(w,h);
92 
93  xdist = (FImage::value_type)w; // init x and
94  ydist = (FImage::value_type)h; // y distances with 'large' values
95 
96  SrcImageIterator sy = src_upperleft;
97  DestImageIterator ry = dest_upperleft;
98  FImage::Iterator xdy = xdist.upperLeft();
99  FImage::Iterator ydy = ydist.upperLeft();
100  SrcImageIterator sx = sy;
101  DestImageIterator rx = ry;
102  FImage::Iterator xdx = xdy;
103  FImage::Iterator ydx = ydy;
104 
105  static const Diff2D left(-1, 0);
106  static const Diff2D right(1, 0);
107  static const Diff2D top(0, -1);
108  static const Diff2D bottom(0, 1);
109 
110  int x,y;
111  if(sa(sx) != background) // first pixel
112  {
113  *xdx = 0.0;
114  *ydx = 0.0;
115  da.set(0.0, rx);
116  }
117  else
118  {
119  da.set(norm(*xdx, *ydx), rx);
120  }
121 
122 
123  for(x=1, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x;
124  x<w;
125  ++x, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x) // first row left to right
126  {
127  if(sa(sx) != background)
128  {
129  *xdx = 0.0;
130  *ydx = 0.0;
131  da.set(0.0, rx);
132  }
133  else
134  {
135  *xdx = xdx[left] + 1.0f; // propagate x and
136  *ydx = ydx[left]; // y components of distance from left pixel
137  da.set(norm(*xdx, *ydx), rx); // calculate distance from x and y components
138  }
139  }
140  for(x=w-2, xdx.x -= 2, ydx.x -= 2, sx.x -= 2, rx.x -= 2;
141  x>=0;
142  --x, --xdx.x, --ydx.x, --sx.x, --rx.x) // first row right to left
143  {
144  float d = norm(xdx[right] + 1.0f, ydx[right]);
145 
146  if(da(rx) < d) continue;
147 
148  *xdx = xdx[right] + 1.0f;
149  *ydx = ydx[right];
150  da.set(d, rx);
151  }
152  for(y=1, ++xdy.y, ++ydy.y, ++sy.y, ++ry.y;
153  y<h;
154  ++y, ++xdy.y, ++ydy.y, ++sy.y, ++ry.y) // top to bottom
155  {
156  sx = sy;
157  rx = ry;
158  xdx = xdy;
159  ydx = ydy;
160 
161  if(sa(sx) != background) // first pixel of current row
162  {
163  *xdx = 0.0f;
164  *ydx = 0.0f;
165  da.set(0.0, rx);
166  }
167  else
168  {
169  *xdx = xdx[top];
170  *ydx = ydx[top] + 1.0f;
171  da.set(norm(*xdx, *ydx), rx);
172  }
173 
174  for(x=1, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x;
175  x<w;
176  ++x, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x) // current row left to right
177  {
178  if(sa(sx) != background)
179  {
180  *xdx = 0.0f;
181  *ydx = 0.0f;
182  da.set(0.0, rx);
183  }
184  else
185  {
186  float d1 = norm(xdx[left] + 1.0f, ydx[left]);
187  float d2 = norm(xdx[top], ydx[top] + 1.0f);
188 
189  if(d1 < d2)
190  {
191  *xdx = xdx[left] + 1.0f;
192  *ydx = ydx[left];
193  da.set(d1, rx);
194  }
195  else
196  {
197  *xdx = xdx[top];
198  *ydx = ydx[top] + 1.0f;
199  da.set(d2, rx);
200  }
201  }
202  }
203  for(x=w-2, xdx.x -= 2, ydx.x -= 2, sx.x -= 2, rx.x -= 2;
204  x>=0;
205  --x, --xdx.x, --ydx.x, --sx.x, --rx.x) // current row right to left
206  {
207  float d1 = norm(xdx[right] + 1.0f, ydx[right]);
208 
209  if(da(rx) < d1) continue;
210 
211  *xdx = xdx[right] + 1.0f;
212  *ydx = ydx[right];
213  da.set(d1, rx);
214  }
215  }
216  for(y=h-2, xdy.y -= 2, ydy.y -= 2, sy.y -= 2, ry.y -= 2;
217  y>=0;
218  --y, --xdy.y, --ydy.y, --sy.y, --ry.y) // bottom to top
219  {
220  sx = sy;
221  rx = ry;
222  xdx = xdy;
223  ydx = ydy;
224 
225  float d = norm(xdx[bottom], ydx[bottom] + 1.0f);
226  if(d < da(rx)) // first pixel of current row
227  {
228  *xdx = xdx[bottom];
229  *ydx = ydx[bottom] + 1.0f;
230  da.set(d, rx);
231  }
232 
233  for(x=1, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x;
234  x<w;
235  ++x, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x) // current row left to right
236  {
237  float d1 = norm(xdx[left] + 1.0f, ydx[left]);
238  float d2 = norm(xdx[bottom], ydx[bottom] + 1.0f);
239 
240  if(d1 < d2)
241  {
242  if(da(rx) < d1) continue;
243  *xdx = xdx[left] + 1.0f;
244  *ydx = ydx[left];
245  da.set(d1, rx);
246  }
247  else
248  {
249  if(da(rx) < d2) continue;
250  *xdx = xdx[bottom];
251  *ydx = ydx[bottom] + 1.0f;
252  da.set(d2, rx);
253  }
254  }
255  for(x=w-2, xdx.x -= 2, ydx.x -= 2, sx.x -= 2, rx.x -= 2;
256  x>=0;
257  --x, --xdx.x, --ydx.x, --sx.x, --rx.x) // current row right to left
258  {
259  float d1 = norm(xdx[right] + 1.0f, ydx[right]);
260 
261  if(da(rx) < d1) continue;
262  *xdx = xdx[right] + 1.0f;
263  *ydx = ydx[right];
264  da.set(d1, rx);
265  }
266  }
267 }
268 
269 /********************************************************/
270 /* */
271 /* distanceTransform */
272 /* */
273 /********************************************************/
274 
275 /** \addtogroup DistanceTransform Distance Transform
276  Perform a distance transform using either the Euclidean, Manhattan,
277  or chessboard metrics.
278 
279  See also: \ref MultiArrayDistanceTransform "multi-dimensional distance transforms"
280 */
281 //@{
282 
283 /** For all background pixels, calculate the distance to
284  the nearest object or contour. The label of the pixels to be considered
285  background in the source image is passed in the parameter 'background'.
286  Source pixels with other labels will be considered objects. In the
287  destination image, all pixels corresponding to background will be assigned
288  the their distance value, all pixels corresponding to objects will be
289  assigned 0.
290 
291  The parameter 'norm' gives the distance norm to be used:
292 
293  <ul>
294 
295  <li> norm == 0: use chessboard distance (L-infinity norm)
296  <li> norm == 1: use Manhattan distance (L1 norm)
297  <li> norm == 2: use Euclidean distance (L2 norm)
298 
299  </ul>
300 
301  If you use the L2 norm, the destination pixels must be real valued to give
302  correct results.
303 
304  <b> Declarations:</b>
305 
306  pass arguments explicitly:
307  \code
308  namespace vigra {
309  template <class SrcImageIterator, class SrcAccessor,
310  class DestImageIterator, class DestAccessor,
311  class ValueType>
312  void distanceTransform(SrcImageIterator src_upperleft,
313  SrcImageIterator src_lowerright, SrcAccessor sa,
314  DestImageIterator dest_upperleft, DestAccessor da,
315  ValueType background, int norm)
316  }
317  \endcode
318 
319 
320  use argument objects in conjunction with \ref ArgumentObjectFactories :
321  \code
322  namespace vigra {
323  template <class SrcImageIterator, class SrcAccessor,
324  class DestImageIterator, class DestAccessor,
325  class ValueType>
326  void distanceTransform(
327  triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
328  pair<DestImageIterator, DestAccessor> dest,
329  ValueType background, int norm)
330  }
331  \endcode
332 
333  <b> Usage:</b>
334 
335  <b>\#include</b> <vigra/distancetransform.hxx><br>
336  Namespace: vigra
337 
338 
339  \code
340 
341  vigra::BImage src(w,h), edges(w,h);
342  vigra::FImage distance(w, h);
343 
344  // empty edge image
345  edges = 0;
346  ...
347 
348  // detect edges in src image (edges will be marked 1, background 0)
349  vigra::differenceOfExponentialEdgeImage(srcImageRange(src), destImage(edges),
350  0.8, 4.0);
351 
352  // find distance of all pixels from nearest edge
353  vigra::distanceTransform(srcImageRange(edges), destImage(distance),
354  0, 2);
355  // ^ background label ^ norm (Euclidean)
356  \endcode
357 
358  <b> Required Interface:</b>
359 
360  \code
361 
362  SrcImageIterator src_upperleft, src_lowerright;
363  DestImageIterator dest_upperleft;
364 
365  SrcAccessor sa;
366  DestAccessor da;
367 
368  ValueType background;
369  float distance;
370 
371  sa(src_upperleft) != background;
372  da(dest_upperleft) < distance;
373  da.set(distance, dest_upperleft);
374 
375  \endcode
376 */
377 doxygen_overloaded_function(template <...> void distanceTransform)
378 
379 template <class SrcImageIterator, class SrcAccessor,
380  class DestImageIterator, class DestAccessor,
381  class ValueType>
382 inline void
383 distanceTransform(SrcImageIterator src_upperleft,
384  SrcImageIterator src_lowerright, SrcAccessor sa,
385  DestImageIterator dest_upperleft, DestAccessor da,
386  ValueType background, int norm)
387 {
388  if(norm == 1)
389  {
390  internalDistanceTransform(src_upperleft, src_lowerright, sa,
391  dest_upperleft, da, background,
392  InternalDistanceTransformL1NormFunctor());
393  }
394  else if(norm == 2)
395  {
396  internalDistanceTransform(src_upperleft, src_lowerright, sa,
397  dest_upperleft, da, background,
398  InternalDistanceTransformL2NormFunctor());
399  }
400  else
401  {
402  internalDistanceTransform(src_upperleft, src_lowerright, sa,
403  dest_upperleft, da, background,
404  InternalDistanceTransformLInifinityNormFunctor());
405  }
406 }
407 
408 template <class SrcImageIterator, class SrcAccessor,
409  class DestImageIterator, class DestAccessor,
410  class ValueType>
411 inline void
413  triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
414  pair<DestImageIterator, DestAccessor> dest,
415  ValueType background, int norm)
416 {
417  distanceTransform(src.first, src.second, src.third,
418  dest.first, dest.second, background, norm);
419 }
420 
421 //@}
422 
423 } // namespace vigra
424 
425 #endif // VIGRA_DISTANCETRANSFORM_HXX
426 

© 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.9.0 (Sat Oct 5 2013)