diff --git a/rtengine/EdgePreserveLab.cc b/rtengine/EdgePreserveLab.cc index 6226fa0cd..53900bae0 100644 --- a/rtengine/EdgePreserveLab.cc +++ b/rtengine/EdgePreserveLab.cc @@ -1,4 +1,9 @@ #include "EdgePreserveLab.h" +#include "boxblur.h" +#include + +#define MAX(a,b) ((a)<(b)?(b):(a)) +#define MIN(a,b) ((a)>(b)?(b):(a)) //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -50,28 +55,33 @@ float *EdgePreserveLab::CreateBlur(float *Source, float LScale, float abScale, f unsigned int x, y, i; unsigned int w1 = w - 1, h1 = h - 1; - float eps = 0.02f; + float eps = 0.0001f; float scL = powf(100.0f,LScale); float scab = powf(200.0f,abScale); + + float * var = new float[w*h]; + boxvar(g, var, 1, 1, w, h); for(y = 0; y != h1; y++){ float *rg = &g[w*y]; for(x = 0; x != w1; x++){ //Estimate the central difference gradient in the center of a four pixel square. (gx, gy) is actually 2*gradient. - float gx = (fabs((rg[x + 1] - rg[x]) + (rg[x + w + 1] - rg[x + w]))); + /*float gx = (fabs((rg[x + 1] - rg[x]) + (rg[x + w + 1] - rg[x + w]))); float gy = (fabs((rg[x + w] - rg[x]) + (rg[x + w + 1] - rg[x + 1]))); - - float Lave = 0;//0.25*((rg[x + 1] + rg[x]) + (rg[x + w + 1] + rg[x + w])); - + + //TODO: combine this with gx, gy if not needing separate quantities float hx = (fabs((rg[x + 1 + n] - rg[x + n]) + (rg[x + w + 1 + n] - rg[x + w + n])) + \ fabs((rg[x + 1 + 2*n] - rg[x + 2*n]) + (rg[x + w + 1 + 2*n] - rg[x + w + 2*n]))); float hy = (fabs((rg[x + w + n] - rg[x + n]) + (rg[x + w + 1 + n] - rg[x + 1 + n])) + \ fabs((rg[x + w + 2*n] - rg[x + 2*n]) + (rg[x + w + 1 + 2*n] - rg[x + 1 + 2*n]))); - + */ + //float gradtot = (gx+gy+hx+hy); + //gradhisto[MAX(0,MIN(32767,(int)gradtot))] ++; //Apply power to the magnitude of the gradient to get the edge stopping function. - a[x + w*y] = scL*expf(-100.0f*(gx + gy + hx + hy)/(EdgeStoppingLuma)); - //a[x + w*y] = LScale*100.0f*expf(-100.0f*SQR(gx + gy)/SQR(EdgeStoppingLuma));///(0.1+rg[x]); + //a[x + w*y] = scL*expf(-100.0f*(gx + gy + hx + hy)/(EdgeStoppingLuma)); + //a[x + w*y] = scL*expf(-var[y*w+x]/SQR(0.02*EdgeStoppingLuma));///(0.1+rg[x]); + a[x + w*y] = scL*expf(-50.0f*sqrt(var[y*w+x])/(EdgeStoppingLuma+eps));///(0.1+rg[x]); //b[x + w*y] = (scab)*expf(-20.0f*(gx + gy + Lave*(hx + hy))/(EdgeStoppingChroma)); //b[x + w*y] = (scab)*expf(-400.0f*SQR(gx + gy + Lave*(hx + hy))/SQR(EdgeStoppingChroma));; diff --git a/rtengine/FTblockDNchroma.cc b/rtengine/FTblockDNchroma.cc index 07d062d80..36b276a8d 100644 --- a/rtengine/FTblockDNchroma.cc +++ b/rtengine/FTblockDNchroma.cc @@ -34,6 +34,7 @@ #include "LUT.h" #include "array2D.h" #include "iccmatrices.h" +#include "boxblur.h" #ifdef _OPENMP #include @@ -51,9 +52,9 @@ //#define ULIM(x,y,z) ((y) < (z) ? LIM(x,y,z) : LIM(x,z,y)) //#define CLIP(x) LIM(x,0,65535) -#define TS 256 // Tile size -#define offset TS/2 // shift between tiles -#define fTS ((TS/2+1)) // shift between Fourier tiles +#define TS 128 // Tile size +#define offset (TS/4) // shift between tiles +#define fTS ((TS/2+1)) // second dimension of Fourier tiles //#define eps 0.01f/(TS*TS) //tolerance namespace rtengine { @@ -79,301 +80,261 @@ namespace rtengine { //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -void ImProcFunctions::RGB_InputTransf(Imagefloat * src, LabImage * dst, LabImage * blur, const procparams::DirPyrDenoiseParams & dnparams, const procparams::DefringeParams & defringe) { - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - // gamma transform input channel data - float gam = dnparams.gamma; - float gamthresh = 0.03; - float gamslope = exp(log((double)gamthresh)/gam)/gamthresh; - float gam3slope = exp(log((double)gamthresh)/3.0f)/gamthresh; - - LUTf gamcurve(65536,0); - - for (int i=0; i<65536; i++) { - gamcurve[i] = (gamma((double)i/65535.0, gam, gamthresh, gamslope, 1.0, 0.0));// * 32768.0f; - } + void ImProcFunctions::RGB_InputTransf(Imagefloat * src, LabImage * dst, const procparams::DirPyrDenoiseParams & dnparams, const procparams::DefringeParams & defringe) { + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + // gamma transform input channel data + float gam = dnparams.gamma; + float gamthresh = 0.03; + float gamslope = exp(log((double)gamthresh)/gam)/gamthresh; + float gam3slope = exp(log((double)gamthresh)/3.0f)/gamthresh; + + LUTf gamcurve(65536,0); + + for (int i=0; i<65536; i++) { + gamcurve[i] = (gamma((double)i/65535.0, gam, gamthresh, gamslope, 1.0, 0.0)) * 32768.0f; + } + + //srand((unsigned)time(0));//test with random data + #ifdef _OPENMP #pragma omp parallel for #endif - for (int i=0; iheight; i++) - for (int j=0; jwidth; j++) { - - float X = xyz_prophoto[0][0]*src->r[i][j] + xyz_prophoto[0][1]*src->g[i][j] + xyz_prophoto[0][2]*src->b[i][j]; - float Y = xyz_prophoto[1][0]*src->r[i][j] + xyz_prophoto[1][1]*src->g[i][j] + xyz_prophoto[1][2]*src->b[i][j]; - float Z = xyz_prophoto[2][0]*src->r[i][j] + xyz_prophoto[2][1]*src->g[i][j] + xyz_prophoto[2][2]*src->b[i][j]; - - X = gamcurve[X]; - Y = gamcurve[Y]; - Z = gamcurve[Z]; - - dst->L[i][j] = Y; - dst->a[i][j] = 0.2f*(X-Y); - dst->b[i][j] = 0.2f*(Y-Z); - - } - - //pre-blur L channel to mitigate artifacts - /*for (int i=1; iheight-1; i++) - for (int j=1; jwidth-1; j++) { - dst->L[i][j] = 0.25*(dst->L[i][j]+0.5*(dst->L[i-1][j]+dst->L[i+1][j]+dst->L[i][j-1]+dst->L[i][j+1]) - +0.25*(dst->L[i-1][j-1]+dst->L[i-1][j+1]+dst->L[i+1][j-1]+dst->L[i+1][j+1])); - - }*/ + for (int i=0; iheight; i++) + for (int j=0; jwidth; j++) { + + float X = xyz_prophoto[0][0]*src->r[i][j] + xyz_prophoto[0][1]*src->g[i][j] + xyz_prophoto[0][2]*src->b[i][j]; + float Y = xyz_prophoto[1][0]*src->r[i][j] + xyz_prophoto[1][1]*src->g[i][j] + xyz_prophoto[1][2]*src->b[i][j]; + float Z = xyz_prophoto[2][0]*src->r[i][j] + xyz_prophoto[2][1]*src->g[i][j] + xyz_prophoto[2][2]*src->b[i][j]; + + X = X<65535.0f ? gamcurve[X] : (gamma((double)X/65535.0, gam, gamthresh, gamslope, 1.0, 0.0)*32768.0f); + Y = Y<65535.0f ? gamcurve[Y] : (gamma((double)Y/65535.0, gam, gamthresh, gamslope, 1.0, 0.0)*32768.0f); + Z = Z<65535.0f ? gamcurve[Z] : (gamma((double)Z/65535.0, gam, gamthresh, gamslope, 1.0, 0.0)*32768.0f); + + dst->L[i][j] = Y; + dst->a[i][j] = 0.2f*(X-Y); + dst->b[i][j] = 0.2f*(Y-Z); + + //Y = 0.05+0.1*((float)rand()/(float)RAND_MAX);//test with random data + //dst->L[i][j] = gamcurve[65535.0f*Y]; + + } - //cplx_wavelet_decomposition Ldecomp(dst->data, dst->W, dst->H, 2 /*maxlvl*/); - //Ldecomp.reconstruct(dst->data); - - impulse_nr (dst, 50.0f/20.0f); - //PF_correct_RT(dst, dst, defringe.radius, defringe.threshold); - - - //blur the image - //dirpyr_ab(dst, blur, dnparams);//use dirpyr here if using it to blur all channels - - EPDParams *p = (EPDParams *)(¶ms->edgePreservingDecompositionUI); - - float EdgeStop = (pow(0.15f, 1.0f/dnparams.gamma)/0.15f) * (float)dnparams.Lamt/10.0f;//compensate edgestopping power according to gamma - EdgePreserveLab epd = EdgePreserveLab(dst->W, dst->H); - blur->data = epd.CreateIteratedBlur(dst->data /*Source*/, (float)p->Scale /*LScale*/, (float)p->EdgeStopping /*abScale*/, \ - EdgeStop /*(float)dnparams.Lamt/10.0f*/ /*EdgeStoppingLuma*/, (float)dnparams.chroma/25.0f/*p->EdgeStopping*/ /*EdgeStoppingChroma*/, \ - 5/*Iterates*/, 0/*Reweightings*/, blur->data/*Blur*/); - - //impulsedenoise (blur); - -/* -#ifdef _OPENMP -#pragma omp parallel -#endif - { - AlignedBuffer* buffer = new AlignedBuffer (MAX(dst->W,dst->H)); - //gaussHorizontal (src->L, tmp1->L, buffer, src->W, src->H, radius, multiThread); - //gaussHorizontal (src->b, tmp1->b, buffer, src->W, src->H, radius, multiThread); - //gaussVertical (tmp1->a, tmp1->a, buffer, src->W, src->H, radius, multiThread); - //gaussVertical (tmp1->b, tmp1->b, buffer, src->W, src->H, radius, multiThread); - gaussHorizontal (dst->L, blur->L, buffer, dst->W, dst->H, (double)dnparams.Lamt/25.0f, multiThread); - gaussVertical (blur->L, blur->L, buffer, dst->W, dst->H, (double)dnparams.Lamt/25.0f, multiThread); - //gaussHorizontal (blur->L, blur->L, buffer, dst->W, dst->H, (double)dnparams.Lamt/25.0f, multiThread); - //gaussVertical (blur->L, blur->L, buffer, dst->W, dst->H, (double)dnparams.Lamt/25.0f, multiThread); - - delete buffer; } -*/ - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - // gamma transform input channel data - /*gam = dnparams.gamma/3.0f; - gamthresh = 0.03; - gamslope = exp(log((double)gamthresh)/gam)/gamthresh; - gam3slope = exp(log((double)gamthresh)/3.0f)/gamthresh; - - for (int i=0; i<65536; i++) { - gamcurve[i] = (gamma((double)i/65535.0, gam, gamthresh, gamslope, 1.0, 0.0)) * 65535.0f; - }*/ + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + void ImProcFunctions::RGB_OutputTransf(LabImage * src, Imagefloat * dst, const procparams::DirPyrDenoiseParams & dnparams) { + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + // gamma transform output channel data + float gam = dnparams.gamma; + float gamthresh = 0.03; + float gamslope = exp(log((double)gamthresh)/gam)/gamthresh; + float igam = 1/gam; + float igamthresh = gamthresh*gamslope; + float igamslope = 1/gamslope; + + float gamL = dnparams.gamma/3.0f; + float gamslopeL = exp(log((double)gamthresh)/gamL)/gamthresh; + float igamL = 1/gamL; + float igamthreshL = gamthresh*gamslopeL; + float igamslopeL = 1/gamslopeL; + + LUTf gamcurve(65536,0); + LUTf igamcurve(65536,0); + LUTf igamcurveL(65536,0); + + for (int i=0; i<65536; i++) { + + igamcurve[i] = (gamma((float)i/32768.0f, igam, igamthresh, igamslope, 1.0, 0.0) * 65535.0f); + + /*float L = (gamma((float)i/65535.0f, igamL, igamthreshL, igamslopeL, 1.0, 0.0) * 200.0f); + float fy = (0.00862069 * L) + 0.137932; // (L+16)/116 + float Y = f2xyz(fy); + igamcurveL[i] = (gamma(Y, gam, gamthresh, gamslope, 1.0, 0.0)) * 32768.0f;*/ + } + #ifdef _OPENMP #pragma omp parallel for #endif - for (int i=0; iheight; i++) { - for (int j=0; jwidth; j++) { - //float wt = expf(-100.0f*fabs(dst->L[i][j]-blur->L[i][j])/((float)dnparams.luma)); - //blur->L[i][j] = wt*dst->L[i][j] + (1-wt)*blur->L[i][j]; - //blur->L[i][j] = dst->L[i][j]; - blur->a[i][j] /*= 32768.0f*0.2f*(blur->a[i][j]-blur->L[i][j]);/*/*= 32768.0f; - blur->b[i][j] /*= 32768.0f*0.2f*(blur->L[i][j]-blur->b[i][j]);/*/*= 32768.0f; - blur->L[i][j] *= 32768.0f;//= gamcurve[32768.0f*blur->L[i][j]]; - - dst->L[i][j] *= 32768.0f;//= gamcurve[32768.0f*dst->L[i][j]]; - dst->a[i][j] *= 32768.0f; - dst->b[i][j] *= 32768.0f; - - } - } - - //dirpyr_ab(blur, blur, dnparams);//use dirpyr here if using it to blur ab channels only - //dirpyrLab_denoise(blur, blur, dnparams);//use dirpyr here if using it to blur ab channels only - -} - -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -void ImProcFunctions::RGB_OutputTransf(LabImage * src, Imagefloat * dst, const procparams::DirPyrDenoiseParams & dnparams) { - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - // gamma transform output channel data - float gam = dnparams.gamma; - float gamthresh = 0.03; - float gamslope = exp(log((double)gamthresh)/gam)/gamthresh; - float igam = 1/gam; - float igamthresh = gamthresh*gamslope; - float igamslope = 1/gamslope; - - float gamL = dnparams.gamma/3.0f; - float gamslopeL = exp(log((double)gamthresh)/gamL)/gamthresh; - float igamL = 1/gamL; - float igamthreshL = gamthresh*gamslopeL; - float igamslopeL = 1/gamslopeL; - - LUTf gamcurve(65536,0); - LUTf igamcurve(65536,0); - LUTf igamcurveL(65536,0); - - for (int i=0; i<65536; i++) { - - igamcurve[i] = (gamma((float)i/32768.0f, igam, igamthresh, igamslope, 1.0, 0.0) * 65535.0f); - - /*float L = (gamma((float)i/65535.0f, igamL, igamthreshL, igamslopeL, 1.0, 0.0) * 200.0f); - float fy = (0.00862069 * L) + 0.137932; // (L+16)/116 - float Y = f2xyz(fy); - igamcurveL[i] = (gamma(Y, gam, gamthresh, gamslope, 1.0, 0.0)) * 32768.0f;*/ - } - -#ifdef _OPENMP -#pragma omp parallel for -#endif - for (int i=0; iH; i++) { - float X,Y,Z; - for (int j=0; jW; j++) { - //input normalized to (0,1) - //Y = igamcurveL[ src->L[i][j] ]; - Y = src->L[i][j]; - X = (5.0f*(src->a[i][j])) + Y; - Z = Y - (5.0f*(src->b[i][j])); - - X = igamcurve[X]; - Y = igamcurve[Y]; - Z = igamcurve[Z]; - - dst->r[i][j] = prophoto_xyz[0][0]*X + prophoto_xyz[0][1]*Y + prophoto_xyz[0][2]*Z; - dst->g[i][j] = prophoto_xyz[1][0]*X + prophoto_xyz[1][1]*Y + prophoto_xyz[1][2]*Z; - dst->b[i][j] = prophoto_xyz[2][0]*X + prophoto_xyz[2][1]*Y + prophoto_xyz[2][2]*Z; - - } - } -} - -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -void ImProcFunctions::RGB_denoise(Imagefloat * src, Imagefloat * dst, /*int Roffset,*/ const procparams::DirPyrDenoiseParams & dnparams, const procparams::DefringeParams & defringe) -{ - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - /*if (plistener) { - plistener->setProgressStr ("Block FT Luma Denoise..."); - plistener->setProgress (0.0); - }*/ - - volatile double progress = 0.0; - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - const short int height=src->height, width=src->width; - const short int hfh=(height+1)/2, hfw=(width+1)/2; - //const short int hfh=(height+1), hfw=(width+1); - - const int blkrad=1; - float noisevar_L = SQR(dnparams.luma * TS * 40.0f); - float noisevar_ab = SQR(dnparams.chroma * TS * 200.0f); - - //int dxr=Roffset&1, dyr=(Roffset&2)/2, dxb=(1-dxr), dyb=(1-dyr); - //int rdx, rdy, bdx, bdy; - - // calculation for tiling - const int numblox_W = ceil(((float)(width))/(offset))+2*blkrad; - const int numblox_H = ceil(((float)(height))/(offset))+2*blkrad; - //const int nrtiles = numblox_W*numblox_H; - // end of tiling calc - - //const float eps = 1.0f; - - - - array2D tilemask_in(TS,TS); - array2D tilemask_out(TS,TS); - - array2D totwt(width,height,ARRAY2D_CLEAR_DATA); - - for (int i=0; iTS/2 ? i-TS : i)); - float vmask = (i1<8 ? SQR(sin(M_PI*i1/16.0f)) : 1.0f); - for (int j=0; jTS/2 ? j-TS : j)); - //tilemask_in[i][j] = vmask * (j1<4 ? SQR(sin(M_PI*(float)j1/8.0f)) : 1.0f); - tilemask_in[i][j] = (exp(-(SQR(i-TS/2-0.5f)+SQR(j-TS/2-0.5f))/(2*SQR(TS/4.0f))) * \ - vmask * (j1<8 ? SQR(sin(M_PI*(float)j1/16.0f)) : 1.0f)); - - tilemask_out[i][j] = (SQR(MAX(0.0f, sin(M_PI*(float)(i-4.0f)/(TS-9) ))) * \ - SQR(MAX(0.0f, sin(M_PI*(float)(j-4.0f)/(TS-9) )))); - //tilemask_out[i][j] = exp(-(SQR(i-TS/2-0.5f)+SQR(j-TS/2-0.5f))/(SQR(TS/4.0f))); - + for (int i=0; iH; i++) { + float X,Y,Z; + for (int j=0; jW; j++) { + //input normalized to (0,1) + //Y = igamcurveL[ src->L[i][j] ]; + Y = src->L[i][j]; + X = (5.0f*(src->a[i][j])) + Y; + Z = Y - (5.0f*(src->b[i][j])); + + X = X<32768.0f ? igamcurve[X] : (gamma((float)X/32768.0f, igam, igamthresh, igamslope, 1.0, 0.0) * 65535.0f); + Y = Y<32768.0f ? igamcurve[Y] : (gamma((float)Y/32768.0f, igam, igamthresh, igamslope, 1.0, 0.0) * 65535.0f); + Z = Z<32768.0f ? igamcurve[Z] : (gamma((float)Z/32768.0f, igam, igamthresh, igamslope, 1.0, 0.0) * 65535.0f); + + //Y = 65535.0f*(0.05+0.1*((float)rand()/(float)RAND_MAX));//test with random data + + dst->r[i][j] = prophoto_xyz[0][0]*X + prophoto_xyz[0][1]*Y + prophoto_xyz[0][2]*Z; + dst->g[i][j] = prophoto_xyz[1][0]*X + prophoto_xyz[1][1]*Y + prophoto_xyz[1][2]*Z; + dst->b[i][j] = prophoto_xyz[2][0]*X + prophoto_xyz[2][1]*Y + prophoto_xyz[2][2]*Z; + + } } } //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - LabImage * labin = new LabImage(width,height); - LabImage * labblur = new LabImage(width,height); - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - // transform RGB input to ersatz Lab - - RGB_InputTransf(src, labin, labblur, dnparams, defringe); - - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - // initialize FT and output data structures - LabImage * labdn = new LabImage(width,height); - float ** Lblox = new float *[8] ; - float ** RLblox = new float *[8] ; - float ** BLblox = new float *[8] ; - - fftwf_complex ** fLblox = new fftwf_complex *[8] ; //for FT - //float ** fLblox = new float *[8] ; //for DCT - - fftwf_complex ** fRLblox = new fftwf_complex *[8] ; - fftwf_complex ** fBLblox = new fftwf_complex *[8] ; - - for( int i = 0 ; i < 8 ; i++ ) { - Lblox[i] = (float *) fftwf_malloc (numblox_W*TS*TS * sizeof (float)); - //RLblox[i] = (float *) fftwf_malloc (numblox_W*TS*TS * sizeof (float)); - //BLblox[i] = (float *) fftwf_malloc (numblox_W*TS*TS * sizeof (float)); + void ImProcFunctions::RGB_denoise(Imagefloat * src, Imagefloat * dst, /*int Roffset,*/ const procparams::DirPyrDenoiseParams & dnparams, const procparams::DefringeParams & defringe) + { - fLblox[i] = (fftwf_complex *) fftwf_malloc (numblox_W*TS*fTS * sizeof (fftwf_complex)); //for FT - //fLblox[i] = (float *) fftwf_malloc (numblox_W*TS*TS * sizeof (float)); //for DCT + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + /*if (plistener) { + plistener->setProgressStr ("Block FT Luma Denoise..."); + plistener->setProgress (0.0); + }*/ + + volatile double progress = 0.0; + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + const short int height=src->height, width=src->width; + const short int hfh=(height+1)/2, hfw=(width+1)/2; + //const short int hfh=(height+1), hfw=(width+1); + + const int blkrad=1; + float noisevar_L = SQR(dnparams.luma * TS * 10.0f); + float noisevar_ab = SQR(dnparams.chroma * TS * 100.0f); + + //int dxr=Roffset&1, dyr=(Roffset&2)/2, dxb=(1-dxr), dyb=(1-dyr); + //int rdx, rdy, bdx, bdy; + + // calculation for tiling + const int numblox_W = ceil(((float)(width))/(offset))+2*blkrad; + const int numblox_H = ceil(((float)(height))/(offset))+2*blkrad; + //const int nrtiles = numblox_W*numblox_H; + // end of tiling calc + + //const float eps = 1.0f; + + + + array2D tilemask_in(TS,TS); + array2D tilemask_out(TS,TS); + + array2D totwt(width,height,ARRAY2D_CLEAR_DATA); + + for (int i=0; iTS/2 ? i-TS : i)); + float vmask = (i1<8 ? SQR(sin(M_PI*i1/16.0f)) : 1.0f); + float vmask2 = (i1<32 ? SQR(sin(M_PI*i1/64.0f)) : 1.0f); + for (int j=0; jTS/2 ? j-TS : j)); + //tilemask_in[i][j] = vmask * (j1<4 ? SQR(sin(M_PI*(float)j1/8.0f)) : 1.0f); + tilemask_in[i][j] = (exp(-(SQR(i-TS/2-0.5f)+SQR(j-TS/2-0.5f))/(2*SQR(TS/4.0f))) * \ + vmask * (j1<8 ? SQR(sin(M_PI*(float)j1/16.0f)) : 1.0f)); + + //tilemask_out[i][j] = (SQR(MAX(0.0f, sin(M_PI*(float)(i-8.0f)/(TS-17) ))) * \ + SQR(MAX(0.0f, sin(M_PI*(float)(j-8.0f)/(TS-17) )))); + tilemask_out[i][j] = vmask2 * (j1<32 ? SQR(sin(M_PI*(float)j1/64.0f)) : 1.0f); + + //tilemask_out[i][j] = exp(-(SQR(i-TS/2-0.5f)+SQR(j-TS/2-0.5f))/(SQR(TS/4.0f))); + + } + } + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + LabImage * labin = new LabImage(width,height); + LabImage * labblur = new LabImage(width,height); + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + // transform RGB input to ersatz Lab + + RGB_InputTransf(src, labin, dnparams, defringe); + + wavelet_decomposition Ldecomp(labin->data, labin->W, labin->H, 5 );//last arg is num levels + WaveletDenoise(Ldecomp, SQR((float)dnparams.Lamt/25.0f)); + Ldecomp.reconstruct(labblur->data); + + + //impulse_nr (dst, 50.0f/20.0f); + //PF_correct_RT(dst, dst, defringe.radius, defringe.threshold); + + + int datalen = labin->W*labin->H; + + for (int i=datalen; i<3*datalen; i++) { + labblur->data[i]=labin->data[i]; + } + /* + wavelet_decomposition adecomp(labin->data+datalen, labin->W, labin->H, 5 );//last arg is num levels + WaveletDenoise(adecomp, SQR((float)dnparams.chroma/5.0f)); + adecomp.reconstruct(labblur->data+datalen); + + wavelet_decomposition bdecomp(labin->data+2*datalen, labin->W, labin->H, 5 );//last arg is num levels + WaveletDenoise(bdecomp, SQR((float)dnparams.chroma/5.0f)); + bdecomp.reconstruct(labblur->data+2*datalen); + */ + + //dirpyr_ab(labin, labblur, dnparams);//use dirpyr here if using it to blur ab channels only + //dirpyrLab_denoise(labin, labblur, dnparams);//use dirpyr here if using it to blur ab channels only - //fRLblox[i] = (fftwf_complex *) fftwf_malloc (numblox_W*TS*fTS * sizeof (fftwf_complex)); - //fBLblox[i] = (fftwf_complex *) fftwf_malloc (numblox_W*TS*fTS * sizeof (fftwf_complex)); - } - - //make a plan for FFTW - fftwf_plan plan_forward_blox, plan_backward_blox; - - int nfwd[2]={TS,TS}; - - //for FT: - plan_forward_blox = fftwf_plan_many_dft_r2c(2, nfwd, numblox_W, Lblox[0], NULL, 1, TS*TS, \ - fLblox[0], NULL, 1, TS*fTS, FFTW_ESTIMATE ); - plan_backward_blox = fftwf_plan_many_dft_c2r(2, nfwd, numblox_W, fLblox[0], NULL, 1, TS*fTS, \ - Lblox[0], NULL, 1, TS*TS, FFTW_ESTIMATE ); - - //for DCT: - /*const fftw_r2r_kind fwdkind[2] = {FFTW_REDFT10, FFTW_REDFT10}; - const fftw_r2r_kind bwdkind[2] = {FFTW_REDFT01, FFTW_REDFT01}; - - plan_forward_blox = fftwf_plan_many_r2r(2, nfwd, numblox_W, Lblox[0], NULL, 1, TS*TS, \ - fLblox[0], NULL, 1, TS*TS, fwdkind, FFTW_ESTIMATE ); - plan_backward_blox = fftwf_plan_many_r2r(2, nfwd, numblox_W, fLblox[0], NULL, 1, TS*TS, \ - Lblox[0], NULL, 1, TS*TS, bwdkind, FFTW_ESTIMATE );*/ - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + // initialize FT and output data structures + LabImage * labdn = new LabImage(width,height); + for (int i=0; i<3*labdn->W*labdn->H; i++) { + labdn->data[i] = 0.0f; + } + float ** Lblox = new float *[8] ; + float ** RLblox = new float *[8] ; + float ** BLblox = new float *[8] ; + + fftwf_complex ** fLblox = new fftwf_complex *[8] ; //for FT + //float ** fLblox = new float *[8] ; //for DCT + + fftwf_complex ** fRLblox = new fftwf_complex *[8] ; + fftwf_complex ** fBLblox = new fftwf_complex *[8] ; + + for( int i = 0 ; i < 8 ; i++ ) { + Lblox[i] = (float *) fftwf_malloc (numblox_W*TS*TS * sizeof (float)); + //RLblox[i] = (float *) fftwf_malloc (numblox_W*TS*TS * sizeof (float)); + //BLblox[i] = (float *) fftwf_malloc (numblox_W*TS*TS * sizeof (float)); + + fLblox[i] = (fftwf_complex *) fftwf_malloc (numblox_W*TS*fTS * sizeof (fftwf_complex)); //for FT + //fLblox[i] = (float *) fftwf_malloc (numblox_W*TS*TS * sizeof (float)); //for DCT + + //fRLblox[i] = (fftwf_complex *) fftwf_malloc (numblox_W*TS*fTS * sizeof (fftwf_complex)); + //fBLblox[i] = (fftwf_complex *) fftwf_malloc (numblox_W*TS*fTS * sizeof (fftwf_complex)); + } + + //make a plan for FFTW + fftwf_plan plan_forward_blox, plan_backward_blox; + + int nfwd[2]={TS,TS}; + + //for FT: + plan_forward_blox = fftwf_plan_many_dft_r2c(2, nfwd, numblox_W, Lblox[0], NULL, 1, TS*TS, \ + fLblox[0], NULL, 1, TS*fTS, FFTW_ESTIMATE ); + plan_backward_blox = fftwf_plan_many_dft_c2r(2, nfwd, numblox_W, fLblox[0], NULL, 1, TS*fTS, \ + Lblox[0], NULL, 1, TS*TS, FFTW_ESTIMATE ); + + //for DCT: + //const fftw_r2r_kind fwdkind[2] = {FFTW_REDFT10, FFTW_REDFT10}; + //const fftw_r2r_kind bwdkind[2] = {FFTW_REDFT01, FFTW_REDFT01}; + + //plan_forward_blox = fftwf_plan_many_r2r(2, nfwd, numblox_W, Lblox[0], NULL, 1, TS*TS, \ + fLblox[0], NULL, 1, TS*TS, fwdkind, FFTW_ESTIMATE ); + //plan_backward_blox = fftwf_plan_many_r2r(2, nfwd, numblox_W, fLblox[0], NULL, 1, TS*TS, \ + Lblox[0], NULL, 1, TS*TS, bwdkind, FFTW_ESTIMATE ); + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% // Main algorithm: Tile loop #pragma omp parallel for schedule(dynamic) //int vblock=0, hblock=0; @@ -395,14 +356,14 @@ void ImProcFunctions::RGB_denoise(Imagefloat * src, Imagefloat * dst, /*int Roff int indx = (hblk)*TS;//index of block in malloc // load Lab high pass data - + for (int i=imin; iL[top+i][left+j]-labblur->L[top+i][left+j]);// luma data //RLblox[vblkmod][(indx + i)*TS+j] = tilemask_in[i][j]*(labin->a[top+i][left+j]-labblur->a[top+i][left+j]);// high pass chroma data //BLblox[vblkmod][(indx + i)*TS+j] = tilemask_in[i][j]*(labin->b[top+i][left+j]-labblur->b[top+i][left+j]);// high pass chroma data - + totwt[top+i][left+j] += tilemask_in[i][j]*tilemask_out[i][j]; } @@ -459,7 +420,7 @@ void ImProcFunctions::RGB_denoise(Imagefloat * src, Imagefloat * dst, /*int Roff } } } - + if (jmaxL[i][j]-labblur->L[i][j]; - float hpdn = labdn->L[i][j]/totwt[i][j];//note that labdn initially stores the denoised hipass data - float hgratio = 0;//MIN(1.0f,(abs(hpdn)+eps)/(abs(hporig)+eps)); - - labdn->L[i][j] = labblur->L[i][j] + (hgratio*hporig+(1-hgratio)*hpdn); - //labdn->L[i][j] = -(hporig)+0.25; - - hporig = labin->a[i][j]-labblur->a[i][j]; - hpdn = labdn->a[i][j]/totwt[i][j]; - labdn->a[i][j] = labblur->a[i][j] ;//+ (hgratio*hporig+(1-hgratio)*hpdn); - //labdn->a[i][j] = (hporig); - - - hporig = labin->b[i][j]-labblur->b[i][j]; - hpdn = labdn->b[i][j]/totwt[i][j]; - labdn->b[i][j] = labblur->b[i][j] ;//+ (hgratio*hporig+(1-hgratio)*hpdn); - //labdn->b[i][j] = (hporig); - - - //labdn->a[i][j] = labin->a[i][j]; - //labdn->b[i][j] = labin->b[i][j]; + // clean up + //#pragma omp single nowait + fftwf_destroy_plan( plan_forward_blox ); + //#pragma omp single nowait + fftwf_destroy_plan( plan_backward_blox ); + + for( int i = 0 ; i < 8 ; i++ ) { + fftwf_free ( Lblox[i]); + //fftwf_free (RLblox[i]); + //fftwf_free (BLblox[i]); + fftwf_free ( fLblox[i]); + //fftwf_free (fRLblox[i]); + //fftwf_free (fBLblox[i]); } - } - - //dirpyr_ab(labdn, labdn, dnparams);//use dirpyr here if using it to blur ab channels only - - dirpyrdenoise(labdn);//denoise ab channels using ImProcFns denoise (stripped to ab channels only) - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - // transform denoised "Lab" to output RGB - - RGB_OutputTransf(labdn, dst, dnparams); - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - delete labin; - delete labdn; - delete labblur; - -}//end of main RB_denoise - -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - -void ImProcFunctions::RGBtile_denoise (fftwf_complex ** fLblox, fftwf_complex ** fRLblox, fftwf_complex ** fBLblox, \ - int vblproc, int hblproc, int blkrad, int numblox_H, int numblox_W, \ - float noisevar_L, float noisevar_ab ) //for FT -//void ImProcFunctions::RGBtile_denoise (float ** fLblox, fftwf_complex ** fRLblox, fftwf_complex ** fBLblox, \ - int vblproc, int hblproc, int blkrad, int numblox_H, int numblox_W, \ - float noisevar_L, float noisevar_ab ) //for DCT -{ - int vblprocmod=vblproc%8; - - float eps = 0.01f/(TS*TS); //tolerance - - float RLblockvar=eps, BLblockvar=eps, Lblockvar=eps; - - for (int i=TS/4; i<3*TS/4; i++) - for (int j=TS/4; jL[i][j]/totwt[i][j];//note that labdn initially stores the denoised hipass data + + labdn->L[i][j] = labblur->L[i][j] + hpdn; + //labdn->L[i][j] = -(hporig)+0.25; + + //hpdn = labdn->a[i][j]/totwt[i][j]; + //labdn->a[i][j] = labblur->a[i][j] ;//+ hpdn; + //labdn->a[i][j] = (hporig); + + + //hpdn = labdn->b[i][j]/totwt[i][j]; + //labdn->b[i][j] = labblur->b[i][j] ;//+ hpdn; + //labdn->b[i][j] = (hporig); + + + labdn->a[i][j] = labblur->a[i][j]; + labdn->b[i][j] = labblur->b[i][j]; } - float nbrvar = (nbrsqave/9.0f)-SQR(nbrave/9.0f); - - float Lwsq = eps+SQR( Lcoeffre)+SQR( Lcoeffim); //for FT - //float Lwsq = eps+SQR( Lcoeff); //for DCT - - //float RLwsq = eps+SQR(RLcoeffre)+SQR(RLcoeffim); - //float BLwsq = eps+SQR(BLcoeffre)+SQR(BLcoeffim); - - //wsqave += Lwsq; - float Lfactor = (4*Lblockvar)/(eps+(Lwsq+nbrvar)+2*Lblockvar); - //float Lfactor = expf(-Lwsq/(9*Lblockvar)); - //float Lfactor = 1;//(2* Lblockvar)/(eps+ Lwsq+ Lblockvar); - //float RLfactor = 1;//(2*RLblockvar)/(eps+RLwsq+RLblockvar); - //float BLfactor = 1;//(2*BLblockvar)/(eps+BLwsq+BLblockvar); - Lwsq = MAX(0.0f, Lwsq-0.25*Lblockvar); - float Lshrinkfactor = Lwsq/(Lwsq + noisevar_L * Lfactor); - //float RLshrinkfactor = RLwsq/(RLwsq+noisevar_ab*RLfactor); - //float BLshrinkfactor = BLwsq/(BLwsq+noisevar_ab*BLfactor); - - //float shrinkfactor = (wsqW*labdn->H; + //cplx_wavelet_decomposition adecomp(labdn->data+numpix, labdn->W, labdn->H, 5 );//last arg is num levels + //WaveletDenoise(adecomp, SQR((float)dnparams.chroma*100.0f)); + //adecomp.reconstruct(labdn->data+numpix); + //cplx_wavelet_decomposition bdecomp(labdn->data+2*numpix, labdn->W, labdn->H, 5 );//last arg is num levels + //WaveletDenoise(bdecomp, SQR((float)dnparams.chroma*100.0f)); + //bdecomp.reconstruct(labdn->data+2*numpix); + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + // transform denoised "Lab" to output RGB + + RGB_OutputTransf(labdn, dst, dnparams); + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + delete labin; + delete labdn; + delete labblur; + + }//end of main RB_denoise - //add row of tiles to output image - for (int hblk=0; hblk < numblox_W; hblk++) { - int left = (hblk-blkrad)*offset; - int bottom = MIN( top+TS,height); - int right = MIN(left+TS, width); - int imin = MAX(0,-top); - int jmin = MAX(0,-left); - int imax = bottom - top; - int jmax = right - left; + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + void ImProcFunctions::RGBtile_denoise (fftwf_complex ** fLblox, fftwf_complex ** fRLblox, fftwf_complex ** fBLblox, \ + int vblproc, int hblproc, int blkrad, int numblox_H, int numblox_W, \ + float noisevar_L, float noisevar_ab ) //for FT + //void ImProcFunctions::RGBtile_denoise (float ** fLblox, fftwf_complex ** fRLblox, fftwf_complex ** fBLblox, \ + int vblproc, int hblproc, int blkrad, int numblox_H, int numblox_W, \ + float noisevar_L, float noisevar_ab ) //for DCT + { + int vblprocmod=vblproc%8; - int indx = hblk*TS; + const float eps = 0.01f/(TS*TS); //tolerance + const float cutoffsq = 8.0f;//frequency cutoff - for (int i=imin; iL[top+i][left+j] += tilemask_out[i][j]*bloxrow_L[(indx + i)*TS+j]/(TS*TS); - //labdn->a[top+i][left+j] += tilemask_out[i][j]*bloxrow_a[(indx + i)*TS+j]/(TS*TS); - //labdn->b[top+i][left+j] += tilemask_out[i][j]*bloxrow_b[(indx + i)*TS+j]/(TS*TS); + + float RLblockvar=eps, BLblockvar=eps, Lblockvar=eps; + + for (int i=TS/4; i<3*TS/4; i++) + for (int j=TS/4; jL[top+i][left+j] += tilemask_out[i][j]*bloxrow_L[(indx + i)*TS+j]/(TS*TS); + //labdn->a[top+i][left+j] += tilemask_out[i][j]*bloxrow_a[(indx + i)*TS+j]/(TS*TS); + //labdn->b[top+i][left+j] += tilemask_out[i][j]*bloxrow_b[(indx + i)*TS+j]/(TS*TS); + + } + } + } -} - #undef TS #undef fTS #undef offset -#undef eps - + //#undef eps + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //experimental dirpyr low-pass filter @@ -813,8 +783,10 @@ void ImProcFunctions::dirpyr_ab(LabImage * data_fine, LabImage * data_coarse, co //delete dirpyrlo[1]; } - -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + void ImProcFunctions::dirpyr_ablevel(LabImage * data_fine, LabImage * data_coarse, int width, int height, LUTf & rangefn_L, LUTf & rangefn_ab, int level, int scale) { @@ -897,6 +869,10 @@ void ImProcFunctions::dirpyr_ablevel(LabImage * data_fine, LabImage * data_coars } } + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + /* void ImProcFunctions::FixImpulse_ab(LabImage * src, LabImage * dst, double radius, int thresh) { @@ -985,6 +961,318 @@ void ImProcFunctions::FixImpulse_ab(LabImage * src, LabImage * dst, double radiu } */ + + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + void ImProcFunctions::WaveletDenoise(cplx_wavelet_decomposition &DualTreeCoeffs, float noisevar ) + { + int maxlvl = DualTreeCoeffs.maxlevel(); + int rad_stage1[8] = {3,2,1,1,1,1,1,1}; + int rad_stage2[8] = {2,1,1,1,1,1,1,1}; + + for (int lvl=0; lvl wiener1(Wlvl,Hlvl); + for (int m=0; m<2; m++) { + float ** ReCoeffs = DualTreeCoeffs.level_coeffs(lvl,0+m); + float ** ImCoeffs = DualTreeCoeffs.level_coeffs(lvl,2+m); + float ** ReParents = DualTreeCoeffs.level_coeffs(lvl+1,0+m); + float ** ImParents = DualTreeCoeffs.level_coeffs(lvl+1,2+m); + int ParentPadding = DualTreeCoeffs.level_pad(lvl+1,0+m); + for (int dir=1; dir<4; dir++) { + //FirstStageWiener (ReCoeffs[dir],ImCoeffs[dir],wiener1,Wlvl,Hlvl,rad_stage1[lvl], noisevar); + //SecondStageWiener(ReCoeffs[dir],ImCoeffs[dir],wiener1,Wlvl,Hlvl,rad_stage2[lvl], noisevar); + + BiShrink(ReCoeffs[dir], ImCoeffs[dir], ReParents[dir], ImParents[dir], Wlvl, Hlvl, lvl, ParentPadding, noisevar); + } + } + } + + } + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + void ImProcFunctions::WaveletDenoise(wavelet_decomposition &WaveletCoeffs, float noisevar ) + { + int maxlvl = WaveletCoeffs.maxlevel(); + int rad_stage1[8] = {3,2,1,1,1,1,1,1}; + int rad_stage2[8] = {2,1,1,1,1,1,1,1}; + + for (int lvl=0; lvl wiener1(Wlvl,Hlvl); + int ParentPadding; + float ** WavParents; + float ** WavCoeffs = WaveletCoeffs.level_coeffs(lvl); + if (lvl* buffer = new AlignedBuffer (W*H); + float* temp = buffer->data; + + if (rad==0) { + for (int row=0; row void copy_out(A ** a, B * b, size_t datalen) -{ +{// for complex wavelet decomposition for (size_t j=0; j (0.25*(a[0][j]+a[1][j]+a[2][j]+a[3][j])); } } + +template +void copy_out(A * a, B * b, size_t datalen) +{// for standard wavelet decomposition + for (size_t j=0; j (a[j]); + } +} // %%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -70,7 +78,7 @@ private: float *testfilt_anal; float *testfilt_synth; - cplx_wavelet_level * dual_tree_coeffs[maxlevels][4];//m_c in old code + wavelet_level * dual_tree[maxlevels][4]; public: @@ -78,11 +86,35 @@ public: cplx_wavelet_decomposition(E * src, int width, int height, int maxlvl); ~cplx_wavelet_decomposition(); + + internal_type ** level_coeffs(int level, int branch) const + { + return dual_tree[level][branch]->subbands(); + } + + int level_W(int level, int branch) const + { + return dual_tree[level][branch]->width(); + } + int level_H(int level, int branch) const + { + return dual_tree[level][branch]->height(); + } + + int level_pad(int level, int branch) const + { + return dual_tree[level][branch]->padding(); + } + + int maxlevel() const + { + return lvltot; + } + template void reconstruct(E * dst); - }; @@ -98,20 +130,34 @@ public: //initialize wavelet filters - first_lev_len = Kingsbury_len; - first_lev_offset = Kingsbury_offset; + first_lev_len = FSFarras_len; + first_lev_offset = FSFarras_offset; first_lev_anal = new float[4*first_lev_len]; first_lev_synth = new float[4*first_lev_len]; for (int n=0; n<2; n++) { for (int m=0; m<2; m++) { for (int i=0; i(src, padding, m_w, m_h, first_lev_anal+first_lev_len*2*n, \ + lvltot=0; + float padding = 0;//1<<(maxlvl-1); + dual_tree[0][2*n+m] = new wavelet_level(src, lvltot, padding, m_w, m_h, first_lev_anal+first_lev_len*2*n, \ first_lev_anal+first_lev_len*2*m, first_lev_len, first_lev_offset); - lvltot=1; while(lvltot < maxlvl) { - dual_tree_coeffs[lvltot][2*n+m] = new cplx_wavelet_level(dual_tree_coeffs[lvltot-1][2*n+m]->lopass()/*lopass*/, 0/*no padding*/, \ - dual_tree_coeffs[lvltot-1][2*n+m]->width(), \ - dual_tree_coeffs[lvltot-1][2*n+m]->height(), \ - wavfilt_anal+wavfilt_len*2*n, wavfilt_anal+wavfilt_len*2*m, wavfilt_len, wavfilt_offset); lvltot++; + dual_tree[lvltot][2*n+m] = new wavelet_level(dual_tree[lvltot-1][2*n+m]->lopass()/*lopass*/, lvltot, 0/*no padding*/, \ + dual_tree[lvltot-1][2*n+m]->width(), \ + dual_tree[lvltot-1][2*n+m]->height(), \ + wavfilt_anal+wavfilt_len*2*n, wavfilt_anal+wavfilt_len*2*m, \ + wavfilt_len, wavfilt_offset); } } } //rotate detail coefficients + float coeffave[5][4][3]; + float root2 = sqrt(2); for (int lvl=0; lvlwidth(); - int Hlvl = dual_tree_coeffs[lvl][0]->height(); - for (int i=0; iwavcoeffs[m][i] + dual_tree_coeffs[lvl][3]->wavcoeffs[m][i])/root2; - dual_tree_coeffs[lvl][3]->wavcoeffs[m][i] = (dual_tree_coeffs[lvl][0]->wavcoeffs[m][i] - dual_tree_coeffs[lvl][3]->wavcoeffs[m][i])/root2; - dual_tree_coeffs[lvl][0]->wavcoeffs[m][i] = wavtmp; + int Wlvl = dual_tree[lvl][0]->width(); + int Hlvl = dual_tree[lvl][0]->height(); + for (int n=0; n<4; n++) + for (int m=1; m<4; m++) + coeffave[lvl][n][m-1]=0; + + for (int m=1; m<4; m++) {//detail coefficients only + for (int i=0; iwavcoeffs[m][i] + dual_tree_coeffs[lvl][2]->wavcoeffs[m][i])/root2; - dual_tree_coeffs[lvl][2]->wavcoeffs[m][i] = (dual_tree_coeffs[lvl][1]->wavcoeffs[m][i] - dual_tree_coeffs[lvl][2]->wavcoeffs[m][i])/root2; - dual_tree_coeffs[lvl][1]->wavcoeffs[m][i] = wavtmp; + float wavtmp = (dual_tree[lvl][0]->wavcoeffs[m][i] + dual_tree[lvl][3]->wavcoeffs[m][i])/root2; + dual_tree[lvl][3]->wavcoeffs[m][i] = (dual_tree[lvl][0]->wavcoeffs[m][i] - dual_tree[lvl][3]->wavcoeffs[m][i])/root2; + dual_tree[lvl][0]->wavcoeffs[m][i] = wavtmp; + + wavtmp = (dual_tree[lvl][1]->wavcoeffs[m][i] + dual_tree[lvl][2]->wavcoeffs[m][i])/root2; + dual_tree[lvl][2]->wavcoeffs[m][i] = (dual_tree[lvl][1]->wavcoeffs[m][i] - dual_tree[lvl][2]->wavcoeffs[m][i])/root2; + dual_tree[lvl][1]->wavcoeffs[m][i] = wavtmp; + + for (int n=0; n<4; n++) coeffave[lvl][n][m-1] += fabs(dual_tree[lvl][n]->wavcoeffs[m][i]); } } + for (int n=0; n<4; n++) + for (int i=0; i<3; i++) + coeffave[lvl][n][i] /= Wlvl*Hlvl; } - + } /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ - - /* function y=reconstruct(w,J,Fsf,sf) */ - + template void cplx_wavelet_decomposition::reconstruct(E * dst) { @@ -183,22 +245,21 @@ public: //rotate detail coefficients float root2 = sqrt(2); for (int lvl=0; lvlwidth(); - int Hlvl = dual_tree_coeffs[lvl][0]->height(); + int Wlvl = dual_tree[lvl][0]->width(); + int Hlvl = dual_tree[lvl][0]->height(); for (int i=0; iwavcoeffs[m][i] + dual_tree_coeffs[lvl][3]->wavcoeffs[m][i])/root2; - dual_tree_coeffs[lvl][3]->wavcoeffs[m][i] = (dual_tree_coeffs[lvl][0]->wavcoeffs[m][i] - dual_tree_coeffs[lvl][3]->wavcoeffs[m][i])/root2; - dual_tree_coeffs[lvl][0]->wavcoeffs[m][i] = wavtmp; + float wavtmp = (dual_tree[lvl][0]->wavcoeffs[m][i] + dual_tree[lvl][3]->wavcoeffs[m][i])/root2; + dual_tree[lvl][3]->wavcoeffs[m][i] = (dual_tree[lvl][0]->wavcoeffs[m][i] - dual_tree[lvl][3]->wavcoeffs[m][i])/root2; + dual_tree[lvl][0]->wavcoeffs[m][i] = wavtmp; - wavtmp = (dual_tree_coeffs[lvl][1]->wavcoeffs[m][i] + dual_tree_coeffs[lvl][2]->wavcoeffs[m][i])/root2; - dual_tree_coeffs[lvl][2]->wavcoeffs[m][i] = (dual_tree_coeffs[lvl][1]->wavcoeffs[m][i] - dual_tree_coeffs[lvl][2]->wavcoeffs[m][i])/root2; - dual_tree_coeffs[lvl][1]->wavcoeffs[m][i] = wavtmp; + wavtmp = (dual_tree[lvl][1]->wavcoeffs[m][i] + dual_tree[lvl][2]->wavcoeffs[m][i])/root2; + dual_tree[lvl][2]->wavcoeffs[m][i] = (dual_tree[lvl][1]->wavcoeffs[m][i] - dual_tree[lvl][2]->wavcoeffs[m][i])/root2; + dual_tree[lvl][1]->wavcoeffs[m][i] = wavtmp; } } } - //y = ConstantArray[0, {vsizetmp, hsizetmp}]; internal_type ** tmp = new internal_type *[4]; for (int i=0; i<4; i++) { tmp[i] = new internal_type[m_w*m_h]; @@ -206,12 +267,14 @@ public: for (int n=0; n<2; n++) { for (int m=0; m<2; m++) { + int skip=1<<(lvltot-1); for (int lvl=lvltot-1; lvl>0; lvl--) { - dual_tree_coeffs[lvl][2*n+m]->reconstruct_level(dual_tree_coeffs[lvl-1][2*n+m]->wavcoeffs[0], wavfilt_synth+wavfilt_len*2*n, \ - wavfilt_synth+wavfilt_len*2*m, wavfilt_len, wavfilt_offset); + dual_tree[lvl][2*n+m]->reconstruct_level(dual_tree[lvl-1][2*n+m]->wavcoeffs[0], wavfilt_synth+wavfilt_len*2*n, \ + wavfilt_synth+wavfilt_len*2*m, wavfilt_len, wavfilt_offset, skip); + skip /=2; } - dual_tree_coeffs[0][2*n+m]->reconstruct_level(tmp[2*n+m], first_lev_synth+first_lev_len*2*n, - first_lev_synth+first_lev_len*2*m, first_lev_len, first_lev_offset); + dual_tree[0][2*n+m]->reconstruct_level(tmp[2*n+m], first_lev_synth+first_lev_len*2*n, + first_lev_synth+first_lev_len*2*m, first_lev_len, first_lev_offset, skip); } } @@ -226,9 +289,156 @@ public: } + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + class wavelet_decomposition + { + public: + + typedef float internal_type; + + private: + + static const int maxlevels = 8;//should be greater than any conceivable order of decimation + + int lvltot; + size_t m_w, m_h;//dimensions + size_t m_w1, m_h1; + + int wavfilt_len, wavfilt_offset; + float *wavfilt_anal; + float *wavfilt_synth; + + int testfilt_len, testfilt_offset; + float *testfilt_anal; + float *testfilt_synth; + + wavelet_level * wavelet_decomp[maxlevels]; + + public: + + template + wavelet_decomposition(E * src, int width, int height, int maxlvl); + + ~wavelet_decomposition(); + + internal_type ** level_coeffs(int level) const + { + return wavelet_decomp[level]->subbands(); + } + + int level_W(int level) const + { + return wavelet_decomp[level]->width(); + } + + int level_H(int level) const + { + return wavelet_decomp[level]->height(); + } + + int level_pad(int level) const + { + return wavelet_decomp[level]->padding(); + } + + int maxlevel() const + { + return lvltot; + } + + template + void reconstruct(E * dst); + + }; + + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + template + wavelet_decomposition::wavelet_decomposition(E * src, int width, int height, int maxlvl) + : lvltot(0), m_w(width), m_h(height), m_w1(0), m_h1(0) + { + m_w1 = width; + m_h1 = height; + + //initialize wavelet filters + + + wavfilt_len = Haar_len; + wavfilt_offset = Haar_offset; + wavfilt_anal = new float[2*wavfilt_len]; + wavfilt_synth = new float[2*wavfilt_len]; + + for (int n=0; n<2; n++) { + for (int i=0; i(src, lvltot/*level*/, padding/*padding*/, m_w, m_h, \ + wavfilt_anal, wavfilt_anal, wavfilt_len, wavfilt_offset); + while(lvltot < maxlvl) { + lvltot++; + wavelet_decomp[lvltot] = new wavelet_level(wavelet_decomp[lvltot-1]->lopass()/*lopass*/, lvltot/*level*/, 0/*no padding*/, \ + wavelet_decomp[lvltot-1]->width(), wavelet_decomp[lvltot-1]->height(), \ + wavfilt_anal, wavfilt_anal, wavfilt_len, wavfilt_offset); + } + + + } + + /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ + /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ + + + template + void wavelet_decomposition::reconstruct(E * dst) { + + // data structure is wavcoeffs[scale][channel={lo,hi1,hi2,hi3}][pixel_array] + + int skip=1<<(lvltot-1); + for (int lvl=lvltot-1; lvl>0; lvl--) { + wavelet_decomp[lvl]->reconstruct_level(wavelet_decomp[lvl-1]->wavcoeffs[0], wavfilt_synth, wavfilt_synth, wavfilt_len, wavfilt_offset, skip); + skip /=2; + } + + internal_type * tmp = new internal_type[m_w*m_h]; + wavelet_decomp[0]->reconstruct_level(tmp, wavfilt_synth, wavfilt_synth, wavfilt_len, wavfilt_offset, skip); + + copy_out(tmp,dst,m_w*m_h); + + delete[] tmp; + + + + } + + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + }; diff --git a/rtengine/cplx_wavelet_filter_coeffs.h b/rtengine/cplx_wavelet_filter_coeffs.h index 724d276cd..0bc7a9786 100644 --- a/rtengine/cplx_wavelet_filter_coeffs.h +++ b/rtengine/cplx_wavelet_filter_coeffs.h @@ -32,11 +32,10 @@ namespace rtengine { // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -/*CplxWavelet AntonB = { - 12,//length of filter - 6,//offset +/*const int AntonB_len = 12;//length of filter +const int AntonB_offset = 6;//offset - {//analysis filter +const float AntonB_anal[2][2][12] = {//analysis filter {{0, -0.08838834764832, 0.08838834764832, 0.69587998903400, 0.69587998903400, 0.08838834764832, -0.08838834764832, 0.01122679215254, 0.01122679215254, 0}, {0, 0, 0, 0.04563588155712, -0.02877176311425, -0.29563588155712 , @@ -44,9 +43,9 @@ namespace rtengine { {{0 , 0 , 0.02674875741081, -0.01686411844287, -0.07822326652899, 0.26686411844288, 0.60294901823636, 0.26686411844287, -0.07822326652899, -0.01686411844287, 0.02674875741081, 0}, {0 , 0 , 0, 0 , 0.04563588155712, -0.02877176311425, - -0.29563588155712 , 0.55754352622850, -0.29563588155713, -0.02877176311425, 0.04563588155712 , 0}} }, + -0.29563588155712 , 0.55754352622850, -0.29563588155713, -0.02877176311425, 0.04563588155712 , 0}} }; - {//synthesis filter +const float AntonB_synth[2][2][12] = {//synthesis filter {{0 , 0 , 0, -0.04563588155712, -0.02877176311425, 0.29563588155712, 0.55754352622850, 0.29563588155713, -0.02877176311425, -0.04563588155712, 0, 0}, {0, 0.02674875741081, 0.01686411844287, -0.07822326652899, -0.26686411844288 , 0.60294901823636, @@ -54,8 +53,8 @@ namespace rtengine { {{0 , 0, -0.04563588155712, -0.02877176311425, 0.29563588155712 , 0.55754352622850 , 0.29563588155713, -0.02877176311425, -0.04563588155712, 0, 0 , 0}, {0.02674875741081 , 0.01686411844287, -0.07822326652899, -0.26686411844288 , 0.60294901823636, -0.26686411844287, - -0.07822326652899, 0.01686411844287 , 0.02674875741081 , 0 , 0, 0}} } -};*/ + -0.07822326652899, 0.01686411844287 , 0.02674875741081 , 0 , 0, 0}} };*/ + /* for (int i=0; i<4; i++) for (int n=0; n<12; n++) { @@ -66,7 +65,7 @@ namespace rtengine { const int FSFarras_len=10;//length of filter -const int FSFarras_offset=5;//offset +const int FSFarras_offset=4;//offset const float FSFarras_anal[2][2][10] = {//analysis filter {{0, -0.08838834764832, 0.08838834764832, 0.69587998903400, 0.69587998903400, 0.08838834764832, -0.08838834764832, 0.01122679215254 , 0.01122679215254, 0}, @@ -92,7 +91,7 @@ const float FSFarras_anal[2][2][10] = {//analysis filter % Image Proc.(ICIP),2000 */ const int Kingsbury_len=10;//length of filter -const int Kingsbury_offset=5;//offset +const int Kingsbury_offset=4;//offset const float Kingsbury_anal[2][2][10] = {//analysis filter {{0.03516384000000, 0, -0.08832942000000, 0.23389032000000, 0.76027237000000, 0.58751830000000, 0, -0.11430184000000 , 0, 0}, @@ -103,6 +102,15 @@ const float Kingsbury_anal[2][2][10] = {//analysis filter //synthesis filter is the reverse (see cplx_wavelet_dec.h) /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ + + const int Haar_len=2;//length of filter + const int Haar_offset=1;//offset + + const float Haar_anal[2][2] = {{0.5,0.5}, {0.5,-0.5}};//analysis filter + + //synthesis filter is the reverse (see cplx_wavelet_dec.h) + + /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ }; diff --git a/rtengine/cplx_wavelet_level.h b/rtengine/cplx_wavelet_level.h index 5ef5a28dc..dc2aaca58 100644 --- a/rtengine/cplx_wavelet_level.h +++ b/rtengine/cplx_wavelet_level.h @@ -25,17 +25,19 @@ #include #include "array2D.h" +#include "gauss.h" namespace rtengine { #define MAX(a,b) ((a) > (b) ? (a) : (b)) #define MIN(a,b) ((a) > (b) ? (b) : (a)) - +#define SQR(x) ((x)*(x)) + ////////////////////////////////////////////////////////////////////////////// template - class cplx_wavelet_level + class wavelet_level { // full size size_t m_w, m_h; @@ -46,6 +48,12 @@ namespace rtengine { // size of padded border size_t m_pad; + // level of decomposition + int lvl; + + // spacing of filter taps + size_t skip; + // array of pointers to lines of coeffs // actually is a single contiguous data array pointed by m_coeffs[0] //T ** m_coeffs; @@ -67,32 +75,32 @@ namespace rtengine { //void dwt_2d(size_t w, size_t h); //void idwt_2d(size_t w, size_t h, int alpha); - void AnalysisFilter (T * src, T * dstLo, T * dstHi, float *filterLo, float *filterHi, + void AnalysisFilter (T * srcbuffer, T * dstLo, T * dstHi, float *filterLo, float *filterHi, int taps, int offset, int pitch, int srclen); void SynthesisFilter (T * srcLo, T * srcHi, T * dst, T *bufferLo, T *bufferHi, float *filterLo, float *filterHi, int taps, int offset, int pitch, int dstlen); + void imp_nr (T* src, int width, int height, double thresh); + + public: T ** wavcoeffs; template - cplx_wavelet_level(E * src, int padding, size_t w, size_t h, float *filterV, float *filterH, int len, int offset) - : m_w(w), m_h(h), m_w2((w+1+2*padding)/2), m_h2((h+1+2*padding)/2), m_pad(padding), wavcoeffs(NULL) + wavelet_level(E * src, int level, int padding, size_t w, size_t h, float *filterV, float *filterH, int len, int offset) + : m_w(w), m_h(h), m_w2(w), m_h2(h), m_pad(padding), wavcoeffs(NULL), lvl(level), skip(1< - void decompose_level(E *src, float *filterV, float *filterH, int len, int offset); + size_t padding() const + { + return m_pad/skip; + } template - void reconstruct_level(E *dst, float *filterV, float *filterH, int len, int offset); + void decompose_level(E *src, float *filterV, float *filterH, int len, int offset, int skip); + + template + void reconstruct_level(E *dst, float *filterV, float *filterH, int len, int offset, int skip); }; @@ -130,7 +143,7 @@ namespace rtengine { template - T ** cplx_wavelet_level::create(size_t n) + T ** wavelet_level::create(size_t n) { T * data = new T[4*n]; T ** subbands = new T*[4]; @@ -145,7 +158,7 @@ namespace rtengine { template - void cplx_wavelet_level::destroy(T ** subbands) + void wavelet_level::destroy(T ** subbands) { if(subbands) { @@ -157,22 +170,52 @@ namespace rtengine { // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% template template - void cplx_wavelet_level::loadbuffer(E * src, E * dst, int pitch, int srclen) + void wavelet_level::loadbuffer(E * src, E * dst, int pitch, int srclen) { E * tmp = dst + m_pad; - memset(dst, 0, (srclen+2*m_pad)*sizeof(E)); - - for(size_t i = 0, j = 0; i - void cplx_wavelet_level::AnalysisFilter (T * src, T * dstLo, T * dstHi, float *filterLo, float *filterHi, + void wavelet_level::AnalysisFilter (T * srcbuffer, T * dstLo, T * dstHi, float *filterLo, float *filterHi, int taps, int offset, int pitch, int srclen) { /* Basic convolution code @@ -195,26 +238,49 @@ namespace rtengine { * treatment of mirror BC's is implemented. * */ - - // calculate coefficients - for(int i = 0; i < (srclen); i+=2) { + //input data is 'skip' rows and cosetlen=srclen/skip columns (which includes padding at either and) + + /*int cosetlen = srclen/skip; + + for (size_t coset=0; cosettaps && itaps && iskip*taps && i - void cplx_wavelet_level::SynthesisFilter (T * srcLo, T * srcHi, T * dst, T *bufferLo, T *bufferHi, - float *filterLo, float *filterHi, int taps, int offset, int pitch, int dstlen) { + void wavelet_level::SynthesisFilter (T * srcLo, T * srcHi, T * dst, T *bufferLo, T *bufferHi, float *filterLo, + float *filterHi, int taps, int offset, int pitch, int dstlen) { /* Basic convolution code * Applies an FIR filter 'filter' with 'len' taps, @@ -239,34 +305,66 @@ namespace rtengine { - // calculate coefficients - - int srclen=(dstlen+1+2*m_pad)/2; - for (int i=0; itaps && i<(cosetlen-taps)) {//bulk + for (int j=0, l=-shift; jtaps && i<(srclen-taps)) {//bulk - for (int j=begin, l=0; jskip*taps && i<(srclen-skip*taps)) {//bulk + for (int j=0, l=-skip*shift; j65535.0f) { + float xxx=tot; + float yyy=1.0f; + } } } @@ -274,26 +372,28 @@ namespace rtengine { // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% template template - void cplx_wavelet_level::decompose_level(E *src, float *filterV, float *filterH, int taps, int offset) { + void wavelet_level::decompose_level(E *src, float *filterV, float *filterH, int taps, int offset, int skip) { T *tmpLo = new T[m_w*m_h2]; T *tmpHi = new T[m_w*m_h2]; - T *buffer = new T[MAX(m_w,m_h)+2*m_pad]; + T *buffer = new T[MAX(m_w2,m_h2)]; /* filter along columns */ for (int j=0; j template - void cplx_wavelet_level::reconstruct_level(E *dst, float *filterV, float *filterH, int taps, int offset) { + void wavelet_level::reconstruct_level(E *dst, float *filterV, float *filterH, int taps, int offset, int skip) { - - //int hfw = (W+1)/2; - //int hfh = (H+1)/2; T *tmpLo = new T[m_w*m_h2]; T *tmpHi = new T[m_w*m_h2]; - int buflen = MAX(m_w,m_h); + int buflen = MAX(m_w2,m_h2); float *bufferLo = new float[buflen]; float *bufferHi = new float[buflen]; /* filter along rows */ for (int i=0; i + void wavelet_level::imp_nr (T* src, int width, int height, double thresh) { + + + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + // impulse noise removal + // local variables + + float hpfabs, hfnbrave; + const float eps = 0.01; + + // buffer for the lowpass image + float * lpf = new float[width*height]; + // buffer for the highpass image + float * impish = new float[width*height]; + + //The cleaning algorithm starts here + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + // modified bilateral filter for lowpass image, omitting input pixel; or Gaussian blur + /* + static float eps = 1.0; + float wtdsum[3], dirwt, norm; + int i1, j1; + + AlignedBuffer* buffer = new AlignedBuffer (MAX(width,height)); + + gaussHorizontal (src, lpf, buffer, width, height, MAX(2.0,thresh-1.0), false); + gaussVertical (lpf, lpf, buffer, width, height, MAX(2.0,thresh-1.0), false); + + delete buffer; + */ + + boxblur(src, lpf, 2, 2, width, height); + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + float impthr = MAX(1.0,5.5-thresh); + + for (int i=0; i < height; i++) + for (int j=0; j < width; j++) { + hpfabs = fabs(src[i*width+j]-lpf[i*width+j]); + //block average of high pass data + for (int i1=MAX(0,i-2), hfnbrave=0; i1<=MIN(i+2,height-1); i1++ ) + for (int j1=MAX(0,j-2); j1<=MIN(j+2,width-1); j1++ ) { + hfnbrave += fabs(src[i1*width+j1]-lpf[i1*width+j1]); + } + hfnbrave = (hfnbrave-hpfabs)/24; + hpfabs>(hfnbrave*impthr) ? impish[i*width+j]=1 : impish[i*width+j]=0; + + }//now impulsive values have been identified + + for (int i=0; i < height; i++) + for (int j=0; j < width; j++) { + if (!impish[i*width+j]) continue; + float norm=0.0; + float wtdsum=0.0; + for (int i1=MAX(0,i-2), hfnbrave=0; i1<=MIN(i+2,height-1); i1++ ) + for (int j1=MAX(0,j-2); j1<=MIN(j+2,width-1); j1++ ) { + if (i1==i && j1==j) continue; + if (impish[i1*width+j1]) continue; + float dirwt = 1/(SQR(src[i1*width+j1]-src[i*width+j])+eps);//use more sophisticated rangefn??? + wtdsum += dirwt*src[i1*width+j1]; + norm += dirwt; + } + //wtdsum /= norm; + if (norm) { + src[i*width+j]=wtdsum/norm;//low pass filter + } + + }//now impulsive values have been corrected + + delete [] lpf; + delete [] impish; + + } + }; diff --git a/rtengine/gauss.h b/rtengine/gauss.h index 19ff31179..8f78e9cee 100644 --- a/rtengine/gauss.h +++ b/rtengine/gauss.h @@ -233,9 +233,11 @@ template void gaussDerivH (T** src, T** dst, AlignedBuffer* buf double* temp = buffer->data; for (int j=1; j void gaussDerivH (T** src, T** dst, AlignedBuffer* buf for (int i=0; idata; - double src0 = 0.5*(src[i][1]-src[i][0]); + double src0 = (src[i][1]-src[i][0]); temp2[0] = B * src0 + b1*src0 + b2*src0 + b3*src0; temp2[1] = B * 0.5*(src[i][2]-src[i][0]) + b1*temp2[0] + b2*src0 + b3*src0; @@ -282,7 +284,7 @@ template void gaussDerivH (T** src, T** dst, AlignedBuffer* buf for (int j=3; j void gaussDerivV (T** src, T** dst, AlignedBuffer* buf #ifdef _OPENMP #pragma omp for #endif - for (int i=0; idata; - for (int j = 1; j @@ -558,10 +559,68 @@ void ImProcFunctions::colorCurve (LabImage* lold, LabImage* lnew) { void ImProcFunctions::impulsedenoise (LabImage* lab) { if (params->impulseDenoise.enabled && lab->W>=8 && lab->H>=8) + { + //impulse_nr (lab, (float)params->impulseDenoise.thresh/10.0 );//20 is normal - //impulse_nr (lab, (float)params->impulseDenoise.thresh/20.0 ); - { cplx_wavelet_decomposition Ldecomp(lab->data, lab->W, lab->H, 1 /*maxlvl*/); - Ldecomp.reconstruct(lab->data);} + for (int i=0; iW*lab->H; i++) { + lab->data[i] *= lab->data[i]/32768.0f; + } + wavelet_decomposition Ldecomp(lab->data, lab->W, lab->H, 5 );//last arg is num levels + //WaveletDenoise(Ldecomp, SQR((float)params->impulseDenoise.thresh*25.0f)); + WaveletDenoise(Ldecomp, SQR((float)params->impulseDenoise.thresh/25.0f)); + + LabImage* labtmp = new LabImage (lab->W,lab->H); + + int lvl = (params->impulseDenoise.thresh>>4)&7; + int branch = (params->impulseDenoise.thresh>>2)&1;//2*re_im + dir + int subband = params->impulseDenoise.thresh&3;//orientation for detail subbands + float noisevar = SQR((float)params->defringe.threshold * 10.0f); + /*for (int i=0; iW*lab->H; i++) { + //float recoeff = Ldecomp.level_coeffs(lvl,branch)[subband][i]/(2<data[i] = sqrt(SQR(recoeff)+SQR(imcoeff)) * (subband != 0 ? 2*shrink : 0.707); + lab->data[i] = Ldecomp.level_coeffs(lvl,branch)[subband][i] + (subband != 0 ? 10000 : 0); + }*/ + //for (int i=0; iW*lab->H; i++) { + // Ldecomp.level_coeffs(4)[0][i] = 0; + //} + Ldecomp.reconstruct(labtmp->data); + + //double radius = (int)(params->impulseDenoise.thresh/10) ; + //boxvar(lab->data, lab->data, radius, radius, lab->W, lab->H); + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + /*int w=lab->W; + int h=lab->H; + for(int y = 0; y != h-1; y++){ + float *rg = &lab->data[w*y]; + for(int x = 0; x != w-1; x++){ + float gx = (fabs((rg[x + 1] - rg[x]) + (rg[x + w + 1] - rg[x + w]))); + float gy = (fabs((rg[x + w] - rg[x]) + (rg[x + w + 1] - rg[x + 1]))); + lab->data[w*y+x] = gx+gy;//sqrt(lab->data[i]/32768.0f)*32768.0f; + } + }*/ + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + //double radius = (double)params->impulseDenoise.thresh/40.0 ; + //AlignedBuffer* buffer = new AlignedBuffer (MAX(lab->W,lab->H)); + //gaussDerivH (lab->L, lab->L, buffer, lab->W, lab->H, radius, false/*multiThread*/); + //gaussVertical(lab->L, lab->L, buffer, lab->W, lab->H, radius, false/*multiThread*/); + //delete buffer; + + //impulse_nr (labtmp, 50.0f/20.0f); + + for (int i=0; iW*lab->H; i++) { + //lab->data[i] = 4*(labtmp->data[i]-lab->data[i])+10000; + lab->data[i] = sqrt(MAX(0,labtmp->data[i]/32768.0f))*32768.0f; + } + delete labtmp; + + impulse_nr (lab, 50.0f/20.0f); + + } } void ImProcFunctions::defringe (LabImage* lab) { @@ -646,7 +705,6 @@ fclose(f);*/ void ImProcFunctions::getAutoExp (LUTu & histogram, int histcompr, double defgain, double clip, \ double& expcomp, int& bright, int& contr, int& black, int& hlcompr, int& hlcomprthresh) { - double corr = 1;//pow(2.0, defgain);//defgain may be redundant legacy of superceded code??? float scale = 65536.0; float midgray=0.15;//0.18445f;//middle gray in linear gamma = 0.18445*65535 @@ -677,15 +735,13 @@ fclose(f);*/ octile[count] += histogram[i]; if (octile[count]>sum/8 || (count==7 && octile[count]>sum/16)) { octile[count]=log(1+i)/log(2); - count++;// = MIN(count+1,7); + count++; } } if (i2 ? (octile[i+1]-octile[3]) : (octile[3]-octile[i]))); } - ospread /= 5; + ospread /= 5;//average width of octiles // compute clipping points based on the original histograms (linear, without exp comp.) int clipped = 0; @@ -743,13 +799,9 @@ fclose(f);*/ //sets the mean or median at middle gray, and the amount that sets the estimated top //of the histogram at or near clipping. - float expcomp1 = log(/*(median/ave)*//*(hidev/lodev)*/midgray*scale/(ave-shc+midgray*shc))/log(2); + float expcomp1 = log(midgray*scale/(ave-shc+midgray*shc))/log(2); float expcomp2 = 0.5*( (15.5f-histcompr-(2*octile[7]-octile[6])) + log(scale/rawmax)/log(2) ); - /*expcomp = (expcomp1*fabs(expcomp2)+expcomp2*fabs(expcomp1))/(fabs(expcomp1)+fabs(expcomp2)); - if (expcomp<0) { - MIN(0.0f,MAX(expcomp1,expcomp2)); - }*/ expcomp = 0.5 * (expcomp1 + expcomp2); float gain = exp(expcomp*log(2)); @@ -765,7 +817,7 @@ fclose(f);*/ float shoulder = ((scale/MAX(1,gain))*(hlcomprthresh/200.0))+0.1; //this is a series approximation of the actual formula for comp, //which is a transcendental equation - float comp = (gain*((float)whiteclip)/scale - 1)*2;//*(1-shoulder/scale); + float comp = (gain*((float)whiteclip)/scale - 1)*2; hlcompr=(int)(100*comp/(MAX(0,expcomp) + 1.0)); hlcompr = MAX(0,MIN(100,hlcompr)); @@ -778,52 +830,12 @@ fclose(f);*/ bright = (midgray-midtmp)*15.0/(0.10833-0.0833*midtmp); } - bright = 0.25*/*(median/ave)*(hidev/lodev)*/MAX(0,bright); + bright = 0.25*MAX(0,bright); //compute contrast that spreads the average spacing of octiles contr = 50.0*(1.1-ospread); contr = MAX(0,MIN(100,contr)); - //diagnostics - //printf ("**************** AUTO LEVELS ****************\n"); - //printf ("gain1= %f gain2= %f gain= %f\n",expcomp1,expcomp2,gain); - //printf ("median: %i average: %f median/average: %f\n",median,ave, median/ave); - //printf ("average: %f\n",ave); - //printf ("median/average: %f\n",median/ave); - //printf ("lodev: %f hidev: %f hidev/lodev: %f\n",lodev,hidev,hidev/lodev); - //printf ("lodev: %f\n",lodev); - //printf ("hidev: %f\n",hidev); - //printf ("rawmax= %d whiteclip= %d gain= %f\n",rawmax,whiteclip,gain); - - //printf ("octiles: %f %f %f %f %f %f %f %f\n",octile[0],octile[1],octile[2],octile[3],octile[4],octile[5],octile[6],octile[7]); - //printf ("ospread= %f\n",ospread); - - - /* - // %%%%%%%%%% LEGACY AUTOEXPOSURE CODE %%%%%%%%%%%%% - // black point selection is based on the linear result (yielding better visual results) - black = (int)(shc * corr); - // compute the white point of the exp. compensated gamma corrected image - double whiteclipg = (int)(CurveFactory::gamma2 (whiteclip * corr / 65536.0) * 65536.0); - - // compute average intensity of the exp compensated, gamma corrected image - double gavg = 0; - for (int i=0; i<65536>>histcompr; i++) - gavg += histogram[i] * CurveFactory::gamma2((int)(corr*(i<10.0) expcomp = 10.0; diff --git a/rtengine/improcfun.h b/rtengine/improcfun.h index 859662ba5..4ff124add 100644 --- a/rtengine/improcfun.h +++ b/rtengine/improcfun.h @@ -29,6 +29,7 @@ #include "labimage.h" #include "LUT.h" #include +#include "cplx_wavelet_dec.h" namespace rtengine { @@ -144,7 +145,7 @@ namespace rtengine { //void L_denoise(Imagefloat * src, LabImage * dst, const procparams::DirPyrDenoiseParams & dnparams);//Emil's FT denoise //void tile_denoise (fftwf_complex ** fLblox, int vblproc, int hblproc, \ int blkrad, int numblox_H, int numblox_W, float noisevar ); - void RGB_InputTransf(Imagefloat * src, LabImage * dst, LabImage * blur, const procparams::DirPyrDenoiseParams & dnparams, const procparams::DefringeParams & defringe); + void RGB_InputTransf(Imagefloat * src, LabImage * dst, const procparams::DirPyrDenoiseParams & dnparams, const procparams::DefringeParams & defringe); void RGB_OutputTransf(LabImage * src, Imagefloat * dst, const procparams::DirPyrDenoiseParams & dnparams); void output_tile_row (float *Lbloxrow, float ** Lhipassdn, float ** tilemask, int height, int width, int top, int blkrad ); void RGB_denoise(Imagefloat * src, Imagefloat * dst, const procparams::DirPyrDenoiseParams & dnparams, const procparams::DefringeParams & defringe); @@ -157,6 +158,17 @@ namespace rtengine { void dirpyr_ab(LabImage * data_fine, LabImage * data_coarse, const procparams::DirPyrDenoiseParams & dnparams); void dirpyr_ablevel(LabImage * data_fine, LabImage * data_coarse, int width, int height, \ LUTf &rangefn_L, LUTf &rangefn_ab, int level, int scale); + void ImStats(float* src, float* dst, int H, int W, int box ); + void WaveletDenoise(cplx_wavelet_decomposition &DualTreeCoeffs, float noisevar ); + void WaveletDenoise(wavelet_decomposition &WaveletCoeffs, float noisevar ); + void BiShrink(float * ReCoeffs, float * ImCoeffs, float * ReParents, float * ImParents, \ + int W, int H, int level, int padding, float noisevar); + void BiShrink(float ** WavCoeffs, float ** WavParents, int W, int H, int level, int padding, float noisevar); + void FirstStageWiener(float* ReCoeffs, float* ImCoeffs, float* wiener1, int W, int H, int rad, float noisevar); + void SecondStageWiener(float* ReCoeffs, float* ImCoeffs, float* wiener1, int W, int H, int rad, float noisevar); + void QCoeffs (float* srcre, float* srcim, float* wiener1, float* dst, int rad, int W, int H); + float UniversalThresh(float * HH_Coeffs, int datalen); + // pyramid equalizer void dirpyr_equalizer (float ** src, float ** dst, int srcwidth, int srcheight, const double * mult );//Emil's directional pyramid equalizer diff --git a/rtengine/impulse_denoise.h b/rtengine/impulse_denoise.h index 0d3c8fe7c..52ad84593 100644 --- a/rtengine/impulse_denoise.h +++ b/rtengine/impulse_denoise.h @@ -37,9 +37,7 @@ void ImProcFunctions::impulse_nr (LabImage* lab, double thresh) { int width = lab->W; int height = lab->H; - - float hpfabs, hfnbrave; - + // buffer for the lowpass image float ** lpf = new float *[height]; // buffer for the highpass image @@ -77,10 +75,10 @@ void ImProcFunctions::impulse_nr (LabImage* lab, double thresh) { for (int i=0; i < height; i++) for (int j=0; j < width; j++) { - - hpfabs = fabs(lab->L[i][j]-lpf[i][j]); + float hfnbrave = 0; + float hpfabs = fabs(lab->L[i][j]-lpf[i][j]); //block average of high pass data - for (i1=MAX(0,i-2), hfnbrave=0; i1<=MIN(i+2,height-1); i1++ ) + for (i1=MAX(0,i-2); i1<=MIN(i+2,height-1); i1++ ) for (j1=MAX(0,j-2); j1<=MIN(j+2,width-1); j1++ ) { hfnbrave += fabs(lab->L[i1][j1]-lpf[i1][j1]); } @@ -94,7 +92,7 @@ void ImProcFunctions::impulse_nr (LabImage* lab, double thresh) { if (!impish[i][j]) continue; norm=0.0; wtdsum[0]=wtdsum[1]=wtdsum[2]=0.0; - for (i1=MAX(0,i-2), hfnbrave=0; i1<=MIN(i+2,height-1); i1++ ) + for (i1=MAX(0,i-2); i1<=MIN(i+2,height-1); i1++ ) for (j1=MAX(0,j-2); j1<=MIN(j+2,width-1); j1++ ) { if (i1==i && j1==j) continue; if (impish[i1][j1]) continue;