[ 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. Coupe, 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 std::fill(average_.begin(),average_.end(),RealPromotePixelType(0.0));
554 RealPromoteScalarType totalweight = 0.0;
555
556 if(smoothPolicy_.usePixel(meanImage_[xyz],varImage_[xyz])){
557
558 RealPromoteScalarType wmax = 0.0;
559 const Coordinate start = xyz- Coordinate(param_.searchRadius_);
560 const Coordinate end = xyz+ Coordinate(param_.searchRadius_);
561
562 if(DIM==2){
563 for (nxyz[1] = start[1]; nxyz[1] <= end[1]; nxyz[1]++)
564 for (nxyz[0] = start[0]; nxyz[0] <= end[0]; nxyz[0]++){
565 //nxyz = xyz + nxyz;
566 if(equal(nxyz,xyz))
567 continue;
568 this->processSinglePair<ALWAYS_INSIDE>(xyz,nxyz,wmax,totalweight);
569 }
570 }
571 else if(DIM==3){
572 for (nxyz[2] = start[2]; nxyz[2] <= end[2]; nxyz[2]++)
573 for (nxyz[1] = start[1]; nxyz[1] <= end[1]; nxyz[1]++)
574 for (nxyz[0] = start[0]; nxyz[0] <= end[0]; nxyz[0]++){
575 if(equal(nxyz,xyz))
576 continue;
577 this->processSinglePair<ALWAYS_INSIDE>(xyz,nxyz,wmax,totalweight);
578 }
579 }
580 else if(DIM==4){
581 for (nxyz[3] = start[3]; nxyz[3] <= end[3]; nxyz[3]++)
582 for (nxyz[2] = start[2]; nxyz[2] <= end[2]; nxyz[2]++)
583 for (nxyz[1] = start[1]; nxyz[1] <= end[1]; nxyz[1]++)
584 for (nxyz[0] = start[0]; nxyz[0] <= end[0]; nxyz[0]++){
585 if(equal(nxyz,xyz))
586 continue;
587 this->processSinglePair<ALWAYS_INSIDE>(xyz,nxyz,wmax,totalweight);
588 }
589 }
590
591
592 if (wmax == 0.0){
593 wmax = 1.0;
594 }
595 // give pixel xyz as much weight as
596 // the maximum weighted other patch
597 this->patchExtractAndAcc<ALWAYS_INSIDE>(xyz,wmax);
598 totalweight += wmax;
599
600 // this if seems total useless to me
601 if (totalweight != 0.0){
602 this->patchAccMeanToEstimate<ALWAYS_INSIDE>(xyz,totalweight);
603 }
604
605 }
606 else{
607 const double wmax = 1.0;
608 this->patchExtractAndAcc<ALWAYS_INSIDE>(xyz,wmax);
609 totalweight += wmax;
610 this->patchAccMeanToEstimate<ALWAYS_INSIDE>(xyz,totalweight);
611 }
612}
613
614
615
616 template<int DIM,class PIXEL_TYPE_IN, class SMOOTH_POLICY>
617 template<bool ALWAYS_INSIDE>
618 inline void BlockWiseNonLocalMeanThreadObject<DIM, PIXEL_TYPE_IN, SMOOTH_POLICY>::processSinglePair(
619 const Coordinate & xyz,
620 const Coordinate & nxyz,
621 RealPromoteScalarType & wmax,
622 RealPromoteScalarType & totalweight
623){
624
625 if(BorderHelper<DIM,ALWAYS_INSIDE>::isInside(nxyz,inImage_)){
626 //if(ALWAYS_INSIDE || inImage_.isInside(nxyz)){
627 if(smoothPolicy_.usePixel(meanImage_[nxyz],varImage_[nxyz])){
628 // here we check if to patches fit to each other
629 // one patch is around xyz
630 // other patch is arround nxyz
631 if(smoothPolicy_.usePixelPair(meanImage_[xyz],varImage_[xyz],meanImage_[nxyz],varImage_[nxyz])){
632 const RealPromoteScalarType distance =this->patchDistance<ALWAYS_INSIDE>(xyz,nxyz);
633 const RealPromoteScalarType w = smoothPolicy_.distanceToWeight(meanImage_[xyz],varImage_[xyz],distance);
634 wmax = std::max(w,wmax);
635 this->patchExtractAndAcc<ALWAYS_INSIDE>(nxyz,w);
636 totalweight+= w;
637 }
638 }
639 }
640}
641
642
643
644 template<int DIM,class PIXEL_TYPE_IN, class SMOOTH_POLICY>
645 template<bool ALWAYS_INSIDE>
646 inline typename BlockWiseNonLocalMeanThreadObject<DIM, PIXEL_TYPE_IN, SMOOTH_POLICY>::RealPromoteScalarType
647BlockWiseNonLocalMeanThreadObject<DIM,PIXEL_TYPE_IN,SMOOTH_POLICY>::patchDistance(
648 const Coordinate & pA,
649 const Coordinate & pB
650){
651
652 // TODO : use a acculator like think to make this more beautiful ?
653 const int f = param_.patchRadius_;
654 Coordinate offset(SkipInitialization),nPa(SkipInitialization),nPb(SkipInitialization);
655 int acu = 0;
656 RealPromoteScalarType distancetotal = 0;
657 int c =0 ;
658 //this->mirrorIfIsOutsidePoint<ALWAYS_INSIDE>(nPa);
659 // this->mirrorIfIsOutsidePoint<ALWAYS_INSIDE>(nPb);
660 #define VIGRA_NLM_IN_LOOP_CODE \
661 nPa = pA+offset; \
662 nPb = pB+offset; \
663 BorderHelper<DIM,ALWAYS_INSIDE>::mirrorIfIsOutsidePoint(nPa,inImage_); \
664 BorderHelper<DIM,ALWAYS_INSIDE>::mirrorIfIsOutsidePoint(nPb,inImage_); \
665 const RealPromoteScalarType gaussWeight = gaussWeight_[c]; \
666 const RealPromotePixelType vA = inImage_[nPa]; \
667 const RealPromotePixelType vB = inImage_[nPb]; \
668 distancetotal += gaussWeight*vigra::sizeDividedSquaredNorm(vA-vB); \
669 ++acu;
670
671 if(DIM==2){
672 for (offset[1] = -f; offset[1] <= f; ++offset[1])
673 for (offset[0] = -f; offset[0] <= f; ++offset[0],++c){
674 VIGRA_NLM_IN_LOOP_CODE;
675 }
676 }
677 else if(DIM==3){
678 for (offset[2] = -f; offset[2] <= f; ++offset[2])
679 for (offset[1] = -f; offset[1] <= f; ++offset[1])
680 for (offset[0] = -f; offset[0] <= f; ++offset[0],++c){
681 VIGRA_NLM_IN_LOOP_CODE;
682 }
683 }
684 else if(DIM==4){
685 for (offset[3] = -f; offset[3] <= f; ++offset[3])
686 for (offset[2] = -f; offset[2] <= f; ++offset[2])
687 for (offset[1] = -f; offset[1] <= f; ++offset[1])
688 for (offset[0] = -f; offset[0] <= f; ++offset[0],++c){
689 VIGRA_NLM_IN_LOOP_CODE;
690 }
691 }
692 #undef VIGRA_NLM_IN_LOOP_CODE
693 return distancetotal / acu;
694}
695
696
697
698 template<int DIM,class PIXEL_TYPE_IN, class SMOOTH_POLICY>
699 template<bool ALWAYS_INSIDE>
700 inline void
701BlockWiseNonLocalMeanThreadObject<DIM,PIXEL_TYPE_IN,SMOOTH_POLICY>::patchExtractAndAcc(
702 const Coordinate & xyz,
703 const RealPromoteScalarType weight
704){
705 Coordinate xyzPos(SkipInitialization),abc(SkipInitialization);
706 Coordinate nhSize3(param_.patchRadius_);
707 const int ns = 2 * param_.patchRadius_ + 1;
708 int count = 0;
709
710 // todo: remove abc vector
711
712 #define VIGRA_NLM_IN_LOOP_CODE \
713 xyzPos = xyz + abc - nhSize3; \
714 if(BorderHelper<DIM,ALWAYS_INSIDE>::isOutside(xyzPos,inImage_)) \
715 average_[count] += inImage_[xyz]* weight; \
716 else \
717 average_[count] += inImage_[xyzPos]* weight; \
718 count++
719
720 if(DIM==2){
721 for (abc[1] = 0; abc[1] < ns; abc[1]++)
722 for (abc[0] = 0; abc[0] < ns; abc[0]++){
723 VIGRA_NLM_IN_LOOP_CODE;
724 }
725 }
726 else if(DIM==3){
727 for (abc[2] = 0; abc[2] < ns; abc[2]++)
728 for (abc[1] = 0; abc[1] < ns; abc[1]++)
729 for (abc[0] = 0; abc[0] < ns; abc[0]++){
730 VIGRA_NLM_IN_LOOP_CODE;
731 }
732 }
733 else if(DIM==4){
734 for (abc[3] = 0; abc[3] < ns; abc[3]++)
735 for (abc[2] = 0; abc[2] < ns; abc[2]++)
736 for (abc[1] = 0; abc[1] < ns; abc[1]++)
737 for (abc[0] = 0; abc[0] < ns; abc[0]++){
738 VIGRA_NLM_IN_LOOP_CODE;
739 }
740 }
741
742 #undef VIGRA_NLM_IN_LOOP_CODE
743}
744
745
746 template<int DIM,class PIXEL_TYPE_IN, class SMOOTH_POLICY>
747 template<bool ALWAYS_INSIDE>
748 inline void
749BlockWiseNonLocalMeanThreadObject<DIM,PIXEL_TYPE_IN,SMOOTH_POLICY>::patchAccMeanToEstimate(
750 const Coordinate & xyz,
751 const RealPromoteScalarType globalSum
752){
753 Coordinate abc(SkipInitialization),xyzPos(SkipInitialization),nhSize(param_.patchRadius_);
754 int count = 0 ;
755 const int ns = 2 * param_.patchRadius_ + 1;
756
757 #define VIGRA_NLM_IN_LOOP_CODE \
758 xyzPos = xyz + abc - nhSize; \
759 if(BorderHelper<DIM,ALWAYS_INSIDE>::isInside(xyzPos,inImage_)){ \
760 estimateMutexPtr_->lock(); \
761 RealPromotePixelType value = estimageImage_[xyzPos]; \
762 const RealPromoteScalarType gw = gaussWeight_[count]; \
763 RealPromotePixelType tmp =(average_[count] / globalSum); \
764 tmp*=gw; \
765 value +=tmp; \
766 estimageImage_[xyzPos] = value; \
767 labelImage_[xyzPos]+=gw; \
768 estimateMutexPtr_->unlock(); \
769 } \
770 count++
771
772 if(DIM==2){
773 for (abc[1] = 0; abc[1] < ns; abc[1]++)
774 for (abc[0] = 0; abc[0] < ns; abc[0]++){
775 VIGRA_NLM_IN_LOOP_CODE;
776 }
777 }
778 if(DIM==3){
779 for (abc[2] = 0; abc[2] < ns; abc[2]++)
780 for (abc[1] = 0; abc[1] < ns; abc[1]++)
781 for (abc[0] = 0; abc[0] < ns; abc[0]++){
782 VIGRA_NLM_IN_LOOP_CODE;
783 }
784 }
785 if(DIM==4){
786 for (abc[3] = 0; abc[3] < ns; abc[3]++)
787 for (abc[2] = 0; abc[2] < ns; abc[2]++)
788 for (abc[1] = 0; abc[1] < ns; abc[1]++)
789 for (abc[0] = 0; abc[0] < ns; abc[0]++){
790 VIGRA_NLM_IN_LOOP_CODE;
791 }
792 }
793 #undef VIGRA_NLM_IN_LOOP_CODE
794}
795
796
797
798 template<int DIM,class PIXEL_TYPE_IN, class SMOOTH_POLICY,class PIXEL_TYPE_OUT>
799 inline void gaussianMeanAndVariance(
800 const vigra::MultiArrayView<DIM,PIXEL_TYPE_IN> & inArray,
801 const double sigma,
802 vigra::MultiArrayView<DIM,PIXEL_TYPE_OUT> & meanArray,
803 vigra::MultiArrayView<DIM,PIXEL_TYPE_OUT> & varArray,
804 vigra::MultiArrayView<DIM,PIXEL_TYPE_OUT> & tmpArray
805){
806
807 // compute mean and variance
808 vigra::gaussianSmoothMultiArray(inArray, meanArray, sigma);
809 // square raw data (use estimate Image to store temp. results)
810 for(int scanOrderIndex=0;scanOrderIndex<inArray.size();++scanOrderIndex){
811 tmpArray[scanOrderIndex]=vigra::pow(inArray[scanOrderIndex],2);
812 }
813 vigra::gaussianSmoothMultiArray(tmpArray,varArray, sigma);
814 for(int scanOrderIndex=0;scanOrderIndex<inArray.size();++scanOrderIndex){
815 PIXEL_TYPE_OUT var = varArray[scanOrderIndex] - vigra::pow(meanArray[scanOrderIndex],2);
816 var = clipLower(var);
817 //makeNegtiveValuesZero(var); // callbyref
818 varArray[scanOrderIndex] = var;
819 }
820}
821
822 template<int DIM,class PIXEL_TYPE_IN, class SMOOTH_POLICY,class PIXEL_TYPE_OUT>
823 inline void gaussianMeanAndVariance(
824 const vigra::MultiArrayView<DIM,PIXEL_TYPE_IN> & inArray,
825 const double sigma,
826 vigra::MultiArrayView<DIM,PIXEL_TYPE_OUT> & meanArray,
827 vigra::MultiArrayView<DIM,PIXEL_TYPE_OUT> & varArray
828){
829 vigra::MultiArray<DIM,PIXEL_TYPE_OUT> tmpArray(inArray.shape());
830 gaussianMeanAndVariance<DIM,PIXEL_TYPE_IN,PIXEL_TYPE_OUT>(inArray,sigma,meanArray,varArray,tmpArray);
831}
832
833 namespace detail_non_local_means{
834
835 template<int DIM, class PIXEL_TYPE_IN,class PIXEL_TYPE_OUT,class SMOOTH_POLICY>
836 void nonLocalMean1Run(
837 const vigra::MultiArrayView<DIM,PIXEL_TYPE_IN> & image,
838 const SMOOTH_POLICY & smoothPolicy,
839 const NonLocalMeanParameter param,
840 vigra::MultiArrayView<DIM,PIXEL_TYPE_OUT> outImage
841){
842
843 typedef PIXEL_TYPE_IN PixelTypeIn;
844 typedef typename vigra::NumericTraits<PixelTypeIn>::RealPromote RealPromotePixelType;
845 typedef typename vigra::NumericTraits<RealPromotePixelType>::ValueType RealPromoteScalarType;
846 typedef SMOOTH_POLICY SmoothPolicyType;
847
848 typedef BlockWiseNonLocalMeanThreadObject<DIM,PixelTypeIn,SmoothPolicyType> ThreadObjectType;
849
850
851 // inspect parameter
852 vigra_precondition(param.stepSize_>=1,"NonLocalMean Parameter: \"stepSize>=1\" violated");
853 vigra_precondition(param.searchRadius_>=1, "NonLocalMean Parameter: \"searchRadius >=1\" violated");
854 vigra_precondition(param.patchRadius_>=1,"NonLocalMean Parameter: \"searchRadius >=1\" violated");
855 vigra_precondition(param.stepSize_-1<=param.patchRadius_,"NonLocalMean Parameter: \"stepSize -1 <= patchRadius\" violated");
856
857 // allocate arrays
858 vigra::MultiArray<DIM,RealPromotePixelType> meanImage(image.shape());
859 vigra::MultiArray<DIM,RealPromotePixelType> varImage(image.shape());
860 vigra::MultiArray<DIM,RealPromotePixelType> estimageImage(image.shape());
861 vigra::MultiArray<DIM,RealPromoteScalarType> labelImage(image.shape());
862
863 // compute mean and variance
864 // last argument is a "buffer" since within "gaussianMeanAndVariance" another array is needed
865 // ==> to avoid an unnecessary allocation we use the estimageImage as a buffer
866 //gaussianMeanAndVariance<DIM,PixelTypeIn,RealPromotePixelType>(image,param.sigmaMean_,meanImage,varImage,estimageImage);
867 gaussianMeanAndVariance<DIM,PixelTypeIn,RealPromotePixelType>(image,param.sigmaMean_,meanImage,varImage);
868
869 // initialize
870 labelImage = RealPromoteScalarType(0.0);
871 estimageImage = RealPromotePixelType(0.0);
872
873 ///////////////////////////////////////////////////////////////
874 { // MULTI THREAD CODE STARTS HERE
875
876
877
878 typedef threading::thread ThreadType;
879 typedef threading::mutex MutexType;
880
881 MutexType estimateMutex;
882 //typedef boost::thread ThreadType;
883
884 const size_t nThreads = param.nThreads_;
885 MultiArray<1,int> progress = MultiArray<1,int>(typename MultiArray<1,int>::difference_type(nThreads));
886
887 // allocate all thread objects
888 // each thread object works on a portion of the data
889 std::vector<ThreadObjectType> threadObjects(nThreads,
890 ThreadObjectType(image, meanImage, varImage, estimageImage, labelImage,
891 smoothPolicy, param, nThreads, estimateMutex,progress)
892 );
893
894 // thread ptr
895 std::vector<ThreadType *> threadPtrs(nThreads);
896 for(size_t i=0; i<nThreads; ++i){
897 ThreadObjectType & threadObj = threadObjects[i];
898 threadObj.setThreadIndex(i);
899 typename ThreadObjectType::RangeType lastAxisRange;
900 lastAxisRange[0]=(i * image.shape(DIM-1)) / nThreads;
901 lastAxisRange[1]=((i+1) * image.shape(DIM-1)) / nThreads;
902 threadObj.setRange(lastAxisRange);
903 // this will start the threads and cal operator()
904 // of the threadObjects
905 threadPtrs[i] = new ThreadType(threadObjects[i]);
906 }
907 for(size_t i=0; i<nThreads; ++i)
908 threadPtrs[i]->join();
909 for(size_t i=0; i<nThreads; ++i)
910 delete threadPtrs[i];
911
912 } // MULTI THREAD CODE ENDS HERE
913 ///////////////////////////////////////////////////////////////
914
915 // normalize estimates by the number of labels
916 // and write that in output
917 for(int scanOrderIndex=0; scanOrderIndex<labelImage.size(); ++scanOrderIndex){
918 if (labelImage[scanOrderIndex] <= RealPromoteScalarType(0.00001))
919 outImage[scanOrderIndex]=image[scanOrderIndex];
920 else
921 outImage[scanOrderIndex]=estimageImage[scanOrderIndex] / labelImage[scanOrderIndex];
922 }
923}
924
925}
926
927 template<int DIM, class PIXEL_TYPE_IN,class PIXEL_TYPE_OUT,class SMOOTH_POLICY>
928 void nonLocalMean(
929 const vigra::MultiArrayView<DIM,PIXEL_TYPE_IN> & image,
930 const SMOOTH_POLICY & smoothPolicy,
931 const NonLocalMeanParameter param,
932 vigra::MultiArrayView<DIM,PIXEL_TYPE_OUT> outImage
933){
934 detail_non_local_means::nonLocalMean1Run<DIM,PIXEL_TYPE_IN,PIXEL_TYPE_OUT,SMOOTH_POLICY>(image,smoothPolicy,param,outImage);
935 if(param.iterations_>1){
936
937 vigra::MultiArray<DIM,PIXEL_TYPE_OUT> tmp(outImage.shape());
938 for(size_t i=0;i<static_cast<size_t>(param.iterations_-1);++i){
939 tmp=outImage;
940 detail_non_local_means::nonLocalMean1Run<DIM,PIXEL_TYPE_OUT,PIXEL_TYPE_OUT,SMOOTH_POLICY>(tmp,smoothPolicy,param,outImage);
941 }
942 }
943}
944
945
946
947
948} // end namespace vigra
949
950
951#endif
vigra::MultiArray::difference_type
view_type::difference_type difference_type
Definition multi_array.hxx:2524
vigra::RGBValue
Class for a single RGB value.
Definition rgbvalue.hxx:128
vigra::RGBValue::RGBValue
RGBValue()
Definition rgbvalue.hxx:209
vigra::TinyVectorBase::size
size_type size() const
Definition tinyvector.hxx:913
vigra::TinyVectorBase::end
iterator end()
Definition tinyvector.hxx:864
vigra
vigra::norm
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition fftw3.hxx:1037
vigra::gaussianSmoothMultiArray
void gaussianSmoothMultiArray(...)
Isotropic Gaussian smoothing of a multi-dimensional arrays.
vigra::labelImage
unsigned int labelImage(...)
Find the connected components of a segmented image.
vigra::sum
NumericTraits< V >::Promote sum(TinyVectorBase< V, SIZE, D1, D2 > const &l)
sum of the vector's elements
Definition tinyvector.hxx:2073
vigra::clipLower
TinyVector< V, SIZE > clipLower(TinyVector< V, SIZE > const &t)
Clip negative values.
Definition tinyvector.hxx:2268
vigra::squaredNorm
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition fftw3.hxx:1044

© 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.12.2