Openholo  v1.0
Open Source Digital Holographic Library
ophWaveAberration.cpp
Go to the documentation of this file.
1 /*M///////////////////////////////////////////////////////////////////////////////////////
2 //
3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
4 //
5 // By downloading, copying, installing or using the software you agree to this license.
6 // If you do not agree to this license, do not download, install, copy or use the software.
7 //
8 //
9 // License Agreement
10 // For Open Source Digital Holographic Library
11 //
12 // Openholo library is free software;
13 // you can redistribute it and/or modify it under the terms of the BSD 2-Clause license.
14 //
15 // Copyright (C) 2017-2024, Korea Electronics Technology Institute. All rights reserved.
16 // E-mail : contact.openholo@gmail.com
17 // Web : http://www.openholo.org
18 //
19 // Redistribution and use in source and binary forms, with or without modification,
20 // are permitted provided that the following conditions are met:
21 //
22 // 1. Redistribution's of source code must retain the above copyright notice,
23 // this list of conditions and the following disclaimer.
24 //
25 // 2. Redistribution's in binary form must reproduce the above copyright notice,
26 // this list of conditions and the following disclaimer in the documentation
27 // and/or other materials provided with the distribution.
28 //
29 // This software is provided by the copyright holders and contributors "as is" and
30 // any express or implied warranties, including, but not limited to, the implied
31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
32 // In no event shall the copyright holder or contributors be liable for any direct,
33 // indirect, incidental, special, exemplary, or consequential damages
34 // (including, but not limited to, procurement of substitute goods or services;
35 // loss of use, data, or profits; or business interruption) however caused
36 // and on any theory of liability, whether in contract, strict liability,
37 // or tort (including negligence or otherwise) arising in any way out of
38 // the use of this software, even if advised of the possibility of such damage.
39 //
40 // This software contains opensource software released under GNU Generic Public License,
41 // NVDIA Software License Agreement, or CUDA supplement to Software License Agreement.
42 // Check whether software you use contains licensed software.
43 //
44 //M*/
45 
46 #include "ophWaveAberration.h"
47 
48 
49 
50 inline double ophWaveAberration::factorial(double x)
51 {
52  if (x == 0)
53  return (1);
54  else
55  return (x == 1 ? x : x * factorial(x - 1));
56 }
57 
58 
59 ophWaveAberration::ophWaveAberration() : nOrder(0), mFrequency(0)
60 {
61 
62  cout << "ophWaveAberration Constructor" << endl;
63  uint wavelength_num = 1;
64 
65  complex_H = new Complex<Real>*[wavelength_num];
66  context_.wave_length = new Real[wavelength_num];
67 }
68 
70 {
71  cout << "ophWaveAberration Destructor" << endl;
72 }
73 
74 
75 double** ophWaveAberration::calculateZernikePolynomial(double n, double m, vector<double> x, vector<double> y, double d)
76 {
77  vector<double>::size_type x_max = x.size();
78  vector<double>::size_type y_max = y.size();
79  double radius = d / 2;
80  double N;
81  double r;
82  double theta;
83  double co ;
84  double si ;
85 
86  double **Z = new double*[x_max];
87  double **A = new double*[x_max];
88  for(int i = 0; i < (int)x_max; i++)
89  {
90  A[i] = new double[y_max];
91  Z[i] = new double[y_max];
92 
93  // memset(A[i], 0, y_max*sizeof(double));
94  }
95 
96  for(int ix = 0; ix < (int)x_max; ix++)
97  {
98  for(int iy = 0; iy < (int)y_max; iy++)
99  {
100  A[ix][iy] = (sqrt(pow(x[ix],2) + pow(y[iy],2)) <= radius);
101  };
102  }
103  // Start : Calculate Zernike polynomial
104 
105  N = sqrt(2 * (n + 1) / (1 + (m == 0))); // Calculate Normalization term
106 
107  if (n == 0)
108  {
109  for(int i=0; i<(int)x_max; i++)
110  memcpy(Z[i], A[i],y_max*sizeof(double));
111  }
112  else
113  {
114  for(int i = 0; i<(int)x_max; i++)
115  memset(Z[i],0, y_max*sizeof(double));
116 
117  for(int ix = 0; ix < (int)x_max; ix++)
118  {
119  for(int iy = 0; iy < (int)y_max; iy++)
120  {
121  r = sqrt(pow(x[ix], 2) + pow(y[iy],2));
122 
123  if (((x[ix] >= 0) && (y[iy] >= 0)) || ((x[ix] >= 0) & (y[iy] < 0)))
124  theta = atan(y[iy] / (x[ix] + 1e-30));
125  else
126  theta = M_PI + atan(y[iy] / (x[ix] + 1e-30));
127 
128  for(int s = 0; s <= (n - abs(m)) / 2; s++)
129  {
130  Z[ix][iy] = Z[ix][iy] + pow((-1),s)*factorial(n - s)*pow((r/radius),(n - 2 * s)) /
131  (factorial(s)*factorial((n + abs(m))/2 - s)*factorial((n - abs(m)) / 2 - s));
132  }
133  co = cos(m*theta);
134  si = sin(m*theta);
135  Z[ix][iy] = A[ix][iy]*N*Z[ix][iy]*((m >= 0)*co - (m < 0)*si);
136  }
137  }
138  }
139  // End : Calculate Zernike polynomial
140  for (int i=0; i < x_max; i++)
141  {
142  delete[] A[i];
143  }
144  delete[] A;
145 
146  return Z;
147 }
148 
149 
150 void ophWaveAberration::imresize(double **X, int Nx, int Ny, int nx, int ny, double **Y)
151 {
152  int fx, fy;
153  double x, y, tx, tx1, ty, ty1, scale_x, scale_y;
154 
155  scale_x = (double)Nx / (double)nx;
156  scale_y = (double)Ny / (double)ny;
157 
158  for (int i = 0; i < nx; i++)
159  {
160  x = (double)i * scale_x;
161 
162  fx = (int)floor(x);
163  tx = x - fx;
164  tx1 = double(1.0) - tx;
165  for (int j = 0; j < ny; j++)
166  {
167  y = (double)j * scale_y;
168  fy = (int)floor(y);
169  ty = y - fy;
170  ty1 = double(1.0) - ty;
171 
172  Y[i][j] = X[fx][fy] * (tx1*ty1) + X[fx][fy + 1] * (tx1*ty) + X[fx + 1][fy] * (tx*ty1) + X[fx + 1][fy + 1] * (tx*ty);
173  }
174  }
175 }
176 
177 
178 
180 {
181  auto start_time = CUR_TIME;
182  const oph::Complex<Real> j(0,1);
183 
184  double wave_lambda = context_.wave_length[0]; // wavelength
185  int z_max = sizeof(zernikeCoefficent)/sizeof(zernikeCoefficent[0]);
186  double *ZC;
187  ZC = zernikeCoefficent;
188 
189 
190  double n, m;
191  double dxa = context_.pixel_pitch[_X]; // Sampling interval in x axis of exit pupil
192  double dya = context_.pixel_pitch[_Y]; // Sampling interval in y axis of exit pupil
193  unsigned int xr = context_.pixel_number[_X];
194  unsigned int yr = context_.pixel_number[_Y]; // Resolution in x, y axis of exit pupil
195 
196  double DE = max(dxa*xr, dya*yr); // Diameter of exit pupil
197  double scale = 1.3;
198 
199  DE = DE * scale;
200 
201  vector<double> xn;
202  vector<double> yn;
203 
204  double max_xn = floor(DE/dxa+1);
205  double max_yn = floor(DE/dya+1);
206 
207  xn.reserve((int)max_xn);
208  for (int i = 0; i < (int)max_xn; i++)
209  {
210  xn.push_back(-DE / 2 + dxa*i);
211  } // x axis coordinate of exit pupil
212 
213  yn.reserve((int)max_yn);
214  for (int i = 0; i < max_yn; i++)
215  {
216  yn.push_back(-DE / 2 + dya*i);
217  }// y axis coordinate of exit pupil
218 
219  double d = DE;
220 
221  vector<double>::size_type length_xn = xn.size();
222  vector<double>::size_type length_yn = yn.size();
223 
224  double **W = new double*[(int)length_xn];
225  double **Temp_W = new double*[(int)length_xn];
226 
227  for (int i = 0; i < (int)length_xn; i++)
228  {
229  W[i] = new double[length_yn];
230  Temp_W[i] = new double[length_yn];
231  }
232 
233  for (int i = 0; i < (int)length_xn; i++)
234  {
235  memset(W[i], 0, length_yn*sizeof(double));
236  memset(Temp_W[i], 0, length_yn * sizeof(double));
237  }
238 
239 
240  // Start : Wavefront Aberration Generation
241  for (int i = 0; i <z_max; i++)
242  {
243  if (ZC[i] != 0)
244  {
245  n = ceil((-3 + sqrt(9 + 8 * i)) / 2); // order of the radial polynomial term
246  m = 2 * i - n * (n + 2); // frequency of the sinusoidal component
247 
248  Temp_W = calculateZernikePolynomial(n, m, xn, yn, d);
249 
250  for(int ii = 0; ii < length_xn; ii++)
251  {
252  for (int jj = 0; jj < length_yn; jj++)
253  {
254  W[ii][jj] = W[ii][jj] + ZC[i] * Temp_W[ii][jj];
255  }
256  }
257  }
258  }
259  // End : Wavefront Aberration Generation
260 
261 
262  for (int i = 0; i < (int)length_xn; i++)
263  {
264  memset(Temp_W[i], 0, length_yn * sizeof(double));
265  }
266 
267  int min_xnn, max_xnn;
268  int min_ynn, max_ynn;
269 
270  min_xnn = (int)round(length_xn / 2 - xr / 2);
271  max_xnn = (int)round(length_xn / 2 + xr / 2 + 1);
272  min_ynn = (int)round(length_yn / 2 - yr / 2);
273  max_ynn = (int)round(length_yn / 2 + yr / 2 + 1);
274 
275  int length_xnn, length_ynn;
276  length_xnn = max_xnn - min_xnn;
277  length_ynn = max_ynn - min_ynn;
278 
279  double **WT = new double*[length_xnn];
280  for (int i = 0; i < length_xnn; i++)
281  {
282  WT[i] = new double[length_ynn];
283  memset(WT[i], 0, length_ynn * sizeof(double));
284  }
285 
286  for (int i = 0; i < length_xnn; i++)
287  {
288  for (int j = 0; j < length_ynn; j++)
289  {
290  WT[i][j] = W[min_xnn+i][min_ynn+j];
291  }
292  }
293 
294  double **WS = new double*[(int)xr];
295  for (int i = 0; i < (int)xr; i++)
296  {
297  WS[i] = new double[yr];
298  memset(WS[i], 0, yr * sizeof(double));
299  }
300 
301  imresize(WT, length_xnn, length_ynn, xr, yr, WS);
302 
303 
304  oph::Complex<Real> **WD = new oph::Complex<Real>*[xr];
305 
306  for(int i = 0; i < (int)xr; i++)
307  WD[i] = new oph::Complex<Real>[yr];
308 
309  for(int ii = 0; ii < (int)xr; ii ++ )
310  {
311  for (int jj = 0; jj < (int)yr; jj++)
312  {
313 
314  WD[ii][jj]= exp(-j * (oph::Complex<Real>)2 * M_PI*WS[ii][jj] / wave_lambda); // Wave Aberration Complex Field
315  }
316  }
317  //WD[x][y]
318 
319  for (int i = 0; i < (int)length_xn; i++)
320  {
321  delete [] W[i];
322  delete [] Temp_W[i];
323  }
324  delete[] W;
325  delete[] Temp_W;
326 
327  for (int i = 0; i < (int)xr; i++)
328  {
329  delete[] WS[i];
330  }
331  delete[] WS;
332 
333  for (int i = 0; i < (int)length_xnn; i++)
334  {
335  delete[] WT[i];
336  }
337  delete[] WT;
338 
339  complex_W = WD;
340 
341  for (int x = 0; x < xr; x++) {
342  for (int y = 0; y < yr; y++) {
343  complex_H[0][x + y * xr] = complex_W[x][y];
344  }
345  }
346 
347 // return WD;
348 
349  auto end_time = CUR_TIME;
350 
351  auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
352 
353  LOG("Implement time : %.5lf sec\n", during_time);
354 }
355 
356 
357 void ophWaveAberration::Free2D(oph::Complex<Real> ** doublePtr)
358 {
359  for (int i = 0; i < (int)context_.pixel_number[_X]; i++)
360  {
361  delete[] doublePtr[i];
362  }
363 }
364 
366 {
367  this->Free2D(complex_W);
368  std::cout << " ophFree" << std::endl;
369 }
370 
371 
372 bool ophWaveAberration::readConfig(const char* fname)
373 {
374  LOG("Reading....%s...\n", fname);
375 
376 
377  /*XML parsing*/
378  tinyxml2::XMLDocument xml_doc;
379  tinyxml2::XMLNode *xml_node;
380  tinyxml2::XMLElement *xml_element;
381  const tinyxml2::XMLAttribute *xml_attribute;
382 
383 
384  if (!checkExtension(fname, ".xml"))
385  {
386  LOG("file's extension is not 'xml'\n");
387  return false;
388  }
389  auto ret = xml_doc.LoadFile(fname);
390  if (ret != tinyxml2::XML_SUCCESS)
391  {
392  LOG("Failed to load file \"%s\"\n", fname);
393  return false;
394  }
395 
396  xml_node = xml_doc.FirstChild();
397  xml_element = xml_node->FirstChildElement("Wavelength");
398  if (!xml_element || tinyxml2::XML_SUCCESS != xml_element->QueryDoubleText(&context_.wave_length[0]))
399  return false;
400 
401 
402  xml_element = xml_node->FirstChildElement("PixelPitchHor");
403  if (!xml_element || tinyxml2::XML_SUCCESS != xml_element->QueryDoubleText(&context_.pixel_pitch[_X]))
404  return false;
405 
406  xml_element = xml_node->FirstChildElement("PixelPitchVer");
407  if (!xml_element || tinyxml2::XML_SUCCESS != xml_element->QueryDoubleText(&context_.pixel_pitch[_Y]))
408  return false;
409 
410  xml_element = xml_node->FirstChildElement("ResolutionHor");
411  if (!xml_element || tinyxml2::XML_SUCCESS != xml_element->QueryIntText(&context_.pixel_number[_X]))
412  return false;
413 
414  xml_element = xml_node->FirstChildElement("ResolutionVer");
415  if (!xml_element || tinyxml2::XML_SUCCESS != xml_element->QueryIntText(&context_.pixel_number[_Y]))
416  return false;
417 
418  xml_element = xml_node->FirstChildElement("ZernikeCoeff");
419  xml_attribute = xml_element->FirstAttribute();
420 
421  for(int i=0; i< 45; i++)
422  {
423  if (!xml_attribute || tinyxml2::XML_SUCCESS != xml_attribute->QueryDoubleValue(&zernikeCoefficent[i]))
424  return false;
425  xml_attribute=xml_attribute->Next();
426 
427  }
428 
429  pixelPitchX = context_.pixel_pitch[_X];
430  pixelPitchY = context_.pixel_pitch[_Y];
431 
434 
435  waveLength = *context_.wave_length;
436 
440 
441  cout << "Wavelength: " << context_.wave_length[0] << endl;
442  cout << "PixelPitch(Horizontal): " << context_.pixel_pitch[_X] << endl;
443  cout << "PixelPitch(Vertical): " << context_.pixel_pitch[_Y] << endl;
444  cout << "Resolution(Horizontal): " << context_.pixel_number[_X] << endl;
445  cout << "Resolution(Vertical): " << context_.pixel_number[_Y] << endl;
446  cout << "Zernike Coefficient: " << endl;
447  for(int i=0; i<45; i++)
448  {
449  if (i!=0 && (i+1)%5 == 0)
450  cout << "z["<<i<<"]="<< zernikeCoefficent[i]<<endl;
451  else
452  cout << "z[" << i << "]=" << zernikeCoefficent[i] <<" ";
453  zernikeCoefficent[i] = zernikeCoefficent[i] * context_.wave_length[0];
454  }
455  int xr = context_.pixel_number[_X];
456  int yr = context_.pixel_number[_Y];
457 
458  complex_H[0] = new Complex<Real>[xr * yr];
459 
460  return true;
461 
462 }
463 
464 void ophWaveAberration::saveAberration(const char* fname)
465 {
466  ofstream fout(fname, ios_base::out | ios_base::binary);
467  fout.write((char *)complex_W, context_.pixel_number[_X] * context_.pixel_number[_Y] * sizeof(oph::Complex<Real>));
468  fout.close();
469 }
470 
471 void ophWaveAberration::readAberration(const char* fname)
472 {
473 
474  complex_W = new oph::Complex<Real>*[context_.pixel_number[_X]];
475  for (int i = 0; i < (int)context_.pixel_number[_X]; i++)
476  complex_W[i] = new oph::Complex<Real>[context_.pixel_number[_Y]];
477 
478  ifstream fin(fname, ios_base::in | ios_base::binary);
480  fin.close();
481 }
482 
483 bool ophWaveAberration::loadAsOhc(const char * fname)
484 {
485  if (!Openholo::loadAsOhc(fname)) {
486  LOG("Failed load file");
487  return false;
488  }
489 
490  pixelPitchX = context_.pixel_pitch[_X];
491  pixelPitchY = context_.pixel_pitch[_Y];
492 
493  int xr = resolutionX = context_.pixel_number[_X];
494  int yr = resolutionY = context_.pixel_number[_Y];
495 
496  waveLength = context_.wave_length[0];
497 #if 0
498  int nColor = 1;
499  for (int c = 0; c < nColor; c++) {
500  for (int x = 0; x < xr; x++) {
501  for (int y = 0; y < yr; y++) {
502  complex_W[x][y] = complex_H[c][x + y * xr];
503  }
504  }
505  }
506 #else
507  for (int x = 0; x < xr; x++) {
508  for (int y = 0; y < yr; y++) {
509  complex_W[x][y] = complex_H[0][x + y * xr];
510  }
511  }
512 #endif
513  return true;
514 }
void abs(const oph::Complex< T > &src, oph::Complex< T > &dst)
Definition: function.h:112
Real factorial(double x)
Factorial using recursion.
ophWaveAberration()
Constructor.
oph::Complex< Real > ** complex_W
double pointer of the 2D data array of a wave aberration
Real * wave_length
Definition: Openholo.h:69
~ophWaveAberration()
Destructor.
void setPixelNumberOHC(const ivec2 pixel_number)
getter/setter for OHC file read and write
Definition: Openholo.h:311
virtual bool loadAsOhc(const char *fname)
Function to read OHC file.
double ** calculateZernikePolynomial(double n, double m, vector< double > x, vector< double > y, double d)
Calculates Zernike polynomial.
void accumulateZernikePolynomial()
Sums up the calculated Zernike polynomials.
XMLError QueryIntText(int *ival) const
Definition: tinyxml2.cpp:1638
float Real
Definition: typedef.h:55
#define CUR_TIME
Definition: function.h:58
void setPixelPitchOHC(const vec2 pixel_pitch)
Definition: Openholo.h:314
bool readConfig(const char *fname)
read configuration from a configration file
bool checkExtension(const char *fname, const char *ext)
Functions for extension checking.
Definition: Openholo.cpp:80
void ophFree(void)
Pure virtual function for override in child classes.
const XMLNode * FirstChild() const
Get the first child node, or null if none exists.
Definition: tinyxml2.h:761
#define _Y
Definition: define.h:84
void setWavelengthOHC(const Real wavelength, const LenUnit wavelength_unit)
Definition: Openholo.h:317
#define _X
Definition: define.h:80
void readAberration(const char *fname)
reads the 2D data array of a wave aberration from a file
oph::ivec2 pixel_number
Definition: Openholo.h:63
void Free2D(oph::Complex< Real > **doublePtr)
deletes 2D memory array using double pointer
void saveAberration(const char *fname)
saves the 2D data array of a wave aberration into a file
XMLError QueryDoubleValue(double *value) const
See QueryIntValue.
Definition: tinyxml2.cpp:1425
virtual bool loadAsOhc(const char *fname)
Function to read OHC file.
Definition: Openholo.cpp:249
XMLError LoadFile(const char *filename)
Definition: tinyxml2.cpp:2157
XMLError QueryDoubleText(double *dval) const
See QueryIntText()
Definition: tinyxml2.cpp:1690
OphConfig context_
Definition: Openholo.h:297
Complex< Real > ** complex_H
Definition: Openholo.h:298
const XMLAttribute * Next() const
The next attribute in the list.
Definition: tinyxml2.h:1146
const XMLElement * FirstChildElement(const char *name=0) const
Definition: tinyxml2.cpp:940
unsigned int uint
Definition: typedef.h:62
void imresize(double **X, int Nx, int Ny, int nx, int ny, double **Y)
Resizes 2D data using bilinear interpolation.
#define M_PI
Definition: define.h:52
const XMLAttribute * FirstAttribute() const
Return the first attribute in the list.
Definition: tinyxml2.h:1471
oph::vec2 pixel_pitch
Definition: Openholo.h:64