Start median rework (#3346)
- Added basic algorithms to `rtengine/median.h` - Converted first occurrences
This commit is contained in:
Flössie
@@ -22,7 +22,7 @@
|
||||
//
|
||||
////////////////////////////////////////////////////////////////
|
||||
|
||||
#include <math.h>
|
||||
#include <cmath>
|
||||
#include <fftw3.h>
|
||||
#include "../rtgui/threadutils.h"
|
||||
#include "rtengine.h"
|
||||
@@ -36,6 +36,7 @@
|
||||
#include "sleef.c"
|
||||
#include "opthelper.h"
|
||||
#include "cplx_wavelet_dec.h"
|
||||
#include "median.h"
|
||||
|
||||
#ifdef _OPENMP
|
||||
#include <omp.h>
|
||||
@@ -48,6 +49,8 @@
|
||||
|
||||
#define epsilon 0.001f/(TS*TS) //tolerance
|
||||
|
||||
#define PIX_SORT(a,b) { if ((a)>(b)) {temp=(a);(a)=(b);(b)=temp;} }
|
||||
|
||||
#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]) ;\
|
||||
@@ -292,7 +295,6 @@ void ImProcFunctions::Median_Denoise( float **src, float **dst, const int width,
|
||||
medianOut[i][j] = medianIn[i][j];
|
||||
}
|
||||
} else if(medianType == MED_3X3STRONG) {
|
||||
float pp[9], temp;
|
||||
int j;
|
||||
|
||||
for (j = 0; j < border; j++) {
|
||||
@@ -300,7 +302,7 @@ void ImProcFunctions::Median_Denoise( float **src, float **dst, const int width,
|
||||
}
|
||||
|
||||
for (; j < width - border; j++) {
|
||||
med3(medianIn[i][j] , medianIn[i - 1][j], medianIn[i + 1][j] , medianIn[i][j + 1], medianIn[i][j - 1], medianIn[i - 1][j - 1], medianIn[i - 1][j + 1], medianIn[i + 1][j - 1], medianIn[i + 1][j + 1], medianOut[i][j]);
|
||||
medianOut[i][j] = median(medianIn[i][j] , medianIn[i - 1][j], medianIn[i + 1][j] , medianIn[i][j + 1], medianIn[i][j - 1], medianIn[i - 1][j - 1], medianIn[i - 1][j + 1], medianIn[i + 1][j - 1], medianIn[i + 1][j + 1]);
|
||||
}
|
||||
|
||||
for(; j < width; j++) {
|
||||
@@ -1859,7 +1861,7 @@ SSEFUNCTION void ImProcFunctions::RGB_denoise(int kall, Imagefloat * src, Imagef
|
||||
}
|
||||
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
|
||||
tm[i][j] = median(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)); //3x3
|
||||
}
|
||||
}
|
||||
} else {
|
||||
@@ -1922,7 +1924,7 @@ SSEFUNCTION void ImProcFunctions::RGB_denoise(int kall, Imagefloat * src, Imagef
|
||||
}
|
||||
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]);
|
||||
tm[i][j] = median(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));
|
||||
}
|
||||
}
|
||||
} else {
|
||||
@@ -1987,7 +1989,7 @@ SSEFUNCTION void ImProcFunctions::RGB_denoise(int kall, Imagefloat * src, Imagef
|
||||
}
|
||||
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]);
|
||||
tm[i][j] = median(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));
|
||||
}
|
||||
}
|
||||
} else {
|
||||
|
||||
@@ -31,6 +31,7 @@
|
||||
#include "../rtgui/myflatcurve.h"
|
||||
#include "rt_math.h"
|
||||
#include "opthelper.h"
|
||||
#include "median.h"
|
||||
|
||||
#ifdef _OPENMP
|
||||
#include <omp.h>
|
||||
@@ -766,7 +767,6 @@ SSEFUNCTION void ImProcFunctions::Badpixelscam(CieImage * src, CieImage * dst, d
|
||||
#pragma omp parallel
|
||||
{
|
||||
int ip, in, jp, jn;
|
||||
float pp[9], temp;
|
||||
#pragma omp for nowait //nowait because next loop inside this parallel region is independent on this one
|
||||
|
||||
for (int i = 0; i < height; i++) {
|
||||
@@ -795,7 +795,7 @@ SSEFUNCTION void ImProcFunctions::Badpixelscam(CieImage * src, CieImage * dst, d
|
||||
jn = j + 2;
|
||||
}
|
||||
|
||||
med3(sraa[ip][jp], sraa[ip][j], sraa[ip][jn], sraa[i][jp], sraa[i][j], sraa[i][jn], sraa[in][jp], sraa[in][j], sraa[in][jn], tmaa[i][j]);
|
||||
tmaa[i][j] = median(sraa[ip][jp], sraa[ip][j], sraa[ip][jn], sraa[i][jp], sraa[i][j], sraa[i][jn], sraa[in][jp], sraa[in][j], sraa[in][jn]);
|
||||
}
|
||||
}
|
||||
|
||||
@@ -827,7 +827,7 @@ SSEFUNCTION void ImProcFunctions::Badpixelscam(CieImage * src, CieImage * dst, d
|
||||
jn = j + 2;
|
||||
}
|
||||
|
||||
med3(srbb[ip][jp], srbb[ip][j], srbb[ip][jn], srbb[i][jp], srbb[i][j], srbb[i][jn], srbb[in][jp], srbb[in][j], srbb[in][jn], tmbb[i][j]);
|
||||
tmbb[i][j] = median(srbb[ip][jp], srbb[ip][j], srbb[ip][jn], srbb[i][jp], srbb[i][j], srbb[i][jn], srbb[in][jp], srbb[in][j], srbb[in][jn]);
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -1374,7 +1374,6 @@ SSEFUNCTION void ImProcFunctions::BadpixelsLab(LabImage * src, LabImage * dst, d
|
||||
#pragma omp parallel
|
||||
{
|
||||
int ip, in, jp, jn;
|
||||
float pp[9], temp;
|
||||
#pragma omp for nowait //nowait because next loop inside this parallel region is independent on this one
|
||||
|
||||
for (int i = 0; i < height; i++) {
|
||||
@@ -1403,7 +1402,7 @@ SSEFUNCTION void ImProcFunctions::BadpixelsLab(LabImage * src, LabImage * dst, d
|
||||
jn = j + 2;
|
||||
}
|
||||
|
||||
med3(sraa[ip][jp], sraa[ip][j], sraa[ip][jn], sraa[i][jp], sraa[i][j], sraa[i][jn], sraa[in][jp], sraa[in][j], sraa[in][jn], tmaa[i][j]);
|
||||
tmaa[i][j] = median(sraa[ip][jp], sraa[ip][j], sraa[ip][jn], sraa[i][jp], sraa[i][j], sraa[i][jn], sraa[in][jp], sraa[in][j], sraa[in][jn]);
|
||||
}
|
||||
}
|
||||
|
||||
@@ -1435,7 +1434,7 @@ SSEFUNCTION void ImProcFunctions::BadpixelsLab(LabImage * src, LabImage * dst, d
|
||||
jn = j + 2;
|
||||
}
|
||||
|
||||
med3(srbb[ip][jp], srbb[ip][j], srbb[ip][jn], srbb[i][jp], srbb[i][j], srbb[i][jn], srbb[in][jp], srbb[in][j], srbb[in][jn], tmbb[i][j]);
|
||||
tmbb[i][j] = median(srbb[ip][jp], srbb[ip][j], srbb[ip][jn], srbb[i][jp], srbb[i][j], srbb[i][jn], srbb[in][jp], srbb[in][j], srbb[in][jn], tmbb[i][j]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@@ -1,751 +0,0 @@
|
||||
/*
|
||||
* This file is part of RawTherapee.
|
||||
*
|
||||
* 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/>.
|
||||
*
|
||||
* <20> 2010 Emil Martinec <ejmartin@uchicago.edu>
|
||||
*
|
||||
*/
|
||||
|
||||
#include <cstddef>
|
||||
#include <cmath>
|
||||
#include "curves.h"
|
||||
#include "labimage.h"
|
||||
#include "improcfun.h"
|
||||
#include "array2D.h"
|
||||
#include "rt_math.h"
|
||||
|
||||
#ifdef _OPENMP
|
||||
#include <omp.h>
|
||||
#endif
|
||||
|
||||
#define CLIPC(a) ((a)>-32000?((a)<32000?(a):32000):-32000)
|
||||
|
||||
#define DIRWT_L(i1,j1,i,j) ( rangefn_L[(data_fine->L[i1][j1]-data_fine->L[i][j]+32768)] )
|
||||
|
||||
#define DIRWT_AB(i1,j1,i,j) ( rangefn_ab[(data_fine->a[i1][j1]-data_fine->a[i][j]+32768)] * \
|
||||
rangefn_ab[(data_fine->L[i1][j1]-data_fine->L[i][j]+32768)] * \
|
||||
rangefn_ab[(data_fine->b[i1][j1]-data_fine->b[i][j]+32768)] )
|
||||
|
||||
//#define NRWT_L(a) (nrwt_l[a] )
|
||||
|
||||
#define NRWT_AB (nrwt_ab[(hipass[1]+32768)] * nrwt_ab[(hipass[2]+32768)])
|
||||
|
||||
|
||||
#define med3(a,b,c) (a<b ? (b<c ? b : (a<c ? c : a)) : (a<c ? a : (b<c ? c : b)))
|
||||
|
||||
#define hmf(a11,a12,a13,a21,a22,a23,a31,a32,a33) \
|
||||
(med3(a22,med3(a22,med3(a12,a22,a32),med3(a21,a22,a23)), \
|
||||
med3(a22,med3(a11,a22,a33),med3(a13,a22,a31))) )
|
||||
|
||||
#define PIX_SORT(a,b) { if ((a)>(b)) {temp=(a);(a)=(b);(b)=temp;} }
|
||||
|
||||
#define med3x3(a0,a1,a2,a3,a4,a5,a6,a7,a8,median) { \
|
||||
p[0]=a0; p[1]=a1; p[2]=a2; p[3]=a3; p[4]=a4; p[5]=a5; p[6]=a6; p[7]=a7; p[8]=a8; \
|
||||
PIX_SORT(p[1],p[2]); PIX_SORT(p[4],p[5]); PIX_SORT(p[7],p[8]); \
|
||||
PIX_SORT(p[0],p[1]); PIX_SORT(p[3],p[4]); PIX_SORT(p[6],p[7]); \
|
||||
PIX_SORT(p[1],p[2]); PIX_SORT(p[4],p[5]); PIX_SORT(p[7],p[8]); \
|
||||
PIX_SORT(p[0],p[3]); PIX_SORT(p[5],p[8]); PIX_SORT(p[4],p[7]); \
|
||||
PIX_SORT(p[3],p[6]); PIX_SORT(p[1],p[4]); PIX_SORT(p[2],p[5]); \
|
||||
PIX_SORT(p[4],p[7]); PIX_SORT(p[4],p[2]); PIX_SORT(p[6],p[4]); \
|
||||
PIX_SORT(p[4],p[2]); median=p[4];} //a4 is the median
|
||||
|
||||
|
||||
namespace rtengine
|
||||
{
|
||||
|
||||
static const int maxlevel = 4;
|
||||
|
||||
//sequence of scales
|
||||
//static const int scales[8] = {1,2,4,8,16,32,64,128};
|
||||
//sequence of pitches
|
||||
//static const int pitches[8] = {1,1,1,1,1,1,1,1};
|
||||
|
||||
//sequence of scales
|
||||
//static const int scales[8] = {1,1,1,1,1,1,1,1};
|
||||
//sequence of pitches
|
||||
//static const int pitches[8] = {2,2,2,2,2,2,2,2};
|
||||
|
||||
//sequence of scales
|
||||
//static const int scales[8] = {1,1,2,2,4,4,8,8};
|
||||
//sequence of pitches
|
||||
//static const int pitches[8] = {2,1,2,1,2,1,2,1};
|
||||
|
||||
//sequence of scales
|
||||
static const int scales[8] = {1, 1, 2, 4, 8, 16, 32, 64};
|
||||
//sequence of pitches
|
||||
static const int pitches[8] = {2, 1, 1, 1, 1, 1, 1, 1};
|
||||
|
||||
//pitch is spacing of subsampling
|
||||
//scale is spacing of directional averaging weights
|
||||
//example 1: no subsampling at any level -- pitch=1, scale=2^n
|
||||
//example 2: subsampling by 2 every level -- pitch=2, scale=1 at each level
|
||||
//example 3: no subsampling at first level, subsampling by 2 thereafter --
|
||||
// pitch =1, scale=1 at first level; pitch=2, scale=2 thereafter
|
||||
|
||||
|
||||
|
||||
|
||||
void ImProcFunctions :: dirpyrLab_denoise(LabImage * src, LabImage * dst, const procparams::DirPyrDenoiseParams & dnparams )
|
||||
{
|
||||
float gam = dnparams.gamma / 3.0;
|
||||
//float gam = 2.0;//min(3.0, 0.1*fabs(c[4])/3.0+0.001);
|
||||
float gamthresh = 0.03;
|
||||
float gamslope = exp(log((double)gamthresh) / gam) / gamthresh;
|
||||
|
||||
LUTf gamcurve(65536, 0);
|
||||
|
||||
//DiagonalCurve* lumacurve = new DiagonalCurve (dnparams.lumcurve, CURVES_MIN_POLY_POINTS);
|
||||
//DiagonalCurve* chromacurve = new DiagonalCurve (dnparams.chromcurve, CURVES_MIN_POLY_POINTS);
|
||||
//LUTf Lcurve(65536);
|
||||
//LUTf abcurve(65536);
|
||||
for (int i = 0; i < 65536; i++) {
|
||||
int g = (int)(CurveFactory::gamma((double)i / 65535.0, gam, gamthresh, gamslope, 1.0, 0.0) * 65535.0);
|
||||
gamcurve[i] = CLIP(g);
|
||||
/*float val = (float)i/65535.0;
|
||||
float Lval = (2*(lumacurve->getVal(val)));
|
||||
float abval = (2*(chromacurve->getVal(val)));
|
||||
|
||||
Lcurve[i] = SQR(Lval);
|
||||
abcurve[i] = SQR(abval);
|
||||
if (i % 1000 ==0) printf("%d Lmult=%f abmult=%f \n",i,Lcurve[i],abcurve[i]);*/
|
||||
}
|
||||
|
||||
//delete lumacurve;
|
||||
//delete chromacurve;
|
||||
|
||||
|
||||
|
||||
//#pragma omp parallel for if (multiThread)
|
||||
for (int i = 0; i < src->H; i++) {
|
||||
for (int j = 0; j < src->W; j++) {
|
||||
//src->L[i][j] = CurveFactory::flinterp(gamcurve,src->L[i][j]);
|
||||
src->L[i][j] = gamcurve[src->L[i][j]];
|
||||
}
|
||||
}
|
||||
|
||||
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
|
||||
|
||||
LUTf rangefn_L(65536);
|
||||
LUTf nrwt_l(65536);
|
||||
|
||||
LUTf rangefn_ab(65536);
|
||||
LUTf nrwt_ab(65536);
|
||||
|
||||
//set up NR weight functions
|
||||
|
||||
//gamma correction for chroma in shadows
|
||||
float nrwtl_norm = ((CurveFactory::gamma((double)65535.0 / 65535.0, gam, gamthresh, gamslope, 1.0, 0.0)) -
|
||||
(CurveFactory::gamma((double)75535.0 / 65535.0, gam, gamthresh, gamslope, 1.0, 0.0)));
|
||||
|
||||
for (int i = 0; i < 65536; i++) {
|
||||
nrwt_l[i] = ((CurveFactory::gamma((double)i / 65535.0, gam, gamthresh, gamslope, 1.0, 0.0) -
|
||||
CurveFactory::gamma((double)(i + 10000) / 65535.0, gam, gamthresh, gamslope, 1.0, 0.0)) ) / nrwtl_norm;
|
||||
//if (i % 100 ==0) printf("%d %f \n",i,nrwt_l[i]);
|
||||
}
|
||||
|
||||
float tonefactor = nrwt_l[32768];
|
||||
|
||||
float noise_L = 10.0 * dnparams.luma;
|
||||
float noisevar_L = SQR(noise_L);
|
||||
|
||||
float noise_ab = 100.0 * dnparams.chroma;
|
||||
float noisevar_ab = SQR(noise_ab);
|
||||
|
||||
|
||||
//set up range functions
|
||||
for (int i = 0; i < 65536; i++) {
|
||||
rangefn_L[i] = (( exp(-(double)fabs(i - 32768) * tonefactor / (1.0 + noise_L)) * (1.0 + noisevar_L) / ((double)(i - 32768) * (double)(i - 32768) + noisevar_L + 1.0)));
|
||||
}
|
||||
|
||||
for (int i = 0; i < 65536; i++) {
|
||||
rangefn_ab[i] = (( exp(-(double)fabs(i - 32768) * tonefactor / (1.0 + 3 * noise_ab)) * (1.0 + noisevar_ab) / ((double)(i - 32768) * (double)(i - 32768) + noisevar_ab + 1.0)));
|
||||
}
|
||||
|
||||
|
||||
for (int i = 0; i < 65536; i++) {
|
||||
nrwt_ab[i] = ((1.0 + abs(i - 32768) / (1.0 + 8 * noise_ab)) * exp(-(double)fabs(i - 32768) / (1.0 + 8 * noise_ab) ) );
|
||||
}
|
||||
|
||||
|
||||
//for (int i=0; i<65536; i+=100) printf("%d %d \n",i,gamcurve[i]);
|
||||
|
||||
|
||||
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
|
||||
|
||||
int level;
|
||||
|
||||
LabImage * dirpyrLablo[maxlevel];
|
||||
int w = (int)((src->W - 1) / pitches[0]) + 1;
|
||||
int h = (int)((src->H - 1) / pitches[0]) + 1;
|
||||
dirpyrLablo[0] = new LabImage(w, h);
|
||||
|
||||
for (level = 1; level < maxlevel; level++) {
|
||||
w = (int)((w - 1) / pitches[level]) + 1;
|
||||
h = (int)((h - 1) / pitches[level]) + 1;
|
||||
dirpyrLablo[level] = new LabImage(w, h);
|
||||
};
|
||||
|
||||
|
||||
//////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
// c[0] = luma = noise_L
|
||||
// c[1] = chroma = noise_ab
|
||||
// c[2] decrease of noise var with scale
|
||||
// c[3] radius of domain blur at each level
|
||||
// c[4] shadow smoothing
|
||||
// c[5] edge preservation
|
||||
|
||||
level = 0;
|
||||
|
||||
int scale = scales[level];
|
||||
|
||||
int pitch = pitches[level];
|
||||
|
||||
//int thresh = 10 * c[8];
|
||||
//impulse_nr (src, src, m_w1, m_h1, thresh, noisevar);
|
||||
|
||||
dirpyr(src, dirpyrLablo[0], 0, rangefn_L, rangefn_ab, pitch, scale, dnparams.luma, dnparams.chroma );
|
||||
|
||||
level = 1;
|
||||
|
||||
while(level < maxlevel) {
|
||||
scale = scales[level];
|
||||
pitch = pitches[level];
|
||||
|
||||
dirpyr(dirpyrLablo[level - 1], dirpyrLablo[level], level, rangefn_L, rangefn_ab, pitch, scale, dnparams.luma, dnparams.chroma );
|
||||
|
||||
level ++;
|
||||
}
|
||||
|
||||
|
||||
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
|
||||
|
||||
for(int level = maxlevel - 1; level > 0; level--) {
|
||||
|
||||
int scale = scales[level];
|
||||
int pitch = pitches[level];
|
||||
idirpyr(dirpyrLablo[level], dirpyrLablo[level - 1], level, rangefn_L, nrwt_l, nrwt_ab, pitch, scale, dnparams.luma, dnparams.chroma/*, Lcurve, abcurve*/ );
|
||||
}
|
||||
|
||||
|
||||
scale = scales[0];
|
||||
pitch = pitches[0];
|
||||
|
||||
// freeing as much memory as possible since the next call to idirpyr will need lots
|
||||
for(int i = 1; i < maxlevel; i++) {
|
||||
delete dirpyrLablo[i];
|
||||
}
|
||||
|
||||
idirpyr(dirpyrLablo[0], dst, 0, rangefn_L, nrwt_l, nrwt_ab, pitch, scale, dnparams.luma, dnparams.chroma/*, Lcurve, abcurve*/ );
|
||||
|
||||
// freeing the last bunch of memory
|
||||
delete dirpyrLablo[0];
|
||||
|
||||
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
|
||||
float igam = 1 / gam;
|
||||
float igamthresh = gamthresh * gamslope;
|
||||
float igamslope = 1 / gamslope;
|
||||
|
||||
for (int i = 0; i < 65536; i++) {
|
||||
gamcurve[i] = (CurveFactory::gamma((float)i / 65535.0, igam, igamthresh, igamslope, 1.0, 0.0) * 65535.0);
|
||||
}
|
||||
|
||||
|
||||
if (dnparams.luma > 0) {
|
||||
for (int i = 0; i < dst->H; i++)
|
||||
for (int j = 0; j < dst->W; j++) {
|
||||
dst->L[i][j] = gamcurve[dst->L[i][j]];
|
||||
}
|
||||
} else {
|
||||
for (int i = 0; i < dst->H; i++)
|
||||
for (int j = 0; j < dst->W; j++) {
|
||||
dst->L[i][j] = gamcurve[src->L[i][j]];
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
void ImProcFunctions::dirpyr(LabImage* data_fine, LabImage* data_coarse, int level,
|
||||
LUTf & rangefn_L, LUTf & rangefn_ab, int pitch, int scale,
|
||||
const int luma, const int chroma )
|
||||
{
|
||||
|
||||
//pitch is spacing of subsampling
|
||||
//scale is spacing of directional averaging weights
|
||||
//example 1: no subsampling at any level -- pitch=1, scale=2^n
|
||||
//example 2: subsampling by 2 every level -- pitch=2, scale=1 at each level
|
||||
//example 3: no subsampling at first level, subsampling by 2 thereafter --
|
||||
// pitch =1, scale=1 at first level; pitch=2, scale=2 thereafter
|
||||
|
||||
|
||||
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
// calculate weights, compute directionally weighted average
|
||||
|
||||
int width = data_fine->W;
|
||||
int height = data_fine->H;
|
||||
|
||||
//generate domain kernel
|
||||
int halfwin = 3;//min(ceil(2*sig),3);
|
||||
int scalewin = halfwin * scale;
|
||||
|
||||
#ifdef _OPENMP
|
||||
#pragma omp parallel for
|
||||
#endif
|
||||
|
||||
for(int i = 0; i < height; i += pitch ) {
|
||||
int i1 = i / pitch;
|
||||
|
||||
for(int j = 0, j1 = 0; j < width; j += pitch, j1++) {
|
||||
float dirwt_l, dirwt_ab, norm_l, norm_ab;
|
||||
//float lops,aops,bops;
|
||||
float Lout, aout, bout;
|
||||
norm_l = norm_ab = 0;//if we do want to include the input pixel in the sum
|
||||
Lout = 0;
|
||||
aout = 0;
|
||||
bout = 0;
|
||||
|
||||
for(int inbr = (i - scalewin); inbr <= (i + scalewin); inbr += scale) {
|
||||
if (inbr < 0 || inbr > height - 1) {
|
||||
continue;
|
||||
}
|
||||
|
||||
for (int jnbr = (j - scalewin); jnbr <= (j + scalewin); jnbr += scale) {
|
||||
if (jnbr < 0 || jnbr > width - 1) {
|
||||
continue;
|
||||
}
|
||||
|
||||
dirwt_l = DIRWT_L(inbr, jnbr, i, j);
|
||||
dirwt_ab = DIRWT_AB(inbr, jnbr, i, j);
|
||||
Lout += dirwt_l * data_fine->L[inbr][jnbr];
|
||||
aout += dirwt_ab * data_fine->a[inbr][jnbr];
|
||||
bout += dirwt_ab * data_fine->b[inbr][jnbr];
|
||||
norm_l += dirwt_l;
|
||||
norm_ab += dirwt_ab;
|
||||
}
|
||||
}
|
||||
|
||||
//lops = Lout/norm;//diagnostic
|
||||
//aops = aout/normab;//diagnostic
|
||||
//bops = bout/normab;//diagnostic
|
||||
|
||||
data_coarse->L[i1][j1] = Lout / norm_l; //low pass filter
|
||||
data_coarse->a[i1][j1] = aout / norm_ab;
|
||||
data_coarse->b[i1][j1] = bout / norm_ab;
|
||||
|
||||
|
||||
/*if (level<2 && i>0 && i<height-1 && j>0 && j<width-1) {
|
||||
Lhmf = hmf(data_fine->L[i-1][j-1], data_fine->L[i-1][j], data_fine->L[i-1][j+1], \
|
||||
data_fine->L[i][j-1], data_fine->L[i][j], data_fine->L[i][j+1], \
|
||||
data_fine->L[i+1][j-1], data_fine->L[i+1][j], data_fine->L[i+1][j+1]);
|
||||
//med3x3(data_fine->L[i-1][j-1], data_fine->L[i-1][j], data_fine->L[i-1][j+1], \
|
||||
data_fine->L[i][j-1], data_fine->L[i][j], data_fine->L[i][j+1], \
|
||||
data_fine->L[i+1][j-1], data_fine->L[i+1][j], data_fine->L[i+1][j+1],Lmed);
|
||||
|
||||
data_coarse->L[i1][j1] = Lhmf;
|
||||
}*/
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
}
|
||||
|
||||
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
void ImProcFunctions::idirpyr(LabImage* data_coarse, LabImage* data_fine, int level, LUTf &rangefn_L, LUTf & nrwt_l, LUTf & nrwt_ab,
|
||||
int pitch, int scale, const int luma, const int chroma/*, LUTf & Lcurve, LUTf & abcurve*/ )
|
||||
{
|
||||
|
||||
int width = data_fine->W;
|
||||
int height = data_fine->H;
|
||||
|
||||
array2D<float> nrfactorL (width, height);
|
||||
|
||||
//float eps = 0.0;
|
||||
|
||||
// c[0] noise_L
|
||||
// c[1] noise_ab (relative to noise_L)
|
||||
// c[2] decrease of noise var with scale
|
||||
// c[3] radius of domain blur at each level
|
||||
// c[4] shadow smoothing
|
||||
|
||||
float noisevar_L = 4 * SQR(25.0 * luma);
|
||||
float noisevar_ab = 2 * SQR(100.0 * chroma);
|
||||
float scalefactor = 1.0 / pow(2.0, (level + 1) * 2); //change the last 2 to 1 for longer tail of higher scale NR
|
||||
|
||||
noisevar_L *= scalefactor;
|
||||
|
||||
// for coarsest level, take non-subsampled lopass image and subtract from lopass_fine to generate hipass image
|
||||
|
||||
// denoise hipass image, add back into lopass_fine to generate denoised image at fine scale
|
||||
|
||||
// now iterate:
|
||||
// (1) take denoised image at level n, expand and smooth using gradient weights from lopass image at level n-1
|
||||
// the result is the smoothed image at level n-1
|
||||
// (2) subtract smoothed image at level n-1 from lopass image at level n-1 to make hipass image at level n-1
|
||||
// (3) denoise the hipass image at level n-1
|
||||
// (4) add the denoised image at level n-1 to the smoothed image at level n-1 to make the denoised image at level n-1
|
||||
|
||||
// note that the coarsest level amounts to skipping step (1) and doing (2,3,4).
|
||||
// in other words, skip step one if pitch=1
|
||||
|
||||
// step (1)
|
||||
|
||||
if (pitch == 1) {
|
||||
|
||||
// step (1-2-3-4)
|
||||
|
||||
#ifdef _OPENMP
|
||||
#pragma omp parallel
|
||||
#endif
|
||||
{
|
||||
|
||||
#ifdef _OPENMP
|
||||
#pragma omp for
|
||||
#endif
|
||||
|
||||
for(int i = 0; i < height; i++)
|
||||
for(int j = 0; j < width; j++)
|
||||
{
|
||||
float hipass[3], hpffluct[3], tonefactor, nrfactor;
|
||||
|
||||
tonefactor = (nrwt_l[data_coarse->L[i][j]]);
|
||||
|
||||
hipass[1] = data_fine->a[i][j] - data_coarse->a[i][j];
|
||||
hipass[2] = data_fine->b[i][j] - data_coarse->b[i][j];
|
||||
|
||||
//Wiener filter
|
||||
//luma
|
||||
if (level < 2) {
|
||||
hipass[0] = data_fine->L[i][j] - data_coarse->L[i][j];
|
||||
hpffluct[0] = SQR(hipass[0]) + SQR(hipass[1]) + SQR(hipass[2]) + 0.001;
|
||||
nrfactorL[i][j] = (1.0 + hpffluct[0]) / (1.0 + hpffluct[0] + noisevar_L /* * Lcurve[data_coarse->L[i][j]]*/);
|
||||
//hipass[0] *= hpffluct[0]/(hpffluct[0]+noisevar_L);
|
||||
//data_fine->L[i][j] = CLIP(hipass[0]+data_coarse->L[i][j]);
|
||||
}
|
||||
|
||||
//chroma
|
||||
//hipass[1] = data_fine->a[i][j]-data_coarse->a[i][j];
|
||||
//hipass[2] = data_fine->b[i][j]-data_coarse->b[i][j];
|
||||
hpffluct[1] = SQR(hipass[1] * tonefactor) + 0.001;
|
||||
hpffluct[2] = SQR(hipass[2] * tonefactor) + 0.001;
|
||||
nrfactor = (hpffluct[1] + hpffluct[2]) / ((hpffluct[1] + hpffluct[2]) + noisevar_ab * NRWT_AB);
|
||||
|
||||
hipass[1] *= nrfactor;
|
||||
hipass[2] *= nrfactor;
|
||||
|
||||
data_fine->a[i][j] = hipass[1] + data_coarse->a[i][j];
|
||||
data_fine->b[i][j] = hipass[2] + data_coarse->b[i][j];
|
||||
}
|
||||
|
||||
if (level < 2)
|
||||
{
|
||||
#ifdef _OPENMP
|
||||
#pragma omp for
|
||||
#endif
|
||||
|
||||
for(int i = 0; i < height; i++)
|
||||
for(int j = 0; j < width; j++) {
|
||||
|
||||
float dirwt_l, norm_l;
|
||||
float nrfctrave = 0;
|
||||
norm_l = 0;//if we do want to include the input pixel in the sum
|
||||
|
||||
for(int inbr = max(0, i - 1); inbr <= min(height - 1, i + 1); inbr++) {
|
||||
for (int jnbr = max(0, j - 1); jnbr <= min(width - 1, j + 1); jnbr++) {
|
||||
dirwt_l = DIRWT_L(inbr, jnbr, i, j);
|
||||
nrfctrave += dirwt_l * nrfactorL[inbr][jnbr];
|
||||
norm_l += dirwt_l;
|
||||
}
|
||||
}
|
||||
|
||||
nrfctrave /= norm_l;
|
||||
//nrfctrave = nrfactorL[i][j];
|
||||
//nrfctrave=1;
|
||||
|
||||
float hipass[3];
|
||||
|
||||
//luma
|
||||
|
||||
/*if (i>0 && i<height-1 && j>0 && j<width-1) {
|
||||
med3x3(nrfactorL[i-1][j-1], nrfactorL[i-1][j], nrfactorL[i-1][j+1], \
|
||||
nrfactorL[i][j-1], nrfactorL[i][j], nrfactorL[i][j+1], \
|
||||
nrfactorL[i+1][j-1], nrfactorL[i+1][j], nrfactorL[i+1][j+1], median);
|
||||
//median = hmf(nrfactorL[i-1][j-1], nrfactorL[i-1][j], nrfactorL[i-1][j+1], \
|
||||
nrfactorL[i][j-1], nrfactorL[i][j], nrfactorL[i][j+1], \
|
||||
nrfactorL[i+1][j-1], nrfactorL[i+1][j], nrfactorL[i+1][j+1]);
|
||||
//median = nrfactorL[i][j];
|
||||
} else {
|
||||
median = nrfactorL[i][j];
|
||||
}*/
|
||||
hipass[0] = nrfctrave * (data_fine->L[i][j] - data_coarse->L[i][j]);
|
||||
//hipass[0] = median*(data_fine->L[i][j]-data_coarse->L[i][j]);
|
||||
//hipass[0] = nrfactorL[i][j]*(data_fine->L[i][j]-data_coarse->L[i][j]);
|
||||
data_fine->L[i][j] = CLIP(hipass[0] + data_coarse->L[i][j]);
|
||||
|
||||
//chroma
|
||||
//hipass[1] = nrfactorab[i][j]*(data_fine->a[i][j]-data_coarse->a[i][j]);
|
||||
//hipass[2] = nrfactorab[i][j]*(data_fine->b[i][j]-data_coarse->b[i][j]);
|
||||
|
||||
//data_fine->a[i][j] = hipass[1]+data_coarse->a[i][j];
|
||||
//data_fine->b[i][j] = hipass[2]+data_coarse->b[i][j];
|
||||
}
|
||||
}//end of luminance correction
|
||||
|
||||
|
||||
|
||||
}//end of pitch=1
|
||||
|
||||
} else {//pitch>1
|
||||
|
||||
LabImage* smooth;
|
||||
|
||||
smooth = new LabImage(width, height);
|
||||
#ifdef _OPENMP
|
||||
#pragma omp parallel
|
||||
#endif
|
||||
|
||||
{
|
||||
|
||||
#ifdef _OPENMP
|
||||
#pragma omp for
|
||||
#endif
|
||||
|
||||
for(int i = 0; i < height; i += pitch) {
|
||||
int ix = i / pitch;
|
||||
|
||||
for(int j = 0, jx = 0; j < width; j += pitch, jx++) {
|
||||
|
||||
//copy common pixels
|
||||
smooth->L[i][j] = data_coarse->L[ix][jx];
|
||||
smooth->a[i][j] = data_coarse->a[ix][jx];
|
||||
smooth->b[i][j] = data_coarse->b[ix][jx];
|
||||
}
|
||||
}
|
||||
|
||||
//if (pitch>1) {//pitch=2; step (1) expand coarse image, fill in missing data
|
||||
#ifdef _OPENMP
|
||||
#pragma omp for
|
||||
#endif
|
||||
|
||||
for(int i = 0; i < height - 1; i += 2)
|
||||
for(int j = 0; j < width - 1; j += 2) {
|
||||
//do midpoint first
|
||||
double norm = 0.0, wtdsum[3] = {0.0, 0.0, 0.0};
|
||||
|
||||
//wtdsum[0]=wtdsum[1]=wtdsum[2]=0.0;
|
||||
for(int ix = i; ix < min(height, i + 3); ix += 2)
|
||||
for (int jx = j; jx < min(width, j + 3); jx += 2) {
|
||||
wtdsum[0] += smooth->L[ix][jx];
|
||||
wtdsum[1] += smooth->a[ix][jx];
|
||||
wtdsum[2] += smooth->b[ix][jx];
|
||||
norm++;
|
||||
}
|
||||
|
||||
norm = 1 / norm;
|
||||
smooth->L[i + 1][j + 1] = wtdsum[0] * norm;
|
||||
smooth->a[i + 1][j + 1] = wtdsum[1] * norm;
|
||||
smooth->b[i + 1][j + 1] = wtdsum[2] * norm;
|
||||
}
|
||||
|
||||
#ifdef _OPENMP
|
||||
#pragma omp for
|
||||
#endif
|
||||
|
||||
for(int i = 0; i < height - 1; i += 2)
|
||||
for(int j = 0; j < width - 1; j += 2) {
|
||||
//now right neighbor
|
||||
if (j + 1 == width) {
|
||||
continue;
|
||||
}
|
||||
|
||||
double norm = 0.0, wtdsum[3] = {0.0, 0.0, 0.0};
|
||||
|
||||
for (int jx = j; jx < min(width, j + 3); jx += 2) {
|
||||
wtdsum[0] += smooth->L[i][jx];
|
||||
wtdsum[1] += smooth->a[i][jx];
|
||||
wtdsum[2] += smooth->b[i][jx];
|
||||
norm++;
|
||||
}
|
||||
|
||||
for (int ix = max(0, i - 1); ix < min(height, i + 2); ix += 2) {
|
||||
wtdsum[0] += smooth->L[ix][j + 1];
|
||||
wtdsum[1] += smooth->a[ix][j + 1];
|
||||
wtdsum[2] += smooth->b[ix][j + 1];
|
||||
norm++;
|
||||
}
|
||||
|
||||
norm = 1 / norm;
|
||||
smooth->L[i][j + 1] = wtdsum[0] * norm;
|
||||
smooth->a[i][j + 1] = wtdsum[1] * norm;
|
||||
smooth->b[i][j + 1] = wtdsum[2] * norm;
|
||||
|
||||
//now down neighbor
|
||||
if (i + 1 == height) {
|
||||
continue;
|
||||
}
|
||||
|
||||
norm = 0.0;
|
||||
wtdsum[0] = wtdsum[1] = wtdsum[2] = 0.0;
|
||||
|
||||
for (int ix = i; ix < min(height, i + 3); ix += 2) {
|
||||
wtdsum[0] += smooth->L[ix][j];
|
||||
wtdsum[1] += smooth->a[ix][j];
|
||||
wtdsum[2] += smooth->b[ix][j];
|
||||
norm++;
|
||||
}
|
||||
|
||||
for (int jx = max(0, j - 1); jx < min(width, j + 2); jx += 2) {
|
||||
wtdsum[0] += smooth->L[i + 1][jx];
|
||||
wtdsum[1] += smooth->a[i + 1][jx];
|
||||
wtdsum[2] += smooth->b[i + 1][jx];
|
||||
norm++;
|
||||
}
|
||||
|
||||
norm = 1 / norm;
|
||||
smooth->L[i + 1][j] = wtdsum[0] * norm;
|
||||
smooth->a[i + 1][j] = wtdsum[1] * norm;
|
||||
smooth->b[i + 1][j] = wtdsum[2] * norm;
|
||||
|
||||
}
|
||||
|
||||
#ifdef _OPENMP
|
||||
#pragma omp for
|
||||
#endif
|
||||
|
||||
// step (2-3-4)
|
||||
for( int i = 0; i < height; i++)
|
||||
for(int j = 0; j < width; j++) {
|
||||
|
||||
float tonefactor = (nrwt_l[smooth->L[i][j]]);
|
||||
//double wtdsum[3], norm;
|
||||
float hipass[3], hpffluct[3], nrfactor;
|
||||
|
||||
hipass[1] = data_fine->a[i][j] - smooth->a[i][j];
|
||||
hipass[2] = data_fine->b[i][j] - smooth->b[i][j];
|
||||
|
||||
//Wiener filter
|
||||
//luma
|
||||
if (level < 2) {
|
||||
hipass[0] = data_fine->L[i][j] - smooth->L[i][j];
|
||||
hpffluct[0] = SQR(hipass[0]) + SQR(hipass[1]) + SQR(hipass[2]) + 0.001;
|
||||
nrfactorL[i][j] = (1.0 + hpffluct[0]) / (1.0 + hpffluct[0] + noisevar_L /* * Lcurve[smooth->L[i][j]]*/);
|
||||
//hipass[0] *= hpffluct[0]/(hpffluct[0]+noisevar_L);
|
||||
//data_fine->L[i][j] = CLIP(hipass[0]+smooth->L[i][j]);
|
||||
}
|
||||
|
||||
//chroma
|
||||
//hipass[1] = data_fine->a[i][j]-smooth->a[i][j];
|
||||
//hipass[2] = data_fine->b[i][j]-smooth->b[i][j];
|
||||
hpffluct[1] = SQR(hipass[1] * tonefactor) + 0.001;
|
||||
hpffluct[2] = SQR(hipass[2] * tonefactor) + 0.001;
|
||||
nrfactor = (hpffluct[1] + hpffluct[2]) / ((hpffluct[1] + hpffluct[2]) + noisevar_ab * NRWT_AB /* * abcurve[smooth->L[i][j]]*/);
|
||||
|
||||
hipass[1] *= nrfactor;
|
||||
hipass[2] *= nrfactor;
|
||||
|
||||
data_fine->a[i][j] = hipass[1] + smooth->a[i][j];
|
||||
data_fine->b[i][j] = hipass[2] + smooth->b[i][j];
|
||||
}
|
||||
|
||||
|
||||
if (level < 2) {
|
||||
#ifdef _OPENMP
|
||||
#pragma omp for
|
||||
#endif
|
||||
|
||||
for(int i = 0; i < height; i++)
|
||||
for(int j = 0; j < width; j++) {
|
||||
|
||||
float dirwt_l, norm_l;
|
||||
float nrfctrave = 0;
|
||||
norm_l = 0;//if we do want to include the input pixel in the sum
|
||||
|
||||
for(int inbr = (i - pitch); inbr <= (i + pitch); inbr += pitch) {
|
||||
if (inbr < 0 || inbr > height - 1) {
|
||||
continue;
|
||||
}
|
||||
|
||||
for (int jnbr = (j - pitch); jnbr <= (j + pitch); jnbr += pitch) {
|
||||
if (jnbr < 0 || jnbr > width - 1) {
|
||||
continue;
|
||||
}
|
||||
|
||||
dirwt_l = DIRWT_L(inbr, jnbr, i, j);
|
||||
nrfctrave += dirwt_l * nrfactorL[inbr][jnbr];
|
||||
norm_l += dirwt_l;
|
||||
}
|
||||
}
|
||||
|
||||
nrfctrave /= norm_l;
|
||||
//nrfctrave = nrfactorL[i][j];
|
||||
//nrfctrave=1;
|
||||
|
||||
|
||||
float hipass[3];
|
||||
|
||||
//luma
|
||||
|
||||
/*if (i>0 && i<height-1 && j>0 && j<width-1) {
|
||||
//med3x3(nrfactorL[i-1][j-1], nrfactorL[i-1][j], nrfactorL[i-1][j+1], \
|
||||
nrfactorL[i][j-1], nrfactorL[i][j], nrfactorL[i][j+1], \
|
||||
nrfactorL[i+1][j-1], nrfactorL[i+1][j], nrfactorL[i+1][j+1], median);
|
||||
median = hmf(nrfactorL[i-1][j-1], nrfactorL[i-1][j], nrfactorL[i-1][j+1], \
|
||||
nrfactorL[i][j-1], nrfactorL[i][j], nrfactorL[i][j+1], \
|
||||
nrfactorL[i+1][j-1], nrfactorL[i+1][j], nrfactorL[i+1][j+1]);
|
||||
} else {
|
||||
median = nrfactorL[i][j];
|
||||
}*/
|
||||
hipass[0] = nrfctrave * (data_fine->L[i][j] - smooth->L[i][j]);
|
||||
//hipass[0] = median*(data_fine->L[i][j]-smooth->L[i][j]);
|
||||
//hipass[0] = nrfactorL[i][j]*(data_fine->L[i][j]-data_coarse->L[i][j]);
|
||||
data_fine->L[i][j] = CLIP(hipass[0] + smooth->L[i][j]);
|
||||
|
||||
|
||||
//chroma
|
||||
//hipass[1] = nrfactorab[i][j]*(data_fine->a[i][j]-data_coarse->a[i][j]);
|
||||
//hipass[2] = nrfactorab[i][j]*(data_fine->b[i][j]-data_coarse->b[i][j]);
|
||||
|
||||
//data_fine->a[i][j] = hipass[1]+data_coarse->a[i][j];
|
||||
//data_fine->b[i][j] = hipass[2]+data_coarse->b[i][j];
|
||||
}
|
||||
}//end of luminance correction
|
||||
|
||||
|
||||
} // end parallel
|
||||
delete smooth;
|
||||
}//end of pitch>1
|
||||
|
||||
}
|
||||
|
||||
|
||||
#undef DIRWT_L
|
||||
#undef DIRWT_AB
|
||||
|
||||
//#undef NRWT_L
|
||||
#undef NRWT_AB
|
||||
|
||||
}
|
||||
|
||||
@@ -35,20 +35,6 @@
|
||||
#include "cplx_wavelet_dec.h"
|
||||
#include "pipettebuffer.h"
|
||||
|
||||
#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
|
||||
|
||||
|
||||
|
||||
namespace rtengine
|
||||
{
|
||||
|
||||
@@ -292,7 +278,6 @@ public:
|
||||
|
||||
// pyramid denoise
|
||||
procparams::DirPyrDenoiseParams dnparams;
|
||||
void dirpyrLab_denoise(LabImage * src, LabImage * dst, const procparams::DirPyrDenoiseParams & dnparams );//Emil's directional pyramid denoise
|
||||
void dirpyr (LabImage* data_fine, LabImage* data_coarse, int level, LUTf &rangefn_L, LUTf &rangefn_ab,
|
||||
int pitch, int scale, const int luma, int chroma );
|
||||
void idirpyr (LabImage* data_coarse, LabImage* data_fine, int level, LUTf &rangefn_L, LUTf & nrwt_l, LUTf & nrwt_ab,
|
||||
|
||||
@@ -36,30 +36,20 @@
|
||||
|
||||
*/
|
||||
|
||||
#include <stdlib.h>
|
||||
#include <stdio.h>
|
||||
#include <math.h>
|
||||
#include <string.h>
|
||||
#include <cstdlib>
|
||||
#include <cstdio>
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
#include "rtengine.h"
|
||||
#include "gauss.h"
|
||||
#include "rawimagesource.h"
|
||||
#include "improcfun.h"
|
||||
#include "opthelper.h"
|
||||
#include "median.h"
|
||||
#include "StopWatch.h"
|
||||
|
||||
#define clipretinex( val, minv, maxv ) (( val = (val < minv ? minv : val ) ) > maxv ? maxv : val )
|
||||
|
||||
#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
|
||||
|
||||
|
||||
namespace
|
||||
{
|
||||
void retinex_scales( float* scales, int nscales, int mode, int s, float high)
|
||||
@@ -629,10 +619,8 @@ void RawImageSource::MSR(float** luminance, float** originalLuminance, float **e
|
||||
#endif
|
||||
|
||||
for (int i = borderL; i < hei - borderL; i++) {
|
||||
float pp[9], temp;
|
||||
|
||||
for (int j = borderL; j < wid - borderL; j++) {
|
||||
med3(luminance[i][j], luminance[i - 1][j], luminance[i + 1][j], luminance[i][j + 1], luminance[i][j - 1], luminance[i - 1][j - 1], luminance[i - 1][j + 1], luminance[i + 1][j - 1], luminance[i + 1][j + 1], tmL[i][j]); //3x3
|
||||
tmL[i][j] = median(luminance[i][j], luminance[i - 1][j], luminance[i + 1][j], luminance[i][j + 1], luminance[i][j - 1], luminance[i - 1][j - 1], luminance[i - 1][j + 1], luminance[i + 1][j - 1], luminance[i + 1][j + 1]); //3x3
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@@ -16,26 +16,100 @@
|
||||
* You should have received a copy of the GNU General Public License
|
||||
* along with RawTherapee. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
#include "rt_math.h"
|
||||
|
||||
#pragma once
|
||||
|
||||
#include <array>
|
||||
#include <algorithm>
|
||||
|
||||
#include "opthelper.h"
|
||||
|
||||
template<typename T, std::size_t N>
|
||||
T median(std::array<T, N> array)
|
||||
{
|
||||
const typename std::array<T, N>::iterator middle = array.begin() + array.size() / 2;
|
||||
std::nth_element(array.begin(), middle, array.end());
|
||||
|
||||
return
|
||||
array.size() % 2
|
||||
? ((*middle + *std::min_element(middle + 1, array.end())) / static_cast<T>(2))
|
||||
: *middle;
|
||||
}
|
||||
|
||||
template<typename T, typename... ARGS>
|
||||
T median(T arg, ARGS... args)
|
||||
{
|
||||
return median(std::array<T, sizeof...(args) + 1>{std::move(arg), std::move(args)...});
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
inline T median3(T a, T b, T c)
|
||||
{
|
||||
return std::max(std::min(a, b), std::min(c, std::max(a, b)));
|
||||
}
|
||||
|
||||
template<>
|
||||
inline vfloat median3(vfloat a, vfloat b, vfloat c)
|
||||
{
|
||||
return vmaxf(vminf(a, b), vminf(c, vmaxf(a, b)));
|
||||
}
|
||||
|
||||
// See http://stackoverflow.com/questions/480960/code-to-calculate-median-of-five-in-c-sharp
|
||||
template<typename T>
|
||||
inline T median5(T a, T b, T c, T d, T e)
|
||||
{
|
||||
if (b < a) {
|
||||
std::swap(a, b);
|
||||
}
|
||||
if (d < c) {
|
||||
std::swap(c, d);
|
||||
}
|
||||
|
||||
if (c < a) {
|
||||
std::swap(b, d);
|
||||
c = a;
|
||||
}
|
||||
|
||||
a = e;
|
||||
|
||||
if (b < a) {
|
||||
std::swap(a, b);
|
||||
}
|
||||
|
||||
if (a < c) {
|
||||
std::swap(b, d);
|
||||
a = c;
|
||||
}
|
||||
|
||||
return std::min(a, d);
|
||||
}
|
||||
|
||||
template<>
|
||||
inline vfloat median5(vfloat a, vfloat b, vfloat c, vfloat d, vfloat e)
|
||||
{
|
||||
const vfloat f = vmaxf(vminf(a, b), vminf(c, d));
|
||||
const vfloat g = vminf(vmaxf(a, b), vmaxf(c, d));
|
||||
return median3(e, f, g);
|
||||
}
|
||||
|
||||
// middle 4 of 6 elements,
|
||||
#define MIDDLE4OF6(s0,s1,s2,s3,s4,s5,d0,d1,d2,d3,d4,d5,temp) \
|
||||
{\
|
||||
d1 = min(s1,s2);\
|
||||
d2 = max(s1,s2);\
|
||||
d0 = min(s0,d2);\
|
||||
d2 = max(s0,d2);\
|
||||
temp = min(d0,d1);\
|
||||
d1 = max(d0,d1);\
|
||||
d1 = std::min(s1,s2);\
|
||||
d2 = std::max(s1,s2);\
|
||||
d0 = std::min(s0,d2);\
|
||||
d2 = std::max(s0,d2);\
|
||||
temp = std::min(d0,d1);\
|
||||
d1 = std::max(d0,d1);\
|
||||
d0 = temp;\
|
||||
d4 = min(s4,s5);\
|
||||
d5 = max(s4,s5);\
|
||||
d3 = min(s3,d5);\
|
||||
d5 = max(s3,d5);\
|
||||
temp = min(d3,d4);\
|
||||
d4 = max(d3,d4);\
|
||||
d3 = max(d0,temp);\
|
||||
d2 = min(d2,d5);\
|
||||
d4 = std::min(s4,s5);\
|
||||
d5 = std::max(s4,s5);\
|
||||
d3 = std::min(s3,d5);\
|
||||
d5 = std::max(s3,d5);\
|
||||
temp = std::min(d3,d4);\
|
||||
d4 = std::max(d3,d4);\
|
||||
d3 = std::max(d0,temp);\
|
||||
d2 = std::min(d2,d5);\
|
||||
}
|
||||
|
||||
// middle 4 of 6 elements, vectorized
|
||||
@@ -61,27 +135,27 @@ d2 = vminf(d2,d5);\
|
||||
|
||||
#define MEDIAN7(s0,s1,s2,s3,s4,s5,s6,t0,t1,t2,t3,t4,t5,t6,median) \
|
||||
{\
|
||||
t0 = min(s0,s5);\
|
||||
t5 = max(s0,s5);\
|
||||
t3 = max(t0,s3);\
|
||||
t0 = min(t0,s3);\
|
||||
t1 = min(s1,s6);\
|
||||
t6 = max(s1,s6);\
|
||||
t2 = min(s2,s4);\
|
||||
t4 = max(s2,s4);\
|
||||
t1 = max(t0,t1);\
|
||||
median = min(t3,t5);\
|
||||
t5 = max(t3,t5);\
|
||||
t0 = std::min(s0,s5);\
|
||||
t5 = std::max(s0,s5);\
|
||||
t3 = std::max(t0,s3);\
|
||||
t0 = std::min(t0,s3);\
|
||||
t1 = std::min(s1,s6);\
|
||||
t6 = std::max(s1,s6);\
|
||||
t2 = std::min(s2,s4);\
|
||||
t4 = std::max(s2,s4);\
|
||||
t1 = std::max(t0,t1);\
|
||||
median = std::min(t3,t5);\
|
||||
t5 = std::max(t3,t5);\
|
||||
t3 = median;\
|
||||
median = min(t2,t6);\
|
||||
t6 = max(t2,t6);\
|
||||
t3 = max(median,t3);\
|
||||
t3 = min(t3,t6);\
|
||||
t4 = min(t4,t5);\
|
||||
median = min(t1,t4);\
|
||||
t4 = max(t1,t4);\
|
||||
t3 = max(median,t3);\
|
||||
median = min(t3,t4);\
|
||||
median = std::min(t2,t6);\
|
||||
t6 = std::max(t2,t6);\
|
||||
t3 = std::max(median,t3);\
|
||||
t3 = std::min(t3,t6);\
|
||||
t4 = std::min(t4,t5);\
|
||||
median = std::min(t1,t4);\
|
||||
t4 = std::max(t1,t4);\
|
||||
t3 = std::max(median,t3);\
|
||||
median = std::min(t3,t4);\
|
||||
}
|
||||
|
||||
#define VMEDIAN7(s0,s1,s2,s3,s4,s5,s6,t0,t1,t2,t3,t4,t5,t6,median) \
|
||||
|
||||
@@ -1,76 +0,0 @@
|
||||
/*
|
||||
* This file is part of RawTherapee.
|
||||
*
|
||||
* Copyright (c) 2004-2010 Gabor Horvath <hgabor@rawtherapee.com>
|
||||
*
|
||||
* 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/>.
|
||||
*/
|
||||
#define MINMAX3(a,b,c,min,max) \
|
||||
{ \
|
||||
if ((a)<(b)) { \
|
||||
if ((b)<(c)) { \
|
||||
(min) = (a); \
|
||||
(max) = (c); \
|
||||
} \
|
||||
else { \
|
||||
(max) = (b); \
|
||||
if ((a)<(c)) \
|
||||
(min) = (a); \
|
||||
else \
|
||||
(min) = (c); \
|
||||
} \
|
||||
} else { \
|
||||
if ((b)>(c)) { \
|
||||
(min) = (c); \
|
||||
(max) = (a); \
|
||||
} \
|
||||
else { \
|
||||
(min) = (b); \
|
||||
if ((a)>(c)) \
|
||||
(max) = (a); \
|
||||
else \
|
||||
(max) = (c); \
|
||||
} \
|
||||
} \
|
||||
}
|
||||
|
||||
#define MIN3(a,b,c,min) \
|
||||
{ \
|
||||
if ((a)<(b)) { \
|
||||
if ((a)<(c)) \
|
||||
(min) = (a); \
|
||||
else \
|
||||
(min) = (c); \
|
||||
} else { \
|
||||
if ((b)>(c)) \
|
||||
(min) = (c); \
|
||||
else \
|
||||
(min) = (b); \
|
||||
} \
|
||||
}
|
||||
|
||||
#define MAX3(a,b,c,min) \
|
||||
{ \
|
||||
if ((a)>(b)) { \
|
||||
if ((a)>(c)) \
|
||||
(max) = (a); \
|
||||
else \
|
||||
(max) = (c); \
|
||||
} else { \
|
||||
if ((b)<(c)) \
|
||||
(max) = (c); \
|
||||
else \
|
||||
(max) = (b); \
|
||||
} \
|
||||
}
|
||||
@@ -33,6 +33,7 @@
|
||||
#include "dcp.h"
|
||||
#include "rt_math.h"
|
||||
#include "improcfun.h"
|
||||
#include "median.h"
|
||||
#ifdef _OPENMP
|
||||
#include <omp.h>
|
||||
#endif
|
||||
@@ -433,13 +434,6 @@ PIX_SORT(p[3],p[6]); PIX_SORT(p[1],p[4]); PIX_SORT(p[2],p[5]); \
|
||||
PIX_SORT(p[4],p[7]); PIX_SORT(p[4],p[2]); PIX_SORT(p[6],p[4]); \
|
||||
PIX_SORT(p[4],p[2]); median=p[4];} //a4 is the median
|
||||
|
||||
#define med5(a0,a1,a2,a3,a4,median) { \
|
||||
p[0]=a0; p[1]=a1; p[2]=a2; p[3]=a3; p[4]=a4; \
|
||||
PIX_SORT(p[0],p[1]) ; PIX_SORT(p[3],p[4]) ; PIX_SORT(p[0],p[3]) ; \
|
||||
PIX_SORT(p[1],p[4]) ; PIX_SORT(p[1],p[2]) ; PIX_SORT(p[2],p[3]) ; \
|
||||
PIX_SORT(p[1],p[2]) ; median=p[2] ;}
|
||||
|
||||
|
||||
RawImageSource::RawImageSource ()
|
||||
: ImageSource()
|
||||
, plistener(NULL)
|
||||
@@ -3018,7 +3012,7 @@ SSEFUNCTION void RawImageSource::cfaboxblur(RawImage *riFlatFile, float* cfablur
|
||||
|
||||
for (int i = 0; i < H; i++) {
|
||||
int iprev, inext, jprev, jnext;
|
||||
int p[5], temp, median;
|
||||
int med;
|
||||
|
||||
if (i < 2) {
|
||||
iprev = i + 2;
|
||||
@@ -3048,12 +3042,11 @@ SSEFUNCTION void RawImageSource::cfaboxblur(RawImage *riFlatFile, float* cfablur
|
||||
//med3x3(riFlatFile->data[iprev][jprev], riFlatFile->data[iprev][j], riFlatFile->data[iprev][jnext],
|
||||
// riFlatFile->data[i][jprev], riFlatFile->data[i][j], riFlatFile->data[i][jnext],
|
||||
// riFlatFile->data[inext][jprev], riFlatFile->data[inext][j], riFlatFile->data[inext][jnext], cfatmp[i*W+j]);
|
||||
med5(riFlatFile->data[iprev][j], riFlatFile->data[i][jprev], riFlatFile->data[i][j],
|
||||
riFlatFile->data[i][jnext], riFlatFile->data[inext][j], median);
|
||||
med = median(riFlatFile->data[iprev][j], riFlatFile->data[i][jprev], riFlatFile->data[i][j], riFlatFile->data[i][jnext], riFlatFile->data[inext][j]);
|
||||
|
||||
// if (riFlatFile->data[i][j]>hotdeadthresh*median || median>hotdeadthresh*riFlatFile->data[i][j]) {
|
||||
if (((int)riFlatFile->data[i][j] << 1) > median || (median << 1) > riFlatFile->data[i][j]) {
|
||||
cfatmp[i * W + j] = median;
|
||||
if (((int)riFlatFile->data[i][j] << 1) > med || (med << 1) > riFlatFile->data[i][j]) {
|
||||
cfatmp[i * W + j] = med;
|
||||
} else {
|
||||
cfatmp[i * W + j] = riFlatFile->data[i][j];
|
||||
}
|
||||
|
||||
@@ -11,10 +11,10 @@ static const float MAXVALF = float(MAXVAL); // float version of MAXVAL
|
||||
static const double MAXVALD = double(MAXVAL); // double version of MAXVAL
|
||||
|
||||
template <typename _Tp>
|
||||
inline const _Tp SQR (_Tp x)
|
||||
inline _Tp SQR (_Tp x)
|
||||
{
|
||||
// return std::pow(x,2); Slower than:
|
||||
return (x * x);
|
||||
return x * x;
|
||||
}
|
||||
|
||||
template<typename _Tp>
|
||||
@@ -31,13 +31,13 @@ inline const _Tp& max(const _Tp& a, const _Tp& b)
|
||||
|
||||
|
||||
template<typename _Tp>
|
||||
inline const _Tp LIM(const _Tp& a, const _Tp& b, const _Tp& c)
|
||||
inline const _Tp& LIM(const _Tp& a, const _Tp& b, const _Tp& c)
|
||||
{
|
||||
return std::max(b, std::min(a, c));
|
||||
}
|
||||
|
||||
template<typename _Tp>
|
||||
inline const _Tp LIM01(const _Tp& a)
|
||||
inline _Tp LIM01(const _Tp& a)
|
||||
{
|
||||
return std::max(_Tp(0), std::min(a, _Tp(1)));
|
||||
}
|
||||
@@ -49,7 +49,7 @@ inline const _Tp ULIM(const _Tp& a, const _Tp& b, const _Tp& c)
|
||||
}
|
||||
|
||||
template<typename _Tp>
|
||||
inline const _Tp CLIP(const _Tp& a)
|
||||
inline _Tp CLIP(const _Tp& a)
|
||||
{
|
||||
return LIM(a, static_cast<_Tp>(0), static_cast<_Tp>(MAXVAL));
|
||||
}
|
||||
@@ -80,7 +80,7 @@ inline const _Tp& max(const _Tp& a, const _Tp& b, const _Tp& c, const _Tp& d)
|
||||
}
|
||||
|
||||
template<typename _Tp>
|
||||
inline const _Tp intp(const _Tp a, const _Tp b, const _Tp c)
|
||||
inline _Tp intp(_Tp a, _Tp b, _Tp c)
|
||||
{
|
||||
// calculate a * b + (1 - a) * c
|
||||
// following is valid:
|
||||
@@ -90,19 +90,19 @@ inline const _Tp intp(const _Tp a, const _Tp b, const _Tp c)
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
T norm1(const T& x, const T& y)
|
||||
inline T norm1(const T& x, const T& y)
|
||||
{
|
||||
return std::abs(x) + std::abs(y);
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
T norm2(const T& x, const T& y)
|
||||
inline T norm2(const T& x, const T& y)
|
||||
{
|
||||
return std::sqrt(x * x + y * y);
|
||||
}
|
||||
|
||||
template< typename T >
|
||||
T norminf(const T& x, const T& y)
|
||||
inline T norminf(const T& x, const T& y)
|
||||
{
|
||||
return std::max(std::abs(x), std::abs(y));
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user