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