Improvements for autoexposure and 'blend' highlight reconstruction.

This commit is contained in:
Emil Martinec
2011-11-27 21:43:09 -06:00
parent 9a37da3eb9
commit ba7dddd663
2 changed files with 94 additions and 62 deletions

View File

@@ -621,7 +621,7 @@ void ImProcFunctions::colorCurve (LabImage* lold, LabImage* lnew) {
double corr = 1;//pow(2.0, defgain);//defgain may be redundant legacy of superceded code???
float scale = 65536.0;
float midgray=0.18445f;//middle gray in linear gamma = 0.18445*65535
float midgray=0.15;//0.18445f;//middle gray in linear gamma = 0.18445*65535
int imax=65536>>histcompr;
@@ -630,7 +630,7 @@ void ImProcFunctions::colorCurve (LabImage* lold, LabImage* lnew) {
//find average luminance
for (int i=0; i<imax; i++) {
sum += histogram[i];
ave += histogram[i] * i;//log((float)(i+1))/log(2.0);
ave += histogram[i] * i;
}
ave /= (sum);
@@ -646,24 +646,29 @@ void ImProcFunctions::colorCurve (LabImage* lold, LabImage* lnew) {
float octile[8]={0,0,0,0,0,0,0,0},ospread=0;
count=0;
for (int i=0; i<imax; i++) {
octile[count] += histogram[i];
if (octile[count]>sum/8) {
octile[count]=log(1+i)/log(2);
count = MIN(count+1,7);
if (count<8) {
octile[count] += histogram[i];
if (octile[count]>sum/8 || (count==7 && octile[count]>sum/16)) {
octile[count]=log(1+i)/log(2);
count++;// = MIN(count+1,7);
}
}
if (i<median) {
if (i<ave) {
//lodev += SQR(ave-i)*histogram[i];
lodev += (log(median+1)-log(i+1))*histogram[i];
lodev += (log(ave+1)-log(i+1))*histogram[i];
losum += histogram[i];
} else {
//hidev += SQR(i-ave)*histogram[i];
hidev += (log(i+1)-log(median+1))*histogram[i];
hidev += (log(i+1)-log(ave+1))*histogram[i];
hisum += histogram[i];
}
}
lodev = (lodev/(log(2)*losum));
hidev = (hidev/(log(2)*hisum));
if (octile[7]>log(imax+1)/log2(2)) {
octile[7]=1.5*octile[6]-0.5*octile[5];
}
// compute weighted average separation of octiles
// for future use in contrast setting
for (int i=1; i<6; i++) {
@@ -680,7 +685,7 @@ void ImProcFunctions::colorCurve (LabImage* lold, LabImage* lnew) {
}
//compute clipped white point
int clippable = (int)(sum * clip);
int clippable = (int)(sum * clip );
clipped = 0;
int whiteclip = (imax) - 1;
while (whiteclip>1 && histogram[whiteclip]+clipped <= clippable) {
@@ -702,85 +707,102 @@ void ImProcFunctions::colorCurve (LabImage* lold, LabImage* lnew) {
ave = ave*(1<<histcompr);
median <<= histcompr;
shc <<= histcompr;
//prevent division by 0
if (ave==0) ave=1;
if (lodev==0) lodev=1;
//compute exposure compensation as geometric mean of the amount that
//sets the median at middle gray, and the amount that sets the clipped white point
//at 65535. Apply a correction factor of median/ave so that high/low key images
//are treated properly
float gain = /*corr**/midgray*scale/(median-shc+midgray*shc);
//sets the mean or median at middle gray, and the amount that sets the estimated top
//of the histogram at or near clipping.
gain = (median/ave)*sqrt(gain*scale/rawmax);
float expcomp1 = log(/*(median/ave)*//*(hidev/lodev)*/midgray*scale/(ave-shc+midgray*shc))/log(2);
float expcomp2 = 0.5*( (15.5f-histcompr-(2*octile[7]-octile[6])) + log(scale/rawmax)/log(2) );
/*expcomp = (expcomp1*fabs(expcomp2)+expcomp2*fabs(expcomp1))/(fabs(expcomp1)+fabs(expcomp2));
if (expcomp<0) {
MIN(0.0f,MAX(expcomp1,expcomp2));
}*/
expcomp = 0.5 * (expcomp1 + expcomp2);
float gain = exp(expcomp*log(2));
gain = /*(median/ave)*/sqrt(gain*scale/rawmax);
black = shc*gain;
expcomp = log(gain)/log(2.0);//convert to stops
black = shc*gain;
//now tune hlcompr to bring back rawmax to 65535
hlcomprthresh = 33;
float shoulder = ((scale/MAX(1,gain))*(hlcomprthresh/200.0))+0.1;
//this is a series approximation of the actual formula for comp,
//which is a transcendental equation
float comp = (gain*whiteclip/scale - 1)*2*(1-shoulder/scale);
float comp = (gain*((float)whiteclip)/scale - 1)*2;//*(1-shoulder/scale);
hlcompr=(int)(100*comp/(MAX(0,expcomp) + 1.0));
hlcompr = MAX(0,MIN(100,hlcompr));
//now find brightness if gain didn't bring ave to midgray using
//the envelope of the actual 'control cage' brightness curve for simplicity
float midtmp = gain*(median)/scale;
float midtmp = gain*sqrt(median*ave)/scale;
if (midtmp<0.1) {
bright = (midgray-midtmp)*15.0/(midtmp);
} else {
bright = (midgray-midtmp)*15.0/(0.10833-0.0833*midtmp);
}
bright = 0.25*(median/ave)*(hidev/lodev)*MAX(0,bright);
bright = 0.25*/*(median/ave)*(hidev/lodev)*/MAX(0,bright);
//compute contrast that spreads the average spacing of octiles
contr = 50.0*(1.1-ospread);
contr = MAX(0,MIN(100,contr));
//diagnostics
/*printf ("**************** AUTO LEVELS ****************\n");
printf ("median: %i\n",median);
printf ("average: %f\n",ave);
printf ("median/average: %f\n",median/ave);
printf ("hidev/lodev: %f\n",hidev/lodev);
printf ("lodev: %f\n",lodev);
printf ("hidev: %f\n",hidev);
printf ("octiles: %f %f %f %f %f %f %f\n",octile[0],octile[1],octile[2],octile[3],octile[4],octile[5],octile[6]);
printf ("ospread= %f\n",ospread);*/
printf ("**************** AUTO LEVELS ****************\n");
printf ("gain1= %f gain2= %f gain= %f\n",expcomp1,expcomp2,gain);
printf ("median: %i average: %f median/average: %f\n",median,ave, median/ave);
//printf ("average: %f\n",ave);
//printf ("median/average: %f\n",median/ave);
printf ("lodev: %f hidev: %f hidev/lodev: %f\n",lodev,hidev,hidev/lodev);
//printf ("lodev: %f\n",lodev);
//printf ("hidev: %f\n",hidev);
printf ("rawmax= %d whiteclip= %d gain= %f\n",rawmax,whiteclip,gain);
printf ("octiles: %f %f %f %f %f %f %f %f\n",octile[0],octile[1],octile[2],octile[3],octile[4],octile[5],octile[6],octile[7]);
printf ("ospread= %f\n",ospread);
/*
// %%%%%%%%%% LEGACY AUTOEXPOSURE CODE %%%%%%%%%%%%%
// black point selection is based on the linear result (yielding better visual results)
black = (int)(shc * corr);
// compute the white point of the exp. compensated gamma corrected image
double whiteclipg = (int)(CurveFactory::gamma2 (whiteclip * corr / 65536.0) * 65536.0);
// %%%%%%%%%% LEGACY AUTOEXPOSURE CODE %%%%%%%%%%%%%
// black point selection is based on the linear result (yielding better visual results)
black = (int)(shc * corr);
// compute the white point of the exp. compensated gamma corrected image
double whiteclipg = (int)(CurveFactory::gamma2 (whiteclip * corr / 65536.0) * 65536.0);
// compute average intensity of the exp compensated, gamma corrected image
double gavg = 0;
for (int i=0; i<65536>>histcompr; i++)
gavg += histogram[i] * CurveFactory::gamma2((int)(corr*(i<<histcompr)<65535 ? corr*(i<<histcompr) : 65535)) / sum;
if (black < gavg) {
int maxwhiteclip = (gavg - black) * 4 / 3 + black; // dont let whiteclip be such large that the histogram average goes above 3/4
//double mavg = 65536.0 / (whiteclipg-black) * (gavg - black);
if (whiteclipg < maxwhiteclip)
whiteclipg = maxwhiteclip;
}
whiteclipg = CurveFactory::igamma2 ((float)(whiteclipg/65535.0)) * 65535.0; //need to inverse gamma transform to get correct exposure compensation parameter
black = (int)((65535*black)/whiteclipg);
expcomp = log(65535.0 / (whiteclipg)) / log(2.0);
if (expcomp<0.0) expcomp = 0.0;*/
// compute average intensity of the exp compensated, gamma corrected image
double gavg = 0;
for (int i=0; i<65536>>histcompr; i++)
gavg += histogram[i] * CurveFactory::gamma2((int)(corr*(i<<histcompr)<65535 ? corr*(i<<histcompr) : 65535)) / sum;
if (black < gavg) {
int maxwhiteclip = (gavg - black) * 4 / 3 + black; // dont let whiteclip be such large that the histogram average goes above 3/4
//double mavg = 65536.0 / (whiteclipg-black) * (gavg - black);
if (whiteclipg < maxwhiteclip)
whiteclipg = maxwhiteclip;
}
whiteclipg = CurveFactory::igamma2 ((float)(whiteclipg/65535.0)) * 65535.0; //need to inverse gamma transform to get correct exposure compensation parameter
black = (int)((65535*black)/whiteclipg);
expcomp = log(65535.0 / (whiteclipg)) / log(2.0);
if (expcomp<0.0) expcomp = 0.0;*/
if (expcomp>10.0) expcomp=10.0;
if (expcomp<-5.0) expcomp = -5.0;
if (expcomp>10.0) expcomp = 10.0;
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#include "calc_distort.h"

View File

@@ -1921,10 +1921,11 @@ void RawImageSource::HLRecovery_blend(float* rin, float* gin, float* bin, int wi
#define SQR(x) ((x)*(x))
float minpt=MIN(MIN(hlmax[0],hlmax[1]),hlmax[2]);
//float maxpt=MAX(MAX(hlmax[0],hlmax[1]),hlmax[2]);
//float medpt=hlmax[0]+hlmax[1]+hlmax[2]-minpt-maxpt;
float maxave=(hlmax[0]+hlmax[1]+hlmax[2])/3;
float minpt=MIN(MIN(hlmax[0],hlmax[1]),hlmax[2]);//min of the raw clip points
//float maxpt=MAX(MAX(hlmax[0],hlmax[1]),hlmax[2]);//max of the raw clip points
//float medpt=hlmax[0]+hlmax[1]+hlmax[2]-minpt-maxpt;//median of the raw clip points
float maxave=(hlmax[0]+hlmax[1]+hlmax[2])/3;//ave of the raw clip points
//some thresholds:
const float clipthresh = 0.95;
const float fixthresh = 0.5;
const float satthresh = 0.5;
@@ -1940,7 +1941,7 @@ void RawImageSource::HLRecovery_blend(float* rin, float* gin, float* bin, int wi
#pragma omp parallel for
for (int col=0; col<width; col++) {
float rgb[ColorCount], cam[2][ColorCount], lab[2][ColorCount], sum[2], chratio;
float rgb[ColorCount], cam[2][ColorCount], lab[2][ColorCount], sum[2], chratio, lratio=0;
float L,C,H,Lfrac;
// Copy input pixel to rgb so it's easier to access in loops
@@ -1956,6 +1957,7 @@ void RawImageSource::HLRecovery_blend(float* rin, float* gin, float* bin, int wi
// Initialize cam with raw input [0] and potentially clipped input [1]
FOREACHCOLOR {
lratio += MIN(rgb[c],clip[c]);
cam[0][c] = rgb[c];
cam[1][c] = MIN(cam[0][c],maxval);
}
@@ -1972,13 +1974,13 @@ void RawImageSource::HLRecovery_blend(float* rin, float* gin, float* bin, int wi
for (int c=1; c < ColorCount; c++)
sum[i] += SQR(lab[i][c]);
}
chratio = sqrt(sqrt(sum[1]/sum[0]));
chratio = (sqrt(sum[1]/sum[0]));
// Apply ratio to lightness in lab space
// Apply ratio to lightness in LCH space
for (int c=1; c < ColorCount; c++)
lab[0][c] *= chratio;
// Transform back from lab to RGB
// Transform back from LCH to RGB
FOREACHCOLOR {
cam[0][c]=0;
for (int j=0; j < ColorCount; j++) {
@@ -2001,6 +2003,14 @@ void RawImageSource::HLRecovery_blend(float* rin, float* gin, float* bin, int wi
bin[col]= MIN(maxave,bfrac*rgb[2]+(1-bfrac)*bin[col]);
}
lratio /= (rin[col]+gin[col]+bin[col]);
L = (rin[col]+gin[col]+bin[col])/3;
C = lratio * 1.732050808 * (rin[col] - gin[col]);
H = lratio * (2 * bin[col] - rin[col] - gin[col]);
rin[col] = L - H / 6.0 + C / 3.464101615;
gin[col] = L - H / 6.0 - C / 3.464101615;
bin[col] = L + H / 3.0;
if ((L=(rin[col]+gin[col]+bin[col])/3) > desatpt) {
Lfrac = MAX(0,(maxave-L)/(maxave-desatpt));
C = Lfrac * 1.732050808 * (rin[col] - gin[col]);