From 80bc1a6afc411191712548ad479b05c733749bc1 Mon Sep 17 00:00:00 2001 From: Emil Martinec Date: Mon, 19 Mar 2012 07:44:42 -0500 Subject: [PATCH] Removing dirpyrLab_denoise.cc (deprecated code); implementing fix in issue 1084 (now moved to color.cc). --- rtengine/color.cc | 12 +- rtengine/dirpyrLab_denoise.cc | 708 ---------------------------------- 2 files changed, 6 insertions(+), 714 deletions(-) delete mode 100644 rtengine/dirpyrLab_denoise.cc diff --git a/rtengine/color.cc b/rtengine/color.cc index 3701d47db..877a567f0 100644 --- a/rtengine/color.cc +++ b/rtengine/color.cc @@ -157,12 +157,12 @@ void Color::hsv2rgb (float h, float s, float v, float &r, float &g, float &b) { float r1,g1,b1; - if (i==0) {r1 = v; g1 = t; b1 = p;} - else if (i==1) {r1 = q; g1 = v; b1 = p;} - else if (i==2) {r1 = p; g1 = v; b1 = t;} - else if (i==3) {r1 = p; g1 = q; b1 = v;} - else if (i==4) {r1 = t; g1 = p; b1 = v;} - else if (i==5) {r1 = v; g1 = p; b1 = q;} + if (i==1) {r1 = q; g1 = v; b1 = p;} + else if (i==2) {r1 = p; g1 = v; b1 = t;} + else if (i==3) {r1 = p; g1 = q; b1 = v;} + else if (i==4) {r1 = t; g1 = p; b1 = v;} + else if (i==5) {r1 = v; g1 = p; b1 = q;} + else /*i==(0|6)*/ {r1 = v; g1 = t; b1 = p;} r = ((r1)*65535.0); g = ((g1)*65535.0); diff --git a/rtengine/dirpyrLab_denoise.cc b/rtengine/dirpyrLab_denoise.cc deleted file mode 100644 index dd79c60b8..000000000 --- a/rtengine/dirpyrLab_denoise.cc +++ /dev/null @@ -1,708 +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 "rtengine.h" -#include -#include -#include "curves.h" -#include "labimage.h" -#include "improcfun.h" -#include "array2D.h" - -#ifdef _OPENMP -#include -#endif - -#define SQR(x) ((x)*(x)) -#define CLIPTO(a,b,c) ((a)>(b)?((a)<(c)?(a):(c)):(b)) -#define CLIPC(a) ((a)>-32000?((a)<32000?(a):32000):-32000) -#define CLIP(a) (CLIPTO(a,0,65535)) - - -#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)(Color::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; iH; i++) { - for (int j=0; jW; j++) { - //src->L[i][j] = Color::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 = ((Color::gamma((double)65535.0/65535.0, gam, gamthresh, gamslope, 1.0, 0.0)) - - (Color::gamma((double)75535.0/65535.0, gam, gamthresh, gamslope, 1.0, 0.0))); - for (int i=0; i<65536; i++) { - nrwt_l[i] = ((Color::gamma((double)i/65535.0, gam, gamthresh, gamslope, 1.0, 0.0) - - Color::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 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] = (Color::gamma((float)i/65535.0, igam, igamthresh, igamslope, 1.0, 0.0) * 65535.0); - } - - - if (dnparams.luma>0) { - for (int i=0; iH; i++) - for (int j=0; jW; j++) { - dst->L[i][j] = gamcurve[dst->L[i][j]]; - } - } else { - for (int i=0; iH; i++) - for (int j=0; jW; 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; ixL[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; jxL[i][jx]; - wtdsum[1] += smooth->a[i][jx]; - wtdsum[2] += smooth->b[i][jx]; - norm++; - } - for (int ix=MAX(0,i-1); ixL[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; ixL[ix][j]; - wtdsum[1] += smooth->a[ix][j]; - wtdsum[2] += smooth->b[ix][j]; - norm++; - } - for (int jx=MAX(0,j-1); jxL[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 - -} -