52 #ifndef VIGRA_NON_LOCAL_MEAN 53 #define VIGRA_NON_LOCAL_MEAN 60 #include "multi_array.hxx" 61 #include "multi_convolution.hxx" 63 #include "threading.hxx" 64 #include "gaussians.hxx" 68 struct NonLocalMeanParameter{
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 80 sigmaSpatial_(sigmaSpatial),
81 searchRadius_(searchRadius),
82 patchRadius_(patchRadius),
83 sigmaMean_(sigmaMean),
85 iterations_(iterations),
103 class RatioPolicyParameter{
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
112 meanRatio_(meanRatio),
123 template<
class PIXEL_TYPE_IN>
127 typedef RatioPolicyParameter ParameterType;
128 typedef typename NumericTraits<PIXEL_TYPE_IN>::RealPromote PixelType;
129 typedef typename NumericTraits<PIXEL_TYPE_IN>::ValueType ValueType;
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_){
140 bool usePixel(
const PixelType & meanA,
const PixelType & varA)
const{
141 return sum(meanA) > epsilon_ &&
sum(varA) > epsilon_;
146 const PixelType & meanA,
const PixelType & varA,
147 const PixelType & meanB,
const PixelType & varB
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_));
155 ValueType distanceToWeight(
const PixelType & meanA,
const PixelType & varA,
const ValueType distance){
156 return exp(-distance /sigmaSquared_);
160 ValueType meanRatio_;
163 ValueType sigmaSquared_;
171 class 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
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)
202 template<
int DIM,
bool ALWAYS_INSIDE>
206 struct BorderHelper<DIM,true>{
207 template<
class COORD,
class IMAGE>
208 static bool isInside(
const COORD & ,
const IMAGE & ){
211 template<
class COORD,
class IMAGE>
212 static bool isOutside(
const COORD & ,
const IMAGE & ){
216 template<
class COORD,
class IMAGE>
217 static void mirrorIfIsOutsidePoint(COORD & ,IMAGE & ){
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);
228 template<
class COORD,
class IMAGE>
229 static bool isOutside(
const COORD & c,
const IMAGE & img){
230 return img.isOutside(c);
233 template<
class COORD,
class IMAGE>
234 static void mirrorIfIsOutsidePoint(COORD & coord,
const IMAGE & img){
235 for(
int c=0;c<DIM;++c){
237 coord[c]=-1*coord[c];
238 else if(coord[c]>= img.shape(c))
239 coord[c] = 2 * img.shape(c) - coord[c] - 1;
249 template<
class PIXEL_TYPE_IN>
253 typedef NormPolicyParameter ParameterType;
254 typedef typename NumericTraits<PIXEL_TYPE_IN>::RealPromote PixelType;
255 typedef typename NumericTraits<PIXEL_TYPE_IN>::ValueType ValueType;
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_){
266 bool usePixel(
const PixelType & meanA,
const PixelType & varA)
const{
267 return sum(varA)>epsilon_;
272 const PixelType & meanA,
const PixelType & varA,
273 const PixelType & meanB,
const PixelType & varB
277 const ValueType v = mean(varA/varB);
279 return (m < meanDist_ && v > varRatio_ && v < (1.0 / varRatio_));
282 ValueType distanceToWeight(
const PixelType & meanA,
const PixelType & varA,
const ValueType distance){
283 return exp(-distance /sigmaSquared_);
290 ValueType sigmaSquared_;
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;
303 typedef NonLocalMeanParameter ParameterType;
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;
316 typedef TinyVector<int,2> RangeType;
319 typedef threading::thread ThreadType;
320 typedef threading::mutex MutexType;
322 BlockWiseNonLocalMeanThreadObject(
323 const InArrayView & inImage,
324 MeanArrayView & meanImage,
325 VarArrayView & varImage,
326 EstimateArrayView & estimageImage,
328 const SmoothPolicyType & smoothPolicy,
329 const ParameterType & param,
330 const size_t nThreads,
331 MutexType & estimateMutex,
332 MultiArray<1,int> & progress
336 meanImage_(meanImage),
338 estimageImage_(estimageImage),
339 labelImage_(labelImage),
340 smoothPolicy_(smoothPolicy),
345 estimateMutexPtr_(&estimateMutex),
347 average_(
std::pow( (double)(2*param.patchRadius_+1), DIM) ),
348 gaussWeight_(
std::pow( (double)(2*param.patchRadius_+1), DIM) ),
349 shape_(inImage.shape()),
353 for(
int dim=0;dim<DIM;++dim)
354 totalSize_*=(shape_[dim]/param.stepSize_);
357 void setRange(
const RangeType & lastAxisRange){
358 lastAxisRange_=lastAxisRange;
360 void setThreadIndex(
const size_t threadIndex){
361 threadIndex_=threadIndex;
369 template<
bool ALWAYS_INSIDE>
370 void processSinglePixel(
const Coordinate & xyz);
372 template<
bool ALWAYS_INSIDE>
373 void processSinglePair(
const Coordinate & xyz,
const Coordinate & nxyz,RealPromoteScalarType & wmax,RealPromoteScalarType & totalweight);
375 template<
bool ALWAYS_INSIDE>
376 RealPromoteScalarType patchDistance(
const Coordinate & xyz,
const Coordinate & nxyz);
378 template<
bool ALWAYS_INSIDE>
379 void patchExtractAndAcc(
const Coordinate & xyz,
const RealPromoteScalarType weight);
381 template<
bool ALWAYS_INSIDE>
382 void patchAccMeanToEstimate(
const Coordinate & xyz,
const RealPromoteScalarType globalSum);
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);
392 void initalizeGauss();
394 void progressPrinter(
const int counter);
398 InArrayView inImage_;
399 MeanArrayView meanImage_;
400 VarArrayView varImage_;
401 EstimateArrayView estimageImage_;
402 LabelArrayView labelImage_;
405 SmoothPolicyType smoothPolicy_;
408 ParameterType param_;
411 RangeType lastAxisRange_;
414 MutexType * estimateMutexPtr_;
415 MultiArrayView<1,int> progress_;
419 BlockAverageVectorType average_;
420 BlockGaussWeightVectorType gaussWeight_;
426 template<
int DIM,
class PIXEL_TYPE_IN,
class SMOOTH_POLICY>
427 inline void BlockWiseNonLocalMeanThreadObject<DIM, PIXEL_TYPE_IN, SMOOTH_POLICY>::initalizeGauss(){
429 const int pr = param_.patchRadius_;
430 Gaussian<RealPromoteScalarType> gaussian(param_.sigmaSpatial_);
432 RealPromoteScalarType
sum = RealPromoteScalarType(0.0);
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);
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);
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);
464 for(
size_t i=0;i<gaussWeight_.size();++i){
465 gaussWeight_[i]/=
sum;
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){
476 for(
size_t ti=0;ti<nThreads_;++ti)
481 std::cout<<
"\rprogress "<<std::setw(10)<<pr<<
" %%"<<std::flush;
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_;
492 this->initalizeGauss();
496 if(param_.verbose_ && threadIndex_==nThreads_-1){
497 std::cout<<
"progress";
501 for (xyz[1] = start; xyz[1] < end; xyz[1] += stepSize)
502 for (xyz[0] = 0; xyz[0] < shape_[0]; xyz[0] += stepSize){
504 if(isAlwaysInside(xyz))
505 this->processSinglePixel<true>(xyz);
507 this->processSinglePixel<false>(xyz);
509 this->progressPrinter(c);
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);
520 this->processSinglePixel<false>(xyz);
522 this->progressPrinter(c);
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);
534 this->processSinglePixel<false>(xyz);
536 this->progressPrinter(c);
540 if(param_.verbose_ && threadIndex_==nThreads_-1){
541 std::cout<<
"\rprogress "<<std::setw(10)<<
"100"<<
" %%"<<
"\n";
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
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;
557 if(smoothPolicy_.usePixel(meanImage_[xyz],varImage_[xyz])){
559 RealPromoteScalarType wmax = 0.0;
560 const Coordinate start = xyz- Coordinate(param_.searchRadius_);
561 const Coordinate end = xyz+ Coordinate(param_.searchRadius_);
564 for (nxyz[1] = start[1]; nxyz[1] <= end[1]; nxyz[1]++)
565 for (nxyz[0] = start[0]; nxyz[0] <= end[0]; nxyz[0]++){
569 this->processSinglePair<ALWAYS_INSIDE>(xyz,nxyz,wmax,totalweight);
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]++){
578 this->processSinglePair<ALWAYS_INSIDE>(xyz,nxyz,wmax,totalweight);
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]++){
588 this->processSinglePair<ALWAYS_INSIDE>(xyz,nxyz,wmax,totalweight);
598 this->patchExtractAndAcc<ALWAYS_INSIDE>(xyz,wmax);
602 if (totalweight != 0.0){
603 this->patchAccMeanToEstimate<ALWAYS_INSIDE>(xyz,totalweight);
608 const double wmax = 1.0;
609 this->patchExtractAndAcc<ALWAYS_INSIDE>(xyz,wmax);
611 this->patchAccMeanToEstimate<ALWAYS_INSIDE>(xyz,totalweight);
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
626 if(BorderHelper<DIM,ALWAYS_INSIDE>::isInside(nxyz,inImage_)){
628 if(smoothPolicy_.usePixel(meanImage_[nxyz],varImage_[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);
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
654 const int f = param_.patchRadius_;
655 Coordinate offset(SkipInitialization),nPa(SkipInitialization),nPb(SkipInitialization);
657 RealPromoteScalarType distancetotal = 0;
661 #define VIGRA_NLM_IN_LOOP_CODE \ 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); \ 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;
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;
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;
693 #undef VIGRA_NLM_IN_LOOP_CODE 694 return distancetotal / acu;
699 template<
int DIM,
class PIXEL_TYPE_IN,
class SMOOTH_POLICY>
700 template<
bool ALWAYS_INSIDE>
702 BlockWiseNonLocalMeanThreadObject<DIM,PIXEL_TYPE_IN,SMOOTH_POLICY>::patchExtractAndAcc(
703 const Coordinate & xyz,
704 const RealPromoteScalarType weight
706 Coordinate xyzPos(SkipInitialization),abc(SkipInitialization);
707 Coordinate nhSize3(param_.patchRadius_);
708 const int ns = 2 * param_.patchRadius_ + 1;
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; \ 718 average_[count] += inImage_[xyzPos]* weight; \ 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;
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;
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;
743 #undef VIGRA_NLM_IN_LOOP_CODE 747 template<
int DIM,
class PIXEL_TYPE_IN,
class SMOOTH_POLICY>
748 template<
bool ALWAYS_INSIDE>
750 BlockWiseNonLocalMeanThreadObject<DIM,PIXEL_TYPE_IN,SMOOTH_POLICY>::patchAccMeanToEstimate(
751 const Coordinate & xyz,
752 const RealPromoteScalarType globalSum
754 Coordinate abc(SkipInitialization),xyzPos(SkipInitialization),nhSize(param_.patchRadius_);
756 const int ns = 2 * param_.patchRadius_ + 1;
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); \ 767 estimageImage_[xyzPos] = value; \ 768 labelImage_[xyzPos]+=gw; \ 769 estimateMutexPtr_->unlock(); \ 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;
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;
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;
794 #undef VIGRA_NLM_IN_LOOP_CODE 799 template<
int DIM,
class PIXEL_TYPE_IN,
class SMOOTH_POLICY,
class PIXEL_TYPE_OUT>
800 inline void gaussianMeanAndVariance(
811 for(
int scanOrderIndex=0;scanOrderIndex<inArray.
size();++scanOrderIndex){
812 tmpArray[scanOrderIndex]=
vigra::pow(inArray[scanOrderIndex],2);
815 for(
int scanOrderIndex=0;scanOrderIndex<inArray.
size();++scanOrderIndex){
816 PIXEL_TYPE_OUT var = varArray[scanOrderIndex] -
vigra::pow(meanArray[scanOrderIndex],2);
819 varArray[scanOrderIndex] = var;
823 template<
int DIM,
class PIXEL_TYPE_IN,
class SMOOTH_POLICY,
class PIXEL_TYPE_OUT>
824 inline void gaussianMeanAndVariance(
831 gaussianMeanAndVariance<DIM,PIXEL_TYPE_IN,PIXEL_TYPE_OUT>(inArray,sigma,meanArray,varArray,tmpArray);
834 namespace detail_non_local_means{
836 template<
int DIM,
class PIXEL_TYPE_IN,
class PIXEL_TYPE_OUT,
class SMOOTH_POLICY>
837 void nonLocalMean1Run(
839 const SMOOTH_POLICY & smoothPolicy,
840 const NonLocalMeanParameter param,
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;
849 typedef BlockWiseNonLocalMeanThreadObject<DIM,PixelTypeIn,SmoothPolicyType> ThreadObjectType;
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");
868 gaussianMeanAndVariance<DIM,PixelTypeIn,RealPromotePixelType>(image,param.sigmaMean_,meanImage,varImage);
871 labelImage = RealPromoteScalarType(0.0);
872 estimageImage = RealPromotePixelType(0.0);
879 typedef threading::thread ThreadType;
880 typedef threading::mutex MutexType;
882 MutexType estimateMutex;
885 const size_t nThreads = param.nThreads_;
890 std::vector<ThreadObjectType> threadObjects(nThreads,
891 ThreadObjectType(image, meanImage, varImage, estimageImage, labelImage,
892 smoothPolicy, param, nThreads, estimateMutex,progress)
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);
906 threadPtrs[i] =
new ThreadType(threadObjects[i]);
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];
918 for(
int scanOrderIndex=0; scanOrderIndex<labelImage.size(); ++scanOrderIndex){
919 if (labelImage[scanOrderIndex] <= RealPromoteScalarType(0.00001))
920 outImage[scanOrderIndex]=image[scanOrderIndex];
922 outImage[scanOrderIndex]=estimageImage[scanOrderIndex] / labelImage[scanOrderIndex];
928 template<
int DIM,
class PIXEL_TYPE_IN,
class PIXEL_TYPE_OUT,
class SMOOTH_POLICY>
931 const SMOOTH_POLICY & smoothPolicy,
932 const NonLocalMeanParameter param,
935 detail_non_local_means::nonLocalMean1Run<DIM,PIXEL_TYPE_IN,PIXEL_TYPE_OUT,SMOOTH_POLICY>(image,smoothPolicy,param,outImage);
936 if(param.iterations_>1){
939 for(
size_t i=0;i<param.iterations_-1;++i){
941 detail_non_local_means::nonLocalMean1Run<DIM,PIXEL_TYPE_OUT,PIXEL_TYPE_OUT,SMOOTH_POLICY>(tmp,smoothPolicy,param,outImage);
unsigned int labelImage(...)
Find the connected components of a segmented image.
const difference_type & shape() const
Definition: multi_array.hxx:1596
linalg::TemporaryMatrix< T > exp(MultiArrayView< 2, T, C > const &v)
void gaussianSmoothMultiArray(...)
Isotropic Gaussian smoothing of a multi-dimensional arrays.
Definition: array_vector.hxx:954
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition: fftw3.hxx:1044
Main MultiArray class containing the memory management.
Definition: multi_array.hxx:2422
Definition: accessor.hxx:43
view_type::difference_type difference_type
Definition: multi_array.hxx:2470
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition: fftw3.hxx:1037
difference_type_1 size() const
Definition: multi_array.hxx:1589
NumericTraits< V >::Promote sum(TinyVectorBase< V, SIZE, D1, D2 > const &l)
sum of the vector's elements
Definition: tinyvector.hxx:2073
TinyVector< V, SIZE > pow(TinyVectorBase< V, SIZE, D1, D2 > const &v, E exponent)
Definition: tinyvector.hxx:2036
Base class for, and view to, vigra::MultiArray.
Definition: multi_array.hxx:652
TinyVector< V, SIZE > clipLower(TinyVector< V, SIZE > const &t)
Clip negative values.
Definition: tinyvector.hxx:2268