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

Wavelet levels Improvments - issue2702

Browse Source
This commit is contained in:
jdc
2015-03-31 08:42:41 +02:00
parent 4f2e23d16f
commit 5ee2bf9431
22 changed files with 2632 additions and 557 deletions
Download Patch File Download Diff File Expand all files Collapse all files

167
rtengine/FTblockDN.cc
View File

@@ -13,7 +13,7 @@
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// 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.
//
@@ -379,8 +379,8 @@ SSEFUNCTION void ImProcFunctions::RGB_denoise(int kall, Imagefloat * src, Imagef
//#endif
if (dnparams.luma==0 && dnparams.chroma==0 && !dnparams.median && !noiseLCurve && !noiseCCurve) {
//nothing to do; copy src to dst or do nothing in case src == dst
if(src != dst)
src->copyData(dst);
if(src != dst)
src->copyData(dst);
if(calclum) {
delete calclum;
calclum = NULL;
@@ -419,11 +419,11 @@ SSEFUNCTION void ImProcFunctions::RGB_denoise(int kall, Imagefloat * src, Imagef
// init luma noisevarL
const float noiseluma=(float) dnparams.luma;
const float noisevarL = (useNoiseLCurve && (denoiseMethodRgb || !isRAW)) ? (float) (SQR(((noiseluma+1.0)/125.0)*(10.+ (noiseluma+1.0)/25.0))) : (float) (SQR((noiseluma/125.0)*(1.0+ noiseluma/25.0)));
const bool denoiseLuminance = (noisevarL > 0.00001f);
const bool denoiseLuminance = (noisevarL > 0.00001f);
//printf("NL=%f \n",noisevarL);
if(useNoiseLCurve || useNoiseCCurve) {
int hei=calclum->height;
int wid=calclum->width;
int wid=calclum->width;
TMatrix wprofi = iccStore->workingSpaceMatrix (params->icm.working);
const float wpi[3][3] = {
@@ -479,7 +479,7 @@ SSEFUNCTION void ImProcFunctions::RGB_denoise(int kall, Imagefloat * src, Imagef
}
delete calclum;
calclum = NULL;
}
}
const short int imheight=src->height, imwidth=src->width;
@@ -550,7 +550,7 @@ SSEFUNCTION void ImProcFunctions::RGB_denoise(int kall, Imagefloat * src, Imagef
array2D<float> tilemask_in(TS,TS);
array2D<float> tilemask_out(TS,TS);
if(denoiseLuminance) {
const int border = MAX(2,TS/16);
const int border = MAX(2,TS/16);
for (int i=0; i<TS; i++) {
float i1 = abs((i>TS/2 ? i-TS+1 : i));
float vmask = (i1<border ? SQR(sin((M_PI*i1)/(2*border))) : 1.0f);
@@ -575,7 +575,7 @@ if(settings->leveldnti ==1) {
}
int numTries = 0;
if(ponder)
printf("Tiled denoise processing caused by Automatic Multizone mode\n");
printf("Tiled denoise processing caused by Automatic Multizone mode\n");
bool memoryAllocationFailed = false;
do {
@@ -593,16 +593,16 @@ do {
if(numtiles == 1)
dsttmp = dst;
else {
dsttmp = new Imagefloat(imwidth,imheight);
#ifdef _OPENMP
#pragma omp parallel for
#endif
for(int i=0;i<imheight;i++)
for(int j=0;j<imwidth;j++) {
dsttmp->r(i,j) = 0.f;
dsttmp->g(i,j) = 0.f;
dsttmp->b(i,j) = 0.f;
}
dsttmp = new Imagefloat(imwidth,imheight);
#ifdef _OPENMP
#pragma omp parallel for
#endif
for(int i=0;i<imheight;i++)
for(int j=0;j<imwidth;j++) {
dsttmp->r(i,j) = 0.f;
dsttmp->g(i,j) = 0.f;
dsttmp->b(i,j) = 0.f;
}
}
//now we have tile dimensions, overlaps
@@ -619,7 +619,7 @@ do {
// these are needed only for creation of the plans and will be freed before entering the parallel loop
fftwf_plan plan_forward_blox[2];
fftwf_plan plan_backward_blox[2];
if(denoiseLuminance) {
float *Lbloxtmp = (float*) fftwf_malloc(max_numblox_W*TS*TS*sizeof (float));
float *fLbloxtmp = (float*) fftwf_malloc(max_numblox_W*TS*TS*sizeof (float));
@@ -636,7 +636,7 @@ do {
plan_forward_blox[1] = fftwf_plan_many_r2r(2, nfwd, min_numblox_W, Lbloxtmp, NULL, 1, TS*TS, fLbloxtmp, NULL, 1, TS*TS, fwdkind, FFTW_MEASURE || FFTW_DESTROY_INPUT );
plan_backward_blox[1] = fftwf_plan_many_r2r(2, nfwd, min_numblox_W, fLbloxtmp, NULL, 1, TS*TS, Lbloxtmp, NULL, 1, TS*TS, bwdkind, FFTW_MEASURE || FFTW_DESTROY_INPUT );
fftwf_free ( Lbloxtmp );
fftwf_free ( fLbloxtmp );
fftwf_free ( fLbloxtmp );
}
#ifndef _OPENMP
@@ -712,7 +712,7 @@ do {
int tilebottom = MIN(imheight,tiletop+tileheight);
int width = tileright-tileleft;
int height = tilebottom-tiletop;
int width2 = (width+1)/2;
int width2 = (width+1)/2;
float realred, realblue;
float interm_med =(float) dnparams.chroma/10.0;
float intermred, intermblue;
@@ -725,7 +725,7 @@ do {
const float noisevarab_r = SQR(realred);
const float noisevarab_b = SQR(realblue);
//input L channel
//input L channel
array2D<float> *Lin;
//wavelet denoised image
LabImage * labdn = new LabImage(width,height);
@@ -886,18 +886,18 @@ do {
//now perform basic wavelet denoise
//arguments 4 and 5 of wavelet decomposition are max number of wavelet decomposition levels;
//and whether to subsample the image after wavelet filtering. Subsampling is coded as
//binary 1 or 0 for each level, eg subsampling = 0 means no subsampling, 1 means subsample
//binary 1 or 0 for each level, eg subsampling = 0 means no subsampling, 1 means subsample
//the first level only, 7 means subsample the first three levels, etc.
//actual implementation only works with subsampling set to 1
float interm_medT = (float) dnparams.chroma/10.0;
float interm_medT = (float) dnparams.chroma/10.0;
bool execwavelet = true;
if(!denoiseLuminance && interm_medT < 0.05f && dnparams.median && (dnparams.methodmed=="Lab" || dnparams.methodmed=="Lonly"))
execwavelet=false;//do not exec wavelet if sliders luminance and chroma are very small and median need
//we considered user don't want wavelet
//we considered user don't want wavelet
if(settings->leveldnautsimpl==1 && dnparams.Cmethod!="MAN")
execwavelet=true;
if(settings->leveldnautsimpl==0 && dnparams.C2method!="MANU")
execwavelet=true;
execwavelet=true;
if(execwavelet) {//gain time if user choose only median sliders L <=1 slider chrom master < 1
wavelet_decomposition* Ldecomp;
wavelet_decomposition* adecomp;
@@ -911,7 +911,7 @@ do {
else levwav=8;//maximum ==> I have increase Maxlevel in cplx_wavelet_dec.h from 8 to 9
if(nrQuality == QUALITY_HIGH)
levwav += settings->nrwavlevel;//increase level for enhanced mode
if(levwav>8) levwav=8;
if(levwav>8) levwav=8;
int minsizetile=min(tilewidth, tileheight);
int maxlev2=8;
if(minsizetile < 256) maxlev2 = 7;
@@ -925,7 +925,7 @@ do {
memoryAllocationFailed = true;
}
float madL[8][3];
if(!memoryAllocationFailed) {
if(!memoryAllocationFailed) {
// precalculate madL, because it's used in adecomp and bdecomp
int maxlvl = Ldecomp->maxlevel();
#ifdef _OPENMP
@@ -978,7 +978,7 @@ do {
adecomp->reconstruct(labdn->a[0]);
}
delete adecomp;
if(!memoryAllocationFailed) {
if(!memoryAllocationFailed) {
wavelet_decomposition* bdecomp = new wavelet_decomposition (labdn->b[0], labdn->W, labdn->H, levwav, 1, 1, max(1,denoiseNestedLevels));
if(bdecomp->memoryAllocationFailed) {
memoryAllocationFailed = true;
@@ -1007,30 +1007,31 @@ do {
bdecomp->reconstruct(labdn->b[0]);
}
delete bdecomp;
if(!memoryAllocationFailed) {
if(denoiseLuminance) {
if(!memoryAllocationFailed) {
if(denoiseLuminance) {
int edge=0;
if(nrQuality==QUALITY_STANDARD) {
if(!WaveletDenoiseAllL(*Ldecomp, noisevarlum, madL))//enhance mode
if(!WaveletDenoiseAllL(*Ldecomp, noisevarlum, madL, NULL, edge))//enhance mode
memoryAllocationFailed = true;
} else /*if(nrQuality==QUALITY_HIGH)*/ {
if(!WaveletDenoiseAll_BiShrinkL(*Ldecomp, noisevarlum, madL))//enhance mode
memoryAllocationFailed = true;
if(!memoryAllocationFailed)
if(!WaveletDenoiseAllL(*Ldecomp, noisevarlum, madL))
if(!WaveletDenoiseAllL(*Ldecomp, noisevarlum, madL, NULL, edge))
memoryAllocationFailed = true;
}
if(!memoryAllocationFailed) {
// copy labdn->L to Lin before it gets modified by reconstruction
Lin = new array2D<float>(width,height);
}
if(!memoryAllocationFailed) {
// copy labdn->L to Lin before it gets modified by reconstruction
Lin = new array2D<float>(width,height);
#ifdef _OPENMP
#pragma omp parallel for num_threads(denoiseNestedLevels) if(denoiseNestedLevels>1)
#endif
for(int i=0;i<height;i++)
for(int j=0;j<width;j++)
(*Lin)[i][j] = labdn->L[i][j];
Ldecomp->reconstruct(labdn->L[0]);
}
for(int i=0;i<height;i++)
for(int j=0;j<width;j++)
(*Lin)[i][j] = labdn->L[i][j];
Ldecomp->reconstruct(labdn->L[0]);
}
}
}
}
@@ -1122,7 +1123,7 @@ do {
// Main detail recovery algorithm: Block loop
//DCT block data storage
if(denoiseLuminance /*&& execwavelet*/) {
if(denoiseLuminance /*&& execwavelet*/) {
//residual between input and denoised L channel
array2D<float> Ldetail(width,height,ARRAY2D_CLEAR_DATA);
//pixel weight
@@ -1236,7 +1237,7 @@ do {
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#ifdef _OPENMP
#pragma omp parallel for num_threads(denoiseNestedLevels) if(denoiseNestedLevels>1)
@@ -1246,7 +1247,7 @@ do {
//may want to include masking threshold for large hipass data to preserve edges/detail
labdn->L[i][j] += Ldetail[i][j]/totwt[i][j]; //note that labdn initially stores the denoised hipass data
}
}
}
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// transform denoised "Lab" to output RGB
@@ -1278,10 +1279,10 @@ do {
}
//convert back to RGB and write to destination array
if (isRAW) {
if(!denoiseMethodRgb) {//Lab mode
if(!denoiseMethodRgb) {//Lab mode
realred /= 100.f;
realblue /= 100.f;
#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic,16) num_threads(denoiseNestedLevels) if(denoiseNestedLevels>1)
#endif
@@ -1297,7 +1298,7 @@ do {
float c_h=SQR(a)+SQR(b);
if(c_h>9000000.f){
a *= 1.f + qhighFactor*realred;
b *= 1.f + qhighFactor*realblue;
b *= 1.f + qhighFactor*realblue;
}
//convert XYZ
float X,Y,Z;
@@ -1326,7 +1327,7 @@ do {
dsttmp->b(i,j) += factor*b_;
}
}
}
}
} else {//RGB mode
#ifdef _OPENMP
#pragma omp parallel for num_threads(denoiseNestedLevels) if(denoiseNestedLevels>1)
@@ -1407,8 +1408,8 @@ do {
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
}
delete labdn;
if(denoiseLuminance)
delete labdn;
if(denoiseLuminance)
delete Lin;
}//end of tile row
@@ -1419,43 +1420,43 @@ do {
}
}
if(denoiseLuminance) {
if(denoiseLuminance) {
for(int i=0;i<denoiseNestedLevels*numthreads;i++) {
fftwf_free(LbloxArray[i]);
fftwf_free(fLbloxArray[i]);
}
}
}
#ifdef _OPENMP
omp_set_nested(oldNested);
#endif
//copy denoised image to output
if(numtiles>1) {
if(!memoryAllocationFailed)
dsttmp->copyData(dst);
if(!memoryAllocationFailed)
dsttmp->copyData(dst);
else if(dst != src)
src->copyData(dst);
src->copyData(dst);
delete dsttmp;
}
if (!isRAW && !memoryAllocationFailed) {//restore original image gamma
#ifdef _OPENMP
#pragma omp parallel for
#endif
for(int i=0;i<dst->height;i++)
for(int j=0;j<dst->width;j++) {
dst->r(i,j) = Color::gammatab_srgb[ dst->r(i,j) ];
dst->g(i,j) = Color::gammatab_srgb[ dst->g(i,j) ];
dst->b(i,j) = Color::gammatab_srgb[ dst->b(i,j) ];
}
}
if (!isRAW && !memoryAllocationFailed) {//restore original image gamma
#ifdef _OPENMP
#pragma omp parallel for
#endif
for(int i=0;i<dst->height;i++)
for(int j=0;j<dst->width;j++) {
dst->r(i,j) = Color::gammatab_srgb[ dst->r(i,j) ];
dst->g(i,j) = Color::gammatab_srgb[ dst->g(i,j) ];
dst->b(i,j) = Color::gammatab_srgb[ dst->b(i,j) ];
}
}
if(denoiseLuminance) {
// destroy the plans
fftwf_destroy_plan( plan_forward_blox[0] );
fftwf_destroy_plan( plan_backward_blox[0] );
fftwf_destroy_plan( plan_forward_blox[1] );
fftwf_destroy_plan( plan_backward_blox[1] );
fftwf_cleanup();
fftwf_cleanup();
}
} while(memoryAllocationFailed && numTries < 2 && (options.rgbDenoiseThreadLimit == 0) && !ponder);
if(memoryAllocationFailed)
@@ -1765,7 +1766,7 @@ SSEFUNCTION void ImProcFunctions::RGBtile_denoise (float * fLblox, int hblproc,
int histo[256] ALIGNED64 = {0};
//calculate histogram of absolute values of wavelet coeffs
for (int i=0; i<datalen; i++) {
for (int i=0; i<datalen; i++) {
histo[min(255,abs((int)DataList[i]))]++;
}
@@ -1877,7 +1878,9 @@ SSEFUNCTION bool ImProcFunctions::WaveletDenoiseAll_BiShrinkL(wavelet_decomposit
float ** WavCoeffs_L = WaveletCoeffs_L.level_coeffs(lvl);
if (lvl==maxlvl-1) {
ShrinkAllL(WaveletCoeffs_L, buffer, lvl, dir, noisevarlum, madL[lvl] );
float vari[3];
int edge=0;
ShrinkAllL(WaveletCoeffs_L, buffer, lvl, dir, noisevarlum, madL[lvl], NULL, edge );
} else {
//simple wavelet shrinkage
float * sfave = buffer[0]+32;
@@ -2069,11 +2072,12 @@ SSEFUNCTION bool ImProcFunctions::WaveletDenoiseAll_BiShrinkAB(wavelet_decomposi
}
bool ImProcFunctions::WaveletDenoiseAllL(wavelet_decomposition &WaveletCoeffs_L, float *noisevarlum, float madL[8][3])//mod JD
bool ImProcFunctions::WaveletDenoiseAllL(wavelet_decomposition &WaveletCoeffs_L, float *noisevarlum, float madL[8][3], float * vari, int edge)//mod JD
{
int maxlvl = WaveletCoeffs_L.maxlevel();
if(edge==1) maxlvl=3;//for refine denoise edge wavelet
int maxWL = 0, maxHL = 0;
for (int lvl=0; lvl<maxlvl; lvl++) {
if(WaveletCoeffs_L.level_W(lvl) > maxWL)
@@ -2081,7 +2085,7 @@ SSEFUNCTION bool ImProcFunctions::WaveletDenoiseAll_BiShrinkAB(wavelet_decomposi
if(WaveletCoeffs_L.level_H(lvl) > maxHL)
maxHL = WaveletCoeffs_L.level_H(lvl);
}
bool memoryAllocationFailed = false;
bool memoryAllocationFailed = false;
#ifdef _OPENMP
#pragma omp parallel num_threads(denoiseNestedLevels) if(denoiseNestedLevels>1)
#endif
@@ -2100,7 +2104,7 @@ SSEFUNCTION bool ImProcFunctions::WaveletDenoiseAll_BiShrinkAB(wavelet_decomposi
#endif
for (int lvl=0; lvl<maxlvl; lvl++) {
for (int dir=1; dir<4; dir++) {
ShrinkAllL(WaveletCoeffs_L, buffer, lvl, dir, noisevarlum, madL[lvl]);
ShrinkAllL(WaveletCoeffs_L, buffer, lvl, dir, noisevarlum, madL[lvl], vari, edge);
}
}
}
@@ -2124,7 +2128,7 @@ SSEFUNCTION bool ImProcFunctions::WaveletDenoiseAll_BiShrinkAB(wavelet_decomposi
maxWL = WaveletCoeffs_L.level_W(lvl);
if(WaveletCoeffs_L.level_H(lvl) > maxHL)
maxHL = WaveletCoeffs_L.level_H(lvl);
}
}
bool memoryAllocationFailed = false;
#ifdef _OPENMP
#pragma omp parallel num_threads(denoiseNestedLevels) if(denoiseNestedLevels>1)
@@ -2160,7 +2164,7 @@ SSEFUNCTION bool ImProcFunctions::WaveletDenoiseAll_BiShrinkAB(wavelet_decomposi
SSEFUNCTION void ImProcFunctions::ShrinkAllL(wavelet_decomposition &WaveletCoeffs_L, float **buffer, int level, int dir,
float *noisevarlum, float * madL )
float *noisevarlum, float * madL, float * vari, int edge )
{
//simple wavelet shrinkage
@@ -2176,7 +2180,12 @@ SSEFUNCTION void ImProcFunctions::ShrinkAllL(wavelet_decomposition &WaveletCoeff
float ** WavCoeffs_L = WaveletCoeffs_L.level_coeffs(level);
float mad_L = madL[dir-1] ;
if(edge==1) {
noisevarlum = blurBuffer; // we need one buffer, but fortunately we don't have to allocate a new one because we can use blurBuffer
for (int i=0; i<W_L*H_L; i++) {
noisevarlum[i]=vari[level];
}
}
float levelFactor = mad_L * 5.f / (float)(level+1);
#ifdef __SSE2__
__m128 magv;