https://doi.org/10.5201/ipol.2015.108
gain_evaluation.cpp
/*gain_evaluation.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/>.
*/
/**
* @file gain_evaluation.cpp
* @brief This files constains the code that implements
* equations (21), (22), (23), (24) and (25) that are used
* to compute the gain of the flutter shutter method
* with respect to the optimal snapshot, in terms of RMSE.
* @author Yohann Tendero <tendero@math.ucla.edu>
*/
#include <math.h>
#include "gain_evaluation.h"
#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
#ifndef SQRT_TWOPI
/**
*\sqrt{2\pi} definition.
*/
#define SQRT_TWOPI 2.506628274631000
#endif
/**
* RIEMANN SUM PRECISION
*/
#define N_RESOLUTION 1000; // for Rieman sum
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////// MSE_snapshot(velocity) //////////////////////
//////////////////////////////////////EQUATION (21) ////////////////////////////
////////////////////////////////////////////////////////////////////////////////
/**
double MSE_snapshot_of_velocity(double velocity,double tilde_deltat)
* @brief This function computes the MSE of a snapshot
* for a given velocity. This is equation (21) up to a \frac 1 \pi
* factor:
* \int_{0}^{\pi} \frac{\|u\|_{L^1(\R)}}
* {T \mathrm{sinc}^2\left( \frac{\xi v T}{2 \pi}\right)}d\xi
* @param double velocity : the velocity;
* @param double tilde_deltat: the exposure time of the snapshot;
* @return The MSE of the snapshot at velcity velocity.
* See equation (21). The \frac{\|u\|_{L^1(\R)}}{\pi}\int_{0}^{\pi}
* multiplicative factor is irrelevant, see equation (23).
*/
double MSE_snapshot_of_velocity(double velocity, double tilde_deltat)
{
if (ABS(velocity)>0)
{
///Initalization.
double a=0;
double b=M_PI; // a,b integral bounds of (21).
// the following is for the Riemann sum.
double DELTA_XI=(b-a)/N_RESOLUTION;
double xi=0;
double tn=0;
for(xi = a; xi < b; xi += DELTA_XI)
{
/* Riemann sum*/
tn += MSE_snapshot_of_velocity_integrand(velocity,
tilde_deltat, xi)*DELTA_XI;
}
return(tn);
}
else // velocity=0, see remark 4 page 10.
{
return(M_PI/tilde_deltat);
}
}
/**
double MSE_snapshot_of_velocity_integrand(double velocity,
double tilde_deltat, double xi)
* @brief This function computes the integrand of equation (21):
* \frac{1}{T \mathrm{sinc}^2\left( \frac{\xi v T}{2 \pi}\right)}
* @param double velocity : the velocity;
* @param double tilde_deltat: the exposure time of the snapshot;
* @param doule xi: the frequency;
* @return The integrand of equation (21).
*/
double MSE_snapshot_of_velocity_integrand(double velocity,
double tilde_deltat, double xi)
{
double xivdt_2=xi*velocity*tilde_deltat/2;
double sinc=1;
if (ABS(xivdt_2)>0)
{
sinc=sin(xivdt_2)/xivdt_2;
}
// (else) sinc=1 since sinc(0)=1
// (see the variable initialization).
return(1/(tilde_deltat*sinc*sinc));
// this is the integrand in (21).
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////// MSE_flutter_(velocity) //////////////////////
//////////////////////////////////////EQUATION (22) ////////////////////////////
////////////////////////////////////////////////////////////////////////////////
/**
double MSE_flutter_of_velocity(double* code, int code_length,double velocity,
double tilde_deltat)
* @brief This function computes the MSE of flutter shutter
* for a given velocity. This is equation (22) up to a \frac 1 \pi
* factor:
* \int_{0}^{\pi}
* \frac{\|\alpha \|^2_{L^2(\R)} \|u\|_{L^1(\R)}}{|\hat \alpha(\xi v)|^2}d\xi
* @param double* code : array of doubles where code[k] contains the \alpha_k
* @param int code_length : length of the code L;
* @param double velocity : the velocity;
* @param double deltat : the temporal sampling of the
* flutter shutter function;
* @return The MSE of the flutter at velcity velocity.
* See equation (22). The \frac{\|u\|_{L^1(\R)}}{\pi}\int_{0}^{\pi}
* multiplicative factor is irrelevant, see equation (23).
*/
double MSE_flutter_of_velocity(double* code, int code_length, double velocity,
double deltat)
{
if (ABS(velocity)>0)
{
///Initalization.
double a=0;
double b=M_PI; // a,b integral bounds of (21).
// the following is for the Riemann sum.
double DELTA_XI=(b-a)/N_RESOLUTION;
double xi=0;
double tn=0;
for(xi = a; xi < b; xi += DELTA_XI)
{
/* Riemann sum*/
tn += MSE_flutter_of_velocity_integrand(code, code_length, velocity,
deltat, xi)*DELTA_XI;
}
return(tn);
}
else // velocity=0, see remark 4 page 10.
{
return(M_PI*MSE_flutter_of_velocity_integrand(code, code_length, 0,
deltat, 0));
}
}
/**
double MSE_flutter_of_velocity_integrand(double* code, int code_length,
double velocity, double tilde_deltat, double xi)
* @brief This function computes the integrand of equation (21):
* \frac{\|\alpha \|^2_{L^2(\R)} \|u\|_{L^1(\R)}}{|\hat \alpha(\xi v)|^2}d\xi
* @param double* code : array of doubles where code[k] contains the \alpha_k
* @param int code_length : length of the code L;
* @param double velocity : the velocity;
* @param double deltat : the temporal sampling of the
* flutter shutter function;
* @param doule xi: the frequency;
* @return The integrand of equation (21).
*/
double MSE_flutter_of_velocity_integrand(double* code, int code_length,
double velocity, double deltat, double xi)
{
///Step1 compute |\hat \alpha( \xi v)|
double modulus_hat_alpha_of_vxi=abs_hat_alpha(code, code_length,
xi*velocity,deltat);
double l2_norm_code_squared=0;
///Step2 : compute ||\alpha||_{l2}^2
for (int k=0; k<code_length; k++)
{
l2_norm_code_squared=l2_norm_code_squared+ code[k]*code[k]*deltat;
}
return(l2_norm_code_squared/(modulus_hat_alpha_of_vxi*
modulus_hat_alpha_of_vxi));
}
/**
* @fn abs_hat_alpha(double* code, int code_length, double xi, double deltat)
* @brief Given a code compute the modulus of the Fourier transform of
* the Flutter-Shutter function, see equation (2):
* \hat \alpha(\xi)=
* \Delta t \mathrm{sinc}\left(\frac{\xi \Delta t}{2\pi}\right)
* e^{-\frac{i \xi \Delta t}{2}} \sum_{k=0}^{L-1} \alpha_k e^{-ik \xi \Delta t}
* @param double* code : array of doubles where code[k]
* contains the \alpha_k;
* @param int code_length : length of the code;
* @param double xi : frequency;
* @param double deltat : the temporal sampling of the
* flutter shutter function;
* @return |\hat \alpha (\xi)| see equation (2).
*/
double abs_hat_alpha(double* code, int code_length, double xi, double deltat)
{
///Gives back the modulus of the Fourier transform, see equation (2)
//initialization
double re_abs_hat_alpha=0; //real part
double im_abs_hat_alpha=0; //imaginary part
//Main loop
for (int k=0; k<code_length; k++)
{
im_abs_hat_alpha =im_abs_hat_alpha+code[k]*sin(-xi*deltat*(k+0.5));
re_abs_hat_alpha =re_abs_hat_alpha+code[k]*cos(-xi*deltat*(k+0.5));
}
if (ABS(xi)>0) //avoiding 0/0
{
im_abs_hat_alpha =im_abs_hat_alpha*sin(xi*deltat/2)
/(xi*deltat/2);//ELSE =1;
re_abs_hat_alpha =re_abs_hat_alpha*sin(xi*deltat/2)
/(xi*deltat/2); //ELSE =1;
}
return(deltat*pow(im_abs_hat_alpha*im_abs_hat_alpha+
re_abs_hat_alpha*re_abs_hat_alpha,0.5));
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////// G(v) computation ////////////////////////////
///////////////////////////////////equation (23) ///////////////////////////////
////////////////////////////////////////////////////////////////////////////////
/**
* @fn gain_in_terms_of_RMSE_at_velocity_v(double* code, int code_length,
double tilde_deltat, double velocity, double deltat)
* @brief This function implements equation (23):
* G(v)&=\sqrt{
*\frac{\int_{0}^{\pi} \frac{1}
* {T^* \mathrm{sinc}^2\left( \frac{\xi v T^*}{2 \pi}\right)}d\xi}
* { \int_{0}^{\pi} \frac{\|\alpha \|^2_{L^2(\R)}}
* {|\hat \alpha(\xi v)|^2}d\xi }}
* @param double* code : array of doubles where code[k] contains
* the \alpha_k
* @param int code_length : length of the code L;
* @param double tilde_deltat : the \deltat of the snapshot
* (best snapshot on average);
* @param double velocity : velocity;
* @param double deltat : the temporal sampling of the
* flutter shutter function;
* @return G(v), see equation (23).
*/
double gain_in_terms_of_RMSE_at_velocity_v(double* code, int code_length,
double tilde_deltat, double velocity, double deltat)
{
return(pow(MSE_snapshot_of_velocity(velocity, tilde_deltat)/
MSE_flutter_of_velocity(code, code_length, velocity,
deltat),0.5));// This is equation (23).
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////// Computation of the average gain \mu /////////
/////////////////////////////////// equation (24) //////////////////////////////
////////////////////////////////////////////////////////////////////////////////
/**
double probability_density(double velocity, int flag_motion_model,
double motion_model_parameter)
* @brief Given a velocity velocity, the motion model and its parameters
* the function computes the probability.
* @param double velocity : the velocity;
* @param int flag_motion_model : motion model selection
* -0 for Gaussian motion model
* -1 For uniform motion model
* @param double motion_model_parameter : motion model parameter;
* @return \rho(velocity).
*/
double probability_density(double velocity, int flag_motion_model,
double motion_model_parameter)
{
if (flag_motion_model==0) //Gaussian
{
// Recall that for the Gaussian case motion_model_parameter is \sigma
return((1/(motion_model_parameter*SQRT_TWOPI)
*exp(-velocity*velocity/
(2*motion_model_parameter*motion_model_parameter))));
}
else // Uniform
// for the Uniform motion model motion_model_parameter is v_{max}
{
return(1/(2*motion_model_parameter)) ;
}
return(0);
}
/**
* @fn average_gain_in_terms_of_RMSE(double* code, int code_length,
double tilde_deltat_star,
int flag_motion_model,
double motion_model_parameter,double deltat,
int num_plot_points,
double *G_of_v,
double *average_gain_mu)
* @brief This function implements equation (24):
* \mu:=\int_\R G(v) \rho(v) dv
* @param double* code : array of doubles where code[k] contains the \alpha_k
* @param int code_length : length of the code L;
* @param double tilde_deltat_star : the optimal exposure time for the snapshot
* @param int flag_motion_model : motion model selection
* -0 for truncated Gaussian of std-dev s;
* -1 For uniform of range s;
* @param double motion_model_parameter : motion model parameter;
* @param double deltat : the \deltat of the numerical flutter shutter;
* @param int num_plot_point : the lengt of output array for plot;
* @param double *G_of_v : will contain G(v), see equation (23);
* @param double *_average_gain_mu : will contain equation (24);
* @return void.
*/
void average_gain_in_terms_of_RMSE(double* code, int code_length,
double tilde_deltat_star,
int flag_motion_model,
double motion_model_parameter,double deltat,
int num_plot_points,
double *G_of_v,
double *average_gain_mu)
{
///Initialization.
// Recall that for the Gaussian case motion_model_parameter is \sigma
// for the Uniform motion model motion_model_parameter is v_{max}
double velocity_max=motion_model_parameter;
if (flag_motion_model==0) velocity_max=4*motion_model_parameter;
//at this point velocity_max contains v_{max} for both motion models.
average_gain_mu[0]=0;
int i=0;
double sum_probabilities=0;
double probability_of_v=0;
///Computation of equation (24).
if (flag_motion_model==1) //uniform motion model
{
probability_of_v=probability_density(0, 1,
motion_model_parameter);
for (double velocity=-velocity_max; velocity<velocity_max;
velocity=velocity+2*velocity_max/num_plot_points)
{
G_of_v[i]=gain_in_terms_of_RMSE_at_velocity_v(code,
code_length, tilde_deltat_star, velocity,
deltat);
// G_of_v[0] is the gain at v=-v_{max}, etc.
sum_probabilities=sum_probabilities+probability_of_v;
// for normalization so that \sum probability_of_v=1.
average_gain_mu[0]=average_gain_mu[0]// 0 in init
+G_of_v[i]
*probability_of_v;
i++;
}
}
else // Gaussian motion model
{
for (double velocity=-velocity_max; velocity<velocity_max;
velocity=velocity+2*velocity_max/num_plot_points)
{
probability_of_v=probability_density(velocity, 0,
motion_model_parameter);
G_of_v[i]=gain_in_terms_of_RMSE_at_velocity_v(code,
code_length, tilde_deltat_star, velocity,
deltat);
// G_of_v[0] is the gain at v=-v_{max}, etc.
sum_probabilities=sum_probabilities+probability_of_v;
// for normalization so that \sum probability_of_v=1.
average_gain_mu[0]=average_gain_mu[0]// 0 in init
+G_of_v[i]
*probability_of_v;
i++;
}
}
average_gain_mu[0]=average_gain_mu[0]/sum_probabilities;
// to ensure that the mass is = 1.
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////// Computation of the dtd-dev //////////////////
/////////////////////////////////// equation (25) //////////////////////////////
////////////////////////////////////////////////////////////////////////////////
/**
* @fn std_dev_gain_in_terms_of_RMSE(int flag_motion_model,
double motion_model_parameter,
int num_plot_points,
double *G_of_v,
double *average_gain_mu,
double *std_dev_gain_sigma)
* @brief This function implements equation (25):
* \sigma:=\sqrt{\int_\R \left|G(v)-\int_\R G(u) \rho(u) du\right|^2 \rho(v) dv}.
* @param double* code : array of doubles where code[k] contains the \alpha_k
* @param int code_length : length of the code L;
* @param double tilde_deltat_star : the optimal exposure time for the snapshot
* @param int flag_motion_model : motion model selection
* -0 for truncated Gaussian of std-dev s;
* -1 For uniform of range s;
* @param double motion_model_parameter : motion model parameter;
* @param double deltat : the \deltat of the numerical flutter shutter;
* @param int num_plot_point : the lengt of output array for plot;
* @param double *G_of_v : contains G(v), see equation (23);
* @param double *_average_gain_mu : contains equation (24);
* @param double *std_dev_gain_sigma : will contains equation (25).
* @return void.
*/
void std_dev_gain_in_terms_of_RMSE(int flag_motion_model,
double motion_model_parameter,
int num_plot_points,
double *G_of_v,
double *average_gain_mu,
double *std_dev_gain_sigma)
{
///Initialization.
// Recall that for the Gaussian case motion_model_parameter is \sigma
// for the Uniform motion model motion_model_parameter is v_{max}
double velocity_max=motion_model_parameter;
if (flag_motion_model==0) velocity_max=4*motion_model_parameter;
//at this point velocity_max contains v_{max} for both motion models.
std_dev_gain_sigma[0]=0;
int i=0;
double sum_probabilities=0;
double probability_of_v=0;
///Computation of equation (25), knowing equations (23) and (24)
if (flag_motion_model==1) //uniform motion model
{
probability_of_v=1/(2*velocity_max) ;
for (double velocity=-velocity_max; velocity<velocity_max;
velocity=velocity+2*velocity_max/num_plot_points)
{
// G_of_v[0] is the gain at v=-v_{max}, etc.
sum_probabilities=sum_probabilities+probability_of_v;
// for normalization so that \sum probability_of_v=1.
std_dev_gain_sigma[0]=std_dev_gain_sigma[0]+
(G_of_v[i]-
average_gain_mu[0])*
(G_of_v[i]-
average_gain_mu[0])*
probability_of_v;
i++;
}
}
else // Gaussian motion model
{
for (double velocity=-velocity_max; velocity<velocity_max;
velocity=velocity+2*velocity_max/num_plot_points)
{
probability_of_v=probability_density(velocity, 0,
motion_model_parameter);
// G_of_v[0] is the gain at v=-v_{max}, etc.
sum_probabilities=sum_probabilities+probability_of_v;
// for normalization so that \sum probability_of_v=1.
std_dev_gain_sigma[0]=std_dev_gain_sigma[0]+
(G_of_v[i]-
average_gain_mu[0])*
(G_of_v[i]-
average_gain_mu[0])*
probability_of_v;
i++;
}
}
std_dev_gain_sigma[0]=pow(std_dev_gain_sigma[0]/
sum_probabilities,.5);
}
