[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

non_local_mean.hxx |

1 /************************************************************************/
2 /* */
3 /* Copyright 2014 by Thorsten Beier and Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 /* RE-IMPLEMENTAION OF THE WORK OF: */
36 /* Pierrick Coupe - pierrick.coupe@gmail.com */
37 /* Jose V. Manjon - jmanjon@fis.upv.es */
38 /* Brain Imaging Center, Montreal Neurological Institute. */
39 /* Mc Gill University */
40 /* */
41 /************************************************************************/
42 /* The ONLM filter is described in: */
43 /* */
44 /* P. Coup���, P. Yger, S. Prima, P. Hellier, C. Kervrann, C. Barillot. */
45 /* An Optimized Blockwise Non Local Means Denoising Filter */
46 /* for 3D Magnetic Resonance Images */
47 /* . IEEE Transactions on Medical Imaging, 27(4):425-441, */
48 /* Avril 2008 */
49 /************************************************************************/
50
51
52 #ifndef VIGRA_NON_LOCAL_MEAN
53 #define VIGRA_NON_LOCAL_MEAN
54
55
56 /*std*/
57 #include <iomanip>
58
59 /*vigra*/
60 #include "multi_array.hxx"
61 #include "multi_convolution.hxx"
62 #include "error.hxx"
63 #include "threading.hxx"
64 #include "gaussians.hxx"
65
66 namespace vigra{
67
68 struct NonLocalMeanParameter{
69
70 NonLocalMeanParameter(
71 const double sigmaSpatial = 2.0,
72 const int searchRadius = 3,
73 const int patchRadius = 1,
74 const double sigmaMean = 1.0,
75 const int stepSize = 2,
76 const int iterations=1,
77 const int nThreads = 8,
78 const bool verbose = true
79 ):
80 sigmaSpatial_(sigmaSpatial),
81 searchRadius_(searchRadius),
82 patchRadius_(patchRadius),
83 sigmaMean_(sigmaMean),
84 stepSize_(stepSize),
85 iterations_(iterations),
86 nThreads_(nThreads),
87 verbose_(verbose){
88 }
89 double sigmaSpatial_;
90 int searchRadius_;
91 int patchRadius_;
92 double sigmaMean_;
93 int stepSize_;
94 int iterations_;
95 int nThreads_;
96 bool verbose_;
97 };
98
99
100 // this has no template since this should
101 // be the same class for different dimensions
102 // and pixel types
103 class RatioPolicyParameter{
104 public:
105 RatioPolicyParameter(
106 const double sigma = 5.0,
107 const double meanRatio = 0.95,
108 const double varRatio = 0.5,
109 const double epsilon = 0.00001
110 ):
111 sigma_(sigma),
112 meanRatio_(meanRatio),
113 varRatio_(varRatio),
114 epsilon_(epsilon){
115 }
116 double sigma_;
117 double meanRatio_;
118 double varRatio_;
119 double epsilon_;
120 };
121
122
123 template<class PIXEL_TYPE_IN>
124 class RatioPolicy{
125
126 public:
127 typedef RatioPolicyParameter ParameterType;
128 typedef typename NumericTraits<PIXEL_TYPE_IN>::RealPromote PixelType;
129 typedef typename NumericTraits<PIXEL_TYPE_IN>::ValueType ValueType;
130
131
132 RatioPolicy(const ParameterType & param)
133 : meanRatio_(static_cast<ValueType>(param.meanRatio_)),
134 varRatio_(static_cast<ValueType>(param.varRatio_)),
135 epsilon_(static_cast<ValueType>(param.epsilon_)),
136 sigmaSquared_(param.sigma_*param.sigma_){
137
138 }
139
140 bool usePixel(const PixelType & meanA, const PixelType & varA)const{
141 return sum(meanA) > epsilon_ && sum(varA) > epsilon_;
142 }
143
144
145 bool usePixelPair(
146 const PixelType & meanA, const PixelType & varA,
147 const PixelType & meanB, const PixelType & varB
148 )const{
149 // Compute mean ratio of mean and variance
150 const ValueType m = mean(meanA/meanB);
151 const ValueType v = mean(varA/varB);
152 return (m > meanRatio_ && m < (1.0 / meanRatio_) && v > varRatio_ && v < (1.0 / varRatio_));
153 }
154
155 ValueType distanceToWeight(const PixelType & meanA, const PixelType & varA, const ValueType distance){
156 return exp(-distance /sigmaSquared_);
157 }
158
159 private:
160 ValueType meanRatio_;
161 ValueType varRatio_;
162 ValueType epsilon_;
163 ValueType sigmaSquared_;
164 };
165
166
167
168 // this has no template since this should
169 // be the same class for different dimensions
170 // and pixel types
171 class NormPolicyParameter{
172 public:
173 NormPolicyParameter(
174 const double sigma = 5.0,
175 const double meanDist = 0.95,
176 const double varRatio = 0.5,
177 const double epsilon = 0.00001
178 ):
179 sigma_(sigma),
180 meanDist_(meanDist),
181 varRatio_(varRatio),
182 epsilon_(epsilon){
183 }
184 double sigma_;
185 double meanDist_;
186 double varRatio_;
187 double epsilon_;
188 };
189
190
191
192 template<class V,int SIZE>
193 inline bool equal(const TinyVector<V,SIZE> & a,const TinyVector<V,SIZE> b){
194 for(int i=0;i<SIZE;++i)
195 if(a[i]!=b[i])
196 return false;
197 return true;
198 }
199
200
201
202 template<int DIM, bool ALWAYS_INSIDE>
203 struct BorderHelper;
204
205 template< int DIM>
206 struct BorderHelper<DIM,true>{
207 template<class COORD,class IMAGE>
208 static bool isInside(const COORD & ,const IMAGE & ){
209 return true;
210 }
211 template<class COORD,class IMAGE>
212 static bool isOutside(const COORD & ,const IMAGE & ){
213 return false;
214 }
215
216 template<class COORD,class IMAGE>
217 static void mirrorIfIsOutsidePoint(COORD & ,IMAGE & ){
218 }
219
220 };
221
222 template< int DIM>
223 struct BorderHelper<DIM,false>{
224 template<class COORD,class IMAGE>
225 static bool isInside(const COORD & c,const IMAGE & img){
226 return img.isInside(c);
227 }
228 template<class COORD,class IMAGE>
229 static bool isOutside(const COORD & c,const IMAGE & img){
230 return img.isOutside(c);
231 }
232
233 template<class COORD,class IMAGE>
234 static void mirrorIfIsOutsidePoint(COORD & coord,const IMAGE & img){
235 for(int c=0;c<DIM;++c){
236 if(coord[c]<0)
237 coord[c]=-1*coord[c];
238 else if(coord[c]>= img.shape(c))
239 coord[c] = 2 * img.shape(c) - coord[c] - 1;
240 }
241 }
242 };
243
244
245
246
247
248
249 template<class PIXEL_TYPE_IN>
250 class NormPolicy{
251
252 public:
253 typedef NormPolicyParameter ParameterType;
254 typedef typename NumericTraits<PIXEL_TYPE_IN>::RealPromote PixelType;
255 typedef typename NumericTraits<PIXEL_TYPE_IN>::ValueType ValueType;
256
257
258 NormPolicy(const ParameterType & param)
259 : meanDist_(static_cast<ValueType>(param.meanDist_)),
260 varRatio_(static_cast<ValueType>(param.varRatio_)),
261 epsilon_(static_cast<ValueType>(param.epsilon_)),
262 sigmaSquared_(param.sigma_*param.sigma_){
263
264 }
265
266 bool usePixel(const PixelType & meanA, const PixelType & varA)const{
267 return sum(varA)>epsilon_;
268 }
269
270
271 bool usePixelPair(
272 const PixelType & meanA, const PixelType & varA,
273 const PixelType & meanB, const PixelType & varB
274 )const{
275 // Compute mean ratio of mean and variance
276 const ValueType m = squaredNorm(meanA-meanB);
277 const ValueType v = mean(varA/varB);
278 //std::cout<<"norms "<<m<<" v "<<v<<"\n";
279 return (m < meanDist_ && v > varRatio_ && v < (1.0 / varRatio_));
280 }
281
282 ValueType distanceToWeight(const PixelType & meanA, const PixelType & varA, const ValueType distance){
283 return exp(-distance /sigmaSquared_);
284 }
285
286 private:
287 ValueType meanDist_;
288 ValueType varRatio_;
289 ValueType epsilon_;
290 ValueType sigmaSquared_;
291 };
292
293
294 // UnstridedArrayTag
295
296 template<int DIM, class PIXEL_TYPE_IN, class SMOOTH_POLICY>
297 class BlockWiseNonLocalMeanThreadObject{
298 typedef PIXEL_TYPE_IN PixelTypeIn;
299 typedef typename NumericTraits<PixelTypeIn>::RealPromote RealPromotePixelType;
300 typedef typename NumericTraits<RealPromotePixelType>::ValueType RealPromoteScalarType;
301
302 typedef typename MultiArray<DIM,int>::difference_type Coordinate;
303 typedef NonLocalMeanParameter ParameterType;
304
305 public:
306 typedef void result_type;
307 typedef MultiArrayView<DIM,PixelTypeIn> InArrayView;
308 typedef MultiArrayView<DIM,RealPromotePixelType> MeanArrayView;
309 typedef MultiArrayView<DIM,RealPromotePixelType> VarArrayView;
310 typedef MultiArrayView<DIM,RealPromotePixelType> EstimateArrayView;
311 typedef MultiArrayView<DIM,RealPromoteScalarType> LabelArrayView;
312 typedef std::vector<RealPromotePixelType> BlockAverageVectorType;
313 typedef std::vector<RealPromoteScalarType> BlockGaussWeightVectorType;
314 typedef SMOOTH_POLICY SmoothPolicyType;
315 // range type
316 typedef TinyVector<int,2> RangeType;
317 //typedef boost::mutex MutexType;
318
319 typedef threading::thread ThreadType;
320 typedef threading::mutex MutexType;
321
322 BlockWiseNonLocalMeanThreadObject(
323 const InArrayView & inImage,
324 MeanArrayView & meanImage,
325 VarArrayView & varImage,
326 EstimateArrayView & estimageImage,
327 LabelArrayView & labelImage,
328 const SmoothPolicyType & smoothPolicy,
329 const ParameterType & param,
330 const size_t nThreads,
331 MutexType & estimateMutex,
332 MultiArray<1,int> & progress
333 )
334 :
335 inImage_(inImage),
336 meanImage_(meanImage),
337 varImage_(varImage),
338 estimageImage_(estimageImage),
339 labelImage_(labelImage),
340 smoothPolicy_(smoothPolicy),
341 param_(param),
342 lastAxisRange_(),
343 threadIndex_(),
344 nThreads_(nThreads),
345 estimateMutexPtr_(&estimateMutex),
346 progress_(progress),
347 average_(std::pow( (double)(2*param.patchRadius_+1), DIM) ),
348 gaussWeight_(std::pow( (double)(2*param.patchRadius_+1), DIM) ),
349 shape_(inImage.shape()),
350 totalSize_()
351 {
352 totalSize_ = 1;
353 for(int dim=0;dim<DIM;++dim)
354 totalSize_*=(shape_[dim]/param.stepSize_);
355 }
356
357 void setRange(const RangeType & lastAxisRange){
358 lastAxisRange_=lastAxisRange;
359 }
360 void setThreadIndex(const size_t threadIndex){
361 threadIndex_=threadIndex;
362 }
363
364
365 void operator()();
366
367 private:
368
369 template<bool ALWAYS_INSIDE>
370 void processSinglePixel(const Coordinate & xyz);
371
372 template<bool ALWAYS_INSIDE>
373 void processSinglePair( const Coordinate & xyz,const Coordinate & nxyz,RealPromoteScalarType & wmax,RealPromoteScalarType & totalweight);
374
375 template<bool ALWAYS_INSIDE>
376 RealPromoteScalarType patchDistance(const Coordinate & xyz,const Coordinate & nxyz);
377
378 template<bool ALWAYS_INSIDE>
379 void patchExtractAndAcc(const Coordinate & xyz,const RealPromoteScalarType weight);
380
381 template<bool ALWAYS_INSIDE>
382 void patchAccMeanToEstimate(const Coordinate & xyz,const RealPromoteScalarType globalSum);
383
384
385 bool isAlwaysInside(const Coordinate & coord)const{
386 const Coordinate r = (Coordinate(param_.searchRadius_) + Coordinate(param_.patchRadius_) +1 );
387 const Coordinate test1 = coord - r;
388 const Coordinate test2 = coord + r;
389 return inImage_.isInside(test1) && inImage_.isInside(test2);
390 }
391
392 void initalizeGauss();
393
394 void progressPrinter(const int counter);
395
396
397 // array views
398 InArrayView inImage_;
399 MeanArrayView meanImage_;
400 VarArrayView varImage_;
401 EstimateArrayView estimageImage_;
402 LabelArrayView labelImage_;
403
404 // policy object
405 SmoothPolicyType smoothPolicy_;
406
407 // param obj.
408 ParameterType param_;
409
410 // thread related;
411 RangeType lastAxisRange_;
412 size_t threadIndex_;
413 size_t nThreads_;
414 MutexType * estimateMutexPtr_;
415 MultiArrayView<1,int> progress_;
416
417
418 // computations
419 BlockAverageVectorType average_;
420 BlockGaussWeightVectorType gaussWeight_;
421 Coordinate shape_;
422 size_t totalSize_;
423 };
424
425
426 template<int DIM,class PIXEL_TYPE_IN, class SMOOTH_POLICY>
427 inline void BlockWiseNonLocalMeanThreadObject<DIM, PIXEL_TYPE_IN, SMOOTH_POLICY>::initalizeGauss(){
428 Coordinate xyz;
429 const int pr = param_.patchRadius_;
430 Gaussian<RealPromoteScalarType> gaussian(param_.sigmaSpatial_);
431 int c=0;
432 RealPromoteScalarType sum = RealPromoteScalarType(0.0);
433 if(DIM==2){
434 for (xyz[1] = -pr; xyz[1] <=pr; ++xyz[1])
435 for (xyz[0] = -pr; xyz[0] <=pr; ++xyz[0],++c){
436 const RealPromoteScalarType distance = norm(xyz);
437 const RealPromoteScalarType w =gaussian(distance);
438 sum+=w;
439 gaussWeight_[c]=w;
440 }
441 }
442 if(DIM==3){
443 for (xyz[2] = -pr; xyz[2] <=pr; ++xyz[2])
444 for (xyz[1] = -pr; xyz[1] <=pr; ++xyz[1])
445 for (xyz[0] = -pr; xyz[0] <=pr; ++xyz[0],++c){
446 const RealPromoteScalarType distance = norm(xyz);
447 const RealPromoteScalarType w =gaussian(distance);
448 sum+=w;
449 gaussWeight_[c]=w;
450 }
451 }
452 if(DIM==4){
453 for (xyz[3] = -pr; xyz[3] <=pr; ++xyz[3])
454 for (xyz[2] = -pr; xyz[2] <=pr; ++xyz[2])
455 for (xyz[1] = -pr; xyz[1] <=pr; ++xyz[1])
456 for (xyz[0] = -pr; xyz[0] <=pr; ++xyz[0],++c){
457 const RealPromoteScalarType distance = norm(xyz);
458 const RealPromoteScalarType w =gaussian(distance);
459 sum+=w;
460 gaussWeight_[c]=w;
461 }
462 }
463 // normalize
464 for(size_t i=0;i<gaussWeight_.size();++i){
465 gaussWeight_[i]/=sum;
466 }
467 }
468
469
470 template<int DIM,class PIXEL_TYPE_IN, class SMOOTH_POLICY>
471 inline void BlockWiseNonLocalMeanThreadObject<DIM, PIXEL_TYPE_IN, SMOOTH_POLICY>::progressPrinter(const int counter){
472 progress_[threadIndex_] = counter;
473 if(threadIndex_==nThreads_-1){
474 if(counter%100 == 0){
475 int c=0;
476 for(size_t ti=0;ti<nThreads_;++ti)
477 c+=progress_[ti];
478 double pr=c;
479 pr/=totalSize_;
480 pr*=100.0;
481 std::cout<<"\rprogress "<<std::setw(10)<<pr<<" %%"<<std::flush;
482 }
483 }
484 }
485
486 template<int DIM,class PIXEL_TYPE_IN, class SMOOTH_POLICY>
487 void BlockWiseNonLocalMeanThreadObject<DIM, PIXEL_TYPE_IN, SMOOTH_POLICY>::operator()(){
488 const int start = lastAxisRange_[0];
489 const int end = lastAxisRange_[1];
490 const int stepSize = param_.stepSize_;
491
492 this->initalizeGauss();
493
494 Coordinate xyz;
495 int c=0;
496 if(param_.verbose_ && threadIndex_==nThreads_-1){
497 std::cout<<"progress";
498 }
499
500 if(DIM==2){
501 for (xyz[1] = start; xyz[1] < end; xyz[1] += stepSize)
502 for (xyz[0] = 0; xyz[0] < shape_[0]; xyz[0] += stepSize){
503
504 if(isAlwaysInside(xyz))
505 this->processSinglePixel<true>(xyz);
506 else
507 this->processSinglePixel<false>(xyz);
508 if(param_.verbose_)
509 this->progressPrinter(c);
510 ++c;
511 }
512 }
513 if(DIM==3){
514 for (xyz[2] = start; xyz[2] < end; xyz[2] += stepSize)
515 for (xyz[1] = 0; xyz[1] < shape_[1]; xyz[1] += stepSize)
516 for (xyz[0] = 0; xyz[0] < shape_[0]; xyz[0] += stepSize){
517 if(isAlwaysInside(xyz))
518 this->processSinglePixel<true>(xyz);
519 else
520 this->processSinglePixel<false>(xyz);
521 if(param_.verbose_)
522 this->progressPrinter(c);
523 ++c;
524 }
525 }
526 if(DIM==4){
527 for (xyz[3] = start; xyz[3] < end; xyz[3] += stepSize)
528 for (xyz[2] = 0; xyz[2] < shape_[2]; xyz[2] += stepSize)
529 for (xyz[1] = 0; xyz[1] < shape_[1]; xyz[1] += stepSize)
530 for (xyz[0] = 0; xyz[0] < shape_[0]; xyz[0] += stepSize){
531 if(isAlwaysInside(xyz))
532 this->processSinglePixel<true>(xyz);
533 else
534 this->processSinglePixel<false>(xyz);
535 if(param_.verbose_)
536 this->progressPrinter(c);
537 ++c;
538 }
539 }
540 if(param_.verbose_ && threadIndex_==nThreads_-1){
541 std::cout<<"\rprogress "<<std::setw(10)<<"100"<<" %%"<<"\n";
542 }
543 }
544
545
546
547 template<int DIM,class PIXEL_TYPE_IN, class SMOOTH_POLICY>
548 template<bool ALWAYS_INSIDE>
549 inline void BlockWiseNonLocalMeanThreadObject<DIM, PIXEL_TYPE_IN, SMOOTH_POLICY>::processSinglePixel(
550 const Coordinate & xyz
551 ){
552 Coordinate nxyz(SkipInitialization);
553 const int searchRadius = param_.searchRadius_;
554 std::fill(average_.begin(),average_.end(),RealPromotePixelType(0.0));
555 RealPromoteScalarType totalweight = 0.0;
556
557 if(smoothPolicy_.usePixel(meanImage_[xyz],varImage_[xyz])){
558
559 RealPromoteScalarType wmax = 0.0;
560 const Coordinate start = xyz- Coordinate(param_.searchRadius_);
561 const Coordinate end = xyz+ Coordinate(param_.searchRadius_);
562
563 if(DIM==2){
564 for (nxyz[1] = start[1]; nxyz[1] <= end[1]; nxyz[1]++)
565 for (nxyz[0] = start[0]; nxyz[0] <= end[0]; nxyz[0]++){
566 //nxyz = xyz + nxyz;
567 if(equal(nxyz,xyz))
568 continue;
569 this->processSinglePair<ALWAYS_INSIDE>(xyz,nxyz,wmax,totalweight);
570 }
571 }
572 else if(DIM==3){
573 for (nxyz[2] = start[2]; nxyz[2] <= end[2]; nxyz[2]++)
574 for (nxyz[1] = start[1]; nxyz[1] <= end[1]; nxyz[1]++)
575 for (nxyz[0] = start[0]; nxyz[0] <= end[0]; nxyz[0]++){
576 if(equal(nxyz,xyz))
577 continue;
578 this->processSinglePair<ALWAYS_INSIDE>(xyz,nxyz,wmax,totalweight);
579 }
580 }
581 else if(DIM==4){
582 for (nxyz[3] = start[3]; nxyz[3] <= end[3]; nxyz[3]++)
583 for (nxyz[2] = start[2]; nxyz[2] <= end[2]; nxyz[2]++)
584 for (nxyz[1] = start[1]; nxyz[1] <= end[1]; nxyz[1]++)
585 for (nxyz[0] = start[0]; nxyz[0] <= end[0]; nxyz[0]++){
586 if(equal(nxyz,xyz))
587 continue;
588 this->processSinglePair<ALWAYS_INSIDE>(xyz,nxyz,wmax,totalweight);
589 }
590 }
591
592
593 if (wmax == 0.0){
594 wmax = 1.0;
595 }
596 // give pixel xyz as much weight as
597 // the maximum weighted other patch
598 this->patchExtractAndAcc<ALWAYS_INSIDE>(xyz,wmax);
599 totalweight += wmax;
600
601 // this if seems total useless to me
602 if (totalweight != 0.0){
603 this->patchAccMeanToEstimate<ALWAYS_INSIDE>(xyz,totalweight);
604 }
605
606 }
607 else{
608 const double wmax = 1.0;
609 this->patchExtractAndAcc<ALWAYS_INSIDE>(xyz,wmax);
610 totalweight += wmax;
611 this->patchAccMeanToEstimate<ALWAYS_INSIDE>(xyz,totalweight);
612 }
613 }
614
615
616
617 template<int DIM,class PIXEL_TYPE_IN, class SMOOTH_POLICY>
618 template<bool ALWAYS_INSIDE>
619 inline void BlockWiseNonLocalMeanThreadObject<DIM, PIXEL_TYPE_IN, SMOOTH_POLICY>::processSinglePair(
620 const Coordinate & xyz,
621 const Coordinate & nxyz,
622 RealPromoteScalarType & wmax,
623 RealPromoteScalarType & totalweight
624 ){
625
626 if(BorderHelper<DIM,ALWAYS_INSIDE>::isInside(nxyz,inImage_)){
627 //if(ALWAYS_INSIDE || inImage_.isInside(nxyz)){
628 if(smoothPolicy_.usePixel(meanImage_[nxyz],varImage_[nxyz])){
629 // here we check if to patches fit to each other
630 // one patch is around xyz
631 // other patch is arround nxyz
632 if(smoothPolicy_.usePixelPair(meanImage_[xyz],varImage_[xyz],meanImage_[nxyz],varImage_[nxyz])){
633 const RealPromoteScalarType distance =this->patchDistance<ALWAYS_INSIDE>(xyz,nxyz);
634 const RealPromoteScalarType w = smoothPolicy_.distanceToWeight(meanImage_[xyz],varImage_[xyz],distance);
635 wmax = std::max(w,wmax);
636 this->patchExtractAndAcc<ALWAYS_INSIDE>(nxyz,w);
637 totalweight+= w;
638 }
639 }
640 }
641 }
642
643
644
645 template<int DIM,class PIXEL_TYPE_IN, class SMOOTH_POLICY>
646 template<bool ALWAYS_INSIDE>
647 inline typename BlockWiseNonLocalMeanThreadObject<DIM, PIXEL_TYPE_IN, SMOOTH_POLICY>::RealPromoteScalarType
648 BlockWiseNonLocalMeanThreadObject<DIM,PIXEL_TYPE_IN,SMOOTH_POLICY>::patchDistance(
649 const Coordinate & pA,
650 const Coordinate & pB
651 ){
652
653 // TODO : use a acculator like think to make this more beautiful ?
654 const int f = param_.patchRadius_;
655 Coordinate offset(SkipInitialization),nPa(SkipInitialization),nPb(SkipInitialization);
656 int acu = 0;
657 RealPromoteScalarType distancetotal = 0;
658 int c =0 ;
659 //this->mirrorIfIsOutsidePoint<ALWAYS_INSIDE>(nPa);
660 // this->mirrorIfIsOutsidePoint<ALWAYS_INSIDE>(nPb);
661 #define VIGRA_NLM_IN_LOOP_CODE \
662 nPa = pA+offset; \
663 nPb = pB+offset; \
664 BorderHelper<DIM,ALWAYS_INSIDE>::mirrorIfIsOutsidePoint(nPa,inImage_); \
665 BorderHelper<DIM,ALWAYS_INSIDE>::mirrorIfIsOutsidePoint(nPb,inImage_); \
666 const RealPromoteScalarType gaussWeight = gaussWeight_[c]; \
667 const RealPromotePixelType vA = inImage_[nPa]; \
668 const RealPromotePixelType vB = inImage_[nPb]; \
669 distancetotal += gaussWeight*vigra::sizeDividedSquaredNorm(vA-vB); \
670 ++acu;
671
672 if(DIM==2){
673 for (offset[1] = -f; offset[1] <= f; ++offset[1])
674 for (offset[0] = -f; offset[0] <= f; ++offset[0],++c){
675 VIGRA_NLM_IN_LOOP_CODE;
676 }
677 }
678 else if(DIM==3){
679 for (offset[2] = -f; offset[2] <= f; ++offset[2])
680 for (offset[1] = -f; offset[1] <= f; ++offset[1])
681 for (offset[0] = -f; offset[0] <= f; ++offset[0],++c){
682 VIGRA_NLM_IN_LOOP_CODE;
683 }
684 }
685 else if(DIM==4){
686 for (offset[3] = -f; offset[3] <= f; ++offset[3])
687 for (offset[2] = -f; offset[2] <= f; ++offset[2])
688 for (offset[1] = -f; offset[1] <= f; ++offset[1])
689 for (offset[0] = -f; offset[0] <= f; ++offset[0],++c){
690 VIGRA_NLM_IN_LOOP_CODE;
691 }
692 }
693 #undef VIGRA_NLM_IN_LOOP_CODE
694 return distancetotal / acu;
695 }
696
697
698
699 template<int DIM,class PIXEL_TYPE_IN, class SMOOTH_POLICY>
700 template<bool ALWAYS_INSIDE>
701 inline void
702 BlockWiseNonLocalMeanThreadObject<DIM,PIXEL_TYPE_IN,SMOOTH_POLICY>::patchExtractAndAcc(
703 const Coordinate & xyz,
704 const RealPromoteScalarType weight
705 ){
706 Coordinate xyzPos(SkipInitialization),abc(SkipInitialization);
707 Coordinate nhSize3(param_.patchRadius_);
708 const int ns = 2 * param_.patchRadius_ + 1;
709 int count = 0;
710
711 // todo: remove abc vector
712
713 #define VIGRA_NLM_IN_LOOP_CODE \
714 xyzPos = xyz + abc - nhSize3; \
715 if(BorderHelper<DIM,ALWAYS_INSIDE>::isOutside(xyzPos,inImage_)) \
716 average_[count] += inImage_[xyz]* weight; \
717 else \
718 average_[count] += inImage_[xyzPos]* weight; \
719 count++
720
721 if(DIM==2){
722 for (abc[1] = 0; abc[1] < ns; abc[1]++)
723 for (abc[0] = 0; abc[0] < ns; abc[0]++){
724 VIGRA_NLM_IN_LOOP_CODE;
725 }
726 }
727 else if(DIM==3){
728 for (abc[2] = 0; abc[2] < ns; abc[2]++)
729 for (abc[1] = 0; abc[1] < ns; abc[1]++)
730 for (abc[0] = 0; abc[0] < ns; abc[0]++){
731 VIGRA_NLM_IN_LOOP_CODE;
732 }
733 }
734 else if(DIM==4){
735 for (abc[3] = 0; abc[3] < ns; abc[3]++)
736 for (abc[2] = 0; abc[2] < ns; abc[2]++)
737 for (abc[1] = 0; abc[1] < ns; abc[1]++)
738 for (abc[0] = 0; abc[0] < ns; abc[0]++){
739 VIGRA_NLM_IN_LOOP_CODE;
740 }
741 }
742
743 #undef VIGRA_NLM_IN_LOOP_CODE
744 }
745
746
747 template<int DIM,class PIXEL_TYPE_IN, class SMOOTH_POLICY>
748 template<bool ALWAYS_INSIDE>
749 inline void
750 BlockWiseNonLocalMeanThreadObject<DIM,PIXEL_TYPE_IN,SMOOTH_POLICY>::patchAccMeanToEstimate(
751 const Coordinate & xyz,
752 const RealPromoteScalarType globalSum
753 ){
754 Coordinate abc(SkipInitialization),xyzPos(SkipInitialization),nhSize(param_.patchRadius_);
755 int count = 0 ;
756 const int ns = 2 * param_.patchRadius_ + 1;
757
758 #define VIGRA_NLM_IN_LOOP_CODE \
759 xyzPos = xyz + abc - nhSize; \
760 if(BorderHelper<DIM,ALWAYS_INSIDE>::isInside(xyzPos,inImage_)){ \
761 estimateMutexPtr_->lock(); \
762 RealPromotePixelType value = estimageImage_[xyzPos]; \
763 const RealPromoteScalarType gw = gaussWeight_[count]; \
764 RealPromotePixelType tmp =(average_[count] / globalSum); \
765 tmp*=gw; \
766 value +=tmp; \
767 estimageImage_[xyzPos] = value; \
768 labelImage_[xyzPos]+=gw; \
769 estimateMutexPtr_->unlock(); \
770 } \
771 count++
772
773 if(DIM==2){
774 for (abc[1] = 0; abc[1] < ns; abc[1]++)
775 for (abc[0] = 0; abc[0] < ns; abc[0]++){
776 VIGRA_NLM_IN_LOOP_CODE;
777 }
778 }
779 if(DIM==3){
780 for (abc[2] = 0; abc[2] < ns; abc[2]++)
781 for (abc[1] = 0; abc[1] < ns; abc[1]++)
782 for (abc[0] = 0; abc[0] < ns; abc[0]++){
783 VIGRA_NLM_IN_LOOP_CODE;
784 }
785 }
786 if(DIM==4){
787 for (abc[3] = 0; abc[3] < ns; abc[3]++)
788 for (abc[2] = 0; abc[2] < ns; abc[2]++)
789 for (abc[1] = 0; abc[1] < ns; abc[1]++)
790 for (abc[0] = 0; abc[0] < ns; abc[0]++){
791 VIGRA_NLM_IN_LOOP_CODE;
792 }
793 }
794 #undef VIGRA_NLM_IN_LOOP_CODE
795 }
796
797
798
799 template<int DIM,class PIXEL_TYPE_IN, class SMOOTH_POLICY,class PIXEL_TYPE_OUT>
800 inline void gaussianMeanAndVariance(
801 const vigra::MultiArrayView<DIM,PIXEL_TYPE_IN> & inArray,
802 const double sigma,
803 vigra::MultiArrayView<DIM,PIXEL_TYPE_OUT> & meanArray,
804 vigra::MultiArrayView<DIM,PIXEL_TYPE_OUT> & varArray,
805 vigra::MultiArrayView<DIM,PIXEL_TYPE_OUT> & tmpArray
806 ){
807
808 // compute mean and variance
809 vigra::gaussianSmoothMultiArray(inArray, meanArray, sigma);
810 // square raw data (use estimate Image to store temp. results)
811 for(int scanOrderIndex=0;scanOrderIndex<inArray.size();++scanOrderIndex){
812 tmpArray[scanOrderIndex]=vigra::pow(inArray[scanOrderIndex],2);
813 }
814 vigra::gaussianSmoothMultiArray(tmpArray,varArray, sigma);
815 for(int scanOrderIndex=0;scanOrderIndex<inArray.size();++scanOrderIndex){
816 PIXEL_TYPE_OUT var = varArray[scanOrderIndex] - vigra::pow(meanArray[scanOrderIndex],2);
817 var = clipLower(var);
818 //makeNegtiveValuesZero(var); // callbyref
819 varArray[scanOrderIndex] = var;
820 }
821 }
822
823 template<int DIM,class PIXEL_TYPE_IN, class SMOOTH_POLICY,class PIXEL_TYPE_OUT>
824 inline void gaussianMeanAndVariance(
825 const vigra::MultiArrayView<DIM,PIXEL_TYPE_IN> & inArray,
826 const double sigma,
827 vigra::MultiArrayView<DIM,PIXEL_TYPE_OUT> & meanArray,
828 vigra::MultiArrayView<DIM,PIXEL_TYPE_OUT> & varArray
829 ){
830 vigra::MultiArray<DIM,PIXEL_TYPE_OUT> tmpArray(inArray.shape());
831 gaussianMeanAndVariance<DIM,PIXEL_TYPE_IN,PIXEL_TYPE_OUT>(inArray,sigma,meanArray,varArray,tmpArray);
832 }
833
834 namespace detail_non_local_means{
835
836 template<int DIM, class PIXEL_TYPE_IN,class PIXEL_TYPE_OUT,class SMOOTH_POLICY>
837 void nonLocalMean1Run(
838 const vigra::MultiArrayView<DIM,PIXEL_TYPE_IN> & image,
839 const SMOOTH_POLICY & smoothPolicy,
840 const NonLocalMeanParameter param,
841 vigra::MultiArrayView<DIM,PIXEL_TYPE_OUT> outImage
842 ){
843
844 typedef PIXEL_TYPE_IN PixelTypeIn;
845 typedef typename vigra::NumericTraits<PixelTypeIn>::RealPromote RealPromotePixelType;
846 typedef typename vigra::NumericTraits<RealPromotePixelType>::ValueType RealPromoteScalarType;
847 typedef SMOOTH_POLICY SmoothPolicyType;
848
849 typedef BlockWiseNonLocalMeanThreadObject<DIM,PixelTypeIn,SmoothPolicyType> ThreadObjectType;
850
851
852 // inspect parameter
853 vigra_precondition(param.stepSize_>=1,"NonLocalMean Parameter: \"stepSize>=1\" violated");
854 vigra_precondition(param.searchRadius_>=1, "NonLocalMean Parameter: \"searchRadius >=1\" violated");
855 vigra_precondition(param.patchRadius_>=1,"NonLocalMean Parameter: \"searchRadius >=1\" violated");
856 vigra_precondition(param.stepSize_-1<=param.patchRadius_,"NonLocalMean Parameter: \"stepSize -1 <= patchRadius\" violated");
857
858 // allocate arrays
859 vigra::MultiArray<DIM,RealPromotePixelType> meanImage(image.shape());
860 vigra::MultiArray<DIM,RealPromotePixelType> varImage(image.shape());
861 vigra::MultiArray<DIM,RealPromotePixelType> estimageImage(image.shape());
862 vigra::MultiArray<DIM,RealPromoteScalarType> labelImage(image.shape());
863
864 // compute mean and variance
865 // last argument is a "buffer" since within "gaussianMeanAndVariance" another array is needed
866 // ==> to avoid an unnecessary allocation we use the estimageImage as a buffer
867 //gaussianMeanAndVariance<DIM,PixelTypeIn,RealPromotePixelType>(image,param.sigmaMean_,meanImage,varImage,estimageImage);
868 gaussianMeanAndVariance<DIM,PixelTypeIn,RealPromotePixelType>(image,param.sigmaMean_,meanImage,varImage);
869
870 // initialize
871 labelImage = RealPromoteScalarType(0.0);
872 estimageImage = RealPromotePixelType(0.0);
873
874 ///////////////////////////////////////////////////////////////
875 { // MULTI THREAD CODE STARTS HERE
876
877
878
879 typedef threading::thread ThreadType;
880 typedef threading::mutex MutexType;
881
882 MutexType estimateMutex;
883 //typedef boost::thread ThreadType;
884
885 const size_t nThreads = param.nThreads_;
886 MultiArray<1,int> progress = MultiArray<1,int>(typename MultiArray<1,int>::difference_type(nThreads));
887
888 // allocate all thread objects
889 // each thread object works on a portion of the data
890 std::vector<ThreadObjectType> threadObjects(nThreads,
891 ThreadObjectType(image, meanImage, varImage, estimageImage, labelImage,
892 smoothPolicy, param, nThreads, estimateMutex,progress)
893 );
894
895 // thread ptr
896 std::vector<ThreadType *> threadPtrs(nThreads);
897 for(size_t i=0; i<nThreads; ++i){
898 ThreadObjectType & threadObj = threadObjects[i];
899 threadObj.setThreadIndex(i);
900 typename ThreadObjectType::RangeType lastAxisRange;
901 lastAxisRange[0]=(i * image.shape(DIM-1)) / nThreads;
902 lastAxisRange[1]=((i+1) * image.shape(DIM-1)) / nThreads;
903 threadObj.setRange(lastAxisRange);
904 // this will start the threads and cal operator()
905 // of the threadObjects
906 threadPtrs[i] = new ThreadType(threadObjects[i]);
907 }
908 for(size_t i=0; i<nThreads; ++i)
909 threadPtrs[i]->join();
910 for(size_t i=0; i<nThreads; ++i)
911 delete threadPtrs[i];
912
913 } // MULTI THREAD CODE ENDS HERE
914 ///////////////////////////////////////////////////////////////
915
916 // normalize estimates by the number of labels
917 // and write that in output
918 for(int scanOrderIndex=0; scanOrderIndex<labelImage.size(); ++scanOrderIndex){
919 if (labelImage[scanOrderIndex] <= RealPromoteScalarType(0.00001))
920 outImage[scanOrderIndex]=image[scanOrderIndex];
921 else
922 outImage[scanOrderIndex]=estimageImage[scanOrderIndex] / labelImage[scanOrderIndex];
923 }
924 }
925
926 }
927
928 template<int DIM, class PIXEL_TYPE_IN,class PIXEL_TYPE_OUT,class SMOOTH_POLICY>
929 void nonLocalMean(
930 const vigra::MultiArrayView<DIM,PIXEL_TYPE_IN> & image,
931 const SMOOTH_POLICY & smoothPolicy,
932 const NonLocalMeanParameter param,
933 vigra::MultiArrayView<DIM,PIXEL_TYPE_OUT> outImage
934 ){
935 detail_non_local_means::nonLocalMean1Run<DIM,PIXEL_TYPE_IN,PIXEL_TYPE_OUT,SMOOTH_POLICY>(image,smoothPolicy,param,outImage);
936 if(param.iterations_>1){
937
938 vigra::MultiArray<DIM,PIXEL_TYPE_OUT> tmp(outImage.shape());
939 for(size_t i=0;i<param.iterations_-1;++i){
940 tmp=outImage;
941 detail_non_local_means::nonLocalMean1Run<DIM,PIXEL_TYPE_OUT,PIXEL_TYPE_OUT,SMOOTH_POLICY>(tmp,smoothPolicy,param,outImage);
942 }
943 }
944 }
945
946
947
948
949 } // end namespace vigra
950
951
952 #endif
vigra::labelImage
unsigned int labelImage(...)
Find the connected components of a segmented image.
vigra::linalg::exp
linalg::TemporaryMatrix< T > exp(MultiArrayView< 2, T, C > const &v)
vigra::gaussianSmoothMultiArray
void gaussianSmoothMultiArray(...)
Isotropic Gaussian smoothing of a multi-dimensional arrays.
std
Definition: array_vector.hxx:954
vigra::squaredNorm
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition: fftw3.hxx:1044
vigra::MultiArray
Main MultiArray class containing the memory management.
Definition: multi_array.hxx:2422
vigra::MultiArrayView::size
difference_type_1 size() const
Definition: multi_array.hxx:1589
vigra
Definition: accessor.hxx:43
vigra::MultiArray< DIM, int >::difference_type
view_type::difference_type difference_type
Definition: multi_array.hxx:2470
vigra::norm
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition: fftw3.hxx:1037
vigra::MultiArrayView::shape
const difference_type & shape() const
Definition: multi_array.hxx:1596
vigra::sum
NumericTraits< V >::Promote sum(TinyVectorBase< V, SIZE, D1, D2 > const &l)
sum of the vector&#39;s elements
Definition: tinyvector.hxx:2073
vigra::pow
TinyVector< V, SIZE > pow(TinyVectorBase< V, SIZE, D1, D2 > const &v, E exponent)
Definition: tinyvector.hxx:2036
vigra::MultiArrayView
Base class for, and view to, vigra::MultiArray.
Definition: multi_array.hxx:652
vigra::clipLower
TinyVector< V, SIZE > clipLower(TinyVector< V, SIZE > const &t)
Clip negative values.
Definition: tinyvector.hxx:2268

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python vigra 1.11.0