liz/rawTherapee
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Adding boxblur.h.

Browse Source
This commit is contained in:
Emil Martinec
2012-02-12 09:58:30 -06:00
parent c492509175
commit 17d9483495
6 changed files with 983 additions and 482 deletions
Download Patch File Download Diff File Expand all files Collapse all files

427
rtengine/FTblockDNchroma.cc
View File

@@ -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; j<jmax; j++) {
Lblox[vblkmod][(indx + i)*TS+j] = tilemask_in[i][j]*(labin->L[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; i<imin; i++)
for (int j=jmin; j<jmax; j++) {
Lblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->L[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; i<imin; i++)
for (int j=0; j<jmin; j++) {
Lblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->L[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; i<TS; i++)
for (int j=jmin; j<jmax; j++) {
Lblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->L[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 (jmax<TS) {
for (int i=imax; i<TS; i++)
for (int j=jmax; j<TS; j++) {
Lblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->L[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; j<jmin; j++)
for (int i=imin; i<imax; i++) {
Lblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->L[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 (imax<TS) {
for (int j=0; j<jmin; j++)
for (int i=imax; i<TS; i++) {
Lblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->L[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; j<TS; j++)
for (int i=imin; i<imax; i++) {
Lblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->L[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; i<imin; i++)
for (int j=jmax; j<TS; j++) {
Lblox[vblkmod][(indx + i)*TS+j]=tilemask_in[i][j]*(labin->L[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 (vblk<blkrad) continue;
@@ -481,8 +463,8 @@ namespace rtengine {
fftwf_execute_dft_c2r(plan_backward_blox,fLblox[vblprocmod],Lblox[vblprocmod]); //for FT
//fftwf_execute_r2r(plan_backward_blox,fLblox[vblprocmod],Lblox[vblprocmod]); //for DCT
//fftwf_execute_dft_c2r(plan_backward_blox,fRLblox[vblprocmod],RLblox[vblprocmod]);
//fftwf_execute_dft_c2r(plan_backward_blox,fBLblox[vblprocmod],BLblox[vblprocmod]);
fftwf_execute_dft_c2r(plan_backward_blox,fRLblox[vblprocmod],RLblox[vblprocmod]);
fftwf_execute_dft_c2r(plan_backward_blox,fBLblox[vblprocmod],BLblox[vblprocmod]);
int topproc = (vblproc-blkrad)*offset;
@@ -499,8 +481,8 @@ namespace rtengine {
fftwf_execute_dft_c2r(plan_backward_blox,fLblox[vblprocmod],Lblox[vblprocmod]); //for FT
//fftwf_execute_r2r(plan_backward_blox,fLblox[vblprocmod],Lblox[vblprocmod]); //for DCT
//fftwf_execute_dft_c2r(plan_backward_blox,fRLblox[vblprocmod],RLblox[vblprocmod]);
//fftwf_execute_dft_c2r(plan_backward_blox,fBLblox[vblprocmod],BLblox[vblprocmod]);
fftwf_execute_dft_c2r(plan_backward_blox,fRLblox[vblprocmod],RLblox[vblprocmod]);
fftwf_execute_dft_c2r(plan_backward_blox,fBLblox[vblprocmod],BLblox[vblprocmod]);
RGBoutput_tile_row (Lblox[vblprocmod], RLblox[vblprocmod], BLblox[vblprocmod], labdn, \
tilemask_out, height, width, topproc, blkrad );
@@ -522,18 +504,18 @@ namespace rtengine {
for( int i = 0 ; i < 8 ; i++ ) {
fftwf_free ( Lblox[i]);
//fftwf_free (RLblox[i]);
//fftwf_free (BLblox[i]);
fftwf_free (RLblox[i]);
fftwf_free (BLblox[i]);
fftwf_free ( fLblox[i]);
//fftwf_free (fRLblox[i]);
//fftwf_free (fBLblox[i]);
fftwf_free (fRLblox[i]);
fftwf_free (fBLblox[i]);
}
delete[] Lblox;
//delete[] RLblox;
//delete[] BLblox;
delete[] RLblox;
delete[] BLblox;
delete[] fLblox;
//delete[] fRLblox;
//delete[] fBLblox;
delete[] fRLblox;
delete[] fBLblox;
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -545,18 +527,14 @@ namespace rtengine {
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->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 = (wsq<noisevar ? 0 : 1);//hard threshold
@@ -685,11 +663,11 @@ namespace rtengine {
fLblox[vblprocmod][(hblproc*TS+i)*fTS+j][1] *= Lshrinkfactor;
//fLblox[vblprocmod][(hblproc*TS+i)*TS+j] *= Lshrinkfactor; //for DCT
//fRLblox[vblprocmod][(hblproc*TS+i)*fTS+j][0] *= Lshrinkfactor;
//fRLblox[vblprocmod][(hblproc*TS+i)*fTS+j][1] *= Lshrinkfactor;
fRLblox[vblprocmod][(hblproc*TS+i)*fTS+j][0] *= RLshrinkfactor;
fRLblox[vblprocmod][(hblproc*TS+i)*fTS+j][1] *= RLshrinkfactor;
//fBLblox[vblprocmod][(hblproc*TS+i)*fTS+j][0] *= Lshrinkfactor;
//fBLblox[vblprocmod][(hblproc*TS+i)*fTS+j][1] *= Lshrinkfactor;
fBLblox[vblprocmod][(hblproc*TS+i)*fTS+j][0] *= BLshrinkfactor;
fBLblox[vblprocmod][(hblproc*TS+i)*fTS+j][1] *= BLshrinkfactor;
}//end of block denoise
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -722,8 +700,8 @@ namespace rtengine {
for (int j=jmin; j<jmax; j++) {
labdn->L[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<<level);
//boxdev(WavCoeffs[3],sigma,rad,rad,W,H);//box blur detail coeffs to estimate local variance
const float root3 = sqrt(3);
const int Wpar = (W+2*pad);
//const int Hpar = (H+1+2*pad)/2;
const float eps = 0.01f;
printf("level=%d ",level);
for (int dir=1; dir<4; dir++) {
@@ -1082,12 +1055,7 @@ void ImProcFunctions::FixImpulse_ab(LabImage * src, LabImage * dst, double radiu
for (int i=0; i<H; i++) {
for (int j=0; j<W; j++) {
//float thresh = root3 * noisevar/sqrt(MAX(sigma[i*W+j]-noisevar, eps));
//int parentloc = ((i)-pad)*Wpar+(j)-pad;
int coeffloc = i*W+j;
//float mag = (SQR(WavCoeffs[dir][coeffloc]) + SQR(WavParents[dir][parentloc]) );
//float shrinkfactor = MAX(0,mag-thresh);
//shrinkfactor /= (shrinkfactor+thresh+eps);
float mag = SQR(WavCoeffs[dir][coeffloc]);
float shrinkfactor = mag/(mag+noisevar*mad*exp(-mag/(3*noisevar*mad))+eps);
//float shrinkfactor = mag/(mag+noisevar*SQR(sigma[coeffloc])+eps);
@@ -1097,15 +1065,17 @@ void ImProcFunctions::FixImpulse_ab(LabImage * src, LabImage * dst, double radiu
}
}
boxblur(sigma, sigma, 1, 1, W, H);
boxblur(sigma, sigma, 1, 1, W, H);//increase smoothness by locally averaging shrinkage
for (int i=0; i<W*H; i++) {
float mag = SQR(WavCoeffs[dir][i]);
float sf = mag/(mag+noisevar*mad+eps);
//float sf1 = mag/(mag+4*noisevar*mad+eps);
//use smoothed shrinkage unless local shrinkage is much less
WavCoeffs[dir][i] *= (SQR(sigma[i])+SQR(sf))/(sigma[i]+sf+eps);
//the following is for testing
//float wdn = WavCoeffs[dir][i]*sf;
//float sf1 = mag/(mag+4*noisevar*mad+eps);
//float wdn1 = WavCoeffs[dir][i]*sf1;
//WavCoeffs[dir][i] = wdn-wdn1;
@@ -1115,132 +1085,6 @@ void ImProcFunctions::FixImpulse_ab(LabImage * src, LabImage * dst, double radiu
delete[] sigma;
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
void ImProcFunctions::FirstStageWiener(float* ReCoeffs, float* ImCoeffs, float* wiener1, int W, int H, int rad, float noisevar) {
//wiener filter
float * sigmare = new float[W*H];
float * sigmaim = new float[W*H];
//box blur detail coeffs to estimate local variance
boxvar(ReCoeffs,sigmare,rad,rad,W,H);
boxvar(ImCoeffs,sigmaim,rad,rad,W,H);
const float eps = 0.01f;
for (int i=0; i<H; i++) {
for (int j=0; j<W; j++) {
//classic Wiener filter
float var = sigmare[i*W+j]+sigmaim[i*W+j];
wiener1[i*W+j] = MAX(var-noisevar, eps)/(var+eps);
}
}
delete[] sigmare;
delete[] sigmaim;
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
void ImProcFunctions::SecondStageWiener(float* ReCoeffs, float* ImCoeffs, float* wiener1, int W, int H, int rad, float noisevar) {
//wiener filter
float * Q = new float[W*H];
float * C = new float[W*H];
//compute matched filter and neighborhood wiener ave
QCoeffs(ReCoeffs, ImCoeffs, wiener1, Q, rad, W, H);
boxblur(wiener1, C, rad, rad, W, H);
const float eps = 0.01f;
for (int i=0; i<H; i++) {
for (int j=0; j<W; j++) {
//classic Wiener filter
float wiener2 = MAX(C[i*W+j]*noisevar,Q[i*W+j]-C[i*W+j]*noisevar)/(Q[i*W+j]+(1-C[i*W+j])*noisevar+eps);
ReCoeffs[i*W+j] *= wiener2;
ImCoeffs[i*W+j] *= wiener2;
}
}
delete[] Q;
delete[] C;
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
void ImProcFunctions::QCoeffs (float* srcre, float* srcim, float* wiener1, float* dst, int rad, int W, int H) {
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//box blur image; box size is (2*rad+1)x(2*rad+1)
AlignedBuffer<float>* buffer = new AlignedBuffer<float> (W*H);
float* temp = buffer->data;
if (rad==0) {
for (int row=0; row<H; row++)
for (int col=0; col<H; col++) {
temp[row*H+col] = wiener1[row*W+col]*(SQR(srcre[row*W+col])+SQR(srcim[row*W+col]));
}
} else {
//horizontal blur
for (int row = 0; row < H; row++) {
int len = rad + 1;
temp[row*W+0] = wiener1[row*W+0]*(SQR(srcre[row*W+0])+SQR(srcim[row*W+0]))/len;
for (int j=1; j<=rad; j++) {
temp[row*W+0] += wiener1[row*W+j]*(SQR(srcre[row*W+j])+SQR(srcim[row*W+j]))/len;
}
for (int col=1; col<=rad; col++) {
temp[row*W+col] = (temp[row*W+col-1]*len +
wiener1[row*W+col+rad]*(SQR(srcre[row*W+col+rad])+SQR(srcim[row*W+col+rad])))/(len+1);
len ++;
}
for (int col = rad+1; col < W-rad; col++) {
temp[row*W+col] = temp[row*W+col-1] + (wiener1[row*W+col+rad]*(SQR(srcre[row*W+col+rad])+SQR(srcim[row*W+col+rad])) - \
wiener1[row*W+col-rad-1]*(SQR(srcre[row*W+col-rad-1])+SQR(srcim[row*W+col-rad-1])))/len;
}
for (int col=W-rad; col<W; col++) {
temp[row*W+col] = (temp[row*W+col-1]*len - wiener1[row*W+col-rad-1]*(SQR(srcre[row*W+col-rad-1])+SQR(srcim[row*W+col-rad-1])))/(len-1);
len --;
}
}
}
if (rad==0) {
for (int row=0; row<H; row++)
for (int col=0; col<H; col++) {
dst[row*W+col] = temp[row*W+col];
}
} else {
//vertical blur
for (int col = 0; col < W; col++) {
int len = rad + 1;
dst[0*W+col] = temp[0*W+col]/len;
for (int i=1; i<=rad; i++) {
dst[0*W+col] += temp[i*W+col]/len;
}
for (int row=1; row<=rad; row++) {
dst[row*W+col] = (dst[(row-1)*W+col]*len + temp[(row+rad)*W+col])/(len+1);
len ++;
}
for (int row = rad+1; row < H-rad; row++) {
dst[row*W+col] = dst[(row-1)*W+col] + (temp[(row+rad)*W+col] - temp[(row-rad-1)*W+col])/len;
}
for (int row=H-rad; row<H; row++) {
dst[row*W+col] = (dst[(row-1)*W+col]*len - temp[(row-rad-1)*W+col])/(len-1);
len --;
}
}
}
delete buffer;
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -1274,5 +1118,104 @@ void ImProcFunctions::FixImpulse_ab(LabImage * src, LabImage * dst, double radiu
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
void ImProcFunctions::WaveletDenoiseAll(wavelet_decomposition &WaveletCoeffs_L, wavelet_decomposition &WaveletCoeffs_a,
wavelet_decomposition &WaveletCoeffs_b, float noisevar_L, float noisevar_ab )
{
int maxlvl = WaveletCoeffs_L.maxlevel();
for (int lvl=0; lvl<maxlvl; lvl++) {
int Wlvl_L = WaveletCoeffs_L.level_W(lvl);
int Hlvl_L = WaveletCoeffs_L.level_H(lvl);
int Wlvl_ab = WaveletCoeffs_a.level_W(lvl);
int Hlvl_ab = WaveletCoeffs_a.level_H(lvl);
float skip_L = WaveletCoeffs_L.level_stride(lvl);
float skip_ab = WaveletCoeffs_a.level_stride(lvl);
float ** WavCoeffs_L = WaveletCoeffs_L.level_coeffs(lvl);
float ** WavCoeffs_a = WaveletCoeffs_a.level_coeffs(lvl);
float ** WavCoeffs_b = WaveletCoeffs_b.level_coeffs(lvl);
ShrinkAll(WavCoeffs_L, WavCoeffs_a, WavCoeffs_b, lvl, Wlvl_L, Hlvl_L, Wlvl_ab, Hlvl_ab,
skip_L, skip_ab, noisevar_L, noisevar_ab);
}
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
void ImProcFunctions::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)
{
//simple wavelet shrinkage
const float eps = 0.01f;
float * sigma = new float[W_L*H_L];
printf("\n level=%d \n",level);
for (int dir=1; dir<4; dir++) {
float mad_L = SQR(UniversalThresh(WavCoeffs_L[dir], W_L*H_L));//*6*/(level+1);
float mad_a = SQR(UniversalThresh(WavCoeffs_a[dir], W_ab*H_ab));//*6*/(level+1);
float mad_b = SQR(UniversalThresh(WavCoeffs_b[dir], W_ab*H_ab));//*6*/(level+1);
printf(" dir=%d mad_L=%f mad_a=%f mad_b=%f \n",dir,sqrt(mad_L),sqrt(mad_a),sqrt(mad_b));
for (int i=0; i<H_ab; i++) {
for (int j=0; j<W_ab; j++) {
int coeffloc_ab = i*W_ab+j;
int coeffloc_L = ((i*skip_L)/skip_ab)*W_L + ((j*skip_L)/skip_ab);
float mag_L = fabs(WavCoeffs_L[dir][coeffloc_L ])+eps;
float mag_a = SQR(WavCoeffs_a[dir][coeffloc_ab])+eps;
float mag_b = SQR(WavCoeffs_b[dir][coeffloc_ab])+eps;
float edgefactor = exp(-mag_L/(sqrt(noisevar_L)*mad_L)) * exp(-mag_a/(3*noisevar_ab*mad_a)) * exp(-mag_b/(3*noisevar_ab*mad_b));
WavCoeffs_a[dir][coeffloc_ab] *= mag_a/(mag_a + noisevar_ab*mad_a*edgefactor + eps);
WavCoeffs_b[dir][coeffloc_ab] *= mag_b/(mag_b + noisevar_ab*mad_b*edgefactor + eps);
}
}
for (int i=0; i<W_L*H_L; i++) {
float mag = SQR(WavCoeffs_L[dir][i]);
float shrinkfactor = mag/(mag+noisevar_L*mad_L*exp(-mag/(3*noisevar_L*mad_L))+eps);
//float shrinkfactor = mag/(mag+noisevar*SQR(sigma[coeffloc])+eps);
//WavCoeffs[dir][coeffloc] *= shrinkfactor;
sigma[i] = shrinkfactor;
}
boxblur(sigma, sigma, 1, 1, W_L, H_L);//increase smoothness by locally averaging shrinkage
for (int i=0; i<W_L*H_L; i++) {
float mag = SQR(WavCoeffs_L[dir][i]);
float sf = mag/(mag+noisevar_L*mad_L+eps);
//use smoothed shrinkage unless local shrinkage is much less
WavCoeffs_L[dir][i] *= (SQR(sigma[i])+SQR(sf))/(sigma[i]+sf+eps);
//the following is for testing
//float wdn = WavCoeffs[dir][i]*sf;
//float sf1 = mag/(mag+4*noisevar*mad+eps);
//float wdn1 = WavCoeffs[dir][i]*sf1;
//WavCoeffs[dir][i] = wdn-wdn1;
}//now luminance coeffs are denoised
}
delete[] sigma;
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
};

464
rtengine/boxblur.h Normal file
View File

@@ -0,0 +1,464 @@
/*
* This file is part of RawTherapee.
*
* Copyright © 2010 Emil Martinec <ejmartin@uchicago.edu>
*
* 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 <http://www.gnu.org/licenses/>.
*/
#ifndef _BOXBLUR_H_
#define _BOXBLUR_H_
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "alignedbuffer.h"
#ifdef _OPENMP
#include <omp.h>
#endif
#define SQR(x) ((x)*(x))
// classical filtering if the support window is small:
template<class T, class A> void boxblur (T** src, A** dst, int radx, int rady, int W, int H) {
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//box blur image; box range = (radx,rady)
AlignedBuffer<float>* buffer = new AlignedBuffer<float> (W*H);
float* temp = buffer->data;
if (radx==0) {
for (int row=0; row<H; row++)
for (int col=0; col<H; col++) {
temp[row*H+col] = (float)src[row][col];
}
} else {
//horizontal blur
for (int row = 0; row < H; row++) {
int len = radx + 1;
temp[row*W+0] = (float)src[row][0]/len;
for (int j=1; j<=radx; j++) {
temp[row*W+0] += (float)src[row][j]/len;
}
for (int col=1; col<=radx; col++) {
temp[row*W+col] = (temp[row*W+col-1]*len + (float)src[row][col+radx])/(len+1);
len ++;
}
for (int col = radx+1; col < W-radx; col++) {
temp[row*W+col] = temp[row*W+col-1] + ((float)(src[row][col+radx] - src[row][col-radx-1]))/len;
}
for (int col=W-radx; col<W; col++) {
temp[row*W+col] = (temp[row*W+col-1]*len - src[row][col-radx-1])/(len-1);
len --;
}
}
}
if (rady==0) {
for (int row=0; row<H; row++)
for (int col=0; col<H; col++) {
dst[row][col] = temp[row*W+col];
}
} else {
//vertical blur
for (int col = 0; col < W; col++) {
int len = rady + 1;
dst[0][col] = temp[0*W+col]/len;
for (int i=1; i<=rady; i++) {
dst[0][col] += temp[i*W+col]/len;
}
for (int row=1; row<=rady; row++) {
dst[row][col] = (dst[(row-1)][col]*len + temp[(row+rady)*W+col])/(len+1);
len ++;
}
for (int row = rady+1; row < H-rady; row++) {
dst[row][col] = dst[(row-1)][col] + (temp[(row+rady)*W+col] - temp[(row-rady-1)*W+col])/len;
}
for (int row=H-rady; row<H; row++) {
dst[row][col] = (dst[(row-1)][col]*len - temp[(row-rady-1)*W+col])/(len-1);
len --;
}
}
}
delete buffer;
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
template<class T, class A> 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<float>* buffer = new AlignedBuffer<float> (W*H);
float* temp = buffer->data;
if (radx==0) {
for (int row=0; row<H; row++)
for (int col=0; col<H; col++) {
temp[row*H+col] = src[row*W+col];
}
} else {
//horizontal blur
for (int row = 0; row < H; row++) {
int len = radx + 1;
temp[row*W+0] = (float)src[row*W+0]/len;
for (int j=1; j<=radx; j++) {
temp[row*W+0] += (float)src[row*W+j]/len;
}
for (int col=1; col<=radx; col++) {
temp[row*W+col] = (temp[row*W+col-1]*len + src[row*W+col+radx])/(len+1);
len ++;
}
for (int col = radx+1; col < W-radx; col++) {
temp[row*W+col] = temp[row*W+col-1] + ((float)(src[row*W+col+radx] - src[row*W+col-radx-1]))/len;
}
for (int col=W-radx; col<W; col++) {
temp[row*W+col] = (temp[row*W+col-1]*len - src[row*W+col-radx-1])/(len-1);
len --;
}
}
}
if (rady==0) {
for (int row=0; row<H; row++)
for (int col=0; col<H; col++) {
dst[row*W+col] = temp[row*W+col];
}
} else {
//vertical blur
for (int col = 0; col < W; col++) {
int len = rady + 1;
dst[0*W+col] = temp[0*W+col]/len;
for (int i=1; i<=rady; i++) {
dst[0*W+col] += temp[i*W+col]/len;
}
for (int row=1; row<=rady; row++) {
dst[row*W+col] = (dst[(row-1)*W+col]*len + temp[(row+rady)*W+col])/(len+1);
len ++;
}
for (int row = rady+1; row < H-rady; row++) {
dst[row*W+col] = dst[(row-1)*W+col] + (temp[(row+rady)*W+col] - temp[(row-rady-1)*W+col])/len;
}
for (int row=H-rady; row<H; row++) {
dst[row*W+col] = (dst[(row-1)*W+col]*len - temp[(row-rady-1)*W+col])/(len-1);
len --;
}
}
}
delete buffer;
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
template<typename T> void boxvar (T* src, T* dst, int radx, int rady, int W, int H) {
AlignedBuffer<float>* buffer1 = new AlignedBuffer<float> (W*H);
AlignedBuffer<float>* buffer2 = new AlignedBuffer<float> (W*H);
float* tempave = buffer1->data;
float* tempsqave = buffer2->data;
//float *tempave2 = new float[H];
AlignedBuffer<float>* buffer3 = new AlignedBuffer<float> (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<W; col++) {
tempave[row*W+col] = (tempave[row*W+col-1]*len - src[row*W+col-radx-1])/(len-1);
tempsqave[row*W+col] = (tempsqave[row*W+col-1]*len - SQR(src[row*W+col-radx-1]))/(len-1);
len --;
}
}
//vertical blur
#ifdef _OPENMP
#pragma omp parallel for
#endif
for (int col = 0; col < W; col++) {
int len = rady + 1;
tempave2[0] = tempave[0*W+col]/len;
dst[0*W+col] = tempsqave[0*W+col]/len;
for (int i=1; i<=rady; i++) {
tempave2[0] += tempave[i*W+col]/len;
dst[0*W+col] += tempsqave[i*W+col]/len;
}
for (int row=1; row<=rady; row++) {
tempave2[row] = (tempave2[(row-1)]*len + tempave[(row+rady)*W+col])/(len+1);
dst[row*W+col] = (dst[(row-1)*W+col]*len + tempsqave[(row+rady)*W+col])/(len+1);
len ++;
}
for (int row = rady+1; row < H-rady; row++) {
tempave2[row] = tempave2[(row-1)] + (tempave[(row+rady)*W+col] - tempave[(row-rady-1)*W+col])/len;
dst[row*W+col] = dst[(row-1)*W+col] + (tempsqave[(row+rady)*W+col] - tempsqave[(row-rady-1)*W+col])/len;
}
for (int row=H-rady; row<H; row++) {
tempave2[row] = (tempave2[(row-1)]*len - tempave[(row-rady-1)*W+col])/(len-1);
dst[row*W+col] = (dst[(row-1)*W+col]*len - tempsqave[(row-rady-1)*W+col])/(len-1);
len --;
}
//now finish off
for (int row=0; row<H; row++) {
dst[row*W+col] = fabs(dst[row*W+col] - SQR(tempave2[row]));
//image_ave += src[row*W+col];
}
}
//image_ave /= (W*H);
delete buffer1;
delete buffer2;
delete buffer3;
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
template<typename T> 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<float>* buffer1 = new AlignedBuffer<float> (W*H);
float* temp = buffer1->data;
AlignedBuffer<float>* buffer2 = new AlignedBuffer<float> (W*H);
float* tempave = buffer2->data;
if (radx==0) {
for (int row=0; row<H; row++)
for (int col=0; col<H; col++) {
temp[row*H+col] = src[row*W+col];
}
} else {
//horizontal blur
for (int row = 0; row < H; row++) {
int len = radx + 1;
temp[row*W+0] = (float)src[row*W+0]/len;
for (int j=1; j<=radx; j++) {
temp[row*W+0] += (float)src[row*W+j]/len;
}
for (int col=1; col<=radx; col++) {
temp[row*W+col] = (temp[row*W+col-1]*len + src[row*W+col+radx])/(len+1);
len ++;
}
for (int col = radx+1; col < W-radx; col++) {
temp[row*W+col] = temp[row*W+col-1] + ((float)(src[row*W+col+radx] - src[row*W+col-radx-1]))/len;
}
for (int col=W-radx; col<W; col++) {
temp[row*W+col] = (temp[row*W+col-1]*len - src[row*W+col-radx-1])/(len-1);
len --;
}
}
}
if (rady==0) {
for (int row=0; row<H; row++)
for (int col=0; col<H; col++) {
tempave[row*W+col] = temp[row*W+col];
}
} else {
//vertical blur
for (int col = 0; col < W; col++) {
int len = rady + 1;
tempave[0*W+col] = temp[0*W+col]/len;
for (int i=1; i<=rady; i++) {
tempave[0*W+col] += temp[i*W+col]/len;
}
for (int row=1; row<=rady; row++) {
tempave[row*W+col] = (tempave[(row-1)*W+col]*len + temp[(row+rady)*W+col])/(len+1);
len ++;
}
for (int row = rady+1; row < H-rady; row++) {
tempave[row*W+col] = tempave[(row-1)*W+col] + (temp[(row+rady)*W+col] - temp[(row-rady-1)*W+col])/len;
}
for (int row=H-rady; row<H; row++) {
tempave[row*W+col] = (tempave[(row-1)*W+col]*len - temp[(row-rady-1)*W+col])/(len-1);
len --;
}
}
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//box blur absolute deviation
if (radx==0) {
for (int row=0; row<H; row++)
for (int col=0; col<H; col++) {
temp[row*H+col] = fabs(src[row*W+col]-tempave[row*W+col]);
}
} else {
//horizontal blur
for (int row = 0; row < H; row++) {
int len = radx + 1;
temp[row*W+0] = fabs(src[row*W+0]-tempave[row*W+0])/len;
for (int j=1; j<=radx; j++) {
temp[row*W+0] += fabs(src[row*W+j]-tempave[row*W+j])/len;
}
for (int col=1; col<=radx; col++) {
temp[row*W+col] = (temp[row*W+col-1]*len + fabs(src[row*W+col+radx]-tempave[row*W+col+radx]))/(len+1);
len ++;
}
for (int col = radx+1; col < W-radx; col++) {
temp[row*W+col] = temp[row*W+col-1] + (fabs(src[row*W+col+radx]-tempave[row*W+col+radx]) - \
fabs(src[row*W+col-radx-1]-tempave[row*W+col-radx-1]))/len;
}
for (int col=W-radx; col<W; col++) {
temp[row*W+col] = (temp[row*W+col-1]*len - fabs(src[row*W+col-radx-1]-tempave[row*W+col-radx-1]))/(len-1);
len --;
}
}
}
if (rady==0) {
for (int row=0; row<H; row++)
for (int col=0; col<H; col++) {
dst[row*W+col] = temp[row*W+col];
}
} else {
//vertical blur
for (int col = 0; col < W; col++) {
int len = rady + 1;
dst[0*W+col] = temp[0*W+col]/len;
for (int i=1; i<=rady; i++) {
dst[0*W+col] += temp[i*W+col]/len;
}
for (int row=1; row<=rady; row++) {
dst[row*W+col] = (dst[(row-1)*W+col]*len + temp[(row+rady)*W+col])/(len+1);
len ++;
}
for (int row = rady+1; row < H-rady; row++) {
dst[row*W+col] = dst[(row-1)*W+col] + (temp[(row+rady)*W+col] - temp[(row-rady-1)*W+col])/len;
}
for (int row=H-rady; row<H; row++) {
dst[row*W+col] = (dst[(row-1)*W+col]*len - temp[(row-rady-1)*W+col])/(len-1);
len --;
}
}
}
delete buffer1;
delete buffer2;
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
template<class T, class A> 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<float>* buffer = new AlignedBuffer<float> (W*H);
float* temp = buffer->data;
if (radx==0) {
for (int row=0; row<H; row++)
for (int col=0; col<H; col++) {
temp[row*H+col] = SQR(src[row*W+col]);
}
} else {
//horizontal blur
for (int row = 0; row < H; row++) {
int len = radx + 1;
temp[row*W+0] = SQR((float)src[row*W+0])/len;
for (int j=1; j<=radx; j++) {
temp[row*W+0] += SQR((float)src[row*W+j])/len;
}
for (int col=1; col<=radx; col++) {
temp[row*W+col] = (temp[row*W+col-1]*len + SQR(src[row*W+col+radx]))/(len+1);
len ++;
}
for (int col = radx+1; col < W-radx; col++) {
temp[row*W+col] = temp[row*W+col-1] + ((float)(SQR(src[row*W+col+radx]) - SQR(src[row*W+col-radx-1])))/len;
}
for (int col=W-radx; col<W; col++) {
temp[row*W+col] = (temp[row*W+col-1]*len - SQR(src[row*W+col-radx-1]))/(len-1);
len --;
}
}
}
if (rady==0) {
for (int row=0; row<H; row++)
for (int col=0; col<H; col++) {
dst[row*W+col] = temp[row*W+col];
}
} else {
//vertical blur
for (int col = 0; col < W; col++) {
int len = rady + 1;
dst[0*W+col] = temp[0*W+col]/len;
for (int i=1; i<=rady; i++) {
dst[0*W+col] += temp[i*W+col]/len;
}
for (int row=1; row<=rady; row++) {
dst[row*W+col] = (dst[(row-1)*W+col]*len + temp[(row+rady)*W+col])/(len+1);
len ++;
}
for (int row = rady+1; row < H-rady; row++) {
dst[row*W+col] = dst[(row-1)*W+col] + (temp[(row+rady)*W+col] - temp[(row-rady-1)*W+col])/len;
}
for (int row=H-rady; row<H; row++) {
dst[row*W+col] = (dst[(row-1)*W+col]*len - temp[(row-rady-1)*W+col])/(len-1);
len --;
}
}
}
delete buffer;
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#endif /* _BOXBLUR_H_ */

49
rtengine/cplx_wavelet_dec.h
View File

@@ -62,9 +62,8 @@ private:
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 first_lev_len, first_lev_offset;
float *first_lev_anal;
@@ -83,7 +82,7 @@ private:
public:
template<typename E>
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<typename E>
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<internal_type>(src, lvltot, padding, m_w, m_h, first_lev_anal+first_lev_len*2*n, \
dual_tree[0][2*n+m] = new wavelet_level<internal_type>(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<internal_type>(dual_tree[lvltot-1][2*n+m]->lopass()/*lopass*/, lvltot, 0/*no padding*/, \
dual_tree[lvltot][2*n+m] = new wavelet_level<internal_type>(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<typename E>
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<typename E>
void reconstruct(E * dst);
@@ -361,15 +367,12 @@ public:
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
template<typename E>
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<internal_type>(src, lvltot/*level*/, padding/*padding*/, m_w, m_h, \
wavelet_decomp[lvltot] = new wavelet_level<internal_type>(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<internal_type>(wavelet_decomp[lvltot-1]->lopass()/*lopass*/, lvltot/*level*/, 0/*no padding*/, \
wavelet_decomp[lvltot] = new wavelet_level<internal_type>(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);

474
rtengine/cplx_wavelet_level.h
View File

@@ -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<typename E>
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)
wavelet_level(E * src, int level, int subsamp, 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), subsamp_out((subsamp>>level)&1)
{
m_w2 = (w+2*skip*padding);
m_h2 = (h+2*skip*padding);
if (subsamp) {
skip = 1;
for (int n=0; n<level; n++) {
skip *= 2-((subsamp>>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<typename E>
void decompose_level(E *src, float *filterV, float *filterH, int len, int offset, int skip);
@@ -168,39 +195,14 @@ namespace rtengine {
}
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
template<typename T> template<typename E>
void wavelet_level<T>::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; i<srclen; i++, j += pitch)
{
int coset = i%skip;
int indx = i/skip;
tmp[coset*cosetlen + indx] = src[j];
}
//even up last row/column if srclen is not a multiple of 'skip'
for (size_t i=srclen; i<srclen+(srclen%skip); i++) {
tmp[i] = tmp[i-skip];
}
// extend each coset mirror-like by padding amount 'm_pad'
for (size_t coset=0; coset<skip*cosetlen; coset+=cosetlen) {
for (size_t i=1; i<=MIN(cosetlen-1,m_pad); i++) {
tmp[coset-i] = tmp[coset+i];
tmp[coset+cosetlen+i-1] = tmp[coset+cosetlen-i-1];
}
}*/
//create padded buffer from src data
for (size_t i = 0, j = 0; i<srclen; i++, j += pitch)
{
@@ -229,40 +231,11 @@ namespace rtengine {
/* 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 '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.
* aligning the 'offset' element of the filter with
* the input pixel, and skipping 'skip' pixels between taps
*
*/
//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; coset<srclen; coset+=cosetlen) {
for (size_t i = 0; i < (cosetlen); i++) {
float lo=0,hi=0;
if (i>taps && i<cosetlen-taps) {//bulk
for (int j=0, l=-offset; j<taps; j++, l++) {
lo += filterLo[j] * src[i-l];//lopass channel
hi += filterHi[j] * src[i-l];//hipass channel
}
} else {//boundary
for (int j=0; j<taps; j++) {
int arg = MAX(0,MIN(i+(offset-j),srclen-1));//clamped BC's
lo += filterLo[j] * src[arg];//lopass channel
hi += filterHi[j] * src[arg];//hipass channel
}
}
dstLo[(pitch*(coset+i))] = lo;
dstHi[(pitch*(coset+i))] = hi;
}
}*/
for (size_t i = 0; i < (srclen); i++) {
float lo=0,hi=0;
@@ -286,57 +259,20 @@ namespace rtengine {
}
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
template<typename T>
void wavelet_level<T>::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; i<srclen; i++, j+=pitch) {
int indx = (i%skip)*cosetlen + i/skip;
bufferLo[indx]=srcLo[j];
bufferHi[indx]=srcHi[j];
}
for (size_t coset=0; coset<srclen; coset+=cosetlen) {
for (size_t i = m_pad; i < (cosetlen-m_pad); i++) {
float tot=0;
if (i>taps && i<(cosetlen-taps)) {//bulk
for (int j=0, l=-shift; j<taps; j++, l++) {
tot += (filterLo[j] * bufferLo[i-l] + filterHi[j] * bufferHi[i-l]);
}
} else {//boundary
if (coset+i-m_pad == srclen) return;
for (int j=0, l=-shift; j<taps; j++, l++) {
int arg = MAX(0,MIN((i-l),srclen-1));//clamped BC's
tot += (filterLo[j] * bufferLo[arg] + filterHi[j] * bufferHi[arg]);
}
}
dst[pitch*(coset+i-m_pad)] = tot;
}
}*/
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// load into buffer
int srclen = (dstlen==m_w ? m_w2 : m_h2);//length of row/col in src (coarser level)
@@ -346,51 +282,114 @@ namespace rtengine {
bufferHi[i]=srcHi[j];
}
int shift=(taps-offset-1);
int shift=skip*(taps-offset-1);
for(size_t i = m_pad; i < (dstlen+m_pad); i++) {
float tot=0;
if (i>skip*taps && i<(srclen-skip*taps)) {//bulk
for (int j=0, l=-skip*shift; j<taps; j++, l+=skip) {
for (int j=0, l=-shift; j<taps; j++, l+=skip) {
tot += (filterLo[j] * bufferLo[i-l] + filterHi[j] * bufferHi[i-l]);
}
} else {//boundary
for (int j=0, l=-skip*shift; j<taps; j++, l+=skip) {
for (int j=0, l=-shift; j<taps; j++, l+=skip) {
int arg = MAX(0,MIN((i-l),srclen-1));//clamped BC's
tot += (filterLo[j] * bufferLo[arg] + filterHi[j] * bufferHi[arg]);
}
}
dst[pitch*(i-m_pad)] = tot;
if (tot<0.0f || tot>65535.0f) {
float xxx=tot;
float yyy=1.0f;
}
}
}
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
template<typename T>
void wavelet_level<T>::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<typename T>
void wavelet_level<T>::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<srclen; i++, j+=pitch) {
bufferLo[i]=srcLo[j];
bufferHi[i]=srcHi[j];
}
for(size_t i = m_pad; i < (dstlen+m_pad-skip); i++) {
dst[pitch*(i-m_pad)] = 0.5*(bufferLo[i] - bufferHi[i] + bufferLo[i+skip] + bufferHi[i+skip]);
}
for(size_t i = (dstlen+m_pad-skip); i < (dstlen+m_pad); i++) {
dst[pitch*(i-m_pad)] = 0.5*(bufferLo[i] - bufferHi[i] + bufferLo[i-skip] + bufferHi[i-skip]);
}
}
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
template<typename T> template<typename E>
void wavelet_level<T>::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<m_w; j++) {
loadbuffer(src+j, buffer, m_w/*pitch*/, m_h/*srclen*/);//pad a column of data and load it to buffer
AnalysisFilter (buffer, tmpLo+j, tmpHi+j, filterV, filterV+taps, taps, offset, m_w/*output_pitch*/, m_h/*srclen*/);
if (subsamp_out) {
AnalysisFilterSubsamp (buffer, tmpLo+j, tmpHi+j, filterV, filterV+taps, taps, offset, m_w/*output_pitch*/, m_h/*srclen*/);
} else {
AnalysisFilter (buffer, tmpLo+j, tmpHi+j, filterV, filterV+taps, taps, offset, m_w/*output_pitch*/, m_h/*srclen*/);
}
}
/* filter along rows */
for (int i=0; i<m_h2; i++) {
loadbuffer(tmpLo+i*m_w, buffer, 1/*pitch*/, m_w/*srclen*/);//pad a row of data and load it to buffer
AnalysisFilter (buffer, wavcoeffs[0]+i*m_w2, wavcoeffs[1]+i*m_w2, filterH, filterH+taps, taps, offset, 1/*output_pitch*/, m_w/*srclen*/);
if (subsamp_out) {
AnalysisFilterSubsamp (buffer, wavcoeffs[0]+i*m_w2, wavcoeffs[1]+i*m_w2, filterH, filterH+taps, taps, offset, 1/*output_pitch*/, m_w/*srclen*/);
} else {
AnalysisFilter (buffer, wavcoeffs[0]+i*m_w2, wavcoeffs[1]+i*m_w2, filterH, filterH+taps, taps, offset, 1/*output_pitch*/, m_w/*srclen*/);
}
loadbuffer(tmpHi+i*m_w, buffer, 1/*pitch*/, m_w/*srclen*/);
AnalysisFilter (buffer, wavcoeffs[2]+i*m_w2, wavcoeffs[3]+i*m_w2, filterH, filterH+taps, taps, offset, 1/*output_pitch*/, m_w/*srclen*/);
if (subsamp_out) {
AnalysisFilterSubsamp (buffer, wavcoeffs[2]+i*m_w2, wavcoeffs[3]+i*m_w2, filterH, filterH+taps, taps, offset, 1/*output_pitch*/, m_w/*srclen*/);
} else {
AnalysisFilter (buffer, wavcoeffs[2]+i*m_w2, wavcoeffs[3]+i*m_w2, filterH, filterH+taps, taps, offset, 1/*output_pitch*/, m_w/*srclen*/);
}
}
//imp_nr (wavcoeffs[0], m_w2, m_h2, 50.0f/20.0f);
@@ -400,8 +399,9 @@ namespace rtengine {
delete[] buffer;
}
/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
template<typename T> template<typename E>
void wavelet_level<T>::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<m_h2; i++) {
SynthesisFilter (wavcoeffs[0]+i*m_w2, wavcoeffs[1]+i*m_w2, tmpLo+i*m_w, bufferLo, bufferHi,
filterH, filterH+taps, taps, offset, 1/*pitch*/, m_w/*dstlen*/);
SynthesisFilter (wavcoeffs[2]+i*m_w2, wavcoeffs[3]+i*m_w2, tmpHi+i*m_w, bufferLo, bufferHi,
filterH, filterH+taps, taps, offset, 1/*pitch*/, m_w/*dstlen*/);
if (subsamp_out) {
SynthesisFilterSubsamp (wavcoeffs[0]+i*m_w2, wavcoeffs[1]+i*m_w2, tmpLo+i*m_w, bufferLo, bufferHi,
filterH, filterH+taps, taps, offset, 1/*pitch*/, m_w/*dstlen*/);
SynthesisFilterSubsamp (wavcoeffs[2]+i*m_w2, wavcoeffs[3]+i*m_w2, tmpHi+i*m_w, bufferLo, bufferHi,
filterH, filterH+taps, taps, offset, 1/*pitch*/, m_w/*dstlen*/);
} else {
SynthesisFilter (wavcoeffs[0]+i*m_w2, wavcoeffs[1]+i*m_w2, tmpLo+i*m_w, bufferLo, bufferHi,
filterH, filterH+taps, taps, offset, 1/*pitch*/, m_w/*dstlen*/);
SynthesisFilter (wavcoeffs[2]+i*m_w2, wavcoeffs[3]+i*m_w2, tmpHi+i*m_w, bufferLo, bufferHi,
filterH, filterH+taps, taps, offset, 1/*pitch*/, m_w/*dstlen*/);
}
}
/* filter along columns */
for (int j=0; j<m_w; j++) {
SynthesisFilter (tmpLo+j, tmpHi+j, dst+j, bufferLo, bufferHi,
filterV, filterV+taps, taps, offset, m_w/*pitch*/, m_h/*dstlen*/);
if (subsamp_out) {
SynthesisFilterSubsamp (tmpLo+j, tmpHi+j, dst+j, bufferLo, bufferHi,
filterV, filterV+taps, taps, offset, m_w/*pitch*/, m_h/*dstlen*/);
} else {
SynthesisFilter (tmpLo+j, tmpHi+j, dst+j, bufferLo, bufferHi,
filterV, filterV+taps, taps, offset, m_w/*pitch*/, m_h/*dstlen*/);
}
}
@@ -435,85 +447,153 @@ namespace rtengine {
}
/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
template<typename T>
void wavelet_level<T>::imp_nr (T* src, int width, int height, double thresh) {
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
template<typename T>
void wavelet_level<T>::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<srclen-skip*taps) {//bulk
for (int j=0, l=-skip*offset; j<taps; j++, l+=skip) {
lo += filterLo[j] * srcbuffer[i-l];//lopass channel
hi += filterHi[j] * srcbuffer[i-l];//hipass channel
}
} else {//boundary
for (int j=0; j<taps; j++) {
int arg = MAX(0,MIN(i+skip*(offset-j),srclen-1));//clamped BC's
lo += filterLo[j] * srcbuffer[arg];//lopass channel
hi += filterHi[j] * srcbuffer[arg];//hipass channel
}
}
dstLo[(pitch*(i/2))] = lo;
dstHi[(pitch*(i/2))] = hi;
}
// 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<double>* buffer = new AlignedBuffer<double> (MAX(width,height));
gaussHorizontal<float> (src, lpf, buffer, width, height, MAX(2.0,thresh-1.0), false);
gaussVertical<float> (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<typename T>
void wavelet_level<T>::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; i<srclen; i++, j+=pitch) {
bufferLo[i]=srcLo[j];
bufferHi[i]=srcHi[j];
}
int shift=skip*(taps-offset-1);//align filter with data
for(size_t i = m_pad; i < (dstlen+m_pad); i++) {
float tot=0;
int i_src = (i+shift)/2;
int begin = (i+shift)%2;
if (i>skip*taps && i<(srclen-skip*taps)) {//bulk
for (int j=begin, l=0; j<taps; j+=2, l+=skip) {
tot += 2*((filterLo[j] * bufferLo[i_src-l] + filterHi[j] * bufferHi[i_src-l]));
}
} else {//boundary
for (int j=begin, l=0; j<taps; j+=2, l+=skip) {
int arg = MAX(0,MIN((i_src-l),srclen-1));//clamped BC's
tot += 2*((filterLo[j] * bufferLo[arg] + filterHi[j] * bufferHi[arg]));
}
}
dst[pitch*(i-m_pad)] = tot;
}
}
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
template<typename T>
void wavelet_level<T>::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<typename T>
void wavelet_level<T>::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];
}
}
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
};

41
rtengine/improcfun.cc
View File

@@ -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; i<lab->W*lab->H; i++) {
int datalen = lab->W*lab->H;
for (int i=0; i<datalen; i++) {
lab->data[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; i<lab->W*lab->H; i++) {
//float recoeff = Ldecomp.level_coeffs(lvl,branch)[subband][i]/(2<<lvl);
//float imcoeff = Ldecomp.level_coeffs(lvl,branch+2)[subband][i]/(2<<lvl);
//float shrink = (SQR(recoeff)+SQR(imcoeff))/(SQR(recoeff)+SQR(imcoeff)+noisevar);
//lab->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<<lvl);
//float imcoeff = Ldecomp.level_coeffs(lvl,branch+2)[subband][i]/(2<<lvl);
//float shrink = (SQR(recoeff)+SQR(imcoeff))/(SQR(recoeff)+SQR(imcoeff)+noisevar);
//lab->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; i<lab->W*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; i<lab->W*lab->H; i++) {
for (int i=0; i<datalen; 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);
//impulse_nr (lab, 50.0f/20.0f);
}
}

10
rtengine/improcfun.h
View File

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