CWO++library  0.32
 All Classes Functions Variables Groups
cwo.h
1 
2 #ifndef _CWO_H
3 #define _CWO_H
4 
5 
6 
7 #include "cwo_lib.h"
8 #include <math.h>
9 #include <stdarg.h>
10 
11 #include <string>
12 #include <memory>
13 
14 #include <assert.h>
15 
16 template<class T>
17 struct cwoVect2{
18  T x, y;
19  cwoVect2(){};
20  cwoVect2(T tx, T ty){ x = tx; y = ty; };
21  void Show(){ std::cout << "(" << x << "," << y << ")" << std::endl; }
22 
23  cwoVect2<T>& operator=(const cwoVect2<T> &a){
24  x = a.x; y = a.y;
25  return *this;
26  };
27  cwoVect2<T>& operator=(const T a){
28  x = a; y = 0.0f;
29  return *this;
30  };
31  cwoVect2<T>& operator+=(const cwoVect2<T> &a){
32  x += a.x; y += a.y;
33  return *this;
34  };
35  cwoVect2<T>& operator+=(const T a){
36  x += a;
37  return *this;
38  };
39  cwoVect2<T>& operator-=(const cwoVect2<T> &a){
40  x -= a.x; y -= a.y;
41  return *this;
42  };
43  cwoVect2<T>& operator-=(const T a){
44  x -= a;
45  return *this;
46  };
47  cwoVect2<T>& operator*=(const cwoVect2<T> &a){ //Hadamard product
48  x *= a.x; y *= a.y;
49  return *this;
50  }
51  cwoVect2<T>& operator*=(const T a){
52  x *= a; y *= a;
53  return *this;
54  };
55 
56  cwoVect2<T>& operator/=(const cwoVect2<T> &a){ //Hadamard product
57  x /= a.x; y /= a.y;
58  return *this;
59  };
60 
61  cwoVect2<T>& operator/=(const T a){
62  x /= a; y /= a;
63  return *this;
64  };
65  const cwoVect2<T> operator+(const cwoVect2<T> &a) const{
66  cwoVect2<T> tmp;
67  tmp.x = x + a.x;
68  tmp.y = y + a.y;
69  return tmp;
70  }
71  const cwoVect2<T> operator+(const T a) const{
72  cwoVect2<T> tmp;
73  tmp.x = x + a;
74  tmp.y = y;
75  return tmp;
76  }
77  const cwoVect2<T> operator-(const cwoVect2<T> &a) const{
78  cwoVect2<T> tmp;
79  tmp.x = x - a.x;
80  tmp.y = y - a.y;
81  return tmp;
82  }
83  const cwoVect2<T> operator-(const T a) const{
84  cwoVect2<T> tmp;
85  tmp.x = x - a;
86  tmp.y = y;
87  return tmp;
88  }
89  const cwoVect2<T> operator*(const cwoVect2<T> &a) const{ //Hadamard product
90  cwoVect2<T> tmp;
91  tmp.x = x * a.x;
92  tmp.y = y * a.y;
93  return tmp;
94  }
95  const cwoVect2<T> operator*(const T a) const{
96  cwoVect2<T> tmp;
97  tmp.x = x * a;
98  tmp.y = y * a;
99  return tmp;
100  }
101  const cwoVect2<T> operator/(const cwoVect2<T> &a) const{ //Hadamard product
102  cwoVect2<T> tmp;
103  tmp.x = x / a.x;
104  tmp.y = y / a.y;
105  return tmp;
106  }
107  const cwoVect2<T> operator/(const T a) const{
108  cwoVect2<T> tmp;
109  tmp.x = x / a;
110  tmp.y = y / a;
111  return tmp;
112  }
113 
114  cwoVect2<T> operator()(const T t1, const T t2) const {
115  cwoVect2<T> tmp;
116  tmp.x = t1;
117  tmp.y = t2;
118  return tmp;
119  }
120 
121 
122 
123 };
124 
125 
127 using cwoint2 = cwoVect2<int>;
128 using cwofloat2 = cwoVect2<float>;
129 
130 
131 
133 class CWO;
134 
135 namespace cwo{
136 
137  void SetThreads(int N);
138  void Crop(CWO *src, CWO *dst, int x1, int y1, int Sx, int Sy);
139  CWO* Expand(CWO *src, CWO *dst);
140 
141  template<class T> T& Expand(T &src, T &dst)
142  {
143  long long int sNx = src.GetNx();
144  long long int sNy = src.GetNy();
145  long long int dNx = dst.GetNx();
146  long long int dNy = dst.GetNy();
147  cwoComplex *p_src = src.GetBuffer();
148  cwoComplex *p_dst = dst.GetBuffer();
149 
150  dst.SetCtx(src.GetCtx());
151  dst.SetSize(dNx, dNy);
152 
153  src.__Expand(
154  p_src, 0, 0, sNx, sNy,
155  p_dst, abs((dNx - sNx) / 2), abs((dNy - sNy) / 2), dNx, dNy,
156  CWO_C2C);
157 
158  /* if (dNx >= sNx && dNy >= sNy){
159  src.__Expand(
160  p_src, 0, 0, sNx, sNy,
161  p_dst, (dNx - sNx) / 2, (dNy - sNy) / 2, dNx, dNy,
162  CWO_C2C);
163  }
164  else{
165  src.__Expand(
166  p_src, (sNx - dNx) / 2, (sNy - dNy) / 2, sNx, sNy,
167  p_dst, 0, 0, dNx, dNy,
168  CWO_C2C);
169  }*/
170 
171  return dst;
172  }
173 
174  void Resize(CWO *src, CWO *dst, float scale_x, float scale_y, int flag);
175 
176 
177  CWO* DiffractAngular(CWO *src, CWO *dst, float z);
178  CWO* DiffractARSSFresnel(CWO *src, CWO *dst, float z);
179 
180  template<class T> T& DiffractAngular(T &src, T &dst, float z){
181  src.SetPropDist(z);
182  dst.SetCtx(src.GetCtx());
183  long long int Nx = src.GetNx();
184  long long int Ny = src.GetNy();
185  cwoComplex *p_src = src.GetBuffer();
186  cwoComplex *p_dst = dst.GetBuffer();
187  dst.__FFT(p_src, p_dst, CWO_C2C);
188  dst.__AngularProp(p_dst, CWO_PROP_MUL_FFT_SHIFT);
189  dst.__IFFT(p_dst, p_dst);
190  dst.__Mul(p_dst, cwo::MakeCplx(1.0f / (Nx*Ny), 0), p_dst);
191  return dst;
192  }
193 
194  template<class T> T& DiffractARSSFresnel(T &src, T &dst, float z){
195  src.SetPropDist(z);
196  dst.SetCtx(src.GetCtx());
197 
198  long long int Nx = src.GetNx();
199  long long int Ny = src.GetNx();
200  cwoComplex *p_src = src.GetBuffer();
201  cwoComplex *p_dst = dst.GetBuffer();
202 
203  double wl = src.GetWaveLength();
204  double spx = src.GetSrcPx();
205  double spy = src.GetSrcPy();
206  double dpx = src.GetDstPx();
207  double dpy = src.GetDstPy();
208  double ox = src.GetSrcOx() + src.GetDstOx();
209  double oy = src.GetSrcOy() + src.GetDstOy();
210  double sx1 = spx / dpx;
211  double sy1 = spy / dpy;
212  double sx2 = dpx / spx;
213  double sy2 = dpy / spy;
214 
215  int cz1x = (int)fabs((wl*fabs(z) / dpx - fabs(2 * sx1*ox)) / (2 * fabs(1 - sx1)*dpx));
216  int cz1y = (int)fabs((wl*fabs(z) / dpy - fabs(2 * sy1*oy)) / (2 * fabs(1 - sy1)*dpy));
217  int u1x = (int)fabs((wl*fabs(z) / dpx - fabs(2 * sx1*ox)) / (2 * (sx1*sx1 - sx1)*dpx));
218  int u1y = (int)fabs((wl*fabs(z) / dpy - fabs(2 * sy1*oy)) / (2 * (sy1*sy1 - sy1)*dpy));
219  int h1x = (int)(wl*fabs(z) / (2 * sx1*dpx*dpx));
220  int h1y = (int)(wl*fabs(z) / (2 * sy1*dpy*dpy));
221 
222  int cz2x = (int)fabs((wl*fabs(z) / spx - fabs(2 * sx2*ox)) / (2 * (sx2*sx2 - sx2)*spx));
223  int cz2y = (int)fabs((wl*fabs(z) / spy - fabs(2 * sy2*oy)) / (2 * (sy2*sy2 - sy2)*spy));
224  int u2x = (int)fabs((wl*fabs(z) / spx - fabs(2 * sx2*ox)) / (2 * fabs(1 - sx2)*spx));
225  int u2y = (int)fabs((wl*fabs(z) / spy - fabs(2 * sy2*oy)) / (2 * fabs(1 - sy2)*spy));
226  int h2x = (int)(wl*fabs(z) / (2 * sx2*spx*spx));
227  int h2y = (int)(wl*fabs(z) / (2 * sy2*spy*spy));
228 
229  src.__ARSSFresnelAperture(p_src);
230  src.__FFT(p_src, p_src, CWO_C2C);
231 
232  src.__ARSSFresnelProp(p_dst);
233  src.__FFTShift(p_dst);
234  src.__RectFillOutside(
235  p_dst, Nx / 2 - h1x, Ny / 2 - h1y,
236  h1x * 2, h1y * 2, cwo::MakeCplx(0, 0));
237  src.__FFTShift(p_dst);
238  src.__FFT(p_dst, p_dst, CWO_C2C);
239 
240  src.__Mul(p_src, p_dst, p_dst);
241  src.__IFFT(p_dst, p_dst);
242  src.__ARSSFresnelCoeff(p_dst);//include the multiplication of 1/(Nx*Ny)
243 
244  //clear padding area for next propagation
245  dst.Rect(dst.GetNx() / 4, dst.GetNy() / 4, dst.GetNx() / 2, dst.GetNy() / 2, cwo::MakeCplx(0, 0), CWO_FILL_OUTSIDE);
246 
247  return dst;
248  }
249 
250 
251  template<class T> T& Diffract(T &src, T &dst, float z, int type = CWO_ANGULAR){
252  switch (type){
253  case CWO_ANGULAR: DiffractAngular(src, dst, z); break;
254  case CWO_ARSS_FRESNEL: DiffractARSSFresnel(src, dst, z); break;
255  }
256  return dst;
257  }
258 
259 
260 
263  //@param affine_mat : pointer to affine matrix (3x3 = 9 float components)
264  //@param deg : rotating angle (in degree)
265  //*/
266  //template<typename T>
267  //void Save(const std::string fname, const T &fld,int div_x, int div_y);
268 
269 
270 
271  std::string ExtractFileName(const std::string &path);
272  std::string ExtractPath(const std::string &path);
273  std::string ExtractExt(const std::string &path);
274 
275 
276 
277 
278 
280 
283  // void DiffractDiv(CWO *u, cwoFloat2 )
284 
285 
286 
287  //int GetAliasingFreeApproxSphWave(float z, float p, float wl){
288  // int r = (int)((wl*fabs(z) / (2.0*p)) / p);
289  // return r;
290  //}
291 
292 
294 
298  void AffineRotX(float *affine_mat, float deg);
300 
304  void AffineRotY(float *affine_mat, float deg);
306 
310  void AffineRotZ(float *affine_mat, float deg);
311 
312 
313  inline cwoComplex MakeCplx(float re, float im){ cwoComplex tmp; CWO_RE(tmp) = re; CWO_IM(tmp) = im; return tmp; };
314  inline cwoFloat2 MakeFloat2(float x, float y){ cwoFloat2 tmp; tmp.x = x; tmp.y = y; return tmp; };
315  inline cwoFloat3 MakeFloat3(float x, float y, float z){ cwoFloat3 tmp; tmp.x = x; tmp.y = y; tmp.z = z; return tmp; };
316  template<class Type> inline cwoVect2<Type> MakeVect2(Type x, Type y) { cwoVect2<Type> tmp; tmp.x = x; tmp.y = y; return tmp; };
317 
318 
319 
321 
324  float GetTime();
325 
326 
327 }
328 
329 
330 
331 
332 
334 
337 class CWO {
338 public:
339  cwoCtx ctx;//current context
340  //cwoCtx prev_ctx; //previous context
341 
342 //protected:
343 public:
344  long long int rnd_seed;
345 
346 public:
347  //****************
348  //For Diffraction
349  //****************
350  void *p_field;//complex amplitude field
351 
352  //****************
353  //other
354  //****************
355  void InitParams();
356  void InitBuffers();
357  void FreeBuffers();
358 
359 
360 public:
361  virtual void* __Malloc(size_t size);
362  virtual void __Free(void **a);
363  virtual void __Memcpy(void *dst, void *src, size_t size);
364  virtual void __Memset(void *p, int c, size_t size);
365  //virtual void __Memset(void *p, cwoComplex c, unsigned int size);
366 
367  virtual void __Expand(
368  void *src, int sx, int sy, int srcNx, int srcNy,
369  void *dst, int dx, int dy, int dstNx, int dstNy,
370  int type);
371 
372  virtual void __ExpandMirror(
373  cwoComplex* src, int srcNx, int srcNy,
374  cwoComplex* dst, int dstNx, int dstNy,
375  int type);
376 
377  virtual void __ShiftedFresnelAperture(cwoComplex *a);
378  virtual void __ShiftedFresnelProp(cwoComplex *a);
379  virtual void __ShiftedFresnelCoeff(cwoComplex *a);
380 
381  virtual void __ARSSFresnelAperture(cwoComplex *a);
382  virtual void __ARSSFresnelProp(cwoComplex *a);
383  virtual void __ARSSFresnelCoeff(cwoComplex *a);
384 
385  virtual void __FresnelConvProp(cwoComplex *a);
386  virtual void __FresnelConvCoeff(cwoComplex *a, float const_val=1.0f);
387 
388  virtual void __FresnelAnalysisTransfer(cwoComplex *a, cwoComplex *b);
389 
390  virtual void __AngularProp(cwoComplex *a, int flag);
391  void __AngularLim(float *fx_c, float *fx_w, float *fy_c, float *fy_w);
392 // virtual void __ShiftedAngularProp(cwoComplex *a);
393 
394  virtual void __HuygensProp(cwoComplex *a);
395 
396  virtual void __FresnelFourierProp(cwoComplex *a);
397  virtual void __FresnelFourierCoeff(cwoComplex *a);
398 
399  virtual void __FresnelDblAperture(cwoComplex *a, float z1);
400  virtual void __FresnelDblFourierDomain(cwoComplex *a, float z1, float z2, cwoInt4 *zp);
401  virtual void __FresnelDblCoeff(cwoComplex *a, float z1, float z2);
402 
403  virtual void __FFT(void *src, void *dst, int type);
404  virtual void __IFFT(void *src, void *dst);
405  virtual void __FFTShift(void *src);
406 
407  void __Add(cwoComplex *a, float b, cwoComplex *c);//c=a+b
408  virtual void __Add( cwoComplex *a, const cwoComplex &b, cwoComplex *c);//c=a+b
409  virtual void __Add( cwoComplex *a, cwoComplex *b, cwoComplex *c);//c=a+b
410 
411  void __Sub( cwoComplex *a, float b, cwoComplex *c);//c=a-b
412  virtual void __Sub( cwoComplex *a, const cwoComplex &b, cwoComplex *c);//c=a-b
413  virtual void __Sub( cwoComplex *a, cwoComplex *b, cwoComplex *c);//c=a-b
414 
415  void __Mul( cwoComplex *a, float b, cwoComplex *c);//c=a*b
416  virtual void __Mul( cwoComplex *a, const cwoComplex &b, cwoComplex *c);//c=a*b
417  virtual void __Mul( cwoComplex *a, cwoComplex *b, cwoComplex *c);//c=a*b
418 
419  void __Div( cwoComplex *a, float b, cwoComplex *c);//c=a/b
420  virtual void __Div( cwoComplex *a, const cwoComplex &b, cwoComplex *c);//c=a/b
421  virtual void __Div( cwoComplex *a, cwoComplex *b, cwoComplex *c);//c=a/b
422 
423  void __Rep(cwoComplex *a, float b, cwoComplex *c); //Reciprocal number c=b/a
424  virtual void __Rep(cwoComplex *a, const cwoComplex &b, cwoComplex *c);//Reciprocal number c=b/a
425  virtual void __Rep(cwoComplex *a, cwoComplex *b, cwoComplex *c);//Reciprocal number c=b/a
426 
427  virtual void __Re(cwoComplex *a , cwoComplex *b);
428  virtual void __Im(cwoComplex *a , cwoComplex *b);
429  virtual void __Conj(cwoComplex *a);
430 
431 
432  virtual void __Intensity(cwoComplex *a, cwoComplex *b); //OK
433  virtual void __Amp(cwoComplex *a , cwoComplex *b);//OK
434  virtual void __Phase(cwoComplex *a , cwoComplex *b, float offset);//OK
435  virtual void __Arg(cwoComplex *a , cwoComplex *b, float scale, float offset);//OK
436  virtual void __Real2Complex(float *src, cwoComplex *dst);
437  virtual void __Phase2Complex(float *src, cwoComplex *dst);
438  virtual void __Arg2Cplx(cwoComplex *src, cwoComplex *dst, float scale, float offset);
439  virtual void __Polar(float *amp, float *ph, cwoComplex *c);
440  virtual void __ReIm(cwoComplex *re, cwoComplex *im, cwoComplex *c);
441  virtual void __RectFillInside(cwoComplex *p, int m, int n, int Sx, int Sy, cwoComplex a);
442  virtual void __RectFillOutside(cwoComplex *p, int m, int n, int Sx, int Sy, cwoComplex a);
443 
444  //virtual void __CopyFloat(float *src, int x1, int y1, int sNx, int sNy, float *dst, int x2, int y2, int dNx, int dNy, int Sx, int Sy);
445  //virtual void __CopyComplex(cwoComplex *src, int x1, int y1, int sNx, int sNy, cwoComplex *dst, int x2, int y2, int dNx, int dNy, int Sx, int Sy);
446 
447  virtual void __FloatToChar(char *dst, float *src, int N);
448  virtual void __CharToFloat(float *dst, char *src, int N);
449 
450  void FromChar(const char *src, int Nx, int Ny);
451  void FromFloat(const float *src, int Nx, int Ny);
452 
453  //void SetFieldType(int type);
454  void SetSize(int Nx, int Ny);
455  void SetPropDist(float z);
456  float GetPropDist();
457 
460  //@return if true, the calculation size was changed, otherwise false
461  //*/
463 
464 public:
465 
466  //###########################################
470  CWO();
471 
473 
479  virtual ~CWO();
480 
481  CWO(CWO &tmp);
482  CWO(int Nx, int Ny);
483 
485 
486  //###########################################
490  /*CWO& operator=(CWO& const tmp);
491  CWO& operator=(float a);
492  CWO& operator=(const cwoComplex &a);*/
493  CWO& operator=(CWO& tmp);
494  CWO& operator=(float a);
495  CWO& operator=(const cwoComplex &a);
496 
497  //CWO(CWO &&a) {
498  // //printf("Move Constructor \n");
499  // if (this != &a) {
500  // p_field = a.p_field;
501  // a.p_field = nullptr;
502  // }
503  // ctx = a.ctx;
504  //
505  //};
506  //CWO& operator=(CWO&& a) {
507  // if (this != &a) {
508  // delete p_field;
509  // p_field = a.p_field;
510  // a.p_field = nullptr;
511  // }
512  // ctx = a.ctx;
513  // printf("move assignment\n");
514  // return *this;
515  //}
516 
517 //#define __NOW_DEBUGGIBNG_OPERATOR
518 #ifdef __NOW_DEBUGGIBNG_OPERATOR
519  template<class T> friend CWO operator+(CWO& const lhs, T& const rhs){ printf("ope1 + \n"); CWO nrv = lhs; nrv += rhs; return nrv; };
520  template<class T> friend CWO&& operator+(CWO&& lhs, T& const rhs){ printf("ope2 + \n"); lhs += rhs; return std::move(lhs); };
521  template<class T> friend CWO&& operator+(const CWO& lhs, T&& rhs){ printf("ope3 + \n"); rhs += lhs; return std::move(rhs); };
522  template<class T> friend CWO&& operator+(CWO&& lhs, T&& rhs){ printf("ope4 + \n"); lhs += std::move(rhs); return std::move(lhs); };
523 
524  template<class T> friend CWO operator-(CWO& const lhs, T& const rhs){ printf("ope1 - \n"); CWO nrv = lhs; nrv -= rhs; return nrv; };
525  template<class T> friend CWO&& operator-(CWO&& lhs, T& const rhs){ printf("ope2 - \n"); lhs -= rhs; return std::move(lhs); };
526  template<class T> friend CWO&& operator-(const T& lhs, CWO&& rhs){ printf("ope3 - \n"); rhs -= lhs; return std::move(rhs); };
527  template<class T> friend CWO&& operator-(CWO&& lhs, T&& rhs){ printf("ope4 - \n"); lhs -= std::move(rhs); return std::move(lhs); };
528 
529  template<class T> friend CWO operator*(CWO& const lhs, T& const rhs){ printf("ope1 * \n"); CWO nrv = lhs; nrv *= rhs; return nrv; };
530  template<class T> friend CWO&& operator*(CWO&& lhs, T& const rhs){ printf("ope2 * \n"); lhs *= rhs; return std::move(lhs); };
531  template<class T> friend CWO&& operator*(const T& lhs, CWO&& rhs){ printf("ope3 * \n"); rhs *= lhs; return std::move(rhs); };
532  template<class T> friend CWO&& operator*(CWO&& lhs, T&& rhs){ printf("ope4 * \n"); lhs *= std::move(rhs); return std::move(lhs); };
533 
534  template<class T> friend CWO operator/(CWO& const lhs, T& const rhs){ printf("ope1 / \n"); CWO nrv = lhs; nrv /= rhs; return nrv; };
535  template<class T> friend CWO&& operator/(CWO&& lhs, T& const rhs){ printf("ope2 / \n"); lhs /= rhs; return std::move(lhs); };
536  template<class T> friend CWO&& operator/(const T& lhs, CWO&& rhs){ printf("ope3 / \n"); rhs /= lhs; return std::move(rhs); };
537  template<class T> friend CWO&& operator/(CWO&& lhs, T&& rhs){ printf("ope4 / \n"); lhs /= std::move(rhs); return std::move(lhs); };
538 #else
539  template<class T=CWO> T operator+(T &a) {
540  T tmp;
541  tmp.Create(GetNx(), GetNy());
542  assert(tmp.GetBuffer() != nullptr);
543  (tmp).ctx = ctx;
544  __Add(GetBuffer(), a.GetBuffer(), tmp.GetBuffer());
545  return tmp;
546  }
547  template<class T=CWO> T operator+(float a) {
548  T tmp;
549  tmp.Create(GetNx(), GetNy());
550  assert(tmp.GetBuffer() != nullptr);
551  (tmp).ctx = ctx;
552  cwoComplex cplx;
553  CWO_RE(cplx) = a; CWO_IM(cplx) = 0.0f;
554  __Add(GetBuffer(), cplx, tmp.GetBuffer());
555  return tmp;
556  }
557  template<class T=CWO> T operator+(const cwoComplex &a) {
558  T tmp;
559  tmp.Create(GetNx(), GetNy());
560  assert(tmp.GetBuffer() != nullptr);
561  (tmp).ctx = ctx;
562  __Add(GetBuffer(), a, tmp.GetBuffer());
563  return (tmp);
564  }
565  template<class T> friend CWO operator+(const T &a, CWO& c) { CWO tmp = c; tmp += a; return tmp; };
566 
567  template<class T = CWO> T operator-(T &a) {
568  T tmp;
569  tmp.Create(GetNx(), GetNy());
570  assert(tmp.GetBuffer() != nullptr);
571  (tmp).ctx = ctx;
572  __Sub(GetBuffer(), a.GetBuffer(), tmp.GetBuffer());
573  return tmp;
574  }
575  template<class T = CWO> T operator-(float a) {
576  T tmp;
577  tmp.Create(GetNx(), GetNy());
578  assert(tmp.GetBuffer() != nullptr);
579  (tmp).ctx = ctx;
580  cwoComplex cplx;
581  CWO_RE(cplx) = a; CWO_IM(cplx) = 0.0f;
582  __Sub(GetBuffer(), cplx, tmp.GetBuffer());
583  return tmp;
584  }
585  template<class T = CWO> T operator-(const cwoComplex &a) {
586  T tmp;
587  tmp.Create(GetNx(), GetNy());
588  assert(tmp.GetBuffer() != nullptr);
589  (tmp).ctx = ctx;
590  __Sub(GetBuffer(), a, tmp.GetBuffer());
591  return (tmp);
592  }
593  template<class T> friend CWO operator-(const T &a, CWO& c) { CWO tmp(c.GetNx(),c.GetNy()); tmp = a; tmp -= c; return tmp; };
594 
595  template<class T = CWO> T operator*(T &a) {
596  T tmp;
597  tmp.Create(GetNx(), GetNy());
598  assert(tmp.GetBuffer() != nullptr);
599  (tmp).ctx = ctx;
600  __Mul(GetBuffer(), a.GetBuffer(), tmp.GetBuffer());
601  return tmp;
602  }
603  template<class T = CWO> T operator*(float a) {
604  T tmp;
605  tmp.Create(GetNx(), GetNy());
606  assert(tmp.GetBuffer() != nullptr);
607  (tmp).ctx = ctx;
608  cwoComplex cplx;
609  CWO_RE(cplx) = a; CWO_IM(cplx) = 0.0f;
610  __Mul(GetBuffer(), cplx, tmp.GetBuffer());
611  return tmp;
612  }
613  template<class T = CWO> T operator*(const cwoComplex &a) {
614  T tmp;
615  tmp.Create(GetNx(), GetNy());
616  assert(tmp.GetBuffer() != nullptr);
617  (tmp).ctx = ctx;
618  __Mul(GetBuffer(), a, tmp.GetBuffer());
619  return (tmp);
620  }
621  template<class T> friend CWO operator*(const T &a, CWO& c) { CWO tmp = c; tmp *= a; return tmp; };
622 
623  template<class T = CWO> T operator/(T &a) {
624  T tmp;
625  tmp.Create(GetNx(), GetNy());
626  assert(tmp.GetBuffer() != nullptr);
627  (tmp).ctx = ctx;
628  __Div(GetBuffer(), a.GetBuffer(), tmp.GetBuffer());
629  return tmp;
630  }
631  template<class T = CWO> T operator/(float a) {
632  T tmp;
633  tmp.Create(GetNx(), GetNy());
634  assert(tmp.GetBuffer() != nullptr);
635  (tmp).ctx = ctx;
636  cwoComplex cplx;
637  CWO_RE(cplx) = a; CWO_IM(cplx) = 0.0f;
638  __Div(GetBuffer(), cplx, tmp.GetBuffer());
639  return tmp;
640  }
641  template<class T = CWO> T operator/(const cwoComplex &a) {
642  T tmp;
643  tmp.Create(GetNx(), GetNy());
644  assert(tmp.GetBuffer() != nullptr);
645  (tmp).ctx = ctx;
646  __Div(GetBuffer(), a, tmp.GetBuffer());
647  return (tmp);
648  }
649  template<class T> friend CWO operator/(const T &a, CWO& c) {
650  CWO tmp;
651  tmp.Create(c.GetNx(), c.GetNy());
652  assert(tmp.GetBuffer() != nullptr);
653  (tmp).ctx = c.ctx;
654  c.__Rep(c.GetBuffer(), a, tmp.GetBuffer());
655  return tmp;
656  };
657 
658  //CWO operator-(CWO &a);
659  //CWO operator-(float a);
660  //CWO operator-(cwoComplex a);
661  //template<class T> friend CWO operator-(const T a, CWO& const c) { CWO tmp = c; tmp -= a; return tmp; };
662 
663  //CWO operator*(CWO &a);
664  //CWO operator*(float a);
665  //CWO operator*(cwoComplex a);
666  //template<class T> friend CWO operator*(const T a, CWO& const c) { CWO tmp = c; tmp *= a; return tmp; };
667 #endif
668 
669  //CWO operator/(CWO &a);
670  //CWO operator/(float a);
671  //CWO operator/(cwoComplex a);
672 
673 
674  CWO& operator+=(CWO& const a);
675  CWO& operator+=(float a);
676  CWO& operator+=(cwoComplex a);
677  CWO& operator-=(CWO &a);
678  CWO& operator-=(float a);
679  CWO& operator-=(cwoComplex a);
680  CWO& operator*=(CWO &a);
681  CWO& operator*=(float a);
682  CWO& operator*=(cwoComplex a);
683  CWO& operator/=(CWO &a);
684  CWO& operator/=(float a);
685  CWO& operator/=(cwoComplex a);
686 
687 
688 
690 
691  //###########################################
695 
696 
698 
701  virtual void Send(CWO &a);
703 
707  virtual void Recv(CWO &a);
708 
710 
712 
716  int Create(int Nx, int Ny=1);
717 
719 
721  virtual void Destroy();
722 
723 
724 // void Attach(CWO &a);
725 
726  void Attach(cwoCtx *ctx, cwoComplex *p) {
727  __Free(&p_field);
728  SetCtx(*ctx);
729  p_field = p;
730  }
731 
732  //###########################################
736 
738 
740  void SetTimer();
741 
743 
746  double EndTimer();
747 
749 
750 
751 
752  //###########################################
756 
758 
762  void SetCtx(const cwoCtx& c){
763  ctx = c;
764  };
765 
767 
771  cwoCtx& GetCtx(){return ctx;};
772 
774 
777  void SetCalcType(int type){
778  ctx.calc_type = type;
779  }
780 
782 
786  void SetPitch(float px, float py){
787  SetSrcPitch(px, py);
788  SetDstPitch(px, py);
789  }
790 
791  void SetPitch(float p){
792  SetSrcPitch(p, p);
793  SetDstPitch(p, p);
794  }
795 
796  void SetPitch(const cwoFloat2 &p){
797  SetPitch(p.x, p.y);
798  }
799 
800  void SetPitch(const cwoVect2<float> &p){
801  SetPitch(p.x, p.y);
802  }
803 
805 
809  void SetSrcPitch(float px, float py){
810  ctx.src_px = px;
811  ctx.src_py = py;
812  }
813  void SetSrcPitch(float p){
814  SetSrcPitch(p, p);
815  }
816 
817  void SetSrcPitch(const cwoFloat2 &p){
818  SetSrcPitch(p.x, p.y);
819  }
820 
821  void SetSrcPitch(const cwoVect2<float> &p){
822  SetSrcPitch(p.x, p.y);
823  }
824 
826 
830  void SetDstPitch(float px, float py){
831  ctx.dst_px = px;
832  ctx.dst_py = py;
833  }
834 
835  void SetDstPitch(float p){
836  SetDstPitch(p, p);
837  }
838 
839  void SetDstPitch(const cwoFloat2 &p){
840  SetDstPitch(p.x, p.y);
841  }
842 
843  void SetDstPitch(const cwoVect2<float> &p){
844  SetDstPitch(p.x, p.y);
845  }
846 
848 
851  void SetWaveLength(float w){
852  ctx.wave_length = w;
853  }
854 
856 
860  void SetOffset(float x, float y){
861  SetSrcOffset(x, y);
862  SetDstOffset(0.0f, 0.0f);
863  }
864 
865  void SetOffset(const cwoFloat2 &p){
866  SetOffset(p.x, p.y);
867  }
868 
869  void SetOffset(const cwoVect2<float> &p){
870  SetOffset(p.x, p.y);
871  }
872 
874 
879  void SetSrcOffset(float x, float y){
880  ctx.src_ox = x;
881  ctx.src_oy = y;
882  }
883 
884  void SetSrcOffset(const cwoFloat2 &p){
885  SetSrcOffset(p.x, p.y);
886  }
887 
888  void SetSrcOffset(const cwoVect2<float> &p){
889  SetSrcOffset(p.x, p.y);
890  }
891 
892 
893 
895 
900  void SetDstOffset(float x, float y){
901  ctx.dst_ox = x;
902  ctx.dst_oy = y;
903  }
904 
905  void SetDstOffset(const cwoFloat2 &p){
906  SetDstOffset(p.x, p.y);
907  }
908 
909  void SetDstOffset(const cwoVect2<float> &p){
910  SetDstOffset(p.x, p.y);
911  }
912 
914 
917  //void SetThreads(int Nx, int Ny=1);
918  void SetThreads(int Nx, int Ny=1);
919 
920 
922 
925  int GetThreads();
926 
928 
931  size_t GetNx(){
932  return ctx.Nx;
933  }
935 
938  size_t GetNy(){
939  return ctx.Ny;
940  }
942 
945  cwoVect2<size_t> GetShape(){
946  return cwo::MakeVect2(GetNx(), GetNy());
947  }
948 
949 
951 
954  float GetWaveLength(){
955  return ctx.wave_length;
956  }
957 
959 
962  float GetWaveNum(){
963  return 2.0f*CWO_PI / GetWaveLength();
964  }
966 
969  float GetDistance(){
970  return ctx.z;
971  }
973 
976  float GetPx(){
977  return ctx.src_px;
978  }
980 
983  float GetPy(){
984  return ctx.src_py;
985  }
986 
987 
988  //cwoFloat2 GetPitch(){
989  // return cwo::MakeFloat2(GetPx(), GetPy());
990  //}
991  //
992  cwoVect2<float> GetPitch(){
993  //cwoVect2<float> ret(GetPx(), GetPy())
994  return cwo::MakeVect2(GetPx(), GetPy());
995  }
996 
998 
1001  float GetSrcPx(){
1002  return GetPx();
1003  }
1005 
1008  float GetSrcPy(){
1009  return GetPy();
1010  }
1011 
1012  cwoVect2<float> GetSrcPitch(){
1013  return cwo::MakeVect2(GetSrcPx(), GetSrcPy());
1014  }
1015 
1016 
1017 
1019 
1022  float GetDstPx(){
1023  return ctx.dst_px;
1024  }
1026 
1029  float GetDstPy(){
1030  return ctx.dst_py;
1031  }
1032 
1033  cwoVect2<float> GettDstPitch(){
1034  return cwo::MakeVect2(GetDstPx(), GetDstPy());
1035  }
1036 
1038 
1041  float GetOx(){
1042  return GetSrcOx();
1043  }
1045 
1048  float GetOy(){
1049  return GetSrcOy();
1050  }
1051 
1052  cwoVect2<float> GetOffset(){
1053  return cwo::MakeVect2(GetOx(), GetOy());
1054  }
1055 
1057 
1060  float GetSrcOx(){
1061  return ctx.src_ox;
1062  }
1064 
1067  float GetSrcOy(){
1068  return ctx.src_oy;
1069  }
1070 
1071  cwoVect2<float> GetSrcOffset(){
1072  return cwo::MakeVect2(GetSrcOx(), GetSrcOy());
1073  }
1074 
1076 
1079  float GetDstOx(){
1080  return ctx.dst_ox;
1081  }
1083 
1086  float GetDstOy(){
1087  return ctx.dst_oy;
1088  }
1089 
1090  cwoVect2<float> GetDstOffset(){
1091  return cwo::MakeVect2(GetDstOx(), GetDstOy());
1092  }
1093 
1094 
1096 
1099  float GetLx(){ return GetSrcLx(); }
1100 
1102 
1105  float GetLy(){ return GetSrcLy(); }
1106 
1107  cwoVect2<float> GetLength(){
1108  return cwo::MakeVect2(GetLx(), GetLy());
1109  }
1110 
1112 
1115  float GetSrcLx(){
1116  return GetSrcPx()*GetNx();
1117  }
1118 
1120 
1123  float GetSrcLy(){
1124  return GetSrcPy()*GetNy();
1125  }
1126 
1127  cwoVect2<float> GetSrcLength(){
1128  return cwo::MakeVect2(GetSrcLx(), GetSrcLy());
1129  }
1130 
1132 
1135  float GetDstLx(){
1136  return GetDstPx()*GetNx();
1137  }
1138 
1140 
1143  float GetDstLy(){
1144  return GetDstPy()*GetNy();
1145  }
1146 
1147  cwoVect2<float> GetDstLength(){
1148  return cwo::MakeVect2(GetDstLx(), GetDstLy());
1149  }
1150 
1151  cwoVect2<int> GetCenterPix() {
1152  int x = GetNx() / 2;
1153  int y = GetNx() / 2;
1154  return cwo::MakeVect2(x,y);
1155  }
1156  //cwoVect2<float> GetCenter() {
1157  // int x = GetNx() / 2;
1158  // int y = GetNx() / 2;
1159  // return cwo::MakeVect2(x, y);
1160  //}
1161 
1162 
1163 
1165 
1168  int GetAliasingFreeApproxSphWave(float z, float p, float wl){
1169  int r = (int)(fabs(z)*asin(wl / (2.0*p)) / p);
1170  return r;
1171  }
1172 
1173 
1174 
1176  // /*!
1177  // @return Current field-type. (e.g. CWO_FLD_COMPLEX)
1178  // */
1179  //int GetFieldType();
1181  // /*!
1182  // @return Current field-type as string. (e.g. string "CWO_FLD_COMPLEX")
1183  // */
1184  //char *GetFieldName();
1185 
1187  void ShowParam();
1188 
1189 
1190 
1192 
1193  //###########################################
1197  int CheckExt(const char* fname, const char* ext);
1198 
1200 
1206  int Load(const std::string& fname, int c = CWO_GREY);
1207 
1209 
1216  int Load(const std::string& fname_amp, const std::string& fname_pha, int c = CWO_GREY);
1217 
1218 
1220 
1228  int Save(const std::string& fname, int bmp_8_24 = 24);
1229  int Save(const std::string& fname, CWO *r, CWO *g = nullptr, CWO *b = nullptr);
1230 
1231  int SaveMonosToColor(const std::string& fname, const std::string& r_name, const std::string& g_name, const std::string& b_name);
1232  int SaveAsImage(const std::string& fname, float i1, float i2, float o1, float o2, int flag = CWO_SAVE_AS_RE, int bmp_8_24 = 24);
1233  int SaveAsImage(const std::string& fname, int flag = CWO_SAVE_AS_RE, int bmp_8_24 = 24);
1234  virtual int SaveAsImage(cwoComplex *p, int Nx, int Ny, const std::string& fname, int flag = CWO_SAVE_AS_RE, int bmp_8_24 = 24);
1235 
1236 
1237  int SaveLineAsText(const std::string& fname, int flag);
1238  int SaveLineAsText(const std::string& fname, int flag, int x1, int y1, int x2, int y2);
1239 
1241 
1242 
1243  //###########################################
1247 
1249 
1252  cwoComplex* GetBuffer(){
1253  return (cwoComplex *)p_field;
1254  }
1255 
1256  int CmpCtx(cwoCtx &a, cwoCtx &b);
1257 
1259 
1262  size_t GetMemSizeCplx();
1264 
1267  size_t GetMemSizeFloat();
1269 
1272  size_t GetMemSizeChar();
1273 
1275 
1279  int GetDfrExpand();
1280  void SetDfrExpand(int flag = 0);
1281 
1282 
1283  void Diffract(float d, int type = CWO_ANGULAR, cwoInt4 *zp = nullptr, cwoComplex *knl_mask = nullptr, CWO *numap = nullptr);
1284  //void Diffract(float d, int type, CWO *tmp_buf1,CWO *tmp_buf2);
1285  void Diffract(float d, int type, CWO *tmp_buf);
1286  void Diffract(int type, CWO *kernel[]);
1287 
1288  void DiffractConv(float d, int type, cwoComplex *ape, cwoComplex *prop, cwoInt4 *zp = nullptr, cwoComplex *knl_mask = nullptr, CWO *numap = nullptr, int prop_calc = 1);
1289  //void DiffractConvImplicit(float d, int type, cwoComplex *ape, cwoComplex *prop, int prop_calc=1);
1290  void DiffractFourier(float d, int type, cwoInt4 *zp=NULL);
1291  void DiffractDirect(float d, CWO *snumap, CWO *dnumap);
1292 
1294 
1299  void DiffractAffine(float *affine_mat, int interpolation = CWO_INTER_LINEAR);
1300 
1301 
1302 
1303 
1304  //void Diffract3D(CWO &a, float d, int type);
1305  //void Diffract(float d, int type, int Dx, int Dy, char *dir = nullptr);
1306 
1307  void DiffractDepth(float z, float dz, int num,
1308  const char* path, const char *name, const char *ext, int type=CWO_ANGULAR);
1309 
1310 
1311  void AngularProp(CWO *a);
1312  void ShiftedAngularProp(CWO *a);
1313 
1314  void FresnelConvProp(CWO *a);
1315  void FresnelConvCoeff(CWO *a);
1316 
1317  void ShiftedFresnelProp(CWO *a);
1318  void ShiftedFresnelAperture(CWO *a);
1319  void ShiftedFresnelCoeff(CWO *a);
1320 
1321  void ARSSFresnelProp(CWO *a);
1322  void ARSSFresnelAperture(CWO *a);
1323  void ARSSFresnelCoeff(CWO *a);
1324 
1325  void FresnelFourierProp(CWO *a);
1326  void FresnelFourierCoeff(CWO *a);
1327 
1328  void FresnelDblAperture(CWO *a,float z1);
1329  void FresnelDblFourierDomain(CWO *a, float z1, float z2, cwoInt4 *zp);
1330  void FresnelDblCoeff(CWO *a, float z1, float z2);
1331 
1332  void FresnelDupProp(CWO *a, float M);
1333 
1334  //CWO* PreDiffract(float z, int type, CWO *a);
1335  void PreDiffract(float z, int type, CWO *a[]);
1336 
1337  void ReleaseTmpBuffer();
1338 
1339  void FresnelCoeff(float z1, float z2, float wl1, float wl2);
1340 
1341  //void AngularProp(float z, int iNx, int iNy);
1342 
1343 
1344  float fresnel_c(float x);
1345  float fresnel_s(float x);
1346  void FresnelInt(
1347  float z, int x1, int y1, int x2, int y2);
1349 
1350  //###########################################
1355 
1356  virtual void __AddSphericalWave(cwoComplex *p, float x, float y, float z, float px, float py, float a, float ph=0.0f);
1357  virtual void __MulSphericalWave(cwoComplex *p, float x, float y, float z, float px, float py, float a);
1358 
1361 
1368  void AddSphericalWave(float x, float y, float z, float a, float ph=0.0f);
1369 
1371 
1379  void AddSphericalWave(float x, float y, float z, float px, float py, float a, float ph=1.0f);
1381 
1391 
1392 
1399  void MulSphericalWave(float x, float y, float z, float a = 1.0f);
1400 
1401  void MulSphericalWave(float x, float y, float z, float px, float py, float a=1.0f);
1402 
1403 
1405 
1414 
1415 
1424  //virtual void __AddApproxSphWave(cwoComplex *p, float x, float y, float z, float px, float py, float a);
1425  virtual void __AddApproxSphWave(cwoComplex *p, float x, float y, float phi, float zx, float zy, float px, float py, float a);
1426 
1427  void AddApproxSphWave(float x, float y, float z, float a=1.0f);
1428  void AddApproxSphWave(float x, float y, float z, float px, float py, float a=1.0f);
1429  void AddApproxSphWave(float x, float y, float phi, float zx, float zy, float px, float py, float a=1.0f);
1430 
1432 
1440  //virtual void __MulApproxSphWave(cwoComplex *p, float x, float y, float phi, float px, float py, float a);
1441  virtual void __MulApproxSphWave(cwoComplex *p, float x, float y, float phi, float zx, float zy, float px, float py, float a);
1442 
1443  void MulApproxSphWave(float x, float y, float z, float a=1.0f);
1444  void MulApproxSphWave(float x, float y, float z, float px, float py, float a=1.0f);
1445  void MulApproxSphWave(float x, float y, float phi, float zx, float zy, float px, float py, float a=1.0f);
1446 
1447 
1449 
1456  void AddPlanarWave(float kx, float ky, float px, float py, float a=1.0f);
1458 
1463  void AddPlanarWave(float kx, float ky, float a=1.0f);
1464 
1466 
1475  void MulPlanarWave(float kx, float ky, float px, float py, float a=1.0f);
1477 
1483  void MulPlanarWave(float kx, float ky, float a=1.0f);
1484 
1486 
1487  //###########################################
1491 
1493 
1498  void Re();
1499 
1501 
1506  void Im();
1507 
1509 
1514  void Conj();
1515 
1517 
1522  void Intensity();
1523 
1525 
1530  void Amp();
1531 
1533 
1540  void Phase(float offset=0.0f);
1541 
1542 
1544 
1569  void Arg(float scale=1.0f, float offset=0.0f);
1570 
1572 // void Cplx();
1574 // int Cplx(CWO &amp, CWO &ph);
1575 
1576 
1577  void Expi();
1578 
1579 
1580 
1583 
1598  void Arg2Cplx(float scale=1.0, float offset=0.0);
1599 
1600 
1601  void Cart2Polar();
1602  void Polar2Cart();
1603 
1604 
1605 
1607 
1621  void ReIm(CWO &re, CWO &im);
1622 
1624 
1625 void Char(); //conver to current data type to char
1626 void Float(); //conver to current data type to float
1627 
1628 void Quant(float alpha);
1629 
1630 void ToDouble(double *a);
1631 void FromDouble(double *a);
1632 
1633 
1635 
1637  virtual void SqrtReal();
1639 
1641  virtual void SqrtCplx();
1642 
1643  void Div(float a);
1644 
1645 // cwoComplex Mul(cwoComplex a, cwoComplex b);
1646 // cwoComplex Conj(cwoComplex a);
1647 
1648  //
1649  cwoComplex Polar(float amp, float arg);
1650 
1651 
1652  //###########################################
1656 
1657  void SetPixel(int x, int y, float a);
1658  void SetPixel(const cwoVect2<int> &v, float a){ SetPixel(v.x, v.y, a); };
1659  void SetPixel(int x, int y, float amp, float ph);
1660  void SetPixel(const cwoVect2<int> &v, float amp, float ph){ SetPixel(v.x, v.y, amp, ph); };
1661  void SetPixel(int x, int y, cwoComplex &a);
1662  void SetPixel(const cwoVect2<int> &v, cwoComplex &a){ SetPixel(v.x, v.y, a); };
1663  void SetPixel(int x, int y, CWO &a);
1664  void SetPixel(const cwoVect2<int> &v, CWO &a){ SetPixel(v.x, v.y, a); };
1665 
1666  void AddPixel(int x, int y, cwoComplex a);
1667  void AddPixel(const cwoVect2<int> &v, cwoComplex a){ AddPixel(v.x, v.y, a); };
1668  void AddPixel(int x, int y, CWO &a);
1669  void AddPixel(const cwoVect2<int> &v, CWO &a){ AddPixel(v.x, v.y, a); };
1670 
1671  void MulPixel(int x, int y, float a);
1672  void MulPixel(const cwoVect2<int> &v, float a){ MulPixel(v.x, v.y, a); };
1673  void MulPixel(int x, int y, cwoComplex a);
1674  void MulPixel(const cwoVect2<int> &v, cwoComplex a){ MulPixel(v.x, v.y, a); };
1675 
1676  void GetPixel(int x, int y, float &a);
1677  void GetPixel(const cwoVect2<int> &v, float &a){ GetPixel(v.x, v.y, a); };
1678  void GetPixel(int x, int y, cwoComplex &a);
1679  void GetPixel(const cwoVect2<int> &v, cwoComplex &a){ GetPixel(v.x, v.y, a); };
1680  cwoComplex GetPixel(int m, int n);
1681  cwoComplex GetPixel(const cwoVect2<int> &v){ return GetPixel(v.x, v.y); };
1682  void GetPixel(int x, int y, CWO &a);
1683  void GetPixel(const cwoVect2<int> &v, CWO &a){ GetPixel(v.x, v.y, a); };
1684 
1685  void GetLineH(int n, int N, CWO &a);
1686  void GetLineV(int m, int N, CWO &a);
1687  void GetLine(int m1, int n1, int m2, int n2, CWO &a);
1688  void SetLineH(int n, CWO &a);
1689  void SetLineV(int m, CWO &a);
1690  void SetLine(int m1, int n1, int m2, int n2, CWO &a);
1691 
1692  virtual void __Copy(
1693  cwoComplex *src, int x1, int y1, int sNx, int sNy,
1694  cwoComplex *dst, int x2, int y2, int dNx, int dNy,
1695  int Sx, int Sy);
1696 
1697  void Copy(CWO &a, int x1, int y1, int x2,int y2, int Sx, int Sy);
1698  //CWO Extract(int x, int y, int Nx, int Ny);
1699 
1700  void Clear(int c=0);
1701  virtual void Fill(cwoComplex a);
1702 
1703  template <class T> void FlipH(T *a);
1704  template <class T> void FlipV(T *a);
1705  void Flip(int mode=0);
1706 
1707  void Crop(int x1, int y1, int Sx, int Sy);
1708  void Crop(const cwoVect2<int> &v, const cwoVect2<int> &s){ Crop(v.x, v.y, s.x, s.y); };
1709  void Crop(int x1, int y1, int Sx, int Sy, CWO &a);
1710  void Crop(const cwoVect2<int> &v, const cwoVect2<int> &s, CWO &a){ Crop(v.x, v.y, s.x, s.y, a); };
1711 
1712  void ShiftX(int s, int flag=0);
1713  void ShiftY(int s, int flag=0);
1714 
1715  void Transpose();
1716  virtual void __Transpose(cwoComplex *pi, cwoComplex *po);
1717 
1718  void JoinH( CWO &a);
1719  void JoinV( CWO &a);
1720 
1722 
1728  void Rect(int m, int n, int Sx, int Sy, cwoComplex &a, int flag=CWO_FILL_INSIDE);
1729  void Rect(const cwoVect2<int> &v, const cwoVect2<int> &s, cwoComplex &a, int flag = CWO_FILL_INSIDE)
1730  {
1731  Rect(v.x, v.y, s.x, s.y, a, flag);
1732  };
1733 
1735 
1742  void Circ(int m, int n, int r, cwoComplex &a, int flag=CWO_FILL_INSIDE);
1743  void Circ(const cwoVect2<int> &v, int r, cwoComplex &a, int flag = CWO_FILL_INSIDE)
1744  {
1745  Circ(v.x, v.y, r, a, flag);
1746  };
1747 
1748  virtual void __CircFillInside(cwoComplex *p, int m, int n, int r, cwoComplex &a);
1749  virtual void __CircFillOutside(cwoComplex *p, int m, int n, int r, cwoComplex &a);
1750 
1752 
1758  void MulCirc(int m, int n, int r, cwoComplex &a);
1759  void MulCirc(const cwoVect2<int> &v, int r, cwoComplex &a){
1760  MulCirc(v.x, v.y, r, a);
1761  }
1762 
1763 
1764 
1767 
1773  void Convolve(CWO &a, float *ker, int Sx, int Sy);
1774 
1775  void ConvolveFFT(CWO *a, CWO *b, CWO *c);
1776 
1777  //
1778  virtual void Expand(int Nx, int Ny);
1779  void ExpandMirror();
1780  void ExpandTwice(CWO *src);
1781  void ExpandHalf(CWO *dst);
1782 
1783 
1784 
1785 
1787 
1801  int ScaleReal(float lim=1.0);
1802 
1804 
1819  int ScaleReal(float i1, float i2, float o1, float o2);
1820 
1821  virtual int __ScaleReal(float i1, float i2, float o1, float o2);
1822 
1824 
1828  int ScaleCplx(float lim=1.0);
1829 
1831 
1841  int ScaleCplx(float i1, float i2, float o1, float o2);
1842 
1843  virtual int __ScaleCplx(float i1, float i2, float o1, float o2);
1844 
1846 
1847 
1848  //###########################################
1852 
1853 
1861  virtual void SetRandSeed(long long int s);
1862  long long int GetRandSeed();
1863 
1865 
1870  float GetRandVal();
1871 
1872 
1873  cwoComplex GetRandComplex();
1874  virtual void __RandPhase(cwoComplex *a, float max, float min);
1875  virtual void __MulRandPhase(cwoComplex *a, float max, float min);
1876 
1877  virtual void RandReal(float max = 1.0f, float min = 0.0f);
1878  void RandPhase(float max=CWO_PI, float min=-CWO_PI);
1879  void SetRandPhase(float max=CWO_PI, float min=-CWO_PI);
1880  void MulRandPhase(float max=CWO_PI, float min=-CWO_PI);
1882 
1883 
1884 
1885 
1886  void __Hanning(cwoComplex *a, int m, int n, int Wx, int Wy);
1887  void __Hamming(cwoComplex *a, int m, int n, int Wx, int Wy);
1888  void __WindowCos(cwoComplex *a, int m, int n, int W, int Wa);
1889  void __WindowTukey(cwoComplex *a, int m, int n, float r);
1890 
1892 
1896  virtual void __WindowFlatCos(cwoComplex *p, int L, int L0, float k);
1897  void __WindowFlatCosCirc(cwoComplex *p, int L, int L0, float k);
1898 
1899  void __WindowSigmoid(cwoComplex *p, int m, int n, float a, float Tmax, float Tmin);
1900 
1901  //cwoComplex inter_nearest(cwoComplex *old, double x, double y);
1902  //cwoComplex inter_linear(cwoComplex *old, float x, float y);
1903  //cwoComplex inter_cubic(cwoComplex *old, float x, float y);
1904 
1905  /*
1906  virtual void __InterNearest(
1907  cwoComplex *p_new, int newNx, int newNy,
1908  cwoComplex *p_old, int oldNx, int oldNy,
1909  float mx, float my,
1910  int xx, int yy);
1911 
1912  virtual void __InterLinear(
1913  cwoComplex *p_new, int newNx, int newNy,
1914  cwoComplex *p_old, int oldNx, int oldNy,
1915  float mx, float my,
1916  int xx, int yy);
1917 
1918  virtual void __InterCubic(
1919  cwoComplex *p_new, int newNx, int newNy,
1920  cwoComplex *p_old, int oldNx, int oldNy,
1921  float mx, float my,
1922  int xx, int yy);
1923 */
1924  virtual void __ResizeNearest(
1925  cwoComplex *p_new, int newNx, int newNy,
1926  cwoComplex *p_old, int oldNx, int oldNy);
1927 
1928  virtual void __ResizeLinear(
1929  cwoComplex *p_new, int newNx, int newNy,
1930  cwoComplex *p_old, int oldNx, int oldNy);
1931 
1932  virtual void __ResizeCubic(
1933  cwoComplex *p_new, int newNx, int newNy,
1934  cwoComplex *p_old, int oldNx, int oldNy);
1935 
1936  virtual void __ResizeLanczos(
1937  cwoComplex *p_new, int newNx, int newNy,
1938  cwoComplex *p_old, int oldNx, int oldNy);
1939 
1940  virtual void __ResizeNearestScale(cwoComplex *p_src, cwoComplex *p_dst, float scale_x, float scale_y);
1941  virtual void __ResizeLinearScale(cwoComplex *p_src, cwoComplex *p_dst, float scale_x, float scale_y);
1942  virtual void __ResizeCubicScale(cwoComplex *p_src, cwoComplex *p_dst, float scale_x, float scale_y);
1943  virtual void __ResizeLanczosScale(cwoComplex *p_src, cwoComplex *p_dst, float scale_x, float scale_y);
1944 
1945 
1946 
1947  void AddNoiseWhite(float max=1.0);
1948  void MulNoiseWhite(float max=1.0);
1949  void AddNoiseGaussian(float mu, float sigma);
1950  void MulNoiseGaussian(float mu, float sigma);
1951 
1952 
1953  //###########################################
1957 
1958 
1960 
1973  void Resize(int dNx, int dNy, int flag=CWO_INTER_LINEAR);
1974  void Resize(const cwoVect2<int> &s, int flag = CWO_INTER_LINEAR){
1975  Resize(s.x, s.y, flag);
1976  }
1977 
1978  void Resize(CWO *dst, float scale_x, float scale_y, int flag = CWO_INTER_LINEAR);
1979 
1980  void Rotate(float deg);
1981 
1982  virtual void AffineAngularSpectrum(float *mat_affine, float px, float py, int flag=CWO_INTER_LINEAR);
1983  void AffineAngularSpectrum(float *mat_affine, const cwoVect2<float> &pitch, int flag = CWO_INTER_LINEAR){
1984  AffineAngularSpectrum(mat_affine, pitch.x, pitch.y, flag);
1985  }
1986 
1987  virtual void ErrorDiffusion(CWO *a = nullptr, int flag = CWO_ED_FS);
1988  virtual void ErrorDiffusionSegmented(CWO *output, int flag = CWO_ED_FS);
1989 
1990 
1991  void RGB2YCbCr(CWO *rgb, CWO *ycbcr);
1992  void YCbCr2RGB(CWO *rgb, CWO *ycbcr);
1993 
1994  void MulMatrix(CWO *a, cwoComplex *b , CWO *c);
1996 
1997  //###########################################
2001  void FourierShift(float m, float n);
2002 
2003  int FFT(int flag=0);
2004  virtual void FFTShift();
2005 
2006  int FFT(int dx, int dy, int flag = 0);
2007  virtual void FFTShift(int dx, int dy);
2008 
2009  virtual void __ScaledFFTCoeff(cwoComplex *p, float sx, float sy, cwoComplex alpha);
2010  virtual void __ScaledFFTKernel(cwoComplex *p, float sx, float sy, cwoComplex alpha);
2011  //void ScaledFFT(int flag=0); //!< Scaled FFT
2012 
2013  virtual void __NUFFT_T1(cwoComplex *p_fld, cwoFloat2 *p_x, int R=2, int Msp=12);
2014  virtual void __NUFFT_T2(cwoComplex *p_fld, cwoFloat2 *p_x, int R=2, int Msp=12);
2015 
2016  //void NUFFT_T1(int R=2, int Msp=12);//!< Non-uniform FFT (Type1)
2017  //void NUFFT_T2(int R=2, int Msp=12);//!< Non-uniform FFT (Type2)
2018 
2019  void NUFFT1(CWO *map, int R=2, int Msp=12);
2020  void NUFFT2(CWO *map, int R=2, int Msp=12);
2021 
2022  void ConvertSamplingMap(int type);
2023  void SamplingMapScaleOnly(int Nx, int Ny, float R, float sgn);
2025 
2026  virtual void Log(float base=10.0f, float eps=1.0f);
2027 
2028  virtual void Gamma(float g);//OK
2029  virtual void Threshold(float max, float min=0.0);//OK
2030  virtual void Threshold(float max_th, float max_val, float min_th, float min_val);
2031  void ThresholdAmp(float max_amp, float min_amp);
2032  void ThresholdToZero(float th);
2033  void ThresholdAbsToZero(float th);
2034  void ThresholdHist(float var);
2035  virtual void Binary(float th=0.0, float max=1.0, float min=0.0);
2036 
2037  virtual void __PickupFloat(float *src, float *pix_p, float pix);
2038  virtual void __PickupCplx(cwoComplex *src, cwoComplex *pix_p, float pix);
2039 
2041 
2046  void Pickup(CWO *a, float pix);
2047 
2048 
2050 
2057  int POC(CWO &ref);
2058 
2059 
2060  //###########################################
2064 
2065  unsigned int Count(float a);
2066  unsigned int CountNonZero();
2067 
2068 
2070 
2077  virtual float Average();
2078 
2079  float Average(int m, int n, int Sx, int Sy);
2080 
2082 
2085  virtual float Variance();
2086 
2087  float Variance(int m, int n, int Sx, int Sy);
2088 
2090 
2093  float StdDev(int m, int n, int Sx, int Sy);
2094 
2095 
2096 
2097  virtual void __MaxMin(cwoComplex *a, float *max, float *min, int *max_x = nullptr, int *max_y = nullptr, int *min_x = nullptr, int *min_y = nullptr);//OK
2098 
2100 
2103  float Max(int *m=NULL, int *n=NULL);
2105 
2108  float Min(int *m = NULL, int *n = NULL);
2109 
2111 
2127  int MaxMin(float *max, float *min, int *max_x = nullptr, int *max_y = nullptr, int *min_x = nullptr, int *min_y = nullptr);//OK
2128 
2130 
2137  void VarianceMap(int sx, int sy);
2138 
2140 
2153  float Histogram(int *hist, int N);
2154 
2156 
2159  virtual cwoComplex TotalSum();
2160 
2161  /*
2162  //
2163  float Variance(float ave);
2164  float SMDx();*/
2165 
2167 
2168 
2169 
2170 
2171  //###########################################
2175 
2177 
2186  virtual float MSE(CWO &ref);
2187 
2189 
2198  float SNR(CWO &ref);
2199 
2201 
2210  float PSNR(CWO &ref);
2211 
2213 
2222  float SSIM(CWO &ref);
2223 
2224  float SpeckleContrast(int m, int n, int Sx, int Sy);
2225 
2227 
2228  //###########################################
2232 
2233  CWO LaplacianPhase();
2234  void PoisonSolverFFT();
2235  void PhaseUnwrap(int type=0);
2236 
2238 
2239 
2240 
2241  //****************
2242  //For PLS
2243  //****************
2244  virtual cwoObjPoint* GetPointBuffer(){ return nullptr; };
2245  virtual int GetPointNum(){return 0;};
2246  virtual void SetPointNum(int num){};
2247  virtual void ScalePoint(float lim){};
2248  void PLS(int flag);
2249 
2250 // virtual void __PLS_Fresnel(){};
2251 // virtual void __PLS_CGH_Fresnel(){};
2252 
2253  virtual void __PLS_Huygens(float ph=0.0f){};
2254  virtual void __PLS_Fresnel(float ph=0.0f){};
2255  virtual void __PLS_CGH_Fresnel(float ph=0.0f){};
2256 
2257 
2258 
2259 
2260 
2261 
2263 //Test code
2265 
2266 //void SamplingMapX(cwoFloat2 *p, int Nx, int Ny, int quadrant);
2267 //void SamplingMapY(cwoFloat2 *p, int Nx, int Ny, int quadrant);
2268 
2269 
2270 //void NUDFT();
2271 
2272 //void __ScaledAngularProp(cwoComplex *a, float px, float py);
2273 
2274 
2275 
2276 
2277 virtual void __InvFx2Fy2(cwoComplex *a);
2278 
2279 //###########################################
2283 float z0;
2284 float pz;
2285 //float fresnel_c_tbl(float x, float z, float z0, float pz);
2286 //float fresnel_s_tbl(float x, float z, float z0, float pz);
2287 //void PrepareTblFresnelApproxInt(int nx, float z0, float pz, int nz);
2288 //void FresnelApproxIntTbl(
2289 // float z, int x1, int y1, int sx1, int sy1);
2290 //
2291 //void FresnelPolyAperture();
2292 
2293 void ParticleField(cwoFloat3 pos, float radius, float amp=1.0f, float init_ph=0.0f);
2294 
2295 
2296 
2297 };
2298 
2299 
2300 
2301 
2302 
2303 //inline cwoComplex cwoCplxPolar(float amp, float ph);
2304 //inline cwoComplex cwoCplxCart(float r, float i);
2305 inline cwoComplex cwoCplxPolar(float amp, float ph)
2306 {
2307  cwoComplex tmp;
2308  CWO_RE(tmp)=amp*cos(ph);
2309  CWO_IM(tmp)=amp*sin(ph);
2310  return tmp;
2311 }
2312 
2313 inline cwoComplex cwoCplxCart(float r, float i)
2314 {
2315  cwoComplex tmp;
2316  CWO_RE(tmp)=r;
2317  CWO_IM(tmp)=i;
2318  return tmp;
2319 }
2320 
2321 
2322 inline double cwoCplx2Arg(cwoComplex a)
2323 {
2324  double arg = atan2(CWO_IM(a), CWO_RE(a));
2325  if (arg < 0.0f) arg += CWO_2PI;
2326  if (arg >= CWO_2PI) arg -= CWO_2PI;
2327  return arg;
2328 }
2329 
2330 
2331 char* cwoFname(const char *arg, ...);
2332 
2333 
2334 
2335 
2336 
2338 
2341 class cwoVect{
2342 public:
2343  float x, y, z;
2344 public:
2345  cwoVect(){};
2346 
2347  cwoVect(float tx, float ty, float tz) {
2348  x = tx; y = ty; z = tz;
2349  };
2350 
2351  void Set(float tx, float ty, float tz){
2352  x = tx; y = ty; z = tz;
2353  }
2354 
2355  cwoVect operator=(cwoVect &a){
2356  cwoVect c;
2357  x =a.x;
2358  y =a.y;
2359  z =a.z;
2360  return *this;
2361  };
2362 
2363  cwoVect operator+(cwoVect &a){ //vector+vector
2364  cwoVect c;
2365  c.x = x+a.x;
2366  c.y = y+a.y;
2367  c.z = z+a.z;
2368  return c;
2369  }
2370 
2371  cwoVect operator-(cwoVect &a){
2372  cwoVect c;
2373  c.x = x - a.x;
2374  c.y = y - a.y;
2375  c.z = z - a.z;
2376  return c;
2377  };
2378 
2379  cwoVect operator*(float a){
2380  cwoVect c;
2381  c.x = x * a;
2382  c.y = y * a;
2383  c.z = z * a;
2384  return c;
2385  };
2386 
2387  cwoVect operator/(float a){
2388  cwoVect c;
2389  if(a == 0.0f) return *this;
2390  c.x = x / a;
2391  c.y = y / a;
2392  c.z = z / a;
2393  return c;
2394  };
2395 
2396  cwoVect& operator+=(cwoVect &a) {
2397  x += a.x;
2398  y += a.y;
2399  z += a.z;
2400  return *this;
2401  };
2402 
2403  cwoVect& operator-=(cwoVect &a) {
2404  x -= a.x;
2405  y -= a.y;
2406  z -= a.z;
2407  return *this;
2408  };
2409 
2410  cwoVect& operator*=(float a){
2411  x *= a;
2412  y *= a;
2413  z *= a;
2414  return *this;
2415  };
2416 
2417  cwoVect& operator/=(float a){
2418  if(a == 0.0f) return *this;
2419  x /= a;
2420  y /= a;
2421  z /= a;
2422  return *this;
2423  };
2424 
2425  float Length(){
2426  double tx = (double)x;
2427  double ty = (double)y;
2428  double tz = (double)z;
2429  return (float)sqrt(tx*tx + ty*ty + tz*tz);
2430  }
2431 
2432  cwoVect Normalize(){
2433  cwoVect tmp;
2434  tmp=*this;
2435 
2436  float len=Length();
2437  len = 1.0f / len;
2438  tmp.x *= len;
2439  tmp.y *= len;
2440  tmp.z *= len;
2441  return tmp;
2442  }
2443 
2444  float Dot(cwoVect &a){
2445  return (x*a.x+y*a.y+z*a.z);
2446  }
2447 
2448  cwoVect Cross(cwoVect &a){
2449  cwoVect vec;
2450  double x1, y1, z1, x2, y2, z2;
2451  x1 = (double)x;
2452  y1 = (double)y;
2453  z1 = (double)z;
2454  x2 = (double)a.x;
2455  y2 = (double)a.y;
2456  z2 = (double)a.z;
2457  vec.x = (float)(y1 * z2 - z1 * y2);
2458  vec.y = (float)(z1 * x2 - x1 * z2);
2459  vec.z = (float)(x1 * y2 - y1 * x2);
2460  return vec;
2461  }
2462 
2463 };
2464 
2465 
2466 
2467 
2469 
2472 class cwoTbl{
2473  float *tbl_sin;
2474  float *tbl_cos;
2475  CWO *tbl_wrp;
2476  int __Nz;
2477 
2479 
2483  void SetNz(int Nz);
2484 public:
2485 
2486  cwoTbl();
2487  ~cwoTbl();
2488 
2490 
2493  int GetNz();
2494 
2496 
2499  void MakeSin(int N);
2500 
2502 
2505  void MakeCos(int N);
2506 
2508 
2516  void MakeWRPTbl(float z, int Nz, float pz, float wl, float px, float py);
2517 
2519 
2521  void ClipWRPTbl();
2522 
2523 
2525 
2529  CWO *GetWRPTbl(int idx);
2530 
2531 
2532 };
2533 
2534 
2535 template<class T>
2536 class cwoMat{
2537 private:
2538  int dim;
2539 public:
2540  T *a;
2541 };
2542 
2543 
2545 
2548 template<typename T>
2549 class cwoRGB{
2550 private:
2551  //float wavelength_r, wavelength_g, wavelength_b;
2552  T cwo_r, cwo_g, cwob;
2553 public:
2554 /* float GetWaveLengthR(){ return cwo_r.GetWaveLength(); };
2555  float GetWaveLengthG(){ return cwo_r.GetWaveLength(); };
2556  float GetWaveLengthB(){ return wavelength_b; };
2557  void SetWaveLengthR(float wl){ wavelength_r=wl; };
2558  void SetWaveLengthG(float wl){ wavelength_g = wl; };
2559  void SetWaveLengthB(float wl){ wavelength_b = wl; };*/
2560 
2561  T* GetCWOR(){ return &cwo_r; };
2562  T* GetCWOG(){ return &cwo_g; };
2563  T* GetCWOB(){ return &cwo_b; };
2564 
2565 };
2566 
2567 
2568 
2569 template<typename T>
2570 class cwoView{
2571  float wavelength;
2572  cwoFloat3 gaze_point;
2573  cwoFloat3 eye_point;
2574  cwoFloat3 src_point;
2575  float threshold_hist_val;
2576  int outNx, outNy;
2577  float pupil_radius;
2578  int type_holo_eye;
2579  int type_eye_gaze;
2580  std::string out_path;
2581  std::string out_name;
2582 
2583 public:
2584  cwoView(){
2585  wavelength = 633e-9;
2586  gaze_point.x = 0.0f;
2587  gaze_point.y = 0.0f;
2588  gaze_point.z = 0.0f;
2589  eye_point.x = 0.0f;
2590  eye_point.x = 0.0f;
2591  eye_point.x = 0.0f;
2592  src_point.x = 0.0f;
2593  src_point.x = 0.0f;
2594  src_point.x = 0.0f;
2595  threshold_hist_val = 5.0f;
2596  outNx = 1024;
2597  outNy = 1024;
2598  pupil_radius = 2e-3;
2599  type_holo_eye = CWO_SHIFTED_ANGULAR;
2600  type_eye_gaze = CWO_SHIFTED_ANGULAR;
2601 
2602  out_path = "./";
2603  out_name = "view_z=%02d_y=%02d_x=%02d.jpg";
2604 
2605  }
2606  void SetWaveLength(float wl){
2607  wavelength = wl;
2608  }
2609  float GetWaveLength(){
2610  return wavelength;
2611  }
2612 
2613  void SetGazePoint(const cwoFloat3 &a){
2614  gaze_point = a;
2615  }
2616  cwoFloat3 GetGazePoint(){
2617  return gaze_point;
2618  }
2619  void SetEyePoint(const cwoFloat3 &a)
2620  {
2621  eye_point = a;
2622  }
2623  cwoFloat3 GetEyePoint()
2624  {
2625  return eye_point;
2626  }
2627  void SetSrcPoint(const cwoFloat3 &a)
2628  {
2629  src_point = a;
2630  }
2631  cwoFloat3 GetSrcPoint()
2632  {
2633  return src_point;
2634  }
2635 
2636  void SetThresholdHist(const float hist){
2637  threshold_hist_val = hist;
2638  }
2639  void SetOutImageSize(const int Nx, const int Ny){
2640  outNx = Nx;
2641  outNy = Ny;
2642  }
2643  int GetOutImageNx(){
2644  return outNx;
2645  }
2646  int GetOutImageNy(){
2647  return outNy;
2648  }
2649 
2650  float GetThresholdHist(){
2651  return threshold_hist_val;
2652  }
2653 
2654  void SetPupilRadius(const float r){
2655  pupil_radius = r;
2656  }
2657  float GetPupilRadius(){
2658  return pupil_radius;
2659  }
2660  void SetTypeSrcToEye(const int type){
2661  type_holo_eye = type;
2662  }
2663  void SetTypeEyeToGaze(const int type){
2664  type_eye_gaze = type;
2665  }
2666  int GetTypeSrcToEye(){
2667  return type_holo_eye;
2668  }
2669  int GetTypeEyeToGaze(){
2670  return type_eye_gaze;
2671  }
2672 
2673  void SetOutPath(const std::string &path){
2674  out_path = path;
2675  }
2676  void SetOutName(const std::string &name){
2677  out_name = name;
2678  }
2679  std::string &GetOutPath(){
2680  return out_path;
2681  }
2682  std::string &GetOutName(){
2683  return out_name;
2684  }
2685 
2687 
2693  void SrcToEye(
2694  T* src, T* dst,
2695  const cwoFloat3 &sp, const cwoFloat3 &ep,
2696  const cwoFloat2 &p_src, const cwoFloat2 &p_eye){
2697 
2698  (*dst) = (*src);
2699 
2700 
2701 
2702  dst->SetSrcOffset(0, 0);
2703  dst->SetDstOffset(ep.x - sp.x, ep.y - sp.y);
2704  //dst->SetSrcOffset(sp.x, sp.y);
2705  //dst->SetDstOffset(ep.x , ep.y );
2706  dst->SetSrcPitch(p_src);
2707  dst->SetDstPitch(p_eye);
2708  //dst->ScaleCplx();
2709  dst->Diffract(ep.z, GetTypeSrcToEye());
2710 
2711 
2712  float pr = GetPupilRadius();
2713  int pri = (int)(pr / p_eye.x);
2714 
2715  printf("pupil_r %f pupil_r(int) %d\n", pr, pri);
2716 
2717  dst->Circ(dst->GetNx() / 2, dst->GetNy() / 2, pri, cwoCplxCart(0, 0), CWO_FILL_OUTSIDE);
2718 
2719  }
2720 
2721 
2723 
2729  void EyeToDst(
2730  T* src, T* dst,
2731  const cwoFloat3 &ep, const cwoFloat3 &gp,
2732  const cwoFloat2 &p_eye, const cwoFloat2 &p_dst){
2733 
2734  (*dst) = (*src);
2735 
2736  double vx = ep.x - gp.x;
2737  double vy = ep.y - gp.y;
2738  double vz = ep.z - gp.z;
2739  double R = sqrt(vx*vx + vy*vy + vz*vz);
2740 
2741  dst->SetOffset(0, 0);
2742  dst->SetDstOffset(-vx, -vy);
2743  dst->SetSrcPitch(p_eye);
2744  dst->SetDstPitch(p_dst);
2745  dst->Diffract(abs(vz), GetTypeEyeToGaze());
2746 
2747  double radx = atan(vx / vz);
2748  double rady = atan(vy / vz);
2749  double wn = (double)dst->GetWaveNum();
2750 
2751  float aff[9];
2752  dst->SetOffset(0, 0);
2753  dst->SetPitch(p_eye);
2754  cwo::AffineRotY(aff, CWO_DEG(+radx));
2755 
2756  dst->Expand(dst->GetNx() * 2, dst->GetNy() * 2);
2757  dst->DiffractAffine(aff);
2758  dst->Expand(dst->GetNx() / 2, dst->GetNy() / 2);
2759  }
2760 
2761 
2762 
2764 
2771  void Eye(T* src, T* dst, const cwoFloat3 &ep, const cwoFloat3 &gp, const cwoFloat3 &sp){
2772 
2773  (*dst) = (*src);
2774 
2775  double vx = ep.x - gp.x;
2776  double vy = ep.y - gp.y;
2777  double vz = ep.z - gp.z;
2778  double R = sqrt(vx*vx + vy*vy + vz*vz);
2779 
2780  double px = src->GetPx();
2781  double py = src->GetPy();
2782 
2783  //printf("vx vy vz %f %f %f \n", vx, vy, vz);
2784 
2785  double radx = atan(vx / vz);
2786  double rady = atan(vy / vz);
2787  double wn = (double)dst->GetWaveNum();
2788 
2789  //printf("rad_x rad_y %f %f\n", CWO_DEG(radx), CWO_DEG(rady));
2790 
2791  dst->SetSrcOffset(0, 0);
2792  dst->SetDstOffset(ep.x - sp.x, ep.y - sp.y);
2793  dst->SetSrcPitch(px, py);
2794  dst->SetDstPitch(px, py);
2795  //dst->ScaleCplx();
2796  dst->Diffract(ep.z, GetTypeSrcToEye());
2797 
2798  //dst->ShowParam();
2799 
2800  //dst->ScaleCplx();
2801 
2802  float pr = GetPupilRadius();
2803  int pri = (int)(pr / dst->GetPx());
2804  printf("pupil_r %f pupil_r(int) %d\n", pr, pri);
2805  dst->Circ(dst->GetNx() / 2, dst->GetNy() / 2, pri, cwoCplxCart(0, 0), CWO_FILL_OUTSIDE);
2806 
2807  float aff[9];
2808  dst->SetOffset(0, 0);
2809  dst->SetPitch(px, py);
2810  cwo::AffineRotY(aff, CWO_DEG(+radx));
2811  //dst->DiffractTilt(aff);
2812 
2813  dst->SetOffset(0, 0);
2814  dst->SetDstOffset(-vx, -vy);
2815  dst->Diffract(abs(vz), GetTypeEyeToGaze());
2816 
2817  }
2818 
2819 
2820 
2821  void ViewXY(
2822  T *src, const float dx, const float dy,
2823  const int Nvx, const int Nvy = 1, char *fname_out = nullptr)
2824  {
2825  T *tmp = (T*)new T;
2826  cwoFloat3 cur_eye_pos = GetEyePoint();
2827  cwoFloat3 cur_gaz_pos = GetGazePoint();
2828  cwoFloat3 cur_src_pos = GetGazePoint();
2829 
2830  for (int i = 0; i<Nvy; ++i){
2831  for (int j = 0; j<Nvx; ++j){
2832 
2833  // SetEyePoint(flt3);
2834 
2835  Eye(src, tmp, cur_eye_pos, cur_gaz_pos, cur_src_pos);
2836 
2837  //tmp->ShowParam();
2838 
2839  int srcNx = src->GetNx();
2840  int srcNy = src->GetNy();
2841  int outNx = GetOutImageNx();
2842  int outNy = GetOutImageNy();
2843 
2844  if (srcNx != outNx || srcNy != outNy) tmp->Resize(outNx, outNy);
2845 
2846  tmp->ScaleCplx();
2847  tmp->Intensity();
2848  float th = GetThresholdHist();
2849  if (th >= 0.0f) tmp->ThresholdHist(th);
2850  if (fname_out == nullptr)
2851  tmp->SaveAsImage(cwoFname("z:\\view_y=%02d_x=%02d.jpg", i, j));
2852  else
2853  tmp->SaveAsImage(cwoFname(fname_out, i));
2854 
2855  cur_eye_pos.x += dx;
2856 
2857  }
2858  cwoFloat3 tmp = GetEyePoint();
2859  cur_eye_pos.x = tmp.x;
2860  cur_eye_pos.y += dy;
2861 
2862  }
2863 
2864  delete tmp;
2865  }
2866 
2868 
2871  void ViewXYDiv(
2872  char *fname_in, const float dx, const float dy,
2873  const int Nvx, const int Nvy = 1, char *fname_out = nullptr)
2874  {
2875 
2876  T *src = (T*)new T;
2877  T *tmp = (T*)new T;
2878  T *tmp2 = (T*)new T;
2879  cwoFloat3 cur_eye_pos = GetEyePoint();
2880  cwoFloat3 cur_gaz_pos = GetGazePoint();
2881  cwoFloat3 cur_src_pos = GetSrcPoint();
2882 
2883  for (int i = 0; i<Nvy; ++i){
2884  for (int j = 0; j < Nvx; ++j){
2885 
2886  // SetEyePoint(flt3);
2887 
2888  FILE *fp = fopen(fname_in, "r");
2889  double src_pos_x, src_pos_y;
2890  char fname[1024];
2891 
2892  tmp->Clear();
2893  while (fscanf(fp, "%lf %lf %s", &src_pos_x, &src_pos_y, fname) != EOF){
2894  src->Load(fname);
2895 
2896  src->SetPitch(1e-6);//
2897 
2898  cur_src_pos.x = src_pos_x;
2899  cur_src_pos.y = src_pos_y;
2900 
2901  printf("%e %e %s \n", src_pos_x, src_pos_y, fname);
2902 
2903  Eye(src, tmp2, cur_eye_pos, cur_gaz_pos, cur_src_pos);
2904 
2905  if (tmp->GetNx() != src->GetNx() || tmp->GetNy() != src->GetNy()){
2906  tmp->Create(src->GetNx(), src->GetNy());
2907  }
2908 
2909  //tmp2->SaveAsImage(cwoFname("%s.jpg",fname), CWO_SAVE_AS_INTENSITY);
2910 
2911  (*tmp) += (*tmp2);
2912  }
2913  fclose(fp);
2914 
2915  //tmp->ShowParam();
2916 
2917  int srcNx = src->GetNx();
2918  int srcNy = src->GetNy();
2919  int outNx = GetOutImageNx();
2920  int outNy = GetOutImageNy();
2921 
2922  if (srcNx != outNx || srcNy != outNy)
2923  tmp->Resize(outNx, outNy);
2924 
2925  tmp->ScaleCplx();
2926  tmp->Intensity();
2927  float th = GetThresholdHist();
2928  if (th >= 0.0f) tmp->ThresholdHist(th);
2929  if (fname_out == nullptr)
2930  tmp->SaveAsImage(cwoFname("z:\\view_y=%02d_x=%02d.jpg", i, j));
2931  else
2932  tmp->SaveAsImage(cwoFname(fname_out, i));
2933 
2934  cur_eye_pos.x += dx;
2935 
2936  }
2937  cwoFloat3 tmp = GetEyePoint();
2938  cur_eye_pos.x = tmp.x;
2939  cur_eye_pos.y += dy;
2940 
2941  }
2942 
2943  delete tmp;
2944  delete tmp2;
2945  delete src;
2946 
2947  }
2948 
2949 
2950 
2951  void ViewXYDiv2(
2952  const std::string &fname_in,
2953  const cwoFloat3 &eye_dv, const cwoFloat3 &gaz_dv,
2954  const cwoFloat3 &Nv,
2955  const cwoFloat2 &p_src, const cwoFloat2 &p_eye, const cwoFloat2 &p_dst){
2956 
2957  T *src = (T*)new T;
2958  T *tmp = (T*)new T;
2959  T *tmp2 = (T*)new T;
2960 
2961  cwoFloat3 cur_eye_pos = GetEyePoint();
2962  cwoFloat3 cur_gaz_pos = GetGazePoint();
2963  cwoFloat3 cur_src_pos = GetSrcPoint();
2964 
2965  cwoFloat3 tmp_flt3;
2966 
2967  int Nvx = Nv.x;
2968  int Nvy = Nv.y;
2969  int Nvz = Nv.z;
2970 
2971  for (int k = 0; k < Nvz; ++k){
2972  for (int i = 0; i < Nvy; ++i){
2973  for (int j = 0; j < Nvx; ++j){
2974 
2975  // SetEyePoint(flt3);
2976 
2977  FILE *fp = fopen(fname_in.c_str(), "r");
2978 
2979  double src_pos_x, src_pos_y;
2980  char fname[1024];
2981 
2982  tmp->Clear();
2983  while (fscanf(fp, "%lf %lf %s", &src_pos_x, &src_pos_y, fname) != EOF){
2984  src->Load(fname);
2985 
2986  cur_src_pos.x = src_pos_x;
2987  cur_src_pos.y = src_pos_y;
2988 
2989  printf("%e %e %s \n", src_pos_x, src_pos_y, fname);
2990 
2991  src->SetWaveLength(GetWaveLength());
2992  tmp->SetWaveLength(GetWaveLength());
2993  tmp2->SetWaveLength(GetWaveLength());
2994 
2995  SrcToEye(src, tmp2, cur_src_pos, cur_eye_pos, p_src, p_eye);
2996 
2997  if (tmp->GetNx() != src->GetNx() || tmp->GetNy() != src->GetNy()){
2998  tmp->Create(src->GetNx(), src->GetNy());
2999  }
3000 
3001  (*tmp) += (*tmp2);
3002  }
3003  fclose(fp);
3004 
3005  (*src) = (*tmp);
3006  EyeToDst(src, tmp, cur_eye_pos, cur_gaz_pos, p_eye, p_dst);
3007 
3008  int srcNx = src->GetNx();
3009  int srcNy = src->GetNy();
3010  int outNx = GetOutImageNx();
3011  int outNy = GetOutImageNy();
3012 
3013  if (srcNx != outNx || srcNy != outNy)
3014  tmp->Resize(outNx, outNy);
3015 
3016  tmp->ScaleCplx();
3017  tmp->Intensity();
3018  float th = GetThresholdHist();
3019  if (th >= 0.0f) tmp->ThresholdHist(th);
3020 
3021  /*if (fname_out == nullptr)
3022  tmp->SaveAsImage(cwoFname("z:\\view_z=%02d_y=%02d_x=%02d.jpg",k, i, j));
3023  else
3024  tmp->SaveAsImage(cwoFname(fname_out, i));
3025  */
3026  tmp->SaveAsImage(cwoFname((GetOutPath() + GetOutName()).c_str(), k, i, j));
3027 
3028 
3029  cur_eye_pos.x += eye_dv.x;
3030  cur_gaz_pos.x += gaz_dv.x;
3031 
3032  }
3033 
3034  tmp_flt3 = GetEyePoint();
3035  cur_eye_pos.x = tmp_flt3.x;
3036  cur_eye_pos.y += eye_dv.y;
3037 
3038  tmp_flt3 = GetGazePoint();
3039  cur_gaz_pos.x = tmp_flt3.x;
3040  cur_gaz_pos.y += gaz_dv.y;
3041 
3042  }
3043 
3044  tmp_flt3 = GetEyePoint();
3045  cur_eye_pos.x = tmp_flt3.x;
3046  cur_eye_pos.y = tmp_flt3.y;
3047  cur_eye_pos.z += eye_dv.z;
3048 
3049  tmp_flt3 = GetGazePoint();
3050  cur_gaz_pos.x = tmp_flt3.x;
3051  cur_gaz_pos.y = tmp_flt3.y;
3052  cur_gaz_pos.z += gaz_dv.z;
3053 
3054  }
3055 
3056  delete tmp;
3057  delete tmp2;
3058  delete src;
3059 
3060  }
3061 
3062  void Divide(const std::string &fname_in, int Nx, int Ny, float px, float py){
3063  //std::unique_ptr<CWO> c_org;
3064  //std::unique_ptr<CWO> c_new;
3065 
3066  CWO *c_org = (CWO*)new CWO;
3067  CWO *c_new = (CWO*)new CWO;
3068 
3069  c_org->Load(fname_in);
3070  int org_Nx = c_org->GetNx();
3071  int org_Ny = c_org->GetNy();
3072 
3073  int new_Nx = int(ceil((float)org_Nx/Nx)*Nx);
3074  int new_Ny = int(ceil((float)org_Ny/Ny)*Ny);
3075 
3076 
3077  c_org->Expand(new_Nx, new_Ny);
3078  c_new->Create(Nx, Ny);
3079 
3080 
3081  std::string fname_out = "../../data/wide_reconst/hol_div_%02d%02d.cwo";
3082  std::string fname_txt = "../../data/wide_reconst/div_out.txt";
3083 
3084 
3085  FILE *fp = fopen(fname_txt.c_str(), "wt");
3086 
3087  for (int i = 0; i < new_Ny; i += Ny){
3088  for (int j = 0; j < new_Nx; j += Nx){
3089  c_org->GetPixel(j, i, *c_new);
3090 
3091  int idx_x = int((float)j / Nx);
3092  int idx_y = int((float)i / Ny);
3093 
3094  std::string name = cwoFname(fname_out.c_str(), idx_y, idx_x);
3095  c_new->Save(name);
3096  c_new->SaveAsImage(name+".jpg");
3097 
3098  float center_x = CWO_X(((j - new_Nx / 2) + Nx / 2)*px);
3099  float center_y = CWO_Y(((i - new_Ny / 2) + Ny / 2)*py);
3100 
3101  std::string txt =
3102  std::to_string(center_x) + " " + std::to_string(center_y) + " " + name + "\n";
3103 
3104  fwrite(txt.c_str(),1, strlen(txt.c_str()), fp);
3105 
3106  }
3107  }
3108 
3109  fclose(fp);
3110 
3111  delete c_org;
3112  delete c_new;
3113 
3114  }
3115 
3116 
3117 };
3118 
3119 
3120 
3121 #endif
3122 
3123 
3124 
3125 
Generated on Mon Apr 2 2018 13:35:29 for CWO++library by   doxygen 1.8.1