Openholo  beta 2.0 Open Source Digital Holographic Library
Get Parameter-AT

## Functions

double ophSig::sigGetParamAT ()
Extraction of distance parameter using axis transfomation . More...

## Detailed Description

This module is related method which extraction of distance parameter using axis transformation method.

# Introduction

• OSH is reconstructed digitally by convolving the complex conjugate of Fresnel zone plate (FZP) with the hologram. This is the same as conventional digital holography, in which the FZP’s distance parameter is set according to the depth location of the objects. However, since the depth location of the objects is unknown, digital reconstructions with different distance parameters are required until we get a focused image.
• Unfortunately, this searching process is time consuming and requires a manual process. Recently several numerical techniques that extract the distance parameter directly from the hologram without reconstructions have been proposed [1-3]; however, these involve a search algorithm [2,3] or a tracking process [1]. Most recently, an auto-focusing technique based on the Wigner distribution analysis of a hologram has been proposed [4].
• In the proposed technique, we extract the distance parameter directly from the hologram without any searching or tracking process. However, a manual process that measures the slope of the Wigner distribution output is required in order to determine the distance parameter. Therefore, we propose to extract the distance parameter directly from the hologram using axis transformation without any manual processes.

# Algorithm

• First, the complex hologram is filtered by a Gaussian low-pass filter with transfer function. The Gaussian low-pass filtered hologram is given by the following equation:

$H_{com}^{lp}\left(x,y\right)=F^{-1}\left\{F\left[H_{com}\left(x,y\right)\right]\times G\left(k_x,k_y\right)\right\} =\int_{z_0-\frac{1}{2}\delta z}^{z_0+\frac{1}{2}\delta z}I_o\left(x,y,z \right)\otimes\frac{j}{\lambda z} \exp\left[\left(\frac{-\pi}{a_{lp} z^2}-j\frac{\pi}{\lambda z}\right)\left(x^2+y^2\right)\right]dz \qquad (1) G\left(k_x,k_y\right)=exp\left\{-\pi\left[\frac{\lambda}{2\pi{NA}_g}\right]^2\left({k_x}^2+{k_y}^2\right)\right\}$

• Where $$F{.}$$ represents the Fourier transform operator, $$a_{lp}\left(z\right)=NA_gNA/\sqrt{NA^2+NA_g^2}\times z$$ determines the radius of the Gaussian low pass filtered hologram. Hence, the Reyleigh range of the Gaussian low-pass filtered hologram evolves into the following equation:

$\Delta z=2{\lambda z}^2/{{(a}_{lp}z}^2\ \pi)=2\lambda/\pi\times({NA}^2+{NA}_g^2)/{({NA}_gNA)}^2$

• In Gaussian low pass filtering, we set $$NA_g$$ such that the Rayleigh range of the FZP is larger than the depth range of the object, i.e. $$\Delta z\geq\delta z$$. The radius of the scanning beam pattern is approximately constant within the depth range of the object. This makes the phase of the Gaussian low pass filtered hologram stationary within the depth range of the specimen [4, 5]; thus, we can extract the FZP which contains the information only about the distance parameter from the hologram.
• Second, a real-only spectrum hologram in the frequency domain is synthesized and is given by the following:

$H_{r-only}^{lp}\left(k_x,k_y\right)=Re\left\{F\left[Re\left(H_{com}^{lp}\left(x,y\right)\right)\right]\right\} +jRe\left\{F\left[Im\left(H_{com}^{lp}\left(x,y\right)\right)\right]\right\} \qquad (2)$

• After projecting the real-only spectrum hologram onto the ky direction, we filter the square of the projected real-only spectrum hologram by a power-fringe-adjusted filter [6-9]. The power-fringe-adjusted filtered output is given by:

$H_{FZP}\left(k_x\right)=\frac{\left[\int{H_{r-only}^{lp}\left(k_x,k_y\right)dk_y}\right]^2}{\left|\int{H_{r-only}^{lp} \left(k_x,k_y\right)dk_y}\right|^2+\delta}\approx exp\left(j\frac{\lambda z_0}{2\pi}k_x^2\right) \qquad (3)$

• Eq. (3), the filtered output is a chirping signal whose chirping rate is determined by the distance parameter. Note that the filtered output is the one-dimensional FZP with the distance parameter.
• Third, the real part of Eq. (3) is extracted and is given by the following:\ .\

$H_{FZP}^{Re}\left(k_x\right)=Re\left\lfloor H_{FZP}\left(k_x\right)\right\rfloor \approx cos\left(\frac{\lambda z_0}{2\pi}k_x^2\right) \qquad (4)$

• $$H_{FZP}^{Re}\left(k_x\right)$$ shown as fig.1-a.
• Fourth, the transformation from original frequency axis to new-frequency axis using interpolation.

$k_x^{new}=k_x^2\ , H_{FZP}^{Re}\left(k_x^{new}\right)\approx cos\left(\frac{\lambda z_0}{2\pi}k_x^{new}\right) \qquad (5)$

• Note that this sinusoidal signal has a single frequency and the frequency of the signal is directly proportional to the distance parameter. Hence, the inverse Fourier transformation of Eq.(5) expresses the delta function pair in the new spatial axis:

$h_{FZP}^{Re}\left(x^{new}\right) =\mathbf{F}^{-\mathbf{1}}\left\{H_{FZP}^{Re}\left(k_x^{new}\right)\right\} \approx\frac{1}{2}\delta\left(x^{new}-\frac{\lambda z_0}{2\pi}\right)+\frac{1}{2}\delta\left(x^{new}+\frac{\lambda z_0}{2\pi}\right) \qquad (6)$

• Note that the location of the delta function pair gives the distance parameter. This can be extracted directly by detecting the location of the maximum value of Eq. (6).
Figure 1. Flowchart.

# Reference

• [1] P. Ferraro, G. Coppola, S. D. Nicola, A. Finizio, and G. Pierattini, “Digital holographic microscope with automatic focus tracking by detecting sample displacement in real time,” Opt. Lett. 28, 1257-1259 (2003).
• [2] M. Liebling and M. Unser, “Autofocus for digital Fresnel holograms by use of a Fresnelet-sparsity criterion,” J. Opt. Soc. Am. A 21, 2424-2430 (2004).
• [3] P. Langehanenberg, B. Kemper, D. Dirksen, and G. von Bally, “Autofocusing in digital holographic phase contrast microscopy on pure phase objects for live cell imaging,” Appl. Opt. 47, D176 (2008).
• [4] T. Kim and T.-C. Poon, “Auto-focusing in optical scanning holography,” Appl. Opt. 48, H153-H159 (2009).
• [5] T. Kim, Y. S. Kim, W. S. Kim, and T.-C. Poon, “Algorithm for converting full-parallax holograms to horizontal parallax -only holograms,” Opt. Lett. 34, 1231-1233 (2009).
• [6] T. Kim and T.-C. Poon, “Extraction of 3-D location of matched 3-D object using power fringe-adjusted filtering and Wigner analysis,” Opt. Eng. 38, 2176-2183 (1999).
• [7] T. Kim and T.-C. Poon, “Experiments of depth detection and image recovery of a remote target using a complex hologram,” Opt. Eng. 43, 1851-1855 (2004).
• [8] T. Kim, T.-C. Poon, and G. Indebetouw, ‘‘Depth detection and image recovery in remote sensing by optical scanning holography,’’ Opt. Eng. 41, 1331-1338 (2002).
• [9] P. Klysubun, G. Indebetouw, T. Kim, and T.-C. Poon, ‘‘Accuracy of three-dimensional remote target location using scanning holographic correlation,’’ Opt. Comm. 184, 357 -366 (2000).

## ◆ sigGetParamAT()

 double ophSig::sigGetParamAT ( )

Extraction of distance parameter using axis transfomation .

Returns
result distance

Definition at line 891 of file ophSig.cpp.

891  {
892
893  int i = 0, j = 0;
894  Real max = 0; double index = 0;
895  float NA_g = (float)0.025;
896  int nx = context_.pixel_number[_X];
897  int ny = context_.pixel_number[_Y];
898
899  OphComplexField Flr(nx, ny);
900  OphComplexField Fli(nx, ny);
901  OphComplexField Hsyn(nx, ny);
902  OphComplexField Hsyn_copy1(nx, ny);
903  OphComplexField Hsyn_copy2(nx, ny);
904  OphRealField Hsyn_copy3(nx, ny);
905
906  OphComplexField Fo(nx, ny);
907  OphComplexField Fon, yn, Ab_yn;
908
909  OphRealField Ab_yn_half;
910  OphRealField G(nx, ny);
911  Real r = 1, c = 1;
912  vector<Real> t, tn;
913
914  for (i = 0; i < nx; i++)
915  {
916  for (j = 0; j < ny; j++)
917  {
918
919  r = (2 * M_PI*(i) / _height - M_PI*(nx - 1) / _height);
920  c = (2 * M_PI*(j) / _width - M_PI*(ny - 1) / _width);
921  G(i, j) = std::exp(-M_PI * (context_.wave_length[0] / (2 * M_PI * NA_g)) * (context_.wave_length[0] / (2 * M_PI * NA_g)) * (c * c + r * r));
922  Flr(i, j)._Val[0] = ComplexH[0](i, j)._Val[0];
923  Fli(i, j)._Val[0] = ComplexH[0](i, j)._Val[1];
924  Flr(i, j)._Val[1] = 0;
925  Fli(i, j)._Val[1] = 0;
926  }
927  }
928
929  fft2(Flr, Flr);
930  fft2(Fli, Fli);
931
932  int xshift = nx / 2;
933  int yshift = ny / 2;
934
935  for (i = 0; i < nx; i++)
936  {
937  int ii = (i + xshift) % nx;
938  for (j = 0; j < ny; j++)
939  {
940  Hsyn(i, j)._Val[_RE] = Flr(i, j)._Val[_RE] * G(i, j);
941  Hsyn(i, j)._Val[_IM] = Fli(i, j)._Val[_RE] * G(i, j);
942  Hsyn_copy1(i, j) = Hsyn(i, j);
943  Hsyn_copy2(i, j) = Hsyn_copy1(i, j) * Hsyn(i, j);
944  Hsyn_copy3(i, j) = pow(sqrt(Hsyn(i, j)._Val[_RE] * Hsyn(i, j)._Val[_RE] + Hsyn(i, j)._Val[_IM] * Hsyn(i, j)._Val[_IM]), 2) + pow(10, -300);
945  int jj = (j + yshift) % ny;
946  Fo(ii, jj)._Val[_RE] = Hsyn_copy2(i, j)._Val[0] / Hsyn_copy3(i, j);
947  Fo(ii, jj)._Val[_IM] = Hsyn_copy2(i, j)._Val[1] / Hsyn_copy3(i, j);
948  }
949  }
950
951  t = linspace(0, 1, nx / 2 + 1);
952  tn.resize(t.size());
953  Fon.resize(1, t.size());
954
955  for (int i = 0; i < tn.size(); i++)
956  {
957  tn.at(i) = pow(t.at(i), 0.5);
958  Fon(0, i)._Val[0] = Fo(nx / 2 - 1, nx / 2 - 1 + i)._Val[0];
959  Fon(0, i)._Val[1] = 0;
960  }
961
962  yn.resize(1, static_cast<int>(tn.size()));
963  linInterp(t, Fon, tn, yn);
964  fft1(yn, yn);
965  Ab_yn.resize(yn.size[_X], yn.size[_Y]);
966  absMat(yn, Ab_yn);
967  Ab_yn_half.resize(1, nx / 4 + 1);
968
969  for (int i = 0; i < nx / 4 + 1; i++)
970  {
971  Ab_yn_half(0, i) = Ab_yn(0, nx / 4 + i - 1)._Val[_RE];
972  if (i == 0) max = Ab_yn_half(0, 0);
973  else
974  {
975  if (Ab_yn_half(0, i) > max)
976  {
977  max = Ab_yn_half(0, i);
978  index = i;
979  }
980  }
981  }
982  index = -(((index + 1) - 120) / 10) / 140 + 0.1;
983
984  return index;
985 }
void fft1(matrix< Complex< T >> &src, matrix< Complex< T >> &dst, int sign=OPH_FORWARD, uint flag=OPH_ESTIMATE)
Function for Fast Fourier transform 1D.
Definition: ophSig.cpp:190
void fft2(matrix< Complex< T >> &src, matrix< Complex< T >> &dst, int sign=OPH_FORWARD, uint flag=OPH_ESTIMATE)
Function for Fast Fourier transform 2D.
Definition: ophSig.cpp:226
vector< Real > linspace(double first, double last, int len)
Generate linearly spaced vector.
Definition: ophSig.cpp:272
void absMat(matrix< Complex< T >> &src, matrix< T > &dst)
Function for extracts Complex magnitude value.
Definition: ophSig.cpp:59
void linInterp(vector< T > &X, matrix< Complex< T >> &in, vector< T > &Xq, matrix< Complex< T >> &out)
Function for Linear interpolation 1D.
Definition: ophSig.cpp:295