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