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
-
-}
-