Fine feature retention of magnetic resonance images using intensity inhomogeneity weighted non-local filtering

Abstract

Electronic noise in the presence of intensity inhomogeneity severely inhibits the robustness of many traditional filtering algorithms. We propose an adaptive filter that utilises the intensity inhomogeneity bias field often observed in magnetic resonance imaging as a map for directing the magnitude of smoothing. This adaptive filter allows for controlled filtering of images where high signal-to-noise (SNR) regions are smoothed less, while low SNR regions are smoothed more intensely. We validated our method with a phantom and MRIs of the lumbar spine. Results showed that the proposed algorithm improved filtering of noise in the presence of intensity inhomogeneity.

Keywords::

• medical imaging and visualisation
• computer aided diagnosis, surgery, therapy and treatment
• imaging and visualisation in biomedical engineering
• imaging and visualisation in biomechanics
• multiscale imaging and visualisation

Introduction

A common form of corruption encountered in magnetic resonance imaging (MRI) of the lumbar spine is spatial inhomogeneity of the signal intensity, or intensity inhomogeneity. Intensity inhomogeneity presents difficult obstacles for many modern noise filters, primarily due to an assumption of a globally homogenous signal-to-noise ratio (SNR). In the presence of local intensity perturbations, static filtering can lead to over- or under-filtering that spatially mimics the local shifts in SNR. Therefore, when confronted with intensity inhomogeneity, it may often be more desirable to utilise adaptive filtering that is sensitive to local changes in signal strength.

Intensity inhomogeneity correction (IIC) has been the principle focus of researchers for decades (Henkelman 1985; Haselgrove and Prammer 1986; McVeigh et al. 1986); the predominant source of corruption manifests when the receiver coil suffers from a steep decline in signal sensitivity in the direction of increasing tissue depth. A plethora of correction methods have been proposed in the literature, and Vovk et al. (2007) provides a review and classification of these methods. Two IIC philosophies that have garnered recent attention are information minimisation (Vovk et al. 2004, 2006) and histogram matching (Styner et al. 2000; Shattuck et al. 2001). A hybridisation of these techniques was proposed by Salvado et al. (2006), which merges philosophies from both in what they term local entropy minimisation with bicubic splines (LEMS). We utilised LEMS to realise the proposed filter and provide a more detailed description of the method in Section 3.2.

In addition to intensity inhomogeneity, MRI is also confounded by corruption of electronic noise. MR noise arises during image reformulation of the raw signal by inverse discrete Fourier transform of the k-space. It has been shown that the noise component of the raw signal is characterised by Rician distribution as it arises from two independent Gaussian random variables (originating from the real and imaginary channels) (Bernstein et al. 1989; Macovski 1996; McRobbie et al. 2006). Similar to IIC, numerous methods have been developed over the years for MRI noise correction (NC). A review by Milanfar (2013) provides an elucidation of some modern image filters. Our analysis focuses on the non-local means (NLM) filter, which arises from a photometric kernel (Buades et al. 2005; Awate and Whitaker 2006).

The developments of IIC and NC methods have largely progressed as independent entities tailored to one artefact or the other. Conversely, our analysis focuses on developing a scheme that successfully addresses smoothing when both IIC and NC are required to enhance the retention of fine features. Various techniques have been proposed in the literature to combat this scenario, such as combining a filtering parameter with an anisotropic filter (Samsonov and Johnson 2004); another study developed a variational approach that utilised half-quadratic optimisation (Geman and Reynolds 1992); others devised a method by adaptively fitting the extracted bias field through simultaneous body and surface coil images (Fan et al. 2003). Kervrann and Boulanger (2006) developed a spatially adaptive filter, similar to an NLM kernel that optimises the weighting function by iteratively growing the search window while adaptively updating the filtering parameter.

A common measurement for describing intensity and noise corruption of an MR image is:

(1)

where at the i-th pixel, y (i) is the observed signal, x (i) is the signal of interest, b (i) is the bias field and n (i) is noise. Our objective is to develop a consistent and reliable technique for extracting the true signal x (i) from the noisy and biased signal y (i).

Solving for (1) elicits the question of which correction method should be performed first to achieve the most optimal results. Previous studies have concluded that IIC should precede NC, since IIC can disrupt the homogeneity of the embedded image noise (Montillo et al. 2003; Palumbo et al. 2011). However, these studies employed static denoising methods and, thus, did not investigate the relevance of adaptive filtering.

The aims of this study were twofold: first was to modify a NLM filter to adaptively update its kernel weight with values determined by estimates of the underlying bias noise. Using this technique, the degree of filtering is adaptively modulated by local SNR values, where low SNR regions experience more filtering and high SNR regions incur less filtering. The second aim was to determine the order in which LEMS and adaptive NLM should be performed. We evaluated the performance of both aims with a synthetic phantom, the popular image of Lena and MR images of the lumbar spine.

Materials

Test images

A 512 × 512 synthetic phantom was generated from a fractal image of repeating circles based on a quaternary form of an Apollonian gasket (Figure 1(a)). The image is spatially symmetric and has highly controlled features and sizes, ideal for analysing spatial intensity inhomogeneity. The circles are an analogous representation of fatty infiltrations and sparsely located adipose tissue dispersed within muscle. Each circle was randomly seeded with an 8-bit grey value of 50 (muscle) or 150 (adipose), which roughly approximates the values of T₂ imaging.

Figure 1 An Appolonian phantom (a) and Lena (g) corrupted with a realistic MR intensity inhomogeneity bias field (b),(h) and 9% Rician noise (c),(i). A zoomed section of the original images (d),(j), with the added inhomogeneity (e),(k) and 9% Rician noise (f),(l).

In addition, we also tested with the popular image of Lena (Figure 1(g)). We wanted to explore the robustness of bias adaptive non-local means (baNLM) when confronted with more complex textures and contrasts than the two-toned phantom. Lena provided a suitable platform for testing such features. Throughout the duration of this manuscript, references to this image are termed ‘Lena’.

Synthetic corruption

A realistic bias field was added to the phantom. This was performed by element-wise multiplication of the bias image and synthetic phantom or Lena (Figure 1(b), (e), (h), and (k)). Next, Rician noise was added to the bias-corrupted images (Figure 1(c), (f), (i), and (l)) at five Rician noise levels: 1%, 3%, 5%, 7%, and 9%. Results were averaged over 10 iterations at each of the five Rician noise levels. The same phantom and image of Lena were used throughout testing.

MR images

The lumbar paraspinal muscles of 26 human subjects were evaluated for this study. The participants were volunteers recruited from the staff, student and patient populations at Vanderbilt University Medical Center (VUMC) and the University of Tennessee (UT), Knoxville. Institutional Review Board (IRB) approval was obtained as well as informed consent for all subjects participating in this study (UT IRB #7393; VUMC IRB #070883).

Axial T₂-weighted images of the target lumbar muscles were acquired on a 1.5 T MR^† unit (^†Achieva; Phillips Medical Systems, Amsterdam, Netherlands) with a SENSE-spine^† phased array coil. The images were captured using a turbo spin echo scanning sequence with: 512 × 512 resolution, 4 s repetition time, 120 ms echo time, 4 mm slice thickness, 90° flip angle and a 225 × 225 mm² field of view. A stack of 50 images was acquired for each subject spanning the lumbar spine from approximately T12 to S1.

Three bilateral paraspinal lumbar muscles were used for testing: the psoas, erector spinae and multifidi (Figure 2). The psoas muscles lie lateral to the vertebral body and anterior to the transverse processes; the erector spinae (superficial and deep) lie in the groove of the vertebral column and transverse processes; and the multifidi sit medially to the erector spinae and in the groove of the spinous and transverse processes (Figure 2). T₂ weighting displays pixel values of muscle as dark grey to black, fat/adipose as white and bone as light grey (Figure 2).

Figure 2 Sample MR slice at L4 with bulk muscle boundaries highlighted. The highlighted muscles include the bilaterally occurring: psoas (bilateral to the vertebral body), erector spinae (posterior to psoas muscles) and multifidi (bilateral to the spinous process and posterior to the transverse processes). (Note: The image has been brightened and contrasted to enhance relevant tissue boundaries).

Semi-automatic segmentation

Bulk muscle measurements were performed by manually segmenting the outer fascial boundaries of each of the six muscles of an image at L1, L3 and L5. Gold standard manual segmentation was performed by NVB (four years segmentation experience), which was validated by ASM (novice segmentation experience) by independently segmenting the same images and comparing with a matched pairs t-test. Intermuscular adipose tissue was further segmented by fuzzy c-means (FCM) 2-cluster classification (Dunn 1973; Bezdek 1981).

Summary of operations

In summary, we tested two different approaches for intensity inhomogeneity and noise correction using both NLM and baNLM (Figure 3).

Figure 3 Two varying methodological paths used for realising the bias adaptive NLM filter. The scheme (top path) passes the intensity corrected image to baNLM, while (bottom path) passes the original image to baNLM and uses LEMS both before filtering, for bias extraction, and after filtering, for intensity correction.

Methods

Notation convention

In all, we tested four different configurations of the corresponding IIC and NC methods. We denote them using the following convention: NLM as , bias-adaptive NLM as , and LEMS as . The combined orders of operation are tabulated in Table 1.

Table 1 Notation for orders of operation.

Orders of operation	Notation
Bias-adaptive NLM → LEMS
NLM → LEMS
LEMS → bias-adaptive NLM
LEMS → NLM

Local entropy minimisation by bicubic splines

LEMS is an optimisation algorithm that estimates the global background bias field through estimates of the entropy in local sub-regions. The entropy, S, is estimated by , where PDF_X (·) is the probability density function of X and N is the window containing the local sub-region. The values of PDF_X (N) are approximated by the histogram of the pixels in the current region divided by the number of pixels in the current region. LEMS optimises the entropy in a piecewise fashion by calculating the entropy in local square neighbourhoods (or knots). The knots are then ranked based on their SNR. At the first knot, the bias field is optimised by minimising the entropy in the local region. Next, the second highest SNR knot is optimised for both the current and previous regions; the third knot is then optimised for the first, second and third regions and so on, until all knots have been optimised. Concurrently, the local bias fields are stitched together with a bicubic spline that smoothes the gradient of the global bias field, which is b (i) from (1).

NLM filtering

Given a noisy image , the NLM algorithm estimates the true image by a weighted average of all pixels in the image, , where represents the family of weights that depend on the similarity between pixels i and j. This measures the grey scale similarity between the vectors y (N_i) and y (N_j), where N_k denotes a square neighbourhood (or the similarity window) centred at the k-th pixel. The similarity is a measurement of the photometric distance, , where a>0 is the standard deviation of the Gaussian kernel. The kernel is defined as:

(2)

where Z(i) is the normalising constant

(3)

The parameter h dictates the strength of filtering by controlling the exponential decay of the weighting function. Larger values of h flatten the decay of the weighting function, thus causing greater filtering. Smaller values of h sharpen the decay, forcing the weights to zero and causing less filtering. The parameterisation of h is the principle point of interest in our analysis of the connection between image denoising and intensity inhomogeneity.

The coupled algorithm: baNLM filtering

By varying the constant h in (2) and (3), the level of local blurring can be accurately controlled. In applications of the NLM filter, h is generally parameterised as a function of the global standard deviation of the noise. This is, however, impractical for lumbar MRIs where: (1) the noise variance is not explicitly known and (2) there is poor spatial homogeneity of the SNR. The inhomogeneous SNR properties that the bias field introduces to the image can be illustrated by corrupting a uniformly noisy signal with a realistic bias field, and then correcting for the intensity inhomogeneity with LEMS (Figure 4).

Figure 4 A realistic MRI inhomogeneous intensity field (a) corrupted with 9% Rician noise (b). Intensity inhomogeneity corrected for with LEMS (c), whereby local SNR values are calculated (d). As expected, the intensity of the SNR image (d) mimics that of the initial intensity field (a).

If we then calculate local SNR values of the corrected image, we find a globally inhomogeneous SNR pattern that mimics that of the original bias field (Figure 4(a),(d)). Since NLM assumes globally homogenous SNR, we conclude that utilising a static value for h may result in over- or under-filtering at certain regions of the image, depending on local SNR values. In contrast, it would be more desirable to parameterise h that more accurately reflect local SNR values. Since the extracted bias field provides us with this estimation, we propose an adaptive parameterisation of the NLM kernel by modulating h with local values of b. This provides the framework for NLM filtering with adaptive properties that are guided by estimations of the SNR.

In the context of the NLM kernel, the extracted bias field provides relative estimations of local SNR values that are used to dictate filter weighting (Figure 5). Figure 5 illustrates the correlation between the probability of the kernel and local SNR values. The green similarity window indicates an approximate mean of the SNR where average filtering approximately occurs. The orange similarity window spans the brightest region of the image and correlates to the highest SNR (Figure 4(d)). The high SNR (or low noise standard deviation) is indicative of higher probability weighting from (2) and (3), and thus produces a sharper filtering curve (Figure 5 (left), orange curve). Conversely, the blue similarity window reflects a region with low SNR where heavier weighting of (2) and (3) generates a flatter curve to reduce the noise with greater variance (low SNR) (Figure 5 (left), blue curve).

Figure 5 Three similarity windows aligned at various coordinates along the bias image (right). The mean pixel intensity contained in each window reflects the probability density (left) for parameterising the adaptive kernel function.

Since h is inversely proportional to the intensity of the extracted bias image, the compliment of the bias is needed to achieve the desired spatially sensitive blurring effects. The compliment of the bias image is found by:

(4)

where is the bias compliment and MAX (·) is the maximum grey value allowed by the bit depth of the image.

Maintaining the notation developed in Section 3.3 for the similarity window and centre pixel, we define the adaptive filtering parameter as the mean of within the corresponding similarity window about the current centre pixel i:

(5)

where β (i) is the adaptive weighting parameter, l is the current pixel in the similarity window, and m is the size of the similarity window. The term β (i) serves to modulate the filtering intensity by dictating the local SNR about the similarity window. In this view, β (i) is an instantaneous scalar multiplier to h, which when combined with (2) is

(6)

and with (3) as

(7)

As shown in Equations 456(7), the baNLM kernel weights utilise values of the bias estimate to modulate the degree of filtering. In our analysis of the two test images, we optimised h for NLM and then initialised baNLM with the same value of h to ensure an unbiased evaluation of each filters' performance.

Performance measures

Evaluation of the baNLM algorithm was performed by four quantitative performance measures. We first used the root mean square error (RMSE), which measures the difference between pixel values of a reference image and the denoised (corrected) image:

(8)

where x (i) is the i-th pixel of the reference or corrected image and N_R is the number of pixels in the reference or corrected image. Lower values of RMSE indicate better performance.

The second measure was Pratt's Figure of Merit (FOM) (Pratt 1978). FOM is a quantitative measure of edge preservation and enhancement, and is calculated from the Canny edge map of the denoised and referenced images:

(9)

where N is MAX{N_corrected, N_ideal}, N_corrected is the number of detected edge pixels in the corrected image and N_ideal is the number of detected edge pixels in the reference image, d_i is the distance between the ith edge pixel of the corrected image's Canny map and the next closest edge pixel in the reference map, and γ is a scaling constant defined to be 1/9 as indicated by Pratt (1978). The index of FOM ∈  [0, 1], where higher values indicate better performance.

The third performance measure used was the structural similarity index (SSIM), which is a quantitative measure of the structural similarities between the reference and denoised images (Wang et al. 2004):

(10)

Lastly, semi-automatic segmentation of the MR dataset was evaluated against the manually segmented gold standard by way of the Dice coefficient. This is a comparative metric used for testing the similarity between two images (Dice 1945), calculated by:

(11)

where x_FCM is the semi-automatically segmented image, x_GS is the manually segmented gold standard and N is the number of pixels in one image. Calculations of Dice ∈ [0, 1], where higher values indicate more accurate segmentation.

Results

Parameterisation of h

Brute force optimisation of h at each noise level is tabulated in Table 2. They reflect values of h that returned the best performances for and . They were likewise used to parameterise the two adaptive filters to ensure equivalent testing.

Table 2 Average performance of the four filters on the synthetic phantom and Lena with various levels of noise.

		Synthetic phantom				Lena
Noise	Metric									h
1%	SSIM	0.7877	0.7828	0.7888	0.7892^*^*	0.8211	0.7555	0.8248^*,^*^*	0.7523
	FOM	0.4988	0.6881	0.5554	0.7313^*,^*^*	0.9295	0.8625	0.9289	0.8605	0.0002
	RMSE	0.2767	0.2821	0.2769	0.2625^*,^*^*	0.1024	0.1581	0.1023^*,^*^*	0.1581
3%	SSIM	0.8094	0.8323	0.8202	0.8424^*,^*^*	0.6877	0.6167	0.6937^*,^*^*	0.6168
	FOM	0.4534	0.8062	0.5334	0.8199^*,^*^*	0.8232^*,^*^*	0.6731	0.8174	0.6827	0.003
	RMSE	0.2427	0.2280	0.2378	0.2168^*,^*^*	0.1053	0.1565	0.1047^*,^*^*	0.1583
5%	SSIM	0.8618	0.8720	0.8761^*,^*^*	0.8641	0.6009	0.4995	0.6075^*,^*^*	0.5023
	FOM	0.4500	0.6013	0.5347	0.6354^*,^*^*	0.6596	0.5622	0.6630^*,^*^*	0.5691	0.02
	RMSE	0.1827	0.1531	0.1616	0.1469^*,^*^*	0.1063	0.1339	0.1061^*,^*^*	0.1368
7%	SSIM	0.8651	0.8681	0.8754^*,^*^*	0.8596	0.5414	0.4521	0.5456^*,^*^*	0.4547
	FOM	0.3739	0.5060	0.4440	0.5620^*,^*^*	0.6916^*,^*^*	0.5407	0.6848	0.5441	0.03
	RMSE	0.1468	0.1243	0.1270	0.1235	0.1100	0.1318	0.1096^*,^*^*	0.1353
9%	SSIM	0.8403	0.8410	0.8434^*,^*^*	0.8389	0.4811	0.3952	0.4840^*,^*^*	0.3985
	FOM	0.3764	0.4142	0.4132	0.4486^*,^*^*	0.6573	0.4599	0.6821^*,^*^*	0.4785	0.05
	RMSE	0.1388	0.1276	0.1374	0.1285	0.1177	0.1277	0.1168^*,^*^*	0.1307

Notes: Statistically significant (p < 0.01) performance values for a matched pairs t-test are indicated with the appropriate symbol (footnote). Values in bold indicate the filter that obtained the best performance and that a test of significance was performed (see footnotes).

^*Statistically different (p < 0.01) than the next closest filter for the same noise level and performance measure.

^**Statistically different (p < 0.01) than the result of the complementing filter (e.g. the complement of is , and the complement of is ).

Phantom

Filtering results of the synthetic phantom showed that baNLM generally outperforms NLM, particularly at lower levels of noise (Table 2). We found that achieved the better results than the three other filters for all FOM and RMSE measures (with the exception of 9% RMSE) and SSIM at 1% and 3% noise, while the results for SSIM at 5%, 7% and 9% noise were better for . Nine per cent RMSE was the only test where a non-bias adaptive filter () did not achieve the best score.

Lena

The performance measures collected on the image of Lena showed that the and procedures consistently outperformed their counterparts, and (Table 2). That is, in contrast to the phantom, the IIC-then-NC sequence produced superior results for every test and condition. Additionally, we observed that outperformed in all measures except for FOM at 1%, 3% and 5% (1% not significantly).

Magnetic resonance imaging

Spatially varying filter strength was achieved in the real MR data (Figure 6). As hypothesised, less smoothing was performed in the posterior muscle region (erector spinae and multifidus), while more was in the anterior region (psoas). Results are reported for the MRI data as determined by the accuracy of the FCM tissue classification (e.g. Figure 7) and quantised by the Dice index. Figure 7 illustrates two examples of the classification accuracy for muscle tissue. Relatively high accuracy by both filters (Figure 7(a)), though baNLM had approximately half as many mutually exclusive misclassifications as NLM (Figure 7(d)). Even greater disparity was observed between the filters in Image 2, where baNLM had approximately seven times less misclassified pixels as NLM (Figure 7(d)). Additionally, the classifications of baNLM of both test Images 1 and 2 were more accurate than NLM (Table 3). baNLM outperformed NLM in the erector spinae and multifidus muscles based on mean paired differences of the Dice coefficient (Table 4). The left erector spinae at L5 was the exception, as no significant differences were observed for both filter comparisons. No differences were observed at the psoas muscle, except at L1. Using the total segmented area, the baNLM algorithm returned higher Dice index values than NLM at all three vertebral levels and both configurations.

Figure 6 The original MR image (a) and compliment of the extracted bias field (e). Shown are the unbiased image corrected with just LEMS (b) and the same unbiased image filtered with (f). Zoomed regions of the unbiased image are shown in (c) and (d), while (g) and (h) are the same zoomed regions for . Greater filtering was observed at (g) relative to (c), while less filtering occurred at (h) relative to (d).

Figure 7 Automatic segmentation accuracy of muscle tissue for two different MR images (a),(b). Each image was processed independently with baNLM () and NLM (), and then classified with FCMs clustering. The classification results of the images from both filters are overlayed on a single image and coded based on correspondence of the classification results. The right-hand bars of (d) depicts regions where baNLM and NLM were mutually correct, while the left-hand bars of (d) indicates where they were mutually incorrect (b),(d). The right-hand bars of (c) indicates that baNLM was correct and NLM was incorrect, while the left-hand bars of (c) indicates that NLM was correct and baNLM was incorrect (a),(c). Incorrectly classified pixels were identified as both false-positives and false-negatives, while correctly identified pixels were associated with only true-positive muscle tissue. The classifications from both images show that baNLM (left-hand bars of (c)) had significantly less mutually exclusive misclassifications than NLM (right-hand bars of (c)) in both images (c).

Table 3 Classification accuracy percentages of test Images 1 and 2 (Figure 7).

Image	Filter	TP	TN	FP	FN	Accuracy
1	baNLM	88.53	6.43	2.92	2.12	94.96
	NLM	89.55	4.39	4.96	1.09	93.94
2	baNLM	76.66	7.19	1.23	14.92	83.85
	NLM	70.12	7.53	0.89	21.47	77.64

Note: Reported are the percentages of pixels classified as: true positive (TP), true negative (TN), false positive (FP), false negative (FN) and percent accuracy.

Table 4 Dice values for the mean differences over 26 subjects at each vertebral level and muscle group for baNLM against NLM.

		Psoas		Erector		Multifidus
Mean difference	Level	Right	Left	Right	Left	Right	Left	Full Seg.
	L1	0.0095	0.0112	0.0077	0.0091	0.0164	0.0148	0.0096
	L3	0.0011	0.0012	0.0044	0.0053	0.0126	0.0123	0.0051
	L5	0.0046	− 0.0012	0.0139	0.0500	0.0080	0.0082	0.0059
	L1	0.0240	0.0112	0.0141	0.0123	0.0242	0.0195	0.0145
	L3	0.0020	0.0004	0.0055	0.0060	0.0143	0.0141	0.0057
	L5	− 0.0011	− 0.0009	0.0240	0.0534	0.0085	0.0116	0.0063

Notes: Full Seg. refers to the collective segmentation of all six muscle groups. Values in bold indicate a significant difference (p < 0.05) for a paired t-test.

Lastly, residual images are a subtraction of the filtered image from unfiltered. They provide a qualitative metric for investigating a filter's ability to perform noise reduction without losing fine features important to the image. The residual of a perfect filter would yield a purely noisy image that contained no discernable features of the original image. An example MRI residual of baNLM (Figure 8(a)) clearly shows better fine features retention than the same image filtered with NLM (Figure 8(b)). Performance was markedly improved in the posterior region at the muscle–fat interface. The anterior region, however, appears to have similar qualitative properties. These results further confirm the findings from Table 4 that showed more accurate segmentation in the posterior muscle, while the segmentation for the anterior muscles was not significantly different.

Figure 8 Residual images of an example MR slice that was filtered by baNLM (a) and the same image filtered by NLM (b). Visual inspection reveals less discernable anatomical features in baNLM than NLM.

Discussion

The fundamental premise behind most kernel-based image filters is a globally homogeneous SNR. However, as is shown in this work (Figure 4), the presence of intensity inhomogeneity necessarily results in an incongruous SNR after performing intensity correction. This can severely compromise the performance of many filters, as well as other post-processing intensity-based techniques. The dilemma is particularly relevant to MRI of the spine, which often exhibits a drastic drop in intensity moving away from the spine due to the use of a one-sided surface coil. We have shown that it is superior to smooth images of this sort with a filter that is spatially sensitive to local changes in the relative intensity of the image.

Results for the order of operations on the phantom indicated that out performs at low levels of noise (1% and 3%). Yet, it remains unclear as to why at higher levels of noise performed better with the FOM and RMSE measure, while performed better for SSIM. Nonetheless, the dynamic version of NLM repeatedly outperformed its static counterpart. Lastly, we note that outperforming was in direct contrast to the findings reported in the studies mentioned earlier (Montillo et al. 2003; Palumbo et al. 2011).

Conversely, the performance measures gathered on Lena showed the opposite of those gathered on the phantom where a NC-then-IIC sequence was clearly more beneficial to restoration of Lena. The discrepancies observed between the two signals may be due to the content of the tested images; that is, the difference between a two-toned phantom and an image (e.g. photograph) that contains more complex shades, contrasts and textures. It appears from our data that NC-then-IIC may perform better when subjected to smoothed, gradient images, whereas IIC-then-NC is better at delineating between sharp edges with little gradient and contrast.

Given the observed differences between the two sequences, we determined that MR images of the lumbar spine, where fascial boundaries are abrupt, correspond more with the phantom image. Therefore, the NC-then-IIC procedures were utilised for testing our proposed adaptive kernel on the test MRIs.

Results of the real MR datasets showed that baNLM performed better than NLM for the erector spinae and multifidus muscles and equivalently for the psoas muscles. The muscle-specific results are not overly surprising considering the spatial location of the muscle groups. The likely scenario was that since the psoas muscles generally sit in the central/anterior region of the image, the bias values associated with the filtering of the psoas muscles were likely initialised to a similar value to h in NLM. This implies that both filters would have been performing the same (or similar) degree of filtering in the local regions of the psoas muscles, which is supported by the residuals of Figure 8. Conversely, the posterior muscles, which sit in a high SNR region, clearly did not require as much filtering as was performed by NLM. It could be argued that it would be possible to ‘tune’ the NLM filter to optimise the accuracy in the posterior region; however, it would likely result in significant over-filtering of the anterior region. This, however, further illustrates the need for spatially sensitive filtering when in the presence of image intensity inhomogeneity.

Our proposal of a bias adaptive kernel opens the door for a slew of other experiments and adaptations that could be performed. For example, we used the mean value of the similarity window around the centre pixel, but it could be worthwhile to investigate the effect of combining the mean value across the moving and centre windows; or perhaps how baNLM performs better with the mean value of the search window; or using the centre pixel value of the bias image. We could also adapt our technique with a search window growth optimisation method (Kervrann and Boulanger 2006). Numerous other possibilities exist that could be explored, one or more of which may prove to perform better than the approach that we have outlined in this manuscript.

We used a Rician noise model to perform our experiments since it is commonly associated with MRI denoising and is widely utilised by the medical image processing community for MRI analysis. However, it is also relevant to note that the noise model associated with MRI SENSES-Encoding is not actually known. This is primarily due to interpolation that occurs between adjacent images. It would certainly be viable to explore the performance of baNLM with other noise models.

Finally, though we have employed LEMS and NLM to develop the fundamental framework behind our proposed filtering algorithm, we have generalised it in a way that is compatible in any kernel-based filtering algorithm. The only criterion for utilising the ideas behind this bias adaptive filter is an analogous kernel parameter to h that controls the filtering strength. Similarly, any IIC algorithm could be employed since the bias image is calculated directly from the processed and unprocessed images.

Conclusion

We have developed a technique for adaptively guiding the NLM filter in a way that utilises intensity inhomogeneity to control the degree of filtering. The premise lies with the notion that intensity inhomogeneity inherently shifts local SNR values, which violates a primary assumption used by most modern image filters. By taking advantage of the correlation between latent bias field and local SNR values, we steer the kernel weights with an extracted image of the bias field to accurately adjust the degree of filtering in local regions of the image. This configuration allows us to more aptly satisfy the assumption that the strength of filtering should be proportional to local SNR values.

Conflict of interest disclosure statement

The contents are solely the responsibility of the authors and do not necessarily represent the official views of NIAMS.

Acknowledgements

The authors would like to thank the following from the University of Tennessee, Knoxville, for providing valuable assistance in this work: A. Ali, L. Bowers, C. Carr, E. E. Abdel Fatah, J. Mitchell, A. Sharma, and G. To.

Additional information

Funding

This publication was made possible by the grant from the National Institute of Arthritis and Musculoskeletal and Skin Diseases (NIAMS) at the National Institutes of Health [grant number 5R01AR055882].