A Clustering Approach to Improve IntraVoxel Incoherent Motion Maps from DW-MRI Using Conditional Auto-Regressive Bayesian Model

: The Intra-Voxel Incoherent Motion (IVIM) model allows to estimate water diffusion and perfusion-related coefﬁcients in biological tissues using diffusion weighted MR images. Among the available approaches to ﬁt the IVIM bi-exponential decay, a segmented Bayesian algorithm with a Conditional Auto-Regressive (CAR) prior spatial regularization has been recently proposed to produce more reliable coefﬁcient estimation. However, the CAR spatial regularization can generate inaccurate coefﬁcient estimation, especially at the interfaces between different tissues. To overcome this problem, the segmented CAR model was coupled in this work with a k -means clustering approach, to separate different tissues and exclude voxels from other regions in the CAR prior speciﬁcation. The proposed approach was compared with the original Bayesian CAR method without clustering and with a state-of-the-art Bayesian approach without CAR. The approaches were tested and compared on simulated images by calculating the estimation error and the coefﬁcient of variation ( CV ). Furthermore, the proposed method was applied to some illustrative real images of oncologic patients. On simulated images, the proposed innovation reduced the average error of 47%, 21% and 58% for D , f and D ∗ , respectively, compared to the state-of-the-art Bayesian approach, and of 48% and 34% for D and f , respectively, compared to the original CAR, while it achieved the same error for D ∗ . The clustering approach was also able to consistently reduce the CV for each coefﬁcient. On real images, the novel approach did not alter the IVIM maps obtained by the original CAR method, with the advantage of reducing their typical blotchy appearance at the boundaries. The proposed approach represents a valuable improvement over the state-of-the-art Bayesian CAR method and provides more reliable IVIM coefﬁcient estimation, and is less sensitive to bias and inconsistency at tissue/tissue and tissue/background interfaces.


Introduction
Intravoxel Incoherent Motion (IVIM) Magnetic Resonance Imaging (MRI) is a type of Diffusion-Weighted MRI (DW-MRI) acquisition, firstly proposed by Le Bihan et al. [1], for estimating diffusion and perfusion tissue properties.In the IVIM mathematical representation, a bi-exponential function of the diffusion signal decay at different b-values is used to simultaneously characterize tissue diffusivity and perfusion within capillaries.Specifically, the IVIM model allows us to estimate the diffusion (D) and pseudo-diffusion coefficients (D * ), along with the perfusion volume fraction ( f ) within the tissue.
IVIM has gained particular interest in oncologic applications, especially in the last decade, thanks to the advancements in MRI hardware, pulse design and post-processing methods [2].In fact, recent works have proposed the adoption of the IVIM model to improve tumor diagnosis [3], prognosis [4] and response to treatment [5].Furthermore, IVIM has also been used to characterize radiation-induced modifications in normal tissues after radiotherapy [6].Nevertheless, the use of the IVIM technique is not confined to oncology, but it has significant potential in clinical and preclinical applications for neurological diseases [7] and other physio-pathological conditions [8,9].
However, IVIM coefficient estimation presents several issues related to the presence of noise over the image, the correct model identification and the amount of tissue perfusion [10].Currently, the most popular techniques to fit IVIM data are (i) the non-linear least-square (LSQ) fitting, based on the Levenberg-Marquardt algorithm [11], and (ii) the Bayesian approach, which provides the posterior probability densities of the coefficients based on prior knowledge [12].
The Bayesian method showed superior performance in reducing errors and variability [12], and also allowed the inclusion of a spatial dependence between voxels, to obtain more regular maps.This is the case of the spatial homogeneity prior [13] and the Conditional Auto-Regressive (CAR) prior specification, which takes into account information from the neighbor voxels [14,15].In particular, it was reported that the CAR specification on a segmented Bayesian method improved the estimation accuracy compared to both LSQ and traditional Bayes approaches, especially for D * [15]; however, a limitation of the CAR method lies in the choice of the neighborhood, as it attributes the same weight to all surrounding voxels without considering the continuity of the structure.This leads to a slight decrease in the quality of the IVIM coefficient maps, especially at the tissue/tissue and tissue/background interfaces, characterized by a blotchy appearance.
The goal of this work was to overcome this limitation by coupling the CAR prior with a k-means clustering approach, which distinguishes the different structures based on the diffusion signal decay along the b-values, and uses this information to set different weights to the contributions of the different neighbor voxels.This technological improvement was aimed to obtain more regular and reliable IVIM parametric maps.The voxel-based estimation of IVIM coefficients can be used in the clinical practice for the characterization of microstructural properties of diffusion and perfusion within subregions of organs/structures.This information can be useful, for example, to define prediction models for treatment response in oncology or to assess changes occurring in physiological and pathological conditions in different applications; therefore, the availability of high quality and accurate maps represents a paramount condition in clinical practice and improved quality maps may have, prospectively, an impact on the quality of care The proposed method was compared to the original CAR approach and a state-of-theart Bayesian algorithm in terms of accuracy and map homogeneity.First, it was applied to a set of simulated phantoms for which the true value of the IVIM coefficients was known.Then, it was applied to some illustrative real images of oncologic patients, to evaluate the estimated maps in a real clinical scenario.

IVIM Fitting Procedure
The decay of the signal intensity SI(i, j, b) in function of b is described in any voxel (i, j) by a bi-exponential function, according to the IVIM model [16]: where D(i, j), D * (i, j) and f (i, j) are the diffusion coefficient, the pseudo-diffusion coeffi- cient and the perfusion volume fraction in voxel (i, j), respectively.The goal of the fitting is to estimate these three coefficients in each voxel (i, j) of a region of interest, based on a set of signal intensity observations SI obs ijb in an image domain I at different b-values.

Clustering Approach
In the proposed approach, a clustering technique is employed to determine the cluster C ij associated with each voxel (i, j) ∈ I.Then, a Bayesian CAR fitting is performed, where the CAR prior specification depends on the clustering solution C ij , (i, j) ∈ I .
The k-means clustering approach is chosen for its simplicity and low computational cost, and because it is one of the most widely adopted clustering methods for medical image segmentation.It is employed by assuming as the distance between any pair of voxels the L1 distance between the vectors of their signal intensities over b ∈ B.Moreover, the number of clusters N C is set equal to 9, which is reasonable from a theoretical viewpoint, as it represents the combinations of two levels (low or high) for the three coefficients to be estimated plus an additional case for the background (2 3 + 1 = 9); therefore, it gives enough flexibility to recognize these macroscopically different behaviors.
The results of the clustering procedure is the cluster C ij ∈ [1, N C ] of each voxel (i, j) ∈ I, where all voxels belonging to the same cluster have the same C ij value.
The clustering was implemented in MATLAB (R2020a, The MathWorks Inc., Natick, MA, USA), with the function "kmeans", imposing 50 replicates with alternative initial cluster centroid positions.

Bayesian CAR Fitting
This fitting procedure is taken from [15], where a segmented approach was considered by dividing the estimation procedure in two steps.
In the first step, the coefficients D(i, j) and f (i, j), where f denotes an intermediate estimation of f , are estimated with the observations at b = 0 and b ∈ B : b ≥ b, with b = 200 s/mm 2 .A simplified version of the IVIM model is used, which only includes the second exponential decay: In the second step, the coefficients D * (i, j) and f (i, j) are estimated using all observa- tions b ∈ B. The complete IVIM model (1) is used while fixing coefficients D(i, j) to the values estimated in the first step (the expected value of the marginal posterior density).
To perform the Bayesian fitting, (2) and (1) are processed to derive the likelihood function for the first and second step, respectively, as reported in the Supplementary Materials.This likelihood also depends on an additional parameter to be estimated, i.e., the standard deviation σ lik of the conditioned densities that determine it.
Then, a prior density for the unknown parameter is assumed.A priori independence between σ lik , each D(i, j) and each f (i, j) is assumed for the first step, and between σ lik , each D * (i, j) and each f (i, j) for the second step.
Separately for each λ (where λ generically denotes each coefficient D, D * , f or f ), a CAR specification is included in the prior by considering the intrinsic CAR formulation of Leroux et al. [17].More specifically, the following distribution for a coefficient in a voxel given the values in the other voxels is taken: where λ c i,j = λ \ {λ(i, j)}, and w α,β i,j denotes the spatial neighborhood matrix.We assume a dependency (w α,β i,j = 0) only for the voxels (α, β) bordering on (i, j).Their values are assigned based on the clustering solution C ij , (i, j) ∈ I .In particular, w α,β i,j = 1 if the bordering voxel (α, β) belongs to the same cluster of voxel (i, j) (i.e., C αβ = C ij ), or w α,β i,j = 0.1 if the bordering voxel (α, β) does not belong to the same cluster of voxel (i, j) (i.e., C αβ = C ij ); therefore: In addition, to smooth the autoregressive component, the actual prior density of each λ consists of a mixture model with two components: the above conditional autoregressive specification (3) and a truncated Gaussian marginal density N (λ min ,λ max ) µ λ , σ 2 0λ , where λ min and λ max denote the minimum and maximum admissible values, respectively.Therefore, the overall prior for each coefficient λ is a mixture between (3) and the Gaussian component, with weight 0.75 for the CAR specification (3) and weight 0.25 for the Gaussian component.The hyperparameters are detailed in Section 1 of the Supplementary Materials.
Finally, the prior of the standard deviation parameter σ lik is an independent Gamma density to respect its positivity, i.e., σ lik ∼ Gamma , 1 .The scale parameter is set equal to 1 to give a large density, while the mean value is set to half of the average SI obs ij0 over (α, β) ∈ I.
The approach was coded in R with package STAN [18], which directly provides the posterior marginal densities of each parameter.
The whole procedure is sketched in Figure 1.The clustering processes the signal intensity SI(i, j) to determine the clusters C ij , which are used in the next steps to estimate the weights in the CAR prior specification.The input of both steps of the Bayesian CAR fitting are thus the output of the clustering, together with the signal intensity SI(i, j) itself.The first step produces the coefficients D(i, j) and f (i, j), which are passed as input to the second step.The second step produces the coefficients D * (i, j) and f (i, j).

Evaluation 2.2.1. Simulated Data
We applied our method on simulated DW-MRI signals with different noise levels to test the performance of the proposed approach.To increase the complexity of the simu-lated dataset, a Shepp-Logan phantom (matrix size 64 × 64) with three regions of interest (ROIs) along with background noise was used as a template for our simulations.Each ROI was characterized by an IVIM coefficient triplet randomly generated from uniform distributions and a specific intensity at b = 0.More specifically, D was sampled in the interval 0.0005 mm 2 s , 0.002 mm 2 s , D * in the interval 0.005 mm 2 s , 0.1 mm 2 s and f in the interval [0.025, 0.4], which include typical values of tumor and normal tissues.Images were gener- ated using (1) for two different b-values sets, i.e., 0, 25, 50, 75, 100, 150, 300, 500, 800 s mm 2 and 0, 25, 50, 75, 100, 300, 600, 1000 s mm 2 .Rician noise was finally added to create the final DW-MRI images with different Signalto-Noise Ratios (SNRs), equal to 10, 25, 50, 100, 150 and 200.Five numerical phantoms for each combination of b-value set and SNR were generated, for a total of 60 simulated cases.
Numerical phantoms were generated using a MATLAB custom-made script.

Real Data
DW-MRI images of five Head-and-Neck (HN) patients affected by oropharyngeal tumor (three males and two females) and four pelvic patients affected by rectal tumor (two males and two females) were considered as examples of IVIM acquisitions for tumor diagnosis.They were collected from retrospective studies approved by the local ethics committee and already published in previous works [14,15,19].
DW-MRI images were acquired using a 1.5 T system (Optima MR 450w, GE Healthcare, Milwaukee, WI, USA) and obtained by single-shot spin-echo echo-planar imaging with: acquisition matrix equal to 128 × 128; in-plane resolution equal to 1.094 × 1.094 mm 2 ; TR/TE equal to 4500 ms/72 ms for HN and 3500 ms/77 ms for pelvis; slice thickness equal to 4 mm for HN and 5 mm for pelvis; same b-values as in the simulated images.The estimated average SNR of the acquired images was equal to 80 and 50 for the HN and pelvis, respectively; in this way, we were able to test the performance of the methods at different SNR values also in the real setting.
Three different tissues were considered in the two datasets: masticatory muscle, parotid gland and tumor for HN; internal obturator muscle, prostate and rectal tumor for pelvis.An average of 4-5 slices containing the three structures of interest were selected for each patient, excluding basal and apical slices to avoid partial volume effects, for a total of 39 two-dimensional images (21 for HN and 18 for pelvis).In the two female pelvic patients, one structure (the prostate) was not considered.

Comparison with State-of-the-Art Methods
Maps estimated with the proposed approach were compared to the CAR approach without clustering, already published in [15], and to a state-of-the-art Bayesian approach without spatial regularization [20].As for the CAR without clustering, it follows the description reported in Section 2.1.2but neglecting the clustering solution.Indeed, the weights w α,β i,j of the CAR prior specification are set equal to 1 for all voxels (α, β) bordering on (i, j), and 0 for the others: As for the standard Bayesian approach, we adopted the implementation provided by Gustafsson et al. [20], considering the default priors and the hyperparameters proposed by the authors.The IVIM signal was fitted using the bi-exponential decay in (1) and considering the whole set of b-values.The three IVIM coefficients and each signal SI(i, j, 0) at b = 0 were estimated under the following priors [20]: a uniform prior for f and each SI(i, j, 0); a lognormal prior for D and D * .Moreover, all prior distributions were truncated to the following intervals: 0, 3 • 10 −3 mm 2 s for D; 0, 0.15 mm 2 s for D * ; [0, 1] for f ; [0, 1] for SI(i, j, 0) being the simulated images in this range.The parameter estimates were given in terms of the mean value from the respective marginalized posterior distribution.
Following the procedures already adopted in previous works [10,14,15], the accuracy of each method on simulated images was evaluated in a voxel in terms of the error with respect to the true value: where λ est is the estimated value of the coefficient and λ true the corresponding true value.
The median of the errors over the voxels of a ROI was finally calculated for each coefficient.Moreover, the noisiness of the map was evaluated in terms of the coefficient of variation (CV), calculated as the ratio between the standard deviation StD(λ est ) and the mean value Mean(λ est ) of the estimated coefficients in a ROI with constant λ true : The average value of the CV over the ROIs was then considered.
A statistical analysis was performed to evaluate the differences among the three approaches.A Kruskal-Wallis test was applied to ε and CV estimated with the three methods for each SNR value, in order to identify significant differences among groups.For those groups presenting significant differences, a Tuckey post hoc test was next performed to determine which couples were significant.Significance was set to p < 0.05, and Bonferroni correction for multiple comparisons was applied.
IVIM maps estimated on real cases were qualitatively evaluated and mean values of IVIM coefficients were calculated within each structure.

Simulated Data
The qualitative comparison of the estimated maps among the three different methods revealed that the CAR clustering approach was able to improve their quality with respect to the Bayesian CAR method without clustering and the state-of-the-art Bayesian approach without CAR (Figure 2).
The Bayesian estimation without CAR generated maps with high contrast between different structures, but maintained a certain level of noise throughout the image, especially at low SNR values.On the contrary, the CAR approach without clustering was able to generate homogeneous values within the ROIs, but at the same time it lost the boundaries between the regions and introduced blurring and blotchy artifacts, especially at low SNR values and for D * .Instead, the proposed approach with clustering better highlighted the contrast between structures at the interfaces, and preserved a high homogeneity within the regions; therefore, it retained the advantages of both other methods.
The error ε maps reported in Figure 3 make it evident that the high errors at the boundaries of the structures in the CAR maps were almost completely eliminated with the clustering approach, especially for high SNR values.In addition, it can be seen from the maps that, in both CAR approaches, ε values were highly dependent on the IVIM coefficient values (i.e., some specific triplets were more difficult to be estimated), whilst the Bayesian method without CAR seemed less dependent.In particular, it showed for D the same error level throughout the whole image.Figure 4 shows for each SNR the barplots of ε and CV averaged on the simulated images with the corresponding significance, while the detailed results for each specific simulation are reported in the Supplementary Materials Tables S1-S4.Similar trends of ε and CV can be observed for the three methods, with higher values at low SNR that decrease at high SNR, as expected.The accuracy was generally higher under the proposed approach with clustering as regards the estimation of D and f : ε on D presented an averaged reduction of 46% compared to the other approaches, which was significant for SNR = 10, SNR = 100 and SNR = 200; ε on f reduced of 21% compared to the Bayesian approach without CAR and of 34% compared to the CAR without clustering, although this was not statistically significant.As for D * , similar errors were obtained for both CAR approaches, lower than for Bayesian approach without CAR (non-significant reduction of about 58%).CV was always the lowest for the proposed approach with the clustering, for each IVIM coefficient and SNR value.Specifically, the introduction of clustering reduced the CV of about 25-55% compared to traditional Bayes, and of about 30% compared to the original CAR approach.Furthermore, for D * maps at low SNR, the two CAR approaches presented similar CVs, and a significant reduction was found only between the original Bayesian approach without CAR and the CAR with clustering.On the contrary, for high SNR the original CAR without clustering and the traditional Bayesian approach were almost comparable, whilst the introduction of clustering significantly reduced the CV compared to both of them; therefore, the CAR with clustering was able to improve not only the quality of the maps, especially on the boundaries of the structures, but also the estimation of the IVIM coefficients with respect to the other methods.

Real Data
The IVIM maps estimated on an illustrative slice of a HN and a pelvic patient using the three different methods are shown in Figure 5 and Figure 6, respectively.The maps estimated with the clustering were almost similar to those estimated with the CAR method without clustering.The most relevant differences lie in the more defined boundaries of contrasted structures (e.g., the external boundary in the HN case, and the bladder wall and cortical bone of the femoral head in the pelvic case) and in the reduced blurred effects.D, f and D * values estimated by the three methods and averaged over the HN and the pelvic patients are reported in Table 1.As expected, the variations between the original CAR and that with clustering were limited, whilst the Bayesian approach without clustering showed different values, especially for D * .

Discussion and Conclusions
In this work, we have proposed an improvement to the Bayesian CAR fitting of the IVIM model [15] by embedding a k-means clustering procedure to specify the weights in the CAR prior specification.The aim was to improve the estimation of the parametric maps for the IVIM coefficients.In fact, in the previous work [15], it was highlighted that the main limitation of the CAR approach was in the way the neighbor voxels were considered in the fitting procedure, since all the eight surrounding voxels were included with the same weight, thus neglecting the heterogeneity of the surrounding tissues.
The simplest solution to this issue from a theoretical point of view is to also estimate all weights of all voxels by means of the Bayesian approach; however, this approach proved to be ineffective due to the huge number of parameters to be estimated and the too high computational effort.In fact, this is not pursued in the literature, while only examples are available in other applications in which the weight of the CAR component is tuned as a whole without distinguishing the contributions of the specific elements in the neighborhood [21].
To overcome this issue, the weights have been provided by embedding a k-means clustering approach in the procedure, based on a distance function that considers the whole set of images at different b-values.This allows us to group the voxels with similar exponential decay and to give them higher weights.It is worth noticing that a simpler clustering approach based only on the anatomical information brought by the image at b = 0 was not sufficient to reduce the blurring effects, since it is not able to capture tissuespecific differences in the diffusion signal decay.This also shows that a manual clustering is rather impossible for this problem, as it should require to analyze "as a whole" the entire set of images at different b-values, which is a very difficult task for a human operator.In fact, no work in the literature considers this manual procedure on diffusion signal dynamic.
As the k-means requires to enter the number N C of clusters as input, we performed a preliminary sensitivity analysis on a subset of data.We found that the clusters usually remain quite invariant within the region of interest, while additional clusters are identified in the background.In the real images, we also found that for N C > 9 different clusters C ij are assigned to voxels in regions disjoint from each other, which belong to the same cluster for N C = 9, while the shape of the connected regions in the clusters does not change.Since the CAR prior is not affected by this difference, we have chosen to limit the number of clusters to 9. On the contrary, on real images, a lower number of clusters may not be able to satisfactorily reduce the blotchy effect produced by the original CAR approach.
Results on simulated data have confirmed that the improved approached based on the clustering is able to face the CAR limitations, showing less noisy IVIM maps and sharper boundaries definition at the interfaces between structures, associated with reduced errors at the boundaries.In addition, comparable results were obtained from simulations generated at different b-values, considering both HN and pelvic setups.
The application to some illustrative real oncologic cases has highlighted that the introduction of the clustering did not alter the IVIM maps estimation obtained by the original CAR method, with the advantage of reducing the typical blotchy appearance at the boundaries, which usually affects the CAR approach, especially at low SNRs and for the estimation of D * .
The proposed method was able to generate IVIM maps with lower errors and higher quality compared to the original CAR and also to the standard Bayesian approach.Compared to other approaches that adopted a spatial homogeneity prior [13], the results with the CAR approach were similar in terms of map quality and accuracy, as already discussed in [15].Even if a fair direct quantitative comparison between these works is not possible, it seems that the introduction of clusters in the CAR prior may give more regularized maps than those presented in [13], at least for simulated images.
Recently, approaches based on deep neural networks (DNNs) have been introduced to estimate IVIM parametric maps [22][23][24].In these works, it was reported a higher quality of the maps compared to those estimated using traditional least squared and Bayesian methods.Specifically, Kaandorp et al. [23] stated that an optimized DNN can improve IVIM estimation for tissues with low perfusion and on images with low SNR, and generate less noisy maps with more robust estimation over neighbor voxels.In our work, we have not compared our method with DNN-based approaches; however, we are confident that our results are in line with those reported in [22][23][24].In fact, thanks to the introduction of a spatial dependency through the CAR prior specification, the continuity between adjacent voxels that share the same diffusion and perfusion properties was guaranteed.Furthermore, our estimation accuracy was high even at very low SNR, probably comparable with that obtained by DNNs (a direct comparison is not possible, as different setups and simulations were adopted).
The major limitation of the proposed method is that, although the clustering approach is computationally efficient, the overall approach suffers from high computational load and time because of the Bayesian estimation.Compared to the previously proposed CAR [15], the computational time is not increased by the presence of different weights, but it is anyway high.A further limitation lies in the reduced number of b-values samplings that were considered in the tests.Though, in our results, the performance improvement seemed independent of the b-values, a future work might be focused on evaluating the impact of b-values choice on the performance.Finally, the proposed method was evaluated on a limited number of real cases, looking only at two body districts and at a reduced number of structures; therefore, its application to a larger dataset would allow a better evaluation of its impact from a clinical perspective.
In conclusion, this paper described a novel IVIM fitting approach, based on embedding the k-means clustering approach in the Bayesian CAR method.The proposed method provided more accurate and robust parametric maps with respect to traditional voxelwise approaches, even at the tissue/tissue and tissue/background interfaces.The tool was proven to be effective in both simulated and real cases, and represents a valuable alternative to fit the bi-exponential IVIM model in biomedical applications.

Figure 1 .
Figure 1.Block diagram of the IVIM fitting procedure with clustering.

Figure 2 .
Figure 2. Parametric maps of D (a), f (b) and D * (c), estimated from a simulated image using the three methods (columns) at different SNR values (rows).The true maps are also reported in the top left corner.

Figure 3 .
Figure 3. Error ε maps of D (a), f (b), D * (c), computed from a simulated image using the three methods (columns) at different SNR values (rows).

Figure 5 .
Figure 5. Parametric maps of D (first row), f (second row) and D * (third row), estimated from a selected slice of a HN patient using the three methods (columns).The tumor, the parotid gland and the masticatory muscle are delineated in red, green and blue, respectively.Red arrows indicate improvement in boundary definition, yellow arrows indicate reduced blurred effects.

Figure 6 .
Figure 6.Parametric maps of D (first row), f (second row) and D * (third row), estimated from a selected slice of a pelvic patient using the three methods (columns).The tumor, the prostate and the internal obturator muscle are delineated in red, green and blue, respectively.Red arrows indicate improvement in boundary definition, yellow arrows indicate reduced blurred effects.

Table 1 .
Values of D, f and D * estimated from the real images in each structure of interest by the three methods, reported in terms of mean and standard deviation over the subjects.