[ 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 */ 00167 //@{ 00168 00169 /********************************************************/ 00170 /* */ 00171 /* seededRegionGrowing3D */ 00172 /* */ 00173 /********************************************************/ 00174 00175 /** \brief Three-dimensional Region Segmentation by means of Seeded Region Growing. 00176 00177 This algorithm implements seeded region growing as described in 00178 00179 The seed image is a partly segmented multi-dimensional array which contains uniquely 00180 labeled regions (the seeds) and unlabeled voxels (the candidates, label 0). 00181 Seed regions can be as large as you wish and as small as one voxel. If 00182 there are no candidates, the algorithm will simply copy the seed array 00183 into the output array. Otherwise it will aggregate the candidates into 00184 the existing regions so that a cost function is minimized. 00185 Candidates are taken from the neighborhood of the already assigned pixels, 00186 where the type of neighborhood is determined by parameter <tt>neighborhood</tt> 00187 which can take the values <tt>NeighborCode3DSix()</tt> (the default) 00188 or <tt>NeighborCode3DTwentySix()</tt>. The algorithm basically works as follows 00189 (illustrated for 6-neighborhood, but 26-neighborhood works in the same way): 00190 00191 <ol> 00192 00193 <li> Find all candidate pixels that are 6-adjacent to a seed region. 00194 Calculate the cost for aggregating each candidate into its adjacent region 00195 and put the candidates into a priority queue. 00196 00197 <li> While( priority queue is not empty) 00198 00199 <ol> 00200 00201 <li> Take the candidate with least cost from the queue. If it has not 00202 already been merged, merge it with it's adjacent region. 00203 00204 <li> Put all candidates that are 4-adjacent to the pixel just processed 00205 into the priority queue. 00206 00207 </ol> 00208 00209 </ol> 00210 00211 <tt>SRGType</tt> can take the following values: 00212 00213 <DL> 00214 <DT><tt>CompleteGrow</tt> <DD> produce a complete tesselation of the volume (default). 00215 <DT><tt>KeepContours</tt> <DD> keep a 1-voxel wide unlabeled contour between all regions. 00216 <DT><tt>StopAtThreshold</tt> <DD> stop when the boundary indicator values exceed the 00217 threshold given by parameter <tt>max_cost</tt>. 00218 <DT><tt>KeepContours | StopAtThreshold</tt> <DD> keep 1-voxel wide contour and stop at given <tt>max_cost</tt>. 00219 </DL> 00220 00221 The cost is determined jointly by the source array and the 00222 region statistics functor. The source array contains feature values for each 00223 pixel which will be used by the region statistics functor to calculate and 00224 update statistics for each region and to calculate the cost for each 00225 candidate. The <TT>RegionStatisticsArray</TT> must be compatible to the 00226 \ref ArrayOfRegionStatistics functor and contains an <em> array</em> of 00227 statistics objects for each region. The indices must correspond to the 00228 labels of the seed regions. The statistics for the initial regions must have 00229 been calculated prior to calling <TT>seededRegionGrowing3D()</TT> 00230 00231 For each candidate 00232 <TT>x</TT> that is adjacent to region <TT>i</TT>, the algorithm will call 00233 <TT>stats[i].cost(as(x))</TT> to get the cost (where <TT>x</TT> is a <TT>SrcImageIterator</TT> 00234 and <TT>as</TT> is 00235 the SrcAccessor). When a candidate has been merged with a region, the 00236 statistics are updated by calling <TT>stats[i].operator()(as(x))</TT>. Since 00237 the <TT>RegionStatisticsArray</TT> is passed by reference, this will overwrite 00238 the original statistics. 00239 00240 If a candidate could be merged into more than one regions with identical 00241 cost, the algorithm will favour the nearest region. If <tt>StopAtThreshold</tt> is active, 00242 and the cost of the current candidate at any point in the algorithm exceeds the optional 00243 <tt>max_cost</tt> value (which defaults to <tt>NumericTraits<double>::max()</tt>), 00244 region growing is aborted, and all voxels not yet assigned to a region remain unlabeled. 00245 00246 In some cases, the cost only depends on the feature value of the current 00247 voxel. Then the update operation will simply be a no-op, and the <TT>cost()</TT> 00248 function returns its argument. This behavior is implemented by the 00249 \ref SeedRgDirectValueFunctor. 00250 00251 <b> Declarations:</b> 00252 00253 pass arguments explicitly: 00254 \code 00255 namespace vigra { 00256 template <class SrcImageIterator, class Shape, class SrcAccessor, 00257 class SeedImageIterator, class SeedAccessor, 00258 class DestImageIterator, class DestAccessor, 00259 class RegionStatisticsArray, class Neighborhood> 00260 void 00261 seededRegionGrowing3D(SrcImageIterator srcul, Shape shape, SrcAccessor as, 00262 SeedImageIterator seedsul, SeedAccessor aseeds, 00263 DestImageIterator destul, DestAccessor ad, 00264 RegionStatisticsArray & stats, 00265 SRGType srgType = CompleteGrow, 00266 Neighborhood neighborhood = NeighborCode3DSix(), 00267 double max_cost = NumericTraits<double>::max()); 00268 } 00269 \endcode 00270 00271 use argument objects in conjunction with \ref ArgumentObjectFactories : 00272 \code 00273 namespace vigra { 00274 template <class SrcImageIterator, class Shape, class SrcAccessor, 00275 class SeedImageIterator, class SeedAccessor, 00276 class DestImageIterator, class DestAccessor, 00277 class RegionStatisticsArray, class Neighborhood> 00278 void 00279 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> src, 00280 pair<SeedImageIterator, SeedAccessor> seeds, 00281 pair<DestImageIterator, DestAccessor> dest, 00282 RegionStatisticsArray & stats, 00283 SRGType srgType = CompleteGrow, 00284 Neighborhood neighborhood = NeighborCode3DSix(), 00285 double max_cost = NumericTraits<double>::max()); 00286 } 00287 \endcode 00288 00289 */ 00290 doxygen_overloaded_function(template <...> void seededRegionGrowing3D) 00291 00292 template <class SrcImageIterator, class Diff_type, class SrcAccessor, 00293 class SeedImageIterator, class SeedAccessor, 00294 class DestImageIterator, class DestAccessor, 00295 class RegionStatisticsArray, class Neighborhood> 00296 void 00297 seededRegionGrowing3D(SrcImageIterator srcul, Diff_type shape, SrcAccessor as, 00298 SeedImageIterator seedsul, SeedAccessor aseeds, 00299 DestImageIterator destul, DestAccessor ad, 00300 RegionStatisticsArray & stats, 00301 SRGType srgType, 00302 Neighborhood, 00303 double max_cost) 00304 { 00305 SrcImageIterator srclr = srcul + shape; 00306 //int w = srclr.x - srcul.x; 00307 int w = shape[0]; 00308 //int h = srclr.y - srcul.y; 00309 int h = shape[1]; 00310 //int d = srclr.z - srcul.z; 00311 int d = shape[2]; 00312 int count = 0; 00313 00314 SrcImageIterator isy = srcul, isx = srcul, isz = srcul; // iterators for the src image 00315 00316 typedef typename RegionStatisticsArray::value_type RegionStatistics; 00317 typedef typename PromoteTraits<typename RegionStatistics::cost_type, double>::Promote CostType; 00318 typedef detail::SeedRgVoxel<CostType, Diff_type> Voxel; 00319 00320 typename Voxel::Allocator allocator; 00321 00322 typedef std::priority_queue< Voxel *, 00323 std::vector<Voxel *>, 00324 typename Voxel::Compare > SeedRgVoxelHeap; 00325 typedef MultiArray<3, int> IVolume; 00326 00327 // copy seed image in an image with border 00328 Diff_type regionshape = shape + Diff_type(2,2,2); 00329 IVolume regions(regionshape); 00330 MultiIterator<3,int> ir = regions.traverser_begin(); 00331 ir = ir + Diff_type(1,1,1); 00332 00333 //IVolume::Iterator iry, irx, irz; 00334 MultiIterator<3,int> iry, irx, irz; 00335 00336 //initImageBorder(destImageRange(regions), 1, SRGWatershedLabel); 00337 initMultiArrayBorder(destMultiArrayRange(regions), 1, SRGWatershedLabel); 00338 00339 copyMultiArray(seedsul, Diff_type(w,h,d), aseeds, ir, AccessorTraits<int>::default_accessor()); 00340 00341 // allocate and init memory for the results 00342 00343 SeedRgVoxelHeap pheap; 00344 int cneighbor; 00345 00346 #if 0 00347 static const Diff_type dist[] = { Diff_type(-1, 0, 0), Diff_type( 0,-1, 0), 00348 Diff_type( 1, 0, 0), Diff_type( 0, 1, 0), 00349 Diff_type( 0, 0,-1), Diff_type( 0, 0, 1) }; 00350 #endif 00351 00352 typedef typename Neighborhood::Direction Direction; 00353 int directionCount = Neighborhood::DirectionCount; 00354 00355 Diff_type pos(0,0,0); 00356 00357 for(isz=srcul, irz=ir, pos[2]=0; pos[2]<d; 00358 pos[2]++, isz.dim2()++, irz.dim2()++) 00359 { 00360 //std::cerr << "Z = " << pos[2] << std::endl; 00361 00362 for(isy=isz, iry=irz, pos[1]=0; pos[1]<h; 00363 pos[1]++, isy.dim1()++, iry.dim1()++) 00364 { 00365 //std::cerr << "Y = " << pos[1] << std::endl; 00366 00367 for(isx=isy, irx=iry, pos[0]=0; pos[0]<w; 00368 pos[0]++, isx.dim0()++, irx.dim0()++) 00369 { 00370 //std::cerr << "X = " << pos[0] << std::endl; 00371 00372 if(*irx == 0) 00373 { 00374 // find candidate pixels for growing and fill heap 00375 for(int i=0; i<directionCount; i++) 00376 { 00377 cneighbor = *(irx + Neighborhood::diff((Direction)i)); 00378 if(cneighbor > 0) 00379 { 00380 CostType cost = stats[cneighbor].cost(as(isx)); 00381 00382 Voxel * voxel = 00383 allocator.create(pos, pos+Neighborhood::diff((Direction)i), cost, count++, cneighbor); 00384 pheap.push(voxel); 00385 } 00386 } 00387 } 00388 } 00389 } 00390 } 00391 00392 // perform region growing 00393 while(pheap.size() != 0) 00394 { 00395 Voxel * voxel = pheap.top(); 00396 pheap.pop(); 00397 00398 Diff_type pos = voxel->location_; 00399 Diff_type nearest = voxel->nearest_; 00400 int lab = voxel->label_; 00401 CostType cost = voxel->cost_; 00402 00403 allocator.dismiss(voxel); 00404 00405 if((srgType & StopAtThreshold) != 0 && cost > max_cost) 00406 break; 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 // free temporary memory 00451 while(pheap.size() != 0) 00452 { 00453 allocator.dismiss(pheap.top()); 00454 pheap.pop(); 00455 } 00456 00457 // write result 00458 transformMultiArray(ir, Diff_type(w,h,d), AccessorTraits<int>::default_accessor(), 00459 destul, ad, detail::UnlabelWatersheds()); 00460 } 00461 00462 template <class SrcImageIterator, class Diff_type, class SrcAccessor, 00463 class SeedImageIterator, class SeedAccessor, 00464 class DestImageIterator, class DestAccessor, 00465 class RegionStatisticsArray, class Neighborhood > 00466 inline void 00467 seededRegionGrowing3D(SrcImageIterator srcul, Diff_type shape, SrcAccessor as, 00468 SeedImageIterator seedsul, SeedAccessor aseeds, 00469 DestImageIterator destul, DestAccessor ad, 00470 RegionStatisticsArray & stats, SRGType srgType, Neighborhood n) 00471 { 00472 seededRegionGrowing3D( srcul, shape, as, seedsul, aseeds, 00473 destul, ad, stats, srgType, n, NumericTraits<double>::max()); 00474 } 00475 00476 template <class SrcImageIterator, class Diff_type, class SrcAccessor, 00477 class SeedImageIterator, class SeedAccessor, 00478 class DestImageIterator, class DestAccessor, 00479 class RegionStatisticsArray > 00480 inline void 00481 seededRegionGrowing3D(SrcImageIterator srcul, Diff_type shape, SrcAccessor as, 00482 SeedImageIterator seedsul, SeedAccessor aseeds, 00483 DestImageIterator destul, DestAccessor ad, 00484 RegionStatisticsArray & stats, SRGType srgType) 00485 { 00486 seededRegionGrowing3D( srcul, shape, as, seedsul, aseeds, 00487 destul, ad, stats, srgType, NeighborCode3DSix()); 00488 } 00489 00490 template <class SrcImageIterator, class Diff_type, class SrcAccessor, 00491 class SeedImageIterator, class SeedAccessor, 00492 class DestImageIterator, class DestAccessor, 00493 class RegionStatisticsArray > 00494 inline void 00495 seededRegionGrowing3D(SrcImageIterator srcul, Diff_type shape, SrcAccessor as, 00496 SeedImageIterator seedsul, SeedAccessor aseeds, 00497 DestImageIterator destul, DestAccessor ad, 00498 RegionStatisticsArray & stats) 00499 { 00500 seededRegionGrowing3D( srcul, shape, as, seedsul, aseeds, destul, ad, 00501 stats, CompleteGrow); 00502 } 00503 00504 template <class SrcImageIterator, class Shape, class SrcAccessor, 00505 class SeedImageIterator, class SeedAccessor, 00506 class DestImageIterator, class DestAccessor, 00507 class RegionStatisticsArray, class Neighborhood> 00508 inline void 00509 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1, 00510 pair<SeedImageIterator, SeedAccessor> img3, 00511 pair<DestImageIterator, DestAccessor> img4, 00512 RegionStatisticsArray & stats, 00513 SRGType srgType, Neighborhood n, double max_cost) 00514 { 00515 seededRegionGrowing3D(img1.first, img1.second, img1.third, 00516 img3.first, img3.second, 00517 img4.first, img4.second, 00518 stats, srgType, n, max_cost); 00519 } 00520 00521 template <class SrcImageIterator, class Shape, class SrcAccessor, 00522 class SeedImageIterator, class SeedAccessor, 00523 class DestImageIterator, class DestAccessor, 00524 class RegionStatisticsArray, class Neighborhood> 00525 inline void 00526 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1, 00527 pair<SeedImageIterator, SeedAccessor> img3, 00528 pair<DestImageIterator, DestAccessor> img4, 00529 RegionStatisticsArray & stats, 00530 SRGType srgType, Neighborhood n) 00531 { 00532 seededRegionGrowing3D(img1.first, img1.second, img1.third, 00533 img3.first, img3.second, 00534 img4.first, img4.second, 00535 stats, srgType, n, NumericTraits<double>::max()); 00536 } 00537 00538 template <class SrcImageIterator, class Shape, class SrcAccessor, 00539 class SeedImageIterator, class SeedAccessor, 00540 class DestImageIterator, class DestAccessor, 00541 class RegionStatisticsArray> 00542 inline void 00543 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1, 00544 pair<SeedImageIterator, SeedAccessor> img3, 00545 pair<DestImageIterator, DestAccessor> img4, 00546 RegionStatisticsArray & stats, SRGType srgType) 00547 { 00548 seededRegionGrowing3D(img1.first, img1.second, img1.third, 00549 img3.first, img3.second, 00550 img4.first, img4.second, 00551 stats, srgType, NeighborCode3DSix()); 00552 } 00553 00554 template <class SrcImageIterator, class Shape, class SrcAccessor, 00555 class SeedImageIterator, class SeedAccessor, 00556 class DestImageIterator, class DestAccessor, 00557 class RegionStatisticsArray> 00558 inline void 00559 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1, 00560 pair<SeedImageIterator, SeedAccessor> img3, 00561 pair<DestImageIterator, DestAccessor> img4, 00562 RegionStatisticsArray & stats) 00563 { 00564 seededRegionGrowing3D(img1.first, img1.second, img1.third, 00565 img3.first, img3.second, 00566 img4.first, img4.second, 00567 stats); 00568 } 00569 00570 } // namespace vigra 00571 00572 #endif // VIGRA_SEEDEDREGIONGROWING_HXX 00573
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|