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

Additional Median filter for Denoise issue2423

Browse Source
This commit is contained in:
jdc
2014-06-25 13:36:08 +02:00
parent 3d0a164941
commit 7636502a7c
10 changed files with 528 additions and 129 deletions
Download Patch File Download Diff File Expand all files Collapse all files

521
rtengine/FTblockDN.cc
View File

@@ -40,7 +40,7 @@
#include "boxblur.h"
#include "rt_math.h"
#include "mytime.h"
#include "sleef.c"
#include "sleef.c"
#include "opthelper.h"
#ifdef _OPENMP
@@ -62,6 +62,65 @@
#define epsilon 0.001f/(TS*TS) //tolerance
#define PIX_SORT(a,b) { if ((a)>(b)) {temp=(a);(a)=(b);(b)=temp;} }
#define med3(a0,a1,a2,a3,a4,a5,a6,a7,a8,median) { \
pp[0]=a0; pp[1]=a1; pp[2]=a2; pp[3]=a3; pp[4]=a4; pp[5]=a5; pp[6]=a6; pp[7]=a7; pp[8]=a8; \
PIX_SORT(pp[1],pp[2]); PIX_SORT(pp[4],pp[5]); PIX_SORT(pp[7],pp[8]); \
PIX_SORT(pp[0],pp[1]); PIX_SORT(pp[3],pp[4]); PIX_SORT(pp[6],pp[7]); \
PIX_SORT(pp[1],pp[2]); PIX_SORT(pp[4],pp[5]); PIX_SORT(pp[7],pp[8]); \
PIX_SORT(pp[0],pp[3]); PIX_SORT(pp[5],pp[8]); PIX_SORT(pp[4],pp[7]); \
PIX_SORT(pp[3],pp[6]); PIX_SORT(pp[1],pp[4]); PIX_SORT(pp[2],pp[5]); \
PIX_SORT(pp[4],pp[7]); PIX_SORT(pp[4],pp[2]); PIX_SORT(pp[6],pp[4]); \
PIX_SORT(pp[4],pp[2]); median=pp[4];} //pp4 = median
#define med2(a0,a1,a2,a3,a4,median) { \
pp[0]=a0; pp[1]=a1; pp[2]=a2; pp[3]=a3; pp[4]=a4; \
PIX_SORT(pp[0],pp[1]) ; PIX_SORT(pp[3],pp[4]) ; PIX_SORT(pp[0],pp[3]) ;\
PIX_SORT(pp[1],pp[4]) ; PIX_SORT(pp[1],pp[2]) ; PIX_SORT(pp[2],pp[3]) ;\
PIX_SORT(pp[1],pp[2]) ; median=pp[2] ;}
#define med5(a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20,a21,a22,a23,a24,median) { \
pp[0]=a0; pp[1]=a1; pp[2]=a2; pp[3]=a3; pp[4]=a4; pp[5]=a5; pp[6]=a6; pp[7]=a7; pp[8]=a8; pp[9]=a9; pp[10]=a10; pp[11]=a11; pp[12]=a12; \
pp[13]=a13; pp[14]=a14; pp[15]=a15; pp[16]=a16; pp[17]=a17; pp[18]=a18; pp[19]=a19; pp[20]=a20; pp[21]=a21; pp[22]=a22; pp[23]=a23; pp[24]=a24; \
PIX_SORT(pp[0], pp[1]) ; PIX_SORT(pp[3], pp[4]) ; PIX_SORT(pp[2], pp[4]) ;\
PIX_SORT(pp[2], pp[3]) ; PIX_SORT(pp[6], pp[7]) ; PIX_SORT(pp[5], pp[7]) ;\
PIX_SORT(pp[5], pp[6]) ; PIX_SORT(pp[9], pp[10]) ; PIX_SORT(pp[8], pp[10]) ;\
PIX_SORT(pp[8], pp[9]) ; PIX_SORT(pp[12], pp[13]) ; PIX_SORT(pp[11], pp[13]) ;\
PIX_SORT(pp[11], pp[12]) ; PIX_SORT(pp[15], pp[16]) ; PIX_SORT(pp[14], pp[16]) ;\
PIX_SORT(pp[14], pp[15]) ; PIX_SORT(pp[18], pp[19]) ; PIX_SORT(pp[17], pp[19]) ;\
PIX_SORT(pp[17], pp[18]) ; PIX_SORT(pp[21], pp[22]) ; PIX_SORT(pp[20], pp[22]) ;\
PIX_SORT(pp[20], pp[21]) ; PIX_SORT(pp[23], pp[24]) ; PIX_SORT(pp[2], pp[5]) ;\
PIX_SORT(pp[3], pp[6]) ; PIX_SORT(pp[0], pp[6]) ; PIX_SORT(pp[0], pp[3]) ;\
PIX_SORT(pp[4], pp[7]) ; PIX_SORT(pp[1], pp[7]) ; PIX_SORT(pp[1], pp[4]) ;\
PIX_SORT(pp[11], pp[14]) ; PIX_SORT(pp[8], pp[14]) ; PIX_SORT(pp[8], pp[11]) ;\
PIX_SORT(pp[12], pp[15]) ; PIX_SORT(pp[9], pp[15]) ; PIX_SORT(pp[9], pp[12]) ;\
PIX_SORT(pp[13], pp[16]) ; PIX_SORT(pp[10], pp[16]) ; PIX_SORT(pp[10], pp[13]) ;\
PIX_SORT(pp[20], pp[23]) ; PIX_SORT(pp[17], pp[23]) ; PIX_SORT(pp[17], pp[20]) ;\
PIX_SORT(pp[21], pp[24]) ; PIX_SORT(pp[18], pp[24]) ; PIX_SORT(pp[18], pp[21]) ;\
PIX_SORT(pp[19], pp[22]) ; PIX_SORT(pp[8], pp[17]) ; PIX_SORT(pp[9], pp[18]) ;\
PIX_SORT(pp[0], pp[18]) ; PIX_SORT(pp[0], pp[9]) ; PIX_SORT(pp[10], pp[19]) ;\
PIX_SORT(pp[1], pp[19]) ; PIX_SORT(pp[1], pp[10]) ; PIX_SORT(pp[11], pp[20]) ;\
PIX_SORT(pp[2], pp[20]) ; PIX_SORT(pp[2], pp[11]) ; PIX_SORT(pp[12], pp[21]) ;\
PIX_SORT(pp[3], pp[21]) ; PIX_SORT(pp[3], pp[12]) ; PIX_SORT(pp[13], pp[22]) ;\
PIX_SORT(pp[4], pp[22]) ; PIX_SORT(pp[4], pp[13]) ; PIX_SORT(pp[14], pp[23]) ;\
PIX_SORT(pp[5], pp[23]) ; PIX_SORT(pp[5], pp[14]) ; PIX_SORT(pp[15], pp[24]) ;\
PIX_SORT(pp[6], pp[24]) ; PIX_SORT(pp[6], pp[15]) ; PIX_SORT(pp[7], pp[16]) ;\
PIX_SORT(pp[7], pp[19]) ; PIX_SORT(pp[13], pp[21]) ; PIX_SORT(pp[15], pp[23]) ;\
PIX_SORT(pp[7], pp[13]) ; PIX_SORT(pp[7], pp[15]) ; PIX_SORT(pp[1], pp[9]) ;\
PIX_SORT(pp[3], pp[11]) ; PIX_SORT(pp[5], pp[17]) ; PIX_SORT(pp[11], pp[17]) ;\
PIX_SORT(pp[9], pp[17]) ; PIX_SORT(pp[4], pp[10]) ; PIX_SORT(pp[6], pp[12]) ;\
PIX_SORT(pp[7], pp[14]) ; PIX_SORT(pp[4], pp[6]) ; PIX_SORT(pp[4], pp[7]) ;\
PIX_SORT(pp[12], pp[14]) ; PIX_SORT(pp[10], pp[14]) ; PIX_SORT(pp[6], pp[7]) ;\
PIX_SORT(pp[10], pp[12]) ; PIX_SORT(pp[6], pp[10]) ; PIX_SORT(pp[6], pp[17]) ;\
PIX_SORT(pp[12], pp[17]) ; PIX_SORT(pp[7], pp[17]) ; PIX_SORT(pp[7], pp[10]) ;\
PIX_SORT(pp[12], pp[18]) ; PIX_SORT(pp[7], pp[12]) ; PIX_SORT(pp[10], pp[18]) ;\
PIX_SORT(pp[12], pp[20]) ; PIX_SORT(pp[10], pp[20]) ; PIX_SORT(pp[10], pp[12]) ;\
median=pp[12];}
#define ELEM_FLOAT_SWAP(a,b) { register float t=(a);(a)=(b);(b)=t; }
namespace rtengine {
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -83,14 +142,58 @@ namespace rtengine {
extern const Settings* settings;
// Median calculation using quicksort
float fq_sort2(float arr[], int n)
{
int low, high ;
int median;
int middle, ll, hh;
low = 0 ; high = n-1 ; median = (low + high) / 2;
for (;;) {
if (high <= low)
return arr[median] ;
if (high == low + 1) {
if (arr[low] > arr[high])
ELEM_FLOAT_SWAP(arr[low], arr[high]) ;
return arr[median] ;
}
middle = (low + high) / 2;
if (arr[middle] > arr[high]) ELEM_FLOAT_SWAP(arr[middle], arr[high]) ;
if (arr[low] > arr[high]) ELEM_FLOAT_SWAP(arr[low], arr[high]) ;
if (arr[middle] > arr[low]) ELEM_FLOAT_SWAP(arr[middle], arr[low]) ;
ELEM_FLOAT_SWAP(arr[middle], arr[low+1]) ;
ll = low + 1;
hh = high;
for (;;) {
do ll++; while (arr[low] > arr[ll]) ;
do hh--; while (arr[hh] > arr[low]) ;
if (hh < ll)
break;
ELEM_FLOAT_SWAP(arr[ll], arr[hh]) ;
}
ELEM_FLOAT_SWAP(arr[low], arr[hh]) ;
if (hh <= median)
low = ll;
if (hh >= median)
high = hh - 1;
}
}
void ImProcFunctions::RGB_denoise(Imagefloat * src, Imagefloat * dst, bool isRAW, const procparams::DirPyrDenoiseParams & dnparams, const procparams::DefringeParams & defringe, const double expcomp)
{
{
//#ifdef _DEBUG
// MyTime t1e,t2e;
// t1e.set();
MyTime t1e,t2e;
t1e.set();
//#endif
static MyMutex FftwMutex;
@@ -108,13 +211,16 @@ namespace rtengine {
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
const short int imheight=src->height, imwidth=src->width;
if (dnparams.luma==0 && dnparams.chroma==0) {
//nothing to do; copy src to dst
memcpy(dst->data,src->data,dst->width*dst->height*3*sizeof(float));
const short int imheight=src->height, imwidth=src->width;
if (dnparams.luma==0 && dnparams.chroma==0 && !dnparams.median) {
//nothing to do; copy src to dst or do nothing in case src == dst
if(src != dst)
memcpy(dst->data,src->data,dst->width*dst->height*3*sizeof(float));
return;
}
if (dnparams.luma!=0 || dnparams.chroma!=0) {
perf=false;
if(dnparams.dmethod=="RGB") perf=true;//RGB mode
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -165,6 +271,7 @@ namespace rtengine {
float incr=1.f;
float noisevar_Ldetail = SQR((float)(SQR(100.-dnparams.Ldetail) + 50.*(100.-dnparams.Ldetail)) * TS * 0.5f * incr);
bool enhance_denoise = dnparams.enhance;
// bool median_denoise = dnparams.median;
int gamlab = settings->denoiselabgamma;//gamma lab essentialy for Luminance detail
if(gamlab > 2) gamlab=2;
if(settings->verbose) printf("Denoise Lab=%i\n",gamlab);
@@ -173,7 +280,7 @@ namespace rtengine {
array2D<float> tilemask_out(TS,TS);
const int border = MAX(2,TS/16);
#ifdef _OPENMP
#pragma omp parallel for
#endif
@@ -294,9 +401,9 @@ namespace rtengine {
{
Lblox = (float*) fftwf_malloc(max_numblox_W*TS*TS*sizeof(float));
fLblox = (float*) fftwf_malloc(max_numblox_W*TS*TS*sizeof(float));
}
float * nbrwt = new float[TS*TS];
}
float * nbrwt = new float[TS*TS];
float * blurbuffer = new float[TS*TS];
#ifdef _OPENMP
#pragma omp for schedule(dynamic) collapse(2)
@@ -321,7 +428,7 @@ namespace rtengine {
//TODO: implement using AlignedBufferMP
//fill tile from image; convert RGB to "luma/chroma"
if (isRAW) {//image is raw; use channel differences for chroma channels
if(!perf){//lab mode
if(!perf){//lab mode
//modification Jacques feb 2013
for (int i=tiletop/*, i1=0*/; i<tilebottom; i++/*, i1++*/) {
int i1 = i - tiletop;
@@ -768,18 +875,188 @@ namespace rtengine {
}
delete dsttmp;
// 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();
}
//median 3x3 in complement
if(dnparams.median) {
int wid=dst->width, hei=dst->height;
float** tm;
tm = new float*[hei];
for (int i=0; i<hei; i++)
tm[i] = new float[wid];
Imagefloat *source;
if (dnparams.luma==0 && dnparams.chroma==0)
source = dst;
else
source = src;
int methmed=0;
int border = 1;
if(dnparams.medmethod=="soft")
methmed=0;
else if(dnparams.medmethod=="33")
methmed=1;
else if(dnparams.medmethod=="55") {
methmed = 3;
border = 2;
}
else if(dnparams.medmethod=="55soft") {
methmed = 2;
border = 2;
}
#pragma omp parallel
{
if(methmed < 2) {
#pragma omp for
for (int i=1; i<hei-1; i++) {
float pp[9],temp;
if(methmed == 0)
for (int j=1; j<wid-1; j++) {
med2(source->r(i,j),source->r(i-1,j),source->r(i+1,j),source->r(i,j+1),source->r(i,j-1),tm[i][j]);//3x3 soft
}
else
for (int j=1; j<wid-1; j++) {
med3(source->r(i,j),source->r(i-1,j),source->r(i+1,j),source->r(i,j+1),source->r(i,j-1),source->r(i-1,j-1),source->r(i-1,j+1),source->r(i+1,j-1),source->r(i+1,j+1),tm[i][j]);//3x3
}
}
} else {
#pragma omp for
for (int i=2; i<hei-2; i++) {
float pp[25],temp;
if(methmed == 3)
for (int j=2; j<wid-2; j++) {
med5(source->r(i,j),source->r(i-1,j),source->r(i+1,j),source->r(i,j+1),source->r(i,j-1),source->r(i-1,j-1),source->r(i-1,j+1),source->r(i+1,j-1),source->r(i+1,j+1),
source->r(i-2,j),source->r(i+2,j),source->r(i,j+2),source->r(i,j-2),source->r(i-2,j-2),source->r(i-2,j+2),source->r(i+2,j-2),source->r(i+2,j+2),
source->r(i-2,j+1),source->r(i+2,j+1),source->r(i-1,j+2),source->r(i-1,j-2),source->r(i-2,j-1),source->r(i+2,j-1),source->r(i+1,j+2),source->r(i+1,j-2),
tm[i][j]);//5x5
}
else
for (int j=2; j<wid-2; j++) {
pp[0]=source->r(i,j);pp[1]=source->r(i-1,j); pp[2]=source->r(i+1,j);pp[3]=source->r(i,j+1);pp[4]=source->r(i,j-1);pp[5]=source->r(i-1,j-1);pp[6]=source->r(i-1,j+1);
pp[7]=source->r(i+1,j-1);pp[8]=source->r(i+1,j+1);pp[9]=source->r(i+2,j);pp[10]=source->r(i-2,j);pp[11]=source->r(i,j+2);pp[12]=source->r(i,j-2);
fq_sort2(pp,13);
tm[i][j]=pp[6];//5x5 soft
}
}
}
#ifdef _OPENMP
#pragma omp for nowait
#endif
for(int i = border; i < hei-border; i++ ) {
for(int j = border; j < wid-border; j++) {
dst->r(i,j) = tm[i][j];
}
}
if(methmed < 2) {
#pragma omp for
for (int i=1; i<hei-1; i++) {
float pp[9],temp;
if(methmed == 0)
for (int j=1; j<wid-1; j++) {
med2(source->b(i,j),source->b(i-1,j),source->b(i+1,j),source->b(i,j+1),source->b(i,j-1),tm[i][j]);
}
else
for (int j=1; j<wid-1; j++) {
med3(source->b(i,j),source->b(i-1,j),source->b(i+1,j),source->b(i,j+1),source->b(i,j-1),source->b(i-1,j-1),source->b(i-1,j+1),source->b(i+1,j-1),source->b(i+1,j+1),tm[i][j]);
}
}
} else {
#pragma omp for
for (int i=2; i<hei-2; i++) {
float pp[25],temp;
if(methmed == 3)
for (int j=2; j<wid-2; j++) {
med5(source->b(i,j),source->b(i-1,j),source->b(i+1,j),source->b(i,j+1),source->b(i,j-1),source->b(i-1,j-1),source->b(i-1,j+1),source->b(i+1,j-1),source->b(i+1,j+1),
source->b(i-2,j),source->b(i+2,j),source->b(i,j+2),source->b(i,j-2),source->b(i-2,j-2),source->b(i-2,j+2),source->b(i+2,j-2),source->b(i+2,j+2),
source->b(i-2,j+1),source->b(i+2,j+1),source->b(i-1,j+2),source->b(i-1,j-2),source->b(i-2,j-1),source->b(i+2,j-1),source->b(i+1,j+2),source->b(i+1,j-2),
tm[i][j]);//5x5
}
else
for (int j=2; j<wid-2; j++) {
pp[0]=source->b(i,j);pp[1]=source->b(i-1,j); pp[2]=source->b(i+1,j);pp[3]=source->b(i,j+1);pp[4]=source->b(i,j-1);pp[5]=source->b(i-1,j-1);pp[6]=source->b(i-1,j+1);
pp[7]=source->b(i+1,j-1);pp[8]=source->b(i+1,j+1);pp[9]=source->b(i+2,j);pp[10]=source->b(i-2,j);pp[11]=source->b(i,j+2);pp[12]=source->b(i,j-2);
fq_sort2(pp,13);
tm[i][j]=pp[6];//5x5 soft
}
}
}
#ifdef _OPENMP
#pragma omp for nowait
#endif
for(int i = border; i < hei-border; i++ ) {
for(int j = border; j < wid-border; j++) {
dst->b(i,j) = tm[i][j];
}
}
if(methmed < 2) {
#pragma omp for
for (int i=1; i<hei-1; i++) {
float pp[9],temp;
if(methmed == 0)
for (int j=1; j<wid-1; j++) {
med2(source->g(i,j),source->g(i-1,j),source->g(i+1,j),source->g(i,j+1),source->g(i,j-1),tm[i][j]);
}
else
for (int j=1; j<wid-1; j++) {
med3(source->g(i,j),source->g(i-1,j),source->g(i+1,j),source->g(i,j+1),source->g(i,j-1),source->g(i-1,j-1),source->g(i-1,j+1),source->g(i+1,j-1),source->g(i+1,j+1),tm[i][j]);
}
}
} else {
#pragma omp for
for (int i=2; i<hei-2; i++) {
float pp[25],temp;
if(methmed == 3)
for (int j=2; j<wid-2; j++) {
med5(source->g(i,j),source->g(i-1,j),source->g(i+1,j),source->g(i,j+1),source->g(i,j-1),source->g(i-1,j-1),source->g(i-1,j+1),source->g(i+1,j-1),source->g(i+1,j+1),
source->g(i-2,j),source->g(i+2,j),source->g(i,j+2),source->g(i,j-2),source->g(i-2,j-2),source->g(i-2,j+2),source->g(i+2,j-2),source->g(i+2,j+2),
source->g(i-2,j+1),source->g(i+2,j+1),source->g(i-1,j+2),source->g(i-1,j-2),source->g(i-2,j-1),source->g(i+2,j-1),source->g(i+1,j+2),source->g(i+1,j-2),
tm[i][j]);//5x5
}
else
for (int j=2; j<wid-2; j++) {
pp[0]=source->g(i,j);pp[1]=source->g(i-1,j); pp[2]=source->g(i+1,j);pp[3]=source->g(i,j+1);pp[4]=source->g(i,j-1);pp[5]=source->g(i-1,j-1);pp[6]=source->g(i-1,j+1);
pp[7]=source->g(i+1,j-1);pp[8]=source->g(i+1,j+1);pp[9]=source->g(i+2,j);pp[10]=source->g(i-2,j);pp[11]=source->g(i,j+2);pp[12]=source->g(i,j-2);
fq_sort2(pp,13);
tm[i][j]=pp[6];//5x5 soft
}
}
}
#ifdef _OPENMP
#pragma omp for
#endif
for(int i = border; i < hei-border; i++ ) {
for(int j = border; j < wid-border; j++) {
dst->g(i,j) = tm[i][j];
}
}
}
for (int i=0; i<hei; i++)
delete [] tm[i];
delete [] tm;
}
//end median
//#ifdef _DEBUG
// if (settings->verbose) {
// t2e.set();
// printf("Denoise performed in %d usec:\n", t2e.etime(t1e));
// }
if (settings->verbose) {
t2e.set();
printf("Denoise performed in %d usec:\n", t2e.etime(t1e));
}
//#endif
}//end of main RGB_denoise
@@ -890,13 +1167,13 @@ SSEFUNCTION void ImProcFunctions::RGBtile_denoise (float * fLblox, int hblproc,
return (( (median-1) + (datalen/2-count_)/((float)(count-count_)) )/0.6745);
}
float ImProcFunctions::Mad( float * RESTRICT DataList, int datalen, int * RESTRICT histo) {
//computes Median Absolute Deviation
//DataList values should mostly have abs val < 65535
for (int i=0; i<65536; i++)
for (int i=0; i<65536; i++)
histo[i]=0;
//calculate histogram of absolute values of HH wavelet coeffs
@@ -916,7 +1193,7 @@ SSEFUNCTION void ImProcFunctions::RGBtile_denoise (float * fLblox, int hblproc,
// interpolate
return (( (median-1) + (datalen/2-count_)/((float)(count-count_)) )/0.6745);
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -951,7 +1228,7 @@ SSEFUNCTION void ImProcFunctions::WaveletDenoiseAll_BiShrink(wavelet_decompositi
madL[lvl][dir-1] = SQR(Mad(WavCoeffs_L[dir], Wlvl_L*Hlvl_L, madHisto));
mada[lvl][dir-1] = SQR(Mad(WavCoeffs_a[dir], Wlvl_ab*Hlvl_ab, madHisto));
madb[lvl][dir-1] = SQR(Mad(WavCoeffs_b[dir], Wlvl_ab*Hlvl_ab, madHisto));
}
}
}
delete [] madHisto;
@@ -982,8 +1259,8 @@ SSEFUNCTION void ImProcFunctions::WaveletDenoiseAll_BiShrink(wavelet_decompositi
//float ** WavPars_b = WaveletCoeffs_b.level_coeffs(lvl+1);
//simple wavelet shrinkage
float * sfave = new float[Wlvl_L*Hlvl_L];
float * WavCoeffsLtemp = new float[Hlvl_ab*Wlvl_ab];
float * sfave = new float[Wlvl_L*Hlvl_L];
float * WavCoeffsLtemp = new float[Hlvl_ab*Wlvl_ab];
// array2D<float> edge(Wlvl_L,Hlvl_L);
@@ -992,8 +1269,8 @@ SSEFUNCTION void ImProcFunctions::WaveletDenoiseAll_BiShrink(wavelet_decompositi
for (int dir=1; dir<4; dir++) {
float mad_L = madL[lvl][dir-1];
float mad_a = noisevar_abr*mada[lvl][dir-1];
float mad_b = noisevar_abb*madb[lvl][dir-1];
float mad_b = noisevar_abb*madb[lvl][dir-1];
//float mad_Lpar = madL[lvl+1][dir-1];
//float mad_apar = mada[lvl+1][dir-1];
@@ -1003,37 +1280,37 @@ SSEFUNCTION void ImProcFunctions::WaveletDenoiseAll_BiShrink(wavelet_decompositi
// float skip_L_ratio = WaveletCoeffs_L.level_stride(lvl+1)/skip_L;
if (noisevar_abr>0.01f || noisevar_abb>0.01f) {
for(int i=0;i<Hlvl_ab;i++)
for(int j=0;j<Wlvl_ab;j++)
for(int i=0;i<Hlvl_ab;i++)
for(int j=0;j<Wlvl_ab;j++)
WavCoeffsLtemp[i*Wlvl_ab+j] = WavCoeffs_L[dir][((i*skip_L)/skip_ab)*Wlvl_L + ((j*skip_L)/skip_ab)];
#ifdef __SSE2__
int j;
__m128 onev = _mm_set1_ps(1.f);
__m128 mad_Lm9v = _mm_set1_ps(mad_L * 9.f);
__m128 mad_av = _mm_set1_ps(mad_a);
__m128 mad_bv = _mm_set1_ps(mad_b);
__m128 epsv = _mm_set1_ps(eps);
__m128 mag_Lv, mag_av, mag_bv;
__m128 tempav, tempbv;
for (int i=0; i<Hlvl_ab; i++) {
int coeffloc_ab = i*Wlvl_ab;
for (j=0; j<Wlvl_ab-3; j+=4, coeffloc_ab+=4) {
tempav = LVFU(WavCoeffs_a[dir][coeffloc_ab]);
#ifdef __SSE2__
int j;
__m128 onev = _mm_set1_ps(1.f);
__m128 mad_Lm9v = _mm_set1_ps(mad_L * 9.f);
__m128 mad_av = _mm_set1_ps(mad_a);
__m128 mad_bv = _mm_set1_ps(mad_b);
__m128 epsv = _mm_set1_ps(eps);
__m128 mag_Lv, mag_av, mag_bv;
__m128 tempav, tempbv;
for (int i=0; i<Hlvl_ab; i++) {
int coeffloc_ab = i*Wlvl_ab;
for (j=0; j<Wlvl_ab-3; j+=4, coeffloc_ab+=4) {
tempav = LVFU(WavCoeffs_a[dir][coeffloc_ab]);
tempbv = LVFU(WavCoeffs_b[dir][coeffloc_ab]);
mag_Lv = LVFU(WavCoeffsLtemp[coeffloc_ab]);
mag_Lv = LVFU(WavCoeffsLtemp[coeffloc_ab]);
mag_av = SQRV(tempav)+epsv;
mag_bv = SQRV(tempbv)+epsv;
mag_Lv = SQRV(mag_Lv) + epsv;
mag_Lv = SQRV(mag_Lv) + epsv;
_mm_storeu_ps(&WavCoeffs_a[dir][coeffloc_ab], tempav * SQRV((onev-xexpf(-(mag_av/mad_av)-(mag_Lv/mad_Lm9v)))));
_mm_storeu_ps(&WavCoeffs_b[dir][coeffloc_ab], tempbv * SQRV((onev-xexpf(-(mag_bv/mad_bv)-(mag_Lv/mad_Lm9v)))));
}
}
for (; j<Wlvl_ab; j++,coeffloc_ab++) {
float mag_L = SQR(WavCoeffsLtemp[coeffloc_ab ])+eps;
float mag_a = SQR(WavCoeffs_a[dir][coeffloc_ab])+eps;
float mag_b = SQR(WavCoeffs_b[dir][coeffloc_ab])+eps;
WavCoeffs_a[dir][coeffloc_ab] *= SQR(1.f-xexpf(-(mag_a/mad_a)-(mag_L/(9.f*mad_L)))/*satfactor_a*/);
WavCoeffs_b[dir][coeffloc_ab] *= SQR(1.f-xexpf(-(mag_b/mad_b)-(mag_L/(9.f*mad_L)))/*satfactor_b*/);
}
}
}//now chrominance coefficients are denoised
#else
//printf(" dir=%d mad_L=%f mad_a=%f mad_b=%f \n",dir,sqrt(mad_L),sqrt(mad_a),sqrt(mad_b));
@@ -1066,24 +1343,24 @@ SSEFUNCTION void ImProcFunctions::WaveletDenoiseAll_BiShrink(wavelet_decompositi
WavCoeffs_b[dir][coeffloc_ab] *= SQR(1.f-xexpf(-(mag_b/mad_b)-(mag_L/(9.f*mad_L)))/*satfactor_b*/);
}
}//now chrominance coefficients are denoised
}//now chrominance coefficients are denoised
#endif
}
if (noisevar_L>0.01f) {
mad_L *= noisevar_L*5.f/(lvl+1);
#ifdef __SSE2__
int j;
__m128 mad_Lv = _mm_set1_ps(mad_L);
__m128 mad_Lm9v = _mm_set1_ps(mad_L * 9.f);
__m128 epsv = _mm_set1_ps(eps);
__m128 mag_Lv;
for (int i=0; i<Hlvl_L; i++) {
if (noisevar_L>0.01f) {
mad_L *= noisevar_L*5.f/(lvl+1);
#ifdef __SSE2__
int j;
__m128 mad_Lv = _mm_set1_ps(mad_L);
__m128 mad_Lm9v = _mm_set1_ps(mad_L * 9.f);
__m128 epsv = _mm_set1_ps(eps);
__m128 mag_Lv;
for (int i=0; i<Hlvl_L; i++) {
int coeffloc_L = i*Wlvl_L;
for (j=0; j<Wlvl_L-3; j+=4,coeffloc_L+=4) {
mag_Lv = SQRV(LVFU(WavCoeffs_L[dir][coeffloc_L]));
_mm_storeu_ps(&sfave[coeffloc_L], mag_Lv / ( mag_Lv + mad_Lv * xexpf(-mag_Lv/mad_Lm9v )+ epsv));
}
}
for (; j<Wlvl_L; j++, coeffloc_L++) {
float mag_L = SQR(WavCoeffs_L[dir][coeffloc_L]);
sfave[coeffloc_L] = mag_L/(mag_L+mad_L*xexpf(-mag_L/(9.f*mad_L))+eps);
@@ -1111,23 +1388,23 @@ SSEFUNCTION void ImProcFunctions::WaveletDenoiseAll_BiShrink(wavelet_decompositi
//gaussVertical<float> (edge, edge, buffer, Wlvl_L, Hlvl_L, 1<<(lvl+1), false);
boxblur(sfave, sfave, lvl+2, lvl+2, Wlvl_L, Hlvl_L);//increase smoothness by locally averaging shrinkage
#ifdef __SSE2__
__m128 tempLv;
__m128 tempL2v;
__m128 sf_Lv;
#ifdef __SSE2__
__m128 tempLv;
__m128 tempL2v;
__m128 sf_Lv;
for (int i=0; i<Hlvl_L; i++) {
int coeffloc_L = i*Wlvl_L;
for (j=0; j<Wlvl_L-3; j+=4,coeffloc_L+=4) {
tempLv = LVFU(WavCoeffs_L[dir][coeffloc_L]);
mag_Lv = SQRV(tempLv);
for (j=0; j<Wlvl_L-3; j+=4,coeffloc_L+=4) {
tempLv = LVFU(WavCoeffs_L[dir][coeffloc_L]);
mag_Lv = SQRV(tempLv);
tempL2v = LVFU(sfave[coeffloc_L]);
sf_Lv = mag_Lv/(mag_Lv + mad_Lv*xexpf(-mag_Lv/mad_Lm9v)+epsv);
sf_Lv = mag_Lv/(mag_Lv + mad_Lv*xexpf(-mag_Lv/mad_Lm9v)+epsv);
_mm_storeu_ps(&WavCoeffs_L[dir][coeffloc_L], tempLv * (SQRV(tempL2v)+SQRV(sf_Lv))/(tempL2v+sf_Lv+epsv));
//use smoothed shrinkage unless local shrinkage is much less
}//now luminance coeffs are denoised
}//now luminance coeffs are denoised
for (; j<Wlvl_L; j++,coeffloc_L++) {
float mag_L = SQR(WavCoeffs_L[dir][coeffloc_L]);
float sf_L = mag_L/(mag_L + mad_L*xexpf(-mag_L/(9.f*mad_L))+eps);
@@ -1137,8 +1414,8 @@ SSEFUNCTION void ImProcFunctions::WaveletDenoiseAll_BiShrink(wavelet_decompositi
}//now luminance coeffs are denoised
}
#else
#else
for (int i=0; i<Hlvl_L; i++)
for (int j=0; j<Wlvl_L; j++) {
@@ -1158,7 +1435,7 @@ SSEFUNCTION void ImProcFunctions::WaveletDenoiseAll_BiShrink(wavelet_decompositi
#endif
}
}
}
delete[] WavCoeffsLtemp;
delete[] sfave;
@@ -1193,7 +1470,7 @@ SSEFUNCTION void ImProcFunctions::WaveletDenoiseAll_BiShrink(wavelet_decompositi
float ** WavCoeffs_b = WaveletCoeffs_b.level_coeffs(lvl);
// printf("Hab : %d\n", Hlvl_ab);
// printf("Wab : %d\n", Wlvl_ab);
// printf("Wab : %d\n", Wlvl_ab);
ShrinkAll(WavCoeffs_L, WavCoeffs_a, WavCoeffs_b, lvl, Wlvl_L, Hlvl_L, Wlvl_ab, Hlvl_ab,
skip_L, skip_ab, noisevar_L, noisevar_abr, noisevar_abb, noi);
@@ -1204,7 +1481,7 @@ SSEFUNCTION void ImProcFunctions::WaveletDenoiseAll_BiShrink(wavelet_decompositi
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SSEFUNCTION 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_abr, float noisevar_abb, LabImage * noi,
int W_L, int H_L, int W_ab, int H_ab,int skip_L, int skip_ab, float noisevar_L, float noisevar_abr, float noisevar_abb, LabImage * noi,
float * madaa, float * madab, float * madaL, bool madCalculated)
{
@@ -1212,25 +1489,25 @@ SSEFUNCTION void ImProcFunctions::ShrinkAll(float ** WavCoeffs_L, float ** WavCo
const float eps = 0.01f;
float * sfave = new float[W_L*H_L];
float * sfavea = new float[W_L*H_L];
float * sfaveb = new float[W_L*H_L];
float * WavCoeffsLtemp = new float[H_ab*W_ab];
float * sfaveb = new float[W_L*H_L];
float * WavCoeffsLtemp = new float[H_ab*W_ab];
// int max;
//printf("\n level=%d \n",level);
for (int dir=1; dir<4; dir++) {
float mada, madb, madL;
if(madCalculated) {
mada = madaa[dir-1];
madb = madab[dir-1];
madL = madaL[dir-1] ;
for (int dir=1; dir<4; dir++) {
float mada, madb, madL;
if(madCalculated) {
mada = madaa[dir-1];
madb = madab[dir-1];
madL = madaL[dir-1] ;
} else {
int * madHisto = new int[65536];
mada = SQR(Mad(WavCoeffs_a[dir], W_ab*H_ab, madHisto));
madb = SQR(Mad(WavCoeffs_b[dir], W_ab*H_ab, madHisto));
madL = SQR(Mad(WavCoeffs_L[dir], W_L*H_L, madHisto));
delete [] madHisto;
delete [] madHisto;
}
// printf(" dir=%d mad_L=%f mad_a=%f mad_b=%f skip_ab=%i \n",dir,sqrt(madL),sqrt(mada),sqrt(madb), skip_ab);
@@ -1239,39 +1516,39 @@ SSEFUNCTION void ImProcFunctions::ShrinkAll(float ** WavCoeffs_L, float ** WavCo
float mad_a = mada*noisevar_abr; // noisevar_abr between 0..2.25=default 100=middle value ==> 582=max
float mad_b = madb*noisevar_abb;
if (noisevar_abr>0.01f || noisevar_abb>0.01f) {
for(int i=0;i<H_ab;i++)
for(int j=0;j<W_ab;j++)
if (noisevar_abr>0.01f || noisevar_abb>0.01f) {
for(int i=0;i<H_ab;i++)
for(int j=0;j<W_ab;j++)
WavCoeffsLtemp[i*W_ab+j] = WavCoeffs_L[dir][((i*skip_L)/skip_ab)*W_L + ((j*skip_L)/skip_ab)];
#ifdef __SSE2__
int j;
__m128 onev = _mm_set1_ps(1.f);
__m128 madLm9v = _mm_set1_ps(madL * 9.f);
__m128 mad_av = _mm_set1_ps(mad_a);
__m128 mad_bv = _mm_set1_ps(mad_b);
__m128 epsv = _mm_set1_ps(eps);
__m128 mag_Lv, mag_av, mag_bv;
for (int i=0; i<H_ab; i++) {
int coeffloc_ab = i*W_ab;
#ifdef __SSE2__
int j;
__m128 onev = _mm_set1_ps(1.f);
__m128 madLm9v = _mm_set1_ps(madL * 9.f);
__m128 mad_av = _mm_set1_ps(mad_a);
__m128 mad_bv = _mm_set1_ps(mad_b);
__m128 epsv = _mm_set1_ps(eps);
__m128 mag_Lv, mag_av, mag_bv;
for (int i=0; i<H_ab; i++) {
int coeffloc_ab = i*W_ab;
for (j=0; j<W_ab-3; j+=4, coeffloc_ab+=4) {
mag_Lv = LVFU(WavCoeffsLtemp[coeffloc_ab]);
mag_Lv = LVFU(WavCoeffsLtemp[coeffloc_ab]);
mag_av = SQRV(LVFU(WavCoeffs_a[dir][coeffloc_ab]))+epsv;
mag_bv = SQRV(LVFU(WavCoeffs_b[dir][coeffloc_ab]))+epsv;
mag_Lv = SQRV(mag_Lv) + epsv;
mag_Lv = SQRV(mag_Lv) + epsv;
_mm_storeu_ps(&sfavea[coeffloc_ab], (onev-xexpf(-(mag_av/mad_av)-(mag_Lv/madLm9v))));
_mm_storeu_ps(&sfaveb[coeffloc_ab], (onev-xexpf(-(mag_bv/mad_bv)-(mag_Lv/madLm9v))));
}
}
for (; j<W_ab; j++,coeffloc_ab++) {
float mag_L = SQR(WavCoeffsLtemp[coeffloc_ab])+eps;
float mag_a = SQR(WavCoeffs_a[dir][coeffloc_ab])+eps;
float mag_b = SQR(WavCoeffs_b[dir][coeffloc_ab])+eps;
sfavea[coeffloc_ab] = (1.f-xexpf(-(mag_a/mad_a)-(mag_L/(9.f*madL))));
sfaveb[coeffloc_ab] = (1.f-xexpf(-(mag_b/mad_b)-(mag_L/(9.f*madL))));
}
}
}//now chrominance coefficients are denoised
#else
for (int i=0; i<H_ab; i++) {
#else
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);
@@ -1283,32 +1560,32 @@ SSEFUNCTION void ImProcFunctions::ShrinkAll(float ** WavCoeffs_L, float ** WavCo
// 'firm' threshold of chroma coefficients
//WavCoeffs_a[dir][coeffloc_ab] *= (1-exp(-(mag_a/mad_a)-(mag_L/(9*madL))));//(coeff_a>2*thresh_a ? 1 : (coeff_a<thresh_a ? 0 : (coeff_a/thresh_a - 1)));
//WavCoeffs_b[dir][coeffloc_ab] *= (1-exp(-(mag_b/mad_b)-(mag_L/(9*madL))));//(coeff_b>2*thresh_b ? 1 : (coeff_b<thresh_b ? 0 : (coeff_b/thresh_b - 1)));
}
}
}//now chrominance coefficients are denoised
#endif
boxblur(sfavea, sfavea, level+2, level+2, W_ab, H_ab);//increase smoothness by locally averaging shrinkage
boxblur(sfaveb, sfaveb, level+2, level+2, W_ab, H_ab);//increase smoothness by locally averaging shrinkage
#ifdef __SSE2__
__m128 sfav, sfbv;
__m128 sfaveav, sfavebv;
for (int i=0; i<H_ab; i++) {
int coeffloc_ab = i*W_ab;
for (j=0; j<W_ab-3; j+=4,coeffloc_ab+=4) {
mag_Lv = LVFU(WavCoeffsLtemp[coeffloc_ab]);
#ifdef __SSE2__
__m128 sfav, sfbv;
__m128 sfaveav, sfavebv;
for (int i=0; i<H_ab; i++) {
int coeffloc_ab = i*W_ab;
for (j=0; j<W_ab-3; j+=4,coeffloc_ab+=4) {
mag_Lv = LVFU(WavCoeffsLtemp[coeffloc_ab]);
mag_av = SQRV(LVFU(WavCoeffs_a[dir][coeffloc_ab]))+epsv;
mag_bv = SQRV(LVFU(WavCoeffs_b[dir][coeffloc_ab]))+epsv;
mag_Lv = SQRV(mag_Lv) + epsv;
sfav = (onev - xexpf(-(mag_av/mad_av)-(mag_Lv/madLm9v)));
mag_Lv = SQRV(mag_Lv) + epsv;
sfav = (onev - xexpf(-(mag_av/mad_av)-(mag_Lv/madLm9v)));
sfbv = (onev - xexpf(-(mag_bv/mad_bv)-(mag_Lv/madLm9v)));
sfaveav = LVFU(sfavea[coeffloc_ab]);
sfavebv = LVFU(sfaveb[coeffloc_ab]);
sfaveav = LVFU(sfavea[coeffloc_ab]);
sfavebv = LVFU(sfaveb[coeffloc_ab]);
//use smoothed shrinkage unless local shrinkage is much less
_mm_storeu_ps( &WavCoeffs_a[dir][coeffloc_ab], LVFU(WavCoeffs_a[dir][coeffloc_ab]) * (SQRV(sfaveav)+SQRV(sfav))/(sfaveav+sfav+epsv));
_mm_storeu_ps( &WavCoeffs_a[dir][coeffloc_ab], LVFU(WavCoeffs_a[dir][coeffloc_ab]) * (SQRV(sfaveav)+SQRV(sfav))/(sfaveav+sfav+epsv));
_mm_storeu_ps( &WavCoeffs_b[dir][coeffloc_ab], LVFU(WavCoeffs_b[dir][coeffloc_ab]) * (SQRV(sfavebv)+SQRV(sfbv))/(sfavebv+sfbv+epsv));
}//now chrominance coefficients are denoised
}//now chrominance coefficients are denoised
for (; j<W_ab; j++,coeffloc_ab++) {
//modification Jacques feb 2013
float mag_L = SQR(WavCoeffsLtemp[coeffloc_ab])+eps;
@@ -1323,8 +1600,8 @@ SSEFUNCTION void ImProcFunctions::ShrinkAll(float ** WavCoeffs_L, float ** WavCo
WavCoeffs_a[dir][coeffloc_ab] *= (SQR(sfavea[coeffloc_ab])+SQR(sfa))/(sfavea[coeffloc_ab]+sfa+eps);
WavCoeffs_b[dir][coeffloc_ab] *= (SQR(sfaveb[coeffloc_ab])+SQR(sfb))/(sfaveb[coeffloc_ab]+sfb+eps);
}//now chrominance coefficients are denoised
}
}//now chrominance coefficients are denoised
}
#else
for (int i=0; i<H_ab; i++) {
for (int j=0; j<W_ab; j++) {
@@ -1344,8 +1621,8 @@ SSEFUNCTION void ImProcFunctions::ShrinkAll(float ** WavCoeffs_L, float ** WavCo
WavCoeffs_a[dir][coeffloc_ab] *= (SQR(sfavea[coeffloc_ab])+SQR(sfa))/(sfavea[coeffloc_ab]+sfa+eps);
WavCoeffs_b[dir][coeffloc_ab] *= (SQR(sfaveb[coeffloc_ab])+SQR(sfb))/(sfaveb[coeffloc_ab]+sfb+eps);
}//now chrominance coefficients are denoised
}
}//now chrominance coefficients are denoised
}
#endif
}
@@ -1403,10 +1680,10 @@ SSEFUNCTION void ImProcFunctions::ShrinkAll(float ** WavCoeffs_L, float ** WavCo
}
}
}
delete[] sfave;
delete[] sfavea;
delete[] sfaveb;
delete[] sfaveb;
delete[] WavCoeffsLtemp;
}