## 2 Materials and Methods
### IVIM Fitting Procedure
#### 2.1.2 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 \(\tilde{f}(i,j)\), where \(\tilde{f}\) denotes an intermediate estimation of \(f\), are estimated with the observations at \(b=0\) and \(b\in B\) : \(b\geq\tilde{b}\), with \(\tilde{b}=200\) s/mm\({}^{2}\). A simplified version of the IVIM model is used, which only includes the second exponential decay:  
\[SI(i,j,b)=SI(i,j,0)\big{[}1-\tilde{f}(i,j)\big{]}e^{-b\ D(i,j)}\ \ \ \ \ \forall\ (i,j)\in\mathcal{I},b\in B:b\geq\tilde{b} \tag{2}\]  
In the second step, the coefficients \(D^{*}(i,j)\) and \(f(i,j)\) are estimated using all observations \(b\in 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 \(\sigma_{lik}\) of the conditioned densities that determine it. Then, a prior density for the unknown parameter is assumed. A priori independence between \(\sigma_{lik}\), each \(D(i,j)\) and each \(\tilde{f}(i,j)\) is assumed for the first step, and between \(\sigma_{lik}\), each \(D^{*}(i,j)\) and each \(f(i,j)\) for the second step. Separately for each \(\lambda\) (where \(\lambda\) generically denotes each coefficient \(D\), \(D^{*}\), \(f\) or \(\tilde{f}\)), a CAR specification is included in the prior by considering the intrinsic CAR formulation of [PERSON] et al. [17]. More specifically, the following distribution for a coefficient in a voxel given the values in the other voxels is taken:  
\[\lambda(i,j)|\lambda^{\mathsf{c}}_{i,j}\sim\mathcal{N}\left(\frac{\sum_{( \alpha,\beta)\in\mathcal{I}}w^{\alpha,\beta}_{i,j}\lambda(\alpha,\beta)}{ \sum_{(\alpha,\beta)\in\mathcal{I}}w^{\alpha,\beta}_{i,j}},\frac{ \sigma^{2}_{\lambda}}{\sum_{(\alpha,\beta)\in\mathcal{I}}w^{ \alpha,\beta}_{i,j}}\right) \tag{3}\]  
where \(\lambda^{\mathsf{c}}_{i,j}=\boldsymbol{\lambda}\setminus\{\lambda(i,j)\}\), and \(w^{\alpha,\beta}_{i,j}\) denotes the spatial neighborhood matrix. We assume a dependency (\(w^{\alpha,\beta}_{i,j}\
eq 0\)) only for the voxels \((\alpha,\beta)\) bordering on \((i,j)\). Their values are assigned based on the clustering solution \(\big{\{}\mathcal{C}_{ij^{\prime}},(i,j)\in\mathcal{I}\big{\}}\). In particular, \(w^{\alpha,\beta}_{i,j}=1\) if the bordering voxel \((\alpha,\beta)\) belongs to the same cluster of voxel \((i,j)\) (i.e., \(\mathcal{C}_{\alpha\beta}=\mathcal{C}_{ij}\)),or \(w_{i,j}^{\alpha,\beta}=0.1\) if the bordering voxel \((\alpha,\beta)\) does not belong to the same cluster of voxel \((i,j)\) (i.e., \(\mathcal{C}_{\alpha\beta}\
eq\mathcal{C}_{ij}\)); therefore:  
\[w_{i,j}^{\alpha,\beta}= \begin{cases}1&(\alpha,\beta)=\{(i-1,j-1),(i-1,j),(i-1,j+1 ),\\ &(i,j-1),(i,j+1),(i+1,j-1),\\ &(i+1,j),(i+1,j+1)\}\\ &\text{and }\mathcal{C}_{\alpha\beta}=\mathcal{C}_{ij}\\ 0.1&(\alpha,\beta)=\{(i-1,j-1),(i-1,j),(i-1,j+1),\\ &(i,j-1),(i,j+1),(i+1,j-1),\\ &(i+1,j),(i+1,j+1)\}\\ &\text{and }\mathcal{C}_{\alpha\beta}\
eq\mathcal{C}_{ij}\\ 0&\text{otherwise}\end{cases}  \tag{4}\]  
In addition, to smooth the autoregressive component, the actual prior density of each \(\lambda\) consists of a mixture model with two components: the above conditional autoregressive specification (3) and a truncated Gaussian marginal density \(\mathcal{N}_{(\lambda_{min},\lambda_{max})}(\mu_{\lambda},\sigma_{0\lambda}^{ 2})\), where \(\lambda_{min}\) and \(\lambda_{max}\) denote the minimum and maximum admissible values, respectively.