Deformable Model Acquisition and Validation by Jochen Lang M.Sc, York University, Toronto, Canada, 1996 B.Eng.(Hons), University of Plymouth, Uk, 1992 Diplom-Ingenieur (FH), Fachhochschule Ulm, Germany, 1992 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy in THE FACULTY OF GRADUATE STUDIES (Department of Computer Science) We accept this thesis as conforming to the required standard The University of British Columbia October 2001 © Jochen Lang, 2001 In presenting t h i s thesis i n p a r t i a l f u l f i l m e n t of the requirements for an advanced degree at the U n i v e r s i t y of B r i t i s h Columbia, I agree that the Library s h a l l make i t f r e e l y a v a i l a b l e for reference and study. I further agree that permission for extensive copying of t h i s thesis for s c h o l a r l y purposes may be granted by the head of my department or by his or her representatives. It i s understood that copying or p u b l i c a t i o n of t h i s thesis for f i n a n c i a l gain s h a l l not be allowed without my written permission. Department of Aeu -Severe e The U n i v e r s i t y of B r i t i s h Columbia Vancouver, Canada Q c . 4 . a Date A b s t r a c t Objects deform in response to contact forces exerted on them. The deformation depends on material properties, the geometry of the object and external forces. This thesis develops a robotic system for automatically acquiring observations of a deforming object and for estimating a model of the deformation from these obser-vations. Models of deformable objects are in wide-spread use in simulation, computer graphics and virtual reality. Deformation, impact and fitting simulation aid man-ufacturing. In computer graphics deformable objects are designed and animated. Medical simulators incorporating physical models of organs and tissue are a signifi-cant emerging virtual reality application. The material properties of deformable models are often assigned based on mechanical (and other) testing of material samples. Material samples do not repre-sent commonly simulated objects well if there is a high variability in the material of the object, e.g., due to mixed material, unknown material, or material with im-perfections. In contrast to material sampling, this thesis develops a method to scan the deformation behavior of a complete object. The scanning is analogous to the scanning of the visual appearance of an object. The scan captures the individual response of a complete object to contact forces. The result of the scan is not only a deformable model but also data which serves to validate the model. This validation step gives a qualitative and sometimes quantitative assessment of the suitability of a model. ii C o n t e n t s Abstract ii Contents iii Acknowledgements vii 1 Introduction 1 1.1 Problem Statement and Overview 5 1.1.1 Hypothesis 5 1.1.2 Approach 5 1.1.3 Computational Deformation Model . . . 6 1.1.4 Overview 7 1.2 Related Work 7 1.2.1 Deformable Models 8 1.2.2 Inverse Problems in Elasticity 12 1.2.3 Theoretical Results 12 1.2.4 Elastic Parameter Estimation 13 1.3 Summary of Related Work 14 2 Elastostatic Object Behavior 15 iii 2.1 Discrete Green's Functions Matrix 16 2.2 Deformation Models in Continuum Mechanics 18 2.2.1 Strain Tensor 18 2.2.2 Stress Tensor 20 2.2.3 Constitutive Equation 21 2.2.4 Navier's Equations 22 2.3 Solution Methods 22 2.3.1 Finite Element Model 23 2.3.2 Boundary Element Model 25 2.4 Other Physics Based Deformable Models 28 2.4.1 Particle Systems 28 2.4.2 Deformable Superquadrics 30 2.5 Summary 31 3 The A C M E Facility 32 3.1 Overview 34 3.1.1 A C M E Devices 35 3.1.2 Control Architecture 37 3.1.3 Programming A C M E 38 3.1.4 Sensor Model 40 3.2 A C M E Deformation Experiment 41 3.3 Geometric Modeling with A C M E 44 3.4 Summary 45 4 Deformation Measurement 47 4.1 Robust Flow Estimation 49 iv 4.1.1 Robust Range-Flow by Additional Imagery 51 4.1.2 Implementation Details 54 4.1.3 Results for Rigid Body Motion 56 4.1.4 Examples of Object Deformation 62 4.2 Estimation of Vertex Displacement 62 4.2.1 Displacement Results 63 4.3 Summary 66 5 Estimation of Discrete Green's Functions 68 5.1 Set-Up of Estimation 69 5.2 Local Behavior 71 5.2.1 Least Squares 72 5.2.2 Measurement Design Matrix 73 5.2.3 Regularization . . 76 5.3 Global Behavior 81 5.3.1 Regularization 82 5.3.2 Outliers 84 5.3.3 Hole Filling .- 86 5.4 Mesh Refinement 87 5.5 Results 89 5.5.1 Plush Toy . . . . 90 5.5.2 Anatomic Soft-Tissue Human Wrist Model 94 6 Inverse Problem Solution 102 6.1 Boundary Element Method for an Inverse Elastic Problem 103 6.1.1 Elastic Parameter Estimation 103 v 6.1.2 Optimization . . 105 6.1.3 Results for Half Ball 106 6.2 Approximate Model 108 6.3 Summary 112 7 Conclusions and Future Work 113 7.1 Conclusions 113 7.2 Future Work 115 Appendix A Linear Triangular Elements 117 Appendix B Geometric Modeling with A C M E 120 Bibliography 125 v i A c k n o w l e d g e m e n t s I would like to thank my thesis supervisors Dr. Robert Woodham and Dr. Dinesh Pai. I am grateful for Bob's thoughtful consideration and questions, as well as his encouragement. I like to thank Dinesh for his enthusiasm and support in exploring the ideas of this thesis. Also, I would like to thank my committee member Dr. Uri Ascher for his valuable comments. I am thankful for many helpful discussions with Doug James on deformable models and haptic rendering. Dr. John Lloyd's diligent work on the University of British Columbia (UBC) Active Measurement facility ( A C M E ) enabled the exper-imental work of this thesis. This thesis also profits from the work of all the other A C M E group (former) members: Rod Barman, Timothy Edmunds, Julian Fong, Stewart Kingdon, Paul Kry, Josh Richmond, and Som Yau. I am also thankful to Don Murray at Point Grey Research (PGR) for help with the P G R Triclops Stereo Vision SDK. Furthermore, I like to acknowledge Dr. Gerhard Roth at the National Research Council, Canada for making his volumetric mesh creation software avail-able to me. Jesse Hoey, Paul K r y and John Lloyd provided helpful comments on a draft of (parts of) this thesis. I am grateful for the financial support of my doctoral studies by the Uni-versity of British Columbia and by the National Sciences and Engineering Research Council of Canada. I would like to acknowledge IRIS Precarn Centres of Excellence and their support of the reality-based modeling project ( A C M E ) at U B C . J O C H E N L A N G The University of British Columbia October 2001 vii C h a p t e r 1 I n t r o d u c t i o n Interactive simulation of physical objects is desirable for numerous applications in the areas of v i r tua l and augmented reality. T h e goal in interactive simulation is to ensure v i r tua l objects resemble their real counterparts as closely as possible but at the same t ime maintain interactivity. T h e consequence of this trade-off is the need to employ an approximate model of the physical behavior of objects. In general, the realism mathematical models of objects can achieve depends on the proper choice of the physical parameters of the model. How should these physical parameters be chosen? T h i s thesis argues that realistic simulated behavior requires parameters which are based on measurements and that these measurements can be acquired efficiently. T h e physical behavior explored is the deformation of solid objects due to contact forces. Th i s thesis shows how to efficiently scan this behavior i n a robotic measurement facility very similar to the scanning of object geometry in a shape scanner. Simulations gain realism if the effect of forces on an object is represented realistically. Forces govern many behaviors of real world objects. In particular, forces deform an object depending on its geometric and material properties. The 1 simulation of realistic v i r tua l deformations necessitates the inclusion of forces, ob-ject geometry and realistic material properties. Realist ic deformation of v i r tua l objects is required in medical t ra ining simulators [31], in haptic interfaces to v i r tua l worlds [113], and in design of robotic tasks involving soft materials [59, 60]. M e d -ical simulators allow students to practice [123], and help to predict and plan some surgical procedures [67]. Med ica l robotics can also profit from user interfaces wi th haptic feedback [112]. Force s imulat ion plays a very important role in the creation and animation of v i r tua l worlds [40]. Forces are a very natural way for a designer to interact w i th v i r tua l objects. Real is t ic force responses are intuit ive because they copy the behavior of our surrounding physical world. A realistic force response is a powerful method to manipulate complex surfaces while meeting other constraints. A n example is the editing of a model of a human face, the skin is modified by the user but forces due to bone, muscle and tissue structure resist the changes and ensure the modification is in accordance wi th physical l imits [126]. Another example is the Virtual Clay system of M c D o n n e l l , Q i n and Wloda rczyk [79] where forces dur ing sculpt ing are simulated in real t ime allowing the creation of v i r tua l sculptures. M a n y computer graphics techniques have evolved over the years to model deformations [40]. These techniques vary from full continuum models for elastic materials, through mass-spring particle systems, to very ad-hoc methods. B u t even the most elaborate deformation model w i l l be unrealistic if the simulated material does not match the real physical object. T h e measurement of these material prop-erties has so far only been done in a "material oriented" way, e.g., in physical and mechanical testing laboratories. A s a,result, many computer graphics methods do not model material properties for deformation simulation of physical objects. Fur-2 thermore, simulation results with a specific computational model are commonly not validated with respect to physical object behavior [40, 31]. The key for future increase in realism is observations. These observations of physical phenomena need to be sufficient to develop a computational model for objects to be simulated. The combination of a measurement approach, a compu-tational model and a robotic measurement system is required. This combination creates a scanner for deformation behavior of objects which is the goal of this the-sis. Scanning deformation and other physical interaction behavior has recently been described by Pai et al. [95]. Scanning the geometry of objects is common practice today. Shape infor-mation for a wide class of everyday objects can be acquired with scanners typically based on optical triangulation. Another physical property beside shape which can be scanned is the reflectance of an object represented by its bidirectional reflectance distribution function (BRDF) [116, 20]. < Given the ubiquitous graphics rendering capabilities of today's computers, reality-based animations of three-dimensional ob-jects are commonplace. Up to very recently, the physical behavior of objects due to contact forces could not be scanned. This is part of the reason why these behaviors are rarely rendered. Scanning physical object behavior due to contact forces may be understood as an inverse problem. In an inverse problem a model is determined based on obser-vations. (Given a parametric model that applies, the process of determining model parameters forms the field of parameter estimation.) Inverse problems related to de-formable models are being researched in the mechanics of solids and structures [66]. In mechanics, computational methods are often more elaborate than in computer graphics since accuracy of simulation results is critical in many applications (e.g., in 3 structural analysis in c iv i l engineering). Surprisingly, the verification and validation of computat ional methods in inverse problems in mechanics are often neglected as well . T h e performance of computat ional methods for inverse problems is commonly evaluated by comparison wi th either an analytic solution or by comparison wi th the same method using finer discretization of the model domain (or boundary) . Direct experimental validation provides a more meaningful characterization of a method. T h e scanning method of this thesis makes experimental val idat ion simple. However, scanning of object deformation does not need to be approached as an inverse problem. Deformation behavior can be directly observed and rendered. T h e interior stress and strain, often the quantities of interest i n mechanics, are not directly necessary to render deformation behavior. The purpose of the deformable model in computer graphics is to allow physically realistic rendering of the defor-mat ion of a real object under force. The model has to encode the response of the object as a whole to a force. W h a t interior structure causes that response is irrele-vant. However, the scanning of object deformation is also the first step in an inverse problem solution. Th i s thesis also contains a .method to estimate physical model parameters based on scanned behavior. T h e model assumption in this case is that the object behaves as a homogeneous linear isotropic elastostatic body. Next , in Section 1.1, a concise problem statement is given which may guide the reader through the following chapters. Section 1.2 reviews work related to this thesis. 4 1.1 P r o b l e m S t a t e m e n t a n d O v e r v i e w 1.1.1 Hypothesis The central hypothesis is: A framework which acquires and validates a computational model of object deformations by displacement and force measurements enhances the utility of the model. This utility is in the rapid availability of the model, and in the knowledge of the applicability and limitations of the model available from the comparison between model behavior and actual object behavior. 1.1.2 Approach This thesis develops a measurement apparatus, a measurement method and an esti-mation procedure to fit a computational deformable model to an object's measured deformation behavior. A result of the fitting process is a model and a residual which characterizes the ability of the model to explain the measurements. The model is able to synthesize (novel) object behavior which allows one to analyze and verify model behavior. The approach taken in this thesis may be characterized as an analysis by synthesis framework [130]. This thesis develops and utilizes the U B C active measurement facility (ACME) for object deformation measurements. A C M E is a robotic facility which automates object measurement tasks. Object deformation is characterized by surface displace-ments and contact forces, and potential topological changes in the case of inelastic deformation. A measurement facility for object deformation has to be capable of applying and measuring contact forces, as well as measuring the displacement at the contact itself and over the free surface of the object. The measurement method of this thesis is a direct method in which the 5 behavior of the object to be modeled is directly measured. • The forces over the surface of the object and the displacement of the free surface of the object are observed. The force and displacement are measured by a robotic probe while it actively deforms the object. The displacement of the free surface is observed visually with a trinocular stereo vision system. The three-dimensional displacement vectors over the surface are obtained with simultaneous (geometric) stereo and optical flow. The measurement of the initial (undeformed) shape of an object is also based on stereo range data. The estimation procedures developed in this thesis are applicable with dif-ferent prior knowledge. A procedure with very weak assumptions about the object behavior is developed, as is a method with very restrictive assumptions about the behavior. Robustness is emphasized during the estimation process. Outliers typical of many visual techniques are handled by measurements which overconstrain the solution and by regularization. 1.1.3 Computational Deformation Model In this thesis the deformable model is a discrete linear elastostatic boundary model. The boundary model is based on a triangular surface tessellation of the object. The model linearly relates the traction on surface elements to the displacement of surface elements. The model coefficients depend on the object shape and its elastic material properties. The model is adequate to render small deformations1 which an object undergoes if a force is exerted on it. In the case of a homogeneous, isotropic and linearly elastic object the above linear boundary model corresponds to a discretized 1Here, small means the small deformation strain tensor for linear-elastic material is employed neglecting second order effects (see Section 2.2.1). Also, the initial tessellation of the objects surface is unchanged during the deformation limiting the approximation of the discrete representation of object geometry. 6 continuum mechanics model. In this case, the coefficients of. the model can be derived with a triangular surface tessellation in a finite element formulation of the object's elastostatics (see Section 2.3). The efficient haptic and visual rendering of this type of boundary model is developed in the thesis work of James [56]. 1.1.4 Overview Work related to this thesis can be found in Section 1.2. The reader unfamiliar with elasticity may find helpful explanations regarding the cited work in Chapter 2. Chapter 2 describes elastostatic object behavior and the deformable model of this thesis in detail. Chapter 3 summarizes the capabilities of the A C M E robotic mea-surement system and gives details of features specific to deformation modeling. The visual deformation measurement approach is introduced in Chapter 4. This Chapter includes results and ground truth comparison for the method. The next two chap-ters describe the estimation method: Chapter 5 explains the procedure for acquiring estimates for the linear boundary model: Chapter 6 describes the inverse problem solution to elastic parameter estimation for objects made of isotropic homogeneous linear elastostatic material. This thesis concludes with Chapter 7. Possibilities for continuation of this work are discussed. 1.2 R e l a t e d W o r k Related work to thesis has been done in a number of fields with varying goals. In computer graphics the efficient representation of deformation behavior of solids has been surveyed by Gibson [40], while Section 1.2.1 reviews deformation models in interactive simulation. The focus of the review is on how parameters for the em-ployed model are chosen. Section 1.2.2 introduces work on related inverse problems 7 in elasticity. Theoretical results have been obtained for the inverse problem of elastic parameter estimation and are summarized in Section 1.2.3. Section 1.2.4 discusses numerical and experimental methods in elastic parameter estimation. Finally, Sec-tion 1.3 gives a brief summary of related work. 1.2.1 Deformable Models This section focuses on deformation models for physically realistic interactive sim-ulation. At the Forschungszentrum Karlsruhe in Germany, the development of medi-cal training simulators has been ongoing since the early 1990's [67]. The simulation is based on KISMET (Kinematic Simulation, Monitoring and Off-Line Programming Environment for Telerobotics). The Karlsruhe Endoscopic Surgery Trainer employs two types of deformable models: particle systems [125] and linear elastostatic finite element models with and without condensation [18]. The stiffness of living tissue is estimated based on ultrasound measurements and compared to stress-strain curves of animal samples obtained in tension test machines and with a compression test [78]. The stress-strain relation is fit to a third order polynomial a = Ee(l + P\e + P2£2) where a is the stress, e the strain, E is a scaling parameter (Young's modulus) and P\,P2 are non-linear shape factors. This relationship is simulated with a particle system. Szekely et al. report [122] a Laparoscopic Surgery Simulator (LASSO) at the ETH Zurich, Switzerland, which includes deformable models of human organs. The organs are modeled as hyper-elastic material. The simulation is based on finite elements computed on a specifically designed parallel computer with 64 DEC Alpha processors. The material parameters are estimated based on measurements with an 8 aspiration tube; a novel instrument which applies a tensile load to a small area of the tissue by suction and observes the displacement visually [129]. The material law employed in simulation and estimation is based on a modification of the Mooney-Rivlin strain energy function W = /i(J\ — 3) + a(J2 — 3) + J 3 — l ) 2 where J\, J2 and J3 are the reduced invariants of the right Cauchy-Green deformation tensor and H, a and K are the material parameters. The estimation is based on a finite element •model of the area surrounding the probe without the inclusion of bulk modulus K. The non-linear least squares error is minimized based on displacement and tension force measured over time. However, due to instabilities, a is arbitrarily set to zero and only \x is estimated. Within the Sharp project at INRIA, Rhone-Alpes, C. Laugier and his group modeled human tissue for medical simulators [27, 29]. The force interaction between an echo-graphic sensor head and the human thigh is modeled. The human thigh is represented with a particle system including linear and non-linear springs. The parameters of the springs are derived from measurements with a robotic probe [27]. The model is optimized to match the measured local force response within 5%. The authors noted the stiff behavior of the human thigh and the difficulty of modeling .this behavior with a particle system. A haptic interface to the deformable model has been presented in [26]. A human liver is also modeled however no measurements of the behavior are reported [29]. The Epidaure project at INRIA, Sophia-Antipolis, is also modeling deformable objects for medical simulation. Bro-Nielsen and Cotin [18] introduced finite element modeling with condensation into interactive simulation. The model applies to lin-ear elastostatic simulation. Their model representation is strongly related to the presentation employed in this thesis as detailed in Section 2.3.1. More recently, 9 Cotin, Delingette and Ayache'[25] utilized a tensor-mass model solved with finite elements for soft tissue in a dynamic simulation. This model is combined with an elastostatic model in a hybrid approach. The topology of the model can be changed in real-time by removing complete tetrahedra. Very recently Picinbono, Delingette and Ayache [99] have adapted a dynamic model with a non-linear large deformation strain tensor. The model is solved with finite elements and mass-lumping with an external incompressibility constraint. A comparison of the behavior of a physical brain sample in a compression test to a simulation is reported by Cotin, Delingette and Ayache [24]. Keeve et al. [64] employ deformation models of facial tissue to predict the outcome of cranofacial surgery. Simulation results are compared with the actual outcome of the operation., Two types of simulation are reported: one based on particle systems with piecewise linear springs; the other based on finite element models with a non-linear strain-stress relationship. The finite element solution is not interactive. Basdogan et al. [7] model the deformation behavior of soft tissue with free form surfaces. Later, De and Srinivasan [28] substituted a physical boundary model. They model the finger-pad as a membrane in plane strain filled with incompressible fluid. The model is linear but varies over time. Simulation results are compared to experimental measurements. Deformable models in VESTA (Virtual Environments for Surgical Training and Augmentation) at the University of California2 are computed with finite ele-ments with mass lumping [131] as in [136]. The non-linear material law is based on Mooney-Rivlin and a Neo-Hookean model. The computation is sped up with a 2The project is a joint cooperation between the University of California, San Francisco, Berkeley and Santa Barbara. 10 geometrically adaptive scheme. The group is also investigating soft tissue behavior during a surgical procedure (needle driving) and under tension [19]. Koch et al. [65] report a craniofacial surgery simulation. This system models tissue behavior with a surface finite element method attached to a volumetric particle system. The model is based on individual bone, fat and soft tissue distribution deduced from a CT scan of the head. The simulation result is compared qualitatively to the result of the surgery. The interactive surgical cutting simulator of Biesler, Maiwald and Gross [11] employs particle systems only. Haptic manipulation of soft objects and tissue based on the finite element method with the small deformation linear strain tensor is reported in [50] and [88]. Debunne at al. [30] describe a space and time adaptive scheme for dynamic simu-lation of soft objects. The solution method employed is finite elements with mass lumping as in the off-line simulation of O'Brien and Hodgins [91]. Other interactive deformation models can be found in [104, 36, 15, 23]. In the context of haptic simulation the capabilities of the human sensory system can provide some guidance to how realistic simulation needs to be. In [12] the capability of the human haptic system is summarized. The stiffness necessary for an object to be perceived rigid is at least 25N/mm. The JND (just noticable difference) of compliance is about 5 — 15% if the displacement is fixed during the trial and the object surface is rigid. The maximum force that can be exerted controllably by a human finger pad is about 100iV with 0.04iV or 1% resolution, whichever is higher. 11 1.2.2 Inverse Problems in Elas t ic i ty Inverse problems aim to discover the cause of a physical phenomenon from its man-ifestation. Problems in elasticity are usually labeled inverse when the unknowns of the problem are not displacements or strains, and not traction or forces [66]. The identification of cracks, cavities or flaws based on some deformation of the object (e.g., [54, 9, 90, 121]) form one set of inverse problems. The estimation of elastic co-efficients given undeformed and deformed object state as well as the corresponding deformation forces (e.g., [117, 39, 107, 45, 97]) is a problem of parameter estima-tion. The boundary conditions may be unknown on part of the boundary. The estimation of these boundary conditions from, e.g., other over-specified boundaries (e.g., [10]) or interior measurements is a problem occurring in non-destructive test-ing. Another set of inverse problems deals with inverse geometry, the calculation of the undeformed object shape given the forces and the deformed shape of the object (e.g., [44, 43]). Applications of inverse elasticity problems can be found, amongst others, in medical imaging [37], non-destructive testing (e.g., [9]), compos-ites (e.g., [46]), civil engineering (e.g., [4]) and earth sciences (e.g., [41]). 1.2.3 Theoret ical Results Related theoretical research in elastic parameter estimation considers the question: do boundary conditions provide enough information to deduce the elastic properties of a material? Various properties of the map of traction and displacement boundary conditions over an inhomogeneous medium have been proven. In an inhomogeneous medium the elastic constants are varying over the medium, i.e., they become func-tions over the domain. Akamatsu, Nakamuru and Steinberg [1] prove that a map of complete boundary conditions determines the Lame functions A, \JL in the case of a 12 two-dimensional, isotropic, linear-elastic, inhomogeneous material (some conditions on the Lame functions and the smoothness of the elastic domain apply). A similar result has later been derived for higher dimensions [85]. Furthermore, Nakamura and Uhlmann [86] also showed that the functions are unique. Additionally, a generic anisotropic linear-elastic tensor in two dimensions can be determined in two dimen-sions [87]. Unfortunately, this result can not be extended to higher dimensions [84]. These results may give some guidance in designing.a measurement procedure and apparatus. They also support the idea that modeling of anisotropic, inhomo-geneous deformable object from boundary measurements is better approached as a direct scanning problem rather than an inverse problem. 1.2.4 Elastic Parameter Estimation Solution strategies for inverse problems of elasticity may be classified based on the need for observations inside a domain of the elastic continuum. Interior observa-tions require a more elaborate experimental setup than observations of the bound-ary alone. Parameter estimation using domain methods is applied to beam struc-tures [114, 4] and for two-dimensional domains [107, 46] and if special measurement techniques (like ultrasound) allow for domain interior measurements [41, 61, 103]. (The solution procedure for all the domain methods in [117, 114, 4, 107, 46, 41, 61, 97] is based on finite elements, except for the finite difference approach of [103].) If the inverse problem involves a three-dimensional continuum and no interior measure-ments are available, the problem is more difficult. The boundary element method has been applied in three-dimensions even if interior measurements are not avail-able [9]. Elastic coefficients can be estimated with at least two different inverse ap-13 proaches (as defined by S. Kubo [66]): a direct method, or an optimization method based on either force error minimization or displacement error minimization [4]. One may isolate the elastic parameters for an homogeneous isotropic material [66, 107] resulting in a direct method. A direct method may also be obtained by re-arranging the boundary element expression [71] (see Section 6). Displacement or strain error minimization are applied in [39, 114, 4] and a sensitivity analysis of each respective method is provided in [39, 115, 5]. In a case-study of a beam structure [5], Banan, Banan and Hjelmstad find that force error minimization leads to very similar results as displacement error minimization. 1.3 S u m m a r y o f R e l a t e d W o r k The trend in interactive simulation of deformable models has been towards phys-ically more accurate models. Finite element solutions to problems with physical constitutive relations of the continuum are often employed. The computational re-quirements limit the size and the complexity of the stress-strain relationship which can successfully be modeled at interactive update rates. Dynamic models with large deformation are especially challenging. The need for an appropriate choice of elas-tic parameters for these increasingly realistic models has been recognized. Elastic parameter estimation is researched in various fields, however no generally accepted solution for three-dimensional continua exists. The scanning of deformation behav-ior developed in this thesis represents an alternative approach to acquire realistic deformation models. 14 C h a p t e r 2 E l a s t o s t a t i c O b j e c t B e h a v i o r Deformation of solids has been studied extensively in physics, civil engineering, me-chanical engineering and biomechanics. The theory of elasticity provides analytical tools to describe the deformation of elastic solids. Linear elasticity theory forms the basis of empirical methods which describe the tension, shear and compression of engineering solids. Section 2.1 states the model of linear static elasticity which this thesis em-ploys. The model is a discrete Green's functions matrix which is based on continuum mechanics. Section 2.2 contains a brief introduction into linear elastostatic models in continuum mechanics which the reader familiar with continuum mechanics may wish to skip. In Section 2.3, solution strategies of the boundary value problem associated with elastic solids are discussed to the extent that they are the basis of this thesis. Section 2.3 details how either the Finite Element Method or the Boundary Element Method may be employed to derive a discrete Green's functions matrix. Again, the reader familiar with these methods may wish to skip these sections. Finally, some common modeling approaches used in computer vision, computer graphics and haptics are summarized in Section 2.4. 15 2.1 D i s c r e t e G r e e n ' s F u n c t i o n s M a t r i x The elasticity of a solid determines its deformation when external forces act on it. In static equilibrium, the net effect of the external forces is zero and the solid's boundary remains stationary (Section 2.2.4 gives a mathematical description). The domain of the problem is a solid 0 with its boundary F. The boundary has two regions of boundary conditions, To and T i , such that r 0 U T i =T and To n F\ = 0. On the boundary To displacements are prescribed and on the boundary T i tractions are prescribed, (see Figure 2.1), i.e., To is fixed while T i is free. A traction p is a force f normalized over the surface area A it is affecting, i.e., p = f/A. (Note, that instead of tractions, forces may be incorporated directly since the undeformed object surface is a constant. Details can be found in Section 2.3.) The Green's functions relate a field of displacement vectors to a field of traction vectors on the boundary of the elastic solid. The boundary in the model of this thesis is discretized into a triangular mesh with k = 0... n vertices (see Figure 2.1). The displacement and traction vector at each vertex are denoted by and p ,^ respectively. The basis functions of the displacement and traction fields are hat functions </> centered at a vertex. The hat function's magnitude is one at the vertex, linearly decays'to zero at surrounding vertices, and is zero elsewhere on the surface (see Appendix A for a mathematical definition). The displacement and traction vectors are combined into block vectors u = [ i v f , . . . , u ^ ] T and p = [pj,..., pT]T, respectively. For a given boundary configu-ration, the block vectors u and p can be rearranged such that all prescribed values are collected in a block vector v. Entries = if the vertex k is on the sur-face To, while = if the vertex k is on the surface I V The complementary 16 Figure 2.1: Discrete Boundary Value Problem associated with an Elastic Solid displacement and traction vectors are entered into a block vector v. The matrix relating the prescribed values v and v is the discrete Green's functions matrix is S (see Equation 2.1). Equation 2.1 is central to this thesis. v = 3v (2.1) In simulation, a discrete Green's functions matrix S and the input vector of prescribed tractions and displacements v allows one to calculate the resulting de-formed shape of the solid as well as the reaction forces. In estimation, the Green's functions matrix S is to be determined from over-specified displacement and trac-tion boundary conditions. The block vectors v and v are measured or known in the estimation of the Green's functions matrix S. This measurement and estimation problem is the topic of the remainder of this thesis. The next sections explore in detail the assumption in the derivation of the Green's functions matrix. In Sec-tion 2.2, the basic definitions of elasticity are reviewed concluding with the Navier's equations. Section 2.3 shows how the Navier's equations can be solved resulting in the discrete Green's functions matrix. 17 2 . 2 D e f o r m a t i o n M o d e l s i n C o n t i n u u m M e c h a n i c s Mathematical descriptions in continuum mechanics are based on some definition of the material stress state and the corresponding strain of the material. In general, a stress is a force density acting at an internal section of material area while strain is the relative change in dimension of a material. Various mathematical definitions of stresses and strains are possible, see Sections 2.2.1 and 2.2.2. The strain of a material is caused by its stress state (and potentially its stress history). The constitutive equations of the material relate stress and strain. Constitutive equations together with an equilibrium condition lead to Navier's equations for an elastic solid. The following brief Sections 2.2.1, 2.2.2 and 2.2.3 give more concise definitions. 2.2.1 S t r a i n T e n s o r The deformation of a solid is the change in the location of an element relative to other elements of the solid: For the definition of a strain tensor, define a displacement vector u from the position x of a point on a body in undeformed state to the same point y on the body in the deformed state. The displacement is u = y — x. Then, we can define a Jacobian from undeformed to deformed state as dx\ dx2 9x3 1 + j = du2 i i du2 du2 dx\ dx2 dx% du3 du3 i , du3 dxi 0x2 9x3 The above Jacobian may be employed to define a strain tensor. Two points on the body P and Q are a squared distance (<5/o)2 = (dx)Tdx apart if the body is undeformed. In the deformed state, they are a distance (SI)2 = (dy)Tdy apart. Since y = x + u it follows that dy = dx + du/dx dx = (I + du/dx)dx. However 18 the Jacobian is defined as J = I + J = / + <9u/<9x and it follows that the change in squared distance between point P and Q is (dl)2 - (5l0)2 = (dyfdy - ( d x ) T d x = (dx)TJTJdx - ( d x ) T c i x = ( d x ) T [ ( J T + J ) + J r J ] d x 3 3 - (dx)T2edx = Y l Y l d x i 2 £ i J d x J (2-2) i = i j = i Equa t ion 2.2 defines the symmetric strain tensor Sij and because of its sym-metry may be wri t ten as a strain vector e = [e\ti £2.2 £3,3 £1,2 £1,3 £ 2 , 3 ] T ' • Equa t ion 2.2 can be expressed as e = D u where D is a differential operator. In engineering applications, the strain is usually assumed to be small , i.e., the Jacobian product JTJ 0 and the strain tensor simplifies em ~ \\JT + J ] . The differential operator D reduces then to three tensile and three shear strains. 2 ^ 0 0 0 24-OX2 0 0 0 0x3 d 9x2 d 9x\ 0 9 9x3 • 0 9 3xi 0 a 9x3 9 9x2 (2.3) T h e Green-Lagrange strain tensor of Equa t ion (2.2) is independent of the path of deformation, which is adequate for linear elasticity. Note however, that if rotational components are significant, JT J Z$> 0, the approximation no longer holds. The approximation is also not suitable for plastic deformation which depend on the path of deformation. Rate-of-deformation and spin tensor are employed in these cases [135]. 19 2.2.2 Stress Tensor a, a a, Figure 2.2: Stress on Infinitesimal Volume Element The stress tensor relates forces aligned along the axes of the coordinate sys-tem to forces on a arbitrary oriented stress plane inside the material. In order to derive the stress tensor, consider the diagram in Figure 2.2. Define a normal sur-face force as f = pSS with the surface traction p acting on the area SS (shaded triangle in Figure 2.2). The unit normal of the stress plane is fi and the surface traction is therefore p = an. The remaining sides of the elementary tetrahedron are aligned with the coordinate axis. These remaining sides have reaction forces acting on them to counterbalance the acting force f. These reaction forces f; are in general not aligned with the normal of their respective sides. The stress on these sides is o\ = fi/SSi. However, simple geometry leads to 552 = nrfs because of the alignment with the respective coordinate axis. Newton's second law for the infinites-imal tetrahedron (i.e., its volume and therefore forces due to its mass tend to zero) becomes f! + f2 + f3 = f. This yields Cauchy's stress formula (Eqn. 2.4) and the stress tensor Oij. Pj = arij = o~ijrij (2.4) 20 T h e stress tensor is symmetric which allows for a vector notation a = °~2,2 03 ,3 <7it2 ° 2 , 3 ] T - T h e Cauchy stress tensor does not incorporate a stress history nor does it allow for large 1 deformations. Other definitions of stress tensors in use are the second Pio la-Kirchhof f pseudo stress tensor, the J a u m a n n stress rate and the Green-Naghdi rate [135] which are less restrictive than the Cauchy stress tensor. T h e stress and strain tensor are related by the constitutive equations for a material which are discussed next. 2.2.3 Constitutive Equation Const i tut ive equations relate stresses to strains and are material dependent. T h e constitutive equation of linear-elastic material is given by the generalized Hooke's law: E q u a t i o n 2.5 contains a m a x i m u m of 36 independent material parameters due to the symmetry of stress and strain tensor. However, for homogeneous isotropic materials the m a t r i x of elastic coefficients has only two independent variables; e.g., the Lame constants A and /j,, or alternatively the Young's modulus E — yi(2p + 3X)/(fi + A) and the Poisson's ratio v = A/(2(/i + A)), or the Shear modulus 7 = /i. and the Poisson's ratio v\ the latter is used in this thesis. E q u a t i o n 2.6 is the tensor notation of the spring constant of the generalized Hooke's law w i t h 5ij being the Kronecker delta (Sij = 1 if i = j;0 otherwise). For comparison, E q u a t i o n 2.7 expresses the spring constant i n m a t r i x notation (a. = C e). 3 3 k=\ 1=1 (2.5) Cijki = z 7z-Sij5ki + 7(8ikO~ji + Sudjk) (2.6) Here, large implies geometric non-linearities can not be ignored. 21 c = 2 7(l-i,) \-1v 271/ \-2u 271/ \-2v 0 0 0 l-2u 2 7 ( l -« / ) \-2v 2-yv \-2w 0 0 0 2yu \-2v 2-yu l-2i/ 2 7 ( l -» / ) l-2i/ 0 0 0 0 0 0 2 7 0 0 0 ' 0 0 0 2 7 0 0 0 0 0 0 2 7 (2.7) 2.2.4 Navier's Equations In general, the forces acting on an infinitesimal material cubic element under de-formation are stresses and volume-dependent body forces bj, as well as acceleration forces p-iij. The dynamic equilibrium is given by f'"> : Ox, • k> and the static equilibrium then is given by dcji dxi + bj = 0 ( 2 -Equations 2.8, 2.2 and 2.5 can be manipulated to derive the Navier's equations. For a homogeneous isotropic material, Navier's equations are given by Equation 2.9. di dxudxi + 7 duk 2v dxkdxj + bi = 0 (2.9) 2.3 Solution Methods Navier's equations (Equation 2.9) can be numerically solved in at least two ways: a volumetric approach based on energy conservation and a boundary approach based on the application of Green's theorem. The potential energy of the deformation over 22 the discretized domain of the solid can be numerically evaluated wi th the F in i te Element M e t h o d ( F E M ) . T h e result of the F E M is a block matr ix for the solid's stiffness which relates the applied forces to the solid's deformation. T h e method of condensation yields a discrete Green's functions mat r ix from this stiffness matr ix . T h i s is briefly outl ined in Section 2.3.1. App l i ca t ion of Green's theorem to the elastic equi l ibr ium equation leads to a boundary integral equation, which can be evaluated by a weighted residual statement. The discretization of the residual statement is the starting point of the direct Boundary Element M e t h o d . T h e B E M solution directly leads to the mat r ix of Green's functions. Th i s approach is described in some detail i n Section 2.3.2. 2.3.1 Finite Element Model The v i r tua l work of elastic forces can be derived from Newton's second law. In the following body forces are ignored for clarity but they could be included tr ivial ly. Equat ion 2.8 relates the forces on an element of an elastic solid. Integration of the equation of equi l ibr ium over the domain of the deformable body yields Employ ing the Cauchy stress formula and divergence theorem one can arrive (see T h e v i r tua l work of the elastic forces is the volume integral of the product of the stress tensor and the strain tensor. [118, pp. 247-249]) at 5W 23 = / eTCTdedQ = / ( D u ) T C D d u d f t (2.10) Equa t ion 2.10 can be made discrete by employing finite volume elements throughout the domain Q of the solid. T h e displacements and forces are approx-imated wi th polynomia l t r ia l functions over the elements. One choice of a finite element is a tetrahedral element wi th four hat-functions centered at the vertices of the tetrahedron (see Figure 2.3). Note that the boundary of the solid is represented as a tr iangular mesh and the displacement field and force field over the boundary is represented wi th hat-shaped basis functions. T h i s is identical to the Boundary Element M o d e l described in Section 2.3.2 (with the exception that it is customary to employ tractions instead of forces in the B E M ) . 0 Interior Vertices O Exterior Vertices Q Exterior Vertices Figure 2.3: Discrete Domain and Boundary (2D) (a) F E M (b) B E M The F in i te Element M o d e l for an elastostatic solid has the form of Equa -t ion 2.11. The displacement block vector fi is related to the force block vector f by the stiffness block mat r ix K . K u = f (2.11) T h e displacement and force block vector can be reordered into block vectors associated wi th the interior of the solid uo. and f^ and into vectors associated wi th 24 the boundary up and fr and it follows (2.12) The displacement of the boundary ur may be calculated from by condensation of Equation 2.12 as K r , r K r , n u r fr K n , r fQ K r,r ur = fr - K r . n K ^ f n • (2.13) In the scenario of this thesis where internal forces do not occur and internal displace-ments are irrelevant the condensed Equation 2.13 becomes K u p = fr- Here, K can be reordered into a discrete Green's functions matrix S by inverting the elements of block columns of K associated with unknown displacements. The Finite Element Method with condensation has been applied in real-time simulation by Bro-Nielsen and Cotin [18]. 2.3.2 Boundary Element Model The boundary element method (BEM) is a discretization of the boundary integral equation for an elastic solid. The boundary integral equation can be derived with the method of virtual work. This derivation will be quickly outlined in this section. Next, the discretization with the linear triangular boundary elements is employed to derive a discrete Green's functions matrix. Equation 2.8 for the static equilibrium holds over the complete body. As mentioned before, in the course of this thesis, we will ignore body forces b inside the solid 0 . Multiplying Equation 2.8 with a virtual displacement field u* yields Equation 2.14. The right-hand-side boundary terms of Equation 2.14 are the work of the virtual displacement on the boundary. Note, that the virtual displacement plays a role as a weighting term. 25 / ^ * u * d f t = f (pk - pk)u*kdT + / (uk-uk)p*kdY (2.14) JU OXj JY1 JTO Integrating Equa t ion 2.14 by parts twice, as well as employing the Cauchy stress formula 2.4 yields udtt = - / pku*kdT - / pku*kdT + / ukpkdT + / ukp*kdT. Jr-i Jr0 JTi Jr0 In dxj Given the fundamental solution for da*k/dxj + A = 0 where A is an unit impulse, we arrive at the boundary integral equation 2.15. lvu*j^,x)Pl(x)dr(x) = c , V ( 0 " , ( 0 + JrPUZ,x)u3(x)dr(x) (2.15) Equa t ion 2.15 is the boundary integral equation wri t ten as a weighted resid-ual statement over the complete surface T of the elastic body. The location of the load application is £, while the response location is x a n d the distance between them is r. Coordinate i is associated wi th the point load, while j is the coordinate of the response. T h e weights of the equations, the terms p*j(£,x) a n d u*j(£,x)> a r e Kelv in ' s fundamental solution for the 3-D Navier 's equations for an infinite body wi th the same material properties subjected to an impulse point load. T h e fun-damental solution is Equat ion 2.16 and Equa t ion 2.17 where the surface normal is denoted by n and 5ij is the Kronecker delta. Pijfox) 1 167r(l — u)jr - 1 8TT(1 - v)r2 (3 - Au)5ij + dr dr dxi dxj (1 - 2v)8ij + 3 dr dr \ dr dx\ dxj J dn dr dr - { l - 2 u ) { ^ - d x - n ' ] (2.16) (2.17) 26 More details on the derivation of the boundary integral equation can be found in many textbooks including [16, 33, 17]. T h e integrals in Equa t ion 2.15 can be numerically evaluated over linear triangular boundary elements w i th Gaus-sian quadrature. In the following the traction vector field and the displacement fundamental solution tensor are wri t ten as p * ( £ , x) a n d u * ( £ , x) respectively. T h e left-hand side of Equa t ion 2.15 at a mesh node k is ( G P ) f c = ^u*(a,x)p(x)^r(x). T h e right-hand side of Equa t ion 2.15 at a mesh node k is ( H u ) f c = ^p*(&,x)u(x)dT(x). T h e traction vector field p ( £ ) over the boundary is approximated wi th shape func-tions p(x) = JliLi <f>i{x)Pi where p^ are the traction values at a node. Individual 3 x 3 block elements G f c m of the fundamental solution mat r ix are Gkm = J u*^k,x)4>m(x)dT(x)- (2-18) Indiv idual 3 x 3 block elements H f c m of the fundamental solution matr ix are H f c m = jV(£fc,x)<M^r(x). (2.19) The integrals are numerically evaluated by Gaussian quadrature. The com-putat ion is performed by a change to local barycentric coordinates for each element. A p p e n d i x A contains details of these steps. The discretized boundary integral equa-tion is the matr ix equation H u = G p . T h e equation has to be rearranged for the specific boundary condition at hand. The prescribed traction and displacements are collected in the block vector v while the right-hand side is the complementary set of tractions and displacements v . Exchanging elements in vectors p and u to yield 27 v and v requires exchanging block columns in matrices G and H to yield G and H . T h e modified equation H v = G v is solved wi th the Green's functions mat r ix S as follows: G ^ H v = S v = v. This completes the derivation of the Green's functions matr ix . Next ap-proximate physics-based models are compared to the formulation of this thesis of Equa t ion 2.1. 2 . 4 O t h e r P h y s i c s B a s e d D e f o r m a b l e M o d e l s Physics based deformable models have been employed extensively in the computer vision, computer graphics and medical imaging communities. However, the term physics based modeling is not employed consistently. In this Section, we briefly describe the central equations in two of the most common approximate models. In Section 2.4.1 the set-up of particle systems is described, while in Section 2.4.2 the model of deformation in deformable superquadrics is summarized. 2.4.1 Particle Systems The modeling of elastic solids by particle systems has been applied in computer graphics since 1981 [100]. T h e particle system approach divides the solid into in -d iv idua l mass points inter-connected by springs (see Figure 2.4). Usual ly the mass points are arranged in a regular lattice borrowing from the atomic structure of crys-tals. The springs are typical ly linear and usually a viscous damping term is added. A linear elastostatic solid modeled as a particle system gives rise to an equa-t ion Ku = f [40]. The stiffness mat r ix K is composed of the indiv idual spring 28 Figure 2.4: Mass-Spr ing or Par t ic le M o d e l stiffnesses and their rest lengths, u is the displacements of the mass points from their rest state and f is the external force vector. The stiffness mat r ix is banded determined by the connectivity of the lattice. Various researchers attempt to set the spring constants to physical properties. Lee, Terzopoulos and Waters set the spring constants of their facial models according to typical human anatomy [73]. T h e actual constants are arbi trary while the material dis t r ibut ion (e.g., muscle, facial tissue, epidermis etc.) is determined by the typica l human anatomy. K o c h et al . set the spring constants according to the material dis t r ibut ion of a ind iv idua l human face based on a C T scan [65]. A g a i n the actual constants are arbitrary. Reznik and Laugier tune the spring in their model of a rubber t ip of a robotic finger according to the Young 's modulus of rubber [106]. Identification from observation for dynamic systems is described in [60]. In [77] a mass-spring cloth model is fit to observation of hanging cloth behavior. D ' A u l i g n a c , Cavusoglu and Laugier [27] model the elastic behavior of the human thigh wi th linear and non-linear springs from observations. T h e observations consist of force and displacement measurements in response to pressure probing wi th a robotic arm. A general difficulty wi th mass-spring models is the lack of a suitable general-purpose tuning algori thm for the springs [40]. T h e springs have no physical connec-t ion to the behavior of an elastic solid. Mass-spring models do not model various 29 physical behavior of solids well including incompressibility and large stiffness. How-ever, mass-spring models can be rendered fast if the springs are'reasonably soft and have been utilized in computer graphics with success. 2.4.2 Deformable Superquadrics Terzopoulos formulated regularization terms based on physical analogy of a mem-brane under tension and a thin plate under tension [124]. Metaxas and Terzopoulos extended this work into deformable superquadrics and introduced it to the com-puter vision and computer graphics fields and later to the medical area [82, 81]. The deformable superquadrics of Metaxas and Terzopoulos are based on Lagrangian dynamics. Their stiffness matrix K combines global and local stiffness coefficients. The position of a point on the surface of the model relative to the center depends on the global deformation g and the local deformation 1, i.e., u = g + 1. This needs to be expressed in generalized coordinates in order to be utilized in the Lagrangian dynamics formulation. The generalized coordinates chosen by Metaxas and Ter-zopoulos for the global deformation are the parameters of a superquadric. Also, the center of mass of the superquadric and the rotation angles serve directly as generalized coordinates. The generalized coordinates of the local deformation d are not independent of other generalized coordinates which introduces a time-varying Jacobian. Metaxas and Terzopoulos represent the local deformation of a body with finite elements over the surface of the solid. They use bilinear quadrilaterals with four nodes and right-angled triangles with three nodes. All of these nodes are on the surface of the solid. In particular, this local finite element model is not 30 a discretization of the physics of a deformable solid. Therefore, the deformable superquadrics approach is not compatible w i th this thesis: 2.5 Summary T h i s chapter details the deformation model employed by this thesis. T h e model is a discrete cont inuum mechanics model of linear elastostatic object behavior, ex-pressed i n terms of a Green's functions mat r ix which can be derived by the Bound-ary Element M e t h o d for homogeneous isotropic linear-elastic material . However, the Green's functions mat r ix representation does not depend on the B E M ; in particular, F E M w i t h condensation leads to similar formulation. The definitions and assump-tions i n cont inuum mechanics for a linear elastostatic solid have been summarized in this chapter. L imi ta t ions of the model are due to geometric approximation, i.e., first-order deformation and material approximation, i.e., linear elasticity. T h e remainder of this thesis utilizes Equa t ion 2.1 in the estimation task. The estimation task is finding the discrete Green's functions matr ix based on observation of vertex displacements and tractions at a vertex. 31 C h a p t e r 3 T h e A C M E F a c i l i t y T h e Act ive Measurement Faci l i ty ( A C M E ) is an integrated robotic facili ty designed to acquire measurements of interactions wi th objects. A t the core of A C M E is a contact manipulat ion system ( C M S ) which executes contact interactions wi th an object to be measured. These interactions can be recorded wi th various sensors including a force-torque sensor, a t r inocular stereo system, a microphone and a high quali ty 3 - C C D color video camera. Sensors which are employed to measure the response of an object from a distance are combined in a field measurement system ( F M S ) . A C M E is implemented as a robotic facility because of the abi l i ty of robots to perform tasks wi th high repeatability and consistent accuracy. T h i s is extremely important if a large number of registered measurements of an object are to be collected. T h e test station is a posit ioning robot to move the object to be measured wi th precision. The C M S is buil t around a robot arm positioned by a motion stage. A gantry robot places the F M S relative to an object. A l l these robots are commercial devices but their assembly, control and programming are unique to A C M E . A s a result A C M E is a tele-robotic system wi th a total of fifteen degrees of freedom. 32 Figure 3.1: A C M E Facil i ty Overview T h i s thesis employs A C M E for the acquisition of measurements dur ing active object deformation. W h i l e the design and implementation of A C M E is a group ef-fort, some features are especially designed to enable the acquisit ion and verification of deformable models. These feature are discussed in some detail in this chapter, while only an overview is provided about A C M E in general. T h e remainder of this chapter is organized as follows. In Section 3.1 a system overview is given including the actuators and sensors of A C M E (Section 3.1.1), the control architecture (Sec-t ion 3.1.2), the programming of A C M E (Section 3.1.3) and the sensor programming model (Section 3.1.4). T h e structure of an A C M E experiment direct ing A C M E to actively deform an object and record the deformation is discussed in Section 3.2. T h e geometric model ing capabil i ty of A C M E is summarized wi th examples uti l ized in this thesis in Section 3.3. 33 3.1 O v e r v i e w The design goal of A C M E is to make it efficient to build realiiy-based models [92], i.e., interactive computational models of real, physical objects based on actual measure-ments. The physical layout of the facility with its variety of sensors and actuators is shown in Figure 3.1. The grouping of these sensors and actuators (ACME de-vices) is outlined in Section 3.1.1. The control of these devices is distributed over a number of computers and layers. This allows one to employ commercial con-trollers when possible but adding advanced features as required. In particular, the control architecture provides the foundation for the high-level control of the sys-tem which makes programming it simple. The architecture is based on tele-robotic control and is described further in Section 3.1.2. The simplicity in programming arises from an environment which integrates a kinematic simulator with the tele-robotic paradigm in a number of ways. On the one hand the simulator serves as a debugging tool for A C M E experiments. An A C M E experiment is the top-level user program written in Java1. On the other hand the simulator is available during execution to verify motion requests. This allows for a simple generate-and-test ap-proach to motion planning. Furthermore, the user programming of A C M E is aided by a client program which can execute independently from the facility or act as a terminal for ACME. The client in stand alone mode helps considerably in sharing A C M E and eases program development. Section 3.1.3 contains more information on the programming environment, while Section 3.2 discusses the A C M E deformation experiment. The interested reader is also directed to [94, 93] for a more complete description of ACME. 1 Java is a trademark of Sun Microsystems Inc., MountainView, Ca., USA 34 Figure 3.2: Measurement of contact force and posi t ion w i t h C M S 3.1.1 A C M E Devices T h e actuators and sensors w i th in A C M E are collectively referred to as devices. T h e y are conceptually d iv ided into the following subsystems (see F igure 3.1 for an overview): • The Contact Manipulation System ( C M S ) is based on a P u m a 260 robot a rm wi th a wrist mounted 6-axis force/torque sensor ( A T I M i n i 40). T h e P u m a 260 has been used to apply sustained contact forces of a least ION to 15N. The A T I M i n i has a m a x i m u m force l imit of 40/V to 120vV (depending on the axis) and a m a x i m u m torque of 2Nm w i th a resolution of 0.047V to 0.12A r and l.ONmm, respectively. A 1.267m(4/t.) long linear motion stage moves the robot a rm to and from its application area and increases the robot arm's workspace. T h e control of the robot arm has been modified to enable the application of forces to a test object. T h e robot a r m can work i n impedance control mode or execute guarded moves moni tor ing impact forces at its t ip . The C M S is capable of measuring object properties such as friction and local stiffness and of apply ing a force impulse to generate contact sounds. Figure 3.2 shows the C M S acquiring force data while in contact wi th a test object. 35 Figure 3.3: A C M E F M S Sub-System • The Field Measurement System ( F M S ) holds sensors for measuring the light and sound fields around the test object (see Figure 3.3). T h e unit is positioned wi th a gantry robot (a luminum structural elements by "Item Industrietechnik und Maschinenbau m b H " powered by Parker Hanni f in Corp . brushless D C motors) and oriented by a two degree-of-freedom pan-and-ti l t uni t (Sagebrush Model-20). T w o vision sensors are part of the sub-system: a h igh quali ty 3 - C C D color video camera wi th computer controlled zoom and focus (Sony D X C 950) and a t r inocular color stereo vision system (Point Grey Research Color Triclops) . A microphone is also attached to the unit . • T h e 3 - D O F Test Station presents the test object to the sensing equipment. The test-object can be mounted r igidly to the test station which positions the object in the plane. The test-station consists of x- and y-translat ional motion stages, as well as a rotation stage. The stages (by Daedal) possess a linear and rotary resolution of ± 0 . 0 0 6 3 5 m m and ± 1 0 arc-min, respectively. • Miscellaneous devices include a time-of-flight laser range finder (Acu i ty A c c u -Range 3000 L I R ) and various lights which can be switched under computer control. 36 A C M E Client Robot A r m Force Sensor Motion Stage Figure 3.4: A C M E Server: Overview of Cont ro l 3.1.2 Con t ro l Archi tec ture The software and hardware control architecture of A C M E is briefly outl ined in this Section. A n overview of the architecture is shown in Figure 3.4. T h e control is client-server based, i.e., an A C M E client w i th which a user interacts and an A C M E server which controls the robotic subsystems introduced above. The A C M E server system is dis tr ibuted between several computers; separating low-level closed-loop control, trajectory generation, data acquisition, and networking. The four layers of the A C M E server software from top to bo t tom are the 37 user's Experiment, the Java A C M E device classes, the native bindings to the devices and finally the low-level device controller. For an example of the top-most layer, the user's Experiment, see Section 3.2. The remaining layers are described below. The layer of high level Java objects called Devices grouped into Actuators and Sensors is next. These Actuators are commanded to move by Motion objects. This is important, since it decouples the functionality of a specific actuator from a motion specification. Motions are grouped into a MotionPlan which is executed as a whole. This design feature allows a MotionPlan to be verified, actuators to be synchronized and run in a common but separate thread. The Actuator objects match the physical actuators of A C M E described in Section 3.1.1. The programming of Sensors is separately described in Section 3.1.4. The next lower layer of software is the native bindings to devices. The native bindings make code exclusive to the A C M E server and is neither available nor necessary to the A C M E client. These layers can be replaced by a simulator. The actuator native bindings are through the Robot Control C Library (RCCL) [75]. RCCL generates real-time trajectories for each actuator device from a MotionPlan. The set-points of these trajectories are sent to the actual controllers of the various actuators. The Puma is controlled by Unimation's control unit while the rest of the devices are controlled by an 8 axis motion control card (Precision Micro Dynamics MC-8) except for the RS-232 interface to the stepper motor of the FMS's pan-tilt head. 3.1.3 Programming A C M E The A C M E client is the stand alone user interface to ACME. The client is the pro-gram development environment for A C M E and is composed of a three-dimensional 38 View Figure 3.5: A C M E 3-D Display A collision event produced by the simulation is rendered wi th the col l iding boxes highlighted. display of the facility, a simple editor and access to the Java compiler . T h e client has access to the complete A C M E Java code which includes a kinematic simulator of the facility. Th i s allows the user to develop, debug and verify a l l Experiments without access to the A C M E server. Once a compiled Experiment is ready to be executed on the A C M E server, the A C M E client connects to the server over the In-ternet, uploads the compiled Experiment and receives any data requested. Because an Experiment is any Java object which adheres to the Experiment interface, the user has complete freedom to implement any desired routines. T h e great advan-tage of this model-based tele-programming approach [76] is that a l l control-loops are closed solely on the server side rul ing out communicat ion delays between server and client. D a t a processing routines can also be run on the server which allows data to be reduced to its essential size before being transmitted to the client. T h e simulator models the entire A C M E facility w i th an enhanced Java3D 39 scene graph. The scene graph contains the geometric, posit ional , and appearance information used in simulation on the client, run-time checking on the server and graphic display on either server or client. A n important feature of this scene graph is a simple on-line collision detection using oriented bounding boxes [42]. Bound-ing boxes are placed around pr incipal system components and can be tailored for different experiments. Figure 3.5 shows an example of a collision event rendered on the graphic display of the client. 3.1.4 Sensor Model A C M E is a robotic system designed to collect data from sensors. T h e programming abstraction of these sensors is an important design feature of A C M E . T h e high-level description of sensors is related to the description of actuators discussed in Section 3.1.1. Sensors in A C M E are widely distr ibuted on the computers that are the A C M E server. A l l actuators provide a position sensor functionality which al-lows the recording of an actuator's posit ion at 100Hz. Th i s functionality is t ightly coupled wi th the R C C L layer of the server. Th i s is also true of the force sensor of the robotic arm. The video digitizer of the color 3 - C C D camera and its serial command interface on the other hand are quite separate from A C M E server. The same is true for t r inocular stereo-head, the laser range-finder and the microphone. Nevertheless, al l these devices are A C M E Sensor objects. Th i s is achieved by adding another client-server layer. A n ind iv idua l SensorServer executes on the computer hosting the sensor hardware, while the A C M E server acts as a client to these in-d iv idua l SensorServers. T h i s design not only unifies the access to sensors from an A C M E Experiment but it also enables asynchronous real-time data acquisition. A n Experiment specifies SensorTrigger events which are uploaded and checked on the 40 server of the sensor. After being started, the server of the sensor continuously ac-quires data checks the trigger condit ion and buffers data to be sent to the A C M E server and eventually to the A C M E client. In this fashion, communicat ion delays are eliminated from the data acquisition loop. For example, the Triclops server incorporates the control of the video dig-itizer, the Point Grey Research Color Triclops control software and access to the Point Grey Research Triclops stereo vision S D K (software development k i t ) . T h i s allows the user to specify the stereo processing to be executed on the Triclops server. Trigger events can take full advantage of the stereo library, e.g., by monitor ing depth in the stereo imagery, or implement whatever image processing is required. A s the trigger activates, the SensorServer may proceed to send raw images or images which are the result of processing in various formats. 3.2 A C M E Deformation Experiment A n A C M E deformation experiment as employed i n this thesis applies a contact force to an object wi th the C M S while monitoring the global deformation wi th the F M S . Th i s section details the methods uti l ized in p lanning the contact wi th the object, executing the probing and acquiring measurement data. T h e structure of the A C M E DeformExperiment is described here in detail since it is the Experiment used in this thesis to acquire measurements and it also serves as an example of an A C M E Ex-periment. The Experiment does assume that a geometric model of the object under test is available. T h e next Section 3.3 describes how the models for this thesis are acquired. Also not part of the DeformExperiment are the measurement processing routines. Th is separation into data acquisition and analysis was convenient for the research methods explored in this thesis. T h e processing of measurements is de-41 ScanPoistion Experiment TestObject * DeformExperiment Pokelt A FmsPositioning CmsPositioning TsPositioning PokeGoalPlanner GuardedMoveinZ SinglePoke extends Class ( Interface ) \ Abstract Class" Figure 3.6: Structure of an A C M E Deformation Exper iment scribed in the following chapters, Chapter 4 to Chapter 6. T h e overall structure of a deformation experiment is depicted in Figure 3.6. The DeformExperiment takes the specification probing, which is object-specific, from a data-base object which implements TestObject. Examples of speci-fication include the m a x i m u m force exerted on the object while probing, the max-i m u m probing depth, the number of steps for the deformation measurement (see Chapter 4) and various files and directories to store data in . The other data-base object implements the ScanPosition interface which contains pre-planned viewpoints and calibration information for the trinocular stereo-head wi th in the F M S . Based on these specifications the DeformExperiment positions the F M S , initializes the PokeGoalPlanner and starts the probing process. A simplified execution flow of the DeformExperiment is shown in Figure 3.6. T h e probing process executes as follows. T h e PokeGoalPlanner generates 42 1. Posi t ion TestStation, FMS and CMS 2. Initialize PokeGoalPlanner 3. Loop over al l positions specified in ScanPosition • Loop over al l specified vertices to be probed — Request probe specification from planner — W h i l e probing not successful * Simulate probing * Execute GuardedMovelnZ * Execute SinglePoke * If not successful change configuration — Process and save a l l da ta 4. Move FMS and CMS to park posit ion Figure 3.7: Simplif ied Execut ion F low of A C M E Deformation Exper iment the next posit ion and orientation for the probe based on the mesh of the object under test and the specification in TestObject. It checks l imits in order to rule out goals which would br ing the robot-arm too close to the TestStation. It also does ray-intersection w i t h the object in order to verify that the probe is contacting the object at the expected location. T h e n the DeformExperiment calls the s imulat ion wi th the necessary motion to reach the goal in a default configuration. O n success it executes this motion. O n failure the configuration is changed wi th a heuristic, either changing the orientation of the object by rotat ing the TestStation or modifying the distance of the robot from the object by moving it wi th the long linear stage. A successfully simulated motion is executed in a three-step process. F i rs t an approximate position away from the object is reached. Second, a GuardedMovelnZ is executed to find the actual posit ion of the surface of the object. A l l measurements are made relative to this posit ion reducing the influence of posit ioning errors of the robot arm. T h i r d a SinglePoke is executed which contacts the surface, applies force, records force and 43 displacements, and triggers the image capture by the stereo-head. T h i s experiment is run in s imulat ion when it is applied to an object for the first t ime. The simulation performs collision checking and the rendering of the mo-tion provides a visual confirmation of the experiment. If this performs satisfactory, A C M E acquires the deformation data autonomously. 3.3 G e o m e t r i c M o d e l i n g w i t h A C M E Geometric modeling of objects in A C M E is necessary for various tasks including contact planning, object modeling and registration. T h i s thesis relies on A C M E being able to acquire geometric models of the object under test. F i r s t , a surface mesh of the object under test is required as discretization for the discrete Green's functions mat r ix described in Section 2.1. A surface mesh is also required for the deformation data acquisition described in Section 3.2 above. Furthermore, the sur-face mesh is also uti l ized in the deformation measurement described in Section 4. Final ly , the estimation of the discrete Green's functions mat r ix in Chapter 5 and Chapter 6 relies on the discrete surface representation as well . It is evident that geometry is essential to this thesis. Despite this, no novel method for in i t ia l model acquisition is proposed. Unfortunately, geometry acquisition techniques are s t i l l far from being fully automated and general. Dur ing the research for this thesis distance measurement w i th a time-of-flight laser range-finder has been investigated [70], as well as a novel stereo-vision matching technique [69]. In the end, a collection of techniques based on data obtained from the tr inocular stereo-head have been found to be effective. The techniques succeed in obtaining a closed (or water-tight) mesh wi th subdivision hierarchy (see Figure 3.8) for the objects in this thesis. T h e process is only semi-automated and some mesh edit ing is necessary to arrive at a water-tight 44 multi-resolution surface mesh. A detailed description of the techniques can be found in A p p e n d i x B . 3.4 S u m m a r y T h e A C M E facility is a robotic measurement facility which enables the execution of experiments in which a probe contacts an object. T h e probe is mounted on a robotic arm equipped wi th a force sensor. T h i s contact manipulat ion system ( C M S ) is able to exert a force onto the surface of an object at a specified location. T h e tr inocular stereo-head wi th in the field measurement system ( F M S ) can be placed to observe the object dur ing deformation from various viewpoints. A robot such as A C M E can execute repetitive tasks wi th high repeatability which is important dur ing exhaustive contact sampling of an object's surface. The user interface of A C M E combined wi th its tele-robotic design make programming contact interactions convenient. A further feature of A C M E necessary for this thesis is the abi l i ty to produce three-dimensional point-clouds registered in an object reference frame. These point-clouds are combined in a mesh generation step which produces a hierarchy of meshes wi th subdivis ion connectivity which can also be texture-mapped for better viewing. 45 (a) Raw Scan Mesh (b) Level 0 (149v, 294f) J O B (c) Level 1 (590v, 1176 f) (d) Level 2 (2354 v, 4704 f) (e) Level 3 (9410v,18816f) (f) Level 2 Mesh Texture-Mapped Figure 3.8: Mesh Reconstruction The raw scan mesh 3.8(a) wi th a total of 15076 vertices(v) and 29095 faces(f) including various external stray polygons. The subdivis ion hierarchy (3.8(b), 3.8(c), 3.8(d) and 3.8(e)) is generated using displaced subdivis ion surface style piercing of the raw scan mesh. 46 C h a p t e r 4 D e f o r m a t i o n M e a s u r e m e n t The deformation measurement task is to estimate how the vertices of the unde-formed triangular mesh are displaced during applicat ion of a load on the object. Displacement measurement has been approached i n various ways. Shape estimation and tracking of deformable objects has been intensively studied in computer vision (e.g., [127, 53, 34, 49, 62, 74]). In mechanical engineering, optical measurement of strain has been developed (e.g., [83, 46, 80]), as well as commercialized (e.g., G O M m b H : Gesellschaft fur Optische Messtechnik, Braunschweig, G e r m a n y 1 ) . In mechanical engineering, the strains considered are usual ly small and often only the two-dimensional result is of interest. However, while strain can be measured d i -rectly visually, it may also be approached in other ways. For example, in s tructural mechanics, strain gauges give local information at cr i t ical points of the structure. Th i s chapter describes the range-flow approach to deformation measurement of this thesis. Range-flow is the recovery of scene changes from three-dimensional images. It clearly solves the measurement task at hand since three-dimensional information is recovered. lSee http://www.gom.com 47 T h e task in range-flow is to match depth information over time. T h e range-flow constraint equation [51, 132] for the depth of a surface Z(x,y,t) is dZ , dZ , ,„ dZ , n ——ax + —— ay — dZ + -7—at = 0. ox oy at T h e range-flow constraint equation is the three-dimensional equivalent of the opt ical flow equation [52]. Similar to optical flow [6] this constraint may be ut i l ized by area-based matching or directly wi th a differential technique. Th i s thesis employs area-based matching. Es t imat ing range-flow is a harder problem than optical flow because of the increased dimensionality and the necessity for quali ty range data. Nevertheless, there has been progress in range-flow estimation [132, 120], sometimes also referred to as scene flow [128]. T h e approach of this thesis is to calculate range-flow based on a combination of simultaneous 2D-opt ical flow wi th stereo-depth information. Simultaneous stereo depth and opt ical flow goes back to at least Ba l l a rd and K i m b a l l [3]. Recently this approach has been followed successfully in [128, 133, 119]. However, in this the-sis, the addit ional imagery is uti l ized to reduce the error in the very noise-sensitive optical flow. Th i s results in a purely local method and no global regularization is employed. T h e high quali ty calibration of the commercial t r inocular stereo sys-tem [105] in A C M E lets this approach succeed. T h e choice of stereo vision and opt ical flow requires objects w i th sufficient visual texture. Sufficient texture is texture which leads to an unique peak in the area-based correlation at the correct dispari ty in the stereo matching and at the correct location in the search area for optical flow. T h i s condition can be met since A C M E is a (robotics) laboratory environment and most objects either have sufficient visual texture or can have visual texture applied to them tr ivial ly. In Section 4.1 methods to robustly estimate flow are briefly reviewed and in 48 Section 4.1.1 the approach of this thesis to robust range-flow estimation is outl ined. T h e computations involved in this method are detailed in Section 4.1.2. Results for r igid body mot ion compared to ground t ru th data are given in Section 4.1.3, while examples for deformable objects are shown in Section 4.1.4. Range-flow results are not direct measurements of the vertex displacement u^, needed for the estimation task defined in Section 5.1. These displacements need to be estimated based on range-flow results. Th i s estimation is described i n Section 4.2 and results are given in Section 4.2.1. A summary in Section 4.3 concludes this chapter. 4 . 1 R o b u s t F l o w E s t i m a t i o n This section gives a selective review of range-flow and optical flow estimation in light of the deformation measurement task at hand. Mos t optical flow techniques can be characterized either as a differential technique or as a region-based matching tech-nique [6]. Differential techniques use the opt ical flow constraint equation which is underdetermined (one equation in two unknowns). L o c a l averaging or a regulariza-t ion term need to be introduced in order to find a solution. For example, the original technique of H o r n and Schunck [52] uses a global smoothness constraint. Area-based matching techniques, on the other hand, impl i c i t ly average the flow wi th in the area being matched and no further constraints are s t r ic t ly necessary. Nevertheless, the results typical ly are noisy and need further processing [6]. Th i s processing may in-clude introducing a pr ior i constraints on the solution v i a regularization (e.g., [68]), fitting results to a parametric model of motion (e.g., [55]), or filtering flow results based on a rel iabil i ty measure (e.g., [2]). M u l t i p l e approaches are often used in combination. Since the displacement to be measured is between the undeformed and de-49 formed object states, only two images are necessary. However, accurate derivatives can only be calculated from a sequence of images (unless the signal is highly over-sampled or the intensity structure is nearly linear [6]).. Lack of accurate derivatives is identified a major source of error in differential approaches [6]. Area-based tech-niques do not require estimates of the t ime derivative of intensity and depth in opt ical flow and range-flow respectively. Therefore, an area-based optical flow ap-proach is implemented to avoid having to deal w i th a sequence of images. Stereo-range data necessitates at least two cameras but most systems today use either a set of stereo camera pairs or at least three cameras [63, 89, 105]. C o n -straints from addit ional imagery are also exploited in the computat ion of optical flow in this thesis. The integration of results from all stereo cameras allows the rejection of outliers, especially if they are generated by viewpoint dependent effects or if they are the result of ambiguous matching results. Specular highlights (see e.g., [13]) are a viewpoint dependent effect that produces noise in the flow estimate. Another source of noise is the lack of texture in the matched area. T h i s can be identified using multiple cameras since the peak in the correlation surface w i l l vary between cameras. The peak w i l l also be inconsistent at motion boundaries asso-ciated wi th depth discontinuities. The peak in the correlation surface w i l l be in different locations since the percentage of area in the matching window support ing the different motions wi l l vary between cameras. M u l t i p l e cameras also help to as-sign a confidence measure for optical flow results which has been identified as a hard problem [6]. T h e addit ional data from mult iple views does not prevent the ut i l izat ion of other approaches to compute a robust optical flow result. In particular, the ro-bust optical flow techniques of [13] (multiple motions wi th in an aperture) or the 50 regularization method based on the sum of squared differerences (SSD) of [68] (min-imizing the effect of discretization errors in the flow) would be compatible with the approach here. However, for the purposes of this thesis, the result obtained from • multiple views is sufficiently robust to serve as input to the vertex displacement estimation stage of Section 4.2. The surface deformation measurement technique utilizing a trinocular stereo vision system of this thesis is described in Section 4.1.1. 4.1.1 Robust Range-Flow by Additional Imagery . Robust range-flow by additional imagery tracks visual surface texture in three spatial dimensions, i.e., it calculates the range-flow of the deformable object's surface. The stereo-disparity between views serves as range information as well as to register the flow calculation. Stereo and optical flow results are both calculated with an area-based matching technique using all three cameras of the stereo-system. An overview of the range-flow calculation is shown in Figure 4.1. In the following, the approach is described and design choices in its implementation clarified. Section 4.1.2 explains implementation details. In the algorithm, a geometric model of the undeformed surface of the object under test is either assumed known or acquired beforehand. This geometric model and the tracking capabilities of the system enable the segmentation of the three-dimensional range-flow into flow at the surface of the object and other flow. With the object undeformed, a stereo-triplet is captured. The disparity of the images and a range image is calculated. The range-image is taken from a known position and orientation in the A C M E system. Based on this spatial information the closest triangular surface patch for each pixel of the range-image is found. Pixels which are within a threshold of the object's surface are accepted as surface pixels. Optionally, 51 T = 0 Figure 4.1: Range-flow processing: overview at T=0 (initialization) and T = l Range-flow shown is at T=4; at T = l there is no deformation since the probe i only just touching the object's surface. 52 the in i t i a l registration may be improved wi th an iterative closest point ( ICP) [8] type algori thm. However, in our set-up this is not necessary since the in i t ia l guess is sufficient for segmenting the range image. A l l image row, column, dispari ty triplets associated wi th surface pixels are input into the tracker. Next , opt ical flow from the ini t ia l izat ion to the first t ime step in the reference camera is calculated for al l pixels which are tracked. In the approach of this thesis, area-based matching uses the normalized cross-correlation. T h e normalized cross-correlation has been selected because of reports of its good performance in area-based matching [35] and because it can be calculated from correlation scores which are. linear w i th respect to image intensity (see Equa t ion 4.1). T h e opt ical flow in the non-reference camera images is calculated for surface pixels in the coordinate frame of the reference camera. T h i s is achieved by using the sub-pixel dispari ty information as a look-up table (see Section 4.1.2 for more details). Each surface pixel then has an optical flow result for each camera. If al l results are in accordance, the flow result is accepted and the th i rd dimension is added based on the dispari ty and stereo tr iangulat ion. A l l points successfully tracked from ini t ia l izat ion to the first t ime step are considered surface points for the following time step. Before discussing the implementation details,, some tracking issues need to be discussed. It is convenient to record the ini t ia l izat ion followed by two time steps, even though only one is s tr ict ly necessary. T h e ini t ial izat ion is done without the probe touching the object, while in the image triplet taken after the first time step the probe is just touching the surface and in the final image triplet the object is completely deformed. T h e separate ini t ia l izat ion allows the range-based segmenta-t ion to succeed in image areas where the force probe otherwise would be very close to the surface. The first time-step on the other hand allows for the correction of 53 row + disparity Top Camera •a-correlation window column + disparity -D Left Cam Oi 0 columns column - a Reference Camera Figure 4.2: Triclops Stereo Image Triple t any in i t ia l offset introduced by the force probe making contact w i th the object. The short sequence of three image triplets means the tracker does not need to have new points to be tracked added in each step. It just tracks a l l the in i t ia l points through the complete sequence. However, note that after the displacement of vertices has been estimated the deformed mesh could be ut i l ized to add addit ional points into the tracker. 4.1.2 Implementation Details In this section, the implementat ion details of the range-flow calculation are described starting wi th the input images. The images from the stereo cameras are rectified such that the epipolar lines are aligned wi th the image axes (see Figure 4.2). The result of the stereo processing wi th the Point Grey Research Triclops stereo vision S D K is a triplet of row, column and sub-pixel dispari ty at each valid pixel in the reference image. A pixel is valid if the stereo processing produces a result. Va l id pixels which are deemed to be on the object's surface are tracked in the reference image. T h e area-based matching correlates a square window centered 54 at a given pixel wi th a square search area in the next image of the sequence from the same camera. T h e mask size of the correlation window and the search area both influence the result of correlation. T h e influence of these parameters and of the image size on the results and computat ion t ime are demonstrated i n Section 4.1.3. T h e optical flow correlation scores of the non-reference cameras are combined wi th the ones of the reference camera wi th the aid of stereo disparities. The pixels of the start and end of the flow vector in the reference image are mapped wi th sub-pixel disparities into the stereo images. T h e correlation between these locations is calculated wi th Equat ion 4.1. Equa t ion 4.1 is the normalized cross correlation NCC of two image windows at column sub-pixel location using bilinear image intensity interpolation. The NCC is specified wi th pixel a at location row ra and column ca plus 0 < p < 1 and pixel b at locat ion rb and q, plus 0 < q < 1. The NCC is calculated from the non-normalized cross-correlation CC and the self correlation SC. Correlat ion of image intensities I C O r r is computed over a window wi th size 2 * n + 1 where n a positive integer. (4.1) CC(ra,ca,rb\cb,p,q) (1 -p)(l- q)ICOrr{r + (1 -p)qlc0rr(r c fc+l) +p(l - q)ICorr{r Cb) +PqIcorr(r SC(ra,c( P) ( ( ! - P ) ( l ~ P)Icorr(ra,Ca,rt +2.0(1 -p)PIcorr(r ca+l) +ppICorr{ra, c a + i , r a , C a + l ) 55 n n /——n k=—n The result of the optical flow calculation is three correlation scores for each goal location within the search window provided a disparity value also exists. In other words, three correlation surfaces registered with goal locations in the refer-ence camera are found. These surfaces have missing values for locations where no disparity value was found in the stereo processing. The maximum peak of each surface is the flow estimate for a given camera. A voting scheme accepts the result only if all cameras are in good accordance, i.e., the peak of the correlation surface of all stereo images is in the 8-neighborhood of the maximum correlation score of the reference image. This voting scheme simply rejects outliers rather than try to find the correct flow. This strategy is in contrast to the stereo approach of sum of the sum of squared differences (SSSD) developed by Kanade et al. [63]. Summing the correlation scores helps to better define a peak but it is not as robust as voting. It would be worth trying to find the combined peak to sub-pixel accuracy after the voting step utilizing the better definition of the peak with the robustness of the voting. However, pixel accuracy is sufficient for this thesis. 4.1.3 Results for Rigid Body Motion In the following, the range-flow algorithm is evaluated by tracking a moving rigid test object. This choice was made because of the difficulty of obtaining ground-truth for surface deformation. The evaluation is done with both an object translating as well as rotating. In Figure 4.3 the estimated optical and range-flow for translation and rotation of the test object are shown. For simplicity, this section provides comparisons of results for different choices of parameters for the translation only. Note also, these comparisons are intended for relative comparison and to provide 56 (a) Optical Flow Density (b) Range-Flow Distribution (c) Optical Flow Density (d) Range-Flow Distribution Figure 4.3: R i g i d Translat ion and Rota t ion of a Test Object Translat ion by 5 m m is shown in 4.3(a) and 4.3(b), rotat ion by 2.0 ° deg. around the object's z-axis (up in the image) is shown in 4.3(c) and 4.3(d). Brighter range-flow corresponds to larger motion. an indicat ion of the expected error. One can expect results of slightly lesser quali ty during deformation or rotation. Th i s degradation is expected due to a significant change in surface normal of the observed surface relative to the cameras and light sources. Mot ions which w i l l cause optical flow to vary over a correlation window w i l l make the matching task more difficult and introduce error i n locat ing the match. Note, however, that the observed translation is not fronto-parallel to the image plane, i.e., the correct flow is changing over a correlation window. T h e following tables compare the range-flow result to the correct t ranslat ional 57 motion. T h e correct motion is obtained from, reading encoder values of the high-precision linear posit ioning devices in A C M E . The table shows the result for tracking the object over eight t ime steps plus one t ime step for ini t ia l izat ion. Surface points in the ini t ia l izat ion t ime step are tracked, no points are added to be tracked at later t ime steps. Therefore, the number of points tracked decreases w i t h each t ime step. T h e motion between steps is 2 .5mm. T h e tables report the true and estimated values for the orientation and length of the 3D range-flow wi th respect to the ini t ia l izat ion step. T h e reference camera of the stereo-head is 81.0cm away from the center of the object coordinate frame (which is a typical distance in the scenario of this thesis to observe deformation from). Furthermore, al l results are for the test-object in Figure 4.3. The first two tables, Table 4.1 and Table 4.2 show results for two different translations. T h e correlation window size is 7 x 7 pixels w i th a search range of 13 x 13 pixels and the image size is 320 x 240. T h e reported 3D-orientation is the angle between the z-axis (optical axis) and the flow vectors. T h e variat ion in the angle decreases wi th greater translation while the variation in vector length increases wi th greater translation. The decrease in orientational variation indicates the influence of the discretization error in the opt ical flow since larger flow vectors can be more precisely oriented in the discrete image plane. The increase in variat ion in the length of the flow vector is due to a reduction in the density of the flow-field combined wi th larger absolute flow values. The method produces results for translation to about ^ mill imeter accuracy which is sufficient for object deformation measurement. The motion to the left-rear in the reference image (Table 4.1) produces similar results to the motion to the rear-right (Table 4.2). In the following, the translat ion to the left-rear of Table 4.1 is evaluated varying the mask size, image size and search area 58 Step Angle [deg] Length [cm] # M e a n a M e a n a Encoder C o m m a n d 1 67.0081 10.2836 0.2678 0.0340 0.2355 0.2500 2 64.7748 5.7253 0.5396 0.0402 0.4864 0.5000 3 64.8529 4.1890 0.8113 0.0464 0.7374 0.7500 4 64.8575 3.5898 1.0849 0.0571 0.9911 1.0000 5 64.9783 2.8942 1.3567 0.0626 1.2378 1.2500 6 65.0915 2.6112 1.6290 0.0697 1.4867 1.5000 7 65.1626 2.3755 1.9019 0.0770 1.7366 1.7500 8 65.2725 1.9633 2.1747 0.0827 1.9929 2.0000 Table 4.1: Range-Flow Results for Translat ion M o t i o n to the left-rear Step Angle [deg] Length [cm] # M e a n a M e a n a Encoder C o m m a n d •1 48.8813 7.4717 0.2367 0.0334 0.2230 0.2.500 2 46.2486 4.5264 0.4926 0.0461 0.4750 0.5000 3 46.0249 3.3112 0.7404 0.0518 0.7146 0.7500 4 45.4467 2.6273 0.9997 0.0507 0.9729 1.0000 5 45.6313 2.6246 1.2471 0.0606 1.2161 1.2500 6 45.4051 2.5060 1.5086 0.0780 1.4705 1.5000 7 45.4808 2.4960 1.7614 0.0889 1.7192 1.7500 8 45.2843 2.3798 2.0266 0.0937 1.9736 2.0000 Table 4.2: Range-Flow Results for Translat ion M o t i o n to the rear-right size. T h e first comparison is intended to quantify the influence of mask size and image size. The search area is kept constant relative to the image size, i.e., the search area is 7 x 7 for the 120 x 160, 13 x 13 for the 320 x 240 and 25 x 25 for the 640 x 480 image. The mask size varies for all images over 3 x 3 , 5 x 5 and 7 x 7 . For easier comparison all results are combined into a table per t ime step and time step 1, 4 and 8 are selected and shown in Table 4.3, Table 4.4 and Table 4.5, respectively. T h e selection of a mask size, image size and search area makes various trade-offs. T h e larger the mask size the more distinct is the peak in the correlation window, 59 Image Size M a s k Size Angle [deg] Length [cm] Mean a M e a n a 160, 120 3, 3 65.3530 20.5544 0.2816 0.1378 5, 5 64.2933 20.3311 0.2615 0.0611 7, 7 64.1025 19.6148 0.2630 0.0863 320, 240 3, 3 69.3456 21.9100 0.3711 0.2286 5, 5 64.9249 14.4413 0.2741 0.0959 7, 7 64.8249 13.1080 0.2761 0.0789 640, 480 3, 3 65.7337 15.5447 0.3072 0.1716 5, 5 63.1974 9.7911 0.2469 0.0460 7, 7 63.4371 9.8475 0.2446 0.0327 Table 4.3: Influence of M a s k and Image Size on F low Result The search area is kept constant relative to the image size. T h e commanded motion is 0.25cm while the encoder value reads 0.2355cm (Step # 1 of Table 4.1). but averaging wi th in a larger mask may introduce a larger error since the apparent motion in the image over the mask may change. Addi t iona l ly , the larger the search area the more likely is the inclusion of the correct flow location but the greater the chance of a mismatch. Furthermore, a larger image improves the localization of the peak but it also requires larger masks to make the correlation result distinct. It is also of interest to look at the stabil i ty of the correlation peak wi th respect to the size of the search area. The more distinct the intensity in the correlation window is over the search area, the more unique the peak in the correlation surface w i l l be. The time complexity of the flow computat ion is 0(sizemaskSizesearch,sizeirnage) assuming the number of tracked pixels is proport ional to the size of the image. For example, the difference between the largest image wi th the largest mask size and the smallest image wi th the smallest mask size is 7 x 7 x 25 x 25 x 640 x 480 n n n ~ 300 3 x 3 x 7 x 7 x 120 x 160 Therefore, it is crucial to be conservative in chosen mask size, search area size and 60 Image Size M a s k Size Angle [deg] Length [cm] M e a n a Mean a 160, 120 3, 3 62.8964 11.9992 0.9983 0.1387 5, 5 62.6035 10.9469 0.9989 0.1115 7, 7 62.5492 10.8381 1.0017 0.1182 320, 240 3, 3 63.4279 10.9815 1.0120 0.1350 5, 5 63.3405 5.2712 1.0713 0.0874 7, 7 63.0772 4.7647 1.0966 0.0946 640, 480 3, 3 62.1532 5.3999 0.9891 0.0714 5, 5 62.2516 3.1899 1.0076 0.0514 7, 7 62.1506 2.8523 1.0060 0.0488 Table 4.4: Influence of M a s k and Image Size on F low Result T h e search area is kept constant relative to the image size. T h e commanded motion is 1.0000cm while the encoder value reads 0.9911cm (Step # 4 of Table 4.1). Image Size M a s k Size Angle [deg] Length [cm] Mean a Mean a 160, 120 3, 3 63.9707 7.5256 1.9500 0.1391 5, 5 62.7593 7.8290 1.9893 0.1479 7, 7 ' 62.7484 7.8810 1.9908 0.1429 320, 240 3, 3 61.3715 3.7981 1.9008 0.1522 5, 5 63.3863 3.1567 2.1400 0.1286 7, 7 63.1105 3.2557 2.1945 0.1330 640, 480 3, 3 65.4546 2.2675 1.9501 0.1205 5, 5 62.4504 1.9104 2.0245 0.0673 7, 7 62.1912 1.9717 2.0256 0.0793 Table 4.5: Influence of M a s k and Image Size on F low Result The search area is kept constant relative to the image size. T h e commanded motion is 2.0000cm while the encoder value reads 1.9929cm (Step # 8 of Table 4.1). 61 image size. Also, note, that standard methods of speeding up the calculation by running sums will be less effective here since the evaluated pixels do not necessarily form a closed region due to invalid stereo pixels. 4.1.4 Examples of Object Deformation Figure 4.4 shows the range-flow and vertex displacement estimated for a deforming object, a plush toy tiger. The intensity coding of the flow result shows the result as quantitatively correct. 4.2 Estimation of Vertex Displacement The final step in the deformation measurement is the estimation of vertex displace-ments given the range-flow of surface points. This is related to the mask-size issue in the range-flow calculation itself. The best resolution would be achieved by select-ing the nearest flow vector to a vertex but it would be extremely sensitive to noise. Selecting a larger aperture allows one to filter noisy individual flow estimates but reduces the resolution. One may fit a local model of the displacement, in particular, the underlying linear shape functions of the neighboring triangular patches. This would couple the estimation since the one-ring neighborhood2 of a vertex has to be taken into account. Another choice is to fit a locally rigid model to the motion. This would require fitting a rotation matrix which is computationally expensive. The rotation model is also in contrast to the linear modeling approach of this thesis, i.e., the linear shape functions for displacements. Another choice is an affine motion model. The approach followed in this thesis is one of outlier rejection followed by 2The one-ring neighborhood of a vertex are all the vertices connected to this vertex by an edge. 62 weighted averaging. In the outlier rejection stage, flow vectors associated with a vertex are filtered. All flow vectors which have the start point in the undeformed object state associated with triangles at the vertex are fed into the filter. The filtering is accomplished by finding the median of the flow values in Cartesian space. The residual of individual flows is calculated and sorted. The median of the residuals is taken as an outlier rejection threshold. As a result half minus one flow vectors are discarded. The remaining flow vectors are weighed by 1/(1 + d) based on the distance d from the vertex and averaged. In the following, the nearest neighbor approach is compared with the robust weighted averaging approach of this thesis. 4.2.1 D i s p l a c e m e n t R e s u l t s Table 4.6 summarizes the results for vertex displacement estimation, again for the same translation test as in Section 4.1.3 and for a rotation test conducted in the same way. The rotation in each step is 2.0° to a total of 16.0°. The vertex displacements shown are estimated by the filtered weighted averaging method. The error reflects both the error in the vertex displacement as well as an error in the initial registration of the surface. The error in the orientation is also due to the missing sub-pixel optical flow estimation. This can be seen from the improvement in orientation as the motion vectors become larger (especially in the orientation test). The following comparison is between four different approaches; estimating vertex displacement by the nearest flow vector to the vertex (method NF), the median of flow vectors within the triangles joining at a vertex (method MF), the weighed 1/(1 + d) average of flow vectors after median filtering (method MFWA) 63 (a) 2D Optical Flow (b) 3D Range-Flow (c) Displacement Figure 4.4: Deformation Measurement The 3D range-flow shows the dis t r ibut ion of the estimated flow vectors w i th brighter range-flow corresponding to larger motions. T h e vertex displacement in 4.4(c) also shows larger movement being brighter. D a r k surface is either stationary or no estimate is available. Translation in X Translation in Y Rotation around Z Angle [deg] Angle [deg] Angle [deg] Length [ cm] Error a Error a Error a Mean a Command 10.9142 23.2285 7.4402 15.0430 23.5574 47.8135 0.2998 0.0518 0.3053 8.6968 20.3784 5.8115 11.7741 24.8044 58.0988 0.5928 0.0884 0.6089 7.9037 18.9947 5.2960 10.6937 19.1940 42.1753 0.8739 0.1468 0.9052 8.0567 21.4754 5.0309 10.2249 17.5430 41.4193 1.1370 0.2130 1.2070 8.3435 23.5471 4.8527 9.9227 13.7172 37.8691 1.4744 0.2350 1.5336 8.0368 24.2161 4.8800 9.9940 12.2357 36.0586 1.7621 0.3434 1.8643 8.1894 26.5517 4.9821 10.1555 10.0026 32.9437 1.9929 0.4729 2.1297 7.9706 28.0945 4.9916 10.1864 8.2385 31.2005 2.1469 0.6727 2.3233 Table 4.6: Vertex Displacement Results for a Translat ion and Rota t ion T h e results are for the test-object in Figure 4.3. Vertex displacements are the median-filtered weighted average 3-D flow wi th in the triangles jo in ing at a vertex (see text for details). T h e displacements are transformed into object coordinates. The errors are the combined error over al l vertices and are a result of both registration error and flow estimation error. 64 Step Method Orientation [deg] Magnitude [cm] Error a Mean a 1 N F 14.6050 18.8691 0.2388 0.0384 M F 3.5504 4.4913 0.2373 0.0236 M F W A 3.0627 3.9388 0.2372 0.0191 M F W A S 3.1004 3.9460 0.2373 0.0191 4 N F 4.6907 6.6571 •1.0030 0.0572 M F 1.8624 2.3026 . 1.0141 0.0245 M F W A 1.7360 2.1340 1.0107 0.0192 M F W A S 1.7168 2.1106 1.0107 0.0192 8 N F 2.6801 3.0260 2.0122 0.0737 M F 1.9392 2.6460 2.0325 0.0528 M F W A 1.8062 2.4697 2.0299 0.0484 M F W A S 1.8042 2.4669 2.0301 0.0488 Table 4.7: Compar ison of Estimates of Vertex Displacements for a Translat ion The results are for the test-object in Figure 4.3. Methods are nearest flow vector to the vertex (method N F ) , the median of flow vectors wi th in the triangles jo in ing at a vertex (method M F ) , the weighed 1/(1 + d) average of flow vectors after median filtering (method M F W A ) and median filtering followed by the weighed squared 1/(1 + d)2 average of flow vectors (method M F W A D ) . See the text for an explanation of the methods. and median filtering followed by the weighed squared 1/(1 + d)2 average of flow vectors (method M F W A D ) . Table 4.7 compares these methods for the translation in the y-direction (to the rear-right in the image) extensively used before. The vertex displacements are estimated from range-flow calculated in a 640 x 480 image wi th a mask of 5 x 5 and a search area of 25 x 25. Table 4.8 compares the methods NF, MF, MFWA and MFWAD for a rotation around the z-axis described above. T h e vertex displacements are again estimated from range-flow calculated in a 640 x 480 image wi th a mask of 5 x 5 and a search area of 25 x 25. T h e rotation penalizes averaging more since the variat ion in flow vectors is much higher over the set of triangles jo ining at a vertex. 65 Step Method Orientation [deg] Magnitude [cm] Error a Mean a Command 1 NF 24.2158 28.0852 0.2923 0.0695 0.2826 MF 20.8808 22.5736 0.2939 0.0299 0.2732 MFWA 20.7130 22.1168 0.2939 0.0281 0.2736 MFWAS 20.5988 22.0136 0.2939 0.0282 0.2738 4 NF 15.4796 16.1011 1.1682 0.1535 1.0821 MF 13.3114 18.3754 1.1531 0.0911 1.0928 MFWA 13.8109 18.7523 1.1531 0.0861 1.0932 MFWAS 13.7860 18.6873 1.1531 0.0839 1.0953 8 NF 7.4156 9.8216 2.3229 0.2536 2.2000 MF 5.5896 9.6690 2.2861 0.1968 2.1748 MFWA 5.4954 9.5279 2.2861 0.1771 2.1867 MFWAS 5.4283 9.3702 2.2861 0.1710 2.1893 Table 4.8: Comparison of Estimates of Vertex Displacements for a Rotation The results are for the test-object in Figure 4.3. Methods are nearest flow vector to the vertex (method NF) , the average of flow vectors within the triangles joining at a vertex (method A V G ) , the filtered average of flow vectors (method FAVG) and the filtered weighted average of flow vectors (method F W A V G ) . See the text for an explanation of the methods. 4.3 Summary This chapter details the vertex displacement measurement approach of this thesis. The choices in implementing the method are discussed. The method implemented is a robust range-flow measurement technique based on simultaneous optical flow and stereo-range data from trinocular imagery. Robustness is achieved by addi-tional imagery. The approach can be characterized as a purely local bottom-up approach. The measurement of the vertex displacements is independent of the de-formation model. This independence is convenient because it allows a comparison of deformation models independent of displacement measurement. Even more im-portantly, since the displacement of each vertex is estimated independently, partial measurements of deformation is trivial. Partial measurement is paramount for all but the simplest deformable objects. Furthermore, covering the surface of a deform-ing object nearly always requires multiple views, each view covering different sets of 66 vertices. Different approaches to deformation measurement are possible. In fact, many researchers have suggested the use of deformation models to regularize various com-puter vision tasks (e.g., [127, 98, 21, 22]). A n approach employing the deformation model as a constraint in the vertex estimation is possible (e.g., [96]). 67 C h a p t e r 5 Estimation of Discrete Green's Functions Est imat ion of discrete Green's functions based on Equa t ion 2.1 can be viewed as a linear estimation problem. T h i s view disregards the underlying structure of the Green's functions mat r ix and is explored in this chapter. T h i s is in contrast to Chapter 6 where the boundary element method derivation of the Green's functions mat r ix of Section 2.3.2 is exploited to solve an inverse problem. Th i s chapter can be understood as either the final step in the modeling process or as precursor to an inverse problem solution (see Chapter 6). Firs t , the local behavior of an elastostatic object around the contact point is discussed. The local behavior is dominated by the compliant motion of the object because of the probe. Typica l ly , the local deformation is significant but the deforma-t ion decays quickly w i t h the distance from the contact point. (Global deformation does not decay quickly w i th large scale motion of the object, e.g., due to articula-t ion of the object as w i t h the example in Section 5.5.2). T h e overall response of the 68 object to the contact force is more varied than the local response across different objects. The topology of the object as well as (in)homogeneity of the material of the object have great influence. The measurement set-up briefly summarized in Section 5.1 is the starting point for the estimation procedure. The local behavior estimation is described in detail in Section 5.2, while the global behavior estimation is discussed in Section 5.3. This chapter concludes with two examples: estimation of the discrete Green's func-tions matrix of a plush toy in Section 5.5.1 and an anatomic soft-tissue human wrist model1 in Section 5.5.2. 5.1 Set-Up of Estimation Equation 2.1 combines the prescribed displacement and traction vectors into block vector v and the complementary unknown displacement and traction vectors into block vector v. Vector v consists of tractions for the free surface and displacements for the fixed surface (see Section. 2.1 for details). The displacements for the fixed surface are zero and all tractions for the free surface are zero except for the contact area. This structure of the problem is important and is exploited later on. The set-up of the measurement and estimation tasks must match the above mathematical description. The model is fixed to a support mounted on the Test Station. The support surface of the object corresponds to the fixed surface I V The rest of the surface is free and corresponds to the free surface I V The Contact Measurement System xThe anatomic soft-tissue human wrist model has a bone structure made of hard plas-tic with a cancellous inner core. The bone structure is surrounded by soft foam with a plastic skin layer. The object is made by Sawbones, Pacific Research Laboratories, Inc. (http://www. sawbones, com). 69 (CMS) probes surface locations on the free surface corresponding to vertices of the object's mesh at the desired reconstruction resolution (see Section 3.3). If the contact area of the probe of the CMS is smaller than the size of the triangles in the mesh, it is reasonable to assume that a single vertex is affected by the probe. This is true for all the models in this thesis. As a result the structure of the block vector v is all zero, except for a single traction vector at the contact point at vertex k. This leads to the structure illustrated in Equation 5.1. never excited unobservable 0 Pfc 0 ui Ufc not observed (5.1) The key behavior of the model is characterized by the displacement response of the m vertices of the free surface F\ due to forces applied to it. This displacement corresponds to the measurable portion of the v block vector. The rows of v which correspond to vertices of the fixed surface To are the tractions of the support surface. These tractions are not measured in ACME. Therefore, the corresponding rows in the Green's functions matrix H cannot be recovered. These same vertices correspond to columns in S which can not be actively excited since their displacement is fixed to zero. The position of the contact probe measures the displacement of the contacted vertex uk, while the force-torque sensor measures the contact force. The force at a probed vertex has to be converted into surface traction for the estimation. The observed force has to be the surface integral of tractions over the affected area. 70 The surface is represented as a tr iangular mesh and the shape functions for the tractions are linear. The components of the force vector fk are the sum over al l integrated tract ion shape-functions on elements connected to the probed vertex, i.e., fk = Ylj \ {Pk + Pk ~^Pk) 1^1- T h e tract ion over triangle j is a-linear combination of tractions plj .. .plj at its vertices 1 . . . 3; while \G\ is the Jacobian from local to global coordinates (see A p p e n d i x A for more details). T h e estimation problem is: given a set of block vectors v and a set of corre-sponding block vectors v find the m x m observable submatr ix of 3 . T h i s problem can be broken up into m separate estimation problems for each of the m columns due the specific structure of Equa t ion 5.1. Each of these m estimation problems per column consists of finding a 3 x 3 block element of matr ix S given a set of observed displacement vectors u , and the corresponding traction vectors p^- T h e problem when i — k is the estimation of the compliance of the object at location k. This is discussed next in Section 5.2. 5.2 Local Behavior T h e local behavior or the self-compliance of a vertex is of special importance. T h e local stiffness is H^ 1 , the inverse of the self-compliance matr ix . T h i s local stiffness is employed in haptic simulation to map displacements of the tool t ip into a force response [58, 57]. The correct force response can be easily calculated at the high haptic rendering rates required (typically « 1000Hz). A l l that is required is to mul t ip ly the local stiffness at the contact point(s) wi th the displacement vector(s) due to the tool . For a single contact vertex, the equation is 71 For this approach, it is assumed that the self-compliance elements of the Green's functions matr ix are invertible. T h e compliance mat r ix w i l l be invertible as long as the response is not r igid or infinitely stiff in any direction. T h i s local model is exact w i th respect to the global model and is automatical ly synchronized wi th the overall deformation response. T h e global deformation response can be updated independently at any time (for more details and a comparison to other approaches see [57]). The self-compliance estimation in the formulation of Section 5.2.1 is a linear least-squares problem given a set of displacements at the contact point and the corresponding applied traction. 5.2.1 Least Squares A vertex location is probed I times. These I probes are executed such that the I directions of the vertex displacement form a spanning set in three dimensions. The force response of the object is recorded giving I t raction vectors. T h e standard least squares problem to be solved for HjTfc is then PfcPfc • • • Pfc In the standard notation for least squares, this becomes Ax = b. T h e normal equations of the problem are then A T A x = A T b . T h e solution of this system can be found by singular value decomposition of the measurement design matr ix A = USV T . Here, the 3 x I ma t r ix U and the 3 x 3 mat r ix V have orthonormal column vectors [ui, • • •, uj] and [vi, • • •, V3 ] , respectively. The 3 x 3 diagonal mat r ix X contains the ordered singular values o\ > 02 > 0-3. The normal equation may then be expressed as V E 2 V T a ; = VSU T 6. 72 r-^T _ 1 2 "fcfc — UfcUfc Reordering this equation leads to where the diagonal matr ix X - 1 contains the ordered singular values l/o~i < I / 0 2 < I / 1 7 3 . Least squares problems wi th a measurement design mat r ix A w i th a poor condit ion number cannot be solved wi th the standard procedure above. Condi t ion numbers of the design mat r ix A which occur based on measured traction vectors are considered next. 5.2.2 Measurement Design Matrix In theory, in a numerical solution of a discrete elastostatic boundary value prob-lem, rank deficient local stiffness matrices (HT^1) are rarely encountered if at al l . For isotropic materials and most geometries local stiffness matrices are diagonally dominant, i.e., if a force is exerted on an object, the objects deforms locally in the direction of the applied force. However, measuring the response of an object is dif-ferent since the design mat r ix A = [pjfcPfc • ' ' PL] n a s to be well conditioned. The tract ion measurements plk are pol luted wi th measurement noise and modeling error. A s mentioned previously, it is ensured that the probe directions span 3-space in setting up the least squares problem. T h i s guarantees that the local compliance model is excited in all directions. Nevertheless, the measurement of the traction response of the object may not span 3-space. One probe direction is normal to the mesh surface of the object, while addit ional directions are offset from the normal. T h i s offset is l imited by the fact that the contact between the probe and the model is not allowed to slip for a measurement of the static force response at a given contact point. In A C M E a probe t ip w i th high friction (grinding stone tool tip) is used in order to use an offset as large as possible. B u t this setup necessarily leads to greater 73 (a) Vertex 38 (b) Vertex 55 Figure 5.1: Probe T i p at Vertex Loca t ion Vertex Ol 0-2 0-3 0-1/0-2 cri/o-3 38 1.6426 0.4317 0.3940 3.8047 4.1696 55 2.5786 0.2638 0.1762 9.7753 14.6327 Vertex VVan \]Var<2 \fVar$ \\pT<aT TTT\\ I K ^ k k - U 112 38 2.3513 2.0164 1.6076 0.4094 3.9993 55 4.1377 3.8120 3.8823 0.6973 6.8947 Table 5.1: Loca l Force Response Locat ion of vertices is shown in Figure 5.1. Response at vertex 38 is well conditioned while the response at 55 is poor ly condit ioned. certainty in measurement of the force response in the direction of the surface normal than in any off-normal direction. For the remainder of this section measurements of the tract ion response of a plush toy tiger (see Figure 3.8) at two locations are considered to il lustrate the discussion. One location is vertex 38 of the base mesh on the back of the tiger as shown in Figure 5.1(a), while the other location is vertex 55 of the base mesh on the head of the tiger as shown in Figure 5.1(b). The tract ion response on the location on the back spans 3-space well. Th i s is apparent from the ratio of singular values listed in Table 5.1. T h e tract ion response at the measurement location on the head 71 of the tiger is very different. Table 5.1 shows that the ratio of singular values for vertex 55 is larger than for vertex 38. T h e condit ion number at vertex 55 is 14.6327 compared to 4.1696 at vertex 38. A- condition number of 14.6327 does not cause any numerical difficulties but here it signals a measurement problem. T h e larger condit ion number is caused by a traction response wi th a dominant direction. T h e dominant direction is due to the head, of the tiger apparently hinging around its neck. The dominant motion is around the axes of this apparent hinge. T h e directional change of the probing does not change the direction of the force response significantly enough. Here, not significantly enough means that the measured force response of the normal direction is small compared to the resolution of the force sensor (see Section 3.1.1) and the probe localization uncertainty. T h e probe localization error contributes to the force error since the normal direction of the probe is changing over the surface. The singular values <7i of the measurement design mat r ix determine how sensitive the least squares solution is to errors in the associated right singular vectors V i . Th i s can be seen by rewri t ing Equa t ion 5.2.1 as . i = i ' Assuming the error is equally distr ibuted in the components of the right-hand-side measurement vector 6, the factor l/<7i acts as an amplification of the error. A design matr ix wi th singular values differing by more than one order of magnitude, as in the case of vertex 55 on the head of the plush toy tiger, causes the signal to noise ratio to decrease by more than one order of magnitude. In Section 5.2.3 regularization is used to deal wi th poorly conditioned traction responses. 75 5.2.3 Regularization Regularizat ion transforms a parameter estimation problem into a constrained opti-mizat ion problem. T h e goal of regularization is to find an estimate which minimizes the residual of a parameter fit (just as in least squares) but constrained by a dis-tance measure from an a priori known solution. Th i s section discusses the trade-off between choosing a solution wi th low residual and choosing a solution which fits the a priori constraint in the context of estimating the local compliance mat r ix S^ f c. The direct regularization methods of truncated singular value decomposition ( T S V D ) and damped singular value decomposition ( D S V D ) are compared wi th the non-regularized solution. The T S V D (truncated at t) solves the problem min||x||2 subject to min | |A t x — t»|j2-T h i s solution is calculated based on the S V D (Equat ion 5.2.2) as \ - r u b where fj is a filter factor and f; = 1 V i < t and f; = 0 V i > t. T h e D S V D employs filter factors based on a regularization threshold A as f. = _2i_-Oi + A Beside comparing the residual of fit and constraint error, another method for choosing the regularization threshold A is the generalized cross validation ( G C V ) . Below, the example of vertex location 55 shown in Figure 5.1(a) is evaluated. Some of the results in this section are obtained wi th the Ma t l ab package Regularization Tools published by P . C . Hansen [48]. Consider the variance in the elements of the estimated compliance matr ix . The variance Varj in the rows j of the transposed compliance matr ix EjTfc are [102, 76 pp. 676-681] Var3 = £ & i= i ^ T h e variance wi th in a row of the mat r ix H^, is identical because of the fact that estimating & T k is three similar least squares problems wi th three parameters each. T h e reciprocal of a singular value 1/CTJ scales the contribution of the correspond-ing basis vector to the variance of a parameter as well as to the estimate of the parameter. These variances describe the uncertainty in the compliance mat r ix in centimetres[cm] per Newton[N}/centimeter2 [cm2]. The variances of the three rows are listed in Table 5.2 for the i l lustrative example of vertex 55 (Figure 5.1(a)). T h e T S V D w i l l weigh the different right-hand-side singular vectors v either by the inverse of the singular value or zero in Equa t ion 5.2.3. If the inverse of a singular values is large, the uncertainty in the associated column of right-hand-side singular vector is large. Th i s leads to large variances V'arj. T h e determination of large corresponds to the choice of the regularization parameter. In the case of estimating the 3x3 local compliance matr ix , the choice can have three possible outcomes, i = l , i = 2 o r i = 3 (no regularization). The three possible outcomes are listed in Table 5.2. A s expected, the norm of the solution decreases wi th the increase in regularization. T h e regularization parameter can be chosen from Table 5.2 given an expected error (residual of fit). In the estimation problem at hand, one may be able to deduce an expected error by comparing all the ind iv idua l estimation problems. However, the error is dependent on location and on the particulars of a measurement and set-t ing a general threshold may be difficult. One approach to deal w i th this difficulty is to estimate the threshold from the data itself. The generalized cross validation does exactly that. The idea is to find a regularized solution which is capable of 77 Vertex 55 t = 3 t = 2 t = l \/Var2 \JVar3 4.1377 4.1377 4.1377 3.8120 3.8120 0.6516 3.8823 0.5722 0.2809 | | P r S j k - l F | | 2 0.6973 1.3721 1.6926 Il3fc/cll2 6.8947 3.8992 2.5143 ox = 2.5786 a2 = 0.2638 a3 = 0.1762 Table 5.2: Regular izat ion wi th the T S V D The table shows a comparison of the variances Vari of the rows of the local compliance matr ix SjQ., residual of fit | | P r S ^ f c — U r | | 2 and the norm of the solution ||Sjfc||2. The number of non-zero singular values is t. 0.015 Figure 5.2: G C V for T S V D M i n i m u m of the generalized cross validation ( G C V ) for regularization wi th the truncated singular value decomposition ( T S V D ) . T h e G C V for each row of H^fc with k = 55 is shown separately, as well as their sum (solid here). T h e G C V is discrete at t = 1 and t = 2 (t is the number of non-zero singular values); however function values are connected for clarity. 78 predict ing an arbi trary measurement if this measurement is dropped from the esti-mat ion. Addi t ional ly , the choice of regularization parameter should be independent of a orthogonal transform of the measurement data. Th i s leads to a G C V [48] to be minimized . T h e G C V is calculated by G (traced " P T (E?=i f ^ ) ) 2 T h e G C V for the regularization w i t h the T S V D is evaluated at two locations and the smaller value in the example is at t = 2 basis vectors (see Figure 5.2). A t first, the T S V D may seem like a good fit to the estimation problem at hand. However, the problem w i t h the solution obtained is that the estimated com-pliance mat r ix itself is rank deficient. However, rank deficiency w i l l cause difficulty i n the haptic simulation because of the use of the inverse of the local compliance mat r ix (i.e., the local stiffness matr ix) . T h e local compliance mat r ix w i l l be com-pletely r igid, i.e., non-compliant i n the direction of the zeroed basis vector. The T S V D regularization has a discontinuous filter function, i.e., it is either 1 or 0. The damped singular value decomposition ( D S V D ) replaces these filter factors by the smoother function of Equa t ion 5.2.3. T h e D S V D solution for the example is listed in Table 5.3 for regularization parameters A resulting in the same residual as for the T S V D of Table 5.2. The G C V curve in case of the D S V D is increasing monotonical ly (see F i g -ure 5.3). In this thesis the regularization parameter is chosen heuristically. The heuristic is based on the empir ical observation that the transi t ion between poorly and well-conditioned tract ion responses occurs at about uxjoj, ~ 10. A choice of regularization parameter of A = CTI/10 — 0 3 and A = 0 if 0 1 / I O — 0-3 < 0 suggests itself. Th i s leads to a regularization parameter of A = 0.08164 for the poorly condi-tioned force response at vertex 55 and to A = 0 for the well-conditioned response at 79 Vertex 55 A = 0 A = 0.32165 A = 0.51677 A = 0.08164 y/Vari 4.1377 1.4809 1.0735 2.8340 \JVar2 3.8120 1.7081 1.2802 2.9026 3.8823 1.4000 1.0181 2.6656 0.6973 1.3721 1.6926 0.8489 6.8948 . 2.9310 2.4457 4.8591 0-1 = 2.5786 cr2 = 0.2638 0-3 = 0.1762 Table 5.3: Regular izat ion wi th the D S V D T h e table shows a comparison of the variances Vari of the rows of the local stiffness matr ix S^ f c, residual of fit | | P T S ^ f e — UT||2 and the no rm of the solution ||S f^c||2. 0 0.2 0.4 0.6 X Figure 5.3: G C V for D S V D M i n i m u m (on the left) of the generalized cross val idat ion ( G C V ) for regularization wi th the damped singular value decomposition ( D S V D ) . T h e G C V for each row of E^f, w i th k — 55 is shown separately, as well as their sum (solid line). vertex 38. The heuristic achieves the goal of regularizing poorly conditioned trac-t ion responses but fitting the data opt imal ly for well conditioned traction matrices. A l l the results for the local object behavior in Section 5.5 are calculated wi th this approach. In general, the D S V D solution has the advantage of having full rank as long as the design matr ix A has full rank. T h i s avoids difficulties in haptic simulation of the stiffness matrices. Furthermore, all the left hand singular vectors u; are damped, especially al l off-normal vectors w i th their corresponding smal l singular values. Th i s 80 stabilizes the solution in the face of noise in the off-normal directions. 5.3 Global Behavior The measurement of the global behavior of the model is quite different from the mea-surement of local behavior of Section 5.2 above. T h i s is not obvious at first since the structure of the problem is identical. T h e design mat r ix A of the ind iv idua l least squares problems for estimating ajk when j ^ k is identical to the design mat r ix of the self compliance mat r ix H^fc. However, the right hand side b of the problem, the measurement of the vertex displacement is very different. T h e measured surface dis-placement in case of the local behavior is the result of the kinematics of the contact probe. For the global behavior the displacement is measured by the deformation measurement technique of Chapter 4. Th i s results in the following main differences which need to be addressed: increased noise level, addi t ional noise in the form of outliers, and incomplete measurements. Outl iers can, on occasion, s t i l l be observed in the displacement measurement despite the filtering process described in Section 4.2 (see Section 4.2 for a description of shortcomings of the filtering process). Th i s should not be a surprise since most visual measurement techniques need to deal w i th outlier noise. Methods to deal w i th these outliers are discussed in Section 5.3.2. Incomplete measurements arise from lack of observed object surfaces due to occlusion by the probe, self-occlusion of the object and incomplete coverage of the object surface by the visual sensor. The problem of incomplete measurements manifests itself bo th in a reduction of the number of observations per estimate at a vertex and in a complete lack of observations for some vertices. Fewer measurements at a single vertex can be dealt wi th by increasing regularization wi th the same 81 methods as in Section 5.2.3. In Section 5.3.1 the application of the techniques of Section 5.2.3 to global behavior is discussed. A t vertices where there is a complete lack of observation, there is no choice but to somehow infer the behavior from neighboring vertices or to infer the behavior from the behavior of the vertex itself due to probing at a different but close-by location. Considering neighboring vertices corresponds to interpolation wi th in a column of the block matr ix S, while similar probes correspond to interpolation wi th in a row of the block matr ix H. A l s o note that obtaining more measurements can reduce the number of vertices w i t h complete lack of knowledge but is not able to deal w i t h systematic failures. Systematic failures are the consistent failure of the deformation measurement technique for a given vertex due to lack of visual texture or occlusion, either by the probe or by the object itself for al l observations. G loba l regularization schemes are discussed in Chapter 6 but see also Sec-tion 5.4. T h i s section only considers interpolation schemes wi th in a block column of the Green's functions matr ix S. Here, scattered data reconstruction is used to fill in elements for each column individually. T h e interpolation of missing displacements is achieved by solving Laplace's equation over the set of unestimated vertices. T h i s hole filling approach is described in Section 5.3.3. 5.3.1 Regularization T h e off-diagonal matr ix elements Sjk w i th j ^ k are estimated ut i l iz ing the damped singular value decomposition ( D S V D ) . T h e method is identical to the estimation of the diagonal elements E^. described in Section 5.2.3, however the regularization parameter A is chosen differently. Addi t ional ly , the estimation is not attempted if the set of displacement measurements is not the result of probe directions which 82 span 3-space. Th i s occurs as a result of a reduction in the. number of observations for a given vertex. In this case all deformation measurements at that vertex are disregarded and the hole filling process is left to fill them in . The estimation process is also integrated wi th the outlier rejection process of Section 5.3.2. T h e choice of regularization parameter in the D S V D for vertices w i th sufficient observation is discussed below. T h e noise in the deformation measurement effects the direction of the dis-placement much more than the magnitude of the displacement (as discussed in Section 4.2). Regularizat ion wi th the D S V D is again well suited for this problem. Smal l variations in direction should not be magnified by small singular values of the.design mat r ix A. T h e D S V D dampens these small singular values. The only difference then, compared to the local behavior, is in the choice of regularization parameter. A large regularization parameter filters the singular values more. Since the noise in the visual displacement measurement is larger than the noise in the measurement of the displacement wi th the C M S , a larger regularization parameter employed for the estimation of the global behavior than for the estimation of the local behavior. The heuristic is for the global behavior A = aj/5 — 0 3 or A = 0 if (7 i / 5 — 0 3 < 0 compared to A = eri/10 — 0 3 . However, the amount of regularization may st i l l not be sufficient at some locations w i th large errors. These locations need to be identified and the regularization increased correspondingly. Another heuristic is employed in these cases. A way of looking for the right amount of regularization at these problem location is to make use of the fact that the magnitudes of traction and displacement are well defined. The local behavior defines a ratio between observed m a x i m u m compliance max{||I2/I(pi112) over al l measurements I at vertex k and the norm 83 of the solution | |H^||2 (the largest singular value o f the solution matr ix) . In the case of the il lustrative example of vertex 55 on the head of the tiger, the max imum observed compliance is max( | j u.^ . j 12 /11 112) =3.1172. Th i s gives a ratio of R A T L 0 K = 711 i 1 ! ? ^ ! ! 1 1 2 i i n = L 5 5 8 8 ' max(| |u T F C | |2 / | |maxp ' A . | |2) T h i s ratio is also calculated for a l l other vertices j ^ k in the column, i.e., ratio j = l l H J f c l l 2 / ( m a : r ( l l u j l l 2 / | |p i l | 2 ) ) - T h e solution S j f c is found w i t h the D S V D wi th A = <7i/5 — 03. T h e estimated solution at a vertex j is suspicious if the ratio is higher than for the local compliance. In such a case, a larger regularization parameter may produce a more desirable solution. Therefore, the D S V D is repeated wi th the regularization parameter A = o i / 5 * ration/ratioj — 0-3. T h i s heuristic achieves the goal of rejecting solutions wi th norms larger than expected but allowing arbitrary large norms if very compliant behavior is observed. 5.3.2 Outliers T h e goal of the outlier rejection process is to reject erroneous observations of the deformation. The problem is to select a subset of measurements which does not contain any spoiled measurements. A number of selection methods (or positive-breakdown methods) exist wi th in robust estimation techniques, e.g., random sample consensus ( R A N S A C ) [14], least median of squares ( L M S ) [110] and least t r immed squares (LTS) [110]. R A N S A C selects P number of samples out of M measurements where P is the number of parameters of the model to estimate. Given a pr ior i knowledge of the max imum measurement error, inlier points (i.e., points which are not outliers) are selected for a given estimate. If at least a m i n i m u m number of points are wi th in the expected measurement error, the estimate is accepted. Final ly , a least 84 squares fit to the inlier points only is computed and this is the estimate produced by R A N S A C . R A N S A C is guaranteed to succeed if the m i n i m u m number of points H for acceptance is greater than ( M + P + l ) / 2 . L M S proceeds identically except it does not use an expected error to classify inlier points. Instead all possible P sets are evaluated and the estimate wi th least median of squared residual is selected. A l l points w i th a squared residual smaller than or equal to the least median of squared residual are considered inlier points. T h e final least squares step and the number of inlier points necessary for L M S to succeed are identical to R A N S A C . F ina l ly , L T S is an improvement over L M S . In L M S the evaluation of a t r ia l set P depends on the errors in the data points P, i.e., the median of squared residuals of a data set may have been considerably lower if the in i t ia l estimate would have been improved by a least squares step. Therefore, L T S improves each subset of points H by improving the fit w i th repeated least squares [111]. The least squares iteration has converged when the H set has become stable. L T S finds the H set wi th the min imum least squares residual. T h e L T S process again is robust even if ( M — P — l ) / 2 measurements are spoiled. T h i s thesis employs L T S implemented very s imilar ly to the description in [111]. One modification of the algori thm is in the selection of the in i t i a l P-set. P = 3 for the estimation of S j f c since each column is represented by the same parameter (as discussed in Section 5.2.3). There is no point to selecting an in i t i a l subset P which is rank-deficient. Therefore, bins are employed to group tract ion measurements to-gether which are obtained wi th the same orientation of the displacement vector. B i n n i n g has been suggested, e.g., in [134]. The other modification is the use of regularization by D S V D inside the least squares iteration steps of L T S . 85 5.3.3 Hole Filling T h e hole filling process interpolates data wi th in the block column k of the Green's functions matr ix Sjk. T h i s is possible since the observable part of the Green's functions mat r ix contains compliance matrices only (see Equa t ion 5.1). A column of the compliance mat r ix can be viewed as the displacement due to a unit t raction applied along the axis of the object coordinate system. A s a consequence, the three columns of the observable part of Sjk for j = 1,..., n a l l represent compatible values. If the displacement field is smooth, the columns are smooth. Th i s smoothness assumption is typical ly violated at the contact point. Here the displacement field has a cusp. Note that K e l v i n ' s fundamental solution for infinite linear-elastic isotropic homogeneous materials decays smoothly wi th 1/r where r is the distance from the contact point. T h e hole filling process is never necessary at the contact point since A C M E ' s contact measurement system guarantees that observations are made. Smooth hole filling can be achieved by solving Laplace 's equation over ver-tices j w i th unestimated Sjk subject to estimated values S;fc at vertices I and S t t,fc = 0 for vertices w which are fixed. Th i s may be computed efficiently in several ways, however time-stepping the diffusion equation sufficiently long is effective here. T h e discrete Laplac ian L at time t and vertex j is calculated using data from the vertex's one-ring neighborhood 2 , Nj where e3i is the edge joining vertex j and i. Exp l i c i t forward Euler discretization is computed as 2The one-ring neighborhood of a vertex is all the vertices joined to that vertex by an Bjle(t + At) = Ejk(t) + At L(S]k{t)) edge. 86 subject to a time-step restriction on A t [32]. T h e process is started wi th the un-estimated Sjk(t = 0) — 0. Constraints can be introduced wi th a weighting factor 0 < 8 < 1 and the Euler time step becomes S3k(t + At) = 8 3jk(0) + (l-8)(3jk(t) + AtL(Ejk(t))). T h i s simple process is quite effective for fill ing a l imited number of holes, even though explicit methods may require a large number of smal l t ime steps. Solving Laplace's equation has also been applied successfully to smooth displacement fields and in geometric mesh fairing [32]. The result of this interpolation process is shown in Figure 5.7. 5 . 4 M e s h R e f i n e m e n t The subdivis ion mesh hierarchy can be exploited to improve rendering quality and reduce measurement and estimation time. T h e object can be contacted at vertices of a coarse resolution mesh but the global behavior may be estimated at any mesh of level / given sufficient spatial resolution of the visual measurement of the defor-mation. A s a result the estimated Green's functions matr ix S at the resolution I of the deformation measurements lacks columns which correspond to odd vertices inserted at that subdivision level. In order to fill in these columns, interpolation wi th in a row is necessary. In this thesis, the missing columns are interpolated for rendering purposes only. T h e interpolation is linear except for the one-ring neighborhood of the local response. Here a special shifting is done to preserve the diagonally dominant struc-ture of the overall Green's functions matr ix and prevent linear dependence columns of S (see Equa t ion 5.2). It is stressed that global inverse problem solutions are by 87 far preferable to such a crude interpolation. However, for rendering the estimated Green's function for the examples presented in Section 5.5, the described approach produced good results. Note as well that this scheme may be ut i l ized for adaptive sampling of the surface. Assume the mesh is probed at mesh resolution Z and the discrete Green's functions matr ix is estimated. Tak ing the columns of the even ver-tices, i.e., corresponding to vertices of mesh level I — 1, the response at odd vertices may be found by interpolation. If the measured response at odd vertices is close to the response at these vertices found by interpolation, sufficient spatial resolution is achieved. (5.2) T h e interpolation is executed on subdivision level Z. Level Z has been sub-divided once from the level Z — 1 and all observable block columns of S are available at level Z — 1. Then , the jth odd vertex on level Z has a response S y inferred if both even vertices k and m of its parent edge have responses 3^ and 3;m, respectively (see Equa t ion 5.2). If so, S j j is l inearly interpolated from these two responses. T h e local behavior, 3jj and S n j - when vertex n is a one-ring neighbor of j, is handled 88 differently. Th i s local behavior is dominated by a peak at vertex location j and linear interpolation would introduce a large error. Instead, these local responses are computed as the weighted average of parent responses which are diagonally shifted (see Equa t ion 5.2, Snj is handled exactly as Sjj)- For example, Sjj is l inearly interpolated from H f^c and S m m and not from Skj nor 3 m j as for far field responses. 5 . 5 R e s u l t s In this section, estimation results for two objects are presented, one a plush toy tiger shown in Figure 3.8, the other a anatomic soft tissue human wrist model . Presen-tat ion of the results on paper is somewhat l imi t ing . V i s u a l and haptic interaction wi th the models provides a far clearer picture of the results 3 . Quantit ies which can be presented here include the residual of fit for the local behavior estimation, i.e., for the diagonal of the Green's functions matr ix . The dis t r ibut ion of the max imum eigenvalue of the compliances (i.e., the diagonal of Green's functions matr ix) are shown to give a relative sense of compliance of the model. For the global behavior, examples of the columns of Green's functions matrices, are shown graphically. Th i s is meaningful since these columns present the response of the model to an applied unit t raction in a given direction at a vertex location. Furthermore, the effect of the regularization process, the outlier rejection by L T S and the hole fi l l ing by Lap la -cian smoothing are presented for the plush toy. Final ly , the effect of different mesh resolutions is presented. 3 Haptic interaction with the plush toy. tiger has been successfully demonstrated to au-diences at several scientific conferences (e.g., 11th IRIS Precarn Annual Conference, June 2001, Ottawa, Canada; 8th Int. Conf. of Computer Vision, July 2001, Vancouver, Canada) 89 '0 200 Time [10ms] 400 0. 0 200 Time [10ms] 400 .(a) Probe at Vertex 38 of Figure 5.1(a) (b) Probe at Vertex 55 of Figure 5.1(b) Figure 5.4: Measured Magni tude of Force and Displacement 5.5.1 Plush Toy Firs t , results for the local behavior of the model are presented. Figure 5.4 shows the magnitude of force and displacement at the probe t ip dur ing deformation. The tip executes a spatial linear motion while in contact wi th the object. A n ideal linear elastostatic object sensed wi th a perfect sensor would cause the force profile to be a mult iple of the displacement. T h e plots in Figure 5:4 show approximately such a behavior, however there are also differences. The force does not remain constant after the probe stops its motion but decays exponentially to a lower value. Th i s decay corresponds to viscoelastic object behavior. T h e plush toy's fill ing material resettles to accommodate the deformation to some degree. T h e overall magnitude of the effect is dominated by the magnitude due to linear elastic behavior. T h e force profile is not zero at the start since before the recording of the profile, the probe approaches the surface unt i l the force reaches a threshold. T h e displacement 90 is offset by a corresponding amount not shown in the graphs. D u r i n g measurement the force profile is l inearly interpolated to zero, as well as the displacement. The force sensor causes some high frequency noise in the force profile. The result of the estimation of the discrete Green's functions is shown in Figure 5.5. Figure 5.5(a) shows the max imum singular value of the diagonal of the Green's functions matr ix . T h e diagonal describes the compliance of the solid. The figure shows that the tiger is maximal ly compliant at the rear of the head. In general the head is more compliant than the back. Th i s corresponds well to the inspection of the plush toy by hand. Figure 5.5(b) shows the residual of fit | | P T E j T f c - X J X | | 2 for the local compliance matrices. T h e residual is uniformly low on the back of the tiger. T h e estimates for the head of tiger have increasingly larger residual further away from the neck. Furthermore, some problem locations are on the ears of the tiger and on the ta i l . T h e dis tr ibut ion of the residual is clearly different from the dis tr ibut ion of the largest singular value. T h i s indicates that the residual depends on the location over the surface and is not just a percentage of the m a x i m u m compliance. T h e results for the global behavior are i l lustrated wi th the two Green's func-tions block columns corresponding to vertex 38 and 55. T h e location of vertex 38 and 55 are shown in Figure 5.1(a) and (b), respectively. Figure 5.6 shows the ef-fect of regularization on the y-column of the Green's functions matr ix . The x- and z-column show similar behavior and are not shown. T h e behavior without regu-larization shows large Green's functions magnitudes wi th an uneven dis tr ibut ion. Regularizat ion ensures that noise in the observations do not cause solutions wi th excessive norms. Th i s can be observed in the reduction of the Green's function y-vectors. The estimates for vertex 38 are satisfactory without regularization, but improve wi th regularization. Notice that on the head of the tiger the movement 91 /A (a) Solution Norm | | E £ f c | | 2 (b) Residual of Fi t | | P T 3 j [ f e - U T | | 2 Figure 5.5: Loca l Compliance The figure is vertex colored, red indicates a large value while blue indicates the smallest value. Vertices not probed due to reachability constraints result in S?k not being estimated and are shown in black. The result shown is obtained wi th a max imum of eight probes per vertex location and observations from two different viewpoints. 92 (c) Probe at Vertex 38 of Fig- (d) Probe at Vertex 55 of ure 5.1(a) Figure 5.1(b) Figure 5.6: Regularizat ion of a Green's Functions Es t ima t ion Est imated y-columns of the discrete Green's function without regularization are shown in 5.6(a) and 5.6(b). T h e result wi th regularization (as in Section 5.3.1) is shown in 5.6(c) and 5.6(d) T h e y-column of the observable por t ion of a block column of the discrete Green's function corresponds to the displacements of the free vertices due to a unit t raction in the y-direction (see Equa t ion 5.1). The displacements are color-coded where red indicates a large value while blue indicates the smallest value. 93 gets damped while it persists around the probe location. The estimates for vertex 55 are in clear need of regularization as discussed in Section 5.2. Observe how there are large displacements in stray directions on the head of the tiger before regular-ization. Aga in , regularization dampens these directions but preserves the observed displacements near the probed vertex. The behavior for vertex 38 is more local than for vertex 55. Figure 5.7 shows the effect of Laplac ian smoothing described in Section 5.3.3. T h e black areas corresponding to unestimated displacements due to lack of infor-mation are smoothly filled in . Figure 5.8 shows some subtle improvements over Figure 5.7 due to addi t ional smoothing of existing visual observations using the Laplac ian and by u t i l iz ing L T S during the estimation. Unfortunately, this is hard to show on paper. T h e major improvement is in the smoothness of the direction of the Green's functions over the surface. Nevertheless, the perimeter of the area of large displacement (shown in red) is better defined after smoothing. T h e L T S step brings some visible improvement for vertex 55, especially on the far side of the tiger's head. T h e difference for vertex 38 is hardly visible. Figure 5.9 shows the same measurements wi th estimation done on a finer mesh. The finer mesh level is better able to resolve the changes in the Green's functions. A s a result the magni-tude of the y-column of the Green's function changes more smoothly which is visible as smoother color transitions in Figure 5.9. 5.5.2 Anatomic Soft-Tissue Human Wrist Model A s before wi th the plush-toy, the results for the local behavior of the soft-tissue wrist are discussed. T w o magnitude plots of forces and displacement during a probing are shown in Figure 5.11 for the locations indicated in Figure 5.10. The force profiles 94 (a) Probe at Vertex 38 of Figure 5.1(a) (b) Probe at Vertex 55 of Figure 5.1(b) :) Probe at Vertex 38 of Figure 5.1(a) (d) Probe at Vertex 55 of Figure 5.1(b) Figure 5.7: Hole fill ing Discrete Green's functions without hole filling are shown in 5.7(a) and 5.7(b), while the result of hole filling by Laplacian smoothing ( Section 5.3.3) is shown in 5.7(c) and 5.7(d). The magnitude of displacements are shown w i t h vertex coloring, where red indicates a large value while blue indicates the smallest value. Black vertices are not estimated. 95 (a) Probe at Vertex 38 of Figure 5.1(a) (b) Probe at Vertex 55 of Figure 5.1(b) (c) Probe at Vertex 38 of Figure 5.1(a) (d) Probe at Vertex 55 of Figure 5.1(b) Figure 5.8: Hole fill ing Results wi th Lap lac ian smoothing (Section 5.3.3) ut i l ized to smooth estimated values is shown in 5.8(a) and 5.8(b). The results are obtained w i t h 25 iterations wi th At = 0.2, 8 = 1 for the probed vertex and all fixed vertices, while 8 = 0.75 for vertices estimated based on visual information). T h e results for Least T r i m m e d Squares (LTS) as described in Section 5.3.2 are shown in 5.8(c) and 5.8(d). T h e magnitude of displacements are shown wi th vertex coloring, where red indicates a large value while blue indicates the smallest value. Black vertices are not estimated. 96 Figure 5.9: Discrete Green's Functions for Different M e s h Resolut ion T h e result for the base mesh of the subdivision hierarchy are shown in 5.9(a) and 5.9(b). T h e results for the next finer mesh resolution (see Section 5.4) are shown in 5.9(c) and 5.9(d). T h e estimation for both meshes is based on the same set of probes and images. Aga in , the magnitude of displacements are shown w i t h vertex coloring, where red indicates a large value while blue indicates the smallest value. Black vertices are not estimated. 97 (a) Vertex 4 (b) Vertex 98 Figure 5.10: Probe T i p at Vertex Loca t ion have similar structure to those of the plush toy. Dur ing the motion of the probe, the force is linear w i t h displacement. However, as the mot ion stops the displacement is maintained wi th an exponentially decreasing force. T h i s is a smal l viscoelastic effect. The stiffness at the two locations shown in 5.10 is quite different. Vertex 98, close to the mount of the object, is quite stiff while vertex 4, on the object's hand, is quite soft (notice the different displacement levels in Figure 5.11). T h e softness on the hand results from the hinging of the hand combined w i t h the longer lever relative to the mount ing surface. T h e varying stiffness of the soft-tissue wrist is also clearly shown in in F i g -ure 5.12(a) which is the plot of the distr ibution of the m a x i m u m singular values of the compliance matrices (the diagonal of the Green's function matr ix) . Figure 5.12(b) shows the residual of fit HP^Sj^ — U T | |2 for the local compliance matrices. T h e residual of fit is quite low and uniform over the mesh w i t h the exception of the thumb. The object's internal bone structure is simplified in the fingers. The struc-ture in the fingers contains no joints, does not extent the full length of the finger, 98 Time [1 Oms] Time [10ms] (a) Probe at Vertex 4 of Figure 5.10(a) (b) Probe at Vertex 98 of Figure 5.10(b) Figure 5.11: Magni tude of Force and Displacement and is fractured towards the finger's t ip . T h e global behavior of the stiff vertex 98 and the compliant vertex 4 is also very different. In Figure 5.13 the corresponding y-column of the discrete Green's functions mat r ix is i l lustrated. The Green's function for vertex 98 near the mount ing surface of the a rm is localized in the v ic in i ty of the vertex. However, the Green's function for vertex 4 on the hand of the human wrist model has no noticeable local effect. Tractions at this vertex 4 deform the complete human wrist model . T h e deformation is larger on the hand than on the lower arm. Th i s is a result of the lever effect due to the mounting of the object and the object's capabil i ty to hinge around its wrist . T h e physical wrist acts as a complex rotational joint causing points on the surface of the hand to describe an approximately circular trajectory. Th i s is approximated as a straight-line motion by the model . 99 (a) Solution Norm | | H L | | 2 (b) Residual of Fi t | |P T 3fc f c - U T | | 2 Figure 5.12: Loca l Compliance T h e figure is vertex colored, red indicates a large value while blue indicates the smallest value. Vertices not probed due to reachability constraints result i n EjTfc not being estimated and are shown in black. The result shown is obtained wi th a max imum of 24 probes per vertex location and observations from six different viewpoints. 100 (a) Probe at Vertex 4 of Figure 5.10(a) i f f mm (b) Probe at Vertex 98 of Figure 5.10(b) a i m l i f t (c) Probe at Vertex 4 of Figure 5.10(a) (d) Probe at Vertex 98 of Figure 5.10(b) Figure 5.13: Discrete Green's Functions for Different Mesh Resolut ion Base mesh level 0 of the subdivision hierarchy (5.13(a) and 5.13(b)) and level 1 (5.13(c) and 5.13(d)) are shown. T h e estimation for bo th meshes is based on the same observations (see Section 5.4). Aga in , the magnitude of displacements are shown wi th vertex coloring, where red indicates a large value while blue indicates the smallest value. Black vertices are not estimated. 101 C h a p t e r 6 Inverse Problem Solution T h e discrete Green's functions mat r ix estimation has been approached in Chapter 5 as a linear estimation problem. T h i s ignored the underlying structure of the Green's functions matr ix . T h e derivation given in Section 2.3.2 is based on a discretization of the boundary integral equation. Therefore, the solution by weighted residuals employed Ke lv in ' s fundamental solution for homogeneous, isotropic, linear-elastic materials. Th i s leads to a specific structure of the Green's functions mat r ix which is exploited below in Section 6.1 for an inverse problem solution to the estimation problem. T h e approach simplifies the estimation task significantly when applicable, however, it is also very restrictive. The assumption of an object made entirely out of a single homogeneous material is l imi t ing and would be a poor approximation for the examples of the plush toy tiger and the anatomic wrist model in Section 5.5. The plush toy is made out of cloth wi th a filling and most l ikely the head is separate from the body inside the cloth. T h e anatomic soft-tissue wrist model has an artificial bone structure which is in itself inhomogeneous. In addit ion, the soft-tissue is simulated wi th foam wi th a separate skin layer. The assumption of isotropy is also violated 102 since the toy-tiger seems to have a preferred hinging direction at the intersection between head and body, while the human wrist model behaves similar to the human wrist joint. (A linear elastic solid model also does in general not fit cloth.) The question then is: can an approximate homogeneous linear elastic model help with the estimation task? This is addressed in Section 6.2. 6 . 1 B o u n d a r y E l e m e n t M e t h o d f o r a n I n v e r s e E l a s t i c P r o b l e m The boundary element method is a convenient way to derive the discrete Green's functions matrix of the elastic behavior of a solid. The derivation is dependent on a fundamental solution of the material behavior in case of an infinite medium with the same material properties. Kelvin's fundamental solution for the 3-D Navier's equations describes the behavior of an infinite isotropic, homogeneous, linear-elastic material subjected to a concentrated unit point load. Kelvin's fundamental solution is the basis for the parameter estimation technique below, which can be considered a type 4b) of Kubo's [66] inverse problem classification. 6.1.1 Elastic Parameter Estimation The weights of Equation 2.15 are the terms p*j(£,x) a n d u*j(£,x)- These weights are calculated by Kelvin's fundamental solution of Equation 2.17 and Equation 2.16. The key observation is the fundamental solution can be separated into terms only de-pendent on geometry and terms only dependent on material constants. Equation 6.1 and Equation 6.2 group the factors depending on geometry and on material. Pij&x) = 1 - 2u ikn3 J- di dj dn 1 — v 8irr2 dr dr dr (6.1) 1 - v 8-7rr2 103 U ^ ' X j 7 ( l - i / ) 1 6 7 r r 7 ( 1 - 1 / ) Wirr' X ' T h e overall mat r ix equation of the boundary element formulation H u = G p is Equa t ion 6.3. Separating the geometric and material properties terms as in Equat ion 6.1 and Equa t ion 6.2 leads to Equa t ion 6.3 (Section 2.3.2 and A p p e n d i x A contain more detail). V ^ H A + T - ^ - H B + C 1 — V 1 — V u = - G A H jz T G B 7(1 - v) 7(1 - u) (6.3) w i th the diagonal matr ix Cu = ~ ^ij) 3 = 1 1 2i/ 1 «.4,, + . -//„. 1 - 1 / ^ 1 - 1 / Equa t ion 6.3 and Equat ion 2.1 form the basis for the inverse problem solu-t ion. A s described in Section 2.3.2, a given boundary configuration enables one to reorder H u = G p into the Green's functions matr ix Equat ion 2.1. Equa t ion 6.3 can be reordered in the same way. However, the separation between mater ial and geometric properties is lost, since the matrices H and G depend nonlinearly on the Poisson's ratio v. The boundary configuration employed in the following is the same as de-scribed in Section 5.1. The support ing surface is fixed, i.e., a displacement of zero is described and the free surface has its t ract ion prescribed. A s before, these values are combined in the block vector v . T h e boundary configuration remains fixed as different vertices on the free surface are probed. Th i s enables the combination of several probings in a common estimation step since the Green's functions mat r ix does not change. For each probing, some part of the right-hand side vector v is observed. Aga in , the observation at the probing location is the local measurement made wi th the robotic arm (see Section 3) while the global deformation is observed 104 by range-flow (as described in Section 4). A n opt imizat ion procedure to find shear modulus 7 and Poisson's ratio v from = v is described in the following Section 6.1.2. 6.1.2 Optimization T h e shear modulus 7 and Poisson's ratio v are found by an iterative opt imizat ion. F i rs t , the surface integrals i n Equa t ion 6.3 need to be calculated. T h i s is a very expensive computat ion (0(TNM) w i t h T the number of triangles, N the number of vertices and M the number of integration points wi th in an element), however, it needs to be done only once for a given triangular mesh. Note, that this is the same as for rendering the elastic response of an analyt ical ly derived model of an object [58]. T h e objective function in the min imiza t ion is the error function n F ( 7 , ^ ) = E ^ ( H ( 7 ^ ) v - v ) 2 . (6.4) A n element of the selection vector a in Equa t ion 6.4 is one if observations of vertex displacement are made and zero otherwise. The objective function of Equa t ion 6.4 is evaluated start ing from an in i t ia l guess for the unknowns 7 and v. In each iteration Equa t ion 6.3 needs to be evaluated w i t h the new values of 7 and v. T h e matrices H ^ , H g , C , Ga and G g remain unchanged and the fundamental solution matrices G and H are found wi th Equat ion 6.3. In the next step the Green's function matr ix 3 ( 7 , v) has to be calculated from the fundamental solution matrices G and H . T h i s calculation involves a matr ix inverse in the size of S (see Section 2.3.2). Th i s mat r ix inverse enables one to arrive at rows for which al l measurements are available because the complete described values v are known. The formulation is also t r iv ia l ly 105 expanded to mult iple probings in the same boundary configuration. Equa t ion 6.5 is the objective function for m measurements in the same boundary configuration. n m F ( l , V) = Yl [ Q 0 • • • arn]ij (S (7 , v) [Vi • • • V m ] - [ V i • • • V m ] ) ? . (6.5) i=0 j=0 T h e M a t l a b function fmins, an implementation of the Nelder -Mead simplex method, has l i t t le difficulty solving the minimiza t ion task at hand. However, note, that simpler methods should be sufficient because of the low dimensionali ty of the problem, the l imited range of possible values for Poisson's ratio v = 0.0 • • • 0.5 (in practice) and the shear modulus 7 acting like a scale factor. 6.1.3 Results for Half Ball T h i s section presents the results of the above described inverse problem solution for a test object: half of a soft ( N e r f ™ ) bal l (see Figure 4.3). T h e estimation results for the half bal l are summarized in Table 6.1. The linear elastic model is a good approximation for the bal l and the model fitting is successful. T h e estimated parameters v and 7, when used in simulat ion, produce results which closely resemble the recorded range-flow (see Figure 6.1). T h e results in Table 6.1 are from measured displacements by range-flow and displacement and force at a probed vertex. Table 6.1 also contains a comparison wi th a destructive test. The test is a simple compression test of a cyl indr ica l mate-rial sample. The material sample is compressed measuring force and the reduction in length of the cylinder. In practice, this works well for metal samples in common testing machines, however for soft elastic material it is far less precise. For the ball test object, the object behavior is greatly influenced by a membrane-like layer over 106 Node No . V 7 E 1 0.1183 0.7809 1.746 2 0.1878 0.6308 1.499 3 0.2227 0.6039 1.477 4 0.2516 0.6555 1.641 • 5 0.1857 0.6977 1.655 6 0.2009 0.6696 1.608 8 0.2089 0.6550 1.584 9 0.2007 0.6646 1.596 11 0.1868 0.7175 1.703 12 0.2082 0.7585 1.833 14 0.1997 0.7993 1.918 15 0.2084 0.7451 1.801 17 0.1192 0.7666 1.716 18 0.1366 0.7562 1.719 20 0.1400 0.7480 -1.705 22 0.1264 0.7121 1.604 26 0.1442 0.9383 2.147 A L L the above 0.1730 0.7283 1.709 Test: - - K, 1.75 Table 6.1: Mate r i a l Properties Es t imat ion Results T h e table lists the estimation results of the material properties (shear modulus 7, Poisson's ratio v and the corresponding Young's modulus E) for the half-ball. Each row in the table contains the result for a separate estimate based on the force probe location at the specified node(s). Us ing data from all probes in a single estimation results in values for the material properties very close to an independent compression test for the material . 107 its surface, i.e., the cylinder cut does need to have that membrane intact. T h e cy l in-drical sample of the ba l l is compressed in a mi l l ing machine. T h e mi l l ing machine ensures that the cylinder is compressed in the vertical direction only. T h e force is measured wi th an electronic balance (Ohaus CT 6000-S) under the sample. For a compressed cyl indr ica l material piece the Young's modulus is E = -j^, where F is the compression force, I and <5; is the height and the change i n height of the cylinder, respectively, and A the area of its circular end (see e.g., [101]). For homogeneous, isotropic, linear-elastic material , the relation between shear modulus 7 and Young's modulus E is 7 = ^(i+v) • T h i s allows the comparison in Table 6.1, which shows an excellent correspondence between the material-based compression test and the object-based inverse problem method. The variat ion of the estimates between sin-gle probes is mostly l ikely due to the misalignment of the probe wi th respect to the surface. Combin ing the measurement in one estimation step averages out the noise nicely. A visual comparison between the measured displacement and the simulated displacement for a probing is shown in Figure 6.1. T h e magnitude of surface displace-ment shows very good correspondence between the s imulat ion and the measurement on the surface area not occluded by the probe. 6.2 Approximate Model T h e above parameter estimation approach seems to work well if the object fits the assumptions of the model. A s previously stated, since these assumption are fairly strong, relaxing the assumptions is desirable. One approach to overcome the l imi ta t ion of homogeneity would be to t ry to estimate the location of interior boundaries between materials at the same time as the material properties. T h i s 108 I Figure 6.1: Compar ison of Measured and Simulated Displacement F r o m left to right: 3D range-flow, estimated displacement, simulated Displacement. T h e 3D range-flow shows the distr ibution of the estimated flow vectors w i th brighter range-flow corresponding to larger motions. T h e node displacement figure also shows larger movement being brighter. Dark surface is either stat ionary or no estimate is available. T h e estimated node displacement forms a par t ia l r ing around the location of the probe; the location itself and the other part of the r ing are occluded by the probe. obviously would be a formidable task but st i l l quite restrictive. Addi t iona l ly , this interior interface between materials is not required for the s imulat ion, only the Green's functions mat r ix is necessary. Instead of est imating the location of an interface, the response of the object to probing at vertex k can be modeled approximately wi th an apparent shear modulus 7k and an apparent Poisson's rat ion v\.. Th i s allows for varying stiffness over the object surface and variat ion in how global the object responds at a locat ion on the surface. It is l imi ted to object responses which remain continuous over the object's surface and are isotropic. The l imita t ions of the homogeneous isotropic linear elastostatic boundary element model are i l lustrated next w i t h the example of the anatomic soft-tissue human wrist model . A n apparent shear modulus 7^ and an apparent Poisson ratio vk are estimated for vertices k = 4 and k = 98 which are shown i n F igure 5.10. T h e estimates are based on 8 probes per vertex and only the local response is ut i l ized. T h e estimates are given in the following table. 109 Vertex No . • Ik . 4 98 0.04658 ' 0.0 19.28 1.157 T h i s approximate model is s t i l l too restrictive for the wrist. N o m i n i m u m can be found for the response at vertex 98 and instead the l imi t of vgg, = 0 is selected (This is found wi th constraint opt imizat ion wi th the Ma t l ab function fmincon). T h e apparent shear modu l i 7 £ are drastically different at the two locations. T h e lower a rm deforms smoothly when probed at Vertex 98 and the object's response is quite local wi th a cusp at the contact location (see Figure 6.2). T h e response has also no preferred direction. However, homogeneous isotropic material cannot produce this response given the boundary configuration: T h i s response is caused by the soft tissue of the object being compressed between the probe and the object's bone structure. Figures 6.2(a) to 6.2(b) compare the result of the param-eter estimation wi th the response estimated by the Green's function estimation of Chapter 5. Figure 6.2(c) shows the magnitude of the difference between the two methods for the y-column of the Green's functions block column Sj,98-However, some of the response of the human wrist model above the wrist joint can be approximated. The bone structure is not r igid enough to restrict the deformation and the complete object bends. Figures 6.2(d) and 6.2(e) show how both estimation method lead to a global bending of the object. T h e magnitude of the difference between the estimated displacements, shown in 6.2(f), is quite even. However, for deformation perpendicular to the axis of the joint greater differences can be expected. 110 (a) Approx. Inverse Problem (b) Green's Functions Estimation (c) Dif. be-tween 6.2(a) and 6.2(b) (d) Approx. Inverse Problem (e) Green's Functions Estimation (f) Dif. be-tween 6.2(d) and 6.2(e) Figure 6.2: Inverse P rob lem Solution and Discrete Green's Functions Es t ima t ion T h e Green's functions estimation method of Chapter 5 can represent deformations which are outside the domain of homogeneous isotropic linear elastostatic models. Figures 6.2(a), 6.2(b) and 6.2(c) show the behavior at vertex 98 of F igure 5.10(a), while 6.2(d), 6.2(e) and 6.2(f) show the comparison for vertex 4 of F igure 5.10(b). The magnitude of displacements are shown wi th vertex coloring, where red indicates a large value while blue indicates the smallest value. B lack vertices are not estimated. I l l 6.3 Summary T h i s section shows how the elasticity of an object can be estimated by solving an inverse problem for the elastic constants of the object. The method is restricted to objects made of isotropic, homogeneous, linear-elastic material . For objects which fit this assumption, the model works well as shown for the test object. T h e material restriction may be relaxed if the model is understood to be local and approximate only. A local approximate model based on linear elasticity may be employed by as-signing different elastic constants at each location on the surface. In fact, the elastic constants become elastic functions over the surface. However, internal structure and art iculat ion cannot be modeled in this way. B u t the discrete Green's function estimation of Chapter 5 s t i l l provides a reasonable approximation in the case of the anatomic soft tissue human wrist model . A n explicit model of articulated soft bodies is perhaps more promising and is the subject of future work. 112 C h a p t e r 7 C o n c l u s i o n s a n d F u t u r e W o r k 7.1 Conclusions This thesis presents a method to scan physical deformation behavior of objects. A model of the deformation behavior is estimated based on measurements of an intact existing object. A method to acquire these measurements including a measurement of the global deformation response is described. The global deformation response is observed w i t h a computer vision technique. The robotic measurement facility A C M E is extended to enable the data acquisition. The same techniques also allow one to validate an estimated deformable model. A C M E serves to acquire a geometric model of an object, to apply contact forces to the model and to record the deformation of the object. A C M E allows the acquisit ion of a large amount of deformation data in an automatic fashion, i.e., A C M E can scan physical deformation behavior of objects. The data acquisition dur ing the scan is based on two devices: a force sensor equipped robotic a rm and a t r inocular stereo-vision system. The robotic arm senses the posit ion and forces at the contact point wi th the object. T h e deformation of the object surface is sensed 113 visually. T h e technique developed for the visual sensing of the surface displacement is based on a robust combination of tr inocular optical flow and stereo. T h e addit ional images in the calculation of the opt ical flow enable a consistency check of the result. T h e combination of optical flow and stereo yields three-dimensional range flow. T h e registration wi th a separately acquired surface mesh of the object leads to estimates of the three dimensional vertex displacements during deformation. T h e acquired data can be uti l ized to estimate and validate deformable models. Th i s thesis develops and implements two estimation approaches. One ap-proach may be characterized as a data driven approach wi th the assumption of linear elastostatic behavior of object. Th i s allows one to fit a deformable model to a wide range of objects but requires a large amount of data. T h i s requirement of extensive data acquisition is mit igated by the use of the measurement robot A C M E . T h e other estimation approach makes the stronger a priori assumptions of homo-geneous linear elastostatic object behavior. T h i s estimation method requires l i t t le data because of this strong assumption. Nevertheless, it describes certain types of objects well; which are made of a single suitable material for instance an object made of closed-cell foam. Mode l ing assumptions in between these two extremes are not investigated but are suggested. The thesis contains an example of model fitting for a plush toy which has been demonstrated successfully to audiences at several scientific conferences (11th IRIS Precarn A n n u a l Conference, June 2001, Ot-tawa, Canada; 8th Int. Conf. of Computer V i s ion , Ju ly 2001, Vancouver, Canada) . Another example discussed is an anatomic soft tissue human wrist model . A de-formation model has also been estimated for a simple homogeneous foam object. To the best of the author's knowledge this thesis presents the first attempt to scan 114 deformation behavior of general objects. 7.2 Future Work One goal for future work could be to embed the presented techniques into a proba-bil ist ic framework based on a measurement theory. T h i s thesis characterizes mea-surements by a comparison w i t h ground-truth in the case of the displacement mea-surement in Chapter 4 and wi th residual of fit in the case of the estimation methods in Chapters 5 and 6. T h i s allows one to assess how well a model explains the measurements. In order to be able acquire a measurement and make a quantita-tive statement if this measurement is explained by a given model, a probabilist ic framework wi th in a measurement theory is necessary. A measurement theory estab-lishes the relations between observations and a physical model . The theory has to model the measurement process, reasoning about the influence of various physical phenomena - some in accordance wi th the model and some not captured by the model . The latter are generally referred to as noise. Bayesian estimation theory is one framework for quantitative statements about model fit and model selection. The other suggestions for future work concern est imating dynamic models of deformation. The trend in deformable models for medical s imulation is towards dynamic non-linear finite element models, despite the fact that current computers can render only very smal l and very soft models in real-time. Dynamic model fitting wi th a large deformation strain tensor would require the t racking of objects during the deformation process at least at twice the max imum occurring frequency. The techniques of this thesis s t i l l apply but the image acquisit ion wi th the Triclops stereo-vision system would most l ikely be too slow for the deformation transient. Higher speed camera systems and strobe lights could solve this problem. 115 Another option would be fitting articulated deformable models. T h i s problem is usually addressed wi th a pr ior i known kinematic model . T h e model of art iculat ion could be fit while min imiz ing the deformation of ind iv idua l l inks. In humans the area around the joints exhibit large deformation (e.g., finger joints) which would be one area to address. 116 A p p e n d i x A Linear Triangular Elements This A p p e n d i x gives the necessary notation to be able to evaluate integrals based on the fundamental solutions in Equations 2.16, 2.17, 6.1 and 6.2. T h e integrals are evaluated wi th linear triangular elements, i.e., the shape functions over a tr iangu-lar boundary element are linear in the values at the vertices of the triangle. T h e fundamental solutions contain the distance r between load location and response location. A component of the distance vector Ti given the vertices xl... x3 of the triangular element in anti-clockwise order can be calculated as rl = x \ - mx] - r]2x2 - (1 - m - r)2)x3 The distance is r = \JYA=\ ri, while the part ia l derivative of the distance is = rf and | ^ = "^v 1 - A component of the unit normal is ??,; = gi/\G\ (see Equa -t ion A . 7 ) . Next the mappings from local coordinates of an element to global coordinates and vice versa are given. Equat ion A . l maps local homogeneous surface coordinates rji (i.e., 773 = 1 — 771—772) of a triangular shape element to global Cartesian coordinates Xi. The vertices of the triangular element are labeled x\ ... x\ in anti-clockwise order. 117 Xl = 1-771 -I- x\f]2 + X3773 ( A - l ) Equat ions A . 2 , A . 3 and A . 4 map global to local coordinates. Vi m ^ 2 ^ 3 ™3™2 1 2 ' 1 2 1 _2 ~f" X X 2 2 3 — arxf + x1x\ ry 1 /yi 2 I ™3 ™2 1 2 ' 1 2 I /-*» 1 /*^ 3 rp2> ™3 x l x 2 — x\x\ ~\~ Xrj ™3™ 1 ™ 1 ™3 1 2 ' 1 2 x x 2 + a^Xr, i i — x x\ + x 2 x 3 rv» 1 2 1 3 2 1 2 ' 1 2 /T/-2 ™3 x l x 2 3 1 — X-yX^ ~\~ X-^X2 1 2 ' 1 2 X X 2 + xxx\ 2 2 — orxf + x 2 x | /-p 1 »^--2 i ™3~.2 1 2 1 2 1 /y 1 /y>3 i xj-^2 2 3 — x 1 x 2 3 \ — x^x2 ~j~ X ^ X 2 (A.2) (A-3) (A-4) T h e displacement and traction vector over an element are expressed in Equa-tions A . 5 and A . 6 , respectively, using identical interpolation functions as the above shape functions, i.e., (fi = rji. Uj (A.5> (A.6) The change to local coordinates requires the Jacobian in Equa t ion A . 7 . 9i 92 9i dF — \G\drjidrj2 with | G | = (A .7) i=l 3x2 9x3 8x2 dx3 dm dr)2 dm dm 9x3 <9xi dx3 dxi <9T?I <9?72 dm dm dxi 8x2 dxi dx2 dm dr]2 dr/2 dm 118 W i t h the local coordinates of Equa t ion A . l , the partials in Equa t ion A . 7 become: After the mappings between local and global coordinates have been estab-lished, it remains to evaluate the integrals. T h e fundamental solution integrals in Equations 2.18, 2.19 and Equa t ion 6.3 are evaluated numerically by Gaussian quadrature. Employ ing L integration steps per element weighted by a factor wi leads from Equat ion 2.19 to A.8, (The other integrals are converted identically.) H f c m = / p*(Sk,x)<p"m{x)dr(x) « E \ G \ i w i ( p t m <Pm(vi,m)i) (A.8) 1=1 T h e integration points and the corresponding weights for the Gaussian quadra-ture can be found in many references, e.g., [33]. Care has to be taken wi th the singular integrals for the self-effects (r —> 0 for k = m). L 119 A p p e n d i x B G e o m e t r i c M o d e l i n g w i t h A C M E T h i s appendix details the in i t ia l geometric model acquisition briefly introduced in Section 3.3. T h e main shape sensor in A C M E is the tr inocular stereo-head wi th in the F M S . T h e sensor produces large amounts of range data in close to real t ime em-ploying its commercial stereo vision S D K . However, range data from stereo is noto-riously noisy. Another feature of stereo ranging is its dependence on image features for matching between the stereo cameras. Add i t i ona l features can be projected onto a surface wi th a light or laser system. Color Triclops produces data of adequate but modest resolution (for a typical viewpoint in A C M E a depth resolution of 2-3 m m is achieved). Geometric object modeling requires that the entire surface of the object is covered, which requires varying viewpoints. Registrat ion of viewpoints in A C M E is achieved by visual calibration of the F M S wi th respect to an optical target, as well 120 as by exploit ing the accuracy of motions of the TestStation and the repeatability of mot ion commands for the F M S . T h e registration of the viewpoint relative to the target is accomplished wi th a variant of an Iterative Closest Point ( ICP) [8] algori thm. The starting point for the minimiza t ion is derived from the kinematics of A C M E . Range images from different viewpoints are combined in a voxel-based volumetric approach w i t h reconstruction software provided by the Nat ional Research Counc i l of Canada [109]. In order for this approach to succeed the data not only has to be in a common coordinate frame but also needs to be noise-filtered. Furthermore, normals for the point-cloud have to be estimated. The range data is filtered by calculating range w i t h variable mask sizes and combining them in a vot ing scheme. A n example of a point-cloud before and after filtering is shown in Figure B . l . The threshold for the filtering is set based on experience and is far from universal. However, since the stereo-head is quite common throughout the machine vis ion research community, the processing strategy and actual setting of thresholds may be of interest to some and are reported here: The filtering is achieved by processing a stereo image triplet twice. One stereo range image is calculated to maximize the accuracy of the result by using a small mask, while the other image is calculated to minimize noise wi th a large mask and tight validation settings. One stereo range image is calculated wi th a small stereo correlation mask of five in a 640 x 480 image. Addi t iona l ly , the Triclops stereo vision S D K ' s uniqueness threshold is kept nearly ineffective wi th a setting of 2.0 (The uniqueness threshold is intended to filter out matches which are not consistent between the right-left and the right-up camera pair) . T h e second stereo range image is calculated wi th a much larger correlation mask size of 13 and the tight validation setting of an uniqueness threshold of 1.0 combined wi th the connected 121 (a) Before Filtering. (b) After Filtering. (c) Before Filtering. (d) After Filtering. Figure B . l : Po in t -Cloud F i l t e r ing T w o point-clouds from Triclops filtered by variable mask size processing and plane fitting. Figures B . l ( a ) and B . l ( b ) shows the calibration target; Figures B . l ( c ) and B . l ( d ) the anatonomic soft-tissue human wrist model . (See Text for details). 122 region threshold of 300 pixels wi th less than 0.75 dispari ty difference. The images are combined in a plane fitting step simultaneously estimating normals and val idat ing a point. A l l points are rejected which are not wi th in 1cm between the two range images. T h e plane fit for each point in the noisy range image is accomplished by singular value decomposition ( S V D ) . A neighborhood of 5 x 5 is used when at least half (12) pixels wi th in the neighborhood have a val id depth. T h e plane fit is characterized by the ratio between smallest singular value over the sum of the two larger singular values. F i t t i n g is done i n metric x — y — z space and a fit wi th a singular value ratio of not more than 0.005 is accepted. The angles between the viewing axes to the point and the estimated normal is checked and the point are kept if the angle is less than 80 degrees. The filtered range images are input in the mesh generation software and a high-resolution mesh is generated. The resolution of the triangular mesh depends on the sampling density of the object's surface. T h e produced mesh contains typical ly some triangles due to support fixtures, some due to noise, and is not watertight. A multi-resolution subdivision mesh is buil t from this high-resolution raw mesh. In the process, holes are filled, triangles edited and the method for mesh simplification and subsequent subdivision refinement eliminates irregularities of the mesh. T h e simplification is done wi th Q S l i m [38]. The remaining process is essentially a dis-placement subdivis ion surface [72] or normal mesh [47] construction process. The implementation of this subdivision construction process is due to D . James and is described wi th more detail in [95]. The result for a plush toy tiger is shown in Figure 3.8. The stereo system also provides calibrated color imagery. These images are used for very simple texture mapping of the object. T h e vertices of the subdivision 123 model are projected into the color images wi th the pinhole camera model . A triangle of the mesh is assumed to be visible from a camera if al l its vertices and its centroid are visible. V i s i b i l i t y is calculated by ray intersection wi th the tr iangular mesh. T h e intersection test is accelerated by an octree space divis ion. The triangle is typical ly visible from a number of views but w i th different foreshortening and resolution. Out of the images where the triangle is visible, the.one wi th the highest viewing quali ty measure is selected for texture mapping. T h e viewing quali ty measure used is the product of image area covered by the triangle times the cosine between view direction and triangle normal . T h e relative global brightness of al l these color images is also adjusted. The produced texture map (see Figure 3.8(f)) contains many discontinuities due to inaccuracies of the mesh, misregistration and varying local brightness. T h e texture maps could be readily improved by blending of the local textures over the entire map el iminat ing abrupt brightness changes at edges and min imiz ing the effect of geometric errors, as in [108]. 124 B i b l i o g r a p h y [1] M . Akamatsu , G . Nakamura, and S. Steinberg. Identification of Lame coeffi-cients from boundary observations. Inverse Problems, 7(3):335-354, 1991. [2] P . Anandan . A computat ional framework and an algori thm for the measure-ment of visual motion. International Journal of Computer Vision, 2:283-310, 1989. [3] D . H . Ba l l a rd and O . A . K i m b a l l . R i g i d body motion from depth and optic flow. Computer Vision. Graphics and Image Processing, 22(1):95-115, 1983. [4] M . R . Banan , M . R . Banan , and K . D . Hjelmstad. Parameter est imation of structures form static response. I. Computa t iona l aspects. Journal of Struc-tural Engineering, 120(ll) :3243-3258, 1994. . [5] M . R . Banan , M . R . Banan , and K . D . Hjelmstad. Parameter estimation of structures form static response. II. Numer ica l simulation studies. Journal of Structural Engineering, 120(11):3259-3283, 1994. [6] J . L . Bar ron , D . J . Fleet, and S.S. Beauchemin. Performance of optical flow techniques. IJCV, 12(l):43-77, 1994. 125 [7] C . Basdogan, C . - H . Ho, M . A . Srinivasan, S .D. Smal l , and S .L . Dawson. Force interactions in laparascopic simulations: haptic rendering of soft tissues. In Medicine Meets Virtual Reality, pages 385-391, 1998. [8] P . J . Besl and N . D . M c K a y . A method for registration of 3-d shapes. IEEE Transactions on Pattern Recognition and Machine Intelligence, 14(2):239-256, Feb 1992. [9] L . M . Bezerra and S. Saigal. A boundary-element formulation for the inverse elastostatics problem ( IESP) of flaw detection. International Journal for Nu-' merical Methods in Engineering, 36(13):2189-2202, 1993. [10] L . M . Bezerra and S. Saigal. Inverse boundary traction reconstruction wi th the B E M . International Journal of Solids and Structures, 32(10):1417-1431, 1995. [11] D . Bielser, V . A . M a i w a l d , and M . H . Gross. Interactive cuts through 3-dimensional soft tissue. In P. Brunet and R . Scopigno, editors, Eurographics, Computer Graphics Forum, 18(3), pages C 3 1 - C 3 8 , Mi l ano , Italy, 1999. [12] J . Biggs and M . A . Srinivasan. Handbook of Virtual Environment Technology, chapter Hapt ic interfaces, K . M . Stanney (editor). Lawrence E r l b a u m Asso-ciates, in press. [13] M . J . Black and P. Anandan . T h e robust estimation of mult iple motions: parametric and piecewise-smooth flow fields. Computer Vision and Image Understanding, 63(1):75-104, 1996. 126 [14] R . C . Bolles and M . A . Fischler. A R A N S A C - b a s e d approach to model fitting and its application to finding cylinders in range data. In Int. Joint Conf. on Art. Intelligence, pages 637-643, Vancouver, Canada , 1981. [15] S. D i B o n a and O . Salvett i . A deformation model for biotissues behavior simulation. In International Conference on Image Processing, volume 2, pages 443-446, Vancouver, Canada, Sept 2000. [16] C A . Brebbia . Boundary element techniques in engineering. Pentech Press, London, U k , 1978. [17] C A . Brebbia , J . C . F . Telles, and L . C . Wrobe l . Boundary element techniques: theory and applications in engineering. Springer-Verlag, Ber l in , 1984. [18] M . Bro-Nielsen and S. C o t i n . Real-t ime volumetric deformable models for surgery simulation using finite elements and condensation. In J . Rossignac and F . X . Si l l ion , editors, Eurographics, Computer Graphics Forum, 15(3), pages 57-66, Poitiers, France, 1996. [19] I. Brouwer, J . Us t in , L . Bentley, A . Sherman, N . Dhruv , and F. Tendick. Measur ing in vivo animal soft tissue properties for haptic modeling in surgical simulation. In Medicine Meets Virtual Reality, pages 69-74, Jan 2001. [20] R . Carceroni and K . Kutu lakos . Mul t i -v i ew scene capture by surfel sampling: from video streams to non-rigid 3d motion, shape and reflectance. In Inter-national Conference on Computer Vision, volume 2, pages 60-67, 2001. [21] S. Chaudhur i and S. Chatterjee. Es t imat ion of mot ion parameters for a de-formable object from range data. In Computer Vision and Pattern Recogni-tion, pages 291^295. I E E E , 1989. [22] S.S. Chen and M . A . Penna. Shape and motion of nonrigid bodies. Computer Vision, Graphics and Image Processing, 36(2/3):175-207, 1986. [23] L F . Cos ta and R . Balaniuk. L e m - an approach for real-time physically based soft tissue simulat ion. In International Conference on Robotics and Automa-tion, pages 2337-2343, Seoul, South Korea , M a y 2001. [24] S. C o t i n , H . Delingette, and N . Ayache. Real- t ime elastic deformations of soft tissues for surgery simualt ion. IEEE Trans, on Visulization and Computer Graphics, 5( l ) :62-73, 1999. [25] S. C o t i n , H . Delingette, and N . Ayache. A hybr id elastic model allowing real-time cut t ing, deformations and force-feedback for surgery t raining and simulat ion. The Visual Computer, 16(8):437-452, 2000. [26] D . d 'Aul ignac , R . Balaniuk, and C . Laugier. A haptic interface for a vir-tual exam of the human thigh. In International Conference on Robotics and Automation, pages 2452-2457, San Francisco, U S A , A p r i l 2000. [27] D . d 'Aul ignac , M . C . Cavusoglu, and C . Laugier. Mode l ing the dynamics of the human thigh for a realistic echographic simulator w i th force feedback. In Int. Conf. on Medical Image Computer-Assisted Intervention, pages 1191-1198, Cambridge, U k , Sept 1999. [28] S. De and M . A . Srinivasan. T h i n walled models for haptic and graphical rendering of soft tissues in surgical simulations. In Medicine Meets Virtual . Reality, pages 94-99, 1999. 128 [29]: F . B o u x de Casson, D . d 'Aul ignac , and C . Laugier. A n interactive model of the human liver. In 7th International Symposium on Experimental Robotics, Honolulu , U S A , December 2000. [30] G . Debunne, M . Desbrun, M . - P . C a n i , and A . H . Bar r . Dynamic real-time deformations using space &; t ime adaptive sampling. In Computer Graphics, Annual Conference Series, pages 31-36, Los Angles, A u g 2001. A C M S I G -G R A P H . [31] H . Delingette. Toward realistic soft-tissue model ing in medical s imulat ion. Proceedings of the IEEE, 86(3):512-523, 1998. [32] M . Desbrun, M . Meyer, P. Schroder, and A . H . Bar r . Implici t fairing of ir-regular meshes using diffusion and curvature flow. In Computer Graphics, Annual Conference Series, pages 317-324, Los Angles, U S A , A u g 1999. A C M S I G G R A P H . [33] J . Dominguez. Boundary Elements in Dynamics. Computa t iona l Mechanics Publ icat ions , Southhamptom, U k , 1993. [34] J .S. Duncan , R . L . Owen, L . H . Staib, and P. Anandan . Measurement of non-rigid motion using contour shape descriptors. In Computer Vision and Pattern Recognition, pages 318-324. I E E E , 1991. [35] P. Fua . A parallel stereo algori thm that produces dense depth maps and preserves image features. Machine Vision and Applications, 6:35-49, 1993. [36] F . Ganovell i , P. Cignoni , C . Montan i , and R Scopigno. A multiresolution model for soft objects support ing interactive cuts and lacerations. In M . Gross 129 and F . R . A . Hopgood, editors, Eurographics, Computer Graphics Forum, 19(3), pages C271-C282, Interlaken, Switzerland, 2000. [37] L . Gao , K . J . Parker, R . M . Lerner, and S .F . Levinson. Imaging of the elastic properties of tissue - a review. Ultrasound in Medicine and Biology, 22(8) :959-977, 1996. [38] M . Gar l and and P.S. Heckbert. Surface simplification using quadric error . metrics. In Computer Graphics, Annual Conference Series, pages 209-216, Los Angles, U S A , A u g 1997. A C M . [39] J . C . Ge l in and O . Ghoua t i . Inverse method for material parameters estimation in the inelastic range. Computational-Mechanics, 16(3):143-150, 1995. [40] S . F . F . Gibson and B . M i r t i c h . A survey of deformable modeling in computer graphics. Technical Report TR-97-19, Mi t sub i sh i Electr ic Research Labora -tory, Cambridge, M A , U S A , Nov. 1997. [41] C . G ioda . Indirect identification of the average elastic characteristics of rock masses. In Structural foundations on rock, volume 1, pages 65-75, 1980. [42] S. Gottschalk, M . L i n , and D . Manocha . OBB- t ree : A hierarchical structure for rapid interference detection. In Computer Graphics, Annual Conference Series, pages 171-180, New Orleans, U S A , A u g 1996. A C M S I G G R A P H . [43] S. Govindjee and P . A . M i h a l i c . Computa t iona l methods for inverse finite elas-tostatics. Computer Methods in Applied Mechanics and Engineering, 136(1-2):47-57, 1996. 130 [44] S. Govindjee and P . A . M i h a l i c . Computa t iona l methods for inverse deforma-tions in quasi-incompressible finite elasticity. Int. Journal for Numer. Meth. in Engineering, 43:821-838, 1998. [45] M . Grediac, C . Toukourou, and A . Vau t r in . Inverse problem in mechanics of structures: a new approach based on displacement field processing. In M . Tanaka and H . D . B u i , editors, IUATM Symposium, Inverse Problems in Engineering Mechanics, Be r l in , Germany, 1993. Springer-Verlag. [46] M . Grediac and A . Vau t r in . A new method for determination of bending rigidities of th in anisotropic plates. Transactions of the ASME, pages 964-968, 1990. [47] I. Guskov, K . Vid imce , W . Sweldens, and P. Schroder. N o r m a l meshes. In Computer Graphics. Annual Conference Series, pages 95-102, New Orleans, U S A , J u l 2000. A C M S I G G R A P H . [48] P . C . Hansen. Regularizat ion tools: A Ma t l ab package for analysis and solution of discrete ill-posed problems. Numerical Algorithms, 6(1/2):1—35, 1994. [49] J . D . He lm, M . A . Sutton, and S.R. Ne i l . Three-dimensional image correlation for surface displacement measurements. In S. F . E l - H a k i m , editor, Videomet-rics III, volume 2350, pages 32-45. S P I E , 1994. [50] K . H i ro t a and T . Kaneko. Representation of soft objects in v i r tua l environ-ment. In 8th Int. Conf. on Artificial Reality and Tele-Existence (ICAT), pages 59-62, Tokyo, Jp, Dec 1998. [51] B . K . P . H o r n and J . G . Harr is . R i g i d body motion from range image sequences. CVGIP: Image Understanding, 53(1):1-13, 1991. [52] B . K . P . H o r n and B . G . Schunk. Determining optical flow. Artificial Intelli-gence, 17(l-3):185-203, 1981. [53] B . Horowitz and A . Pent land. Recovery of non-rigid motion and structure. In Computer Vision and Pattern Recognition, pages 325-330. I E E E , 1991. [54] S . -C. Hsieh and T . M u r a . Nondestructive cavity identification i n structures. International Journal of Solids and Structures, 30(12):1579-1587, 1993. [55] T . S . Huang and S.D. Bloste in . Robust algorithms for motion estimation based on two sequential stereo image pairs. In Computer Vision and Pattern Recog-nition, pages 518-523. I E E E , 1985. [56] D . James. Multiresolution Green's Function Methods for Interactive Simula-tion of Large-scale Elastostatic Objects and Other Physical Systems in Equilib-rium. P h D thesis, Univers i ty of B r i t i s h Columbia , expected November 2001. [57] D . James and D . K . P a i . A unified treatment of elastostatic contact simulation for real t ime haptics. Haptics-e, 2(1), 2001. [58] D . L . James and D . K . P a i . Ar tDefo accurate real t ime deformable objects. In Computer Graphics, Annual Conference Series, pages 65-72, Los Angles, U S A , A u g 1999. A C M S I G G R A P H : [59] A . Joukhadar, C . B a r d , and C . Laugier. P lann ing dextrous operations using physical models. In International Conference on Robotics and Automation, pages 748-753. I E E E , 1994. [60] A . Joukhadar, F . Garat , and C . Laugier. Parameter identification for dynamic simulation. In International Conference on Robotics and Automation, pages 1928-1933, Albuquerque, N M , U S A , 1997. I E E E . [61] F . K a l l e l and M . Ber t rand . Tissue elasticity reconstruction using linear per-turbat ion method. IEEE' Transactions on Medical Imaging, 15(3):299-313, 1996. [62] C . Kambhamet tu , D . B . Goldgof, and M . He. Determinat ion of motion param-eters and estimation of point correspondences in small nonrigid deformations. In Computer Vision and Pattern Recognition, pages 943-946. I E E E , 1994. [63] T . Kanade, A . Yoshida , K . Oda , H . Kano , and M . Tanaka. A stereo machine for video-rate dense depth mapping and its new applicat ion. In Computer Vision and Pattern Recognition, pages 196-202. I E E E , 1996. [64] E . Keeve, S. G i r o d , R . K i k i n s , and B . G i r o d . Deformable modeling of facial tissue for craniofacial surgery simulation. Computer Aided Surgery, 3(5):228-238, 1998. [65] R . M . K o c h , M . H . Gross, F . R . Carls , D . F . von Bi i r en , G . Fankhauser, and Y . I . H . Par ish . Simulat ing facial surgery using finite element models. In Com-puter Graphics, Annual Conference Series, pages 421-428, New Orleans, U S A , A u g 1996. A C M S I G G R A P H . [66] S. K u b o . Inverse problems related to the mechanics and fracture of solids and structures. JSME International Journal, Series 1, 31(2):157-166, 1988. [67] U . Ki ihnapfe l , H . K . Cakmak , and H . Maafi . Endoscopic surgery t raining using vi r tua l reality and deformable tissue simaulation. Computers & Graphics, 24:671-682, 2000. [68] S . -H. L a i and B . C . Vemur i . Reliable and efficient computat ion of optical flow. International Journal of Computer Vision, 29(2):87-105, 1998. 133 [69] J . Lang . Contour stereo matching by correlation in a constraint satisfaction approach. In Vision Interface, pages 19-26, 1998. [70] J . L a n g and D . K . P a i . Bayesian estimation of distance and surface normal w i th a time-of-flight laser rangefinder. In International Conference on 3-D . Digital Imaging and Modeling. I E E E , October 1999. [71] J . L a n g and D . K . P a i . Es t imat ion of elastic constants from 3D range-flow. In 3rd International Conference on 3-D Digital Imaging and Modeling, pages 331-338, Quebec Ci ty , Canada, 2001. I E E E . [72] A . W . F . Lee, H . More ton , and H . Hoppe. Displaced subdivision surfaces. In Computer Graphics, Annual Conference Series, pages 85-94, New Orleans, U S A , J u l 2000. A C M S I G G R A P H . [73] Y . Lee, D . Terzopoulos, and K . Waters. Realist ic modeling of facial animation. In Computer Graphics, Annual Conference Series, pages 55-62, Los Angles, U S A , A u g 1995. A C M S I G G R A P H . [74] W . - H . L iao , S .J . Aggarwal , and J . K . Aggarwal . The reconstruction of dynamic 3D structure of biological objects using stereo microscope images. Machine Vision and Applications, 9(4):166-178, 1997. [75] J . L l o y d and V . Hayward . Trajectory generation in m u l t i - R C C L . Journal of Robotic Systems, 10(3):369-390, 1993. [76] J . E . L l o y d , J .S. Beis, D . K . Pa i , and D . G . Lowe. Model-based telerobotics wi th vision. In International Conference on Robotics and Automation, pages 1297-1304. I E E E , 1997. 134 [77] J . Louchet, X . Provot, and D Crochemore. Evolu t ionary identification of cloth animation models. In Computer Animation and Simulation Workshop, Eurographics, pages 245-252, Maastr icht , Netherlands, Sept 1995. [78] H . Maafi and U . Ki ihnapfe l . Noninvasive measurement of elastic properties of l iv ing tissue. In 13th Int. Congress on Comp. Assisted Radiology (CARS), pages 865-870, 1999. [79] K . T . M c D o n n e l l , H . Q i n , and R . A . Wloda rczyk . V i r t u a l clay: A real-time sculpt ing system wi th haptic toolkits . In Proc. Symposium on Interactive 3D Graphics, pages 179-190. A C M , M a r 2001. [80] T . Merz , D . Paulus, and H . Niemann. L ine segmentation for interferograms of continously deforming objects. In F . S . C h a u ; C . T . L i , editor, Int. Conference on Experimental Mechanics: Advances and Applications, volume 2921, pages 325-330. S P I E , 1997. [81] D . Metaxas . Physics-Based Deformable Models - Applications to Computer Vision. Graphics and Medical Imaging. K luwer Academic Publisher , Boston, 1996. [82] D . Metaxas and D . Terzopoulos. Shape and nonrigid motion estimation through physics-based synthesis. IEEE Transactions on Pattern Recognition and Machine Intelligence, 15(6):580-591, 1993. [83] S . M . Metwal l i , A . R . Ragab, A . H . Karnel , and A . A . Saheb. Determinat ion of plastic stress-strain behavior by digital-image-processing technique. Proc. -of the Society for Experimental Mechanics, Inc., 44:414-422, 1987. 135 [84] G . Nakamura and K . Tanuma. A nonuniqueness theorem for an inverse boundary value problem in elasticity. SIAM Journal on Applied Mathematics, 56(2):602-610, 1996. [85] G . Nakamura and G . U h l m a n n . Identification of Lame parameters by bound-ary measurements. American Journal of Mathematics, 115(5):1161-1187, 1993. [86] G . Nakamura and G . U h l m a n n . G loba l uniqueness for an inverse boundary-problem arising in elasticity. Inventiones Mathematicae, 118(3):457-474, 1994. [87] G . Nakamura and G . U h l m a n n . Inverse problems at the boundary for an elastic medium. SIAM Journal on Mathematical Analysis, 26(2):263-279, 1995. [88] J . - C . Nebel . Soft tissue model l ing from 3D scanned data. In Deform, Geneva, Switzerland, Nov. 2000. [89] J . - C . Nebel , F . J . Rodr iguez-Miguel , and W . P. Cockshott . Stroboscobic stereo rangefinder. In 3rd International Conference on 3-D Digital Imaging and Modeling, pages 59-64, Quebec Ci ty , Canada, 2001. I E E E . [90] B . Novotny. Some aspects of inverse problem analysis for a layered half-space. Computational Mechanics, 14(3):266-276, 1994. [91] J . F . O ' B r i e n and J . K . Hodgins. Graphica l modeling and animat ion of bri t t le fracture. In Computer Graphics, Annual Conference Series, pages 137-146, Los Angles, U S A , A u g 1999. A C M S I G G R A P H . [92] D . K . Pa i , editor. Reality-based modeling and applications in reverse engineer-ing, computer graphics and VR, Workshop. I E E E , 2000. 136 [93] D . K . Pa i , J . Lang , J . E . L l o y d , a n d . J . L . R ichmond . Reality-based modeling wi th A C M E : A progress report. In D . Rus and S. Singh, editors, 7th Inter-national Symposium on Experimental Robotics, LNCIS 271, pages 121-130, Honolu lu , U S A , December 2000. Springer-Verlag. [94] D . K . P a i , J . Lang , J . E . L l o y d , and R . J . Woodham. A C M E , a te lerobot ic active measurement facility. In P. Corke and J . Trevelyan, editors, 6th International Symposium on Experimental Robotics, LNCIS 250, pages 391-400, Sydney, Aus t ra l i a , M a r c h 1999. Springer-Verlag. [95] D . K . Pa i , K . van den Doel , D . L . James, J . Lang , J . E . L l o y d , J . L . R ichmond , and S . H . Y a u . Scanning physical interaction behavior of 3D objects. In Com-puter Graphics, Annual Conference Series, pages 87-96, Los Angles, A u g 2001. A C M S I G G R A P H . [96] X . Papademetris. Estimation of 3D left ventricular deformation from medical images using biomechanical models. P h D thesis, Yale University, M a y 2000. [97] H . W . Park, S. Shin , and H . S . Lee. Determinat ion of an op t imal regularization factor in system identification w i t h Tikhonov regularization for linear elas-tic continua. International Journal for Numerical Methods in Engineering, 51:1211-1230, 2001. [98] A . Pent land. Au tomat i c extraction of deformable part models. International Journal of Computer Vision, 4(2):107-126, 1990. [99] G . Pic inbono, H . Delingette, and N . Ayache. Non-linear and anisotropic elastic soft tissue models for medical s imulat ion. In International Conference on Robotics and Automation, pages 1371-1376, Seoul, South Korea , M a y 2001. 137 [100] S . M . P ia t t and N . I . Badler . A n i m a t i n g facial expression. In Computer Graph-ics, Annual Conference Series, pages 245-252, Dallas, U S A , A u g 1981. A C M S I G G R A P H . [101] E . P . Popov. Engineering mechanics of solids. Prentice H a l l , Englewood Cliffs, N J , U S A , 1st edition, 1990. [102] W . H . Press, S .A . Teukolsky, W . T . Vetterl ing, and B . P . Flannery. Numerical recipes in C: the art of scientific computing. Cambridge Univers i ty Press, 2 edition, 1992. [103] K . R . Raghavan and A . E . Yagle. Inverse and opt imal drive problems in elas-t ici ty imaging of soft tissue. In Medical Imaging Conference, pages 1886-1890. I E E E , 1993. [104] R . Ramanathan and D . Metaxas. Dynamic deformable models for enhanced haptic rendering in v i r tua l environments. In Virtual Reality, pages 31-35. I E E E , 2000. [105] Point Grey Research. Triclops on-line manual, h t tp : / /www.p tgrey .com/ . [106] D . Reznik and C . Laugier . Dynamic simulation and v i r tua l control of a de-formable fingertip. In International Conference on Robotics and Automation, pages 1669-1674. I E E E , 1996. [107] Zantout R . N . and Zheng Y . F . Geodesies: A tool for solving material properties inverse problems. In Int. Conference on Industrial Technology, pages 391-394. I E E E , 1994. 138 [108] C . Rocchin i , P. Cignoni , C . Montan i , and P . Scopigno. M u l t i p l e texture stitch-ing and blending on 3d objects. In 10th Workshop on Rendering, Eurographics, pages 173-180, Granada , Spain, 1999. [109] G . R o t h and E . Wibowoo . A n efficient volumetric method for bui ld ing closed triangular meshes from 3-d image and point data. In Proc. Graphics Interface, pages 173-180, 1997. [110] P . J . Rousseeuw. Least median of squares regression. Journal of the American Statistical Association, 79(388):871-880, 1984. [ I l l ] P . J . Rousseeuw and K . V a n Driessen. Compu t ing L T S regression for large data sets. Technical report, Univers i ty of Antwerp, 1999. [112] M . A . Sagar, D . Bul l ivant , G . D . Mal l inson , P . J . Hunter , and I . W . Hunter. A v i r tua l environment and model of the eye for surgical s imulat ion. In Computer Graphics, Annual Conference Series, pages 205-212, Orlando, U S A , J u l 1994. A C M S I G G R A P H . [113] J . K . Salisbury and M . A . Srinivasan. Phantom-based haptic interaction wi th v i r tua l objects. IEEE Computer Graphics & Applications, 17(5):6-10, 1997. [114] M . Sanayei and M . J . Saletnik; Parameter estimation of structures form static strain measurements. I. Formulat ion. Journal of Structural Engineering, 122(5):555-562, 1996. [115] M . Sanayei and M . J . Saletnik. Parameter estimation of structures form static strain measurements. II. Er ro r sensitivity analysis. Journal of Structural En-gineering, 122(5):555-562, 1996. 139 , [116] Y . Sato, M . Wheeler, and K . Ikeuchi. Object shape and reflectance modeling from observation. In Computer Graphics. Annual Conference Series, pages. 379-387, Los Angles, U S A , A u g 1997. A C M . ; [117] D . S . Schnur and N . Zabaras. A n inverse method for determining elastic-. . material properties and a material interface. International Journal for Nu-merical Methods in Engineering, 33(10):2039-2057, 1992. [118] A . A . Shabana. Dynamics of Multibody Systems. John W i l e y & Sons, New York , U S A , 1989. [119] H . Spies, B . Jahne, and J . L . Bar ron . Dense range flow from depth and intensity data. In International Conference on Pattern Recognition, volume 1, pages 131-134, Barcelona, Spain, 2000. [120] H . Spies, B . Jahne, and J . L . Bar ron . Regularised range flow. In European Con-ference on Computer Vision, L C N S 1843/44, pages 785-799, D u b l i n , Ireland, 2000. I E E E , Springer-Verlag, Ber l in . [121] G . E . Stavroulakis and H . Antes. Nondestructive elastostatic identification of unilateral cracks through bem and neural networks. Computational-Mechanics, 20(5):439-451, 1997. [122] G . Szekely, C . Brechbiihler, J . D u a l , R . Enzler , J . Hug , R . Hutter , N . Ironmon-ger, M . Kauer , V . Meier, P. Niederer, A . Rhomberg, P . Schmid, G . Schweitzer, M . Thaler , V . Vuskovic, G . Troster, U . Haller , and M . Bajka. V i r t u a l reality-based simulat ion of endoscopic surgery. Presence, 9(3):310-333, 2000. 140 [123] F . Tendick, M . Downes, T . Goktek in , M . C . Cavusoglu, D . Feygin, X . W u , • R . E y a l , M . Hegarty, and L . W . Way. A v i r tua l environment testbed for train-ing laparoscopic surgical skills. Presence, 9(3):236-255, 2000. [124] D . Terzopoulos. Regularizat ion of inverse visual problems involving disconti-nuities. IEEE Transactions on Pattern Recognition and Machine Intelligence, 8(4):413-424, 1986. [125]- D . Terzopoulos and K . Fleischer. Mode l ing inelastic deformation: viscoelas-ticity, plasticity, fracture. In Computer Graphics, Annual Conference Series, pages 269-278, At lan ta , U S A , A u g 1988. A C M S I G G R A P H . [126] D . Terzopoulos and K . Waters. Physical ly-based facial modeling, analysis and animation. Journal of Visualization and Computer Animation, 1:73-80, 1990. [127] D . Terzopoulos, A . P . W i t k i n , and M . Kaas . Constraints on deformable models: Recovering 3D shape and nonrigid motion. Artificial Intelligence, 36(1):91-123, 1988. [128] S. Vedula , S. Baker. P. Rander, R . Col l ins , and T . Kanade. Three-dimensional -.'. scene flow. In International Conference on Computer Vision, pages 722-729, 1999. [129] V . Vuskovic, M . Kauer , G . Szekely, and M . Reidy. Realist ic force feedback for v i r tua l reality based diagnostic surgery simulators. In International Con-ference on Robotics and Automation, pages 1592-1598, San Francisco, U S A , A p r i l 2000. 141 [130] R . J . Woodham, E . Catanzar i t i , and A . K . Mackwor th . Analys is by syntheis ) . in computat ional vision wi th applicat ion to remote sensing. Computational Intelligence, l(2):71-79, 1985. [131] X . W u , M . S . Downes, T . Gotekin , and F . Tendick. Adapt ive nonlinear finite elements for deformable body simulat ion using dynamic progressive meshes. In A . Chalmers and T . - M . Rhyne , editors, Eurographics, Computer Graphics Forum, 20(3), pages C349-C358, Manchester, U k , 2001. [132] M . Yamamoto , P . Boulanger, J . - A . Bera ld in , and M . R ioux . Direct estimation of range flow on deformable shape from a video rate range camera. IEEE Transactions on Pattern Recognition and Machine Intelligence, 15(l):82-89, 1993. [133] Y . Zhang and C . Kambhamet tu . Integrated 3d scene flow and structure re-covery from mult iview image sequences. In Computer Vision and Pattern Recognition, volume 2, pages 674-681, 2000. [134] Z . Zhang. Parameter estimation techniques: A tutorial w i th application to conic fitting. Image and Vision Computing, 15(l):59-76, 1997. [135] Z . - H . Zhong. Finite Element Procedures for Contact-Impact Problems. Oxford Univers i ty Press, Oxford, U k , 1993. [136] Y . Zhuang and J . Canny. Hapt ic interactions wi th global deformations. In International Conference on Robotics and Automation, pages 2428-2433, San Francisco, U S A , A p r i l 2000. 142
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Deformable model acquisition and validation
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Deformable model acquisition and validation Lang, Jochen 2001
pdf
Page Metadata
Item Metadata
Title | Deformable model acquisition and validation |
Creator |
Lang, Jochen |
Date Issued | 2001 |
Description | Objects deform in response to contact forces exerted on them. The deformation depends on material properties, the geometry of the object and external forces. This thesis develops a robotic system for automatically acquiring observations of a deforming object and for estimating a model of the deformation from these observations. Models of deformable objects are in wide-spread use in simulation, computer graphics and virtual reality. Deformation, impact and fitting simulation aid manufacturing. In computer graphics deformable objects are designed and animated. Medical simulators incorporating physical models of organs and tissue are a significant emerging virtual reality application. The material properties of deformable models are often assigned based on mechanical (and other) testing of material samples. Material samples do not represent commonly simulated objects well if there is a high variability in the material of the object, e.g., due to mixed material, unknown material, or material with imperfections. In contrast to material sampling, this thesis develops a method to scan the deformation behavior of a complete object. The scanning is analogous to the scanning of the visual appearance of an object. The scan captures the individual response of a complete object to contact forces. The result of the scan is not only a deformable model but also data which serves to validate the model. This validation step gives a qualitative and sometimes quantitative assessment of the suitability of a model. |
Extent | 6342662 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-10-08 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0051715 |
URI | http://hdl.handle.net/2429/13753 |
Degree |
Doctor of Philosophy - PhD |
Program |
Computer Science |
Affiliation |
Science, Faculty of Computer Science, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2001-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2001-714896.pdf [ 6.05MB ]
- Metadata
- JSON: 831-1.0051715.json
- JSON-LD: 831-1.0051715-ld.json
- RDF/XML (Pretty): 831-1.0051715-rdf.xml
- RDF/JSON: 831-1.0051715-rdf.json
- Turtle: 831-1.0051715-turtle.txt
- N-Triples: 831-1.0051715-rdf-ntriples.txt
- Original Record: 831-1.0051715-source.json
- Full Text
- 831-1.0051715-fulltext.txt
- Citation
- 831-1.0051715.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0051715/manifest