"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

Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)