A comparative study between different intensity standardization methods has found that more sophisticated methods that consider all images simultaneously either directly or at a pixel level improve performance [29]. Another study suggested that it is more appropriate to first perform intensity inhomogeneity correction and then to perform intensity standardization [30]. In WB-MRI, intensity standardization and intensity homogenization are intertwined. Thus, there has been some work on intensity standardization between different parts of an image and between different WB-MRI scans. A simple method identifies peak and trough features in the histograms of different coronal sections or stations. It then fits a continuum for these histogram features along the sections and stations [31]. However, different organs have different histograms, and organs can also contain pathologies that violate this continuity assumption. In a more comprehensive method, the joint intensity histogram of pairs of acquired images is non-rigidly registered to the joint histogram of reference images of the same contrasts to offer a vector field mapping the intensities to their corrected values [32, 33]. Other methods are correcting based on overlapping regions of single-image blocks [34] or of their extension for multispectral image blocks [35]. This method is extended even further to start with a bias correction for the image of each coil separately using N4 [36, 37]. It then continues to register the joint histograms of consecutive stations using their overlap as well as to register to a median station histogram to perform intensity standardization for the whole image [37]. Beyond smooth intensity nonuniformities, the method in this work needs to and does represent the abrupt changes in the intensity inhomogeneity field corresponding to the junctions between consecutive coils. The formation and preservation of discontinuities of the inhomogeneities under smoothing has been achieved with inhomogeneous diffusion also known as anisotropic diffusion. There has been a significant amount of work on these types of methods that preserve discontinuities across regions or equivalently smooth only within regions [38, 39, 40, 41, 42]. An early, but yet extensive, method for piecewise image smoothing and simultaneous noise estimation is based on the Minimum Description Length (MDL) principle formulation [41, 42]. The method allows the separation of an image into piecewise smooth regions [41, 42]. In effect, these methods are similar to the result of the solution of appropriate Euler-Lagrange partial differential equations [41, 42]. Along the same lines, methods for solving geometric partial differential equations over the image domain have also been proposed [39, 40]. Other early methods can separate the image only into piecewise constant components [38]. Image piecewise constancy has also been enforced with the \(L_{1}\) norm constraint [43, 44, 45]. Some of the methods for anisotropic diffusion have been applied directly to biomedical MRI images for denoising [46, 47]. Also, image piecewise smoothness has been combined with smoothness of the nonuniformity, e.g., retinex [11, 12]. However, direct anisotropic diffusion of an image only assumes high-frequency spatial noise and incorrectly interprets shades and the discontinuities they can create as being due to image structure or anatomy [42, 48]. Hence, in biomedical MRI, anisotropic diffusion methods interpret regions with different shading and abrupt transitions between them as the result of tissue differences. The proposed method addresses this limitation by applying the anisotropic diffusion to the shading nonuniformity artifact instead of applying it directly to the image. It thus obtains a piecewise smooth nonuniformity field. The proposed postacquisition method performs a joint restoration of two images of the same anatomic region with different contrasts. The method uses a non-parametric Bayesian coring formulation for intensity restoration.