diff --git a/rtengine/FTblockDNchroma.cc b/rtengine/FTblockDNchroma.cc index 36b276a8d..81d904641 100644 --- a/rtengine/FTblockDNchroma.cc +++ b/rtengine/FTblockDNchroma.cc @@ -85,10 +85,9 @@ namespace rtengine { //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% // gamma transform input channel data - float gam = dnparams.gamma; + float gam = dnparams.gamma-0.5; 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); @@ -113,7 +112,7 @@ namespace rtengine { 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->a[i][j] = 0.5f*(X-Y); dst->b[i][j] = 0.2f*(Y-Z); //Y = 0.05+0.1*((float)rand()/(float)RAND_MAX);//test with random data @@ -131,31 +130,17 @@ namespace rtengine { //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% // gamma transform output channel data - float gam = dnparams.gamma; + float gam = dnparams.gamma-0.5; 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 @@ -167,7 +152,7 @@ namespace rtengine { //input normalized to (0,1) //Y = igamcurveL[ src->L[i][j] ]; Y = src->L[i][j]; - X = (5.0f*(src->a[i][j])) + Y; + X = (2.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); @@ -257,30 +242,27 @@ namespace rtengine { 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); + int datalen = labin->W * labin->H; + wavelet_decomposition Ldecomp(labin->data, labin->W, labin->H, 5/*maxlevels*/, 0/*subsampling*/ ); + wavelet_decomposition adecomp(labin->data+datalen, labin->W, labin->H, 5, 1 );//last args are maxlevels, subsampling + wavelet_decomposition bdecomp(labin->data+2*datalen, labin->W, labin->H, 5, 1 );//last args are maxlevels, subsampling + + float noisevarL = SQR(dnparams.Lamt/25.0f);//TODO: clean up naming confusion about params + float noisevarab = SQR(dnparams.chroma/25.0f); + + //WaveletDenoise(Ldecomp, SQR((float)dnparams.Lamt/25.0f)); + WaveletDenoiseAll(Ldecomp, adecomp, bdecomp, noisevarL, noisevarab); + + Ldecomp.reconstruct(labblur->data); + adecomp.reconstruct(labblur->data+datalen); + bdecomp.reconstruct(labblur->data+2*datalen); + //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 @@ -303,14 +285,14 @@ namespace rtengine { 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)); + 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)); + 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 @@ -361,8 +343,8 @@ namespace rtengine { for (int j=jmin; jL[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 + 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]; } @@ -374,15 +356,15 @@ namespace rtengine { for (int i=0; iL[2*imin-i-1][left+j]-labblur->L[2*imin-i-1][left+j]); - //RLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->a[2*imin-i-1][left+j]-labblur->a[2*imin-i-1][left+j]); - //BLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->b[2*imin-i-1][left+j]-labblur->b[2*imin-i-1][left+j]); + RLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->a[2*imin-i-1][left+j]-labblur->a[2*imin-i-1][left+j]); + BLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->b[2*imin-i-1][left+j]-labblur->b[2*imin-i-1][left+j]); } if (jmin>0) { for (int i=0; iL[2*imin-i-1][2*jmin-j-1]-labblur->L[2*imin-i-1][2*jmin-j-1]); - //RLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->a[2*imin-i-1][2*jmin-j-1]-labblur->a[2*imin-i-1][2*jmin-j-1]); - //BLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->b[2*imin-i-1][2*jmin-j-1]-labblur->b[2*imin-i-1][2*jmin-j-1]); + RLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->a[2*imin-i-1][2*jmin-j-1]-labblur->a[2*imin-i-1][2*jmin-j-1]); + BLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->b[2*imin-i-1][2*jmin-j-1]-labblur->b[2*imin-i-1][2*jmin-j-1]); } } } @@ -391,15 +373,15 @@ namespace rtengine { for (int i=imax; iL[height+imax-i-1][left+j]-labblur->L[height+imax-i-1][left+j]); - //RLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->a[height+imax-i-1][left+j]-labblur->a[height+imax-i-1][left+j]); - //BLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->b[height+imax-i-1][left+j]-labblur->b[height+imax-i-1][left+j]); + RLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->a[height+imax-i-1][left+j]-labblur->a[height+imax-i-1][left+j]); + BLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->b[height+imax-i-1][left+j]-labblur->b[height+imax-i-1][left+j]); } if (jmaxL[height+imax-i-1][width+jmax-j-1]-labblur->L[height+imax-i-1][width+jmax-j-1]); - //RLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->a[height+imax-i-1][width+jmax-j-1]-labblur->a[height+imax-i-1][width+jmax-j-1]); - //BLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->b[height+imax-i-1][width+jmax-j-1]-labblur->b[height+imax-i-1][width+jmax-j-1]); + RLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->a[height+imax-i-1][width+jmax-j-1]-labblur->a[height+imax-i-1][width+jmax-j-1]); + BLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->b[height+imax-i-1][width+jmax-j-1]-labblur->b[height+imax-i-1][width+jmax-j-1]); } } } @@ -408,15 +390,15 @@ namespace rtengine { for (int j=0; jL[top+i][2*jmin-j-1]-labblur->L[top+i][2*jmin-j-1]); - //RLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->a[top+i][2*jmin-j-1]-labblur->a[top+i][2*jmin-j-1]); - //BLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->b[top+i][2*jmin-j-1]-labblur->b[top+i][2*jmin-j-1]); + RLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->a[top+i][2*jmin-j-1]-labblur->a[top+i][2*jmin-j-1]); + BLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->b[top+i][2*jmin-j-1]-labblur->b[top+i][2*jmin-j-1]); } if (imaxL[height+imax-i-1][2*jmin-j-1]-labblur->L[height+imax-i-1][2*jmin-j-1]); - //RLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->a[height+imax-i-1][2*jmin-j-1]-labblur->a[height+imax-i-1][2*jmin-j-1]); - //BLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->b[height+imax-i-1][2*jmin-j-1]-labblur->b[height+imax-i-1][2*jmin-j-1]); + RLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->a[height+imax-i-1][2*jmin-j-1]-labblur->a[height+imax-i-1][2*jmin-j-1]); + BLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->b[height+imax-i-1][2*jmin-j-1]-labblur->b[height+imax-i-1][2*jmin-j-1]); } } } @@ -426,21 +408,21 @@ namespace rtengine { for (int j=jmax; jL[top+i][width+jmax-j-1]-labblur->L[top+i][width+jmax-j-1]); - //RLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->a[top+i][width+jmax-j-1]-labblur->a[top+i][width+jmax-j-1]); - //BLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->b[top+i][width+jmax-j-1]-labblur->b[top+i][width+jmax-j-1]); + RLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->a[top+i][width+jmax-j-1]-labblur->a[top+i][width+jmax-j-1]); + BLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->b[top+i][width+jmax-j-1]-labblur->b[top+i][width+jmax-j-1]); } if (imin>0) { for (int i=0; iL[2*imin-i-1][width+jmax-j-1]-labblur->L[2*imin-i-1][width+jmax-j-1]); - //RLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->a[2*imin-i-1][width+jmax-j-1]-labblur->a[2*imin-i-1][width+jmax-j-1]); - //BLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->b[2*imin-i-1][width+jmax-j-1]-labblur->b[2*imin-i-1][width+jmax-j-1]); + RLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->a[2*imin-i-1][width+jmax-j-1]-labblur->a[2*imin-i-1][width+jmax-j-1]); + BLblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->b[2*imin-i-1][width+jmax-j-1]-labblur->b[2*imin-i-1][width+jmax-j-1]); } } } //end of tile padding //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - Lblox[vblkmod][(indx + TS/2)*TS+TS/2]=32768.0f;//testing: locate block center + //Lblox[vblkmod][(indx + TS/2)*TS+TS/2]=32768.0f;//testing: locate block centers }//end of filling block row @@ -448,8 +430,8 @@ namespace rtengine { fftwf_execute_dft_r2c(plan_forward_blox,Lblox[vblkmod],fLblox[vblkmod]); // FT an entire row of tiles //fftwf_execute_r2r(plan_forward_blox,Lblox[vblkmod],fLblox[vblkmod]); // DCT an entire row of tiles - //fftwf_execute_dft_r2c(plan_forward_blox,RLblox[vblkmod],fRLblox[vblkmod]);// FT an entire row of tiles - //fftwf_execute_dft_r2c(plan_forward_blox,BLblox[vblkmod],fBLblox[vblkmod]);// FT an entire row of tiles + fftwf_execute_dft_r2c(plan_forward_blox,RLblox[vblkmod],fRLblox[vblkmod]);// FT an entire row of tiles + fftwf_execute_dft_r2c(plan_forward_blox,BLblox[vblkmod],fBLblox[vblkmod]);// FT an entire row of tiles if (vblkL[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->a[i][j]/totwt[i][j]; + labdn->a[i][j] = labblur->a[i][j] + hpdn; + hpdn = labdn->b[i][j]/totwt[i][j]; + labdn->b[i][j] = labblur->b[i][j] + hpdn; - //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]; + //labdn->a[i][j] = labblur->a[i][j]; + //labdn->b[i][j] = labblur->b[i][j]; } } @@ -641,11 +619,11 @@ namespace rtengine { //float Lcoeff = fLblox[vblprocmod][(hblproc*TS+i)*TS+j]; //for DCT - //float RLcoeffre = fRLblox[vblprocmod][(hblproc*TS+i)*fTS+j][0]; - //float RLcoeffim = fRLblox[vblprocmod][(hblproc*TS+i)*fTS+j][1]; + float RLcoeffre = fRLblox[vblprocmod][(hblproc*TS+i)*fTS+j][0]; + float RLcoeffim = fRLblox[vblprocmod][(hblproc*TS+i)*fTS+j][1]; - //float BLcoeffre = fBLblox[vblprocmod][(hblproc*TS+i)*fTS+j][0]; - //float BLcoeffim = fBLblox[vblprocmod][(hblproc*TS+i)*fTS+j][1]; + float BLcoeffre = fBLblox[vblprocmod][(hblproc*TS+i)*fTS+j][0]; + float BLcoeffim = fBLblox[vblprocmod][(hblproc*TS+i)*fTS+j][1]; /*double nbrave=0, nbrsqave=0, coeffsq; int vblnbrmod, hblnbrmod; @@ -664,20 +642,20 @@ namespace rtengine { 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); + 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 freqfactor = 1.0f-MAX((expf(-(SQR(i)+SQR(j))/cutoffsq)),(expf(-(SQR(TS-i)+SQR(j))/cutoffsq))); float Lfactor = 1;//freqfactor;//*(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 RLfactor = 1;//(2*RLblockvar)/(eps+RLwsq+RLblockvar); + float BLfactor = 1;//(2*BLblockvar)/(eps+BLwsq+BLblockvar); + float Lshrinkfactor = SQR(Lwsq/(Lwsq + noisevar_L * Lfactor*exp(-Lwsq/(3*noisevar_L)))); - //float RLshrinkfactor = RLwsq/(RLwsq+noisevar_ab*RLfactor); - //float BLshrinkfactor = BLwsq/(BLwsq+noisevar_ab*BLfactor); + float RLshrinkfactor = RLwsq/(RLwsq+noisevar_ab*RLfactor*exp(-Lwsq/(3*noisevar_L))); + float BLshrinkfactor = BLwsq/(BLwsq+noisevar_ab*BLfactor*exp(-Lwsq/(3*noisevar_L))); //float shrinkfactor = (wsqL[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); + 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); } } @@ -1022,7 +1000,7 @@ void ImProcFunctions::FixImpulse_ab(LabImage * src, LabImage * dst, double radiu float threshsq = noisevar;//*SQR(UniversalThresh(WaveletCoeffs.level_coeffs(lvl)[3], Wlvl*Hlvl)); - BiShrink(WavCoeffs, WavParents, Wlvl, Hlvl, lvl, ParentPadding, threshsq/*noisevar*/); + Shrink(WavCoeffs, Wlvl, Hlvl, lvl, threshsq/*noisevar*/); } } @@ -1065,15 +1043,10 @@ void ImProcFunctions::FixImpulse_ab(LabImage * src, LabImage * dst, double radiu //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - void ImProcFunctions::BiShrink(float ** WavCoeffs, float ** WavParents, int W, int H, int level, int pad, float noisevar) + void ImProcFunctions::Shrink(float ** WavCoeffs, int W, int H, int level, float noisevar) { - //bivariate shrinkage of Sendur & Selesnick + //simple wavelet shrinkage float * sigma = new float[W*H]; - int rad = 256;//3*(1<* buffer = new AlignedBuffer (W*H); - float* temp = buffer->data; - - if (rad==0) { - for (int row=0; row + * + * 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 . + */ +#ifndef _BOXBLUR_H_ +#define _BOXBLUR_H_ + +#include +#include +#include +#include "alignedbuffer.h" + +#ifdef _OPENMP +#include +#endif + +#define SQR(x) ((x)*(x)) + +// classical filtering if the support window is small: + +template void boxblur (T** src, A** dst, int radx, int rady, int W, int H) { + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + //box blur image; box range = (radx,rady) + + AlignedBuffer* buffer = new AlignedBuffer (W*H); + float* temp = buffer->data; + + if (radx==0) { + for (int row=0; row void boxblur (T* src, A* dst, int radx, int rady, int W, int H) { + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + //box blur image; box range = (radx,rady) i.e. box size is (2*radx+1)x(2*rady+1) + + AlignedBuffer* buffer = new AlignedBuffer (W*H); + float* temp = buffer->data; + + if (radx==0) { + for (int row=0; row void boxvar (T* src, T* dst, int radx, int rady, int W, int H) { + + AlignedBuffer* buffer1 = new AlignedBuffer (W*H); + AlignedBuffer* buffer2 = new AlignedBuffer (W*H); + float* tempave = buffer1->data; + float* tempsqave = buffer2->data; + + //float *tempave2 = new float[H]; + AlignedBuffer* buffer3 = new AlignedBuffer (H); + float* tempave2 = buffer3->data; + + //float image_ave = 0; + + //box blur image channel; box size = 2*box+1 + //horizontal blur +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (int row = 0; row < H; row++) { + int len = radx + 1; + tempave[row*W+0] = src[row*W+0]/len; + tempsqave[row*W+0] = SQR(src[row*W+0])/len; + for (int j=1; j<=radx; j++) { + tempave[row*W+0] += src[row*W+j]/len; + tempsqave[row*W+0] += SQR(src[row*W+j])/len; + } + for (int col=1; col<=radx; col++) { + tempave[row*W+col] = (tempave[row*W+col-1]*len + src[row*W+col+radx])/(len+1); + tempsqave[row*W+col] = (tempsqave[row*W+col-1]*len + SQR(src[row*W+col+radx]))/(len+1); + len ++; + } + for (int col = radx+1; col < W-radx; col++) { + tempave[row*W+col] = tempave[row*W+col-1] + (src[row*W+col+radx] - src[row*W+col-radx-1])/len; + tempsqave[row*W+col] = tempsqave[row*W+col-1] + (SQR(src[row*W+col+radx]) - SQR(src[row*W+col-radx-1]))/len; + } + for (int col=W-radx; col void boxdev (T* src, T* dst, int radx, int rady, int W, int H) { + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + //box blur image; box range = (radx,rady) i.e. box size is (2*radx+1)x(2*rady+1) + + AlignedBuffer* buffer1 = new AlignedBuffer (W*H); + float* temp = buffer1->data; + + AlignedBuffer* buffer2 = new AlignedBuffer (W*H); + float* tempave = buffer2->data; + + if (radx==0) { + for (int row=0; row void boxsqblur (T* src, A* dst, int radx, int rady, int W, int H) { + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + //box blur image; box range = (radx,rady) i.e. box size is (2*radx+1)x(2*rady+1) + + AlignedBuffer* buffer = new AlignedBuffer (W*H); + float* temp = buffer->data; + + if (radx==0) { + for (int row=0; row - cplx_wavelet_decomposition(E * src, int width, int height, int maxlvl); + cplx_wavelet_decomposition(E * src, int width, int height, int maxlvl, int subsampling); ~cplx_wavelet_decomposition(); @@ -122,11 +121,9 @@ public: // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% template - cplx_wavelet_decomposition::cplx_wavelet_decomposition(E * src, int width, int height, int maxlvl) - : lvltot(0), m_w(width), m_h(height), m_w1(0), m_h1(0) + cplx_wavelet_decomposition::cplx_wavelet_decomposition(E * src, int width, int height, int maxlvl, int subsampling) + : lvltot(0), subsamp(subsampling), m_w(width), m_h(height) { - m_w1 = width; - m_h1 = height; //initialize wavelet filters @@ -187,11 +184,11 @@ public: for (int m=0; m<2; m++) { 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, \ + dual_tree[0][2*n+m] = new wavelet_level(src, lvltot, subsamp, 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); while(lvltot < maxlvl) { 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][2*n+m] = new wavelet_level(dual_tree[lvltot-1][2*n+m]->lopass()/*lopass*/, lvltot, subsamp, 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, \ @@ -305,9 +302,8 @@ public: static const int maxlevels = 8;//should be greater than any conceivable order of decimation - int lvltot; + int lvltot, subsamp; size_t m_w, m_h;//dimensions - size_t m_w1, m_h1; int wavfilt_len, wavfilt_offset; float *wavfilt_anal; @@ -322,7 +318,7 @@ public: public: template - wavelet_decomposition(E * src, int width, int height, int maxlvl); + wavelet_decomposition(E * src, int width, int height, int maxlvl, int subsampling); ~wavelet_decomposition(); @@ -346,11 +342,21 @@ public: return wavelet_decomp[level]->padding(); } + int level_stride(int level) const + { + return wavelet_decomp[level]->stride(); + } + int maxlevel() const { return lvltot; } + int subsample() const + { + return subsamp; + } + template void reconstruct(E * dst); @@ -361,15 +367,12 @@ public: // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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) + wavelet_decomposition::wavelet_decomposition(E * src, int width, int height, int maxlvl, int subsampling) + : lvltot(0), subsamp(subsampling), m_w(width), m_h(height) { - m_w1 = width; - m_h1 = height; //initialize wavelet filters - wavfilt_len = Haar_len; wavfilt_offset = Haar_offset; wavfilt_anal = new float[2*wavfilt_len]; @@ -396,11 +399,11 @@ public: int padding = 0;//pow(2, maxlvl);//must be a multiple of two lvltot=0; - wavelet_decomp[lvltot] = new wavelet_level(src, lvltot/*level*/, padding/*padding*/, m_w, m_h, \ + wavelet_decomp[lvltot] = new wavelet_level(src, lvltot/*level*/, subsamp, 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] = new wavelet_level(wavelet_decomp[lvltot-1]->lopass()/*lopass*/, lvltot/*level*/, subsamp, 0/*no padding*/, \ wavelet_decomp[lvltot-1]->width(), wavelet_decomp[lvltot-1]->height(), \ wavfilt_anal, wavfilt_anal, wavfilt_len, wavfilt_offset); } @@ -417,15 +420,15 @@ public: // data structure is wavcoeffs[scale][channel={lo,hi1,hi2,hi3}][pixel_array] - int skip=1<<(lvltot-1); + //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; + wavelet_decomp[lvl]->reconstruct_level(wavelet_decomp[lvl-1]->wavcoeffs[0], wavfilt_synth, wavfilt_synth, wavfilt_len, wavfilt_offset, subsamp); + //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); + wavelet_decomp[0]->reconstruct_level(tmp, wavfilt_synth, wavfilt_synth, wavfilt_len, wavfilt_offset, subsamp); copy_out(tmp,dst,m_w*m_h); diff --git a/rtengine/cplx_wavelet_level.h b/rtengine/cplx_wavelet_level.h index dc2aaca58..943873be4 100644 --- a/rtengine/cplx_wavelet_level.h +++ b/rtengine/cplx_wavelet_level.h @@ -29,6 +29,8 @@ namespace rtengine { +#undef MAX +#undef MIN #define MAX(a,b) ((a) > (b) ? (a) : (b)) #define MIN(a,b) ((a) > (b) ? (b) : (a)) #define SQR(x) ((x)*(x)) @@ -51,6 +53,9 @@ namespace rtengine { // level of decomposition int lvl; + // whether to subsample the output + bool subsamp_out; + // spacing of filter taps size_t skip; @@ -80,6 +85,17 @@ namespace rtengine { 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 AnalysisFilterHaar (T * srcbuffer, T * dstLo, T * dstHi, int pitch, int srclen); + void SynthesisFilterHaar (T * srcLo, T * srcHi, T * dst, T *bufferLo, T *bufferHi, int pitch, int dstlen); + + void AnalysisFilterSubsamp (T * srcbuffer, T * dstLo, T * dstHi, float *filterLo, float *filterHi, + int taps, int offset, int pitch, int srclen); + void SynthesisFilterSubsamp (T * srcLo, T * srcHi, T * dst, T *bufferLo, T *bufferHi, + float *filterLo, float *filterHi, int taps, int offset, int pitch, int dstlen); + + void AnalysisFilterSubsampHaar (T * srcbuffer, T * dstLo, T * dstHi, int pitch, int srclen); + void SynthesisFilterSubsampHaar (T * srcLo, T * srcHi, T * dst, T *bufferLo, T *bufferHi, int pitch, int dstlen); + void imp_nr (T* src, int width, int height, double thresh); @@ -88,11 +104,17 @@ namespace rtengine { T ** wavcoeffs; template - 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<>level)&1) { - m_w2 = (w+2*skip*padding); - m_h2 = (h+2*skip*padding); + if (subsamp) { + skip = 1; + for (int n=0; n>n)&1); + } + } + m_w2 = (subsamp_out ? ((w+1+2*skip*padding)/2) : (w+2*skip*padding)); + m_h2 = (subsamp_out ? ((h+1+2*skip*padding)/2) : (h+2*skip*padding)); m_pad= skip*padding; wavcoeffs = create((m_w2)*(m_h2)); @@ -130,6 +152,11 @@ namespace rtengine { return m_pad/skip; } + size_t stride() const + { + return skip; + } + template void decompose_level(E *src, float *filterV, float *filterH, int len, int offset, int skip); @@ -168,39 +195,14 @@ namespace rtengine { } // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + template template void wavelet_level::loadbuffer(E * src, E * dst, int pitch, int srclen) { E * tmp = dst + m_pad; memset(dst, 0, (MAX(m_w2,m_h2))*sizeof(E)); - /*int cosetlen = (srclen+1)/skip; - - //create buffer with 'skip' rows and 'cosetlen' columns from src data - //'skip' is the spacing of taps on the wavelet filter to be applied to src rows/columns - //therefore there are 'skip' cosets of the row/column data, each of length 'cosetlen' - //'pitch' is 1 for rows, W for columns - for (size_t i = 0, j = 0; itaps && i 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, - * aligning the 'offset' element of the filter - * with the input pixel, and skipping 'pitch' pixels - * between taps (eg pitch=1 for horizontal filtering, - * pitch=W for vertical, pitch=W+1,W-1 for diagonals. - * Currently diagonal filtering is not supported - * for the full source array, until a more sophisticated - * treatment of mirror BC's is implemented. + * Applies an FIR filter 'filter' with filter length 'taps', + * aligning the 'offset' element of the filter with + * the input pixel, and skipping 'skip' pixels between taps * */ - - // load into buffer - /* - int srclen=(dstlen+(dstlen%skip)+2*m_pad); //length of row/col in src (coarser level) - int cosetlen = srclen/skip; //length of coset (skip is spacing of taps in filter) - - for (size_t i=0, j=0; itaps && i<(cosetlen-taps)) {//bulk - for (int j=0, l=-shift; jskip*taps && i<(srclen-skip*taps)) {//bulk - for (int j=0, l=-skip*shift; j65535.0f) { - float xxx=tot; - float yyy=1.0f; - } + } } + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + template + void wavelet_level::AnalysisFilterHaar (T * srcbuffer, T * dstLo, T * dstHi, int pitch, int srclen) { + + /* Basic convolution code + * Applies a Haar filter + * + */ + + for(size_t i = 0; i < (srclen - skip); i++) { + dstLo[(pitch*(i))] = 0.5*(srcbuffer[i] + srcbuffer[i+skip]); + dstHi[(pitch*(i))] = 0.5*(srcbuffer[i] - srcbuffer[i+skip]); + } + + for(size_t i = (srclen-skip); i < (srclen); i++) { + dstLo[(pitch*(i))] = 0.5*(srcbuffer[i] + srcbuffer[i-skip]); + dstHi[(pitch*(i))] = 0.5*(srcbuffer[i] - srcbuffer[i-skip]); + } + + } + + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + template + void wavelet_level::SynthesisFilterHaar (T * srcLo, T * srcHi, T * dst, T *bufferLo, T *bufferHi, int pitch, int dstlen) { + + /* Basic convolution code + * Applies a Haar filter + * + */ + + int srclen = (dstlen==m_w ? m_w2 : m_h2);//length of row/col in src (coarser level) + + for (size_t i=0, j=0; i template 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_w2,m_h2)]; + T *buffer = new T[MAX(m_w,m_h)]; /* filter along columns */ for (int j=0; j template void wavelet_level::reconstruct_level(E *dst, float *filterV, float *filterH, int taps, int offset, int skip) { @@ -415,16 +415,28 @@ namespace rtengine { /* filter along rows */ for (int i=0; i - void wavelet_level::imp_nr (T* src, int width, int height, double thresh) { + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + template + void wavelet_level::AnalysisFilterSubsamp (T * srcbuffer, T * dstLo, T * dstHi, float *filterLo, float *filterHi, + int taps, int offset, int pitch, int srclen) { + /* Basic convolution code + * Applies an FIR filter 'filter' with filter length 'taps', + * aligning the 'offset' element of the filter with + * the input pixel, and skipping 'skip' pixels between taps + * Output is subsampled by two + */ - // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - // impulse noise removal - // local variables + // calculate coefficients - float hpfabs, hfnbrave; - const float eps = 0.01; + for(int i = 0; i < (srclen); i+=2) { + float lo=0,hi=0; + if (i>skip*taps && i* 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; - } + + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + template + void wavelet_level::SynthesisFilterSubsamp (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 filter length 'taps', + * aligning the 'offset' element of the filter with + * the input pixel, and skipping 'skip' pixels between taps + * Output is subsampled by two + */ + + + + // calculate coefficients + + int srclen = (dstlen==m_w ? m_w2 : m_h2);//length of row/col in src (coarser level) + + //fill a buffer with a given row/column of data + for (size_t i=0, j=0; iskip*taps && i<(srclen-skip*taps)) {//bulk + for (int j=begin, l=0; j + void wavelet_level::AnalysisFilterSubsampHaar (T * srcbuffer, T * dstLo, T * dstHi, int pitch, int srclen) { + + /* Basic convolution code + * Applies a Haar filter + * Output is subsampled by two + */ + + // calculate coefficients + + for(size_t i = 0; i < (srclen - skip); i+=2) { + dstLo[(pitch*(i/2))] = 0.5*(srcbuffer[i] + srcbuffer[i+skip]); + dstHi[(pitch*(i/2))] = 0.5*(srcbuffer[i] - srcbuffer[i+skip]); + } + + for(size_t i = (srclen-skip)-(srclen-skip)&1; i < (srclen); i+=2) { + dstLo[(pitch*(i/2))] = 0.5*(srcbuffer[i] + srcbuffer[i-skip]); + dstHi[(pitch*(i/2))] = 0.5*(srcbuffer[i] - srcbuffer[i-skip]); + } + + } + + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + template + void wavelet_level::SynthesisFilterSubsampHaar (T * srcLo, T * srcHi, T * dst, T *bufferLo, T *bufferHi, int pitch, int dstlen) { + + /* Basic convolution code + * Applies a Haar filter + * Output was subsampled by two + */ + + + // calculate coefficients + + for(size_t i = m_pad; i < (dstlen+m_pad-skip); i+=2) { + dst[pitch*(i-m_pad)] = bufferLo[i/2]+bufferHi[i/2]; + dst[pitch*(i+skip-m_pad)] = bufferLo[i/2]-bufferHi[i/2]; + } + + if ((dstlen+m_pad)&1) { + dst[pitch*(dstlen-1)] = bufferLo[(dstlen+m_pad-1)/2]+bufferHi[(dstlen+m_pad-1)/2]; + } + + } + + + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + }; diff --git a/rtengine/improcfun.cc b/rtengine/improcfun.cc index d91ddda00..7a43e78d3 100644 --- a/rtengine/improcfun.cc +++ b/rtengine/improcfun.cc @@ -561,14 +561,22 @@ void ImProcFunctions::colorCurve (LabImage* lold, LabImage* lnew) { if (params->impulseDenoise.enabled && lab->W>=8 && lab->H>=8) { //impulse_nr (lab, (float)params->impulseDenoise.thresh/10.0 );//20 is normal - - for (int i=0; iW*lab->H; i++) { + int datalen = lab->W*lab->H; + for (int i=0; idata[i] *= lab->data[i]/32768.0f; + //lab->data[i] = 5000; } - 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)); - + //lab->data[100*lab->W+100] = 20000; + wavelet_decomposition Ldecomp(lab->data, lab->W, lab->H, 5/*maxlevel*/, 1/*subsampling*/ ); + wavelet_decomposition adecomp(lab->data+datalen, lab->W, lab->H, 5, 1 );//last args are maxlevels, subsampling + wavelet_decomposition bdecomp(lab->data+2*datalen, lab->W, lab->H, 5, 1 );//last args are maxlevels, subsampling + + float noisevar_L = SQR((float)params->impulseDenoise.thresh/25.0f); + float noisevar_ab = SQR((float)params->defringe.threshold / 5.0f); + + //WaveletDenoise(Ldecomp, SQR((float)params->impulseDenoise.thresh/25.0f)); + WaveletDenoiseAll(Ldecomp, adecomp, bdecomp, noisevar_L, noisevar_ab); + LabImage* labtmp = new LabImage (lab->W,lab->H); int lvl = (params->impulseDenoise.thresh>>4)&7; @@ -576,17 +584,20 @@ void ImProcFunctions::colorCurve (LabImage* lold, LabImage* lnew) { 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); - }*/ + //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); - + adecomp.reconstruct(lab->data+datalen); + bdecomp.reconstruct(lab->data+2*datalen); + //double radius = (int)(params->impulseDenoise.thresh/10) ; //boxvar(lab->data, lab->data, radius, radius, lab->W, lab->H); @@ -612,13 +623,13 @@ void ImProcFunctions::colorCurve (LabImage* lold, LabImage* lnew) { //impulse_nr (labtmp, 50.0f/20.0f); - for (int i=0; iW*lab->H; i++) { + for (int i=0; idata[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); + //impulse_nr (lab, 50.0f/20.0f); } } diff --git a/rtengine/improcfun.h b/rtengine/improcfun.h index 4ff124add..be5a92e97 100644 --- a/rtengine/improcfun.h +++ b/rtengine/improcfun.h @@ -161,14 +161,14 @@ namespace rtengine { 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 WaveletDenoiseAll(wavelet_decomposition &WaveletCoeffs_L, wavelet_decomposition &WaveletCoeffs_a, + wavelet_decomposition &WaveletCoeffs_b, float noisevar_L, float noisevar_ab ); 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); + void Shrink(float ** WavCoeffs, int W, int H, int level, float noisevar); + void ShrinkAll(float ** WavCoeffs_L, float ** WavCoeffs_a, float ** WavCoeffs_b, int level, \ + int W_L, int H_L, int W_ab, int H_ab, int skip_L, int skip_ab, float noisevar_L, float noisevar_ab); 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