OCT Image Restoration Using Non-Local Deep Image Prior

Abstract

In recent years, convolutional neural networks (CNN) have been widely used in image denoising for their high performance. One difficulty in applying the CNN to medical image denoising such as speckle reduction in the optical coherence tomography (OCT) image is that a large amount of high-quality data is required for training, which is an inherent limitation for OCT despeckling. Recently, deep image prior (DIP) networks have been proposed for image restoration without pre-training since the CNN structures have the intrinsic ability to capture the low-level statistics of a single image. However, the DIP has difficulty finding a good balance between maintaining details and suppressing speckle noise. Inspired by DIP, in this paper, a sorted non-local statics which measures the signal autocorrelation in the differences between the constructed image and the input image is proposed for OCT image restoration. By adding the sorted non-local statics as a regularization loss in the DIP learning, more low-level image statistics are captured by CNN networks in the process of OCT image restoration. The experimental results demonstrate the superior performance of the proposed method over other state-of-the-art despeckling methods, in terms of objective metrics and visual quality.

1. Introduction

Medical images are an important application of digital imaging combined with biology to visually express important medical information. Optical coherence tomography (OCT) is a rapidly evolving imaging model with high resolution, dynamic range, and non-invasive and non-destructive detection, suitable for detecting a variety of biological tissues [1,2]. However, OCT is susceptible to speckle noise during acquisition, resulting in poor imaging results and preventing the physician from making accurate diagnoses [3]. Therefore, OCT image denoising is very important for subsequent medical applications.

In the past few decades, a variety of methods have been developed for reducing the speckle noise in OCT images. These algorithms can be categorized as transformation-based, spatial filtering, and sparse representations. In the transform domain, there are mainly the wavelet transform [4,5], the curvelet transform [6], and the contourlet transform [7]. In the transform domain, the coefficients representing the image are assumed to be distinguishable from the speckle noise by using shrinkage or threshold. However, when the decomposed sub-bands come to recomposition, the transform domain algorithms are prone to producing some artifacts. Spatial filtering techniques are mainly represented by non-local means (NLM) methods [8,9,10,11,12], which have the advantage of the self-similarity of natural images by comparing patches with non-local neighborhoods [8]. Based on the double Gaussian anisotropic kernels [9], the uncorrupted probability of each pixel (PNLM) [10], or the weighted maximum likelihood estimation [11], the non-local means is able to reduce the speckle noise more effectively. Non-local means-based methods take advantage of the self-similarity of images in a non-local manner; they are usually more robust compared with local-based filtering. However, these techniques tend to over-smooth the image details and leave a few artifacts in the despeckling results. Sparse representation has been widely applied to signal processing, such as multi-focus image fusion [13] and OCT image despeckling [14,15]. Sparse representation techniques in OCT adopt a customized scanning pattern in which a few B-scans are captured with high signal-to-noise ratio (SNR) which can be used to improve the quality of low-SNR B-scans. Based on the high-SNR images, multi-scale sparsity-based tomographic denoising (MSBTD) [14] trains the sparse representation dictionaries to improve the quality of neighboring low-SNR B-scans. Furthermore, in sparsity-based simultaneous denoising and interpolation (SBSDI) [15], the dictionaries are improved by constructing from previously collected datasets high-SNR images from the target imaging subject. These sparse representation techniques can preserve most image details, but they tend to leave some noise in the despeckling result.

In recent years, with the advent of deep convolutional neural networks, the ability to exploit the spatial correlations and extract the features at multiple resolutions has been greatly improved by using a hierarchical network structure. Current advancements in neural networks have shown their great applicability in supervised and unsupervised signal preprocessing and classification. Many phases of the biosignal process have been augmented with the use of deep neural network-based methods, such as electron transport proteins identifying [16,17,18] and OCT despeckling [19,20]. Deep convolutional neural networks have become a popular tool for image denoising [21,22,23]. The denoising CNN (DnCNN) [21] method showed that the residual learning and batch normalization were particularly useful for image denoising. The denoising performance of CNN can be improved by symmetric convolutional-deconvolutional layers [22] and the prior observation model [23]. However, to train their numerous parameters, these networks generally require large amounts of high-quality data, which is an inherent limitation for OCT despeckling. To tackle this issue, frame averaging has been used for training the network [19,20], but the results were not satisfactory since the ground truth images were obtained only by frame averaging and still contained speckle noise.

Deep image prior (DIP) [24] shows that the structure of the generator network is sufficient to capture a large number of low-level image statistics—that is, before any learning has shown that CNN can produce a denoised image without a “clean” dataset or a training dataset. The main reason is that compared with the noise, the convolution layer can easily recover the signal which has more self-similarity. However, as the reconstruction loss of DIP only measures the L1 or L2 loss between the network output and the input image, it is not suitable for multiple speckle noise suppressions.

In this paper, a non-local deep image prior network is proposed for OCT image restoration by adding a sorted non-local statics as autocorrelation loss in the reconstruction function. The sorted non-local statics measures the signal autocorrelation by calculating the correlation between each patch with its non-local neighbors in the differences between the constructed image and the input image, and then sorting these correlation values to select the most correlated patches. By minimizing the signal autocorrelation loss in the DIP learning, more non-local similarity image statistics are captured by CNN in the process of OCT image restoration. It allows the recovery of most of the OCT signals.

In summary, the two main contributions in this paper are as follows:

(1) A sorted non-local statics is proposed which measures the signal autocorrelation by calculating and sorting the correlation between each patch with its non-local neighbors in the differences between the constructed image and the input image.

(2) The sorted non-local statics is used as an autocorrelation loss in the deep image prior learning framework to get rid of the speckle noise in the OCT image.

2. Non-Local Deep Image Prior

2.1. Deep Image Prior Model

In this section, we briefly review the deep image prior model [24].

Deep convolutional networks are one of the most commonly-used techniques in image denoising. Traditional deep convolutional networks need to train a large amount of image information, and the training parameters of the deep convolutional network can be determined by:

$$\widehat{\theta }=\underset{\theta }{\arg \min }{\displaystyle \sum _i||x_{\mathrm{label}}^i}-f(\theta |z^i)||$$

(1)

where $\widehat{\theta }$ represents the training parameters, $x_{\mathrm{label}}^i$ represents the $i$th original image, $z^i$ represents the $i$th input noisy image, and $f$ represents the deep convolutional network. However, it is difficult to prepare a large amount of high-quality data for the OCT image despeckling. The DIP model indicates that a large amount of training data is not needed, and the deep convolution network itself has the ability to learn the low-dimensional information of noise images. In the DIP [24], the goal of the image restoration process is to recover the original image $\widehat{x}$ from the corrupted image $x$. Such problems are often formulated as an optimization task:

$$\widehat{x}=\underset{x}{\arg \min }(E(\widehat{x};x)+R(\widehat{x}))$$

(2)

where $E(\widehat{x};x)$ is a task-dependent data item and $R(\widehat{x})$ is a regularizer. Similar to the adaptive dictionary learning, $R(\widehat{x})$ can be replaced by the implicit prior captured by the deep convolution network, so the DIP model can be written as:

$$\widehat{\theta }=\underset{\theta }{\arg \min }E(f(\theta |z);x)$$

(3)

where the network input $z$ of the DIP network is random noise. The parameters $\widehat{\theta }$ which are represented $\widehat{x}$ can be trained by using an optimizer such as gradient descent. In the DIP denoising model, $E(\widehat{x};x)$ is the reconstruction loss, which measures the error between the constructed image and the input image:

$$E(\widehat{x};x)=\left|\right|\widehat{x}-x|{|}^2.$$

(4)

Therefore, Equation (3) becomes:

$$\widehat{\theta }=\underset{\theta }{\arg \min }||f(\theta |z)-x|{|}^2.$$

(5)

Figure 1 is the OCT despeckling results of the DIP. Figure 1a is the original noisy OCT image, and Figure 1b–d is the despeckling results after 500, 1000 and 1500 iterations, respectively. It can be seen that when the DIP is trained to reconstruct the OCT image, it generates along the way an intermediate clean version of the image, before overfitting the noise. It can be seen from Figure 1c that the clean version of the deep image prior allows the recovery of most of the image details while containing some speckle noise. In other words, based on the reconstruction loss in Equation (4), the DIP model has difficulty finding a good balance between maintaining details and suppressing speckle noise.

2.2. Non-Local Deep Image Prior Model

The convolutional neural networks (CNNs) have the intrinsic ability to solve inverse problems such as denoising without any pre-training. The main reason is that compared with the noise, the convolution layer is better able to recover the signal which has more self-similarity. However, for OCT despeckling, the DIP model has difficulty finding a good balance between maintaining details and suppressing speckle noise. One reason is that the L2 loss in Equation (4) does not distinguish between the differences originating from signal and noise.

The differences between the constructed image and the input image can be defined as:

$$d=\left(\widehat{x}-x\right).$$

(6)

As shown in Figure 2c, in the difference image, the signal shows stronger self-similarity than the noise. A sorted non-local statics is proposed in this section to measure the signal autocorrelation. Minimizing the signal autocorrelation of the difference image can make the network parameters converge more quickly with the OCT signals.

We apply block matching to measure the non-local self-similarity in the difference image by calculating the correlation between the reference patch and each patch in the neighborhood ${\mathbb{C}}_i$. In this paper, the autocorrelation coefficient is applied as the metric:

$$c\left(d_i,d_j\right)={\sigma }_{d_i}{}_{d_j}/{\sigma }_{d_i}{\sigma }_{d_j}j\in {\mathbb{C}}_i.$$

(7)

where $d_i$ is the reference patch in difference image centered at $i$, ${\sigma }_{d_i}{}_{d_j}$ is the covariance of $d_i$ and $d_j$, and ${\sigma }_{d_i}$ and ${\sigma }_{d_j}$ denote the standard deviation of $d_i$ and $d_j$, respectively.

Note that the edge details in the difference image are not correlated with all their non-local neighborhood. To solve this problem, a sorted non-local statistic is proposed by sorting and selecting the mth largest autocorrelation coefficient. Based on the non-local statistic, the uncorrelation of difference image is obtained as follows:

$$u_l({\widehat{x}}_i)=1-{\displaystyle \sum _{m=1}^l{c}_m\left(d_i,d_j\right)}/lj\in {\mathbb{C}}_i,$$

(8)

where $c_m\left(d_i,d_j\right)$ is the mth largest one among $c\left(d_i,d_j\right)$. It can be seen in Figure 2d, the statistic provides a measure of how independent a patch is to its most similar non-local neighbors in the difference image. The signal uncorrelation of the differences is obtained as follows:

$$u_{\mathrm{s}}({\widehat{x}}_i)={\widehat{x}}_i\cdot u_l({\widehat{x}}_i).$$

(9)

The optimization loss function is therefore:

$$\begin{array}{cc}\hfill E(\widehat{x};x)& =L_{\mathrm{S}}+L_{\mathrm{C}}\hfill \\ & =\Vert \widehat{x}-x\Vert +\alpha \Vert 1-u_{\mathrm{s}}(\widehat{\mathit{x}})\Vert \hfill \end{array}$$

(10)

where the $L_{\mathrm{C}}$ maximizes the signal uncorrelation or minimizes the signal autocorrelation in the difference image. The weight $\alpha$ is used to balance two terms in the loss function.

Therefore, the parameters of the neural network can be trained by:

$$\widehat{\theta }=\underset{\theta }{\arg \min }\left(\Vert \widehat{x}-x\Vert +\alpha \Vert 1-u_{\mathrm{s}}(\widehat{\mathit{x}})\Vert \right)$$

(11)

Figure 2e,f shows the denoising results of proposed non-local deep image prior (NLM-DIP) after 450 and 900 iterations, respectively. Compared with the DIP results in Figure 1c,d, it can be seen that the NLM-DIP results are obviously clearer and maintain more details.

2.3. Network Structure and the NLM-DIP Algorithm

It can be seen from Figure 3 that the network structure employed in this work is based on the same U-Net type [25] 5-layer “hourglass” architecture with skip-connections as DIP. The purpose of using an hourglass network is to repeatedly obtain information contained in images of different scales. The architecture consists of an encoding path (upper side) and a decoding path (lower side). Each scale of the encoding path consists of two 5 × 5 convolution layers, the first convolution layer with stride 2 for down-sampling. Each convolution layer is followed by a batch normalization (BN) and a leaky rectified linear unit (LReLU). Each scale of the decoding path consists of a 5 × 5 convolution layer, a 1 × 1 convolution layer and an up-sampling layer respectively, each convolution layer followed by the BN and the LReLU. As shown in Figure 3, the numbers of convolution kernels in the 5-layer network are 16, 32, 64, 128, and 128, respectively. The parameters of the neural network can be trained by the proposed reconstruction loss in Equation (11), as shown in Figure 3.

The proposed NLM-DIP algorithm is summarized as follows:

Algorithm NLM-DIP

Input: Maximum iteration number MaxIt, network initialization ${\theta }^0$, noisy image $x$ and the network input $z$(random noise). For $i=1$ to MaxIt do   ${\widehat{x}}^i=f({\theta }^i|z)$   Calculate the signal uncorrelation of the differences using (9).   Calculate the reconstruction loss using (10).   Train the network using (11). End For Return $\widehat{x}=f({\theta }^{\mathrm{MaxIt}}|z)$.

In each iteration, the input $z$ are perturbed with an additive normal noise which is used as a noise-based regularization. The proposed sorted non-local static selects $l=20$ most correlated patches ($3\times 3$) in the $9\times 9$ neighborhood ${\mathbb{C}}_i$. The $\alpha$ in Equation (10) is used to balance the L1 loss and the proposed signal autocorrelation loss in the reconstruction function. In this paper, $\alpha$ is set to 0.85 according to the experimental results. We set a learning rate of 0.001 and 900 iterations for the network training. Currently the code of the proposed NLM-DIP is not optimized and the network training was done in GPU (NVIDIA GTX1060). The neural network was implemented using PyTorch. Because the NLM-DIP needs to calculate and sort the non-local correlation coefficients in each iteration, it takes 0.40 s per iteration for a 900 × 450 OCT image, while the DIP takes 0.18 s per iteration.

3. Experimental Results

In this section, the proposed NLM-DIP is compared with TNode [12], SBSDI [15], PNLM [10], and DIP [24] on two types of retinal SDOCT image datasets [15]. The synthetic image dataset contains 18 spectral domain OCT (SDOCT) noise images generated from high resolution images that were then subsampled; each image was obtained by an 840 nm wavelength SDOCT imaging system from Bioptigen with an axial resolution of ~4.5 μm per pixel in the tissue. The real image dataset contains 39 SDOCT noise images acquired at a natively low sample rate.

3.1. OCT Despeckling Results

In order to facilitate the subjective comparison of the denoising effect, two images were randomly selected from the synthetic and the real image datasets, respectively. The experimental results are shown in Figure 4 and Figure 5. Figure 4a and Figure 5a are the original noisy OCT images. Figure 4b and Figure 5b are the TNode results which contain some noise. Figure 4c and Figure 5c show the SBSDI results. It can be seen that SBSDI allowed the recovery of most of the image information while containing some speckle noise. On the contrary, the PNLM results shown in Figure 4d and Figure 5d removed most speckle noise while losing some details. By contrast, the results of DIP and NLM-DIP are more natural-looking. Compared with DIP, NLM-DIP generally output smoother surfaces in homogeneous regions and maintained layer texture and boundary clarity without introducing artifacts.

In addition, some interesting areas in the field in the OCT image have been selected in Figure 4 and Figure 5 which are marked with green boxes, and the partial area is enlarged. The results are shown in Figure 6 and Figure 7. The TNode results show that the images containing noise and some local details were not clear. As can be seen from Figure 6c and Figure 7c, the SBSDI algorithm blurred the local layer structure of the image, and the denoised result still contained a slight noise. The PNLM algorithm performed well in low-intensity regions, but the layer structure of the denoised result was unclear. The results suggest the DIP algorithm can be used for OCT denoising, but the layer structure of the image was ambiguous. The NLM-DIP algorithm efficiently removed the speckle noise. In addition, it maintained the details and the layer structure of the image, which was beneficial to further layer segmentation or clinical diagnosis and analysis.

3.2. Image Quality Metrics

Considering that the experimental images in this chapter have been processed and compared with the clinical images of medical experiments, there is a lack of original images, and it is inconvenient to perform signal-to-noise ratio measurements. In order to make objective quantitative comparisons of various denoising methods, this section uses the average contrast to noise ratio (CNR) and the average equivalent number of looks (ENL) as evaluation measures. CNR measures the contrast level between the selected regions from the object in an image and the background region in the image. It can be calculated as:

$${\mathrm{CNR}}_m=10\log \frac{|u_m-u_b|}{\sqrt{0.5({\sigma }_m^2+{\sigma }_b^2)}},$$

(12)

where $u_m$ and ${\sigma }_m$ are the mean and standard deviation of the $m$th interest foreground region (e.g., the areas in solid rectangles inFigure 4a, $u_b$ and ${\sigma }_b$ are the mean and standard deviation of the pixels in the background region (e.g., the dashed rectangle area in Figure 4a). ENL measures smoothness in areas that should have a homogeneous appearance but are corrupted by speckle:

$${\mathrm{ENL}}_m=u_m^2/{\sigma }_m^2.$$

(13)

To quantitatively compare different despeckling algorithms, we averaged the CNR and ENL images in the two datasets. It can be seen from

Table 1. Metrics comparison of different despeckling algorithms on the synthetic dataset.

	TNode	SBSDI	PNLM	DIP	NLM-DIP
CNR(dB)	8.67	10.72	10.25	10.32	11.17
ENL	544	1267	1420	1012	1468

and

Table 2. Metrics comparison of different despeckling algorithms on the real dataset.

	TNode	SBSDI	PNLM	DIP	NLM-DIP
CNR(dB)	9.92	11.28	11.12	11.14	11.65
ENL	596	740	793	595	1306

that the NLM-DIP algorithm provided the best despeckling performance in terms of CNR and ENL, compared with the other algorithm.

4. Conclusions

Current advancements in neural networks show their great applicability for unsupervised signal preprocessing. In this paper, a non-local deep image prior network was proposed for OCT image restoration. By adding a sorted non-local statics as regularization loss in the DIP learning, more low-level image statistics were captured by CNN networks in the process of OCT image restoration. The experimental results show that the proposed NLM-DIP allowed the recovery of most of the OCT signal while getting rid of the speckle noise.