[ 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 /* 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_watersheds3D_HXX 00038 #define VIGRA_watersheds3D_HXX 00039 00040 #include "voxelneighborhood.hxx" 00041 #include "multi_array.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 #if 0 00085 NeighborhoodTraverser<SrcIterator, Neighborhood3D> c(xs), cend(c); 00086 #endif /* #if 0 */ 00087 00088 NeighborhoodCirculator<SrcIterator, Neighborhood3D> c(xs), cend(c); 00089 00090 do { 00091 if(sa(c) < v) 00092 { 00093 v = sa(c); 00094 o = c.directionBit(); 00095 } 00096 else if(sa(c) == my_v && my_v == v) 00097 { 00098 o = o | c.directionBit(); 00099 } 00100 } 00101 while(++c != cend); 00102 } 00103 else 00104 { 00105 #if 0 00106 RestrictedNeighborhoodTraverser<SrcIterator, Neighborhood3D> c(xs, atBorder), cend(c); 00107 #endif /* #if 0 */ 00108 00109 RestrictedNeighborhoodCirculator<SrcIterator, Neighborhood3D> c(xs, atBorder), cend(c); 00110 do { 00111 if(sa(c) < v) 00112 { 00113 v = sa(c); 00114 o = c.directionBit(); 00115 } 00116 else if(sa(c) == my_v && my_v == v) 00117 { 00118 o = o | c.directionBit(); 00119 } 00120 } 00121 while(++c != cend); 00122 } 00123 if (o==0) local_min_count++; 00124 da.set(o, xd); 00125 }//end x-iteration 00126 }//end y-iteration 00127 }//end z-iteration 00128 return local_min_count; 00129 } 00130 00131 template <class SrcIterator, class SrcAccessor,class SrcShape, 00132 class DestIterator, class DestAccessor, 00133 class Neighborhood3D> 00134 unsigned int watershedLabeling3D( SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa, 00135 DestIterator d_Iter, DestAccessor da, 00136 Neighborhood3D) 00137 { 00138 //basically needed for iteration and border-checks 00139 int w = srcShape[0], h = srcShape[1], d = srcShape[2]; 00140 int x,y,z, i; 00141 00142 //declare and define Iterators for all three dims at src 00143 SrcIterator zs = s_Iter; 00144 SrcIterator ys(zs); 00145 SrcIterator xs(ys); 00146 00147 // temporary image to store region labels 00148 typedef vigra::MultiArray<3,int> LabelVolume; 00149 typedef LabelVolume::traverser LabelTraverser; 00150 00151 LabelVolume labelvolume(srcShape); 00152 00153 //Declare traversers for all three dims at target 00154 LabelTraverser zt = labelvolume.traverser_begin(); 00155 LabelTraverser yt(zt); 00156 LabelTraverser xt(yt); 00157 00158 // Kovalevsky's clever idea to use 00159 // image iterator and scan order iterator simultaneously 00160 // memory order indicates label 00161 LabelVolume::iterator label = labelvolume.begin(); 00162 00163 // initialize the neighborhood traversers 00164 00165 #if 0 00166 NeighborOffsetTraverser<Neighborhood3D> nc(Neighborhood3D::CausalFirst); 00167 NeighborOffsetTraverser<Neighborhood3D> nce(Neighborhood3D::CausalLast); 00168 #endif /* #if 0 */ 00169 00170 NeighborOffsetCirculator<Neighborhood3D> nc(Neighborhood3D::CausalFirst); 00171 NeighborOffsetCirculator<Neighborhood3D> nce(Neighborhood3D::CausalLast); 00172 ++nce; 00173 // pass 1: scan image from upper left front to lower right back 00174 // to find connected components 00175 00176 // Each component will be represented by a tree of pixels. Each 00177 // pixel contains the scan order address of its parent in the 00178 // tree. In order for pass 2 to work correctly, the parent must 00179 // always have a smaller scan order address than the child. 00180 // Therefore, we can merge trees only at their roots, because the 00181 // root of the combined tree must have the smallest scan order 00182 // address among all the tree's pixels/ nodes. The root of each 00183 // tree is distinguished by pointing to itself (it contains its 00184 // own scan order address). This condition is enforced whenever a 00185 // new region is found or two regions are merged 00186 i=0; 00187 for(z = 0; z != d; ++z, ++zs.dim2(), ++zt.dim2()) 00188 { 00189 ys = zs; 00190 yt = zt; 00191 00192 for(y = 0; y != h; ++y, ++ys.dim1(), ++yt.dim1()) 00193 { 00194 xs = ys; 00195 xt = yt; 00196 00197 for(x = 0; x != w; ++x, ++xs.dim0(), ++xt.dim0(), ++i) 00198 { 00199 *xt = i; // default: new region 00200 00201 //queck whether there is a special borde threatment to be used or not 00202 AtVolumeBorder atBorder = isAtVolumeBorderCausal(x,y,z,w,h,z); 00203 00204 //We are not at the border! 00205 if(atBorder == NotAtBorder) 00206 { 00207 00208 #if 0 00209 nc = NeighborOffsetTraverser<Neighborhood3D>(Neighborhood3D::CausalFirst); 00210 #endif /* #if 0 */ 00211 00212 nc = NeighborOffsetCirculator<Neighborhood3D>(Neighborhood3D::CausalFirst); 00213 00214 do 00215 { 00216 // Direction of NTraversr Neighbor's direction bit is pointing 00217 // = Direction of voxel towards us? 00218 if((*xs & nc.directionBit()) || (xs[*nc] & nc.oppositeDirectionBit())) 00219 { 00220 int neighborLabel = xt[*nc]; 00221 00222 // find the root label of a label tree 00223 while(neighborLabel != label[neighborLabel]) 00224 { 00225 neighborLabel = label[neighborLabel]; 00226 } 00227 00228 if(neighborLabel < *xt) // always keep the smallest among the possible neighbor labels 00229 { 00230 label[*xt] = neighborLabel; 00231 *xt = neighborLabel; 00232 } 00233 else 00234 { 00235 label[neighborLabel] = *xt; 00236 } 00237 } 00238 ++nc; 00239 }while(nc!=nce); 00240 } 00241 //we are at a border - handle this!! 00242 else 00243 { 00244 00245 #if 0 00246 nc = NeighborOffsetTraverser<Neighborhood3D>(Neighborhood3D::nearBorderDirectionsCausal(atBorder,0)); 00247 #endif /* #if 0 */ 00248 00249 nc = NeighborOffsetCirculator<Neighborhood3D>(Neighborhood3D::nearBorderDirectionsCausal(atBorder,0)); 00250 int j=0; 00251 while(nc.direction() != Neighborhood3D::Error) 00252 { 00253 // Direction of NTraversr Neighbor's direction bit is pointing 00254 // = Direction of voxel towards us? 00255 if((*xs & nc.directionBit()) || (xs[*nc] & nc.oppositeDirectionBit())) 00256 { 00257 int neighborLabel = xt[*nc]; 00258 00259 // find the root label of a label tree 00260 while(neighborLabel != label[neighborLabel]) 00261 { 00262 neighborLabel = label[neighborLabel]; 00263 } 00264 00265 if(neighborLabel < *xt) // always keep the smallest among the possible neighbor labels 00266 { 00267 label[*xt] = neighborLabel; 00268 *xt = neighborLabel; 00269 } 00270 else 00271 { 00272 label[neighborLabel] = *xt; 00273 } 00274 } 00275 nc.turnTo(Neighborhood3D::nearBorderDirectionsCausal(atBorder,++j)); 00276 } 00277 } 00278 } 00279 } 00280 } 00281 00282 // pass 2: assign one label to each region (tree) 00283 // so that labels form a consecutive sequence 1, 2, ... 00284 DestIterator zd = d_Iter; 00285 00286 unsigned int count = 0; 00287 00288 i= 0; 00289 00290 for(z=0; z != d; ++z, ++zd.dim2()) 00291 { 00292 DestIterator yd(zd); 00293 00294 for(y=0; y != h; ++y, ++yd.dim1()) 00295 { 00296 DestIterator xd(yd); 00297 00298 for(x = 0; x != w; ++x, ++xd.dim0(), ++i) 00299 { 00300 if(label[i] == i) 00301 { 00302 label[i] = count++; 00303 } 00304 else 00305 { 00306 label[i] = label[label[i]]; // compress trees 00307 } 00308 da.set(label[i]+1, xd); 00309 } 00310 } 00311 } 00312 return count; 00313 } 00314 00315 00316 /** \addtogroup SeededRegionGrowing Region Segmentation Algorithms 00317 Region growing, watersheds, and voronoi tesselation 00318 */ 00319 //@{ 00320 00321 /********************************************************/ 00322 /* */ 00323 /* watersheds3D */ 00324 /* */ 00325 /********************************************************/ 00326 00327 /** \brief Region Segmentation by means of the watershed algorithm. 00328 00329 <b> Declarations:</b> 00330 00331 pass arguments explicitly: 00332 \code 00333 namespace vigra { 00334 template <class SrcIterator, class SrcAccessor,class SrcShape, 00335 class DestIterator, class DestAccessor, 00336 class Neighborhood3D> 00337 unsigned int watersheds3D(SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa, 00338 DestIterator d_Iter, DestAccessor da, 00339 Neighborhood3D neighborhood3D); 00340 } 00341 \endcode 00342 00343 use argument objects in conjunction with \ref ArgumentObjectFactories : 00344 \code 00345 namespace vigra { 00346 template <class SrcIterator, class SrcAccessor,class SrcShape, 00347 class DestIterator, class DestAccessor, 00348 class Neighborhood3D> 00349 unsigned int watersheds3D(triple<SrcIterator, SrcShape, SrcAccessor> src, 00350 pair<DestIterator, DestAccessor> dest, 00351 Neighborhood3D neighborhood3D); 00352 } 00353 \endcode 00354 00355 use with 3D-Six-Neighborhood: 00356 \code 00357 namespace vigra { 00358 00359 template <class SrcIterator, class SrcAccessor,class SrcShape, 00360 class DestIterator, class DestAccessor> 00361 unsigned int watersheds3DSix(triple<SrcIterator, SrcShape, SrcAccessor> src, 00362 pair<DestIterator, DestAccessor> dest); 00363 00364 } 00365 \endcode 00366 00367 use with 3D-TwentySix-Neighborhood: 00368 \code 00369 namespace vigra { 00370 00371 template <class SrcIterator, class SrcAccessor,class SrcShape, 00372 class DestIterator, class DestAccessor> 00373 unsigned int watersheds3DTwentySix(triple<SrcIterator, SrcShape, SrcAccessor> src, 00374 pair<DestIterator, DestAccessor> dest); 00375 00376 } 00377 \endcode 00378 00379 This function implements the union-find version of the watershed algorithms 00380 as described in 00381 00382 J. Roerdink, R. Meijster: "<em>The watershed transform: definitions, algorithms, 00383 and parallelization stretegies</em>", Fundamenta Informaticae, 41:187-228, 2000 00384 00385 The source volume is a boundary indicator such as the gradient magnitude 00386 of the trace of the \ref boundaryTensor(). Local minima of the boundary indicator 00387 are used as region seeds, and all other voxels are recursively assigned to the same 00388 region as their lowest neighbor. Pass \ref vigra::NeighborCode3DSix or 00389 \ref vigra::NeighborCode3DTwentySix to determine the neighborhood where voxel values 00390 are compared. The voxel type of the input volume must be <tt>LessThanComparable</tt>. 00391 The function uses accessors. 00392 00393 ...probably soon in VIGRA: 00394 Note that VIGRA provides an alternative implementaion of the watershed transform via 00395 \ref seededRegionGrowing3D(). It is slower, but handles plateaus better 00396 and allows to keep a one pixel wide boundary between regions. 00397 00398 <b> Usage:</b> 00399 00400 <b>\#include</b> <<a href="watersheds3D_8hxx-source.html">vigra/watersheds3D.hxx</a>><br> 00401 Namespace: vigra 00402 00403 Example: watersheds3D of the gradient magnitude. 00404 00405 \code 00406 typedef vigra::MultiArray<3,int> IntVolume; 00407 typedef vigra::MultiArray<3,double> DVolume; 00408 DVolume src(DVolume::difference_type(w,h,d)); 00409 IntVolume dest(IntVolume::difference_type(w,h,d)); 00410 00411 float gauss=1; 00412 00413 vigra::MultiArray<3, vigra::TinyVector<float,3> > temp(IntVolume::difference_type(w,h,d)); 00414 vigra::gaussianGradientMultiArray(srcMultiArrayRange(vol),destMultiArray(temp),gauss); 00415 00416 IntVolume::iterator temp_iter=temp.begin(); 00417 for(DVolume::iterator iter=src.begin(); iter!=src.end(); ++iter, ++temp_iter) 00418 *iter = norm(*temp_iter); 00419 00420 // find 6-connected regions 00421 int max_region_label = vigra::watersheds3DSix(srcMultiArrayRange(src), destMultiArray(dest)); 00422 00423 // find 26-connected regions 00424 max_region_label = vigra::watersheds3DTwentySix(srcMultiArrayRange(src), destMultiArray(dest)); 00425 00426 \endcode 00427 00428 <b> Required Interface:</b> 00429 00430 \code 00431 SrcIterator src_begin; 00432 SrcShape src_shape; 00433 DestIterator dest_begin; 00434 00435 SrcAccessor src_accessor; 00436 DestAccessor dest_accessor; 00437 00438 // compare src values 00439 src_accessor(src_begin) <= src_accessor(src_begin) 00440 00441 // set result 00442 int label; 00443 dest_accessor.set(label, dest_begin); 00444 \endcode 00445 */ 00446 doxygen_overloaded_function(template <...> unsigned int watersheds3D) 00447 00448 template <class SrcIterator, class SrcAccessor, class SrcShape, 00449 class DestIterator, class DestAccessor, 00450 class Neighborhood3D> 00451 unsigned int watersheds3D( SrcIterator s_Iter, SrcShape srcShape, SrcAccessor sa, 00452 DestIterator d_Iter, DestAccessor da, Neighborhood3D neighborhood3D) 00453 { 00454 //create temporary volume to store the DAG of directions to minima 00455 if ((int)Neighborhood3D::DirectionCount>7){ //If we have 3D-TwentySix Neighborhood 00456 00457 vigra::MultiArray<3,int> orientationVolume(srcShape); 00458 00459 int local_min_count= preparewatersheds3D( s_Iter, srcShape, sa, 00460 destMultiArray(orientationVolume).first, destMultiArray(orientationVolume).second, 00461 neighborhood3D); 00462 00463 return watershedLabeling3D( srcMultiArray(orientationVolume).first, srcShape, srcMultiArray(orientationVolume).second, 00464 d_Iter, da, 00465 neighborhood3D); 00466 } 00467 else{ 00468 00469 vigra::MultiArray<3,unsigned char> orientationVolume(srcShape); 00470 00471 int local_min_count= preparewatersheds3D( s_Iter, srcShape, sa, 00472 destMultiArray(orientationVolume).first, destMultiArray(orientationVolume).second, 00473 neighborhood3D); 00474 00475 return watershedLabeling3D( srcMultiArray(orientationVolume).first, srcShape, srcMultiArray(orientationVolume).second, 00476 d_Iter, da, 00477 neighborhood3D); 00478 } 00479 } 00480 00481 template <class SrcIterator, class SrcShape, class SrcAccessor, 00482 class DestIterator, class DestAccessor> 00483 inline unsigned int watersheds3DSix( vigra::triple<SrcIterator, SrcShape, SrcAccessor> src, 00484 vigra::pair<DestIterator, DestAccessor> dest) 00485 { 00486 return watersheds3D(src.first, src.second, src.third, dest.first, dest.second, NeighborCode3DSix()); 00487 } 00488 00489 template <class SrcIterator, class SrcShape, class SrcAccessor, 00490 class DestIterator, class DestAccessor> 00491 inline unsigned int watersheds3DTwentySix( vigra::triple<SrcIterator, SrcShape, SrcAccessor> src, 00492 vigra::pair<DestIterator, DestAccessor> dest) 00493 { 00494 return watersheds3D(src.first, src.second, src.third, dest.first, dest.second, NeighborCode3DTwentySix()); 00495 } 00496 00497 }//namespace vigra 00498 00499 #endif //VIGRA_watersheds3D_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|