2.1.

Wave Atoms Transform

In this section, we briefly describe a recently developed half multiscale and multidirectional representation called the wave atoms transform.³⁸ Wave atoms can be seen as a variant of 2-D wavelet packets and obey the parabolic scaling of curvelet wavelength $\sim$ ${(\text{diameter})}^2$. Wave atoms can be adapted to arbitrary local directions of a pattern and also can sparsely represent anisotropic patterns aligned with the axes. Moreover, the warped oscillatory functions and oriented textures in the wave atoms have been proven to have a dramatically sparser expansion compared to some other fixed standard representations, such as Gabor filters, wavelets, and curvelets, which make the wave atoms transform a better option for image characteristic extraction. Wave atoms interpolate precisely between Gabor filters and directional wavelets. Moreover, wave atoms come either as an orthonormal basis or a tight frame of directional wave packets and are particularly well suited for representing oscillatory patterns and textures in the images.³⁹

Let the wave atoms be ${\phi }_{\mu }(x)$ with the subscript $\mu =(j,\mathit{m},\mathit{n})=(j,m_1,m_2,n_1,n_2)$, where the five quantities $j$, $m_1$, $m_2$, $n_1$, and $n_2$ all are integer values. A point ($x_{\mu }$, ${\omega }_{\mu }$) is indexed in phase space by

Here $M>0$, the oscillations within the envelope of a wave atom in $x$ have a wavelength $\sim 2^{-2j}$. Here, the subscript $j$ indexes different “dyadic coronae,” whereas the additional subscript $M$ labels the different wave numbers $w_{\mu }$ within each dyadic corona.

Let $g$ be a real-valued continuous function with support included in the domain of [$-7\pi /6$, $5\pi /6$] and such that for $|\omega |\le \pi /3$, $g{(\pi /2-\omega )}^2+g{(\pi /2+\omega )}^2=1$ and $g(-\pi /2-2\omega )=g(\pi /2+\omega )$. Define $\nu$ as the inverse Fourier transform (FT) $v(t)={(2\pi )}^{-1}\int g(\omega ){\mathrm{e}}^{i\omega t}\mathrm{d}\omega$ and

The FT of ${\psi }_m^0$ is given by³⁸

Here, ${\tau }_m={(-1)}^m$, ${\alpha }_m=(\pi /2)(m+0.5)$ and $\sum _m{|{\mathrm{\Psi }}_m^0(\omega )|}^2=1$. The function $g$ is an appropriate real-valued, $C^{\infty }$ bump function, compactly supported on an interval of length $2\pi$. Writing basis functions as ${\psi }_{m,n}^j(x)={\psi }_m^j(x-2^{-j}n)=2^{0.5j}{\psi }_m^0(2^{j}x-n)$, the coefficients of the transform can be obtained by

Assuming that the function $u$ is discretized at $x_{\varphi }=\varphi h$, $h=1/N$ and $\mathrm{\Phi }=1,2,\dots N$, then we have,

This algorithm can be implemented by the following steps: (1) fast Fourier transform (FFT) of $u(x_{\varphi })$, (2) wrap the product $\overline{{\mathrm{\Psi }}_m^j}\tilde{u}$ by periodicity inside the interval [$-2^j\pi$, $2^j\pi$] for each ($j$, $m$), and (3) perform an inverse FFT.

By individually taking products of 1-D wave packets, 2-D orthonormal basis functions with four bumps can be formed in frequency.^38,39 The 2-D extension can be formed by the products

Here, $H$ denotes the Hilbert transform operator. The combinations ${\phi }_{\mu }^{(1)}=0.5({\phi }_{\mu }^{+}+{\phi }_{\mu }^{-})$ and ${\phi }_{\mu }^{(2)}=0.5({\phi }_{\mu }^{+}-{\phi }_{\mu }^{-})$ from the wave atoms frame are denoted jointly as $\{{\phi }_{\mu }\}=\{{\phi }_{\mu }^{(1)},{\phi }_{\mu }^{(2)}\}$, and we have

The coefficients of 2-D wave atoms transform can be given by³⁸

2.2.

Wave Atoms Cycle Spinning

Orthogonal wave atoms transform is not translation invariant. This means that visual distortions, such as Gibbs-like phenomenon or ringing effects, will be produced in the image discontinuous point neighborhood area (e.g., edges and textures) due to lack of translation invariant in the translation processing. Fortunately, any denoiser can be turned into a translation invariant denoiser by performing cycle spinning.⁴¹ The denoiser is applied to several shifted copies of the image, then the resulting denoised image is shifted back to the original position and the results are averaged. As the wave atoms transform is not a translation invariant, this approach will result in different estimates of the original image with statistically different noises, which are reduced by linear averaging on all denoised results. Denoting the 2-D circular shift by $S_{i,j}$, the forward and inverse wave atoms transform by WA and ${\mathrm{WA}}^{-1}$, and the threshold operator by $T(.)$, therefore the cycle spinning will be performed as⁴¹

2.3.

Algorithm Implementation

As described in the above section, in order to inhibit visual distortions, such as ringing effects and Gibbs-like phenomenon, the wave atoms based image despeckling algorithm is proposed by combining cycle spinning technology in the process of a hard threshold denoising theme. The denoising procedure of wave atoms cyclic spinning shrinkage is similar to the steps in the traditional wavelet transform method and includes seven steps. The seven steps are: logarithm transformation, cycle spinning operation, forward wave atom transform, thresholding set, inverse wave atoms transform, reverse cycle spinning, and the exponential transformations to convert the denoising images from the logarithmic scale to the linear scale. The procedure for speckle reduction in OCT imaging based on wave atoms cycle spinning shrinkage is shown in Fig. 1.

Fig. 1

The basic work frame of speckle reduction for OCT images based on wave atoms transform.

The specific algorithm steps are as follows:

3.1.

Quality Metrics

For a quantitative comparison of the various thresholding filters, such as wavelet²⁹ and curvelet,³⁶ we used several blind speckle denoising metrics to quantify the despeckled image quality. The metrics included SNR, CNR, and ENL.²⁹ The CNR measures the contrast between an area of the image feature and an area of background noise, while the ENL measures the smoothness of homogeneous regions of interests (ROIs) that are corrupted by speckle noise. The ENL is a good metric for the quantification of noise reduction in homogeneous areas, and it increases with noise reduction. The ENL in an intensity image is a measurement of the statistical speckle fluctuations. Thus, ENL gives essentially an idea about the smoothness in regions of the images that are supposed to have a homogeneous appearance but are affected by noise.

The definitions for SNR, CNR, and ENL are defined as²⁹

Moreover, we also use the cross-correlation (XCOR) value as an indicator of similarity between different speckle fields, which measures the similarity between the images before and after denoising. Its definition is given by³⁶

3.2.

Results and Discussions

This section gives a detailed qualitative and quantitative analysis of the proposed OCT speckle reduction algorithm based on wave atoms transform. The image data are acquired with a spectrometer-based 890 nm Fourier domain OCT system.⁴² Briefly, the spectrometer-based Fourier domain OCT uses a super luminescent diode light source, which has a central wavelength of 890 nm and full width at half maximum bandwidth of 150 nm. A modified scanning head from a commercial Zeiss Stratus OCT was used. The optical power on the human eye was set at 650 μW. The charge-coupled device (CCD) integration time was set at 50 μs. The system sensitivity was measured to be about 100 dB at around zero imaging depth. The 6 dB sensitivity roll-off distance was found to be at an imaging depth of 1.6 mm. The imaging process includes background signal subtraction, linear interpolation to convert data from the linear wavelength space to the linear wavenumber space, and numerical dispersion. The lateral and axial resolutions were measured to be 20.0 and 3.5 μm, respectively, in the air.

We applied our algorithm to the acquired OCT image of a human retinal image with $512\times 1536\text{pixels}$ ($\mathrm{axial}\text{transverse}$). Here, the retinal image was first reduced to $512\times 512\text{pixels}$ by averaging groups of three adjacent axial scans, as shown in Fig. 2(a). As shown in Fig. 2(a), six ROIs, including three nonhomogeneous (labeled 1, 3, and 5) and three homogeneous (labeled 2, 4, and 6) areas, respectively, were selected in the retinal image and marked with solid rectangles in red. Meanwhile, a biological tissue-free area in the top region of the retinal image was selected as a “background” region, which was marked with a dashed rectangle in red and used for calculating SNR of the images. A region in the black rectangle window in Fig. 2(a) was enlarged 1.8 times and used for comparative analysis. Meanwhile, wavelet²⁹ and curvelet³⁵ based thresholding filterings were also applied to the same OCT image. For objective and equivalent comparisons, cycle spinning technology was also applied to remove Gibbs artifacts in the wavelet filtering and curvelet filtering processes. Note that we performed the wavelet- and curvelet-based thresholdings with the best results by obtaining the most appropriate threshold by a trial and error method in the image processing. The threshold is chosen to be 3.6 and 1.1 times of the noise variances, which are obtained from the median estimator of the highest sub-band of the transform for wavelet- and curvelet-based methods, respectively.^35,36 Comparing wavelet- and curvelet-based methods with the proposed wave atoms-based methods for a similar cross-correlation XCOR value, the threshold value in the wave atoms despeckling process is chosen to be $K=1.0$. In all cases, the images were processed with MATLAB on an Intel dual-core 2.2 GHz laptop.

Fig. 2

Denoised results of a human retinal image with $512\times 512\text{pixels}$. (a) Original image, (b) denoised image with wavelet, (c) curvelet, and (d) proposed wave atoms shrinkage, respectively. The areas in solid rectangles (marked in red) are the selected regions of interests (ROIs) and the dashed red rectangle area is the selected background region used for quantitative comparison of the performance of all speckle reduction algorithms based on wavelet, curvelet, and the proposed wave atom transforms applied to the original image. The red elliptical area points to structures with clearer visibility in the denoised OCT images; white scale bars represent 0.80 mm.

Figure 2 shows the original and despeckled OCT images of a human retina. The unprocessed image, shown in Fig. 2(a), has a grainy appearance due to the presence of speckle noise. Figure 2(b) shows the denoised image by the traditional wavelet method and the same OCT image after despeckling by curvelet is shown in Fig. 2(c). Figure 2(d) shows the denoised image using the proposed wave atoms transform method. It can be seen from the results that all three methods can remove most speckle noises and improve the visualization of small morphological features, such as outer plexiform layer (OPL) and retinal pigment epithelium, which are not shown clearly in the original raw OCT image. Meanwhile, it also can be seen that the processed OCT image, as shown in Fig. 2(d) with the proposed wave atoms method, shows much better visual effects and clearer detailed morphological features indicated, for example, by the red elliptical area in Fig. 2(d), than with the wavelet and curvelet methods. However, it needs to be pointed out that, like other image despeckling methods, the proposed wave atoms method also will blur or injure some feature details and edge information in the processed image. For example, there are some discontinuities in the OPL layer in the despeckled image compared to the original image. Therefore, it is necessary to keep a balance between obtaining the highest level of despeckling in the processed image and maintaining minimum damage to the original image in the practice applications.

For further comparative analysis, enlarged views of the region in the solid rectangle marked in black in Fig. 2(a) are presented in Fig. 3. It can be seen from Fig. 3 that there is almost no residual speckle noise pattern in the denoised images with wavelet, curvelet, and wave atoms shrinkages. However, there appears to be significant visual aberration (ringing effects and artificial effect in some areas in the despeckled image) and blurring of features in the denoised images by wavelet- and curvelet-domain filtering methods compared to the proposed wave atoms transform method. Namely, the proposed wave atoms shrinkage method shows much better visual effects and clearer detailed morphological features in the denoised images, as indicated, for example, by the two black arrows in the Fig. 3 (right image), the OPL and external limiting membrane (ELM). For the wavelet method, as shown in Fig. 3 (left), the ELM layer is almost invisible and the OPL layer cannot be resolved easily in the wavelet-based despeckled image, while in the curvelet-based despeckled image, as shown in Fig. 3 (middle), we can see that the ELM layer is distorted and discontinuous. Moreover, the OPL layer edge is distorted with mottling, and it cannot be resolved easily. In the wave atoms-based despeckled image, as shown in Fig. 3 (right), it shows a nicely continuous OPL layer, and the ELM layer is very clear except for where there are a few minor discontinuities as indicated by the white arrow.

Fig. 3

Enlarged copies of the solid black rectangle region in OCT retinal images [as shown in Fig. 2(a)] for close comparison of the performance of three image processing algorithms. The black arrows point to layer structures with clearer visibility and the white arrow points to the discontinuous layer structure in the denoised OCT images. The white scale bars represent 0.42 mm.

The reason for this performance difference is that the wavelet transform is only optimal in representing 1-D singularities but not optimal in representing 2-D image contours or layers because of poor orientation selectivity and anisotropy. Therefore, it cannot effectively represent images consisting of different regions separated by boundaries that often appear in OCT images. The curvelet transform gives an optimal sparse representation because it obeys a more precise parabolic balance between oscillations and support size, and it concentrates the energy of an object with an arbitrary discontinuity curve in just a few coefficients.³⁴ Meanwhile, wave atoms can be adapted to arbitrary local directions of a pattern and also can sparsely represent anisotropic patterns aligned with the axes. Moreover, the warped oscillatory functions and oriented textures in wave atoms have been proven to have a dramatically sparser expansion compared to wavelets and curvelets, and wave atoms transform is particularly well suited for representing oscillatory patterns and textures in the images³⁹ As a result, the wave atoms transform would generate relatively larger transformed coefficients for such continuous and weak layer features (like ELM), while the wavelet and curvelet coefficients for the weak ELM signals are small and are easily attenuated along with speckle noise. The proposed wave atoms method, therefore, is a highly promising preprocessing approach that can enhance retinal OCT image quality for further quantitative analysis, such as segmentation for retina layer measurement and analysis.

To quantitatively analyze the performance of the proposed algorithm, we also calculated the image quality metrics, such as SNR, CNR, ENL, and XOCR. The CNR values were averaged over six ROIs ($m=1\to 6$) within the solid red rectangle in Fig. 2(a). Similarly, the ENL values were averaged over the three ROIs, $m=2$, 4, 6, in Fig. 2(a). The selected background region in the dashed rectangle was used to calculate the background noise level. Note that the CNR and ENL were computed using a logarithmic scale while SNR and XCOR were calculated using a linear scale. Table 1 gives the results of the qualitative metrics for the original and the filtered images. It clearly shows that the denoised results with curvelet and the proposed wave atoms transforms are better than the traditional wavelet thresholding method. Note that we performed the wavelet- and curvelet-based thresholdings with the best results by obtaining the most appropriate threshold by a trial and error method in the image processing. For similar XCOR values, for example, we took the XCOR values of 0.925, 0.924, and 0.926 for the wavelet-, curvelet-, and wave atoms-based methods, respectively. As shown in Table 1, our proposed wave atoms shrinkage method ($K=1.0$) made further improvements to the SNR, 8.43 dB and 3.12 dB, compared to the wavelet- and curvelet-based methods, which also greatly improving the image quality of the CNR and ENL.

Table 1

Image quality metrics for OCT images.

In speckle reduction with wave atoms shrinkage, the parameter $K$ controls the degree of speckle reduction, and there is a trade-off between speckle reduction degree and edge preservation in practical applications. Figure 4 shows the trend of SNR improvement and cross-correlation XCOR with a different threshold factor, $K$, in the denoising process. As shown in Fig. 4, the XCOR decreases with the increment of threshold factor $K$. However, SNR increases with the increment of $K$ until $K$ reaches 1.1. The XCOR decreases sharply initially at small $K$ values, and then it does not change significantly for $K$ between 0.8 and 1.2. Increasing the threshold further would eventually lead to the loss of image features and decrease the XCOR. Moreover, Fig. 4 also shows acceptable and stable despeckled results (SNR and XCOR) for a wide available range of threshold $K$ (from 0.7 to 1.2) in the wave atoms despeckling process.

Fig. 4

Measurement of SNR improvement and cross-correlation XCOR as a function of different threshold $K$ values (from 0 to 1.2) in an OCT retinal image by wave atoms despeckling algorithm. The algorithm improves the maximum SNR of 58.70 dB at $K=1.1$, while the cross-correlation XCOR between the original and despeckled images is 0.923.