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

non_local_mean.hxx VIGRA

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
647 BlockWiseNonLocalMeanThreadObject<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
701 BlockWiseNonLocalMeanThreadObject<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
749 BlockWiseNonLocalMeanThreadObject<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::labelImage
unsigned int labelImage(...)
Find the connected components of a segmented image.
vigra::MultiArrayView::shape
const difference_type & shape() const
Definition: multi_array.hxx:1648
vigra::linalg::exp
linalg::TemporaryMatrix< T > exp(MultiArrayView< 2, T, C > const &v)
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:2474
vigra::MultiArray< DIM, int >::difference_type
view_type::difference_type difference_type
Definition: multi_array.hxx:2522
vigra::norm
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition: fftw3.hxx:1037
vigra::MultiArrayView::size
difference_type_1 size() const
Definition: multi_array.hxx:1641
vigra::sum
NumericTraits< V >::Promote sum(TinyVectorBase< V, SIZE, D1, D2 > const &l)
sum of the vector'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:704
vigra::clipLower
TinyVector< V, SIZE > clipLower(TinyVector< V, SIZE > const &t)
Clip negative values.
Definition: tinyvector.hxx:2268
vigra::gaussianSmoothMultiArray
void gaussianSmoothMultiArray(...)
Isotropic Gaussian smoothing of a multi-dimensional arrays.

© 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.1 (Fri May 19 2017)