Class-Specific Interferometric Phase Estimation Using Patch-Based Importance Sampling

Interferometric phase (InPhase) estimation, that is, the denoising of modulo- $2\pi $ phase images from sinusoidal $2\pi $ -periodic and noisy observations, is a challenging inverse problem with wide applications in many coherent imaging techniques. This paper introduces a novel approach to InPhase restoration based on an external data set and importance sampling. In the proposed method, a class-specific data set of clean patches is clustered using a mixture of circular symmetric Gaussian (csMoG) distributions. For each noisy patch, a “home-cluster”, i.e., the closest cluster in the external data set, is identified. An InPhase estimator, termed as Shift-invariant Importance Sampling (SIS) estimator, is developed using the principles of importance sampling. The SIS estimator uses samples from the home-cluster to perform the denoising operation. Both the clustering mechanism and the estimation technique are developed for complex-valued signals by taking into account patch shift invariance, which is an important property for an efficient InPhase denoiser. The effectiveness of the proposed algorithm is shown using experiments conducted on a semi-real InPhase data set constructed using human face images and medical imaging applications involving real magnetic resonance imaging (MRI) data. It is observed that, in most of the experiments, the SIS estimator shows better results compared to the state-of-the-art algorithms, yielding a minimum improvement of 1 dB in peak signal-to-noise ratio (PSNR) for low to high noise levels.


I. INTRODUCTION
Coherent imaging techniques, especially interferometry, plays an important role in many present-day technologies. Interferometry is a family of techniques that measure the phase differences of two or more waves (electromagnetic or acoustic) to infer physical parameters, such as small displacements, refractive index, topography of an irregular surface, etc. It has many applications in remote sensing [1], [2], medical diagnostic [3]- [5], surveillance [6], [7], weather forecasting [8], [9], and photography [10]. Some of the very relevant technologies in which phase imaging is a crucial part include interferometric synthetic aperture radar and sonar (InSAR/InSAS) [1], [2], [11]- [14], magnetic resonance imaging (MRI) [15], [16], optical The associate editor coordinating the review of this manuscript and approving it for publication was Yu Zhang.
Often, in a coherent imaging system, a physical quantity of interest is coded in a phase image. However, inference of phase is a challenging inverse problem due to the degradation in the sensing mechanism. Let z ∈ C be a complex-valued observation in a noiseless coherent imaging system which is related to the phase φ ∈ R via where j = √ −1. A major issue with this observation model is that the measured signal z depends only on the principal (wrapped) value of the absolute phase φ. Such sensing mechanism provides only a wrapped version of φ, usually defined in the interval [−π, π). We denote the wrapping operation as where W is the wrapping operator that performs modulo 2π wrapping, given by The wrapped phase is termed ''interferometric phase'', which is abbreviated as InPhase (φ 2π ) in this paper. Inference of the original phase (φ) from the observed InPhase (φ 2π ) is known as phase unwrapping (PU).
In addition to the above-mentioned non-linear distortion, the measured phase in a practical acquisition system is further corrupted by the noise inherent to the acquisition mechanism and the electronic equipment. Let a noisy observation be denoted as where n ∈ C is a complex-valued perturbation, often termed noise. Model (4) captures the essential features of most interferometric phase estimation problems, including magnetic resonance imaging (MRI) [13], [14], [19]- [21]. The observed noisy InPhase can be written as φ noisy 2π = W φ noisy = angle (z) .
Here, the observed InPhase φ noisy 2π is not only wrapped but also noisy. Estimation of the noiseless wrapped phase (φ 2π ) from noisy wrapped observation (φ noisy 2π ) is termed as phase denoising (PD) or InPhase estimation.
Though there have been a few attempts to address PD and PU jointly [14], [22]- [24], the common strategy is to follow a two-step approach in which PD is followed by PU. This paper focuses on PD, the first step of the two-step approach. We emphasize that PD is very different from natural image denoising due to the wrapping discontinuities present in InPhase images.

A. RELATED WORK
One of the simplest and most straightforward strategies in PD is to assume that the phase is constant in a small local neighbourhood. Local polynomial approximation (LPA) [25] and PEARLS [26] are two representative algorithms from this class. LPA approximates the phase in a rectangular window using a zero-order polynomial. The denoising performance of such algorithms depends crucially on the size of the window. LPA's performance is considerably affected by its fixed-size windows. PEARLS is another polynomial approximation algorithm that provides an adaptive-window framework to address LPA's limitations. It uses a zero-order LPA to estimate the window size selection and a first-order LPA to perform the filtering operations. However, PEARLS fails to deal with sharp discontinuities since the first-order polynomials are not adequate to model such structures.
Many recent image denoising methodologies benefit from the study of self-similarity and sparsity, often using patch-based techniques [27]- [29]. Natural images exhibit non-local self-similarity, meaning that they contain many similar patches at different locations. Among the selfsimilarity-based algorithms, BM3D [30] is a famous representative example, in which image denoising is accomplished through collaborative filtering of groups of self-similar patches. Although BM3D is not directly applicable for PD, there are a few variants of it proposed exclusively for InPhase applications [31]- [33]; among them, CD-BM3D and ImRe-BM3D [32] are state-of-the-art. The CD-BM3D algorithm performs the collaborative filtering through a 3D high order SVD (HOSVD); on the other hand, ImRe-BM3D incorporates a 4D HOSVD-based filter and treats the real and imaginary parts of the complex-valued signal as a pair. One major disadvantage of BM3D-based algorithms is that they do not incorporate a data-adaptive transform, which limits their performance for images with fine details, singularities, or sharp and curved edges. Although some data-adaptive versions of BM3D have been proposed [34], they are not developed for complex-valued signals.
SpInPhase [21] and NL-MoGInPhase [35] represent another type of recent PD algorithms that exploit self-similarity in complex-valued images. SpInPhase is a dictionary-based denoising method that utilizes sparse coding in the complex domain. In NL-MoGInPhase, PD is accomplished by modelling the complex-valued phase patches using mixtures of Gaussian (MoG) densities. The data-adaptive capabilities of these algorithms obtained by learning a dictionary or a MoG from the noisy data make its performance better compared to CD-BM3D and ImRe-BM3D [32].
Although the data-driven representations like SpInPhase and NL-MoGInPhase yield very powerful PD algorithms, the fact that, in many applications, the phase is locally smooth makes fixed representations in the Fourier domain still powerful tools in PD. Such locally smooth phase images can be approximated by a few windowed Fourier transform (WFT) coefficients. Windowed Fourier filtering (WFF) [36] is a relevant work in this category. Although WFF yields stateof-the-art PD performance in many InPhase applications, it is limited by its fixed-resolution due to the lack of an adaptive window. SURE-fuse WFF [37] is a recent PD algorithm, based on a multi-resolution WFF analysis to overcome the fixed-resolution limitations of the original WFF. In that algorithm, the WFF estimates having different resolutions are fused in a linear, pixel-wise, and data-adaptive manner by minimizing an unbiased estimate of the mean square error derived using Stein's lemma. However, SURE-fuse WFF suffers from high computational complexity since it requires the computation and fusing of WFF estimates having different window sizes.
One common drawback of most of the PD algorithms discussed above is that they do not have an algorithmic structure to make use of prior knowledge from an external database. In many practical applications, the image to be denoised is known to belong to a certain class. In such a situation, a large data set of clean images, termed ''external data set'' [38]- [40], from the same class, may be available, which can be used to tailor a statistical prior to enhance the denoising performance. Further development of this idea can be observed in [41]- [43], in which the principles of importance sampling (IS) are exploited to approximate the intractable minimum mean square error (MMSE) integral, which is a usual roadblock in image denoising. An important contribution of those works is that they bring flexibility in image denoising irrespective of the noise statistics.
Recently proposed class-specific image denoising algorithms [41]- [43] were mainly developed for real-valued images, and they do not adapt well to PD, in which the underlying signals are complex-valued. Although SpInPhase [21] and NL-MoGInPhase [35] can be modified to exploit a class-specific external data set [44], they were not originally developed to learn the dictionary (or MoG) from external noise-free data. This limitation motivates the need for a PD algorithm that can make use of clean external data sets, often available in many practical applications. Figure 1 illustrates a motivational application of class-specific PD, aiming to reduce MRI scanning time, which is an active research topic in medical imaging [45], [46]. In Fig. 1, a noisy MRI interferogram, possibly obtained from a quick scanning [47], of a person is denoised using clean MRI data from other people. A similar approach is applicable in other InPhase scenarios; for example, a noisy InSAR image of a mountain could be denoised using the InSAR data of other mountains. To the best of our knowledge, such class-specific strategies have not been explored in PD, and this paper is the first to advance in this direction.

B. PAPER CONTRIBUTIONS
In this paper, we investigate how to exploit an external data set to develop a class-specific PD algorithm. Inspired by [41]- [43], we propose a new estimation technique for InPhase images, which involves two major parts: • A shift-invariant clustering mechanism developed using csMoG densities to cluster external class-specific InPhase patches • A patch-based shift-invariant InPhase estimator developed using IS technique to denoise the noisy InPhase images using the clustered external patches. The remainder of this paper is organized as follows. Section II provides a brief introduction to importance sampling necessary to the following algorithm derivation. In Section III, the class-specific InPhase estimator is developed. The experiments and results are provided in Section IV and the paper concludes in Section V.

II. INTRODUCTION TO IMPORTANCE SAMPLING
Let us consider a random variable X distributed as X ∼ p X and a function f that is zero outside a set A. In many applications, it can happen that the set A is ''small'' or it is in the tail of the distribution p X as shown in Fig. 2. In such cases, the probability P(A) = A p X (x) dx is small. Now, to compute an estimate of any statistical characteristic of f (x), say the expectation E p X [f (X)], a plain Monte Carlo sampling from the distribution p X is not an effective method. This is because a plain sampling from p X could fail to provide even one sample from the set A. This is a typical problem that appears in many fields, namely, in high energy physics, Bayesian inference, rare event simulation for finance, insurance etc. [48,Chapter 9].
An intuitive approach is to sample from a different distribution, say q, that yields more samples from the ''important'' region, that is the set A in the above example. Such a sampling strategy is termed importance sampling (IS). Let us consider the expectation formula where D ⊆ R d is the domain of p X . 1 Let q be another positive probability density function defined on R d as shown in Fig. 2. IS is based on rewriting (6) by making q the sampling distribution instead of p X : where E q denotes the expectation operator under the density q. Here p X is a nominal density and q is an importance density. A Monte-Carlo approximation based on (7) yields the IS estimate of µ = E p X [f (X)] given by where the samples By the strong law of large numbers, P lim n→∞ µ q = µ = 1 [48,Chapter 9]. In the following section, we propose a new PD algorithm by adopting the IS strategy.

III. CLASS-SPECIFIC InPhase ESTIMATION USING PATCH-BASED IMPORTANCE SAMPLING
Our proposal is based on the assumption that an external clean database of the same class of the InPhase image to be denoised is available. Following the standard procedure in a patch-based approach, patches are extracted from both the noisy image and the external database. In this work, we extract square shaped patches of size √ m × √ m from the images in an overlapped manner with unite stride. The readers may refer to [21], [44] for further details of extraction and aggregation of patches. Let {z i } i=N p i=1 be the patches extracted from a noisy image z image ∈ C N 1 ×N 2 , where N 1 × N 2 is the pixel-size of the noisy image and N p is the total number of patches extracted from it. A noisy patch z i , using vectorized form of (4), can be written as where z i , n i ∈ C m and φ i ∈ R m . The noise n i is assumed to be i.i.d. zero-mean circular Gaussian [49] with variance σ 2 , i.e., n i ∼ CN (0, σ 2 I m ). All the operations in (9) are to be understood as component-wise. Let the set of external complex-valued patches be denoted as . We propose an InPhase estimation algorithm, which comprises the following steps: i) the external patches are clustered using a suitable clustering algorithm; ii) for each noisy patch, an external ''home-cluster'' is identified as the importance region; iii) InPhase estimation is accomplished by an estimator, developed using IS strategy, which makes use of the samples drawn from the home-cluster. These three steps are discussed in the following subsections.

A. CLUSTERING THE EXTERNAL DATABASE
The mixture of Gaussian (MoG) has been identified as a good model for multivariate complex-valued InPhase patches [35]. In this work, clustering is accomplished by fitting a mixture of Gaussian (MoG) to the external patches. Following [35], we assume circular symmetry (cs) to the Gaussian components of the MoG, which is a reasonable assumption that yields state-of-the-art PD performance [35]. The cs Gaussian expression for a signal y ∈ C m is given by where := E yy H , E yy T = 0, and E [y] = 0. The cs assumption has the following advantages: i) it simplifies the expression for Gaussian density involving complex-valued variables; ii) the cs Gaussian has a ''shift-invariance'' property, which is desirable in the context of InPhase patch modelling. To explain shift-invariance, let us consider a patch x i = e jψ i and another patch x i = e j(ψ i +1θ) , where 1 ∈ R m is a vector with all its elements equal to 1. Here x i is obtained by applying a common phase shift θ to the components of x i , as shown in Fig. 3. An efficient representation should treat x i and x i as similar patches irrespective of the common phase shift θ. This, in turn, facilitates the exploitation of the self-similarity between x i and x i in the InPhase estimation step. We term this property as shift-invariance and it is straightforward to verify that the expression (10) is shift-invariant. 2 In this work, we use circular symmetric MoG (csMoG), i.e., MoG whose Gaussian components are circularly symmetric, to model the InPhase patches. The csMoG density is given by where K is the number of components of the mixture model, k and α k are the covariance matrix and the weight of the k th Gaussian component, respectively, with 0 ≤ α k ≤ 1 and K k=1 α k = 1. The parameters of the csMoG are learned using the classical expectation maximization (EM) algorithm run on the external patches. The components of the learned csMOG are associated with K different clusters.
The cluster labels of the external patches are identified using the ''posterior weights'' (γ ik ) learned using the EM algorithm. We remark that γ ik indicates the ''responsibility'' of a particular MoG component k in the generation of the patch x i and is defined as After the EM algorithm converges, the cluster label of an external patch x i , denoted as L (x i ), is assigned to the cluster that has maximum value of γ ik , i.e., The pseudo code for the clustering is provided in Algorithm 1, whose main part is the EM algorithm. The algorithm can easily be derived by following a standard EM derivation adapted to csMoG density (see [50,Chapter 9] for further details).

B. IDENTIFICATION OF HOME-CLUSTER FOR THE NOISY PATCHES
Section III-A discussed the clustering of the external database. Our strategy is to make use of these learned clusters in the denoising process. For each noisy patch z i , a homecluster is identified and n samples from the home-cluster are randomly chosen to denoise z i . Each cluster is related to one Gaussian component of the MoG. In order to identify the home-cluster for a patch z i , the probability density values of z i under all the Gaussian components are calculated and the one that has maximum density value is considered as the
3 M-step: for k = 1, 2, . . . , K , Log Likelihood Evaluation: Convergence Check: where k is the covariance matrix estimated from the clean external patches using Algorithm 1. In (14), the covariance matrix corresponding to the noisy data is obtained by adding the noise variance σ 2 to k .

C. InPhase ESTIMATION WITH SHIFT-INVARIANT PATCHES
Let the distribution of the clean and noisy patches be ∼ p z , respectively. Given the prior p ψ , our objective is to estimate φ i 2π := W φ i , where W is the wrapping operator defined in (3). We adopt the following minimum risk Bayesian criterion: Therefore As discussed in Section III-A, for InPhase images, there can be many patches that are just phase-shifted versions of each other. These patches should be identified as self-similar, apart from a phase shift, to improve the efficiency of the estimator. This aspect was considered in Section III-A while developing the shift-invariant clustering. We now derive a shift-invariant InPhase estimator in this section. Let θ ∈ R be a common phase shift of a patch φ ∈ R m and φ 0 ∈ R m be the remaining part, i.e., We assume that φ 0 and θ are independent random variables, i.e., We use the same risk function as in (16); for a given patch φ, dropping the patch index i, where φ is an estimate of φ. Now (19) is rearranged by using (17) and (18) as follows: Since the model (9) assumes n ∼ CN (0, σ 2 I m ), the probability p z|ψ 0 ,θ (z|φ 0 , θ) in (20) can be written as Also, we denote the angle of the noisy patch as: Using (21) and (22), the risk function (20) is further expanded. In the following derivations, we use three constants c 1 , c 2 , and c 3 to account for the constants independent of φ.
In (23), we have assumed that p θ is approximately constant in any interval of length 2π; the constants accumulated from every 2π period is accommodated in c 2 . The required estimate is defined as φ k is computed in Appendix by solving the equation where I 1 is the modified Bessel function of the first order, and Computing φ k 2π for all values of k and stacking them leads to the vector solution as follows: We remark that the expression (28) is independent of the shift θ. Assuming that the class-specific external samples ψ 1 0 , ψ 2 0 , . . . , ψ n 0 ∼ p ψ 0 are available, the estimate of φ 2π can be computed as In (29), ζ i and b i depend on the external patch ψ i 0 and the noisy patch z. In order to understand the shift-invariance property of the estimator, let us substitute a shifted patch ψ i = ψ i 0 + θ1 instead of ψ i 0 in (29). In this case, it can easily be verified that the corresponding ζ i is ζ i − θ. Hence the corresponding exponent term of the estimator becomes  (29) with ψ i and the estimator automatically handles the shift contained in ψ i . By applying this shift-invariance property, (29) can be modified as VOLUME 8, 2020 which allows direct sampling from the external database without manipulating the shift, i.e., sampling from the density p ψ .

D. SHIFT-INVARIANT InPhase ESTIMATION USING IMPORTANCE SAMPLING
The estimator (30) assumes that the external data samples x i = e jψ i i=n i=1 are available to denoise a patch z. Let us consider a large external database x i = e jψ i i=N ext i=1 having N ext patches, where N ext n. As discussed in Section II, a random sampling of n patches from a large database may fail to yield even one patch with I 1 (b i ) ≈ 0 in (30). In order to tackle this issue, we use the IS strategy.
To develop an IS strategy, we make use of the learned clusters of the external patches. From (11), p ψ e jψ i can be written as Now, we choose an importance density (q ψ ), based on the home-cluster of the noisy patches computed using (14). Let k * = L (z) be the home-cluster of patch z.
Using (31) and (32), is evaluated as The RHS of (34) is evaluated using (12), which yields Inserting (35) in (33) leads to We emphasize that the expectation operation in (36) is w.r.t. the importance density q ψ . The IS estimator of φ 2π can be written as where α k * and γ ik * are easily computed from Algorithm 1 and expression (14) as follows: The estimator (37) is termed as Shift-invariant Importance Sampling (SIS) estimator. The pseudo-code of the proposed SIS estimation algorithm by summarizing the discussions in Sections III-A to III-C is shown in Algorithm 2.

SIS estimation using
where e jψ l 's are the drawn from the home-cluster k * = L (z i ) .

IV. EXPERIMENTS AND RESULTS
In this section, we illustrate the effectiveness of the proposed SIS estimator. Two types of data sets are used: i) semi-real InPhase data constructed from human face images and ii) real MRI interferograms.

A. EXPERIMENTS USING SEMI-REAL DATA SET
Let us consider a complex-valued image e jφ image ∈ C N 1 ×N 2 .
Here, we consider human face images (from FEI face data set [51]) as φ image . We would like to emphasize that the data considered here are not the phase images collected by means of any coherent imaging technique; they are the real images of the human faces captured using normal cameras. All the images considered here are of size 90 × 65. The pixel values are scaled to the range 0 to 9 to ensure a maximum number of wrapping equal to 1. Face images from 4 different persons, shown in Figs. 4a to 4d, are considered as the absolute phase φ image . The InPhase images, i.e., W φ image , are constructed from these four face images using the wrapping operator defined in (3). These images are shown in Figs. 4e to 4h. Since these InPhase images are synthetically created from real images, we term them as ''semi-real'' data. For the external database, 81 face images of the same size and the same range of pixel values as that of the test images, are considered. The external data set and the corresponding InPhase images are shown in Fig. 5a and Fig. 5b respectively. The following are the settings of the parameters: patch size 10×10, number of clusters K = 50, number of patch samples used for the SIS estimator in (37) n = 1000. Observations are generated as per the model (9), for low to high levels of noise corresponding to σ ∈ {0.5, 0.6, 0.7, 0.8, 0.9, 1}. Observations of all these noise levels are denoised separately  for each of the four test images in Figs. 4a to 4d. We choose 4 different algorithms to compare the performance of the SIS estimator: i) CD-BM3D [32], ii) WFF [36], iii) SpInPhase [21], and iv) SpInPhaseExt [21]. Although the first three algorithms do not make use of an external data set, they are selected for the comparison since they are state-of-the-art in InPhase denoising. However, we chose a fourth algorithm by modifying SpInPhase [21], which we term SpInPhaseExt. The SpInPhaseExt learns a complexvalued sparse dictionary from the external data set and uses it for denoising the noisy images. To the best of our knowledge, there are no other recent algorithms developed for classspecific InPhase estimation.
To evaluate the quality of the estimates, given that our main objective is InPhase denoising, we adopt peak signal-to-noise ratio (PSNR), defined as: [dB], (40) where φ image is the true unwrapped phase image, φ image 2π is the estimated InPhase, N is the number of pixels, and W is the wrapping operator defined in (3). In Figs. 6a to 6e, the images estimated from the noisy face-1 (Fig. 4i, σ = 0.9) using various algorithms are shown. From these estimates, it can be observed that the ability of SIS in preserving the minute details is much better, when compared to other estimates. For instance, the fine details associated with the eye region of the other estimates are either lost or partially preserved, whereas SIS retains these details. This claim is well supported by the PSNR value of SIS, which shows more than 2 dB improvements over the others. Also, the estimation error obtained by computing the difference between the true and the estimated InPhase is shown Figs. 6f to 6j. The estimation error of SIS is less structured, i.e., the features of the human face are less visible in SIS error image compared to the errors of the other two estimators. This, in turn, indicates that the SIS estimate is of better quality.
Next, we compare the performance of all algorithms using all the test images (Figs. 4a to 4d) by considering six different noise levels. The PSNR values are plotted in Fig. 7. For the low-noise case (σ = 0.5), SIS and SpInPhase show almost equal performance. This might be due to the fact that, at lownoise, the SpInPhase is able to learn a good quality dictionary. But as the noise level increases, the performance of SpIn-Phase degrades in comparison with SIS. When the noise is very high (σ = 1), SIS shows at least 2 dB improvement in PSNR compared to all other algorithms, for all the four test cases. One remarkable observation from Fig. 7 is that the SpInPhaseExt, which exploits the external data sets by learning dictionaries from them, does not perform as good as its original version SpInPhase. This indicates that the stateof-the-art dictionary learning based technique for PD is not well designed to exploit a class-specific external data set, highlighting the significance of SIS, since such data sets are often available in many practical applications. From this comparative study, it can be concluded that the SIS estimator outperforms all other previous methods, especially when the noise level is high.
In Fig. 8, we test the performance of the SIS estimator by varying the number of clusters (MoG components). The noisy InPhase data shown in Figs. 4i to 4l (σ = 0.9) are denoised using SIS estimators with cluster numbers K ∈ {10, 20, 30, . . . , 100}. The PSNR for each K averaged over all the four test cases are shown in Fig. 8. It can be observed that the PSNR increases from K = 10 to K = 50, and thereafter, it remains almost the same with a very slight increase. Similar observations were obtained with the InPhase data for other noise levels. We note that the number of clusters for the proposed SIS estimator is heuristically chosen to be K = 50, and further research on this aspect is left as future work.
In Table 1, we have also considered an additional algorithm termed ''SIS full sampling'' estimator as a reference algorithm for the time comparisons. SIS full sampling estimator is the SIS estimator where the importance sampling step (clustering) is discarded and a full sampling is considered. This means that for the estimation of each noisy patch, all the external patches (nearly half a million patches) are used, resulting in a huge run time of 2306 seconds. From Table 1, it can be verified that the SIS estimator yields very close results (≈ 0.5 dB difference) compared to the full sampling in a very short run time of 45.02 seconds. This clearly illustrates the role of importance sampling in SIS estimation. We would like to remark that although the WFF and SpInPhaseExt algorithms are much faster than the SIS algorithm, their PSNR values are considerably low. VOLUME 8, 2020

B. EXPERIMENTS USING REAL MRI INTERFEROGRAMS
Next, we present a challenging experiment in which a noisy medical image of a particular person is restored using clean images from completely different persons. Restoring medical images from a large external database of previously scanned images is a challenging task, which may have useful SpInPhase [21], SpExt: SpInPhaseExt [21], CD-BM3D [32], and WFF [36].
applications. For instance, in an MRI scan, the noise level depends on the scanning time. A patient may not be able to have very long scanning sessions due to medical reasons. In such a case, if the noisy images from a quick scan could be denoised with the help of an external database, it will be a great leap. Here, we open the door to such research possibilities.
The experiments presented here are conducted using MRI phase images from four different persons. These images are the interferograms of the head region obtained through the scanning along the front, side, and top orientations (see Fig. 9). The scanning was carried out on a 1.5 T GE Signa clinical scanner operating within Western General Hospital (WGH), University of Edinburgh [52]. A particular scanning orientation, say front view, is taken as a specific class. The front view interferograms of persons 1, 2, and 3 are used to construct an external data set. Patches from this external database are used to denoise interferograms of the fourth person's front view scan. The experiment is repeated for the other two scanning orientations, i.e., side view and top view. Figure 10 illustrates the restoration of a very noisy interferogram (σ = 0.9), where the estimates have PSNR values around 26 dB. The same parameter settings as described in Section IV-A are used.
The images estimated using different algorithms from a high-noise observation (σ = 0.9) are shown in Fig. 12.
We emphasize that the MRI phase images here considered are very challenging to denoise due to the presence of abrupt discontinuities. The WFF, being a Fourier-transform-based PD algorithm, usually performs better for images with smooth variations, since such images have better sparse representations in the Fourier domain. Also, CD-BM3D does not perform well for images with sharp discontinuities. It can be observed from Fig. 12 that WFF and CD-BM3D estimates are oversmoothed, and hence many minute details are lost. The visual quality of SpInPhaseExt estimate is better compared to WFF and CD-BM3D. However, SIS estimates are the best among all the estimates, approximately having at least 1 dB better PSNR.
The same experiments are repeated for the noise levels σ ∈ {0.5, 0.6, 0.7, 0.8, 0.9, 1} and the PSNR values are plotted in Fig. 11. As expected, the performance of WFF and CD-BM3D are not very competitive. Although SpIn-Phase shows good performance for very low noise levels, the PSNR degrades quickly as the dictionary representation of non-smooth high-noisy images is a difficult task. SpIn-PhaseExt is the only algorithm that shows competitive results with SIS, but at medium to high noise levels, SIS clearly outperforms SpInPhaseExt.
It is to be noted that, in these experiments, the external database is constructed using the images from just 3 persons. A larger external database with more images will definitely improve the restoration ability of the SIS estimator.
Limitations and future works: The present work does not discuss an optimum strategy for clustering the external database. The number of clusters (K ) associated with the proposed SIS estimator was decided based on a heuristic tuning, as shown in Fig. 8. However, the nature of the InPhase data has a vital role in the clustering strategy. As future work, we suggest further research on optimum data-driven clustering strategies. Also, the proposed algorithm uses fixed-sized square patches irrespective of the nature of data. A generalization of the current work by incorporating size-and shape-adaptive patches is expected to bring improvements in the algorithm performance, to which we will devote future work.

V. CONCLUSION
A novel patch-based method for restoring class-specific InPhase images has been proposed. The new algorithm uses patches from a clean external data set. The main part of the algorithm is an estimator, termed SIS estimator, which is developed based on importance sampling. Also, a clustering mechanism is developed to cluster the external data set by fitting a mixture of complex Gaussians. For each noisy patch, the SIS estimator identifies a home-cluster, draws samples from it, and performs the InPhase estimation. Both the clustering and the estimation technique are developed for complex-valued patches by considering the shift-invariance property.
The effectiveness of the SIS estimator is illustrated through a challenging medical imaging application in which the MRI interferogram of a particular person is denoised using the scan images from different persons. Also, a comparative study among the state-of-the-art PD algorithms is done using semi-real human face InPhase images, which concludes that SIS is a robust PD algorithm under a wide range of noise levels, yielding at least 1dB better PSNR than the previous state of the arts.