USFDC Home  USF Electronic Theses and Dissertations   RSS 
Material Information
Subjects
Notes
Record Information

Full Text 
xml version 1.0 encoding UTF8 standalone no
record xmlns http:www.loc.govMARC21slim xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.loc.govstandardsmarcxmlschemaMARC21slim.xsd leader nam Ka controlfield tag 001 001670319 003 fts 005 20051216093243.0 006 med 007 cr mnuuuuuu 008 051109s2005 flu sbm s000 0 eng d datafield ind1 8 ind2 024 subfield code a E14SFE0001177 035 (OCoLC)62248641 SFE0001177 040 FHM c FHM 049 FHMM 090 TK7885 (Online) 1 100 Zhang, Yong. 0 245 Robust algorithms for property recovery in motion modeling, medical imaging and biometrics h [electronic resource] / by Yong Zhang. 260 [Tampa, Fla.] : b University of South Florida, 2005. 500 Includes vita. 502 Thesis (Ph.D.)University of South Florida, 2005. 504 Includes bibliographical references. 516 Text (Electronic thesis) in PDF format. 538 System requirements: World Wide Web browser and PDF reader. Mode of access: World Wide Web. Title from PDF of title page. Document formatted into pages; contains 134 pages. 520 ABSTRACT: The past two decades has witnessed growing interest in physics based techniques in computer vision, computer graphics and medical imaging. The main advantage of a physical model is its mathematical rigor and physical soundness, which makes it an ideal tool to study complex nonrigid motion. However, since a model based on continuum mechanics is computationally demanding, an idealized framework is often adopted where physical motion parameters are significantly simplified, which inevitably affects the accuracy and reliability of modeling results. In this study, a new modeling approach is developed that features the reconstruction of actual material properties such as the Young's modulus and the Poisson's ratio. Justified by the constitutive law and mathematical considerations, the Young's modulus is identified as a unique physical motion parameter.By imposing an adaptive smoothness constraint, the Young's modulus helps preserve the local characteristics (discontinuity) of an object's deformation, a role similar to the weighting coefficient in the study of edgepreserving visual surface reconstruction. The contribution of this work is fourfold: (1) two recovery algorithms are developed to solve the inverse elastic problem: A deterministic algorithm that is based on the GaussNewton method and the general cross validation, and a stochastic algorithm that is based on the constrained genetic evolution; (2) a new modeling approach is proposed that has the ability to recover nonrigid motion in terms of the physical parameters. The use of recovered parameters can be implemented within a boundarydriven motion synthesis scheme; (3) A sensitivity method is proposed to evaluate the impact of different parameters. 590 Adviser: Dr. Dmitry B. Goldgof. Coadviser: Dr. Sudeep Sarkar 653 Nonrigid deformation. Elasticity. Computer vision. Face recognition. 690 Dissertations, Academic z USF x Computer Science and Engineering Doctoral. 773 t USF Electronic Theses and Dissertations. 4 856 u http://digital.lib.usf.edu/?e14.1177 PAGE 1 Robust Algorithms for Prop ert y Reco v ery in Motion Mo deling, Medical Imaging and Biometrics b y Y ong Zhang A dissertation submitted in partial fulllmen t of the requiremen ts for the degree of Do ctor of Philosoph y in Computer Science and Engineering Departmen t of Computer Science and Engineering College of Engineering Univ ersit y of South Florida CoMa jor Professor: Dmitry B. Goldgof, Ph.D. CoMa jor Professor: Sudeep Sark ar, Ph.D. Don R. Hilb elink, Ph.D. An thon y J. Llew ellyn, Ph.D. Leonid V. Tsap, Ph.D. Date of Appro v al: Ma y 3, 2005 Keyw ords: Nonrigid Deformation, Elasticit y Computer Vision, F ace Recognition. c r Cop yrigh t 2005, Y ong Zhang PAGE 2 DEDICA TION T o m y paren ts PAGE 3 A CKNO WLEDGMENTS I w ould lik e to thank m y t w o ma jor advisors, Dr. Goldgof and Dr. Sark ar, for their academic guidance and y ears of supp ort. I w as strongly inruenced b y their \vision" of what Computer Vision should b e as a relativ ely new discipline. The one thing that I appreciate the most is that they ga v e me great freedom to pursue m y researc h in terest and pro vide opp ortunities so that I can w ork on a wide range of topics. I couldn't ha v e c hosen b etter sup ervisors. I am v ery fortunate to ha v e the opp ortunit y to w ork with Dr. Hall. Without his insigh t and exp ertise, the t w o pap ers on the constrained genetic approac h w ould not b e p ossible. I am grateful to Dr. Hilb elink, Dr. Llew ellyn and Dr. Tsap for taking the time from their busy sc hedule to b e the committee mem b er and pro viding v ery helpful commen ts on the draft. I also thank the Departmen t of Computer Science and Engineering at the Univ ersit y of South Florida for pro viding the lab facilit y that allo w ed me to complete m y dissertation. Most imp ortan t of all, I wish to thank m y family for their unfailing supp ort and contin ually helping me in ev ery w a y p ossible. PAGE 4 T ABLE OF CONTENTS LIST OF T ABLES iii LIST OF FIGURES iv ABSTRA CT v CHAPTER 1 INTR ODUCTION 1 1.1 Ov erview 1 1.2 F undamen tal Issues and Related W ork 2 1.2.1 Prop ert y Reco v ery and In v erse Problem 2 1.2.2 Computational Issues and Algorithms 3 1.2.3 Nonrigid Motion Mo deling 6 1.3 Con tributions 9 CHAPTER 2 RECO VER Y ALGORITHMS 10 2.1 F orw ard Mo del: Linear Elastic Deformation 10 2.2 Reco v ery Algorithm: Deterministic Approac h 12 2.2.1 Tikhono v Regularization 12 2.2.2 Baseline Solution Strategy 13 2.2.3 GCV Solution Strategy 17 2.2.4 Uniqueness of Solution 19 2.3 Reco v ery Algorithm: Constrained Genetic Approac h 22 2.3.1 Genetic Enco ding 22 2.3.2 Rank Based Constrain t 24 2.3.3 Balancing Fitness and P enalt y b y Sto c hastic Ranking 28 2.3.4 Genetic Op erators 29 2.4 Deterministic or Sto c hastic Approac h 32 2.4.1 Exp erimen ts with GA, CGA and GNM 32 2.4.2 Exp erimen ts with CGA and GNM 37 CHAPTER 3 APPLICA TION IN NONRIGID MOTION MODELING 45 3.1 Motion P arameters in Ph ysical Mo del 45 3.1.1 Ph ysical Motion P arameters 45 3.1.2 Ph ysical and Mathematical Justications 46 3.2 Syn thesis with Reco v ered P arameter 52 3.2.1 Three P ersp ectiv es 52 3.2.2 Diric hletT yp e Motion Syn thesis with Reco v ered Prop ert y 55 3.3 Exp erimen ts 58 i PAGE 5 3.3.1 Reco v ery of Ph ysical Motion P arameter 58 3.3.2 Motion Syn thesis with Reco v ered P arameter 63 CHAPTER 4 APPLICA TION IN BURN SCAR ASSESSMENT 68 4.1 Mo del Construction 69 4.1.1 Extraction of Natural P oin t F eatures 69 4.1.2 Adaptiv e Meshing 71 4.2 Estmate Scar Elasticit y 72 4.2.1 Range Scanner Setting 72 4.2.2 Boundary Condition Sp ecication 73 4.3 Exp erimen tal Results 74 4.3.1 Data Set 74 4.3.2 Quan tify Scar Elasticit y 74 4.3.3 Relativ e Elastic Index 75 CHAPTER 5 APPLICA TION IN F A CE RECOGNITION 78 5.1 Hyp othesis: F ace Anatom y and Muscle Biomec hanics 79 5.2 Computational Metho d 81 5.3 Exp erimen tal Results 84 CHAPTER 6 SENSITIVITY ANAL YSIS 87 6.1 Computational Metho d 88 6.1.1 Primary Problem 89 6.1.2 P erformance Measure and Dimensional Sensitivit y 90 6.1.3 Adjoin t State Metho d 92 6.1.4 Deriv ativ e Computation 93 6.1.5 Normalized Sensitivit y 94 6.2 Exp erimen ts with Syn thetic Mo del 95 6.2.1 Mo del Conguration and P erformance Measure 95 6.2.2 Normalized Sensitivit y Distribution 97 6.2.3 Dimensional Sensitivit y Distribution 99 6.2.4 Relativ e Error Distribution 99 6.3 Exp erimen ts with Burn Scar 101 6.3.1 Mo del Conguration 101 6.3.2 Normalized Sensitivit y 102 6.3.3 Dimensional Sensitivit y 102 6.4 Exp erimen ts with Deformed F ace Mo del 104 6.4.1 F ace Mo del Conguration 104 6.4.2 Normalized Sensitivit y 104 6.4.3 Dimensional Sensitivit y 105 CHAPTER 7 CONCLUSIONS AND DISCUSSIONS 109 REFERENCES 113 ABOUT THE A UTHOR End P age ii PAGE 6 LIST OF T ABLES T able 1. General Ob ject Reconstruction Diagram. 1 T able 2. In v ersion Algorithms for Elastic Prop ert y Reconstruction. 5 T able 3. Link T able for Main taining Original Spatial Connections. 24 T able 4. CGA Exp erimen ts with Dieren t P opulation Sizes. 39 T able 5. P erformance of CGA and GNM with Dieren t Initial Conditions. 40 T able 6. Adaptiv e Smo othing in Motion Analysis and Surface Reconstruction. 51 T able 7. Three P ersp ectiv es for Motion Syn thesis. 55 T able 8. Syn thesized Motion along CrossSection Line. 65 T able 9. The P erformance of P oin t F eature Detector. 71 T able 10. The Estimated Relativ e Scar Elasticit y (kP a). 74 T able 11. Ph ysician's Rating and Relativ e Elastic Index (REI). 77 T able 12. F ace Exp erimen ts. 84 iii PAGE 7 LIST OF FIGURES Figure 1. Genetic Co ding of the Y oung's Mo dulus in a Finite Elemen t Mo del. 23 Figure 2. The Rank Based Sc heme for Incorp orating Qualitativ e Prior Information. 27 Figure 3. Syn thetic Mo del and the Qualitativ e Prior Kno wledge. 33 Figure 4. Reconstructed Y oung's Mo dulus b y GA, CGA and GNM. 34 Figure 5. Reconstructed Y oung's Mo dulus b y the CGA with Biased Prior. 36 Figure 6. Conguration of the F orw ard Mo del. 37 Figure 7. Con v ergence Beha vior of CGA with Dieren t P opulation Sizes. 39 Figure 8. Con v ergence P athes of CGA and GNM with Dieren t Starting P oin ts. 44 Figure 9. Results of GNM with V arious Initial Conditions. 44 Figure 10. Congurations of Spring Before and After Stretc hing. 52 Figure 11. Adaptiv e Con trol of Motion Field b y the Y oung's Mo dulus. 53 Figure 12. Conguration of a Syn thetic Spring. 58 Figure 13. Tw o F orw ard Sim ulations with Dieren t but Prop ortional Elasticities. 59 Figure 14. In v erse Solutions with Dieren t Regularization P arameters. 61 Figure 15. Norms against Regularization P arameter and SemiCon v ergence. 62 Figure 16. P erformances of Reco v ery Algorithms. 62 Figure 17. Image Sequence of a Bandage for Motion Syn thesis. 64 Figure 18. Finite Elemen t Mo del and the Reco v ered Elasticit y 64 Figure 19. Comparison of Tw o Exp erimen ts. 67 Figure 20. P osition Shift Bet w een the Corresp onding F eatures. 70 Figure 21. Adaptiv e Meshing Using P oin t F eatures. 72 iv PAGE 8 Figure 22. The Geometry of K2T Range Scanner for Burn Scar Study 73 Figure 23. Estimated Elasticit y and Strain. 76 Figure 24. Correlation Bet w een Ph ysician's Rating and REIs. 77 Figure 25. Strain P attern of a Mo died F ace. 83 Figure 26. F our Images Acquired for Eac h Sub ject. 83 Figure 27. R OCs for Regular and Mo died F aces. 85 Figure 28. Syn thetic Mo del Conguration. 96 Figure 29. Deformed Mo dels. 97 Figure 30. Normalized Sensitivit y Con tour Maps. 98 Figure 31. Dimensional Sensitivit y Con tour Maps. 98 Figure 32. Relativ e Error Distributions. 100 Figure 33. Deformation of Burn Scar. 100 Figure 34. Normalized Sensitivit y Maps of Scar. 101 Figure 35. Dimensional Sensitivit y Maps of Scar. 103 Figure 36. F ace Deformation with Expression. 103 Figure 37. Normalized Sensitivit y Maps of F ace Mo del. 105 Figure 38. Dimensional Sensitivit y Maps of F ace Mo del. 107 v PAGE 9 R OBUST ALGORITHMS F OR PR OPER TY RECO VER Y IN MOTION MODELING, MEDICAL IMA GING AND BIOMETRICS Y ong Zhang ABSTRA CT The past t w o decades has witnessed gro wing in terest in ph ysics based tec hniques in computer vision, computer graphics and medical imaging. The main adv an tage of a ph ysical mo del is its mathematical rigor and ph ysical soundness, whic h mak es it an ideal to ol to study complex nonrigid motion. Ho w ev er, since a mo del based on con tin uum mec hanics is computationally demanding, an idealized framew ork is often adopted where ph ysical motion parameters are signican tly simplied, whic h inevitably aects the accuracy and reliabilit y of mo deling results. In this study a new mo deling approac h is dev elop ed that features the reconstruction of actual material prop erties suc h as the Y oung's mo dulus and the P oisson's ratio. Justied b y the constitutiv e la w and mathematical considerations, the Y oung's mo dulus is iden tied as a unique ph ysical motion parameter. By imp osing an adaptiv e smo othness constrain t, the Y oung's mo dulus helps preserv e the lo cal c haracteristics (discon tin uit y) of an ob ject's deformation, a role similar to the w eigh ting co ecien t in the study of edgepreserving visual surface reconstruction. The con tribution of this w ork is fourfold: (1) t w o reco v ery algorithms are dev elop ed to solv e the in v erse elastic problem: A deterministic algorithm that is based on the GaussNewton metho d and the general cross v alidation, and a sto c hastic algorithm that is based on the constrained genetic ev olution; (2) a new mo deling approac h is prop osed that has the abilit y to reco v er nonrigid motion in terms of the ph ysical parameters. The use of reco v ered parameters can b e implemen ted within a b oundarydriv en motion syn thesis sc heme; (3) A vi PAGE 10 sensitivit y metho d is prop osed to ev aluate the impact of dieren t parameters. The metho d uses the adjoin t state equation and hence is suitable for large scale mo dels. (4) the prop osed mo deling approac h has b een applied to burn scar assessmen t and face recognition. vii PAGE 11 CHAPTER 1 INTR ODUCTION 1.1 Ov erview One of the fundamen tal issues in computer vision and image understanding is to extract basic attributes of an ob ject from 2D images. Those attributes can b e either directly utilized or further pro cessed to facilitate higher lev el tasks suc h as classication, recognition and syn thesis etc. The commonly studied attributes and asso ciated reconstruction metho ds are listed in T able 1. Muc h progress has b een made in reconstructing 3D shap e and optical prop erties (surface rerectance) under the common theme of \shap e from X". Ho w ev er, little has b een done to reco v er the material prop erties that are hidden b eneath the visible surface suc h as the elasticit y and the viscosit y It is in teresting to note that, ev en in the famous Marr's diagram, only the geometrical primitiv es (surface and shap e) are emphasized. On the other hand, h uman vision can easily detect the subtleties in the w a y an ob ject mo v es and deforms, and further relate it to the underlying mec hanical prop erties. F or example, simply lo oking at the expressions of a bab y's face and a wrinkled face can lead us to b eliev e (no matter ho w implicitly suc h a p erception is formed) that the bab y's skin is more elastic, whic h is probably T able 1. General Ob ject Reconstruction Diagram. A ttributes Examples Reconstruction Metho ds geometric 3D shap e/structure shap e from X surface rerectance and color radiometric mo dels kinematic deformation matrix ane motion mo dels dynamic force and b oundary conditions ph ysical motion mo dels mec hanical elasticit y and viscosit y prop ert y estimation 1 PAGE 12 one of the reasons that w e nd the bab y's face more app ealing. As another example, b y observing the w a y clouds mo v e, sea w a v es propagate and a rag raps in the wind, a p erson can dev elop a sense that the ob jects ha v e dieren t material prop erties, ev en though he or she ma y not b e able to giv e a quan titativ e description using the language of con tin uum mec hanics. One of the ob jectiv es of this researc h is to dev elop a metho d that can reco v er material prop erties from images, in particular, the Y oung's mo dulus of elastic ob jects. In addition to the aforemen tioned theoretical in terest in material p erception, reconstructing elastic prop ert y is also of in terest to researc hers in man y practical elds. F or example, iden tifying elasticit y abnormalities in soft tissue has high p oten tial in early cancer detection [37 94 63 98 34 ], b ecause pathological pro cesses often cause prop ert y c hanges in the diseased tissues that amplify themselv es as hard inclusions. Kno wing elastic prop erties of deformable ob jects is also critical for ph ysicsbased mo deling suc h as motion trac king, realistic animation, scar assessmen t, visualization and surgical planning [74 127 135 82 99 35], b ecause the qualit y of elastic prop erties could ha v e a strong impact on the p erformance of ph ysical mo dels. Reco v ering elastic prop erties has also found applications in structural damage iden tication, geoph ysical exploration, w afer engineering, rob otics design and comp osite material c haracterization [22 7, 38 67 45 ]. In this study applications of the elastic prop ert y reco v ery will b e demonstrated in three areas: (1) nonrigid motion mo deling; (2) quan titativ e burn scar assessmen t; (3) face recognition (biometrics). 1.2 F undamen tal Issues and Related W ork 1.2.1 Prop ert y Reco v ery and In v erse Problem Lik e man y vision problems, inferring the elastic prop erties from the observ ed deformation is an illp osed in v erse problem. A problem can b e regarded as either a forw ard problem or an in v erse problem dep ending on its causeeect alignmen t. A problem that describ es the eect based on the giv en cause is a forw ard problem, while a problem that is aimed at understanding the cause of a consequence is an in v erse problem. In man y disciplines, the 2 PAGE 13 forw ard problem is to nd out the b eha vior of a system resp onding to the external excitation (cause). F or example, w e can mo del the nonrigid deformation of an elastic ob ject caused b y the surface force. On the other hand, the in v erse problem is to infer the quan tities that directly or indirectly caused the eect, suc h as the force and the elastic prop erties. A ph ysical la w can sp ecify the cause and the corresp onding eect with mathematical equations. It has long b een found that in v erse problems are c haracterized b y their unstable solutions, while the forw ard problems are usually w ell b eha v ed n umerically By pro viding a formal mathematical denition of the w ellp osedness and the illp osedness, Hadamard [53 ] laid the foundation of in v erse problem theory and pa v ed the road of nding a meaningful solution of man y in v erse problems. In the Hadamard sense, a w ellp osed problem m ust ha v e a solution that satises the follo wing three conditions: existence, uniqueness and stabilit y (con tin uous dep endence of the solution on the input data). A problem is illp osed if one of those conditions failed. Most of the forw ard problems are w ellp osed and the corresp onding in v erse problems are illp osed. The dicult y in solving an illp osed in v erse problem is often related to the nonuniqueness and the instabilit y of the solution. In practice, an in v erse problem is often transformed in to an optimization problem. Man y algorithms ha v e b een dev elop ed to reduce the solution space and to smo oth the solution path in order to restore the uniqueness and the stabilit y of the n umerical solution. The basic idea is to utilize a priori kno wledge ab out the ph ysical system. F or example, deterministic regularization is a commonly used metho d for parameter estimation [30 ]. Statistical metho ds in the Ba y esian framew ork ha v e also b een in v estigated [33 124 ]. In the con text of computer vision, P oggio et al. [106 ] presen t a comprehensiv e review of the in v erse problem theory and its implication to image understanding and visual p erception. 1.2.2 Computational Issues and Algorithms The regularized functional of a nonlinear illp osed in v erse problem is dicult to solv e, b ecause its con v exit y can only b e guaran teed v ery lo cally In the deterministic domain, the use of gradien tbased metho ds is common [30 ]. The gradien tbased metho ds are sensitiv e to 3 PAGE 14 the starting p oin t and ma y stuc k in the lo cal extrema. Therefore, solving in v erse problems with the sto c hastic algorithms has receiv ed m uc h atten tion recen tly The h ybrid algorithm that utilizes the strengths of b oth the deterministic and the sto c hastic metho ds has also b een prop osed [22 ]. The related w ork is summerized in T able 2. The list only includes the pap ers that deal with the computational asp ect of solving an in v erse elastic problem. F or studies on the direct measuremen t and the dev elopmen t of sp ecial imaging mo dalities suc h as the magnetic resonance elastograph y the ultrasonic elastograph y and the optical coherence tomograph y readers are refered to [42 31 34 95 71 ]. It is apparen t that most of the studies c hose the deterministic algorithms that use either direct or iterativ e in v ersions. Kallel and Bertrand [68 ] used the NewtonRaphson metho d com bined with a nite elemen t mo del to estimate the Y oung's mo dulus of syn thetic tissues. Do yley et al [27 ] studied a mo died iterativ e NewtonRaphson metho d. V an Houten et al [139 ] prop osed a m ultiresolution metho d (subzone) to impro v e b oth the eciency and robustness. Other deterministic algorithms suc h as the lev el set metho d [6 ], the adjoin t state metho d [90 ], the steep est descen t searc h [135 ] and the iterativ e regularization [29 ] ha v e also b een rep orted. Due to the lac k of quan titativ e comparison data for the t w o approac hes (deterministic vs. sto c hastic), it is often a dicult task to design an optimal solution strategy that is b oth robust and ecien t. In this study b oth the deterministc and the sto c hastic approac hes will b e in v estigated. Sp ecically a NewtonRaphson metho d based on the general cross v alidation and a constrained genetic algorithm will b e discussed. Sev eral comparativ e exp erimen ts w ere also p erformed to examine the robustness and the con v ergence rate of the t w o algorithms. Both linear and nonlinear constitutiv e mo dels ha v e b een studied. The linear mo del is computationally attractiv e and usually a go o d appro ximation for man y ob jects. In this study a linear elastic mo del will b e used. 4 PAGE 15 T able 2. In v ersion Algorithms for Elastic Prop ert y Reconstruction. Authors Date Algorithms Mo dels Exp erimen ts (Materials) Agly amo v et al [3] 2004 Deterministic Linear Throm b osis Ameur et al [6] 2004 Deterministic Linear Numerical Arns et al [7 ] 2002 Deterministic Linear P orous media Barb one et al [12 ] 2004 Deterministic Linear Numerical Chiroiu et al [21 ] 2000 Genetic Nonlinear Crystal lattice Chiwiaco wsky et al [22 ] 2004 Genetic Linear Numerical Do yley et al [27 ] 2000 Deterministic Linear Numerical Engl et al [29 ] 2003 Deterministic Linear P olymer Eskin et al [32 ] 2002 Deterministic Linear Numerical F ranca et al [38 ] 2004 Deterministic Linear Silicon w afer F uruk a w a et al [43 ] 1997 Genetic Nonlinear CrMo steel Garb o czi et al [45 ] 1995 Deterministic Linear Numerical Hori et al [60 ] 2003 Deterministic Linear Numerical Huang et al [61 ] 2001 Deterministic Nonlinear Numerical Ji et al [65 ] 2004 Deterministic Linear Numerical John et al [66 ] 2003 Genetic Linear P elvis b one Joukhadar et al [67 ] 1997 Genetic Linear Numerical Kallel et al [68 ] 1996 Deterministic Linear Numerical Kirkpatric k et al [71 ] 2002 Deterministic Linear Tissues Maurice et al [79 ] 2004 Deterministic Linear V essel w alls Miga [83 ] 2003 Deterministic Linear Numerical Muth upillai et al [87 ] 1995 Deterministic Linear Gel Nak am ura et al [88 ] 1993 Deterministic Linear Numerical Ob erai et al [90 ] 2003 Deterministic Linear Numerical O'Donnell et al [91 ] 2004 Deterministic Nonlinear Heart tissue Oliphan t et al [92 ] 2001 Deterministic Linear Phan tom P ellotBarak at et al [100 ] 2004 Deterministic Linear Phan tom Plew es et al [105 ] 2000 Deterministic Linear Numerical Ragha v an et al [108 ] 1994 Deterministic Linear Numerical Samani et al [113 ] 2001 Deterministic Linear Breast tissues Sumi et al [122 ] 1998 Deterministic Linear Numerical Tsap et al [135 ] 1998 Deterministic Nonlinear Skin V an Houten et al [139 ] 1999 Deterministic Linear Numerical Zh u et al [149 ] 2003 Deterministic Linear Numerical Pap ers that use the hybrid appr o ach c ombining the deterministic and genetic algorithms ar e classie d as Genetic. Pap ers that c onsider e d ge ometric al or material line arities to simplify the forwar d mo del ar e classie d as Line ar. No mor e than two p ap ers ar e sele cte d fr om the same r ese ar ch gr oup (authors) unless the algorithms use d ar e signic antly dier ent. 5 PAGE 16 1.2.3 Nonrigid Motion Mo deling Muc h of this study is conducted in the con text of nonrigid motion mo deling. Nonrigid motion is more dicult y to analyze than its rigid coun terpart, b ecause it deals with a m uc h wider range of ob jects and more complex deformation patterns. It is unlik ely that a single solution sc heme can b e found that is applicable to all t yp es of nonrigid motion. Therefore, the follo wing strategy has b een adopted: (1) classify nonrigid motion in to dieren t t yp es and deal with eac h t yp e individually; (2) use a motion mo del to facilitate the analysis b ecause a go o d mo del can signican tly reduce the reco v ery and trac king complexities. Huang [62 ] classied nonrigid motion in to three basic t yp es: Articulated motion, elastic motion and ruid motion. Kam bhamettu et al. [69 ] further extended the classication in to more detailed subgroups based on the geometrical c haracteristics of deformation. Aggarw al et al [2 ] surv ey ed the existing nonrigid motion mo deling tec hniques, particularly those based on geometrical and ph ysical mo dels. A k ey issue in the mo delbased approac h is to dene and reco v er motion parameters. Motion parameters pla y an imp ortan t role of carrying v arious constrain ts that w e wish to imp ose on the motion. The go v erning equations of a motion mo del are c haracterized b y lo cal descriptions of the degrees of freedom of deforming ob jects. Theoretically this lo calit y prop ert y enables us to ac hiev e a mo deling accuracy close to an analytical solution with an innitely small discretization unit, suc h as a nite elemen t or a nite dierence patc h. Unfortunately suc h a solution strategy quic kly b ecomes in tractable b ecause the computational cost asso ciated with v ery ne mesh is prohibitiv ely exp ensiv e. Therefore, a mo del that can describ e the motion on b oth lo cal and global scales is strongly desirable. P osing global constrain ts on motion mo del through appropriate motion parameters is a commonly used strategy An extreme example is the motion mo del of a rigid b o dy where lo cal motion is automatically global b ecause of the rigidit y constrain t. The ideal motion parameters that can p ose global constrain ts should p osses follo wing attributes: (1) the n um b er of parameters is manageable to reduce the computational com6 PAGE 17 plexit y; (2) parameters remain constan t during deformation; (3) parameters are w ell dened (either ph ysically or kinematically) so that their v alues can b e v eried. Because of the in timate relationship b et w een an ob ject's shap e and motion, global constrain ts are often form ulated as v arious shap e functions. A constrained deformable sup erquadric function has b een used for sim ultaneous shap e reconstruction and ob ject recognition [126 80 ]. Along the same line, a h yp erquadric function com bined with a kinematic mo del w as prop osed to reco v er the nonrigid motion without the need of corresp ondence [75 ]. The dra wbac k of using a global parametric shap e function is that the parameters are temp oral functions and ha v e to b e con tin uously up dated during trac king. A truncated series of vibration mo des has b een used to break do wn the nonrigid motion in to dieren t comp onen ts [101 ]. The small set of lo w frequency mo dal parameters allo ws for a stable estimation and represen tation of motion eld. A similar frequencybased metho d w as emplo y ed to analyze the nonrigid motion in 4D v olume images [89 ]. The main limitation of those approac hes is that the mo dal parameters do not ha v e the clear ph ysical meanings and their v alues are dicult y to v erify The aforemen tioned parameters are based primarily on the kinematic c haracteristics of motion. The adv an tage of kinematic parameters is their geometrical in tuitiv eness and generalit y The main disadv an tage is the lac k of ph ysical descriptions of motion. F or instance, forces that cause the motion as w ell as the ob ject's in ternal resp onse are not accoun ted in the kinematic parameters. Ob jects b eha v e dieren tly under the same loading condition. An in tuitiv e explanation is that an ob ject's capabilit y to resist b eing deformed is related to its unique material constitution. An ob ject's capabilit y to resist external forces can b e view ed as a sp ecial in ternal constrain t that is automatically imp osed on its motion equations. This t yp e of constrain t can b e quan tied through a constitutiv e equation, whic h is in turn determined b y an ob ject's in trinsic material prop erties. Imp osing motion constrain ts based on material c haracteristics through the constitutiv e equation is ph ysically plausible and is consisten t with the practice of con tin uum mec hanics. 7 PAGE 18 Tsap et al [136 ] ha v e utilized the reco v ered elasticit y to assist nonrigid motion mo deling. The reco v ery algorithm is based on an one dimensional gradien t decen t metho d. The motion analysis is carried out within a m ultila y er mesh renemen t framew ork. The metho d has also b een applied to burn scar study [135 ]. Mo deling nonrigid motion using a full scale ph ysical mo del has b een rep orted in the study of ruid motion [82 ], all though no attempt has b een made to reco v er the material prop ert y of ruid suc h as the viscosit y In this study based on thorough analysis of the c haracteristics of motion equation and mathematical justication of prop ert y con trolled adaptiv e motion smo othing, the Y oung's mo dulus is iden tied as an unique ph ysical motion parameter to p ose con tin uit y constrain ts on the deformation of elastic ob jects. As a ph ysical motion parameter, the Y oung's mo dulus has the follo wing desirable attributes: 1. As a distributed parameter, the Y oung's mo dulus can b e assigned to eac h elemen t of a nite elemen t mo del. A t the rst glance, it seems an exp ensiv e option for a mo del of man y elemen ts. Ho w ev er, biological ob jects usually ha v e prop ert y abnormalities in a few isolated areas, and the bac kground can b e regarded as homogeneous. Therefore, an ob ject's elastic prop ert y can b e adequately represen ted b y a few distinct v alues. 2. The Y oung's mo dulus is a function of temp erature in thermo dynamics. F ortunately isothermal condition is v alid in most cases. The v ariation of Y oung's mo dulus of biological materials within the normal range of b o dy temp erature is negligible and th us can b e viewd as a c onstant This feature distinguishes the Y oung's mo dulus from kinematic parameters. In long sequence trac king and medical image registration, the Y oung's mo dulus only needs to b e reco v ered once. 3. The Y oung's mo dulus has a clear ph ysical in terpretation in elasticit y theory and its v alue can b e v eried. The ground truth Y oung's mo dulus of biological tissues can b e gathered from standard tensile test [41 ]. 8 PAGE 19 1.3 Con tributions In summary the main con tributions of this w ork are: 1. In the con text of nonrigid motion analysis, the Y oung's mo dulus is iden tied as a unique ph ysical motion parameter, rather than a regular v ariable in a forw ard nite elemen t mo del. The uniqueness of using the Y oung's mo dulus as a motion parameter lies in the fact that it essen tially con trols the distribution of motion eld with a mec hanism that is similar to the regularized adaptiv e smo othing. More imp ortan tly the Y oung's mo dulus p oses an implicit con tin uit y constrain t on the motion inside an ob ject, whic h allo ws us to syn thesize the in terior motion from b oundary observ ations b y solving a w ellp osed Diric hlet t yp e b oundary v alue problem. 2. Tw o new reco v ery algorithms are dev elop ed to solv e the illp osed nonlinear in v erse elastic problem, one is a deterministic approac h based on the GaussNewton minimization and the general cross v alidation selection (GCV), and the other is a sto c hastic approac h based on the constrained ev olutionary computation. The coupling of GaussNewton iteration and GCV searc h impro v es b oth the eciency and the robustness of the algorithm. The constrained genetic algorithm pro v es to b e v aluable in the case where a go o d initial condition is not a v ailable, esp ecially for large scale mo dels. 3. A sensitivit y metho d is dev elop ed to impro v e the mo deling accuracy The metho d allo ws us to quan tify the impact of dieren t parameters using the sensitivit y maps. The most vulnerable areas in a mo del can b e iden tied with the aid of sensitivit y analysis. The algorithm uses the adjoin t state metho d and therefore is suitable for large scale nite elemen t mo dels. 4. The prop osed mo deling approac h has b een applied to burn scar assessmen t and face recognition. The use of natural features enhances the applicabilit y of the curren t scar rating metho d. The exp erimen t using elastic strain pattern of prole images indicates that the prop osed new biometrics holds great p oten tial in face recognition. 9 PAGE 20 CHAPTER 2 RECO VER Y ALGORITHMS 2.1 F orw ard Mo del: Linear Elastic Deformation The dynamics of a deformable b o dy dened in a b ounded domain (n R n ; @ n = 1 + 2 ) under the isothermal condition is go v erned b y the follo wing partial dieren tial equation with the Cauc h y data C d [41 132 24 ]: @ 2 u @ t 2 = r T + F ; in n ; (1) C d = f u j 1 = u ; @ u @ n 2 = g g ; on @ n = 1 + 2 : (2) where u ( u; v ; w ) denotes the displacemen t v ector of a p oin t at ( x; y ; z ), F ( f x ; f y ; f z ) is the b o dy force, is the mass densit y is the stress tensor, T denotes transp ose, r is the div ergence op erator with resp ect to a tensor, n is the out w ard unit normal on the b oundary and ( u ; g ) are the Diric hlet and Neumann data on the b oundary ( 1 ; 2 ), resp ectiv ely The motion equation of elastic b o dy is simplied with t w o linear conditions. A Cauc h y strain tensor (geometrically linear) is used, = 1 2 [ r u + ( r u ) T ] ; (3) where the gradien t op erator r is dened with resp ect to the displacemen t v ector: r u = 2 6 6 6 6 6 4 @ u @ x @ u @ y @ u @ z @ v @ x @ v @ y @ v @ z @ w @ x @ w @ y @ w @ z 3 7 7 7 7 7 5 ; (4) 10 PAGE 21 and consequen tly strain tensor can also b e expressed in 3D Cartesian co ordinate as: = 2 6 6 6 6 6 4 @ u @ x 1 2 ( @ u @ y + @ v @ x ) 1 2 ( @ u @ z + @ w @ x ) 1 2 ( @ v @ x + @ u @ y ) @ v @ y 1 2 ( @ v @ z + @ w @ y ) 1 2 ( @ w @ x + @ u @ z ) 1 2 ( @ w @ y + @ v @ z ) @ w @ z 3 7 7 7 7 7 5 : (5) Using the linear constitutiv e la w (material linear, the generalized Ho ok e's la w) and isotropic prop ert y the follo wing stressstrain relationship can b e deriv ed: = ( tr ) I + 2 = ( r u ) I + r u + ( r u ) T ; (6) where I is the iden tit y matrix, tr denotes trace, and are the Lam e constan ts. The symmetry of stress tensor ( = T ) is automatically deriv ed from the symmetry of strain tensor ( = T ). With ab o v e linear conditions, the go v erning motion equations can rewritten in terms of displacemen t that describ es the deformation of an inhomogeneous, isotropic and linear elastic ob ject: @ 2 u @ t 2 = r [ ( r u ) I + r u + ( r u ) T ] + F ; in n ; (7) C d = f u j 1 = u ; @ u @ n 2 = g g ; on @ n = 1 + 2 : (8) Material prop erties more commonly kno wn in the engineering literature suc h as the Y oung's mo dulus ( E ) and the P oisson's ratio ( ) are related to the Lam e constan ts through the follo wing transformations and can b e directly substituted in to the motion equation: = E 2(1 + ) ; = E (1 + )(1 2 ) : (9) Since b oth the Y oung's mo dulus and P oisson's ratio are considered as spatial functions, the motion equations (7, 8) will not b e further simplied in to the Na vier's equation. 11 PAGE 22 2.2 Reco v ery Algorithm: Deterministic Approac h The in v ersion of a forw ard equation is often illp osed and situation gets w orse with a highly nonlinear system of elasticit y reconstruction. Ev en if the forw ard problem can b e simplied through linearization, the corresp onding in v erse problem is often nonlinear [30 ]. V arious regularization metho ds ha v e to b e used to alleviate the illp osedness, among whic h the Tikhono v functional is a p opular one [131 ]. 2.2.1 Tikhono v Regularization Considering the static case of elastic deformation: r [ ( r u ) I + r u + ( r u ) T ] + F = 0 ; (10) a dieren tial op erator B for displacemen t u with resp ect to E can b e dened: B ( E ) = r [ ( r ( )) I + r ( ) + ( r ( )) T ] : (11) Rearranging equation (10) using B ( E ), a nonlinear op erator equation can b e deriv ed: B ( E ) u + F = 0 ; (12) A ( E ) = F B ( E ) 1 = u : (13) Assuming that the nonlinear op erator A : X Y is con tin uously F r ec het dieren tiable in the Hilb ert space ( X Y ), an in v erse problem of estimating the Y oung's mo dulus from noisy data is obtained: A ( E ) = u ; (14) k u t u k ; (15) where u t is the true displacemen t, and u is the corrupted data with a noise lev el of This in v erse problem is v ery lik ely illp osed in the Hadamard sense [53 ]. The existence of an in v erse solution is the most dicult to pro v e. F ortunately for practical problems, 12 PAGE 23 although the existence of an in v erse solution cannot b e guaran teed with noisy data, at least the ph ysical realit y itself oers an solution that can b e appro ximated. In the case of the Y oung's mo dulus reco v ery w e kno w that there m ust b e sp ecic Y oung's mo dulus v alues asso ciated with eac h real ob ject. If an in v erse system is insolv able, it usually can b e found out that either the mo del conguration or b oundary sp ecication do es not mak e ph ysical sense. Therefore, w e are more concerned ab out obtaining a n umerical solution that is unique and stable, assuming that suc h a solution do es exist. The reco v ery problem is form ulated as a Tikhono v functional in its v ariational form: E = min E R n T ( E ) ; (16) T ( E ) = k D A ( E ) u k 2 + k S ( E E ) k 2 ; (17) where D is a pro jection op erator that maps the con tin uous degrees of freedom to the discrete measuremen t co ordinate, E denotes the prior kno wledge, is the regularization parameter and S is the smo othness matrix (a discretized v ersion of the gradien t or Laplacian op erators). The solution of Tikhono v functional represen ts a compromise b et w een the delities to observ ation data and prior kno wledge. It can also b e view ed as a regular leastsquare solution, sub ject to the constrain t of a small norm. 2.2.2 Baseline Solution Strategy The minimizer of the Tikhono v functional can b e found b y either a deterministic gradien t metho d or a sto c hastic global searc h. The GaussNewton metho d is used to solv e the nonlinear illp osed in v erse problem of (16). Since the con v exit y of a nonlinear Tikhono v functional can only b e guaran teed lo cally It is exp ected that the algorithm starts with a p oin t not v ery far from the true solution. In the follo wing discussions, P is used to denote the generic parameter to b e estimated, whic h could b e the Y oung's mo dulus, the P oisson's ratio or an y other quan tities of in terest. u is the data v ector that consists of b oth measured displacemen t and stress. The in ter13 PAGE 24 p olation matrix D is drop ed assuming that the outputs from forw ard mo del ha v e b een computed on the measuremen t co ordinate. The GaussNewton strategy is to linearize the nonlinear system and then construct an iterativ e gradien tbased searc h algorithm. Considering the Tikhono v functional of a nonlinear op erator A ( P ): T ( P ) = k A ( P ) u k 2 + k S ( P P ) k 2 ; (18) (18) is rewritten in terms of the linearized expression around a lo cal p oin t P k : L ( P ) = k A ( P k ) + J ( P k )( P P k ) u k 2 + k S ( P P ) k 2 ; (19) where J ( P ) is the F r ec het deriv ativ e and L denotes the linearized functional. L ( P ) is then minimized through the v anishing gradien t condition: g ( P ) = r L ( P ) = 0 ; (20) J T ( P k )( A ( P k ) + J ( P k )( P P k ) u ) + S T S ( P P ) = 0 : (21) Replacing P with P k + 1 in (21), the iterativ e xed p oin t form ula of estimating P k + 1 is obtained: P k + 1 = P k + 4 P ; (22) 4 P = [ J T ( P k ) J ( P k ) + S T S ] 1 [ J T ( P k )( u A ( P k )) S T S ( P k P )] : (23) It can b e seen that (23) is the normal equation of the follo wing linear system: 2 6 4 J ( P k ) p S 3 7 5 4 P = 2 6 4 u A ( P k ) p S ( P k P ) 3 7 5 : (24) The Conjugate Gradien t metho d is an ecien t solv er for the large scale linear system. During eac h GaussNewton iteration for minimizing the Tikhono v functional, the Conjugate Gradien t metho d is used to solv e (24). The computed incremen t 4 P will then b e used to up date the Tikhono v solution P k + 1 14 PAGE 25 The baseline algorithm for solving an illp osed elastic problem is giv en as follo ws: Baseline Algorithm: 1. select the initial solution 2. For i = 1 to maximum tests of regularization parameter 3. compute the initial Tikhonov functional with P(i) and beta(i) 4. For k = 1 to maximum Gauss_Newton iterations 5. compute the Frechet derivative J(P) 6. compute the gradient g(P) 7. if converge criteria is satisfied 8. P(i) = P(k) 9. exit Gauss_Newton iteration 10. else 11. compute the increment with the Conjugate Gradient: dP 12. update solution: P(k) = P(k1) + gamma dP 13. compute new Tikhonov functional with P(k) 14. if (new Tikhonov < old Tikhonov) 15. accept new solution: P(k) 16. update the Tikhonov functional with P(k) 17. else 18. choose a smaller step size: gamma 19. go to step 12 20. End for 21. if P(i) satisfies the discrepancy principle 22. exit 23. else 24. choose beta(i+1) using the norms and interpolation of beta(i), beta(i1), ... 25. End for The baseline algorithm emplo ys a strategy that tries out man y p ossible v alues of regularization parameter ( i ) and c ho oses the b est one that satises certain criteria suc h as the discrepancy principle [86 ]: k A ( P ( i )) u k : (25) 15 PAGE 26 An in tuitiv e in terpretation of the discrepancy principle is that it is meaningless to pursue a solution whose residual is b elo w the b ound of noise lev el. If the Gaussian data noise is assumed with N (0 ; ): u = A ( P true ) + ; (26) and solution P ( ) is close to the true solution P true then the exp ectation of residual norm can b e appro ximated as follo w: fk A ( P ( )) u k 2 g fk k 2 g = N 2 ; (27) with N as the dimension of data v ector. So, the outer lo op of the baseline algorithm should b e stopp ed as so on as the discrepancy rule is satised: j k A ( P ( )) u k 2 N 2 j < tol er ance (28) If the GaussNewton solution with a giv en do es not satisfy the discrepancy criteria, a new regularization parameter has to b e c hosen and tested. This implemen tation sc heme is relativ ely exp ensiv e, b ecause p oten tially man y dieren t v alues ha v e to b e tested, and for eac h the nonlinear Tikhono v functional (18) has to b e solv ed. If the algorithm can start with a i that is close to the optimal v alue and the residual norm is monotonic around i the in terp olation and brac k eting tec hniques can b e applied to searc h for the next i +1 using the residual and solution norms of previous v alues ( i i 1 ...), un til the optimal one is found [9 ]. The main problem with the metho ds based on the discrepancy rule is that information ab out data noise ma y not b e a v ailable in practice. Therefore, w e ha v e to lo ok in to more heuristicallybased metho ds suc h as the generalized cross v alidation (GCV) [140 50 ]. 16 PAGE 27 2.2.3 GCV Solution Strategy GCV is designed based on the lea v eoneout principle and has b een successfully applied in linear settings. Considering a linear system and its Tikhono v functional in standard form: G x = y ; (29) T ( x ) = k G x y k 2 + k x k 2 ; (30) the ab o v e functional can b e rewritted with the k th elemen t of data missing: T ( x k ) = N X i 6 = k k ( G x ) i y i k 2 + k x k 2 ; (31) where x k is the solution of minimizing T ( x k ), and N is the dimension of y The philosoph y b ehind the cross v alidation is that if an elemen t of data is left out, the corresp onding solution of minimizing (31) should b e able to predict the missing elemen t reasonably w ell. In other w ords, if the used in minimization of (31) is a go o d c hoice, the dierence b et w een the predicted and observ ed k th datum should b e small. So, a cross v alidation function can b e constructed b y summing up the squared distance for all elemen ts: C V ( ) = N X k =1 (( G x k ) k y k ) 2 ; (32) and the optimal regularization parameter can b e found b y minimizing the CV function: = min C V ( ). T o a v oid solving N equations in v olv ed in C V ( ), the follo wing form ula can b e used that has b een thoroughly studied b y W ahaba in the con text of tting spline mo del [140 ]: C V ( ) = N X k =1 j ( G x k ) k y k j 2 j 1 Q k k ( ) j 2 ; (33) where Q ( ) = G ( G T G + I ) 1 G T and x ( ) = ( G T G + I ) 1 G T y the solution of T ( x ) = 0. 17 PAGE 28 T o remo v e its undesired prop ert y of dep endence on orthogonal data transformation, (33) w as extended to the generalized cross v alidation (GCV) [50 ]: GC V ( ) = k G x ( ) y k 2 ( tr ace [ I Q ( )]) 2 ; (34) and GC V ( ) can b e minimized through either direct in v ersion or iterativ e pro cedures [51 ]. F or the nonlinear in v erse problem (14), it is noticed that its linearized Tikhono v functional: L ( P ) = k A ( P k ) + J ( P k )( P k + 1 P k ) u k 2 + k S ( P k + 1 P ) k 2 ; (35) is a regularization of the follo wing linear system: J V = d (36) with J = J ( P k ), V = P k +1 P and d = J P k J P A ( P k ) + u So, L ( P ) is rewritten as: L ( P ) = k J V d k 2 + k S V k 2 : (37) Using the same form ula as in (34), the GCV for L ( P ) is obtained as: GC V ( ) = k J V ( ) d k 2 ( tr ace [ I Q 0 ( )]) 2 ; (38) with Q 0 ( ) = J ( J T J + S T S ) 1 J T and V ( ) = ( J T J + S T S ) 1 J T d No w that the regularization parameter is estimated for a linearized functional (35) without the requiremen t of prior kno wledge of data noise, a sc heme is needed to solv e the original nonlinear Tikhono v functional (18) utilizing the GCV metho d. If w e simply place the GCV pro cedure at the end of the outer lo op of baseline algorithm, w e w ould end up with the same situation of solving a nonlinear function man y times. A more ecien t approac h is to incorp orate the GCV inside the GaussNewton iteration so that and P can b e solv ed sim ultaneously [52 ]. The basic idea is that, for eac h GaussNewton iteration, a new reg18 PAGE 29 ularization parameter corresp onding to L ( P ) in (35) will b e estimated using the GCV, and ev en tually b oth the regularization parameter and solution will con v erge. The mo died algorithm is listed in the follo wing pseudo co de: GaussNewton with GCV: 1. select the initial solution 2. For k = 1 to maximum Gauss_Newton iterations 3. compute the Frechet derivative J(P) 4. compute the gradient g(P) 5. compute beta(k) using the GCV 6. compute the increment with the Conjugate Gradient: dP 7. update solution: P(k) = P(k1) + gamma dP 8. compute new Tikhonov functional with P(k), beta(k) 9. compute old Tikhonov functional with P(k1), beta(k) 10. if (new Tikhonov < old Tikhonov) 11. accept new solution: P(k) 12. else 13. choose a smaller step size: gamma 14. go to step 7 15. if both beta(k) and and P(k) converge 16. exit 17. End for 2.2.4 Uniqueness of Solution So far, the stabilit y of in v erse solution has b een the primary concern. Situation b ecomes more complex with the non uniqueness issue, b ecause ev en in the ideal case with zero data error, certain form ulation of in v erse problem can lead to the existence of m ultiple solutions. Let us consider a one dimensional linear elastic problem with the Diric hlet condition: d dx E ( x ) du dx = d dx [ E ( x ) ( x )] = 0 ; x 2 [ x 1 ; x 2 ] : (39) 19 PAGE 30 u ( x 1 ) = u 1 ; u ( x 2 ) = u 2 : (40) where = xx = du dx = xx = E ( x ) = E ( x ) du dx u ( x ) is kno wn in the domain, implying that ( x ) is also kno wn. This is a t ypical heat diusion problem and E ( x ) cannot b e determined uniquely from u ( x ) or ( x ), b ecause in tegration of (39) giv es a solution of E ( x ) with an arbitrary constan t C : E ( x ) = C ( x ) : (41) A unique solution of E ( x ) can b e obtained if w e ha v e either the Neumann data or a direct measuremen t of E ( x ) in the domain. F or example, giv en the follo wing Neumann condition: E ( x ) du dx x 1 = g ; (42) C = g and a solution can b e obtained ev erywhere from ( x ): E ( x ) = g ( x ) : (43) Similarly giv en a measuremen t E ( x a ), x a 2 [ x 1 ; x 2 ], a unique solution with C = E ( x a ) ( x a ) is obtained, E ( x ) = E ( x a ) ( x a ) ( x ) : (44) F or a 2D elastic ob ject, the equilibrium and constitutiv e equations are as follo ws: @ xx @ x + @ xy @ y = 0 ; @ y y @ y + @ xy @ x = 0 ; (45) xx = 2 xx + ( xx + y y ) ; y y = 2 y y + ( xx + y y ) ; xy = 2 xy : (46) Com bining t w o equations in (45) and substitute ij b y ij a single go v erning equation is obtained: @ @ x [2 ( xx + xy ) + ( xx + y y )] + @ @ y [ 2 ( y y + xy ) + ( xx + y y )] = 0 : (47) 20 PAGE 31 (47) is rewritten b y replacing ( ; ) with ( E ) using (9) and obtain an equation with E as the only unkno wn (a constan t P oisson ratio is used in the range of 0.490 0.495): @ @ x [ E f 1 ( x; y )] + @ @ y [ E f 2 ( x; y )] = 0 ; (48) f 1 ( x; y ) = xx + xy 1 + + ( xx + xy ) (1 + )(1 2 ) ; f 2 ( x; y ) = y y + xy 1 + + ( xx + xy ) (1 + )(1 2 ) : (49) After dieren tiation, (48) is further transformed to: f 1 ( x; y ) @ E @ x + f 2 ( x; y ) @ E @ y = f 3 ( x; y ) E : (50) f 3 ( x; y ) = @ f 1 ( x; y ) @ x + @ f 2 ( x; y ) @ y : (51) (50) is a t ypical rst order linear partial dieren tial equation with resp ect to E T o obtain a unique solution of (50), the follo wing Cauc h y problem [109 ] has to b e solv ed: f 1 ( x; y ) @ E @ x + f 2 ( x; y ) @ E @ y = f 3 ( x; y ) E ; n 2 R 2 : (52) C = f x = h 1 ( t ) ; y = h 2 ( t ) g ; n : (53) E = h 3 ( t ) ; on C : (54) where C is the initial curv e in n, dened b y t w o giv en functions ( h 1 h 2 ). F unction h 3 on curv e C denes the Cauc h y data, and t is a realv alued parameter. The Cauc h y theorem requires that C is not c haracteristic with resp ect to (52), i.e. C should not b e tangen t to an y c haracteristic curv e of (52). The Cauc h y theorem ensures the uniqueness of solution only in a lo cal neigh b orho o d that con tains C with the Cauc h y data dened on it. This lo cal area ma y or ma y not corresp ond to the whole mo deling domain, dep ending up on the partial dieren tial equation, initial curv e and the Cauc h y data. If the Cauc h y data can b e sp ecied along the domain b oundary ( @ n), E ( x; y ) ma y b e determined uniquely pro vided that the sp ecied b oundary is not 21 PAGE 32 c haracteristic and all in tegral curv es of (52) pass through the mo deling domain. A recen t study b y Barb one and Bam b er [12 ] sho ws in teresting results on sp ecial congurations of 2D elastic mo dels. Readers are refered to [138 88 ] for indepth discussions on the uniqueness of in v erse solution through Diric hlettoNeumann mapping. The k ey p oin t is that, to obtain a unique elasticit y estimation, the kno wledge of either stress or elasticit y itself is needed in the mo deling domain (lik ely along the b oundary). Displacemen t data along will not b e sucien t to guaran tee a unique solution. If displacemen t and material prop erties are measurable on the b oundary stress/traction data can b e deriv ed and added to the data v ector ( u ) in the reco v ery algorithm. 2.3 Reco v ery Algorithm: Constrained Genetic Approac h Since the deterministic metho ds are sensitiv e to the initial conditions, the sto c hastic metho ds ha v e also b een studied recen tly In this section, a constrained genetic algorithm (CGA) will b e discussed that fo cus on the follo wing issues: (1) enco ding a 2D/3D spatial v ariable suc h as E ( x ) in a 1D genetic computational unit; (2) expressing and incorp orating qualitativ e prior kno wledge in the genetic algorithm as a rank based constrain t; (3) balancing con tributions from the measuremen t data whic h has uncertain ties and the prior kno wledge through sto c hastic ranking. 2.3.1 Genetic Enco ding Giv en the nite elemen t mo del of a deformed ob ject, its Y oung's mo dulus can b e in terpreted as a c hromosome in a CGA. The Y oung's mo dulus v alue of eac h elemen t is enco ded as a gene in the c hromosome through a onetoone mapping function (Figure 1). As a result, if the nite elemen t mo del has N elemen ts, the corresp onding c hromosome w ould ha v e N genes. Eac h c hromosome in the p opulation p o ol represen ts a p ossible Y oung' mo dulus distribution. If dynamic meshing (m ultigrid) is used in the nite elemen t mo del, more sophisticated enco ding sc hemes ha v e to b e considered that allo w the size and shap e of c hromosomes to c hange adaptiv ely during the ev olution [48 49 115 ]. 22 PAGE 33 E1 E2 E3 E4 E5 E6 E7 E8 E9 g1 g2 g3 g4 g5 g6 g7 g8 g9 Finite Element Mesh E(i) = Young's modulus of the ith element Chromosome g1 = E1 . g9 = E9Figure 1. Genetic Co ding of the Y oung's Mo dulus in a Finite Elemen t Mo del. This onetoone mapping function is straigh tforw ard to implemen t and w orks w ell for a 1D nite elemen t mo del. Ho w ev er, information ab out the spatial connections among the neigh b oring elemen ts in a 2D/3D nite elemen t mesh is completely lost in this enco ding sc heme. Computations that rely on spatial information cannot b e p erformed prop erly F or example, in man y studies, our in terest is to iden tify only a few isolated areas of abnormal Y oung's mo dulus v alues from the bac kground of relativ ely uniform Y oung's mo dulus distribution, and the information ab out genes' spatial distribution is needed in a CGA to accomplish the task ecien tly More imp ortan tly the smo othing eect represen ted b y the deriv ativ e op erators in the deterministic regularization cannot b e realized in a rankbased CGA without suc h spatial information. T o o v ercome those shortcomings of onetoone mapping, a mec hanism is designed to remem b er the original spatial connections among neigh b oring elemen ts. An auxiliary link table is created for all c hromosomes. T able 3 sho ws an example that uses a 4neigh b or connection for the 2D quadrilateral mesh in Figure 1. Similarly a 3neigh b or connection can b e considered for a triangle mesh. This link table will b e used in the constrained m utation op eration. Because the Y oung's mo dulus has 23 PAGE 34 T able 3. Link T able for Main taining Original Spatial Connections. Gene ID V alues 4Neigh b or Connections g1 E1 g2 = E2, g4 = E4 g2 E2 g1 = E1, g3 = E3, g5 = E5 g3 E3 g2 = E2, g6 = E6 g4 E4 g1 = E1, g5 = E5, g7 = E7 g5 E5 g2 = E2, g4 = E4, g6 = E6, g8 = E8 g6 E6 g3 = E3, g5 = E5, g9 = E9 g7 E7 g4 = E4, g8 = E8 g8 E8 g5 = E5, g7 = E7, g9 = E9 g9 E9 g6 = E6, g8 = E8 con tin uous v alues that can v ary in the range of sev eral orders of magnitude, a realv alued (double) enco ding approac h is used [36 ]. 2.3.2 Rank Based Constrain t In the sto c hastic framew ork of genetic computation, the illp osedness of a nonlinear in v erse problem implies man y lo cal plateaus in the landscap e of the admissible solution space. By ha ving a div erse p opulation p o ol, genetic algorithms can explore a m uc h wider solution space than the gradien t metho ds and th us ha v e a b etter c hance to nd the optimal solution b y escaping the lo cal minima. The illp osed nature of the in v erse problem also sho ws up as highly unstable solutions that are ph ysically meaningless. T o o v ercome those n umerical difculties asso ciated with the illp osedness, v arious constrain ts that represen t prior kno wledge m ust b e imp osed. In a CGA, constrain ts can tak e the form of p enalt y functions [23 ], whic h are equiv alen t to the regularization stabilizers or the preconditioners in the deterministic metho ds [30 ]. The ob jectiv e function to b e minimized b y a CGA is form ulated as a com bination of the tness function and the p enalt y function: obj ( E ) = k Q ( E ) b k 2 n + P ( E ) (55) 24 PAGE 35 where Q ( E ) denotes the output from the forw ard mo del, b is the data v ector, n is the n um b er of no des on whic h the measuremen t is made, P ( E ) is the p enalt y function, and is the w eigh t co ecien t. In studies of spline and surface tting, a smo othness constrain t imp osed on the spline and the surface is often form ulated as a quadratic in tegral functional suc h as the Sob olev norm. Similarly the Y oung's mo dulus as a distributed parameter can b e regarded as a spatial function that p ossesses certain degrees of con tin uit y and the smo othness constrain ts can b e realized as: P ( E ) = Z R jr E ( x ) j 2 dx (56) where r could b e an y deriv ativ e op erator suc h as the gradien t or the Laplacian. Kno wledge ab out the relativ e Y oung's mo dulus distribution can b e obtained from either an exp ert's visual assessmen t or from lo w lev el image cues suc h as in tensit y color and texture. This t yp e of prior kno wledge is often expressed qualitativ ely and is not suited for p enalt y functions that are comp osed of con tin uous dieren tial op erators. T o utilize the qualitativ e prior kno wledge in a genetic algorithm, An alternativ e rankbased metho d is prop osed to compute the p enalt y function. In eac h p ossible solution (c hromosome), elemen ts (genes) will b e rank ed based on their relativ e Y oung's mo dulus v alues and their p ositions will b e recorded in a sorted rank table. Similarly the qualitativ e prior kno wledge is transformed in to another rank table. F or eac h elemen t, the dierence of its rank ed p ositions in the t w o rank tables is computed. The rank discrepancy of all elemen ts is then summed up to represen t the distance b et w een the solution and the prior kno wledge: P ( E ) = m X i =1 k r i R i k 2 (57) where r i is the rank p osition of elemen t i in the rank table for the solution, R i is the rank p osition of elemen t i in the rank table based on prior kno wledge, and m is the n um b er of elemen ts. 25 PAGE 36 Figure 2 illustrates the ranking sc heme with a simple t w odimensional mo del of four elemen ts. As the qualitativ e prior kno wledge, the Y oung's mo dulus v alue of eac h elemen t is lab eled as \high", \mid" and \lo w". This qualitativ e information is then transformed in to a rank table where elemen ts are sorted in descending order from \high" to \lo w". If n elemen ts ( n > 1) ha v e the same lab el, they all can ha v e n p oten tial rank p ositions, whic h will b e determined b y their coun terparts in the solution rank table. F or instance, b oth elemen t (1) and elemen t (4) are lab eled as \lo w", therefore their rank p osition can b e either 3 or 4. In Solution 1, the ranks of elemen t (2) and elemen t (3) in the solution table matc h exactly with their ranks in the table of prior kno wledge. Ho w ev er, for elemen t (1), its rank in the solution table is 3, while its rank in prior kno wledge table is 3 or 4. In the case of m ultiple ranks, the v alue that is closest to its coun terpart in the solution table is selected, whic h is 3 for elemen t (1). Similarly for elemen t (4), 4 is selected from its m ultiple prior rank of 3 or 4. The nal p enalt y v alue for Solution 1 is zero ( P ( E ) = 0). In con trast, Solution 2 has a Y oung's mo dulus distribution that is quite dieren t from the sp ecied prior kno wledge and th us a sh ued rank table, whic h leads to a higher p enalt y v alue of 6. One p oten tial problem with this rankbased p enalt y function is that t w o solutions of the same rank table are p enalized equally ev en if their absolute Y oung's mo dulus v alues are quite dieren t. As demonstrated in Figure 2, Solution 3 has a Y oung's mo dulus distribution that is 10 times higher than that of Solution 2. But their solution rank tables are exactly the same and therefore receiv e the same amoun t of p enalt y of 6. This problem is related to the nonuniqueness of the in v erse solution and can b e resolv ed b y incorp orating b oth displacemen t and force in the data v ector b In other w ords, Solution 2 and Solution 3 will ha v e dieren t tness functions b ecause of their dieren t b and Q ( E ) v alues. This am biguit y can also b e partially resolv ed b y in tro ducing another p enalt y function that sp ecies a range (b oth the upp er b ound and the lo w er b ound) within whic h the Y oung's mo dulus of an elemen t is allo w ed to v ary This rankbased approac h has the adv an tage that it is in trinsically piecewise and helps preserv e the parameter discon tin uit y (although a strong smo othness constrain t can still b e 26 PAGE 37 Prior Knowledge low high mid low (1) (2) (3) (4) (2) (3) (1) (4) high mid low low 1 2 3,4 3,4 Prior Rank Table Element ID E (kPa) Prior Rank Element ID E (kPa) Solu. Rank Prior Rank (2) (3) (1) (4) 40 16 8 3 1 2 3 4 1 2 3 4 Solution 1 8 40 3 16 (1) (2) (3) (4) p(E) = 33 + 11 + 22 + 44 = 0 2 2 2 2 Penalty of Solution 1 Solution 2 4 22 6 35 (1) (2) (3) (4) Element ID E (kPa) Solu. Rank Prior Rank (4) (2) (3) (1) 35 22 6 4 1 2 3 4 3 1 2 4 p(E) = 44 + 21 + 32 + 13 = 6 2 2 2 2 Penalty of Solution 2 Element ID E (kPa) Solu. Rank Prior Rank (4) (2) (3) (1) 350 220 60 40 1 2 3 4 3 1 2 4 Solution 3 (1) (2) (3) (4) 40 220 60 350 p(E) = 44 + 21 + 32 + 13 = 6 2 2 2 2 Penalty of Solution 3Figure 2. The Rank Based Sc heme for Incorp orating Qualitativ e Prior Information. The numb ers in the p ar entheses denote element ID and the numb ers at the lowerright c orner of elements r epr esent the Y oung's mo dulus. 27 PAGE 38 imp osed in the areas of little Y oung's mo dulus v ariation). This approac h is particularly suitable for studies that aim at iden tifying and quan tifying prop ert y abnormalities. In our previous studies on burn scar assessmen t [135 ], qualitativ e prior kno wledge w as collected from ph ysicians who isolate and rate the scars based on a relativ e rating scale. Automatic metho ds for extracting the information directly from image in tensit y texture and color can also b e considered. 2.3.3 Balancing Fitness and P enalt y b y Sto c hastic Ranking In the constrained ob jectiv e function (55), the tness and p enalt y terms are computed on dieren t quan tities. The tness is measured as dierence of displacemen t (meter), while the p enalt y is based on the dierence of rank orders (unitless). An optimal w eigh t co ecien t is needed to balance their con tributions. If is to o small, the data noise will not b e p enalized enough and the resulting solution b ecomes unstable. If is to o large, solution will b e forced in to a smo othed prior space and most of the data signals, i.e. the information used to infer the Y oung's mo dulus, will b e lost due to o v erp enalization. In the deterministic domain, sev eral c hoices are a v ailable for determining the optimal regularization parameter. If the noise lev el of observ ation data is kno wn, metho ds based on the discrepancy principle [86 ] suc h as the Miller metho d [84 ] can b e considered. In case that information ab out the data noise is not a v ailable, heuristic metho ds suc h as the generalized cross v alidation [140 ] or the Lcurv e [55 ] metho d are commonly used. Ho w ev er, those metho ds are deterministic and not suited for handling rankbased qualitativ e data. In a recen t study on the constrained ev olutionary optimization, Runarsson and Y ao [110 ] presen ted a sto c hastic metho d to strik e a balance b et w een the tness and p enalt y functions, without the need of computing the w eigh t co ecien t explicitly Determination of an optimal is related to the dynamic and adaptiv e ranking of individual c hromosome in a p opulation. Ranking is based on the relativ e dominance of either the tness function or the p enalt y function b et w een t w o adjacen t individuals. The balance of dominance is ac hiev ed b y in tro ducing a bubblesortlik e dynamic ranking pro cedure for an individual 28 PAGE 39 to win a comparison. It is found that this sto c hastic metho d is w ell suited for handling qualitativ e prior kno wledge in elasticit y reconstruction. Readers are refered to [110 ] for detailed discussions on the metho d and related implemen tation issues. 2.3.4 Genetic Op erators T o minimize the ob jectiv e function (55), Sev eral imp ortan t genetic op erators are sp ecied. Mutation: A standard Gaussian m utation op erator is used, g j = g j + N (0 ; 1) g j (58) where g j is the v alue of gene j b efore the m utation, g j is the v alue of gene j after the m utation. N (0 ; 1) is a random Gaussian n um b er (mean = 0, standard deviation = 1) and is the m utation step size. The dynamic con trol of m utation step size is determined b y a predened deca y rate as: ( k + 1) = ( k ) ; (59) where k is the generation coun ter. is usually set in the range of 0.99 1.0. A smaller v alue of ma y help sp eed up the con v ergence rate but has the risk of premature con v ergence. Other than Gaussian m utation, other t yp es of realco ded m utation op erators are also considered suc h as uniform m utation and step m utation. Those m utation approac hes did not sho w an y adv an tage o v er Gaussian m utation in terms of b oth solution accuracy and computation eciency Exp erimen ts w ere conducted with the metho d of dynamically setting the m utation probabilit y ( P m ) based on p opulation statistics and no signican t impro v emen t w as found in terms of the solution accuracy One of the purp oses of ha ving a relativ ely high m utation rate is to main tain the p opulation div ersit y and prev en t premature con v ergence. T o utilize the spatial connectivit y among nite elemen ts as recorded in the link table (T able 3), eac h gene will b e compared with its neigh b oring genes after the m utation. If the gene has a Y oung's mo dulus v alue that is m uc h higher/lo w er than the maxm um/minim um of its neigh b oring 29 PAGE 40 genes, the gene will b e assigned a mean v alue of its neigh b ors. This pro cess has the eect of partially restoring the smo othness prop ert y of the quan titativ e regularization functional that is missed in the rankbased p enalt y function (57). Cr ossover: It is also found that the onep oin t crosso v er op erator and the m ultiplep oin t crosso v er op erator p erformed equally w ell, at least for this particular in v erse problem setting. The crosso v er probabilit y ( P c ) is xed to 0.7. The onep oin t crosso v er function is implemen ted in a traditional fashion: c hildren are generated b y joining t w o paren ts at a randomly selected crosso v er p osition and then sw apping eac h sides. Par ent Sele ction and R eplac ement: The paren t selection op erator is implemen ted as tournamen t selection (k=2). Exp erimen ts w ere conducted with a wide range of replacemen t ratios (0.05 1.0), i.e., the p ercen tage of paren t p opulation to b e replaced b y new c hromosomes. Giv en a p opulation of size S and a replacemen t ratio of r the n um b er of paren ts to b e replaced is: N = r S Smaller replacemen t ratios ( < 0.3) did not yield satisfying results b ecause of the lac k of con tributions from new c hromosomes to the p opulation div ersit y Ho w ev er, no signican t dierence w as observ ed with replacemen t ratios ranging from 0.5 to 1.0. On the other hand, the higher the replacemen t ratio, the longer the sim ulation time. In [110 ], it w as sho wn that the sto c hastic ranking probabilit y ( P f ) in the range of (0.45 0.475) ga v e the b est result for 13 b enc hmark functions. Tw o P f v alues (0.45 and 0.475) w ere tested in exp erimen ts, and no noticeable impro v emen t or degradation w as observ ed in the algorithm's p erformance. So a P f v alue of 0.45 is used in all of the follo wing comparison studies. An elitism strategy is also enabled during the replacemen t op eration, in whic h some elite mem b ers of the old generation are c hosen to surviv e to the next generation without comp etition (p oten tially b eing replaced b y a b etter ospring). Giv en the genetic parameters sp ecied ab o v e, trial tests with dieren t elite ratios (1% 15%) sho w ed that 3% ga v e a sligh tly b etter result on a v erage (although v ery marginal). So an elite ratio of 3% is c hosen in all the exp erimen ts. The pseudo co de for the constrained genetic algorithm is as follo ws: 30 PAGE 41 1 S /* p opulation size */ 2 r = 0 : 8 /* r eplac ement r atio */ 3 = 0 : 997 /* mutation de c ay r ate */ 4 = 0 : 2 /* initial mutation step size */ 5 P c = 0 : 2 /* onep oint cr ossover pr ob ability */ 6 P f = 0 : 45 /* sto chastic r anking pr ob ability */ 7 for i = 1 to p opulation S(0) 8 c hromosome (i) = a v alue from a Gaussian distribution 9 end for 10 for k = 1 to maxim um generation 11 for i = 1 to p opulation S (k) 12 E (i) = c hromosome (i) 13 call forw ard nite elemen t mo del with E (i) 14 compute tness and p enalt y functions for c hromosome (i) 15 end 16 rank the whole p opulation S (k) with sto c hastic ranking ( P f ) 17 pic k a set of paren ts from p opulation S (k) through tournamen t selection 18 use onep oin t crosso v er to generate N ospring with a probabilit y of P c 19 for i = 1 to N ospring 20 Mutate c hromosome(i): g = g + N (0 ; 1) g 21 E (i) = c hromosome (i) 22 call forw ard nite elemen t mo del with E (i) 23 compute tness and p enalt y functions for c hromosome (i) 24 end 25 pic k N paren ts from p opulation S (k) through tournamen t selection 26 rank N pic k ed paren ts and N ospring with sto c hastic ranking ( P f ) 27 for i = 1 to N 28 for j = 1 to N 29 if rank of paren t (i) rank of ospring (j) 30 then 31 replace paren t (i) with ospring (j) 32 remo v e ospring (j) 33 end 34 end 35 a new p opulation S (k+1) 36 up date m utation step size: (k+1) = (k) 37 end 31 PAGE 42 2.4 Deterministic or Sto c hastic Approac h 2.4.1 Exp erimen ts with GA, CGA and GNM A syn thetic n umerical mo del is used to compare the ecacy of three reco v ery algorithms: a regular genetic algorithm (GA), the prop osed constrained genetic algorithm (CGA), and the deterministic GaussNewton metho ds (GNM). The forw ard mo del is a t w o dimensional thin shell (5 cm b y 5 cm) and is discretized to 61 no des and 100 triangle elemen ts (Figure 3). A small square at the upp erleft corner of the mo del is the abnormal area and has a higher Y oung's mo dulus v alue of 250 kP a. The rest of bac kground elemen ts ha v e a uniform Y oung's mo dulus of 50 kP a. The elemen t t yp e is 3no de triangular shell without outofplane deformation ( w comp onen t in z direction). A linear sc heme is used in the in terp olation (shap e) functions: u ( x; y ) = H i u i + H j u j + H k u k ; (60) v ( x; y ) = H i v i + H j v j + H k v k ; (61) H i ( x; y ) = 1 2 A [( x j y k x k y j ) + ( y j y k ) x + ( x k x j ) y ] ; (62) H j ( x; y ) = 1 2 A [( x k y i x i y k ) + ( y k y i ) x + ( x i x k ) y ] ; (63) H k ( x; y ) = 1 2 A [( x i y j x j y i ) + ( y i y j ) x + ( x j x i ) y ] ; (64) where H denotes the in terp olation function, A is the area of triangle elemen t (see Figure 3 (c)), ( i; j; k ) are used to index three no dal p ositions and ( x; y ) represen t the lo cation in a global 2D co ordinate. Boundary conditions are sp ecied as follo ws: F orces (concen trated loads) are applied to the no des of top b oundary (1.2N on eac h no de). Displacemen t constrain ts (Diric hlet t yp e) are xed to zero on the no des of b ottom b oundary (see Figure 3 (a)). Those b oundary conditions are used to generate noisefree displacemen t data on eac h no de. As a result, the data v ector b includes forces of top b oundary no des as w ell as displacemen ts of all no des other than those on the b ottom b oundary 32 PAGE 43 250 250 250 250 50 (1.2 N) Fixed Boundary Nodal Forces High High High High Low (x i y i )(ui vi) (uj vj)(x j y j ) (x k, y k )(uk vk) (u v)(x,y) (a) F orw ard mo del (b) Prior kno wledge (c) T riangle elemen t Figure 3. Syn thetic Mo del and the Qualitativ e Prior Kno wledge. In (a), the Y oung's mo dulus of the b ackgr ound ar e a is 50 kPa. In (b), the b ackgr ound ar e a is lab ele d as \low". (c) il lustr ates the line ar interp olation scheme of triangle element use d in e quations (31) to (35). ( u; v ) ar e displac ement c omp onents in a 2D c o or dinate, ( x; y ) denotes the lo c ation of a p oint and ( i; j; k ) ar e use d to index thr e e no des. Three algorithms (GA, CGA and GNM) w ere tested with the noisefree data to study to what degree the Y oung's mo dulus in the abnormal area can b e repro duced. The data v ector includes b oth displacemen t and force to ensure the uniqueness of the in v erse solution. As prior kno wledge, four elemen ts in the upp erleft square w ere iden tied as a p oten tial abnormal area. These four elemen ts w ere lab eled as \high" and other elemen ts in the bac kground w ere lab eled as \lo w" (Figure 3). 3% white noise w as then added to the noisefree data. Again, the GA, CGA and GNM w ere studied using the noisy data to assess the p erformance. F or this small 2D mo del, w e found that a p opulation size of 50 is large enough for the GA and CGA to con v erge to a go o d solution with less than 300 generations. The p erformance of the three algorithms w as ev aluated based on the shap e and absolute Y oung's mo dulus v alue of the reconstructed abnormal area. F or the GA and CGA exp erimen ts, w e ran 30 sim ulations for eac h case and therefore the conclusions are dra wn from the a v erage of 30 sim ulations. The reconstructed Y oung's mo dulus using the three algorithms are sho wn in Figure 4. Giv en the noisefree data, all algorithms w ere able to nd a nearoptimal solution with the abnormal area w ell dened in terms of b oth the geometry and absolute v alues. The results from the CGA and GNM are almost iden tical. Giv en the noise corrupted data, b oth 33 PAGE 44 (a) GNM with 0% data noise (b) GNM with 3% data noise (c) CGA with 0% data noise (d) CGA with 3% data noise (e) GA with 0% data noise (f ) GA with 3% data noise Figure 4. Reconstructed Y oung's Mo dulus b y GA, CGA and GNM. The gur es ar e plotte d using gr eysc ale, with black and white indic ating high and low Y oung's mo dulus values, r esp e ctively. The gur es wer e gener ate d using the GMV (the Gener al Mesh Viewer) visualization softwar e [47 ] that may add c ertain de gr e es of smo othing ee ct on the Y oung's mo dulus distribution (the Y oung's mo dulus value inside an element is c onstant). 34 PAGE 45 GNM and CGA still successfully iden tied the abnormal area, although the b oundaries w ere blurred due to the smo othing eect. Ho w ev er, the GA failed to nd an acceptable solution as evidenced b y the m ultiple abnormal areas (Figure 4 (f )) that do not exist in the forw ard mo del. A Gaussian distribution with a mean of 50 kP a w as used as the initial condition for all the exp erimen ts. In the exp erimen ts of GNM, the regularization parameter w as set to 1.93E8, and the smo othness matrix w as the zero order nite dierence matrix, i.e. the iden tit y matrix. If necessary higher order matrix can b e used. F or example, for an ev enly spaced solution v ector of 4 elemen ts, the rst order nite dierence matrix is: 2 6 6 6 6 6 6 6 6 4 1 1 0 0 1 2 1 0 0 1 2 1 0 0 1 1 3 7 7 7 7 7 7 7 7 5 : (65) Since the higher order matrices ha v e a strong tendency to o v ersmo oth the v ector, w e used the iden tit y matrix in all the exp erimen ts. The sim ulation w as stopp ed when either a maxim um GaussNewton iteration n um b er is reac hed (100 for this exp erimen t) or a stopping criteria is met. The stopping rule is dened as follo w: conv er g ence = k g r ad ( E ) k = k ( E k ) k < tol er ance (66) where ( E k ) is the linearized Tikhono v functional with resp ect to the curren t solution E k and g r ad ( E ) is the gradien t op erator of Tikhono v functional: g r ad ( E ) = D T ( E k )( Q ( E k ) + D ( E k )( E E k ) b ) + W T W ( E E pr ior ) : (67) The con v ergence tolerance w as set to 1.9e10 for all the exp erimen ts with GNM, including those in the next section. This tolerance v alue is relativ ely stringen t. A go o d initial 35 PAGE 46 high high high high (a) Biased prior kno wledge (d) Result of CGA with 3% data noise Figure 5. Reconstructed Y oung's Mo dulus b y the CGA with Biased Prior. The same noisy data (3% noise) was use d. The prior lab els four elements as \high", which ar e inc onsistent with the true solution. The nal solution is cle arly inruenc e d by the biase d prior. solution guess usually allo ws the algorithm to con v erge with less than 20 iteration steps, while a p o or initial p oin t results in the maxim um iteration n um b er to b e reac hed (div erge). One imp ortan t issue in using the rankbased CGA for the Y oung's mo dulus reconstruction is the qualit y of prior kno wledge. F or example, if the prior kno wledge used in the algorithm do es not corresp ond to the true Y oung's mo dulus distribution v ery w ell, ho w m uc h will the nal solution b e aected? An exp erimen t w as carried out with sligh tly biased prior kno wledge. The mo del conguration, initial condition and noisy data w ere the same as those in the previous exp erimen ts. The biased prior and the reconstructed Y oung's mo dulus are sho wn in Figure 5. A rectangle of four elemen ts w as iden tied as the p oten tial abnormal area and w as lab eled "high" incorrectly but still close to the true abnormalit y square. It can b e seen that the shap e of the reconstructed abnormal area is inruenced b y the biased prior, although the CGA still managed to nd the true abnormalit y square. The prior kno wledge has a tendency to attract the in v erse solution to w ards itself. In applications, w e do not exp ect the prior kno wledge to b e sev erely biased, but its sideeect should not b e ignored and m ust b e corrected as m uc h as p ossible. It is w orth noting that there are t w o basic t yp es of prior in v olv ed in the Y oung's mo dulus reconstruction with the CGA. One is the smo othness prior that is implicit in the quadratic 36 PAGE 47 Nodal Force (0.5 N) Fixed Boundary 400 kPa50 kPa i j k ls = 1 s t t = 1 t = 1 s = 1 (a) F orw ard mo del (b) Quadrilateral elemen t Figure 6. Conguration of the F orw ard Mo del. numb er of elements = 900, numb er of no des = 961, numb er of elements in the abnormal ar e a = 60. (b) il lustr ates the interp olation scheme of quadrilater al element use d in the e quations (39) to (44). ( u; v ) ar e displac ement c omp onents, ( x; y ) denotes the lo c ation in glob al c o or dinate system, ( s; t ) r epr esents the lo c al c o or dinates within an element, and ( i; j; k ; l ) ar e use d to index four no des. p enalt y term, and the other one is explicitly sp ecied as the rank ed prior table. The former is a global constrain t while the later is lo cal in nature. Therefore, the inruence of a bad lo cal prior cannot b e remedied b y adjusting the global regularization parameter or sto c hastic ranking in CGA. More study is needed to nd more sophisticated and eectiv e metho ds so that the biased prior kno wledge can b e handled prop erly 2.4.2 Exp erimen ts with CGA and GNM There ha v e b een debates on whether the genetic algorithms ha v e an y adv an tage o v er the traditional gradien t decen t metho ds, b ecause the computational eciency of the genetic algorithms is often a concern, esp ecially for largescale problems. T o understand ho w v arious factors migh t aect the p erformance of a CGA in solving in v erse problems, w e carried out another set of exp erimen ts using b oth the CGA and GNM. 37 PAGE 48 The forw ard mo del is similar to the one that w e used in the previous exp erimen ts, except a m uc h larger parameter space of 900 elemen ts (Figure 6), and hence it is more dicult to nd an optimal solution of the corresp onding in v erse problem. The mo del has a size of 10 b y 10 (cm) and is meshed with quadrilateral elemen ts. 60 elemen ts in the cen ter of the mo del are the abnormal area with a Y oung's mo dulus v alue of 400 kP a, and the bac kground elemen ts ha v e a lo w er Y oung's mo dulus v alue of 50 kP a. The elemen t t yp e used is 4no des quadrilateral shell. Isoparametric mapping and the Lagrange form ula are used in the in terp olation (shap e) functions: u ( x; y ) = H i u i + H j u j + H k u k + H l u l ; (68) v ( x; y ) = H i v i + H j v j + H k v k + H l v l ; (69) H i ( s; t ) = 1 4 (1 s )(1 t ) ; (70) H j ( s; t ) = 1 4 (1 + s )(1 t ) ; (71) H k ( s; t ) = 1 4 (1 + s )(1 + t ) ; (72) H l ( s; t ) = 1 4 (1 s )(1 + t ) ; (73) where H denotes the in terp olation function, ( i; j; k ; l ) are used to index four no dal p ositions, ( s; t ) denote the lo cal co ordinates of quadrilateral elemen t, and ( x; y ) represen t the global co ordinates (see Figure 6 (b)). The b oundary conditions w ere sp ecied as follo ws: The mo del w as compressed b y forces exerted on the top b oundary (0.5N on eac h no de) while the b ottom b oundary w as xed (Diric hlet t yp e with all displacemen t comp onen ts set to zero). First, an appropriate p opulation size needs to b e determined. A series of exp erimen ts w as conducted with dieren t p opulation sizes (10, 20, 50, 100, 500, 1000). The CGA sim ulation w as stopp ed whenev er the ob jectiv e function reac hed a preset minim um (0.03 in all exp erimen ts). Eac h sim ulation w as started with a go o d initial guess to ensure that the solutions obtained are equally mec hanically meaningful (close to the true solution). F or 38 PAGE 49 T able 4. CGA Exp erimen ts with Dieren t P opulation Sizes. P opulation Size 10 20 50 100 500 1000 Generations 1563 825 637 559 430 195 Time (min) 17.6 24.5 44.2 82.3 286.7 339.5 Results are rep orted as the a v erage of 30 runs for eac h p opulation size. All sim ulations w ere started with a go o d initial p oin t and stopp ed when the ob jectiv e function reac hed a preset v alue. Figure 7. Con v ergence Beha vior of CGA with Dieren t P opulation Sizes. Each p oint on the curve is plotte d using the aver age value of 30 simulations. eac h c hosen p opulation size, w e ran 30 sim ulations and the a v erage results are listed in T able 4 and plotted in Figure 7. It is clear that the CGA with a larger p opulation size can reac h a sp ecied minim um with man y less generations than the CGA with a smaller p opulation size, but requires m uc h more sim ulation time. In other w ords, for the CGA with a larger p opulation size, eac h generation tak es more time due to the computational cost asso ciated with its large p o ol. Although the CGA with a smaller p opulation size uses m uc h less execution time, this b enet is often oset b y the concern of premature con v ergence due to a less div ersied p o ol that is common to smaller p opulations. T aking the ab o v e factors in to accoun t, a p opulation size of 100 w as c hosen that can main tain enough p opulation div ersit y y et still has a \reasonable" computational cost. 39 PAGE 50 T able 5. P erformance of CGA and GNM with Dieren t Initial Conditions. GNM CGA Initial Num. Con v erg. Av erage Av erage Num. Con v erg. Av erage Av erage Condition of Ratio Iterat. Con v erg. of Ratio Generat. Con v erg. N ( ; ) Sim ul. Steps Time Sim ul. Coun t Time (min) (min) N (50, 10) 100 95% 8 9 30 100% 315 45 N (100,50) 100 82% 17 20 30 97% 549 76 N (200,50) 100 44% 25 33 30 93% 834 137 N (300,50) 100 17% 34 39 30 91% 1494 251 N (400,50) 100 3.8% 52 61 30 88% 3483 582 N (500,50) 100 4.1% 48 56 30 84% 7347 1249 F or e ach of the 6 initial c ondition distributions (r ow), we r an 100 GNM simulations and 30 CGA simulations fr om which the c onver genc e r atio, the aver age iter ation steps (gener ation c ount) and the aver age c onver genc e time wer e c ompute d. The c onver genc e r atio was c ompute d as the numb er of c onver ge d simulations divide d by the total numb er of simulations. Only the r esults of c onver ge d simulations wer e use d in the c omputation of the aver age iter ation steps (gener ation c ount) and the aver age c onver genc e time. Note that the true solution is 400 kPa in the abnormal ar e a, and 50 kPa in the b ackgr ound, with a me an of 73 kPa. The b eha vior of a CGA that has M c hromosomes in its p o ol can b e roughly view ed as b eing equiv alen t to M lo cal gradien t metho ds that are run sim ultaneously The dierence is that the c hromosomes in a CGA ha v e certain degrees of information sharing among themselv es through the crosso v er op eration, while the gradien t metho ds running in parallel are completely indep enden t of eac h other. It w ould b e in teresting to compare the p erformance of a CGA of M c hromosomes with that of M indep enden t gradien t metho ds. The gradien t metho ds are kno wn to b e sensitiv e to the initial conditions. A guess is used that ob eys the normal distribution la w, N ( ; ). In other w ords, the initial Y oung's mo dulus v alues assigned to 900 elemen ts should ha v e a sample mean of and a standard deviation of The exp erimen ts used 6 suc h initial conditions (unit = kP a): N (50 ; 10), N (100 ; 50), N (200 ; 50), N (300 ; 50), N (400 ; 50), N (500 ; 50), in the order of gradually deviating from the true solution (T able 5). All of the exp erimen ts w ere p erformed on a Sun Sparc ultra 5 mac hine (248 MHz, 2560 Mb). 40 PAGE 51 The exp erimen ts are designed as follo ws: for eac h of the 6 initial condition distributions, w e ran 30 CGA sim ulations and 100 GNM sim ulations and rep ort the results in T able 5. So, eac h ro w in T able 5 represen ts the a v erage p erformance of 30 CGA sim ulations and 100 GNM sim ulations, all with their initial guesses generated b y the same normal distribution function. A GNM sim ulation w as stopp ed if the con v ergence criteria (37) w as satised with a sp ecied tolerance v alue (1.9e10) or a preset maxim um iteration n um b er w as reac hed (400 in this exp erimen t). Similarly a CGA sim ulation w as stopp ed if the ob jectiv e function v alue w as b elo w a threshold (0.05) or the c hange of ob jectiv e function b et w een t w o consecutiv e generations w as b elo w a threshold (0.0005). The t w o algorithms w ere then ev aluated based on three p erformance indexes: Con v ergence ratio, iteration steps (generation coun t) and con v ergence time. The c onver genc e r atio is dened as the n um b er of con v erged sim ulations divided b y the total n um b er of sim ulations. F or the GNM, the con v ergence ratio dropp ed rapidly as the initial condition mo v ed a w a y from the true solution (the third column in T able 5). F or example, starting with t w o relativ ely bad initial conditions ( N (400 ; 50) and N (500 ; 50)), most of the GNM sim ulations failed to con v erge (or con v erged to a lo cal minim um far a w a y from the true solution). In con trast, the con v ergence ratio of the CGA is more stable. The p erformance deterioration of the CGA caused b y a p o or initialization is m uc h less drastic (the sev en th column in T able 5). The gain in the con v ergence ratio of a CGA comes with a price of longer c onver genc e time A single CGA sim ulation to ok m uc h more time to con v erge than a single GNM simulation did (the fth and nin th columns in T able 5). On the other hand, if w e view the outcome of a single CGA of 100 c hromosomes as b eing equal to that of 100 indep enden t GNMs, then running 100 indep enden t GNMs seems more exp ensiv e than running a single CGA. Ho w ev er, this comparison could b e misleading. It should b e stressed that the b ottlenec k of a GNM sim ulation is the computation of the Jacobian matrix. In our curren t implemen tation, a nite dierence appro ximation metho d is used, whic h implies that for a problem of m parameters, at least m 2 calls of the forw ard mo del w ere needed for eac h 41 PAGE 52 iteration step. This accoun ts for ab out 92% of the total GNM sim ulation time. If a more ecien t metho d suc h as the adjoin t state metho d is used [90 ], the sim ulation time of a GNM can b e signican tly reduced, and the adv an tage of a single CGA against 100 indep enden t GNMs will lik ely disapp ear. F urthermore, if w e are dealing with a complex 3D problem, a larger p opulation size migh t b e needed, whic h will also increase the computational cost of the CGA. The iter ation steps and gener ation c ount sho w similar trends in that they increased gradually as the initial condition got w orse. But they are often implemen tationdep enden t and th us are less informativ e ab out the algorithm's b eha vior than the con v ergence ratio and con v ergence time. F rom the ab o v e exp erimen tal results and analysis, sev eral observ ations can b e made: 1. The CGA is more robust than the GNM in the sense that the CGA is less sensitiv e to the initial conditions. 2. A single GNM is more ecien t than a single CGA in terms of the con v ergence time. If a go o d initialization is a v ailable, the GNM is a more attractiv e c hoice. 3. A b etter solution strategy is to com bine the GNM and the CGA in a w a y that b oth their strengths can b e fully utilized. F or instance, the CGA can b e used to nd a relativ ely go o d initial appro ximation, up on whic h the GNM can b e launc hed to sp eed up the searc h. T o illustrate the con v ergence b eha viors of the CGA and the GNM, the minimization steps of four example runs w ere plotted in Figure 8. Since it is dicult to visualize the solution landscap e in high dimensions, the space w as simplied in the follo wing w a y: F or eac h solution v ector, t w o a v eraged v alues w ere computed, one for the abnormal area and one for the bac kground: E a = I X i =1 E i ; E b = J X j =1 E i ; where I and J are the n um b er of elemen ts in the abnormal area and the bac kground, resp ectiv ely Using the t w o a v eraged v alues as new parameters, eac h solution and the 42 PAGE 53 corresp onding ob jectiv e function can b e plotted in a 2D con tour map. Eac h p oin t in the 2D map represen ts a p ossible solution. The true solution sho wn as a blac k circle is lo cated at the global minim um of solution space (ob jectiv e function = 0). As sho wn in Figure 8 (a,c), if starting with a p oin t close to the true solution, the GNM to ok only a few steps to reac h the true solution while the CGA could tak e a couple of h undred generations. Giv en a p o or initial guess, Figure 8 (d) depicts the con v ergence path of a successful CGA sim ulation, while Figure 8 (b) sho ws a GNM sim ulation that failed to con v erge to the true solution (stuc k in a lo cal minim um). It should b e emphasized that Figure 8 is only used for illustration. The actual solution space of 900 parameters is far more complex than the 2D con tour map. The are n umerous lo cal extrema that cannot b e seen in this extremely simplied t w o parameter space. The results of 100 GNM sim ulations are plotted in Figure 9. The 100 examples are randomly pic k ed from the total of 600 sim ulations (w e had 100 sim ulations for eac h of the six initial conditions). The plot helps us visualize the inruence of initialization on the b eha vior of lo cal gradien t metho ds. It is apparen t that as the initial condition deviated a w a y from the true solution, more and more GNM sim ulations failed to con v erge (trapp ed in lo cal minim um). Based on the ab o v e comparison results, the deterministic GaussNewton metho d is used in all of the follo wing application exp erimen ts, due to the consideration of its computational eciency 43 PAGE 54 0 100 200 300 400 500 600 0 100 200 300 400 500 600 Ea (kPa)Eb (kPa) 500 300 150 80 40 20 1 10 5 0 100 200 300 400 500 600 0 100 200 300 400 500 600 Ea (kPa)Eb (kPa) 500 300 150 80 40 20 10 5 1 (a) GNM with initial guess: (b) GNM with initial guess: E a =48 kP a, E b =58 kP a E a =489 kP a, E b =469 kP a 0 100 200 300 400 500 600 0 100 200 300 400 500 600 Ea (kPa)Eb (kPa) 500 300 150 80 40 20 10 5 1 0 100 200 300 400 500 600 0 100 200 300 400 500 600 Ea (kPa)Eb (kPa) 500 300 150 80 40 20 10 5 1 (c) CGA with initial guess: (d) CGA with initial guess: E a =62 kP a, E b =50 kP a E a =467 kP a, E b =515 kP a Figure 8. Con v ergence P athes of CGA and GNM with Dieren t Starting P oin ts. The c ontours (dashe d lines) r epr esent the values of obje ctive function. The black cir cle indic ates the p osition of the true solution. Note that this is a VER Y simplie d view of the solution sp ac e. The actual solution landsc ap e in higher dimensions is much mor e c omplex and is char acterize d by numer ous lo c al minimums. 0 100 200 300 400 500 600 0 100 200 300 400 500 600 Ea (kPa)Eb (kPa) 500 300 150 80 40 20 10 5 1 Figure 9. Results of GNM with V arious Initial Conditions. The cir cles and l le d triangles r epr esent the c ases that c onver ge d and diver ge d, r esp e ctively. 44 PAGE 55 CHAPTER 3 APPLICA TION IN NONRIGID MOTION MODELING 3.1 Motion P arameters in Ph ysical Mo del Let's rst examine the motion equations and iden tify the ph ysical motion parameters that can b e reco v ered and used for motion analysis. Justications will then b e pro vided for c ho osing certain parameters based on b oth ph ysical and mathematical considerations. 3.1.1 Ph ysical Motion P arameters In motion equations (7, 8), four t yp es of quan tit y can b e iden tied: Displacemen t ( u ), material prop erties ( E ), b oundary conditions ( u ; g ) and b o dy force ( F ). Displacemen t v ector represen ts the solution of motion equations and certainly is not a candidate of motion parameter. Bo dy forces (mainly the gra vit y) can b e considered as kno wn constan ts in most applications and will not b e treated as motion parameters. By the classical elasticit y theory [120 76 ], w e can seek solutions of u throughout the mo deling domain b y sp ecifying appropriate material prop erties and b oundary conditions. Therefore, material prop erties and b oundary conditions are the p ossible c hoices of motion parameter. Direct measuremen t of actual forces (Neumann data g ) is usually not a v ailable. But displacemen ts (Diric hlet data u ) can b e estimated from image sequences and robust algorithms ha v e b een dev elop ed using either feature corresp ondence or optical ro w. With the displacemen t data, ph ysical mo delbased motion analysis can b e c haracterized as a Diric hlet b oundary v alue problem. Material prop erties are therefore c hosen as ph ysical motion parameters. Although all motion parameters can b e rev o v ered sim ultaneously a sequen tial approac h seems more practical in whic h all parameters are set to constan t, except the one to b e esti45 PAGE 56 mated. F or man y biological materials, P oisson's ratio can b e appro ximated as a constan t. F or example, soft tissues can b e view ed as nearincompressible, and a v alue in the range of 0.490 0.495 is computationally safe. Mass densit y can b e acquired from literature or direct measuring. Note that mass densit y is only needed in a transien t sim ulation, and is canceled out with the temp oral deriv ativ e in static case. The fo cus is on reco v ering and using the Y oung's mo dulus as the only ph ysical parameter to facilitate motion mo deling. 3.1.2 Ph ysical and Mathematical Justications Motion lies in the core of mec hanics and has b een studied from b oth kinematic and dynamic p ersp ectiv es. In computer vision and graphics, motion analysis is pursued along the same routes: b oth kinematic and ph ysical mo dels are used. Ho w ev er, the parameters (constrain ts) in tro duced to vision and graphic mo dels are often based on geometrical considerations and ma y not b e ph ysically realistic. the c hoice of ph ysical motion parameters is motiv ated b y b oth mathematical con v enience and ph ysical soundness. Constr aint, Constitutive Equation and Material Pr op erty: Field equations (including motion equation) deriv ed from the la w of conserv ation (mass, momen tum, energy and entrop y) ha v e t w o imp ortan t features: (1) they are lo cal, and (2) they are v alid for all t yp es of medium irresp ectiv e of their material constitutions. Mathematically this results in an underdetermined system whic h requires extra constrain ts to ensure a unique solution. Although an y form of constrain t can b e applied, almost all constitutiv e equations in con tin uum mec hanics are p ostulated based on exp erimen tal observ ations and mathematical idealization. F rom a ph ysical p oin t of view, deformation is a result of ob ject's resp onse to external forces, so it is natural to form ulate the constitutiv e equation as a relationship b et w een stress and deformation. The Ho ok e's la w for elastic solids and the Stok e's la w for viscous ruids are w ell kno wn examples. Constitutiv e equations are deriv ed from macroscopic observ ations, and ha v e supp orts from fundamen tal constitutiv e axioms as w ell as microscopic statistical dynamics. In tuitiv ely the strong in ternal constrain t (resistance to b eing deformed) of solids is related to the strong in teratomic b ounds, while the fact that gases exhibit little resis46 PAGE 57 tance to external pressure is a consequence of their w eak in termolecular b ounds. Rigorous treatmen t of the constitutiv e theory is not in the scop e of this pap er. What is emphasized is that imp osing motion constrain ts through the constitutiv e equation has a rm ro ot in fundamen tal ph ysics and is w ell justied [133 24]. The stronger the constrain t the simpler the motion equation. F or example, equation (7) deriv ed with t w o linear constrain ts (3, 6) is m uc h easier to handle than the original eld equation (1). Note that (3) is a kinematic description while (6) is a constitutiv e equation. With an ev en stronger constrain t of = 0, (7) b ecomes the motion equation of rigid solids where the div ergence of stress disapp ears (more precisely the rigidit y condition is G = 0, with G as the Green strain tensor: G = 1 2 [ r u + ( r u ) T + r u ( r u ) T ]). On the other hand, constrain t m ust b e signican tly relaxed for ruid motion. Unlik e elastic solids that can remem b er and restore the original conguration after stress is remo v ed, most ruids are memoryless with resp ect to their initial conguration and th us ruid motion is often studied with Euler descriptions (v elo cit y and strainrate tensor) rather than Lagrangian descriptions (displacemen t and Cauc h yGreen strain tensor). Because of the w eak constrain t, ruid motion equations are more dicult to handle. A daptive Smo othing of Motion Field: Another motiv ation for studying material properties as ph ysical motion parameters is that they pla y an imp ortan t role in the regularized motion syn thesis. If motion eld is view ed as an unkno wn function of space, u ( x ), x R n one of the goals of motion syn thesis is to estimate a relativ ely dense distribution of u from a sparse observ ation u Lik e man y vision problems [106 ], motion syn thesis is lik ely illp osed and m ust b e regularized: u = ar g :min fk u u k 2 + R ( u ) g ; (74) where u is the motion eld to b e syn thesized, u is the observ ed sparse and noisy motion data, is the regularization parameter, and R ( u ) is a regularization (smo othness) function. If discon tin uit y exists in the motion eld, an adaptiv e regularization sc heme is needed to a v oid o v ersmo othing. Adaptiv e smo othing has b een in tensiv ely studied in surface recon47 PAGE 58 struction [46 16 128 129 ], and the basic algorithm can b e represen ted as an optimization of the follo wing energy function: W ( f ) = k H f b k 2 + Z n w ( x ) ( D f ( x )) d x ; (75) where f ( x ) ; x R 2 is a piecewise smo oth surface, H is a linear op erator, b is the observ ation data, w ( x ) is a w eigh ting function, D f ( x ) is a deriv ativ e of f and is an energy function dened on D f ( x ). Assuming that motion eld can also b e appro ximated as a piecewise smo oth function, motion syn thesis problem (74) can b e form ulated as a minimization of the follo wing functional: W ( u ) = k u u k 2 + Z n w ( x ) ( D u ( x )) d x : (76) On the other hand, motion of a con tin uum can also b e studied with the v ariational approac h. Without losing generalit y a static case of linear elastic b o dy in the in tegral form is considered: W ( u ) = Z 1 1 ( u u ) + Z n 2 ( 1 2 ( r u + ( r u ) T ) + Z n [ 1 8 ( tr ( r u + ( r u ) T )) 2 + 1 4 ( r u + ( r u ) T ) ( r u + ( r u ) T )] + Z n F u + Z 2 F s u ; (13b) where F s is surface force, 1 and 2 are Lagrangian m ultipliers, n denotes ob ject's conguration with a b oundary split in to t w o sections: 1 for Diric hlet condition and 2 for Neumann condition. Using v ariational calculus ( W = 0), motion equations (7, 8) can b e obtained as the Euler equations of (13b). This deriv ation is also kno wn as the HellingerReissner v ariational principle in elasticit y theory [58 ]. Comparing (76) and (13b), it can b e seen that the residual norm k u u k 2 is equiv alen t to the Diric hlet condition R 1 1 ( u u ), and the smo othness norm is related to the strain 48 PAGE 59 energy the third in tegral in (13b). It is in teresting to note that material prop erties, ( x ) and ( x ), pla y a role similar to w ( x ) in (76) as the adaptiv e regularization con troller, although the strain energy is more in v olv ed than man y energy functions used in surface reconstruction. Material prop erties b eing adaptiv e regularization con troller is more ob vious in dieren tial equations. In the study of m ultiscale image description, P erona and Malik [102 ] prop osed an adaptiv e smo othing sc heme that utilizes the anisotropic diusion equation: @ f ( x ; t ) @ t = r ( c r f ) ; (77) c ( x ; t ) = R ( r f ) : (78) where f ( x ; t ) is image in tensit y x is image co ordinate, t is a scale v ariable, and c is the diusion co ecien t that con trols the degree of smo othness. Note that c itself is a function of in tensit y gradien t. F or comparison, the motion equation (7) is rewritten in terms of the Y oung's mo dulus ( E ): @ 2 u @ t 2 = r [ E (1 + )(1 2 ) ( r u ) I + E 2(1 + ) ( r u + ( r u ) T )] + F : (79) It is not dicult y to see the similarit y b et w een (77) and (79). The Ho ok e's la w (6) implies that, under the same loading condition, E is a monotonically decreasing function of deformation, whic h is exactly the desired prop ert y of co ecien t c in diusiv e adaptiv e smo othing (see Figure 4 in [102 ] and related discussions). In other w ords, a spatially v arying elastic prop ert y E ( x ) in motion equation (79) has the eect of adaptiv ely smo othing the motion eld ( u ) with a mec hanism that c haracterizes the nonlinear diusion equation (77). It should b e p oin ted out, though, unlik e the motion equation that has b een w ell established in the elasticit y theory the adoption of nonlinear diusion equation in image smo othing is motiv ated b y its mathematical con v enience, alb eit its mark ed success. F or indepth discussions on the nonlinear diusion equation with resp ect to adaptiv e smo othing, readers are refered to [5 141 ]. 49 PAGE 60 Sev eral observ ations can b e made from the ab o v e comparisons: 1. Both E and c are monotonically decreasing functions of the deriv ativ es of their corresp onding unkno wn v ariables ( r u ; r f ), a necessary condition to a v oid o v ersmo othing. 2. Similar to the surface reconstruction where small c ( x ) or w ( x ) indicate the existence of edges (for example, c ( x ) = 0 or w ( x ) = 0), in motion analysis, small E ( x ) in an elastic b o dy corresp onds to the discon tin uit y in motion eld u (see the follo wing spring example). 3. The Y oung's mo dulus can b e view ed as a regularization con troller in motion syn thesis. 4. The Y oung's mo dulus is an actual ph ysical quan tit y w elldened in elastic dynamics. In surface reconstruction or image restoration, c is an abstract co ecien t that is p ostulated for computational con v enience, and is not considered or in tended to b e used as a "true" ph ysical quan tit y ev en though c do es ha v e a ph ysical denition of conductivit y in thermo dynamics and h ydro dynamics. 5. As a consequence of b eing a ph ysical quan tit y the Y oung's mo dulus is unique for eac h ob ject and can b e used as a constan t (under the isothermal condition). On the other hand, man y dieren t c functions ha v e b een prop osed to ac hiev e adaptiv e smo othing. The comparativ e study is summarized in T able 1. Adaptiv e smo othing will b e revisited in the discussion of motion syn thesis. T o further illustrate the adaptiv e smo othing eect of E on motion eld, a tension analysis of a simple spring is presen ted. The spring is 19 cm long and consists of 3 sections of dieren t materials, with the middle section of a m uc h smaller Y oung's mo dulus (Figure 10). Numerical computation w as p erformed with a nite elemen t mo del that discretizes the spring in to 19 elemen ts (1 cm eac h) and 20 no des. The spring w as xed at one end and stretc hed b y 5 cm at the other end. The deformation of eac h no de is measured with resp ect to its prestretc hed p osition. The distributions of displacemen t u (x) and the Y oung's mo dulus E ( x ) are plotted in Figure 11. The middle section of small E causes a large "jump" 50 PAGE 61 T able 6. Adaptiv e Smo othing in Motion Analysis and Surface Reconstruction. Ph ysicsbased Motion Syn thesis Surface/image Reconstruction Unkno wn Displacemen t, u ( x ) surface or image in tensit y f ( x ) V ariational R 1 1 ( u u ) + (Energy) R n [ 1 8 E (1+ )(1 2 ) ( tr ( r u + ( r u ) T )) 2 + 1 4 E 2(1+ ) ( r u + ( r u ) T ) ( r u + ( r u ) T )] + R n 2 ( 1 2 ( r u + ( r u ) T ) + R n F u + R 2 su k H f b k 2 + R n w ( x ) ( D f ( x )) d x Smo othness Con troller E ( x ) w ( x ) PDE @ 2 u @ t 2 = F + r [ E (1+ )(1 2 ) ( r u ) I ] + r [ E 2(1+ ) ( r u + ( r u ) T )] @ f ( x ;t ) @ t = r ( c r f ) Smo othness Con troller E ( x ) / 1 r u c = g ( r f ), c / 1 r f Comparison (1) E is a measurable ph ysical quan tit y (1) c is a co ecien t without the ( E c ) corresp onding ph ysical denition in surface/image reconstruction. (2) As an in trinsic material prop ert y (2) V arious forms of c can b e E can b e used as a xed constan t. constructed as functions of the deriv ativ es of f ( x ). 51 PAGE 62 Figure 10. Congurations of Spring Before and After Stretc hing. in the motion eld, i.e. a motion discon tin uit y while t w o sections of large E ha v e relativ ely small and smo oth displacemen ts. This has the same eect of setting c ( x ) or w ( x ) to zero at edges of a surface/image. It should b e emphasized again that this adaptiv e motion smo othing phenomenon exhibited b y an elastic b o dy is determined b y its ph ysical prop ert y that cannot b e c hanged arbitrarily F or example, setting E = 0 do es not mak e ph ysical sense in the con text of solid mec hanics. Instead, adaptiv e smo othing of motion eld is determined b y the relativ e prop ert y dierence (heterogeneit y) that could b e in the range of sev eral orders of magnitude [28 ]. Based on biomec hanical tests [42 ], the Y oung's mo dulus of t ypical materials found in biological tissues are in the follo wing ranges (kP a): b one ( 10 7 ), collegian ( 10 6 ), elastin ( 10 3 ) and blo o d v essel ( 10 2 ). 3.2 Syn thesis with Reco v ered P arameter 3.2.1 Three P ersp ectiv es Motion syn thesis can b e studied from sev eral p ersp ectiv es: (1) regularized motion reconstruction, (2) Ba y esian inference, and (3) Diric hlet t yp e ph ysical mo deling. F rom a data in terp olation p oin t of view, motion syn thesis can b e view ed as a regularized reconstruction problem where a smo oth motion eld u ( x ) is estimated from sparse and noisy data u ( x ) b y minimizing function (76). V arious w ( x ) and functions can b e designed to 52 PAGE 63 (a) Distribution of the Y oung's mo dulus ( E ) (b) Motion eld discon tin uit y Figure 11. Adaptiv e Con trol of Motion Field b y the Y oung's Mo dulus. Note the sp atial c orr esp ondenc e b etwe en the se ction of smal l E and the disc ontinuity in motion eld. E is plotte d on lo garithm sc ale. Displac ement of e ach no de is c ompute d with r esp e ct to its original p osition: u = x x 0 a v oid o v ersmo othing of motion eld. Sh ulman and Aloimonos [118 ] prop osed a regularized approac h to in terp olate a smo oth motion eld with discon tin uit y preserv ed. A unique feature of their approac h is that the co ecien t w ( x ) is considered as appro ximately a constan t, although it is not a ph ysical quan tit y Motion syn thesis can also b e view ed as an inference pro cess in the Ba y esian framew ork. With motion observ ation b eing appro ximated b y a measuremen t op erator G : u = G u + ; (80) the p osterior probabilit y of u giv en the noisy data u can b e computed as: p ( u = u ) = p ( u = u ) p ( u ) p ( u ) : (81) Regularized motion syn thesis can b e treated as a sp ecial case of Ba y esian inference. Assuming that the lik eliho o d and prior v ariables in (81) are Mark o v random elds (MRF), b y the HammersleyCliord theorem [13 ], their probabilities ha v e a Gibbs distribution: p ( u = u ) = 1 Z 1 e W 1 ( u = u ) ; p ( u ) = 1 Z 2 e W 2 ( u ) ; (82) 53 PAGE 64 where Z i is a normalizing constan t kno wn as the partition function, and W i is the energy function. By treating the evidence p ( u ) as a constan t, (81) is rewritten as: p ( u = u ) = 1 p ( u ) Z 1 Z 2 exp [ W 1 ( u = u ) W 2 ( u )] : (83) So, the p osterior v ariable is also a MRF. Since exp onen tial function is monotonic, Maximum A Posterior (MAP) solution of a MRF is equiv alen t to the minimizer of its energy function: arg. max u [ p ( u = u )] = arg. min u [ W 1 ( u = u ) + W 2 ( u )] : (84) If data is kv arian t and noise is i.i.d additiv e normal, N ( = 0 ; = 2 1 I ), then: p ( ) = 1 (2 ) k = 2 j j 1 = 2 exp[ 1 2 ( ) T 1 ( )] = 1 (2 ) k = 2 k 1 exp[ 1 2 2 1 2 ] ; (85) whic h implies that the lik eliho o d densit y is prop ortional to the exp onen tial function of residual norm: p ( u = u ) / exp[ 1 2 2 1 k u G u k 2 ] ; W 1 ( u = u ) = 1 2 2 1 k u G u k 2 : (86) Similarly with a h yp othesis of motion eld as: u N ( = u m ; = 2 2 I ), the prior densit y is prop ortional to the exp onen tial function of smo othness term in the regularization approac h: p ( u ) / exp[ 1 2 2 2 k u u m k 2 ] ; W 2 ( u ) = 1 2 2 2 k u u m k 2 : (87) With W 1 and W 2 in (86) and (87), (84) b ecomes: arg. max u [ p ( u = u )] = arg. min u [ 1 2 2 1 k u G u k 2 + 1 2 2 2 k u u m k 2 ] : (88) Therefore, Ba y esian inference is reduced to the classical quadratic regularization form, as a consequence of MRF and Gaussian condition on b oth prior distribution and data noise (lik eliho o d). Readers are refered to [124 ] for a thorough co v erage on the statistical tec hniques for solving in v erse problems. 54 PAGE 65 T able 7. Three P ersp ectiv es for Motion Syn thesis. P ersp ectiv es Comp onen t 1 Comp onen t 2 Regularized Reconstruction Residual Norm Regularization Norm k u u k 2 R n w ( x ) ( D u ( x )) d x Ba y esian Inference Lik eliho o d Solution Prior p ( u = u ) = p ( u = u ) p ( u ) p ( u ) N (0 ; 2 1 I ), u N ( u m ; 2 I ), u = G u + p ( u = u ) / exp[ 1 2 2 1 k u G u k 2 ] p ( u ) / exp[ 1 2 2 2 k u u m k 2 ] Ph ysical Motion Equation Diric hlet Condition Constitutiv e Constrain t (Elasticit y Theory) b y Material Prop erties u = u ; on @ n r + F = 0, = ( r u ) I + r u + ( r u ) T In previous sections, the analogy b et w een the regularized motion reconstruction (13) and elastic dynamics in v ariational form (13b) has b een discussed, particularly the role of Y oung's mo dulus as an adaptiv e smo othness con troller of motion eld, the equiv alence of residual norm to Diric hlet condition, and the corresp ondence b et w een regularization term and strain energy Motion equations in PDE forms (7,8) can b e deriv ed from the v ariational form b y the HellingerReissner principle [58 ]. Three p ersp ectiv es are summarized in T able 2, with an observ ation of follo wing equiv alencies: (1) residual norm lik eliho o d Diric hlet condition; (2) regularization norm prior strain energy (constitutiv e constrain t). 3.2.2 Diric hletT yp e Motion Syn thesis with Reco v ered Prop ert y T o b e consisten t with our original motiv ation of analyzing nonrigid motion using ph ysical motion parameters, motion syn thesis is approac hed b y solving a Diric hlet b oundary v alue problem. This approac h has the follo wing adv an tages: (1) Diric hlet problem is w ellp osed and pro ofs can b e readily found in computational ph ysics and applied mathematics [24 73 ]; (2) giv en material prop erties and b oundary measuremen t, motion can b e syn thesized throughout the mo deling domain; (3) syn thesis algorithm can b e implemen ted with nite elemen t metho d. The coupling of trac king and syn thesis can b e realized in a mo del whose 55 PAGE 66 trac king mo dule is resp onsible for monitoring and measuring b oundary deformation while the syn thesis mo dule computes motions inside the ob ject. Reconstructing a smo oth function from noisy data with discon tin uit y preserv ed is a dicult y problem, b ecause the lo cation of discon tin uit y is unkno wn a priori In surface reconstruction, Geman and Geman [46 ] presen ted a \line pro cess" mo del b y treating b oth surface function and discon tin uit y as dual Mark o v random elds. Blak e and Zisserman [16 ] prop osed a ph ysicsbased approac h b y form ulating surface reconstruction as \w eak spring" or \w eak thin plate" mo dels. Since then, t w o general solution strategies ha v e b ee adopted: (1) reconstruct surface and discon tin uit y sim ultaneously; (2) reconstruct surface with the discon tin uities predened. The main dicult y of reconstructing surface and discon tin uit y sim ultaneously is that the energy function to b e minimized is noncon v ex. V arious deterministic and sto c hastic algorithms ha v e b een prop osed, suc h as graduated noncon v exit y [16 ], sim ulated annealing [46 ], lo cal v alidation [129 ] and others. If discon tin uit y can b e predetermined, as suggested b y T erzop oulos [129 ], discon tin uit ypreserving reconstruction problem w ould b e m uc h easier to handle n umerically In other w ords, the prior kno wledge ab out discon tin uit y w ould b e obtained from other cues, rather than b eing inferred directly from noisy data. Unfortunately in surface reconstruction, suc h prior kno wledge of discon tin uit y is rarely a v ailable, although recen t progress in Gibbs prior learning [150 ] sheds some ligh ts on the p ossibilit y of obtaining generic structural information from natural scenes. In ph ysicsbased motion syn thesis, if the Y oung's mo dulus is kno wn, the second reconstruction strategy b ecomes a feasible option, b ecause the Y oung's mo dulus is an adaptiv e smo othing con troller with its lo w v alues indicating the presence of motion discon tin uit y Homogeneous prop ert y is often assumed in motion mo deling. The v alidit y of using homogeneous prop ert y is taskdep enden t. In applications that are dominated b y global motion trac king, errors in lo cal motion caused b y the use of homogeneous prop ert y migh t b e acceptable. On the other hand, in applications where detailed lo cal deformation has signican t 56 PAGE 67 implications, prop ert y heterogeneit y m ust b e tak e in to accoun t. A simple homogeneous prop ert y could lead to large mo deling errors (see examples in Section 5). This is esp ecially true in medical imaging, where lo cal deformations could pro vide in v aluable diagnostic information. Implemen tation of the Diric hlett yp e motion syn thesis can b e realized b y discretizing the motion equation with nite elemen t metho d. Deriv ation of nite elemen t form ulation using the Ra yleighRitz or Galerkin metho ds can b e found in man y textb o oks. The algebraic equation that assem bles all elemen ts is: Ku = F g (89) where K is the stiness matrix that includes material prop erties, u is no dal displacemen t solution, and F g is the generalized force. In practice, Diric hlet condition can b e implemen ted using the p enalt y metho d. F or example, to sp ecify a displacemen t at no de i u i = y a p enalt y term can b e added to (89), k u i = k y By setting k to a v ery large n um b er with resp ect to the corresp onding diagonal elemen t in the stiness matrix, k k ii Diric hlet condition can b e satised n umerically [151 ]. Motion syn thesis is approac hed with Diric hlet data b ecause it is more accessible than Neumann data, alb eit the later is more in tuitiv e in terms of driving the deformation. If necessary force can b e deriv ed from displacemen t with appropriate mapping functions. By constructing a p oten tial function of in tensit y gradien t and using a suitable pro jection, a virtual force can b e deriv ed from image measuremen t [126 ]. T o obtain actual b oundary force, transformation similar to Diric hlettoNeumann mapping can b e used: E ; ( u ) = n f E (1 + )(1 2 ) ( tr ( u )) I + E (1 + ) ( u ) g @ n (90) where E ; ( u ) is the Diric hlettoNeumann map, n is the unit outer normal to @ n, and u is the sp ecied b oundary displacemen t. The ph ysical in terpretation of ( u ) is that the b oundary displacemen t is transformed in to b oundary traction. This deriv ation of actual force requires the kno wledge of material prop erties, at least on the b oundary 57 PAGE 68 E = 2500 kPa E = 500 kPa E = 500 kPa U 1 = 0 (fixed) U 20 = 5 cm 19 cmFigure 12. Conguration of a Syn thetic Spring. The syn thesis sc heme can b e extended to transien t cases b y solving an initialb oundary v alue problem. In motion analysis, a mo del with reco v ered ph ysical parameter has the adv an tage that the parameter can b e used as a constan t without the need of up date in the subsequen t trac king and syn thesis. Only feature that need to b e trac k ed is ob ject's b oundary whic h is usually the most salien t image cue. 3.3 Exp erimen ts Tw o exp erimen ts are presen ted to demonstrate the prop osed motion analysis metho d: (1) a 1D syn thetic spring is used to illustrate the reco v ery algorithms; (2) an image sequence of a 2D elastic ob ject is used to demonstrate ho w the reco v ered heterogeneous Y oung's mo dulus can b e utilized to syn thesize in terior motion from b oundary observ ations. 3.3.1 Reco v ery of Ph ysical Motion P arameter The adv an tage of exp erimen ting with syn thetic mo del is that its w ell con trolled nature (data error and true solution are kno wn and adjustable) allo ws us to examine all asp ects of the reco v ery algorithm with v arious test scenario. The information con tained in a simple spring exp erimen t is surprisingly ric h, whic h helps us to gain insigh ts on the b eha vior of algorithms and to optimize the implemen tation. F orwar d Mo del: The mo del conguration is as follo ws (Figure 12): (1) spring has a length of 19 cm, discretized in to 19 elemen ts of equal length; (2) elemen t is of uniaxial tensioncompression t yp e without b ending eect. Eac h elemen t has a crosssection area of 0.1 cm 2 and eac h no de has one degree of freedom; (3) material prop erties are isotropic and linear. A P oisson ratio of 0.495 is assigned to all elemen ts. Prop ert y heterogeneit y in terms of Y oung's mo dulus is as follo w: all elemen ts ha v e a v alue of 500 kP a, except the v e 58 PAGE 69 Figure 13. Tw o F orw ard Sim ulations with Dieren t but Prop ortional Elasticities. elemen ts in the middle that ha v e a higher v alue of 2500 kP a. This t yp e of conguration is of in terest in medical imaging b ecause diseased tissues often sho w up as hard inclusions surrounded b y soft tissues; (4) one end of spring is xed and the other end is pulled b y 5 cm, whic h is equiv alen t to a force of 1.6667 NT; (5) spring is in equilibrium after deformation (static sim ulation). Inverse Pr oblem: Tw o t yp es of data are considered for reco v ering spring's Y oung's mo dulus: (1) data that includes only displacemen ts, and (2) data that includes b oth displacemen ts and b oundary forces. Ideal displacemen t data is generated b y running the forw ard mo del with prop erties and b oundary conditions as sp ecied ab o v e. Noisy data is obtained b y adding 1% white noise to the ideal data. Solution Uniqueness and Data T yp es: A few examples are used to demonstrate n umerically that the uniqueness of in v erse solution is determined b y the data typ e (Diric hlet or Neumann). Tw o forw ard sim ulations are conducted using mo dels that are iden tical except that the second mo del has an elasticit y 10 times higher than the rst one. Giv en the same Diric hlet conditions (xed at no de 1 and pulled at no de 20 b y 5 cm), t w o mo dels generated the same deformation on eac h no de (Figure 13). This exercise indicates that eac h displacemen t data set is related to an innite n um b er of p ossible elasticit y distributions as long as they are prop ortional to eac h other: E i = bE j with b b eing a real constan t. Therefore, a unique elasticit y solution cannot b e obtained from displacemen t data only 59 PAGE 70 Tw o in v erse exp erimen ts are used to sho w the non uniqueness of elasticit y solution asso ciated with displacemen t data. F or eac h exp erimen t, the same ideal data w as used, u = ( u 1 ; u 2 ; :::; u 20 ), but dieren t initial guesses (800 kP a and 1600 kP a, b oth uniform). The reco v ered elasticit y is prop ortional to the true solution, but with dieren t absolute v alues. The exp erimen ts imply that the solution space of in v erse problem form ulated with only displacemen t data has man y minima of equal magnitude. The exp erimen ts also indicate that only r elative elasticit y can b e reco v ered with displacemen t data. F or applications where the primary in terest is to iden tify lo cal prop ert y abnormalities, reco v ering a relativ e elasticit y distribution is certainly an attractiv e alternativ e. With a data set that includes b oth displacemen t and force, u = ( u 1 ; u 2 ; :::; u 20 ; F 20 ), a unique elasticit y solution w as obtained regardless of the initial v alues b eing used (not to o far a w a y from the true solution). Since actual force measuremen t is rarely a v ailable in practice, surface traction can b e deriv ed using (90), pro vided that the Y oung's mo dulus on the b oundary is kno wn. F or the spring mo del, 5 cm displacemen t at no de 20 ( E = 500 kP a) is equiv alen t to a pulling force of 1.6667 Newton. Solution Stability and Data Quality: In addition to the uniqueness issue, data noise causes more serious problem of solution instabilit y In the next exp erimen t, a data set that w as corrupted b y 1% white noise w as used, u = ( u 1 ; u 2 ; :::; u 20 ; F 20 ). The results are plotted in Figure 14. The elasticit y solution obtained without regularization ( = 0) sho ws a t ypical \zigzag" pattern, indicating high instabilit y of in v erse solution. Small p erturbation errors in the data are dramatically amplied in solutions. The illp osed nature of in v erse problems is often o v erlo ok ed and solution is computed b y naiv e least square approac h, i.e. = 0. This example demonstrates that the solution obtained with a minimal residual norm is ph ysically meaningless. In other w ords, the regular least square metho d w ould \o v ert" solution to the noisy data to a degree that noise propagation completely dominates the solution prole. On the other hand, if the residual norm is p enalized to o m uc h with a large inevitably the solution will b e o v ersmo othed (Figure 14 (b)). The k ey is therefore to balance the 60 PAGE 71 (a) T rue solution (b) = 1 : 0 E 10 (c) = 1 : 0 E 14 (d) = 1 : 6 E 15 (e) = 1 : 0 E 17 (f ) = 0 Figure 14. In v erse Solutions with Dieren t Regularization P arameters. con tributions from data and prior with an optimal regularization parameter (Figure 14 (d)). I plot the solution norm k E k residual norm k A ( E ) u k and relativ e error norm k E E tr ue k = k E tr ue k against regularization parameter in Figure 15. The t ypical semicon v ergence prop ert y of regularization algorithm is clearly visible and the optimal is iden tied to b e 1.6E15. The small v alue of is mainly due to the large magnitude disparit y b et w een solution and data, whic h can b e scaled prop erly if n umerical underro w is a concern. Performanc e: The baseline and GCVbased algorithms are applied to the same in v erse problem to assess their p erformance. The mo del conguration of spring and b oundary condition are the same as those in previous exp erimen ts. T o ensure con v ergence of GaussNewton iteration in GCVbased metho d, the iteration w as started with a larger and enforcing a co oling sequence using i +1 = max ( + r i ), with b eing the curren tly estimated v alue from GCV function and r < 1 : 0. The initial solution is a uniform v alue of 800 (kP a). 61 PAGE 72 (a) V ariation of norms (b) Semicon v ergence prop ert y Figure 15. Norms against Regularization P arameter and SemiCon v ergence. (a) (b) (c) Figure 16. P erformances of Reco v ery Algorithms. (a) R e c over e d elasticity maps using b aseline and GCVb ase d algorithms. (b) Conver genc e of r e gularization p ar ameter in GCVb ase d algorithm. (c) Conver genc e of solution in GCVb ase d algorithm. E t denotes true solution. The results are sho wn in Figure 16. Both algorithms are successful in terms of the qualit y of reco v ered elasticit y with the heterogeneit y zone in the middle of spring clearly iden tied. Both algorithms con v erged within 15 iterations, but the baseline approac h to ok considerable more time b ecause of the large n um b er of nonlinear equation (18) that has to b e solv ed. The con v erged from the GCVbased algorithm is sligh tly smaller (0.47E15) than the v alue c hosen b y the baseline algorithm (1.57E15). 62 PAGE 73 3.3.2 Motion Syn thesis with Reco v ered P arameter Data A c quisition and Par ameter R e c overy: The ob ject used in the exp erimen t is an elastic bandage with a dimension of 11.6 b y 7.4 cm (Figure 17). The bandage is designed for w ound care and made of material that has a Y oung's mo dulus close to that of articial skin (p olytetraruoro eth ylene). T o create prop ert y heterogeneit y less elastic tap es w ere attac hed at b oth ends of bandage. Measured Y oung's mo dulus of bandage and tap e from tensile test are 160 and 7900 (kP a), resp ectiv ely As a result, the nal elastic ob ject has three sections: an easytodeform section in the middle and t w o lessdeformable sections on b oth ends. A 7x10 rectangle grid w as mark ed on the bandage to facilitate corresp ondence matc hing. An image sequence w as acquired with Minolta range camera while the bandage w as gradually stretc hed laterally F rame 6 and frame 7 are the same except that, in frame 7, the cen tral p ortion of bandage w as temp orally blo c k ed b y another ob ject (Figure 17 (d)). F rame 7 will b e used to demonstrate the b oundarybased motion syn thesis sc heme in the presence of o cclusion. A nite elemen t mo del w as built directly on top of the grid, with 70 no des and 54 quadrilateral elemen ts (Figure 18 (a)). F rame 1 and frame 4 w ere used to reco v er the Y oung's mo dulus. Displacemen ts w ere computed with the aid of grid corresp ondence and assigned to eac h no de of mo del. T o reco v er the Y oung's mo dulus with a higher resolution, the mesh w as rened through bilinear in terp olation (Figure 18 (b)). The Y oung's mo dulus on b oundary elemen ts w ere kno wn and used to deriv e b oundary stress b y (90). The results are plotted in Figure 18 (d). The o v erall distribution of reco v ered elasticit y matc hes v ery w ell with measured elasticit y with the cen tral heterogeneit y section clearly iden tied. Signic anc e of R e c over e d Heter o gene ous Pr op erty: Motions in the in terior of a deformable ob ject can b e syn thesized from b oundary observ ations b y solving a Diric hlett yp e forw ard problem, pro vided that material prop erties are kno wn. In this case, the use of homogeneous prop ert y is not v alid, b ecause it automatically leads to a simple smo oth motion eld. 63 PAGE 74 (a) F rame 1 (b) F rame 4 (c) F rame 6 (d) F rame 7 with o cclusion Figure 17. Image Sequence of a Bandage for Motion Syn thesis. (a) Original nite elemen t mesh (b) Rened nite elemen t mesh 0 5 10 15 20 0 5 10 15 20 25 30 0 2 4 6 8 10 x 10 6 0 5 10 15 20 0 5 10 15 20 25 30 0 2 4 6 8 10 x 10 6 (c) Measured elasticit y distribution (d) Reco v ered elasticit y distribution Figure 18. Finite Elemen t Mo del and the Reco v ered Elasticit y 64 PAGE 75 T able 8. Syn thesized Motion along CrossSection Line. No de Measured Syn thesized Motion Error Syn thesized Motion Error ID motion (m) (homog. elasticit y) (%) (reco v ered elasticit y) (%) 31 0.0121 0.0121 0 0.0121 0 32 0.0118 0.0089 25 0.0120 2 33 0.0108 0.0058 46 0.0123 14 34 0.0099 0.0021 79 0.0080 19 35 0.0009 0.0010 211 0.0012 33 36 0.0056 0.0042 25 0.0065 16 37 0.0132 0.0076 42 0.0155 17 38 0.0189 0.0109 42 0.0181 4 39 0.0191 0.0146 24 0.0180 6 40 0.0186 0.0186 0 0.0186 0 Finite elemen t mo del is used to syn thesize in terior motion of bandage in frame 7 that is partially blo c k ed. Displacemen ts on bandage b oundary w ere measured b et w een frame 1 and frame 7 and then assigned to the b oundary no des of mo del. F or comparison, t w o exp erimen ts w ere p erformed, one with an assumed homogeneous Y oung's mo dulus of 7900 (kP a) and another one with the reco v ered Y oung's mo dulus (Figure 18 (d)). Syn thesized motion at the cen ter of eac h elemen t w as plotted as small squares in frame 6 (Figure 19 (a,b)). Note that F rame 6 is the same as frame 7 except the o cclusion. There exist signican t mismatc hes b et w een true p ositions and syn thesized p ositions using homogeneous Y oung's mo dulus (Figure 19 (a)), esp ecially in the middle section. Mo del with reco v ered Y oung's mo dulus reduces syn thesis error signican tly (Figure 19 (b)). T o visualize the adaptiv e con trol of Y oung's mo dulus on motion eld, a crosssection is dra wn in the cen ter of bandage along the direction of stretc hing for b oth displacemen t and elasticit y distributions (Figure 19 (c,d,e,f )). As exp ected, homogenous elasticit y causes an unrealistic smo oth transition in motion eld while reco v ered heterogeneous elasticit y helps preserv e the observ ed motion discon tin uit y The results are also listed in T able 3. This example also demonstrates the adv an tage of the prop osed syn thesis sc heme o v er other approac hes that require trac king all features explicitly whic h could b e defeated b y the existence of o cclusion. 65 PAGE 76 This motion syn thesis sc heme using the reco v ered ob jectsp ecic heterogeneous prop ert y holds great p oten tials in medical imaging. Ba jcsy et al [11 ] presen ted an elastic registration metho d based on a simplied Na vier's equation that w as widely adopted in medical image registration and surgery planning [99 54 ]. Metho d based on b oundary mapping w as also prop osed for nonrigid registration [25 ]. Man y researc hers recognize that it is essen tial to incorp orate anatomical structure with spatially v ariable prop erties in to biomec hanical mo dels. This syn thesis sc heme can also b e used in realistic visualization and animation, esp ecially those use ph ysicsbased mo deling tec hniques [74 ]. Vivid sim ulation of liv e ob jects that in v olv es subtle deformation of soft tissues, suc h as facial expression, is only p ossible if detailed kno wledge of tissue prop erties is incorp orated in to the mo del. 66 PAGE 77 (a) Syn thesized motion (b) Syn thesized motion (homogeneous elasticit y) (reco v ered heterogeneous elasticit y) (c) Crosssection of elasticit y (d) Crosssection of elasticit y (homogeneous) (reco v ered heterogeneous) (e) Crosssection of syn thesized motion (f ) Crosssection of syn thesized motion (homogeneous elasticit y) (reco v ered heterogeneous elasticit y) Figure 19. Comparison of Tw o Exp erimen ts. Synthesize d motion using the assume d homo gene ous elasticity (a,c,e) and the r e c over e d heter o gene ous elasticity (b,d,f ). Me asur e d displac ement is shown as solid line in (e,f ). Synthesize d displac ements ar e depicte d as cir cles and triangles in (e,f ). 67 PAGE 78 CHAPTER 4 APPLICA TION IN BURN SCAR ASSESSMENT Eac h y ear more than one million p eople suer burn injuries in Canada and the US [111 ]. Accurate rating of scar condition is needed in order to design an eectiv e treatmen t plan. Scar rating in clinical settings is done using the V ancouv er rating scale or its v arian ts b y whic h exp erts assess the v ascularit y pigmen tation, pliabilit y and size of scars [121 145 ]. These rating metho ds suer from their sub jectiv e nature and lo w consistency among rates. There is a strong need for a quan titativ e and ob jectiv e scar rating metho d based on the bioph ysical prop erties of skin tissue [107 ]. In this application, a ph ysical mo delbased scar rating metho d is presen ted whic h fo cuses on estimating the relativ e elasticit y of scars. The con tributions of the prop osed metho d are: (1) nonin v asiv e natural image features are used to generate an adaptiv e mesh and to build a nite elemen t mo del. In previous w orks, articial mark ers (ink stamps) on the skin to facilitate mo del construction had b een used. The use of natural features enables us to quan tify the elastic prop ert y without the need of an y mark ers, and hence greatly enhances the applicabilit y of the prop osed rating metho d; (2) a robust pro cedure of quan tifying the scar elasticit y is established using the regularization metho d so that the noisy data can b e handled prop erly Tw o imp ortan t issues in scar elasticit y estimation are the illp osedness and the computational cost. The stabilit y of in v erse solution can b e partially restored through the regularization metho ds. The follo wing approac hes are exp erimen ted in the con text of burn scar assessmen t: (1) reduce the parameter space b y p osing stronger constrain ts; (2) reduce the parameter space b y using an adaptiv e mesh. 68 PAGE 79 4.1 Mo del Construction The mec hanical b eha vior of skin is determined b y the presence of collagen b ers, elastin b ers and lubricating ground substance [77 42 130 ]. Burn scars tend to ha v e random organization of collagen b ers during tissue regeneration. In the previous studies of burn scar ev aluation, it has b een sho wn that the use of elastic and isotropic mo del is adequate for this particular task [135 107 ]. 4.1.1 Extraction of Natural P oin t F eatures The ShiT omasi detector is used to extract salien t p oin ts in scar images [117 114 ]. F or eac h pixel p in the image, a 2 2 matrix ( g ) is computed within a windo w ( w ): g = 2 6 4 P w @ I @ x @ I @ x P w @ I @ x @ I @ y P w @ I @ y @ I @ x P w @ I @ y @ I @ y 3 7 5 (91) where I is the image in tensit y and ( x; y ) represen t ro w and column. The rst deriv ativ es is computed b y con v olving the in tensit y with the deriv ativ e of Gaussian lter ( G ): @ I @ x = @ G @ x I @ I @ y = @ G @ y I The co ecien ts of gradien t matrix is then computed: @ I @ x @ I @ x ; @ I @ y @ I @ y ; @ I @ x @ I @ y ; @ I @ y @ I @ x The co ecien ts of all pixels inside w are summed up to pro duce the matrix g of pixel p The eigen v alues of matrix g is then computed: ( 1 ; 2). Giv en a threshold T satisfaction of the condition: min( 1 ; 2) > T indicates that the windo w con tains a corner/p oin t with t w o strong edges along the eigen v ector directions. A v alue of T = 1 is used for all of the scar exp erimen ts. With this v alue, 400800 features can b e ensured in the scar images, whic h suces the need of building a relativ ely dense nite elemen t mesh. Another threshold h is used, the minim um distance b et w een t w o adjacen t p oin ts, to con trol feature distribution. Since scars ha v e lo w er in tensit y than normal skins, the v alue of h is c hanged adaptiv ely based on the in tensit y v ariation. An area is rst selected that con tains only scars and skins. The selected area is then equalized to highligh t the in tensit y con trast b et w een scars and skins. F or eac h feature p oin t, its h is computed based on the a v erage in tensit y in w scaled linearly b y the o v erall in tensit y v ariation in the mo deling 69 PAGE 80 (a) F eature in frame 1 (b) Corresp ondence in frame 2 Figure 20. P osition Shift Bet w een the Corresp onding F eatures. Note that fe atur es in (a) and (b) ar e c orr esp onding p airs. The white squar e indic ates the c ompute d fe atur e p osition by the fe atur e dete ctor. The shift b etwe en (a) and (b) is 3 pixels. area: h = H min + ( I w I min ) H max H min I max I min ; where h is the minim um distance b et w een a feature and its neigh b ors, I w is the a v erage in tensit y in w ( I max I min ) are the maxim um and minim um in tensities in mo deling area after equalization, and ( H max H min ) are usersp ecied maxim um and minim um h that corresp ond to ( I max I min ). After computing h for eac h feature, all features are sorted based on the n um b er of their neigh b ors that are in conrict with their h v alue. By "conrict" w e mean that the distance b et w een a feature and its neigh b ors is less than its h The feature is then iterativ ely remo v ed that has the most conricts un til no conrict exists. As a result, the nal distribution has more features in scars than in skins. All images ha v e the resolution of 640 480 pixels. F eatures are extracted using a windo w size of 9 9. The minim um distance ( H min ) is in the range of 1520 pixels, and H max is set to b e 1.3 2.5 times larger than H min The c hoice of H min dep ends on the size of scars. If H min is to o large, there ma y not b e enough features in scar areas and vise v ersa. The same argumen t applies to the c hoice of H max Those heuristic parameters m ust b e tuned for a sp ecic setting dep ending up on the qualit y of images and size of scars. T o quan tify the precision of feature extraction, 218 corresp onding features w ere examined in scar images. F or eac h patien t, t w o frames w ere tak en (b efore and after deformation). F eatures w ere then extracted from t w o frames. There exist small shifts b et w een the com70 PAGE 81 T able 9. The P erformance of P oin t F eature Detector. Min. Max. Av erage feature p osition shift (pixel) 0 4.5 2.1 actual displacemen t (pixel) 11.9 36.2 24.8 error (shift/displacemen t) 0.0% 13.3% 7.2% puted p ositions of corresp onding features in t w o frames. An example is sho wn in Figure 20. The white square indicates the computed p osition b y the feature detector. The corresp ondence pair has a shift of 3 pixels. Note that the feature patterns c hanged sligh tly b et w een t w o frames due to the v ariations in ligh ting, pro jection and skin deformation. The actual displacemen ts due to skin deformation b et w een t w o corresp ondences are 26 pixels. Since our in terest is to obtain go o d displacemen t data, the extraction error is dened as the ratio of the computed p osition shift to the actual displacemen t. F or the example sho wn in Figure 20, the extraction errors is 11.5%. F or 218 feature pairs examined, the results are summarized in T able 9. With curren t imaging setting and feature extraction metho d, the a v erage error in tro duced in feature extraction is less than 10%. 4.1.2 Adaptiv e Meshing Giv en a set of p oin ts distributed on ob ject's surface, an adaptiv e mesh based on the Delauna y principle [19 ] can b e generated. The follo wing meshing pro cedure is used: (1) extract p oin t features from images; (2) select the p oin ts inside the region of in terest (R OI) for whic h a ph ysical mo del will b e built; (3) link the p oin ts on the b oundary of R OI to form a closed p olygon or surface; (4) generate an adaptiv e Delauna y mesh using feature p oin ts as its no des; (5) rene the mesh. Renemen t is p erformed to insert new no des in to the area of in terest as a qualit y assurance pro cedure. These guidelines are follo w ed in the no de insertion: (1) new no des are added at the region of high curv ature to ensure accurate surface represen tation; (2) new no des should sub divide the thin elemen t in to regular elemen ts to reduce mo deling errors; (3) new no des are added at regions of prop ert y discon tin uities whic h could cause mo deling errors. Figure 21 sho ws an example of adaptiv e triangle mesh. 71 PAGE 82 (a) (b) (c) Figure 21. Adaptiv e Meshing Using P oin t F eatures. (a) Sc ar image. (b) Point fe atur es and R OI. (c) A daptive triangle mesh with r enement at the b oundary. 4.2 Estmate Scar Elasticit y 4.2.1 Range Scanner Setting A K2T range scanner is used to obtain 3D displacemen t b et w een corresp onding features. The setting of K2T scanner is illustrated in Figure 22 (a). The K2T system consists of a CCD camera and a structured ligh t pro jector, and computes the depth from images of strip ed ligh t patterns. An example of structured ligh t patterns is sho wn in Figure 22 (b). In scar study the eectiv e imaging area is ab out 3035 cm cub e. The distance b et w een patien ts and the range camera is ab out 100130 cm. A set of images (24 frames) w as tak en while patien t's skin w as stretc hed. Since the scar condition of most patien ts did not allo w the use of con tact devices, patien t's skin w as pulled gen tly b y hands to a v oid pain and further damage. Therefore, it w as not p ossible to measure the forces applied to patien t's skin. 3D displacemen t ( u ) obtained from in tensit y and range images is used to compute the ob jectiv e function in an \outputleastsquared" form: k F ( E ) u k 2 where F ( E ) denotes the forw ard mo del. Our in terest is to determine the r elative elasticit y of scars. Go o d information can b e obtained ab out the geometry of scars from image cues (in tensit y color and texture). T o reduce the computational complexit y sev eral simplications are made: (1) the geometry of scars is kno wn; (2) the elasticit y of scar is higher than that of normal skin; (3) the elasticit y 72 PAGE 83 100130 cm100 cmLight Projector CCD Object (a) K2T imaging geometry (b) Ligh t strip e pattern Figure 22. The Geometry of K2T Range Scanner for Burn Scar Study of bac kground normal skin is kno wn. The minimization of Tikhono v functional is then transformed in to a onedimensional searc h problem with t w o steps: 1. Determine the regularization parameter using the Lcurv e metho d [55 ]. 2. Change scar elasticit y ( E s ) un til the minimizer of the simplied functional is found: T ( E s ) = k F ( E s ) u k 2 + k W ( E s E s ) k 2 : (92) In Lcurv e metho d, the log of the regularized solution norm l og jj E jj is plotted against the log of the residual norm l og jj F ( E ) u jj for a range of The optimal is c hosen to b e the one that corresp onds to the corner of the Lcurv e. The in tuitiv e in terpretation of the Lcurv e metho d is that the solution norm will go to the maxim um and the residual norm will go to the minim um as approac hes zero, and vice v ersa. 4.2.2 Boundary Condition Sp ecication Displacemen ts are sp ecied at the b oundary no des that are also feature p oin ts. The feature p oin ts inside the mo deling domain are designated as the con trolling no des, on whic h the measured displacemen ts and the sim ulated displacemen ts b y forw ard mo del are used to calculate the residual norm. Since only displacemen t is a v ailable, the elasticit y of normal skins is used as a reference to determine the relativ e elasticit y of scars. The rep orted Y oung's 73 PAGE 84 T able 10. The Estimated Relativ e Scar Elasticit y (kP a). P atien ts Scar elasticit y Scar elasticit y Absolute Relativ e (Natural features) (Arti. mark ers) error error 970407 53 46 7 13.2% 970416 12 14 2 16.6% 970425 41 37 4 9.7% 970922 19 18 1 5.3% Err ors ar e c ompute d b etwe en the r esults using natur al fe atur es and articial markers. The elasticity of normal skin use d in exp eriments is 5 kPa. mo dulus v alues of soft tissues are in the range of 1100 kP a [42 130 ]. In scar assessmen t, a Y oung's mo dulus of 5 kP a is used for normal skin. It is w orth noting that the b o dy force should not b e confused with the b oundary conditions in normal sense suc h as surface force/traction. The v ariation of b o dy force (mainly gra vit y) is a signican t factor in the sim ulation of largesized ob jects suc h as the sea w ater in a ba y or the planet earth. But for small ob jects suc h as scars, b o dy force can b e treated as a constan t. 4.3 Exp erimen tal Results 4.3.1 Data Set The data set includes images of four patien ts that w ere tak en at dieren t healing stages. The b oundaries b et w een scars and normal skins can b e clearly iden tied in those images. 4.3.2 Quan tify Scar Elasticit y Tw o indep enden t exp erimen ts are carried out, one using articial mark ers and one using natural features, to determine if the natural feature is a viable option for replacing the articial mark er. The results for all patien ts are listed in T able 10. Although direct measuremen ts of scar's Y oung's mo dulus are not a v ailable for those patien ts, the fact that the dierences b et w een the estimated Y oung's mo dulus using articial mark ers and using natural features are v ery small suggests a high lev el of consistency The results of patien t74 PAGE 85 970922 are plotted in Figure 23. As exp ected, exp erimen t using articial mark ers sho ws a sligh tly b etter p erformance than that using natural features, as indicated b y the minimization curv es (Figure 23 (f )). But the disparit y b et w een t w o curv es is v ery small. More imp ortan tly the minim um p oin ts of t w o curv es are almost iden tical (18 kP a and 19 kP a). The strain distributions also matc h w ell with scar distribution. The errors in measured displacemen ts mainly came from t w o sources: (1) feature extraction caused b y illumination and texture v ariations; (2) corresp ondence mismatc h. Since the corresp ondence w as man ually established, the disparit y b et w een the mo deling results of using natural features and using articial mark ers is mainly caused b y the feature extraction error. As analyzed b efore, the feature extraction error is less than 10% on a v erage, whic h is an acceptable v alue for burn scar assessmen t. 4.3.3 Relativ e Elastic Index Curren t clinical scar assessmen t metho ds are based on a relativ e rating scale. This study will in v estigate if there exists a p ositiv e correlation b et w een the clinical rating and the estimated elasticit y ratio of scars to normal skins, so that a standard can b e established that is quan titativ e and comparable to the one used b y ph ysicians. The Relativ e Elastic Index (REI) is dened as the ratio of a v erage Y oung's mo dulus of scars to that of normal skins: R E I = P n i =1 E i =n P m j =1 E j =m ; (93) where n and m are the n um b ers of elemen ts inside scars and normal skins, resp ectiv ely T able 11 lists ph ysician's ratings, as w ell as REIs using b oth articial mark ers and natural features. In all exp erimen ts, patien t's skin w as stretc hed along t w o directions (horizon tal and v ertical). So, eac h patien t w as studied t wice. In Figure 24, the REIs using natural features and articial mark ers are plotted against the ph ysician's ratings. Tw o imp ortan t observ ations can b e made: (1) there exists a go o d correlation b et w een the ph ysician's ratings and the REIs using either natural features or articial mark ers; (2) the REIs computed using natural feature are consisten t with those using articial mark ers, esp ecially for less 75 PAGE 86 (a) (d) (b) (e) (c) (f ) Figure 23. Estimated Elasticit y and Strain. (a) Sc ar image. (b) After b eing str etche d. (c) R ange image. (d) Str ain using natur al fe atur es. (e) Str ain using marker. (f ) Minimization curves. In (d) and (e), dark c olor indic ates low str ain value of sc ars. 76 PAGE 87 T able 11. Ph ysician's Rating and Relativ e Elastic Index (REI). Ph ysician's REI using REI using rating natural features articial mark ers scar970407(hori) 4.5 9.5 8.0 scar970407(v ert) 4.5 9.1 7.3 scar970416(hori) 3.0 2.8 3.3 scar970416(v ert) 3.0 4.1 3.8 scar970425(hori) 3.3 6.3 5.0 scar970425(v ert) 3.3 7.4 5.6 scar970922(hori) 2.0 2.2 2.3 scar970922(v ert) 2.0 2.8 2.0 normalskin 1.0 1.0 1.0 F or e ach p atient, the skin was str etche d in two dir e ctions, horizontal ly and vertic al ly. The REI of normal skin (1.0) is dene d as the b aseline. Figure 24. Correlation Bet w een Ph ysician's Rating and REIs. damaged scars. The sligh t deviation from the linear monotonic relationship is probably caused b y the complex scar patterns, whic h aect b oth feature extraction and matc hing. This problem can b e solv ed b y calibrating the mo del against the ground truth, a p ossibilit y that is curren tly under in v estigation. It should b e p oin ted out that, due to the limited n um b er of patien ts, no statistical conclusion w as dra wn at this p oin t ab out the correlation b et w een the ph ysicians' rating and the REI. 77 PAGE 88 CHAPTER 5 APPLICA TION IN F A CE RECOGNITION During the past a couple of y ears, biometrics researc h has receiv ed considerable atten tion due to its high p oten tial in securit y related applications. There exists a wide v ariet y of biometrics tec hniques, some are relativ ely mature while others are still in their infancy Eac h biometrics tec hnique has its pros and cons, and it is not p ossible to nd a single one that can solv e all practical problems. Therefore, there is alw a ys a need for new biometrics. Other than ngerprin t, face recognition is probably the most natural (and hence a p opular) biometrics b ecause w e ha v e dev elop ed the abilit y to recognize faces automatically with no conscious eort. Curren t face recognition metho ds rely on visible photometric or geometric attributes that are presen t in in tensit y images. Based on large amoun t of researc h and b enc hmark studies [20 40 103 ], it has b een recognized that those metho ds suer from problems asso ciated with follo wing factors: (1) illumination and p ose v ariation; (2) mak eup, hairs and glasses; (3) plastic surgery; (4) face deformation during expression (dynamic face analysis in video sequence). F uture face recognition metho ds m ust address those dicult issues. A new class of features (or biometrics) is prop osed that is deriv ed from the computed strain pattern exhibited during face expression. The metho d has sev eral adv an tages: 1. Elastic strain pattern is directly related to the material prop ert y of underlying facial m uscles. Our h yp othesis is that, if the anatomical structure of an individual face (geometry distribution and strength of b ones and m uscles) is unique, then this anatomical uniqueness should b e rerected in the elastic strain pattern. Strain pattern could b e less sensitiv e to the illumination and p ose c hanges as w ell as camourage using mak eup. 78 PAGE 89 2. The computation of elastic strain map requires at least t w o frames that capture the face deformation during expression. Most face recognition metho ds use static images only and dynamic face expression has b een considered as an adv erse factor that ma y cause p erformance degradation [142 ]. Ho w ev er, a recen t study b y Y aco ob and Da vis [143 ] indicates that deformed face (smiling face) is more recognizable and actually increases the iden tication rate. This study will go one step further b ey ond the visible cues in face expression to reco v er the elastic strain pattern that migh t help rev eal the underlying anatomical individualit y 3. Anatom ybased ph ysical mo del has b een widely used in realistic facial animation and surgery sim ulation [97 127 74 ]. Essa and P en tland [116 ] prop osed a ph ysical mo del to represen t facial motion and distinguish face expressions. They used the mo del to estimate visual m uscle activ ation and to generate motionenergy templates for eac h expression. Ho w ev er, w e are not a w are of an y study that uses ph ysical mo del to compute elastic strain pattern for the purp ose of face recognition. A nite elemen t mo del based on con tin uum mec hanics enables us to compute the strain map accurately 5.1 Hyp othesis: F ace Anatom y and Muscle Biomec hanics Biometrics refers to the iden tication of an individual based on distinctiv e ph ysiological or b eha vioral c haracteristics. In face recognition, the h yp othesis is that eac h individual has a unique visible face pattern in terms of shap e, color or texture that can b e utilized for recognition. This uniqueness in visible pattern is certainly related to, and probably determined b y the uniqueness of the underlying anatomical structure and biomec hanical comp ositions, whic h could also b e useful for recognition, if they are measurable. Ma jor anatomical units of h uman face are: b ones (skull), m uscles, skin, blo o d v essels and nerv es [39 ]. The use of skull measuremen t in iden tication (craniofacial analysis) and its forensic implications is w ell do cumen ted [64 ]. But it is doubtful that a practical biometrics can b e deriv ed due to the limitation of sp ecial imaging mo dalit y (Xra y is needed for skull measuremen t). It is also dicult to capture facial nerv e patterns with curren t 79 PAGE 90 imaging tec hnologies. The blo o d v essel patterns, ho w ev er, has b een utilized in b oth facial thermograph y and irisretina scans with the aid of infrared camera. F ace expression is con trolled b y m uscles and emotional states that trigger and c hange the m uscle mo v emen ts, whic h could b e quan tied b y the elastic strain pattern. Elastic strain pattern computed from the observ ed face deformation not only rerects an individual's emotional state, but, more imp ortan tly also rev eals the in trinsic m uscle prop erties, and therefore can b e utilized for recognition. This unique strain pattern asso ciated with an individual can remain unc hanged for a long p erio d of time, although it is reasonable to exp ect some v ariations caused b y aging (loss of elastin b ers and m uscle elasticit y), injuries and plastic op erations. An ideal face mo del to compute strain pattern should incorp orate all anatomical details. Ho w ev er, suc h a fullscale mo del w ould b e to o exp ensiv e to b e realistic for face recognition. T o strik e a balance b et w een mo deling accuracy and computational eciency w e ha v e to mak e a compromise on what to mo del and ho w to mo del. Although a whole face con tains more information ab out an individual, it is more practical to mo del a p ortion of face whose deformation is dominated b y one or a few ma jor m uscles. A section is c hosen that is b et w een the c heek b one and ja w line (side view) and co v ered b y the masseter m uscle. The masseter m uscle is a large, thic k and roughly rectangular plate that is resp onsible for ja w action. Its deformation will b e mo deled b et w een t w o p ositions, namely the closed mouth and the op en mouth. Masseter m uscle is of striated t yp e. Due to the lo cation just b eneath the skin tissue, its con traction has an immediate eect on skin motion. Therefore, w e do not treat skin as a separate functional la y er. Instead, w e mo del the deformation of m uscle and skin together as an in tegrated mec hanical unit. This singlela y er mo del satises requiremen ts of face recognition. But t w ola y er mo del ma y b e considered, esp ecially for older p eople where the motion of m uscle and aged skin are less sync hronized (wrinkles). 80 PAGE 91 5.2 Computational Metho d There are four ma jor steps in computing elastic strain pattern: (1) feature extraction and motion measuremen t; (2) nite elemen t construction; (3) strain computation; (4) con v ersion of strain map to in tensit y image. F e atur e Extr action and Motion Me asur ement: First, salien t p oin ts will b e extracted from in tensit y images and then construct a p olygonal surface using the Delauna y principle. F eature extraction is p erformed using the ShiT omasi detector [117 ]. Due to the nonrigid nature of face deformation, w e presen tly establish corresp ondence b et w een t w o frames man ually to ensure the qualit y of displacemen t v ector. The displacemen t data will b e used in the nite elemen t mo del to sp ecify the Diric hlet b oundary condition. This man ual pro cessing is not new in early face recognition researc h where ey es and other facial features are man ually lo cated (FERET test [103 ]). F uture w ork will use an automated corresp ondence matc hing metho d. Delaunay Meshing and Mo del Construction: The Delauna y triangulation is used to generate an adaptiv e triangle mesh with a set of p oin ts that are randomly distributed on face images. A simple meshing strategy is designed: (1) select p oin ts that denes the b oundary of the region to b e mo deled; (2) link those b oundary p oin ts to form a p olygon; (3) select the p oin ts that are inside the p olygon; (4) generate a triangle mesh that is adaptiv e to all the selected p oin ts. T o ensure the mesh qualit y a lo cal mesh renemen t pro cedure is used that detects badshap ed elemen ts and impro v es the initial mesh accordingly: (1) no de will b e added at the region of high curv ature to increase the accuracy of surface represen tation; (2) new no de will sub divide the long and thin elemen t in to more regularshap ed elemen t. F orwar d Mo deling and Str ain Computation: The deformation of face tissues can b e describ ed b y a motion equation: r [ ( r u ) I + G r u + G ( r u ) T ] + f b = @ 2 u @ t 2 ; (94) 81 PAGE 92 where G and are the Lam e constan ts. Equation (94) is solv ed n umerically after b eing discretized o v er the Delauna y triangle mesh using the nite elemen t metho d in v ariational form ulation [151 ]. Since the primary in terest is in static deformation only the nal forw ard mo del in the discrete matrix form b ecomes: Ku = F (95) where K is stiness matrix and F is the generalized force. T o compute the strain pattern of a deformed face, w e ha v e to supplemen t appropriate b oundary conditions. The Diric hlet condition is sp ecied using the measured displacemen t data on the feature p oin ts. As a result, it b ecomes an o v ersp ecied b oundary v alue problem of the rst kind and can b e readily solv ed using an iterativ e n umerical solv er. Str ain Conversion and PCA A nalysis: The standard principle comp onen t analysis (PCA) [137 ] is used for similarit y computation. The strain maps is normalized in to in tensit y images using a simple linear transformation: e x e min e max e min = I x I min I max I min (96) where ( e max ; e min ) are the maxim um and minim um strain v alues for all sub jects, e x is the strain v alue to b e con v erted, I x is the con v erted in tensit y v alue. ( I max ; I min ) are set to 255 and 0, resp ectiv ely Because strain v alue spans o v er a wide range, information ma y b e lost with this linear con v ersion. More sophisticated con v ersion metho d will b e considered in future in v estigations. Before PCA analysis, all con v erted strain in tensit y images are scaled based on the geometry of facial landmarks (nose, ear, mouth and ja w line). A rectangle region of the size of 150 b y 160 pixels (original face images are 640 b y 480 pixels) in the cen ter of strain in tensit y image w as then c hopp ed out. This rectangle strain image w as then used in PCA analysis. An example strain pattern is sho wn in Figure 25. The implemen tation of PCA algorithm used in this study can b e found in [14 ]. 82 PAGE 93 (a) Close mouth (b) Op en mouth (a) Original strain map (b) Con v erted strain in tensit y Figure 25. Strain P attern of a Mo died F ace. The fac e is attache d with a less str etchable tap e, which c orr esp onds to smal l str ain (low intensity). (a) (b) (c) (d) Figure 26. F our Images Acquired for Eac h Sub ject. (a) close mouth + bright light, (b) op en mouth + bright light, (c) close mouth + low light, (d) op en mouth + low light. 83 PAGE 94 T able 12. F ace Exp erimen ts. Gallery Prob e Exp. 1 27 (RF,BL) 27 (RF,LL) Exp. 2 20 (RF,BL) + 7 (MF,BL) 20 (RF,LL) + 7 (MF,LL) RF: Regular F ace. MF: Mo died F ace. BL: Brigh t Ligh t. LL: Lo w Ligh t. 5.3 Exp erimen tal Results Tw o sideview images (in tensit y + range) w ere acquired for eac h sub ject (op enmouth, closemouth) with a Minolta Vivid900 range scanner. Eac h sub ject w as imaged under t w o illumination conditions: brigh t and lo w ligh t conditions. As a result, eac h sub ject has 4 images: op enmouth + brigh tligh t, closemouth + brigh tligh t, op enmouth + lo wligh t and closemouth + lo wligh t (see Figure 26). The complete data set con tains 27 sub jects (108 images). T o study the ecacy of the prop osed metho d in the presence of mo died faces (tissue prop erties are c hanged due to surgery trauma or burn, either in ten tionally or acciden tally), a transparen t rectangular tap e w as attac hed on the face of 7 sub jects. The tap e is less stretc hable and th us has the same eect of mo died faces and results in abnormal strain patterns. Tw o exp erimen ts w ere p erformed with dieren t gallery and prob e comp ositions (see T able 1). The results are presen ted with the standard receiv er op erating c haracteristic (R OC) curv es that is used to quan tify iden tication p erformance. The normalization pro cedure as rep orted in FR VT2002 [40 ] w as used in the R OC computation. The purp ose of Exp erimen t1 is to in v estigate whether elastic strain pattern has enough discrimination p o w er for recognition of regular faces. Exp erimen t2 is to understand whether the mo died faces are more dicult to recognize or more distinguishable due to their unusual strain pattern. The R OC curv es of t w o exp erimen ts are sho wn in Figure 27. F or the gallery and prob e sets that include only regular faces, a v erication rate of 67.4% w as observ ed at a false alarm rate of 5%. Although the results w ere obtained with a relativ ely small data set, the 84 PAGE 95 Figure 27. R OCs for Regular and Mo died F aces. p erformance is still quite promising. The data sets do not con tain fron tal views that are commonly used in face recognition test. As exp ected, the exp erimen t with gallery and prob e sets that con tain mo died faces sho ws a b etter p erformance. A t a false alarm rate of 5%, w e w ere able to ac hiev e a v erication rate of 78.2%. 6 of 7 mo died faces w ere correctly iden tied in Exp erimen t2 while only 3 of 7 regular faces of the same sub jects w ere iden tied in Exp erimen t1. The t w o exp erimen ts suggests that a p erson who c hanged his/her app earance b y plastic surgery or other approac hes actually has a b etter c hance to b e detected b ecause medical surgery cause prop ert y c hanges of facial tissues, whic h is hard to detect using the metho ds that rely on visible cues only The rationale for the new biometrics is that the underlying anatomical structure is unique for eac h individual, and this ph ysiological in v arian t can b e explored for face recognition. The prop osed metho d has the adv an tage that recognition reac hes b ey ond the visible faces that are usually seen and used b y other face recognition metho ds, whic h can b e confused b y v arious camourage using plastic surgery or mak eup. Sev eral issues that need to b e addressed in future in v estigations are: (1) a larger data set is needed to thoroughly ev aluate the p erformance; (2) corresp ondence matc hing needs 85 PAGE 96 to b e automated. Optical ro w and motion measuremen t in video sequence also can b e considered; (4) using strain pattern on the whole face rather than a p ortion of face; (5) a comparativ e study with other metho ds is necessary 86 PAGE 97 CHAPTER 6 SENSITIVITY ANAL YSIS In comparison to geometrical and massspring mo dels, ph ysical mo dels that are based on con tin uum mec hanics ha v e high computational complexit y As a result, v arious assumptions are often made to simplify the mo del and its parameters. F or example, a commonly used assumption is that the material prop erties of an ob ject are isotropic and homogeneous. Ho w ev er, results from large amoun t of biomec hanical tests [42 ] indicate that the mec hanical b eha vior of man y biological materials, esp ecially soft tissues, cannot b e accurately describ ed b y suc h a simplied mo del. Certain t yp es of m uscles (suc h as sk eletal m uscle) are c haracterized b y strong anisotropic b eha viors. More imp ortan tly the prop ert y heterogeneit y of sev eral orders of magnitude is also common in h uman organs [28 ]. Recognizing the inadequacy of simplied ph ysical mo dels, researc hers ha v e started to in v estigate to what degree the v arious assumptions, esp ecially those ab out the material prop erties and the b oundary conditions, ma y aect the mo del's p erformance b y means of sensitivit y analysis. F or example, Altero vitz et al [4 ] studied the inruence of b oth the ph ysiciancon trolled parameters and the in trinsic material parameters on the accuracy of needle insertion sim ulation. In the study of mo delbased breast cancer diagnosis, T anner el al [123 ] compared the results of biomec hanical mo dels with dieren t settings of b oundary condition and material prop erties. Ho w ev er, those studies w ere done on the caseb ycase basis using ad ho c comparison metho ds, and therefore the conclusions cannot b e readily generalized to other domains. Moreo v er, the exp erimen ts w ere p erformed based on the assumption that the mo del has homogeneous material prop erties, whic h implies that the solution of a Diric hlet t yp e problem could b e indep enden t of the in ternal prop ert y v ariation, and hence the subsequen t sensitivit y analysis results ma y not b e v alid. Another limitation 87 PAGE 98 of their metho ds is that the sensitivit y data is incapable of pro viding a complete picture of the mo del's spatial resp onse to the parameter v ariation on eac h p oin t of the mo del. In this section, a lo cal gradien tbased computation metho d is discussed that can b e used to conduct a systematic and comprehensiv e sensitivit y analysis of an y motion mo del. Sp ecically ph ysicsbased nonrigid motion mo deling will b enet from suc h a sensitivit y analysis in the follo wing asp ects: 1. the prop osed metho d allo ws us to compare the relativ e imp ortance of dieren t parameters using the normalized sensitivit y data and quan tify the impact of v arious assumptions on mo del's p erformance using the dimensional sensitivit y data; 2. the algorithm is designed based on the adjoin t state metho d, whic h signican tly reduces the computational cost and is suited for handling large scale n umerical mo dels; 3. the sensitivit y con tour map enables us to iden tify the vulnerable areas of a mo del that are most aected b y a p o or assumption, so that further impro v emen t can b e made. 4. as disscussed in the previous sections, the parameter v alue can b e obtained b y either the direct measuremen t or the indirect inference. But those acquisition pro cedures are timeconsuming and exp ensiv e. It w ould b e economic to rst conduct a sensitivit y analysis to iden tify the primary parameters and then to concen trate our eort on the acquisition of those parameters. 6.1 Computational Metho d Sensitivit y analysis is closely related to optimization problems often encoun tered in the traditional mo del calibration tasks, suc h as optimal shap e design, b oundary condition sp ecication, discretization strategy in v estigation, as w ell as material prop ert y assignmen t [72 ]. In this study the fo cus is on assessing the impact of t w o elastic parameters (the Y oung's mo dulus and the P oisson's ratio) on mo del's p erformance, based on the computed sensitivit y information. Without loss of generalit y w e will discuss the deformation of a linear elastic 88 PAGE 99 b o dy and its resp onse to parameter p erturbation. Ho w ev er, the metho dology dev elop ed here can b e readily applied to nonlinear systems. 6.1.1 Primary Problem Giv en the go v erning equations of elastic dynamics, a v ariational metho d can b e used to deriv e the nite elemen t form ulation of the partial dieren tial equations so that a n umerical solution of the forw ard mo del can b e obtained. F or example, using the Galerkin metho d [151 ], the primary problem in the static equilibrium can b e discretized in to a linear nite elemen t system: K u = b; (97) where K is the stiness matrix, u denotes the no dal solution, and b is the load v ector. F or eac h nite elemen t, an in terp olation matrix ( R i ) (the basis function) is used to relate the displacemen ts in the lo cal elemen t co ordinate ( d i ) and the global co ordinate ( u ), whic h ensures that the requiremen t of strain compatibilit y is satised: d i = R i u: (98) After the discretization, the corresp onding constitutiv e equation (2) and the straindisplacemen t equation (3) in the discrete matrix form are: i = K i i ; (99) i = D i u; (100) where K i is the co ecien t of stiness matrix and D i is a matrix deriv ed from the spatial dieren tiation of R i and the com bination of ro w items. Note that K i is comprised of the material prop erties ( E ; ) that will b e studied in the sensitivit y analysis. 89 PAGE 100 With the ab o v e notations, the stiness matrix and the load v ector in (8) that assem ble all the nite elemen ts b ecome: K = M X i =1 Z V D T i K i D i dV ; (101) b = M X i =1 Z V R T i q 1 i dV + M X i =1 Z S R T i q 2 i dS + L X j =1 q 3 j ; (102) where M is the total n um b er of elemen ts, V denotes the size of elemen t, q 1 i is the b o dy force within an elemen t, q 2 i is the surface traction term, q 3 j is the p oin t load, and L is the n um b er of p oin t loads exerted on the b o dy Since the in terest is in the sensitivit y assessmen t of a mec hanical system go v erned b y (8) to the parameter v ariations, (8) is rewritten in terms of the parameter v ector: K ( p ) u = b ( p ) ; (103) where p denotes the generic parameter to b e in v estigated, whic h could b e the material prop erties, the stress or strain distributions, the loading conditions (either the Diric hlet t yp e or the Neumann t yp e), ob ject's geometrical conguration, and ev en the meshing sc heme itself. Here w e concen trate on the sensitivit y analysis of t w o spatial parameters, the Y oung's mo dulus and the P oisson's ratio: p = f p j g = f E j ; j g ; j = 1 ; 2 ; ::: M 6.1.2 P erformance Measure and Dimensional Sensitivit y The simplest sensitivit y analysis metho d that can b e done in nonrigid motion mo deling is based on a T a ylor series describing the relationship b et w een the state v ariable ( u ) and the parameter ( p ): u ( p + 4 p ) = u ( p ) + M X j =1 @ u @ p j 4 p j + 1 2 M X j =1 M X l =1 @ 2 u @ p j @ p l 4 p j 4 p l + ::: (104) 90 PAGE 101 where the rstorder deriv ativ es ( @ u @ p j ) are the co ecien ts of the Jacobian matrix: J = [ @ u i @ p j ] ; ( i = 1 ; 2 ; ::: N ; j = 1 ; 2 ; ::: M ), with N and M as the total n um b er of no des and the total n um b er of elemen ts in a mo del, resp ectiv ely The Jacobian matrix pro vides a lo cal description of the resp onse of the state v ariable to the parameter c hange and is adequate for analyzing a simple system suc h as the massspring mo del. Ho w ev er, for most practical engineering problems, more sophisticated sensitivit y analysis metho d is needed to help us understand ho w a giv en mo del dep ends on certain parameters, so that the most signican t ones can b e iden tied and impro v ed. T o this end, a p erformance function (ob jectiv e function) has to b e dened: H ( u; p ) = Z n f ( u; p ) d n : (105) The function dened in (16) is only a generic in tegral form. The exact form ula of f ( u; p ) is often taskdep enden t, and the use of Euclidean norm is common. Giv en a p erformance function H ( u; p ), the dimensional sensitivity ( S ) can b e dened as its total deriv ativ e with resp ect to the parameter v ector: S = dH ( u; p ) dp = Z n @ f ( u; p ) @ p + @ f ( u; p ) @ u du dp d n ; (106) or in the discrete form: S = [ S 1 ; S 2 ; :::S j ; :::S M ] ; (107) S j = dH ( u; p ) dp j = @ H ( u; p ) @ p j + @ H ( u; p ) @ u du dp j : (108) Without causing an y confusion, from no w on, (18) and (19) will b e represen ted with a single expression: S = dH ( u; p ) dp = @ H ( u; p ) @ p + @ H ( u; p ) @ u du dp : (109) 91 PAGE 102 6.1.3 Adjoin t State Metho d The dimensional sensitivit y can b e computed using the nite dierence appro ximation, the direct dieren tiation or the adjoin t state metho d [18 ]. The adjoin t state metho d is c hosen b ecause it is more ecien t for large scale problems. F or a nite elemen t mo del of mo dest size that are commonly encoun tered in ph ysicsbased nonrigid motion sim ulation, the n um b er of elemen ts could easily reac h a couple of thousands. Therefore, the direct computation of sensitivit y distribution w ould b e to o costly The adjoin t state metho d is based on the principle of v ariational calculus [56 ], and can b e deriv ed from the go v erning partial dieren tial equations in the con tin uous space (15), or from the linear system in the discrete space (14). This study will discuss the adjoin t state approac h follo wing the second route. (14) is rst dieren tiated with resp ect to the parameter p : @ K ( p ) @ p u + K ( p ) du dp = @ b ( p ) @ p ; (110) (21) is then m ultiplied b y an arbitrary v ector (adjoin t state v ariable): T @ K ( p ) @ p u + T K ( p ) du dp = T @ b ( p ) @ p : (111) Com bining (20) and (22), w e ha v e: dH ( u; p ) dp = @ H ( u; p ) @ p + @ H ( u; p ) @ u du dp T K ( p ) du dp T @ K ( p ) @ p u + T @ b ( p ) @ p : (112) Since is arbitrary a is c hosen suc h that the follo wing equation is satised (the second and third terms in the righ t side of (23): @ H ( u; p ) @ u du dp T K ( p ) du dp = 0 : (113) 92 PAGE 103 (24) is equiv alen t to: ( du dp ) T K T ( p ) = ( du dp ) T ( @ H ( u; p ) @ u ) T ; (114) whic h leads to the adjoin t state equation: K ( p ) = ( @ H ( u; p ) @ u ) T (115) In the ab o v e deriv ation, w e to ok adv an tage of the fact that the stiness matrix is symmetric. Also note that the adjoin t problem (26) is in the same form as the primary problem (14), except for the load term, ( @ H ( u;p ) @ u ) T whic h is also kno wn as the pseudoload. Substituting (24) bac k in to (23), it leads to: dH ( u; p ) dp = @ H ( u; p ) @ p + T ( @ b ( p ) @ p @ K ( p ) @ p u ) : (116) As a result, w e only need to solv e the primary problem (14) and the adjoin t problem (26) once, in order to obtain the state solution ( u ) and the adjoin t state solution ( ), and then use (27) to compute the dimensional sensitivit y The adjoin t state form ulation can also b e deriv ed using the Lagrange m ultiplier metho d b y treating (14) as the state equation, (16) as a constrain t function and as the Lagrange m ultiplier [26 ]. 6.1.4 Deriv ativ e Computation The partial deriv ativ e terms in (27) are needed in order to calculate the dimensional sensitivit y If the p erformance function is dened in an explicit form, the second term in (27) can b e computed as: @ H ( u; p ) @ p = Z n @ f ( u; p ) @ p d n : (117) In most applications, f ( u; p ) is dened as a function of the state v ariable only f ( u ), and therefore its partial deriv ativ e with resp ect to the parameter b ecomes zero. 93 PAGE 104 With a nite elemen t form ulation, the partial deriv ativ es of the stiness matrix ( K ) and the load v ector ( b ) with resp ect to the parameter ( p ) can b e expressed as: @ K ( p ) @ p = M X i =1 Z V D T i @ K i ( p ) @ p D i dV ; (118) @ b ( p ) @ p = M X i =1 Z V R T i @ q 1 i ( p ) @ p dV + M X i =1 Z S R T i @ q 2 i ( p ) @ p dS + L X j =1 @ q 3 j ( p ) @ p : (119) Since the in terest is in the sensitivit y v alues with resp ect to the material prop erties only ( p = f E ; v g ), the deriv ativ e computations in v olv ed in (29) and (30) will only cause a sligh t increase of the computational cost. 6.1.5 Normalized Sensitivit y Dimensional sensitivit y as dened in (20) is suitable for studying the dep endence of a mo del on a single parameter (note that the v ector p = [ p 1 ; p 2 ; :::p M ] of a nite elemen t mo del is view ed as a single parameter here). Because eac h ph ysical parameter usually has its o wn unit, a comparison among dieren t parameters suc h as the Y oung's mo dulus and the P oisson's ratio using the dimensional sensitivit y is v ery dicult and p oten tially could b e misleading. In order to remo v e the dimensional eect, a normalize d sensitivity ( S n ) is dened as: S n = p H ( u; p ) dH ( u; p ) dp = p H ( u; p ) @ H ( u; p ) @ p + @ H ( u; p ) @ u du dp (120) The normalized sensitivit y measures the p ercen tage c hange of the p erformance function to the p ercen tage c hange of parameters. Therefore, it allo ws us to ev aluate and rank the imp ortance of dieren t parameters in terms of their impact on the mo del on a common basis. Note that the normalized sensitivit y do es not require an y c hanges of the computational 94 PAGE 105 pro cedure (adjoin t state metho d) dev elop ed for the dimensional sensitivit y If the statistical information ab out the measuremen t data is a v ailable, a normalized sensitivit y function w eigh ted b y the measuremen t condence co ecien t can also b e considered. 6.2 Exp erimen ts with Syn thetic Mo del 6.2.1 Mo del Conguration and P erformance Measure The mo del used in this exp erimen t represen ts a 2D elastic ob ject that has a size of 10x8 ( cm 2 ). The ob ject is discretized in to a nite elemen t mesh with 357 no des and 320 elemen ts (Figure 28). The elemen ts are of quadrilateral thin shell t yp e. All elemen ts are assigned with isotropic prop erties ( E = 50 kP a, = 0.460), except for the 32 elemen ts in the middle of the mo deling domain that ha v e higher v alues ( E = 500 kP a, = 0.495), creating the prop ert y heterogeneit y This heterogeneously distributed prop erties are regarded as the true parameters: p t = f E t ; t g The Diric hlet b oundary conditions are sp ecied as follo ws: all no des on the left side b oundary of the mo del are xed (displacemen t = 0 cm), while the compression loads are exerted along the righ t side b oundary (displacemen t = 3 cm). The no des on the top and b ottom b oundaries are set free. With the true parameters and b oundary conditions as sp ecied ab o v e, the observ ation data ( u ) are generated b y running a forw ard sim ulation. The data v ector is dened on all no des ( u = [ u 1 ; u 2 ; ::: u N ] T ). It should b e noted that, in real applications, the state v ariable ma y b e measured only on parts of the ob ject. Let's rst dene a p erformance function in the form of discrete Euclidean norm: H ( u ) = Z n f ( u ) d n = jj f ( u ) jj 2 = N X i =1 w 2 i [ u i ( p a ) u i ] 2 (121) where w i is a w eigh t co ecien t assigned to the ith no de, p a denotes the assumed parameters, u ( p a ) is the solution that is obtained using the assumed parameters and u represen ts the 95 PAGE 106 E = 50 kPa v = 0.460 E = 500 kPa v = 0.495 Fixed Fixed Fixed Fixed Fixed 10 cm 3 (cm) 3 (cm) 3 (cm) 3 (cm) 3 (cm) 8 cmFigure 28. Syn thetic Mo del Conguration. observ ation data that is generated with the true parameters. In all of the exp erimen ts, the unit w eigh t co ecien t: w i = 1 is used. This p erformance function consists of the state v ariable only So the marginal sensitivit y (the second term in (27)) b ecomes zero. In situations where the parameter measuremen ts are a v ailable, they can also b e incorp orated in to the p erformance function ( H ( u; p )). If stress or strain analysis is of in terest, they can also b e included to construct a Neumann t yp e p erformance function. Since homogeneit y is a commonly used assumption, a mo del is tested that has the homogeneous prop erties: p a = f E a = 50 kP a, a = 0 : 460 g The deformed mo dels with the true heterogeneous prop erties ( p t = f E t ; t g ) and the assumed homogeneous prop erties ( p a = f E a ; a g ) are sho wn in Figure 29. The discrepancy b et w een the t w o deformed shap es is clearly visible. Next, it will b e demonstrated that (1) the normalized sensitivit y can b e used to iden tify the parameter that is most resp onsible for the discrepancy; (2) the dimensional sensitivit y is helpful for determining ho w go o d an assumption is. 96 PAGE 107 (a) (b) Figure 29. Deformed Mo dels. (a) the true heter o gene ous pr op erties and (b) the assume d homo gene ous pr op erties. 6.2.2 Normalized Sensitivit y Distribution The normalized sensitivities with resp ect to the Y oung's mo dulus and the P oisson's ratio are presen ted in Figure 30. The sensitivit y distribution is plotted elemen t b y elemen t as a con tour map ( S n = [ S n 1 ; S n 2 ; :::S nM ] T ). The sensitivit y co ecien t in an elemen t ( S nj ) represen ts the p ercen tage c hange of the p erformance function due to a 1% c hange of the parameter v alue in that elemen t. The dierence b et w een the t w o maps in terms of the magnitude of sensitivit y co ecien t suggests that the mo del is more sensitiv e to the Y oung's mo dulus than to the P oisson's ratio. In other w ords, with this particular mo del conguration, Y oung's mo dulus has a signican t impact on the mo deling accuracy while the eect of the P oisson's ratio is secondary F urther insp ection of Figure 30 (a) rev eals that the mo del is most sensitiv e at the propert y discon tin uities. The geometry of the high sensitivit y area matc hes the shap e of prop ert y heterogeneit y So, the sensitivit y map is also informativ e ab out the lo cation of parameter abnormalities to whic h eort should b e dedicated to correct the original homogeneous assumption. 97 PAGE 108 0.85 0.85 0.80 0.68 0.59 0.43 0.21 0.80 0.68 0.59 0.43 0.21 0.02 0.02 0.04 0.06 0.04 0.06 0.01 0.01 (a) S n of the Y oung's mo dulus (b) S n of the P oisson's ratio Figure 30. Normalized Sensitivit y Con tour Maps. 9.1e10 8.4e10 7.1e10 4.9e10 2.6e10 2.6e10 4.9e10 7.1e10 8.4e10 9.1e10 9.7e12 4.4e12 1.6e12 1.6e12 4.4e12 9.7e12 (a) S with E heter o = 70 kP a (b) S with E heter o = 440 kP a Figure 31. Dimensional Sensitivit y Con tour Maps. 98 PAGE 109 6.2.3 Dimensional Sensitivit y Distribution It should b e p oin ted out that certain form ulations of the p erformance function ma y result in a sensitivit y v alue of innit y F or example, with a p erformance function as dened in (32), the examination of (31) suggests that, as the assumed parameters approac h the true parameters ( p a p t ), the normalized sensitivit y go es to innit y ( S n 1 ). Although it rarely o ccurs in real applications that our rst assumption is as go o d as the true parameters, it is still wise to in terpret the normalized sensitivit y with cautions. In the con text of nonrigid motion mo deling, a t w ostep sensitivit y analysis pro cedure is recommended: (1) the normalized sensitivit y can b e used to quan tify the relativ e imp ortance of parameters in terms of their impact on the mo del's b eha vior; (2) once a parameter is considered as signican t, the dimensional sensitivit y can b e used to assess the assumptions ab out the parameter. Without m uc h loss of rigor, the closer the assumed parameters ( p a ) are to the true parameters ( p t ), the smaller the dimensional sensitivit y is. The exp erimen ts in v estigate t w o assumptions ab out the Y oung's mo dulus inside the area of heterogeneit y (the 32 elemen ts in the cen ter of the mo del): E 1 = 70 kP a and E 2 = 440 kP a. The Y oung's mo dulus in the bac kground elemen ts and the P oisson's ratio are the same as the true v alues. The dimensional sensitivit y con tours using those t w o assumptions are sho wn in Figure 31 (unit = [ m=P a ]). The smaller magnitude of the second map indicates that the assumption ( E 2 = 440 kP a) is m uc h closer to the true v alue ( E t = 500 kP a) and is a b etter assumption than the rst one ( E 1 = 70 kP a). In the study of parameter estimation, the in v erse algorithms also utilize the sensitivit y information, although the regularization is more imp ortan t in solving the illp osed in v erse problems [148 ]. 6.2.4 Relativ e Error Distribution The distribution of relativ e mo deling error can also c haracterize the complex relationship b et w een a mo del's b eha vior and its parameters. The relativ e mo deling error is dened as follo w: er r = j u ( p a ) u j j u j : (122) 99 PAGE 110 0.41 0.36 0.21 0.03 0.08 0.31 0.22 0.10 0.03 0.03 0.0003 0.0003 0.0003 0.0004 0.0008 0.0003 0.0003 0.0004 0.0004 0.0001 (a) er r with E a and t (b) er r with E t and a Figure 32. Relativ e Error Distributions. E t : true Y oung's mo dulus, E a : assume d Y oung's mo dulus, t : true Poisson 's r atio, a : assume d Poisson 's r atio. (a) scar with grid (b) deformed scar Figure 33. Deformation of Burn Scar. T o determine the con tribution to the relativ e mo deling error b y an individual parameter, t w o exp erimen ts w ere carried out, one with the true heterogeneous Y oung's mo dulus and the assumed homogeneous P oisson's ratio, and one with the assumed homogeneous Y oung's mo dulus and the true heterogeneous P oisson's ratio. The results are sho wn in Figure 32. The errors caused b y an incorrect homogeneous assumption of the Y oung's mo dulus is around 10% to 40%, while the errors due to a bad assumption of the P oisson's ratio is less than 0.08%, whic h conrms the conclusion dra wn from the sensitivit y analysis that the Y oung's mo dulus is more signican t than the P oisson's ratio under the conguration of this syn thetic mo del. 100 PAGE 111 0.16 0.14 0.12 0.10 0.06 0.08 0.06 0.04 0.04 0.03 0.01 0.02 0.02 0.02 (a) S n of the Y oung's mo dulus (b) S n of the P oisson's ratio Figure 34. Normalized Sensitivit y Maps of Scar. 6.3 Exp erimen ts with Burn Scar The exp erimen t with syn thetic mo del has implications to man y practical mo deling problems, b ecause the prop ert y v alues used in the exp erimen t are based on the published data. F or example, in needle insertion sim ulation, a Y oung's mo dulus in the range of 40160 kP a w as used for prostate tissue [4 ]. In breast mo deling, the v ariation of E among dieren t tissues can b e up to an order of magnitude, 1088 kP a for the skin tissue and 1184 kP a for glandular tissue [10 123 ]. The P oisson's ratio used in soft tissue mo deling is often in the range of 0.4650.495. Therefore, the mo deling scenario of syn thetic mo del is v ery realistic and the sensitivit y analysis demonstrated with the syn thetic mo del can b e readily extended to real ob jects. In the next exp erimen t, the sensitivit y metho d will b e applied to burn scar. 6.3.1 Mo del Conguration A patien t's arm with a burn scar w as used in the exp erimen t. Using an ink stamp, a rectangle grid w as created on the arm's surface. Tw o images w ere then tak en b efore and after the patien t's arm w as stretc hed (see Figure 33). A nite elemen t mo del with quadrilateral thin shell elemen t w as then constructed that matc hes the rectangle grid (see [135 ] for detailed descriptions of the mo deling pro cedure). A Diric hlet condition w as sp ecied around the b oundary using the measured displacemen t data and the in ternal no des w ere set free. 101 PAGE 112 6.3.2 Normalized Sensitivit y The nite elemen t mo del co v ers b oth the scar (in the middle) and the surrounding normal skins. A homogeneous prop ert y distribution throughout the mo del w as assumed for b oth the Y oung's mo dulus and the P oisson's ratio ( E = 60 kP a, = 0 : 495). The normalized sensitivit y con tours computed with the ab o v e setting are sho wn in Figure 34. Again, the larger v alues in Figure 34 (a) indicate that the scar mo del is more sensitiv e to the c hange of the Y oung's mo dulus than to that of the P oisson's ratio, although the dierence b et w een the sensitivit y v alues of E and is not as large as those observ ed in the exp erimen ts with syn thetic mo del. It can b e seen that the region of high normalized sensitivit y v alues matc hed w ell with the o v erall shap e of scar (in b oth the Y oung's mo dulus map and the P oisson's ratio map). This similarit y b et w een the scar b oundary and the sensitivit y con tours indicates that the mo del is v ery sensitiv e to the prop ert y c hange b et w een scar and normal skin. It has b een observ ed that the Y oung's mo dulus of the burn scars could b e 210 times higher than that of the normal skin (5100 kP a) [135 ]. Based on the normalized sensitivit y map of the P oisson's ratio in Figure 34 (b), it seems reasonable to exp ect that the same t yp e of heterogeneit y of P oisson's ratio exists b et w een the scar and the normal skins, although the magnitude of the heterogeneit y still remains unclear. 6.3.3 Dimensional Sensitivit y A homogeneous elasticit y is certainly not adequate for scar assessmen t study The next example is used to demonstrates ho w to nd a b etter prop ert y v alue based on the dimensional sensitivit y information. The Y oung's mo dulus v alue of scar (8 elemen ts in the middle of the mo del) w as gradually increased from the initial 60 kP a un til no further impro v emen t can b e observ ed in terms of the mo del's dimensional sensitivit y distribution. T o quan tify the ev aluation pro cess, an a v erage dimensional sensitivit y v alue is used: S av e = 1 M M X j =1 S j : (123) 102 PAGE 113 8.1e11 5.5e11 4.5e11 3.2e11 2.0e11 2.6e12 5.8e12 7.2e12 2.0e11 1.2e11 (a) S with E scar = 60 kP a (b) S with E scar = 450 kP a Figure 35. Dimensional Sensitivit y Maps of Scar. (a) Closed mouth (b) Op en mouth Figure 36. F ace Deformation with Expression. Note that the skin c over e d by the tap e has r elatively smal l deformation. As discussed in the previous section, the smaller S av e is, the closer the corresp onding assumption is to the true distribution. The nal Y oung's mo dulus v alue of scar whic h resulted in a minim um S av e is 450 kP a. The dimensional sensitivit y con tour maps of t w o assumptions ( E scar = 60 kP a and E scar = 450 kP a) are sho wn in Figure 35. The smaller sensitivit y co ecien t v alues in the second map (b) clearly indicates that the heterogeneous assumption ( E back g r ound = 60 kP a, E scar = 450 kP a) is a b etter c hoice than the original homogeneous assumption ( E back g r ound = 60 kP a, E scar = 60 kP a). Note that the landscap e of the dimensional sensitivit y con tour in Figure 35 (b) b ecomes more \rat" than that in Figure 35 (a), whic h is a result of using more accurate prop ert y v alues in the mo del. 103 PAGE 114 6.4 Exp erimen ts with Deformed F ace Mo del F ace recognition is usually p erformed using static in tensit y images [103 137 ]. Recen tly the in v estigation of face biometrics has b een extended to the dynamic domain suc h as facial expression [142 143 ], motion eld [116 ], and strain patterns [147 ]. The new biometrics attempt to explore the uniqueness of soft tissue prop erties and the asso ciated elastic patterns exhibited during face expressions. Ho w ev er, no quan titativ e ph ysical link b et w een the h yp othesis of new biometrics and the existence of prop ert y abnormalit y has b een established. In the next exp erimen t, it is sho wn that the sensitivit y analysis ma y pro vide quan titativ e measures that are useful for the design of new face biometrics. 6.4.1 F ace Mo del Conguration Tw o prole images (side view) of a sub ject's face w ere tak en using a Minolta 3D range camera, one with mouth closed and the other one with mouth op en (see Figure 36). T o facilitate corresp ondence matc hing, a set of p oin ts (8x8) w as mark ed on the sub ject's face. A rectangular transparen t tap e w as attac hed on the face surface close to the lo w er ja w. The tap e is less elastic and is used to sim ulate prop ert y c hange of facial tissue resulting from either a plastic surgery or an acciden t. The face mo del consists of 64 no des and 49 quadrilateral elemen ts. F ace deformation measured from t w o images will then b e used b y the nite elemen t mo del to infer an y prop ert y c hange. 6.4.2 Normalized Sensitivit y The normalized sensitivit y w as used to determine whic h parameter, the Y oung's mo dulus or the P oisson's ratio, is more informativ e ab out the existence of face prop ert y c hange. In other w ords, to whic h parameter v ariation the face mo del is more sensitiv e. As in the previous exp erimen ts, a Diric hlet condition w as sp ecied along the mo del b oundary and the in ternal no des w ere set free. Homogeneous prop erties ( E = 50 kP a, = 0 : 495) w ere also assumed. The resulting normalized sensitivit y distributions are plotted in Figure 37. 104 PAGE 115 0.18 0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.06 0.04 0.02 0.01 0.005 0.002 0.005 (a) S n of the Y oung's mo dulus (b) S n of the P oisson's ratio Figure 37. Normalized Sensitivit y Maps of F ace Mo del. In b oth maps, the dense con tour lines that surround the area of high S n v alues clearly outline the lo cation/shap e of the tap e. The a v erage sensitivit y v alue inside the tap e is 0.148 for the Y oung's mo dulus and 0.013 for the P oisson's ratio, a dierence of ab out one order of magnitude. So, the face mo del is more sensitivit y to the c hange of the Y oung's mo dulus than to the c hange of the P oisson's ratio. This observ ation is consisten t with that dra wn from the exp erimen ts using syn thetic mo del and burn scar. The quan titativ e sensitivit y information has follo wing implications: (1) it is easier to detect/estimate the c hange of the Y oung's mo dulus than c hange of P oisson's ratio; (2) The elasticit y pattern is a b etter candidate to design an eectiv e face biometrics. 6.4.3 Dimensional Sensitivit y In the exp erimen t with burn scar, w e sho w that the dimensional sensitivit y can b e used to impro v e the initial v alue b y c hanging the prop ert y v alue inside the abnormal area. The pro cedure is p ossible b ecause w e kno w the exact lo cation of the abnormal area. F or example, the b oundary b et w een the scar and the normal skin is clearly visible on the images. Ho w ev er, in face recognition, the c hange in visible cues (color, texture and in tensit y) ma y not corresp ond to an y underlying prop ert y abnormalities, either due to the relativ ely lo w 105 PAGE 116 image qualit y or the mak eup put on the face. So, w e devise a new pro cedure that utilizes the dimensional sensitivit y itself to guide the searc h for a b etter prop ert y v alue. The exp erimen ts start with an initial homogeneous prop ert y and compute the dimensional sensitivit y Elemen ts are then selected that ha v e higher sensitivit y v alues as the p oten tial abnormal areas. The prop ert y v alue inside an abnormal area will b e c hanged un til its sensitivit y v alue is close to the a v erage of en tire mo del. The ab o v e steps will b e rep eated un til the sensitivit y v alues throughout the mo del are b elo w a usersp ecied threshold. In other w ords, the up dated v alue will get closer to the true prop ert y v alues after eac h step. A few in termediate results and the nal sensitivit y con tours are sho wn in Figure 38. Based on the initial sensitivit y v alues, four elemen ts w ere selected as the abnormal area, all co v ered b y the tap e. The Y oung's mo dulus of the four elemen ts w ere then c hanged to 230 kP a and 810 kP a, resp ectiv ely with the corresp onding sensitivit y maps sho wn in Figure 38 (a,b). Note that the sensitivit y v alues inside the tap ed area are lo w er than that of the surrounding regions. Using the information in Figure 38 (b), eigh t more elemen ts that ha v e higher sensitivit y v alues w ere selected, three b elo w the tap e and v e ab o v e the tap e. The Y oung's mo dulus of the eigh t elemen ts w ere c hanged to 150 kP a and 370 kP a, resp ectiv ely and the nal sensitivit y is giv en in Figure 38 (d). The nal sensitivit y map has m uc h smaller magnitude and sho ws relativ ely uniform landscap e, implying that the nal Y oung's mo dulus v alue is more accurate, and therefore more suited for face recognition using a ph ysical mo del. Mo delbased nonrigid and articulated motion analysis requires careful design and calibration of the mo del b eing used, whether it is a geometrical mo del or a ph ysical mo del. V arious assumptions ab out the mo del and its parameters m ust b e thoroughly ev aluated. A d ho c sensitivit y analysis metho d is impractical for this task, esp ecially when dealing with large scale problems. a systematic sensitivit y analysis metho d is prop osed that is capable of making a quan titativ e and reasoned diagnosis of the mo del's p erformance related to our assumptions ab out the parameters. The prop osed metho d is form ulated using the rstorder 106 PAGE 117 6.0e11 2.0e11 1.4e10 1.0e10 2.0e11 6.0e11 1.0e10 1.4e10 1.8e10 4.0e11 2.0e11 1.0e11 6.0e11 8.0e11 4.0e11 2.0e11 1.0e11 (a) E abnor mal = 230 kP a (b) E abnor mal = 910 kP a in 4 abnormal elemen ts in 4 abnormal elemen ts 4.0e11 3.0e11 1.0e11 5.0e11 2.0e11 3.0e11 4.0e11 3.0e11 2.0e11 1.0e11 4.0e11 3.0e11 0.7e12 2.0e11 1.0e11 4.0e11 2.0e11 1.0e11 (c) E abnor mal = 150 kP a (d) E abnor mal = 370 kP a in 8 abnormal elemen ts in 8 abnormal elemen ts Figure 38. Dimensional Sensitivit y Maps of F ace Mo del. 107 PAGE 118 lo cal gradien t information. The adjoin t state metho d is emplo y ed to reduce the computational cost of large scale n umerical mo dels. A t w ostep pro cedure is recommended: (1) the normalized sensitivit y is suited for iden tifying the most signican t parameters; (2) the dimensional sensitivit y can b e used to assess and impro v e our assumption ab out a particular parameter. Based on the exp erimen ts with syn thetic mo dels, burn scars, and faces, sev eral observ ations can b e made: (1) the mo dels are more sensitiv e to the c hange of the Y oung's mo dulus than to the c hange of the P oisson's ratio; (2) the mo dels are most sensitiv e at the prop ert y discon tin uities (heterogeneit y); (3) sensitivit y map is informativ e ab out the lo cation/geometry of the prop ert y abnormalities; It should b e stressed that the exp erimen ts rep orted here w ere p erformed using a linear elastic mo del. The system resp onse of a nonlinear mo del to the parameter v ariation is more complex. In future in v estigations, the prop osed metho d will b e applied to address the follo wing issues: (1) Ho w sensitiv e an elastic mo del is to an assumed isotropic prop ert y? (2) What is the sensitivit y distribution of an elastic mo del under v arious loading conditions (b oundary conditions)? (3) What is the sensitivit y resp onse of a nonlinear mo del? Those cases will b e studied that are geometrically nonlinear (Green strain tensor) and materially nonlinear (visco elastic and plastic materials). 108 PAGE 119 CHAPTER 7 CONCLUSIONS AND DISCUSSIONS A new ph ysicsbased approac h for mo deling nonrigid motion is prop osed. The no v elt y of the prop osed metho d lies in the fact that the motion parameters are actual ph ysical parameters. The ph ysical parameters ha v e clear ph ysical meanings, and more imp ortan tly they only need to b e reco v ered once. In the trac king exp erimen ts, a sc heme is demonstrated that uses the reco v ered elasticit y and the b oundary corresp ondence to syn thesize the degrees of freedom of the in ternal no des, whic h reduces the trac king complexit y signican tly The prop osed mo deling approac h can b e applied to a wide range of domains. In particular, its application is demonstrated in burn scar assessmen t and face recognition. In burn scar study an metho d based on relativ e elasticit y ratio and natural features is dev elop ed. The use of natural features and adaptiv e mesh not only increases the computational eciency but also allo ws us to w ork with scar images without the need of tagging articial mark ers. Burn scar damage is quan tied b y estimating its Y oung's mo dulus using a simplied regularization metho d. Exp erimen ts sho w that the prop osed metho d is robust when presen ted with noisy data. The p ositiv e correlation b et w een the ph ysician's rating and REIs (using b oth natural features and articial mark ers) suggests that the prop osed approac h can pro vide a quan titativ e and ob jectiv e ev aluation of burn scar damage. In face recognition study a new biometrics is prop osed that is based on the consideration of anatomical and biomec hanical c haracteristics of facial tissues. Elastic strain pattern inferred from facial expression can rev eal an individual's biometrics signature asso ciated with the underlying anatomical structure, and th us has the p oten tial for face recognition. A metho d based on the con tin uum mec hanics in nite elemen t form ulation is emplo y ed to compute the strain pattern. Exp erimen ts sho w v ery promising results. The prop osed 109 PAGE 120 metho d is quite dieren t from other face recognition metho ds, and b oth its adv an tages and limitations, as w ell as future researc h for impro v emen t are discussed. A comparativ e study is conducted to ev aluate the strength and w eakness of t w o reco v ery algorithms. In all of the exp erimen ts, b oth the CGA and the GNM sho w ed their robustness against data noise b y con v erging to a stable and nearoptimal solution, attributed to the constrain t imp osed through the p enalt y function. Because of its sto c hastic nature, the CGA is more demanding than the GNM metho d in terms of the need of computational resources. Although the CGA is a computationally in tensiv e metho d, it has a v ery desirable feature of b eing less sensitiv e to the initial conditions. By main taining a div ersied p opulation p o ol, CGA can explore the solution space on a global scale without b eing caugh t in lo cal extrema, whic h has b een the main obstacle to the gradien tbased metho ds. The prop osed motion analysis approac h are v alid for mo deling man y defomable ob jects. But the follo wing areas are iden tied to whic h future w orks need to b e dedicated. Line ar or Nonline ar Mo del: Linear mo del is conceptually and computationally attractiv e and suces our need in most cases. Ho w ev er, large amoun t of biomec hanical exp erimen ts sho wn that man y liv e ob jects, esp ecially soft tissues, do not ob ey the Ho ok e's la w strictly In addition, Cauc h y strain tensor deriv ed from innitesimal deformation assumption is no longer suitable for ob jects of large deformation. When to abandon linear mo del is a tric ky issue b ecause nonlinear mo del is far more complex. The v alidit y of using linear mo del is taskdep enden t and the c hoice is often made based on mo deler's kno wledge on ob ject's material c haracteristics. One p ossible solution is to deriv e empirical constitutiv e form ulation (pseudoelasticit y) from tensile test and use it in motion equation [42 ]. Another p ossibilit y is to carry out a thorough n umerical exp erimen t to assess the feasibilit y of using linear mo del in a particular task. Unfortunately suc h kind of n umerical exercise is rarely practiced in motion analysis. Isotr opic or A nisotr opic Pr op erty: Similar problem arises in the c hoice of isotropic or anisotropic prop ert y Man y biological materials exhibit anisotropic deformation under stress. Anisotropic prop ert y is determined b y material comp osition and substructure (ori110 PAGE 121 en tation and arrangemen t of constitutional units). Blo o d v essels and sk eletal m uscles are go o d examples. F or a 2D anisotropic ob ject, reco v ering its Y oung's mo dulus is equiv alen t to solving an in v erse problem with t w o spatially v arying unkno wns ( E x ; E y ). This should not p ose a sev ere mathematical c hallenge to the reco v ery algorithm, if the constitutiv e and motion equations are prop erly form ulated b y taking anisotropic prop ert y in to accoun t. A lternative R e gularization Metho ds: In addition to the Tikhono v regularization, there exists a large p o ol of alternativ e regularization metho ds that can b e used for nonrigid motion analysis, noticeably the iterativ e regularization and m ultigrid regularization. As to statistical metho ds, w e refer readers to [124 ] for classical applicationorien ted materials and [33 ] for recen t dev elopmen ts. The fact that the Tikhono v functional of nonlinear problem is dicult y to solv e raises the in terest of in v estigating iterativ e regularization metho d. The fundamen tal theory of iterativ e regularization is based on the "semicon v ergence" prop ert y observ ed in the plot of solution error against the iteration coun t. F or example, in a Landw eb er iteration, solution error gradually decreases with iteration coun t and then blo ws out after passing an optimal iteration n um b er. This selfregularization prop ert y allo ws us to stabilize the solution b y stopping iteration early The iteration coun t pla ys a role that is similar to the regularization parameter. Cho osing an appropriate stopping rule is still a c hallenging issue and the discrepancy principle seems to b e the b est c hoice [15 30 ]. Multigrid metho d is w ell kno wn for its stabilit y and eciency in solving largescale w ellp osed problems. In nite dimensions, m ultigrid is exp ected to yield a stable solution for a discretized in v erse problem, but no theory has b een rigorously established. Studies ha v e sho wn that grid renemen t itself (dynamic meshing or adaptiv e meshing) p osses the propert y of regularization [78 8]. One attractiv e feature of m ultigrid metho d is that parameter space and hence computational complexit y can b e signican tly reduced. Grid renemen t can b e done in either parameter space (Y oung's mo dulus) or state v ariable space (displacemen t or force). If renemen t is carried out in parameter space, the nest grid is usually predened and iteration is implemen ted through the \zonation" tec hnique where the pa111 PAGE 122 rameter in a particular zone is assumed constan t. This \zonation" tec hnique eectiv ely reduces the parameter space and has the eect of stabilizing the solution and accelerating the con v ergence. Applications of m ultigrid metho d in estimating tissue elasticit y and aquifer h ydraulic transmissivit y ha v e b een rep orted in [139 6]. One of the diculties in implemen ting m ultigrid metho d is to c ho ose an appropriate rule to stop the renemen t pro cess. The p ossibilit y of com bining the Tikhono v regularization with m ultigrid metho d is also in v estigated [8]. 112 PAGE 123 REFERENCES [1] J. K. Aggarw al and Q. Cai, \Articulated and elastic nonrigid motin: a review," IEEE workshop on Notion of Nonrigid and A rticulate d obje cts Austin, TX, pp. 214, 1992. [2] J. K. Aggarw al, Q. Cai, W. Liao, and B. Sabata, \Nonrigid Motion Analysis: Articulated and Elastic Motion," Computer Vision and Image Understanding 70(2), pp. 142156, 1998. [3] S. Agly amo v, A. R. Sk o v oro da, J. M. Rubin, M. O'Donnell, S. Y. Emeliano v, \Mo delbased reconstructiv e elasticit y imaging of deep v enous throm b osis," IEEE T r ansactions on Ultr asonics, F err o ele ctrics, and F r e quency Contr ol v ol. 51, no. 5, pp. 521531, 2004. [4] R. Altero vitz, K. Goldb erg, J. P ouliot, R. T asc hereau, and ICho w Hsu, \Needle insertion and radioactiv e seed implan tation in h uman tissues: sim ulation and sensitivit y analysis," Pr o c. of the 2003 IEEE Inter. Conf. on R ob otics and A utomation pp. 17931799, 2003. [5] L. Alv arez, P L. Lions, and J. M. Morel, \Image selectiv e smo othing and edge detection b y nonlinear diusion (I I)," SIAM Journal of Numeric al A nalysis v ol. 29, pp. 845866, 1992. [6] H. B. Ameur, M. Burger, and B. Hac kl, \Lev el set metho ds for geometric in v erse problems in linear elasticit y ," Inverse Pr oblems v ol. 20, pp. 673696, 2004. [7] C. H. Arns, M. A. Knac kstedt; W. V. Pinczewski, E. J. Garb o czi, \Computation of linear elastic prop erties from microtomographic images: metho dology and matc h to theory and exp erimen t," Journal of Ge ophysics v ol. 67, no. 5, pp. 13961405, 2002. [8] U. M. Asc her and E. Hab er, \Grid renemen t and scaling for distributed parameter estimation problems," Inverse Pr oblems v ol. 17, pp. 571590, 2001. [9] K. E. A tkinson, A n Intr o duction to Numeric al A nalysis 2nd Edition, Wiley New Y ork, 1989. [10] F. S. Azar, D. N. Metaxas, and M. D. Sc hnall, \Metho ds for mo deling and predicting mec hanical deformations of the breast under external p erturbations," Me dic al Image A nalysis v ol. 6, pp. 127, 2002. [11] R. Ba jcsy and S. Ko v acic, \Multiresolution elastic matc hing," Computer Vision, Gr aphics and Image Pr o c essing v ol. 46, pp. 121, 1989. 113 PAGE 124 [12] P E. Barb one and J. C. Bam b er, \Quan titativ e elasticit y imaging: what can and cannot b e inferred from strain images," Physics in Me dicine and Biolo gy v ol. 47, pp. 21472164, 2002. [13] J. Besag, \Spatial in teraction and the statistical analysis of lattice systems," Journal of the R oyal Statistic al So ciety B, 36(2), pp. 192236, 1974. [14] R. Bev eridge and B. Drap er, "Ev aluation of face recognition algorithms (V ersion 5.0), h ttp://www.cs.colostate.edu/ev alfacerec /index.h tml. [15] A. Binder, M. Hank e, and O. Sc herzer, \On the Landw eb er iteration for nonlinear illp osed problems," J. Inverse Il lp ose d Pr obl. V ol. 4, pp. 381389, 1996. [16] A. Blak e and A. Zisserman, Visual R e c onstruction Cam bridge, MA: MIT Press, 1987. [17] C. Bruyns and M. Ottensmey er, \Measuremen ts of softtissue mec hanical prop erties to supp ort dev elopmen t of a ph ysically based virtual animal mo del," Pr o c. of the Me dic al Image Computing and ComputerAssiste d Intervention, 5th Int. Conf. (MICCAI) T oky o, Japan, pp. 2528, 2002. [18] D. G. Cacuci, Sensitivity and Unc ertainty A nalysis: The ory (V olume I) A CR C Press, 2003. [19] G. F. Carey Computational Grids : Gener ation, A daptation, and Solution Str ate gies. Series in Computational and Ph ysical Pro cesses in Mec hanics and Thermal Science. T a ylor & F rancis, 1997. [20] R. Chelleppa, C. L. Wilson and S. Sirob ey "Human and mac hine recognition of faces : A surv ey" Pr o c e e dings of the IEEE v ol. 83, no. 5, pp. 705740, 1995. [21] C. Chiroiu, L. Mun tean u, V. Chiroiu, P P Delsan to, and M Scalerandi, \A genetic algorithm for determination of the elastic constan ts of a mono clinic crystal," Inverse Pr oblems v ol. 16, pp. 121132, 2000. [22] L. D. Chiwiaco wsky H. F. C. V elho, and P Gasbarri, \A v ariational approac h for solving an in v erse vibration problem," Inverse Pr oblems, Design and Optimization Symp osium Rio de Janeiro, Brazil, 2004. [23] C. C. A. Co ello, \A surv ey of constrain t handling tec hniques used with ev olutionary algorithms," T ec h. Rep ort LaniaRI9904, Lab oratorio Nacional de Informatica Av anzada, 1999. [24] R. Couran t and D. Hilb ert, Metho ds of Mathematic al Physics v ol. 1, John Wiley and Sons, New Y ork, 1953. [25] C. Da v atzik os, \Nonlinear elastic registration of brain images with tumor pathology using a biomec hanical mo del," IEEE T r ans. Me d. Imag. v ol. 18, no. 7, pp. 580592, 1999. [26] K. Dems and Z. Mroz, "V ariational approac h to sensitivit y analysis in thermoelasticit y J. Thermal Str esses (10), pp. 283306, 1987. 114 PAGE 125 [27] M. M. Do yley P M. Meaney and J. C. Bam b er, \Ev aluation of an iterativ e reconstruction metho d for quan titativ e elastograph y ," IEEE Phys. Me d. Biol. v ol. 45, pp. 15211540, 2000. [28] F. A. Duc k, Physic al Pr op erties of Tissues, A c ompr ehensive R efer enc e Bo ok Academic Press, 1990. [29] H. W. Engl, \Iden tication of parameters in p olymer crystallization, semiconductor mo dels and elasticit y via iterativ e regularization metho ds," Il lp ose and Inverse pr oblems (editor, V. G. Romano v), Brill Academic Publisher, 2003. [30] H. W. Engl, M. Hank e, and A. Neubauer, R e gularization of Inverse Pr oblems Klu w er Academic Publishers, c1996. [31] R. Q. Erk amp, A. R. Sk o v oro da, S. Y. Emeliano v, and M. O'Donnell, \Measuring the nonlinear elastic prop erties of tissuelik e phan toms," IEEE T r ansactions on Ultr asonics, F err o ele ctrics, and F r e quency Contr ol v ol. 51, no. 4, pp. 410419, 2004. [32] G. Eskin and J. Ralston, \On the in v erse b oundary v alue problem for linear isotropic elasticit y ," Inverse Pr oblems v ol. 18, pp. 907921, 2002. [33] S. N. Ev ans and P B. Stark, \In v erse problems as statistics," Inverse Pr oblems v ol. 18, pp. R55R97, 2002. [34] M. F atemi, A. Manduca, and J. F. Greenleaf, \Imaging elastic prop erties of biological tissues b y lo wfrequency harmonic vibration," Pr o c e e dings of the IEEE v ol. 91 no. 10 pp. 15031519, 2003. [35] M. F erran t, A. Naba vi, B. Macq, E. A. Jolesz, R. Kikinis, and S. K. W areld, \Registration of 3D in traop erativ e MR images of the brain using a niteelemen t biomec hanical mo del," IEEE T r ans. Me d. Imag. v ol. 20, no. 12, pp. 13841397, 2001. [36] D. B. F ogel, Evolutionary Computation Piscata w a y NJ. IEEE Press, 1995. [37] J. B. F o wlk es, S. Y. Y emely ano v, J. G. Pip e, P L. Carson, R. S. Adler, A. P Sarv azy an, and A. R. Sk o v oro da, \P ossibilit y of cancer detection b y means of measuremen t of elastic prop erties," R adiolo gy v ol. 185, pp. 206207, 1992. [38] D. R. F ranca and A. Blouin, \Alloptical measuremen t of inplane and outofplane Y oung's mo dulus and P oisson's ratio in silicon w afers b y means of vibration mo des," Me as. Sci. T e chnol. v ol. 15, pp. 859868, 2004. [39] L. A. F ried, A natomy of the He ad, Ne ck, F ac e, and Jaws Henry Kimpton Publishers, London, 2nd ed., 1980. [40] P .J. Phillips, P Grother, R.J Mic heals, D.M. Blac kburn, E T abassi, and J.M. Bone, \F ace Recognition V endor T est (FR VT) 2002: Ov erview and Summary", h ttp://www.frvt.org/FR VT2002/do cumen ts.h tm, 2003. [41] Y. C. F ung, F oundations of Solid Me chanics Pren ticeHall, Englew o o d Clis, NJ. 1965. 115 PAGE 126 [42] Y. C. F ung, Biome chanics; Me chanic al Pr op erties of Living Tissues SpringerV erlag, 2nd ed. 1993. [43] T. F uruk a w a and G. Y aga w a, \Inelastic constitutiv e parameters iden tication using an ev olutionary algorithm with con tin uous individuals," International Journal F or Numeric al Metho ds in Engine ering v ol. 40, pp. 10711090, 1997. [44] M. A. Gaballa, T. E. Ra y a, B. R. Simon, and S. Goldman, \Arterial mec hanics in sp on taneously h yp ertensiv e rats: Mec hanical prop erties, h ydraulic conductivit y and t w ophase(solid/ruid) nite elemen t mo dels" Cir c. R es. 71, pp. 145158, 1992. [45] E. J. Garb o czi and A. R. Da y \An algorithm for computing the eectiv e linear elastic prop erties of heterogeneous materials: 3D results for comp osites with equal phase P oisson ratios," Journal of Physics and Me chanics of Solids v ol. 43, pp. 13491362, 1995. [46] S. Geman and D. Geman, \Sto c hastic relaxation, Gibbs distributions, and the Ba y esian restoration of images," IEEE T r ans. Pattern A nalysis and Machine Intel ligenc e v ol. 6, no. 11, pp. 721741, 1984. [47] Los Alamos National Lab oratory GMV (The Gener al Mesh Viewer) h ttp://wwwxdiv.lanl.go v/X CM/gm v/GMVHome. h tml. [48] D. E. Goldb erg, Genetic A lgorithms in Se ar ch, Optimization and Machine L e arning Reading, MA. AddisonW esley 1989. [49] D. E. Goldb erg, K. Deb, and B. Korb, \Do not w orry b e messy ," Pr o c. 4th Int. Conf. Genetic A lgorithms San Mateo, CA, Morgan Kaufmann, pp. 2430, 1991. [50] G. H. Golub, M. Heath, and G. W ah ba, \Generalized crossv alidation as a metho d for c ho osing a go o d ridge parameter," T e chnometrics v ol. 21, pp. 215223, 1979. [51] G. H. Golub and U. V on Matt, \Generalized crossv alidation for large scale problems," Journal of Computational and Gr aphic al Statistics v ol 6, no. 1 pp. 134, 1997. [52] E. Hab er and D. Olden burg, \A GCV based metho d for nonlinear illp osed problems," Computational Ge oscienc es V ol. 4, No. 1, pp. 4163, 2000. [53] J. Hadamard, \Lectures on Cauc h y's Problem in Linear P artial Dieren tial Equations" Y ale University Pr ess 1923. [54] A. Hagemann, K. Rohr, H. S. Stiel, U. Sp etzger, and J. M. Gilsbac h, \Biomec hanical mo deling of the h uman head for ph ysically based, nonrigid image registration," IEEE T r ans. Me d. Imag. V ol. 18, pp. 875884, 1999. [55] P C. Hansen, \Analysis of discrete illp osed problems b y means of the Lcurv e," SIAM R ev. v ol. 34, pp. 561580, 1992. [56] E. J. Haug, K. K. Choi, and V. Komk o v, Design Sensitivity A nalysis of Structur al Systems Academic Press, 1986. 116 PAGE 127 [57] D. Haumann, \Mo deling the ph ysical b eha vior of rexible ob jects," T opics in Physic al lyBase d Mo deling, A. Barr et al. e ditors, A CM SIGGRAPH'87, Course Notes V ol. 17, A CM, 1987. [58] E. Reissner, \On v ariational principles in elasticit y ," Pr o c e e d. of Symp osia in Applie d Mathematics v ol. 8, pp. 16, MacGra wHill, 1958. [59] C. T. Hsiao, G. Chahine, and N. Gumero v, \Application of a h ybrid genetic/P o w ell algorithm and a b oundary elemen t metho d to electrical imp edance tomograph y ," Journal of Computational Physics v ol. 173, no. 2, pp. 433454, 2001. [60] M. Hori, and K. Oguni, \In v erse analysis metho d for iden tication of lo cal elastic prop erties b y using displacemen t data," Inverse Pr oblems in Engine ering Me chanics IV (editor, M. T anak a), Nagano, Japan, pp. 111119, 2003. [61] CH, Huang, \A nonlinear in v erse vibration problem of estimating the timedep enden t stiness co ecien ts b y conjugate gradien t metho d," International Journal for Numeric al Metho ds in Engine ering v ol. 50, no. 7 pp. 15451558, 2001. [62] T. S. Huang, \Mo deling, analysis, and visualization of nonrigid ob ject motion," Pr oc e e dings of 10th ICPR pp. 361364, 1990. [63] M. F. Insana and J. C. Bam b er (editors), Sp ecial issue on tissue motion and elasticit y imaging, Physics in Me dicine and Biolo gy v ol. 45, no. 6, pp. 14091714, 2000. [64] M. Y. Iscan, and R. P Helmer (Editors) F or ensic A nalysis of the Skul l: Cr aniofacial A nalysis, R e c onstruction, and Identic ation New Y ork, NY. WileyLiss, c1993. [65] L. Ji and J. McLaughlin, \Reco v ery of the Lam e parameter in biological tissues," Inverse Pr oblems v ol. 20, no. 1, pp. 124, 2004. [66] A. John, W. Kus, and P Oran tek, \Material co ecien ts iden tication of b one tissues using ev olutionary algorithms," Inverse Pr oblems in Engine ering Me chanics IV (editor, M. T anak a), Nagano, Japan, pp. 95101, 2003. [67] A. Joukhadar, T. Garat, and C. Laugier, \Constrain tbased iden tication of a dynamic mo del," Pr o c e e dings of Intel ligent R ob ots and Systems, IR OS'97 v ol. 1, pp. 337342, 1997. [68] F. Kallel and M. Bertrand, \Tissue elasticit y reconstruction using linear p erturbation metho d," IEEE T r ans. Me d. Imag. v ol. 15, no. 3, pp. 299313, 1996. [69] C. Kam bhamettu, D. B. Goldgof, D. T erzop oulos and T. S. Huang, \Nonrigid motion analysis," Handb o ok of PRIP: Computer Vision editor, Tza y Y oung, pp. 405430, Academic Press, San Diego, CA, 1994. [70] M. Kass, A. Witkin, and D. T erzop oulos, \Snak es: activ e con tour mo dels," International Journal of Computer Vision v ol. 1, no. 4, pp. 321331, 1988. 117 PAGE 128 [71] S. J. Kirkpatric k and D. D. Duncan, \Optical assessmen t of tissue mec hanics," Handb o ok of Optic al Biome dic al Diagnostics (editor, V alery T uc hin), published b y SPIE, 2002. [72] M. Kleib er and T. Hisada (editors), Design Sensitivity A nalysis A tlan ta T ec hnology Publication, 1992. [73] R. J. Knops and L. E. P a yne, Uniqueness The or ems in Line ar Elasticity SpringerV erlag, New Y ork, 1971. [74] R. M. Ko c h, M. H. Gross, F. R. Carls, D. F. v on Buren, G. F ankhauser, and Y. I. H. P arish, \Sim ulating facial surgery using nite elemen t mo dels," Pr o c. of SIGGRAPH's 96 pp. 421428, 1996. [75] S. Kumar, D. Goldgof, \Reco v ery of Global Nonrigid Motion a Mo del Based Approac h Without P oin t Corresp ondences", Journal of the Optic al So ciety of A meric a A 17(9), pp. 16171626, 2000. [76] L. D. Landau and E. M. Lifshitz, The ory of Elasticity 3rd edition, P ergamon Press, Oxford, 1986. [77] W. F. Larrab ee and J. A. Galt, \A nite elemen t mo del of skin deformation. (iii) the nite elemen t mo del," L aryngosc op e v ol. 96, pp. 413419, 1986. [78] J. Liu, \A m ultiresolution metho d for distributed parameter estimation," SIAM, J. Sci. Comput. v ol. 14, no. 2, pp. 389405, 1993. [79] R. L. Maurice, J. Oha y on, Y. F retign y M. Bertrand, G. Soulez, and G. Cloutier, \Nonin v asiv e v ascular elastograph y: theoretical framew ork," IEEE T r ansactions on Me dic al Imaging v ol. 23, no. 2, pp. 164180, 2004. [80] D. Metaxas, and D. T erzop oulos, \Shap e and nonrigid motion estimation through ph ysicsbased syn thesis," IEEE P AMI v ol. 15, no. 6, pp. 580591, 1993. [81] D. Metaxas, and D. T erzop oulos, \Constrained deformable sup erquadrics and nonrigid motion trac king," Pr o c e e dings of Computer Vision and Pattern R e c o gnition pp. 337343, June 1991. [82] D. Metaxas, Physicsb ase d Deformable Mo dels: Applic ations to Computer Vision, Gr aphics and Me dic al Imaging Klu w er Academic Publishers, 1997. [83] M. I. Miga, \A new approac h to elastograph y using m utual information and nite elemen ts," Phys. Me d. Biol. v ol. 48, pp. 467480, 2003. [84] K. Miller, \Least squares metho ds for illp osed problems with prescrib ed b ound," SIAM J. Math. A nal. v ol. 1, pp. 5274, 1970. [85] Z. Mic halewicz, \A surv ey of constrain t handling tec hniques in ev olutionary computation metho ds," Pr o c. of the 4th A nnual Confer enc e on Evolutionary Pr o gr amming (editors: J. R. McDonnell, R. G. Reynolds, and D. B. F ogel), pp. 135{155, MIT Press, 1995. 118 PAGE 129 [86] V. A. Morozo v, \On the solution of functional equations b y the metho d of regularization," Soviet Math. Dokl. v ol. 7, pp. 414417, 1966. [87] R. Muth upillai, D. J. Lomas, P J. Rossman, J. F. Greenleaf, A. Manduca, and R. L. Ehman, \Magnetic resonance elastograph y b y direct visualization of propagating acoustic strain w a v es," Scienc e v ol. 269, pp. 18541857, 1995. [88] G. Nak am ura and G. Uhlmann, \Iden tication of Lam e parameters b y b oundary measuremen ts," A meric an Journal Of Math. v ol. 115, pp. 11611187, 1993. [89] C. Nastar and N. Ay ac he, \F requencybased nonrigid motion analysis: Application to four dimensional medical images," IEEE T r ansaction on Pattern A nalysis and Machine Intel ligenc e v ol. 18, no. 11, pp. 10671079, 1996. [90] A. A. Ob erai, N. H. Gokhale1, and G. R. F eijo o, \Solution of in v erse problems in elasticit y imaging using the adjoin t metho d," Inverse Pr oblems v ol. 19, pp. 297313, 2003. [91] M. O'Donnell, A. R. Sk o v oro da, \Prosp ects for elasticit y reconstruction in the heart," IEEE T r ans. on Ultr asonics, F err o ele ctrics, and F r e quency Contr ol v ol. 51, no, 3, pp. 322328, 2004. [92] T. E. Oliphan t, A. Manduca, R. L. Ehman, J. F. Greenleaf, \Complexv alued stiness reconstruction for magnetic resonance elastograph y b y algebraic in v ersion of the dieren tial equation," Magnetic R esonanc e Me dcine v ol. 45, no. 2, pp. 299310, 2001. [93] R. Olmi, M. Bini, and S. Priori, \A genetic algorithm approac h to image reconstruction in electrical imp edance tomograph y ," IEEE T r ans. Evolutionary Computation v ol. 4, no. 1, pp. 8388, 2000. [94] J. Ophir, I. Cesp edes, H. P onnek an ti, Y. Y azdi, and X. Li, \Elastograph y: A quan titativ e metho d for measuring the elasticit y of biological tissues," Ultr asonic Imaging v ol. 13, pp. 111134, 1991. [95] J. Ophir, S. K. Alam, B. S. Garra, F. Kallel, E. E. Konofagou, T. Krousk op, C. R. B. Merritt, R. Righetti, R. Souc hon, S. Sriniv asan, and T. V arghese, \Elastograph y: Imaging the elastic prop erties of soft tissues with ultrasound," Journal of Me d. Ultr asound v ol. 29, pp. 155171, 2002. [96] P Oran tek, \Hybrid ev olutionary algorithms in optimization of structures under dynamical loads," Pr o c. of IUT AM Symp. on Evolutionary Metho ds in Me chanics, SMIA v ol. 117, Klu w er, 2004. [97] F. I. P ark e, and K. W aters, Computer F acial A nimation A.K. P eters, W ellesley Massac h usetts, 1997. [98] K. J. P ark er, L. Gao, R. M. Lerner, and S. F. Levinson, \T ec hniques for elastic imaging: a review," IEEE Engine ering in Me dicine and Biolo gy v ol. 15, pp. 5259, 1996. 119 PAGE 130 [99] K. D. P aulsen, M. I. Miga, F. E. Kennedy P J. Ho op es, A. Harto v, and D. W. Rob erts, \A computational mo del for trac king subsurface tissue deformation during sterotactic neurosurgery ," IEEE T r ans. Biome d. Eng. v ol. 46, no. 2, pp. 213225, 1999. [100] C. P ellotBarak at, F. F rouin, M. F. Insana, A. Hermen t, \Ultrasound elastograph y based on m ultiscale estimations of regularized displacemen t elds," IEEE T r ansactions on Me dic al Imaging v ol. 23, no. 2, pp. 153163, 2004. [101] A. P en tland and B. Horo witz, \Reco v ery of Nonrigid Motion and Structure," IEEE T r ansaction on Pattern A nalysis and Machine Intel ligenc e v ol. 13, no. 7, pp. 730742, 1991. [102] P P erona and J. Malik, \Scalespace and edge detection using anisotropic diusion," IEEE T r ansaction on Pattern A nalysis and Machine Intel ligenc e v ol. 12, no. 7, pp. 629639, 1990. [103] P Phillips, H. Mo on, S. Rizvi, and P Rauss, "The FERET Ev aluation Metho dology for F aceRecognition Algorithms" IEEE T r ans. on Pattern A na. and Machine Intel l. 22(10), pp. 10901104, 2000. [104] G. Pip erno, E. H. F rei, and M. Moshitzky \Breast cancer screening b y imp edance measuremen ts" F r ontiers Me d. Biol. Eng., 2(2), pp. 111117, 1990. [105] D. B. Plew es, J. Bishop, A. Samani, J. Sciarretta, \Visualization and quan tization of breast cancer biomec hanical prop erties with magnetic resonance elastograph y ," Physics in Me dicine and Biolo gy v ol. 45, no. 6, pp. 15911610, 2000. [106] T. P oggio, V. T orre and C. Ko c h, \Computational vision and regularization theory ," Natur e v ol. 317, pp. 314319, 1985. [107] P S. P o w ers, S. Sark ar, D. B. Goldgof, C. W. Cruse, and L. V. Tsap, \Scar Assessmen t: Curren t problems and future solutions," Journal of Burn Car e and R ehabilitation v ol. 20, no. 1, Jan. 1999. [108] K. R. Ragha v an and A. Y agle, \F orw ard and in v erse problems in imaging the elasticit y of soft tissue" IEEE T r ans. Nucle ar Scienc e 41(4), pp. 16391647, 1994. [109] H. Rhee, R. Aris and N. R. Am undson, First Or der Partial Dier ential Equations Pren tice Hall, Englew o o d Clis, New Jersey 1986. [110] T. P Runarsson and X. Y ao, \Sto c hastic ranking of constrained ev olutionary optimization," IEEE T r ansactions on Evolutionary Computation v ol. 4, no. 3, pp. 284294, 2000. [111] J. R. Sae, B. Da vis, and P Williams, \Recen t outcomes in the treatmen t of burn injury in the United States: A rep ort from the American Burn Asso ciation P atien t Registry ," Journal of Burn Car e and R ehabilitation v ol. 16, pp. 219232, Ma y 1995. [112] R. Salomon, \Ev olutionary algorithms and gradien t searc h: Similarities and dierences," IEEE T r ansactions on Evolutionary Computation v ol. 2, no. 2, pp. 4555, 1998. 120 PAGE 131 [113] A. Samani, J. Bishop, and D. B. Plew es, \A constrained mo dulus reconstruction tec hnique for breast cancer assessmen t," IEEE T r ans. Me d. Imag. v ol. 20, no. 9, pp. 877885, 2001. [114] C. Sc hmid, R. Mohr, and C. Bauc khage, \Ev aluation of in terest p oin t detectors," International Journal of Computer Vision v ol. 37, no. 2, pp. 151172, 2000. [115] C. G. Shaefer, \The AR GOT strategy: Adaptiv e represen tation genetic optimizer tec hnique," Pr o c. 2nd Int. Conf. Genetic A lgorithms and Their Applic ations Hillsdale, NJ, La wrence Erlbaum, pp. 5055, 1987. [116] M. Shah and R. Jain, MotionBase d R e c o gnition Academic Publishers, Netherlands, 1997. [117] J. Shi and C. T omasi, \Go o d features to trac k," IEEE Confer enc e on Computer Vision and Pattern R e c o gnition (CVPR), Se attle pp. 593600, June 1994. [118] D. Sh ulman and J. Y. Aloimonos, \Nonrigid motion in terpretation: a regularized approac h," Pr o c. R oy. So c. L ondon v ol. B233, pp. 217234, 1988. [119] A. R. Sk o v oro da, S. Y. Emeliano v, M. A. Lubinski, A. P Sarv azy an, and M. O'Donnell, \Theoretical analysis and v erication of ultrasound displacemen t and strain imaging," IEEE T r ans. Ultr asonics, F err o ele ctrics, and F r e quency Contr ol v ol. 41, no. 3, pp. 302313, 1994. [120] I. S. Sok olnik o, Mathematic al The ory of Elasticity 2nd ed. McGra wHill Bo ok Co., New Y ork, 1956. [121] T. Sulliv an, J. Smith, J. Kermo de, E. McIv er, and D. J. Courtemanc he, \Rating the burn scar," Journal of Burn Car e R ehabilatation v ol. 11, pp. 256260, 1990. [122] C. Sumi and K Nak a y ama, \A robust n umerical solution to reconstruct a globally relativ e shear mo dulus distribution from strain measuremen ts" IEEE T r ans. Me d. Imag. 17(3), pp. 419428, 1998. [123] C. T anner, A. Degenhard, J. A. Sc hnab el, A. CastellanoSmith, C. Ha y es, L. I. Sono da, M. O. Leac h, D. R. Hose, D. L. G. Hill, and D. J. Ha wk es, \A metho d for the comparison of biomec hanical breast mo dels," IEEE Workshop on Mathematic al Metho ds in Biome dic al Image A nalysis Kauai, USA, pp. 1118, 2001. [124] A. T aran tola, Inverse Pr oblem The ory: Metho ds for Fitting and Mo del Par ameter Estimation New Y ork, Elsevier, 1987. [125] D. T erzop oulos, and A. Witkin, \Ph ysically based mo dels with rigid and deformable comp onen ts," IEEE Computer Gr aphics and Applic ations 8(6), pp. 4151, 1988. [126] D. T erzop oulos, and D. Metaxas, \Dynamic 3D mo dels with lo cal and global deformations: deformable sup erquadrics," IEEE T r ans. Pattern A nal. and Machine Intel l. 13(7), pp. 703714, 1991. 121 PAGE 132 [127] D. T erzop oulos, and K. W aters, \Analysis and syn thesis of facial image sequences using ph ysical and anatomical mo dels," IEEE T r ans. Pattern A nalysis and Machine Intel ligenc e 15(6), pp. 569579, 1993. [128] D. T erzop oulos, \Regularization of visual problems in v olving discon tin uities," IEEE T r ans. Pattern A nalysis and Machine Intel ligenc e v ol. 8, no. 4, pp. 413424, 1986. [129] D. T erzop oulos, \The computation of visiblesurface represen tations," IEEE T r ans. Pattern A nalysis and Machine Intel ligenc e v ol. 10, no. 4, pp. 417438, 1988. [130] J. G. Thac k er, \The elastic prop erties of h uman skin in vivo ," Ph.D. Dissertation, Univ ersit y of Virginia, Charlottesville, V A, 1976. [131] A. N. Tikhono v, A. Gonc harsky V. Stepano v and A. Y agola, Numeric al Metho ds for the Solution of Il lPose d Pr oblems Klu w er, Dordrec h t, 1995. [132] S. P Timoshenk o and J. N. Go o dier, The ory of Elasticity 3rd ed., McGra wHill, New Y ork, 1970. [133] C. T ruesdell and W. Noll, \The nonlinear eld theories of mec hanics", In: Handbuch der Physik (ed. S. Fl ugge), v ol. I I I/3, SpringerV erlag, Berlin, 1965. [134] R. Y. Tsai, and T. S. Huang, \Uniqueness and estimation of threedimensional motion parameters of rigid ob jects with curv ed surfaces," IEEE T r ansaction on Pattern A nalysis and Machine Intel ligenc e V ol. 6, pp. 1326, 1984. [135] L. V. Tsap, D. B. Goldgof, S. Sark ar, and P S. P o w ers, \A visionbased tec hnique for ob jectiv e assessmen t of burn scars," IEEE T r ans. Me d. Imag. v ol. 17, no. 4, pp. 620633, 1998. [136] L. V. Tsap, D. B. Goldgof and S. Sark ar, \Nonrigid motion analysis based on dynamic renemen t of nite elemen t mo dels" IEEE T r ans. Pattern A nal. and Machine Intel l. 22(5), pp. 526543, 2000. [137] M. T urk, and A. P en tland, \Eigenfaces for recognition," Journal of Co gnitive Neur oscienc e (3)1, pp. 7186, 1991. [138] G. Uhlmann, \Dev elopmen ts in in v erse problem since Calderon's foundational pap er," Harmonic A nalysis and PDE Chicago, IL: Univ ersit y of Chicago Press. [139] E. E. W. V an Houten, K. D. P aulsen, M. I. Miga, F. E. Kennedy and J. B. W ea v er, \An o v erlapping subzone tec hnique for MRbased elastic prop ert y reconstruction," Magnetic R esonanc e in Me dicine v ol. 42, pp. 779786, 1999. [140] G. W ah ba, \Spline mo dels for observ ational data," SIAM CBMSNSF Regional Conference Series in Applied Mathematics, V. 59, Philadelphia, 1990. [141] J. W eic k ert, A nisotr opic Diusion in Image Pr o c essing T eubnerV erlag, Stuttgart, 1998. 122 PAGE 133 [142] Y. Y aco ob, H. Lam, and L. S. Da vis, \Recognizing faces sho wing expressions" International Workshop on A utomatic F ac e and Gestur e R e c o gnition Zuric h, pp. 278283, 1995. [143] Y. Y aco ob, and L. S. Da vis, \Smiling faces are b etter for face recognition" International Confer enc e on F ac e R e c o gnition and Gestur e A nalysis W ashingtonDC, pp. 5257, 2002. [144] F. Y eung, S. F. Levinson, D. F u, and K. J. P ark er, \F eatureAdaptiv e Motion T rac king of Ultrasound Image Sequences Using A deformable Mesh", IEEE T r ans. Me d. Imag. 17(6), pp. 945956, 1998. [145] E. K. Y eong, R. Mann, L. H. Engra v, M. Goldb erg, V. Cain, B. Costa, M. Mo ore, D. Nak am ura, J. Lee, \Impro v ed burn scar assessmen t with use of a new scarrating scale," Journal of Burn Car e and R ehabilitation v ol. 18, no. 4, pp. 353355, 1997. [146] Y. Zhang, D. B. Goldgof, S. Sark ar, and L. V. Tsap, \Mo delbased nonrigid motion analysis using natural feature adaptiv e mesh," Pr o c e e dings of International Confer enc e on Pattern R e c o gnition pp. 839843, 2000. [147] Y. Zhang, S. J. Kundu, D. B. Goldgof, S. Sark ar, and L. V. Tsap, \Elastic face, An anatom ybased biometrics b ey ond visible cue" 17th International Confer enc e on Pattern R e c o gnition, V ol. 2 Cam bridge, UK, pp. 1922, 2004. [148] Y. Zhang, D. B. Goldgof, and S. Sark ar, \T o w ards ph ysicallysound registration using ob jectsp ecic prop erties for regularization," 2nd Inter. Workshop on Biome d. Image R e gistr ation pp. 358366, 2003. [149] Y. Zh u, T. J. Hall, and J. Jiang, \A niteelemen t approac h for Y oung's mo dulus reconstruction," IEEE T r ansactions on Me dic al Imaging v ol. 22, no. 7, pp. 890901, 2003. [150] S. C. Zh u and D. Mumford, \Prior learning and Gibbs reactiondiusion" IEEE T r ans. Pattern A nalysis and Machine Intel ligenc e v ol, 19, no. 11, pp. 123650, 1997. [151] O. C. Zienkiewicz, The Finite Element Metho d 3rd edition, McGra wHill, 1977. 123 PAGE 134 ABOUT THE A UTHOR Y ong Zhang receiv ed B.S. degree in Marine Science from Ocean Univ ersit y of Qingdao, China in 1986 and M.S. degree in Marine Science from the same sc ho ol in 1989. He also receiv ed M.S. degree in Computer Science and Engineering from Univ ersit y of South Florida in 2000. He is curren tly a PhD candidate at Univ ersit y of South Florida. His researc h in terests include computer vision, motion mo deling, image pro cessing, pattern recognition and articial in telligence. 