[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
vigra/watersheds3d.hxx | ![]() |
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 2006-2007 by F. Heinrich, B. Seppke, 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_watersheds3D_HXX 00037 #define VIGRA_watersheds3D_HXX 00038 00039 #include "voxelneighborhood.hxx" 00040 #include "multi_array.hxx" 00041 #include "multi_localminmax.hxx" 00042 #include "labelvolume.hxx" 00043 #include "seededregiongrowing3d.hxx" 00044 #include "watersheds.hxx" 00045 00046 namespace vigra 00047 { 00048 00049 template <class SrcIterator, class SrcAccessor, class SrcShape, 00050 class DestIterator, class DestAccessor, class Neighborhood3D> 00051 int preparewatersheds3D( SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa, 00052 DestIterator d_Iter, DestAccessor da, Neighborhood3D) 00053 { 00054 //basically needed for iteration and border-checks 00055 int w = srcShape[0], h = srcShape[1], d = srcShape[2]; 00056 int x,y,z, local_min_count=0; 00057 00058 //declare and define Iterators for all three dims at src 00059 SrcIterator zs = s_Iter; 00060 SrcIterator ys(zs); 00061 SrcIterator xs(ys); 00062 00063 //Declare Iterators for all three dims at dest 00064 DestIterator zd = d_Iter; 00065 00066 for(z = 0; z != d; ++z, ++zs.dim2(), ++zd.dim2()) 00067 { 00068 ys = zs; 00069 DestIterator yd(zd); 00070 00071 for(y = 0; y != h; ++y, ++ys.dim1(), ++yd.dim1()) 00072 { 00073 xs = ys; 00074 DestIterator xd(yd); 00075 00076 for(x = 0; x != w; ++x, ++xs.dim0(), ++xd.dim0()) 00077 { 00078 AtVolumeBorder atBorder = isAtVolumeBorder(x,y,z,w,h,d); 00079 typename SrcAccessor::value_type v = sa(xs); 00080 // the following choice causes minima to point 00081 // to their lowest neighbor -- would this be better??? 00082 // typename SrcAccessor::value_type v = NumericTraits<typename SrcAccessor::value_type>::max(); 00083 int o = 0; // means center is minimum 00084 typename SrcAccessor::value_type my_v = v; 00085 if(atBorder == NotAtBorder) 00086 { 00087 NeighborhoodCirculator<SrcIterator, Neighborhood3D> c(xs), cend(c); 00088 00089 do { 00090 if(sa(c) < v) 00091 { 00092 v = sa(c); 00093 o = c.directionBit(); 00094 } 00095 else if(sa(c) == my_v && my_v == v) 00096 { 00097 o = o | c.directionBit(); 00098 } 00099 } 00100 while(++c != cend); 00101 } 00102 else 00103 { 00104 RestrictedNeighborhoodCirculator<SrcIterator, Neighborhood3D> c(xs, atBorder), cend(c); 00105 do { 00106 if(sa(c) < v) 00107 { 00108 v = sa(c); 00109 o = c.directionBit(); 00110 } 00111 else if(sa(c) == my_v && my_v == v) 00112 { 00113 o = o | c.directionBit(); 00114 } 00115 } 00116 while(++c != cend); 00117 } 00118 if (o==0) local_min_count++; 00119 da.set(o, xd); 00120 }//end x-iteration 00121 }//end y-iteration 00122 }//end z-iteration 00123 return local_min_count; 00124 } 00125 00126 template <class SrcIterator, class SrcAccessor,class SrcShape, 00127 class DestIterator, class DestAccessor, 00128 class Neighborhood3D> 00129 unsigned int watershedLabeling3D( SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa, 00130 DestIterator d_Iter, DestAccessor da, 00131 Neighborhood3D) 00132 { 00133 typedef typename DestAccessor::value_type LabelType; 00134 00135 //basically needed for iteration and border-checks 00136 int w = srcShape[0], h = srcShape[1], d = srcShape[2]; 00137 int x,y,z; 00138 00139 //declare and define Iterators for all three dims at src 00140 SrcIterator zs = s_Iter; 00141 DestIterator zd = d_Iter; 00142 00143 // temporary image to store region labels 00144 detail::UnionFindArray<LabelType> labels; 00145 00146 // initialize the neighborhood traversers 00147 NeighborOffsetCirculator<Neighborhood3D> nc(Neighborhood3D::CausalFirst); 00148 NeighborOffsetCirculator<Neighborhood3D> nce(Neighborhood3D::CausalLast); 00149 ++nce; 00150 // pass 1: scan image from upper left front to lower right back 00151 // to find connected components 00152 00153 // Each component will be represented by a tree of pixels. Each 00154 // pixel contains the scan order address of its parent in the 00155 // tree. In order for pass 2 to work correctly, the parent must 00156 // always have a smaller scan order address than the child. 00157 // Therefore, we can merge trees only at their roots, because the 00158 // root of the combined tree must have the smallest scan order 00159 // address among all the tree's pixels/ nodes. The root of each 00160 // tree is distinguished by pointing to itself (it contains its 00161 // own scan order address). This condition is enforced whenever a 00162 // new region is found or two regions are merged 00163 for(z = 0; z != d; ++z, ++zs.dim2(), ++zd.dim2()) 00164 { 00165 SrcIterator ys = zs; 00166 DestIterator yd = zd; 00167 00168 for(y = 0; y != h; ++y, ++ys.dim1(), ++yd.dim1()) 00169 { 00170 SrcIterator xs = ys; 00171 DestIterator xd = yd; 00172 00173 for(x = 0; x != w; ++x, ++xs.dim0(), ++xd.dim0()) 00174 { 00175 LabelType currentLabel = labels.nextFreeLabel(); // default: new region 00176 00177 //check whether there is a special border treatment to be used or not 00178 AtVolumeBorder atBorder = isAtVolumeBorderCausal(x,y,z,w,h,d); 00179 00180 //We are not at the border! 00181 if(atBorder == NotAtBorder) 00182 { 00183 00184 nc = NeighborOffsetCirculator<Neighborhood3D>(Neighborhood3D::CausalFirst); 00185 00186 do 00187 { 00188 // Direction of NTraversr Neighbor's direction bit is pointing 00189 // = Direction of voxel towards us? 00190 if((sa(xs) & nc.directionBit()) || (sa(xs,*nc) & nc.oppositeDirectionBit())) 00191 { 00192 currentLabel = labels.makeUnion(da(xd,*nc), currentLabel); 00193 } 00194 ++nc; 00195 }while(nc!=nce); 00196 } 00197 //we are at a border - handle this!! 00198 else 00199 { 00200 nc = NeighborOffsetCirculator<Neighborhood3D>(Neighborhood3D::nearBorderDirectionsCausal(atBorder,0)); 00201 int j=0; 00202 while(nc.direction() != Neighborhood3D::Error) 00203 { 00204 // Direction of NTraversr Neighbor's direction bit is pointing 00205 // = Direction of voxel towards us? 00206 if((sa(xs) & nc.directionBit()) || (sa(xs,*nc) & nc.oppositeDirectionBit())) 00207 { 00208 currentLabel = labels.makeUnion(da(xd,*nc), currentLabel); 00209 } 00210 nc.turnTo(Neighborhood3D::nearBorderDirectionsCausal(atBorder,++j)); 00211 } 00212 } 00213 da.set(labels.finalizeLabel(currentLabel), xd); 00214 } 00215 } 00216 } 00217 00218 unsigned int count = labels.makeContiguous(); 00219 00220 // pass 2: assign one label to each region (tree) 00221 // so that labels form a consecutive sequence 1, 2, ... 00222 zd = d_Iter; 00223 for(z=0; z != d; ++z, ++zd.dim2()) 00224 { 00225 DestIterator yd(zd); 00226 00227 for(y=0; y != h; ++y, ++yd.dim1()) 00228 { 00229 DestIterator xd(yd); 00230 00231 for(x = 0; x != w; ++x, ++xd.dim0()) 00232 { 00233 da.set(labels[da(xd)], xd); 00234 } 00235 } 00236 } 00237 return count; 00238 } 00239 00240 00241 /** \addtogroup SeededRegionGrowing 00242 */ 00243 //@{ 00244 00245 /** \brief Generate seeds for watershed computation and seeded region growing. 00246 00247 The source image is a boundary indicator such as the gradient magnitude 00248 or the trace of the \ref boundaryTensor(). Seeds are generally generated 00249 at locations where the boundaryness (i.e. the likelihood of the point being on the 00250 boundary) is very small. In particular, seeds can be placed by either 00251 looking for local minima (possibly including minimal plateaus) of the boundaryness, 00252 of by looking at level sets (i.e. regions where the boundaryness is below a threshold). 00253 Both methods can also be combined, so that only minima below a threshold are returned. 00254 The particular seeding strategy is specified by the <tt>options</tt> object 00255 (see \ref SeedOptions). 00256 00257 The pixel type of the input image must be <tt>LessThanComparable</tt>. 00258 The pixel type of the output image must be large enough to hold the labels for all seeds. 00259 (typically, you will use <tt>UInt32</tt>). The function will label seeds by consecutive integers 00260 (starting from 1) and returns the largest label it used. 00261 00262 Pass \ref vigra::EightNeighborCode or \ref vigra::FourNeighborCode to determine the 00263 neighborhood where pixel values are compared. 00264 00265 The function uses accessors. 00266 00267 <b> Declarations:</b> 00268 00269 pass arguments explicitly: 00270 \code 00271 namespace vigra { 00272 template <class SrcIterator, class SrcAccessor, 00273 class DestIterator, class DestAccessor, 00274 class Neighborhood = EightNeighborCode> 00275 unsigned int 00276 generateWatershedSeeds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 00277 DestIterator upperleftd, DestAccessor da, 00278 Neighborhood neighborhood = EightNeighborCode(), 00279 SeedOptions const & options = SeedOptions()); 00280 } 00281 \endcode 00282 00283 use argument objects in conjunction with \ref ArgumentObjectFactories : 00284 \code 00285 namespace vigra { 00286 template <class SrcIterator, class SrcAccessor, 00287 class DestIterator, class DestAccessor, 00288 class Neighborhood = EightNeighborCode> 00289 unsigned int 00290 generateWatershedSeeds(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00291 pair<DestIterator, DestAccessor> dest, 00292 Neighborhood neighborhood = EightNeighborCode(), 00293 SeedOptions const & options = SeedOptions()); 00294 } 00295 \endcode 00296 00297 <b> Usage:</b> 00298 00299 <b>\#include</b> <vigra/watersheds.hxx><br> 00300 Namespace: vigra 00301 00302 For detailed examples see watershedsRegionGrowing(). 00303 */ 00304 doxygen_overloaded_function(template <...> unsigned int generateWatershedSeeds3D) 00305 00306 #if 0 00307 template <unsigned int N, class T1, class C1, class T2, class C2> 00308 class Neighborhood> 00309 unsigned int 00310 generateWatershedSeeds3D(MultiArrayView<N, T1, C1> in, MultiArrayView<N, T2, C2> out, 00311 Neighborhood neighborhood, 00312 SeedOptions const & options = SeedOptions()) 00313 { 00314 using namespace functor; 00315 00316 vigra_precondition(in.shape() == out.shape(), 00317 "generateWatershedSeeds3D(): Shape mismatch between input and output."); 00318 00319 vigra_precondition(options.mini != SeedOptions::LevelSets || 00320 options.thresholdIsValid<SrcType>(), 00321 "generateWatershedSeeds3D(): SeedOptions.levelSets() must be specified with threshold."); 00322 00323 MultiArray<N, UInt8> seeds(in.shape()); 00324 00325 if(options.mini == SeedOptions::LevelSets) 00326 { 00327 transformMultiArray(srcMultiArrayRange(in), destMultiArray(seeds), 00328 ifThenElse(Arg1() <= Param(options.thresh), Param(1), Param(0))); 00329 } 00330 else 00331 { 00332 localMinima(in, seeds, 00333 LocalMinmaxOptions().neighborhood(Neighborhood::DirectionCount) 00334 .markWith(1.0) 00335 .threshold(options.thresh) 00336 .allowAtBorder() 00337 .allowPlateaus(options.mini == SeedOptions::ExtendedMinima)); 00338 } 00339 00340 return labelVolumeWithBackground(srcMultiArrayRange(seeds), destMultiArray(out), 00341 neighborhood, 0); 00342 } 00343 00344 template <class SrcIterator, class SrcAccessor, 00345 class DestIterator, class DestAccessor> 00346 inline unsigned int 00347 generateWatershedSeeds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 00348 DestIterator upperleftd, DestAccessor da, 00349 SeedOptions const & options = SeedOptions()) 00350 { 00351 return generateWatershedSeeds(upperlefts, lowerrights, sa, upperleftd, da, 00352 EightNeighborCode(), options); 00353 } 00354 00355 template <class SrcIterator, class SrcAccessor, 00356 class DestIterator, class DestAccessor, 00357 class Neighborhood> 00358 inline unsigned int 00359 generateWatershedSeeds(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00360 pair<DestIterator, DestAccessor> dest, 00361 Neighborhood neighborhood, 00362 SeedOptions const & options = SeedOptions()) 00363 { 00364 return generateWatershedSeeds(src.first, src.second, src.third, 00365 dest.first, dest.second, 00366 neighborhood, options); 00367 } 00368 00369 template <class SrcIterator, class SrcAccessor, 00370 class DestIterator, class DestAccessor> 00371 inline unsigned int 00372 generateWatershedSeeds(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00373 pair<DestIterator, DestAccessor> dest, 00374 SeedOptions const & options = SeedOptions()) 00375 { 00376 return generateWatershedSeeds(src.first, src.second, src.third, 00377 dest.first, dest.second, 00378 EightNeighborCode(), options); 00379 } 00380 00381 #endif 00382 00383 /********************************************************/ 00384 /* */ 00385 /* watersheds3D */ 00386 /* */ 00387 /********************************************************/ 00388 00389 /** \brief Region Segmentation by means of the watershed algorithm. 00390 00391 <b> Declarations:</b> 00392 00393 pass arguments explicitly: 00394 \code 00395 namespace vigra { 00396 template <class SrcIterator, class SrcAccessor,class SrcShape, 00397 class DestIterator, class DestAccessor, 00398 class Neighborhood3D> 00399 unsigned int watersheds3D(SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa, 00400 DestIterator d_Iter, DestAccessor da, 00401 Neighborhood3D neighborhood3D); 00402 } 00403 \endcode 00404 00405 use argument objects in conjunction with \ref ArgumentObjectFactories : 00406 \code 00407 namespace vigra { 00408 template <class SrcIterator, class SrcAccessor,class SrcShape, 00409 class DestIterator, class DestAccessor, 00410 class Neighborhood3D> 00411 unsigned int watersheds3D(triple<SrcIterator, SrcShape, SrcAccessor> src, 00412 pair<DestIterator, DestAccessor> dest, 00413 Neighborhood3D neighborhood3D); 00414 } 00415 \endcode 00416 00417 use with 3D-Six-Neighborhood: 00418 \code 00419 namespace vigra { 00420 00421 template <class SrcIterator, class SrcAccessor,class SrcShape, 00422 class DestIterator, class DestAccessor> 00423 unsigned int watersheds3DSix(triple<SrcIterator, SrcShape, SrcAccessor> src, 00424 pair<DestIterator, DestAccessor> dest); 00425 00426 } 00427 \endcode 00428 00429 use with 3D-TwentySix-Neighborhood: 00430 \code 00431 namespace vigra { 00432 00433 template <class SrcIterator, class SrcAccessor,class SrcShape, 00434 class DestIterator, class DestAccessor> 00435 unsigned int watersheds3DTwentySix(triple<SrcIterator, SrcShape, SrcAccessor> src, 00436 pair<DestIterator, DestAccessor> dest); 00437 00438 } 00439 \endcode 00440 00441 This function implements the union-find version of the watershed algorithms 00442 as described in 00443 00444 J. Roerdink, R. Meijster: "<em>The watershed transform: definitions, algorithms, 00445 and parallelization strategies</em>", Fundamenta Informaticae, 41:187-228, 2000 00446 00447 The source volume is a boundary indicator such as the gradient magnitude 00448 of the trace of the \ref boundaryTensor(). Local minima of the boundary indicator 00449 are used as region seeds, and all other voxels are recursively assigned to the same 00450 region as their lowest neighbor. Pass \ref vigra::NeighborCode3DSix or 00451 \ref vigra::NeighborCode3DTwentySix to determine the neighborhood where voxel values 00452 are compared. The voxel type of the input volume must be <tt>LessThanComparable</tt>. 00453 The function uses accessors. 00454 00455 ...probably soon in VIGRA: 00456 Note that VIGRA provides an alternative implementation of the watershed transform via 00457 \ref seededRegionGrowing3D(). It is slower, but handles plateaus better 00458 and allows to keep a one pixel wide boundary between regions. 00459 00460 <b> Usage:</b> 00461 00462 <b>\#include</b> <vigra/watersheds3D.hxx><br> 00463 Namespace: vigra 00464 00465 Example: watersheds3D of the gradient magnitude. 00466 00467 \code 00468 typedef vigra::MultiArray<3,int> IntVolume; 00469 typedef vigra::MultiArray<3,double> DVolume; 00470 DVolume src(DVolume::difference_type(w,h,d)); 00471 IntVolume dest(IntVolume::difference_type(w,h,d)); 00472 00473 float gauss=1; 00474 00475 vigra::MultiArray<3, vigra::TinyVector<float,3> > temp(IntVolume::difference_type(w,h,d)); 00476 vigra::gaussianGradientMultiArray(srcMultiArrayRange(vol),destMultiArray(temp),gauss); 00477 00478 IntVolume::iterator temp_iter=temp.begin(); 00479 for(DVolume::iterator iter=src.begin(); iter!=src.end(); ++iter, ++temp_iter) 00480 *iter = norm(*temp_iter); 00481 00482 // find 6-connected regions 00483 int max_region_label = vigra::watersheds3DSix(srcMultiArrayRange(src), destMultiArray(dest)); 00484 00485 // find 26-connected regions 00486 max_region_label = vigra::watersheds3DTwentySix(srcMultiArrayRange(src), destMultiArray(dest)); 00487 00488 \endcode 00489 00490 <b> Required Interface:</b> 00491 00492 \code 00493 SrcIterator src_begin; 00494 SrcShape src_shape; 00495 DestIterator dest_begin; 00496 00497 SrcAccessor src_accessor; 00498 DestAccessor dest_accessor; 00499 00500 // compare src values 00501 src_accessor(src_begin) <= src_accessor(src_begin) 00502 00503 // set result 00504 int label; 00505 dest_accessor.set(label, dest_begin); 00506 \endcode 00507 */ 00508 doxygen_overloaded_function(template <...> unsigned int watersheds3D) 00509 00510 template <class SrcIterator, class SrcAccessor, class SrcShape, 00511 class DestIterator, class DestAccessor, 00512 class Neighborhood3D> 00513 unsigned int watersheds3D( SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa, 00514 DestIterator d_Iter, DestAccessor da, Neighborhood3D neighborhood3D) 00515 { 00516 //create temporary volume to store the DAG of directions to minima 00517 if ((int)Neighborhood3D::DirectionCount>7){ //If we have 3D-TwentySix Neighborhood 00518 00519 vigra::MultiArray<3,int> orientationVolume(srcShape); 00520 00521 preparewatersheds3D( s_Iter, srcShape, sa, 00522 destMultiArray(orientationVolume).first, destMultiArray(orientationVolume).second, 00523 neighborhood3D); 00524 00525 return watershedLabeling3D( srcMultiArray(orientationVolume).first, srcShape, srcMultiArray(orientationVolume).second, 00526 d_Iter, da, 00527 neighborhood3D); 00528 } 00529 else{ 00530 00531 vigra::MultiArray<3,unsigned char> orientationVolume(srcShape); 00532 00533 preparewatersheds3D( s_Iter, srcShape, sa, 00534 destMultiArray(orientationVolume).first, destMultiArray(orientationVolume).second, 00535 neighborhood3D); 00536 00537 return watershedLabeling3D( srcMultiArray(orientationVolume).first, srcShape, srcMultiArray(orientationVolume).second, 00538 d_Iter, da, 00539 neighborhood3D); 00540 } 00541 } 00542 00543 template <class SrcIterator, class SrcShape, class SrcAccessor, 00544 class DestIterator, class DestAccessor> 00545 inline unsigned int watersheds3DSix( vigra::triple<SrcIterator, SrcShape, SrcAccessor> src, 00546 vigra::pair<DestIterator, DestAccessor> dest) 00547 { 00548 return watersheds3D(src.first, src.second, src.third, dest.first, dest.second, NeighborCode3DSix()); 00549 } 00550 00551 template <class SrcIterator, class SrcShape, class SrcAccessor, 00552 class DestIterator, class DestAccessor> 00553 inline unsigned int watersheds3DTwentySix( vigra::triple<SrcIterator, SrcShape, SrcAccessor> src, 00554 vigra::pair<DestIterator, DestAccessor> dest) 00555 { 00556 return watersheds3D(src.first, src.second, src.third, dest.first, dest.second, NeighborCode3DTwentySix()); 00557 } 00558 00559 }//namespace vigra 00560 00561 #endif //VIGRA_watersheds3D_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|