OpenPaperOpt/utilities.h

Go to the documentation of this file.

#pragma once
#ifndef UTILITIES_H
#define UTILITIES_H




#include <stdio.h>
#include <string>
#include <vector>
#include <complex>
#include <stdarg.h>
#include <functional>
#include <algorithm>
#include <math.h>
#include "V3.h"
#include <time.h>


#define MIN(a,b)   (((a)<=(b))?(a):(b))

#define MAX(a,b)   (((a)>=(b))?(a):(b))

#define ABS(a)     (((a) < 0) ? -(a) : (a))


#define FLT_QNAN   0xffc00000

#define DATA_MOVE(data,dist)    (data)=(double*)(((char*)(data))+(dist))
#define DATA_ADD(data,dist)     ((double*)(((char*)(data))+(dist)))

#pragma warning(disable:4996)

using namespace std;
//============================================================================
template <class C> class AutoDelete {
public:

        AutoDelete()           {ptr=0;}
        AutoDelete(C* p)       {ptr=p;}
        ~AutoDelete()          {delete ptr;ptr=0;}

        C* operator->()  const {return ptr;}
        C& operator*()   const {return *ptr;}
        C* Ptr()         const {return ptr;}
        operator C*()    const {return Ptr();}
        void clear()           {delete ptr;ptr=0;}
        void set(C* p)         {delete ptr;ptr=p;}
        void operator=(C* p)   {set(p);}

        C*   detatch()         {C* p=ptr;ptr=0;return p;}

private:
        AutoDelete(const AutoDelete& p) {} //Deny copy

        C* ptr;
};

//===============================================================================
template <class C> class AutoDeleteV {
public:
        AutoDeleteV()           {ptr=0;}
        AutoDeleteV(C* p)       {ptr=p;}
        ~AutoDeleteV()          {delete [] ptr;ptr=0;}

        C* operator->()   const {return ptr;}
        C& operator*()    const {return *ptr;}
        C* Ptr()          const {return ptr;}
        operator C*()     const {return Ptr();}

        void clear()            {delete [] ptr;ptr=0;}
        void set(C* p)          {delete [] ptr;ptr=p;}
        void operator=(C* p)    {set(p);}

        C*   detatch()          {C* p=ptr;ptr=0;return p;}

private:
        AutoDeleteV(const AutoDeleteV& p) {} //Deny copy

        C* ptr;
};
//=============================================================================

class StdExcept {
public:
        StdExcept() {}
        StdExcept(const char* format,...){
                char buf[2048];
                va_list(arglist);
                va_start(arglist,format);
                vsnprintf(buf,2048,format,arglist);
                mMsg=buf;
                printf("/n%s/n",buf);
        }
        string GetMessage() const {return mMsg;}
private:
        string mMsg;
};




//==============================================================================

template <int numbins>
class SimpleDist {
public:
        SimpleDist() {mMin=mBinWidth=0;memset(mF,0,sizeof(mF));}
        ~SimpleDist() {}

        int GetDataSize() {return sizeof(*this);}

        void InitFromSamples(double minVal,double maxVal,double* data,int numData,int dataDist);

        double InvF(double y); //Calc F^-1(y), y=[0..1)

private:
        double mMin,mBinWidth;//Start of distributaion and width of bins
        double mF[numbins-1]; //The f-function assume 1 in last
};

template <int numbins>
void SimpleDist<numbins>::InitFromSamples(double minVal,double maxVal,double* data,int numData,int dataDist) 
{
        //Sanitycheck parameters
        if (numData<2) throw StdExcept("SimpleDist::InitFromSamples: numData<2");
        if (minVal>=maxVal) throw StdExcept("SimpleDist::InitFromSamples: minVal(%f)>=maxVal(%f) error",minVal,maxVal);
        if (numbins<2) throw StdExcept("SimpleDist::InitFromSamples: numbins<2"); //Well, will compile away...
        if (dataDist<4) throw StdExcept("SimpleDist::InitFromSamples: invalid dataDist (%d)",dataDist);

        //Uppdate locals
        mMin      = minVal;
        mBinWidth = (maxVal-minVal)/(double)(numbins);

        //Accumulate points
        int numPnts=0;
        int bins[numbins+1]; //+1 to capture case when data==maxVal

        memset(bins,0,sizeof(bins));
        for (int i=numData;i;i--) {
                double val=*data;
                data=(double*)(((char*)data)+dataDist);
                if (val>=minVal && val<=maxVal) {
                        numPnts++;
                        int binNo=(int)((val-minVal)/mBinWidth);
                        if (binNo<0 || binNo>=numbins+1) throw StdExcept("SimpleDist<numbins>::InitFromSamples: Dodamn :)"); //TODO: Remove
                        bins[binNo]++;
                }
        }
        //Check that we caught atleast two points
        if (numPnts<2) throw StdExcept("SimpleDist::InitFromSamples: Less than 2 points in intervall");

        //Calc f-function
        double nPntInv=1.0/numPnts;
        int    pntAcc=0;

        for (int i=0;i<(numbins-1);i++) {
                pntAcc+=bins[i];
                mF[i]=pntAcc*nPntInv;
        }
}

//----------------------------------------------------------------------------------------

template <int numbins>
double SimpleDist<numbins>::InvF(double y) {

        //Binsearch for the bin we'r in
        double* f=upper_bound(mF,mF+(numbins-1),y);

        double l=(f==mF)?0.0:*(f-1);       //Left value
        double r=(f==mF+(numbins-1))?1.0:*f; //Right value

        //If no samples in bin, return midpos
        if (fabs(r-l)<1e-15)
                return ((f-mF)+0.5)*mBinWidth+mMin;

        double pos=mBinWidth*(y-l)/(r-l);  //Ofs in this bin  
        pos+=(f-mF)*mBinWidth+mMin;       //Add bin offset
        return pos;
}


//=================================================================================

typedef SimpleDist<13> ZDist;

//=================================================================================
class SubdivPatch {
public:
        SubdivPatch() {min=0;max=0;p=0;numPnt=0;}
        SubdivPatch(double mi,double ma,double* pp,int np) {min=mi;max=ma;p=pp,numPnt=np;}

        bool operator<(const SubdivPatch& p) const  {return min<p.min;}
        bool operator==(const SubdivPatch& p) const {return min==p.min;}

        double  min,max;
        double* p;
        int    numPnt;
};

//=================================================================================
class SubdivEnv {
public:
        double minVal,maxVal;
        double minDelta;
        int minPoints;
        double maxPatch;
        int dataDist;

        vector <SubdivPatch> patches;
};
//=================================================================================
class fiberData {
public:
        fiberData() {x=0;y=0;z=0;}
        fiberData(double xx,double yy,double zz) {x=xx;y=yy;z=zz;}
        bool operator<(const fiberData& f) const {return 1;}
        bool operator==(const fiberData& f) const {return 1;}

        double* getZPointer() {return &z;}

        bool compare_x(double xx) const {return x<=xx;}
        bool compare_y(double yy) const {return y<=yy;}


        double x,y,z;
};
//================================================================================
typedef vector<fiberData>  fsamples;
typedef fsamples::iterator fsamplesiter;

//Predicate-templates
template<class _Tx,class _Ty>
struct compare_x : binary_function<_Tx, _Ty, bool> {
        bool operator()(const _Tx& _X, const _Ty& _Y) const
        {return _X.compare_x(_Y);}
};
template<class _Tx,class _Ty>
struct compare_y : binary_function<_Tx, _Ty, bool> {
        bool operator()(const _Tx& _X, const _Ty& _Y) const
        {return _X.compare_y(_Y);}
};
//=================================================================================

//== Creaton of local_weight(y)
class TreeNode {
public:
        TreeNode() {value=0;};
        ~TreeNode() {};

        TreeNode* Left()    const {return left;}
        TreeNode* Right()   const {return right;}
        double    Value()   const {return value;}

        void Left(TreeNode* n)    {left=n;}
        void Right(TreeNode* n)   {right=n;}
        void Value(double d)      {value=d;}

        void operator+=(double d) {value+=d;}

private:
        AutoDelete<TreeNode> left;
        AutoDelete<TreeNode> right;
        double               value;
};

//======================================================================================

class FiberSample {
public:
        int   operator<(const FiberSample& f)  const {return d[0]<f.d[0];}
        int   operator==(const FiberSample& f) const {return d[0]==f.d[0] && d[1]==f.d[1] && d[2]==f.d[2];}
        double operator[](int idx)              const {return d[idx];}

        double* PtrToFirst()                    {return &(d[0]);}

        double d[3]; //Len,Wid,Curl
};

//=========================================================================================
namespace p3d{

        static const double gMinBaseLen   = 1000.0/8.0;    //um, stop subdividing base here
        static const double gTolerance    = 0.03f;         //%,  itterate untill within this tollerance
        static const double gLowTolerance = 0.09f;         //%,  itterate untill within this tollerance if <65% cf
        static const double gEnergyMultipler = 1.0f/500.0f;//Angular scaler //TODO: Fit this parameter to samples
        extern int GenerateCurlVector(double len,double rad,double formfac,int numSegs,double* angles,int packetsLow=25,int packetsHig=35);
}
//===========================================================================================
namespace p3d{
#define SIMVOL_FLAG_ANGLE   0x0001 
#define SIMVOL_FLAG_CURL    0x0002 
#define SIMVOL_FLAG_POS     0x0004 
#define SIMVOL_FLAG_POS3D   0x0008 
#define SEGT_ROTZ           0x0001
}



#if _MSC_VER >= 1000
#pragma once
#endif // _MSC_VER >= 1000

//#define NEW_PROP_IN_LAYER 0




#define EPS 2.220446049250313e-016
#define INF 1e48
#define PI 3.1415926535897932384626433832795
const double THICKNESS_LIMIT = 1e-4;//100*EPS;


#include <stdlib.h>
#include <iostream>
#include <complex>
using namespace std;
//Complex number representation
typedef std::complex<double> complex_LC;  



#define R2_IM1 2147483563
#define R2_IM2 2147483399
#define R2_AM (1.0/R2_IM1)
#define R2_IMM1 (R2_IM1-1)
#define R2_IA1 40014
#define R2_IA2 40692
#define R2_IQ1 53668
#define R2_IQ2 52774
#define R2_IR1 12211
#define R2_IR2 3791
#define R2_NTAB 32
#define R2_NDIV (1+R2_IMM1/R2_NTAB)
#define R2_EPS 1.2e-7
#define R2_RNMX (1.0-R2_EPS)
//long rand_idnum;
//=-347;

inline int InitialiseGenerateRandom(void){
        extern long rand_idnum;
        srand ( time(NULL) );
        rand_idnum = (long)(-rand()/double(RAND_MAX)*1000);
        cout << "Random generator initialised: " << rand_idnum << endl;
        return(1);
}

inline double GenerateRandom(void){
        extern long rand_idnum;
        // FROM JON JENSEN
        //----------------------------------------------------------------------------------------
        // From Numerical Recipies, slightly modified code, same generator.
        int j;
        long k;
        static long idum2=123456789;
        static long iy=0;
        static long iv[R2_NTAB];
        float  temp;

        if (rand_idnum <= 0) {                            //Initialize.
                if (-(rand_idnum) < 1)                          //Be sure to prevent rand_idnum = 0.
                        rand_idnum=1;
                else
                        rand_idnum = -(rand_idnum);
                idum2=(rand_idnum);
                for (j=R2_NTAB+7;j>=0;j--) {                    //Load the shuffle table (after 8 warm-ups).
                        k=(rand_idnum)/R2_IQ1;
                        rand_idnum=R2_IA1*(rand_idnum-k*R2_IQ1)-k*R2_IR1;
                        if (rand_idnum < 0)
                                rand_idnum += R2_IM1;
                        if (j < R2_NTAB)
                                iv[j] = rand_idnum;
                }
                iy=iv[0];
        }

        k=(rand_idnum)/R2_IQ1;                            //Start here when not initializing.
        rand_idnum=R2_IA1*(rand_idnum-k*R2_IQ1)-k*R2_IR1; //Compute rand_idnum=(IA1*rand_idnum) % IM1 without
        if (rand_idnum < 0)                               //overflows by Schrage's method.
                rand_idnum += R2_IM1;
        k=idum2/R2_IQ2;
        idum2=R2_IA2*(idum2-k*R2_IQ2)-k*R2_IR2;           //Compute idum2=(IA2*rand_idnum) % IM2 likewise.
        if (idum2 < 0)
                idum2 += R2_IM2;
        j=iy/R2_NDIV;                                     //Will be in the range 0..NTAB-1.
        iy=iv[j]-idum2;                                   //Here rand_idnum is shufled, rand_idnum and idum2 are
        //combined to generate output.
        iv[j] = rand_idnum;
        if (iy < 1)
                iy += R2_IMM1;
        if ((temp=(float)R2_AM*iy) > R2_RNMX)                    //Because users don't expect endpoint values.
                return (float)R2_RNMX;
        else
                return temp;
};




inline double square(complex_LC a)
{
        double temp = a.real()*a.real() + a.imag()*a.imag();
        return (temp*temp);
}

// MIN and MAX functions
inline double grMax(double n1,double n2){if (n1 > n2){return n1;} else return n2;}
inline double grMin(double n1, double n2){if (n1 < n2){return n1;} else return n2;}



int CalcQ(const std::complex<double> n1, const std::complex<double> n2,const V3<double> wp_lmn_local, const V3<double> wp_lmn_new,
                  double *Qss, double *Qsp, double *Qps, double *Qpp);

int MicroScattering2(V3<double>& wp_lmn_local, V3<double>& wp_pol_local, double wp_s_local, 
                                         double wp_p_local, complex_LC refr_i,complex_LC refr_t,V2<float> rms,
                                         double lambda, double &reflfactor);

// Implement complex arcsin
complex<double> cmplxasin(complex<double> z);



#endif