Retinal Optical Coherence Tomography Image Denoising Using Modified Soft Thresholding Wavelet Transform

The block diagram of the implemented PSO-WT (Particle Swarm Optimization based modified soft thresholding Wavelet Transform) method is shown in Figure 2.

2.png

Figure 2. Block diagram of proposed PSO based modified soft thresholding for OCT image denoising

The noisy retinal OCT image (z) is input to the log transform block which convert multiplicative speckle noise into additive noise. By using log transforming the image we can use near additive models which are simpler to work [10]. So, followed the similar approaches as in studies [8, 10, 17, 19]. The 2D discrete wavelet transform (DWT) is taken to obtain approximation (low-frequency sub-band) and detailed coefficient (high-frequency sub-band data) representation of the input image data (z) of size N×N. The size of each coefficient matrix is (N/2×N/2), i and j are used for coefficient matrix indexing.

$\mathrm{y}=$ 2D_DWT(z) $=\left(y_{\mathrm{i} . \mathrm{j}}^{L L}, y_{\mathrm{i} . \mathrm{j}}^{L H}, y_{\mathrm{i} . \mathrm{j}}^{H L}, y_{\mathrm{i} . \mathrm{j}}^{H H}\right)$       (1)

The detailed wavelet coefficients are filtered using the modified threshold function given by Cao et al. [19]. It incorporates three hyper-parameters that decide the degree of curvature (β), degree of shrinkage of high frequency coefficient (α), and threshold value (K). The modified wavelet coefficients are obtained using Eqns. (2) and (3):

$\tilde{y}_{i, j=1}\left\{\begin{array}{cc}x_{i, j}+\frac{2}{\pi} \tan ^{-1}\left(\beta \frac{x_{i, j}}{T h}+\operatorname{sgn}\left(y_{i . j}\right) \alpha\right) y_{t h},\left|y_{i, j}\right| \geq y_{t h} \\ \frac{2}{\pi} \tan ^{-1}(\alpha) y_{i, j},|y i, j|<y_{t h}\end{array}\right.$        (2)

x_(i,j) is calculated by adding or subtracting the threshold value from the wavelet coefficient depending on whether y is greater than or less than y_th. The threshold value is calculated by:

$y_{t h}=K \frac{m d n\left(\left|H H_{10 c t}\,\,\right|\right)}{0.6754} \frac{\sqrt{2 \ln M}}{L_d}$       (3)

where, d is the decomposition level of the wavelet transform, L_ddepends on decomposition level. M is the length of the signal.

The quality of denoised images depends on the values set for the three hyper parameters in Eqns. (2) and (3). We have used Particle Swarm Optimization (PSO) algorithm to optimize hyper-parameter values with PSNR as the fitness function. PSO mimics navigation logic followed by a flock of birds. It is used to find the optimum solution while maximizing or minimizing function value. Particles represent a bird’s population [22].

$P_i=\left\{\alpha_i, \beta_i, K_i\right\}$        (4)

Every particle starts with random values of hyper-parameters in different directions. Particles move in the problem space following a set of rules. The image is reconstructed by inverse 2D DWT and inverse log transform for every particle and all the iterations. The fitness value i.e., PSNR of each particle solution is computed between input noisy image and reconstructed image.

$Fitness$ $=P S N R_{p i}$             (5)

If the particles fitness (Fitness_Pi) value is greater than the particle best (Pbest) so far then, particle best value is modified. In case the particles fitness (Fitness_Pi) value is greater than the global fitness (Gbest) value, the global value is modified. Then, the velocity and position (value) of the particles are updated for the next round of iteration.

For modified soft thresholding scheme each particle position (value) is comprises three parameters (α), (β) and k. Hence Eqns. (6) and (7) are used to modify the velocity and value of α for next iteration.

$\begin{array}{r}v_{-} \alpha_i^{i t r+1}=w v_{-} \alpha_i^{i t r}+c 1 r 1\left(\text {Pbest}_{\alpha_i}-\alpha_i^{i t r}\right) \\ +c 2 r 2\left(\text { Gbest_} \alpha_i-\alpha_i^{i t r}\right)\end{array}$        (6)

$\alpha_i^{i t r+1}=\alpha_i^{i t r}+v_{-} \alpha_i^{i t r 1+1}$           (7)

where, c1, c2 are acceleration coefficients, itr is the number of iterations the particles would undergo, and w is inertia weight. The r1 and r2 are random numbers used to set the initial values of tuning parameters and $v_{-} \alpha_i^{i t r+1}$ is the velocity for i^thα particle in next iteration.

The Eqns. (8) and (9) are used to modify the velocity and value of β:

$\begin{aligned} v_{\beta_i^{i t r+1}}\,\,=w v_{\beta_i^{i t r}} & +c 1 r 1\left(\text { Pbest}_{\beta_i}-\beta_i^{i t r}\right)  +c 2 r 2\left(\text { Gbest_} \beta_i-\beta_i^{i t r}\right)\end{aligned}$          (8)

$\beta_i^{i t r+1}=\beta_i^{i t r}+v_{-} \beta_i^{i t r 1+1}$       (9)

and Eqns. (10) and (11) for modifying k for next iteration.

$\begin{aligned} & v_{-} k_i^{i t r+1}=w v_{-} k_i^{i t r}+c 1 r 1\left(\text { Pbest}_{k_i}-k_i^{\mathrm{itr}}\right) \\ &+\mathrm{c} 2 \mathrm{r} 2\left(\text { Gbest_} \mathrm{k}_{\mathrm{i}}-\mathrm{k}_{\mathrm{i}}^{\mathrm{itr}}\right)\end{aligned}$      (10)

$k_i^{i t r+1}=k_i^{i t r}+v_{-} k_i^{i t r 1+1}$       (11)

The PSO algorithm is executed for several iterations. The final optimized solution for three hyper-parameter which maximizes the PSNR is obtained. The optimized values of hyper-parameters α, β, and k using PSO to give better performance may be different for different images, based on image quality, scaling, the orientation of image and intensity distribution, etc. Thus, the proposed system provides automation of hyperparameter tuning for modified soft thresholding wavelet denoising approach, irrespective of image the main steps of algorithm are stated in Algorithm.

Algorithm: Modified soft threshold wavelet transform with PSO for OCT image denoising

The performance of the proposed method for retinal OCT images is evaluated using parameters MSE, PSNR, SSIM, and CNR.

MSE (Mean squre Error): This is a basic parameter used to compare the original noisy image(z) and denoised image (d), and N total number of pixels.

$\mathrm{MSE}=\frac{1}{N} \sum_{j=1}^N\left(Z_j-D_J\right)^2$

PSNR (Peak signal to noise ratio): When SNR is expressed with reference to the maximum power of intensity within the image, it is referred to as PSNR.

$\mathrm{PSNR}=20 \log _{10} \frac{\left(\mathbf{2}^{\mathrm{B}}-\mathbf{1}\right)}{\sqrt{\mathrm{MSE}}}$

where, B is the number of bits per pixel and $2^B-1$ is the maximum value of gray level.

Contrast to Noise Ratio (CNR): This is a patch/window-based parameter. This gives a contrast relationship between background and image features. It specifies the ability to visualize the image features in a noisy environment:

$\mathrm{CNR}=\frac{1}{\mathrm{~F}} \sum_{\mathrm{f}=1}^{\mathrm{f}=\mathrm{F}} \frac{\left|\mu_{\mathrm{f}}-\mu_{\mathrm{b}}\right|}{\sqrt{\left|\sigma_{\mathrm{f}}^2-\sigma_{\mathrm{b}}^2\right|}}$

where, μ_b and $\sigma_{\mathrm{b}}^2$ are mean and the variance of background image, μ_f and $\sigma_f^2$ are mean and the variance of patches of foreground image.

Structural similarity index (SSIM): SSIM gives the structural relationship between the original image and the denoised image. This parameter justifies the human visual system more closely than PSNR. SSIM compares the two images based on luminance, contrast, and structure parameters.

$\operatorname{SSIM}=\frac{\left(2 \mathrm{u}_{\mathrm{z}} u_{\mathrm{d}}+\mathrm{b} 1\right)\left(2 \sigma_{\mathrm{zd}}\,+\,\mathrm{b} 2\right)}{\left(\mathrm{u}_{\mathrm{z}}^2\,+\mathrm{u}_{\mathrm{d}}^2\,+\mathrm{b} 1\right)\left(\sigma_{\mathrm{z}}^2\,+\sigma_{\mathrm{d}}^2+\mathrm{b} 2\right)}$

where, μ_z, σ_z, are the mean, standard deviation of noisy input image, μ_d and σ_d are the mean, standard deviation of denoised image(d). σ_zd is cross-covariance for input noisy image(z) and denoised image(d). b1 and b2 are constants calculated from dynamic range of pixxel values.

In implementation the SSIM function from MATLAB is used with all default parameters values except dynamic range and radius. Dynamic range is calculated from the denoised image, and the radius value is kept at 2.

The modified soft thresholding-based denoising algorithm preserves the finer details of images. Still, the performance depends on tuning parameters that affect artifacts, edge preservation in the reconstructed image, and filtering of speckles carrying information. It is observed that the use of PSO based tuning of parameters improves results marginally as compared to the same methods without PSO. The advantage here is that there is no need to manually find and set the hyper parameter values whenever quality of input image (due to Scaling, orientation of the image, intensity distribution, etc.) is changed.

The proposed method is implemented for denoising any OCT images of different characteristics and getting optimum PSNR (fitness function for PSO) value. In Future work, if any specific parameter is to be improved that can considered as fitness function and its effectiveness on overall performance can be tested. Also, effectiveness of the proposed denoising algorithm can be tested by segmentation or classification tasks using denoised OCT images for disease detection.