[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
vigra/seededregiongrowing3d.hxx | ![]() |
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 Diff_type pos = voxel->location_; 00400 Diff_type nearest = voxel->nearest_; 00401 int lab = voxel->label_; 00402 CostType cost = voxel->cost_; 00403 00404 allocator.dismiss(voxel); 00405 00406 if((srgType & StopAtThreshold) != 0 && cost > max_cost) 00407 break; 00408 00409 irx = ir + pos; 00410 isx = srcul + pos; 00411 00412 if(*irx) // already labelled region / watershed? 00413 continue; 00414 00415 if((srgType & KeepContours) != 0) 00416 { 00417 for(int i=0; i<directionCount; i++) 00418 { 00419 cneighbor = * (irx + Neighborhood::diff((Direction)i)); 00420 if((cneighbor>0) && (cneighbor != lab)) 00421 { 00422 lab = SRGWatershedLabel; 00423 break; 00424 } 00425 } 00426 } 00427 00428 *irx = lab; 00429 00430 if((srgType & KeepContours) == 0 || lab > 0) 00431 { 00432 // update statistics 00433 stats[*irx](as(isx)); 00434 00435 // search neighborhood 00436 // second pass: find new candidate pixels 00437 for(int i=0; i<directionCount; i++) 00438 { 00439 if(*(irx + Neighborhood::diff((Direction)i)) == 0) 00440 { 00441 CostType cost = stats[lab].cost(as(isx, Neighborhood::diff((Direction)i))); 00442 00443 Voxel * new_voxel = 00444 allocator.create(pos+Neighborhood::diff((Direction)i), nearest, cost, count++, lab); 00445 pheap.push(new_voxel); 00446 } 00447 } 00448 } 00449 } 00450 00451 // free temporary memory 00452 while(pheap.size() != 0) 00453 { 00454 allocator.dismiss(pheap.top()); 00455 pheap.pop(); 00456 } 00457 00458 // write result 00459 transformMultiArray(ir, Diff_type(w,h,d), AccessorTraits<int>::default_accessor(), 00460 destul, ad, detail::UnlabelWatersheds()); 00461 } 00462 00463 template <class SrcImageIterator, class Diff_type, class SrcAccessor, 00464 class SeedImageIterator, class SeedAccessor, 00465 class DestImageIterator, class DestAccessor, 00466 class RegionStatisticsArray, class Neighborhood > 00467 inline void 00468 seededRegionGrowing3D(SrcImageIterator srcul, Diff_type shape, SrcAccessor as, 00469 SeedImageIterator seedsul, SeedAccessor aseeds, 00470 DestImageIterator destul, DestAccessor ad, 00471 RegionStatisticsArray & stats, SRGType srgType, Neighborhood n) 00472 { 00473 seededRegionGrowing3D( srcul, shape, as, seedsul, aseeds, 00474 destul, ad, stats, srgType, n, NumericTraits<double>::max()); 00475 } 00476 00477 template <class SrcImageIterator, class Diff_type, class SrcAccessor, 00478 class SeedImageIterator, class SeedAccessor, 00479 class DestImageIterator, class DestAccessor, 00480 class RegionStatisticsArray > 00481 inline void 00482 seededRegionGrowing3D(SrcImageIterator srcul, Diff_type shape, SrcAccessor as, 00483 SeedImageIterator seedsul, SeedAccessor aseeds, 00484 DestImageIterator destul, DestAccessor ad, 00485 RegionStatisticsArray & stats, SRGType srgType) 00486 { 00487 seededRegionGrowing3D( srcul, shape, as, seedsul, aseeds, 00488 destul, ad, stats, srgType, NeighborCode3DSix()); 00489 } 00490 00491 template <class SrcImageIterator, class Diff_type, class SrcAccessor, 00492 class SeedImageIterator, class SeedAccessor, 00493 class DestImageIterator, class DestAccessor, 00494 class RegionStatisticsArray > 00495 inline void 00496 seededRegionGrowing3D(SrcImageIterator srcul, Diff_type shape, SrcAccessor as, 00497 SeedImageIterator seedsul, SeedAccessor aseeds, 00498 DestImageIterator destul, DestAccessor ad, 00499 RegionStatisticsArray & stats) 00500 { 00501 seededRegionGrowing3D( srcul, shape, as, seedsul, aseeds, destul, ad, 00502 stats, CompleteGrow); 00503 } 00504 00505 template <class SrcImageIterator, class Shape, class SrcAccessor, 00506 class SeedImageIterator, class SeedAccessor, 00507 class DestImageIterator, class DestAccessor, 00508 class RegionStatisticsArray, class Neighborhood> 00509 inline void 00510 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1, 00511 pair<SeedImageIterator, SeedAccessor> img3, 00512 pair<DestImageIterator, DestAccessor> img4, 00513 RegionStatisticsArray & stats, 00514 SRGType srgType, Neighborhood n, double max_cost) 00515 { 00516 seededRegionGrowing3D(img1.first, img1.second, img1.third, 00517 img3.first, img3.second, 00518 img4.first, img4.second, 00519 stats, srgType, n, max_cost); 00520 } 00521 00522 template <class SrcImageIterator, class Shape, class SrcAccessor, 00523 class SeedImageIterator, class SeedAccessor, 00524 class DestImageIterator, class DestAccessor, 00525 class RegionStatisticsArray, class Neighborhood> 00526 inline void 00527 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1, 00528 pair<SeedImageIterator, SeedAccessor> img3, 00529 pair<DestImageIterator, DestAccessor> img4, 00530 RegionStatisticsArray & stats, 00531 SRGType srgType, Neighborhood n) 00532 { 00533 seededRegionGrowing3D(img1.first, img1.second, img1.third, 00534 img3.first, img3.second, 00535 img4.first, img4.second, 00536 stats, srgType, n, NumericTraits<double>::max()); 00537 } 00538 00539 template <class SrcImageIterator, class Shape, class SrcAccessor, 00540 class SeedImageIterator, class SeedAccessor, 00541 class DestImageIterator, class DestAccessor, 00542 class RegionStatisticsArray> 00543 inline void 00544 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1, 00545 pair<SeedImageIterator, SeedAccessor> img3, 00546 pair<DestImageIterator, DestAccessor> img4, 00547 RegionStatisticsArray & stats, SRGType srgType) 00548 { 00549 seededRegionGrowing3D(img1.first, img1.second, img1.third, 00550 img3.first, img3.second, 00551 img4.first, img4.second, 00552 stats, srgType, NeighborCode3DSix()); 00553 } 00554 00555 template <class SrcImageIterator, class Shape, class SrcAccessor, 00556 class SeedImageIterator, class SeedAccessor, 00557 class DestImageIterator, class DestAccessor, 00558 class RegionStatisticsArray> 00559 inline void 00560 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1, 00561 pair<SeedImageIterator, SeedAccessor> img3, 00562 pair<DestImageIterator, DestAccessor> img4, 00563 RegionStatisticsArray & stats) 00564 { 00565 seededRegionGrowing3D(img1.first, img1.second, img1.third, 00566 img3.first, img3.second, 00567 img4.first, img4.second, 00568 stats); 00569 } 00570 00571 } // namespace vigra 00572 00573 #endif // VIGRA_SEEDEDREGIONGROWING_HXX 00574
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|