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