domain: spi info

text:
Alberto Badias
email: abadias@unizar.es

Iciar Alfaro

David Gonzalez
email: gonzal@unizar.es

Francisco Chinesta
email: francisco.chinesta@ensam.eu

Elias Cueto
email: ecueto@unizar.es

MORPH-DSLAM: Model Order Reduction for Physics-Based Deformable SLAM

published or not. The documents may come

INTRODUCTION

H UMAN visual sense is able to extract and process an enormous quantity of information, simply by observing a few images. A similar behavior is present in standard cameras, where an image can provide a high quantity of information. A strong development in computer vision has been carried out along last decades to extract relevant information from images. In this work, we want to focus in 3D reconstruction methods since our goal is to provide valuable information about the physical changes that happen in the surrounding space of an agent (a navigation robot or a user, employing augmented reality tools for the latter).

Visual augmented reality is based in the addition of graphical information to interact with reality, so it requires great efforts to fix virtual objects within the real world. This means that we need to know the map of the real world to locate the virtual objects. This process is known as Structurefrom-Motion (SfM), and it consists in the estimation of rigid 3D structures from planar images using different camera view points.

Using the same tools, but with a different goal, we may want a robot to understand the surrounding spatial environment, with the ability of understanding the scene where it is moving. Sometimes not only geometric measurements from a sensor are needed, in some cases we should provide the agent with the ability of semantic perception and understanding of the physical changes that happen. We live in a dynamic world, where changes happen at any time, so robots need to have this comprehension ability to really understand the surrounding space.

Some authors suggest the use of the term Spatial AI [START_REF] Davison | Futuremapping: The computational structure of spatial ai systems[END_REF].

It consists of a 3-layer hierarchical system for the spatial understanding of the environment [START_REF]Novel hardware for spatial ai[END_REF], namely Robust Localization, Dense Mapping and Semantic Understanding. In our work we suggest to add a fourth layer to understand the physical changes that occur around the agent (see Fig. 1).

The physical understanding may be gathered within the semantic understanding, but we consider that both are very relevant features and must be considered separately.

Simulated Reality

Fig. 1: Levels in the Spatial AI approach. Image inspired in the CVPR presentation of A. Davison [START_REF]Novel hardware for spatial ai[END_REF].

The first level of the Spatial AI approach is related to the ability of the agent to locate itself in space and navigate through it in a robust manner. The great development in the last decades has brought new algorithms of Simultaneous Localization and Mapping (SLAM) performing localization tasks very efficiently and robustly [START_REF] Mur-Artal | Orb-slam: a versatile and accurate monocular slam system[END_REF]. Level 2 refers to the dense mapping that reconstructs the environment with a high level of detail [START_REF] Newcombe | DTAM: Dense tracking and mapping in real-time[END_REF]. Level 3 requires the system to identify and recognize the objects around it and to be able to establish semantic relationships [START_REF] Chen | Deeplab: Semantic image segmentation with deep lutional nets, atrous convolution, and fully connected crfs[END_REF]. And finally, at level 4, we place the term defined as Simulated Reality (SR) [START_REF] Stoica | A berkeley view of systems challenges for ai[END_REF], where the system needs to know the physics of the surrounding environment in order to have a deeper understanding of the world. SR is based on continuously simulating the physical phenomena of the environment with which the agent interacts. Learning this kind of concepts is a long and complex process, since we must teach a machine how to interpret the physics.

In the field of computer graphics, recent investigations have focused in the development of simplified physics (see, to name but a few [START_REF] Li | Learning particle dynamics for manipulating rigid bodies, deformable objects, and fluids[END_REF], [START_REF] Kennedy | Autonomous precision pouring from unknown containers[END_REF]) in order to overcome obvious computational limitations. However, these simplified approaches provide results with no known accuracy bounds. Dealing with deformable solids, there are some physics engines to approximate the dynamics, but they hardly comply with basic conservation laws. We consider there is no better approach than really solving the (in)elasticity equations, so we suggest to solve the high-fidelity equations. However, standard solvers are not able to work in real time, and even less on portable devices. Since we use visual cameras, the video frequencies are usually fixed in 30 or 60 frames per second, so this is why we suggest the use of Model Order Reduction (MOR) techniques to work in this range of frequencies. MOR methods are based on the compression of the data to work in a reduced space, but at the same time they must recover the original high-fidelity data in a fast way, fulfilling visualization frequencies, with a minimal loss in accuracy. Moreover, this loss is limited by well established error bounds for the vast majority of MOR techniques. The power of these techniques allows us to simulate realistic behaviors that describe and adapt perfectly to the changes that occur in the environment.

In this paper we propose the use of reduced-order numerical methods to estimate the behavior of real deformable solids. Certainly it is usually thought that solving the equations of solid mechanics can be a heavy process, not suitable for video frequencies. But in recent years some tools have been developed to reduce the complexity of these mechanical models. They allow to approximate the true physics of solids at real-time rates.

Including this introduction section, we have divided this paper into 8 parts. Section 2 includes other relevant works in this field from a computer vision perspective. Section 3 explains our method in a few words. Section 4 shows the mechanical formulation of the problem. Section 5 shows the dimensionality reduction problem and the formulation of some MOR methods. Section 6 explains the inverse problem assimilating the data with a deformable implementation of ORB-SLAM2. Section 7 covers the experimentation phase. And finally, section 8 contains the conclusions of the work.

RELATED WORK

The standard approach to compute rigid SfM is based in bundle adjustment (BA). This technique tries to minimize the reprojection error between observed objects and projected objects (observed in more than one image), allowing the estimation of the 3D points and camera locations [START_REF] Triggs | Bundle adjustment-a modern synthesis[END_REF]. SfM method is normally a post-processing technique giving, in general words, more accurate results as much information is extracted from images. But some applications, such as automatic navigation, need an on-line mapping in real time. This means that the agent needs to compute SfM method in real time to be able to navigate in an unknown world, receiving the name of Simultaneous Localization And Mapping (SLAM), involving tasks of estimating the 3D environment structure at the same time the camera is located in the map.

Many works have been developed in this area, allowing to use different sensors (monocular [START_REF] Davison | Monoslam: Real-time single camera slam[END_REF], stereo [START_REF] Paz | Large-scale 6-dof slam with stereo-in-hand[END_REF], RGB-D [START_REF] Gutierrez-Gomez | Dense rgb-d visual odometry using inverse depth[END_REF] or laser technology [START_REF] Cole | Using laser range data for 3d slam in outdoor environments[END_REF], among others). We can make a division in visual SLAM between direct [START_REF] Engel | Lsd-slam: Large-scale direct monocular slam[END_REF] and featurebased [START_REF] Mur-Artal | Orb-slam: a versatile and accurate monocular slam system[END_REF]. Direct approaches use directly the image pixels minimizing the photometric error, while feature-based systems extract features from the images and use them as keypoints to apply BA techniques. Another classification can be made between filtered [START_REF] Smith | Estimating Uncertain Spatial Relationships in Robotics[END_REF] or keyframe-based [START_REF] Klein | Parallel tracking and mapping for small ar workspaces[END_REF] techniques, resulting the first idea in a filtering technique based in the propagation of probability functions and the second idea in a sparse implementation requiring graph optimizations. The work carried out for the last decades has allowed these techniques to move from only research to the development stage, appearing some companies that offer consumer tools based in SLAM techniques.

Another classification, with special interest for us, can be made regarding the position of the objects in the scene. Standard implementations of SfM (and SLAM) assume a rigid world, but as we said, sometimes this assumption cannot be made. The term non-rigid structure from motion (NRSfM) appeared to take into account this kind of movements, and it assumes the camera can move inside the world but also scene objects can move (or deform), as we can see in the right side of Fig. 2. This is an ill-posed problem as we are not able to build the 3D objects from 2D points triangulation. In fact, a different configuration (non-rigid object deformation) can be observed in each image, making the direct estimation of the geometries difficult. During the last years, many works regarding NRSfM have been developed trying to solve this problem. Some of them used factorization approaches to obtain low-rank representations of the objects from image streams [START_REF] Bregler | Recovering nonrigid 3d shape from image streams[END_REF], where the 3D shape in each frame is a linear combination of a set of basis shapes (previously applied to rigid shapes by [START_REF] Tomasi | Shape and motion from image streams under orthography: a factorization method[END_REF]). Also the non-rigid movements have been computed as a union of subspaces to model complex motions with clustering local subspaces [START_REF] Zhu | Complex nonrigid motion 3d reconstruction by union of subspaces[END_REF]. Some other works computed the trajectory space by a linear combination of basis trajectories instead of basis shapes using generic bases [START_REF] Akhter | Nonrigid structure from motion in trajectory space[END_REF]. Nonlinear dimensionality reduction methods have also been used to model the 3D shape as a non-linear combination of basis shapes, using kernel functions [START_REF] Gotardo | Kernel non-rigid structure from motion[END_REF], allowing a reduced number of bases if the problem is non-linear and the kernel is well defined. Also bayesian implementations have been used [START_REF] Agudo | A scalable, efficient, and accurate solution to non-rigid structure from motion[END_REF].

Priors were introduced to narrow down the search for an optimal solution [START_REF] Torresani | Nonrigid structurefrom-motion: Estimating shape and motion with hierarchical priors[END_REF], and some other works allow a prior less estimation of the structure from motion assuming compressible 3d objects as a block sparse dictionary learning problem [START_REF] Kong | Prior-less compressible structure from motion[END_REF], and also applying dense variational reconstructions requiring GPU acceleration [START_REF] Garg | Dense variational reconstruction of non-rigid surfaces from monocular video[END_REF].

Other works reconstruct the shape of a deformable object using templates previously created [START_REF] Bartoli | Shape-from-template[END_REF], obtaining accurate results, but requiring 2D parameterizations as they work only with the external surface, also with curved geometries [START_REF] Gallardo | Shape-fromtemplate with curves[END_REF]. It means some energy constraints need to be imposed to assure boundary continuity and isometric deformations. Other works estimate the surface normals to help in the deformation tracking [START_REF] Lamarca | Defslam: Tracking and mapping of deforming scenes from monocular sequences[END_REF]. And also simple physics-based implementations have tried to solve the motion of the external surface using physical priors and Kalman filter formulations [START_REF] Agudo | Sequential non-rigid structure from motion using physical priors[END_REF].

There is also another approach that involves the use of the classical theory of finite element methods to estimate and track deformable objects with some advantages like robustness and real forces estimation in large deformations, applied to 1D objects [START_REF] Ilic | Non-linear beam model for tracking large deformations[END_REF] but also to more complex surfaces [START_REF] Metaxas | Constrained deformable superquadrics and nonrigid motion tracking[END_REF], [START_REF] Mcinerney | A finite element model for 3d shape reconstruction and nonrigid motion tracking[END_REF].

To the best of our knowledge, there is no alternative approach, however, in which the physics of the system under consideration is solved rigorously. This is due, undoubtedly, to the high computational cost of solving a high-fidelity model by employing state-of-the-art techniques such as finite elements or finite differences. This motivates our choice of reduced-order models, whose accuracy as a function of the complexity of the model is nowadays well-known and can be controlled by the analyst.

OVERVIEW OF THE PROPOSED METHOD

To overcome the just mentioned difficulties, we propose to track deformable solids from video sequences by solving the (possibly non-linear) elasticity problem. The goal is to create an agent able to understand the physical deformations of the objects to finally show some information to a user (e.g. stress or strain fields, or deformations of internal and hidden areas). We make use of reduced-order and parametric models, with which we solve the inverse problem to estimate the kinematic (displacements and strain fields) and dynamic (forces, stress) states, along with the rest of parameter values. Equivalently, the problem can be viewed as an example of data assimilation, in which we obtain parameter values by measuring displacements in the object external surface. Our approach is based on a two-stage problem. First, in an off-line stage we precompute a reducedorder, multiparametric approximation of the displacement field of the solids. This can be viewed as the estimation of a response surface. Second, an on-line stage is carried out to optimize the parameter values that better approximate the measured displacements in the video sequence under real time constraints.

We do not impose any artificial spatio-temporal restriction or prior, nor ad hoc energy laws, since we are solving the actual physics in a finite element framework, guaranteeing energy conservation-dissipative behaviors do not impose any additional conceptual difficulty to this method-and elastic (or hyperelastic) constitutive laws. Of course, initial and boundary conditions need to be applied in the computation step, but they can also be considered as parameters in the parameterized problem. The proposed method can be applied both to articulated and continuous solids, beating any other existing method in terms of frequency rates (since it only requires an evaluation of the parametric solution, rather than a true simulation), accuracy and robustness (noisy observations are minimized without the need of Bayesian implementations). Our method has no dependency on the acquisition system used (monocular, stereo or RGB-D systems), but we think, however, that one of the biggest advantages appears with the use of monocular sensors, since the introduction of parametric models helps in the minimization process, converging on an admissible solution.

Strictly speaking, we can say that our procedure is not a structure-from-motion method, since we are not building a 3D structure from the motion data. The proposed method is related to the Shape from Template family of methods [START_REF] Bartoli | Shape-from-template[END_REF], in which we register an off-line (reduced-order) parametric model to the camera observations.

PROBLEM FORMULATION

Very often the problem of NRSfM has been formulated respect to the motion of the visible portion of the deforming solid (the domain for which some information is captured by the camera). However, it is well-known that the motion of the visible part of the solid is heavily influenced by the whole geometry, including internal details, as described by the Navier equations of solid mechanics. Therefore, in this paper we formulate the problem as a complete, threedimensional description of the solid [START_REF] Bonet | Nonlinear Continuum Mechanics for Finite Element Analysis[END_REF], thus avoiding other NRSfM formulations focused only in the visible part of the solid whose motion was modeled by shell finite elements [START_REF] Agudo | Sequential non-rigid structure from motion using physical priors[END_REF], [START_REF] Bartoli | Shape-from-template[END_REF]. We are interested in some problems like, for example, augmented reality for mechanical design or surgery guidance, some problems that allow us to assume that a full three-dimensional model is available off-line.

Kinematic description of deforming solids

Neglecting inertial effects and assuming, for simplicity, that there is no change in temperature, we consider a continuous solid Ω, where we identify a point Q, with coordinates X = (X, Y, Z) respect to World (W in Fig. 3). These coordinates are defined in the undeformed configuration ϕ 0 : Ω → R 3 . The object undergoes a deformation so that point Q moves to the new position Q with coordinates X = (X , Y , Z ), respect to W, in the deformed configuration ϕ t : Ω → R 3 . The function that explains the movement of the object is φ Obj : Ω × I → R 3 , where I ∈ R is the time interval where the movement is produced. Therefore, making use of a Lagrangian description of the amount of movement (material, always following the same infinitesimal particle), we can express the deformed configuration of point Q as

X = φ Obj (X, t), Π t φ Obj v u Ψ Cam X Y Z W φ pix Π 0 Q Q T 0 T t ϕ 0 ϕ t ϕ 0,pix ϕ t,pix
q pix q pix T Fig. 3: Kinematic modeling of the proposed formulation.

where t is the instant of time elapsed between both states. The deformation gradient tensor F (a two-point tensor) is the fundamental quantity that explains the deformation produced between two neighboring particles between the initial and the deformed state (assuming continuity in the mapping function). It is defined, with respect to the initial configuration, as

F = ∂X ∂X = ∂φ Obj (X, t) ∂X = ∇φ Obj .
It is convenient to define a strain measurement independent of the type of movement (rigid-solid translation or rotation, relative deformations or the combination of all of them). For independence with respect to to rotations, we use the right Cauchy-Green deformation tensor, defined as

C = F T F .
And finally, using the Green-Lagrange strain tensor (material frame of reference) we can measure the deformations independently of rigid-body motions

E = 1 2 (C -I) = 1 2 (∇U ) T + ∇U + (∇U ) T ∇U , (1)
with U (X) = X -X the material displacement.

Camera-centered kinematic description

Let us assume that we obtain an image of the boundary of an opaque solid at the undeformed configuration, ∂Ω 0 . We consider a standard perspective projection camera [START_REF] Hartley | Multiple view geometry in computer vision[END_REF] (after intrinsic calibration and lens distortion correction) applying a transformation Π 0 : ∂Ω 0 → T , where T ∈ R 2 represents the image space. This camera operator Π 0 maps a point Q from World coordinates to q pix in pixel coordinates x pix = (u, v) making use of the camera intrinsic and extrinsic parameters. The same process occurs with projection Π t : ∂Ω t → T , where points in the deformed configuration ϕ t are projected to the image ϕ t,pix ⊂ R 2 .

Solving the unknown kinematics

Using a monocular implementation, the unknown mapping functions in our problem for any time instant t ∈ I are

•

The displacements of the current configuration in a material description (i.e., with respect to the reference configuration) φ Obj (X, t).

•

The camera projection Π t , but only the extrinsic component T t , since it is the part that can change with the camera movement (and also the redundant function Ψ Cam : R 3 → R 3 , that maps the movement of the camera with respect to its initial position).

The estimation of the mapping Π t for each frame capturing rigid scenes is known, as we said, as the Structurefrom-Motion (SfM) problem and it involves the estimation of the extrinsic parameters T t (and the three-dimensional position of the points of the object, Q t ), by minimizing a reprojection error in L 2 -norm, commonly. If we also add the estimation of φ Obj , where objects can deform at any time instant, the problem becomes more difficult (ill-posed). It is possible that a different deformed configuration appears at every instant, with different camera positions. In this paper we assume that the camera captures both static and deformed points. Rigid or static points are related to the general scene (static objects, blue points in Fig. 4) and are used to estimate the pose of the camera, and non-static points are used to determine the deformations of the nonrigid objects in the scene (red points in Fig. 4). We assume all deformable objects have fixed boundary conditions, so at least one point of the boundary is fixed in the scene. However, we could also assume that the objects, in addition to be deformed, move or rotate around the scene, in which case it would be necessary to estimate the relative position of the camera with respect to the object solving a similar but double problem (PnP problem [START_REF] Lepetit | Epnp: An accurate o (n) solution to the pnp problem[END_REF] and displacements).

The position and orientation (pose) of the camera T t is estimated by using fixed points (image features in our case) along with consensus techniques such as RANSAC [START_REF] Fischler | Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography[END_REF] to identify the subset of fixed points in the scene. To estimate the camera pose, the reprojection error d(x pix , xpix ) is minimized, with d the euclidean distance between the Π i ϕ i Fig. 4: In general, we consider applications where an important part of the scene is static (blue points), while some regions (red points) are non-rigid or deformable. pixel coordinates of a point observed in an image and the pixel coordinates of the point projected on the image (from the 3D point x and the camera projection Π t , parameters that are also optimized). We use ORB-SLAM2 [START_REF] Mur-Artal | An open-source slam system for monocular, stereo, and rgb-d cameras[END_REF] for the estimation of the camera pose and the static mapping.

The next task is to estimate the object deformations from its projection in the images. We use the red points in Fig. 4, or in other words, the outliers. RANSAC techniques applied to SLAM usually neglect those points that do not remain fixed in the scene. But our true interest lies in the information associated with the deformation of the objects. It is located precisely there, in the outliers. Therefore, these points provide the most valuable information about the displacements suffered by the objects, which we have to identify after carrying out a registration of the objects.

Equilibrium in the deformable solid

In previous sections we included a brief summary of the kinematic analysis of the deformable solid, but to fully resolve the mechanical problem it is necessary to take into account the stress analysis to guarantee equilibrium. If we split the solid by an arbitrary plane of unit normal n, a force per unit area appears in any point of the newly created surface,

t(n) = σn; σ = 3 i,j=1
σ ij e i ⊗ e j where t is the stress vector, n is the vector normal to the plane and σ is the Cauchy stress tensor.

The stress magnitude that is conjugate to the gradient of deformation tensor F (their product gives an energy measure) is the so-called first Piola-Kirchhoff stress tensor P , defined as

P = JσF -T ,
where J = det(F ) is the Jacobian of the transformation φ. The stress measure conjugate to tensor E is the second Piola-Kirchhoff stress tensor S defined as S = JF -1 σF -T , so we can write equivalently

P = F S.
In order to guarantee equilibrium in the solid, we must satisfy the following equilibrium equation, in the absence of inertia terms (assuming a quasi-static process):

∇ • P + B = 0 in Ω t , (2)
where B are the applied body forces. These equations need to be supplemented with appropriate boundary conditions,

σ n = t on Γ t , u = ū on Γ u
where t are the traction forces acting along the boundary region Γ t in terms of the unit normal to the boundary, n, and finally ū are the prescribed displacements defined in Γ u . The weak form of Eq. ( 2) is created to relax continuity and derivability restrictions on the functions that approximate the final solution, and is obtained after multiplying by an arbitrary function (admissible variation) and integrating by parts, Ω S : δE dΩ = Γt t δu dΓ.

(3)

Constitutive equations

The last ingredient to describe the whole mechanical process is to define its constitutive nature. In this work we consider hyperelastic material laws instead of linear elasticity, as linear materials may introduce artificial increments in volume when large deformations occur. Hyperelastic materials are those whose stress-strain relationship derives from a strain energy density function W . For these materials, the second Piola-Kirchhoff tensor S can be obtained from the strain energy density function W as

S = ∂W (E) ∂E . (4)
Adding Eq. ( 4) into Eq. ( 3), the problem to solve by using finite element methods [START_REF] Fish | A first course in finite elements[END_REF] is

Ω δE : C : E dΩ = Γt t δu dΓ, ( 5
)
where C is the fourth-order constitutive tensor. Section 7 shows a practical example of hyperelastic materials.

MODEL ORDER REDUCTION (MOR) METHODS

One of the most important aspects of our work is that we are assuming that the movement φ Obj is bounded and can be parameterized and projected onto a low-dimensional manifold. This means that the objects we simulate do not have the infinite number of degrees of freedom typical of continuum mechanics. Due to the laws of continuum mechanics, strong correlations exist between the displacement suffered by each material point. Actually, we can express the solution as a function of a much smaller number of degrees of freedom, what we call dimensions in our reduced space.

Dimensionality Reduction

The computational complexity of the continuum mechanics problem prevents its solution under real-time constraints. Note the non-linear character imposed by the strain measures introduced in Eq. ( 1), and also the common non-linear terms of hyperelastic laws. To guarantee that we comply with video frequencies, a reduction of the complexity of the problem is mandatory. In addition, the parameter dependency (material properties, geometry parameters or boundary conditions, among others) increases the complexity of the problem. Let us assume that the governing equation depends on a vector of parameters µ ∈ P, where P is the set of all possible values of the parameters, and a compact subset of R nparam . The manifold or solution set is M, where all solutions U (µ) remain,

M = ϕ(P) = {U (µ) ∈ V : µ ∈ P ⊂ R nparam },
where ϕ is again the solution map and V is a suitable Hilbert space. The solution map ϕ is defined as

ϕ : P → V , µ → U (µ).
It is not possible to work with the original solution set due to the finite memory of computers, so we have to make a simplification by choosing a suitable discretization technique to obtain the high-fidelity solution set M h , also known as the discrete manifold, usually obtained by finite element techniques. We assume that M h is so close to M that very small differences appear at the discretization points. The discretized solution of the equation is U h (µ) and it belongs to a finite-dimensional subspace V h of dimension N h ,

M h = ϕ h (P) = {U h (µ) ∈ V h : µ ∈ P} ⊂ V h .
After applying the approximation to U (µ), we obtain also the discrete solution map

ϕ h : P → V h , µ → U h (µ)
There are many MOR methods to reduce the dimensionality of the problem with the goal of obtaining a reduced basis with an optimal number of components. The reduced solution thus belongs to a low-dimensional subspace V r ⊂ V h of dimension N r N h where the construction of the basis of V r is generated from a set of r functions. The particular form of constructing these functions can be separated between a posteriori and a priori methods. The classical methods are a posteriori, where some data is coming from an unkown source and statistical tools like Principal Component Analysis [START_REF] Karhunen | Über lineare Methoden in der Wahrscheinlichkeitsrechnung[END_REF], [START_REF] Loève | Probability theory[END_REF] are applied to project the data in a reduced and more efficient space. The Proper Orthogonal Decomposition (POD) method [START_REF] Berkooz | The proper orthogonal decomposition in the analysis of turbulent flows[END_REF] is a common example. But, if a multiparametric problem has to be solved, the curse of dimensionality appears, so a priori methods can be applied to solve the mechanical problem directly in the reduced space, without evaluating the whole set of solutions in the high-dimensionality space.

To better understand the procedure, we detail next the type of model order reduction methods commonly applied to continuum mechanics and how they can be of great help in our work.

MOR in Computational Mechanics

We are looking for a suitable approximation of a space, time and parameter-dependent field U (X, µ, t). Under the term model order reduction we encompass a family of methods that look for an approximation of the type

U (X, t, µ) ≈ nMOR i=1 F i (X) • G i (t) • H i (µ), (6)
where F i (X), G i (t) and H i (µ) represent the so-called modes of the approximation, i.e., n MOR vector-valued functions approximated in a finite element sense that best approximate the unknown field U . The symbol "•" represents the Hadamard or Schur (component-wise) product of vectors.

The simplest way to determine an optimal approximation of the type given by Eq. ( 6), sometimes referred to as affine or separate, is by applying Principal Component Analysis of a set of snapshots of the evolution of the system, called Proper Orthogonal Decomposition in computational mechanics [START_REF] Benner | A survey of projection-based model reduction methods for parametric dynamical systems[END_REF]. Once the modes have been identified, by injecting the approximation given by Eq. ( 6) into the weak form of the differential equation governing the problem, the reduced system can be solved.

Other methods like Reduced Basis (RB) employ directly a set of snapshots of the system, whose parameter value is chosen under well-defined error criterion [START_REF] Quarteroni | Reduced basis methods for partial differential equations: an introduction[END_REF].

On the other hand, Proper Generalized decompositions (PGD) [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF] [START_REF] Cueto | Proper generalized decompositions: an introduction to computer implementation with Matlab[END_REF] do not employ any snapshot of the system. This approach presents several advantages. First of all, it is not necessary to solve the problem for all of the parameter combinations in the high-fidelity space. The problem is solved in the separated space, Eq. ( 6), and it can be computed off-line and stored advantageously in memory in the form of a collection of vectors. These vectors include only the nodal values, at mesh vertices, of the functions F i , G i and H i . At any other point of the model, the value of the field U is obtained by finite element interpolation. Secondly, this separated representation presents important advantages for the computation of inverse problems. The sensitivities (parameter influence in the solution) can be computed straightfoward thanks to this precise form of the approximation, very important to boost the optimization step. A detailed formulation is given in Section 6.

Proper Generalized Decompositions (PGD)

As we said, PGD method does not need any observation of the system to construct the low-rank approximation to the unknown field. The method is based in two stages: an off-line step where the computation of the low-rank approximation given by Eq. ( 6) is performed; and an online phase where the parametric solution is evaluated under real-time constraints. To better explain the basics of PGD, we consider an introductory example. As a toy model problem we consider linear elasticity: the case of a cantilever beam in which the position of the applied load is considered as a parameter of the system.

How PGD works: linear elasticity

Let us assume a vertical and constant force F , applied over the upper surface of a cantilever beam. The unknown field of the problem is the displacement U (X, s), which depends on the particular material position X and on the position of the load, s (see Fig. 5). The load can be applied at any point of the upper part s ∈ Γ ⊂ ∂Ω, with Γ the portion of the boundary Γ t with non-vanishing Neumann conditions. The strong form of the governing equation (linear elasticity) under infinitesimal strain theory [START_REF] Fish | A first course in finite elements[END_REF] is

U (X, s) =    ∇ • σ + B = 0 in Ω, σ n = t on Γ t , U = Ū on Γ U , ( 7
)
where σ is the stress tensor, B is the body force per unit volume (assumed vanishing in this example, for the sake of simplicity), t(s) are the traction forces applied to the solid and Ū are the prescribed displacements (clamped support at point 1). The constitutive equations for elastic materials, also known as Hooke's law, represents the material behavior and relates stress and strain as σ = C:ε, with C the fourth-order constitutive tensor and ε the (small) strain tensor. Assuming a hyperelastic framework (of which linear elasticity is a particular example), this constitutive tensor can be obtained by differentiating the strain energy density function W , Eq. ( 4), that in this case takes the form

W (ε) = 1 2 Eε 2 ,
where parameter E is the Young's modulus. Finally, the strain-displacement (kinematic) equations are

ε = 1 2 [∇U + (∇U ) T ].
Keeping in mind that U depends on X = (X, Y ) and s, the weak form of the governing equation under the PGD formalism can be described as

Γ Ω ∇ s U * : σ dΩ d Γ = Γ Γt U * t dΓ d Γ ( 8
)
where U * ∈ H 1 0 is an arbitrary test function in the appropriate Sobolev space of functions vanishing along the Dirichlet boundary of the solid and ∇ s = 1 2 [∇ + (∇) T ] is the symmetric gradient operator. The solution U (X, s) is assumed to be expressed in separate form,

U (X, s) ≈ nMOR i=1 F i (X) • G i (s), ( 9
)
where n MOR is the number of sums (rank of the tensor decomposition of the solution), F i (X) are the space modes in separate variables and G i (s) are the modes referred to the position of the load. An alternating direction scheme is used to obtain the n MOR functions F i (X) and G i (s) so that, at iteration p we search the p + 1 functions R(X) and S(s) that enrich the solution

U p+1 (X, s) = U p (X, s) + R(X) • S(s).
x Domain (mm) Spatial Amplitude Therefore, the test function is

U * (X, s) = R * (X) • S(s) + R(X) • S * (s).
The load F has been designed as a unitary force acting along the vertical axis Y ,

t = F δ(X -s)e Y , X ∈ Γ t ,
where δ is the Dirac delta function. Function t needs also to be expressed in separate form to comply with the PGD formalism, i.e.,

t Y ≈ m j=1 h j (X) k j (s); t X = 0,
where m is the number of sums to approximate function t Y and h j (X), k j (s) are the functions in space and load position, respectively. To solve the non-linear product of functions R(X) and S(s), the algorithm uses an alternating direction scheme based on the fixed point method at every enrichment iteration. For a detailed description of the problem and a Matlab implementation, the interested reader can consult [START_REF] Cueto | Proper generalized decompositions: an introduction to computer implementation with Matlab[END_REF].

The first modes of F (X) (only its vertical component, since it is a vector-valued displacement field) and G(s) are shown in Fig. 6.(a) and 6.(b), respectively. Finally, the sum of the products of pairs of modes builds the approximate solution, as it is shown in Fig. 6

.(c).

It is important to highlight the fact that the accuracy of the proposed reduced-order method can be controlled by the user. The number of terms in Eq. ( 9), n MOR controls this level of accuracy, with the help of a suitable error estimator, as the ones proposed in [START_REF] Alfaro | An error estimator for real-time simulators based on model order reduction[END_REF] [47] [START_REF] Ammar | An error estimator for separated representations of highly multidimensional models[END_REF], for instance.

A DEFORMABLE IMPLEMENTATION OF ORB-SLAM2

We propose in this section a deformable implementation of the ORB-SLAM2 system [START_REF] Mur-Artal | An open-source slam system for monocular, stereo, and rgb-d cameras[END_REF]. ORB-SLAM2 is a sparse and feature-based SLAM system with outstanding results of accuracy [START_REF] Mur-Artal | Orb-slam: a versatile and accurate monocular slam system[END_REF]. It is used to build a rigid scene and locate the camera in any frame. In our work, after a few seconds of video capturing the static scene, we get the 3D point cloud and apply an automatic and on-line registration against a CAD model of the object, based in RANSAC methods and the ICP (Iterative Closest Point) algorithm [START_REF] Besl | Method for registration of 3-d shapes[END_REF].

The registration step allows us to know which points are deformable and, subsequently, to estimate the displacements of the 3D SLAM points (U SLAM (X, µ)) using the precomputed data from the CAD object

U SLAM (X, µ) = f U Neighs,CAD (X, µ) . ( 10
)
To do that, a mapping function between our undeformed 3D model and the 2D image features need to be established. It means we can afterwards deform the real object, extract image features in any frame and estimate the parameter values µ with respect to the precomputed set of solutions. The goal is to obtain the displacements U SLAM (X, µ) to add to the static set of deformable 3D points initially estimated using the SLAM technique. This solution must maintain the separate multiparametric structure to obtain a specific solution for each possible value of our set of parameters.

Therefore, the displacements of the SLAM points are expressed as a function of the closest registered CAD points. The procedure can be seen in the general overview schema (Fig. 7). Fig. 7: Deformable object tracking and augmentation procedure based in ORB features matching.

(u,v) Q(X,Y,Z) dQ(X,Y,Z,µ) Image projection Q(X,Y,Z) + dQ(X,Y,Z,µ) Deformable Reprojection Minimization Parameter J (µ) = nmeas j=1 ρ u meas pix (x j ) -Πt U MOR (X j , µ) 2 2D Πt(Q + dQ(µ))
In the visual scanned point cloud there are 3D points that remain static (used to fix camera position and orientation) and 3D points that are assigned to the deformable part of the domain in the current observation. Let us focus on the points belonging to the deformable part, Q, with coordinates (X, Y, Z). The position in its deformed configuration is Q , where U = (X -X) = (X , Y , Z ) -(X, Y, Z), but this displacement U is actually estimated from the deformations of the CAD object, which depend on the multiparametric solution projected in separate variables. The point in its deformed configuration (Q ) is observed from the image as q with pixel coordinates (u , v ), which must match the projection q e = Π t (Q + U (µ)) for the optimal set of parameters µ. The reprojection error in the image is therefore defined as d(q , q e ) and the optimal parameter values of our parametric solution that minimize this error for the whole set of observed points are obtained from the next functional

J (µ) = nmeas j=1 ρ u meas pix (x j ) -Π t U MOR (X j , µ) 2 (11)
where n meas are the number of measurements, u meas pix the pixel coordinates measured in the image, U MOR (X j , µ) the 3D displacement predicted by the reduced-order parametric solution, Π t the projection to the image plane at any instant t (frame), and ρ the Huber robust function. It is important to highlight that U MOR (X j , µ) represents the whole solution of our hyperelastic problem in question by storing the position, deformations and stresses of each point, with parametric dependence and represented in the reduced space. However, to simplify the process, U in the functional defined in Eq. ( 11) only represents the deformed position of each point. Using Eq. ( 11) we can estimate the value of the µ parameters that minimize the reprojection error

μ = argmin µ J (µ).
This minimization can be advantageously accomplished by employing the Levenberg-Marquardt algorithm [START_REF] Levenberg | A method for the solution of certain non-linear problems in least squares[END_REF][51], which takes the form

[J T J + λ diag(J T J )] δ = J T U meas pix (x j ) -Π t U MOR (X j , µ) ,
where δ is the perturbation parameter that moves the parameters in the direction of steepest descent and λ an algorithmic parameter that adaptively varies the behavior between the Gradient Descent Method and the Gauss-Newton Method. Finally, J is the Jacobian matrix, defined as

J = ∂U (x, µ) ∂µ . (12)
Here lies precisely the key point that makes the separate representation of the solution so powerful in the parameter optimization process. Recall the general form of PGD, Eq. ( 6), where we can notice that due to its separate form, a separate differentiation can also be performed. We can precompute the derivative function of each parameter, Eq.( 12), offline, since the finite element approximation computes the derivative as a multiplication of matrices. These derivative functions are usually known as sensitivities of the solution [START_REF] Chinesta | Pgd-based computational vademecum for efficient design, optimization and control[END_REF] with respect to each parameter µ k . We can interpret them as the local variation of the function U (X, µ) with respect to the changes in parameter µ k

J k = N i=1 F 1 i (X) • . . . • ∂F k i (µ k ) µ k • . . . • F nparam i (µ nparam ).
It allows us to avoid the exploration of all of the parametric space in search of the derivatives, and directly apply the differentiation on the separated variable vectors, involving a great reduction in the computational task. In sum, we need the partial components of the derivative terms of the

Q I W T C I Π C W T O q O W C ∂u I (x,µ) ∂µ k ∂U O (X,µ) ∂µ k
Fig. 8: Transformations and projection of the precomputed sensitivities from the original system of coordinates to the image reference.

multiparametric solution projected on the image, since we are minimizing the reprojection error (2D), but we precomputed the sensitivities in the original 3D space, so we have to project these sensitivities onto the image plane. A graphic example to understand this process appears in Fig. 8, where

J k = ∂U O (X, µ) ∂µ k (13)
refers to the partial derivative term of the multiparametric 3D solution of the CAD object with respect to the parameter µ k , evaluated at point Q and with the origin of coordinates defined on that object (O), and

j k = ∂u I (x, µ) ∂µ k ( 14
)
represents the same partial derivative projected on the image plane and evaluated on the point q . To obtain the value of this sensitivity in the image plane, j k , we use the chain rule

j k = ∂u I (x, µ) ∂U C (X, µ) ∂U C (X, µ) ∂U W (X, µ) ∂U W (X, µ) ∂U O (X, µ) ∂U O (X, µ) ∂µ k ,
where u I (x, µ) is the position of point q with respect to the image plane, U C (X, µ) with respect to the camera coordinates, U W (X, µ) with respect to the global coordinate axis (World) and U O (X, µ) with respect to the coordinate axis defined in the object. The solution derivative terms are known in the object reference, Eq. ( 13), and their consecutive transformations in the Euclidean space R 3 are carried out by the matrices W T O and W T C , which are assumed to be known since they are estimated by the static part of ORB-SLAM2 and the registration step. However, the projection on the image plane I Π C does not allow such a simple treatment, since the camera conical projection used depends on the location of the 3D point with respect to the image plane. This means that the projection of the derivative terms R 3 → R 2 must take into account the position of the point where the gradients are evaluated. This projection is known as the image Jacobian, and is expressed as

I J Π C = ∂u I (x, µ) ∂U C (X, µ) = 1 Z C f x 0 -X C Z 2 C f x 0 1 Z C f y -Y C Z 2 C f y ,
where (X C , Y C , Z C ) are the 3D coordinates of point Q with respect to the camera axes. Finally, the projection of the multiparametric derivative terms with respect to parameter µ k on the image plane is expressed as

∂u I (x, µ) ∂µ k = I J Π C • W T C -1 • W T O • ∂U O (X, µ) ∂µ k , (15)
where we remember that the sensitivities defined with respect to the object axes are interpolated according to the neighbors in the separate variables precomputed solution.

Regarding to the technical aspects, we use the feature search method developed in the ORB-SLAM2 code to track the deformations. Time requirements for the minimization algorithm are much smaller than the computation time of the SLAM tasks, so a simple parallelization of the whole code in a few threads allows us to work in a frequency of 30 fps using the CPU. We stop the Mapping tasks and work in Localization mode in ORB-SLAM2 once deformations are appearing.

EXPERIMENTS

We provide two examples to support our method. Both have been tested with different video sequences working in real time.

Rubber boot seal

This example shows the behavior of a rubber boot seal, the part that protects the gearshift of a vehicle. The object is deformed forced by the rotation of the gearshift in one direction, from -52 • to 52 • , but the simulation has been only applied within the parametric range θ ∈ [0 • , 52 • ], taking advantage of revolution symmetry. The material behavior applied to the rubber is a nonlinear law, called Neo-Hookean model, with the strain energy density

W = µ 2 (I 1 -3) -µ ln(J) + λ 2 (ln(J)) 2 ,
where I 1 is the first invariant of the right Cauchy-Green deformation tensor

I 1 = λ 2 1 + λ 2 2 + λ 2 3 ,
and where λ i are the eigenvalues. The mesh of the object has 25795 nodes and 18403 linear hexahedral elements. Applying a sparse implementation of PGD [START_REF] Ibá Ñez | A multidimensional data-driven sparse identification technique: the sparse proper generalized decomposition[END_REF], we obtain the parametric solution in the projected space, depending in space X = (X, Y, Z) and one parameter, θ, where the parameter optimization is a very fast process as we precompute the sensitivities. Due to the projection process of the data, a certain error in the approximation is made, as it can be seen in Fig. 9. With a total number of 29 modes, we get a cartesian mean error of 0.18 mm (the normalized error is 0.0041). For this solution, the compression factor is 92.58%, computed as

C(%) = 1 - M P M O • 100,
where M P is the memory storage used by the projected solution and M O is the memory cost necessary to store the original data. In other words, using a reduced space to express our solution we only need a 7.42% of the storage memory, allowing real time evaluations. Since we assume that the geometry of the deformable model is known, but not its external appearance (texture) or its location in the scene, we apply a registration process between the static point cloud and the CAD object. The object to be deformed has geometry of revolution, so we must actually optimize the value of two parameters µ = (θ, α Z ), where parameter θ measures the deformation and α Z the angle of orientation of the axis of rotation, a geometric parameter that we cannot estimate in the registration process (see Fig. 10). This is why, when the deformation occurs, we have to estimate the value of θ in each frame. To favor the minimization process, the value of α Z is minimized throughout the entire sequence, so after a few frames observing the deformed object, the system itself is able to estimate with great degree of certainty the value of α Z , which involves also a better accuracy in the estimation of parameter θ.

Fig. 11 shows some frames extracted from a video sequence with the original object, and the virtual object with the estimated displacements. The central and right frames show the stress field in colors, result of solving the whole mechanical problem. Finally, in Figs. 12 and 13 the 3D error is estimated. These error values have been calculated as the euclidean distances between the surface of the deformed virtual object and the object scanned with an RGB-D camera (Intel RealSense D415). It is important to note that using an accurate model of a complex geometry like the one in this example gives better results than the scanning of an RGB-D system, mainly due to the precision of the camera itself. For example, it can be seen in image 1 of Fig. 12 that the RGB-D camera is not able to accurately estimate the folds of the real object. In any case, the median error using a monocular camera is 1.33 mm (Fig. 13) and the mean error is 1.49 mm. These very accurate values are the result of using a detailed model of the geometry, an adequate description of the material law, and a measurement system as accurate as the non-deformable implementation of ORB-SLAM2.

Paper sheet buckling

The second example simulates the buckling phenomenon of a paper sheet. From a numerical simulation point of view, buckling is a complex problem that is usually solved using linear approximations to obtain the critical loads (buckling modes) that produce the beginning of buckling. For this, an eigenvalue analysis is used, looking for the loads for which the model stiffness matrix becomes singular [START_REF] Novoselac | Linear and nonlinear buckling and post buckling analysis of a bar with the influence of imperfections[END_REF]. However, if a complete analysis of the buckling problem is desired, in order to know in detail the values of the stresses and deformations of the object in question, it is necessary to apply more complex methods such as the Riks method, introducing small imperfections on the geometry [START_REF] Novoselac | Linear and nonlinear buckling and post buckling analysis of a bar with the influence of imperfections[END_REF].

Again, the hypothesis of the use of a reduced model in this example is based on the fact that the number of degrees of freedom of the problem is not high. In other words, although it may appear that the object is arbitrarily deformed, there is a great restriction due to the geometry of the paper, so the number of degrees of freedom of the problem is very small, actually living in a reduced dimension manifold.

The paper sheet material has been considered anisotropic (reinforced with oriented fibers) to emulate the real behavior of the object [START_REF] Schulgasser | On the in-plane elastic constants of paper[END_REF], [START_REF] Yokoyama | Evaluation of in-plane orthotropic elastic constants of paper and paperboard[END_REF]. The model mesh has been discretized in 1333 nodes, 31 in the vertical direction and 43 the horizontal direction. To model the buckling problem we use one single parameter (µ B ) discretized at 801 buckling positions. Applying again the sparse implementation of PGD [START_REF] Ibá Ñez | A multidimensional data-driven sparse identification technique: the sparse proper generalized decomposition[END_REF] it is possible to obtain a projected solution u = (X, Y, µ B ), reaching a compression factor of 98.36%.

The sequence was recorded with an RGB-D camera (StereoLabs Zed Mini) in order to estimate the error made by our monocular method, where Fig. 15 plots the mean error in millimeters between the model and the real object for each frame. Fig. 14 shows some frames extracted from the sequence, where colors represent the Von Mises stress field. Note that this error, that may seem big, is comparable to the best of the state-of-the-art, with one important difference: we are solving the physics of the problem, and therefore are subject to errors in the calibration of material parameters, for instance [START_REF] Lamarca | Defslam: Tracking and mapping of deforming scenes from monocular sequences[END_REF]. In addition, deformations in the considered sequence are much bigger than those in [START_REF] Parashar | Robust isometric nonrigidstructure-from-motion[END_REF], with a maximum vertical displacement value (also corresponding to the maximum error peak) of 16.9 cm. Relative error is thus on the order of 10% at this maximal value.

The fact that both deformations and displacements (as a rigid solid) of the object are allowed makes it also necessary to estimate the deformation parameters and the pose of the object in each frame. This means that the robustness of ORB-SLAM2 in its rigid part cannot be exploited and the errors are slightly higher than those of the previous example, but our method is more general and robust to displacements. The 3D mean error of the whole sequence is 3.98 mm, measured respect to the ground truth data coming from the RGB-D camera and taking into account there is a relevant camera error of 1.2 mm (orange line in Fig. 15). This error is one order of magnitude less than the one reported in state of the art references, such as [START_REF] Lamarca | Defslam: Tracking and mapping of deforming scenes from monocular sequences[END_REF], for instance, that reported errors up to 30-35 mm in the displacement prediction.

CONCLUSIONS

We propose in this paper the application of computational mechanics as a tool for estimating the deformations suffered by the objects perceived by a standard camera. This allows us to augment the video sequence with information regarding the displacements, forces and stresses suffered by the objects. To do this, we rely on the use of models to obtain a parametric solution (using finite element methods) that we use in the on-line phase. This parametric solution is stored in a compressed way to optimize both the cost of storage and its evaluation in real time.

One of the great differences of this work with respect to the bibliography is that we do not have to impose any type of spatio-temporal restriction, nor apply energy conservation laws implemented ad-hoc, as we are solving the real physics of the deformable solids. This also implies that the result is defined for the whole volume of the modeled objects, not only for the external surface, allowing us to capture the real behavior of the objects although partial occlusions appear. It also provides information about how the interior of the objects is changing, even though we are not able to capture them with the camera. It assums a clear advantage for biomedical applications where it is necessary to know, for example, the location of the blood vessels or cancerous bodies inside the organs. Also different material properties can be applied in any point of the solid, giving different results depending on the stiffness, but assuming more expensive previous work.

As it is shown in the experiments section, our method has a great computational efficiency, is CPU-based and quite robust against noisy observations. It is not necessary to precompute sequences to extract trajectories or learning shapes, as we obtain the information directly from the models. But the solution must be obtained for a parametric set of loads, states and deformations enough to be able to work in any sequence.

Finally, this work can use any type of model order reduction method, and applied to any type of camera device, such as monocular, stereo or RGB-D systems.

Fig. 2 :

 2 Fig. 2: a) Rigid SfM allowing a reconstruction of the world from different views. b) Non-rigid SfM implies that both the camera and the scene may take different positions and shapes (time-dependent).

Fig. 5 :

 5 Fig. 5: Cantilever beam problem. A moving load F is parameterized by its position coordinate s.

Fig. 6 :

 6 Fig. 6: Separated solution of the vertical displacement. (a) First five modes depending on space. (b) First five modes depending on the position of the load. (c) Built solution using the sum of the product of modes. Note how the result of the PGD approximation provides with a sort of response surface, without the need of any parameter space sampling strategy. Finally, a slice of the surface with the prescribed load position s 0 builds the solution.

Fig. 9 :θFig. 10 :

 910 Fig.9: Cartesian mean 3D error in displacements due to the projection of the solution in the reduced space respect to the number of modes.

Fig. 11 :

 11 Fig. 11: Some images showing the original and the augmented frames. In the central and right images we can see the stress field that the object is suffering when the deformation appears.

Fig. 12 :Fig. 13 :

 1213 Fig.12: Error of the monocular estimation, measured assuming the ground-truth is coming from a stereo camera. Due to the high complexity of the geometry, the stereo system can not be more accurate than the sum of model and monocular camera in some areas, producing a constant error that can be shown in the undeformed frame 1 .

Fig. 14 :Fig. 15 :

 1415 Fig. 14: Some images showing the original and the augmented frames in the second example. Colors in images represent the real-time stress field that the object is suffering when the deformation appears.

ACKNOWLEDGMENTS

We are grateful to our colleague Prof. J. M. Martínez Montiel, whose help with some technical aspects is much appreciated.

This work has been partially funded by the ESI Group through the ESI Chair at ENSAM Arts et Metiers Institute of Technology, and the project 2019-0060 "Simulated Reality" at the University of Zaragoza. We also thank the support of the Spanish Ministry of Economy and Competitiveness, grant CICYT-DPI2017-85139-C2-1-R and the Regional Government of Aragon and the European Social Fund.

Alberto Badías is currently a postdoctoral researcher at the University of Zaragoza, Spain. He received the B.S. degree in mechanical engineering in 2011, the M.S. degree in industrial engineering (industrial automation and robotics) in 2014, the M.S. degree in biomedical engineering in 2016 and the Ph.D. degree in 2020, all from the University of Zaragoza. He has been working in the area of computer vision and robotics developing 3D reconstructions and new image processing tools, and is working currently in model order reduction methods in applied mechanics and bioengineering.

Iciar Alfaro is an associate professor of continuum mechanics at the university of Zaragoza. Her research has been devoted to the development of innovative numerical methods for complex problems in continuum mechanics. She developed natural element formulations for forming processes, and, more recently, she works in the development of advanced techniques for augmented reality with strong physical meaning.

David Gonz ález is currently a professor at the university of Zaragoza. He obtained his Ph.D. degree in applied mathematics in the same university in 2004, with a thesis on the development of meshless methods for the Lagrangian simulation of free-surface flows. Since then, he has been working on the development of advanced model order reduction techniques, with an eye towards their application in real-time simulation, computational surgery and augmented reality.