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