https://doi.org/10.5201/ipol.2015.108
demo_fluttercode.cpp
/*demo_fluttercode.cpp*/
/*
* Copyright 2014 IPOL Image Processing On Line http://www.ipol.im/
*
* This program 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.
*
* This program 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 this program. If not, see <http://www.gnu.org/licenses/>.
*/
/**
* @mainpage Flutter shutter code optimizer
*
* README.txt:
* @verbinclude README.txt
*/
/**
* @file demo_fluttercode.cpp
* @brief Flutter shutter code optimizer main file,
* that implements algorithm 1.
* @author Yohann Tendero <tendero@math.ucla.edu>
*/
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <math.h>
#include <fstream>
#ifndef ABS
/**
* absolute value definition.
*/
#define ABS(x) (((x) > 0) ? (x) : (-(x)))
#endif
#ifndef M_PI
/**
* M_PI is a POSIX definition for Pi
*/
#define M_PI 3.14159265358979323846
#endif
#include "flutter_optimiser_gaussian.h"
#include "flutter_optimiser_uniform.h"
#include "gain_evaluation.h"
#include "optimal_snapshot.h"
using namespace std;
/**
* @brief main function call
*
INPUTs:
*Motion model : 0: centered Gaussian (0,s) => argv[1]
1: centered Uniform on [-s,s]
* the s parameter.
*N : length of the code
*Deltat : time steps
Outputs
txt file containing the a_k.
EPS file containing the stairs a_k
EPS file containing the FT of the flutter function
EPS file containing the comparative average_gain_in_terms_of_RMSE of the
code and the corresoinding snapshot for velocity_max.
*
*/
int main(int argc, char **argv)
{
//Check arguments : IN OUT;
if (argc < 8)
{
std::cerr << " ***************************************** " << std::endl
<< " ********** Flutter shutter code calculator ****** " << std::endl
<< " EXAMPLE: ./flutter_optimizer 0 0.25 52 5 code.txt " << std::endl
<< " TF_code.txt rendement.txt mean_sdt.txt " << std::endl
<< " Usage: input : motion _model and parameter******** " << std::endl
<< " Code_length, Exposure time factor multiplication*** " << std::endl
<< " code_length exposure_time_factor code_file_name**** " << std::endl
<< "Fourier_transform_file_name efficency_file_name***** " << std::endl
<< " efficency_rist_file_name*************************** " << std::endl
<< " *************************************************** " << std::endl
<< " ************* Yohann Tendero, 2014 ************** " << std::endl
<< " *************************************************** " << std::endl;
return 1;
}
/**
* @brief : Overall description:
* Step 1 : read parameter;
* Step 2 : Given the motion model : compute optimal snapshot on average
* (the optimal exposure time to use in a standard camera to ensure the
* optimal MSE aking the deconvolution into account);
* Step 3 : Given the motion model : compute the optimal flutter shutter
* camera design (the optimal aperture strategy);
* Step 4-5 : store results.
*
**/
////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
//// Step1. Reading arguments, init ////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
int flag_motion_model=atoi(argv[1]);
/// 0 For Gaussian, 1 for uniform, 2 for handcrafted
double motion_model_parameter=ABS(atof(argv[2]));
/// motion model is N(0,\sigma) \sigma=motion_model_parameter;
/// if Gaussian and [-motion_model_parameter,motion_model_parameter]
/// if uniform.
int code_length=atoi(argv[3]); // length of the code L
double exposure_factor=atof(argv[4]);// tune the support of the
// flutter shutter function as a factor of the optimal snapshot
// exposure time.
double* code = new double[code_length]; //to store the flutter shutter code
const int num_plot_points=300; // number of points for plots
double T_star; // to store the optimal exposure time (std camera)
double velocity_max=motion_model_parameter;// maximum velocity.
if (flag_motion_model==0) velocity_max=4*motion_model_parameter;
// storing files names for output.
char* code_file_name=argv[5];
char* TF_code_file_name=argv[6];
char* SNR_comparison_file_name=argv[7];
char* mean_sdt_file_name=argv[8];
////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
//// Step2 optimal snapshot average computation ///////////////////////////
////////////////////////////////////////////////////////////////////////////
T_star=optimal_snapshot(flag_motion_model,
motion_model_parameter);
double deltat=exposure_factor*T_star/code_length;
if (deltat*velocity_max>1)
{
std::cerr << " Deltat too big! convergence not guaranteed"
<< std::endl;
std::cerr << " Increase the code length L." << std::endl;
}
////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
//// Step3 Code computation ///////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
if (flag_motion_model==0) //Gaussian motion model
{
gaussian_optimiser(code,motion_model_parameter,code_length,deltat);
}
if (flag_motion_model==1) //Uniform motion model
{
uniform_optimiser(code,motion_model_parameter,code_length,deltat);
}
//Recall that at this point the code is normalized so that \int\alpha(t)dt=1
////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
//// Step4 : Writing the code file /////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
///WRITE the code to file : code_file_name
ofstream file_code(code_file_name, ios::out | ios::trunc);
// opening the file (erase if previous exists)
if (file_code)
{
for (int k=0; k<code_length; k++)
{
file_code << code[k] <<";" << endl;
}
file_code.close();
}
else
cerr << "Unable to open file !" << endl;
///WRITE the Fourier transfor of the code to file : TF_code_file_name :
///This one contains (columns by columns)
/// \xi for \xi \in [-\2pi velocity_max, 2 pi \velocitymax ]
/// The ideal Fourier Transform
/// The Fourier transform of the code
/// Flutter shutter codes are defined up to a constant renormalization,
/// the following renormalize so \hat \alpha(0)=\sqrt[4] \W(0)=1.
double w_0=0;
if (flag_motion_model==0) w_0=gaussian_w(motion_model_parameter,0);
if (flag_motion_model==1) w_0=uniform_w(motion_model_parameter,0);
// at this point w_0 contains \W(0).
ofstream file_TFcode(TF_code_file_name, ios::out | ios::trunc);
if (file_TFcode)
{
if (flag_motion_model==0) // Gaussian motion model
for (double xi=-2*velocity_max*M_PI; xi<=2*velocity_max*M_PI;
xi=xi+4*velocity_max*M_PI/num_plot_points)
{
file_TFcode << xi
<< " "
<< pow(gaussian_w(motion_model_parameter, xi)/w_0,0.25)
<<" "
<< abs_hat_alpha(code, code_length, xi,deltat) << endl;
}
if (flag_motion_model==1) // Uniform motion model
for (double xi=-2*velocity_max*M_PI; xi<=2*velocity_max*M_PI;
xi=xi+4*velocity_max*M_PI/num_plot_points)
{
file_TFcode << xi
<< " "
<< pow(uniform_w(motion_model_parameter, xi)/w_0,0.25)
<<" "
<< abs_hat_alpha(code, code_length, xi,deltat) << endl;
}
file_TFcode.close();
}
else
cerr << "Unable to open file !" << endl;
////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
//// Step5 : Code comparison with snapshot ////////////////////////////////
////////////////////////////////////////////////////////////////////////////
// Compute the relative efficiency, the relative efficiency on average,
// the std_dev_gain_sigma and store them.
int i=0;
double* average_gain_RMSE = new double[num_plot_points];
double* average_gain_mu = new double[1];
double* std_dev_gain_sigma = new double[1];
if (flag_motion_model==0) velocity_max=4*motion_model_parameter;
average_gain_in_terms_of_RMSE( code, code_length, T_star,
flag_motion_model, motion_model_parameter,
deltat,num_plot_points,
average_gain_RMSE,
average_gain_mu);
std_dev_gain_in_terms_of_RMSE(flag_motion_model,
motion_model_parameter,
num_plot_points,
average_gain_RMSE,
average_gain_mu,
std_dev_gain_sigma);
ofstream file_comparison_file_name(SNR_comparison_file_name,
ios::out | ios::trunc);
// opening the file (erase if previous exists)
if (file_comparison_file_name)
{
double velocity=-velocity_max;
for (i=0; i<num_plot_points; i++)
{
file_comparison_file_name << velocity
<< " "
<< average_gain_RMSE[i]
<< " "
<< average_gain_mu[0]
<< " "
<< probability_density(velocity, flag_motion_model,
motion_model_parameter)
<< " "
<< "1"
<< endl;
velocity=velocity+2*velocity_max/num_plot_points;
}
file_comparison_file_name.close();
}
else
cerr << "Unable to open file !" << endl;
printf("Maximum blur support for optimal snapshot \
(\\Delta t^* |velocity_{max}|) in pixel(s): %f\n",
T_star*velocity_max);
ofstream file_mean_sdt_file_name(mean_sdt_file_name, ios::out | ios::trunc);
// opening the file (erase if previous exists)
if (file_mean_sdt_file_name)
{
file_mean_sdt_file_name
<< "Gain with respect to the optimal snapshot \
in terms of RMSE, summary" << std::endl
<< "*average gain: " <<
average_gain_mu[0] << std::endl
<<"*Standard deviation: " << std_dev_gain_sigma[0] << endl;
file_mean_sdt_file_name.close();
}
else
cerr << "Unable to open file !" << endl;
delete [] code;
delete [] average_gain_RMSE;
delete [] average_gain_mu;
delete [] std_dev_gain_sigma;
return 0;
}
