Structural adaptive anisotropic recursive filter for blind medical image deconvolution

Received Nov 29, 2018 Revised May 1, 2019 Accepted May 15, 2019 Performance of radiographic diagnosis and therapeutic intervention heavily depends on the quality of acquired images. Over decades, a range of preprocessing for image enhancement has been explored. Among the most recent proposals is iterative blinded image deconvolution, which aims to identify the inheritant point spread function, degrading images during acquisition. Thus far, the technique has been known for its poor convergence and stability and was recently superseded by non-negativity and support constraints recursive image filtering. However, the latter requires a priori on intrinsic properties of imaging sensor, e.g., distribution, noise floor and field of view. Most importantly, since homogeneity assumption was implied by deconvolution, recovered degrading function was global, disregarding fidelity of underlying objects. This paper proposes a modified recursive filtering with similar nonnegativity constraints, but also taking into account local anisotropic structure of content. The experiment reported herein demonstrates its superior convergence property, while also preserving crucial image feature.


INTRODUCTION
Recent advances in medical imaging technology has so far enabled high performance computerized radiographic diagnosis and therapeutic intervention [1][2][3][4][5]. More specifically, it has been widely applied, for examples, in patient specific anatomical modeling, lesion extraction and more recently in unsupervised deep learning [6]. Thus far, degradation is one of major impeding factors in their success. Although in practice, it is led by a series of complex processes imaging signal underwent during acquisition, for simplicity, the term is typically characterized by linear deconvolution of blurring kernel and an additive noise [7], as expressed in (1).
where g and f are degraded and (presumably) original images, respectively. In the spatial domain of (x, y)  R 2 , h and n are convolutional kernel and noise, respectively. In analyzing degradation process, h is sometimes referred to as blur filter or, in our context, point spread function (PSF). Reconstructing the original image f, given the degraded g and a priori on (or sometimes, unknown) noise model, n, is however not trivial. Its key element involves estimating a PSF and its respective inverse (h -1 ). This process is called deconvolution and depicted in Figure 1. Depending on assumption of the degrading model, PSF can be estimated by calculating its governing (for example, blurring) parameters [8,9]. However, several other factors, such as out of focus, motion, and geometrical distortion, etc. may have equally contributed to degraded image quality. Identifying these simple model parameters was unable to completely restore such adverse effects. Determining types of degradation, their ranking and interactions, indeed are not trivial, especially without access to imaging modality calibration. Another approach, called blind image deconvolution (BID) [7,10,11], tackled this problem by estimating each element in PSF kernel/ matrix (or its inverse), subject to some criteria, but without prior on degradation sources. BID combines PSF estimation and deconvolution into a single reciprocal process. Estimation of PSF is usually done by iteratively updating vectorized kernel using a gradient descent variant [12], whereby in each cycle its values are varied with respect to pre-defined objective function and constraints. Provided that an optimum exists, upon convergence, resultant inversed PSF are able to closely recover the true non-degraded image. It can be noted that noise model was not incorporated into the inversed PSF and hence would have caused instability in case of low SNR. Other studies thus opted for operations in frequency domain, in which linear property of convolution can be exploited [13,14], i.e., F* (u, v) = G (u, v) H -1 (u, v), where capital letters refer to (estimated) true, degraded images and PSF in frequency domain, respectively. Accordingly, noise may be dismissed by selectively processing only in the lower frequency spectrum. Although iterative variants exist, PSF estimation in frequency domain are of close-form and more efficient, compared to that on spatial one. However, its main drawback was that knowledges about noise properties are prerequisite, without which severe instability could occur. To remedy this adverse effect, Wiener filter [15,16], Wavelet [17] and Curvelet [18] based methods were proposed.
Note may be drawn from the literature that both spatial and frequency domain operations have their pros and cons. While the former can greatly benefit from straightforward yet intuitive constraints imposition, the latter is more efficient, with available fast spatial-frequency conversion algorithms. BID on both domains, called a nonnegativity and support constraints recursive filtering (NAS-RIF) [7,19] was introduced and recently enhanced [20][21][22][23][24][25][26]. Its primary contribution was to overcome instability issue found in conventional BID. As its name suggested, NAS-RIF imposes irreducible, absolutely summable, i.e., ( || h || < ), and being invertible, i.e., h -1 , properties on a PSF, while maintaining the same assumption (i.e., real value and positive definite) on the true image as IBD. The support on the image was defined within a region of interest (ROI). This has made NAS-RIF particularly suitable for medical imaging, where anatomical object is generally acquired in the center of a matrix and surrounded by uniform background. Specifically, NAS-RIF divided the degraded image into two regions, i.e., inside and outside support (Dsup), whose cost functions were determined with different objectives and constraints. A generic NAS-RIF algorithm is summarized in Figure 2. Given a degraded image g (x, y), NAS-RIF recursively determines the optimal inversed PSF, denoted here as u (x, y). At each iteration the elements of u are optimally adjusted subject to a cost function, consisting of negativity penalizes for pixels within and outside support region. Some enhanced proposal suggested additional regularization constraints, e.g., DC gain [21] or more realistic boundary support [22]. Noise issues were elevated by using designated filters [23][24][25]. For better performance nonetheless, deconvolution was performed in the frequency domain.
The summary of NAS-RIF algorithm (see text for detailed explanation): a. Definition  )x,y( : Estimate of the true image at k th iteration  )x,y( : FIR filter parameters of dimension at iteration k  : Cost function, given parameter setting Despite its great stability, shortcoming of NAS-RIF was slow convergence rate. Moreover, the inversed PSF (u) was a compromise between two penalize terms, derived from pixel intensity. This paper, therefore proposes a structural adaptive anisotropic term being introduced in iterative optimization. It was computed, taking into account local orientation pattern of object structure. Its main contribution was not only emphasizing on updates in favor of feature preservation, but also promoting faster convergence as fidelity was enhanced. The remaining of this paper is organized as follow: Section 2 describes the proposed method in more detail. Section 3 reports experimental results and relevant analyses. Finally, section 4 states concluding remark of this study.

RESEARCH METHOD
This paper partly adopted conventional NAS-RIF following the process, depicted in Figure 2a. It iteratively adjusted FIR filter, u, and simultaneously its output, that was an intermediate estimation of the true image, f*. To ensure nonnegativity and support constraints, this image was projected onto a non-linear (NL) space that diminished pixel intensities outside the support region (Dsup) to that of the background (LB). The corresponding non-expansive map, f*NL were then subtracted with its precedent, resulting in error matrix, given in (2) [20].
where , , 0 if , 0 and , ∈ if , 0 and , ∈ if , ∈ Instead of simply minimizing this matrix, which would prematurely bring estimation to halt, the error image was then divided by supporting region. Ideally, object pixels bled outside Dsup due to degrading (e.g., blurring) PSF should be drawn back inside, leaving only the background. This could be achieved by an edgeenhancing FIR. However, exaggerating this adjustment could lead to negative pixels by deconvolution. To maintain the balance of these constraints, basic NAS-RIF defined the cost function to penalize negative pixels inside the Dsup and background discrepancies outside. According to the recent modification [21], irreducibility of FIR was also ensured by regularization. The cost function adopted in this study was thus given in (3). Using this function, it is trivial to prove that its gradient with respect to an FIR element. Given the cost function J and its respective gradient, J, iterative non-linear optimization was done using a conjugate gradient method [26].
where sgn and  were, respectively, a signed function and an empirical factor weighting FIR regularization. Upon convergence, when the difference between successive true image estimations fell within a pre-defined threshold, the resultant FIR filter was then applied to the degraded image, g, producing the final true image restoration.
It was, however, reported in the recent NAS-RIF literature that noise reduction and a priori on underlying pixel distributions are essential determinant in its stability and restoration result. In addition, penalizing cost function, while sufficient for typical photographic images, did not consider structural fidelity in an image, hence undermining anatomical features, crucial for the subsequent analyses. Inspired by intuitive constraints augmentation found in the recent works, this paper therefore incorporated an anisotropic measure into NAS-RIF optimization. Structural anisotropic measure was introduced in [27] and later improved in [28]. In those studies, it was used to orient and adjust the extent of an adaptive FIR filter so that it aligned with underlying pixel orientation pattern. Anisotropic measure within neighborhood surrounding a pixel, p, is given in (4).  In order to avoid further complicating the cost function (3), which would inevitably cause even more local minima, this study instead encouraged feature preservation by adjusting updating step size t according to relative anisotropic strength, as expressed in (5) where  was a typical stepping size taken in each iteration. The benefits of adjusting step size according to total relative anisotropic were two folds. Firstly, while steering away from non-negativity, original NAS-RIF cost function tended to stumble around an overly smoothing kernel. With anisotropic controlled step size, on the other hand, as the bled pixels were gathered inside, implying more pronounced object boundaries, the relative anisotropic measure also increased and so was the confidence in such adjustment. This effectively accelerated NAS-RIF convergence. Secondly, involving anisotropic measure into the optimization also helped lessen the dependency on having to meticulously initialize the supporting region [22]. It is also worth emphasizing here that, anisotropic measure was computed within a neighborhood of specified extent and not from an isolated pixel. It was thus robust against imaging noise [27].

RESULTS AND ANALYSIS
Without loss of generalization, the proposed enhanced NAS-RIF algorithm was examined by applying to both synthetic and medical images corrupted with known degradation. The images were encoded as 2D matrix of grayscale intensities, whose values were stored and processed in floating point format.

Anisotropic strength as image contrast regularization
As pointed out in [20] and subsequent works, trivial all-zero condition could be prevented by imposing a total sum constraint on FIR kernel. We found that it did not, however, rectify a uniform FIR kernel that would bring the image contrast tremendously down to an all-grey. To demonstrate that in addition to structural pattern [27] anisotropic measure is also responsive to such condition (and thus was a viable means of circumvent this problem) relationships between contrast appearances and respective total anisotropic strength are shown in Figure 4. During an early stage of optimization, image contrast could be regularized by anisotropic strength. More specifically, as the FIR proceeded away from trivial all-zeros, the measure helped increasing its confidence by further stepping in that direction.

Visual enhancement
An MR scan of a human brain on a uniform background whose matrix size was 350x350 pixels, was then employed in the next experiment. Comparison between enhancement made by a generic NAS-RIF and the proposed enhancement against an original MR image are illustrated in Figure 5. The results are snapshots at the 80 th iteration. The resultant inversed FIR, u, brought sharper edge and better separation between tissue and the skull. Dynamic range of pixel intensities was much improved, compared to the generic NAS-RIF. It is also worth noted here that, instead of thousands of iterations usually required by a generic NAS-RIF to converge [7,19], the proposed enhanced NAS-RIF gave an estimation with already higher fidelity and contrast at much early cycles. Moreover, no other priors were needed. Figure 6 compares anisotropic strengths during the first 80 th iterations between a generic and the proposed NAS-RIF implementation, and respective enlarged original image estimations. It is evident from the graph that in the proposed implementation, the strength accelerated at faster rate, which well corresponded to much enhanced appearance. It was thus a suitable metric for a NAS-RIF optimization constraint and well conformed to the preliminary hypothesis of this study.

Numerical assessments
To quantitatively elucidate the proposed NAS-RIF scheme, especially in terms of noise immunity, numerical assessment was performed on simulated adulteration To this end, a phantom image was degraded with Gaussian blur and polluted with Rician noise (to emulate what happens in MR acquisition). Peak SNR [29,30], was then computed for the original, adulterated, and generic and proposed NAS-RIF enhanced images. The peak SNR (PSNR) and corresponding restored images at 40 th iteration were listed in Table 1 and shown in Figure 7, respectively.  Although quantitatively and visually, there was only slight improvement in PSNR over the generic NAS-RIF, the proposed method was equally if not better immune to additive noise. To elucidate that the proposed enhanced NAS-RIF could equally well applied to other images, the anisotropic strength and convergence were compared against the baseline one. The experiments were performed on README, Brain Phantom and Real Abdominal CT images. It is evident from the Figure

CONCLUSION
Blind image deconvolution is an ill-posed problem that was designed to restore the true image, undergone degradation by an unknown PSF and possibly by random noise. Variational BID elevates this problem by iteratively estimate the PSF (or its inversed) subject to some pre-defined criteria. NAS-RIF is another well accepted variation BID that imposed non-negativity and supports constraints over a sequence of restored image, during the optimization. Nonetheless, it is prone to noise and had low convergence rate. Many attempts had been made in the literature to address these issues, by suggesting various FIR regularization schemes, selectively filtering the projected image, or accurately defining the object support, etc. This paper put emphasis on quality of the shape and object definition and thus proposing a structural adaptive metric, i.e., anisotropic strength. Its advantages are robustness against noise and intuitively representing characteristics of local orientation pattern. Unlike other recent works, this paper did not augment anisotropic strength into an already complicate NAS-RIF cost function, or else it would have caused minima traps and created another unnecessary expression to be weighted and balanced. On the contrary, it was used simply to adjust the step size in each kernel update. The benefits were two folds; it accelerated convergence as object boundaries became more pronounced and structural appearance. It concisely represented the structural appearance of underlying object and thus lessen the need of precise initial support.
The experimental results reported herein confirmed visually and numerically that the proposed NAS-RIF had much higher convergence rate, offered restoration of better quality, and was equally immune to synthetic noise. It was therefore believed that the proposed method could offer a new direction toward improving the performance of the widely adopted NAS-RIF, especially in the fields of medical imaging, computer aided diagnosis (CAD), and digital anatomy.

ACKNOWLEDGEMENT
This work was supported by Suranaree University of Technology, under One Research One Graduate (OROG) grant scheme.