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