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

Init retinex add tone-mapping

Browse Source
This commit is contained in:
Desmis
2015-11-07 18:47:37 +01:00
parent af0572629e
commit 1892bbf772
10 changed files with 428 additions and 242 deletions
Download Patch File Download Diff File Expand all files Collapse all files

533
rtengine/ipretinex.cc
View File

@@ -217,10 +217,12 @@ void RawImageSource::MSR(float** luminance, float** originalLuminance, float **e
float gain2 = (float) deh.gain / 100.f; //def =1 not use
gain2 = useHslLin ? gain2 * 0.5f : gain2;
float offse = (float) deh.offs; //def = 0 not use
int scal = deh.scal; //def=3
scal = 3;//disabled scal
int iter = deh.iter;
int gradient = deh.scal;
int scal = 3;//disabled scal
int nei = (int) 2.8f * deh.neigh; //def = 220
float vart = (float)deh.vart / 100.f;//variance
float gradvart = (float)deh.grad;
float strength = (float) deh.str / 100.f; // Blend with original L channel data
float limD = (float) deh.limd;
limD = pow(limD, 1.7f);//about 2500 enough
@@ -271,283 +273,342 @@ void RawImageSource::MSR(float** luminance, float** originalLuminance, float **e
if (deh.retinexMethod == "highli" || deh.retinexMethod == "highliplus") {
moderetinex = 3;
}
float aahi = 49.f / 99.f; ////reduce sensibility 50%
float bbhi = 1.f - aahi;
float high;
high = bbhi + aahi * (float) deh.highl;
retinex_scales( RetinexScales, scal, moderetinex, nei, high );
int H_L = height;
int W_L = width;
float *src[H_L] ALIGNED16;
float *srcBuffer = new float[H_L * W_L];
for (int i = 0; i < H_L; i++) {
src[i] = &srcBuffer[i * W_L];
}
#ifdef _OPENMP
#pragma omp parallel for
#endif
for (int i = 0; i < H_L; i++)
for (int j = 0; j < W_L; j++) {
src[i][j] = luminance[i][j] + eps;
luminance[i][j] = 0.f;
for(int it=1; it<iter+1; it++) {//iter nb max of iterations
float grad=1.f;
float sc = 3.f;
if(gradient==0) {
grad=1.f;
sc=3.f;
}
float *out[H_L] ALIGNED16;
float *outBuffer = new float[H_L * W_L];
for (int i = 0; i < H_L; i++) {
out[i] = &outBuffer[i * W_L];
}
const float logBetaGain = xlogf(16384.f);
float pond = logBetaGain / (float) scal;
if(!useHslLin) {
pond /= log(elogt);
}
float *buffer = new float[W_L * H_L];;
#ifdef _OPENMP
#pragma omp parallel
#endif
{
for ( int scale = scal - 1; scale >= 0; scale-- ) {
if(scale == scal - 1) {
gaussianBlur<float> (src, out, W_L, H_L, RetinexScales[scale], buffer);
} else { // reuse result of last iteration
gaussianBlur<float> (out, out, W_L, H_L, sqrtf(SQR(RetinexScales[scale]) - SQR(RetinexScales[scale + 1])), buffer);
else if(gradient==1) {
grad=0.25f*it+0.75f;
sc=-0.5f*it+4.5f;
}
else if(gradient==2) {
grad=0.5f*it+0.5f;
sc=-it+6.f;
}
else if(gradient==3) {
grad=0.666f*it+0.333f;
sc=-it+6.f;
}
else if(gradient==4) {
grad=0.8f*it+0.2f;
sc=-it+6.f;
}
else if(gradient==5) {
grad=0.9f*it+0.1f;
sc=-it+6.f;
}
else if(gradient==-1) {
grad=-0.125f*it+1.125f;
sc=3.f;
}
float varx;
float limdx, ilimdx;
if(gradvart!=0) {
if(gradvart==1) {
varx=vart*(-0.125f*it+1.125f);
limdx=limD*(-0.125f*it+1.125f);
ilimdx=1.f/limdx;
}
#ifdef __SSE2__
vfloat pondv = F2V(pond);
vfloat limMinv = F2V(ilimD);
vfloat limMaxv = F2V(limD);
#endif
#ifdef _OPENMP
#pragma omp for
#endif
for (int i = 0; i < H_L; i++) {
int j = 0;
#ifdef __SSE2__
if(useHslLin) {
for (; j < W_L - 3; j += 4) {
_mm_storeu_ps(&luminance[i][j], LVFU(luminance[i][j]) + pondv * (LIMV(LVFU(src[i][j]) / LVFU(out[i][j]), limMinv, limMaxv) ));
}
} else {
for (; j < W_L - 3; j += 4) {
_mm_storeu_ps(&luminance[i][j], LVFU(luminance[i][j]) + pondv * xlogf(LIMV(LVFU(src[i][j]) / LVFU(out[i][j]), limMinv, limMaxv) ));
}
}
#endif
if(useHslLin) {
for (; j < W_L; j++) {
luminance[i][j] += pond * (LIM(src[i][j] / out[i][j], ilimD, limD));
}
} else {
for (; j < W_L; j++) {
luminance[i][j] += pond * xlogf(LIM(src[i][j] / out[i][j], ilimD, limD)); // /logt ?
}
}
else if(gradvart==2) {
varx=vart*(-0.2f*it+1.2f);
limdx=limD*(-0.2f*it+1.2f);
ilimdx=1.f/limdx;
}
else if(gradvart==-1) {
varx=vart*(0.125f*it+0.875f);
limdx=limD*(0.125f*it+0.875f);
ilimdx=1.f/limdx;
}
else if(gradvart==-2) {
varx=vart*(0.4f*it+0.6f);
limdx=limD*(0.4f*it+0.6f);
ilimdx=1.f/limdx;
}
}
}
delete [] buffer;
delete [] outBuffer;
delete [] srcBuffer;
else {
varx=vart;
limdx=limD;
ilimdx=ilimD;
}
scal=round(sc);
float aahi = 49.f / 99.f; ////reduce sensibility 50%
float bbhi = 1.f - aahi;
float high;
high = bbhi + aahi * (float) deh.highl;
retinex_scales( RetinexScales, scal, moderetinex, nei/grad, high );
mean = 0.f;
stddv = 0.f;
// I call mean_stddv2 instead of mean_stddv ==> logBetaGain
int H_L = height;
int W_L = width;
mean_stddv2( luminance, mean, stddv, W_L, H_L, maxtr, mintr);
// printf("mean=%f std=%f delta=%f maxtr=%f mintr=%f\n", mean, stddv, delta, maxtr, mintr);
float *src[H_L] ALIGNED16;
float *srcBuffer = new float[H_L * W_L];
// mean_stddv( luminance, mean, stddv, W_L, H_L, logBetaGain, maxtr, mintr);
if (dehatransmissionCurve && mean != 0.f && stddv != 0.f) { //if curve
float asig = 0.166666f / stddv;
float bsig = 0.5f - asig * mean;
float amax = 0.333333f / (maxtr - mean - stddv);
float bmax = 1.f - amax * maxtr;
float amin = 0.333333f / (mean - stddv - mintr);
float bmin = -amin * mintr;
for (int i = 0; i < H_L; i++) {
src[i] = &srcBuffer[i * W_L];
}
asig *= 500.f;
bsig *= 500.f;
amax *= 500.f;
bmax *= 500.f;
amin *= 500.f;
bmin *= 500.f;
#ifdef _OPENMP
#pragma omp parallel for
#endif
for (int i = 0; i < H_L; i++)
for (int j = 0; j < W_L; j++) {
src[i][j] = luminance[i][j] + eps;
luminance[i][j] = 0.f;
}
float *out[H_L] ALIGNED16;
float *outBuffer = new float[H_L * W_L];
for (int i = 0; i < H_L; i++) {
out[i] = &outBuffer[i * W_L];
}
const float logBetaGain = xlogf(16384.f);
float pond = logBetaGain / (float) scal;
if(!useHslLin) {
pond /= log(elogt);
}
float *buffer = new float[W_L * H_L];;
#ifdef _OPENMP
#pragma omp parallel
#endif
{
float absciss;
for ( int scale = scal - 1; scale >= 0; scale-- ) {
if(scale == scal - 1) {
gaussianBlur<float> (src, out, W_L, H_L, RetinexScales[scale], buffer);
} else { // reuse result of last iteration
gaussianBlur<float> (out, out, W_L, H_L, sqrtf(SQR(RetinexScales[scale]) - SQR(RetinexScales[scale + 1])), buffer);
}
#ifdef __SSE2__
vfloat pondv = F2V(pond);
vfloat limMinv = F2V(ilimdx);
vfloat limMaxv = F2V(limdx);
#endif
#ifdef _OPENMP
#pragma omp for schedule(dynamic,16)
#pragma omp for
#endif
for (int i = 0; i < H_L; i++ )
for (int j = 0; j < W_L; j++) { //for mintr to maxtr evalate absciss in function of original transmission
if (LIKELY(fabsf(luminance[i][j] - mean) < stddv)) {
absciss = asig * luminance[i][j] + bsig;
} else if (luminance[i][j] >= mean) {
absciss = amax * luminance[i][j] + bmax;
} else { /*if(luminance[i][j] <= mean - stddv)*/
absciss = amin * luminance[i][j] + bmin;
for (int i = 0; i < H_L; i++) {
int j = 0;
#ifdef __SSE2__
if(useHslLin) {
for (; j < W_L - 3; j += 4) {
_mm_storeu_ps(&luminance[i][j], LVFU(luminance[i][j]) + pondv * (LIMV(LVFU(src[i][j]) / LVFU(out[i][j]), limMinv, limMaxv) ));
}
} else {
for (; j < W_L - 3; j += 4) {
_mm_storeu_ps(&luminance[i][j], LVFU(luminance[i][j]) + pondv * xlogf(LIMV(LVFU(src[i][j]) / LVFU(out[i][j]), limMinv, limMaxv) ));
}
}
luminance[i][j] *= (-1.f + 4.f * dehatransmissionCurve[absciss]); //new transmission
}
}
// median filter on transmission ==> reduce artifacts
if (deh.medianmap) {
int wid = W_L;
int hei = H_L;
float *tmL[hei] ALIGNED16;
float *tmLBuffer = new float[wid * hei];
int borderL = 1;
for (int i = 0; i < hei; i++) {
tmL[i] = &tmLBuffer[i * wid];
}
#ifdef _OPENMP
#pragma omp parallel for
#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
if(useHslLin) {
for (; j < W_L; j++) {
luminance[i][j] += pond * (LIM(src[i][j] / out[i][j], ilimdx, limdx));
}
} else {
for (; j < W_L; j++) {
luminance[i][j] += pond * xlogf(LIM(src[i][j] / out[i][j], ilimdx, limdx)); // /logt ?
}
}
}
}
#ifdef _OPENMP
#pragma omp parallel for
#endif
for (int i = borderL; i < hei - borderL; i++ ) {
for (int j = borderL; j < wid - borderL; j++) {
luminance[i][j] = tmL[i][j];
}
}
delete [] tmLBuffer;
}
delete [] buffer;
delete [] outBuffer;
delete [] srcBuffer;
mean = 0.f;
stddv = 0.f;
// I call mean_stddv2 instead of mean_stddv ==> logBetaGain
// mean_stddv( luminance, mean, stddv, W_L, H_L, 1.f, maxtr, mintr);
mean_stddv2( luminance, mean, stddv, W_L, H_L, maxtr, mintr);
// printf("mean=%f std=%f delta=%f maxtr=%f mintr=%f\n", mean, stddv, delta, maxtr, mintr);
}
float epsil = 0.1f;
mini = mean - vart * stddv;
if (mini < mintr) {
mini = mintr + epsil;
}
maxi = mean + vart * stddv;
if (maxi > maxtr) {
maxi = maxtr - epsil;
}
delta = maxi - mini;
// printf("maxi=%f mini=%f mean=%f std=%f delta=%f maxtr=%f mintr=%f\n", maxi, mini, mean, stddv, delta, maxtr, mintr);
if ( !delta ) {
delta = 1.0f;
}
float cdfactor = gain2 * 32768.f / delta;
maxCD = -9999999.f;
minCD = 9999999.f;
// mean_stddv( luminance, mean, stddv, W_L, H_L, logBetaGain, maxtr, mintr);
if (dehatransmissionCurve && mean != 0.f && stddv != 0.f) { //if curve
float asig = 0.166666f / stddv;
float bsig = 0.5f - asig * mean;
float amax = 0.333333f / (maxtr - mean - stddv);
float bmax = 1.f - amax * maxtr;
float amin = 0.333333f / (mean - stddv - mintr);
float bmin = -amin * mintr;
asig *= 500.f;
bsig *= 500.f;
amax *= 500.f;
bmax *= 500.f;
amin *= 500.f;
bmin *= 500.f;
#ifdef _OPENMP
#pragma omp parallel
#pragma omp parallel
#endif
{
float cdmax = -999999.f, cdmin = 999999.f;
{
float absciss;
#ifdef _OPENMP
#pragma omp for
#pragma omp for schedule(dynamic,16)
#endif
for ( int i = 0; i < H_L; i ++ )
for (int j = 0; j < W_L; j++) {
// float cd = cdfactor * ( luminance[i][j] * logBetaGain - mini ) + offse;
float cd = cdfactor * ( luminance[i][j] - mini ) + offse;
if(cd > cdmax) {
cdmax = cd;
}
if(cd < cdmin) {
cdmin = cd;
}
float str = strength;
if(lhutili) { // S=f(H)
{
float HH = exLuminance[i][j];
float valparam;
if(useHsl || useHslLin) {
valparam = float((shcurve->getVal(HH) - 0.5f));
} else {
valparam = float((shcurve->getVal(Color::huelab_to_huehsv2(HH)) - 0.5f));
for (int i = 0; i < H_L; i++ )
for (int j = 0; j < W_L; j++) { //for mintr to maxtr evalate absciss in function of original transmission
if (LIKELY(fabsf(luminance[i][j] - mean) < stddv)) {
absciss = asig * luminance[i][j] + bsig;
} else if (luminance[i][j] >= mean) {
absciss = amax * luminance[i][j] + bmax;
} else { /*if(luminance[i][j] <= mean - stddv)*/
absciss = amin * luminance[i][j] + bmin;
}
str *= (1.f + 2.f * valparam);
luminance[i][j] *= (-1.f + 4.f * dehatransmissionCurve[absciss]); //new transmission
}
}
// median filter on transmission ==> reduce artifacts
if (deh.medianmap && it==1) {//only one time
int wid = W_L;
int hei = H_L;
float *tmL[hei] ALIGNED16;
float *tmLBuffer = new float[wid * hei];
int borderL = 1;
for (int i = 0; i < hei; i++) {
tmL[i] = &tmLBuffer[i * wid];
}
#ifdef _OPENMP
#pragma omp parallel for
#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
}
}
// if(exLuminance[i][j] > 65535.f*hig && higplus) str *= hig;
luminance[i][j] = clipretinex( cd, 0.f, 32768.f ) * str + (1.f - str) * originalLuminance[i][j];
#ifdef _OPENMP
#pragma omp parallel for
#endif
for (int i = borderL; i < hei - borderL; i++ ) {
for (int j = borderL; j < wid - borderL; j++) {
luminance[i][j] = tmL[i][j];
}
}
delete [] tmLBuffer;
}
#ifdef _OPENMP
#pragma omp critical
#endif
{
maxCD = maxCD > cdmax ? maxCD : cdmax;
minCD = minCD < cdmin ? minCD : cdmin;
// I call mean_stddv2 instead of mean_stddv ==> logBetaGain
// mean_stddv( luminance, mean, stddv, W_L, H_L, 1.f, maxtr, mintr);
mean_stddv2( luminance, mean, stddv, W_L, H_L, maxtr, mintr);
}
}
// printf("cdmin=%f cdmax=%f\n",minCD, maxCD);
Tmean = mean;
Tsigma = stddv;
Tmin = mintr;
Tmax = maxtr;
float epsil = 0.1f;
if (shcurve) {
delete shcurve;
}
mini = mean - varx * stddv;
if (mini < mintr) {
mini = mintr + epsil;
}
maxi = mean + varx * stddv;
if (maxi > maxtr) {
maxi = maxtr - epsil;
}
delta = maxi - mini;
// printf("maxi=%f mini=%f mean=%f std=%f delta=%f maxtr=%f mintr=%f\n", maxi, mini, mean, stddv, delta, maxtr, mintr);
if ( !delta ) {
delta = 1.0f;
}
float cdfactor = gain2 * 32768.f / delta;
maxCD = -9999999.f;
minCD = 9999999.f;
#ifdef _OPENMP
#pragma omp parallel
#endif
{
float cdmax = -999999.f, cdmin = 999999.f;
#ifdef _OPENMP
#pragma omp for
#endif
for ( int i = 0; i < H_L; i ++ )
for (int j = 0; j < W_L; j++) {
// float cd = cdfactor * ( luminance[i][j] * logBetaGain - mini ) + offse;
float cd = cdfactor * ( luminance[i][j] - mini ) + offse;
if(cd > cdmax) {
cdmax = cd;
}
if(cd < cdmin) {
cdmin = cd;
}
float str = strength;
if(lhutili) { // S=f(H)
{
float HH = exLuminance[i][j];
float valparam;
if(useHsl || useHslLin) {
valparam = float((shcurve->getVal(HH) - 0.5f));
} else {
valparam = float((shcurve->getVal(Color::huelab_to_huehsv2(HH)) - 0.5f));
}
str *= (1.f + 2.f * valparam);
}
}
// if(exLuminance[i][j] > 65535.f*hig && higplus) str *= hig;
luminance[i][j] = clipretinex( cd, 0.f, 32768.f ) * str + (1.f - str) * originalLuminance[i][j];
}
#ifdef _OPENMP
#pragma omp critical
#endif
{
maxCD = maxCD > cdmax ? maxCD : cdmax;
minCD = minCD < cdmin ? minCD : cdmin;
}
}
// printf("cdmin=%f cdmax=%f\n",minCD, maxCD);
Tmean = mean;
Tsigma = stddv;
Tmin = mintr;
Tmax = maxtr;
if (shcurve) {
delete shcurve;
}
}
}
}