## 3 Methods
### Calculation of Bending Energy Based on Thin-Plate Spline (TPS) Method in the Spatial Domain  
The LMs indicated as yellow points in Figure 2 were placed by an experienced radiologist ([PERSON]) with three years of experience in MR imaging diagnosis of abdomens. Although this manual operation may be replaced by automatic methods, such as Hough transformation used in our preliminary research [26], we planned to use these gold standards in an ideal condition to accurately evaluate the performance of our methods. Figure 2 shows the coordinates of intersecting points of grids  
Figure 1: magnetic resonance (MR) tagging images of a healthy liver (F0) at the maximal inspiratory (frame1, left) phase and the maximal expiratory (Frame 9, right) phase, using 16 mm coronal (**a**); 16 mm sagittal (**b**); 20 mm coronal (**c**); and 20 mm sagittal (**d**) grids. For comparison; (**e,f**) are the images of a cirrhotic liver (F4) using 20 mm coronal and 16 mm sagittal grids, respectively. Note that the grid is deformed over 1 s of forced exhalation, and the deformation of the healthy liver is greater than that of the cirrhotic liver, which reflects the firmness of the parenchyma. determined on (a) the first frame (before warping) and the corresponding points on (b) the last frame (after warping). The reason for setting an LM at the time of inhalation versus the time of expiration is that the difference in the strength of deformation between the normal liver and the abnormal liver conspicuously varies when the BE value is calculated when the volume of deformation of the liver is the largest. These coordinate values are processed to quantify the degree of distortion or deformation of the liver parenchyma. We first applied a TPS-based method referring to a physical analogy involving the bending of a thin sheet of metal, which has been well recognized and used extensively in engineering, to calculate BE values within the liver region. The TPS-based method is widely used for image transformation [25], which requires the setting of the LMs on the image before transformation and the corresponding LM transformation image. In this study, we first set a point \((x_{i},y_{i})\) as an LM on the left of the intersection point of the tag at the first frame (Figure 2a) on the MR tagging image \(I(x,y)\). Then, its corresponding LM is tracked as \((x_{i},y_{i})\) at the ninth frame (Figure 2b). Here, \(1\leq i\leq n\), \(n\) is the number of LMs; \(I(x\prime,y\prime)\) is an image after transformation with the TPS-based method by the mapping function \(f(x,y)=[f_{x}(x,y),f_{y}(x,y)]\), which maps each point \((x_{i},y_{i})\) to its homolog \((x\prime_{i},y\prime_{i})\). To calculate the bending energy \(I_{f}\), the TPS fits the mapping function \(f(x,y)\) between corresponding point-sets by minimizing the following energy function:  
\[E=\iint\limits_{R^{2}}\Bigg{(}\Bigg{(}\frac{\partial^{2}f}{\partial x^{2}} \Bigg{)}^{2}+2\bigg{(}\frac{\partial^{2}f}{\partial x\partial y}\bigg{)}^{2} +\bigg{(}\frac{\partial^{2}f}{\partial y^{2}}\bigg{)}^{2}\Bigg{)}dxdy. \tag{1}\]  
There are many solutions for the function \(E(x,y)\), and the \(I_{f}\) value is selected as the minimal value among these solutions, which satisfies the condition of \(dE=0\). The mapping function  
\[f(x,y)=a_{1}+a_{x}x+a_{y}y+\sum_{i=1}^{n}w_{i}U(|(x,y)-(x_{i},y_{i})|), \tag{2}\]  
minimizes the energy function in Equation (3), where \(U(r)\) can be defined as \(U(r)=|r|\), \(r=|(x,y)-(x_{i},y_{i})|\) is the distance between \((x_{i},y_{i})\) and \((x\prime_{i},y\prime_{i})\), and \(\mathbf{a}\) and \(\mathbf{w}\) are the coefficients between \(I(x,y)\) and \(I(x\prime,y\prime)\), calculated as follows:  
\[\left[ \begin{array}{c}\mathbf{w}\\ \mathbf{a}\end{array} \right]=\left[ \begin{array}{c}\mathbf{v}\\ 0\end{array} \right]\Bigg{[} \begin{array}{cc}\mathbf{K}&\mathbf{P}\\ \mathbf{P}^{T}&0\end{array} \Bigg{]}^{-1}. \tag{3}\]  
Figure 2: Landmarks (LMs) marked in yellow point were manually set by a radiologist.