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