# Real-Time Deformable Models of Non-Linear Tissues by Model Reduction Techniques

Siamak Niroomandi, Icíar Alfaro, Elías Cueto, Francisco Chinesta

► To cite this version:

Siamak Niroomandi, Icíar Alfaro, Elías Cueto, Francisco Chinesta. Real-Time Deformable Models of Non-Linear Tissues by Model Reduction Techniques. 22nd International Conference of Theoretical and Applied Mechanics (ICTAM XXII), 2008, Adelaide, Australia. pp.25-29. hal-00290481

HAL Id: hal-00290481

https://hal.science/hal-00290481v1

Submitted on 16 May 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# REAL-TIME SIMULATION OF NON-LINEAR TISSUES BY MODEL REDUCTION TECHNIQUES

Siamak Niroomandi *, Icíar Alfaro*, Elías Cueto* & Francisco Chinesta**

*Aragón Institute of Engineering Research. University of Zaragoza, Spain

** LMSP, UMR 8106 CNRS-ENSAM- ESEM. Paris, France.

*Summary* In this paper we introduce a new technique for the real-time simulation of non-linear tissue behavior based on a model reduction technique known as Proper Orthogonal (POD) or Karhunen-Loève Decompositions. The technique is based upon the construction of a complete model (using Finite Element modelling or other numerical technique, for instance, but possibly from experimental data) and the extraction and storage of the relevant information in order to construct a model with very few degrees of freedom, but that takes into account the highly non-linear response of most living tissues. We present its application to the simulation of palpation a human cornea and study the limitations and future needs of the proposed technique.

## INTRODUCTION

Real-time surgery simulation has attracted the attention of a wide community of researchers, from computer scientists to mechanical engineers, together with computational geometers, surgeons, etc. The utility of such techniques are obvious, and they include, for instance, surgery planning, training of surgeons in image-guided surgery or minimally-invasive surgery, etc.

Such simulators should provide a physically more or less accurate response such that, with the use of haptic devices, a realistic feedback is transmitted to the surgeon in terms of both visual feedback and force feedback. By “accurate response” we mean that an advanced user should not encounter “unphysical” sensations when handling the simulator. We definitely do not pursue an accurate solution in engineering terms. Following [2], “... the model may be physically correct if it looks right”.

For that to be possible, it is commonly accepted that a minimum bandwidth of 20-60 Hz for visual feedback and 300-1000 Hz for haptic display is necessary. In this paper we focus our attention in the second requirement for the deformable model. One of the main limitations of existing real-time simulation algorithms is that they not take into account the anisotropic and highly non-linear response of virtually all soft tissues.

The technique here presented is based upon existing data on the behavior of the simulated tissues. These data can be obtained after numerical simulations made off-line and stored in memory. But they can be also obtained from physical experiments, for instance. For the work here presented we have chosen the first option, and FE models of the organs being simulated will be considered as “exact” to compare with. From these data we extract the relevant information about the (non-linear) behavior of the tissues, with the help of Karhunen-Loève [3] decompositions and employ it to construct a very fast Galerkin method with very few degrees of freedom.

In order to show the performance of the method, we have chosen to simulate the behavior of the human cornea, although the technique is equally applicable to any other soft tissue. The cornea presents a highly non-linear response, with anisotropic and heterogeneous behavior due to its internal collagen fiber reinforcement. As an accurate enough model we have implemented that employed in [1]. The interested reader is referred to that paper and references therein for further details on the mechanical response of the cornea.

## MODEL REDUCTION TECHNIQUES

In this technique we assume that the evolution of a certain field $T(x, t)$ is known. In practical applications (assume that we have performed off-line some numerical simulations, for instance), this field is expressed in a discrete form which is known at the nodes of a spatial mesh and for some times $t^m$. Thus, we consider that $T(x_i, t^m) = T^m(x_i) \equiv T_i^m$ ($t^m = m \times \Delta t$) are known. We can also write $T^m$ for the vector containing the nodal degrees of freedom at time $t^m$. The main idea of the Karhunen-Loève (KL) decomposition is to obtain the most typical or characteristic structure $\phi(x)$ among these $T^m(x), \forall m$. This is equivalent to obtain a function that maximizes $\alpha$:

$$ \alpha = \frac{\sum_{m=1}^{m=M} \left[ \sum_{i=1}^{i=N} \phi(x_i) T^m(x_i) \right]^2}{\sum_{i=1}^{i=N} (\phi(x_i))^2} \quad (1) $$

where $N$ represents the number of nodes of the complete model and $M$ the number of computed time steps. The maximization leads to:

$$ \sum_{i=1}^{i=N} \left\{ \sum_{j=1}^{j=N} \left[ \sum_{m=1}^{m=M} T^m(x_i) T^m(x_j) \phi(x_j) \right] \tilde{\phi}(x_i) \right\} = \alpha \sum_{i=1}^{i=N} \tilde{\phi}(x_i) \phi(x_i); \forall \tilde{\phi} \quad (2) $$

Defining the vector $\phi$ such that its $i$-th component is $\phi(x_i)$, Eq. (2) takes the following matrix form: $\tilde{\phi}^T c \phi = \alpha \tilde{\phi}^T \phi; \forall \tilde{\phi} \Rightarrow c \phi = \alpha \phi$ where the two-point correlation matrix is given by $c_{ij} = \sum_{m=1}^{m=M} T^m(x_i) T^m(x_j) \Leftrightarrow$
$c = \sum_{m=1}^{m=M} T^m(T^{m})^T$, which is symmetric and positive definite. If we define the matrix $Q$ containing the discrete field history: $Q = [T^1T^2...T^M]$ then it is easy to verify that the matrix $c$ results $c = QQ^T$.

If some direct simulations have been carried out, we can determine $T_i^m, \forall i \in [1, \dots, N]$ and $\forall m \in [1, \dots, M]$, and from these solutions the $n$ eigenvectors related to the $n$-highest eigenvalues that are expected to contain the most important information about the problem solution. For this purpose we solve the eigenvalue problem defined before, retaining all the eigenvalues $\phi_k$ belonging to the interval defined by the highest eigenvalue and that value divided by a large enough value ($10^8$ in our simulations). In practice $n$ is much lower than $N$, and this constitutes the main advantage of the technique. Thus, we can try to use these $n$ eigenfunctions $\phi_k$ for approximating the solution of a problem slightly different to the one that has served to define $T_i^m$. For this purpose we need to define the matrix $B = [\phi_1 \dots \phi_n]$.

Now, if we consider the linear system of equations coming from the discretization of a generic problem, in the form: $GT^m = H^{m-1}$, where the superscript refers to the time step, then, assuming that the unknown vector contains the nodal degrees of freedom, it can be expressed as: $T^m = \sum_{i=1}^{n} \zeta_i^m \phi_i = B\zeta^m$ from which results $GT^m = H^{m-1} \Rightarrow GB\zeta^m = H^{m-1}$, and by multiplying both terms by $B^T$ we obtain $B^TGB\zeta^m = B^TH^{m-1}$, which proves that the final system of equations is of low order, i.e. the dimension of $B^TGB$ is $n \times n$, with $n \ll N$.

## RESULTS. PALPATION OF THE HUMAN CORNEA

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Six first eigenmodes of the problem employed as global basis for the reduced model simulation.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Vertical displacement field for a first position of the load. Complete model (left) vs. reduced model (right). eigenmodes employed in this case was only six, see Fig. 1.</figcaption></figure>

The first test for the proposed technique consists of simulating the palpation of the cornea with a surgical instrument. In order to validate the results, a load was applied to the complete FE model in the central region of the model. The obtained result was compared to the one obtained by employing the model reduction techniques presented before, for the load applied at the same location.

Once the complete model is solved, the most important eigenmodes are extracted from the computed displacements field, together with the initial tangent stiffness matrix. The number of

The displacement field obtained for the complete model is compared to that of the reduced model. We chose different positions of the load and compared the results. For a first location of the load, the obtained vertical displacement is shown in Fig. 2. An error of the order of 20% is obtained, that in our experience is a typical bound, and is on the order of magnitude of other real-time simulation algorithms, that employ linear constitutive equations.

## CONCLUSIONS

In this paper a novel strategy is presented for real-time interactive simulation of non-linear anisotropic tissues. The reduced models are constructed by employing a set of “high quality” global basis functions (as opposed to general-purpose, locally supported FE shape functions) in a Galerkin framework. These functions are constructed after some direct simulations of the organs performed by standard FE or BE techniques, for instance. These simulations (their tangent stiffness matrices) are made off-line and stored in memory prior to beginning with the real-time simulation.

Results obtained showed good accordance with complete model results, and ran at frequencies of around 400-500 Hz, enough for real-time requirements, even for this very rude code prototypes.

## References

[1] V. Alastrué, B. Calvo, E. Peña, and M. Doblaré. Biomechanical modeling of refractive corneal surgery. *Journal of Biomechanical Engineering-Tractions of the ASME*, 128:150–160, 2006.

[2] M. Bro-Nielsen and S. Cotin. Real-time volumetric deformable models for surgery simulation using finite elements and condensation. *Computer Graphics Forum*, 15(3):57–66, 1996.

[3] K. Karhunen. Über lineare Methoden in der wahrscheinlichkeitsrechnung. *Ann. Acad. Sci. Fennicae, ser. Al. Math. Phys.*, 37, 1946.