[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 2003-2007 by Kasim Terzic, Christian-Dennis Rahn */ 00004 /* and Ullrich Koethe */ 00005 /* */ 00006 /* This file is part of the VIGRA computer vision library. */ 00007 /* The VIGRA Website is */ 00008 /* http://hci.iwr.uni-heidelberg.de/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_SEEDEDREGIONGROWING_3D_HXX 00038 #define VIGRA_SEEDEDREGIONGROWING_3D_HXX 00039 00040 #include <vector> 00041 #include <stack> 00042 #include <queue> 00043 #include "utilities.hxx" 00044 #include "stdimage.hxx" 00045 #include "stdimagefunctions.hxx" 00046 #include "seededregiongrowing.hxx" 00047 #include "multi_pointoperators.hxx" 00048 #include "voxelneighborhood.hxx" 00049 00050 namespace vigra { 00051 00052 namespace detail { 00053 00054 template <class COST, class Diff_type> 00055 class SeedRgVoxel 00056 { 00057 public: 00058 Diff_type location_, nearest_; 00059 COST cost_; 00060 int count_; 00061 int label_; 00062 int dist_; 00063 00064 SeedRgVoxel() 00065 //: location_(0,0,0), nearest_(0,0,0), cost_(0), count_(0), label_(0) 00066 { 00067 location_ = Diff_type(0,0,0); 00068 nearest_ = Diff_type(0,0,0); 00069 cost_ = 0; 00070 count_ = 0; 00071 label_ = 0; 00072 } 00073 00074 SeedRgVoxel(Diff_type const & location, Diff_type const & nearest, 00075 COST const & cost, int const & count, int const & label) 00076 : location_(location), nearest_(nearest), 00077 cost_(cost), count_(count), label_(label) 00078 { 00079 int dx = location_[0] - nearest_[0]; 00080 int dy = location_[1] - nearest_[1]; 00081 int dz = location_[2] - nearest_[2]; 00082 dist_ = dx * dx + dy * dy + dz * dz; 00083 } 00084 00085 void set(Diff_type const & location, Diff_type const & nearest, 00086 COST const & cost, int const & count, int const & label) 00087 { 00088 location_ = location; 00089 nearest_ = nearest; 00090 cost_ = cost; 00091 count_ = count; 00092 label_ = label; 00093 00094 int dx = location_[0] - nearest_[0]; 00095 int dy = location_[1] - nearest_[1]; 00096 int dz = location_[2] - nearest_[2]; 00097 dist_ = dx * dx + dy * dy + dz * dz; 00098 } 00099 00100 struct Compare 00101 { 00102 // must implement > since priority_queue looks for largest element 00103 bool operator()(SeedRgVoxel const & l, 00104 SeedRgVoxel const & r) const 00105 { 00106 if(r.cost_ == l.cost_) 00107 { 00108 if(r.dist_ == l.dist_) return r.count_ < l.count_; 00109 00110 return r.dist_ < l.dist_; 00111 } 00112 00113 return r.cost_ < l.cost_; 00114 } 00115 bool operator()(SeedRgVoxel const * l, 00116 SeedRgVoxel const * r) const 00117 { 00118 if(r->cost_ == l->cost_) 00119 { 00120 if(r->dist_ == l->dist_) return r->count_ < l->count_; 00121 00122 return r->dist_ < l->dist_; 00123 } 00124 00125 return r->cost_ < l->cost_; 00126 } 00127 }; 00128 00129 struct Allocator 00130 { 00131 ~Allocator() 00132 { 00133 while(!freelist_.empty()) 00134 { 00135 delete freelist_.top(); 00136 freelist_.pop(); 00137 } 00138 } 00139 00140 SeedRgVoxel * create(Diff_type const & location, Diff_type const & nearest, 00141 COST const & cost, int const & count, int const & label) 00142 { 00143 if(!freelist_.empty()) 00144 { 00145 SeedRgVoxel * res = freelist_.top(); 00146 freelist_.pop(); 00147 res->set(location, nearest, cost, count, label); 00148 return res; 00149 } 00150 00151 return new SeedRgVoxel(location, nearest, cost, count, label); 00152 } 00153 00154 void dismiss(SeedRgVoxel * p) 00155 { 00156 freelist_.push(p); 00157 } 00158 00159 std::stack<SeedRgVoxel<COST,Diff_type> *> freelist_; 00160 }; 00161 }; 00162 00163 } // namespace detail 00164 00165 /** \addtogroup SeededRegionGrowing 00166 Region segmentation and voronoi tesselation 00167 */ 00168 //@{ 00169 00170 /********************************************************/ 00171 /* */ 00172 /* seededRegionGrowing3D */ 00173 /* */ 00174 /********************************************************/ 00175 00176 /** \brief Three-dimensional Region Segmentation by means of Seeded Region Growing. 00177 00178 This algorithm implements seeded region growing as described in 00179 00180 The seed image is a partly segmented multi-dimensional array which contains uniquely 00181 labeled regions (the seeds) and unlabeled voxels (the candidates, label 0). 00182 Seed regions can be as large as you wish and as small as one voxel. If 00183 there are no candidates, the algorithm will simply copy the seed array 00184 into the output array. Otherwise it will aggregate the candidates into 00185 the existing regions so that a cost function is minimized. 00186 Candidates are taken from the neighborhood of the already assigned pixels, 00187 where the type of neighborhood is determined by parameter <tt>neighborhood</tt> 00188 which can take the values <tt>NeighborCode3DSix()</tt> (the default) 00189 or <tt>NeighborCode3DTwentySix()</tt>. The algorithm basically works as follows 00190 (illustrated for 6-neighborhood, but 26-neighborhood works in the same way): 00191 00192 <ol> 00193 00194 <li> Find all candidate pixels that are 6-adjacent to a seed region. 00195 Calculate the cost for aggregating each candidate into its adajacent region 00196 and put the candidates into a priority queue. 00197 00198 <li> While( priority queue is not empty) 00199 00200 <ol> 00201 00202 <li> Take the candidate with least cost from the queue. If it has not 00203 already been merged, merge it with it's adjacent region. 00204 00205 <li> Put all candidates that are 4-adjacent to the pixel just processed 00206 into the priority queue. 00207 00208 </ol> 00209 00210 </ol> 00211 00212 <tt>SRGType</tt> can take the following values: 00213 00214 <DL> 00215 <DT><tt>CompleteGrow</tt> <DD> produce a complete tesselation of the volume (default). 00216 <DT><tt>KeepContours</tt> <DD> keep a 1-voxel wide unlabeled contour between all regions. 00217 <DT><tt>StopAtThreshold</tt> <DD> stop when the boundary indicator values exceed the 00218 threshold given by parameter <tt>max_cost</tt>. 00219 <DT><tt>KeepContours | StopAtThreshold</tt> <DD> keep 1-voxel wide contour and stop at given <tt>max_cost</tt>. 00220 </DL> 00221 00222 The cost is determined jointly by the source array and the 00223 region statistics functor. The source array contains feature values for each 00224 pixel which will be used by the region statistics functor to calculate and 00225 update statistics for each region and to calculate the cost for each 00226 candidate. The <TT>RegionStatisticsArray</TT> must be compatible to the 00227 \ref ArrayOfRegionStatistics functor and contains an <em> array</em> of 00228 statistics objects for each region. The indices must correspond to the 00229 labels of the seed regions. The statistics for the initial regions must have 00230 been calculated prior to calling <TT>seededRegionGrowing3D()</TT> 00231 00232 For each candidate 00233 <TT>x</TT> that is adjacent to region <TT>i</TT>, the algorithm will call 00234 <TT>stats[i].cost(as(x))</TT> to get the cost (where <TT>x</TT> is a <TT>SrcImageIterator</TT> 00235 and <TT>as</TT> is 00236 the SrcAccessor). When a candidate has been merged with a region, the 00237 statistics are updated by calling <TT>stats[i].operator()(as(x))</TT>. Since 00238 the <TT>RegionStatisticsArray</TT> is passed by reference, this will overwrite 00239 the original statistics. 00240 00241 If a candidate could be merged into more than one regions with identical 00242 cost, the algorithm will favour the nearest region. If <tt>StopAtThreshold</tt> is active, 00243 and the cost of the current candidate at any point in the algorithm exceeds the optional 00244 <tt>max_cost</tt> value (which defaults to <tt>NumericTraits<double>::max()</tt>), 00245 region growing is aborted, and all voxels not yet assigned to a region remain unlabeled. 00246 00247 In some cases, the cost only depends on the feature value of the current 00248 voxel. Then the update operation will simply be a no-op, and the <TT>cost()</TT> 00249 function returns its argument. This behavior is implemented by the 00250 \ref SeedRgDirectValueFunctor. 00251 00252 <b> Declarations:</b> 00253 00254 pass arguments explicitly: 00255 \code 00256 namespace vigra { 00257 template <class SrcImageIterator, class Shape, class SrcAccessor, 00258 class SeedImageIterator, class SeedAccessor, 00259 class DestImageIterator, class DestAccessor, 00260 class RegionStatisticsArray, class Neighborhood> 00261 void 00262 seededRegionGrowing3D(SrcImageIterator srcul, Shape shape, SrcAccessor as, 00263 SeedImageIterator seedsul, SeedAccessor aseeds, 00264 DestImageIterator destul, DestAccessor ad, 00265 RegionStatisticsArray & stats, 00266 SRGType srgType = CompleteGrow, 00267 Neighborhood neighborhood = NeighborCode3DSix(), 00268 double max_cost = NumericTraits<double>::max()); 00269 } 00270 \endcode 00271 00272 use argument objects in conjunction with \ref ArgumentObjectFactories : 00273 \code 00274 namespace vigra { 00275 template <class SrcImageIterator, class Shape, class SrcAccessor, 00276 class SeedImageIterator, class SeedAccessor, 00277 class DestImageIterator, class DestAccessor, 00278 class RegionStatisticsArray, class Neighborhood> 00279 void 00280 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> src, 00281 pair<SeedImageIterator, SeedAccessor> seeds, 00282 pair<DestImageIterator, DestAccessor> dest, 00283 RegionStatisticsArray & stats, 00284 SRGType srgType = CompleteGrow, 00285 Neighborhood neighborhood = NeighborCode3DSix(), 00286 double max_cost = NumericTraits<double>::max()); 00287 } 00288 \endcode 00289 00290 */ 00291 doxygen_overloaded_function(template <...> void seededRegionGrowing3D) 00292 00293 template <class SrcImageIterator, class Diff_type, class SrcAccessor, 00294 class SeedImageIterator, class SeedAccessor, 00295 class DestImageIterator, class DestAccessor, 00296 class RegionStatisticsArray, class Neighborhood> 00297 void 00298 seededRegionGrowing3D(SrcImageIterator srcul, Diff_type shape, SrcAccessor as, 00299 SeedImageIterator seedsul, SeedAccessor aseeds, 00300 DestImageIterator destul, DestAccessor ad, 00301 RegionStatisticsArray & stats, 00302 SRGType srgType, 00303 Neighborhood, 00304 double max_cost) 00305 { 00306 SrcImageIterator srclr = srcul + shape; 00307 //int w = srclr.x - srcul.x; 00308 int w = shape[0]; 00309 //int h = srclr.y - srcul.y; 00310 int h = shape[1]; 00311 //int d = srclr.z - srcul.z; 00312 int d = shape[2]; 00313 int count = 0; 00314 00315 SrcImageIterator isy = srcul, isx = srcul, isz = srcul; // iterators for the src image 00316 00317 typedef typename RegionStatisticsArray::value_type RegionStatistics; 00318 typedef typename PromoteTraits<typename RegionStatistics::cost_type, double>::Promote CostType; 00319 typedef detail::SeedRgVoxel<CostType, Diff_type> Voxel; 00320 00321 typename Voxel::Allocator allocator; 00322 00323 typedef std::priority_queue< Voxel *, 00324 std::vector<Voxel *>, 00325 typename Voxel::Compare > SeedRgVoxelHeap; 00326 typedef MultiArray<3, int> IVolume; 00327 00328 // copy seed image in an image with border 00329 Diff_type regionshape = shape + Diff_type(2,2,2); 00330 IVolume regions(regionshape); 00331 MultiIterator<3,int> ir = regions.traverser_begin(); 00332 ir = ir + Diff_type(1,1,1); 00333 00334 //IVolume::Iterator iry, irx, irz; 00335 MultiIterator<3,int> iry, irx, irz; 00336 00337 //initImageBorder(destImageRange(regions), 1, SRGWatershedLabel); 00338 initMultiArrayBorder(destMultiArrayRange(regions), 1, SRGWatershedLabel); 00339 00340 copyMultiArray(seedsul, Diff_type(w,h,d), aseeds, ir, AccessorTraits<int>::default_accessor()); 00341 00342 // allocate and init memory for the results 00343 00344 SeedRgVoxelHeap pheap; 00345 int cneighbor; 00346 00347 #if 0 00348 static const Diff_type dist[] = { Diff_type(-1, 0, 0), Diff_type( 0,-1, 0), 00349 Diff_type( 1, 0, 0), Diff_type( 0, 1, 0), 00350 Diff_type( 0, 0,-1), Diff_type( 0, 0, 1) }; 00351 #endif 00352 00353 typedef typename Neighborhood::Direction Direction; 00354 int directionCount = Neighborhood::DirectionCount; 00355 00356 Diff_type pos(0,0,0); 00357 00358 for(isz=srcul, irz=ir, pos[2]=0; pos[2]<d; 00359 pos[2]++, isz.dim2()++, irz.dim2()++) 00360 { 00361 //std::cerr << "Z = " << pos[2] << std::endl; 00362 00363 for(isy=isz, iry=irz, pos[1]=0; pos[1]<h; 00364 pos[1]++, isy.dim1()++, iry.dim1()++) 00365 { 00366 //std::cerr << "Y = " << pos[1] << std::endl; 00367 00368 for(isx=isy, irx=iry, pos[0]=0; pos[0]<w; 00369 pos[0]++, isx.dim0()++, irx.dim0()++) 00370 { 00371 //std::cerr << "X = " << pos[0] << std::endl; 00372 00373 if(*irx == 0) 00374 { 00375 // find candidate pixels for growing and fill heap 00376 for(int i=0; i<directionCount; i++) 00377 { 00378 cneighbor = *(irx + Neighborhood::diff((Direction)i)); 00379 if(cneighbor > 0) 00380 { 00381 CostType cost = stats[cneighbor].cost(as(isx)); 00382 00383 Voxel * voxel = 00384 allocator.create(pos, pos+Neighborhood::diff((Direction)i), cost, count++, cneighbor); 00385 pheap.push(voxel); 00386 } 00387 } 00388 } 00389 } 00390 } 00391 } 00392 00393 // perform region growing 00394 while(pheap.size() != 0) 00395 { 00396 Voxel * voxel = pheap.top(); 00397 pheap.pop(); 00398 00399 if((srgType & StopAtThreshold) != 0 && voxel->cost_ > max_cost) 00400 break; 00401 00402 Diff_type pos = voxel->location_; 00403 Diff_type nearest = voxel->nearest_; 00404 int lab = voxel->label_; 00405 00406 allocator.dismiss(voxel); 00407 00408 irx = ir + pos; 00409 isx = srcul + pos; 00410 00411 if(*irx) // already labelled region / watershed? 00412 continue; 00413 00414 if((srgType & KeepContours) != 0) 00415 { 00416 for(int i=0; i<directionCount; i++) 00417 { 00418 cneighbor = * (irx + Neighborhood::diff((Direction)i)); 00419 if((cneighbor>0) && (cneighbor != lab)) 00420 { 00421 lab = SRGWatershedLabel; 00422 break; 00423 } 00424 } 00425 } 00426 00427 *irx = lab; 00428 00429 if((srgType & KeepContours) == 0 || lab > 0) 00430 { 00431 // update statistics 00432 stats[*irx](as(isx)); 00433 00434 // search neighborhood 00435 // second pass: find new candidate pixels 00436 for(int i=0; i<directionCount; i++) 00437 { 00438 if(*(irx + Neighborhood::diff((Direction)i)) == 0) 00439 { 00440 CostType cost = stats[lab].cost(as(isx, Neighborhood::diff((Direction)i))); 00441 00442 Voxel * new_voxel = 00443 allocator.create(pos+Neighborhood::diff((Direction)i), nearest, cost, count++, lab); 00444 pheap.push(new_voxel); 00445 } 00446 } 00447 } 00448 } 00449 00450 // write result 00451 transformMultiArray(ir, Diff_type(w,h,d), AccessorTraits<int>::default_accessor(), 00452 destul, ad, detail::UnlabelWatersheds()); 00453 } 00454 00455 template <class SrcImageIterator, class Diff_type, class SrcAccessor, 00456 class SeedImageIterator, class SeedAccessor, 00457 class DestImageIterator, class DestAccessor, 00458 class RegionStatisticsArray, class Neighborhood > 00459 inline void 00460 seededRegionGrowing3D(SrcImageIterator srcul, Diff_type shape, SrcAccessor as, 00461 SeedImageIterator seedsul, SeedAccessor aseeds, 00462 DestImageIterator destul, DestAccessor ad, 00463 RegionStatisticsArray & stats, SRGType srgType, Neighborhood n) 00464 { 00465 seededRegionGrowing3D( srcul, shape, as, seedsul, aseeds, 00466 destul, ad, stats, srgType, n, NumericTraits<double>::max()); 00467 } 00468 00469 template <class SrcImageIterator, class Diff_type, class SrcAccessor, 00470 class SeedImageIterator, class SeedAccessor, 00471 class DestImageIterator, class DestAccessor, 00472 class RegionStatisticsArray > 00473 inline void 00474 seededRegionGrowing3D(SrcImageIterator srcul, Diff_type shape, SrcAccessor as, 00475 SeedImageIterator seedsul, SeedAccessor aseeds, 00476 DestImageIterator destul, DestAccessor ad, 00477 RegionStatisticsArray & stats, SRGType srgType) 00478 { 00479 seededRegionGrowing3D( srcul, shape, as, seedsul, aseeds, 00480 destul, ad, stats, srgType, NeighborCode3DSix()); 00481 } 00482 00483 template <class SrcImageIterator, class Diff_type, class SrcAccessor, 00484 class SeedImageIterator, class SeedAccessor, 00485 class DestImageIterator, class DestAccessor, 00486 class RegionStatisticsArray > 00487 inline void 00488 seededRegionGrowing3D(SrcImageIterator srcul, Diff_type shape, SrcAccessor as, 00489 SeedImageIterator seedsul, SeedAccessor aseeds, 00490 DestImageIterator destul, DestAccessor ad, 00491 RegionStatisticsArray & stats) 00492 { 00493 seededRegionGrowing3D( srcul, shape, as, seedsul, aseeds, destul, ad, 00494 stats, CompleteGrow); 00495 } 00496 00497 template <class SrcImageIterator, class Shape, class SrcAccessor, 00498 class SeedImageIterator, class SeedAccessor, 00499 class DestImageIterator, class DestAccessor, 00500 class RegionStatisticsArray, class Neighborhood> 00501 inline void 00502 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1, 00503 pair<SeedImageIterator, SeedAccessor> img3, 00504 pair<DestImageIterator, DestAccessor> img4, 00505 RegionStatisticsArray & stats, 00506 SRGType srgType, Neighborhood n, double max_cost) 00507 { 00508 seededRegionGrowing3D(img1.first, img1.second, img1.third, 00509 img3.first, img3.second, 00510 img4.first, img4.second, 00511 stats, srgType, n, max_cost); 00512 } 00513 00514 template <class SrcImageIterator, class Shape, class SrcAccessor, 00515 class SeedImageIterator, class SeedAccessor, 00516 class DestImageIterator, class DestAccessor, 00517 class RegionStatisticsArray, class Neighborhood> 00518 inline void 00519 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1, 00520 pair<SeedImageIterator, SeedAccessor> img3, 00521 pair<DestImageIterator, DestAccessor> img4, 00522 RegionStatisticsArray & stats, 00523 SRGType srgType, Neighborhood n) 00524 { 00525 seededRegionGrowing3D(img1.first, img1.second, img1.third, 00526 img3.first, img3.second, 00527 img4.first, img4.second, 00528 stats, srgType, n, NumericTraits<double>::max()); 00529 } 00530 00531 template <class SrcImageIterator, class Shape, class SrcAccessor, 00532 class SeedImageIterator, class SeedAccessor, 00533 class DestImageIterator, class DestAccessor, 00534 class RegionStatisticsArray> 00535 inline void 00536 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1, 00537 pair<SeedImageIterator, SeedAccessor> img3, 00538 pair<DestImageIterator, DestAccessor> img4, 00539 RegionStatisticsArray & stats, SRGType srgType) 00540 { 00541 seededRegionGrowing3D(img1.first, img1.second, img1.third, 00542 img3.first, img3.second, 00543 img4.first, img4.second, 00544 stats, srgType, NeighborCode3DSix()); 00545 } 00546 00547 template <class SrcImageIterator, class Shape, class SrcAccessor, 00548 class SeedImageIterator, class SeedAccessor, 00549 class DestImageIterator, class DestAccessor, 00550 class RegionStatisticsArray> 00551 inline void 00552 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1, 00553 pair<SeedImageIterator, SeedAccessor> img3, 00554 pair<DestImageIterator, DestAccessor> img4, 00555 RegionStatisticsArray & stats) 00556 { 00557 seededRegionGrowing3D(img1.first, img1.second, img1.third, 00558 img3.first, img3.second, 00559 img4.first, img4.second, 00560 stats); 00561 } 00562 00563 } // namespace vigra 00564 00565 #endif // VIGRA_SEEDEDREGIONGROWING_HXX 00566
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|