/*routines.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 routines.cpp
 * @brief Routines for estimating the probability
 * density of flutter shutter codes:
 * Functions implemented
 * - Read file (std::vector<double> fileToVector(const char *name))
 * - Fourier transform of flutter code
 * (double abs_hat_alpha(std::vector<double> code, double xi))
 * - \W' function estimation
 * double w_prime_estimator(double xi, std::vector<double> code,
                         double epsilon)
 *
 * @author Yohann Tendero <tendero@math.ucla.edu>
 */



#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <fstream>
#include <vector>
#include "routines.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




using namespace std;
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
///////////////////////////// READ FLUTTER CODE FILE ///////////////////////////
//////////////////////// STORE THE VALUES IN A VECTOR //////////////////////////
////////////////////////////////////////////////////////////////////////////////
/**
 * @fn fileToVector(const char *name)
 * @brief Given the file name oof a text file containing a flutter code
 * store the values in a vector of doubles
 * @param const char *name: the file name containing the flutter code
 * @return std::vector<double>  code: the vector containing the code values.
 * \alpha_k
 */
std::vector<double> fileToVector(const char *name)
{
    std::vector<double>  code;
    ifstream fin;
    fin.open(name); // open a file
    if (!fin.good())
        cerr << "file not found" << endl; // in case file not found.
    std::ifstream input (name);
    std::string lineData;

    while(getline(input, lineData))
    {
        code.push_back(::atof(lineData.c_str()));
    }

    return code;
}





////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
///////////////////////////// MODULUS OF THE CODE COMPUTATION///////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
/**
 * @fn abs_hat_alpha(double* code, double xi)
 * @brief Given a code compute the modulus of the Fourier transform of
 * the Flutter-Shutter function defined by :
 * \alpha(t)=\sum_{k=0}^{L-1} \alpha_k \mathbb{1}_{[k\Deltat,(k+1)\Deltat[}(t)
 * This corresponds to the Algo 2, step 1.
 * @param std::vector<double> code: vector of doubles where code[k]
 * contains the \alpha_k;
 * @param double xi : frequency;
 * @return |\hat \alpha (\xi)|, see equation (2)
 */
double abs_hat_alpha(std::vector<double> code, double xi)
{
    ///Gives back the modulus of the Fourier transform
    //initialization
    double re_abs_hat_alpha=0; //real part
    double im_abs_hat_alpha=0; //imaginary part
    //Main loop
    //cout << "code.size" << code.size() << "\n" ;
    for (unsigned k=0; k<code.size(); k++)
    { //cout << "code de k" << code[k] << "\n" ;
        im_abs_hat_alpha =im_abs_hat_alpha+code[k]*sin(-xi*(k+0.5));
        re_abs_hat_alpha =re_abs_hat_alpha+code[k]*cos(-xi*(k+0.5));
    }

    if (ABS(xi)>0) //avoiding 0/0
    {
        im_abs_hat_alpha =im_abs_hat_alpha*sin(xi/2)
                          /(xi/2);//ELSE =1;
        re_abs_hat_alpha =re_abs_hat_alpha*sin(xi/2)
                          /(xi/2); //ELSE =1;
    }

    return(pow(im_abs_hat_alpha*im_abs_hat_alpha+
                      re_abs_hat_alpha*re_abs_hat_alpha,0.5));
    ///RETURNS the modulus at xi of the flutter shutter function given its code:
    /// \hat \alpha(\xi)=  sinc(\frac{ \xi}{2\pi})
    /// \sum_{k=0}^{L-1} \alpha_k e^{-i\xi   \frac{2k+1}{2}}
}





////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
///////////////////////////// Function W prime estimator ///////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
/**
 * @fn w_prime_estimator(double xi, std::vector<double> code,
 * double deltat, double epsilon)
 * @brief Given a code and the time step estimate \W'(\xi)
 * by numerical differenciation \frac{\W(\xi+\epsilon)-\W(\xi)}{\epsilon}
 * This corresponds to Algo 2, steps 2-3.
 * @param double xi : frequency;
 * @param std::vector<double> code : vector containing the flutter
 * code coefficients \alpha_k
 * @param double epsilon: precision for numerical differenciation of \W
 * @return \W'(\xi), see algorithm 2, step 2, theorem 3 equation (26)
 */
double w_prime_estimator(double xi, std::vector<double> code,
                         double epsilon)
{
    if (ABS(xi)>0) //avoiding 0/0
    {return((pow(abs_hat_alpha(code, xi+epsilon),4)*
            pow(sin((xi+epsilon)/2)/((xi+epsilon)/2),2)-
            (pow(abs_hat_alpha(code, xi),4))*// this is \frac 1 sinc
            pow(sin(xi/2)/(xi/2),2))/epsilon);}
     return((pow(abs_hat_alpha(code, xi+epsilon),4)-
            (pow(abs_hat_alpha(code, xi),4)))/epsilon);
}