xmipp_funcs.h

Go to the documentation of this file.

/***************************************************************************
*
* Authors:     Carlos Oscar S. Sorzano (coss@cnb.csic.es)
*
* Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
*
* 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 2 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, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
* 02111-1307  USA
*
*  All comments concerning this program package may be sent to the
*  e-mail address 'xmipp@cnb.csic.es'
***************************************************************************/

#ifndef CORE_XMIPP_FUNCS_H
#define CORE_XMIPP_FUNCS_H

#include <complex>
#include "xmipp_macros.h"
#include <vector>

class FileName;

// For timing functions
// Uncomment next line timing functions are giving problems in your system
//#define _NO_TIME
#ifndef _NO_TIME
#include <unistd.h>
#ifndef __MINGW32__
#include <sys/times.h>
#endif
#ifdef _IRIX65
#include <sys/types.h>
#include <time.h>
#endif
#endif






class Tabsinc
{
public:

    double sampl;
    double isampl;
    int xmax;
    int no_elem;
    double* tabulatedsinc;

public:
    Tabsinc(const double dd, const int xx)
    {
        sampl = dd;
        isampl = 1.0/sampl;
        xmax = xx;
        filltable();
    }

    // Destructor
    virtual ~Tabsinc()
    {
        delete tabulatedsinc;
    }

#define TSINCVALUE(Tsinc, x,y) \
     { \
      int TSINCVALUEaux=(int)(x * Tsinc.isampl); \
         y=Tsinc.tabulatedsinc[ABS(TSINCVALUEaux)]; \
     }


    double operator()(double val) const
    {
        int aux=(int)(val * isampl);
        return tabulatedsinc[ABS(aux)];
    }

    void filltable()
    {
        no_elem = (int)(xmax / sampl);
        tabulatedsinc = new double[no_elem];
        tabulatedsinc[0] = 1;
        for (int i = 1; i < no_elem; i++)
        {
            double xx = (double) i * sampl * PI;
            tabulatedsinc[i] = sin(xx) / xx;
        }
    }
};

class KaiserBessel
{
protected:
    double alpha, v, r; 
    int N; 
    int K; 
    double vtable; 
    int ntable;
    std::vector<double> i0table;
    double dtable; 
    double alphar; 
    double fac; 
    double vadjust;
    double facadj; 
    void build_I0table(); 
    double fltb;

public:
    KaiserBessel(double alpha_, int K, double r_,
                 double v_, int N_, double vtable_=0.,
                 int ntable_ = 5999);

    double I0table_maxerror();
    std::vector<double> dump_table() const
    {
        return i0table;
    }

    double sinhwin(double x) const;

    double i0win(double x) const;

    inline double i0win_tab(double x) const
    {
        double xt;
        if(x<0.)
            xt = -x*fltb+0.5;
        else
            xt = x*fltb+0.5;
        return i0table[ (int) xt];
    }

    int get_window_size() const
    {
        return K;
    }
};

#if defined(__APPLE__) || defined(__MINGW32__)

void sin_cos(double angle, double * sine, double * cosine);
#endif

int solve_2nd_degree_eq(double a,
                        double b,
                        double c,
                        double& x1,
                        double& x2,
                        double prec = XMIPP_EQUAL_ACCURACY);

double gaussian1D(double x, double sigma, double mu = 0);

double tstudent1D(double x, double df, double sigma, double mu = 0);

double icdf_gauss(double p);

double cdf_gauss(double x);

double cdf_tstudent(int k, double t);

double cdf_FSnedecor(int d1, int d2, double x);

double icdf_FSnedecor(int d1, int d2, double p);

double gaussian2D(double x,
                  double y,
                  double sigmaX,
                  double sigmaY,
                  double ang,
                  double muX = 0,
                  double muY = 0);



size_t divide_equally(size_t N, size_t size, size_t rank, size_t &first, size_t &last);

size_t divide_equally_group(size_t N, size_t size, size_t myself);

template <class T>
void computeStats(const std::vector<T> &V, double& avg, double& stddev,
                  T& minval, T& maxval)
{
    if (V.size()<= 0)
        return;

    avg = 0;
    stddev = 0;

    minval = maxval = V[0];

    size_t nmax=V.size();
    const T* ptr=&V[0];
    for(size_t n=0; n<nmax; ++n, ++ptr)
    {
        double val=*ptr;
        avg += val;
        stddev += val * val;

        if (val > maxval)
            maxval = val;
        else if (val < minval)
            minval = val;
    }
    avg /= nmax;

    if (nmax > 1)
    {
        stddev = stddev / nmax - avg * avg;
        stddev *= nmax / (nmax - 1);

        // Foreseeing numerical instabilities
        stddev = sqrt(static_cast< double >(ABS(stddev)));
    }
    else
        stddev = 0;
}

template <class T>
void computeAvgStddev(const std::vector<T> &V, double& avg, double& stddev)
{
    if (V.size()<= 0)
        return;

    avg = 0;
    stddev = 0;

    size_t nmax=V.size();
    const T* ptr=&V[0];
    for(size_t n=0; n<nmax; ++n, ++ptr)
    {
        double val=*ptr;
        avg += val;
        stddev += val * val;
    }
    avg /= nmax;

    if (nmax > 1)
    {
        stddev = stddev / nmax - avg * avg;
        stddev *= nmax / (nmax - 1);

        // Foreseeing numerical instabilities
        stddev = sqrt(static_cast< double >(ABS(stddev)));
    }
    else
        stddev = 0;
}

template <class T>
void initConstant(std::vector<T> &V, T &value)
{
    const T* ptr=&V[0];
    size_t nmax=V.size();
    for(size_t n=0; n<nmax; ++n, ++ptr)
        *ptr=value;
}



template <typename T>
void RealImag2Complex(const T* _real,
                      const T* _imag,
                      std::complex< double >* _complex,
                      int length)
{
    T* aux_real = (T*) _real;
    T* aux_imag = (T*) _imag;
    double* aux_complex = (double*) _complex;

    for (int i = 0; i < length; i++)
    {
        *aux_complex++ = (double)(*aux_real++);
        *aux_complex++ = (double)(*aux_imag++);
    }
}

template <typename T>
void AmplPhase2Complex(const T* _ampl,
                       const T* _phase,
                       std::complex< double >* _complex,
                       int length)
{
    T* aux_ampl = (T*) _ampl;
    T* aux_phase = (T*) _phase;
    double* aux_complex = (double*) _complex;

    for (int i = 0; i < length; i++)
    {
        double ampl = (double)(*aux_ampl++);
        double phase = (double)(*aux_phase++);
        *aux_complex++ = ampl * cos(phase);
        *aux_complex++ = ampl * sin(phase);
    }
}

template <typename T>
void Complex2RealImag(const std::complex< double >* _complex,
                      T* _real,
                      T* _imag,
                      int length)
{
    T* aux_real = (T*) _real;
    T* aux_imag = (T*) _imag;
    double* aux_complex = (double*) _complex;

    for (int i = 0; i < length; i++)
    {
        *aux_real++ = (T)(*aux_complex++);
        *aux_imag++ = (T)(*aux_complex++);
    }
}

template <typename T>
void Complex2AmplPhase(const std::complex< double >* _complex,
                       T* _ampl,
                       T* _phase,
                       int length)
{
    T* aux_ampl = (T*) _ampl;
    T* aux_phase = (T*) _phase;
    double* aux_complex = (double*) _complex;

    for (int i = 0; i < length; i++)
    {
        double re = *aux_complex++;
        double im = *aux_complex++;
        *aux_ampl++ = sqrt(re * re + im * im);
        *aux_phase++ = atan2(im, re);
    }
}

void init_random_generator(int seed = -1);

unsigned int randomize_random_generator();

double rnd_unif();

double rnd_unif(double a, double b);

double rnd_student_t(double nu);

double rnd_student_t(double nu, double a, double b);

double rnd_gaus();

double rnd_gaus(double a, double b);

double gaus_within_x0(double x0, double mean = 0, double stddev = 1);

double gaus_outside_x0(double x0, double mean = 0, double stddev = 1);

double gaus_up_to_x0(double x0, double mean = 0, double stddev = 1);

double gaus_from_x0(double x0, double mean = 0, double stddev = 1);

double student_outside_probb(double p, double degrees_of_freedom);

double student_within_t0(double t0, double degrees_of_freedom);

double student_outside_t0(double t0, double degrees_of_freedom);

double student_up_to_t0(double t0, double degrees_of_freedom);

double student_from_t0(double t0, double degrees_of_freedom);

double chi2_up_to_t0(double t0, double degrees_of_freedom);

double chi2_from_t0(double t0, double degrees_of_freedom);

double rnd_log(double a, double b);

#if defined _NO_TIME || defined __MINGW32__
typedef int ProcessorTimeStamp; // Any other kind of data will do
typedef int TimeStamp; // Any other kind of data will do
struct tm* localtime_r (const time_t *clock, struct tm *result);
#else
typedef struct tms ProcessorTimeStamp; // Renaming of the time structure
typedef size_t TimeStamp;              // Timestamp in miliseconds
#endif

void time_config();

#if !defined _NO_TIME && !defined __MINGW32__

void annotate_processor_time(ProcessorTimeStamp* time);
#endif

void annotate_time(TimeStamp* time);

#if !defined _NO_TIME && !defined __MINGW32__

void acum_time(ProcessorTimeStamp* orig, ProcessorTimeStamp* dest);
#endif

double elapsed_time(ProcessorTimeStamp& time, bool _IN_SECS = true);

#if !defined _NO_TIME && !defined __MINGW32__

void print_elapsed_time(ProcessorTimeStamp& time, bool _IN_SECS = true);
#endif

void print_elapsed_time(TimeStamp& time, bool _IN_SECS = true);

class Timer
{
private:
  struct timeval  tv;
  struct tm tm;
  size_t tic_time;

public:
  size_t now();
  size_t tic();
  size_t toc(const char * msg=NULL, bool inSecs=true);
  size_t elapsed();
};

#if !defined _NO_TIME && !defined __MINGW32__

double time_to_go(ProcessorTimeStamp& time, double fraction_done);
#endif

void init_progress_bar(long total);

void progress_bar(long act_time);

char * getCurrentTimeString();


class BaseListener
{
public :
    BaseListener(): verbosity(0), cancel(false)
    {}
    
    virtual ~BaseListener() {}

    virtual void OnInitOperation(unsigned long _est_it) = 0;

    virtual void OnReportOperation(const std::string& _rsOp) = 0;

    virtual void OnProgress(unsigned long _it) = 0;

    virtual const unsigned& getVerbosity() const
    {
        return verbosity;
    }

    virtual unsigned& setVerbosity()
    {
        return verbosity;
    }

    virtual const bool& OnUserCancel() const
    {
        return cancel;
    }

    virtual bool& setCancel()
    {
        return cancel;
    }

private:
    unsigned verbosity;
    bool cancel;
};

class TextualListener: public BaseListener
{
public :

    virtual void OnInitOperation(unsigned long _est_it);

    virtual void OnProgress(unsigned long _it);

    virtual void OnReportOperation(const std::string& _rsOp);
};

void TimeMessage(const std::string &message);

size_t xmippFREAD(void* dest, size_t size, size_t nitems, FILE*& fp,
                  bool reverse = false);

size_t xmippFWRITE(const void* src,
                   size_t size,
                   size_t nitems,
                   FILE*& fp,
                   bool reverse = false);

void mapFile(const FileName &filename, char*&map,size_t &size, int &fileDescriptor, bool readOnly=true);

void unmapFile(char *&map, size_t &size, int& fileDescriptor);

#define little22bigendian(x) swapbytes((unsigned char*)& x,sizeof(x))

void swapbytes(char* v, unsigned long n);

bool IsBigEndian(void);

//
// process_mem_usage(double &, double &) - takes two doubles by reference,
// attempts to read the system-dependent data for a process' virtual memory
// size and resident set size, and return the results in KB.
//
// On failure, returns 0.0, 0.0
void processMemUsage(double& vm_usage, double& resident_set);

bool IsLittleEndian(void);

void printMemoryUsed();

bool compareTwoFiles(const FileName &fn1, const FileName &fn2, size_t offset = 0);



#endif