"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "src/camera/robertson02.cpp" between
pfstools-2.1.0.tgz and pfstools-2.2.0.tgz

About: pfstools are a set of command line programs (and one GUI program) for reading, writing, manipulating, and viewing high-dynamic range (HDR) images and video frames (similar as the netpbm package does for low-dynamic range images).

robertson02.cpp  (pfstools-2.1.0.tgz)	:	robertson02.cpp  (pfstools-2.2.0.tgz)
skipping to change at line 42		skipping to change at line 42
#include <iostream>		#include <iostream>
#include <vector>		#include <vector>
#include <cstdlib>		#include <cstdlib>
#include <math.h>		#include <math.h>
#include <responses.h>		#include <responses.h>
#include <robertson02.h>		#include <robertson02.h>
#include <iostream>		#include <iostream>
		#include <cstring>
using namespace std;		using namespace std;
#define PROG_NAME "robertson02"		#define PROG_NAME "robertson02"
// maximum iterations after algorithm accepts local minima		// maximum iterations after algorithm accepts local minima
#define MAXIT 100		#define MAXIT 100
// maximum accepted error		// maximum accepted error
#define MAX_DELTA 1e-5f		#define MAX_DELTA 1e-5f
extern bool verbose; /* verbose should be declared for each standalone code */		extern bool verbose; /* verbose should be declared for each standalone code */
float normalize_rcurve( float *rcurve, int M );		float normalize_rcurve( float *rcurve, int M );
float corr_vec3( float *a, float *b )		inline float max3( float a[3])
{		{
float ma = (a[0] + a[1] + a[2])/3.f;		float max = (a[0]>a[1]) ? a[0] : a[1];
float mb = (b[0] + b[1] + b[2])/3.f;		return (a[2]>max) ? a[2] : max;
float sa = 1e-5, sb = 1e-5, cov = 0;
for( int cc = 0; cc < 3; cc++ ) {
sa += powf(a[cc] - ma, 2.f);
sb += powf(b[cc] - mb, 2.f);
cov += (a[cc] - ma)*(b[cc] - mb);
}
return cov / sqrt(sa * sb);
}		}
float get_exposure_compensation( const Exposure &ex )		float get_exposure_compensation( const Exposure &ex )
{		{
return ex.exposure_time * ex.aperture*ex.aperture * ex.iso/3.125f / (1.0592f * 11.4f / 3.125f);		return ex.exposure_time * ex.aperture*ex.aperture * ex.iso/3.125f / (1.0592f * 11.4f / 3.125f);
}		}
int robertson02_applyResponseRGB( pfs::Array2D *rgb_out[],		int robertson02_applyResponseRGB( pfs::Array2D *rgb_out[],
const ExposureList *imgs[],		const ExposureList *imgs[],
const float *resp_curve[],		const float *resp_curve[],
const float *weights,		const float *weights,
int M,		int M,
bool deghosting )		const noise_parameters *camera,
		const float deghosting_threshold)
{		{
// number of exposures		// number of exposures
int N = imgs[0]->size();		int N = imgs[0]->size();
// frame size		// frame size
int width = rgb_out[0] -> getCols( );		int width = rgb_out[0] -> getCols( );
int height = rgb_out[0] -> getRows( );		int height = rgb_out[0] -> getRows( );
// number of saturated pixels		// number of saturated pixels
int saturated_pixels = 0;		int saturated_pixels = 0;
		float mmax[3]= {1e-30, 1e-30, 1e-30};
		/*
float pw[M];		float pw[M];
for( int i = 0 ; i < M ; i++ ) {		for( int i = 0 ; i < M ; i++ ) {
pw[i] = weights[i];		pw[i] = weights[i];
}		}
for( int i = 1; i < 10; i++ )		for( int i = 1; i < 10; i++ )
pw[M-i] = 0;		pw[M-i] = 0;
float ew[N]; // Per-exposure weight to account for noise charactetistic		float ew[N]; // Per-exposure weight to account for noise charactetistic
// For each exposure		// For each exposure
for( int i = 0 ; i < N ; i++ ) {		for( int i = 0 ; i < N ; i++ ) {
const Exposure &ex = (*imgs[0])[i];		const Exposure &ex = (*imgs[0])[i];
// TODO: Estimate read-out noise and use for exposure weighting		// TODO: Estimate read-out noise and use for exposure weighting
// TODO: Better blending if difference is weights is high		// TODO: Better blending if difference is weights is high
// This is a noise-optimial reconstruction assuming that the input is affect ed only		// This is a noise-optimial reconstruction assuming that the input is affect ed only
// by the photon noise (no read-out noise), which is most often the case		// by the photon noise (no read-out noise), which is most often the case
ew[i] = ex.exposure_time;		ew[i] = ex.exposure_time;
VERBOSE_STR << "Ex: " << i << " exp. time: " << ex.exposure_time << "\t iso: " << ex.iso		VERBOSE_STR << "Ex: " << i << " exp. time: " << ex.exposure_time << "\t iso: " << ex.iso
<< "\t aperture: " << ex.aperture << "\t weight: " << ew[i] << std::endl;		<< "\t aperture: " << ex.aperture << "\t weight: " << ew[i] << std::endl;
}		}
		*/
float X[3*N];		bool photon_noise_only = false;
bool skip_exp[N];		if( strcmp( camera->camera_name, "Empty" ) == 0 )
float w[3][N];		{
bool saturated_exp[N];		photon_noise_only = true;
		VERBOSE_STR << "Merging using the Poisson Photon Noise Estimator" << endl;
// for( int i = 0; i < M; i++ )		}
// cerr << "i = " << i << "; weight = " << weights[i] << endl;		else
		VERBOSE_STR << "Merging using the Variance Weighted Estimator" << endl;
// For each pixel		// For each pixel
int skipped_deghost = 0;		int skipped_deghost = 0;
for( int j = 0; j < width * height; j++ )		for( int j = 0; j < width * height; j++ )
{		{
bool saturatedRGB[3] = { false, false, false };		bool saturatedRGB[3] = { false, false, false };
float best_mw = -1;		// First compute scene radiances, weights and saturated exposures
int ref = -1;		float X[3][N];
		float w[3][N];
if( deghosting ) {		bool saturated_exp[N];
		fill(saturated_exp, saturated_exp + N, false);
		bool skip_exp[N];
		fill(skip_exp, skip_exp + N, false);
		float t[N];
		float g[N];
// Deghosting - find the exposure that could be used as a reference		for( int i = 0; i < N; i++ )
for( int i = 0 ; i < N ; i++ )		for( int cc = 0; cc < 3; cc++ )
{		{
saturated_exp[i] = false;		const Exposure &ex = (*imgs[cc])[i];
for( int cc=0; cc < 3; cc++ )		t[i] = ex.exposure_time;
{		g[i] = ex.iso / 100;
const Exposure &ex = (*imgs[cc])[i];
float ti = get_exposure_compensation( ex );
int m = (int) (*ex.yi)( j );		int m = (int) (*ex.yi)( j );
X[i*3+cc] = resp_curve[cc][ m ] / ti;		X[cc][i] = resp_curve[cc][ m ] / ( t[i] * g[i] * camera->color_coefficie
w[cc][i] = weights [ m ] * ew[i];		nts[cc] );
		if( photon_noise_only )
		w[cc][i] = t[i];
		else
		{
		float var = resp_curve[cc][ m ] * g[i] * camera->color_coefficients[cc
		]
		+ powf(camera->std_readout * g[i] * camera->color_coeffici
		ents[cc], 2)
		+ powf(camera->std_adc * camera->color_coefficients[cc], 2
		);
		w[cc][i] = powf( t[i] * g[i] * camera->color_coefficients[cc], 2 ) /
		var;
		}
saturated_exp[i] |= ( m <= 1 || m >= M-2 );		// Test for overexposed pixels
		saturated_exp[i] |= resp_curve[cc][m] >= resp_curve[cc][M-1];
		}
// cerr << X[i*3+cc] << ", ";		// Deghosting - find the exposure that could be used as a reference
		// Currently using middle exposure
		// TODO: Perform histogram analysis to pick best reference
		int ref = int(N/2);
		if( deghosting_threshold != -1 )
		{
		// Compute the number of pixels to be skipped due to deghosting
		for( int i = 0; i < N; i++ )
		if(i != ref && !saturated_exp[ref] && !saturated_exp[i])
		{
		// First compute the Pearson correlation coefficient wrt reference
		float ref_mean = 0, ref_var = 1e-5, cur_mean = 0, cur_var = 1e-5, cov
		= 0;
		for( int cc = 0; cc < 3; cc++ )
		{
		ref_mean += X[cc][ref];
		cur_mean += X[cc][i];
		}
		ref_mean /= 3;
		cur_mean /= 3;
		for( int cc = 0; cc < 3; cc++ )
		{
		ref_var += powf(X[cc][ref] - ref_mean, 2);
		cur_var += powf(X[cc][i] - cur_mean, 2);
		cov += (X[cc][ref] - ref_mean)*(X[cc][i] - cur_mean);
		}
}		// Next check whether total deviation from the reference is acceptable
skip_exp[i] = saturated_exp[i];		float std_sum = 0;
		float deviation_from_ref = 0;
		for( int cc = 0; cc < 3; cc++ )
		{
		// Reciprocal of the weight is the exposure compensated variance
		// Square root of this is thus the exposure compensated standard dev
		iation
		// Assuming a gaussian distribution, 95% of the population should li
		e 2 stds away from mean
		std_sum += deghosting_threshold * sqrt(1/w[cc][ref] + 1/w[cc][i]);
		deviation_from_ref += fabs( (float)(X[cc][ref] - X[cc][i]) );
		}
float mw = (w[0][i] + w[1][i] + w[2][i])/3.f;		// Skipping the Pearson coefficient test as reults are not satisfactor
if( mw > best_mw && !saturated_exp[i] ) {		y
best_mw = mw;		// Fixing a better reference might solve this
ref = i;		// float rho = cov / sqrt(ref_var/2 * cur_var/2);
		// skip_exp[i] = rho < 0.9 || deviation_from_ref > std_sum;
		skip_exp[i] = deviation_from_ref > std_sum;
		if( skip_exp[i] )
		skipped_deghost++;
}		}
// cerr << " ; ";		}
}
// Deghosting		/* // Deghosting
if( ref != -1 && !saturated_exp[ref] ) {		if( ref != -1 && !saturated_exp[ref] ) {
const Exposure &ex = (*imgs[0])[ref];		const Exposure &ex = (*imgs[0])[ref];
float ti = get_exposure_compensation( ex );		float ti = get_exposure_compensation( ex );
float Y_ref = (X[ref*3] + X[ref*3+1] + X[ref*3+2])/3.f;		float Y_ref = (X[ref*3] + X[ref*3+1] + X[ref*3+2])/3.f;
const float noise_level = 0.01;		const float noise_level = 0.01;
float var_ref = (noise_level*noise_level)*Y_ref/ti; // variance of the photon noise		float var_ref = (noise_level*noise_level)*Y_ref/ti; // variance of the photon noise
for( int i = 0 ; i < N ; i++ )		for( int i = 0 ; i < N ; i++ )
skipping to change at line 204		skipping to change at line 250
float var_exp = (noise_level*noise_level)*Y_exp / ti;		float var_exp = (noise_level*noise_level)*Y_exp / ti;
float ci = 1.96 * sqrt(var_ref+var_exp);		float ci = 1.96 * sqrt(var_ref+var_exp);
// Skip a frame if color correlation too low or difference to referenc e too high		// Skip a frame if color correlation too low or difference to referenc e too high
skip_exp[i] = ( rho < 0.9 ) || fabs( Y_exp - Y_ref ) > ci;		skip_exp[i] = ( rho < 0.9 ) || fabs( Y_exp - Y_ref ) > ci;
if( skip_exp[i] )		if( skip_exp[i] )
skipped_deghost++;		skipped_deghost++;
}		}
}		}
}
/* int skip_bits = 0;		int skip_bits = 0;
for( int i = 0 ; i < N ; i++ ) {		for( int i = 0 ; i < N ; i++ ) {
if( !skip_exp[i] )		if( !skip_exp[i] )
skip_bits++; // |= (1<<i);		skip_bits++; // |= (1<<i);
}*/		}
/* for( int i = 0 ; i < N ; i++ ) {		for( int i = 0 ; i < N ; i++ ) {
skip_exp[i] = false;		skip_exp[i] = false;
}*/		}
		*/
// cerr << endl;		// cerr << endl;
// For each color channels		// For each color channels
for( int cc=0; cc < 3; cc++ )		for( int cc=0; cc < 3; cc++ )
{		{
// all exposures for each pixel		// all exposures for each pixel
float sum = 0.0f;		float sum = 0.0f;
float div = 0.0f;		float div = 0.0f;
float div2 = 0.0f;
// For each exposure		// For each exposure
for( int i = 0 ; i < N ; i++ )		for( int i = 0 ; i < N ; i++ )
{		{
if( deghosting && (skip_exp[i] || saturated_exp[i]) )		if( saturated_exp[i] || (deghosting_threshold != -1 && skip_exp[i]) )
continue;		continue;
const Exposure &ex = (*imgs[cc])[i];		sum += X[cc][i] * w[cc][i];
int m = (int) (*ex.yi)( j );		div += w[cc][i];
float ti = get_exposure_compensation( ex );
// Testing the noise-optimal approach
/* float iso2 = ex.iso*ex.iso;
float sigma_r2 = 0.00001;
float w = ex.exposure_time*ex.exposure_time*iso2/(iso2*ex.exposure_tim
e*resp_curve[cc][m] + 2*sigma_r2);
sum += w * resp_curve[cc][ m ] / ti;
div += w;*/
// float w = weights [ m ] * ew[i];
// sum += w * resp_curve[cc][ m ] / ti;
// div += w;
// const float noise_level = 0.01;
// float var = (noise_level*noise_level)*resp_curve[cc][ m ]/(ti*ti);
// float w = 1/var;
// float w = log(sqrt(resp_curve[cc][ m ] / ti / var))+3;
// float w = log(sqrt(m));
float w = weights [ m ] * ew[i];
sum += w * resp_curve[cc][ m ] / ti;
div += w;
if( j < 10 && cc == 1 )		// if( j < 10 && cc == 1 )
std::cerr << "j = " << j << " ;i = " << i << " ;w = " << w << std::end		// std::cerr << "j = " << j << " ;i = " << i << " ;w = " << w[cc][i] <<
l;		std::endl;
div2 += weights [ m ];
}		}
/* float mean = sum / div;		if( div >= 1e-15 )
		{
sum = 0;		(*rgb_out[cc])(j) = sum / div;
div = 0;		mmax[cc] = ( mmax[cc] > (*rgb_out[cc])(j) ) ? mmax[cc] : (*rgb_out[cc])(
// For each exposure		j);
for( int i = 0 ; i < N ; i++ )		}
{		else
const Exposure &ex = (*imgs[cc])[i];		{
int m = (int) (*ex.yi)( j );		(*rgb_out[cc])(j) = -1;
float ti = ex.ti;
float ww = 1/(fabs(resp_curve[cc][ m ]/ti - mean) + 1e-7);
sum += weights [ m ] * ew[i] * resp_curve[cc][ m ] * ww / ti;
div += weights [ m ] * ew[i] * ww;
}*/
if( div2 >= 0.0001f )
(*rgb_out[cc])(j) = sum / div;
else {
saturatedRGB[cc] = true;		saturatedRGB[cc] = true;
}		}
// saturatedRGB[cc] = false;
// (*rgb_out[cc])(j) = saturated_exp[ref];
}		}
if( saturatedRGB[0] || saturatedRGB[1] || saturatedRGB[2] ) {		if( saturatedRGB[0] || saturatedRGB[1] || saturatedRGB[2] ) {
saturated_pixels++;		saturated_pixels++;
}		}
for( int cc=0; cc < 3; cc++ ) {		/* for( int cc=0; cc < 3; cc++ ) {
// TODO: Some fancy stuff to deal with distorted colors		// TODO: Some fancy stuff to deal with distorted colors
if( saturatedRGB[cc] ) {		if( saturatedRGB[cc] ) {
// If none of the exposures can give actual pixel value, make		// If none of the exposures can give actual pixel value, make
// the best (clipped) estimate using the longest or the		// the best (clipped) estimate using the longest or the
// shortest exposure;		// shortest exposure;
const Exposure &ex = (*imgs[cc])[0];		const Exposure &ex = (*imgs[cc])[0];
// If pixel value > gray level, use the shortest exposure;		// If pixel value > gray level, use the shortest exposure;
float short_long = ( (*ex.yi)(j) > M/2 ) ? 1.f : -1.f;		float short_long = ( (*ex.yi)(j) > M/2 ) ? 1.f : -1.f;
skipping to change at line 329		skipping to change at line 333
int m = (int) (*ex.yi)( j );		int m = (int) (*ex.yi)( j );
float ti = get_exposure_compensation( ex );;		float ti = get_exposure_compensation( ex );;
if( ti*short_long < best_ti*short_long ) {		if( ti*short_long < best_ti*short_long ) {
best_ti = ti;		best_ti = ti;
best_v = (int)(*ex.yi)(j);		best_v = (int)(*ex.yi)(j);
}		}
}		}
(*rgb_out[cc])(j) = 1/best_ti * resp_curve[cc][best_v];		(*rgb_out[cc])(j) = 1/best_ti * resp_curve[cc][best_v];
}		}
}		}*/
}		}
		// Fill in nan values and normalise
		for( int j = 0; j < width * height; j++ )
		for( int cc=0; cc < 3; cc++ )
		{
		if( (*rgb_out[cc])(j) == -1)
		(*rgb_out[cc])(j) = mmax[cc];
		(*rgb_out[cc])(j) /= max3(mmax);
		}
VERBOSE_STR << "Exposure pixels skipped due to deghosting: " <<		VERBOSE_STR << "Exposure pixels skipped due to deghosting: " <<
(float)skipped_deghost*100.f/(float)(width*height*N) << "%" << std::endl;		(float)skipped_deghost*100.f/(float)(width*height*N) << "%" << std::endl;
return saturated_pixels;		return saturated_pixels;
}		}
int robertson02_applyResponse( pfs::Array2D *xj,		int robertson02_applyResponse( pfs::Array2D *xj,
const ExposureList &imgs,		const ExposureList &imgs,
const float *rcurve,		const float *rcurve,
End of changes. 33 change blocks.
110 lines changed or deleted		133 lines changed or added

---