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