The intensity of an interferogram can be expressed as
\begin{equation} I(\vec{x})=I_0(\vec{x})+I_a(\vec{x})cos[\phi(\vec{x})] \end{equation}
where $(\vec{x})=(x,y)$ is the 2-D coordinates of the individual pixel in the interferogram, $I_0(\vec{x})$ is the background or mean intensity, $I_a(\vec{x})$ is the modulation amplitude, and $\phi(\vec{x})$ is the angular phase information.
The wavelet transform of an interferogram $I(\vec{x})$ is defined as
\begin{equation} W(\vec{b},a,\theta)=a^{-n} \int_{-\infty}^{\infty} I(\vec{x}) \cdot \psi^{*}\left[a^{-1} \cdot r_{-\theta}(\vec{x}-\vec{b})\right]d^2\vec{x} \end{equation}
where $a>0$ is a scale dilation parameter corresponding to the width of the wavelet, $\vec{b}$ is a 2-D translation parameter corresponding to the position of the wavelet, $\theta$ is a rotation parameter, $r_{-\theta}$ is the conventional 2x2 rotation operator matrix, $\psi (\vec{x})$ is the mother wavelet, $\psi^{*}$ denotes the complex conjugate of $\psi$, $n$ is the normalization parameter, and $W(\vec{b},a,\theta)$ is the wavelet transform coefficient.
The 2-D Morlet wavelet can be chosen and modified for interferogram processing. The modified wavelet is expressed as:
\begin{equation}
\psi(\vec{x})=exp\left(-m|\vec{x}|^2+i2\pi x\right)
\end{equation}
where $m$ is a parameter to be determined as described below.
It has been shown that $m$ and $n$ should be set to 2.
It has been proved that the phase at point $\vec{b}$ can be calculated from
\begin{equation}
\phi(\vec{b}) = \Phi = tan^{-1}\left(\frac{Im\left[W(\vec{b},a,\theta)\right]}{Re\left[W(\vec{b},a,\theta)\right]-I_0\frac{\pi}{m} exp\left(-\frac{\pi^2}{m}\right)}\right)
\end{equation}
where $Im$ and $Re$ denote imaginary and real parts of a complex value, respectively. $I_0$ can be set to the mean intensity of the interferogram.
Wavelet transform calculation involves a convolution of two functions. The time-consuming convolution operation can be implemented using the fast Fourier transforms (FFT) described as follows
\begin{equation}
W(\vec{b},a,\theta) = a^{-n}F^{-1}\left\{F\left[I(\vec{x})\right] \cdot F\left[\zeta(\vec{x})\right]\right\}
\end{equation}
where ``$\ast$" denotes convolution and ``$\cdot$" denotes dot product, $F$ and $F^{-1}$ denote the forward and inverse fast Fourier transforms, respectively.
The wavelet transform algorithm can be summarized as follows: for any arbitrary pair $a$ and $\theta$, the wavelet transform coefficients at all points in the interferogram can be calculated by using the aforementioned equation. The pair $(a,\theta)$ that yield the maximum wavelet transform coefficient at each point denote the local fringe period (frequency is $1/a$) and orientation at the same point, and the phase can also be calculated. It should be noted that there is no available theoretical guidance on determining the sampling intervals of parameters $a$ and $\theta$. In practice, a trade-off between finer parameter intervals and shorter computation times needs to be considered.
Once the phase information, fringe frequency and orientation are obtained, they can be directly employed for further interferogram analysis and processing such as determination of fringe order gradients by
\begin{equation}
g_x=a^{-1}cos(\theta), g_y=a^{-1}sin(\theta)
\end{equation}
The phase distributions can also be utilized to filter interferograms through reconstructing fringes. For a 255 gray level interferogram, the filtered interferogram can be obtained, after normalizing background intensities and modulation amplitudes, from
\begin{equation}
I^\prime(\vec{x})=127.5[1+cos\phi(\vec{x})]
\end{equation}
Hi Dr. Wang, Can you comment
Hi Dr. Wang,
Can you comment on the case of a heterodyne frequency. When I(x,y) = I0(x,y) + Ia(x,y) Cos[2πfx + 2πfy + Φ(x,y)]; as in Tomassini's Wavelet transform paper. What would be the equation to recover the phase?
Having Trouble implementing this code
Hi.
I am having trouble coding this. Maybe there is an error in the way I am handling complx numbers or there is a problem in defining x and y. In the past I have been using 1-D CWT methods. When I try to code 2-D CWT, the Absolute Value of W looks like the fringes. In 1-D CWT, abs (W) does not look like fringes. Can you post the code for Wavelet Transform Method or look at my code below?
Thanks.
My IDL code:
ii = dcomplex(0., 1.)
isize = size(image1, /dimensions)
x = dindgen(isize[0]) # replicate(1., isize[1])
y = dindgen(isize[1]) ## replicate(1., isize[0])
SumSq = (x^2. + y^2.)
Fimage = FFT(image1, /double)
BestW = dcomplexarr(isize[0], isize[1])
m = 2.
OneArr = Replicate(DComplex(1,1),isize[0], isize[1])
For t = -!Pi/4, !Pi/4, !Pi/16. Do Begin
For a = 25, 25 Do Begin
Delta = Exp(-m*a^(-2.)*SumSq + ii*2*!Pi*a^(-1.)*(x*Cos(t) + y*Sin(t)))
W = a^(-2.)*FFT(Fimage * FFT(Delta, /double), /Inverse, /double)
Mask = abs(abs(W)) GT abs(abs(BestW))
Mask = Mask + ii*Mask
BestW = BestW * (OneArr - Mask) + W * Mask
EndFor
EndFor
Problem Solved
Sorry, I was using the wrong definition of x and y. The zeros of the x and y grids should be centered in the array.