gwo.h

// Copyright (C) Tomoyoshi Shimobaba 2011-


#ifndef _GWO_H
#define _GWO_H

#include "cwo.h"
#include "cwoCu.h"
#include "gwo_lib.h"

#include <string>
#include <assert.h> 

//#ifndef _GWO_LIB_H
//#define _GWO_LIB_H
//#include "gwo_lib.h"

class GWO : public CWO
{
        int stream_mode;
        float *p_rnd;

public:
        //cwoStream cstream;
        gwoCtx gtx;

        GWO();
        ~GWO();

        GWO(GWO &tmp);
        GWO(int Nx, int Ny);
        
#ifdef __NOW_DEBUGGIBNG_OPERATOR
        template<class T> friend GWO operator+(GWO& const lhs, T& const rhs){ printf("ope1 + \n"); GWO nrv = lhs; nrv += rhs; return nrv; };
        template<class T> friend GWO&& operator+(GWO&& lhs, T& const rhs){ printf("ope2 + \n"); lhs += rhs; return std::move(lhs); };
        template<class T> friend GWO&& operator+(const GWO&  lhs, T&& rhs){ printf("ope3 + \n");        rhs += lhs;     return std::move(rhs); };
        template<class T> friend GWO&& operator+(GWO&& lhs, T&& rhs){ printf("ope4 + \n"); lhs += std::move(rhs); return std::move(lhs); };

        template<class T> friend GWO operator-(GWO& const lhs, T& const rhs){ printf("ope1 - \n"); GWO nrv = lhs; nrv -= rhs; return nrv; };
        template<class T> friend GWO&& operator-(GWO&& lhs, T& const rhs){ printf("ope2 - \n"); lhs -= rhs; return std::move(lhs); };
        template<class T> friend GWO&& operator-(const T&  lhs, GWO&& rhs){ printf("ope3 - \n");        rhs -= lhs;     return std::move(rhs); };
        template<class T> friend GWO&& operator-(GWO&& lhs, T&& rhs){ printf("ope4 - \n"); lhs -= std::move(rhs); return std::move(lhs); };

        template<class T> friend GWO operator*(GWO& const lhs, T& const rhs){ printf("ope1 * \n"); GWO nrv = lhs; nrv *= rhs; return nrv; };
        template<class T> friend GWO&& operator*(GWO&& lhs, T& const rhs){ printf("ope2 * \n"); lhs *= rhs; return std::move(lhs); };
        template<class T> friend GWO&& operator*(const T&  lhs, GWO&& rhs){ printf("ope3 * \n");        rhs *= lhs;     return std::move(rhs); };
        template<class T> friend GWO&& operator*(GWO&& lhs, T&& rhs){ printf("ope4 * \n"); lhs *= std::move(rhs); return std::move(lhs); };

        template<class T> friend GWO operator/(GWO& const lhs, T& const rhs){ printf("ope1 / \n"); GWO nrv = lhs; nrv /= rhs; return nrv; };
        template<class T> friend GWO&& operator/(GWO&& lhs, T& const rhs){ printf("ope2 / \n"); lhs /= rhs; return std::move(lhs); };
        template<class T> friend GWO&& operator/(const T&  lhs, GWO&& rhs){ printf("ope3 / \n");        rhs /= lhs;     return std::move(rhs); };
        template<class T> friend GWO&& operator/(GWO&& lhs, T&& rhs){ printf("ope4 / \n"); lhs /= std::move(rhs); return std::move(lhs); };
#endif

        GWO& operator=(GWO& tmp)
        {
                size_t size = tmp.GetMemSizeCplx();

                assert(tmp.GetBuffer() != nullptr);

                if (GetNx() != tmp.GetNx() || GetNy() != tmp.GetNy() || GetBuffer() == nullptr) {
                        SetSize(tmp.GetNx(), tmp.GetNy());
                        __Free((void**)&p_field);
                        p_field = __Malloc(size);
                        assert(p_field != nullptr);
                }
                else {
                        SetSize(tmp.GetNx(), tmp.GetNy());
                }
                __Memcpy(GetBuffer(), tmp.GetBuffer(), size);

                ctx = tmp.ctx;  //copy context

                return *this;
        }

        GWO& operator=(float a)
        {
                cwoComplex tmp;
                CWO_RE(tmp) = a;
                CWO_IM(tmp) = 0.0f;
                Fill(tmp);
                return *this;
        }

        GWO& operator=(const cwoComplex &a)
        {
                Fill(a);
                return *this;
        }

        template<class T = GWO> T operator+(T &a) {
                T tmp;
                tmp.Create(GetNx(), GetNy());
                assert(tmp.GetBuffer() != nullptr);
                (tmp).ctx = ctx;
                __Add(GetBuffer(), a.GetBuffer(), tmp.GetBuffer());
                return tmp;
        }
        template<class T = GWO>  T operator+(float a) {
                T tmp;
                tmp.Create(GetNx(), GetNy());
                assert(tmp.GetBuffer() != nullptr);
                (tmp).ctx = ctx;
                cwoComplex cplx;
                CWO_RE(cplx) = a; CWO_IM(cplx) = 0.0f;
                __Add(GetBuffer(), cplx, tmp.GetBuffer());
                return tmp;
        }
        template<class T = GWO>  T operator+(const cwoComplex &a) {
                T tmp;
                tmp.Create(GetNx(), GetNy());
                assert(tmp.GetBuffer() != nullptr);
                (tmp).ctx = ctx;
                __Add(GetBuffer(), a, tmp.GetBuffer());
                return (tmp);
        }
        template<class T> friend GWO operator+(const T &a, GWO& c) { GWO tmp = c; tmp += a; return tmp; };

        template<class T = GWO> T operator-(T &a) {
                T tmp;
                tmp.Create(GetNx(), GetNy());
                assert(tmp.GetBuffer() != nullptr);
                (tmp).ctx = ctx;
                __Sub(GetBuffer(), a.GetBuffer(), tmp.GetBuffer());
                return tmp;
        }
        template<class T = GWO>  T operator-(float a) {
                T tmp;
                tmp.Create(GetNx(), GetNy());
                assert(tmp.GetBuffer() != nullptr);
                (tmp).ctx = ctx;
                cwoComplex cplx;
                CWO_RE(cplx) = a; CWO_IM(cplx) = 0.0f;
                __Sub(GetBuffer(), cplx, tmp.GetBuffer());
                return tmp;
        }
        template<class T = GWO>  T operator-(const cwoComplex &a) {
                T tmp;
                tmp.Create(GetNx(), GetNy());
                assert(tmp.GetBuffer() != nullptr);
                (tmp).ctx = ctx;
                __Sub(GetBuffer(), a, tmp.GetBuffer());
                return (tmp);
        }
        template<class T> friend GWO operator-(const T &a, GWO& c) { GWO tmp(c.GetNx(), c.GetNy()); tmp = a; tmp -= c; return tmp; };

        template<class T = GWO> T operator*(T &a) {
                T tmp;
                tmp.Create(GetNx(), GetNy());
                assert(tmp.GetBuffer() != nullptr);
                (tmp).ctx = ctx;
                __Mul(GetBuffer(), a.GetBuffer(), tmp.GetBuffer());
                return tmp;
        }
        template<class T = GWO>  T operator*(float a) {
                T tmp;
                tmp.Create(GetNx(), GetNy());
                assert(tmp.GetBuffer() != nullptr);
                (tmp).ctx = ctx;
                cwoComplex cplx;
                CWO_RE(cplx) = a; CWO_IM(cplx) = 0.0f;
                __Mul(GetBuffer(), cplx, tmp.GetBuffer());
                return tmp;
        }
        template<class T = GWO>  T operator*(const cwoComplex &a) {
                T tmp;
                tmp.Create(GetNx(), GetNy());
                assert(tmp.GetBuffer() != nullptr);
                (tmp).ctx = ctx;
                __Mul(GetBuffer(), a, tmp.GetBuffer());
                return (tmp);
        }
        template<class T> friend GWO operator*(const T &a, GWO& c) { GWO tmp = c; tmp *= a; return tmp; };

        template<class T = GWO> T operator/(T &a) {
                T tmp;
                tmp.Create(GetNx(), GetNy());
                assert(tmp.GetBuffer() != nullptr);
                (tmp).ctx = ctx;
                __Div(GetBuffer(), a.GetBuffer(), tmp.GetBuffer());
                return tmp;
        }
        template<class T = GWO>  T operator/(float a) {
                T tmp;
                tmp.Create(GetNx(), GetNy());
                assert(tmp.GetBuffer() != nullptr);
                (tmp).ctx = ctx;
                cwoComplex cplx;
                CWO_RE(cplx) = a; CWO_IM(cplx) = 0.0f;
                __Div(GetBuffer(), cplx, tmp.GetBuffer());
                return tmp;
        }
        template<class T = GWO> T operator/(const cwoComplex &a) {
                T tmp;
                tmp.Create(GetNx(), GetNy());
                assert(tmp.GetBuffer() != nullptr);
                (tmp).ctx = ctx;
                __Div(GetBuffer(), a, tmp.GetBuffer());
                return (tmp);
        }
        template<class T> friend GWO operator/(const T &a, GWO& c) {
                GWO tmp;
                tmp.Create(c.GetNx(), c.GetNy());
                assert(tmp.GetBuffer() != nullptr);
                (tmp).ctx = c.ctx;
                c.__Rep(c.GetBuffer(), a, tmp.GetBuffer());
                return tmp;
        };

        //CWO& operator=(float a) { return *this; }

        //GWO(GWO&& a) {//!< move constructor
        //      *this = std::move(a);
        //}

        //GWO& operator=(GWO&& a) {//!< move constructor
        //      if (p_field) {
        //              delete p_field;
        //      }
        //      p_field = a.p_field;
        //      a.p_field = nullptr;
        //      return *this;
        //}


        /*friend GWO operator+(GWO& const lhs, GWO& const rhs)
        {
                printf("gpu ope 1\n");
                GWO nrv;
                nrv = lhs;
                nrv += rhs;
                return nrv;
        }

        friend GWO&& operator+(GWO&& lhs, GWO& const rhs)
        {
                printf("gpu ope 2\n");
                lhs += rhs;
                return std::move(lhs);
        }

        friend CWO&& operator+(GWO& const lhs, GWO&& rhs)
        {
                printf("gpu ope 3\n");
                rhs += lhs;
                return std::move(rhs);
        }

        friend GWO&& operator+(GWO&& lhs, GWO&& rhs)
        {
                printf("gpu ope 4\n");
                lhs += std::move(rhs);
                return std::move(lhs);
        }


        friend GWO operator-(GWO& const lhs, GWO& const rhs)
        {
                printf("gpu ope 1\n");
                GWO nrv;
                nrv = lhs;
                nrv -= rhs;
                return nrv;
        }

        friend GWO&& operator-(GWO&& lhs, GWO& const rhs)
        {
                printf("gpu ope 2\n");
                lhs -= rhs;
                return std::move(lhs);
        }

        friend CWO&& operator-(GWO& const lhs, GWO&& rhs)
        {
                printf("gpu ope 3\n");
                rhs -= lhs;
                return std::move(rhs);
        }

        friend GWO&& operator-(GWO&& lhs, GWO&& rhs)
        {
                printf("gpu ope 4\n");
                lhs -= std::move(rhs);
                return std::move(lhs);
        }*/

        void SetDev(int dev);
        void SetThreads(int Nx, int Ny);
        int GetThreadsX();
        int GetThreadsY();

        //void Create(int dev, int Nx, int Ny); 
        //int Create(int Nx, int Ny=1, int Nz=1); 
        void Delete();

        void Send(CWO &a);
        void Recv(CWO &a);

        void CreateStream();
        cwoStream GetStream();
        void DestroyStream();
        void SyncStream();

/*      void SetMemcpyFlag(int i=0);
        int GetMemcpyFlag();
        void SetPagelockFlag(int i=0);
        int GetPagelockFlag();
*/      void SetStreamMode(int mode);
        int GetStreamMode();

        int Load(const std::string& fname, int c=CWO_GREY);
        int Load(const std::string& fname_amp, const std::string& fname_pha, int c = CWO_GREY);
        int Save(const std::string& fname, int bmp_8_24 = 24);
        int Save(const std::string& fname, CWO *r, CWO *g = nullptr, CWO *b = nullptr);
        //int SaveMonosToColor(char* fname, char *r_name, char *g_name, char *b_name);
        int SaveAsImage(const std::string& fname, float i1, float i2, float o1, float o2, int flag = CWO_SAVE_AS_RE);
        int SaveAsImage(const std::string& fname, int flag = CWO_SAVE_AS_RE, CWO *r = nullptr, CWO *g = nullptr, CWO *b = nullptr);
        int SaveAsImage(cwoComplex *p, int Nx, int Ny, const std::string& fname, int flag = CWO_SAVE_AS_RE, int bmp_8_24 = 24);
        
public:
        void* __Malloc(size_t size);
        void __Free(void **a);
        void __Memcpy(void *dst, void *src, size_t size);
        void __Memset(void *p, int c, size_t  size);

        void Fill(cwoComplex pix);

        /*void __Expand(void *src, int srcNx, int srcNy, 
                                                void *dst, int dstNx, int dstNy, int type);*/
        void __Expand(
                void *src, int sx, int sy, int srcNx, int srcNy, 
                void *dst, int dx, int dy, int dstNx, int dstNy,
                int type);


        void __ShiftedFresnelAperture(cwoComplex *a);
        void __ShiftedFresnelProp(cwoComplex *a);
        void __ShiftedFresnelCoeff(cwoComplex *a);

        void __ARSSFresnelAperture(cwoComplex *a);
        void __ARSSFresnelProp(cwoComplex *a);
        void __ARSSFresnelCoeff(cwoComplex *a);

        void __FresnelConvProp(cwoComplex *a);
        void __FresnelConvCoeff(cwoComplex *a, float const_val=1.0f);

        void __AngularProp(cwoComplex *a, int flag);
        //void __ShiftedAngularProp(cwoComplex *a);

        void __HuygensProp(cwoComplex *a);
//      void __HuygensCoeff(cwoComplex *a);

        void __FresnelFourierProp(cwoComplex *a);
        void __FresnelFourierCoeff(cwoComplex *a);
        
        void __FresnelDblAperture(cwoComplex *a, float z1);
        void __FresnelDblFourierDomain(cwoComplex *a, float z1, float z2, cwoInt4 *zp);
        void __FresnelDblCoeff(cwoComplex *a, float z1, float z2);

        void __FFT(void *src, void *dst, int type);
        void __IFFT(void *src, void *dst);
        void __FFTShift(void *src);

        void __NUFFT_T1(cwoComplex *p_fld, cwoFloat2 *p_x, int R=2, int Msp=12);
        void __NUFFT_T2(cwoComplex *p_fld, cwoFloat2 *p_x, int R=2, int Msp=12);

        void __Add(cwoComplex *a, const cwoComplex &b, cwoComplex *c);//c=a+b
        void __Add(cwoComplex *a, cwoComplex *b, cwoComplex *c);//c=a+b
        void __Sub(cwoComplex *a, const cwoComplex &b, cwoComplex *c);//c=a-b
        void __Sub(cwoComplex *a, cwoComplex *b, cwoComplex *c);//c=a-b
        void __Mul(cwoComplex *a, const cwoComplex &b, cwoComplex *c);//c=a*b
        void __Mul(cwoComplex *a, cwoComplex *b, cwoComplex *c);//c=a*b
        void __Div(cwoComplex *a, const cwoComplex &b, cwoComplex *c);//c=a/b
        void __Div(cwoComplex *a, cwoComplex *b, cwoComplex *c);//c=a/b
        void __Rep(cwoComplex *a, float b, cwoComplex *c);//Reciprocal number c=b/a
        void __Rep(cwoComplex *a, const cwoComplex &b, cwoComplex *c);//Reciprocal number c=b/a
        void __Rep(cwoComplex *a, cwoComplex *b, cwoComplex *c);//Reciprocal number c=b/a

        void __AddSphericalWave(cwoComplex *p, float x, float y, float z, float px, float py, float a);
        void __MulSphericalWave(cwoComplex *p, float x, float y, float z, float px, float py, float a);
        void __AddApproxSphWave(cwoComplex *p, float x, float y, float z, float zx, float zy, float px, float py, float a);
        void __MulApproxSphWave(cwoComplex *p, float x, float y, float z, float zx, float zy, float px, float py, float a);

        void __Re(cwoComplex*a , cwoComplex *b);
        void __Im(cwoComplex*a , cwoComplex *b);
        void __Intensity(cwoComplex*a , cwoComplex *b);
        void __Amp(cwoComplex*a , cwoComplex *b);
        void __Phase(cwoComplex*a , cwoComplex *b, float offset);
        void __Arg(cwoComplex *a , cwoComplex *b, float scale, float offset);
        void __Real2Complex(float *src, cwoComplex *dst);
        void __Phase2Complex(float *src, cwoComplex *dst);
        void __Arg2Cplx(cwoComplex *src, cwoComplex *dst, float scale, float offset);
        void __Polar(float *amp, float *ph, cwoComplex *c);
        void __ReIm(cwoComplex *re, cwoComplex *im, cwoComplex *c);
        //      void __Gamma(float *src, float gamma);
        void __RectFillInside(cwoComplex *p, int x, int y, int Sx, int Sy, cwoComplex a);
        void __RectFillOutside(cwoComplex *p, int x, int y, int Sx, int Sy, cwoComplex a);

        void __CircFillInside(cwoComplex *p, int m, int n, int r, const cwoComplex a);
        void __CircFillOutside(cwoComplex *p, int m, int n, int r, const cwoComplex a);

        void __FloatToChar(char *dst, float *src, int N);
        void __CharToFloat(float *dst, char *src, int N);

        void __Copy(
                        cwoComplex *src, int x1, int y1, int sNx, int sNy,
                        cwoComplex *dst, int x2, int y2, int dNx, int dNy, 
                        int Sx, int Sy);

        void FFTShift();

        void SqrtReal();
        void SqrtCplx();

        void Log(float base = 10.0f, float eps = 1.0f);
        void Gamma(float g);
        void Threshold(float max, float min=0.0);
        void __PickupFloat(float *src, float *pix_p, float pix);
        void __PickupCplx(cwoComplex *src, cwoComplex *pix_p, float pix);
        void Binary(float th = 0.0, float max = 1.0, float min = 0.0);


        float Average();
        float Variance();


        void SetRandSeed(long long s);
        void RandReal(float max = 1.0f, float min = 0.0f);
        void __RandPhase(cwoComplex *a, float max, float min);
        void __MulRandPhase(cwoComplex *a, float max, float min);


        void __WindowFlatCos(cwoComplex *p, int L, int L0, float k);


        void __MaxMin(cwoComplex *a, float *max, float *min, int *max_x = nullptr, int *max_y = nullptr, int *min_x = nullptr, int *min_y = nullptr);

        cwoComplex TotalSum();



        int __ScaleReal(float i1, float i2, float o1, float o2);
        int __ScaleCplx(float i1, float i2, float o1, float o2);


        void __ResizeNearest(
                cwoComplex *p_new, int dNx, int dNy, cwoComplex *p_old, int sNx, int sNy);
        void __ResizeLinear(
                cwoComplex *p_new, int dNx, int dNy, cwoComplex *p_old, int sNx, int sNy);
        void __ResizeCubic(
                cwoComplex *p_new, int dNx, int dNy, cwoComplex *p_old, int sNx, int sNy);
        void __ResizeLanczos(
                cwoComplex *p_new, int dNx, int dNy, cwoComplex *p_old, int sNx, int sNy);


        void __ResizeNearestScale(cwoComplex *p_src, cwoComplex *p_dst, float scale_x, float scale_y);
        void __ResizeLinearScale(cwoComplex *p_src, cwoComplex *p_dst, float scale_x, float scale_y);
        void __ResizeCubicScale(cwoComplex *p_src, cwoComplex *p_dst, float scale_x, float scale_y);
        void __ResizeLanczosScale(cwoComplex *p_src, cwoComplex *p_dst, float scale_x, float scale_y);


        void AffineAngularSpectrum(float *mat_affine, float px, float py, int flag);



        void ErrorDiffusion(CWO *a = nullptr, int flag = CWO_ED_FS);
        void ErrorDiffusionSegmented(CWO *a = nullptr, int flag = CWO_ED_FS);

        float MSE(CWO &ref);


//Test code
virtual void __ArbitFresnelDirect(
        cwoComplex *p1, cwoComplex *p2, 
        cwoFloat2 *p_x1, cwoFloat2 *p_x2, 
        float *p_d1, float *p_d2);

virtual void __ArbitFresnelCoeff(
        cwoComplex *p, cwoFloat2 *p_x2, float *p_d2);


        
};


#endif