[ 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,d);
    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.1 (3 Dec 2010)