diff --git a/rtengine/FTblockDN.cc b/rtengine/FTblockDN.cc index 71c65cfb0..9f74d9f2c 100644 --- a/rtengine/FTblockDN.cc +++ b/rtengine/FTblockDN.cc @@ -22,7 +22,7 @@ // //////////////////////////////////////////////////////////////// -#include +#include #include #include "../rtgui/threadutils.h" #include "rtengine.h" @@ -36,6 +36,7 @@ #include "sleef.c" #include "opthelper.h" #include "cplx_wavelet_dec.h" +#include "median.h" #ifdef _OPENMP #include @@ -48,6 +49,8 @@ #define epsilon 0.001f/(TS*TS) //tolerance +#define PIX_SORT(a,b) { if ((a)>(b)) {temp=(a);(a)=(b);(b)=temp;} } + #define med2(a0,a1,a2,a3,a4,median) { \ pp[0]=a0; pp[1]=a1; pp[2]=a2; pp[3]=a3; pp[4]=a4; \ PIX_SORT(pp[0],pp[1]) ; PIX_SORT(pp[3],pp[4]) ; PIX_SORT(pp[0],pp[3]) ;\ @@ -292,7 +295,6 @@ void ImProcFunctions::Median_Denoise( float **src, float **dst, const int width, medianOut[i][j] = medianIn[i][j]; } } else if(medianType == MED_3X3STRONG) { - float pp[9], temp; int j; for (j = 0; j < border; j++) { @@ -300,7 +302,7 @@ void ImProcFunctions::Median_Denoise( float **src, float **dst, const int width, } for (; j < width - border; j++) { - med3(medianIn[i][j] , medianIn[i - 1][j], medianIn[i + 1][j] , medianIn[i][j + 1], medianIn[i][j - 1], medianIn[i - 1][j - 1], medianIn[i - 1][j + 1], medianIn[i + 1][j - 1], medianIn[i + 1][j + 1], medianOut[i][j]); + medianOut[i][j] = median(medianIn[i][j] , medianIn[i - 1][j], medianIn[i + 1][j] , medianIn[i][j + 1], medianIn[i][j - 1], medianIn[i - 1][j - 1], medianIn[i - 1][j + 1], medianIn[i + 1][j - 1], medianIn[i + 1][j + 1]); } for(; j < width; j++) { @@ -1859,7 +1861,7 @@ SSEFUNCTION void ImProcFunctions::RGB_denoise(int kall, Imagefloat * src, Imagef } else for (int j = 1; j < wid - 1; j++) { - med3(source->r(i, j), source->r(i - 1, j), source->r(i + 1, j), source->r(i, j + 1), source->r(i, j - 1), source->r(i - 1, j - 1), source->r(i - 1, j + 1), source->r(i + 1, j - 1), source->r(i + 1, j + 1), tm[i][j]); //3x3 + tm[i][j] = median(source->r(i, j), source->r(i - 1, j), source->r(i + 1, j), source->r(i, j + 1), source->r(i, j - 1), source->r(i - 1, j - 1), source->r(i - 1, j + 1), source->r(i + 1, j - 1), source->r(i + 1, j + 1)); //3x3 } } } else { @@ -1922,7 +1924,7 @@ SSEFUNCTION void ImProcFunctions::RGB_denoise(int kall, Imagefloat * src, Imagef } else for (int j = 1; j < wid - 1; j++) { - med3(source->b(i, j), source->b(i - 1, j), source->b(i + 1, j), source->b(i, j + 1), source->b(i, j - 1), source->b(i - 1, j - 1), source->b(i - 1, j + 1), source->b(i + 1, j - 1), source->b(i + 1, j + 1), tm[i][j]); + tm[i][j] = median(source->b(i, j), source->b(i - 1, j), source->b(i + 1, j), source->b(i, j + 1), source->b(i, j - 1), source->b(i - 1, j - 1), source->b(i - 1, j + 1), source->b(i + 1, j - 1), source->b(i + 1, j + 1)); } } } else { @@ -1987,7 +1989,7 @@ SSEFUNCTION void ImProcFunctions::RGB_denoise(int kall, Imagefloat * src, Imagef } else for (int j = 1; j < wid - 1; j++) { - med3(source->g(i, j), source->g(i - 1, j), source->g(i + 1, j), source->g(i, j + 1), source->g(i, j - 1), source->g(i - 1, j - 1), source->g(i - 1, j + 1), source->g(i + 1, j - 1), source->g(i + 1, j + 1), tm[i][j]); + tm[i][j] = median(source->g(i, j), source->g(i - 1, j), source->g(i + 1, j), source->g(i, j + 1), source->g(i, j - 1), source->g(i - 1, j - 1), source->g(i - 1, j + 1), source->g(i + 1, j - 1), source->g(i + 1, j + 1)); } } } else { diff --git a/rtengine/PF_correct_RT.cc b/rtengine/PF_correct_RT.cc index fce4a14c0..e64c750bb 100644 --- a/rtengine/PF_correct_RT.cc +++ b/rtengine/PF_correct_RT.cc @@ -31,6 +31,7 @@ #include "../rtgui/myflatcurve.h" #include "rt_math.h" #include "opthelper.h" +#include "median.h" #ifdef _OPENMP #include @@ -766,7 +767,6 @@ SSEFUNCTION void ImProcFunctions::Badpixelscam(CieImage * src, CieImage * dst, d #pragma omp parallel { int ip, in, jp, jn; - float pp[9], temp; #pragma omp for nowait //nowait because next loop inside this parallel region is independent on this one for (int i = 0; i < height; i++) { @@ -795,7 +795,7 @@ SSEFUNCTION void ImProcFunctions::Badpixelscam(CieImage * src, CieImage * dst, d jn = j + 2; } - med3(sraa[ip][jp], sraa[ip][j], sraa[ip][jn], sraa[i][jp], sraa[i][j], sraa[i][jn], sraa[in][jp], sraa[in][j], sraa[in][jn], tmaa[i][j]); + tmaa[i][j] = median(sraa[ip][jp], sraa[ip][j], sraa[ip][jn], sraa[i][jp], sraa[i][j], sraa[i][jn], sraa[in][jp], sraa[in][j], sraa[in][jn]); } } @@ -827,7 +827,7 @@ SSEFUNCTION void ImProcFunctions::Badpixelscam(CieImage * src, CieImage * dst, d jn = j + 2; } - med3(srbb[ip][jp], srbb[ip][j], srbb[ip][jn], srbb[i][jp], srbb[i][j], srbb[i][jn], srbb[in][jp], srbb[in][j], srbb[in][jn], tmbb[i][j]); + tmbb[i][j] = median(srbb[ip][jp], srbb[ip][j], srbb[ip][jn], srbb[i][jp], srbb[i][j], srbb[i][jn], srbb[in][jp], srbb[in][j], srbb[in][jn]); } } } @@ -1374,7 +1374,6 @@ SSEFUNCTION void ImProcFunctions::BadpixelsLab(LabImage * src, LabImage * dst, d #pragma omp parallel { int ip, in, jp, jn; - float pp[9], temp; #pragma omp for nowait //nowait because next loop inside this parallel region is independent on this one for (int i = 0; i < height; i++) { @@ -1403,7 +1402,7 @@ SSEFUNCTION void ImProcFunctions::BadpixelsLab(LabImage * src, LabImage * dst, d jn = j + 2; } - med3(sraa[ip][jp], sraa[ip][j], sraa[ip][jn], sraa[i][jp], sraa[i][j], sraa[i][jn], sraa[in][jp], sraa[in][j], sraa[in][jn], tmaa[i][j]); + tmaa[i][j] = median(sraa[ip][jp], sraa[ip][j], sraa[ip][jn], sraa[i][jp], sraa[i][j], sraa[i][jn], sraa[in][jp], sraa[in][j], sraa[in][jn]); } } @@ -1435,7 +1434,7 @@ SSEFUNCTION void ImProcFunctions::BadpixelsLab(LabImage * src, LabImage * dst, d jn = j + 2; } - med3(srbb[ip][jp], srbb[ip][j], srbb[ip][jn], srbb[i][jp], srbb[i][j], srbb[i][jn], srbb[in][jp], srbb[in][j], srbb[in][jn], tmbb[i][j]); + tmbb[i][j] = median(srbb[ip][jp], srbb[ip][j], srbb[ip][jn], srbb[i][jp], srbb[i][j], srbb[i][jn], srbb[in][jp], srbb[in][j], srbb[in][jn], tmbb[i][j]); } } } diff --git a/rtengine/dirpyrLab_denoise.cc b/rtengine/dirpyrLab_denoise.cc deleted file mode 100644 index ee434d75f..000000000 --- a/rtengine/dirpyrLab_denoise.cc +++ /dev/null @@ -1,751 +0,0 @@ -/* - * This file is part of RawTherapee. - * - * RawTherapee is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * RawTherapee is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with RawTherapee. If not, see . - * - * � 2010 Emil Martinec - * - */ - -#include -#include -#include "curves.h" -#include "labimage.h" -#include "improcfun.h" -#include "array2D.h" -#include "rt_math.h" - -#ifdef _OPENMP -#include -#endif - -#define CLIPC(a) ((a)>-32000?((a)<32000?(a):32000):-32000) - -#define DIRWT_L(i1,j1,i,j) ( rangefn_L[(data_fine->L[i1][j1]-data_fine->L[i][j]+32768)] ) - -#define DIRWT_AB(i1,j1,i,j) ( rangefn_ab[(data_fine->a[i1][j1]-data_fine->a[i][j]+32768)] * \ -rangefn_ab[(data_fine->L[i1][j1]-data_fine->L[i][j]+32768)] * \ -rangefn_ab[(data_fine->b[i1][j1]-data_fine->b[i][j]+32768)] ) - -//#define NRWT_L(a) (nrwt_l[a] ) - -#define NRWT_AB (nrwt_ab[(hipass[1]+32768)] * nrwt_ab[(hipass[2]+32768)]) - - -#define med3(a,b,c) (a(b)) {temp=(a);(a)=(b);(b)=temp;} } - -#define med3x3(a0,a1,a2,a3,a4,a5,a6,a7,a8,median) { \ -p[0]=a0; p[1]=a1; p[2]=a2; p[3]=a3; p[4]=a4; p[5]=a5; p[6]=a6; p[7]=a7; p[8]=a8; \ -PIX_SORT(p[1],p[2]); PIX_SORT(p[4],p[5]); PIX_SORT(p[7],p[8]); \ -PIX_SORT(p[0],p[1]); PIX_SORT(p[3],p[4]); PIX_SORT(p[6],p[7]); \ -PIX_SORT(p[1],p[2]); PIX_SORT(p[4],p[5]); PIX_SORT(p[7],p[8]); \ -PIX_SORT(p[0],p[3]); PIX_SORT(p[5],p[8]); PIX_SORT(p[4],p[7]); \ -PIX_SORT(p[3],p[6]); PIX_SORT(p[1],p[4]); PIX_SORT(p[2],p[5]); \ -PIX_SORT(p[4],p[7]); PIX_SORT(p[4],p[2]); PIX_SORT(p[6],p[4]); \ -PIX_SORT(p[4],p[2]); median=p[4];} //a4 is the median - - -namespace rtengine -{ - -static const int maxlevel = 4; - -//sequence of scales -//static const int scales[8] = {1,2,4,8,16,32,64,128}; -//sequence of pitches -//static const int pitches[8] = {1,1,1,1,1,1,1,1}; - -//sequence of scales -//static const int scales[8] = {1,1,1,1,1,1,1,1}; -//sequence of pitches -//static const int pitches[8] = {2,2,2,2,2,2,2,2}; - -//sequence of scales -//static const int scales[8] = {1,1,2,2,4,4,8,8}; -//sequence of pitches -//static const int pitches[8] = {2,1,2,1,2,1,2,1}; - -//sequence of scales -static const int scales[8] = {1, 1, 2, 4, 8, 16, 32, 64}; -//sequence of pitches -static const int pitches[8] = {2, 1, 1, 1, 1, 1, 1, 1}; - -//pitch is spacing of subsampling -//scale is spacing of directional averaging weights -//example 1: no subsampling at any level -- pitch=1, scale=2^n -//example 2: subsampling by 2 every level -- pitch=2, scale=1 at each level -//example 3: no subsampling at first level, subsampling by 2 thereafter -- -// pitch =1, scale=1 at first level; pitch=2, scale=2 thereafter - - - - -void ImProcFunctions :: dirpyrLab_denoise(LabImage * src, LabImage * dst, const procparams::DirPyrDenoiseParams & dnparams ) -{ - float gam = dnparams.gamma / 3.0; - //float gam = 2.0;//min(3.0, 0.1*fabs(c[4])/3.0+0.001); - float gamthresh = 0.03; - float gamslope = exp(log((double)gamthresh) / gam) / gamthresh; - - LUTf gamcurve(65536, 0); - - //DiagonalCurve* lumacurve = new DiagonalCurve (dnparams.lumcurve, CURVES_MIN_POLY_POINTS); - //DiagonalCurve* chromacurve = new DiagonalCurve (dnparams.chromcurve, CURVES_MIN_POLY_POINTS); - //LUTf Lcurve(65536); - //LUTf abcurve(65536); - for (int i = 0; i < 65536; i++) { - int g = (int)(CurveFactory::gamma((double)i / 65535.0, gam, gamthresh, gamslope, 1.0, 0.0) * 65535.0); - gamcurve[i] = CLIP(g); - /*float val = (float)i/65535.0; - float Lval = (2*(lumacurve->getVal(val))); - float abval = (2*(chromacurve->getVal(val))); - - Lcurve[i] = SQR(Lval); - abcurve[i] = SQR(abval); - if (i % 1000 ==0) printf("%d Lmult=%f abmult=%f \n",i,Lcurve[i],abcurve[i]);*/ - } - - //delete lumacurve; - //delete chromacurve; - - - - //#pragma omp parallel for if (multiThread) - for (int i = 0; i < src->H; i++) { - for (int j = 0; j < src->W; j++) { - //src->L[i][j] = CurveFactory::flinterp(gamcurve,src->L[i][j]); - src->L[i][j] = gamcurve[src->L[i][j]]; - } - } - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - - - LUTf rangefn_L(65536); - LUTf nrwt_l(65536); - - LUTf rangefn_ab(65536); - LUTf nrwt_ab(65536); - - //set up NR weight functions - - //gamma correction for chroma in shadows - float nrwtl_norm = ((CurveFactory::gamma((double)65535.0 / 65535.0, gam, gamthresh, gamslope, 1.0, 0.0)) - - (CurveFactory::gamma((double)75535.0 / 65535.0, gam, gamthresh, gamslope, 1.0, 0.0))); - - for (int i = 0; i < 65536; i++) { - nrwt_l[i] = ((CurveFactory::gamma((double)i / 65535.0, gam, gamthresh, gamslope, 1.0, 0.0) - - CurveFactory::gamma((double)(i + 10000) / 65535.0, gam, gamthresh, gamslope, 1.0, 0.0)) ) / nrwtl_norm; - //if (i % 100 ==0) printf("%d %f \n",i,nrwt_l[i]); - } - - float tonefactor = nrwt_l[32768]; - - float noise_L = 10.0 * dnparams.luma; - float noisevar_L = SQR(noise_L); - - float noise_ab = 100.0 * dnparams.chroma; - float noisevar_ab = SQR(noise_ab); - - - //set up range functions - for (int i = 0; i < 65536; i++) { - rangefn_L[i] = (( exp(-(double)fabs(i - 32768) * tonefactor / (1.0 + noise_L)) * (1.0 + noisevar_L) / ((double)(i - 32768) * (double)(i - 32768) + noisevar_L + 1.0))); - } - - for (int i = 0; i < 65536; i++) { - rangefn_ab[i] = (( exp(-(double)fabs(i - 32768) * tonefactor / (1.0 + 3 * noise_ab)) * (1.0 + noisevar_ab) / ((double)(i - 32768) * (double)(i - 32768) + noisevar_ab + 1.0))); - } - - - for (int i = 0; i < 65536; i++) { - nrwt_ab[i] = ((1.0 + abs(i - 32768) / (1.0 + 8 * noise_ab)) * exp(-(double)fabs(i - 32768) / (1.0 + 8 * noise_ab) ) ); - } - - - //for (int i=0; i<65536; i+=100) printf("%d %d \n",i,gamcurve[i]); - - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - - - int level; - - LabImage * dirpyrLablo[maxlevel]; - int w = (int)((src->W - 1) / pitches[0]) + 1; - int h = (int)((src->H - 1) / pitches[0]) + 1; - dirpyrLablo[0] = new LabImage(w, h); - - for (level = 1; level < maxlevel; level++) { - w = (int)((w - 1) / pitches[level]) + 1; - h = (int)((h - 1) / pitches[level]) + 1; - dirpyrLablo[level] = new LabImage(w, h); - }; - - - ////////////////////////////////////////////////////////////////////////////// - - - // c[0] = luma = noise_L - // c[1] = chroma = noise_ab - // c[2] decrease of noise var with scale - // c[3] radius of domain blur at each level - // c[4] shadow smoothing - // c[5] edge preservation - - level = 0; - - int scale = scales[level]; - - int pitch = pitches[level]; - - //int thresh = 10 * c[8]; - //impulse_nr (src, src, m_w1, m_h1, thresh, noisevar); - - dirpyr(src, dirpyrLablo[0], 0, rangefn_L, rangefn_ab, pitch, scale, dnparams.luma, dnparams.chroma ); - - level = 1; - - while(level < maxlevel) { - scale = scales[level]; - pitch = pitches[level]; - - dirpyr(dirpyrLablo[level - 1], dirpyrLablo[level], level, rangefn_L, rangefn_ab, pitch, scale, dnparams.luma, dnparams.chroma ); - - level ++; - } - - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - - - for(int level = maxlevel - 1; level > 0; level--) { - - int scale = scales[level]; - int pitch = pitches[level]; - idirpyr(dirpyrLablo[level], dirpyrLablo[level - 1], level, rangefn_L, nrwt_l, nrwt_ab, pitch, scale, dnparams.luma, dnparams.chroma/*, Lcurve, abcurve*/ ); - } - - - scale = scales[0]; - pitch = pitches[0]; - - // freeing as much memory as possible since the next call to idirpyr will need lots - for(int i = 1; i < maxlevel; i++) { - delete dirpyrLablo[i]; - } - - idirpyr(dirpyrLablo[0], dst, 0, rangefn_L, nrwt_l, nrwt_ab, pitch, scale, dnparams.luma, dnparams.chroma/*, Lcurve, abcurve*/ ); - - // freeing the last bunch of memory - delete dirpyrLablo[0]; - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - - float igam = 1 / gam; - float igamthresh = gamthresh * gamslope; - float igamslope = 1 / gamslope; - - for (int i = 0; i < 65536; i++) { - gamcurve[i] = (CurveFactory::gamma((float)i / 65535.0, igam, igamthresh, igamslope, 1.0, 0.0) * 65535.0); - } - - - if (dnparams.luma > 0) { - for (int i = 0; i < dst->H; i++) - for (int j = 0; j < dst->W; j++) { - dst->L[i][j] = gamcurve[dst->L[i][j]]; - } - } else { - for (int i = 0; i < dst->H; i++) - for (int j = 0; j < dst->W; j++) { - dst->L[i][j] = gamcurve[src->L[i][j]]; - } - } - -} - -void ImProcFunctions::dirpyr(LabImage* data_fine, LabImage* data_coarse, int level, - LUTf & rangefn_L, LUTf & rangefn_ab, int pitch, int scale, - const int luma, const int chroma ) -{ - - //pitch is spacing of subsampling - //scale is spacing of directional averaging weights - //example 1: no subsampling at any level -- pitch=1, scale=2^n - //example 2: subsampling by 2 every level -- pitch=2, scale=1 at each level - //example 3: no subsampling at first level, subsampling by 2 thereafter -- - // pitch =1, scale=1 at first level; pitch=2, scale=2 thereafter - - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - // calculate weights, compute directionally weighted average - - int width = data_fine->W; - int height = data_fine->H; - - //generate domain kernel - int halfwin = 3;//min(ceil(2*sig),3); - int scalewin = halfwin * scale; - -#ifdef _OPENMP - #pragma omp parallel for -#endif - - for(int i = 0; i < height; i += pitch ) { - int i1 = i / pitch; - - for(int j = 0, j1 = 0; j < width; j += pitch, j1++) { - float dirwt_l, dirwt_ab, norm_l, norm_ab; - //float lops,aops,bops; - float Lout, aout, bout; - norm_l = norm_ab = 0;//if we do want to include the input pixel in the sum - Lout = 0; - aout = 0; - bout = 0; - - for(int inbr = (i - scalewin); inbr <= (i + scalewin); inbr += scale) { - if (inbr < 0 || inbr > height - 1) { - continue; - } - - for (int jnbr = (j - scalewin); jnbr <= (j + scalewin); jnbr += scale) { - if (jnbr < 0 || jnbr > width - 1) { - continue; - } - - dirwt_l = DIRWT_L(inbr, jnbr, i, j); - dirwt_ab = DIRWT_AB(inbr, jnbr, i, j); - Lout += dirwt_l * data_fine->L[inbr][jnbr]; - aout += dirwt_ab * data_fine->a[inbr][jnbr]; - bout += dirwt_ab * data_fine->b[inbr][jnbr]; - norm_l += dirwt_l; - norm_ab += dirwt_ab; - } - } - - //lops = Lout/norm;//diagnostic - //aops = aout/normab;//diagnostic - //bops = bout/normab;//diagnostic - - data_coarse->L[i1][j1] = Lout / norm_l; //low pass filter - data_coarse->a[i1][j1] = aout / norm_ab; - data_coarse->b[i1][j1] = bout / norm_ab; - - - /*if (level<2 && i>0 && i0 && jL[i-1][j-1], data_fine->L[i-1][j], data_fine->L[i-1][j+1], \ - data_fine->L[i][j-1], data_fine->L[i][j], data_fine->L[i][j+1], \ - data_fine->L[i+1][j-1], data_fine->L[i+1][j], data_fine->L[i+1][j+1]); - //med3x3(data_fine->L[i-1][j-1], data_fine->L[i-1][j], data_fine->L[i-1][j+1], \ - data_fine->L[i][j-1], data_fine->L[i][j], data_fine->L[i][j+1], \ - data_fine->L[i+1][j-1], data_fine->L[i+1][j], data_fine->L[i+1][j+1],Lmed); - - data_coarse->L[i1][j1] = Lhmf; - }*/ - } - } - - - - -} - -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -void ImProcFunctions::idirpyr(LabImage* data_coarse, LabImage* data_fine, int level, LUTf &rangefn_L, LUTf & nrwt_l, LUTf & nrwt_ab, - int pitch, int scale, const int luma, const int chroma/*, LUTf & Lcurve, LUTf & abcurve*/ ) -{ - - int width = data_fine->W; - int height = data_fine->H; - - array2D nrfactorL (width, height); - - //float eps = 0.0; - - // c[0] noise_L - // c[1] noise_ab (relative to noise_L) - // c[2] decrease of noise var with scale - // c[3] radius of domain blur at each level - // c[4] shadow smoothing - - float noisevar_L = 4 * SQR(25.0 * luma); - float noisevar_ab = 2 * SQR(100.0 * chroma); - float scalefactor = 1.0 / pow(2.0, (level + 1) * 2); //change the last 2 to 1 for longer tail of higher scale NR - - noisevar_L *= scalefactor; - - // for coarsest level, take non-subsampled lopass image and subtract from lopass_fine to generate hipass image - - // denoise hipass image, add back into lopass_fine to generate denoised image at fine scale - - // now iterate: - // (1) take denoised image at level n, expand and smooth using gradient weights from lopass image at level n-1 - // the result is the smoothed image at level n-1 - // (2) subtract smoothed image at level n-1 from lopass image at level n-1 to make hipass image at level n-1 - // (3) denoise the hipass image at level n-1 - // (4) add the denoised image at level n-1 to the smoothed image at level n-1 to make the denoised image at level n-1 - - // note that the coarsest level amounts to skipping step (1) and doing (2,3,4). - // in other words, skip step one if pitch=1 - - // step (1) - - if (pitch == 1) { - - // step (1-2-3-4) - -#ifdef _OPENMP - #pragma omp parallel -#endif - { - -#ifdef _OPENMP - #pragma omp for -#endif - - for(int i = 0; i < height; i++) - for(int j = 0; j < width; j++) - { - float hipass[3], hpffluct[3], tonefactor, nrfactor; - - tonefactor = (nrwt_l[data_coarse->L[i][j]]); - - hipass[1] = data_fine->a[i][j] - data_coarse->a[i][j]; - hipass[2] = data_fine->b[i][j] - data_coarse->b[i][j]; - - //Wiener filter - //luma - if (level < 2) { - hipass[0] = data_fine->L[i][j] - data_coarse->L[i][j]; - hpffluct[0] = SQR(hipass[0]) + SQR(hipass[1]) + SQR(hipass[2]) + 0.001; - nrfactorL[i][j] = (1.0 + hpffluct[0]) / (1.0 + hpffluct[0] + noisevar_L /* * Lcurve[data_coarse->L[i][j]]*/); - //hipass[0] *= hpffluct[0]/(hpffluct[0]+noisevar_L); - //data_fine->L[i][j] = CLIP(hipass[0]+data_coarse->L[i][j]); - } - - //chroma - //hipass[1] = data_fine->a[i][j]-data_coarse->a[i][j]; - //hipass[2] = data_fine->b[i][j]-data_coarse->b[i][j]; - hpffluct[1] = SQR(hipass[1] * tonefactor) + 0.001; - hpffluct[2] = SQR(hipass[2] * tonefactor) + 0.001; - nrfactor = (hpffluct[1] + hpffluct[2]) / ((hpffluct[1] + hpffluct[2]) + noisevar_ab * NRWT_AB); - - hipass[1] *= nrfactor; - hipass[2] *= nrfactor; - - data_fine->a[i][j] = hipass[1] + data_coarse->a[i][j]; - data_fine->b[i][j] = hipass[2] + data_coarse->b[i][j]; - } - - if (level < 2) - { -#ifdef _OPENMP - #pragma omp for -#endif - - for(int i = 0; i < height; i++) - for(int j = 0; j < width; j++) { - - float dirwt_l, norm_l; - float nrfctrave = 0; - norm_l = 0;//if we do want to include the input pixel in the sum - - for(int inbr = max(0, i - 1); inbr <= min(height - 1, i + 1); inbr++) { - for (int jnbr = max(0, j - 1); jnbr <= min(width - 1, j + 1); jnbr++) { - dirwt_l = DIRWT_L(inbr, jnbr, i, j); - nrfctrave += dirwt_l * nrfactorL[inbr][jnbr]; - norm_l += dirwt_l; - } - } - - nrfctrave /= norm_l; - //nrfctrave = nrfactorL[i][j]; - //nrfctrave=1; - - float hipass[3]; - - //luma - - /*if (i>0 && i0 && jL[i][j] - data_coarse->L[i][j]); - //hipass[0] = median*(data_fine->L[i][j]-data_coarse->L[i][j]); - //hipass[0] = nrfactorL[i][j]*(data_fine->L[i][j]-data_coarse->L[i][j]); - data_fine->L[i][j] = CLIP(hipass[0] + data_coarse->L[i][j]); - - //chroma - //hipass[1] = nrfactorab[i][j]*(data_fine->a[i][j]-data_coarse->a[i][j]); - //hipass[2] = nrfactorab[i][j]*(data_fine->b[i][j]-data_coarse->b[i][j]); - - //data_fine->a[i][j] = hipass[1]+data_coarse->a[i][j]; - //data_fine->b[i][j] = hipass[2]+data_coarse->b[i][j]; - } - }//end of luminance correction - - - - }//end of pitch=1 - - } else {//pitch>1 - - LabImage* smooth; - - smooth = new LabImage(width, height); -#ifdef _OPENMP - #pragma omp parallel -#endif - - { - -#ifdef _OPENMP - #pragma omp for -#endif - - for(int i = 0; i < height; i += pitch) { - int ix = i / pitch; - - for(int j = 0, jx = 0; j < width; j += pitch, jx++) { - - //copy common pixels - smooth->L[i][j] = data_coarse->L[ix][jx]; - smooth->a[i][j] = data_coarse->a[ix][jx]; - smooth->b[i][j] = data_coarse->b[ix][jx]; - } - } - - //if (pitch>1) {//pitch=2; step (1) expand coarse image, fill in missing data -#ifdef _OPENMP - #pragma omp for -#endif - - for(int i = 0; i < height - 1; i += 2) - for(int j = 0; j < width - 1; j += 2) { - //do midpoint first - double norm = 0.0, wtdsum[3] = {0.0, 0.0, 0.0}; - - //wtdsum[0]=wtdsum[1]=wtdsum[2]=0.0; - for(int ix = i; ix < min(height, i + 3); ix += 2) - for (int jx = j; jx < min(width, j + 3); jx += 2) { - wtdsum[0] += smooth->L[ix][jx]; - wtdsum[1] += smooth->a[ix][jx]; - wtdsum[2] += smooth->b[ix][jx]; - norm++; - } - - norm = 1 / norm; - smooth->L[i + 1][j + 1] = wtdsum[0] * norm; - smooth->a[i + 1][j + 1] = wtdsum[1] * norm; - smooth->b[i + 1][j + 1] = wtdsum[2] * norm; - } - -#ifdef _OPENMP - #pragma omp for -#endif - - for(int i = 0; i < height - 1; i += 2) - for(int j = 0; j < width - 1; j += 2) { - //now right neighbor - if (j + 1 == width) { - continue; - } - - double norm = 0.0, wtdsum[3] = {0.0, 0.0, 0.0}; - - for (int jx = j; jx < min(width, j + 3); jx += 2) { - wtdsum[0] += smooth->L[i][jx]; - wtdsum[1] += smooth->a[i][jx]; - wtdsum[2] += smooth->b[i][jx]; - norm++; - } - - for (int ix = max(0, i - 1); ix < min(height, i + 2); ix += 2) { - wtdsum[0] += smooth->L[ix][j + 1]; - wtdsum[1] += smooth->a[ix][j + 1]; - wtdsum[2] += smooth->b[ix][j + 1]; - norm++; - } - - norm = 1 / norm; - smooth->L[i][j + 1] = wtdsum[0] * norm; - smooth->a[i][j + 1] = wtdsum[1] * norm; - smooth->b[i][j + 1] = wtdsum[2] * norm; - - //now down neighbor - if (i + 1 == height) { - continue; - } - - norm = 0.0; - wtdsum[0] = wtdsum[1] = wtdsum[2] = 0.0; - - for (int ix = i; ix < min(height, i + 3); ix += 2) { - wtdsum[0] += smooth->L[ix][j]; - wtdsum[1] += smooth->a[ix][j]; - wtdsum[2] += smooth->b[ix][j]; - norm++; - } - - for (int jx = max(0, j - 1); jx < min(width, j + 2); jx += 2) { - wtdsum[0] += smooth->L[i + 1][jx]; - wtdsum[1] += smooth->a[i + 1][jx]; - wtdsum[2] += smooth->b[i + 1][jx]; - norm++; - } - - norm = 1 / norm; - smooth->L[i + 1][j] = wtdsum[0] * norm; - smooth->a[i + 1][j] = wtdsum[1] * norm; - smooth->b[i + 1][j] = wtdsum[2] * norm; - - } - -#ifdef _OPENMP - #pragma omp for -#endif - - // step (2-3-4) - for( int i = 0; i < height; i++) - for(int j = 0; j < width; j++) { - - float tonefactor = (nrwt_l[smooth->L[i][j]]); - //double wtdsum[3], norm; - float hipass[3], hpffluct[3], nrfactor; - - hipass[1] = data_fine->a[i][j] - smooth->a[i][j]; - hipass[2] = data_fine->b[i][j] - smooth->b[i][j]; - - //Wiener filter - //luma - if (level < 2) { - hipass[0] = data_fine->L[i][j] - smooth->L[i][j]; - hpffluct[0] = SQR(hipass[0]) + SQR(hipass[1]) + SQR(hipass[2]) + 0.001; - nrfactorL[i][j] = (1.0 + hpffluct[0]) / (1.0 + hpffluct[0] + noisevar_L /* * Lcurve[smooth->L[i][j]]*/); - //hipass[0] *= hpffluct[0]/(hpffluct[0]+noisevar_L); - //data_fine->L[i][j] = CLIP(hipass[0]+smooth->L[i][j]); - } - - //chroma - //hipass[1] = data_fine->a[i][j]-smooth->a[i][j]; - //hipass[2] = data_fine->b[i][j]-smooth->b[i][j]; - hpffluct[1] = SQR(hipass[1] * tonefactor) + 0.001; - hpffluct[2] = SQR(hipass[2] * tonefactor) + 0.001; - nrfactor = (hpffluct[1] + hpffluct[2]) / ((hpffluct[1] + hpffluct[2]) + noisevar_ab * NRWT_AB /* * abcurve[smooth->L[i][j]]*/); - - hipass[1] *= nrfactor; - hipass[2] *= nrfactor; - - data_fine->a[i][j] = hipass[1] + smooth->a[i][j]; - data_fine->b[i][j] = hipass[2] + smooth->b[i][j]; - } - - - if (level < 2) { -#ifdef _OPENMP - #pragma omp for -#endif - - for(int i = 0; i < height; i++) - for(int j = 0; j < width; j++) { - - float dirwt_l, norm_l; - float nrfctrave = 0; - norm_l = 0;//if we do want to include the input pixel in the sum - - for(int inbr = (i - pitch); inbr <= (i + pitch); inbr += pitch) { - if (inbr < 0 || inbr > height - 1) { - continue; - } - - for (int jnbr = (j - pitch); jnbr <= (j + pitch); jnbr += pitch) { - if (jnbr < 0 || jnbr > width - 1) { - continue; - } - - dirwt_l = DIRWT_L(inbr, jnbr, i, j); - nrfctrave += dirwt_l * nrfactorL[inbr][jnbr]; - norm_l += dirwt_l; - } - } - - nrfctrave /= norm_l; - //nrfctrave = nrfactorL[i][j]; - //nrfctrave=1; - - - float hipass[3]; - - //luma - - /*if (i>0 && i0 && jL[i][j] - smooth->L[i][j]); - //hipass[0] = median*(data_fine->L[i][j]-smooth->L[i][j]); - //hipass[0] = nrfactorL[i][j]*(data_fine->L[i][j]-data_coarse->L[i][j]); - data_fine->L[i][j] = CLIP(hipass[0] + smooth->L[i][j]); - - - //chroma - //hipass[1] = nrfactorab[i][j]*(data_fine->a[i][j]-data_coarse->a[i][j]); - //hipass[2] = nrfactorab[i][j]*(data_fine->b[i][j]-data_coarse->b[i][j]); - - //data_fine->a[i][j] = hipass[1]+data_coarse->a[i][j]; - //data_fine->b[i][j] = hipass[2]+data_coarse->b[i][j]; - } - }//end of luminance correction - - - } // end parallel - delete smooth; - }//end of pitch>1 - -} - - -#undef DIRWT_L -#undef DIRWT_AB - -//#undef NRWT_L -#undef NRWT_AB - -} - diff --git a/rtengine/improcfun.h b/rtengine/improcfun.h index ccc77ab9f..4d03e2747 100644 --- a/rtengine/improcfun.h +++ b/rtengine/improcfun.h @@ -35,20 +35,6 @@ #include "cplx_wavelet_dec.h" #include "pipettebuffer.h" -#define PIX_SORT(a,b) { if ((a)>(b)) {temp=(a);(a)=(b);(b)=temp;} } - -#define med3(a0,a1,a2,a3,a4,a5,a6,a7,a8,median) { \ -pp[0]=a0; pp[1]=a1; pp[2]=a2; pp[3]=a3; pp[4]=a4; pp[5]=a5; pp[6]=a6; pp[7]=a7; pp[8]=a8; \ -PIX_SORT(pp[1],pp[2]); PIX_SORT(pp[4],pp[5]); PIX_SORT(pp[7],pp[8]); \ -PIX_SORT(pp[0],pp[1]); PIX_SORT(pp[3],pp[4]); PIX_SORT(pp[6],pp[7]); \ -PIX_SORT(pp[1],pp[2]); PIX_SORT(pp[4],pp[5]); PIX_SORT(pp[7],pp[8]); \ -PIX_SORT(pp[0],pp[3]); PIX_SORT(pp[5],pp[8]); PIX_SORT(pp[4],pp[7]); \ -PIX_SORT(pp[3],pp[6]); PIX_SORT(pp[1],pp[4]); PIX_SORT(pp[2],pp[5]); \ -PIX_SORT(pp[4],pp[7]); PIX_SORT(pp[4],pp[2]); PIX_SORT(pp[6],pp[4]); \ -PIX_SORT(pp[4],pp[2]); median=pp[4];} //pp4 = median - - - namespace rtengine { @@ -292,7 +278,6 @@ public: // pyramid denoise procparams::DirPyrDenoiseParams dnparams; - void dirpyrLab_denoise(LabImage * src, LabImage * dst, const procparams::DirPyrDenoiseParams & dnparams );//Emil's directional pyramid denoise void dirpyr (LabImage* data_fine, LabImage* data_coarse, int level, LUTf &rangefn_L, LUTf &rangefn_ab, int pitch, int scale, const int luma, int chroma ); void idirpyr (LabImage* data_coarse, LabImage* data_fine, int level, LUTf &rangefn_L, LUTf & nrwt_l, LUTf & nrwt_ab, diff --git a/rtengine/ipretinex.cc b/rtengine/ipretinex.cc index 6773e7a12..06f779875 100644 --- a/rtengine/ipretinex.cc +++ b/rtengine/ipretinex.cc @@ -36,30 +36,20 @@ */ -#include -#include -#include -#include +#include +#include +#include +#include #include "rtengine.h" #include "gauss.h" #include "rawimagesource.h" #include "improcfun.h" #include "opthelper.h" +#include "median.h" #include "StopWatch.h" #define clipretinex( val, minv, maxv ) (( val = (val < minv ? minv : val ) ) > maxv ? maxv : val ) -#define med3(a0,a1,a2,a3,a4,a5,a6,a7,a8,median) { \ -pp[0]=a0; pp[1]=a1; pp[2]=a2; pp[3]=a3; pp[4]=a4; pp[5]=a5; pp[6]=a6; pp[7]=a7; pp[8]=a8; \ -PIX_SORT(pp[1],pp[2]); PIX_SORT(pp[4],pp[5]); PIX_SORT(pp[7],pp[8]); \ -PIX_SORT(pp[0],pp[1]); PIX_SORT(pp[3],pp[4]); PIX_SORT(pp[6],pp[7]); \ -PIX_SORT(pp[1],pp[2]); PIX_SORT(pp[4],pp[5]); PIX_SORT(pp[7],pp[8]); \ -PIX_SORT(pp[0],pp[3]); PIX_SORT(pp[5],pp[8]); PIX_SORT(pp[4],pp[7]); \ -PIX_SORT(pp[3],pp[6]); PIX_SORT(pp[1],pp[4]); PIX_SORT(pp[2],pp[5]); \ -PIX_SORT(pp[4],pp[7]); PIX_SORT(pp[4],pp[2]); PIX_SORT(pp[6],pp[4]); \ -PIX_SORT(pp[4],pp[2]); median=pp[4];} //pp4 = median - - namespace { void retinex_scales( float* scales, int nscales, int mode, int s, float high) @@ -629,10 +619,8 @@ void RawImageSource::MSR(float** luminance, float** originalLuminance, float **e #endif for (int i = borderL; i < hei - borderL; i++) { - float pp[9], temp; - for (int j = borderL; j < wid - borderL; j++) { - med3(luminance[i][j], luminance[i - 1][j], luminance[i + 1][j], luminance[i][j + 1], luminance[i][j - 1], luminance[i - 1][j - 1], luminance[i - 1][j + 1], luminance[i + 1][j - 1], luminance[i + 1][j + 1], tmL[i][j]); //3x3 + tmL[i][j] = median(luminance[i][j], luminance[i - 1][j], luminance[i + 1][j], luminance[i][j + 1], luminance[i][j - 1], luminance[i - 1][j - 1], luminance[i - 1][j + 1], luminance[i + 1][j - 1], luminance[i + 1][j + 1]); //3x3 } } diff --git a/rtengine/median.h b/rtengine/median.h index d5e88d9de..800c3767f 100644 --- a/rtengine/median.h +++ b/rtengine/median.h @@ -16,26 +16,100 @@ * You should have received a copy of the GNU General Public License * along with RawTherapee. If not, see . */ -#include "rt_math.h" + +#pragma once + +#include +#include + +#include "opthelper.h" + +template +T median(std::array array) +{ + const typename std::array::iterator middle = array.begin() + array.size() / 2; + std::nth_element(array.begin(), middle, array.end()); + + return + array.size() % 2 + ? ((*middle + *std::min_element(middle + 1, array.end())) / static_cast(2)) + : *middle; +} + +template +T median(T arg, ARGS... args) +{ + return median(std::array{std::move(arg), std::move(args)...}); +} + +template +inline T median3(T a, T b, T c) +{ + return std::max(std::min(a, b), std::min(c, std::max(a, b))); +} + +template<> +inline vfloat median3(vfloat a, vfloat b, vfloat c) +{ + return vmaxf(vminf(a, b), vminf(c, vmaxf(a, b))); +} + +// See http://stackoverflow.com/questions/480960/code-to-calculate-median-of-five-in-c-sharp +template +inline T median5(T a, T b, T c, T d, T e) +{ + if (b < a) { + std::swap(a, b); + } + if (d < c) { + std::swap(c, d); + } + + if (c < a) { + std::swap(b, d); + c = a; + } + + a = e; + + if (b < a) { + std::swap(a, b); + } + + if (a < c) { + std::swap(b, d); + a = c; + } + + return std::min(a, d); +} + +template<> +inline vfloat median5(vfloat a, vfloat b, vfloat c, vfloat d, vfloat e) +{ + const vfloat f = vmaxf(vminf(a, b), vminf(c, d)); + const vfloat g = vminf(vmaxf(a, b), vmaxf(c, d)); + return median3(e, f, g); +} // middle 4 of 6 elements, #define MIDDLE4OF6(s0,s1,s2,s3,s4,s5,d0,d1,d2,d3,d4,d5,temp) \ {\ -d1 = min(s1,s2);\ -d2 = max(s1,s2);\ -d0 = min(s0,d2);\ -d2 = max(s0,d2);\ -temp = min(d0,d1);\ -d1 = max(d0,d1);\ +d1 = std::min(s1,s2);\ +d2 = std::max(s1,s2);\ +d0 = std::min(s0,d2);\ +d2 = std::max(s0,d2);\ +temp = std::min(d0,d1);\ +d1 = std::max(d0,d1);\ d0 = temp;\ -d4 = min(s4,s5);\ -d5 = max(s4,s5);\ -d3 = min(s3,d5);\ -d5 = max(s3,d5);\ -temp = min(d3,d4);\ -d4 = max(d3,d4);\ -d3 = max(d0,temp);\ -d2 = min(d2,d5);\ +d4 = std::min(s4,s5);\ +d5 = std::max(s4,s5);\ +d3 = std::min(s3,d5);\ +d5 = std::max(s3,d5);\ +temp = std::min(d3,d4);\ +d4 = std::max(d3,d4);\ +d3 = std::max(d0,temp);\ +d2 = std::min(d2,d5);\ } // middle 4 of 6 elements, vectorized @@ -61,27 +135,27 @@ d2 = vminf(d2,d5);\ #define MEDIAN7(s0,s1,s2,s3,s4,s5,s6,t0,t1,t2,t3,t4,t5,t6,median) \ {\ -t0 = min(s0,s5);\ -t5 = max(s0,s5);\ -t3 = max(t0,s3);\ -t0 = min(t0,s3);\ -t1 = min(s1,s6);\ -t6 = max(s1,s6);\ -t2 = min(s2,s4);\ -t4 = max(s2,s4);\ -t1 = max(t0,t1);\ -median = min(t3,t5);\ -t5 = max(t3,t5);\ +t0 = std::min(s0,s5);\ +t5 = std::max(s0,s5);\ +t3 = std::max(t0,s3);\ +t0 = std::min(t0,s3);\ +t1 = std::min(s1,s6);\ +t6 = std::max(s1,s6);\ +t2 = std::min(s2,s4);\ +t4 = std::max(s2,s4);\ +t1 = std::max(t0,t1);\ +median = std::min(t3,t5);\ +t5 = std::max(t3,t5);\ t3 = median;\ -median = min(t2,t6);\ -t6 = max(t2,t6);\ -t3 = max(median,t3);\ -t3 = min(t3,t6);\ -t4 = min(t4,t5);\ -median = min(t1,t4);\ -t4 = max(t1,t4);\ -t3 = max(median,t3);\ -median = min(t3,t4);\ +median = std::min(t2,t6);\ +t6 = std::max(t2,t6);\ +t3 = std::max(median,t3);\ +t3 = std::min(t3,t6);\ +t4 = std::min(t4,t5);\ +median = std::min(t1,t4);\ +t4 = std::max(t1,t4);\ +t3 = std::max(median,t3);\ +median = std::min(t3,t4);\ } #define VMEDIAN7(s0,s1,s2,s3,s4,s5,s6,t0,t1,t2,t3,t4,t5,t6,median) \ diff --git a/rtengine/minmax.h b/rtengine/minmax.h deleted file mode 100644 index 69e939b72..000000000 --- a/rtengine/minmax.h +++ /dev/null @@ -1,76 +0,0 @@ -/* - * This file is part of RawTherapee. - * - * Copyright (c) 2004-2010 Gabor Horvath - * - * RawTherapee is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * RawTherapee is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with RawTherapee. If not, see . - */ -#define MINMAX3(a,b,c,min,max) \ -{ \ -if ((a)<(b)) { \ - if ((b)<(c)) { \ - (min) = (a); \ - (max) = (c); \ - } \ - else { \ - (max) = (b); \ - if ((a)<(c)) \ - (min) = (a); \ - else \ - (min) = (c); \ - } \ -} else { \ - if ((b)>(c)) { \ - (min) = (c); \ - (max) = (a); \ - } \ - else { \ - (min) = (b); \ - if ((a)>(c)) \ - (max) = (a); \ - else \ - (max) = (c); \ - } \ -} \ -} - -#define MIN3(a,b,c,min) \ -{ \ -if ((a)<(b)) { \ - if ((a)<(c)) \ - (min) = (a); \ - else \ - (min) = (c); \ -} else { \ - if ((b)>(c)) \ - (min) = (c); \ - else \ - (min) = (b); \ -} \ -} - -#define MAX3(a,b,c,min) \ -{ \ -if ((a)>(b)) { \ - if ((a)>(c)) \ - (max) = (a); \ - else \ - (max) = (c); \ -} else { \ - if ((b)<(c)) \ - (max) = (c); \ - else \ - (max) = (b); \ -} \ -} diff --git a/rtengine/rawimagesource.cc b/rtengine/rawimagesource.cc index 218ba9fe2..7be671560 100644 --- a/rtengine/rawimagesource.cc +++ b/rtengine/rawimagesource.cc @@ -33,6 +33,7 @@ #include "dcp.h" #include "rt_math.h" #include "improcfun.h" +#include "median.h" #ifdef _OPENMP #include #endif @@ -433,13 +434,6 @@ PIX_SORT(p[3],p[6]); PIX_SORT(p[1],p[4]); PIX_SORT(p[2],p[5]); \ PIX_SORT(p[4],p[7]); PIX_SORT(p[4],p[2]); PIX_SORT(p[6],p[4]); \ PIX_SORT(p[4],p[2]); median=p[4];} //a4 is the median -#define med5(a0,a1,a2,a3,a4,median) { \ -p[0]=a0; p[1]=a1; p[2]=a2; p[3]=a3; p[4]=a4; \ -PIX_SORT(p[0],p[1]) ; PIX_SORT(p[3],p[4]) ; PIX_SORT(p[0],p[3]) ; \ -PIX_SORT(p[1],p[4]) ; PIX_SORT(p[1],p[2]) ; PIX_SORT(p[2],p[3]) ; \ -PIX_SORT(p[1],p[2]) ; median=p[2] ;} - - RawImageSource::RawImageSource () : ImageSource() , plistener(NULL) @@ -3018,7 +3012,7 @@ SSEFUNCTION void RawImageSource::cfaboxblur(RawImage *riFlatFile, float* cfablur for (int i = 0; i < H; i++) { int iprev, inext, jprev, jnext; - int p[5], temp, median; + int med; if (i < 2) { iprev = i + 2; @@ -3048,12 +3042,11 @@ SSEFUNCTION void RawImageSource::cfaboxblur(RawImage *riFlatFile, float* cfablur //med3x3(riFlatFile->data[iprev][jprev], riFlatFile->data[iprev][j], riFlatFile->data[iprev][jnext], // riFlatFile->data[i][jprev], riFlatFile->data[i][j], riFlatFile->data[i][jnext], // riFlatFile->data[inext][jprev], riFlatFile->data[inext][j], riFlatFile->data[inext][jnext], cfatmp[i*W+j]); - med5(riFlatFile->data[iprev][j], riFlatFile->data[i][jprev], riFlatFile->data[i][j], - riFlatFile->data[i][jnext], riFlatFile->data[inext][j], median); + med = median(riFlatFile->data[iprev][j], riFlatFile->data[i][jprev], riFlatFile->data[i][j], riFlatFile->data[i][jnext], riFlatFile->data[inext][j]); // if (riFlatFile->data[i][j]>hotdeadthresh*median || median>hotdeadthresh*riFlatFile->data[i][j]) { - if (((int)riFlatFile->data[i][j] << 1) > median || (median << 1) > riFlatFile->data[i][j]) { - cfatmp[i * W + j] = median; + if (((int)riFlatFile->data[i][j] << 1) > med || (med << 1) > riFlatFile->data[i][j]) { + cfatmp[i * W + j] = med; } else { cfatmp[i * W + j] = riFlatFile->data[i][j]; } diff --git a/rtengine/rt_math.h b/rtengine/rt_math.h index 898f1397c..702fb5360 100644 --- a/rtengine/rt_math.h +++ b/rtengine/rt_math.h @@ -11,10 +11,10 @@ static const float MAXVALF = float(MAXVAL); // float version of MAXVAL static const double MAXVALD = double(MAXVAL); // double version of MAXVAL template -inline const _Tp SQR (_Tp x) +inline _Tp SQR (_Tp x) { // return std::pow(x,2); Slower than: - return (x * x); + return x * x; } template @@ -31,13 +31,13 @@ inline const _Tp& max(const _Tp& a, const _Tp& b) template -inline const _Tp LIM(const _Tp& a, const _Tp& b, const _Tp& c) +inline const _Tp& LIM(const _Tp& a, const _Tp& b, const _Tp& c) { return std::max(b, std::min(a, c)); } template -inline const _Tp LIM01(const _Tp& a) +inline _Tp LIM01(const _Tp& a) { return std::max(_Tp(0), std::min(a, _Tp(1))); } @@ -49,7 +49,7 @@ inline const _Tp ULIM(const _Tp& a, const _Tp& b, const _Tp& c) } template -inline const _Tp CLIP(const _Tp& a) +inline _Tp CLIP(const _Tp& a) { return LIM(a, static_cast<_Tp>(0), static_cast<_Tp>(MAXVAL)); } @@ -80,7 +80,7 @@ inline const _Tp& max(const _Tp& a, const _Tp& b, const _Tp& c, const _Tp& d) } template -inline const _Tp intp(const _Tp a, const _Tp b, const _Tp c) +inline _Tp intp(_Tp a, _Tp b, _Tp c) { // calculate a * b + (1 - a) * c // following is valid: @@ -90,19 +90,19 @@ inline const _Tp intp(const _Tp a, const _Tp b, const _Tp c) } template -T norm1(const T& x, const T& y) +inline T norm1(const T& x, const T& y) { return std::abs(x) + std::abs(y); } template -T norm2(const T& x, const T& y) +inline T norm2(const T& x, const T& y) { return std::sqrt(x * x + y * y); } template< typename T > -T norminf(const T& x, const T& y) +inline T norminf(const T& x, const T& y) { return std::max(std::abs(x), std::abs(y)); }