[ 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 /*       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_watersheds3D_HXX
00038 #define VIGRA_watersheds3D_HXX
00039 
00040 #include "voxelneighborhood.hxx"
00041 #include "multi_array.hxx"
00042 
00043 namespace vigra
00044 {
00045 
00046 template <class SrcIterator, class SrcAccessor, class SrcShape,
00047           class DestIterator, class DestAccessor, class Neighborhood3D>
00048 int preparewatersheds3D( SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
00049                          DestIterator d_Iter, DestAccessor da, Neighborhood3D)
00050 {
00051     //basically needed for iteration and border-checks
00052     int w = srcShape[0], h = srcShape[1], d = srcShape[2];
00053     int x,y,z, local_min_count=0;
00054         
00055     //declare and define Iterators for all three dims at src
00056     SrcIterator zs = s_Iter;
00057     SrcIterator ys(zs);
00058     SrcIterator xs(ys);
00059         
00060     //Declare Iterators for all three dims at dest
00061     DestIterator zd = d_Iter;
00062         
00063     for(z = 0; z != d; ++z, ++zs.dim2(), ++zd.dim2())
00064     {
00065         ys = zs;
00066         DestIterator yd(zd);
00067         
00068         for(y = 0; y != h; ++y, ++ys.dim1(), ++yd.dim1())
00069         {
00070             xs = ys;
00071             DestIterator xd(yd);
00072 
00073             for(x = 0; x != w; ++x, ++xs.dim0(), ++xd.dim0())
00074             {
00075                 AtVolumeBorder atBorder = isAtVolumeBorder(x,y,z,w,h,d);
00076                 typename SrcAccessor::value_type v = sa(xs);
00077                 // the following choice causes minima to point
00078                 // to their lowest neighbor -- would this be better???
00079                 // typename SrcAccessor::value_type v = NumericTraits<typename SrcAccessor::value_type>::max();
00080                 int o = 0; // means center is minimum
00081                 typename SrcAccessor::value_type my_v = v;
00082                 if(atBorder == NotAtBorder)
00083                 {
00084 #if 0
00085                     NeighborhoodTraverser<SrcIterator, Neighborhood3D>  c(xs), cend(c);
00086 #endif /* #if 0 */
00087 
00088                     NeighborhoodCirculator<SrcIterator, Neighborhood3D>  c(xs), cend(c);
00089                     
00090                     do {
00091                         if(sa(c) < v)
00092                         {  
00093                             v = sa(c);
00094                             o = c.directionBit();
00095                         }
00096                         else if(sa(c) == my_v && my_v == v)
00097                         {
00098                             o =  o | c.directionBit();
00099                         }
00100                     }
00101                     while(++c != cend);
00102                 }
00103                 else
00104                 {
00105 #if 0
00106                     RestrictedNeighborhoodTraverser<SrcIterator, Neighborhood3D>  c(xs, atBorder), cend(c);
00107 #endif /* #if 0 */
00108 
00109                     RestrictedNeighborhoodCirculator<SrcIterator, Neighborhood3D>  c(xs, atBorder), cend(c);
00110                     do {
00111                         if(sa(c) < v)
00112                         {  
00113                             v = sa(c);
00114                             o = c.directionBit();
00115                         }
00116                         else if(sa(c) == my_v && my_v == v)
00117                         {
00118                             o =  o | c.directionBit();
00119                         }
00120                     }
00121                     while(++c != cend);
00122                 }
00123                 if (o==0) local_min_count++; 
00124                 da.set(o, xd);
00125             }//end x-iteration
00126         }//end y-iteration
00127     }//end z-iteration
00128     return local_min_count;
00129 }
00130 
00131 template <class SrcIterator, class SrcAccessor,class SrcShape,
00132           class DestIterator, class DestAccessor,
00133           class Neighborhood3D>
00134 unsigned int watershedLabeling3D( SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
00135                                   DestIterator d_Iter, DestAccessor da,
00136                                   Neighborhood3D)
00137 {
00138     //basically needed for iteration and border-checks
00139     int w = srcShape[0], h = srcShape[1], d = srcShape[2];
00140     int x,y,z, i;
00141         
00142     //declare and define Iterators for all three dims at src
00143     SrcIterator zs = s_Iter;
00144     SrcIterator ys(zs);
00145     SrcIterator xs(ys);
00146         
00147     // temporary image to store region labels
00148     typedef vigra::MultiArray<3,int> LabelVolume;
00149     typedef LabelVolume::traverser LabelTraverser;
00150     
00151     LabelVolume labelvolume(srcShape);
00152         
00153     //Declare traversers for all three dims at target
00154     LabelTraverser zt = labelvolume.traverser_begin();
00155     LabelTraverser yt(zt);
00156     LabelTraverser xt(yt);
00157 
00158     // Kovalevsky's clever idea to use
00159     // image iterator and scan order iterator simultaneously
00160     // memory order indicates label
00161     LabelVolume::iterator label = labelvolume.begin();
00162     
00163     // initialize the neighborhood traversers
00164 
00165 #if 0
00166     NeighborOffsetTraverser<Neighborhood3D> nc(Neighborhood3D::CausalFirst);
00167     NeighborOffsetTraverser<Neighborhood3D> nce(Neighborhood3D::CausalLast);
00168 #endif /* #if 0 */
00169 
00170     NeighborOffsetCirculator<Neighborhood3D> nc(Neighborhood3D::CausalFirst);
00171     NeighborOffsetCirculator<Neighborhood3D> nce(Neighborhood3D::CausalLast);
00172     ++nce;
00173     // pass 1: scan image from upper left front to lower right back
00174     // to find connected components
00175 
00176     // Each component will be represented by a tree of pixels. Each
00177     // pixel contains the scan order address of its parent in the
00178     // tree.  In order for pass 2 to work correctly, the parent must
00179     // always have a smaller scan order address than the child.
00180     // Therefore, we can merge trees only at their roots, because the
00181     // root of the combined tree must have the smallest scan order
00182     // address among all the tree's pixels/ nodes.  The root of each
00183     // tree is distinguished by pointing to itself (it contains its
00184     // own scan order address). This condition is enforced whenever a
00185     // new region is found or two regions are merged
00186     i=0;
00187     for(z = 0; z != d; ++z, ++zs.dim2(), ++zt.dim2())
00188     {
00189         ys = zs;
00190         yt = zt;
00191 
00192         for(y = 0; y != h; ++y, ++ys.dim1(), ++yt.dim1())
00193         {
00194             xs = ys;
00195             xt = yt;
00196 
00197             for(x = 0; x != w; ++x, ++xs.dim0(), ++xt.dim0(), ++i)
00198             {
00199                 *xt = i; // default: new region    
00200 
00201                 //queck whether there is a special borde threatment to be used or not
00202                 AtVolumeBorder atBorder = isAtVolumeBorderCausal(x,y,z,w,h,z);
00203                     
00204                 //We are not at the border!
00205                 if(atBorder == NotAtBorder)
00206                 {
00207 
00208 #if 0
00209                     nc = NeighborOffsetTraverser<Neighborhood3D>(Neighborhood3D::CausalFirst);
00210 #endif /* #if 0 */
00211 
00212                     nc = NeighborOffsetCirculator<Neighborhood3D>(Neighborhood3D::CausalFirst);
00213                 
00214                     do
00215                     {            
00216                         //   Direction of NTraversr       Neighbor's direction bit is pointing
00217                         // = Direction of voxel           towards us?
00218                         if((*xs & nc.directionBit()) || (xs[*nc] & nc.oppositeDirectionBit()))
00219                         {
00220                             int neighborLabel = xt[*nc];
00221 
00222                             // find the root label of a label tree
00223                             while(neighborLabel != label[neighborLabel])
00224                             {
00225                                 neighborLabel = label[neighborLabel];
00226                             }
00227 
00228                             if(neighborLabel < *xt) // always keep the smallest among the possible neighbor labels
00229                             {
00230                                 label[*xt] = neighborLabel;
00231                                 *xt = neighborLabel;
00232                             }
00233                             else
00234                             {
00235                                 label[neighborLabel] = *xt;
00236                             }
00237                         }
00238                         ++nc;
00239                     }while(nc!=nce);
00240                 }
00241                 //we are at a border - handle this!!
00242                 else
00243                 {
00244 
00245 #if 0
00246                     nc = NeighborOffsetTraverser<Neighborhood3D>(Neighborhood3D::nearBorderDirectionsCausal(atBorder,0));
00247 #endif /* #if 0 */
00248 
00249                     nc = NeighborOffsetCirculator<Neighborhood3D>(Neighborhood3D::nearBorderDirectionsCausal(atBorder,0));
00250                     int j=0;
00251                     while(nc.direction() != Neighborhood3D::Error)
00252                     {
00253                         //   Direction of NTraversr       Neighbor's direction bit is pointing
00254                         // = Direction of voxel           towards us?
00255                         if((*xs & nc.directionBit()) || (xs[*nc] & nc.oppositeDirectionBit()))
00256                         {
00257                             int neighborLabel = xt[*nc];
00258 
00259                             // find the root label of a label tree
00260                             while(neighborLabel != label[neighborLabel])
00261                             {
00262                                 neighborLabel = label[neighborLabel];
00263                             }
00264 
00265                             if(neighborLabel < *xt) // always keep the smallest among the possible neighbor labels
00266                             {
00267                                 label[*xt] = neighborLabel;
00268                                 *xt = neighborLabel;
00269                             }
00270                             else
00271                             {
00272                                 label[neighborLabel] = *xt;
00273                             }
00274                         }
00275                         nc.turnTo(Neighborhood3D::nearBorderDirectionsCausal(atBorder,++j));
00276                     }
00277                 }
00278             }
00279         }
00280     }
00281 
00282     // pass 2: assign one label to each region (tree)
00283     // so that labels form a consecutive sequence 1, 2, ...
00284     DestIterator zd = d_Iter;
00285 
00286     unsigned int count = 0; 
00287 
00288     i= 0;
00289 
00290     for(z=0; z != d; ++z, ++zd.dim2())
00291     {
00292         DestIterator yd(zd);
00293 
00294         for(y=0; y != h; ++y, ++yd.dim1())
00295         {
00296             DestIterator xd(yd);
00297 
00298             for(x = 0; x != w; ++x, ++xd.dim0(), ++i)
00299             {
00300                 if(label[i] == i)
00301                 {
00302                     label[i] = count++;
00303                 }
00304                 else
00305                 {
00306                     label[i] = label[label[i]]; // compress trees
00307                 }
00308                 da.set(label[i]+1, xd);
00309             }
00310         }
00311     }
00312     return count;
00313 }
00314 
00315 
00316 /** \addtogroup SeededRegionGrowing Region Segmentation Algorithms
00317     Region growing, watersheds, and voronoi tesselation
00318 */
00319 //@{
00320 
00321 /********************************************************/
00322 /*                                                      */
00323 /*                     watersheds3D                     */
00324 /*                                                      */
00325 /********************************************************/
00326 
00327 /** \brief Region Segmentation by means of the watershed algorithm.
00328 
00329     <b> Declarations:</b>
00330 
00331     pass arguments explicitly:
00332     \code
00333     namespace vigra {
00334         template <class SrcIterator, class SrcAccessor,class SrcShape,
00335                   class DestIterator, class DestAccessor,
00336                   class Neighborhood3D>
00337         unsigned int watersheds3D(SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
00338                                   DestIterator d_Iter, DestAccessor da,
00339                                   Neighborhood3D neighborhood3D);
00340     }
00341     \endcode
00342 
00343     use argument objects in conjunction with \ref ArgumentObjectFactories :
00344     \code
00345     namespace vigra {
00346         template <class SrcIterator, class SrcAccessor,class SrcShape,
00347                   class DestIterator, class DestAccessor,
00348                   class Neighborhood3D>
00349         unsigned int watersheds3D(triple<SrcIterator, SrcShape, SrcAccessor> src,
00350                                   pair<DestIterator, DestAccessor> dest,
00351                                   Neighborhood3D neighborhood3D);
00352     }
00353     \endcode
00354 
00355     use with 3D-Six-Neighborhood:
00356     \code
00357     namespace vigra {    
00358     
00359         template <class SrcIterator, class SrcAccessor,class SrcShape,
00360                   class DestIterator, class DestAccessor>
00361         unsigned int watersheds3DSix(triple<SrcIterator, SrcShape, SrcAccessor> src,
00362                                      pair<DestIterator, DestAccessor> dest);
00363                                     
00364     }
00365     \endcode
00366 
00367     use with 3D-TwentySix-Neighborhood:
00368     \code
00369     namespace vigra {    
00370     
00371         template <class SrcIterator, class SrcAccessor,class SrcShape,
00372                   class DestIterator, class DestAccessor>
00373         unsigned int watersheds3DTwentySix(triple<SrcIterator, SrcShape, SrcAccessor> src,
00374                                            pair<DestIterator, DestAccessor> dest);
00375                                     
00376     }
00377     \endcode
00378 
00379     This function implements the union-find version of the watershed algorithms
00380     as described in
00381 
00382     J. Roerdink, R. Meijster: "<em>The watershed transform: definitions, algorithms,
00383     and parallelization stretegies</em>", Fundamenta Informaticae, 41:187-228, 2000
00384 
00385     The source volume is a boundary indicator such as the gradient magnitude
00386     of the trace of the \ref boundaryTensor(). Local minima of the boundary indicator
00387     are used as region seeds, and all other voxels are recursively assigned to the same 
00388     region as their lowest neighbor. Pass \ref vigra::NeighborCode3DSix or 
00389     \ref vigra::NeighborCode3DTwentySix to determine the neighborhood where voxel values 
00390     are compared. The voxel type of the input volume must be <tt>LessThanComparable</tt>.
00391     The function uses accessors. 
00392     
00393     ...probably soon in VIGRA:
00394     Note that VIGRA provides an alternative implementaion of the watershed transform via
00395     \ref seededRegionGrowing3D(). It is slower, but handles plateaus better 
00396     and allows to keep a one pixel wide boundary between regions.
00397     
00398     <b> Usage:</b>
00399 
00400     <b>\#include</b> <<a href="watersheds3D_8hxx-source.html">vigra/watersheds3D.hxx</a>><br>
00401     Namespace: vigra
00402 
00403     Example: watersheds3D of the gradient magnitude.
00404 
00405     \code
00406     typedef vigra::MultiArray<3,int> IntVolume;
00407     typedef vigra::MultiArray<3,double> DVolume;
00408     DVolume src(DVolume::difference_type(w,h,d));
00409     IntVolume dest(IntVolume::difference_type(w,h,d));
00410 
00411     float gauss=1;
00412 
00413     vigra::MultiArray<3, vigra::TinyVector<float,3> > temp(IntVolume::difference_type(w,h,d));
00414     vigra::gaussianGradientMultiArray(srcMultiArrayRange(vol),destMultiArray(temp),gauss);
00415 
00416     IntVolume::iterator temp_iter=temp.begin();
00417     for(DVolume::iterator iter=src.begin(); iter!=src.end(); ++iter, ++temp_iter)
00418         *iter = norm(*temp_iter);
00419     
00420     // find 6-connected regions
00421     int max_region_label = vigra::watersheds3DSix(srcMultiArrayRange(src), destMultiArray(dest));
00422 
00423     // find 26-connected regions
00424     max_region_label = vigra::watersheds3DTwentySix(srcMultiArrayRange(src), destMultiArray(dest));
00425     
00426     \endcode
00427 
00428     <b> Required Interface:</b>
00429 
00430     \code
00431     SrcIterator src_begin;
00432     SrcShape src_shape;
00433     DestIterator dest_begin;
00434 
00435     SrcAccessor src_accessor;
00436     DestAccessor dest_accessor;
00437     
00438     // compare src values
00439     src_accessor(src_begin) <= src_accessor(src_begin)
00440 
00441     // set result
00442     int label;
00443     dest_accessor.set(label, dest_begin);
00444     \endcode
00445 */
00446 doxygen_overloaded_function(template <...> unsigned int watersheds3D)
00447 
00448 template <class SrcIterator, class SrcAccessor, class SrcShape,
00449           class DestIterator, class DestAccessor,
00450           class Neighborhood3D>
00451 unsigned int watersheds3D( SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
00452                            DestIterator d_Iter, DestAccessor da, Neighborhood3D neighborhood3D)
00453 {
00454     //create temporary volume to store the DAG of directions to minima
00455     if ((int)Neighborhood3D::DirectionCount>7){  //If we have 3D-TwentySix Neighborhood
00456         
00457         vigra::MultiArray<3,int> orientationVolume(srcShape);
00458 
00459         int local_min_count= preparewatersheds3D( s_Iter, srcShape, sa, 
00460                                                   destMultiArray(orientationVolume).first, destMultiArray(orientationVolume).second,
00461                                                   neighborhood3D);
00462      
00463         return watershedLabeling3D( srcMultiArray(orientationVolume).first, srcShape, srcMultiArray(orientationVolume).second,
00464                                     d_Iter, da,
00465                                     neighborhood3D);
00466     }
00467     else{
00468                 
00469         vigra::MultiArray<3,unsigned char> orientationVolume(srcShape);
00470 
00471         int local_min_count= preparewatersheds3D( s_Iter, srcShape, sa, 
00472                                                   destMultiArray(orientationVolume).first, destMultiArray(orientationVolume).second,
00473                                                   neighborhood3D);
00474      
00475         return watershedLabeling3D( srcMultiArray(orientationVolume).first, srcShape, srcMultiArray(orientationVolume).second,
00476                                     d_Iter, da,
00477                                     neighborhood3D);
00478     }
00479 }
00480 
00481 template <class SrcIterator, class SrcShape, class SrcAccessor,
00482           class DestIterator, class DestAccessor>
00483 inline unsigned int watersheds3DSix( vigra::triple<SrcIterator, SrcShape, SrcAccessor> src, 
00484                                      vigra::pair<DestIterator, DestAccessor> dest)
00485 {
00486     return watersheds3D(src.first, src.second, src.third, dest.first, dest.second, NeighborCode3DSix());
00487 }
00488 
00489 template <class SrcIterator, class SrcShape, class SrcAccessor,
00490           class DestIterator, class DestAccessor>
00491 inline unsigned int watersheds3DTwentySix( vigra::triple<SrcIterator, SrcShape, SrcAccessor> src, 
00492                                            vigra::pair<DestIterator, DestAccessor> dest)
00493 {
00494     return watersheds3D(src.first, src.second, src.third, dest.first, dest.second, NeighborCode3DTwentySix());
00495 }
00496 
00497 }//namespace vigra
00498 
00499 #endif //VIGRA_watersheds3D_HXX

© 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.6.0 (5 Nov 2009)