## 3 Methods
### Calculation of Bending Energy Based on Thin-Plate Spline (TPS) Method in the Spatial Domain
First, the LMs on the intersection points of grids were marked on the first frame, Frame 1 (**a**), of the MR tagging image before deformation. The corresponding LMs were then tracked and placed on the image after non-rigid transformation in the last frame, Frame 9 (**b**). \(\mathbf{P}\), \(\mathbf{V}\), and \(\mathbf{K}\) are minor determinants of \(L\), where  
\[P=\left[ \begin{array}{ccc}1&x_{1}&y_{1}\\ 1&x_{2}&y_{2}\\...&...&...\\ 1&x_{n}&y_{n}\end{array} \right],V=\left[ \begin{array}{ccc}{x_{1}}^{\prime}&{ x_{2}}^{\prime}&...&{x_{n}}^{\prime}\\ {y_{1}}^{\prime}&{y_{2}}^{\prime}&...&{y_{n}}^{\prime}\end{array} \right],K= \left[ \begin{array}{cccc}0&U(r_{12})&...&U(r_{1n})\\ U(r_{21})&0&...&U(r_{2n})\\...&...&...&...\\ U(r_{n1})&U(r_{n2})&...&0\end{array} \right] \tag{4}\]  
and  
\[\mathrm{L}=\left[ \begin{array}{cc}K&P\\ P^{\mathrm{T}}&O\end{array} \right]. \tag{5}\]  
Based on the selected LMs on two images, the BE value of \(I_{f}\) is finally calculated by  
\[I_{f}=V({L_{n}}^{-1}{KL_{n}}^{-1})V^{T}. \tag{6}\]  
Therefore, the bending energy calculated by Equation (6) is regarded as the value of the non-rigid deformation of the liver. Although the TPS method is a very traditional deformation method, the use of the BE value as a classification feature is rare. Calculation of the Difference in the Power Spectral Values of the Tags by FFT in the Frequency Domain  
Manually placing the LMs is time-consuming and impractical for clinical use. To overcome this disadvantage, our method for measuring non-rigid deformation is modified to the frequency domain processing using a 2D FFT method. Because the main patterns in the MR tagging image are grids that periodically appear at a certain principal frequency \(f_{0}\) (or wavelength) corresponding to the distance between the intervals of the grids; any changes in the interval of the grids will make its frequency shift from \(f_{0}\) to \(f_{0}+\Delta f\). A larger change in the shape of the grids results in a wider distribution of the range of \(\Delta f\), which represents different frequency components. Such a deformation, very similar to the optic phenomenon, appears to be lost in focus at a frequency point of \(f_{0}\) on its 2D power spectrum image. This implies another way of quantifying the non-rigid deformation by applying the 2D FFT method to an MR tagging image and then comparing the difference in the power spectral values in a candidate region around \(f_{0}\). FFT is an efficient algorithm for computing the discrete Fourier transform (DFT); the calculation of the DFT on a limited dataset that greatly decreases the computing iterations and time is optimized. 2D DFT on a 256 \(\times\) 256 MR tagging image can be written as  
\[F(u,v)=\sum_{x=0}^{255}\sum_{y=0}^{255}F(x,y)e^{-j2\pi(xu/256+yv/256)} \tag{7}\]  
where \(F(x,y)\) is the gray value at point \((x,y)\) on an MR tagging image in the spatial domain, and \(F(u,v)\) is the complex value at \((u,v)\) corresponding to the transverse and longitudinal frequency components of \(F(x,y)\) in the frequency domain derived from FFT, where \(x,y,u,v\in(0,1,2,3\cdots 256-1)\). \(j\) is the imaginary unit. The matrix of the power spectrum is calculated as  
\[PWS_{F}(u,v)=\left|F(u,v)\right|^{2}=Fr^{2}(u,v)+Fi^{2}(u,v) \tag{8}\]  
where \(PWS_{F}(u,v)\) is the power spectrum of \(F(u,v)\), and \(Fr(u,v)\) and \(Fi(u,v)\)are the real and imaginary parts of \(F(u,v)\), respectively. The low-frequency components (small gray-filled squares) distribute in the four corners of matrix \(F(u,v)\), and the highest frequencies are located in the middle of the matrix as shown in Figure 3a. To visualize the power spectrum into a gray-scale image for further processing, we re-arrange the four quadrants from the power spectrum matrix \(PWS_{F}(u,v)\) into \(PWS_{I}(u,v)\) as shown in Figure 3b. Because of the symmetries of the spectrum, the entire set of quadrant positions can be diagonally replaced following the direction of the arrows indicated in Figure 3a after replacingthe highest frequencies at the corners of the matrix.