UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Deformable model acquisition and validation Lang, Jochen 2001

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

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

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 p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t of the requirements f o r 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 t h a t the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and study. I f u r t h e r agree t h a t p e r m i s s i o n f o r e x t e n s i v e c o p y i n g of t h i s t h e s i s f o r s c h o l a r l y purposes may be g r a n t e d by the head of my department or by h i s or her r e p r e s e n t a t i v e s . I t i s u n d e r s t o o d t h a t c o p y i n g or p u b l i c a t i o n of t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l not be a l l o w e d without my w r i t t e n p e r m i s s i o n .  Aeu  Department of  -Severe e  The U n i v e r s i t y of B r i t i s h Columbia Vancouver, Canada Date  Qc.4.  a  Abstract 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.  ii  Contents  Abstract  ii  Contents  iii  Acknowledgements  vii  1 Introduction 1.1  1.2  1.3  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  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  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  2.4  2.5  Solution Methods  22  2.3.1  Finite Element Model  23  2.3.2  Boundary Element Model  25  Other Physics Based Deformable Models  28  2.4.1  Particle Systems  28  2.4.2  Deformable Superquadrics  30  Summary  31  3 The A C M E Facility  32  3.1  4  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  Deformation Measurement 4.1  47  Robust Flow Estimation  49  iv  4.2  4.3  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  Estimation of Vertex Displacement  62  4.2.1  63  Displacement Results  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 6.1  102  Boundary Element Method for an Inverse Elastic Problem  103  6.1.1  103  Elastic Parameter Estimation  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 7.1  113  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  vi  Acknowledgements 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. U r i 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 ( U B C ) Active Measurement facility ( A C M E ) enabled the experimental 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 ( P G R ) for help with the P G R Triclops Stereo Vision S D K . Furthermore, I like to acknowledge Dr. Gerhard Roth at the National Research Council, Canada for making his volumetric mesh creation software available 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 University 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 .  JOCHEN  The University of British Columbia October 2001  vii  LANG  Chapter 1  Introduction  Interactive s i m u l a t i o n of physical objects is desirable for numerous applications i n the areas of v i r t u a l and augmented reality. T h e goal i n interactive s i m u l a t i o n is to ensure v i r t u a l objects resemble their real counterparts as closely as possible but at the same t i m e m a i n t a i n interactivity. T h e consequence of this trade-off is the need to employ an approximate m o d e l of the physical behavior of objects. In general, the realism m a t h e m a t i c a l models of objects can achieve depends on the proper choice of the physical parameters of the model. H o w should these physical parameters be chosen?  T h i s thesis argues that realistic simulated behavior requires parameters  w h i c h 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. T h i s thesis shows how to efficiently scan this behavior i n a robotic measurement  facility very similar to the scanning of object geometry i n a shape  scanner. Simulations gain realism if the effect of forces on an object is represented realistically.  Forces govern m a n y behaviors of real world objects.  I n particular,  forces deform an object depending on its geometric and material properties.  1  The  simulation of realistic v i r t u a l deformations necessitates the inclusion of forces, object geometry and realistic material properties.  R e a l i s t i c deformation of v i r t u a l  objects is required i n medical t r a i n i n g simulators [31], i n haptic interfaces to v i r t u a l worlds [113], and i n design of robotic tasks i n v o l v i n g soft materials [59, 60]. M e d ical simulators allow students to practice [123], and help to predict and p l a n some surgical procedures [67]. M e d i c a l robotics can also profit from user interfaces w i t h haptic feedback [112]. Force s i m u l a t i o n plays a very i m p o r t a n t role i n the creation a n d a n i m a t i o n of v i r t u a l worlds [40]. Forces are a very n a t u r a l way for a designer to interact w i t h v i r t u a l objects. R e a l i s t i c force responses are i n t u i t i v e because they copy the behavior of our s u r r o u n d i n g physical w o r l d . A realistic force response is a powerful m e t h o d to manipulate complex surfaces while meeting other constraints. A n example is the editing of a m o d e l of a h u m a n 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 i n accordance w i t h physical limits [126]. A n o t h e r example is the Virtual  Clay  system of M c D o n n e l l , Q i n and W l o d a r c z y k [79] where forces d u r i n g s c u l p t i n g are simulated i n real time allowing the creation of v i r t u a l sculptures. M a n y computer graphics techniques have evolved over the years to m o d e l deformations [40]. These techniques vary from full c o n t i n u u m models for elastic materials, through mass-spring particle systems, to very ad-hoc methods. B u t even the most elaborate deformation m o d e l w i l l be unrealistic if the simulated material does not m a t c h the real physical object. T h e measurement of these m a t e r i a l properties has so far only been done i n a "material oriented" way, e.g., i n physical and mechanical testing laboratories. A s a,result, m a n y computer graphics methods do not m o d e l material properties for deformation simulation of physical objects. F u r -  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 computational model and a robotic measurement system is required. This combination creates a scanner for deformation behavior of objects which is the goal of this thesis. 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 information 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 objects 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 observations. (Given a parametric model that applies, the process of determining model parameters forms the field of parameter estimation.) Inverse problems related to deformable 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  s t r u c t u r a l analysis i n c i v i l engineering). Surprisingly, the verification and validation of c o m p u t a t i o n a l methods i n inverse problems i n mechanics are often neglected as well. T h e performance of c o m p u t a t i o n a l methods for inverse problems is c o m m o n l y evaluated by comparison w i t h either an analytic solution or by comparison w i t h the same m e t h o d using finer discretization of the m o d e l d o m a i n (or b o u n d a r y ) . D i r e c t experimental validation provides a more meaningful characterization of a m e t h o d . T h e scanning m e t h o d of this thesis makes experimental v a l i d a t i o n 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. T h e purpose of the deformable m o d e l i n computer graphics is to allow physically realistic rendering of the deform a t i o n of a real object under force. T h e m o d e l has to encode the response of the object as a whole to a force. W h a t interior structure causes t h a t response is irrelevant. However, the scanning of object deformation is also the first step i n an inverse p r o b l e m solution. T h i s thesis also contains a . m e t h o d to estimate physical model parameters based on scanned behavior. T h e model assumption i n this case is that the object behaves as a homogeneous linear isotropic elastostatic body. N e x t , in Section 1.1, a concise p r o b l e m statement is given w h i c h may guide the reader through the following chapters. Section 1.2 reviews work related to this thesis.  4  1.1 1.1.1  Problem Statement and  Overview  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 estimation 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 ( A C M E ) for object deformation measurements. A C M E is a robotic facility which automates object measurement tasks. Object deformation is characterized by surface displacements 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 deformations which an object 1  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 Here, 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. 1  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 measurement 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 chapters 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  Related  Work  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 employed 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, Section 1.3 gives a brief summary of related work. 1.2.1  Deformable Models  This section focuses on deformation models for physically realistic interactive simulation. At the Forschungszentrum Karlsruhe in Germany, the development of medical 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 E T H 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 MooneyRivlin strain energy function W = /i(J\ — 3) + a(J — 3) + 2  J3 — l ) where J\, J 2  2  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 linear 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 California are computed with finite ele2  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 The project is a joint cooperation between the University of California, San Francisco, Berkeley and Santa Barbara. 2  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 simulation 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 P r o b l e m s i n E l a s t i c i t y  Inverse problems aim to discover the cause of a physical phenomenon from its manifestation. 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 coefficients 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 estimation. 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 testing. 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]), composites (e.g., [46]), civil engineering (e.g., [4]) and earth sciences (e.g., [41]). 1.2.3  T h e o r e t i c a l 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 functions 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 dimensions [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, inhomogeneous 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 observations require a more elaborate experimental setup than observations of the boundary alone. Parameter estimation using domain methods is applied to beam structures [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 measurements are available, the problem is more difficult. The boundary element method has been applied in three-dimensions even if interior measurements are not available [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 of R e l a t e d W o r k  The trend in interactive simulation of deformable models has been towards physically more accurate models. Finite element solutions to problems with physical constitutive relations of the continuum are often employed. The computational requirements 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 elastic 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 behavior developed in this thesis represents an alternative approach to acquire realistic deformation models.  14  Chapter 2  Elastostatic Object Behavior Deformation of solids has been studied extensively in physics, civil engineering, mechanical 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 employs. 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  Discrete Green's Functions  Matrix  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 U T i =T and To n F\ = 0. 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 = [ivf,..., u ^ ]  T  and p = [pj,...,  p], T  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 face To, while  =  =  if the vertex k is on the sur-  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 deformed 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 traction 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 Section 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  Deformation Models in C o n t i n u u m  Mechanics  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  Strain Tensor  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\  1 +j =  du2 dx\  dx2 i  du3 dxi  i  9x3  du2 dx2  du2 dx%  du3 0x2  i  ,  du 9x3 3  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) = (dx) dx apart if the body 2  T  is undeformed. In the deformed state, they are a distance (SI) = (dy) dy apart. 2  T  Since y = x + u it follows that dy = dx + du/dx dx = (I + du/dx)dx. However  18  the J a c o b i a n is defined as J = I + J = / + <9u/<9x and it follows that the change i n squared distance between point P and Q is  (dl)  2  - (5l )  2  0  =  (dyfdy  - ( d x ) d x = (dx) J Jdx  =  (dx) [(J  T  -  (dx) 2edx  T  T  +J)  T  - (dx) cix T  +J J]dx r  3  3  = YlYl  T  T  J  d x i 2 £ i  d x  J  (-) 2  2  i=ij=i  E q u a t i o n 2.2 defines the s y m m e t r i c strain tensor Sij and because of its s y m m e t r y m a y be w r i t t e n as a strain vector e = [e\ i £2.2 £3,3 £1,2 £1,3 2 , 3 ] ' • E q u a t i o n £  T  t  2.2 can be expressed as e = D u where D is a differential operator. I n engineering applications, the strain is usually assumed to be small, i.e., the J a c o b i a n p r o d u c t JJ T  0 and the strain tensor simplifies e  ~ \\J  + J ] . T h e differential operator  T  m  D reduces then to three tensile and three shear strains.  2 ^ 0  0  0  24-  0  OX2  0  0  d  d  0x3  9x2  9x\  0  9 9x3  •0  9 3xi  a  0  9x  3  (2.3)  9 9x2  T h e Green-Lagrange strain tensor of E q u a t i o n (2.2) is independent  of the  p a t h of deformation, which is adequate for linear elasticity. N o t e however, that if rotational components are significant, J  T  J Z$> 0, the a p p r o x i m a t i o n no longer holds.  T h e a p p r o x i m a t i o n is also not suitable for plastic deformation w h i c h depend on the p a t h of deformation.  Rate-of-deformation  cases [135].  19  and spin tensor are employed i n these  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 system 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 surface 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 55 = nrfs because of the 2  alignment with the respective coordinate axis. Newton's second law for the infinitesimal tetrahedron (i.e., its volume and therefore forces due to its mass tend to zero) becomes f! + f + 2  f3  = f.  This yields Cauchy's stress formula (Eqn. 2.4) and the stress tensor Oij.  Pj = arij = o~ijrij 20  (2.4)  T h e stress tensor is s y m m e t r i c which allows for a vector n o t a t i o n a  °~2,2  03,3  <7i 2  °2,3] T  t  =  T h e C a u c h y stress tensor does not incorporate a stress  h i s t o r y nor does it allow for l a r g e  deformations. O t h e r definitions of stress tensors  1  in use are the second P i o l a - K i r c h h o f f pseudo stress tensor, the J a u m a n n stress rate a n d the G r e e n - N a g h d i rate [135]  w h i c h are less restrictive t h a n the C a u c h y stress  tensor. T h e stress a n d s t r a i n tensor are related by the constitutive equations for a m a t e r i a l w h i c h are discussed next.  2.2.3  Constitutive Equation  C o n s t i t u t i v e equations relate stresses to strains a n d are m a t e r i a l dependent.  The  constitutive equation of linear-elastic m a t e r i a l is given by the generalized Hooke's law: 3  3  (2.5)  k=\ 1=1  E q u a t i o n 2.5 contains a m a x i m u m of 36 independent m a t e r i a l parameters due to the s y m m e t r y of stress a n d s t r a i n tensor.  However, for homogeneous isotropic  materials the m a t r i x of elastic coefficients has o n l y two independent variables; e.g., the  Lame  constants  3X)/(fi + A) a n d the  A a n d /j,, or alternatively the Young's m o d u l u s E — yi(2p +  and the  Poisson's ratio v = A/(2(/i + A)),  Poisson's ratio v\  or the  the latter is used i n this thesis.  Shear m o d u l u s Equation  2.6 is the  tensor n o t a t i o n 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  expresses the spring constant i n m a t r i x notation (a. = C  Cijki =  z  e).  7z-Sij5ki + 7(8ikO~ji + Sudjk)  Here, large implies geometric non-linearities can not be ignored.  21  7  (2.6)  2.7  = /i.  2 (l-i,)  271/  7  \-1v  \-2u  c =  2.2.4  0  0  0  2 (l-«/)  2-yv \-2w  0  0  0  2-yu  7  l-2u  271/ \-2v  \-2v  2yu \-2v  l-2i/  2 (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  7  (2.7)  Navier's Equations  In general, the forces acting on an infinitesimal material cubic element under deformation 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  2.3  +  du  7  k  2v  dxkdxj  + bi =  0  (2.9)  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 d o m a i n of the solid can be numerically evaluated w i t h the F i n i t e Element M e t h o d ( F E M ) . T h e result of the F E M is a block m a t r i x for the solid's stiffness w h i c h relates the applied forces to the solid's deformation. T h e m e t h o d of condensation yields a discrete Green's functions m a t r i x from this stiffness m a t r i x . T h i s is briefly o u t l i n e d i n Section 2.3.1. A p p l i c a t i o n of Green's theorem to the elastic e q u i l i b r i u m equation leads to a b o u n d a r y integral equation, w h i c h can be evaluated by a weighted residual statement. T h e discretization of the residual statement is the starting point of the direct B o u n d a r y Element M e t h o d . T h e B E M solution directly leads to the m a t r i x of Green's functions. T h i s approach is described i n some detail i n Section 2.3.2.  2.3.1  Finite Element Model  T h e v i r t u a l work of elastic forces can be derived from N e w t o n ' s second law. I n the following b o d y forces are ignored for clarity but they could be included trivially. E q u a t i o n 2.8 relates the forces on an element of an elastic solid. Integration of the equation of e q u i l i b r i u m over the d o m a i n  of the deformable b o d y yields  E m p l o y i n g the C a u c h y stress formula and divergence theorem one can arrive (see [118, pp. 247-249]) at  T h e v i r t u a l work of the elastic forces is the volume integral of the product of the stress tensor and the strain tensor.  5W 23  =  /  e C dedQ  =  / (Du) CDdudft  T  T  (2.10)  T  E q u a t i o n 2.10 can be made discrete by e m p l o y i n g finite volume elements throughout the d o m a i n Q of the solid. T h e displacements and forces are approxi m a t e d w i t h p o l y n o m i a l t r i a l functions over the elements.  O n e choice of a finite  element is a tetrahedral element w i t h four hat-functions centered at the vertices of the tetrahedron (see F i g u r e 2.3). Note that the b o u n d a r y of the solid is represented as a triangular mesh and the displacement field a n d force field over the b o u n d a r y is represented w i t h hat-shaped basis functions. T h i s is identical to the B o u n d a r y Element M o d e l described i n Section 2.3.2 ( w i t h the exception that it is customary to employ tractions instead of forces i n the B E M ) .  0  Interior Vertices  O  Exterior Vertices  Q  Exterior Vertices  F i g u r e 2.3: Discrete D o m a i n and B o u n d a r y (2D) (a) F E M (b) B E M  T h e F i n i t e E l e m e n t M o d e l for an elastostatic solid has the form of E q u a tion 2.11. T h e displacement block vector fi is related to the force block vector f by the stiffness block m a t r i x K .  Ku = f  (2.11)  T h e displacement and force block vector can be reordered into block vectors associated w i t h the interior of the solid uo. and f^ and into vectors associated w i t h 24  the boundary up and fr and it follows K ,r r  K , n  u  K ,n r  r  fr  (2.12)  fQ  r  The displacement of the boundary u r may be calculated from by condensation of Equation 2.12 as K r,r  u r = fr - K r . n K ^ f n •  (2.13)  In the scenario of this thesis where internal forces do not occur and internal displacements 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  /  f (p - p )u* dT + / (u -u )p* dY  ^* *dft = u  k  JY  OXj  JU  k  k  k  JTO  1  k  (2.14)  k  Integrating E q u a t i o n 2.14 b y parts twice, as well as e m p l o y i n g the C a u c h y stress formula 2.4 yields  In  dxj  udtt = - / p u* dT - / p u* dT + / u p dT + / u p* dT. Jr-i  k  k  k  Jr  k  k  k  JTi  0  G i v e n t h e fundamental solution for da* /dxj k  k  Jr  k  0  + A = 0 where A is an u n i t impulse,  we arrive at the b o u n d a r y integral equation 2.15.  l u* ^, ) (x)dr( ) v  j  x  Pl  = c, (0",(0 +  x  V  J PUZ,x)u (x)dr( ) r  3  (2.15)  x  E q u a t i o n 2.15 is the b o u n d a r y integral equation w r i t t e n as a weighted residu a l statement over the complete surface T of the elastic body. T h e l o c a t i o n of the load a p p l i c a t i o n is £, while the response location is x  a  n  d the distance between  t h e m is r. C o o r d i n a t e i is associated w i t h t h e 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  K e l v i n ' s fundamental solution for the 3-D N a v i e r ' s equations for a n infinite b o d y w i t h the same m a t e r i a l properties subjected to an impulse point load.  T h e fun-  damental solution is E q u a t i o n 2.16 a n d E q u a t i o n 2.17 where the surface n o r m a l is denoted by  n a n d 5ij  is the Kronecker delta.  dr dr (3 - Au)5ij + dxi dxj  1 167r(l — u)jr Pijfox)  -1 8TT(1  -  dr dr \ dr (1 - 2v)8ij + 3 v)r dx\ dxj J dn dr dr - ^-dx- '  (2.16)  2  { l  26  2 u ) {  n  ]  (2.17)  M o r e details on the derivation of the b o u n d a r y integral equation can be found i n m a n y textbooks i n c l u d i n g [16, 33, 17].  T h e integrals i n E q u a t i o n 2.15  can be numerically evaluated over linear triangular b o u n d a r y elements w i t h Gaussian quadrature.  In the following the t r a c t i o n vector field a n d the displacement  fundamental solution tensor are w r i t t e n as p * ( £ , x)  a  n  d u * ( £ , x) respectively. T h e  left-hand side of E q u a t i o n 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 E q u a t i o n 2.15 at a mesh node k is  (Hu)  =  f c  ^p*(&,x)u(x)dT(x).  T h e t r a c t i o n vector field p ( £ ) over the b o u n d a r y is a p p r o x i m a t e d w i t h shape functions  p(x) = JliLi <f>i{x)Pi  3 x 3 block elements G f c  m  where p^ are the traction values at a node. I n d i v i d u a l  of the fundamental solution m a t r i x are  Gkm =  J u*^ , )4>m(x)dT(x)k x  I n d i v i d u a l 3 x 3 block elements H f c  H  f  c  m  =  m  (2-18)  of the fundamental solution m a t r i x are  jV(£ ,x)<M^r(x). fc  (2.19)  T h e integrals are n u m e r i c a l l y evaluated by G a u s s i a n quadrature. T h e comp u t a t i o n 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. T h e discretized b o u n d a r y integral equation is the m a t r i x equation H u = G p . T h e equation has to be rearranged for the specific b o u n d a r y condition at h a n d . T h e prescribed t r a c t i o n and displacements are collected i n the block vector v while the right-hand side is the complementary set of tractions and displacements v . E x c h a n g i n g elements i n vectors p and u to yield 27  v and v requires exchanging block columns i n matrices G and H to yield G and H . T h e modified equation H v = G v is solved w i t h the Green's functions m a t r i x S as follows:  G^Hv  = S v = v.  T h i s completes the derivation of the Green's functions m a t r i x .  N e x t ap-  proximate physics-based models are compared to the formulation of this thesis of E q u a t i o n 2.1.  2.4  Other Physics Based Deformable  Models  P h y s i c s based deformable models have been employed extensively i n the computer vision, computer graphics and medical i m a g i n g communities. physics based m o d e l i n g is not employed consistently.  However, the t e r m  In this Section, we briefly  describe the central equations i n two of the most c o m m o n approximate models. I n Section 2.4.1 the set-up of particle systems is described, while i n Section 2.4.2 the m o d e l of deformation i n deformable superquadrics is s u m m a r i z e d .  2.4.1  Particle Systems  T h e m o d e l i n g of elastic solids by particle systems has been applied i n computer graphics since 1981 [100]. T h e particle system approach divides the solid into i n d i v i d u a l mass points inter-connected by springs (see F i g u r e 2.4). U s u a l l y the mass points are arranged i n a regular lattice b o r r o w i n g from the atomic structure of crystals. T h e springs are t y p i c a l l y linear and usually a viscous d a m p i n g term is added. A linear elastostatic solid modeled as a particle system gives rise to an equat i o n K u = f [40].  T h e stiffness m a t r i x K is composed of the i n d i v i d u a l spring  28  F i g u r e 2.4: M a s s - S p r i n g or P a r t i c l e 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. T h e stiffness m a t r i x is banded determined by the connectivity of the lattice. Various researchers attempt to set the s p r i n g constants to physical properties. Lee, Terzopoulos a n d Waters set the spring constants of their facial models according to typical h u m a n anatomy [73]. T h e actual constants are a r b i t r a r y while the m a t e r i a l d i s t r i b u t i o n (e.g., muscle, facial tissue, epidermis etc.) is determined by the t y p i c a l h u m a n anatomy.  K o c h et al. set the spring constants according to the m a t e r i a l  d i s t r i b u t i o n of a i n d i v i d u a l h u m a n face based on a C T scan [65]. A g a i n the actual constants are arbitrary.  R e z n i k and L a u g i e r tune the spring i n their m o d e l of a  rubber t i p of a r o b o t i c finger according to the Y o u n g ' s m o d u l u s of rubber [106]. Identification from observation for d y n a m i c systems is described i n [60]. In [77] a mass-spring cloth m o d e l is fit to observation of hanging cloth behavior. D ' A u l i g n a c , Cavusoglu and L a u g i e r [27] m o d e l the elastic behavior of the h u m a n t h i g h w i t h linear and non-linear springs from observations.  T h e observations consist of force  and displacement measurements i n response to pressure p r o b i n g w i t h a robotic a r m . A general difficulty w i t h mass-spring models is the lack of a suitable generalpurpose t u n i n g a l g o r i t h m for the springs [40]. T h e springs have no physical connection to the behavior of an elastic solid. M a s s - s p r i n g models do not m o d e l various  29  physical behavior of solids well including incompressibility and large stiffness. However, 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 membrane under tension and a thin plate under tension [124]. Metaxas and Terzopoulos extended this work into deformable superquadrics and introduced it to the computer 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 Terzopoulos 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.  A l l 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 t h this thesis:  2.5  Summary  T h i s chapter details the deformation m o d e l employed by this thesis.  The model  is a discrete c o n t i n u u m mechanics m o d e l of linear elastostatic object behavior, expressed i n terms of a Green's functions m a t r i x w h i c h can be derived by the B o u n d ary E l e m e n t M e t h o d for homogeneous isotropic linear-elastic material. However, the Green's functions m a t r i x representation does not depend on the B E M ; i n particular, F E M w i t h condensation leads to similar formulation. T h e definitions and assumptions i n c o n t i n u u m mechanics for a linear elastostatic solid have been s u m m a r i z e d i n this chapter. L i m i t a t i o n s of the m o d e l are due to geometric a p p r o x i m a t i o n , i.e., first-order  deformation and m a t e r i a l a p p r o x i m a t i o n , i.e., linear elasticity.  T h e remainder of this thesis utilizes E q u a t i o n 2.1 i n the estimation task. T h e estimation task is finding the discrete Green's functions m a t r i x based on observation of vertex displacements and tractions at a vertex.  31  Chapter 3  The A C M E  Facility  T h e A c t i v e Measurement F a c i l i t y ( A C M E ) is an integrated robotic facility designed to acquire measurements of interactions w i t h objects.  A t the core of A C M E is a  contact m a n i p u l a t i o n system ( C M S ) w h i c h executes contact interactions w i t h an object to be measured.  These interactions can be recorded w i t h various sensors  i n c l u d i n g a force-torque sensor, a t r i n o c u l a r stereo system, a microphone and a high quality 3 - C C D color video camera.  Sensors w h i c h are employed to measure  the  response of an object from a distance are combined i n a field measurement system (FMS). A C M E is implemented as a robotic facility because of the a b i l i t y of robots to perform tasks w i t h high repeatability and consistent accuracy. T h i s is extremely i m p o r t a n t if a large number of registered measurements  of an object are to be  collected. T h e test station is a p o s i t i o n i n g robot to move the object to be measured w i t h precision. T h e C M S is b u i l t a r o u n d a robot a r m positioned b y a m o t i o n 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 p r o g r a m m i n g are unique to A C M E . A s a result A C M E is a tele-robotic system w i t h a t o t a l of fifteen degrees of freedom.  32  F i g u r e 3.1: A C M E Facility O v e r v i e w  T h i s thesis employs A C M E for the acquisition of measurements d u r i n g active object deformation. W h i l e the design and implementation of A C M E is a group effort, some features are especially designed to enable the acquisition a n d verification of deformable models. These feature are discussed i n some detail i n this chapter, while o n l y an overview is provided a b o u t 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 i n c l u d i n g the actuators a n d sensors of A C M E (Section 3.1.1), the control architecture (Section 3.1.2), the p r o g r a m m i n g of A C M E (Section 3.1.3) and the sensor p r o g r a m m i n g m o d e l (Section 3.1.4). T h e structure of a n A C M E experiment d i r e c t i n g A C M E to actively deform an object and record the deformation is discussed i n Section 3.2. T h e geometric m o d e l i n g c a p a b i l i t y of A C M E is s u m m a r i z e d w i t h examples utilized in this thesis in Section 3.3.  33  3.1  Overview  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 measurements. 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 devices) 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 controllers when possible but adding advanced features as required. In particular, the control architecture provides the foundation for the high-level control of the system 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 telerobotic 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 Java . On the other hand the simulator is available during 1  execution to verify motion requests. This allows for a simple generate-and-test approach 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 A C M E . 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 A C M E . 1  Java is a trademark of Sun Microsystems Inc., MountainView, Ca., USA  34  F i g u r e 3.2: M e a s u r e m e n t of contact force a n d p o s i t i o n w i t h C M S  3.1.1 The  A C M E Devices  actuators a n d sensors w i t h i n A C M E are collectively referred to as devices.  T h e y are conceptually d i v i d e d into the following subsystems (see F i g u r e 3.1 for a n overview):  • The  Contact Manipulation System  ( C M S ) is based o n a P u m a 260 robot a r m  w i t h a wrist m o u n t e d 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 a p p l y sustained contact forces of a least ION t o  15N.  T h e A T I M i n i has a m a x i m u m force limit of 40/V to 120vV (depending on the axis) and a m a x i m u m torque of 2Nm and  l.ONmm,  w i t h a resolution of 0.047V to 0.12A  respectively. A 1.267m(4/t.) long linear  r  motion stage  the robot a r m to a n d from its application area a n d increases the robot arm's workspace.  T h e control of the robot a r m has been modified to enable the  application of forces to a test object. T h e robot a r m can w o r k i n impedance control mode or execute guarded moves m o n i t o r i n g i m p a c t forces at its t i p . T h e C M S is capable of measuring object properties such as friction a n d local stiffness a n d of a p p l y i n g a force impulse to generate contact sounds. F i g u r e 3.2 shows the C M S a c q u i r i n g force d a t a while i n contact w i t h a test object.  35  moves  F i g u r e 3.3: A C M E F M S S u b - S y s t e m  • The  Field Measurement System  ( F M S ) holds sensors for measuring the light  a n d sound fields a r o u n d the test object (see F i g u r e 3.3). T h e u n i t is positioned w i t h a gantry r o b o t ( a l u m i n u m s t r u c t u r a l elements by "Item Industrietechnik u n d M a s c h i n e n b a u m b H " powered by P a r k e r H a n n i f i n C o r p . brushless D C motors) a n d oriented by a two degree-of-freedom pan-and-tilt u n i t (Sagebrush Model-20).  T w o vision sensors are part of the sub-system:  a high quality  3 - C C D color video camera w i t h computer controlled z o o m a n d focus (Sony D X C 950) and a t r i n o c u l a r color stereo vision system (Point G r e y Research C o l o r T r i c l o p s ) . A microphone is also attached to the unit.  • The 3-DOF  Test Station  presents the test object t o the sensing equipment.  T h e test-object c a n be m o u n t e d rigidly to the test station w h i c h positions the object i n the plane. T h e test-station consists of x - a n d y - t r a n s l a t i o n a l m o t i o n stages, as well as a rotation stage. T h e stages (by Daedal) possess a linear a n d rotary resolution of ± 0 . 0 0 6 3 5 m m and ± 1 0 a r c - m i n , respectively.  • Miscellaneous devices include a time-of-flight laser range finder ( A c u i t y A c c u Range 3000 L I R ) a n d various lights which can be switched under computer control.  36  A C M E Client  Robot A r m Force Sensor Motion Stage  F i g u r e 3.4: A C M E Server: O v e r v i e w of C o n t r o l  3.1.2  Control Architecture  T h e software a n d hardware control architecture of A C M E is briefly o u t l i n e d i n this Section.  A n overview of the architecture is shown i n F i g u r e 3.4.  T h e control is  client-server based, i.e., an A C M E client w i t h w h i c h a user interacts a n d an A C M E server w h i c h controls the robotic subsystems i n t r o d u c e d above. T h e A C M E server system is d i s t r i b u t e d between several computers; separating low-level closed-loop control, trajectory generation, d a t a acquisition, a n d networking. The  four layers of the A C M E server software from top to b o t t o m 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 A C M E . The client is the program development environment for A C M E and is composed of a three-dimensional  38  View  F i g u r e 3.5: A C M E 3-D D i s p l a y A collision event p r o d u c e d by the simulation is rendered w i t h the c o l l i d i n g boxes highlighted.  display of the facility, a simple editor and access to the J a v a c o m p i l e r . T h e client has access to the complete A C M E J a v a code which includes a k i n e m a t i c simulator of the facility.  T h i s allows the user to develop, debug and verify a l l  w i t h o u t access to the A C M E server.  Once a compiled  Experiments  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 Internet, uploads the c o m p i l e d an  Experiment  Experiment  and receives any d a t a requested. Because  is any J a v a object w h i c h adheres to the  Experiment  user has complete freedom to implement any desired routines.  interface, the  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 r u l i n g out c o m m u n i c a t i o n delays between server a n d client. D a t a processing routines can also be r u n on the server w h i c h allows d a t a to be reduced to its essential size before being t r a n s m i t t e d to the client. The  simulator models the entire A C M E facility w i t h an enhanced J a v a 3 D  39  scene graph. T h e scene graph contains the geometric, positional, a n d appearance information used i n s i m u l a t i o n on the client, run-time checking on the server and graphic display on either server or client. A n i m p o r t a n t feature of this scene graph is a simple on-line collision detection using oriented b o u n d i n g boxes [42]. B o u n d ing boxes are placed a r o u n d p r i n c i p a l system components a n d can be tailored for different experiments. F i g u r e 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 d a t a from sensors. T h e p r o g r a m m i n g abstraction of these sensors is an i m p o r t a n t design feature of A C M E . T h e highlevel description of sensors is related to the description of actuators discussed i n Section 3.1.1. Sensors i n A C M E are widely d i s t r i b u t e d on the computers that are the A C M E server. A l l actuators provide a position sensor functionality w h i c h a l lows the recording of an actuator's position at coupled w i t h the R C C L layer of the server.  100Hz. T h i s functionality is tightly T h i s is also true of the force sensor  of the robotic a r m . T h e video digitizer of the color 3 - C C D camera a n d its serial c o m m a n d interface on the other h a n d are quite separate from A C M E server. T h e same is true for t r i n o c u l a r stereo-head, the laser range-finder a n d the microphone. Nevertheless, a l l these devices are A C M E  Sensor objects.  another client-server layer. A n i n d i v i d u a l  T h i s is achieved by a d d i n g  SensorServer  executes on the computer  hosting the sensor hardware, while the A C M E server acts as a client to these i n dividual ACME  SensorServers. Experiment  Experiment  T h i s design not only unifies the access to sensors from an  but it also enables asynchronous real-time d a t a acquisition. A n  specifies  SensorTrigger  40  events w h i c h are uploaded a n d checked on the  server of the sensor. After being started, the server of the sensor continuously acquires d a t a checks the trigger c o n d i t i o n and buffers d a t a to be sent to the A C M E server and eventually to the A C M E client. In this fashion, c o m m u n i c a t i o n delays are eliminated from the d a t a acquisition loop. For example, the itizer, the  Triclops  server incorporates the control of the video dig-  Point Grey Research Color Triclops  Point Grey Research Triclops  control software and access to the  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 m o n i t o r i n g depth in the stereo imagery, or implement whatever image processing is required. A s the trigger activates, the  SensorServer m a y  proceed to send raw images or images w h i c h  are the result of processing i n 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 w i t h the C M S while m o n i t o r i n g the global deformation w i t h the F M S . T h i s section details the methods utilized i n p l a n n i n g the contact w i t h the object, executing the p r o b i n g a n d a c q u i r i n g measurement data. T h e structure of the A C M E  DeformExperiment  is described here i n detail since it is the  Experiment used  i n this  thesis to acquire measurements a n d it also serves as an example of an A C M E Ex-  periment.  The  Experiment  does assume that a geometric m o d e l of the object under  test is available. T h e next Section 3.3 describes how the models for this thesis are acquired. A l s o not part of the  DeformExperiment  are the measurement processing  routines. T h i s separation into d a t a acquisition a n d analysis was convenient for the research methods explored i n this thesis.  41  T h e processing of measurements is de-  Experiment A  TestObject DeformExperiment  FmsPositioning CmsPositioning  *  ScanPoistion TsPositioning Pokelt  PokeGoalPlanner  SinglePoke  GuardedMoveinZ  extends  Class  (  Interface  ) \ Abstract Class"  F i g u r e 3.6: Structure of an A C M E D e f o r m a t i o n E x p e r i m e n t  scribed i n the following chapters, C h a p t e r 4 to C h a p t e r 6. T h e overall structure of a deformation experiment is depicted i n F i g u r e 3.6. The  DeformExperiment  takes the specification probing, w h i c h is object-  specific, from a data-base object w h i c h implements fication  TestObject.  E x a m p l e s of speci-  include the m a x i m u m force exerted on the object while probing, the max-  i m u m p r o b i n g depth, the number of steps for the deformation measurement (see C h a p t e r 4) a n d various files a n d directories to store d a t a i n . T h e other data-base object implements the  ScanPosition interface  w h i c h contains pre-planned viewpoints  and calibration information for the trinocular stereo-head w i t h i n the F M S . B a s e d on these specifications the  PokeGoalPlanner and DeformExperiment  DeformExperiment  positions the F M S , initializes the  starts the probing process. A simplified execution flow of the  is shown i n F i g u r e 3.6.  T h e p r o b i n g process executes as follows.  42  The  PokeGoalPlanner  generates  1. P o s i t i o n  TestStation, FMS a n d CMS  2. Initialize  PokeGoalPlanner  3. L o o p over a l l positions specified i n  ScanPosition  • L o o p over a l l specified vertices to be p r o b e d — Request probe specification from planner — W h i l e p r o b i n g not successful * Simulate p r o b i n g * Execute * Execute  GuardedMovelnZ SinglePoke  * If not successful change configuration — Process and save a l l d a t a 4. M o v e FMS and CMS to park p o s i t i o n F i g u r e 3.7: Simplified E x e c u t i o n F l o w of A C M E Deformation E x p e r i m e n t  the next p o s i t i o n and orientation for the probe based on the mesh of the object under test and the specification i n  TestObject.  It checks l i m i t s i n order to rule out  goals w h i c h w o u l d b r i n g the r o b o t - a r m too close to the  TestStation.  It also does ray-  intersection w i t h the object i n order to verify that the probe is contacting the object at the expected location. T h e n the  DeformExperiment  calls the s i m u l a t i o n w i t h the  necessary m o t i o n to reach the goal i n a default configuration. O n success it executes this m o t i o n . O n failure the configuration is changed w i t h a heuristic, either changing the orientation of the object by r o t a t i n g the  TestStation  or modifying the distance  of the robot from the object by m o v i n g it w i t h the long linear stage. A successfully simulated m o t i o n is executed i n a three-step process. F i r s t an approximate position away from the object is reached. Second, a  GuardedMovelnZ is executed  to find the  actual position of the surface of the object. A l l measurements are made relative to this position reducing the influence of p o s i t i o n i n g errors of the robot a r m . 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 r u n i n s i m u l a t i o n w h e n it is applied to an object for the first time. T h e s i m u l a t i o n performs collision checking and the rendering of the motion provides a visual confirmation of the experiment. If this performs satisfactory, A C M E acquires the deformation d a t a autonomously.  3.3  Geometric Modeling with  ACME  G e o m e t r i c m o d e l i n g of objects i n A C M E is necessary for various tasks i n c l u d i n g contact p l a n n i n g , object m o d e l i n g a n d 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 m a t r i x described i n Section 2.1.  A surface mesh is also required for the  deformation d a t a acquisition described i n Section 3.2 above. F u r t h e r m o r e , the surface mesh is also utilized i n the deformation measurement described i n Section 4. F i n a l l y , the estimation of the discrete Green's functions m a t r i x i n C h a p t e r 5 a n d C h a p t e r 6 relies on the discrete surface representation as well.  It is evident that  geometry is essential to this thesis. Despite this, no novel m e t h o d for i n i t i a l m o d e l acquisition is proposed. Unfortunately, geometry acquisition techniques are still far from being fully automated and general. D u r i n g the research for this thesis distance measurement w i t h a time-of-flight laser range-finder has been investigated [70], as well as a novel stereo-vision m a t c h i n g technique [69]. In the end, a collection of techniques based on d a t a obtained from the t r i n o c u l a r stereo-head have been found to be effective. T h e techniques succeed i n o b t a i n i n g a closed (or water-tight) mesh w i t h s u b d i v i s i o n hierarchy (see F i g u r e 3.8) for the objects i n this thesis. T h e process is only semi-automated a n d some mesh e d i t i n g 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  Summary  T h e A C M E facility is a robotic measurement facility w h i c h enables the execution of experiments i n which a probe contacts an object. T h e probe is m o u n t e d on a robotic arm equipped w i t h a force sensor. T h i s contact m a n i p u l a t i o n system ( C M S ) is able to exert a force onto the surface of an object at a specified location. T h e t r i n o c u l a r stereo-head w i t h i n the field measurement system ( F M S ) can be placed to observe the object d u r i n g deformation from various viewpoints. A robot such as A C M E can execute repetitive tasks w i t h high repeatability which is i m p o r t a n t d u r i n g exhaustive contact s a m p l i n g of an object's surface.  T h e user interface of A C M E combined  w i t h its tele-robotic design make p r o g r a m m i n g contact interactions convenient. A further feature of A C M E necessary for this thesis is the ability to produce threedimensional point-clouds registered i n an object reference frame. These point-clouds are combined i n a mesh generation step w h i c h produces a hierarchy of meshes w i t h s u b d i v i s i o n connectivity w h i c h 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  F i g u r e 3.8: M e s h R e c o n s t r u c t i o n T h e raw scan mesh 3.8(a) w i t h a total of 15076 vertices(v) and 29095 faces(f) i n c l u d i n g various external stray polygons. T h e s u b d i v i s i o n hierarchy (3.8(b), 3.8(c), 3.8(d) and 3.8(e)) is generated using displaced s u b d i v i s i o n surface style piercing of the raw scan mesh. 46  Chapter 4  Deformation Measurement T h e deformation measurement  task is to estimate how the vertices of the unde-  formed triangular mesh are displaced d u r i n g a p p l i c a t i o n 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 i n 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 ) . In 1  mechanical engineering, the strains considered are u s u a l l y s m a l l and often o n l y 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. F o r example, i n s t r u c t u r a l mechanics, strain gauges give local information at critical points of the structure. T h 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  information is recovered.  See http://www.gom.com  l  47  task at h a n d since three-dimensional  T h e task i n range-flow is to m a t c h depth information over time. T h e rangeflow constraint equation [51, 132] for the depth of a surface Z(x,y,t) is  dZ ox  ,  ,„  dZ ,  dZ ,  n  ——ax + —— ay — dZ + -7—at = 0.  oy  at  T h e range-flow constraint equation is the three-dimensional equivalent of the o p t i c a l flow equation [52].  S i m i l a r to o p t i c a l flow [6] this constraint m a y be utilized by  area-based m a t c h i n g or directly w i t h a differential technique. T h i s thesis employs area-based m a t c h i n g . E s t i m a t i n g range-flow is a harder p r o b l e m t h a n o p t i c a l flow because of the increased dimensionality a n d the necessity for quality range data. Nevertheless, there has been progress i n 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 c o m b i n a t i o n of simultaneous 2 D - o p t i c a l flow w i t h stereo-depth information. Simultaneous stereo depth and o p t i c a l flow goes back to at least B a l l a r d and K i m b a l l [3]. R e c e n t l y this approach has been followed successfully i n [128, 133, 119]. However, i n this thesis, the a d d i t i o n a l imagery is utilized to reduce the error i n the very noise-sensitive optical  flow.  T h i s results i n a purely local m e t h o d and no global regularization  is employed. T h e high quality c a l i b r a t i o n of the c o m m e r c i a l t r i n o c u l a r stereo syst e m [105] i n A C M E lets this approach succeed. The  choice of stereo vision and o p t i c a l flow requires objects w i t h sufficient  visual texture.  Sufficient texture is texture w h i c h leads to an unique peak i n the  area-based correlation at the correct d i s p a r i t y i n the stereo m a t c h i n g and at the correct location i n 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 t h e m trivially. In Section 4.1 methods to robustly estimate flow are briefly reviewed and i n 48  Section 4.1.1 the approach of this thesis to robust range-flow estimation is o u t l i n e d . T h e c o m p u t a t i o n s involved i n this m e t h o d are detailed i n Section 4.1.2. R e s u l t s for rigid b o d y m o t i o n compared to g r o u n d t r u t h d a t a are given i n Section 4.1.3, while examples for deformable objects are shown i n Section 4.1.4. Range-flow results are not direct measurements of the vertex displacement u^, needed for the estimation task defined i n Section 5.1. These displacements need to be estimated based on range-flow results. T h i s estimation is described i n Section 4.2 and results are given i n Section 4.2.1. A s u m m a r y i n Section 4.3 concludes this chapter.  4.1  Robust  Flow  Estimation  T h i s section gives a selective review of range-flow a n d optical flow e s t i m a t i o n i n light of the deformation measurement task at h a n d . M o s t optical flow techniques can be characterized either as a differential technique or as a region-based m a t c h i n g technique [6]. Differential techniques use the o p t i c a l flow constraint equation w h i c h is underdetermined (one equation i n two u n k n o w n s ) . L o c a l averaging or a regularization t e r m need to be introduced i n 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 m a t c h i n g techniques, on the other h a n d , i m p l i c i t l y average the flow w i t h i n the area being matched a n d no further constraints are s t r i c t l y necessary. Nevertheless, the results t y p i c a l l y are noisy and need further processing [6]. T h i s processing m a y i n clude i n t r o d u c i n g a p r i o r i constraints on the solution v i a regularization (e.g., [68]), fitting  results to a parametric model of m o t i o n (e.g., [55]), or filtering flow results  based on a reliability measure (e.g., [2]).  M u l t i p l e approaches are often used i n  combination. Since the displacement to be measured is between the undeformed a n d de49  formed object states, only two images are necessary. However, accurate derivatives can only be calculated from a sequence of images (unless the signal is h i g h l y oversampled or the intensity structure is nearly linear [6]).. L a c k of accurate derivatives is identified a major source of error i n differential approaches [6]. Area-based techniques do not require estimates of the time derivative of intensity a n d depth i n o p t i c a l flow a n d range-flow respectively. Therefore, an area-based optical flow approach is implemented to avoid h a v i n g to deal w i t h a sequence of images. Stereo-range d a t a 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 a d d i t i o n a l imagery are also exploited i n the c o m p u t a t i o n of optical flow in this thesis.  T h e integration of results from a l l stereo cameras allows the  rejection of outliers, especially if they are generated by viewpoint dependent effects or if they are the result of ambiguous m a t c h i n g results.  Specular highlights (see  e.g., [13]) are a viewpoint dependent effect that produces noise i n the flow estimate. A n o t h e r source of noise is the lack of texture i n the matched area.  T h i s can be  identified using m u l t i p l e cameras since the peak i n the correlation surface w i l l vary between cameras.  T h e peak w i l l also be inconsistent at m o t i o n boundaries asso-  ciated w i t h depth discontinuities. T h e peak i n the correlation surface w i l l be i n different locations since the percentage of area i n the m a t c h i n g window s u p p o r t i n g the different motions w i l l vary between cameras. M u l t i p l e cameras also help to assign a confidence measure for optical flow results w h i c h has been identified as a h a r d p r o b l e m [6]. T h e a d d i t i o n a l d a t a from m u l t i p l e views does not prevent the u t i l i z a t i o n of other approaches to compute a robust optical flow result.  I n particular, the ro-  bust optical flow techniques of [13] (multiple motions w i t h i n an aperture) or the  50  regularization method based on the sum of squared differerences (SSD) of [68] (minimizing 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 areabased 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 threedimensional 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 i n i t i a l registration m a y be i m p r o v e d w i t h an iterative closest point ( I C P ) [8] type a l g o r i t h m . However, i n our set-up this is not necessary since the i n i t i a l guess is sufficient for segmenting the range image. A l l image row, c o l u m n , d i s p a r i t y triplets associated w i t h surface pixels are i n p u t into the tracker.  N e x t , o p t i c a l flow from  the i n i t i a l i z a t i o n to the first t i m e step i n the reference camera is calculated for a l l pixels w h i c h are tracked. In the approach of this thesis, area-based m a t c h i n g uses the n o r m a l i z e d cross-correlation. T h e normalized cross-correlation has been selected because of reports of its good performance i n area-based m a t c h i n g [35] and because it can be calculated from correlation scores w h i c h are. linear w i t h respect to image intensity (see E q u a t i o n 4.1). T h e o p t i c a l flow i n the non-reference camera images is calculated for surface pixels i n the coordinate frame of the reference camera. T h i s is achieved by using the sub-pixel d i s p a r i t y information as a look-up table (see Section 4.1.2 for more details). E a c h surface pixel then has an o p t i c a l flow result for each camera. If a l l results are i n accordance, the flow result is accepted and the t h i r d dimension is added based on the d i s p a r i t y and stereo t r i a n g u l a t i o n . A l l points successfully tracked from i n i t i a l i z a t i o n to the first time step are considered surface points for the following time step. Before discussing the i m p l e m e n t a t i o n details,, some t r a c k i n g issues need to be discussed. It is convenient to record the i n i t i a l i z a t i o n followed by two time steps, even though only one is s t r i c t l y necessary.  T h e initialization is done w i t h o u t the  probe touching the object, while i n the image triplet taken after the first time step the probe is just touching the surface and i n the final image triplet the object is completely deformed. T h e separate i n i t i a l i z a t i o n allows the range-based segmentation to succeed i n image areas where the force probe otherwise w o u l d be very close to the surface.  T h e first time-step on the other h a n d allows for the correction of  53  Top Camera  •a-  row + disparity  correlation window column + disparity columns  0 Oi  -D  column  - a  Left Cam  Reference Camera  F i g u r e 4.2: Triclops Stereo Image T r i p l e t  any i n i t i a l offset i n t r o d u c e d by the force probe m a k i n g contact w i t h the object. T h e short sequence of three image triplets means the tracker does not need to have new points to be tracked added i n each step. It just tracks a l l the i n i t i a l points t h r o u g h the complete sequence.  However, note that after the displacement of vertices has  been estimated the deformed mesh could be utilized to a d d a d d i t i o n a l points into the tracker.  4.1.2  Implementation Details  In this section, the i m p l e m e n t a t i o n details of the range-flow calculation are described starting w i t h the i n p u t images. T h e images from the stereo cameras are rectified such that the epipolar lines are aligned w i t h the image axes (see F i g u r e 4.2). T h e result of the stereo processing w i t h the Point  Grey Research  Triclops stereo vision  S D K is a triplet of row, c o l u m n a n d sub-pixel d i s p a r i t y at each valid pixel i n the reference image. A pixel is valid if the stereo processing produces a result. V a l i d pixels w h i c h are deemed to be on the object's surface are tracked i n the reference image. T h e area-based m a t c h i n g correlates a square window centered  54  at a given p i x e l w i t h a square search area i n the next image of the sequence from the same camera. T h e mask size of the correlation w i n d o w a n d the search area b o t h influence the result of correlation.  T h e influence of these parameters and of the  image size on the results a n d c o m p u t a t i o n time are demonstrated i n Section 4.1.3. T h e o p t i c a l flow correlation scores of the non-reference cameras are combined w i t h the ones of the reference camera w i t h the aid of stereo disparities. T h e pixels of the start a n d end of the flow vector i n the reference image are m a p p e d w i t h sub-pixel disparities into the stereo images. T h e correlation between these locations is calculated w i t h E q u a t i o n 4.1. E q u a t i o n 4.1 is the n o r m a l i z e d cross correlation NCC  of two image windows at c o l u m n s u b - p i x e l location using bilinear  image intensity interpolation. r  a  and c o l u m n c  is specified w i t h p i x e l a at l o c a t i o n row  plus 0 < p < 1 and p i x e l b at l o c a t i o n r  a  T h e NCC  T h e NCC  a n d q, plus 0 < q < 1.  b  is calculated from the non-normalized cross-correlation CC a n d the self  correlation SC.  C o r r e l a t i o n of image intensities I r r is computed over a w i n d o w C O  w i t h size 2 * n + 1 where n a positive integer.  (4.1) CC(r ,c ,r \c ,p,q) a  a  b  (1 -p)(l-  b  q)I rr{r CO  + (1  -p)ql rr(r  c  c0  +p(l  fc+l)  - q)I orr{r  Cb)  C  +PqIcorr(r  SC(r ,c a  (  P)  ((! - P ) ( l ~  P)Icorr(r ,C ,r a  +2.0(1 -p) I (r P  C  c  a  55  a +  t  a+l)  c  corr  +ppI orr{r ,  a  i,r , C a  a +  l)  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 reference 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  F i g u r e 4.3: R i g i d T r a n s l a t i o n and R o t a t i o n of a Test O b j e c t T r a n s l a t i o n by 5 m m is shown i n 4.3(a) and 4.3(b), r o t a t i o n by 2.0 ° deg. a r o u n d the object's z-axis (up i n the image) is shown i n 4.3(c) and 4.3(d). B r i g h t e r range-flow corresponds to larger m o t i o n .  an i n d i c a t i o n of the expected error. O n e can expect results of slightly lesser quality d u r i n g deformation or rotation. T h i s degradation is expected due to a significant change i n surface n o r m a l of the observed surface relative to the cameras a n d light sources. M o t i o n s w h i c h w i l l cause optical flow to vary over a correlation w i n d o w w i l l make the m a t c h i n g task more difficult a n d introduce error i n l o c a t i n g the m a t c h . Note, however, that the observed translation is not fronto-parallel to the image plane, i.e., the correct flow is changing over a correlation w i n d o w . T h e following tables compare the range-flow result to the correct t r a n s l a t i o n a l  57  m o t i o n . T h e correct m o t i o n is obtained from, reading encoder values of the highprecision linear positioning devices i n A C M E . T h e table shows the result for t r a c k i n g the object over eight time steps plus one time step for i n i t i a l i z a t i o n . Surface points i n the i n i t i a l i z a t i o n time step are tracked, no points are added to be tracked at later time steps. Therefore, the number of points tracked decreases w i t h each time step. T h e m o t i o n between steps is 2 . 5 m m . T h e tables report the true a n d estimated values for the orientation and length of the 3 D range-flow w i t h respect to the i n i t i a l i z a t i o n 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 t y p i c a l distance i n the scenario of this thesis to observe deformation from).  Furthermore, a l l results are for the test-object i n  F i g u r e 4.3. T h e first two tables, Table 4.1 a n d Table 4.2 show results for two different translations. T h e correlation w i n d o w size is 7 x 7 pixels w i t h 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 v a r i a t i o n i n the angle decreases w i t h greater translation while the variation i n vector length increases w i t h greater translation. T h e decrease i n orientational variation indicates the influence of the discretization error i n the o p t i c a l flow since larger flow vectors can be more precisely oriented i n the discrete image plane. T h e increase i n v a r i a t i o n i n the length of the flow vector is due to a reduction i n the density of the flow-field combined w i t h larger absolute flow values. T h e m e t h o d produces results for translation to about ^ m i l l i m e t e r accuracy w h i c h is sufficient for object deformation measurement.  The  m o t i o n to the left-rear i n the reference image (Table 4.1) produces similar results to the m o t i o n to the rear-right (Table 4.2). I n the following, the t r a n s l a t i o n to the left-rear of Table 4.1 is evaluated v a r y i n g the mask size, image size and search area  58  Step  #  A n g l e [deg] Mean  L e n g t h [cm] a  Mean  a  Encoder  Command  0.0340 0.0402  0.2355  0.2500  0.4864 0.7374  0.5000  0.9911  1.0000 1.2500 1.5000 1.7500  1  67.0081  10.2836  0.2678  2 3  64.7748 64.8529  5.7253 4.1890  0.5396 0.8113  4 5 6 7  64.8575 64.9783 65.0915 65.1626 65.2725  3.5898 2.8942 2.6112 2.3755 1.9633  1.0849 1.3567 1.6290 1.9019 2.1747  8  0.0464 0.0571 0.0626 0.0697 0.0770 0.0827  0.7500  1.2378 1.4867 1.7366 1.9929  2.0000  Table 4.1: R a n g e - F l o w Results for T r a n s l a t i o n M o t i o n to the left-rear Step  #  A n g l e [deg] Mean a  •1 2 3 4 5 6 7  48.8813 46.2486 46.0249 45.4467 45.6313  7.4717 4.5264 3.3112 2.6273 2.6246  45.4051 45.4808  8  45.2843  2.5060 2.4960 2.3798  Mean 0.2367 0.4926 0.7404 0.9997 1.2471 1.5086 1.7614 2.0266  L e n g t h [cm] a Encoder 0.0334 0.0461 0.0518 0.0507 0.0606 0.0780 0.0889 0.0937  Command  0.2230 0.4750 0.7146 0.9729 1.2161 1.4705 1.7192 1.9736  0.2.500 0.5000 0.7500 1.0000 1.2500 1.5000 1.7500 2.0000  Table 4.2: R a n g e - F l o w Results for T r a n s l a t i o n 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. T h e 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. T h e mask size varies for all images over 3 x 3 , 5 x 5 and 7 x 7 . For easier comparison a l l results are combined into a table per time step and time step 1, 4 and 8 are selected and shown i n T a b l e 4.3, Table 4.4 and T a b l e 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 i n the correlation window,  59  Image Size  M a s k Size  L e n g t h [cm]  A n g l e [deg] Mean  a  Mean  a  20.5544 20.3311  0.2816 0.2615  0.1378 0.0611  0.2630 0.3711 0.2741 0.2761 0.3072 0.2469 0.2446  0.0863  160, 120  3, 3 5, 5 7, 7  65.3530 64.2933 64.1025  19.6148  320, 240  3, 3 5, 5 7, 7  640, 480  3, 3 5, 5 7, 7  69.3456 64.9249 64.8249 65.7337  21.9100 14.4413 13.1080 15.5447  63.1974 63.4371  9.7911 9.8475  0.2286 0.0959 0.0789 0.1716 0.0460 0.0327  Table 4.3: Influence of M a s k a n d Image Size o n F l o w R e s u l t T h e search area is kept constant relative to the image size. T h e c o m m a n d e d m o t i o n is 0.25cm while the encoder value reads 0.2355cm (Step # 1 of Table 4.1).  but averaging w i t h i n a larger mask m a y introduce a larger error since the apparent m o t i o n i n the image over the mask m a y change. A d d i t i o n a l l y , the larger the search area the more likely is the inclusion of the correct flow l o c a t i o n but the greater the chance of a m i s m a t c h . 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 stability of the correlation peak w i t h respect to the size of the search area.  T h e more distinct the intensity i n the correlation  window is over the search area, the more unique the peak i n the correlation surface w i l l be. T h e time c o m p l e x i t y of the flow c o m p u t a t i o n is  0(size kSize h,sizei mas  searc  assuming the number of tracked pixels is p r o p o r t i o n a l to the size of the image. F o r example, the difference between the largest image w i t h the largest mask size a n d the smallest image w i t h the smallest mask size is 7 x 7 x 25 x 25 x 640 x 480  ~ 300 n  n  n  3 x 3 x 7 x 7 x 120 x 160 Therefore, it is crucial to be conservative i n chosen mask size, search area size a n d 60  )  rnage  Image Size  M a s k Size  A n g l e [deg] Mean  L e n g t h [cm] a  Mean  a  160, 120  3, 3  62.8964  11.9992  0.9983  0.1387  320, 240  5, 5 7, 7 3, 3  62.6035 62.5492 63.4279  0.9989 1.0017 1.0120  0.1115 0.1182 0.1350  5, 7, 3, 5, 7,  63.3405 63.0772 62.1532  10.9469 10.8381 10.9815 5.2712 4.7647 5.3999  0.0874 0.0946 0.0714  62.2516 62.1506  3.1899 2.8523  1.0713 1.0966 0.9891 1.0076 1.0060  640, 480  5 7 3 5 7  0.0514 0.0488  Table 4.4: Influence of M a s k a n d Image Size on F l o w R e s u l t T h e search area is kept constant relative to the image size. T h e c o m m a n d e d m o t i o n is 1.0000cm while the encoder value reads 0.9911cm (Step # 4 of Table 4.1).  Image Size  M a s k Size  160, 120  3, 3  63.9707  5, 5 7, 7 '  62.7593 62.7484  3, 3  61.3715 63.3863  320, 240  5, 5 7, 7 640, 480  3, 3 5, 5 7, 7  A n g l e [deg] Mean a  63.1105 65.4546 62.4504 62.1912  7.5256 7.8290 7.8810 3.7981 3.1567 3.2557 2.2675 1.9104 1.9717  L e n g t h [cm] Mean a 1.9500  0.1391  1.9893 1.9908 1.9008  0.1479 0.1429 0.1522  2.1400  0.1286  2.1945 1.9501 2.0245 2.0256  0.1330 0.1205 0.0673 0.0793  Table 4.5: Influence of M a s k and Image Size on F l o w R e s u l t T h e search area is kept constant relative to the image size. T h e c o m m a n d e d m o t i o n 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 displacements 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 selecting 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 neighborhood of a vertex has to be 2  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 The one-ring neighborhood of a vertex are all the vertices connected to this vertex by an edge. 2  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  Displacement  Results  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  F i g u r e 4.4: Deformation Measurement T h e 3 D range-flow shows the d i s t r i b u t i o n of the estimated flow vectors w i t h brighter range-flow corresponding to larger motions. T h e vertex displacement i n 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 Angle [deg] a Error 10.9142 23.2285 8.6968 20.3784 7.9037 18.9947 8.0567 21.4754 8.3435 23.5471 8.0368 24.2161 8.1894 26.5517 7.9706 28.0945  Translation in Y Angle [deg] a Error 7.4402 15.0430 5.8115 11.7741 5.2960 10.6937 5.0309 10.2249 4.8527 9.9227 4.8800 9.9940 4.9821 10.1555 4.9916 10.1864  Rotation around Z Length [ cm] Angle [deg] a a Command Error Mean 0.3053 23.5574 47.8135 0.2998 0.0518 24.8044 58.0988 0.5928 0.0884 0.6089 19.1940 42.1753 0.8739 0.1468 0.9052 17.5430 41.4193 1.1370 0.2130 1.2070 13.7172 37.8691 1.4744 0.2350 1.5336 12.2357 36.0586 1.7621 0.3434 1.8643 10.0026 32.9437 1.9929 0.4729 2.1297 8.2385 31.2005 2.1469 0.6727 2.3233  Table 4.6: V e r t e x Displacement Results for a T r a n s l a t i o n and R o t a t i o n T h e results are for the test-object i n F i g u r e 4.3. V e r t e x displacements are the median-filtered weighted average 3-D flow w i t h i n the triangles j o i n i n g at a vertex (see text for details). T h e displacements are transformed into object coordinates. T h e errors are the combined error over a l l vertices a n d are a result of b o t h registration error and flow estimation error.  64  Step  Method  1  NF MF MFWA MFWAS NF MF MFWA MFWAS NF MF MFWA MFWAS  4  8  Orientation [deg] Magnitude [cm] Error a Mean a 14.6050 18.8691 0.2388 0.0384 3.5504 4.4913 0.2373 0.0236 3.0627 3.9388 0.2372 0.0191 3.1004 3.9460 0.2373 0.0191 4.6907 6.6571 •1.0030 0.0572 1.8624 2.3026 . 1.0141 0.0245 1.7360 2.1340 1.0107 0.0192 1.7168 2.1106 1.0107 0.0192 2.6801 0.0737 3.0260 2.0122 1.9392 2.6460 2.0325 0.0528 1.8062 2.4697 2.0299 0.0484 1.8042 2.4669 2.0301 0.0488  Table 4.7: C o m p a r i s o n of E s t i m a t e s of V e r t e x Displacements for a T r a n s l a t i o n T h e results are for the test-object i n F i g u r e 4.3. M e t h o d s are nearest flow vector to the vertex (method N F ) , the m e d i a n of flow vectors w i t h i n the triangles j o i n i n g at a vertex (method M F ) , the weighed 1/(1 + d) average of flow vectors after m e d i a n filtering (method M F W A ) a n d m e d i a n filtering followed by the weighed squared 1/(1 + d) average of flow vectors (method M F W A D ) . See the text for an explanation of the methods. 2  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 i n the y-direction (to the rear-right i n the image) extensively used before.  The  vertex displacements are estimated from range-flow calculated i n a 640 x 480 image w i t h 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 r o t a t i o n  around the z-axis described above. T h e vertex displacements are again estimated from range-flow calculated i n a 640 x 480 image w i t h a mask of 5 x 5 a n d a search area of 25 x 25. T h e rotation penalizes averaging more since the v a r i a t i o n i n flow vectors is much higher over the set of triangles j o i n i n g at a vertex.  65  Step  Method  1  NF MF MFWA MFWAS NF MF MFWA MFWAS NF MF MFWA MFWAS  4  8  Orientation [deg] Error a 24.2158 28.0852 20.8808 22.5736 20.7130 22.1168 20.5988 22.0136 15.4796 16.1011 13.3114 18.3754 13.8109 18.7523 13.7860 18.6873 7.4156 9.8216 5.5896 9.6690 5.4954 9.5279 5.4283 9.3702  Magnitude [cm] Mean a Command 0.2923 0.0695 0.2826 0.2939 0.0299 0.2732 0.2939 0.0281 0.2736 0.2939 0.0282 0.2738 1.1682 0.1535 1.0821 1.1531 0.0911 1.0928 1.1531 0.0861 1.0932 1.1531 0.0839 1.0953 2.3229 0.2536 2.2000 2.2861 0.1968 2.1748 2.2861 0.1771 2.1867 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 N F ) , 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 additional imagery. The approach can be characterized as a purely local bottom-up approach. The measurement of the vertex displacements is independent of the deformation model. This independence is convenient because it allows a comparison of deformation models independent of displacement measurement. Even more importantly, 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 deforming object nearly always requires multiple views, each view covering different sets of 66  vertices. Different approaches to deformation measurement are possible. I n fact, m a n y researchers have suggested the use of deformation models to regularize various computer vision tasks (e.g., [127, 98, 21, 22]). A n approach e m p l o y i n g the deformation m o d e l as a constraint i n the vertex estimation is possible (e.g., [96]).  67  C h a p t e r  5  Estimation of Discrete Green's Functions E s t i m a t i o n of discrete Green's functions based on E q u a t i o n 2.1 can be viewed as a linear estimation p r o b l e m . T h i s view disregards the u n d e r l y i n g structure of the Green's functions m a t r i x and is explored i n this chapter.  T h i s is i n contrast  to  C h a p t e r 6 where the b o u n d a r y element m e t h o d derivation of the Green's functions m a t r i x of Section 2.3.2 is exploited to solve an inverse p r o b l e m . T h i s chapter can be understood as either the final step i n the m o d e l i n g process or as precursor to an inverse p r o b l e m solution (see C h a p t e r 6). F i r s t , the local behavior of an elastostatic object around the contact point is discussed. T h e local behavior is dominated by the compliant m o t i o n of the object because of the probe. T y p i c a l l y , the local deformation is significant but the deformation decays quickly w i t h the distance from the contact point. ( G l o b a l deformation does not decay quickly w i t h large scale m o t i o n of the object, e.g., due to articulation of the object as w i t h the example i n 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 functions matrix of a plush toy in Section 5.5.1 and an anatomic soft-tissue human wrist model in Section 5.5.2. 1  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 setup 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 The anatomic soft-tissue human wrist model has a bone structure made of hard plastic 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). x  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.  0  ui  Pfc  Ufc  0  not observed  never excited  unobservable  (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 A C M E . 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 u , while the force-torque sensor measures the contact force. The force at k  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  T h e surface is represented as a triangular mesh and the shape functions for the tractions are linear. T h e components of the force vector f  k  are the s u m over all  integrated t r a c t i o n shape-functions on elements connected to the p r o b e d vertex, i.e., fk = Ylj \ {Pk + Pk ~^Pk) of tractions p  lj  .. .p  lj  1^1-  T h e t r a c t i o n over triangle j is a-linear c o m b i n a t i o n  at its vertices 1 . . . 3; while \G\ is the J a c o b i a n from local to  global coordinates (see A p p e n d i x A for more details). T h e estimation p r o b l e m is: given a set of block vectors v and a set of corres p o n d i n g block vectors v find the m x m observable s u b m a t r i x of 3 . T h i s p r o b l e m can be broken up into m separate estimation problems for each of the m columns due the specific structure of E q u a t i o n 5.1. E a c h of these m estimation problems per c o l u m n consists of finding a 3 x 3 block element  of m a t r i x S given a set of  observed displacement vectors u , and the corresponding t r a c t i o n vectors p^- T h e p r o b l e m w h e n i — k is the estimation of the compliance of the object at l o c a t i o n k. T h i s is discussed next i n 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 ^ , the inverse of the self-compliance m a t r i x . T h i s local stiffness 1  is employed i n haptic simulation to map displacements of the tool t i p into a force response [58, 57]. T h e correct force response can be easily calculated at the high haptic rendering rates required (typically «  1000Hz).  A l l that is required is to  m u l t i p l y the local stiffness at the contact point(s) w i t h the displacement vector(s) due to the t o o l . For a single contact vertex, the equation is  71  For this approach, it is assumed t h a t the self-compliance elements of the Green's functions m a t r i x are invertible. T h e compliance m a t r i x w i l l be invertible as long as the response is not rigid or infinitely stiff i n any direction.  T h i s local m o d e l  is exact w i t h respect to the global m o d e l and is a u t o m a t i c a l l y synchronized w i t h 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]). T h e self-compliance estimation i n the formulation of Section 5.2.1 is a linear least-squares p r o b l e m 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 i n three dimensions. T h e force response of the object is recorded giving I traction vectors. T h e standard least squares p r o b l e m to be solved for  HjT is cf  then  _ 12 "fcfc — fc fc r-^T  PfcPfc • • • Pfc  U  U  In the standard notation for least squares, this becomes Ax equations  = b.  The  normal  A A x = A b . T h e solution of this system T  of the p r o b l e m are then  T  can be found by singular value decomposition of the measurement design m a t r i x A = USV . Here, the 3 x I m a t r i x U and the 3 x 3 m a t r i x V have o r t h o n o r m a l T  c o l u m n vectors [ui, • • •, uj] and [vi, • • •, V3], respectively. T h e 3 x 3 diagonal m a t r i x X contains the ordered singular values o\ > 02 > 0-3. T h e n o r m a l equation may then be expressed as  V E V a ; = VSU 6. 2  T  T  72  R e o r d e r i n g this equation leads to  where the diagonal m a t r i x X I / 1 7 3 .  -  1  contains the ordered singular values l/o~i < I / 0 2 <  Least squares problems w i t h a measurement design m a t r i x A w i t h a poor  c o n d i t i o n number cannot be solved w i t h the standard procedure above. C o n d i t i o n numbers of the design m a t r i x A w h i c h occur based on measured t r a c t i o n vectors are considered next.  5.2.2  Measurement Design Matrix  I n theory, i n a numerical solution of a discrete elastostatic b o u n d a r y value problem, rank deficient local stiffness matrices (HT^ ) are rarely encountered if at a l l . 1  For isotropic materials a n d most geometries local stiffness matrices are diagonally dominant, i.e., if a force is exerted on an object, the objects deforms locally i n the direction of the applied force. However, measuring the response of an object is different since the design m a t r i x A = [pjfcPfc • ' ' PL] t r a c t i o n measurements p  l k  n  a  s  to be well conditioned. T h e  are p o l l u t e d w i t h measurement noise a n d m o d e l i n g error.  A s mentioned previously, it is ensured that the probe directions span 3-space in setting up the least squares p r o b l e m . T h i s guarantees that the local compliance m o d e l is excited i n all directions. Nevertheless, the measurement of the traction response of the object m a y not span 3-space. O n e probe d i r e c t i o n is n o r m a l to the mesh surface of the object, while a d d i t i o n a l directions are offset from the n o r m a l . T h i s offset is l i m i t e d by the fact that the contact between the probe a n d the model is not allowed to slip for a measurement of the static force response at a given contact point. I n A C M E a probe t i p w i t h high friction (grinding stone tool tip) is used i n 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  F i g u r e 5.1: P r o b e T i p at V e r t e x L o c a t i o n Vertex 38  Ol 1.6426  0-2 0.4317  0-3 0.3940  55  2.5786  0.2638  0.1762  Vertex  VVan 2.3513  \]Var<2 2.0164  \fVar$  38 55  4.1377  3.8120  0-1/0-2 3.8047  cri/o-  3  4.1696 14.6327  1.6076  9.7753 \\pT<aT TTT\\ IK ^kk112 0.4094  3.8823  0.6973  6.8947  U  3.9993  Table 5.1: L o c a l Force Response L o c a t i o n of vertices is shown i n F i g u r e 5.1. Response at vertex 38 is well conditioned while the response at 55 is p o o r l y c o n d i t i o n e d .  certainty i n measurement of the force response i n the d i r e c t i o n of the surface n o r m a l t h a n i n any off-normal direction. For the remainder of this section measurements of the t r a c t i o n response of a plush toy tiger (see F i g u r e 3.8) at two locations are considered to illustrate the discussion. O n e location is vertex 38 of the base mesh on the back of the tiger as shown i n F i g u r e 5.1(a), while the other location is vertex 55 of the base mesh on the head of the tiger as shown i n F i g u r e 5.1(b). T h e t r a c t i o n response on the location on the back spans 3-space well. T h i s is apparent from the ratio of singular values listed i n Table 5.1. T h e t r a c t i o n response at the measurement location on the head  71  of the tiger is very different. T a b l e 5.1 shows that the ratio of singular values for vertex 55 is larger t h a n for vertex 38. T h e c o n d i t i o n number at vertex 55 is 14.6327 compared to 4.1696 at vertex 38. A- condition number of 14.6327 does not cause any n u m e r i c a l difficulties but here it signals a measurement problem. T h e larger condition number is caused by a t r a c t i o n response w i t h a dominant direction. T h e dominant direction is due to the head, of the tiger apparently hinging around its neck. T h e dominant m o t i o n is around the axes of this apparent hinge. T h e directional change of the p r o b i n g does not change the direction of the force response significantly enough. Here, not significantly enough means that the measured force response of the n o r m a l 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 n o r m a l direction of the probe is changing over the surface. T h e singular values <7i of the measurement  design m a t r i x determine how  sensitive the least squares solution is to errors i n the associated right singular vectors V i . T h i s can be seen by r e w r i t i n g E q u a t i o n 5.2.1 as  .  i=i  '  A s s u m i n g the error is equally d i s t r i b u t e d i n the components of the right-hand-side measurement vector 6, the factor l/<7i acts as an amplification of the error. A design m a t r i x w i t h singular values differing by more t h a n one order of magnitude, as i n the case of vertex 55 on the head of the p l u s h toy tiger, causes the signal to noise ratio to decrease by more t h a n one order of magnitude. I n Section 5.2.3 regularization is used to deal w i t h p o o r l y conditioned t r a c t i o n responses.  75  5.2.3  Regularization  R e g u l a r i z a t i o n transforms a parameter estimation p r o b l e m into a constrained o p t i m i z a t i o n p r o b l e m . T h e goal of regularization is to find an estimate w h i c h m i n i m i z e s the residual of a parameter fit (just as i n least squares) but constrained by a distance measure from an a priori k n o w n solution. T h i s section discusses the trade-off between choosing a solution w i t h low residual a n d choosing a solution w h i c h fits the a priori constraint i n the context of estimating the local compliance m a t r i x S^ . fc  T h e direct regularization methods of truncated singular value decomposition  ( T S V D ) a n d d a m p e d singular value decomposition ( D S V D ) are compared w i t h the non-regularized solution. T h e T S V D (truncated at t) solves the p r o b l e m  min||x||2 subject to m i n | | A x — t»|j2t  T h i s solution is calculated based on the S V D ( E q u a t i o n 5.2.2) as  \-r  u  b  where fj is a filter factor a n d f; = 1 V i < t a n d 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 c o m p a r i n g the residual of fit a n d constraint error, another m e t h o d for choosing the regularization threshold A is the generalized cross validation ( G C V ) . Below, the example of vertex location 55 shown i n F i g u r e 5.1(a) is evaluated. Some of the results i n this section are obtained w i t h the M a t l a b package  Regularization  Tools published b y P . C . Hansen [48]. Consider the variance i n the elements of the estimated compliance m a t r i x . T h e variance Varj  i n the rows j of the transposed compliance m a t r i x EjT are [102, fc  76  pp. 676-681] Var  3  = £ i=i  & ^  T h e variance w i t h i n a row of the m a t r i x H^, is identical because of the fact that estimating &  T k  is three similar least squares problems w i t h three parameters each.  T h e reciprocal of a singular value 1/CTJ scales the c o n t r i b u t i o n of the corresponding basis vector to the variance of a parameter as well as to the estimate of the parameter.  These variances describe the uncertainty i n the compliance m a t r i x i n  centimetres[cm]  per Newton[N}/centimeter  [cm ]. T h e variances of the three rows  2  2  are listed i n T a b l e 5.2 for the illustrative 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 i n E q u a t i o n 5.2.3. If the inverse of a singular values is large, the uncertainty i n the associated c o l u m n of right-hand-side singular vector is large.  T h i s leads to large variances V'arj.  T h e determination  of large corresponds to the choice of the regularization parameter. estimating the 3 x 3  I n the case of  local compliance m a t r i x , the choice can have three possible  outcomes, i = l , i = 2 o r i = 3 (no regularization). T h e three possible outcomes are listed i n T a b l e 5.2. A s expected, the n o r m of the solution decreases w i t h the increase i n regularization. T h e regularization parameter can be chosen from Table 5.2 given an expected error (residual of fit). In the estimation p r o b l e m at hand, one may be able to deduce an expected error by c o m p a r i n g a l l the i n d i v i d u a l estimation problems. However, the error is dependent on location and on the particulars of a measurement and sett i n g a general threshold may be difficult. O n e approach to deal w i t h this difficulty is to estimate the threshold from the d a t a itself. T h e generalized cross v a l i d a t i o n does exactly that.  T h e idea is to find a regularized solution w h i c h is capable of  77  Vertex 55 \/Var \JVar3 2  ||P Sjk-lF|| r  2  Il fc/cll2 3  o  x  = 2.5786 a  2  t = 3  t = 2  t = l  4.1377 3.8120 3.8823 0.6973  4.1377 3.8120 0.5722 1.3721  4.1377 0.6516 0.2809 1.6926  6.8947  3.8992  2.5143  = 0.2638 a  3  = 0.1762  Table 5.2: R e g u l a r i z a t i o n w i t h the T S V D T h e table shows a comparison of the variances Vari of the rows of the local r compliance m a t r i x SjQ., residual of fit | | P S ^ f c — U | | 2 and the n o r m of the solution ||Sj ||2. T h e number of non-zero singular values is t. r  fc  0.015  F i g u r e 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 w i t h the truncated singular value decomposition ( T S V D ) . T h e G C V for each row of H^ w i t h k = 55 is shown separately, as well as their s u m (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. fc  78  p r e d i c t i n g an a r b i t r a r y measurement if this measurement is d r o p p e d from the estim a t i o n . A d d i t i o n a l l y , the choice of regularization parameter should be independent of a orthogonal transform of the measurement data. T h i s leads to a G C V [48] to be m i n i m i z e d . T h e G C V is calculated by G (traced  "  P (E?=i T  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 i n the example is at t = 2 basis vectors (see F i g u r e 5.2). A t first, the T S V D m a y seem like a good fit to the estimation problem at h a n d . However, the problem w i t h the solution obtained is that the estimated compliance m a t r i x 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 m a t r i x (i.e., the local stiffness m a t r i x ) . T h e local compliance m a t r i x w i l l be completely rigid, 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. T h e d a m p e d singular value decomposition ( D S V D ) replaces these filter factors by the smoother function of E q u a t i o n 5.2.3. T h e D S V D solution for the example is listed i n Table 5.3 for regularization parameters A resulting i n the same residual as for the T S V D of Table 5.2. T h e G C V curve i n case of the D S V D is increasing m o n o t o n i c a l l y (see F i g ure 5.3).  In this thesis the regularization parameter is chosen heuristically.  The  heuristic is based on the e m p i r i c a l observation that the transition between p o o r l y and well-conditioned t r a c t i o n 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. T h i s leads to a regularization parameter of A = 0.08164 for the p o o r l y conditioned 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  2.8340  3.8120 3.8823  1.4809 1.7081  1.0735  \JVar2  1.2802 1.0181  2.9026 2.6656  1.4000  0.6973  1.3721  1.6926  0.8489  6.8948  . 2.9310  2.4457  4.8591  0-1 = 2.5786 cr = 0.2638 0-3= 0.1762 2  Table 5.3: R e g u l a r i z a t i o n w i t h the D S V D T h e table shows a comparison of the variances Vari of the rows of the local stiffness m a t r i x S^ , residual of fit | | P S ^ — U ||2 and the n o r m of the solution ||S^ ||2. T  T  fc  f e  0  fc  0.2  X  0.4  0.6  F i g u r e 5.3: G C V for D S V D M i n i m u m (on the left) of the generalized cross v a l i d a t i o n ( G C V ) for regularization w i t h the d a m p e d singular value decomposition ( D S V D ) . T h e G C V for each row of E^f, w i t h k — 55 is shown separately, as well as their s u m (solid line).  vertex 38. T h e heuristic achieves the goal of regularizing p o o r l y conditioned traction responses but fitting the d a t a o p t i m a l l y for well conditioned t r a c t i o n matrices. A l l the results for the local object behavior i n Section 5.5 are calculated w i t h this approach. In general, the D S V D solution has the advantage of having full rank as long as the design m a t r i x A has full rank. T h i s avoids difficulties i n haptic simulation of the stiffness matrices. Furthermore, all the left h a n d singular vectors u; are damped, especially all off-normal vectors w i t h their corresponding s m a l l singular values. T h i s  80  stabilizes the solution i n the face of noise i n the off-normal directions.  5.3  Global Behavior  T h e measurement of the global behavior of the m o d e l is quite different from the measurement of local behavior of Section 5.2 above. T h i s is not obvious at first since the structure of the p r o b l e m is identical. T h e design m a t r i x A of the i n d i v i d u a l least squares problems for estimating aj  k  when j ^ k is identical to the design m a t r i x of  the self compliance m a t r i x H^ . However, the right h a n d side b of the problem, the fc  measurement of the vertex displacement is very different. T h e measured surface displacement i n 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 C h a p t e r 4. T h i s results i n the following m a i n differences w h i c h need to be addressed:  increased noise level, a d d i t i o n a l noise i n the form of  outliers, and incomplete measurements. Outliers can, on occasion, still be observed i n the displacement  measurement  despite the filtering process described i n Section 4.2 (see Section 4.2 for a description of shortcomings of the filtering process). T h i s should not be a surprise since most visual measurement  techniques need to deal w i t h outlier noise.  M e t h o d s to deal  w i t h these outliers are discussed i n 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. T h e p r o b l e m of incomplete measurements manifests itself b o t h i n a reduction of the number of observations per estimate at a vertex and i n a complete lack of observations for some vertices. Fewer measurements at a single vertex can be dealt w i t h by increasing regularization w i t h the same  81  methods as i n 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 p r o b i n g at a different but close-by location. C o n s i d e r i n g neighboring vertices corresponds to interpolation w i t h i n a c o l u m n of the block m a t r i x S, while similar probes correspond to interpolation w i t h i n a row of the block m a t r i x H. A l s o note that o b t a i n i n g 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 all observations. G l o b a l regularization schemes are discussed i n C h a p t e r 6 but see also Section 5.4. T h i s section only considers interpolation schemes w i t h i n a block c o l u m n of the Green's functions m a t r i x S. Here, scattered d a t a reconstruction is used to fill i n elements for each c o l u m n 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 i n Section 5.3.3.  5.3.1  Regularization  T h e off-diagonal m a t r i x elements Sj  k  w i t h j ^ k are estimated u t i l i z i n g the d a m p e d  singular value decomposition ( D S V D ) . T h e m e t h o d is identical to the estimation of the diagonal elements E^. described i n Section 5.2.3, however the regularization parameter A is chosen differently. A d d i t i o n a l l y , the estimation is not attempted if the set of displacement measurements is not the result of probe directions w h i c h  82  span 3-space. T h i s occurs as a result of a reduction i n the. number of observations for a given vertex.  In this case all deformation measurements at t h a t vertex are  disregarded and the hole filling process is left to fill t h e m i n . T h e estimation process is also integrated w i t h the outlier rejection process of Section 5.3.2.  T h e choice  of regularization parameter i n the D S V D for vertices w i t h sufficient observation is discussed below. T h e noise i n the deformation measurement effects the direction of the displacement much more t h a n the magnitude of the displacement Section 4.2).  (as discussed i n  R e g u l a r i z a t i o n w i t h the D S V D is again well suited for this p r o b l e m .  S m a l l variations i n direction should not be magnified by s m a l l singular values of the.design m a t r i x A.  T h e D S V D dampens these s m a l l singular values. T h e o n l y  difference then, compared to the local behavior, is i n the choice of regularization parameter. A large regularization parameter filters the singular values more. Since the noise i n the visual displacement measurement is larger t h a n the noise i n the measurement of the displacement w i t h the C M S , a larger regularization parameter employed for the estimation of the global behavior than for the estimation of the local behavior. T h e 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 still not be sufficient at some locations w i t h large errors. These locations need to be identified and the regularization increased correspondingly. A n o t h e r heuristic is employed i n these cases. A way of l o o k i n g 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.  T h e local behavior defines a ratio between observed m a x i m u m  compliance max{||I2/I(pi112)  over all measurements I at vertex k a n d the n o r m  83  of the solution ||H^||2 (the largest singular value o f the s o l u t i o n m a t r i x ) . In the case of the illustrative example of vertex 55 on the head of the tiger, the m a x i m u m observed compliance is max( | j u.^. j 12 /11  R  A  T  L  0  K  =  112) = 3 . 1 1 7 2 . T h i s gives a ratio of  ! a x p ' i. | |i2n) m a x 711 ( | | u i |!| ? 2^ / | !| m 1  T  112  FC  =  L  5  5  8  8  '  A  T h i s ratio is also calculated for a l l other vertices j ^ k i n the c o l u m n , i.e., ratio j = ll Jfcll2/( H  m a : r  (ll jll2/||pil|2))u  T h e solution S j  f c  is found w i t h the D S V D w i t h A =  <7i/5 — 03. T h e estimated s o l u t i o n at a vertex j is suspicious if the ratio is higher t h a n for the local compliance. I n such a case, a larger regularization parameter may produce a more desirable s o l u t i o n . Therefore, the D S V D is repeated w i t h the regularization parameter A = o i / 5 * ration/ratioj  — 0-3. T h i s heuristic achieves the  goal of rejecting solutions w i t h norms larger t h a n expected but allowing a r b i t r a r y 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.  T h e p r o b l e m is to select a subset of measurements w h i c h does not  contain any spoiled measurements.  A number of selection methods (or positive-  breakdown methods) exist w i t h i n robust estimation techniques, e.g., r a n d o m sample consensus ( R A N S A C ) [14], least m e d i a n of squares ( L M S ) [110] and least t r i m m e d squares ( L T S ) [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 m o d e l to estimate.  G i v e n a p r i o r i knowledge of  the m a x i m u m measurement error, inlier points (i.e., points w h i c h are not outliers) are selected for a given estimate.  If at least a m i n i m u m number of points are  w i t h i n the expected measurement error, the estimate is accepted. 84  F i n a l l y , a least  squares fit to the inlier points only is c o m p u t e d 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 t h a n ( M + P + l ) / 2 . L M S proceeds identically except it does not use an expected error to classify inlier points. Instead a l l possible P sets are evaluated and the estimate w i t h least median of squared residual is selected. A l l points w i t h a squared residual smaller t h a n 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 i n a l l y , L T S is an improvement over L M S . I n L M S the evaluation of a t r i a l set P depends on the errors i n the d a t a points P, i.e., the m e d i a n of squared residuals of a d a t a set m a y have been considerably lower if the i n i t i a l estimate w o u l d have been improved by a least squares step.  Therefore, L T S improves each subset of  points H by i m p r o v i n g the fit w i t h repeated least squares [111]. T h e least squares iteration has converged when the H set has become stable.  L T S finds the H set  w i t h the m i n i m u m 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 i m i l a r l y to the description i n [111]. One modification of the a l g o r i t h m is i n the selection of the i n i t i a l P-set. P = 3 for the estimation of S j  f c  since each c o l u m n is represented by the same parameter  (as  discussed i n Section 5.2.3). T h e r e is no point to selecting an i n i t i a l subset P w h i c h is rank-deficient. Therefore, bins are employed to group t r a c t i o n measurements together w h i c h are obtained w i t h the same orientation of the displacement vector. B i n n i n g has been suggested, e.g., i n [134].  T h e 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 d a t a w i t h i n the block c o l u m n k of the Green's functions m a t r i x Sj . k  T h i s is possible since the observable part of the Green's  functions m a t r i x contains compliance matrices o n l y (see E q u a t i o n 5.1).  A column  of the compliance m a t r i x can be viewed as the displacement due to a unit t r a c t i o n applied along the axis of the object coordinate system. A s a consequence, the three columns of the observable part of Sj  k  for j = 1,..., n a l l represent compatible values.  If the displacement field is smooth, the columns are s m o o t h .  T h i s smoothness  assumption is t y p i c a l l y 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 s m o o t h l y w i t h 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. S m o o t h hole filling can be achieved by solving L a p l a c e ' s equation over vertices j w i t h unestimated  Sj  subject to estimated values S;fc at vertices I and  k  S ,fc = 0 for vertices w w h i c h are fixed. T h i s may be c o m p u t e d efficiently i n several tt  ways, however time-stepping the diffusion equation sufficiently long is effective here. T h e discrete L a p l a c i a n L at time t and vertex j is calculated using d a t a from the vertex's one-ring n e i g h b o r h o o d , Nj 2  where e i is the edge j o i n i n g vertex j and i. E x p l i c i t forward E u l e r discretization is 3  computed as B (t jle  + At) = E (t) jk  + At  L(S {t)) ]k  The one-ring neighborhood of a vertex is all the vertices joined to that vertex by an edge. 2  86  subject to a time-step restriction on A t [32]. T h e process is started w i t h the u n estimated Sj (t k  = 0) — 0. Constraints can be i n t r o d u c e d w i t h a weighting factor  0 < 8 < 1 a n d the E u l e r time step becomes  S (t 3k  + At) = 8 3 (0) jk  + (l-8)(3 (t) jk  +  AtL(E (t))). jk  T h i s simple process is quite effective for filling a l i m i t e d number of holes, even though explicit methods m a y require a large number of s m a l l time steps. S o l v i n g Laplace's equation has also been applied successfully to s m o o t h displacement fields and i n geometric mesh fairing [32]. T h e result of this interpolation process is shown i n F i g u r e 5.7.  5.4  Mesh  Refinement  T h e s u b d i v i s i o n mesh hierarchy can be exploited to improve rendering quality a n d reduce measurement a n d estimation time. T h e object can be contacted at vertices of a coarse resolution mesh but the global behavior m a y be estimated at any mesh of level / given sufficient spatial resolution of the visual measurement of the deform a t i o n . A s a result the estimated Green's functions m a t r i x S at the resolution I of the deformation measurements lacks columns w h i c h correspond to o d d vertices inserted at that s u b d i v i s i o n level. I n order to fill i n these columns, i n t e r p o l a t i o n w i t h i n a row is necessary. In this thesis, the missing columns are interpolated for rendering purposes only. T h e i n t e r p o l a t i o n is linear except for the one-ring neighborhood of the local response. Here a special shifting is done to preserve the diagonally dominant structure of the overall Green's functions m a t r i x a n d prevent linear dependence columns of S (see E q u a t i o n 5.2). It is stressed that global inverse p r o b l e m solutions are b y 87  far preferable to such a crude interpolation. However, for rendering the estimated Green's function for the examples presented i n Section 5.5, the described approach produced good results. Note as well that this scheme may be utilized for adaptive s a m p l i n g of the surface. A s s u m e the mesh is probed at mesh resolution Z and the discrete Green's functions m a t r i x is estimated. T a k i n g the columns of the even vertices, i.e., corresponding to vertices of mesh level I — 1, the response at o d d vertices may be found by interpolation. If the measured response at o d d 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 s u b d i v i s i o n level Z. Level Z has been subd i v i d e d once from the level Z — 1 and all observable block columns of S are available at level Z — 1. T h e n , the j  th  o d d vertex on level Z has a response S y inferred if b o t h  even vertices k and m of its parent edge have responses 3^ and 3; , m  respectively  (see E q u a t i o n 5.2). If so, S j j is linearly interpolated from these two responses. T h e local behavior, 3jj  and S - when vertex n is a one-ring neighbor of j, is handled n j  88  differently.  T h i s local behavior is d o m i n a t e d 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 w h i c h are diagonally shifted (see E q u a t i o n 5.2, S j n  is handled exactly as Sjj)-  linearly interpolated from H^fc and S  m  m  For example, Sjj nor 3 j  and not from S j k  m  is  as for far field  responses.  5.5  Results  In this section, estimation results for two objects are presented, one a p l u s h toy tiger shown i n F i g u r e 3.8, the other a anatomic soft tissue h u m a n wrist m o d e l . Present a t i o n of the results on paper is somewhat l i m i t i n g . V i s u a l and haptic interaction w i t h the models provides a far clearer picture of the results . Quantities w h i c h can 3  be presented here include the residual of fit for the local behavior estimation, i.e., for the diagonal of the Green's functions m a t r i x . T h e d i s t r i b u t i o n of the m a x i m u m eigenvalue of the compliances (i.e., the diagonal of Green's functions m a t r i x ) 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. T h i s is meaningful since these columns present the response of the m o d e l to an applied unit t r a c t i o n i n a given direction at a vertex location. Furthermore, the effect of the regularization process, the outlier rejection by L T S and the hole filling by L a p l a cian s m o o t h i n g are presented for the p l u s h toy. F i n a l l y , the effect of different mesh resolutions is presented. H a p t i c interaction with the plush toy. tiger has been successfully demonstrated to audiences 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) 3  89  '0  200 Time [10ms]  400  .(a) Probe at Vertex 38 of Figure 5.1(a)  0. 0  200 Time [10ms]  400  (b) Probe at Vertex 55 of Figure 5.1(b)  F i g u r e 5.4: M e a s u r e d M a g n i t u d e of Force a n d Displacement  5.5.1  Plush Toy  F i r s t , results for the local behavior of the model are presented.  F i g u r e 5.4 shows  the magnitude of force a n d displacement at the probe t i p d u r i n g deformation. T h e tip executes a spatial linear m o t i o n w h i l e i n contact w i t h the object. A n ideal linear elastostatic object sensed w i t h a perfect sensor w o u l d cause the force profile to be a m u l t i p l e of the displacement.  T h e plots i n F i g u r e 5:4 show a p p r o x i m a t e l y such  a behavior, however there are also differences. T h e force does not r e m a i n constant after the probe stops its m o t i o n b u t decays exponentially to a lower value. T h i s decay corresponds to viscoelastic object behavior. T h e p l u s h toy's filling 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 u n t i l the force reaches a threshold. T h e displacement  90  is offset by a corresponding amount not shown i n the graphs. D u r i n g measurement the force profile is linearly interpolated to zero, as well as the displacement.  The  force sensor causes some high frequency noise i n the force profile. T h e result of the estimation of the discrete Green's functions is shown i n F i g u r e 5.5. F i g u r e 5.5(a) shows the m a x i m u m singular value of the diagonal of the Green's functions m a t r i x . T h e diagonal describes the compliance of the solid. T h e figure shows that the tiger is m a x i m a l l y compliant at the rear of the head. In general the head is more compliant t h a n the back. T h i s corresponds well to the inspection of the plush toy by h a n d . F i g u r e 5.5(b) shows the residual of fit | | P E j T T  - XJ || X  f c  2  for  the local compliance matrices. T h e residual is u n i f o r m l y 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 p r o b l e m locations are on the ears of the tiger and on the t a i l . T h e d i s t r i b u t i o n of the residual is clearly different from the d i s t r i b u t i o n 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 illustrated w i t h the two Green's functions block columns corresponding to vertex 38 and 55. T h e l o c a t i o n of vertex 38 and 55 are shown i n F i g u r e 5.1(a) and (b), respectively. F i g u r e 5.6 shows the effect of regularization on the y - c o l u m n of the Green's functions m a t r i x . T h e x- and z-column show similar behavior and are not shown. T h e behavior w i t h o u t regularization shows large Green's functions magnitudes w i t h an uneven d i s t r i b u t i o n . R e g u l a r i z a t i o n ensures that noise i n the observations do not cause solutions w i t h excessive norms.  T h i s can be observed i n the reduction of the Green's function  y-vectors. T h e estimates for vertex 38 are satisfactory w i t h o u t regularization, but improve w i t h regularization. N o t i c e that on the head of the tiger the movement  91  /A (a) Solution Norm | | E £ | | 2 fc  (b) Residual of F i t | | P 3 j [ T  - U || T  fe  2  F i g u r e 5.5: L o c a l C o m p l i a n c e T h e figure is vertex colored, red indicates a large value while blue indicates the smallest value. Vertices not p r o b e d due to reachability constraints result i n S? not being estimated a n d are shown i n black. T h e result shown is o b t a i n e d w i t h a m a x i m u m of eight probes per vertex location and observations from two different k  viewpoints.  92  (c) Probe at Vertex 38 of Figure 5.1(a)  (d) Probe at Vertex 55 of Figure 5.1(b)  F i g u r e 5.6: R e g u l a r i z a t i o n of a Green's F u n c t i o n s E s t i m a t i o n E s t i m a t e d y-columns of the discrete Green's function w i t h o u t regularization are shown i n 5.6(a) a n d 5.6(b). T h e result w i t h regularization (as i n Section 5.3.1) is shown i n 5.6(c) a n d 5.6(d) T h e y - c o l u m n of the observable p o r t i o n of a block c o l u m n of the discrete Green's function corresponds to the displacements of the free vertices due to a unit traction i n the y-direction (see E q u a t i o n 5.1). T h e displacements are color-coded where red indicates a large value w h i l e blue indicates the smallest value.  93  gets d a m p e d while it persists around the probe location. T h e estimates for vertex 55 are i n clear need of regularization as discussed i n Section 5.2. Observe how there are large displacements i n stray directions on the head of the tiger before regularization. A g a i n , regularization dampens these directions but preserves the observed displacements near the p r o b e d vertex. T h e behavior for vertex 38 is more local t h a n for vertex 55. F i g u r e 5.7 shows the effect of L a p l a c i a n s m o o t h i n g described i n Section 5.3.3. T h e black areas corresponding to unestimated displacements due to lack of inform a t i o n are s m o o t h l y filled i n . F i g u r e 5.8 shows some subtle improvements over F i g u r e 5.7 due to a d d i t i o n a l s m o o t h i n g of existing visual observations using the L a p l a c i a n and by u t i l i z i n g L T S d u r i n g the estimation. Unfortunately, this is h a r d to show on paper.  T h e major improvement is i n the smoothness of the direction  of the Green's functions over the surface. Nevertheless, the perimeter of the area of large displacement (shown i n 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 h a r d l y visible.  F i g u r e 5.9 shows the  same measurements w i t h estimation done on a finer mesh. T h e finer mesh level is better able to resolve the changes i n the Green's functions. A s a result the magnitude of the y - c o l u m n of the Green's function changes more smoothly w h i c h is visible as smoother color transitions i n F i g u r e 5.9.  5.5.2  Anatomic Soft-Tissue Human Wrist Model  A s before w i t h 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 d u r i n g a p r o b i n g are shown i n F i g u r e 5.11 for the locations indicated i n F i g u r e 5.10. T h e 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)  F i g u r e 5.7: Hole filling Discrete Green's functions w i t h o u t hole filling are shown in 5.7(a) a n d 5.7(b), while the result of hole filling by L a p l a c i a n s m o o t h i n g ( Section 5.3.3) is shown i n 5.7(c) a n d 5.7(d). T h e magnitude of displacements are shown w i t h vertex coloring, where red indicates a large value while blue indicates the smallest value. B l a c k 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)  F i g u r e 5.8: Hole filling Results w i t h L a p l a c i a n s m o o t h i n g (Section 5.3.3) u t i l i z e d to s m o o t h estimated values is shown i n 5.8(a) and 5.8(b). T h e results are obtained w i t h 25 iterations w i t h At = 0.2, 8 = 1 for the probed vertex and a l l fixed vertices, while 8 = 0.75 for vertices estimated based o n visual information). T h e results for Least T r i m m e d Squares ( L T S ) as described i n Section 5.3.2 are shown i n 5.8(c) a n d 5.8(d). T h e magnitude of displacements are shown w i t h vertex coloring, where red indicates a large value w h i l e blue indicates the smallest value. B l a c k vertices are not estimated.  96  F i g u r e 5.9: Discrete Green's Functions for Different M e s h R e s o l u t i o n T h e result for the base mesh of the s u b d i v i s i o n hierarchy are shown i n 5.9(a) a n d 5.9(b). T h e results for the next finer mesh resolution (see S e c t i o n 5.4) are shown i n 5.9(c) a n d 5.9(d). T h e estimation for b o t h meshes is based o n the same set of probes a n d images. A g a i n , the magnitude of displacements are shown w i t h vertex coloring, where red indicates a large value while blue indicates the smallest value. B l a c k vertices are not estimated.  97  (a) Vertex 4  (b) Vertex 98  F i g u r e 5.10: P r o b e T i p at Vertex L o c a t i o n  have similar structure to those of the plush toy. D u r i n g the m o t i o n of the probe, the force is linear w i t h displacement.  However, as the m o t i o n stops the displacement  is maintained w i t h an exponentially decreasing force.  T h i s is a s m a l l viscoelastic  effect. T h e stiffness at the two locations shown i n 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 i n F i g u r e 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 m o u n t i n g surface. T h e v a r y i n g stiffness of the soft-tissue wrist is also clearly shown i n i n F i g ure 5.12(a) w h i c h is the plot of the d i s t r i b u t i o n of the m a x i m u m singular values of the compliance matrices (the diagonal of the Green's function m a t r i x ) . F i g u r e 5.12(b) shows the residual of fit  HP^Sj^ —  U | | 2 for the local compliance matrices. T  The  residual of fit is quite low and uniform over the mesh w i t h the exception of the thumb. T h e object's internal bone structure is simplified i n the fingers. T h e structure i n the fingers contains no joints, does not extent the full length of the  98  finger,  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)  F i g u r e 5.11: M a g n i t u d e of Force and Displacement  and is fractured towards the finger's t i p . T h e global behavior of the stiff vertex 98 and the compliant vertex 4 is also very different.  I n F i g u r e 5.13 the corresponding y - c o l u m n of the discrete Green's  functions m a t r i x is illustrated. T h e Green's function for vertex 98 near the m o u n t i n g surface of the a r m is localized i n the v i c i n i t y of the vertex. However, the Green's function for vertex 4 on the h a n d of the h u m a n wrist m o d e l has no noticeable local effect.  Tractions at this vertex 4 deform the complete h u m a n wrist m o d e l .  The  deformation is larger on the h a n d t h a n on the lower arm. T h i s is a result of the lever effect due to the m o u n t i n g of the object and the object's c a p a b i l i t y to hinge around its wrist. T h e physical wrist acts as a complex rotational joint causing points on the surface of the h a n d to describe an approximately circular trajectory. T h i s is approximated as a straight-line m o t i o n by the m o d e l .  99  (a) Solution Norm | | H L | | 2  (b) Residual of F i t | | P 3 f c T  fc  -  U || T  2  F i g u r e 5.12: L o c a l C o m p l i a n c e 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 EjT  fc  not being estimated and are shown i n black. T h e result shown is o b t a i n e d w i t h a m a x i m u m of 24 probes per vertex location and observations from six different viewpoints.  100  (a) Probe at Vertex 4 of Figure 5.10(a)  (b) Probe at Vertex 98 of Figure 5.10(b)  aim lift  i  f  f  mm (c) Probe at Vertex 4 of Figure 5.10(a)  (d) Probe at Vertex 98 of Figure 5.10(b)  F i g u r e 5.13: Discrete Green's Functions for Different M e s h R e s o l u t i o n Base mesh level 0 of the s u b d i v i s i o n hierarchy (5.13(a) and 5.13(b)) a n d level 1 (5.13(c) and 5.13(d)) are shown. T h e estimation for b o t h meshes is based on the same observations (see Section 5.4). A g a i n , the m a g n i t u d e of displacements are shown w i t h vertex coloring, where red indicates a large value while blue indicates the smallest value. B l a c k vertices are not estimated.  101  C h a p t e r  6  Inverse P r o b l e m Solution T h e discrete Green's functions m a t r i x estimation has been approached i n C h a p t e r 5 as a linear estimation p r o b l e m . T h i s ignored the u n d e r l y i n g structure of the Green's functions m a t r i x . T h e derivation given i n Section 2.3.2 is based on a discretization of the b o u n d a r y integral equation.  Therefore, the solution by weighted residuals  employed K e l v i n ' s fundamental solution for homogeneous,  isotropic, linear-elastic  materials. T h i s leads to a specific structure of the Green's functions m a t r i x w h i c h is exploited below i n 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. T h e assumption of an object made entirely out of a single homogeneous m a t e r i a l is l i m i t i n g and would be a poor a p p r o x i m a t i o n for the examples of the plush toy tiger and the anatomic wrist m o d e l i n Section 5.5. T h e plush toy is made out of cloth w i t h a filling and most likely the head is separate from the b o d y inside the cloth. T h e anatomic soft-tissue wrist m o d e l has an artificial bone structure w h i c h is i n itself inhomogeneous.  In addition, the soft-tissue is simulated  w i t h foam w i t h a separate skin layer. T h e 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  Boundary  Element  Method  for a n Inverse  Elastic  Problem  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 dependent 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 1 - v  ikn  J-  3  8-7rr  103  2  1 —v  dr dr  dr  di dj  dn  8irr  2  (6.1)  U  ^ '  X  j 7  7 ( 1 - 1 / ) Wirr'  (l-i/)167rr  X  '  T h e overall m a t r i x equation of the b o u n d a r y element formulation H u = G p is E q u a t i o n 6.3. Separating the geometric and material properties terms as i n E q u a t i o n 6.1 and E q u a t i o n 6.2 leads to E q u a t i o n 6.3 (Section 2.3.2 and A p p e n d i x A contain more detail).  V ^  1  — VA H  + T-^-HB +  1  —V  C  u =  7(1 - v)  w i t h the diagonal m a t r i x Cu =  ~ 3=1  ^ij)  -GA H 1  jz  TGB  7(1 - u)  (6.3)  2i/  1 «.4,, + . -//„. 1-1/ ^ 1 - 1 /  E q u a t i o n 6.3 and E q u a t i o n 2.1 form the basis for the inverse p r o b l e m solution. A s described i n Section 2.3.2, a given b o u n d a r y configuration enables one to reorder H u = G p into the Green's functions m a t r i x E q u a t i o n 2.1. E q u a t i o n 6.3 can be reordered i n the same way. However, the separation between m a t e r i a l and geometric properties is lost, since the matrices H and G depend nonlinearly on the Poisson's ratio v. T h e b o u n d a r y configuration employed i n the following is the same as described i n Section 5.1. T h e s u p p o r t i n g surface is fixed, i.e., a displacement of zero is described and the free surface has its t r a c t i o n prescribed. A s before, these values are combined i n the block vector v . T h e b o u n d a r y configuration remains fixed as different vertices on the free surface are probed. T h i s enables the c o m b i n a t i o n of several probings i n a c o m m o n estimation step since the Green's functions m a t r i x does not change.  For each probing, some part of the right-hand side vector v is  observed. A g a i n , the observation at the p r o b i n g location is the local measurement made w i t h the robotic a r m (see Section 3) while the global deformation is observed  104  by range-flow (as described i n Section 4). A n o p t i m i z a t i o n procedure to find shear m o d u l u s 7 a n d Poisson's ratio v from  = v is described i n the following  Section 6.1.2.  6.1.2  Optimization  T h e shear m o d u l u s 7 and Poisson's ratio v are found by an iterative o p t i m i z a t i o n . F i r s t , the surface integrals i n E q u a t i o n 6.3 need to be calculated. T h i s is a very expensive c o m p u t a t i o n (0(TNM)  w i t h T the number of triangles, N the number  of vertices a n d M the number of integration points w i t h i n 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 a n a l y t i c a l l y derived m o d e l of an object [58]. T h e objective function i n the m i n i m i z a t i o n is the error function n  F(7,^) = E ^ (  H  (7^)v-v)  2  .  (6.4)  A n element of the selection vector a i n E q u a t i o n 6.4 is one if observations of vertex displacement are made and zero otherwise. T h e objective function of E q u a t i o n 6.4 is evaluated starting from an i n i t i a l guess for the unknowns 7 and v. I n each iteration E q u a t i o n 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 w i t h E q u a t i o n 6.3. I n the next step the Green's function m a t r i x 3 ( 7 , v) has to be calculated from the fundamental solution matrices G and H . T h i s calculation involves a m a t r i x inverse i n the size of S (see Section 2.3.2). T h i s m a t r i x inverse enables one to arrive at rows for which all measurements are available because the complete described values v are k n o w n . T h e formulation is also t r i v i a l l y  105  expanded to m u l t i p l e probings i n the same b o u n d a r y configuration. E q u a t i o n 6.5 is the objective function for m measurements i n the same b o u n d a r y configuration.  n F  (l,  V) =  m Yl  [ 0 • • • rn]ij Q  a  ( S ( 7 , v) [ V i • • • V ] m  [  V i  • • • V ])?. m  (6.5)  i=0 j=0 T h e M a t l a b function fmins, an i m p l e m e n t a t i o n of the N e l d e r - M e a d simplex m e t h o d , has little difficulty solving the m i n i m i z a t i o n task at h a n d . However, note, that simpler methods should be sufficient because of the low d i m e n s i o n a l i t y of the p r o b l e m , the l i m i t e d range of possible values for Poisson's ratio v = 0.0 • • • 0.5 (in practice) and the shear m o d u l u s 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 p r o b l e m solution for a test object:  half of a soft ( N e r f ™ ) ball (see F i g u r e 4.3).  T h e estimation  results for the half ball are s u m m a r i z e d i n Table 6.1. T h e linear elastic m o d e l is a good a p p r o x i m a t i o n for the ball and the m o d e l fitting is successful. T h e estimated parameters v and 7, when used i n s i m u l a t i o n , produce results w h i c h closely resemble the recorded range-flow (see F i g u r e 6.1). T h e results i n Table 6.1 are from measured displacements by range-flow and displacement and force at a p r o b e d vertex.  Table 6.1 also contains a comparison  w i t h a destructive test. T h e test is a simple compression test of a c y l i n d r i c a l material sample. T h e material sample is compressed measuring force and the reduction in length of the cylinder. In practice, this works well for metal samples i n c o m m o n 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  E  1  0.1183  7 0.7809  2 3 4 5 6  0.1878 0.2227 0.2516 0.1857  0.6308 0.6039 0.6555 0.6977  0.2009  0.6696  1.641 • 1.655 1.608  8 9  0.2089 0.2007  11 12 14 15  0.1868 0.2082 0.1997 0.2084  0.6550 0.6646 0.7175 0.7585 0.7993 0.7451  1.584 1.596 1.703 1.833 1.918 1.801  17 18 20 22 26  0.1192 0.1366 0.1400 0.1264 0.1442  0.7666 0.7562 0.7480 0.7121 0.9383  1.716 1.719 -1.705 1.604 2.147  A L L the above  0.1730  0.7283  Test:  -  -  1.746 1.499 1.477  1.709 1.75  K,  Table 6.1: M a t e r i a l Properties E s t i m a t i o n Results T h e table lists the estimation results of the m a t e r i a l properties (shear modulus 7, Poisson's ratio v a n d the corresponding Y o u n g ' s m o d u l u s E) for the half-ball. E a c h row i n the table contains the result for a separate estimate based on the force probe location at the specified node(s). U s i n g d a t a from a l l probes i n a single estimation results i n 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 c y l i n drical sample of the b a l l is compressed i n a m i l l i n g machine. T h e m i l l i n g machine ensures that the cylinder is compressed i n the vertical direction only. T h e force is measured w i t h an electronic balance (Ohaus  CT 6000-S) under the sample. For a  compressed c y l i n d r i c a l m a t e r i a l piece the Y o u n g ' s m o d u l u s is E = -j^, where F is the compression force, I a n d <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 m o d u l u s 7 a n d Y o u n g ' s m o d u l u s E is 7 = ^(i+v) • T h i s allows the comparison i n T a b l e 6.1, w h i c h shows an excellent correspondence between the material-based compression test a n d the object-based inverse p r o b l e m m e t h o d . T h e v a r i a t i o n of the estimates between single probes is mostly likely due to the misalignment of the probe w i t h respect to the surface. C o m b i n i n g the measurement i n one estimation step averages out the noise nicely. A visual comparison between the measured displacement and the simulated displacement for a p r o b i n g is shown i n F i g u r e 6.1. T h e m a g n i t u d e of surface displacement shows very good correspondence between the s i m u l a t i o n 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 m o d e l . A s previously stated, since these assumption are fairly strong, relaxing the assumptions is desirable.  O n e approach to overcome  the l i m i t a t i o n of homogeneity w o u l d be to t r y to estimate the location of interior boundaries between materials at the same time as the m a t e r i a l properties. 108  This  I F i g u r e 6.1: C o m p a r i s o n of M e a s u r e d a n d S i m u l a t e d D i s p l a c e m e n t F r o m left to right: 3 D range-flow, estimated displacement, simulated Displacement. T h e 3 D range-flow shows the d i s t r i b u t i o n of the estimated flow vectors w i t h brighter range-flow corresponding to larger motions. T h e node displacement figure also shows larger movement being brighter. D a r k surface is either s t a t i o n a r y or no estimate is available. T h e estimated node displacement forms a p a r t i a l r i n g a r o u n d the location of the probe; the location itself and the other part of the r i n g are occluded by the probe.  obviously w o u l d be a formidable task but still quite restrictive. A d d i t i o n a l l y , this interior interface between materials is not required for the s i m u l a t i o n , o n l y the Green's functions m a t r i x is necessary. Instead of e s t i m a t i n g the location of an interface, the response of the object to p r o b i n g at vertex k can be modeled approximately w i t h an apparent shear m o d u l u s 7k and an apparent Poisson's r a t i o n v\.. T h i s allows for v a r y i n g stiffness over the object surface a n d v a r i a t i o n i n how global the object responds at a l o c a t i o n on the surface. It is l i m i t e d to object responses which remain continuous over the object's surface and are isotropic. T h e l i m i t a t i o n s of the homogeneous isotropic linear elastostatic b o u n d a r y element m o d e l are i l l u s t r a t e d next w i t h the example of the a n a t o m i c soft-tissue h u m a n wrist m o d e l . A n apparent shear m o d u l u s 7^ a n d an apparent P o i s s o n ratio v  k  are estimated for vertices k = 4 and k = 98 w h i c h are shown i n F i g u r e 5.10. T h e  estimates are based o n 8 probes per vertex and only the local response is utilized. T h e estimates are given i n the following table.  109  Vertex N o . . 4 98  • Ik 0.04658 ' 0.0  19.28 1.157  T h i s approximate m o d e l 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 a n d instead the l i m i t of vgg, = 0 is selected ( T h i s is found w i t h constraint o p t i m i z a t i o n w i t h the M a t l a b function fmincon). T h e apparent shear m o d u l i £ are drastically different at the two locations. 7  T h e lower a r m deforms s m o o t h l y w h e n p r o b e d at V e r t e x 98 a n d the object's response is quite local w i t h a cusp at the contact location (see F i g u r e 6.2).  The  response has also no preferred direction. However, homogeneous isotropic m a t e r i a l cannot produce this response given the b o u n d a r y configuration: T h i s response is caused by the soft tissue of the object being compressed between the probe a n d the object's bone structure. Figures 6.2(a) to 6.2(b) compare the result of the p a r a m eter e s t i m a t i o n w i t h the response estimated by the Green's function e s t i m a t i o n of C h a p t e r 5. F i g u r e 6.2(c) shows the m a g n i t u d e of the difference between the two methods for the y - c o l u m n of the Green's functions block c o l u m n Sj,98However, some of the response of the h u m a n wrist m o d e l above the wrist joint can be a p p r o x i m a t e d . T h e bone structure is not rigid enough to restrict the deformation a n d the complete object bends.  Figures 6.2(d) and 6.2(e) show how  b o t h estimation m e t h o d lead to a global b e n d i n g of the object. T h e m a g n i t u d e of the difference between the estimated displacements, shown i n 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. between 6.2(a) and 6.2(b)  (d) Approx. Inverse Problem  (e) Green's Functions Estimation  (f) Dif. between 6.2(d) and 6.2(e)  F i g u r e 6.2: Inverse P r o b l e m Solution and Discrete Green's F u n c t i o n s E s t i m a t i o n T h e Green's functions estimation m e t h o d of C h a p t e r 5 can represent deformations w h i c h are outside the d o m a i n of homogeneous isotropic linear elastostatic models. Figures 6.2(a), 6.2(b) a n d 6.2(c) show the behavior at vertex 98 of F i g u r e 5.10(a), while 6.2(d), 6.2(e) a n d 6.2(f) show the comparison for vertex 4 of F i g u r e 5.10(b). T h e magnitude of displacements are shown w i t h vertex coloring, where red indicates a large value while blue indicates the smallest value. B l a c k vertices are not estimated.  Ill  6.3  Summary  T h i s section shows how the elasticity of an object can be estimated by solving an inverse p r o b l e m for the elastic constants of the object. T h e m e t h o d is restricted to objects made of isotropic, homogeneous, linear-elastic material. For objects w h i c h fit this assumption, the model works well as shown for the test object. T h e m a t e r i a l restriction may be relaxed if the m o d e l is understood to be local and approximate only. A local approximate m o d e l based on linear elasticity may be employed by assigning different elastic constants at each location on the surface. I n fact, the elastic constants become elastic functions over the surface.  However, internal structure  and a r t i c u l a t i o n cannot be modeled i n this way. B u t the discrete Green's function estimation of C h a p t e r 5 still provides a reasonable a p p r o x i m a t i o n i n the case of the anatomic soft tissue h u m a n wrist m o d e l . A n explicit m o d e l of articulated soft bodies is perhaps more promising and is the subject of future work.  112  Chapter 7  Conclusions and Future W o r k  7.1  Conclusions  T h i s thesis presents a m e t h o d to scan physical deformation behavior of objects. A m o d e l of the deformation behavior is estimated based on measurements of an intact existing object. A method to acquire these measurements i n c l u d i n g a measurement of the global deformation response is described. T h e global deformation is observed w i t h a computer vision technique.  T h e robotic measurement  response facility  A C M E is extended to enable the d a t a acquisition. T h e same techniques also allow one to validate an estimated deformable model. A C M E serves to acquire a geometric m o d e l of an object, to a p p l y contact forces to the m o d e l and to record the deformation of the object.  A C M E allows  the acquisition of a large amount of deformation d a t a i n an automatic fashion, i.e., A C M E can scan physical deformation behavior of objects.  T h e d a t a acquisition  d u r i n g the scan is based on two devices: a force sensor equipped robotic a r m and a t r i n o c u l a r stereo-vision system. T h e robotic a r m senses the p o s i t i o n and forces at the contact point w i t h 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 c o m b i n a t i o n of t r i n o c u l a r optical flow a n d stereo. T h e a d d i t i o n a l images i n the calculation of the o p t i c a l flow enable a consistency check of the result. T h e c o m b i n a t i o n of optical flow a n d stereo yields three-dimensional range flow. T h e registration w i t h a separately acquired surface mesh of the object leads to estimates of the three dimensional vertex displacements d u r i n g deformation. T h e acquired d a t a can be utilized to estimate a n d validate deformable models. T h i s thesis develops and implements two estimation approaches.  O n e ap-  proach m a y be characterized as a d a t a driven approach w i t h the assumption of linear elastostatic behavior of object. T h i s allows one to fit a deformable m o d e l to a wide range of objects but requires a large amount of data. T h i s requirement of extensive d a t a acquisition is m i t i g a t e d 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 m e t h o d requires little d a t a because of this strong assumption. Nevertheless, it describes certain types of objects well; w h i c h are made of a single suitable material for instance an object made of closed-cell foam.  M o d e l i n g assumptions i n between these two extremes  are not investigated but are suggested. fitting  T h e thesis contains an example of m o d e l  for a p l u s h toy w h i c h has been demonstrated successfully to audiences at  several scientific conferences (11th I R I S P r e c a r n A n n u a l Conference, J u n e 2001, O t tawa, C a n a d a ; 8th Int. Conf. of C o m p u t e r V i s i o n , J u l y 2001, Vancouver, C a n a d a ) . A n o t h e r example discussed is an anatomic soft tissue h u m a n wrist m o d e l . A deformation m o d e l has also been estimated for a simple homogeneous foam object. T o 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  O n e goal for future work could be to embed the presented techniques into a probabilistic framework based on a measurement theory. T h i s thesis characterizes measurements by a comparison w i t h g r o u n d - t r u t h i n the case of the displacement measurement i n C h a p t e r 4 a n d w i t h residual of fit i n the case of the estimation methods i n C h a p t e r s 5 and 6. measurements.  T h i s allows one to assess how well a m o d e l explains the  I n order to be able acquire a measurement a n d make a quantita-  tive statement if this measurement is explained by a given m o d e l , a probabilistic framework w i t h i n a measurement theory is necessary. A measurement theory establishes the relations between observations and a physical m o d e l . T h e theory has to m o d e l the measurement process, reasoning about the influence of various physical phenomena - some i n accordance w i t h the m o d e l and some not captured by the m o d e l . T h e latter are generally referred to as noise. Bayesian estimation theory is one framework for quantitative statements about m o d e l fit a n d m o d e l selection. T h e other suggestions for future work concern e s t i m a t i n g d y n a m i c models of deformation. T h e trend i n deformable models for m e d i c a l s i m u l a t i o n is towards d y n a m i c non-linear finite element models, despite the fact that current computers can render only very s m a l l a n d very soft models i n real-time. D y n a m i c model  fitting  w i t h a large deformation s t r a i n tensor would require the t r a c k i n g of objects d u r i n g the deformation process at least at twice the m a x i m u m o c c u r r i n g frequency. T h e techniques of this thesis still a p p l y but the image acquisition w i t h the Triclops stereo-vision system would most likely be too slow for the deformation transient. Higher speed camera systems a n d strobe lights could solve this p r o b l e m .  115  A n o t h e r o p t i o n w o u l d be fitting articulated deformable models. T h i s p r o b l e m is usually addressed w i t h a p r i o r i k n o w n k i n e m a t i c m o d e l . T h e m o d e l of articulation could be fit while m i n i m i z i n g the deformation of i n d i v i d u a l links. In humans the area around the joints exhibit large deformation (e.g., finger joints) w h i c h w o u l d be one area to address.  116  A p p e n d i x  A  Linear Triangular Elements T h i s A p p e n d i x gives the necessary notation to be able to evaluate integrals based on the fundamental solutions i n E q u a t i o n s 2.16, 2.17, 6.1 a n d 6.2. T h e integrals are evaluated w i t h linear triangular elements, i.e., the shape functions over a t r i a n g u lar b o u n d a r y element are linear i n the values at the vertices of the triangle. T h e fundamental solutions contain the distance r between load location a n d response location. A component of the distance vector Ti given the vertices x ... l  x  3  of the  triangular element i n anti-clockwise order can be calculated as  r = x \ - mx] - r] x - (1 - m - r) )x 2  l  2  3  2  T h e distance is r = \JYA=\ i, while the p a r t i a l derivative of the distance is r  and | ^ =  = f r  " ^ v - A component of the unit n o r m a l is ??,; = gi/\G\ (see E q u a 1  tion A . 7 ) . N e x t the mappings from local coordinates of an element to global coordinates and vice versa are given. E q u a t i o n A . l maps local homogeneous surface coordinates rji (i.e., 773 = 1 — 771—772) of a triangular shape element to global Cartesian coordinates Xi. T h e vertices of the triangular element are labeled x\ ... x\ i n anti-clockwise order.  117  X = 1-771 -I- x\f]2 + X 773  (A-l)  3  l  E q u a t i o n s A . 2 , A . 3 and A . 4 m a p global to local coordinates.  Vi  ^2^3 1 2  I  1 /yi 2 ™3 ™2 1 2 1 ' ™ 1 1™3 2 ™3™ 1 2 ' 1 2 ry  m  ™3™2 1 2  '  _2  1 I  /-*»  1 /*^3  X X2  x x2  2 3  — arxf  1  /y 1 /y>3  i  xj-^2  1  2  x  x  2  3  1  2  — orxf 3 \  — x x — x^x ~j~  (A.2)  3  1  /T/-2 ™3 — X-yX^~\~ X - ^ X 2 2 2 + l x2x\ + x x| x  X2  x x\  x\x\ ~\~ Xrj l 2 i i + a^Xr, — x x\ + x x x  3  X  +  —  rp2> ™3  x  1  rv» 1 2 3 2 1 2 ' 1 2 1 2 ' 1 2 /-p 1 ^»--2 i ™3~.2 1 2 1 2  ~f"  (A-3)  2  X^X2  (A-4)  2  T h e displacement and traction vector over an element are expressed i n E q u a tions A . 5 and A . 6 , respectively, using identical interpolation functions as the above shape functions, i.e., (fi = rji.  (A.5>  Uj  (A.6)  T h e change to local coordinates requires the J a c o b i a n i n E q u a t i o n A . 7 .  dF  —  with | G |  =  \G\drjidrj2  (A.7)  i=l 9i 92 9i  3x2 9x  3  dm  8x2 dx  3  dr)  dm  dm  9x3 <9xi  dx  dxi  dm  dm  2  <9T?I <9?72  3  dxi  8x2  dxi dx2  dm  dr]2  dr/2  118  dm  W i t h the local coordinates of E q u a t i o n A . l , the partials i n E q u a t i o n A . 7 become:  After the mappings between local a n d global coordinates have been established, it remains to evaluate the integrals. T h e fundamental solution integrals i n E q u a t i o n s 2.18, 2.19 a n d E q u a t i o n 6.3 are evaluated n u m e r i c a l l y by G a u s s i a n quadrature.  E m p l o y i n g L integration steps  per element weighted by a factor wi leads from E q u a t i o n 2.19 to A . 8 , ( T h e other integrals are converted identically.)  L  H  f c m  = /  p*(Sk,x)<p"m{x)dr(x)  « E\  G  \ i i ( p t m <Pm(vi,m)i) w  (A.8)  1=1  T h e integration points a n d the corresponding weights for the G a u s s i a n quadrature can be found i n many references, e.g., [33].  C a r e has to be taken w i t h the  singular integrals for the self-effects (r —> 0 for k = m).  119  A p p e n d i x  B  Geometric Modeling with A C M E T h i s a p p e n d i x details the i n i t i a l geometric m o d e l acquisition briefly i n t r o d u c e d i n Section 3.3. The  m a i n shape sensor i n A C M E is the t r i n o c u l a r stereo-head w i t h i n the  F M S . T h e sensor produces large amounts of range d a t a i n close to real time employing its commercial stereo vision S D K . However, range d a t a from stereo is notoriously noisy. A n o t h e r feature of stereo ranging is its dependence o n image features for m a t c h i n g between the stereo cameras. A d d i t i o n a l features can be projected onto a surface w i t h a light or laser system.  Color Triclops  produces d a t a of adequate but  modest resolution (for a t y p i c a l viewpoint i n A C M E a depth resolution of 2-3 m m is achieved). Geometric object m o d e l i n g requires that the entire surface of the object is covered, w h i c h requires varying viewpoints. R e g i s t r a t i o n of v i e w p o i n t s i n A C M E is achieved by visual calibration of the F M S w i t h respect to an optical target, as well  120  TestStation  as by e x p l o i t i n g the accuracy of motions of the  and the repeatability  of m o t i o n commands for the F M S . T h e registration of the viewpoint relative to  Iterative Closest Point  the target is accomplished w i t h a variant of an  ( I C P ) [8]  a l g o r i t h m . T h e starting point for the m i n i m i z a t i o n is derived from the kinematics of A C M E . Range images from different viewpoints are combined i n a voxel-based volumetric approach w i t h reconstruction software p r o v i d e d by the N a t i o n a l Research C o u n c i l of C a n a d a [109]. I n order for this approach to succeed the d a t a not only has to be i n a c o m m o n coordinate frame but also needs to be noise-filtered. Furthermore, normals for the p o i n t - c l o u d have to be estimated. T h e range d a t a is filtered b y calculating range w i t h variable mask sizes and c o m b i n i n g t h e m i n a v o t i n g scheme. A n example of a p o i n t - c l o u d before and after filtering  is shown i n F i g u r e B . l . T h e threshold for the  filtering  is set based on  experience and is far from universal. However, since the stereo-head is quite c o m m o n throughout the machine v is io n research community, the processing strategy and actual setting of thresholds m a y be of interest to some and are reported here: T h e filtering is achieved by processing a stereo image triplet twice. One stereo range image is calculated to m a x i m i z e the accuracy of the result by using a small mask, while the other image is calculated to m i n i m i z e noise w i t h a large mask and tight validation settings. O n e stereo range image is calculated w i t h a small stereo correlation mask of five i n a 640 x 480 image. vision S D K ' s (The  uniqueness threshold  uniqueness threshold  A d d i t i o n a l l y , the  Triclops  stereo  is kept nearly ineffective w i t h a setting of 2.0  is intended to filter out matches w h i c h are not consistent  between the right-left and the right-up camera p a i r ) .  T h e second stereo range  image is calculated w i t h a m u c h larger correlation mask size of 13 and the tight validation setting of an  uniqueness threshold  121  of 1.0 combined w i t h the connected  (a) Before Filtering.  (b) After Filtering.  (c) Before Filtering.  (d) After Filtering.  Figure B . l : P o i n t - C l o u d F i l t e r i n g 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 a n a t o n o m i c soft-tissue h u m a n wrist m o d e l . (See Text for details).  122  region threshold of 300 pixels w i t h less t h a n 0.75 d i s p a r i t y difference. T h e images are combined i n a plane fitting step simultaneously estimating normals a n d v a l i d a t i n g a point. A l l points are rejected w h i c h are not w i t h i n 1cm between the two range images. T h e plane fit for each point i n 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 w i t h i n the neighborhood have a valid depth. T h e plane fit is characterized by the ratio between smallest singular value over the s u m of the two larger singular values. F i t t i n g is done i n metric x — y — z space a n d a fit w i t h a singular value ratio of not more t h a n 0.005 is accepted. T h e angles between the viewing axes to the point and the estimated n o r m a l is checked a n d the point are kept if the angle is less t h a n 80 degrees. T h e filtered range images are i n p u t i n the mesh generation software a n d a high-resolution mesh is generated. T h e resolution of the triangular mesh depends on the s a m p l i n g density of the object's surface. T h e produced mesh contains t y p i c a l l y some triangles due to support fixtures, some due to noise, and is not watertight. A multi-resolution s u b d i v i s i o n mesh is built from this high-resolution raw mesh. I n the process, holes are filled, triangles edited a n d the m e t h o d for mesh simplification and subsequent s u b d i v i s i o n refinement eliminates irregularities of the mesh.  The  simplification is done w i t h Q S l i m [38]. T h e r e m a i n i n g process is essentially a displacement s u b d i v i s i o n surface [72] or n o r m a l mesh [47] construction process. T h e implementation of this subdivision construction process is due to D . James a n d is described w i t h more detail i n [95]. T h e result for a plush toy tiger is shown i n F i g u r e 3.8. T h e stereo system also provides calibrated color imagery. These images are used for very simple texture m a p p i n g of the object. T h e vertices of the s u b d i v i s i o n  123  model are projected into the color images w i t h the pinhole camera m o d e l . A triangle of the mesh is assumed to be visible from a camera if a l l its vertices and its centroid are visible. V i s i b i l i t y is calculated by ray intersection w i t h the triangular mesh. T h e intersection test is accelerated by an octree space division. T h e triangle is t y p i c a l l y visible from a number of views but w i t h different foreshortening and resolution. O u t of the images where the triangle is visible, the.one w i t h the highest v i e w i n g quality measure is selected for texture m a p p i n g . T h e viewing quality measure used is the p r o d u c t of image area covered by the triangle times the cosine between view direction and triangle n o r m a l . T h e relative global brightness of a l l these color images is also adjusted.  T h e produced texture m a p (see F i g u r e 3.8(f)) contains m a n y  discontinuities due to inaccuracies of the mesh, misregistration and v a r y i n g local brightness.  T h e texture maps could be readily improved by blending of the local  textures over the entire m a p e l i m i n a t i n g abrupt brightness changes at edges and m i n i m i z i n g the effect of geometric errors, as i n [108].  124  Bibliography [1] M . A k a m a t s u , G . N a k a m u r a , a n d S. Steinberg. Identification of L a m e coefficients from b o u n d a r y observations.  Inverse Problems,  7(3):335-354, 1991.  [2] P . A n a n d a n . A c o m p u t a t i o n a l framework and an a l g o r i t h m for the measurement of visual m o t i o n .  International Journal of Computer Vision,  2:283-310,  1989. [3] D . H . B a l l a r d and O . A . K i m b a l l . flow.  R i g i d b o d y m o t i o n from d e p t h and optic  Computer Vision. Graphics and Image Processing,  [4] M . R . B a n a n , M . R . B a n a n , a n d K . D . Hjelmstad.  Parameter estimation of  structures form static response. I. C o m p u t a t i o n a l aspects.  tural Engineering,  120(ll):3243-3258, 1994.  [5] M . R . B a n a n , M . R . B a n a n , and K . D . Hjelmstad.  Journal of Struc-  . Parameter estimation of  structures form static response. II. N u m e r i c a l simulation studies.  Structural Engineering,  IJCV,  Journal of  120(11):3259-3283, 1994.  [6] J . L . B a r r o n , D . J . Fleet, and S.S. Beauchemin. techniques.  22(1):95-115, 1983.  12(l):43-77, 1994.  125  Performance of optical flow  [7] C . Basdogan, C . - H . H o , M . A . Srinivasan, S . D . S m a l l , and S . L . D a w s o n . Force interactions i n laparascopic simulations: h a p t i c rendering of soft tissues. In  Medicine Meets Virtual Reality, [8] P . J . B e s l and N . D . M c K a y .  pages 385-391, 1998.  A m e t h o d for registration of 3-d shapes.  IEEE  Transactions on Pattern Recognition and Machine Intelligence,  14(2):239-256,  Feb 1992. [9] L . M . B e z e r r a and S. Saigal. A boundary-element formulation for the inverse elastostatics p r o b l e m ( I E S P ) of flaw detection.  merical Methods in Engineering, [10] L . M . B e z e r r a and S. Saigal. the B E M .  International Journal for Nu-'  36(13):2189-2202, 1993.  Inverse b o u n d a r y t r a c t i o n reconstruction w i t h  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 t h r o u g h 3-  dimensional soft tissue. In P . B r u n e t and R . Scopigno, editors,  Eurographics,  C o m p u t e r G r a p h i c s F o r u m , 18(3), pages C 3 1 - C 3 8 , M i l a n o , Italy, 1999. [12] J . Biggs and M . A . Srinivasan.  Handbook of Virtual Environment Technology,  chapter H a p t i c interfaces, K . M . Stanney (editor). Lawrence E r l b a u m Associates, i n press. [13] M . J . B l a c k and P . A n a n d a n .  T h e robust estimation of m u l t i p l e motions:  parametric and piecewise-smooth flow  Understanding,  63(1):75-104, 1996.  126  fields.  Computer Vision and Image  [14] R . C . Bolles and M . A . Fischler. A R A N S A C - b a s e d approach to model and its application to finding cylinders i n range data. In  Art. Intelligence,  Int. Joint Conf. on  pages 637-643, Vancouver, C a n a d a , 1981.  [15] S. D i B o n a and O . Salvetti. simulation. In  fitting  A deformation m o d e l for biotissues behavior  International Conference on Image Processing,  volume 2, pages  443-446, Vancouver, C a n a d a , Sept 2000. [16] C A . B r e b b i a .  Boundary element techniques in engineering.  Pentech Press,  L o n d o n , U k , 1978. [17] C A . B r e b b i a , J . C . F . Telles, and L . C . W r o b e l .  Boundary element techniques:  theory and applications in engineering. [18] M . Bro-Nielsen and S. C o t i n .  Springer-Verlag, B e r l i n , 1984.  Real-time volumetric deformable models for  surgery simulation using finite elements and condensation. and F . X . Sillion, editors,  Eurographics,  In J . Rossignac  C o m p u t e r G r a p h i c s F o r u m , 15(3),  pages 57-66, Poitiers, France, 1996. [19] I. Brouwer, J . U s t i n , L . Bentley, A . Sherman, N . D h r u v , and F. Tendick. M e a s u r i n g i n vivo a n i m a l soft tissue properties for h a p t i c modeling i n surgical simulation. In  Medicine Meets Virtual Reality,  pages 69-74, J a n 2001.  [20] R . Carceroni and K . K u t u l a k o s . M u l t i - v i e w scene capture by surfel sampling: from video streams to non-rigid 3d motion, shape and reflectance. In  national Conference on Computer Vision,  Inter-  volume 2, pages 60-67, 2001.  [21] S. C h a u d h u r i and S. Chatterjee. E s t i m a t i o n of m o t i o n parameters for a deformable object from range data. In  tion,  pages 291^295. I E E E , 1989.  Computer Vision and Pattern Recogni-  [22] S.S. C h e n and M . A . P e n n a . Shape and m o t i o n of n o n r i g i d bodies.  Vision, Graphics and Image Processing,  Computer  36(2/3):175-207, 1986.  [23] L F . C o s t a and R . B a l a n i u k . L e m - an approach for real-time physically based soft tissue simulation. In  tion,  International Conference on Robotics and Automa-  pages 2337-2343, Seoul, S o u t h K o r e a , M a y 2001.  [24] S. C o t i n , H . Delingette, and N . A y a c h e . R e a l - t i m e elastic deformations of soft tissues for surgery simualtion.  Graphics,  IEEE Trans, on Visulization and Computer  5 ( l ) : 6 2 - 7 3 , 1999.  [25] S. C o t i n , H . Delingette, and N . Ayache.  A h y b r i d elastic m o d e l allowing  real-time c u t t i n g , deformations and force-feedback for surgery t r a i n i n g and simulation.  The Visual Computer,  16(8):437-452, 2000.  [26] D . d ' A u l i g n a c , R . B a l a n i u k , and C . Laugier. t u a l exam of the h u m a n t h i g h . In  Automation,  A haptic interface for a vir-  International Conference on Robotics and  pages 2452-2457, San Francisco, U S A , A p r i l 2000.  [27] D . d ' A u l i g n a c , M . C . Cavusoglu, and C . Laugier. M o d e l i n g the d y n a m i c s of the h u m a n t h i g h for a realistic echographic simulator w i t h force feedback. I n Int.  Conf. on Medical Image Computer-Assisted Intervention,  pages 1191-1198,  C a m b r i d g e , U k , Sept 1999. [28] S. D e and M . A . Srinivasan.  T h i n walled models for haptic and graphical  rendering of soft tissues i n surgical simulations. In  . Reality,  pages 94-99, 1999.  128  Medicine Meets Virtual  [29]: F . B o u x de Casson, D . d ' A u l i g n a c , and C . Laugier. A n interactive m o d e l of the h u m a n liver. In 7th  International Symposium on Experimental Robotics,  H o n o l u l u , U S A , December 2000. [30] G . D e b u n n e , M . D e s b r u n , M . - P . C a n i , and A . H . B a r r .  D y n a m i c real-time  Computer Graphics,  deformations using space &; time adaptive s a m p l i n g . In  Annual Conference Series,  pages 31-36, L o s Angles, A u g 2001. A C M S I G -  GRAPH. [31] H . Delingette.  T o w a r d realistic soft-tissue m o d e l i n g i n medical s i m u l a t i o n .  Proceedings of the IEEE,  86(3):512-523, 1998.  [32] M . D e s b r u n , M . Meyer, P . Schroder, and A . H . B a r r . regular meshes using diffusion and curvature flow.  Annual Conference Series,  In  Implicit fairing of ir-  Computer Graphics,  pages 317-324, L o s Angles, U S A , A u g 1999. A C M  SIGGRAPH. [33] J . D o m i n g u e z .  Boundary Elements in Dynamics.  C o m p u t a t i o n a l Mechanics  P u b l i c a t i o n s , S o u t h h a m p t o m , U k , 1993. [34] J . S . D u n c a n , R . L . Owen, L . H . Staib, and P . A n a n d a n . Measurement of nonrigid m o t i o n using contour shape descriptors. In  Recognition, [35] P . F u a .  Computer Vision and Pattern  pages 318-324. I E E E , 1991.  A parallel stereo a l g o r i t h m that produces dense depth maps and  preserves image features.  Machine Vision and Applications,  [36] F . Ganovelli, P . C i g n o n i , C . M o n t a n i , and R Scopigno.  6:35-49, 1993.  A multiresolution  model for soft objects s u p p o r t i n g interactive cuts and lacerations. In M . Gross  129  and F . R . A . H o p g o o d , editors,  Eurographics,  C o m p u t e r G r a p h i c s F o r u m , 19(3),  pages C 2 7 1 - C 2 8 2 , Interlaken, S w i t z e r l a n d , 2000. [37] L . G a o , K . J . Parker, R . M . Lerner, and S . F . L e v i n s o n . Imaging of the elastic properties of tissue - a review.  Ultrasound in Medicine and Biology,  22(8) :959-  977, 1996. [38] M . G a r l a n d and P . S . Heckbert. . metrics.  In  Surface simplification using quadric error  Computer Graphics, Annual Conference Series,  pages 209-216,  Los Angles, U S A , A u g 1997. A C M . [39] J . C . G e l i n and O . G h o u a t i . Inverse m e t h o d for material parameters estimation i n the inelastic range.  Computational-Mechanics,  16(3):143-150, 1995.  [40] S . F . F . G i b s o n and B . M i r t i c h . A survey of deformable m o d e l i n g i n computer graphics. Technical R e p o r t T R - 9 7 - 1 9 , M i t s u b i s h i E l e c t r i c Research L a b o r a tory, C a m b r i d g e , M A , U S A , N o v . 1997. [41] C . G i o d a . Indirect identification of the average elastic characteristics of rock masses. I n  Structural foundations on rock,  volume 1, pages 65-75, 1980.  [42] S. Gottschalk, M . L i n , and D . M a n o c h a . O B B - t r e e : A hierarchical structure for r a p i d interference detection. I n  Series,  Computer Graphics, Annual Conference  pages 171-180, N e w 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 . C o m p u t a t i o n a l methods for inverse finite elastostatics.  Computer Methods in Applied Mechanics and Engineering,  2):47-57, 1996.  130  136(1-  [44] S. Govindjee and P . A . M i h a l i c . C o m p u t a t i o n a l methods for inverse deformations i n quasi-incompressible finite elasticity.  in Engineering,  43:821-838, 1998.  [45] M . Grediac, C . T o u k o u r o u , and A . V a u t r i n . of structures:  Int. Journal for Numer. Meth.  Inverse p r o b l e m i n mechanics  a new approach based on displacement field processing.  M . T a n a k a and H . D . B u i , editors,  IUATM Symposium,  In  Inverse P r o b l e m s i n  E n g i n e e r i n g Mechanics, B e r l i n , Germany, 1993. Springer-Verlag. [46] M . G r e d i a c and A . V a u t r i n .  A new method for determination of bending  rigidities of t h i n anisotropic plates.  Transactions of the ASME,  pages 964-  968, 1990. [47] I. Guskov, K . V i d i m c e , W . Sweldens, and P . Schroder. N o r m a l meshes. In  Computer Graphics. Annual Conference Series,  pages 95-102, N e w Orleans,  U S A , J u l 2000. A C M S I G G R A P H . [48] P . C . Hansen. R e g u l a r i z a t i o n tools: A M a t l a b package for analysis and solution of discrete ill-posed problems.  Numerical Algorithms,  6(1/2):1—35, 1994.  [49] J . D . H e l m , M . A . S u t t o n , and S . R . N e i l . Three-dimensional image correlation for surface displacement measurements. In S. F . E l - H a k i m , editor,  rics III, volume  Videomet-  2350, pages 32-45. S P I E , 1994.  [50] K . H i r o t a and T . K a n e k o . Representation of soft objects i n v i r t u a l environment. In  8th Int. Conf. on Artificial Reality and Tele-Existence (ICAT),  59-62, Tokyo, J p , Dec 1998. [51] B . K . P . H o r n and J . G . H a r r i s . R i g i d b o d y m o t i o n from range image sequences.  CVGIP: Image Understanding,  53(1):1-13, 1991.  pages  [52] B . K . P . H o r n and B . G . Schunk. D e t e r m i n i n g optical flow.  gence,  Artificial Intelli-  17(l-3):185-203, 1981.  [53] B . H o r o w i t z and A . P e n t l a n d . Recovery of non-rigid m o t i o n 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 . H u a n g and S . D . B l o s t e i n . R o b u s t algorithms for m o t i o n estimation based on two sequential stereo image pairs. In  nition,  Computer Vision and Pattern Recog-  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 Equil rium.  P h D thesis, U n i v e r s i t y of B r i t i s h C o l u m b i a , expected November 2001.  [57] D . James and D . K . P a i . A unified treatment of elastostatic contact simulation  Haptics-e,  for real time haptics.  [58] D . L . James and D . K . P a i . In  2(1), 2001.  A r t D e f o accurate real time deformable objects.  Computer Graphics, Annual Conference Series,  U S A , A u g 1999. A C M  pages 65-72, Los Angles,  SIGGRAPH:  [59] A . J o u k h a d a r , C . B a r d , and C . Laugier. P l a n n i n g dextrous operations using physical models. In  International Conference on Robotics and Automation,  pages 748-753. I E E E , 1994. [60] A . J o u k h a d a r , F . G a r a t , and C . Laugier. Parameter identification for d y n a m i c simulation. In  International Conference on Robotics and Automation,  1928-1933, A l b u q u e r q u e , N M , U S A , 1997. I E E E .  pages  [61] F . K a l l e l and M . B e r t r a n d . Tissue elasticity reconstruction using linear perturbation method.  IEEE'  Transactions  on Medical Imaging,  15(3):299-313,  1996. [62] C . K a m b h a m e t t u , D . B . Goldgof, and M . He. D e t e r m i n a t i o n of m o t i o n p a r a m eters and estimation of point correspondences i n small n o n r i g i d deformations. In Computer  Vision  and Pattern  Recognition,  pages 943-946. I E E E , 1994.  [63] T . K a n a d e , A . Y o s h i d a , K . O d a , H . K a n o , and M . T a n a k a . A stereo machine for video-rate dense depth m a p p i n g and its new a p p l i c a t i o n . In Vision  and Pattern  Recognition,  Computer  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 m o d e l i n g of facial tissue for craniofacial surgery s i m u l a t i o n . Computer  Aided Surgery, 3(5):228-  238, 1998. [65] R . M . K o c h , M . H . Gross, F . R . Carls, D . F . von B i i r e n , G . Fankhauser, and Y . I . H . P a r i s h . S i m u l a t i n g facial surgery using finite element models. In Computer Graphics, Annual A u g 1996. A C M  Conference  Series, pages 421-428, N e w Orleans, U S A ,  SIGGRAPH.  [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 . K i i h n a p f e l , H . K . C a k m a k , and H . Maafi. E n d o s c o p i c surgery t r a i n i n g using v i r t u a l reality and deformable tissue s i m a u l a t i o n .  Computers  &  Graphics,  24:671-682, 2000. [68] S . - H . L a i and B . C . V e m u r i . Reliable and efficient c o m p u t a t i o n of o p t i c a l flow. International  Journal  of Computer  Vision, 29(2):87-105, 1998. 133  [69] J . L a n g . C o n t o u r stereo m a t c h i n g by correlation i n a constraint satisfaction approach. In  Vision Interface,  [70] J . L a n g and D . K . P a i .  pages 19-26,  Bayesian estimation of distance and surface n o r m a l  w i t h a time-of-flight laser rangefinder.  . Digital Imaging and Modeling. [71] J . L a n g and D . K . P a i . In 3rd  1998.  In  International Conference on 3-D  I E E E , October  1999.  E s t i m a t i o n of elastic constants from 3 D range-flow.  International Conference on 3-D Digital Imaging and Modeling,  331-338, Quebec City, C a n a d a , 2001.  pages  IEEE.  [72] A . W . F . Lee, H . M o r e t o n , and H . H o p p e . D i s p l a c e d s u b d i v i s i o n surfaces.  Computer Graphics, Annual Conference Series,  In  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. Realistic m o d e l i n g of facial a n i m a t i o n . 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 i a o , S.J. A g g a r w a l , and J . K . A g g a r w a l . T h e reconstruction of d y n a m i c 3 D structure of biological objects using stereo microscope images.  Vision and Applications,  9(4):166-178,  1997.  [75] J . L l o y d and V . H a y w a r d . Trajectory generation in m u l t i - R C C L .  Robotic Systems,  10(3):369-390,  Journal of  1993.  [76] J . E . L l o y d , J.S. Beis, D . K . P a i , and D . G . Lowe. w i t h vision. In  Machine  Model-based telerobotics  International Conference on Robotics and Automation,  1297-1304. I E E E ,  1997.  134  pages  [77] J . Louchet, X . Provot, and D Crochemore. cloth a n i m a t i o n models. In  E v o l u t i o n a r y identification of  Computer Animation and Simulation Workshop,  Eurographics, pages 245-252, M a a s t r i c h t , Netherlands, Sept 1995. [78] H . Maafi and U . K i i h n a p f e l . of l i v i n g tissue. In  Noninvasive measurement of elastic properties  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 . W l o d a r c z y k . s c u l p t i n g system w i t h haptic toolkits. In  Graphics,  V i r t u a l clay: A real-time  Proc. Symposium on Interactive 3D  pages 179-190. A C M , M a r 2001.  [80] T . M e r z , D . P a u l u s , and H . N i e m a n n . L i n e 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 . M e t a x a s .  Physics-Based Deformable Models - Applications to Computer  Vision. Graphics and Medical Imaging.  Kluwer Academic Publisher, Boston,  1996. [82] D . M e t a x a s and D . Terzopoulos. t h r o u g h physics-based synthesis.  Shape and nonrigid m o t i o n estimation  IEEE Transactions on Pattern Recognition  and Machine Intelligence,  15(6):580-591, 1993.  [83] S . M . M e t w a l l i , A . R . R a g a b , A . H . K a r n e l , and A . A . Saheb. D e t e r m i n a t i o n of plastic stress-strain behavior by digital-image-processing technique.  the Society for Experimental Mechanics, Inc.,  135  Proc. -of  44:414-422, 1987.  [84]  G . N a k a m u r a and  K . Tanuma.  A nonuniqueness theorem for  b o u n d a r y value problem i n elasticity. 56(2):602-610,  an  inverse  SIAM Journal on Applied Mathematics,  1996.  [85] G . N a k a m u r a and G . U h l m a n n . Identification of L a m e parameters by b o u n d ary  measurements.  American Journal of Mathematics,  115(5):1161-1187,  1993. [86]  G . N a k a m u r a and G . U h l m a n n . G l o b a l uniqueness for an inverse b o u n d a r y p r o b l e m arising i n elasticity.  [87]  Inventiones Mathematicae,  118(3):457-474,  1994.  G . N a k a m u r a and G . U h l m a n n . Inverse problems at the b o u n d a r y for an elastic medium.  SIAM Journal on Mathematical Analysis,  [88] J . - C . N e b e l . Soft tissue m o d e l l i n g from 3 D scanned data. In  26(2):263-279,  Deform,  1995.  Geneva,  S w i t z e r l a n d , N o v . 2000. [89]  J . - C . N e b e l , F . J . R o d r i g u e z - M i g u e l , and W . P. Cockshott. Stroboscobic stereo  International Conference on 3-D Digital Imaging and  rangefinder.  In 3rd  Modeling,  pages 59-64, Quebec C i t y , C a n a d a , 2001.  IEEE.  [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 . H o d g i n s . G r a p h i c a l modeling and a n i m a t i o n of brittle 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 . P a i , editor.  Reality-based modeling and applications in reverse engineer-  ing, computer graphics and VR, W o r k s h o p . 136  I E E E , 2000.  [93]  D . K . P a i , J . L a n g , J . E . L l o y d , a n d . J . L . R i c h m o n d . Reality-based m o d e l i n g w i t h A C M E : A progress report. In D . R u s and  S. Singh, editors,  7th Inter-  national Symposium on Experimental Robotics, LNCIS 271,  pages 121-130,  H o n o l u l u , U S A , December 2000. Springer-Verlag. [94]  D . K . P a i , J . L a n g , J . E . L l o y d , and R . J . W o o d h a m . A C M E , a t e l e r o b o t i c active measurement facility. In P. C o r k e and J . Trevelyan, editors,  6th International  Symposium on Experimental Robotics, LNCIS 250, A u s t r a l i a , M a r c h 1999. [95]  pages 391-400, Sydney,  Springer-Verlag.  D . K . P a i , K . van den D o e l , D . L . James, J . L a n g , J . E . L l o y d , J . L . R i c h m o n d , and S . H . Y a u . Scanning physical interaction behavior of 3D objects. In  puter Graphics, Annual Conference Series,  Com-  pages 87-96, Los Angles, A u g  2001.  A C M SIGGRAPH. [96]  X . Papademetris.  Estimation of 3D left ventricular deformation from medical  images using biomechanical models. [97]  P h D thesis, Yale University, M a y 2000.  H . W . P a r k , S. S h i n , and H . S . Lee. D e t e r m i n a t i o n of an o p t i m a l regularization factor in system identification w i t h T i k h o n o v regularization for linear elastic continua. 51:1211-1230,  [98]  International Journal for Numerical Methods in Engineering, 2001.  A . P e n t l a n d . A u t o m a t i c extraction of deformable part models.  Journal of Computer Vision, [99]  4(2):107-126,  International  1990.  G . P i c i n b o n o , H . Delingette, and N . Ayache. Non-linear and anisotropic elastic soft tissue models for  medical simulation.  Robotics and Automation,  In  International Conference on  pages 1371-1376, Seoul, S o u t h K o r e a , M a y  137  2001.  [100] S . M . P i a t t and N . I . B a d l e r . A n i m a t i n g facial expression. In  ics, Annual Conference Series,  Computer Graph-  pages 245-252, Dallas, U S A , A u g 1981. A C M  SIGGRAPH. [101] E . P . P o p o v .  Engineering mechanics of solids.  Prentice H a l l , E n g l e w o o d Cliffs,  N J , U S A , 1st edition, 1990. [102] W . H . Press, S . A . Teukolsky, W . T . Vetterling, and B . P . Flannery.  recipes in C: the art of scientific computing.  Numerical  C a m b r i d g e U n i v e r s i t y Press, 2  edition, 1992. [103] K . R . R a g h a v a n and A . E . Yagle. Inverse and o p t i m a l drive problems i n elasticity imaging of soft tissue. I n  Medical Imaging Conference,  pages 1886-1890.  I E E E , 1993. [104] R . R a m a n a t h a n and D . M e t a x a s . D y n a m i c deformable models for enhanced haptic rendering i n v i r t u a l environments.  In  Virtual Reality,  pages 31-35.  I E E E , 2000. [105] P o i n t G r e y Research. Triclops on-line m a n u a l , h t t p : / / w w w . p t g r e y . c o m / . [106] D . R e z n i k and C . Laugier. D y n a m i c simulation and v i r t u a l control of a deformable fingertip. I n  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,  I E E E , 1994.  138  pages 391-394.  [108] C . R o c c h i n i , P . C i g n o n i , C . M o n t a n i , and P . Scopigno. M u l t i p l e texture stitching and blending on 3d objects. In 10th  Workshop on Rendering,  Eurographics,  pages 173-180, G r a n a d a , S p a i n , 1999. [109] G . R o t h and E . W i b o w o o . A n efficient v o l u m e t r i c m e t h o d for b u i l d i n g closed triangular meshes from 3-d image and point data. I n  Proc. Graphics Interface,  pages 173-180, 1997. [110] P . J . Rousseeuw. Least median of squares regression.  Statistical Association,  Journal of the American  79(388):871-880, 1984.  [ I l l ] P . J . Rousseeuw and K . V a n Driessen.  C o m p u t i n g L T S regression for large  d a t a sets. Technical report, U n i v e r s i t y of A n t w e r p , 1999. [112] M . A . Sagar, D . B u l l i v a n t , G . D . M a l l i n s o n , P . J . H u n t e r , and I . W . Hunter. A v i r t u a l environment and m o d e l of the eye for surgical s i m u l a t i o n . In  Graphics, Annual Conference Series,  Computer  pages 205-212, O r l a n d o , U S A , J u l 1994.  A C M SIGGRAPH. [113] J . K . S a l i s b u r y and M . A . Srinivasan. P h a n t o m - b a s e d h a p t i c interaction w i t h v i r t u a l objects.  IEEE Computer Graphics & Applications,  [114] M . Sanayei and M . J . Saletnik;  17(5):6-10, 1997.  Parameter estimation of structures form  static strain measurements. I. F o r m u l a t i o n .  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. E r r o r sensitivity analysis.  gineering,  122(5):555-562, 1996.  139  Journal of Structural En-  , [116] Y . Sato, M . Wheeler, and K . Ikeuchi. O b j e c t shape and reflectance m o d e l i n g from observation.  In  Computer Graphics. Annual Conference Series,  pages.  379-387, L o s Angles, U S A , A u g 1997. A C M . ; [117] D . S . Schnur and N . Zabaras. . .  A n inverse m e t h o d for determining elastic-  material properties and a material interface.  merical Methods in Engineering,  International Journal for Nu-  33(10):2039-2057, 1992.  Dynamics of Multibody Systems.  [118] A . A . Shabana.  J o h n W i l e y & Sons, N e w  Y o r k , U S A , 1989. [119] H . Spies, B . Jahne, and J . L . B a r r o n . Dense range flow from depth and intensity data. I n  International Conference on Pattern Recognition,  volume 1, pages  131-134, B a r c e l o n a , Spain, 2000. [120] H . Spies, B . Jahne, and J . L . B a r r o n . Regularised range flow. I n  ference on Computer Vision,  European Con-  L C N S 1843/44, pages 785-799, D u b l i n , Ireland,  2000. I E E E , Springer-Verlag, B e r l i n . [121] G . E . Stavroulakis and H . Antes.  Nondestructive elastostatic identification  of unilateral cracks through b e m and neural networks.  Mechanics,  Computational-  20(5):439-451, 1997.  [122] G . Szekely, C . Brechbiihler, J . D u a l , R . E n z l e r , J . H u g , R . H u t t e r , N . Ironmonger, M . K a u e r , V . Meier, P. Niederer, A . R h o m b e r g , P . S c h m i d , G . Schweitzer, M . T h a l e r , V . Vuskovic, G . Troster, U . Haller, and M . B a j k a . V i r t u a l realitybased s i m u l a t i o n of endoscopic surgery.  140  Presence,  9(3):310-333, 2000.  [123] F . Tendick, M . Downes, T . G o k t e k i n , M . C . Cavusoglu, D . Feygin, X . W u , • R . E y a l , M . Hegarty, and L . W . Way. A v i r t u a l environment testbed for t r a i n ing laparoscopic surgical skills.  Presence,  9(3):236-255, 2000.  [124] D . Terzopoulos. R e g u l a r i z a t i o n of inverse v i s u a l problems involving discontinuities.  IEEE Transactions on Pattern Recognition and Machine Intelligence,  8(4):413-424, 1986. [125]- D . Terzopoulos and K . Fleischer. M o d e l i n g inelastic deformation: viscoelasticity, plasticity, fracture. I n  Computer Graphics, Annual Conference Series,  pages 269-278, A t l a n t a , U S A , A u g 1988. A C M S I G G R A P H . [126] D . Terzopoulos and K . Waters. P h y s i c a l l y - b a s e d facial modeling, analysis a n d animation.  Journal of Visualization and Computer Animation,  1:73-80, 1990.  [127] D . Terzopoulos, A . P . W i t k i n , and M . K a a s . Constraints on deformable models: Recovering 3 D shape and nonrigid m o t i o n .  Artificial Intelligence,  36(1):91-  123, 1988. [128] S. V e d u l a , S. Baker. P. Rander, R . C o l l i n s , and T . K a n a d e . Three-dimensional -.'. scene flow. I n  International Conference on Computer Vision,  pages 722-729,  1999. [129] V . Vuskovic, M . K a u e r , G . Szekely, and M . Reidy.  Realistic force feedback  for v i r t u a l reality based diagnostic surgery simulators. In  ference on Robotics and Automation, A p r i l 2000.  141  International Con-  pages 1592-1598, San Francisco, U S A ,  [130] R . J . W o o d h a m , E . C a t a n z a r i t i , and A . K . M a c k w o r t h .  A n a l y s i s by syntheis  ) . i n c o m p u t a t i o n a l vision w i t h a p p l i c a t i o n to remote sensing.  Intelligence, l(2):71-79,  Computational  1985.  [131] X . W u , M . S . Downes, T . G o t e k i n , and F . Tendick. A d a p t i v e nonlinear finite elements for deformable b o d y s i m u l a t i o n using d y n a m i c progressive meshes. In A . C h a l m e r s and T . - M . R h y n e , editors,  Eurographics,  Computer Graphics  F o r u m , 20(3), pages C 3 4 9 - C 3 5 8 , Manchester, U k , 2001. [132] M . Y a m a m o t o , P . Boulanger, J . - A . B e r a l d i n , and M . R i o u x . D i r e c t 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 . Z h a n g and C . K a m b h a m e t t u . Integrated 3d scene flow and structure recovery from m u l t i v i e w image sequences.  Recognition,  In  Computer Vision and Pattern  volume 2, pages 674-681, 2000.  [134] Z . Z h a n g . Parameter estimation techniques: A t u t o r i a l w i t h application to conic fitting. [135] Z . - H . Zhong.  Image and Vision Computing,  15(l):59-76, 1997.  Finite Element Procedures for Contact-Impact Problems.  Oxford  U n i v e r s i t y Press, Oxford, U k , 1993. [136] Y . Z h u a n g and J . Canny. H a p t i c interactions w i t h global deformations.  International Conference on Robotics and Automation, Francisco, U S A , A p r i l 2000.  142  In pages 2428-2433, San  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0051715/manifest

Comment

Related Items