Openholo  v1.1
Open Source Digital Holographic Library
Openholo.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 "Openholo.h"
47 
48 #include <windows.h>
49 #include <fileapi.h>
50 #include <omp.h>
51 #include "sys.h"
52 #include "ImgCodecOhc.h"
53 #include "ImgControl.h"
54 
56  : Base()
57  , plan_fwd(nullptr)
58  , plan_bwd(nullptr)
59  , fft_in(nullptr)
60  , fft_out(nullptr)
61  , pnx(1)
62  , pny(1)
63  , pnz(1)
64  , fft_sign(OPH_FORWARD)
65  , OHC_encoder(nullptr)
66  , OHC_decoder(nullptr)
67  , complex_H(nullptr)
68 {
69  context_ = { 0 };
70  fftw_init_threads();
71  fftw_plan_with_nthreads(omp_get_max_threads());
74 }
75 
77 {
78  if (OHC_encoder) {
79  delete OHC_encoder;
80  OHC_encoder = nullptr;
81  }
82  if (OHC_decoder) {
83  delete OHC_decoder;
84  OHC_decoder = nullptr;
85  }
86  fftw_cleanup_threads();
87 }
88 
89 bool Openholo::checkExtension(const char * fname, const char * ext)
90 {
91  string filename(fname);
92  string fext(ext);
93  size_t find = filename.find_last_of(".");
94  size_t len = filename.length();
95 
96  if (find > len)
97  return false;
98 
99  if (!filename.substr(find).compare(fext))
100  return true;
101  else
102  return false;
103 }
104 
105 bool Openholo::saveAsImg(const char * fname, uint8_t bitsperpixel, uchar* src, int width, int height)
106 {
107  LOG("Saving...%s...\n", fname);
108  bool bOK = true;
109  auto start = CUR_TIME;
110  int _width = width, _height = height;
111 
112  int padding = 0;
113  int _byteperline = ((_width * bitsperpixel / 8) + 3) & ~3;
114  int _pixelbytesize = _height * _byteperline;
115  int _filesize = _pixelbytesize;
116  bool hasColorTable = (bitsperpixel <= 8) ? true : false;
117  int _headersize = sizeof(bitmap);
118  int _iColor = (hasColorTable) ? 256 : 0;
119 
120  if (_width % 4 != 0) {
121  padding = 4 - _width % 4;
122  }
123 
124  rgbquad *table = nullptr;
125 
126  if (hasColorTable) {
127  _headersize += _iColor * sizeof(rgbquad);
128  table = new rgbquad[_iColor];
129  memset(table, 0, sizeof(rgbquad) * _iColor);
130  for (int i = 0; i < _iColor; i++) { // for gray-scale
131  table[i].rgbBlue = i;
132  table[i].rgbGreen = i;
133  table[i].rgbRed = i;
134  }
135  }
136 
137  _filesize += _headersize;
138 
139  bool bConvert = _stricmp(PathFindExtensionA(fname) + 1, "bmp") ? true : false;
140 
141  uchar *pBitmap = new uchar[_filesize];
142  memset(pBitmap, 0x00, _filesize);
143 
144  bitmap bitmap;
145  memset(&bitmap, 0, sizeof(bitmap));
146  int iCur = 0;
147 
148  bitmap.fileheader.signature[0] = 'B';
149  bitmap.fileheader.signature[1] = 'M';
150  bitmap.fileheader.filesize = _filesize;
152 
154  bitmap.bitmapinfoheader.width = _width;
155  bitmap.bitmapinfoheader.height = _height;
157  bitmap.bitmapinfoheader.bitsperpixel = bitsperpixel;
159  bitmap.bitmapinfoheader.imagesize = _pixelbytesize;
160  bitmap.bitmapinfoheader.ypixelpermeter = 0;// Y_PIXEL_PER_METER;
161  bitmap.bitmapinfoheader.xpixelpermeter = 0;// X_PIXEL_PER_METER;
163 
164  memcpy(&pBitmap[iCur], &bitmap.fileheader, sizeof(fileheader));
165  iCur += sizeof(fileheader);
166  memcpy(&pBitmap[iCur], &bitmap.bitmapinfoheader, sizeof(bitmapinfoheader));
167  iCur += sizeof(bitmapinfoheader);
168 
169  if (hasColorTable) {
170  memcpy(&pBitmap[iCur], table, sizeof(rgbquad) * _iColor);
171  iCur += sizeof(rgbquad) * _iColor;
172  }
173 
174 
175 
176  if (context_.bRotation) {
177  ImgControl *pControl = ImgControl::getInstance();
178  uchar *pTmp = new uchar[_pixelbytesize];
179 
180  pControl->Rotate(180.0, src, pTmp, _width + padding, _height, _width + padding, _height, bitsperpixel / 8);
181  memcpy(&pBitmap[iCur], pTmp, _pixelbytesize);
182  delete[] pTmp;
183  iCur += _pixelbytesize;
184  }
185  else {
186  if (padding != 0)
187  {
188  for (int i = 0; i < _height; i++)
189  {
190  memcpy(&pBitmap[iCur], &src[_width * i], _width);
191  iCur += _width;
192  memset(&pBitmap[iCur], 0x00, padding);
193  iCur += padding;
194  }
195  }
196  else
197  {
198  memcpy(&pBitmap[iCur], src, _pixelbytesize);
199  iCur += _pixelbytesize;
200  }
201  }
202 
203 
204 
205  if (!bConvert) {
206  if (iCur != _filesize)
207  bOK = false;
208  else {
209  FILE *fp;
210  fopen_s(&fp, fname, "wb");
211  if (fp == nullptr)
212  bOK = false;
213  else {
214  fwrite(pBitmap, 1, _filesize, fp);
215  fclose(fp);
216  }
217  }
218  }
219  else {
220  ImgControl *pControl = ImgControl::getInstance();
221  pControl->Save(fname, pBitmap, _filesize);
222  }
223 
224  if(hasColorTable && table) delete[] table;
225  delete[] pBitmap;
226  auto end = CUR_TIME;
227 
228  auto during = ((std::chrono::duration<Real>)(end - start)).count();
229 
230  LOG("%.5lfsec...done\n", during);
231 
232  return bOK;
233 }
234 
235 
236 uchar * Openholo::loadAsImg(const char * fname)
237 {
238  FILE *infile;
239  fopen_s(&infile, fname, "rb");
240  if (infile == nullptr) { LOG("No such file"); return 0; }
241 
242  // BMP Header Information
243  fileheader hf;
244  bitmapinfoheader hInfo;
245  fread(&hf, sizeof(fileheader), 1, infile);
246  if (hf.signature[0] != 'B' || hf.signature[1] != 'M') { LOG("Not BMP File"); return 0; }
247 
248  fread(&hInfo, sizeof(bitmapinfoheader), 1, infile);
249  fseek(infile, hf.fileoffset_to_pixelarray, SEEK_SET);
250 
251  oph::uchar *img_tmp;
252  if (hInfo.imagesize == 0) {
253  int nSize = (((hInfo.width * hInfo.bitsperpixel / 8) + 3) & ~3) * hInfo.height;
254  img_tmp = new uchar[nSize];
255  fread(img_tmp, sizeof(oph::uchar), nSize, infile);
256  }
257  else {
258  img_tmp = new uchar[hInfo.imagesize];
259  fread(img_tmp, sizeof(oph::uchar), hInfo.imagesize, infile);
260  }
261  fclose(infile);
262 
263  return img_tmp;
264 }
265 
266 bool Openholo::saveAsOhc(const char * fname)
267 {
268  std::string fullname = fname;
269  if (!checkExtension(fname, ".ohc")) fullname.append(".ohc");
270  OHC_encoder->setFileName(fullname.c_str());
271 
272  // Clear vector
274 
275  ohcHeader header;
276  OHC_encoder->getOHCheader(header);
277  auto wavelength_num = header.fieldInfo.wavlenNum;
278 
279  for (uint i = 0; i < wavelength_num; i++)
281 
282  if (!OHC_encoder->save()) return false;
283 
284  return true;
285 }
286 
287 bool Openholo::loadAsOhc(const char * fname)
288 {
289  std::string fullname = fname;
290  if (!checkExtension(fname, ".ohc")) fullname.append(".ohc");
291  OHC_decoder->setFileName(fullname.c_str());
292  if (!OHC_decoder->load()) return false;
293 
297 
298  vector<Real> wavelengthArray;
299  OHC_decoder->getWavelength(wavelengthArray);
300  context_.wave_length = new Real[wavelengthArray.size()];
301  for (int i = 0; i < wavelengthArray.size(); i++)
302  context_.wave_length[i] = wavelengthArray[i];
303 
305 
306  context_.k = (2 * M_PI) / context_.wave_length[0];
309 
310  return true;
311 }
312 
313 bool Openholo::loadAsImgUpSideDown(const char * fname, uchar* dst)
314 {
315  FILE *infile;
316  fopen_s(&infile, fname, "rb");
317  if (infile == nullptr) { LOG("No such file"); return false; }
318 
319  // BMP Header Information
320  fileheader hf;
321  bitmapinfoheader hInfo;
322  fread(&hf, sizeof(fileheader), 1, infile);
323  if (hf.signature[0] != 'B' || hf.signature[1] != 'M') { LOG("Not BMP File"); return false; }
324 
325  fread(&hInfo, sizeof(bitmapinfoheader), 1, infile);
326  fseek(infile, hf.fileoffset_to_pixelarray, SEEK_SET);
327 
328  oph::uchar* img_tmp;
329  if (hInfo.imagesize == 0) {
330  img_tmp = new oph::uchar[hInfo.width*hInfo.height*(hInfo.bitsperpixel / 8)];
331  fread(img_tmp, sizeof(oph::uchar), hInfo.width*hInfo.height*(hInfo.bitsperpixel / 8), infile);
332  }
333  else {
334  img_tmp = new oph::uchar[hInfo.imagesize];
335  fread(img_tmp, sizeof(oph::uchar), hInfo.imagesize, infile);
336  }
337  fclose(infile);
338 
339  // data upside down
340  int bytesperpixel = hInfo.bitsperpixel / 8;
341  //int rowsz = bytesperpixel * hInfo.width;
342 
343  int rowsz = ((hInfo.width * bytesperpixel) + 3) & ~3;
344 
345  for (oph::uint k = 0; k < hInfo.height*rowsz; k++) {
346  int r = k / rowsz;
347  int c = k % rowsz;
348  ((oph::uchar*)dst)[(hInfo.height - r - 1)*rowsz + c] = img_tmp[r*rowsz + c];
349  }
350 
351  delete[] img_tmp;
352  return true;
353 }
354 
355 bool Openholo::getImgSize(int & w, int & h, int & bytesperpixel, const char * fname)
356 {
357 
358  char szExtension[_MAX_EXT] = { 0, };
359  sprintf(szExtension, "%s", PathFindExtension(fname) + 1);
360 
361  if (_stricmp(szExtension, "bmp")) { // not bmp
362  return false;
363  }
364 
365  // BMP
366  FILE *infile;
367  fopen_s(&infile, fname, "rb");
368  if (infile == NULL) { LOG("No Image File"); return false; }
369 
370  // BMP Header Information
371  fileheader hf;
372  bitmapinfoheader hInfo;
373  fread(&hf, sizeof(fileheader), 1, infile);
374  if (hf.signature[0] != 'B' || hf.signature[1] != 'M') return false;
375  fread(&hInfo, sizeof(bitmapinfoheader), 1, infile);
376  //if (hInfo.bitsperpixel != 8) { printf("Bad File Format!!"); return 0; }
377 
378  w = hInfo.width;
379  h = hInfo.height;
380  bytesperpixel = hInfo.bitsperpixel / 8;
381  fclose(infile);
382 
383  return true;
384 }
385 
386 void Openholo::ImageRotation(double rotate, uchar* src, uchar* dst, int w, int h, int channels)
387 {
388  auto begin = CUR_TIME;
389  int channel = channels;
390  int nBytePerLine = ((w * channel) + 3) & ~3;
391 
392  int origX, origY;
393  double radian = rotate * M_PI / 180.0;
394  double cc = cos(radian);
395  double ss = sin(-radian);
396  double centerX = (double)w / 2.0;
397  double centerY = (double)h / 2.0;
398 
399  uchar pixel;
400  uchar R, G, B;
401 
402  int num_threads = 1;
403  for (int y = 0; y < h; y++) {
404  for (int x = 0; x < w; x++) {
405  origX = (int)(centerX + ((double)y - centerY)*ss + ((double)x - centerX)*cc);
406  origY = (int)(centerY + ((double)y - centerY)*cc - ((double)x - centerX)*ss);
407 
408  pixel = 0;
409  R = G = B = 0;
410  if ((origY >= 0 && origY < h) && (origX >= 0 && origX < w)) {
411  int offsetX = origX * channel;
412  int offsetY = origY * nBytePerLine;
413  B = src[offsetY + offsetX + 0];
414  G = src[offsetY + offsetX + 1];
415  R = src[offsetY + offsetX + 2];
416 
417  }
418  if (channel == 1) {
419  dst[y * nBytePerLine + x] = pixel;
420  }
421  else if (channel == 3) {
422  int tmpX = x * 3;
423  int tmpY = y * nBytePerLine;
424  dst[tmpY + tmpX + 0] = B;
425  dst[tmpY + tmpX + 1] = G;
426  dst[tmpY + tmpX + 2] = R;
427  }
428  }
429  }
430  auto end = CUR_TIME;
431  LOG("Image Rotated (%d threads): (%d/%d) (%lf degree) : %lf(s)\n",
432  num_threads,
433  w, h, rotate,
434  ((chrono::duration<Real>)(end - begin)).count());
435 }
436 
437 void Openholo::imgScaleBilinear(uchar* src, uchar* dst, int w, int h, int neww, int newh, int channels)
438 {
439  auto begin = CUR_TIME;
440  int channel = channels;
441  int nBytePerLine = ((w * channel) + 3) & ~3;
442  int nNewBytePerLine = ((neww * channel) + 3) & ~3;
443  int num_threads = 1;
444 #ifdef _OPENMP
445  int y;
446 #pragma omp parallel
447  {
448  num_threads = omp_get_num_threads();
449 #pragma omp for private(y)
450  for (y = 0; y < newh; y++)
451 #else
452  for (int y = 0; y < newh; y++)
453 #endif
454  {
455  int nbppY = y * nNewBytePerLine;
456  for (int x = 0; x < neww; x++)
457  {
458  float gx = (x / (float)neww) * (w - 1);
459  float gy = (y / (float)newh) * (h - 1);
460 
461  int gxi = (int)gx;
462  int gyi = (int)gy;
463 
464  if (channel == 1) {
465  uint32_t a00, a01, a10, a11;
466 
467  a00 = src[gxi + 0 + gyi * nBytePerLine];
468  a01 = src[gxi + 1 + gyi * nBytePerLine];
469  a10 = src[gxi + 0 + (gyi + 1) * nBytePerLine];
470  a11 = src[gxi + 1 + (gyi + 1) * nBytePerLine];
471 
472  float dx = gx - gxi;
473  float dy = gy - gyi;
474 
475  float w1 = (1 - dx) * (1 - dy);
476  float w2 = dx * (1 - dy);
477  float w3 = (1 - dx) * dy;
478  float w4 = dx * dy;
479 
480  dst[x + y * neww] = int(a00 * w1 + a01 * w2 + a10 * w3 + a11 * w4);
481  }
482  else if (channel == 3) {
483  uint32_t b00[3], b01[3], b10[3], b11[3];
484  int srcX = gxi * channel;
485  int dstX = x * channel;
486 
487  b00[0] = src[srcX + 0 + gyi * nBytePerLine];
488  b00[1] = src[srcX + 1 + gyi * nBytePerLine];
489  b00[2] = src[srcX + 2 + gyi * nBytePerLine];
490 
491  b01[0] = src[srcX + 3 + gyi * nBytePerLine];
492  b01[1] = src[srcX + 4 + gyi * nBytePerLine];
493  b01[2] = src[srcX + 5 + gyi * nBytePerLine];
494 
495  b10[0] = src[srcX + 0 + (gyi + 1) * nBytePerLine];
496  b10[1] = src[srcX + 1 + (gyi + 1) * nBytePerLine];
497  b10[2] = src[srcX + 2 + (gyi + 1) * nBytePerLine];
498 
499  b11[0] = src[srcX + 3 + (gyi + 1) * nBytePerLine];
500  b11[1] = src[srcX + 4 + (gyi + 1) * nBytePerLine];
501  b11[2] = src[srcX + 5 + (gyi + 1) * nBytePerLine];
502 
503  float dx = gx - gxi;
504  float dy = gy - gyi;
505 
506  float w1 = (1 - dx) * (1 - dy);
507  float w2 = dx * (1 - dy);
508  float w3 = (1 - dx) * dy;
509  float w4 = dx * dy;
510 
511  dst[dstX + 0 + nbppY] = int(b00[0] * w1 + b01[0] * w2 + b10[0] * w3 + b11[0] * w4);
512  dst[dstX + 1 + nbppY] = int(b00[1] * w1 + b01[1] * w2 + b10[1] * w3 + b11[1] * w4);
513  dst[dstX + 2 + nbppY] = int(b00[2] * w1 + b01[2] * w2 + b10[2] * w3 + b11[2] * w4);
514  }
515  }
516  }
517 #ifdef _OPENMP
518  }
519 #endif
520  auto end = CUR_TIME;
521  LOG("Scaled img size (%d threads): (%d/%d) => (%d/%d) : %lf(s)\n",
522  num_threads,
523  w, h, neww, newh,
524  ((chrono::duration<Real>)(end - begin)).count());
525 }
526 
527 void Openholo::convertToFormatGray8(unsigned char * src, unsigned char * dst, int w, int h, int bytesperpixel)
528 {
529  int idx = 0;
530  unsigned int r = 0, g = 0, b = 0;
531  for (int i = 0; i < w*h*bytesperpixel; i++)
532  {
533  unsigned int blue = src[i + 0];
534  unsigned int green = src[i + 1];
535  unsigned int red = src[i + 2];
536  dst[idx++] = (blue + green + red) / 3;
537  i += bytesperpixel - 1;
538  }
539 }
540 
541 void Openholo::fft1(int n, Complex<Real>* in, int sign, uint flag)
542 {
543  pnx = n;
544  bool bIn = true;
545 
546  if (fft_in == nullptr)
547  fft_in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * n);
548  if (fft_out == nullptr)
549  fft_out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * n);
550 
551  if (!in){
552  in = new Complex<Real>[pnx];
553  bIn = false;
554  }
555 
556  for (int i = 0; i < n; i++) {
557  fft_in[i][_RE] = in[i].real();
558  fft_in[i][_IM] = in[i].imag();
559  }
560 
561  fft_sign = sign;
562 
563  if (!bIn) delete[] in;
564 
565  if (sign == OPH_FORWARD)
566  plan_fwd = fftw_plan_dft_1d(n, fft_in, fft_out, sign, flag);
567  else if (sign == OPH_BACKWARD)
568  plan_bwd = fftw_plan_dft_1d(n, fft_in, fft_out, sign, flag);
569  else {
570  LOG("failed fftw : wrong sign");
571  fftFree();
572  return;
573  }
574 }
575 
576 
577 void Openholo::fft2(oph::ivec2 n, Complex<Real>* in, int sign, uint flag)
578 {
579  pnx = n[_X], pny = n[_Y];
580  int N = pnx * pny;
581 
582  if (fft_in == nullptr)
583  fft_in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);
584  if (fft_out == nullptr)
585  fft_out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);
586 
587  if (in != nullptr)
588  {
589  int i;
590  for (i = 0; i < N; i++) {
591  fft_in[i][_RE] = in[i][_RE];
592  fft_in[i][_IM] = in[i][_IM];
593  }
594  }
595 
596  fft_sign = sign;
597 
598  if (sign == OPH_FORWARD)
599  plan_fwd = fftw_plan_dft_2d(pny, pnx, fft_in, fft_out, sign, flag);
600  else if (sign == OPH_BACKWARD)
601  plan_bwd = fftw_plan_dft_2d(pny, pnx, fft_in, fft_out, sign, flag);
602  else {
603  LOG("failed fftw : wrong sign");
604  fftFree();
605  return;
606  }
607 }
608 
609 void Openholo::fft3(oph::ivec3 n, Complex<Real>* in, int sign, uint flag)
610 {
611  pnx = n[_X], pny = n[_Y], pnz = n[_Z];
612  bool bIn = true;
613  if (fft_in == nullptr)
614  fft_in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * pnx * pny * pnz);
615  if (fft_out == nullptr)
616  fft_out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * pnx * pny * pnz);
617 
618  if (!in) {
619  in = new Complex<Real>[pnx * pny * pnz];
620  bIn = false;
621  }
622 
623  for (int i = 0; i < pnx * pny * pnz; i++) {
624  fft_in[i][_RE] = in[i].real();
625  fft_in[i][_IM] = in[i].imag();
626  }
627 
628  fft_sign = sign;
629 
630  if (!bIn) delete[] in;
631 
632  if (sign == OPH_FORWARD)
633  plan_fwd = fftw_plan_dft_3d(pnz, pny, pnx, fft_in, fft_out, sign, flag);
634  else if (sign == OPH_BACKWARD)
635  plan_bwd = fftw_plan_dft_3d(pnz, pny, pnx, fft_in, fft_out, sign, flag);
636  else {
637  LOG("failed fftw : wrong sign");
638  fftFree();
639  return;
640  }
641 }
642 
643 void Openholo::fftExecute(Complex<Real>* out, bool bReverse)
644 {
645  if (fft_sign == OPH_FORWARD)
646  fftw_execute(plan_fwd);
647  else if (fft_sign == OPH_BACKWARD)
648  fftw_execute(plan_bwd);
649  else {
650  LOG("failed fftw : wrong sign");
651  out = nullptr;
652  fftFree();
653  return;
654  }
655 
656  if (!bReverse) {
657  int i;
658 #ifdef _OPENMP
659 #pragma omp for private(i)
660 #endif
661  for (i = 0; i < pnx * pny * pnz; i++) {
662  out[i][_RE] = fft_out[i][_RE];
663  out[i][_IM] = fft_out[i][_IM];
664  }
665  }
666  else {
667  int div = pnx * pny * pnz;
668  int i;
669 #ifdef _OPENMP
670 #pragma omp for private(i)
671 #endif
672  for (i = 0; i < pnx * pny * pnz; i++) {
673  out[i][_RE] = fft_out[i][_RE] / div;
674  out[i][_IM] = fft_out[i][_IM] / div;
675  }
676 
677  }
678 
679  fftFree();
680 }
681 
683 {
684  if (plan_fwd) {
685  fftw_destroy_plan(plan_fwd);
686  plan_fwd = nullptr;
687  }
688  if (plan_bwd) {
689  fftw_destroy_plan(plan_bwd);
690  plan_fwd = nullptr;
691  }
692  fftw_free(fft_in);
693  fftw_free(fft_out);
694 
695  plan_bwd = nullptr;
696  fft_in = nullptr;
697  fft_out = nullptr;
698 
699  pnx = 1;
700  pny = 1;
701  pnz = 1;
702 }
703 
704 void Openholo::fftwShift(Complex<Real>* src, Complex<Real>* dst, int nx, int ny, int type, bool bNormalized)
705 {
706  const int N = nx * ny;
707 
708  Complex<Real>* tmp = new Complex<Real>[N];
709  memset(tmp, 0., sizeof(Complex<Real>) * N);
710  fftShift(nx, ny, src, tmp);
711 
712  fftw_complex *in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);
713  fftw_complex *out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);
714 
715  memcpy(in, tmp, sizeof(Complex<Real>) * N);
716 
717  fftw_plan plan = nullptr;
718  if (!plan_fwd && !plan_bwd) {
719  plan = fftw_plan_dft_2d(ny, nx, in, out, type, OPH_ESTIMATE);
720  fftw_execute(plan);
721  }
722  else {
723  if (type == OPH_FORWARD)
724  fftw_execute_dft(plan_fwd, in, out);
725  else if (type == OPH_BACKWARD)
726  fftw_execute_dft(plan_bwd, in, out);
727  }
728 
729  int normalF = 1;
730  if (bNormalized) normalF = N;
731  memset(tmp, 0, sizeof(Complex<Real>) * N);
732 
733  int k;
734 #pragma omp parallel for private(k) firstprivate(normalF)
735  for (k = 0; k < N; k++) {
736  tmp[k][_RE] = out[k][_RE] / normalF;
737  tmp[k][_IM] = out[k][_IM] / normalF;
738  }
739  fftw_free(in);
740  fftw_free(out);
741  if (plan)
742  fftw_destroy_plan(plan);
743 
744  memset(dst, 0, sizeof(Complex<Real>) * N);
745  fftShift(nx, ny, tmp, dst);
746  delete[] tmp;
747 }
748 
749 void Openholo::fftShift(int nx, int ny, Complex<Real>* input, Complex<Real>* output)
750 {
751  int hnx = nx / 2;
752  int hny = ny / 2;
753 
754  int i;
755 #ifdef _OPENMP
756 #pragma omp parallel for private(i) firstprivate(hnx, hny)
757 #endif
758  for (i = 0; i < nx; i++)
759  {
760  for (int j = 0; j < ny; j++)
761  {
762  int ti = i - hnx; if (ti < 0) ti += nx;
763  int tj = j - hny; if (tj < 0) tj += ny;
764 
765  output[ti + tj * nx] = input[i + j * nx];
766  }
767  }
768 
769 }
770 
771 void Openholo::setWaveNum(int nNum)
772 {
773  context_.waveNum = nNum;
774  if (context_.wave_length != nullptr) {
775  delete[] context_.wave_length;
776  context_.wave_length = nullptr;
777  }
778 
779  context_.wave_length = new Real[nNum];
780 }
781 
782 
784 {
785  ohcHeader header;
786  OHC_encoder->getOHCheader(header);
787  auto wavelength_num = header.fieldInfo.wavlenNum;
788  for (uint i = 0; i < wavelength_num; i++) {
789  if (complex_H[i]) {
790  delete[] complex_H[i];
791  complex_H[i] = nullptr;
792  }
793  }
794  if (complex_H) {
795  delete[] complex_H;
796  complex_H = nullptr;
797  }
798  if (context_.wave_length) {
799  delete[] context_.wave_length;
800  context_.wave_length = nullptr;
801  }
802  if (OHC_encoder) {
803  delete OHC_encoder;
804  OHC_encoder = nullptr;
805  }
806  if (OHC_decoder) {
807  delete OHC_decoder;
808  OHC_decoder = nullptr;
809  }
810 }
811 
fileheader fileheader
Definition: struct.h:76
uint32_t dibheadersize
Definition: struct.h:57
virtual bool saveAsImg(const char *fname, uint8_t bitsperpixel, uchar *src, int width, int height)
Function for creating image files.
Definition: Openholo.cpp:105
Real k
Definition: Openholo.h:67
#define _IM
Definition: complex.h:57
Definition: Base.h:55
void fftFree(void)
Definition: Openholo.cpp:682
bool Rotate(double rotate, unsigned char *src, unsigned char *dst, int w, int h, int neww, int newh, int ch)
Definition: ImgControl.cpp:195
#define OPH_ESTIMATE
Definition: define.h:76
unsigned char uchar
Definition: typedef.h:64
void fftShift(int nx, int ny, Complex< Real > *input, Complex< Real > *output)
Swap the top-left quadrant of data with the bottom-right , and the top-right quadrant with the bottom...
Definition: Openholo.cpp:749
ohcFieldInfoHeader fieldInfo
void ImageRotation(double rotate, uchar *src, uchar *dst, int w, int h, int channels)
Definition: Openholo.cpp:386
adapt ASM r
void fft3(ivec3 n, Complex< Real > *in, int sign=OPH_FORWARD, uint flag=OPH_ESTIMATE)
Functions for performing fftw 3-dimension operations inside Openholo.
vec2 ss
Definition: Openholo.h:68
#define _Y
Definition: define.h:84
bool checkExtension(const char *fname, const char *ext)
Functions for extension checking.
Definition: Openholo.cpp:89
void getOHCheader(ohcHeader &_Header)
return true
Definition: ophGen.cpp:844
#define _RE
Definition: complex.h:54
vec2 pixel_pitch
Definition: Openholo.h:65
uint32_t imagesize
Definition: struct.h:63
def k(wvl)
Definition: Depthmap.py:16
void getComplexFieldData(OphComplexField &cmplx_field, uint wavelen_idx)
Definition: ImgCodecOhc.h:81
#define OPH_BACKWARD
Definition: define.h:67
void setWaveNum(int nNum)
Definition: Openholo.cpp:771
uint8_t signature[2]
Definition: struct.h:51
Definition: struct.h:69
Definition: struct.h:75
virtual void releaseFldData()
void addComplexFieldData(const OphComplexField &data)
void imgScaleBilinear(uchar *src, uchar *dst, int w, int h, int neww, int newh, int channels=1)
Function for change image size.
Definition: Openholo.cpp:437
#define CUR_TIME
Definition: function.h:58
bool getImgSize(int &w, int &h, int &bytesperpixel, const char *fname)
Function for getting the image size.
Definition: Openholo.cpp:355
uint16_t bitsperpixel
Definition: struct.h:61
void fftExecute(Complex< Real > *out, bool bReverse=false)
Execution functions to be called after fft1, fft2, and fft3.
Definition: Openholo.cpp:643
void fft2(ivec2 n, Complex< Real > *in, int sign=OPH_FORWARD, uint flag=OPH_ESTIMATE)
Functions for performing fftw 2-dimension operations inside Openholo.
void convertToFormatGray8(uchar *src, uchar *dst, int w, int h, int bytesperpixel)
Function for convert image format to gray8.
Definition: Openholo.cpp:527
ss
Definition: 2D-RS.py:25
ImgDecoderOhc * OHC_decoder
Definition: Openholo.h:314
virtual ~Openholo(void)=0
Destructor.
Definition: Openholo.cpp:76
uint waveNum
Definition: Openholo.h:69
bitmapinfoheader bitmapinfoheader
Definition: struct.h:77
uint32_t ypixelpermeter
Definition: struct.h:64
fclose('all')
uint32_t numcolorspallette
Definition: struct.h:66
bool setFileName(const std::string &_fname)
Definition: ImgCodecOhc.cpp:92
ivec2 pixel_number
Definition: Openholo.h:64
parameters h
bool Save(const char *path, BYTE *pBuf, UINT len, int quality=100)
Definition: ImgControl.cpp:22
#define M_PI
Definition: define.h:52
virtual bool loadAsOhc(const char *fname)
Function to read OHC file.
Definition: Openholo.cpp:287
uint32_t filesize
Definition: struct.h:52
Openholo(void)
Constructor.
Definition: Openholo.cpp:55
uint32_t fileoffset_to_pixelarray
Definition: struct.h:54
uint32_t width
Definition: struct.h:58
virtual uchar * loadAsImg(const char *fname)
Function for loading image files.
Definition: Openholo.cpp:236
#define OPH_PLANES
Definition: define.h:145
uint32_t height
Definition: struct.h:59
virtual bool saveAsOhc(const char *fname)
Function to write OHC file
Definition: Openholo.cpp:266
uint16_t planes
Definition: struct.h:60
#define _Z
Definition: define.h:88
Complex< Real > ** complex_H
Definition: Openholo.h:307
float Real
Definition: typedef.h:55
void fftwShift(Complex< Real > *src, Complex< Real > *dst, int nx, int ny, int type, bool bNormalized=false)
Convert data from the spatial domain to the frequency domain using 2D FFT on CPU. ...
Definition: Openholo.cpp:704
#define _X
Definition: define.h:80
ImgEncoderOhc * OHC_encoder
OHC file format Variables for read and write.
Definition: Openholo.h:313
OphConfig context_
Definition: Openholo.h:306
#define OPH_COMPRESSION
Definition: define.h:146
#define OPH_FORWARD
Definition: define.h:66
Real * wave_length
Definition: Openholo.h:70
void fft1(int n, Complex< Real > *in, int sign=OPH_FORWARD, uint flag=OPH_ESTIMATE)
Functions for performing fftw 1-dimension operations inside Openholo.
Definition: Openholo.cpp:541
bool loadAsImgUpSideDown(const char *fname, uchar *dst)
Function for loading image files | Output image data upside down.
Definition: Openholo.cpp:313
void getWavelength(std::vector< double_t > &wavlen_array)
uint32_t compression
Definition: struct.h:62
int w
Definition: 2D-RS.py:16
unsigned int uint
Definition: typedef.h:62
virtual void ophFree(void)
Pure virtual function for override in child classes.
Definition: Openholo.cpp:783
uint32_t xpixelpermeter
Definition: struct.h:65
bool bRotation
Definition: Openholo.h:71