[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

vigra/labelvolume.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_LABELVOLUME_HXX
00039 #define VIGRA_LABELVOLUME_HXX
00040 
00041 
00042 #include "voxelneighborhood.hxx"
00043 #include "multi_array.hxx"
00044 
00045 namespace vigra{
00046 
00047 /** \addtogroup Labeling Connected Components Labeling
00048      The 3-dimensional connected components algorithms may use either 6 or 26 connectivity.
00049      By means of a functor the merge criterium can be defined arbitrarily.
00050 */
00051 //@{
00052 
00053 /********************************************************/
00054 /*                                                      */
00055 /*                        labelVolume                   */
00056 /*                                                      */
00057 /********************************************************/
00058 
00059 /** \brief Find the connected components of a segmented volume.
00060 
00061     <b> Declarations:</b>
00062 
00063     pass arguments explicitly:
00064     \code
00065     namespace vigra {
00066 
00067         template <class SrcIterator, class SrcAccessor,class SrcShape,
00068                   class DestIterator, class DestAccessor,
00069                   class Neighborhood3D>
00070         unsigned int labelVolume(SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
00071                                  DestIterator d_Iter, DestAccessor da,
00072                                  Neighborhood3D neighborhood3D);
00073 
00074         template <class SrcIterator, class SrcAccessor,class SrcShape,
00075                           class DestIterator, class DestAccessor,
00076                           class Neighborhood3D, class EqualityFunctor>
00077         unsigned int labelVolume(SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
00078                                  DestIterator d_Iter, DestAccessor da,
00079                                  Neighborhood3D neighborhood3D, EqualityFunctor equal);
00080 
00081     }
00082     \endcode
00083 
00084     use argument objects in conjunction with \ref ArgumentObjectFactories :
00085     \code
00086     namespace vigra {
00087 
00088         template <class SrcIterator, class SrcAccessor,class SrcShape,
00089                   class DestIterator, class DestAccessor,
00090                   class Neighborhood3D>
00091         unsigned int labelVolume(triple<SrcIterator, SrcShape, SrcAccessor> src,
00092                                  pair<DestIterator, DestAccessor> dest,
00093                                  Neighborhood3D neighborhood3D);
00094 
00095         template <class SrcIterator, class SrcAccessor,class SrcShape,
00096                  class DestIterator, class DestAccessor,
00097                  class Neighborhood3D, class EqualityFunctor>
00098         unsigned int labelVolume(triple<SrcIterator, SrcShape, SrcAccessor> src,
00099                                  pair<DestIterator, DestAccessor> dest,
00100                                  Neighborhood3D neighborhood3D, EqualityFunctor equal);
00101 
00102     }
00103     \endcode
00104     
00105     use with 3D-Six-Neighborhood:
00106     \code
00107     namespace vigra {    
00108     
00109         template <class SrcIterator, class SrcAccessor,class SrcShape,
00110                   class DestIterator, class DestAccessor>
00111         unsigned int labelVolumeSix(triple<SrcIterator, SrcShape, SrcAccessor> src,
00112                                     pair<DestIterator, DestAccessor> dest);
00113                                     
00114     }
00115     \endcode
00116 
00117     Connected components are defined as regions with uniform voxel
00118     values. Thus, <TT>SrcAccessor::value_type</TT> either must be
00119     equality comparable (first form), or an EqualityFunctor must be
00120     provided that realizes the desired predicate (second form). The
00121     destination's value type should be large enough to hold the labels
00122     without overflow. Region numbers will be a consecutive sequence
00123     starting with one and ending with the region number returned by
00124     the function (inclusive).
00125 
00126     Return:  the number of regions found (= largest region label)
00127 
00128     <b> Usage:</b>
00129 
00130     <b>\#include</b> <<a href="labelvolume_8hxx-source.html">vigra/labelvolume.hxx</a>><br>
00131     Namespace: vigra
00132 
00133     \code
00134     typedef vigra::MultiArray<3,int> IntVolume;
00135     IntVolume src(IntVolume::difference_type(w,h,d));
00136     IntVolume dest(IntVolume::difference_type(w,h,d));
00137     
00138     // find 6-connected regions
00139     int max_region_label = vigra::labelVolumeSix(srcMultiArrayRange(src), destMultiArray(dest));
00140 
00141     // find 26-connected regions
00142     int max_region_label = vigra::labelVolume(srcMultiArrayRange(src), destMultiArray(dest), NeighborCode3DTwentySix());
00143     \endcode
00144 
00145     <b> Required Interface:</b>
00146 
00147     \code
00148     SrcIterator src_begin;
00149     SrcShape shape;
00150     DestIterator dest_begin;
00151 
00152     SrcAccessor src_accessor;
00153     DestAccessor dest_accessor;
00154 
00155     SrcAccessor::value_type u = src_accessor(src_begin);
00156 
00157     u == u                      // first form
00158 
00159     EqualityFunctor equal;      // second form
00160     equal(u, u)                 // second form
00161 
00162     int i;
00163     dest_accessor.set(i, dest_begin);
00164     \endcode
00165 
00166 */
00167 doxygen_overloaded_function(template <...> unsigned int labelVolume)
00168 
00169 template <class SrcIterator, class SrcAccessor,class SrcShape,
00170           class DestIterator, class DestAccessor,
00171           class Neighborhood3D>
00172 unsigned int labelVolume(SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
00173                          DestIterator d_Iter, DestAccessor da,
00174                          Neighborhood3D neighborhood3D)
00175 {
00176         return labelVolume(s_Iter, srcShape, sa, d_Iter, da, neighborhood3D, std::equal_to<typename SrcAccessor::value_type>());
00177 }
00178 
00179 template <class SrcIterator, class SrcAccessor,class SrcShape,
00180           class DestIterator, class DestAccessor,
00181           class Neighborhood3D>
00182 unsigned int labelVolume(triple<SrcIterator, SrcShape, SrcAccessor> src,
00183                          pair<DestIterator, DestAccessor> dest,
00184                          Neighborhood3D neighborhood3D)
00185 {
00186     return labelVolume(src.first, src.second, src.third, dest.first, dest.second, neighborhood3D, std::equal_to<typename SrcAccessor::value_type>());
00187 }
00188 
00189 template <class SrcIterator, class SrcAccessor,class SrcShape,
00190           class DestIterator, class DestAccessor,
00191           class Neighborhood3D, class EqualityFunctor>
00192 unsigned int labelVolume(triple<SrcIterator, SrcShape, SrcAccessor> src,
00193                          pair<DestIterator, DestAccessor> dest,
00194                          Neighborhood3D neighborhood3D, EqualityFunctor equal)
00195 {
00196     return labelVolume(src.first, src.second, src.third, dest.first, dest.second, neighborhood3D, equal);
00197 }
00198 
00199 template <class SrcIterator, class SrcAccessor,class SrcShape,
00200           class DestIterator, class DestAccessor,
00201           class Neighborhood3D, class EqualityFunctor>
00202 unsigned int labelVolume(SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
00203                          DestIterator d_Iter, DestAccessor da,
00204                          Neighborhood3D, EqualityFunctor equal)
00205 {
00206 
00207     //basically needed for iteration and border-checks
00208     int w = srcShape[0], h = srcShape[1], d = srcShape[2];
00209     int x,y,z, i;
00210         
00211     //declare and define Iterators for all three dims at src
00212     SrcIterator zs = s_Iter;
00213     SrcIterator ys(zs);
00214     SrcIterator xs(ys);
00215         
00216     // temporary image to store region labels
00217     typedef vigra::MultiArray<3,int> LabelVolume;
00218     typedef LabelVolume::traverser LabelTraverser;
00219     
00220     LabelVolume labelvolume(srcShape);
00221         
00222     //Declare traversers for all three dims at target
00223     LabelTraverser zt = labelvolume.traverser_begin();
00224     LabelTraverser yt(zt);
00225     LabelTraverser xt(yt);
00226 
00227     // Kovalevsky's clever idea to use
00228     // image iterator and scan order iterator simultaneously
00229     // memory order indicates label
00230     LabelVolume::iterator label = labelvolume.begin();
00231     
00232     // initialize the neighborhood traversers
00233 
00234 #if 0
00235     NeighborOffsetTraverser<Neighborhood3D> nc(Neighborhood3D::CausalFirst);
00236     NeighborOffsetTraverser<Neighborhood3D> nce(Neighborhood3D::CausalLast);
00237 #endif /* #if 0 */
00238 
00239     NeighborOffsetCirculator<Neighborhood3D> nc(Neighborhood3D::CausalFirst);
00240     NeighborOffsetCirculator<Neighborhood3D> nce(Neighborhood3D::CausalLast);
00241     ++nce;
00242     // pass 1: scan image from upper left front to lower right back
00243     // to find connected components
00244 
00245     // Each component will be represented by a tree of pixels. Each
00246     // pixel contains the scan order address of its parent in the
00247     // tree.  In order for pass 2 to work correctly, the parent must
00248     // always have a smaller scan order address than the child.
00249     // Therefore, we can merge trees only at their roots, because the
00250     // root of the combined tree must have the smallest scan order
00251     // address among all the tree's pixels/ nodes.  The root of each
00252     // tree is distinguished by pointing to itself (it contains its
00253     // own scan order address). This condition is enforced whenever a
00254     // new region is found or two regions are merged
00255     i=0;
00256     for(z = 0; z != d; ++z, ++zs.dim2(), ++zt.dim2())
00257     {
00258         ys = zs;
00259         yt = zt;
00260 
00261         for(y = 0; y != h; ++y, ++ys.dim1(), ++yt.dim1())
00262         {
00263             xs = ys;
00264             xt = yt;
00265 
00266             for(x = 0; x != w; ++x, ++xs.dim0(), ++xt.dim0(), ++i)
00267             {
00268                 *xt = i; // default: new region    
00269 
00270                 //queck whether there is a special borde threatment to be used or not
00271                 AtVolumeBorder atBorder = isAtVolumeBorderCausal(x,y,z,w,h,z);
00272                     
00273                 //We are not at the border!
00274                 if(atBorder == NotAtBorder)
00275                 {
00276                     
00277 
00278 #if 0
00279                     nc = NeighborOffsetTraverser<Neighborhood3D>(Neighborhood3D::CausalFirst);
00280 #endif /* #if 0 */
00281 
00282                     nc = NeighborOffsetCirculator<Neighborhood3D>(Neighborhood3D::CausalFirst);
00283                 
00284                     do
00285                     {            
00286                         // if colors are equal
00287                         if(equal(*xs, xs[*nc]))
00288                         {
00289                             int neighborLabel = xt[*nc];
00290 
00291                             // find the root label of a label tree
00292                             while(neighborLabel != label[neighborLabel])
00293                             {
00294                                 neighborLabel = label[neighborLabel];
00295                             }
00296 
00297                             if(neighborLabel < *xt) // always keep the smallest among the possible neighbor labels
00298                             {
00299                                 label[*xt] = neighborLabel;
00300                                 *xt = neighborLabel;
00301                             }
00302                             else
00303                             {
00304                                 label[neighborLabel] = *xt;
00305                             }
00306                         }
00307                         ++nc;
00308                     }while(nc!=nce);
00309                 }
00310                 //we are at a border - handle this!!
00311                 else
00312                 {
00313 
00314 #if 0
00315                     nc = NeighborOffsetTraverser<Neighborhood3D>(Neighborhood3D::nearBorderDirectionsCausal(atBorder,0));
00316 #endif /* #if 0 */
00317 
00318                     nc = NeighborOffsetCirculator<Neighborhood3D>(Neighborhood3D::nearBorderDirectionsCausal(atBorder,0));
00319                     int j=0;
00320                     while(nc.direction() != Neighborhood3D::Error)
00321                     {
00322                         //   colors equal???
00323                         if(equal(*xs, xs[*nc]))
00324                         {
00325                             int neighborLabel = xt[*nc];
00326 
00327                             // find the root label of a label tree
00328                             while(neighborLabel != label[neighborLabel])
00329                             {
00330                                 neighborLabel = label[neighborLabel];
00331                             }
00332 
00333                             if(neighborLabel < *xt) // always keep the smallest among the possible neighbor labels
00334                             {
00335                                 label[*xt] = neighborLabel;
00336                                 *xt = neighborLabel;
00337                             }
00338                             else
00339                             {
00340                                 label[neighborLabel] = *xt;
00341                             }
00342                         }
00343                         nc.turnTo(Neighborhood3D::nearBorderDirectionsCausal(atBorder,++j));
00344                     }
00345                 }
00346             }
00347         }
00348     }
00349 
00350     // pass 2: assign one label to each region (tree)
00351     // so that labels form a consecutive sequence 1, 2, ...
00352     DestIterator zd = d_Iter;
00353 
00354     unsigned int count = 0; 
00355 
00356     i= 0;
00357 
00358     for(z=0; z != d; ++z, ++zd.dim2())
00359     {
00360         DestIterator yd(zd);
00361 
00362         for(y=0; y != h; ++y, ++yd.dim1())
00363         {
00364             DestIterator xd(yd);
00365 
00366             for(x = 0; x != w; ++x, ++xd.dim0(), ++i)
00367             {
00368                 if(label[i] == i)
00369                 {
00370                         label[i] = count++;
00371                 }
00372                 else
00373                 {
00374                         label[i] = label[label[i]]; // compress trees
00375                 }
00376                 da.set(label[i]+1, xd);
00377             }
00378         }
00379     }
00380     return count;
00381 }
00382 
00383 /********************************************************/
00384 /*                                                      */
00385 /*                    labelVolumeSix                    */
00386 /*                                                      */
00387 /********************************************************/
00388 
00389 /** \brief Find the connected components of a segmented volume
00390      using the 6-neighborhood.
00391      
00392      See \ref labelVolume() for detailed documentation.
00393 
00394 */
00395 template <class SrcIterator, class SrcAccessor,class SrcShape,
00396           class DestIterator, class DestAccessor>
00397 unsigned int labelVolumeSix(triple<SrcIterator, SrcShape, SrcAccessor> src,
00398                             pair<DestIterator, DestAccessor> dest)
00399 {
00400     return labelVolume(src.first, src.second, src.third, dest.first, dest.second, NeighborCode3DSix(), std::equal_to<typename SrcAccessor::value_type>());
00401 }
00402 
00403 
00404 
00405 
00406 /********************************************************/
00407 /*                                                      */
00408 /*               labelVolumeWithBackground              */
00409 /*                                                      */
00410 /********************************************************/
00411 
00412 /** \brief Find the connected components of a segmented volume,
00413      excluding the background from labeling.
00414 
00415     <b> Declarations:</b>
00416 
00417     pass arguments explicitly:
00418     \code
00419     namespace vigra {
00420 
00421         template <class SrcIterator, class SrcAccessor,class SrcShape,
00422                           class DestIterator, class DestAccessor,
00423                           class Neighborhood3D, class ValueType>
00424         unsigned int labelVolumeWithBackground(    SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
00425                                                           DestIterator d_Iter, DestAccessor da,
00426                                                           Neighborhood3D neighborhood3D, ValueType background_value);
00427 
00428         template <class SrcIterator, class SrcAccessor,class SrcShape,
00429                           class DestIterator, class DestAccessor,
00430                           class Neighborhood3D, class ValueType, class EqualityFunctor>
00431         unsigned int labelVolumeWithBackground(    SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
00432                                                             DestIterator d_Iter, DestAccessor da,
00433                                                           Neighborhood3D neighborhood3D, ValueType background_value,
00434                                                             EqualityFunctor equal);
00435 
00436     }
00437     \endcode
00438 
00439     use argument objects in conjunction with \ref ArgumentObjectFactories :
00440     \code
00441     namespace vigra {
00442 
00443         template <class SrcIterator, class SrcAccessor,class SrcShape,
00444                           class DestIterator, class DestAccessor,
00445                           class Neighborhood3D, class ValueType>
00446         unsigned int labelVolumeWithBackground(    triple<SrcIterator, SrcShape, SrcAccessor> src,
00447                                                           pair<DestIterator, DestAccessor> dest,
00448                                                           Neighborhood3D neighborhood3D, ValueType background_value);
00449 
00450         template <class SrcIterator, class SrcAccessor,class SrcShape,
00451                           class DestIterator, class DestAccessor,
00452                           class Neighborhood3D, class ValueType, class EqualityFunctor>
00453         unsigned int labelVolumeWithBackground(    triple<SrcIterator, SrcShape, SrcAccessor> src,
00454                                                         pair<DestIterator, DestAccessor> dest,
00455                                                         Neighborhood3D neighborhood3D, ValueType background_value,
00456                                                         EqualityFunctor equal);
00457 
00458     }
00459     \endcode
00460 
00461     Connected components are defined as regions with uniform voxel
00462     values. Thus, <TT>SrcAccessor::value_type</TT> either must be
00463     equality comparable (first form), or an EqualityFunctor must be
00464     provided that realizes the desired predicate (second form). All
00465     voxel equal to the given '<TT>background_value</TT>' are ignored
00466     when determining connected components and remain untouched in the
00467     destination volume.
00468 
00469     The destination's value type should be large enough to hold the
00470     labels without overflow. Region numbers will be a consecutive
00471     sequence starting with one and ending with the region number
00472     returned by the function (inclusive).
00473 
00474     Return:  the number of regions found (= largest region label)
00475 
00476     <b> Usage:</b>
00477 
00478     <b>\#include</b> <<a href="labelvolume_8hxx-source.html">vigra/labelvolume.hxx</a>><br>
00479     Namespace: vigra
00480 
00481     \code
00482     typedef vigra::MultiArray<3,int> IntVolume;
00483     IntVolume src(IntVolume::difference_type(w,h,d));
00484     IntVolume dest(IntVolume::difference_type(w,h,d));
00485 
00486     // find 6-connected regions
00487     int max_region_label = vigra::labelVolumeWithBackground(
00488     srcMultiArrayRange(src), destMultiArray(dest), NeighborCode3DSix(), 0);
00489     \endcode
00490 
00491     <b> Required Interface:</b>
00492 
00493     \code
00494     SrcIterator src_begin;
00495     SrcShape shape;
00496     DestIterator dest_begin;
00497 
00498     SrcAccessor src_accessor;
00499     DestAccessor dest_accessor;
00500 
00501     SrcAccessor::value_type u = src_accessor(src_begin);
00502 
00503     u == u                      // first form
00504 
00505     EqualityFunctor equal;      // second form
00506     equal(u, u)                 // second form
00507 
00508     int i;
00509     dest_accessor.set(i, dest_begin);
00510     \endcode
00511 
00512 */
00513 doxygen_overloaded_function(template <...> unsigned int labelVolumeWithBackground)
00514 
00515 template <class SrcIterator, class SrcAccessor,class SrcShape,
00516           class DestIterator, class DestAccessor,
00517           class Neighborhood3D,
00518           class ValueType>
00519 unsigned int labelVolumeWithBackground(SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
00520                                        DestIterator d_Iter, DestAccessor da,
00521                                        Neighborhood3D neighborhood3D, ValueType backgroundValue)
00522 {
00523     return labelVolumeWithBackground(s_Iter, srcShape, sa, d_Iter, da, neighborhood3D, backgroundValue, std::equal_to<typename SrcAccessor::value_type>());
00524 }
00525 
00526 template <class SrcIterator, class SrcAccessor,class SrcShape,
00527           class DestIterator, class DestAccessor,
00528           class Neighborhood3D,
00529           class ValueType>
00530 unsigned int labelVolumeWithBackground(triple<SrcIterator, SrcShape, SrcAccessor> src,
00531                                        pair<DestIterator, DestAccessor> dest,
00532                                        Neighborhood3D neighborhood3D, ValueType backgroundValue)
00533 {
00534     return labelVolumeWithBackground(src.first, src.second, src.third, dest.first, dest.second, neighborhood3D, backgroundValue, std::equal_to<typename SrcAccessor::value_type>());
00535 }
00536 
00537 template <class SrcIterator, class SrcAccessor,class SrcShape,
00538           class DestIterator, class DestAccessor,
00539           class Neighborhood3D,
00540           class ValueType, class EqualityFunctor>
00541 unsigned int labelVolumeWithBackground(triple<SrcIterator, SrcShape, SrcAccessor> src,
00542                                        pair<DestIterator, DestAccessor> dest,
00543                                        Neighborhood3D neighborhood3D, ValueType backgroundValue, EqualityFunctor equal)
00544 {
00545     return labelVolumeWithBackground(src.first, src.second, src.third, dest.first, dest.second, neighborhood3D, backgroundValue, equal);
00546 }
00547  
00548 template <class SrcIterator, class SrcAccessor,class SrcShape,
00549           class DestIterator, class DestAccessor,
00550           class Neighborhood3D,
00551           class ValueType, class EqualityFunctor>
00552 unsigned int labelVolumeWithBackground(SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa,
00553                                        DestIterator d_Iter, DestAccessor da,
00554                                        Neighborhood3D,
00555                                        ValueType backgroundValue, EqualityFunctor equal)
00556 {
00557 
00558     //basically needed for iteration and border-checks
00559     int w = srcShape[0], h = srcShape[1], d = srcShape[2];
00560     int x,y,z, i;
00561         
00562     //declare and define Iterators for all three dims at src
00563     SrcIterator zs = s_Iter;
00564     SrcIterator ys(zs);
00565     SrcIterator xs(ys);
00566         
00567     // temporary image to store region labels
00568     typedef vigra::MultiArray<3,int> LabelVolume;
00569     typedef LabelVolume::traverser LabelTraverser;
00570     
00571     LabelVolume labelvolume(srcShape);
00572         
00573     //Declare traversers for all three dims at target
00574     LabelTraverser zt = labelvolume.traverser_begin();
00575     LabelTraverser yt(zt);
00576     LabelTraverser xt(yt);
00577 
00578     // Kovalevsky's clever idea to use
00579     // image iterator and scan order iterator simultaneously
00580     // memory order indicates label
00581     LabelVolume::iterator label = labelvolume.begin();
00582     
00583     // initialize the neighborhood traversers
00584 
00585 #if 0
00586     NeighborOffsetTraverser<Neighborhood3D> nc(Neighborhood3D::CausalFirst);
00587     NeighborOffsetTraverser<Neighborhood3D> nce(Neighborhood3D::CausalLast);
00588 #endif /* #if 0 */
00589 
00590     NeighborOffsetCirculator<Neighborhood3D> nc(Neighborhood3D::CausalFirst);
00591     NeighborOffsetCirculator<Neighborhood3D> nce(Neighborhood3D::CausalLast);
00592     ++nce;
00593     // pass 1: scan image from upper left front to lower right back
00594     // to find connected components
00595 
00596     // Each component will be represented by a tree of pixels. Each
00597     // pixel contains the scan order address of its parent in the
00598     // tree.  In order for pass 2 to work correctly, the parent must
00599     // always have a smaller scan order address than the child.
00600     // Therefore, we can merge trees only at their roots, because the
00601     // root of the combined tree must have the smallest scan order
00602     // address among all the tree's pixels/ nodes.  The root of each
00603     // tree is distinguished by pointing to itself (it contains its
00604     // own scan order address). This condition is enforced whenever a
00605     // new region is found or two regions are merged
00606     i=0;
00607     int backgroundAddress=-1;
00608     for(z = 0; z != d; ++z, ++zs.dim2(), ++zt.dim2())
00609     {
00610         ys = zs;
00611         yt = zt;
00612 
00613         for(y = 0; y != h; ++y, ++ys.dim1(), ++yt.dim1())
00614         {
00615             xs = ys;
00616             xt = yt;
00617 
00618             for(x = 0; x != w; ++x, ++xs.dim0(), ++xt.dim0(), ++i)
00619             {
00620                 if(equal(*xs, backgroundValue)){
00621                         if(backgroundAddress==-1) backgroundAddress = i;
00622                         label[i]=backgroundAddress;
00623                 }
00624                 else{
00625                     *xt = i; // default: new region    
00626 
00627                     //queck whether there is a special borde threatment to be used or not
00628                     AtVolumeBorder atBorder = isAtVolumeBorderCausal(x,y,z,w,h,z);
00629                         
00630                     //We are not at the border!
00631                     if(atBorder == NotAtBorder)
00632                     {
00633                                 
00634 
00635 #if 0
00636                         nc = NeighborOffsetTraverser<Neighborhood3D>(Neighborhood3D::CausalFirst);
00637 #endif /* #if 0 */
00638 
00639                         nc = NeighborOffsetCirculator<Neighborhood3D>(Neighborhood3D::CausalFirst);
00640                     
00641                         do
00642                         {            
00643                             // if colors are equal
00644                             if(equal(*xs, xs[*nc]))
00645                             {
00646                                 int neighborLabel = xt[*nc];
00647 
00648                                 // find the root label of a label tree
00649                                 while(neighborLabel != label[neighborLabel])
00650                                 {
00651                                     neighborLabel = label[neighborLabel];
00652                                 }
00653 
00654                                 if(neighborLabel < *xt) // always keep the smallest among the possible neighbor labels
00655                                 {
00656                                     label[*xt] = neighborLabel;
00657                                     *xt = neighborLabel;
00658                                 }
00659                                 else
00660                                 {
00661                                     label[neighborLabel] = *xt;
00662                                 }
00663                             }
00664                             ++nc;
00665                         }while(nc!=nce);
00666                     }
00667                     //we are at a border - handle this!!
00668                     else
00669                     {
00670 
00671 #if 0
00672                         nc = NeighborOffsetTraverser<Neighborhood3D>(Neighborhood3D::nearBorderDirectionsCausal(atBorder,0));
00673 #endif /* #if 0 */
00674 
00675                         nc = NeighborOffsetCirculator<Neighborhood3D>(Neighborhood3D::nearBorderDirectionsCausal(atBorder,0));
00676                         int j=0;
00677                         while(nc.direction() != Neighborhood3D::Error)
00678                         {
00679                             //   colors equal???
00680                             if(equal(*xs, xs[*nc]))
00681                             {
00682                                 int neighborLabel = xt[*nc];
00683 
00684                                 // find the root label of a label tree
00685                                 while(neighborLabel != label[neighborLabel])
00686                                 {
00687                                     neighborLabel = label[neighborLabel];
00688                                 }
00689 
00690                                 if(neighborLabel < *xt) // always keep the smallest among the possible neighbor labels
00691                                 {
00692                                     label[*xt] = neighborLabel;
00693                                     *xt = neighborLabel;
00694                                 }
00695                                 else
00696                                 {
00697                                     label[neighborLabel] = *xt;
00698                                 }
00699                             }
00700                             nc.turnTo(Neighborhood3D::nearBorderDirectionsCausal(atBorder,++j));
00701                         }
00702                     }
00703                 }
00704             }
00705         }
00706     }
00707 
00708     // pass 2: assign one label to each region (tree)
00709     // so that labels form a consecutive sequence 1, 2, ...
00710     DestIterator zd = d_Iter;
00711 
00712     unsigned int count = 0; 
00713 
00714     i= 0;
00715 
00716     for(z=0; z != d; ++z, ++zd.dim2())
00717     {
00718         DestIterator yd(zd);
00719 
00720         for(y=0; y != h; ++y, ++yd.dim1())
00721         {
00722             DestIterator xd(yd);
00723 
00724             for(x = 0; x != w; ++x, ++xd.dim0(), ++i)
00725             {
00726                 if(label[i] == backgroundAddress){
00727                     label[i]=0;
00728                     continue;
00729                 }
00730                 
00731                 if(label[i] == i)
00732                 {
00733                     label[i] = count++;
00734                 }
00735                 else
00736                 {
00737                     label[i] = label[label[i]]; // compress trees
00738                 }
00739                 da.set(label[i]+1, xd);
00740             }
00741         }
00742     }
00743     return count;
00744 }
00745 
00746 //@}
00747 
00748 } //end of namespace vigra
00749 
00750 #endif //VIGRA_LABELVOLUME_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)