[ VIGRA Homepage | Class Index | Function Index | File Index | Main Page ]
![]() |
vigra/watersheds.hxx | ![]() |
---|
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 1998-2005 by Ullrich Koethe */ 00004 /* Cognitive Systems Group, University of Hamburg, Germany */ 00005 /* */ 00006 /* This file is part of the VIGRA computer vision library. */ 00007 /* ( Version 1.5.0, Dec 07 2006 ) */ 00008 /* The VIGRA Website is */ 00009 /* http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/ */ 00010 /* Please direct questions, bug reports, and contributions to */ 00011 /* koethe@informatik.uni-hamburg.de or */ 00012 /* vigra@kogs1.informatik.uni-hamburg.de */ 00013 /* */ 00014 /* Permission is hereby granted, free of charge, to any person */ 00015 /* obtaining a copy of this software and associated documentation */ 00016 /* files (the "Software"), to deal in the Software without */ 00017 /* restriction, including without limitation the rights to use, */ 00018 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00019 /* sell copies of the Software, and to permit persons to whom the */ 00020 /* Software is furnished to do so, subject to the following */ 00021 /* conditions: */ 00022 /* */ 00023 /* The above copyright notice and this permission notice shall be */ 00024 /* included in all copies or substantial portions of the */ 00025 /* Software. */ 00026 /* */ 00027 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00028 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00029 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00030 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00031 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00032 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00033 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00034 /* OTHER DEALINGS IN THE SOFTWARE. */ 00035 /* */ 00036 /************************************************************************/ 00037 00038 #ifndef VIGRA_WATERSHEDS_HXX 00039 #define VIGRA_WATERSHEDS_HXX 00040 00041 #include <functional> 00042 #include "mathutil.hxx" 00043 #include "stdimage.hxx" 00044 #include "pixelneighborhood.hxx" 00045 00046 namespace vigra { 00047 00048 template <class SrcIterator, class SrcAccessor, 00049 class DestIterator, class DestAccessor, 00050 class Neighborhood> 00051 unsigned int watershedLabeling(SrcIterator upperlefts, 00052 SrcIterator lowerrights, SrcAccessor sa, 00053 DestIterator upperleftd, DestAccessor da, 00054 Neighborhood neighborhood) 00055 { 00056 int w = lowerrights.x - upperlefts.x; 00057 int h = lowerrights.y - upperlefts.y; 00058 int x,y,i; 00059 00060 SrcIterator ys(upperlefts); 00061 SrcIterator xs(ys); 00062 00063 // temporary image to store region labels 00064 typedef IImage LabelImage; 00065 typedef LabelImage::traverser LabelTraverser; 00066 00067 LabelImage labelimage(w, h); 00068 00069 LabelTraverser yt = labelimage.upperLeft(); 00070 LabelTraverser xt(yt); 00071 00072 // Kovalevsky's clever idea to use 00073 // image iterator and scan order iterator simultaneously 00074 LabelImage::ScanOrderIterator label = labelimage.begin(); 00075 00076 // initialize the neighborhood circulators 00077 NeighborOffsetCirculator<Neighborhood> ncstart(Neighborhood::CausalFirst); 00078 NeighborOffsetCirculator<Neighborhood> ncstartBorder(Neighborhood::North); 00079 NeighborOffsetCirculator<Neighborhood> ncend(Neighborhood::CausalLast); 00080 ++ncend; 00081 NeighborOffsetCirculator<Neighborhood> ncendBorder(Neighborhood::North); 00082 ++ncendBorder; 00083 00084 // pass 1: scan image from upper left to lower right 00085 // to find connected components 00086 00087 // Each component will be represented by a tree of pixels. Each 00088 // pixel contains the scan order address of its parent in the 00089 // tree. In order for pass 2 to work correctly, the parent must 00090 // always have a smaller scan order address than the child. 00091 // Therefore, we can merge trees only at their roots, because the 00092 // root of the combined tree must have the smallest scan order 00093 // address among all the tree's pixels/ nodes. The root of each 00094 // tree is distinguished by pointing to itself (it contains its 00095 // own scan order address). This condition is enforced whenever a 00096 // new region is found or two regions are merged 00097 xs = ys; 00098 xt = yt; 00099 00100 *xt = 0; 00101 00102 ++xs.x; 00103 ++xt.x; 00104 for(x = 1; x != w; ++x, ++xs.x, ++xt.x) 00105 { 00106 if((*xs & Neighborhood::directionBit(Neighborhood::West)) || 00107 (xs[Neighborhood::west()] & Neighborhood::directionBit(Neighborhood::East))) 00108 { 00109 *xt = xt[Neighborhood::west()]; 00110 } 00111 else 00112 { 00113 *xt = x; 00114 } 00115 } 00116 00117 ++ys.y; 00118 ++yt.y; 00119 for(y = 1; y != h; ++y, ++ys.y, ++yt.y) 00120 { 00121 xs = ys; 00122 xt = yt; 00123 00124 for(x = 0; x != w; ++x, ++xs.x, ++xt.x) 00125 { 00126 NeighborOffsetCirculator<Neighborhood> nc(x == w-1 00127 ? ncstartBorder 00128 : ncstart); 00129 NeighborOffsetCirculator<Neighborhood> nce(x == 0 00130 ? ncendBorder 00131 : ncend); 00132 *xt = x + w*y; // default: new region 00133 for(; nc != nce; ++nc) 00134 { 00135 if((*xs & nc.directionBit()) || (xs[*nc] & nc.oppositeDirectionBit())) 00136 { 00137 int neighborLabel = xt[*nc]; 00138 // find the root label of a label tree 00139 while(neighborLabel != label[neighborLabel]) 00140 { 00141 neighborLabel = label[neighborLabel]; 00142 } 00143 if(neighborLabel < *xt) // always keep the smallest among the possible labels 00144 { 00145 label[*xt] = neighborLabel; 00146 *xt = neighborLabel; 00147 } 00148 else 00149 { 00150 label[neighborLabel] = *xt; 00151 } 00152 } 00153 } 00154 } 00155 } 00156 00157 // pass 2: assign one label to each region (tree) 00158 // so that labels form a consecutive sequence 1, 2, ... 00159 DestIterator yd(upperleftd); 00160 00161 unsigned int count = 0; 00162 i = 0; 00163 for(y=0; y != h; ++y, ++yd.y) 00164 { 00165 DestIterator xd(yd); 00166 for(x = 0; x != w; ++x, ++xd.x, ++i) 00167 { 00168 if(label[i] == i) 00169 { 00170 label[i] = ++count; 00171 } 00172 else 00173 { 00174 label[i] = label[label[i]]; 00175 } 00176 da.set(label[i], xd); 00177 } 00178 } 00179 return count; 00180 } 00181 00182 template <class SrcIterator, class SrcAccessor, 00183 class DestIterator, class DestAccessor> 00184 void prepareWatersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 00185 DestIterator upperleftd, DestAccessor da, 00186 FourNeighborCode) 00187 { 00188 int w = lowerrights.x - upperlefts.x; 00189 int h = lowerrights.y - upperlefts.y; 00190 int x,y; 00191 00192 SrcIterator ys(upperlefts); 00193 SrcIterator xs(ys); 00194 00195 DestIterator yd = upperleftd; 00196 00197 for(y = 0; y != h; ++y, ++ys.y, ++yd.y) 00198 { 00199 xs = ys; 00200 DestIterator xd = yd; 00201 00202 for(x = 0; x != w; ++x, ++xs.x, ++xd.x) 00203 { 00204 AtImageBorder atBorder = isAtImageBorder(x,y,w,h); 00205 typename SrcAccessor::value_type v = sa(xs); 00206 // the following choice causes minima to point 00207 // to their lowest neighbor -- would this be better??? 00208 // typename SrcAccessor::value_type v = NumericTraits<typename SrcAccessor::value_type>::max(); 00209 int o = 0; // means center is minimum 00210 if(atBorder == NotAtBorder) 00211 { 00212 NeighborhoodCirculator<SrcIterator, FourNeighborCode> c(xs), cend(c); 00213 do { 00214 if(sa(c) <= v) 00215 { 00216 v = sa(c); 00217 o = c.directionBit(); 00218 } 00219 } 00220 while(++c != cend); 00221 } 00222 else 00223 { 00224 RestrictedNeighborhoodCirculator<SrcIterator, FourNeighborCode> c(xs, atBorder), cend(c); 00225 do { 00226 if(sa(c) <= v) 00227 { 00228 v = sa(c); 00229 o = c.directionBit(); 00230 } 00231 } 00232 while(++c != cend); 00233 } 00234 da.set(o, xd); 00235 } 00236 } 00237 } 00238 00239 template <class SrcIterator, class SrcAccessor, 00240 class DestIterator, class DestAccessor> 00241 void prepareWatersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 00242 DestIterator upperleftd, DestAccessor da, 00243 EightNeighborCode) 00244 { 00245 int w = lowerrights.x - upperlefts.x; 00246 int h = lowerrights.y - upperlefts.y; 00247 int x,y; 00248 00249 SrcIterator ys(upperlefts); 00250 SrcIterator xs(ys); 00251 00252 DestIterator yd = upperleftd; 00253 00254 for(y = 0; y != h; ++y, ++ys.y, ++yd.y) 00255 { 00256 xs = ys; 00257 DestIterator xd = yd; 00258 00259 for(x = 0; x != w; ++x, ++xs.x, ++xd.x) 00260 { 00261 AtImageBorder atBorder = isAtImageBorder(x,y,w,h); 00262 typename SrcAccessor::value_type v = sa(xs); 00263 // the following choice causes minima to point 00264 // to their lowest neighbor -- would this be better??? 00265 // typename SrcAccessor::value_type v = NumericTraits<typename SrcAccessor::value_type>::max(); 00266 int o = 0; // means center is minimum 00267 if(atBorder == NotAtBorder) 00268 { 00269 // handle diagonal and principal neighbors separately 00270 // so that principal neighbors are preferred when there are 00271 // candidates with equal strength 00272 NeighborhoodCirculator<SrcIterator, EightNeighborCode> 00273 c(xs, EightNeighborCode::NorthEast); 00274 for(int i = 0; i < 4; ++i, c += 2) 00275 { 00276 if(sa(c) <= v) 00277 { 00278 v = sa(c); 00279 o = c.directionBit(); 00280 } 00281 } 00282 --c; 00283 for(int i = 0; i < 4; ++i, c += 2) 00284 { 00285 if(sa(c) <= v) 00286 { 00287 v = sa(c); 00288 o = c.directionBit(); 00289 } 00290 } 00291 } 00292 else 00293 { 00294 RestrictedNeighborhoodCirculator<SrcIterator, EightNeighborCode> 00295 c(xs, atBorder), cend(c); 00296 do 00297 { 00298 if(!c.isDiagonal()) 00299 continue; 00300 if(sa(c) <= v) 00301 { 00302 v = sa(c); 00303 o = c.directionBit(); 00304 } 00305 } 00306 while(++c != cend); 00307 do 00308 { 00309 if(c.isDiagonal()) 00310 continue; 00311 if(sa(c) <= v) 00312 { 00313 v = sa(c); 00314 o = c.directionBit(); 00315 } 00316 } 00317 while(++c != cend); 00318 } 00319 da.set(o, xd); 00320 } 00321 } 00322 } 00323 00324 /** \addtogroup SeededRegionGrowing Region Segmentation Algorithms 00325 Region growing, watersheds, and voronoi tesselation 00326 */ 00327 //@{ 00328 00329 /********************************************************/ 00330 /* */ 00331 /* watersheds */ 00332 /* */ 00333 /********************************************************/ 00334 00335 /** \brief Region Segmentation by means of the watershed algorithm. 00336 00337 This function implements the union-find version of the watershed algorithms 00338 as described in 00339 00340 J. Roerdink, R. Meijster: "<em>The watershed transform: definitions, algorithms, 00341 and parallelization stretegies</em>", Fundamenta Informaticae, 41:187-228, 2000 00342 00343 The source image is a boundary indicator such as the gradient magnitude 00344 of the trace of the \ref boundaryTensor(). Local minima of the boundary indicator 00345 are used as region seeds, and all other pixels are recursively assigned to the same 00346 region as their lowest neighbor. Pass \ref vigra::EightNeighborCode or 00347 \ref vigra::FourNeighborCode to determine the neighborhood where pixel values 00348 are compared. The pixel type of the input image must be <tt>LessThanComparable</tt>. 00349 The function uses accessors. 00350 00351 Note that VIGRA provides an alternative implementaion of the watershed transform via 00352 \ref seededRegionGrowing(). It is slower, but handles plateaus better 00353 and allows to keep a one pixel wide boundary between regions. 00354 00355 <b> Declarations:</b> 00356 00357 pass arguments explicitly: 00358 \code 00359 namespace vigra { 00360 template <class SrcIterator, class SrcAccessor, 00361 class DestIterator, class DestAccessor, 00362 class Neighborhood = EightNeighborCode> 00363 unsigned int 00364 watersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 00365 DestIterator upperleftd, DestAccessor da, 00366 Neighborhood neighborhood = EightNeighborCode()) 00367 } 00368 \endcode 00369 00370 use argument objects in conjunction with \ref ArgumentObjectFactories: 00371 \code 00372 namespace vigra { 00373 template <class SrcIterator, class SrcAccessor, 00374 class DestIterator, class DestAccessor, 00375 class Neighborhood = EightNeighborCode> 00376 unsigned int 00377 watersheds(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00378 pair<DestIterator, DestAccessor> dest, 00379 Neighborhood neighborhood = EightNeighborCode()) 00380 } 00381 \endcode 00382 00383 <b> Usage:</b> 00384 00385 <b>\#include</b> "<a href="watersheds_8hxx-source.html">vigra/watersheds.hxx</a>"<br> 00386 Namespace: vigra 00387 00388 Example: watersheds of the gradient magnitude. 00389 00390 \code 00391 vigra::BImage in(w,h); 00392 ... // read input data 00393 00394 vigra::FImage gradx(x,y), grady(x,y), gradMag(x,y); 00395 gaussianGradient(srcImageRange(src), destImage(gradx), destImage(grady), 3.0); 00396 combineTwoImages(srcImageRange(gradx), srcImage(grady), destImage(gradMag), 00397 vigra::MagnitudeFunctor<float>()); 00398 00399 // the pixel type of the destination image must be large enough to hold 00400 // numbers up to 'max_region_label' to prevent overflow 00401 vigra::IImage labeling(x,y); 00402 int max_region_label = watersheds(srcImageRange(gradMag), destImage(labeling)); 00403 00404 \endcode 00405 00406 <b> Required Interface:</b> 00407 00408 \code 00409 SrcImageIterator src_upperleft, src_lowerright; 00410 DestImageIterator dest_upperleft; 00411 00412 SrcAccessor src_accessor; 00413 DestAccessor dest_accessor; 00414 00415 // compare src values 00416 src_accessor(src_upperleft) <= src_accessor(src_upperleft) 00417 00418 // set result 00419 int label; 00420 dest_accessor.set(label, dest_upperleft); 00421 \endcode 00422 */ 00423 template <class SrcIterator, class SrcAccessor, 00424 class DestIterator, class DestAccessor, 00425 class Neighborhood> 00426 unsigned int 00427 watersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 00428 DestIterator upperleftd, DestAccessor da, Neighborhood neighborhood) 00429 { 00430 SImage orientationImage(lowerrights - upperlefts); 00431 SImage::traverser yo = orientationImage.upperLeft(); 00432 00433 prepareWatersheds(upperlefts, lowerrights, sa, 00434 orientationImage.upperLeft(), orientationImage.accessor(), neighborhood); 00435 return watershedLabeling(orientationImage.upperLeft(), orientationImage.lowerRight(), orientationImage.accessor(), 00436 upperleftd, da, neighborhood); 00437 } 00438 00439 template <class SrcIterator, class SrcAccessor, 00440 class DestIterator, class DestAccessor> 00441 inline unsigned int 00442 watersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 00443 DestIterator upperleftd, DestAccessor da) 00444 { 00445 return watersheds(upperlefts, lowerrights, sa, upperleftd, da, EightNeighborCode()); 00446 } 00447 00448 template <class SrcIterator, class SrcAccessor, 00449 class DestIterator, class DestAccessor, 00450 class Neighborhood> 00451 inline unsigned int 00452 watersheds(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00453 pair<DestIterator, DestAccessor> dest, Neighborhood neighborhood) 00454 { 00455 return watersheds(src.first, src.second, src.third, dest.first, dest.second, neighborhood); 00456 } 00457 00458 template <class SrcIterator, class SrcAccessor, 00459 class DestIterator, class DestAccessor> 00460 inline unsigned int 00461 watersheds(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00462 pair<DestIterator, DestAccessor> dest) 00463 { 00464 return watersheds(src.first, src.second, src.third, dest.first, dest.second); 00465 } 00466 00467 //@} 00468 00469 } // namespace vigra 00470 00471 #endif // VIGRA_WATERSHEDS_HXX
© Ullrich Köthe (koethe@informatik.uni-hamburg.de) |
html generated using doxygen and Python
|