Restoration of Bi-Contrast MRI Data for Intensity Uniformity with Bayesian Coring of Co-Occurrence Statistics †

The reconstruction of MRI data assumes a uniform radio-frequency field. However, in practice, the radio-frequency field is inhomogeneous and leads to anatomically inconsequential intensity non-uniformities across an image. An anatomic region can be imaged with multiple contrasts reconstructed independently and be suffering from different non-uniformities. These artifacts can complicate the further automated analysis of the images. A method is presented for the joint intensity uniformity restoration of two such images. The effect of the intensity distortion on the auto-co-occurrence statistics of each image as well as on the joint-co-occurrence statistics of the two images is modeled and used for their non-stationary restoration followed by their back-projection to the images. Several constraints that ensure a stable restoration are also imposed. Moreover, the method considers the inevitable differences between the signal regions of the two images. The method has been evaluated extensively with BrainWeb phantom brain data as well as with brain anatomic data from the Human Connectome Project (HCP) and with data of Parkinson’s disease patients. The performance of the proposed method has been compared with that of the N4ITK tool. The proposed method increases tissues contrast at least 4.62 times more than the N4ITK tool for the BrainWeb images. The dynamic range with the N4ITK method for the same images is increased by up to +29.77%, whereas, for the proposed method, it has a corresponding limited decrease of −1.15%, as expected. The validation has demonstrated the accuracy and stability of the proposed method and hence its ability to reduce the requirements for additional calibration scans.


Introduction
Structural MRI can provide high resolution three-dimensional data sets of organs and tissues within the body [1].The images can be used to examine tissue integrity as well as to diagnose a variety of disorders [2,3].The examination can be qualitative, but it can also be quantitative [3].The extensive quantification and analysis of the data requires further processing with automated methods such as registration and segmentation [4].This further processing is hampered by the presence of imaging artifacts of non-biological origin, s.a.intensity non-uniformities across an image.This artifact mainly stems from the inhomogeneity of the Radio-Frequency (RF) field due to the coil as well as due to its interaction with the subject.The inhomogeneity is more pronounced in high field MRI, which is ≥3.0 Tesla, despite the higher contrast to noise ratio that it can potentially provide.
It is possible to calibrate for the non-uniformity of the radio-frequency field with prospective methods that involve additional acquisitions.A prospective acquisition method uses the frequency responses to parameterized acquisition sequences [5,6].Another method involves the acquisition of a series of proton density weighted images in advance to estimate the inhomogeneity [7].The imaging of physical and geometric phantoms has also been used to approximate the combined non-uniformities of the transmission and the receiver coil(s).The physical correction methods typically involve additional acquisitions that are valid only for particular MRI sequences as well as geometries.Thus, prospective methods may not be sufficient and may reduce contrast.They also increase acquisition time that renders them more vulnerable to increased heat deposition on tissues and to motion artifacts.A general physical formulation for the complicated interaction between the radio-frequency field and the anatomy of the human body is not currently tractable.
Several retrospective restoration methods have also been proposed.They often restore individual images [8].Some are based on a primarily data-driven criterion and assume a distinction between lower spatial frequencies corresponding to the non-uniformities and higher ones corresponding to the anatomy.A direct implementation of this uses a low-pass homomorphic filter [9].In parallel imaging reconstruction, polynomial fitting of the coil data effectively performs low pass filtering that provides the respective sensitivity maps [10].Another primarily data-driven approach represents an image with its spatial intensity derivatives of significant magnitude that are assumed to correspond to tissues boundaries [11,12].The reintegration of the normalized intensity gradients gives piecewise constant regions [11].Piecewise constancy has also been implemented with the Total Variation (TV) [13].However, the TV prior is a Laplacian distribution that is uni-modal and hence favors intensity transitions of uniform magnitude throughout an image.The regions of the same tissue that are spatially disconnected are restored independently.The spatial piecewise smoothness has also been used together with intensity clustering, s.a.fuzzy C-means, over local image windows [14][15][16].These formulations have been extended with level sets that can provide an explicit representation for two [17] or more regions corresponding to tissues [18].
Retrospective restoration has also been performed based on the histogram.A global histogram statistic, the entropy, has been maximized [19].However, the entropy is optimized for aligned mode distributions and that can be limiting [20].Thus, the histogram has been combined with the coring methodology [21] that was originally developed for noise removal [22,23].In the context of MRI inhomogeneity correction, coring has been combined with a spatial smoothness constraint for the field [24,25].However, the histogram methods are sensitive to the level of image noise [4,26].The histogram based methods can also lead to instabilities in the dynamic range and in the image domain [19,24,25].
The intensity correction has been addressed together with other image analysis steps for specific applications.It has been combined with local image denoising for diffusion weighted images [27] and with intensity standardization [4,28].Intensity statistics and the histogram have also been combined with constraints for tissue classification [14,16].The latter have been implemented with tissue intensity priors and with tissue spatial priors to perform simultaneously intensity restoration as well as explicit tissue segmentation [29].Patch-based priors have also been used for intensity correction [30].
A patient imaging protocol typically includes multiple sequences, each of which is reconstructed independently.The resulting images can suffer from different non-uniformities.A variety of retrospective methods has been developed for their joint non-uniformity correction.A data dependent method fuses an MRI image with a Positron Emission Tomography (PET) dataset and processes the resulting image with a low-pass filter [31].Another data dependent method uses a variational formulation that enforces smoothness of the non-uniformity fields and also preserves the differential structure of the images [32].The joint restoration of two images has used the minimization of the entropy of their joint intensity histogram [33].A method for the simultaneous restoration of multiple co-registered images minimizes the sum of the entropies of voxelwise stack vectors over the entire image domain [34].The bi-contrast and multi-contrast restoration methods assume that the valid signal domains of all the images involved are identical.
A retrospective method has been developed in this study that performs a joint intensity restoration of images of two contrasts representing a certain anatomic region of a subject.In general, the retrospective methods can benefit from regularity properties of the anatomy and from the physical properties of the non-uniformities that are valid for a range of MRI contrast mechanisms and geometries.The proposed method uses the Bayesian coring methodology that involves implicit non-parametric priors.The priors are based on the auto-co-occurrence statistics for each image [35,36] as well as on the cross-or joint-co-occurrence statistics of the two images.The effect of the intensity distortions on the co-occurrence statistics is modeled and the statistics are restored.Additional constraints have been imposed to achieve a stable and robust restoration.These include the smoothness of the non-uniformity, the co-occurrence representation, the standardization of the dynamic ranges, and the consideration of the inevitable difference between the signal regions of the two images.
The effectiveness of the method has been demonstrated extensively with BrainWeb brain simulated images [37].It has also been demonstrated with brain anatomic datasets of the Human Connectome Project (HCP) (Data were provided (in part) by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 National Institutes of Health (NIH) Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University.)and with brain anatomic data sets of Parkinson's disease patients.The performance of the method has also been compared with the performance of the N4ITK tool [38,39].The comparison showed the superior performance of the proposed method.Thus, the method has been shown to improve the accuracy, the efficiency, and the stability of the restoration.It also reduces the requirements for additional acquisitions for calibration.

General Bayesian Formulation
The Bayesian formulation involves the model of the distortion that is used both directly to give the likelihood and indirectly to give the prior.In this work, the distortion is multiplicative in image domain and so the statistical distortion filter and the restoration are non-stationary.The prior is represented non-parametrically.The Bayesian estimation is repeated iteratively, t.

Spatial and Statistical Image Representation
The general formulation is introduced for a single anatomic image v(x), with x = (x, y, z), which is corrupted by a multiplicative spatial intensity non-uniformity b(x).This is due to the combined effect of the transmit and the receive MRI radio-frequency inhomogeneities over an underlying latent anatomic image u.The image is also corrupted with additive Rayleigh noise, n [40].The noise is independent and identically distributed.That is, the images model is: where the multiplication • is the voxelwise, Hadamard product.The probability distributions of v, u, and b are defined over the entire image domain.The domain is discrete and the probability distributions are computed voxelwise.The probability distributions of u(x) and b(x), P u (u(x)) and P b (b(x); 1, σ 2 b ), respectively, are assumed to be independent.It is shown in the Appendix A that they give the probability distribution of P v (v) as: where * is the convolution.That is, the latent density of u is convolved with a non-stationary distortion filter over the dynamic range of u.
The image regions with b(x) > 1 become brighter and the image regions with b(x) < 1 become darker.Each of these two ranges corresponds to a mode in the distribution of the non-uniformity representing the distortion of the statistics.Thus, the Point Spread Function (PSF) of the distortion is a bimodal distribution.It is taken to be where G(•) is a Gaussian distribution, 2 is a regularization parameter, and k is a proportionality constant.The proportionality constant k is computed numerically so that the area of the PSF, P b is equal to one.The bimodal shape of P b is shown in Figure 1.

Posterior Expectation for Voxelwise Intensity Restoration
An overview of the Bayesian framework is in Figure 2. It gives the estimate for the posterior expectation of the assumed latent intensity, û = E(u|v) = P (u|v) (u|v)udu.The substitution of Bayes' law expansion for P (u|v) (u|v) leads to: û = E(u|v) = P (u|v) (u|v)udu = P (v|u) (v|u)P u (u)udu This expression involves the probability distribution for latent variable u, P u (u).The latent distribution is estimated from P v (v) with the deconvolution of the distortion P b (v − u|u) from the actual intensity distribution P v (v).The deconvolution provides the prior as a non-parametric distribution, Pu (u).
The PSF for the distortion of the intensity co-occurrences, P b , increases linearly with intensity value.It is non-stationary in its direct domain that are the dynamic ranges.Thus, the deconvolution has to be non-stationary in the direct domain.The deconvolution is performed with the regularized Van Cittert method [41,42].That is, the restoration is an iterative regularized algorithm formulated as where β is the regularization constant and P 0 v is the probability of the initial image.The completion of the iterations provides the estimate Pu (u).The Van Cittert algorithm is formulated for stationary deconvolution filters.In this work, it is extended to accommodate the non-stationary PSF of P b .

Back-Projection of Intensity Restoration to the Images
The restored statistics are estimated and enforced back to the image domain in a voxelwise manner.The likelihood P (v|u) (v|u) in Equation ( 4) is equivalent to the non-stationary PSF in Equation (2) of the intensity statistics.The prior P u (u) in Equation ( 4) can be approximated by Pu (u) = Cu as estimated from Equation ( 5) for the co-occurrence statistics of image u.These two terms are substituted in Equation ( 4) to give The density of P b (v − u; 0, (σ b u) 2 ) for values of u that are far from v tends to zero.Thus, the discretization of Equation ( 6) considers only a finite neighborhood ∆u ∈ N v to give: The size of the distortion filter P b (∆u; 0, (σ b u) 2 ) and of the neighborhood N v are non-stationary and increase linearly with intensity v.
The value of the initial restoration field, b −1 , that is, the inverse of the distortion field, can be derived from û as it is given from Equation (7).This is achieved with: This provides the voxelwise restoraton expression.The restoration is iterative and the updated estimates of the restoration are used at every iteration, t.At t, the estimate of v t provides an updated estimate for ût .In turn, these two, v t and ût , provide  The value of the restoration field for iteration t at location x depends only on the intensity of v t (x).Thus, the restoration can be precomputed and stored in a one-dimensional array of size equal to that of the dynamic range.It can then be indexed to expedite the back-projection of the voxelwise restoration to the image at iteration t.The sequence of steps of the method are shown in Figure 3.The restoration in the proposed methodology is joint for two images.

Methods
The Bayesian formulation is extended in various ways to represent issues related to the specific problem.An extension is for smooth spatial non-uniformities.Another extension is for the joint restoration of bi-contrast data.The identification of the signal domains of the images allows the restoration only over regions where the non-uniformity assumptions are valid.The numerical implementation of the restoration is iterative and the dynamic ranges are enforced to be stable.

Spatial and Statistical Image Representation
The method involves images of two different constrasts, v i (x), where i = 0, 1.The two images are in the same anatomic space and their domains have the same spatial sampling grid.They can be corrupted by different multiplicative spatial intensity non-uniformities, b i , that represent the MRI radio-frequency inhomogeneities over the underlying anatomic images u i , i = 0, 1.Each image is also corrupted with additive Rayleigh noise, n i [40].That is, the model becomes and is identical for both images.
The Taylor series expansion of the non-uniformities, b i (x), around a voxel x 0 gives The first order term provides a linear approximation of the non-uniformity within a spherical neighborhood N ρ of radius ρ = |x − x 0 |.The effects of the second and higher order terms are neglected within N ρ .The statistical representation of images v 0 and v 1 is based on intensities η 0 and η 1 within neighborhood N ρ .Their counts give the co-occurrence statistics as [35,43]: They give the auto-co-occurrences for i ≡ j and the joint-co-occurrences for i = j.The auto-co-occurrences are dominated by their diagonal entries and thus they are weighted down with the sigmoid 1/(1 + e −(k 1 |η 0 −η 1 |+k 2 ) ).The different tissues, or tissues interfaces, of the anatomic images u i are assumed to correspond to distinct modes of the co-occurrence statistics.An example of the co-occurrences of a pair of T 1 and T 2 BrainWeb phantom images [37] is in Figure 4.A median filtering is applied to v i to remove the high frequency noise n i .The objective of the analysis is to separate the remaining two products in v i , i = 0, 1 present in Equation ( 9) to obtain the factors b i and u i , respectively.

Statistical Representation of Intensity Non-Uniformities
The intensity distortions b i and the latent images u i are assumed to be generated by independent random variables.As shown in the Appendix A, the statistics of the products b i • u i correspond to the convolutions of the co-occurrence statistics of the anatomic images C u i u i and C u 0 u 1 with the PSFs of the corresponding distortions.The PSFs are non-stationary to account for the multiplications in the spatial domain.
The effect of the intensity distortion field b i in a neighborhood N ρ around x 0 can be approximated by the effect of the zero order term b(x 0 ) and by the effect of the first order term ∇b(x)| x 0 in Equation ( 10) [35].The zeroth order term scales the auto-co-occurrences radially around the origin of C u i u i and the first order term rotates them around the origin.The point spread functions affecting C u i u i can be more efficiently represented in polar coordinates (r i , φ i ).The standard deviation of the radial scaling is linearly proportional to the radial coordinate σ r i ∝ r i .The standard deviation of the rotation σ φ increases with ρ and is largest along the diagonal and zero along the axes.The application of the PSFs of the intensity distortions to the auto-co-occurrence statistics of the assumed underlying anatomic images u i , which represent the auto-co-occurrences of the distorted images.
The effects of the zero order terms of the non-uniformities b i of the two images on the joint-co-occurrences C u 0 u 1 are also considered.The PSF affecting C u 0 u 1 is represented in Cartesian coordinates.Its sizes are linearly proportional to their distances from the origin σ η i ∝ η i .The relation between η i and r i is The application of the PSF of the intensity distortion to the joint-co-occurrence statistics of the assumed latent anatomic images u 0 and u 1 that correspond to and represents the joint-co-occurrence statistics of the distorted images.

Non-Stationary Restoration of the Co-Occurrence Statistics
As elaborated in Section 3.2 and in the Appendix A, the PSFs of the non-uniformities are non-stationary in co-occurrence space with their standard deviations being linearly proportional to the intensity coordinates.They give the convolutions in Equations ( 12) and (13).Hence, the deconvolutions of these PSF are also non-stationary.They are implemented with the non-stationary Van Cittert method in Equation (5).The restoration of the co-occurrence statistics C v i v i , i = 0, 1, provide their corresponding restored estimates, Cu i u i , i = 0, 1, respectively.The restoration of the joint co-occurrences C v 0 v 1 gives the estimate Cu 0 u 1 .

Back-Projection of the Co-Occurrences Restoration to the Images
The general Equation (7) for the case of the auto-co-occurrence statistics gives the posterior expectation of the intensities in polar coordinates in neighborhoods ∆r i ∈ N r i and ∆φ i ∈ N φ i as where P b ((∆r i , φ i ); 0, (σ b,r i r i ) 2 , σ 2 b,φ i ) is as given in Equation ( 12).This expression for the PSF, P b , also gives Cu i u i from C v i v i .The posterior expectation results from separable filtering.Equation (7) for the case of the joint-co-occurrence statistics gives the posterior expectation of the intensities in Cartesian coordinates in neighborhoods ∆η 0 ∈ N η 0 and ∆η 1 ∈ N η 1 as where P b ((∆η 0 , ∆η 1 ); 0, (σ b,η 0 η 0 ) 2 , (σ b,η 1 η 1 ) 2 ) is given in Equation ( 13) that also gives Cu 0 u 1 (η 0 , η 1 ) from C v 0 v 1 (η 0 , η 1 ).The conditional restoration filtering is also separable.The general Equation ( 8) provides the restoration values from the auto-co-occurrence statistics in Equation ( 14) as R s i,t (r, φ) = ri r i , i = 0, 1.It also provides the restoration value from the joint-co-occurrence statistics in Equation ( 15) as The intensity restoration values are used to construct rectangular restoration matrices with sizes equal to the respective dynamic ranges.They are the auto-co-occurrence restoration matrices and the joint-co-occurrence restoration matrix.The repeated indexing of these matrices in an iteration expedites the back-projection of the restoration to the images.

Estimation of the Cumulative Intensity Restoration
At the first iteration, t = 0, the restoration field is initialized with W i,t (x) = 1, ∀x, and thus, equivalently, the image is initialized with the acquired data v i,t=0 .At subsequent iterations, the intensity co-occurrence statistics index the corresponding auto-restoration and joint-restoration matrices to provide an initial incremental estimate of the restoration The estimates b −1 i,inc,t at iteration t multiply the estimate of the corresponding cumulative restoration, W i,t−1 of the previous iteration t − 1 to give The estimates W (x) are smoothed with a spatial Gaussian filter G x; σ 2 s,i to give the restoration fields These are the estimates of the smoothed cumulative restoration fields.The estimates W i,t , i = 0, 1, at iteration t, are applied to v i,t to provide u i,t = v i,t+1 that are the updated estimates of the underlying latent anatomic images.
The end condition of the iterations involves the standard deviation of W i,t , σ(W i,t ) = The standard deviation decreases with iterations.The iterations stop when the value of σ i,t reaches a minimum for at least one of the two images i.The field W i,t that corresponds to the iteration with the minimum value of σ gives t min = argmin t=0,...t max σ i,t , which provides the restoration field W rest i = W i,t min .The restored images are given by u rest i = W rest i v i .A maximum number of iterations t max is also imposed.

Validity of Image Domains and Stability of the Images' Dynamic Ranges
The dynamic range of an anatomic MRI imaging sequence can be variable.The range is typically from a few hundred to a few thousand intensity levels.The dynamic range does not have a physical significance and may not even have a statistical significance.It is also affected by artifacts.In an MRI image, these can be regions of only noise, artifacts from blood flow, or artifacts resulting from tissues' susceptibilities.The latter two can cause outlying intensities of very low or very high values.
The non-uniformity is a physical characteristic of the signal regions of an image.Thus, the dynamic range of an image is pre-processed to identify its valid sub-range for which the non-uniformity model holds.The limitation of the dynamic ranges result in distributions that are statistically meaningful.The limitation of the intensity ranges also make the filtering for deconvolution more efficient.This standardizes the dynamic ranges to normal values that correspond to the meaningful tissues contrasts intended for by the imaging sequence.It also gives the Region of Interest (ROI) in an image, I ROI,i .The combined ROI of the two images is their union, I ROI = (I ROI,0 ∪ I ROI,1 ).The standardization of the dynamic ranges makes the magnitudes of the distortion PSFs meaningful.
The upper parts of the image dynamic ranges are downsampled linearly.The resulting complete dynamic ranges in MRI data may also be downsampled linearly.
The reference intensity η re f i to identify and standardize the valid range is the intensity value giving a high cumulative percentage, the 90% of the dynamic range of the ROI that still corresponds to tissues.The dynamic range of the noise delimited by the maximum intensity of the noise range or minimum signal intensity, η min i , is represented as a low fraction, 10%, of η re f i .The minimum signal intensity η min i = 0.1η 0.9 i for an image depends on the standard deviation of the Rayleigh noise σ N i .The dynamic range is preserved exactly up to a large upper value, 150% of η ), for i = j, and for i = j, respectively.The same ranges are used for the Van Cittert restorations of these statistics that enable with Equation ( 14) the computation of R s i and with Equation ( 15) the computation of R b i .The joint restoration R b i also considers the complete dynamic range of the complementary image, i = j.The restoration gain factors in R s i and R b i that correspond to entries out of the valid intensity ranges are set to unity.The next step is the back-projection to the images to obtain the initial estimates of the restoration fields b −1 i,inc that use a window N x in Equation ( 16).The back-projection with Equation ( 16) for a rough initial estimate is over the valid part of the ROI of image i and is set to unity outside.The resulting estimate is smoothed as described in Equation ( 17) to give W i using a window depending on σ s .The spatial smoothing considers a weight field that is equal to one in the valid region I ROI,i = 1 and is much less than one in the remaining ROI.As a result, the smooth non-uniformity field from Equation (17) in regions with opposite validity in the two images tends smoothly to unity farther from the valid region and towards the invalid regions of the image with I ROI,i (x) = 0.
The stability of the dynamic ranges along the iterations is imposed with two constraints.The first is in the spatial domain for the outputs of Equation ( 16).The pixel-wise average of gains in the restoration is set to unity with the normalization b The second constraint is in the statistics and uses the restoration field from Equation (17).The reference intensity η re f i,t along the iterations is constrained to be equal to that of the original image η re f i,t=0 .To achieve this, the final restoration fields are rescaled with W i,t ← W i,t × (η re f i,t=0 /η re f i,t ).That is, the reference intensities of v i,t = W i,t v i,0 are the same as those of v i,0 .

Results
The method was tested with BrainWeb brain images of the phantom simulator.It was also tested with real brain datasets of the Human Connectome Project (HCP) and with images of Parkinson's disease patients.

Implementation and Efficiency
The method was implemented in C++ as a command line program.It accepts a configuration file containing the filenames of the two images to be restored and the filenames of the images with the regions of interest that in this case are the corresponding brain masks.The program also accepts two significant command line arguments.These are the standard deviations of the size of the deconvolution filters, σ r i , and the standard deviation of the spatial smoothing filter of the restoration fields, σ s i , i = 0, 1.The remaining parameters of the program are optional.The program was compiled with gcc (version 5.3.0).The gcc compiler was part of a Cygwin environment interfacing Microsoft Windows 10 (Microsoft, Redmond, WA, USA).The Operating System (OS) was installed and the program was executed on an Intel Core i7 processor of 2.60 GHz and 64 bits.This processor has four cores and can accommodate eight threads.The processor is also combined with 8.00 GB of RAM.The C++ program can be used with appropriate compilers in other OSs as well as on different processors to be converted into binary.It is thus a multi-platform implementation.
The parameters were identical for all datasets of the phantom as well as the volunteer image groups.The co-occurrence statistics are of 3D neighborhoods N x of size ρ = 6 pixels with subsampling at regular intervals ∆ρ = 2 pixels along all axes.The dynamic range is limited to restrict the size of the co-occurrence matrices and the cost of their restoration filtering.The low value for the neighborhood size, ρ, allows the setting of a low value for the angular distortion, σ φ i = 4 • .The first variable parameters of the method are σ r 0 = σ r 1 .The parameters σ r i , i = 0, 1 are expressed as a fraction of the dynamic ranges and set to 0.02, also The constants of the sigmoid are k 1 = 6/15 and k 2 = −6.These parameters are used for the deconvolution filters with regularization 2 = 0.01.The parameters of the Van Cittert deconvolution filtering are β = 0.3 with a total of four iterations.The second variable parameters are the standard deviations of the spatial Gaussian filters σ s 0 = σ s 1 that are set to 140 pixels.The Gaussian smoothing filtering is performed separably.The iterative optimization allows using a value of σ r i that is an under-estimate and a value of σ s i that is an over-estimate in each iteration.The parameter controlling the iterations is t max = 20 for all datasets.

Description of the Phantom BrainWeb Brain Images
The first set was from the 1.5 Tesla BrainWeb brain MR simulator.The images were of the most commonly used anatomic brain MRI contrast mechanisms, T 1 weighted, T 1 w, and T 2 weighted, T 2 w.It consists of eight representative pairs of T 1 w and T 2 w images of the BrainWeb simulator [37].Both images in a pair have a resolution of 1.0 × 1.0 × 1.0 mm 3 and a matrix of size 181 × 181 × 217.They are corrupted with simulated non-uniformities of levels b = 0% −20% −40% −60% −80% −100% and a noise of n = 5% as well as an image with b = 40% and a noise of n = 3%.Only the brain region from the entire available head region was used.The brain region of the BrainWeb images is available through the union of the tissue type classification images.

Validation Measures for the Phantom BrainWeb Brain Images
The restorations of all the BrainWeb phantom brain images were evaluated with a measure of the contrast between the intensity statistics of the gray matter and of the white matter tissues.The contrast was quantified with the Coefficient of Joint Variation (CJV) [35,44,45]: where µ GM and µ W M are the mean values of the two tissues intensities, and σ GM and σ W M are their standard deviations.This measure represents the contrast between the intensity distributions of the GM and of the WM tissue regions.The statistics and CJV in Equation ( 18) are considered separately for each of the single contrast images T 1 w and T 2 w to give CJV T 1 and CJV T 2 , respectively.The ratio of the CJV of the restored CJV v t,i to the one of the corresponding original CJV v 0,i image is also computed , for i = T 1 , T 2 .A low value for the CJV v t is desired and thus a value for the ratio CJV ratio v t below unity indicates effective restoration.The BrainWeb brain phantom also makes available the original undistorted templates u i both for the tissue types and for the intensity values.The validation for these datasets is also performed by considering the absolute value of the difference between the restored images and the corresponding underlying undistorted templates di f f v t,i = v t,i − u i .This difference is desired to be kept at a low value.The ratio of the difference from the restored v t,i − u i to the difference from the original is also computed.It is desired to have a low value below unity for improved effectiveness.

Experiments with the Phantom BrainWeb Brain Images
The restorations of the brain part in the BrainWeb phantom images was evaluated with the intensity statistics of the GM and of the WM regions using the CJV given in Equation (18).The ground truth value of the CJV was computed from the regions occupied by each of the two tissues available from the BrainWeb MR simulator [37,46].The CJV was computed for the T 1 w images to give CJV T 1 and for the T 2 w images to give CJV T 2 .Their values are given in Table 1.The Table shows the values for the original corrupted image pairs and the values for the intensity restored images.In parentheses is the ratio CJV ratio t,i , i = T 1 , T 2 .The performance of the original and of the co-occurrence restored images for non-uniformities of very low levels, 0-20%, is comparable.The improvements resulting from the restoration methods are more apparent for higher intensity non-uniformity levels.The image noise increases the CJV.The difference between the corresponding original noise free and non-uniformity free phantom image and the restored ones di f f t,i are also computed and shown in Table 2.In parentheses is the ratio di f f ratio t,i , i = T 1 , T 2 .In all cases, its value is equal to or less than unity that shows that the restoration is effective.The duration of the restoration for a BrainWeb image pair lasts on average approximately 2 h 13 min.The restorations of the Brainweb phantom images with noise n = 5% and with the highest level of non-uniformity, b RF = 100%, are shown in Figure 5.This figure contains sections from the original and from the restored images as well as the co-occurrence statistics.In this example, the cerebellum in both the T 1 w and the T 2 w images becomes brighter and thus its statistics become closer to those of the corresponding mean tissue statistics over the remaining image regions.
The proposed restoration method has been compared with the N4ITK tool provided through Slicer3D (version 4.6.2) [38,39].The default N4ITK parameters in Slicer3D were used to restore the BrainWeb images.Similarly to the proposed method, the brain regions were provided to N4ITK.The N4ITK restoration is very weak for all images.The restoration for low levels of non-uniformity, 0-60%, leads to a significant loss of contrast for both the T 1 weighted images and the T 2 weighted images.This leads to an increase in the CJV and to CJV ratio > 1.In high levels of non-uniformity, 80% and 100%, the restoration method with the N4ITK methodology for the T 1 weighted images decreases contrast less based on the CJV.The restoration of the T 1 w BrainWeb image with non-uniformity of 100% with the proposed method gives CJV ratio i=T 1 ,Co−oc.= 0.63 as shown in Table 1, whereas N4ITK gives CJV ratio i=T 1 ,N4ITK = 2.91.That is, the performance of the proposed method is 4.62 times better than that of N4ITK.The restoration of the T 2 weighted images also remains weak for high levels of non-uniformity.The restoration of the T 2 w BrainWeb image with non-uniformity of 100% with the proposed method gives CJV ratio i=T 2 ,Co−oc.= 0.65 as shown in Table 1, whereas N4ITK gives CJV ratio i=T 2 ,N4ITK = 7.67.That is, the performance of the proposed method for this image is 11.8 times better than that of N4ITK.Thus, the restoration of the T 2 weighted images with the N4ITK tool performs poorly for all levels of non-uniformity.This is due to the low contrast between the distributions of the GM and of the WM in the T 2 weighted BrainWeb images.The results based on di f f show a similar poor performance for N4ITK.Thus, in all cases, the proposed method performs significantly better than N4ITK.
The effects of the restoration methods on the extent of the occupied dynamic ranges of the images are also measured.The percentage of change of the size of the dynamic ranges is measured for every image restoration and for both methods.The average value over the eight images for each contrast and for each method is computed.The value of the average with the proposed method for the T 1 w images is −1.15%, and, for the T 2 w images, it is −9.78%.The non-zero changes are due to the noise removal from the images prior to the standardization of the dynamic ranges.The decrease is also due to the sharpening of the distributions at the high intensity parts of the dynamic ranges.The decrease for the T 2 w images is greater than for the T 1 w images.This is because, in the T 2 w images, the tissues distributions correspond to higher intensity ranges.The above values are expected and demonstrate the stability of the proposed method.The average changes in the dynamic ranges with the N4ITK tool were also measured.The average change with N4ITK for the T 1 w images is +29.77%,and, for the T 2 w images, it is −9.75%.The significant positive increase in the size of the dynamic range for the T 1 w images demonstrates an instability for the N4ITK method.

Description of the Real Images
Anatomic datasets from two studies of the Human Connectome Project (HCP) [47][48][49][50] were used.The first was anatomic data from the Lifespan (LS) pilot study and the second was anatomic data from the Retest study.The data of the HCP LS study was from 27 volunteers from six age groups 4-6, 8-9, 14-15, 25-35, 45-55, and 65-75 years old.The second was data from the HCP Retest study.It was data from 45 aged volunteers, 31 women and 14 men.The volunteers were imaged at 3.0 Tesla (Siemens, Connectom and Prisma).A T 1 weighted 3D structural MPRAGE sequence was acquired sagittally with TR/TE/TI = 2400/2.14/1000ms.A T 2 weighted 3D structural SPACE sequence was acquired sagittally with TR/TE = 3200/565 ms.The voxel resolution of both the T 1 w and the T 2 w images is 0.8 × 0.8 × 0.8 mm 3 .The matrix size of both the T 1 w and the T 2 w images is 208 × 300 × 320.The youngest group of age range 4-6 years old of the HCP LS pilot project was scanned using a customized pediatric head coil and an optimized protocol for that age range.
Another real dataset was of brain images of patients at an advanced stage of Parkinson's disease.The datasets were retrieved retrospectively from the clinical database of University Hospital Goettingen for pre-operative assessment before implantation of electrodes for deep brain stimulation.The use of the patient data was retrospective, fully anonymized, and according to the guidelines of the local ethical committee for clinical research.The patients gave informed consent for all imaging procedures.All pre-operative images were acquired under anaesthesia to reduce motion artifacts and increase anatomical precision.The quality of the images was evaluated by observation and images with artifacts were removed.The data were from 60 patients, 41 men of average age 59.91 ± 7.51 years and 19 women of average age 65.86 ± 7.40 years.
The Parkinson's disease patients were imaged at 3.0 Tesla (Siemens, TrioTim).A standard T 1 weighted 3D structural MPRAGE sequence was acquired sagittally with TR/TE/TI = 2000/2.98/900ms that gave an in-plane resolution of 1.0 × 1.0 mm 2 and a slice thickness of 1.1 mm.The matrix size of the T 1 w images is sufficient to cover the entire head region and of size at least 240 × 256 × 160.
A high-resolution 3D SPACE structural T 2 weighted imaging that was axially planned was also acquired with TR/TE = 1500/355 ms that gave an in-plane resolution of 0.63 × 0.63 mm 2 and a slice thickness of 1.80 mm.The size of the images was 308 × 384 × 80.In the high resolution T 2 w, the caudal cerebellum is not always depicted, depending on the skull size and "brain fitting" in a standard protocol for all Parkinson's disease patients.The primary target in Parkinson's disease diagnostic is the high resolution imaging of the mesencephalon and of the basal ganglia.The in-plane orientation of the two sequences are complementary to enable a more complete clinical evaluation.
The brain regions for the real images were used by the method and were extracted with the BET tool [51].The two brain datasets were placed in the same resolution and in the same sampling grid.The reference was the T 1 w image and the T 2 w image was resampled with a closest neighbor filter [38].The images were smoothed with a median filter of size 3 × 3 × 3 mm 3 along the axes.

Validation Measure for the Real Brain Images
The validation measure uses the entropy of the histogram of the original image, H orig , and the entropy of the histogram of the restored image, H rest .The actual measure is the ratio The entropies are expected values of exponents of densities of intensities.Thus, in Equation ( 19), they are used as exponents of e.A successful restoration decreases the value of the entropy of an image.Thus, both the value of the numerator in Equation ( 19) and the entire ratio, H ratio exp , become negative.The successful restoration also means improved tissues contrasts.

Experiments with Brain Images of the Human Connectome Project
Table 3 gives the statistics for measure H ratio exp for the T 1 w and for the T 2 w images for all 27 volunteers of the Human Connectome Project (HCP) Lifespan (LS) dataset.Table 4 gives the statistics for measure H ratio exp for the T 1 w and for the T 2 w images for all 45 volunteers of the Human Connectome Project (HCP) Retest dataset.The tables show that H ratio exp are indeed negative in all cases.Thus, the restorations are successful for all images.The duration of the restoration for an HCP image pair lasts on average approximately 6 h 53 min.The restoration of representative T 1 w and T 2 w images of an HCP volunteer are shown in Figure 6.The GM tissues in the cortex, subcortical regions, and the cerebellum become of more uniform intensity in both images.The intensities of the WM also become more uniform.The statistics of both the T 1 w and the T 2 w images in Figure 6 show different distributions corresponding to the three tissues, even along the diagonal of the auto-co-occurrences that are minimally involved in the restoration.In the T 1 w statistics, there are separate distributions corresponding to the WM and to the GM.The distribution corresponding to the CSF appears.They also show a sharper distribution corresponding to the border between the WM and the GM.In the T 2 w image, the different distributions of the WM and of the CSF become apparent.The consideration of both images shows their successful restoration.The statistics of the initial T 2 w image have a heavy tail at the high intensity part of the dynamic range from contributions of the GM and of the CSF.Thus, the standardization of that part of the dynamic range to a limited range results in its compression.That compression causes a shift of densities from the lower intensity range to a higher intensity range.

Experiments with the Brain Images of Parkinson's Disease Patients
Table 5 gives the statistics for measure H ratio exp for the T 1 w and the T 2 w images for all 60 patients.The table shows that H ratio exp is negative in all cases.Thus, the restorations are successful for all images.The duration of the restoration for an image pair of a Parkinson's disease patient lasts on average approximately 2 h 12 min.The restoration of representative T 1 w and T 2 w images with extensive non-uniformities of a Parkinson's disease patient are shown in Figure 7.The GM tissues in the cortex, subcortical regions, and the cerebellum become of more uniform intensity both in the T 1 w and in the T 2 w images.The intensity of the WM also becomes more uniform.The statistics of the T 1 w image in Figure 7 show three different distributions corresponding to the three tissues.The diagonal auto-co-occurrences, even though they are minimally involved in the restoration, still consist of three separate distributions.The statistics of the T 2 w image also show distinct distributions for the WM and the GM even though they have a low contrast.They also show a sharper distribution corresponding to the border between the WM and the GM regions.The distribution corresponding to the CSF regions appears.This shows the successful restoration of the T 2 w image.The presence of separate tissues distributions in the statistics of the restored real images gives the entropy based measure a semantic meaning in terms of tissues intensity uniformities and tissues contrasts.
The HCP LS images in Figure 6 and the Parkinson's disease images in Figure 7 show that the sub-cortical regions can have a variable intensity non-uniformity effect.The restoration removes these non-uniformities.The mesencephalon and sub-cortical regions called basal ganglia and substantia nigra are implicated in Parkinson's disease.The restoration of the uniformity in these brain regions can improve the analysis to characterize the appearance of Parkinson's disease in MRI data [52,53].

Discussion
The Bayesian coring formulation is general and non-parametric.In this method, it is used with the auto-and the joint-co-occurrence statistics.These are higher order statistics that favor the dominant distributions of the data and decrease their spread compared to those of the intensity histogram.A sigmoid applied to the auto-co-occurrence statistics removes the dominance of their diagonal and hence increases the sensitivity to the statistics corresponding to the borders between regions.Thus, the statistics improve the discriminability between the distributions of the dominant tissues.In brain images, the extensive interface between the gray matter and the white matter tissues gives rise to one of the dominant joint distributions.
The effect of the spatial non-uniformity field on the co-occurrence statistics is modeled as a bimodal point spread function.One of the modes corresponds to multiplicative factors above unity for the regions that become brighter and the other mode corresponds to multiplicative factors below unity for the regions that become darker.The distortion filter is assumed to be non-stationary.It is used in the Bayesian formulation as both the likelihood and to compute the prior.The bimodal distribution for the distortion can be more robust to the presence of even higher field intensity non-uniformities.
The method also standardizes and preserves the valid dynamic ranges along the iterations.Thus, it offers stability and efficiency.The statistics are global over the entire ROI and this enables the uniform restoration of spatially loosely connected or even disconnected regions of the same tissue such as of the different brain gray matter centers.The method also accommodates the inevitable difference between the signal regions of the two images of different contrasts.
The restoration of MRI images for artifacts similar to the one in this work typically involve the iterative optimization of a parametric energy function.Such methods involve the computation of the histogram, its back-projection to image space, and the smoothing of the entire spatial restoration field.These steps are typically repeated in every iteration for the smoothing coefficient of each spatial basis function such as a b-spline basis.This is very intensive computationally.The method developed in this work is also iterative.However, the histogram computation and the smoothing of the entire 3D spatial field is only done once per iteration.This is significantly more efficient computationally.
The efficiency of the implementation of the methodology can be improved by taking advantage of parallelization within and between its various steps.This can be achieved by combining the C++ implementation with multi-threading.The duration of the restoration also depends on the magnitude of the non-uniformities present in the images of the pair.
The method was tested with images of the BrainWeb brain simulator.It was also tested with brain data from the Human Connectome Project (HCP) as well as data from a database of brain images of elderly patients with advanced Parkinson's disease.Image pairs of T 1 w and T 2 w images have been restored jointly for intensity uniformity.The method restores the images successfully, despite the high intensity distortions that are present in most of them.The Parkinson's disease patients images have a lower tissue contrast and some also contain lesions.The effectiveness of the method for this data demonstrates its efficiency.The joint restoration is mutually beneficial for both images.The proposed restoration method has been compared with the N4ITK tool.The N4ITK showed a very poor performance.The co-occurrence method performs significantly better than N4ITK.In addition, this method preserves and restores the dynamic ranges of the images.The N4ITK tool can expand the dynamic ranges of the T 1 w images that shows its instability.
A simplified model of the physical effects of the non-uniformities is assumed to design this methodology.This can lead to limitations in specific contexts.The method smooths the spatial non-uniformity field with a filter of a finite size.This allows a spatial non-uniformity of a certain magnitude.In images where the spatial non-uniformity has a lower magnitude than that of the allowed non-uniformity, the restoration may result in a loss of contrast between the tissues statistics.A complication is that the MRI physical non-uniformity also depends on the signal properties of the tissues.This is effectively ignored in this non-parametric methodology.Thus, the physical non-uniformity cannot be fully recovered [19].An MRI artifact that is currently ignored is the effect of the radio-frequency non-uniformity on stimulated echoes that can result in non-uniformities of the contrast, such as in the T 2 w images.However, the restoration method in this study, as it emphasizes the intensity transitions, is more robust to the latter artifact.
The design of the methodology also makes several assumptions to improve efficiency.The restoration assumes that the removal of the noise is independent from the removal of the non-uniformity.However, the regions that become darker as a result of the restoration have lower contrast to noise ratio.Thus, the joint consideration of non-uniformity and noise may lead to an improved restoration.The method also requires the field of view for the two images to be provided.This may be more challenging for images with complicated pathology such as brain tumors.In these cases, the restoration methodology may have to be combined with other methods such as segmentation.The imaging of a brain tumor pathology may also benefit from specific imaging sequences beyond only T 1 w and T 2 w.

Conclusions
In conclusion, a method has been developed that uses local intensity co-occurrences across regions for an effective joint restoration of two anatomic MRI images from intensity non-uniformities.The restoration is efficient and robust with respect to contrast mechanisms as well as subject geometries.Thus, it decreases the need for calibration scans and can reduce acquisition time.The method can be used to improve subsequent processing steps such as image segmentation.

Figure 1 .
Figure 1.The bimodal shape of the PSF, P b (•), of the distortion of the statistics.

Figure 2 .
Figure 2. Overview of the non-parametric Bayesian formulation for image restoration.The restoration is repeated iteratively.

Figure 3 .
Figure 3. Overview of the series of steps in the implementation of the Bayesian formulation for image restoration.The restoration is joint for two images and iterative.

Figure 4 .
Figure 4.The auto-co-occurrence statistics and the joint-co-occurrence statistics of a T 1 and a T 2 image of the BrainWeb phantom without non-uniformity, b = 0%, and with noise of n = 5% [37].The densities in the statistics are displayed in logarithmic scale.The individual distributions of the Gray Mater (GM), White Matter (WM), and Cerebrospinal Fluid (CSF) are apparent.
up to an intensity that corresponds to 300% of η re f i to give the maximum intensity η max i and the range (1.5 × η 0.9 i , 3.0 × η 0.9 i ].The intensity ranges of [η min i back-projected to the respective images ROIs to give the valid signal regions.The regions in image I ROI with intensities in ranges [0, η min i ] and with intensities beyond η max i are considered invalid.The auto-co-occurrences C ii in Equation (11) and the joint-co-occurrences C 01 are computed at least over the valid joint ranges ([η min i , η max i ], [η min j , η max j ]

Figure 5 .
Figure 5.The restoration of a T 1 w and a T 2 w BrainWeb image pair with non-uniformity of 100% and noise of 5%.The restoration makes the cerebellum brighter and the statistics sharper.

Figure 6 .
Figure 6.Example restoration of a T 1 w and a T 2 w image pair for the HCP LS dataset.The intensities of the white matter and the gray matter in both the T 1 w image and in the T 2 w image become more uniform.The statistical distributions become sharper.

Figure 7 .
Figure 7. Example restoration of a T 1 w and a T 2 w image pair of a Parkinson's disease patient.The intensities of the white matter in both the T 1 w image and in the T 2 w image become more uniform.The statistical distributions become sharper.

Table 1 .
Validation for BrainWeb phantom T 1 w and T 2 w images with CJV i of GM and WM tissue regions.A low value indicates improved performance.In parentheses is the ratio of the restored to the original, CJV ratio i .Low values and less than unity indicate improved performance.

Table 2 .
Validation for BrainWeb phantom T 1 w and T 2 w images with difference to underlying anatomic images di f f t,i .A low value indicates improved performance.In parentheses is the ratio of after to before the restoration, di f f ratio i .Low values and less than unity indicate improved performance.

Table 3 .
Statistics of H ratio exp for the T 1 w and the T 2 w images for the 27 HCP LS volunteers.The H ratio exp are significantly negative for all images and hence all the restorations are successful.

Table 4 .
Statistics of H ratio exp for the T 1 w and the T 2 w images for the 45 HCP Retest volunteers.The H ratio exp are significantly negative for all images and hence all the restorations are successful.

Table 5 .
Statistics of H ratio exp for the T 1 w and the T 2 w images for the 60 Parkinson's disease patients.The H ratio exp are significantly negative for all images and hence all the restorations are successful.