Openholo  v4.2
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 
60  : nOrder(0)
61  , mFrequency(0)
62  , complex_W(nullptr)
63 {
64 
65  cout << "ophWaveAberration Constructor" << endl;
66  uint wavelength_num = 1;
67 
68  complex_H = new Complex<Real>*[wavelength_num];
69  context_.wave_length = new Real[wavelength_num];
70 }
71 
73 {
74  cout << "ophWaveAberration Destructor" << endl;
75 }
76 
77 
78 double** ophWaveAberration::calculateZernikePolynomial(double n, double m, vector<double> x, vector<double> y, double d)
79 {
80  vector<double>::size_type x_max = x.size();
81  vector<double>::size_type y_max = y.size();
82  double radius = d / 2;
83  double N;
84  double r;
85  double theta;
86  double co ;
87  double si ;
88 
89  double **Z = new double*[x_max];
90  double **A = new double*[x_max];
91  for(int i = 0; i < (int)x_max; i++)
92  {
93  A[i] = new double[y_max];
94  Z[i] = new double[y_max];
95 
96  // memset(A[i], 0, y_max*sizeof(double));
97  }
98 
99  for(int ix = 0; ix < (int)x_max; ix++)
100  {
101  for(int iy = 0; iy < (int)y_max; iy++)
102  {
103  A[ix][iy] = (sqrt(pow(x[ix],2) + pow(y[iy],2)) <= radius);
104  };
105  }
106  // Start : Calculate Zernike polynomial
107 
108  N = sqrt(2 * (n + 1) / (1 + (m == 0))); // Calculate Normalization term
109 
110  if (n == 0)
111  {
112  for(int i=0; i<(int)x_max; i++)
113  memcpy(Z[i], A[i],y_max*sizeof(double));
114  }
115  else
116  {
117  for(int i = 0; i<(int)x_max; i++)
118  memset(Z[i],0, y_max*sizeof(double));
119 
120  for(int ix = 0; ix < (int)x_max; ix++)
121  {
122  for(int iy = 0; iy < (int)y_max; iy++)
123  {
124  r = sqrt(pow(x[ix], 2) + pow(y[iy],2));
125 
126  if (((x[ix] >= 0) && (y[iy] >= 0)) || ((x[ix] >= 0) & (y[iy] < 0)))
127  theta = atan(y[iy] / (x[ix] + 1e-30));
128  else
129  theta = M_PI + atan(y[iy] / (x[ix] + 1e-30));
130 
131  for(int s = 0; s <= (n - abs(m)) / 2; s++)
132  {
133  Z[ix][iy] = Z[ix][iy] + pow((-1),s)*factorial(n - s)*pow((r/radius),(n - 2 * s)) /
134  (factorial(s)*factorial((n + abs(m))/2 - s)*factorial((n - abs(m)) / 2 - s));
135  }
136  co = cos(m*theta);
137  si = sin(m*theta);
138  Z[ix][iy] = A[ix][iy]*N*Z[ix][iy]*((m >= 0)*co - (m < 0)*si);
139  }
140  }
141  }
142  // End : Calculate Zernike polynomial
143  for (size_t i=0; i < x_max; i++)
144  {
145  delete[] A[i];
146  }
147  delete[] A;
148 
149  return Z;
150 }
151 
152 
153 void ophWaveAberration::imresize(double **X, int Nx, int Ny, int nx, int ny, double **Y)
154 {
155  int fx, fy;
156  double x, y, tx, tx1, ty, ty1, scale_x, scale_y;
157 
158  scale_x = (double)Nx / (double)nx;
159  scale_y = (double)Ny / (double)ny;
160 
161  for (int i = 0; i < nx; i++)
162  {
163  x = (double)i * scale_x;
164 
165  fx = (int)floor(x);
166  tx = x - fx;
167  tx1 = double(1.0) - tx;
168  for (int j = 0; j < ny; j++)
169  {
170  y = (double)j * scale_y;
171  fy = (int)floor(y);
172  ty = y - fy;
173  ty1 = double(1.0) - ty;
174 
175  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);
176  }
177  }
178 }
179 
180 
181 
183 {
184  auto start_time = CUR_TIME;
185  const oph::Complex<Real> j(0,1);
186 
187  double wave_lambda = context_.wave_length[0]; // wavelength
188  int z_max = sizeof(zernikeCoefficent)/sizeof(zernikeCoefficent[0]);
189  double *ZC;
190  ZC = zernikeCoefficent;
191 
192 
193  double n, m;
194  double dxa = context_.pixel_pitch[_X]; // Sampling interval in x axis of exit pupil
195  double dya = context_.pixel_pitch[_Y]; // Sampling interval in y axis of exit pupil
196  unsigned int xr = context_.pixel_number[_X];
197  unsigned int yr = context_.pixel_number[_Y]; // Resolution in x, y axis of exit pupil
198 
199  double DE = max(dxa*xr, dya*yr); // Diameter of exit pupil
200  double scale = 1.3;
201 
202  DE = DE * scale;
203 
204  vector<double> xn;
205  vector<double> yn;
206 
207  double max_xn = floor(DE/dxa+1);
208  double max_yn = floor(DE/dya+1);
209 
210  xn.reserve((int)max_xn);
211  for (int i = 0; i < (int)max_xn; i++)
212  {
213  xn.push_back(-DE / 2 + dxa*i);
214  } // x axis coordinate of exit pupil
215 
216  yn.reserve((int)max_yn);
217  for (int i = 0; i < max_yn; i++)
218  {
219  yn.push_back(-DE / 2 + dya*i);
220  }// y axis coordinate of exit pupil
221 
222  double d = DE;
223 
224  vector<double>::size_type length_xn = xn.size();
225  vector<double>::size_type length_yn = yn.size();
226 
227  double **W = new double*[(int)length_xn];
228  double **Temp_W = new double*[(int)length_xn];
229 
230  for (int i = 0; i < (int)length_xn; i++)
231  {
232  W[i] = new double[length_yn];
233  Temp_W[i] = new double[length_yn];
234  }
235 
236  for (int i = 0; i < (int)length_xn; i++)
237  {
238  memset(W[i], 0, length_yn*sizeof(double));
239  memset(Temp_W[i], 0, length_yn * sizeof(double));
240  }
241 
242 
243  // Start : Wavefront Aberration Generation
244  for (int i = 0; i < z_max; i++)
245  {
246  if (ZC[i] != 0)
247  {
248  n = ceil((-3 + sqrt(9 + 8 * i)) / 2); // order of the radial polynomial term
249  m = 2 * i - n * (n + 2); // frequency of the sinusoidal component
250 
251  Temp_W = calculateZernikePolynomial(n, m, xn, yn, d);
252 
253  for(size_t ii = 0; ii < length_xn; ii++)
254  {
255  for (size_t jj = 0; jj < length_yn; jj++)
256  {
257  W[ii][jj] = W[ii][jj] + ZC[i] * Temp_W[ii][jj];
258  }
259  }
260  }
261  }
262  // End : Wavefront Aberration Generation
263 
264 
265  for (int i = 0; i < (int)length_xn; i++)
266  {
267  memset(Temp_W[i], 0, length_yn * sizeof(double));
268  }
269 
270  int min_xnn, max_xnn;
271  int min_ynn, max_ynn;
272 
273  min_xnn = (int)round(length_xn / 2 - xr / 2);
274  max_xnn = (int)round(length_xn / 2 + xr / 2 + 1);
275  min_ynn = (int)round(length_yn / 2 - yr / 2);
276  max_ynn = (int)round(length_yn / 2 + yr / 2 + 1);
277 
278  int length_xnn, length_ynn;
279  length_xnn = max_xnn - min_xnn;
280  length_ynn = max_ynn - min_ynn;
281 
282  double **WT = new double*[length_xnn];
283  for (int i = 0; i < length_xnn; i++)
284  {
285  WT[i] = new double[length_ynn];
286  memset(WT[i], 0, length_ynn * sizeof(double));
287  }
288 
289  for (int i = 0; i < length_xnn; i++)
290  {
291  for (int j = 0; j < length_ynn; j++)
292  {
293  WT[i][j] = W[min_xnn+i][min_ynn+j];
294  }
295  }
296 
297  double **WS = new double*[(int)xr];
298  for (int i = 0; i < (int)xr; i++)
299  {
300  WS[i] = new double[yr];
301  memset(WS[i], 0, yr * sizeof(double));
302  }
303 
304  imresize(WT, length_xnn, length_ynn, xr, yr, WS);
305 
306 
307  oph::Complex<Real> **WD = new oph::Complex<Real>*[xr];
308 
309  for(int i = 0; i < (int)xr; i++)
310  WD[i] = new oph::Complex<Real>[yr];
311 
312  for(int ii = 0; ii < (int)xr; ii ++ )
313  {
314  for (int jj = 0; jj < (int)yr; jj++)
315  {
316 
317  WD[ii][jj]= exp(-j * (oph::Complex<Real>)2 * M_PI*WS[ii][jj] / wave_lambda); // Wave Aberration Complex Field
318  }
319  }
320  //WD[x][y]
321 
322  for (int i = 0; i < (int)length_xn; i++)
323  {
324  delete [] W[i];
325  delete [] Temp_W[i];
326  }
327  delete[] W;
328  delete[] Temp_W;
329 
330  for (int i = 0; i < (int)xr; i++)
331  {
332  delete[] WS[i];
333  }
334  delete[] WS;
335 
336  for (int i = 0; i < (int)length_xnn; i++)
337  {
338  delete[] WT[i];
339  }
340  delete[] WT;
341 
342  complex_W = WD;
343 
344  for (unsigned int x = 0; x < xr; x++) {
345  for (unsigned int y = 0; y < yr; y++) {
346  complex_H[0][x + y * xr] = complex_W[x][y];
347  }
348  }
349 
350 // return WD;
351 
352  auto end_time = CUR_TIME;
353 
354  auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
355 
356  LOG("Implement time : %.5lf sec\n", during_time);
357 }
358 
359 
361 {
362  for (int i = 0; i < (int)context_.pixel_number[_X]; i++)
363  {
364  delete[] doublePtr[i];
365  }
366 }
367 
369 {
370  this->Free2D(complex_W);
371  std::cout << " ophFree" << std::endl;
372 }
373 
374 
375 bool ophWaveAberration::readConfig(const char* fname)
376 {
377  LOG("Reading....%s...\n", fname);
378 
379 
380  /*XML parsing*/
381  tinyxml2::XMLDocument xml_doc;
382  tinyxml2::XMLNode *xml_node;
383  tinyxml2::XMLElement *xml_element;
384  const tinyxml2::XMLAttribute *xml_attribute;
385 
386 
387  if (!checkExtension(fname, ".xml"))
388  {
389  LOG("file's extension is not 'xml'\n");
390  return false;
391  }
392  auto ret = xml_doc.LoadFile(fname);
393  if (ret != tinyxml2::XML_SUCCESS)
394  {
395  LOG("Failed to load file \"%s\"\n", fname);
396  return false;
397  }
398 
399  xml_node = xml_doc.FirstChild();
400  xml_element = xml_node->FirstChildElement("Wavelength");
401  if (!xml_element || tinyxml2::XML_SUCCESS != xml_element->QueryDoubleText(&context_.wave_length[0]))
402  return false;
403 
404 
405  xml_element = xml_node->FirstChildElement("PixelPitchHor");
406  if (!xml_element || tinyxml2::XML_SUCCESS != xml_element->QueryDoubleText(&context_.pixel_pitch[_X]))
407  return false;
408 
409  xml_element = xml_node->FirstChildElement("PixelPitchVer");
410  if (!xml_element || tinyxml2::XML_SUCCESS != xml_element->QueryDoubleText(&context_.pixel_pitch[_Y]))
411  return false;
412 
413  xml_element = xml_node->FirstChildElement("ResolutionHor");
414  if (!xml_element || tinyxml2::XML_SUCCESS != xml_element->QueryIntText(&context_.pixel_number[_X]))
415  return false;
416 
417  xml_element = xml_node->FirstChildElement("ResolutionVer");
418  if (!xml_element || tinyxml2::XML_SUCCESS != xml_element->QueryIntText(&context_.pixel_number[_Y]))
419  return false;
420 
421  xml_element = xml_node->FirstChildElement("ZernikeCoeff");
422  xml_attribute = xml_element->FirstAttribute();
423 
424  for(int i=0; i< 45; i++)
425  {
426  if (!xml_attribute || tinyxml2::XML_SUCCESS != xml_attribute->QueryDoubleValue(&zernikeCoefficent[i]))
427  return false;
428  xml_attribute=xml_attribute->Next();
429 
430  }
431 
432  pixelPitchX = context_.pixel_pitch[_X];
433  pixelPitchY = context_.pixel_pitch[_Y];
434 
437 
438  waveLength = *context_.wave_length;
439 
443 
444  cout << "Wavelength: " << context_.wave_length[0] << endl;
445  cout << "PixelPitch(Horizontal): " << context_.pixel_pitch[_X] << endl;
446  cout << "PixelPitch(Vertical): " << context_.pixel_pitch[_Y] << endl;
447  cout << "Resolution(Horizontal): " << context_.pixel_number[_X] << endl;
448  cout << "Resolution(Vertical): " << context_.pixel_number[_Y] << endl;
449  cout << "Zernike Coefficient: " << endl;
450  for(int i=0; i<45; i++)
451  {
452  if (i!=0 && (i+1)%5 == 0)
453  cout << "z["<<i<<"]="<< zernikeCoefficent[i]<<endl;
454  else
455  cout << "z[" << i << "]=" << zernikeCoefficent[i] <<" ";
456  zernikeCoefficent[i] = zernikeCoefficent[i] * context_.wave_length[0];
457  }
458  int xr = context_.pixel_number[_X];
459  int yr = context_.pixel_number[_Y];
460 
461  complex_H[0] = new Complex<Real>[xr * yr];
462 
463  return true;
464 
465 }
466 
467 void ophWaveAberration::saveAberration(const char* fname)
468 {
469  ofstream fout(fname, ios_base::out | ios_base::binary);
470  fout.write((char *)complex_W, context_.pixel_number[_X] * context_.pixel_number[_Y] * sizeof(oph::Complex<Real>));
471  fout.close();
472 }
473 
474 void ophWaveAberration::readAberration(const char* fname)
475 {
476 
478  for (int i = 0; i < (int)context_.pixel_number[_X]; i++)
480 
481  ifstream fin(fname, ios_base::in | ios_base::binary);
483  fin.close();
484 }
485 
486 bool ophWaveAberration::loadAsOhc(const char * fname)
487 {
488  if (!Openholo::loadAsOhc(fname)) {
489  LOG("Failed load file");
490  return false;
491  }
492 
493  pixelPitchX = context_.pixel_pitch[_X];
494  pixelPitchY = context_.pixel_pitch[_Y];
495 
496  int xr = resolutionX = context_.pixel_number[_X];
497  int yr = resolutionY = context_.pixel_number[_Y];
498 
499  waveLength = context_.wave_length[0];
500  for (int x = 0; x < xr; x++) {
501  for (int y = 0; y < yr; y++) {
502  complex_W[x][y] = complex_H[0][x + y * xr];
503  }
504  }
505 
506  return true;
507 }
void abs(const oph::Complex< T > &src, oph::Complex< T > &dst)
Definition: function.h:113
Real factorial(double x)
Factorial using recursion.
ophWaveAberration()
Constructor.
#define M_PI
Definition: define.h:52
~ophWaveAberration()
Destructor.
void setPixelNumberOHC(const ivec2 pixel_number)
getter/setter for OHC file read and write
Definition: Openholo.h:501
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
void setPixelPitchOHC(const vec2 pixel_pitch)
Definition: Openholo.h:504
bool readConfig(const char *fname)
read configuration from a configration file
float Real
Definition: typedef.h:55
bool checkExtension(const char *fname, const char *ext)
Functions for extension checking.
Definition: Openholo.cpp:86
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
void setWavelengthOHC(const Real wavelength, const LenUnit wavelength_unit)
Definition: Openholo.h:507
vec2 pixel_pitch
Definition: Openholo.h:101
void readAberration(const char *fname)
reads the 2D data array of a wave aberration from a file
#define _X
Definition: define.h:92
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
ivec2 pixel_number
Definition: Openholo.h:99
virtual bool loadAsOhc(const char *fname)
Function to read OHC file.
Definition: Openholo.cpp:280
#define _Y
Definition: define.h:96
XMLError LoadFile(const char *filename)
Definition: tinyxml2.cpp:2150
XMLError QueryDoubleText(double *dval) const
See QueryIntText()
Definition: tinyxml2.cpp:1690
Complex< Real > ** complex_H
Definition: Openholo.h:489
OphConfig context_
Definition: Openholo.h:485
const XMLAttribute * Next() const
The next attribute in the list.
Definition: tinyxml2.h:1146
oph::Complex< Real > ** complex_W
double pointer of the 2D data array of a wave aberration
Real * wave_length
Definition: Openholo.h:106
const XMLElement * FirstChildElement(const char *name=0) const
Definition: tinyxml2.cpp:940
unsigned int uint
Definition: typedef.h:62
#define CUR_TIME
Definition: function.h:58
void imresize(double **X, int Nx, int Ny, int nx, int ny, double **Y)
Resizes 2D data using bilinear interpolation.
const XMLAttribute * FirstAttribute() const
Return the first attribute in the list.
Definition: tinyxml2.h:1471