UFDC Home  myUFDC Home  Help 



Full Text  
GUIDED SURGERY USING RAPID PROTOTYPING PATIENTSPECIFIC GUIDES: A METHODOLOGY TO QUANTIFY MECHANICAL STABILITY AND UNIQUENESS OF FIT By BARBARA GARITA A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2007 2007 Barbara Garita To Donato, Noemi, Esteban, Pablo, Ifi, Ignacio and Natalia. ACKNOWLEDGMENTS First, I must thank my committee for taking the time to work with me during this endeavor. I must also express gratitude to my advisor Dr. Frank Bova for persevering with me throughout the time it took me to complete this research and write this dissertation. Dr. Friedman, Dr. Gilland, Dr. Kim, and Dr. Vemuri have also been generous enough to dedicate their time and expertise to improve my work; I am grateful for their contributions and their kind support. I am grateful as well to Dr. Didier Rajon, for his timely appraisals of my work and useful suggestions throughout these years. I must express my profound gratitude to Dr. Jonathan Earle for invaluable support and friendship, especially during the last year of my dissertation. Jeff Citty has been an immense support as well; working with him and the EFTP program for the last 4 years has given me a true sense of belonging while at the University of Florida. I need to express my gratitude and deep appreciation to Serggio who always stood next to me in good times and bad times. I must acknowledge, as well the many friends and colleagues who assisted, advised and supported my research and writing efforts over the years, my gratitude goes to Ashley, Mafer, Ryanne, April, Luis, Mike, Pam, Atchar and Rolando. Finally I sincerely thank my parents and siblings for encouraging my decision to stay away from home for so many years in order to accomplish these academic dreams. I find the time spent separated from them irreplaceable and I hope I will get the chance to make it all up in the near future. TABLE OF CONTENTS page A C K N O W L E D G M E N T S ..............................................................................................................4 L IS T O F T A B L E S ................................................................................. 7 LIST O F FIG U RE S ................................................................. 8 ABSTRACT ........................................... .. ......... ........... 12 CHAPTER 1 IN T R O D U C T IO N ....................................................................................... .................... 14 2 BACKGROUND AND SIGNIFICANCE............................................................. ........... 19 Im ageGuided Surgery .................. ................ ............... ........ .. .......... .............. 19 FrameBase and Frameless Stereotactic Approaches..........................................20 M odel B asked G uidance........... .... ..... .......... ...... .............. ................ ............. 22 Customized positioning frames for pedicle screw placement..............................22 Customized position frames for intracranial surgery ............................................24 S ig n ifican ce .................................................................................2 7 F racial L andm arks and A ngles ........................................................................ .................. 28 The Finite Elem ent M ethod and ITK ........................................................... ............... 29 Elem ents U sed in the Sim ulation ............................................................................. 30 T h e sp a c e tru ss ................................................................................................... 3 0 The linear tetrahedron ....................... ................ ................... ........ 34 Characteristics of Stiffness M atrices .................................................................... ....35 L in ear S o lv ers................................. ....................................................... ............... 3 5 3 M A T E R IA L S ........................................................................................................... .. 4 7 4 METHODOLOGY FOR MEASURING STABILITY............................... 54 Introduction ................................................................................. 54 System Components and FEM Setup .......................................................................55 E external F forces ............................................................................................... ....... 56 R ig id L ay e r ..........................................................................................5 6 B e d o f S p rin g s ........................................................................................................... 5 8 Program Parameters ............................................................................... 59 In p u t P aram eters ....................................................................................................5 9 S p rin g h eig h t .............................................................5 9 E xternal force m agnitu de ................................................................................... 59 Rigid layer height ..................................... .............. .....................60 D ifferen ce in elasticity ....................................................................................... 6 0 O u tp u t P aram eters ..................................................................................................... 6 0 C criteria ........ ............................................................... 62 T N T R atio .....................................................6 3 Single Point C riterion .................. ................................ .. .. ........ .............. 64 TNT threshold ........................... ................ 65 Consistency of TN T threshold ........................................ ........................... 67 Single Force Criterion .................................... ..... .......... .............. .. 67 5 PROGRAM DESCRIPTION AND VALIDATION.................................... 82 P ro g ram D e scrip tio n ................................................... ............................................. .. 82 V alidation of F E M elem ents................................................................................. .... ..... 85 Linear Spring V alidation ........................................................... ............... 85 Tetrahedron Element Validation .................................. .....................................86 A Sim ple M odel V alidation.................................................. ............................... 86 6 RESULTS AND D ISCU SSION .................................................. ............................... 93 7 CONCLUSION AND FUTURE WORK .................. ........... .....................109 APPENDIX A PRINCIPLE OF VIRTUAL WORK AND ITK FEM MODULE.................... ......... 111 FEM Term inology for Truss E lem ents ................................................................................ Principle of V irtual W ork ................................................................................... 112 ITK FEM and Principles of Virtual W ork ...................................................... .......... 113 D isplacem ent F unctions ................................................................ ... .....................115 Sh ap e F u n action s ...................1...................1...................5......... Strain D isplacem ent Equation ......................................................... .............. 116 L ocal E lem ent Stiffness M atrix ........................................................................... ... 117 G auss Q uadrature ................ ...................................... ........... ......... 118 ITK FEM M odule ............... ................. ........... ................... ....... ... 119 L IST O F R E F E R E N C E S ............................................................................. ..........................123 B IO G R A PH IC A L SK E T C H ......................................................................... .. ...................... 127 6 LIST OF TABLES Table page 41 Relationship between the tetrahedral height and a measure of the rigidity of the rigid lay er ................... ............................................................. ................ 7 0 42 Table that shows the relationship between the difference in elastic modulus between the rigid layer and the bed of springs and the rigidity measure and Cholmod's rcond value......................................................... 70 43 Single point criteria .................. .................. .................. ........... ......... .... 70 51 Displacement in the X, Y and Z direction from the example in Figure 54 tested with the ITKFEM module and ANSYS................................. ................................... 87 61 Measurement of the facial angles for the three anatomical models..............................99 62 Instability score, presented as the percentage of forces that resulted in instability and the total forces applied on a fram e .................................. ............... ............... 99 LIST OF FIGURES Figure pe 11 Components involved in this investigation and the prototypes for the new technique for performing imageguided surgery with patientspecific reference guides................... 18 21 Pretarget localization and presurgical planning is done in a computer workstation.......36 22 Procedures related to framebased image guidance surgery.................... ..................37 23 Frameless stereotactic surgery uses fiducials that are glued to the patient's skin.............37 24 Typical layout when frameless stereotactic surgery is employed.............................. 38 25 Tem plate designs for the spine ................................................ .............................. 38 26 R outer guides for craniotom y ................................................. ............................... 39 27 Biopsy guides ............................................................... ..... ...... ........ 40 28 R apid prototyping m machines ........................................... .................. ............... 41 29 Chart to visualize the significance of the present work ................. ............... ..............42 210 Facial landmarks ........................ ........ .. .... .... ........ .... .... 43 211 Facial angles ................................................... 44 212 Space truss element has six displacements to establish its deformed position ................45 213 Application of a force perpendicular to the axis of the space truss ..............................45 214 Space truss in local coordinates. ............................................................. .....................46 215 The linear tetrahedron has a unique numbering scheme .............................................46 31 The three models used to evaluate the methodology for testing unique fit .....................49 32 Masks built for testing the methodology of fit............ .................. ..................50 33 C ap lik e referen ce fram e......................................................................... .....................5 1 34 Virtually created half sphere test m odel ........................................ ........................ 51 35 Placement of pins for a simple test on the real stability of the large and small masks......52 36 Coronal view of placement of the pins onto the real masks ............................................52 37 The method employed to demonstrate the difference in stability between the large and sm all m asks. ...........................................................................53 41 The methodology for stability of fit involves simulating a rigid layer onto of a deform able bed of springs .......................................................... .. ............... 71 42 Displacement of a mask as it responds to external forces. ..............................................72 43 Triangulation of a customized positioning device proper to STL format........................72 44 D direction of the N main ...................................... ......... ....................... ......... 73 45 N odes point on the rigid fram e ................................................ .............................. 73 46 In a wedge three tetrahedral elem ents are fitted ..................................... .................74 47 Schematic showing how the maximum rigidity fault is calculated .................................75 48 The behavior of a linear spring element under the application of an axial force ..............75 49 Schematic that demonstrates how each spring's displacement was evaluated by identifying the normal (ni) and tangential (ti) displacement magnitudes...........................76 410 Schematic that shows why it is important to distinguish between the springs that are mostly in compression or tension and those that are basically slipping or have a high T N T ra tio .............................................................................. 7 7 411 Bounding volume for a spring to be in compression ......................................................77 412 Twodimensional interpretation of the single point criteria ............................................78 413 Schematic to demonstrate the importance of evaluating the percentage slip of TNT ratio for each spring ......................... ......... .. ........... .. ............78 414 Schematic to shows at the level at which the TNT threshold was defined......................79 415 Empirical determination that the TNT threshold greater than 0.6 correlates with the training data, which says that the large masks are stable .............................................79 416 Empirical determination that the TNT threshold set smaller than 0.8 correlates with the training data, which says that the small masks are unstable.....................................80 417 A case of the large Asian mask, which shows the consistency of an upper TNT thresh old of 0 .7 ............................................................................ 80 418 Schematic showing the relation between the volume of the entire frame (fvol) and the com pression patch volume e (cvol) ......................................................... ............... 81 419 B ox ratio interpretation. ......................... ....................... .......... .... .... 1 51 Flow chart of program for testing the uniqueness of fit.............................................88 52 Original surface is extended in the positive normal direction to form the top surface and the points extended in the negative normal direction to become boundary co n d itio n p o in ts...................................................................... 9 0 53 Examples of al linear spring element or space truss tested in all 3dimensions ...............90 54 Example of two tetrahedral elements used in the validation of the tetrahedral elem ents used in the program ................................................. ............................... 91 55 A corer frame that was used to validate the program ............................................... 92 61 Placem ent of Caucasian m asks onto solid m odel ........................................ .................100 62 Placement of Asian masks onto solid model .............. ..........................................100 63 Placement of infant masks onto solid model .............. ..........................................101 64 Placement of a caplike reference frame onto an upper head model............................101 65 Surface plot color coded by the maximum full rotation angle for the Caucasian m a sk s ................... ........................................................... ................ 10 2 66 Surface plot color coded by the maximum full rotation angle for the Asian masks........102 67 Surface plot color coded by the maximum tilt angle for the infant masks ....................103 68 Surface plot color coded by maximum tilt angle for the cap reference frame design .....103 69 Surface plots color coded by maximum tilt angle for the half sphere test.......................104 610 Normalized histogram of box ratios for each tilt angle (0) for the Caucasian large mask..................... ...................................... 104 611 Normalized histogram of box ratios for each tilt angle (0) for the Caucasian small mask..................... ...................................... 105 612 Normalized histogram of box ratios for each tilt angle (0) for the Asian large mask. ....105 613 Normalized histogram of box ratios for each tilt angle (0) for the Asian small mask.....106 614 Normalized histogram of box ratios for each tilt angle (0) for the infant large mask. ....106 615 Normalized histogram of box ratios for each tilt angle (0) for the infant small mask.....107 616 Normalized histogram of box ratios for each tilt angle (0) for the cap reference frame d e sig n ........ ........ .......................................................................... 10 7 617 Normalized histogram of box ratios for each tilt angle (0) for the half sphere test........ 108 618 P profiles of solid m models ......................................................................... ....................108 A1 Differential element subjected to real stresses and in equilibrium. ...............................121 A2 The greater part of the ITK FEM library ................................................................... 122 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy GUIDED SURGERY USING RAPID PROTOTYPING PATIENTSPECIFIC GUIDES: A METHODOLOGY TO QUANTIFY MECHANICAL STABILITY AND UNIQUENESS OF FIT By Barbara Garita August 2007 Chair: Frank Bova Major: Biomedical Engineering Rapid prototyping technology (3dimensional printing) in conjunction with the virtual manipulation of computergenerated 3D anatomical models rendered from diagnostic images (CT and MRI) is a novel technique for directing imageguided surgery. By incorporating a graphical user interface, the surgeon can plan an intracranial surgical procedure by fabricating a patientspecific reference frame (masklike facial frame) through the use of rapid prototyping technology. Early testing of these patientspecific frames revealed the necessity to include a step prior to fabrication to analyze the design of the frame for stability of fit. In response to this need, the objective in the present investigation is to develop a method that numerically assesses the frame's stability. The approach was to create a finiteelement simulation that models the patient specific reference frame as a rigid layer composed of linear tetrahedron elements under a deformable material simulated by a bed of linear spring elements. The program analyzes the displacement of the rigid layer by simulating external forces distributed around its perimeter. The stability of fit is characterized by the response of the linear springs. Three patients with different facial anatomical features were selected as training data for testing the simulation program. For each case two patientspecific reference frames, one that involves more surface area and surface variations such as contours of the nose bridge and another that consists of less surface area and less surface variations, were built with Rapid prototyping based on the patient's 3D anatomical facial models. In all three cases the methodology correlated a stability score with the stability of the real reference frames. In addition two more cases were evaluated; the first did not incorporate the facial structures and was design to hold onto the top of the head in a caplike manner, the second case was a virtually created half sphere. In both cases it was expected to find no a unique fit, and correspondingly the simulation found both cases to be unstable. CHAPTER 1 INTRODUCTION Starting in 2003, researchers in the radiosurgery and biology lab (RSB) at the University of Florida's McKnight Brain Institute developed the capability to fabricate a patientspecific reference frame that has a potential of becoming a new type of imageguided surgery (IGS) technique [33]. The fabrication of the patientspecific reference frame involves coding a series of algorithms that render a detailed 3dimensional lifelike anatomical model from computer tomography (CT) or magnetic resonance imaging (MRI) [6] datasets (Figure 11). Onto a patientspecific 3dimensional anatomical model a surface is selected (virtually "painted") from which the patientspecific reference frames are created. The ability to create, visualize, and manipulate 3dimensional virtual anatomical models on a computer monitor heightened the enormous potential of utilizing medical image datasets to create customized reference frames for medical diagnostics and treatment, specifically in the area of image guidance surgery [36]. The objective the patientspecific reference frame is to replace existing technology of framedbased (fixing with pins a frame to the skull) and/or frameless (using fiducials and extra tracking equipment in the operation room) stereotactic systems [6]. Stereotactic techniques are used in neurological research and/or surgery to direct a surgical instrument (needle or electrode) into a specific location in the brain [1]. The patientspecific reference frames (or customized guiding template) are designed to tightly fit to a specific patient's facial anatomy contour [6] and allow guidance for neurosurgical applications such as craniotomies, biopsies, and deep brain stimulation [34]. Several clinical trials were initiated to test this new imageguided surgical approach. One of the first involved developing a customized guiding template that fitted onto vertebral bodies and served as a guide for pedicle screw placement. As the templates for vertebral bodies were tested (through a nonsignificant risk IRB approved clinical trial) on actual patients in the operation room (OR), the research team quickly identified a problem that was defined within the context of biomedical engineering: the computer selected surface for the reference frame could not accurately predict the surface that the surgeon was able to present (because of the residual soft tissue), therefore "no unique surface" for the template could be identified. Eventually, the idea of using customized guiding templates for the spine was abandoned and the development of this approach for cranial procedures became the primary focus of the design application. The idea of using customized guiding templates for stereotactic surgery proved to have significant potential [33], but the problem of stability of fit was still apparent. The problem of the stability of fit for the reference frames became of great importance because of the need to add attachments to the reference frame for particular intracranial guidance surgeries (Figure 11D and 11E). As the distance between the reference frame and the attachments increases, the lever arm increases and small movements on the reference frame will results in exaggerated moments on the attachments tips (Figure 11D). For this reason, the accuracy of the patientspecific reference frames depends upon the identification of a broad and conforming patient surface. This surface must allow a reference frame mask (negative from the anatomical model's surface) to be designed in a manner that prevents both rotations and translations. In other words, the mask must fit tightly onto the patient's surface and not slide or wiggle. Hence, the initiative of the present investigation stems from the necessity to solve this stability problem by developing a numerical methodology to evaluate the patientspecific frame's design by testing its stability of fit prior to fabrication with Rapid prototyping technology. The methodology developed for testing stability of fit relies on the creation of a finite element simulation program that models contact interactions between patientspecific frame and the patient's anatomical face model (Figure 11). The approach used is the idealization of a contact mechanics problem by means of the simulating the patientspecific reference frame as a rigid layer composed of linear tetrahedron elements under a deformable material simulated by a bed of linear springs (or truss) elements. The approach will simulate the movement of the rigid tetrahedral layer by the application of external forces onto the frame's surface. As training data for the methodology, three solid lifelike anatomical face models of three patients with notably dissimilar facial profiles were created with rapid prototyping and for each patients two frames one the is known to be stable and the other one known to be unstable were also created. The methodology for testing stability of the frames involves two criteria, both based on the resulting simulated movement of the frames onto the patients' 3dimensional anatomical models. The first criterion is concern with identifying locations of excessive movement and the second criterion provides with a yes or no answer to the stability of the frame on a per force basis. It will be shown that these two criteria are able to predict the degree of instability of the selected surface mask construct; by identifying when the construct is able to resist translations and rotations to typical external force directions. The methodology to test uniqueness of fit presented in this investigation is an essential component to the new technique for imageguided surgery that involves rapid prototyping. The design of the patientspecific reference frame has the imperative requirement to have a unique fit onto the patient's anatomy, this is the major concern of this new imageguided surgery technology [33]. Therefore the work presented here, which numerically evaluates the mechanical stability and uniqueness of fit of the patientspecific reference frame, will overcome this major problem and will allow a surgeon to use the new technology with more confidence and accuracy. In this way, the new IGS technology is a step closer towards being distributed to other the medical communities. In addition, the methodology has been written in the same language and uses the same open source libraries as the main graphical user interface, from which the patient specific reference frames are designed, so it can be easily adapted as a module component. A rr C D E Figure 11. Components involved in this investigation and the prototypes for the new technique for performing imageguided surgery with patientspecific reference guides. A) The graphical user interface that can create 3dimensional patient anatomical model from a series of 2dimensional medical image datasets, and the selection of a contour surface (in red) for the fabrication of a patientspecific reference frames. B) Physical patient's anatomical model created with Rapid prototyping for evaluating the stability of fit. C) Patientspecific reference frames that conformingly fits onto the patient's model. D) Virtual model and reference frame with attachments to guide a biopsy probe. As the moment arm increases, small rotations and translation on the mask will have an exaggerated effect at the level of the attachment tip. E) Physical reference frame and attachments to guide a biopsy probe. CHAPTER 2 BACKGROUND AND SIGNIFICANCE ImageGuided Surgery During the early 20th century, Sir Victor Horsley recognized the need for targeting specific diseased regions in the brain without injuring neighboring healthy structures. He hired an engineer to provide with a method to solve the problem; the engineer came up with the idea of using three imaginary orthogonal spatial planes (horizontal, frontal and sagittal) [23] to specify any location in the brain. Later, during the mid 20th century, Spiegel and Wycis et al. [41] introduced the first stereotactic application. As it was originally envisioned, the stereotactic approach relates to movement in 3dimensional space; therefore, stereotactic surgery uses a 3 dimensional coordinate system to locate small targets inside the brain. Spiegel and Wycis et al. [44], began the work of locating specific regions in the brain using intrinsic reference points and in 1952 [43] they published the first stereotactic map of the brain using the pineal body as the only reference point. However, this localization technique proved unreliable because of anatomical differences between patients and as the surgeons moved further from their reference point they lost accuracy. Seeking to improve the reliability of their technique, Spiegel and Wycis et al. developed the approach of fitting a custom cap to the patient and attaching a head ring to this cap; in this way they established a relation between the brain atlas and a reference head ring [43]. Besides the Spiegel and Wycis custom cap apparatus, during the following year several others investigators developed stereotactic devices for different purposes and with different usage methods [18]. Some systems were based on Cartesian coordinates, others used polar coordinates for placing electrodes at target locations; while some designs moved the patient's head towards a target location; others used interlocking arcs to move the apparatus to get to a specific location [18]. With the advent of computer tomography (CT), imaging of the brain became readily available. The BrownRobertsWells (BRW) system was the first to be designed particularly for use with CT images and was the first to implement a computer base algorithm to calculate coordinates and a trajectory to a target point [22]. As a result, a specific patient's anatomical information obtained through CT imaging was, for the first time, properly visualized and referenced to a head ring [36]. Later, in 1987, Schad et al. [36] introduced multimodality image registration, which increased the accuracy and tissue differentiation for stereotactic targeting. FrameBase and Frameless Stereotactic Approaches Today image guidance procedures can be divided in two categories, according to different referencing methods employed. These are generally termed framebase and frameless stereotactic approaches. The purpose of both is to provide the surgeon with an understanding of the anatomy beneath the visible operative surface. State of the art stereotactic procedures allow the surgeon to physically relate spatial information from a set of diagnostic images (represented by a patient's 3dimensional virtual model), onto a patient's actual, realworld anatomy in a reliable and accurate manner. In general, this is accomplished by means of a discernible reference frame (or landmarks) that are integrated in both the virtual patient's model and the realworld patient and correlates the information amongst them. The creation of a diagnostic plan, for both framebase and frameless stereotactic surgeries, is performed on the virtual, computer generated patient, where image datasets (MR and CT) depict the anatomical information together with the rigid frame landmarks for framebased or skin landmarks for frameless stereotactic surgeries (Figure 21) [6]. Framebase stereotactic surgery is the most accurate technique available [32], but it is not free of significant disadvantages. This represents the most accurate method because it rigidly fixes a reference frame onto the patient's skull, through the use of minimally invasive pins. For example, in the case of a stereotactic framebased biopsy a head ring is applied to the patient prior to imaging (Figure 22A). The diagnostic scan will allow the surgeon to visualize the intracranial anatomy and create a virtual surgical plan on a computer workstation. With a graphical user interface the surgeon will define the intracranial target location and the computer system will mathematically calculates the corresponding coordinates relating the virtual plan to the real world patient anatomy (Figure 21). The coordinates are set on the physical stereotactic system (Figure 22B) (which is directly fixed onto the rigid frame), providing a guide or trajectory for the surgeon [8, 35]. The main shortcomings of framebase stereotactic surgery include the extend of time that the procedure takes (from to presurgical frame application to surgery completion), the specific expertise required for frame application, the high cost of replacement equipment, and the discomfort that the head ring causes on the patient. Indeed, the rigid frame is a discomfort for the patient at the time of imaging and during the preparation of the sterile field, and it gets in the way of the surgical procedures in general. Nevertheless, framebase imageguided surgery provides the gold standard for precision and accuracy for intracranial stereotactic procedures. Frameless stereotactic surgery has been presented as an alternative to framebased procedures. The main difference between framebased and frameless stereotactic systems, is that instead of a head ring, fiducials (landmarks) are used (Figure 23). For cranial surgery these are either MR and CT identifiable markers that are temporally glued onto the patient's skin [2, 51] or anatomical landmarks are defined on 3dimensional rendered images base upon the patient's CT or MR datasets [39]. In both cases, the procedures do not have to be completed in an extended and continuous span of time. Both framebased and frameless procedures permit the same virtual computer planning to be performed by the surgeon; however frameless stereotactic surgery requires instrumenttracking equipment within the operating room (OR) [51]. The main drawback of the frameless approach is the necessity for registration of the fiducials with a reference frame used to track movement between the patient and the surgical tools [6]. The registration requires special cameras (electromagnetic) and additional computer workstations to be present in an already crowded operating room (Figure 24). Although electromagnetic tracking allows for nonline of sight tracking, the surgical team is still burdened with the additional computers in the OR [6], special traceable instruments, the need to minimize ferrous metals from with the operative field and the registration procedures. Model Based Guidance Customized positioning frames for pedicle screw placement As mentioned in the introduction, researchers in the RSB lab initially built customized guiding templates to fit onto vertebral bodies for the placement of pedicle screws. It was during the testing of these spine templates that the problem of uniqueness of fit was first addressed. The following section provides background information on the first implementations of custom guiding templates for the spine and also provides additional information to aid in better understanding of the problem at hand. Stereotactic guidance with the aid of customized guiding templates initiated in 1999 for the placement of pedicle screws [20]. In 1999, Goffin et al. [20] first published a paper that explained about a guide for the posterior transarticular screw fixation of C1C2 (cervical 1 to cervical 2); and in 2001, Goffin et al. [21], Bimbaum et al. [6] and Yoo et al. [50] published papers referring again to template guides for placement of pedicle screws in the spine (Figure 25 AC). While Goffin et al. [20] developed templates for the fixation of the cervical spine, Birnbaum et al. [6] and Yoo et al. [50] tested particularly on the lumbar spine of cadavers and humans. All three research groups used rendered 3dimesional CT information to create the template guides based on the bony surface of the spine. Goffin et al. [20] printed the 3D solid models of the cervical spine with stereolithography and physically defined the correct placement of the pedicle screws, and then he translated the trajectory to a CAD (Computer Aided Design) program. Goffin et al. [21] and Birnbaum et al. [6] used computer controlled milling machines or NC (numerically controlled) machines to create the template guides; Yoo et al. used FDM (Fused Deposition modeling), a machine from Stratasys (Stratasys, Tempe, AZ), that follows the same principle of stereolithography but is composed of a movable nozzles where melted plastic forms the template one layer at a time. Goffin et al. [21] manufactured the device on an acrylate resin and the template had a horseshoe appearance that was positioned on the posterior aspect of C2 (Figure 25 A). He tested on some cadavers and two patients. On the cadavers he initially experienced rotational instability of the template and increased the contact area to include more vertebral support. On one of the two patient cases he was successful; while on the other the template guide did not perform well. Birnbaum et al. [6] created polycarbonate templates with NC machines (Figure 2 5 B), and conducted a cadaveric study comparing computerassisted surgery with individual templates on thirteen lumbar spines. Under imageguided technique, he found two misplaced screws; while all were correctly placed when using the template guides. Yoo's templates were constructed from ABS plastic and expanded over a large contact area comprising both transverse processes, the laminae areas and the spinous process [50]. He tested on dry dissected vertebrae. In all cases the positions of the screws were evaluated with postoperative CT scans. A later paper by Van Cleynenbreugel et al. [47], published in 2002, focuses on the validation of the template guides for the C1C2 transarticular screws (continuing on Goffin's work). As the others had previously done, he used a preoperative CT to base the design of the template and the testing was done on cadavers. During the first part of the study, he ran into some serious instability issues, particularly in the direction of the application of the drilling force; the result was dramatic displacement of the template guide on the order of 2 to 9 mm. Consequently, they altered the design of the template guides to include more support and had significantly better results. The appearance of the templates developed by the three researchers resembled a block that overlaps a broad surface area over the spine. From the experience gained in the RSB Lab, it was found that the templates that overlap a broad area on the spine create the necessity of removing a large amount of muscle tissue to expose the mating bony surface of the spine; this design is not beneficial for the patient and time consuming for the surgeon; in addition residual soft tissue causes the surface conformity between the template and the anatomy to be disrupted. As a consequence, the work on custom guiding templates for the spine was abandoned; however important knowledge was learned from the experience. Figure 15 shows the template designs published by other investigators and the initial designs of the spine template from the RSB lab. Customized position frames for intracranial surgery The first to suggest customized position frames for the head was D'Urso et al. [16]. His approach involved implanting reference sockets on a patient's skull prior to scanning. Then from the diagnostic images he create a physical 3dimensional model of the patient head that included the mounting sockets. The surgical trajectory was planned on the solid model by attaching a reference arc to the mounting sockets. By correspondence, the arc was refitted onto the patient and directed the surgeon to the predetermined location in the brain. The technology developed in the RSB for performing imageguided surgery by use of patientspecific frames has been shown to have superior accuracy and precision when compared to previously documented framebased procedures [33]. Rajon et al. [33] used a simple phantom to measure the point localization accuracy of the customized guiding frames. He found that the technology could position a probe tip with an accuracy of 1.7mm. As mentioned before, the standard in terms of accuracy for imageguided surgery is the framebased approach. The overall system accuracy of a standard frame based procedure was measured by Maciunas el at. [27] and determined to be 2.28mm. Their accuracy measurement is an estimated measure of the imaging accuracy and the mechanical system accuracy extrapolated for the ideal case of zero slice thickness (CT or MR slice thickness) to remove image localization related errors. Today in the RSB lab the researchers have developed customized positioning frames with specific component attachments to perform craniotomies and biopsies (Figure 26 and Figure 2 7). Both of these procedures are incorporated into a graphical user interface, with the aim of presenting the technology as another tool available to the neurosurgeon; one where he or she can preplan a surgical procedure and fabricate his or her own patientspecific frames and surgical guiding attachments. Rapid prototyping and RPD Designer. Rapid prototyping (RP) can be accomplished using a number of approaches. One such approach is 3dimensional printing RP (Figure 28A and B), another is subtracting rapid prototyping (SRP) (Figure 28D). Both RP technologies provide a means to fabricate a physical 3dimensional solid model from information contained in a stereolithography (STL) file. The file consists of planar triangular facets that represent a 3 dimensional surface [29]. In the case of 3dimensional RP, the first step in creating the solid model is the virtual slicing of the STL file by a preprocessing program [49]. For each virtual slice a layer is printed in real space. The bottom layer is built first, followed by additional layers until the last layer is built (Figure 28B). The RP system purchased in the RSB lab is built by ZCorp model 310 (Burlington, MA) and has a build accuracy of 0.3 mm [33]. The system creates the 3dimensional model by translating a thin layer of powder from a reservoir location onto the building chamber. Subsequently, a binder fluid fuses the powder as specified by the slicing of the STL file. Afterwards, the build chamber moves vertically downwards allowing for a new layer of powder to be placed onto the previous layer, and as the binder fluid is placed again it will also glue within layers. When all layers are built the solid model is taken out of the building chamber (Figure 28C) and infiltrated with cyanoacrylate glue to assure durability [33]. The 3dimensional solid, or in the case of this investigation the patientspecific frame, can be gas sterilized to be used within a sterile operative field [33]. The RSB lab purchase another machine to perform rapid prototyping, this machine is a programmable milling machine, and instead of performing 3dimensional printing, it starts with a solid object and removes unwanted material (Figure 28D). Because it removes material, this type of rapid prototyping is called subtractive rapid prototyping (SRP). There are several reasons why the researchers invested in the model MDX650 (Roland DG Corp., Knoxville TN); the main reasons are its speed and the capability of using a variety of materials to fabricate the reference guides. With the subtractive rapid prototyping machine the fabrication time was been cut in a third of the original time, from 11 hours to about 3 or 4 hours [34]. Also the reference frames created with this machine can be sterilized at high a temperature which is much faster than cold sterilization required by the patientspecific reference frames created with the 3 dimensional RP machine. In addition, early accuracy tests demonstrate that the precision of the SRP machine is superior to the 3D printing machine [34]. The creation of the STL file is conducted by a graphical user interface created in the RSB lab, called RPD designer [33]. The RPD designer has capabilities of rendering a solid model from a set of diagnostic images. The patientspecific frame is created by selecting a reference surface (called selected surface) by "painting" the surface of the virtual 3dimensional rendered anatomical model [33] (Figure 11). Significance There is a great potential for the new proposed idea of using diagnostic images and rapid prototyping to create a customized guiding frame for use in imageguided surgery. The most significant benefits of this new technology are rooted on the fact that it will eliminate the need for fixed frame application thereby reducing patient discomfort. It also has the ability to significantly reduce the cost of replacement parts related to framebase stereotactic surgery. In addition, the use of the patientspecific reference frames will eliminate registration processes required by frameless stereotactic surgery and allow for reorganization of an overwhelmingly crowded OR. The new technology will provide the surgeon with the ability to plan and create his or her own surgical guides while reducing patient wait times during surgical planning. In addition, the surgeon will be able to remove and replace the patientspecific reference frame onto the patient as required, and maintaining high accuracy on a range of cranial procedures such as shunt placement, skin incision placement, bone removal during a craniotomy, biopsies, tumor resections and deep brain stimulation [34]. Besides making stereotactic neurosurgery more minimally invasive, the patientspecific reference frame could maintain framebased accuracy [33]. Furthermore, the cost for creating a frame and a whole set of attachments to perform model based imageguided surgery is estimated to be approximately $50.00 per patient in materials cost [6]. Rajon et al. [34] has already created the software that provides the surgeon with the ability to directly build a virtual 3dimensional model from diagnostic image datasets. This program allows selecting a surface on the virtual 3dimensional model for the customized guiding frame and appropriately includes components to effectuate his or her approach towards a patient specific intracranial surgery. The significance of the present investigation is based of the fact that some level of testing is required prior to the fabrication of the customized guiding frame. The results of this testing will provide the clinician with prior information of frameonpatient stability. In addition, it will allow the clinician to effectuate his surgical plan with more confidence and accuracy. Facial Landmarks and Angles The stability of fit of the patientspecific frames (or masks) is greatly benefited by the surface characteristics of the human face. If a reference frame is to be designed based on the superficial surface of the human head, the contours of the human face serve the purpose of anchoring landmarks. The facial characteristics are ideal because of various neighboring irregular and curved surfaces. Besides, the facial contour is in close proximity to the brain, the skin is smooth and in most patients, close to the bony structure, to which the brain is naturally referenced. Therefore, it is important to obtain a basic understanding of the surface anatomy of the face. Figure 210 illustrates the facial anatomical landmarks of importance. On the forehead one can palpate the frontal eminences, which are prominences of the frontal bone that vary among individuals and are often unsymmetrical [10]. Caudal and medial to the frontal eminences is the glabella, a smooth somewhat triangular bony prominence on the frontal lobe located above the nasion, which is the point of intersection between the two nasal bones and the frontal bone. The nasion represents the root of the nose. Lateral to the glabella, bilaterally, one finds the superciliary arches: two arched elevations of the frontal bone that tend to be more prominent on the medial aspect of the face. Some males have notably prominent superciliary arches that may significantly contribute to the stability of the mask. Superciliary aches tend to be small in females and absent in children. Caudal and lateral to the orbit, bilaterally, is the zygomatic bone, which forms the prominence, known commonly as the cheek bone. The prominence of the zygomatic bone can vary significantly among individuals. Finally, the orbital area is also an important area of depression that contains the eye within the bony skull [10]. Furthermore, based on these facial surface anatomy landmarks, one can extrapolate facial angles that are important to our discussion, as they also play an important role in determining the stability of the mask (Fig 2 11). The Finite Element Method and ITK The patientspecific frames are complex in geometry and shape, and the best way to predict its behavior as it response to small external forces is by applying the principles of finite element analysis (FEA). The FEM involves dividing (discretizing) the entire frame into smaller parts or finite elements [5]. By dividing the entire frame into smaller elements, it is possible to obtain an accurate prediction of the displacement of the frame by approximating the displacement of each element to a polynomial function [28]. The patientspecific frames are modeled as rigid frames composed of linear tetrahedral elements under a deformable interface composed of truss elements (or linear springs) (Figure 5 2). The deformable layer, composed of truss elements, is attached at one end to the rigid layer, and on the other end all the truss elements are fixed (cannot moved in any direction x, y and z); this fixing of the elements on one end represent patient's surface. In addition, FEM has been used extensively in medically related mechanical analyses and modeling [12], such as in orthopedics [26], dentistry [48], head injury [25] and organ analysis[39]; and to evaluate stresses in human bones [37], to mention a few. For the present investigation, a linear 3dimensional truss element was added to the local ITK FEM (ITK 1.8, Kitware Inc., Clifton Park, NY) library because it was not available. The ITK FEM module was created for image registration and has three iterative solvers to solve large systems of linear equations; these solvers are slow and not adapted for the cases in the present investigation. A new linear solver, Cholmod [13] was adapted which can solve sparse systems of equations in seconds compared to hours that the available ITK's linear solvers take. Also, the ITK library was carefully studied and manipulated to make the methodology defined in this work progress as fast as possible (Appendix A). The ITK FEM module was chosen for this assignment to provide compatibility to the already built and working graphical user interface (RPD designer). Elements Used in the Simulation Two simple elements were used for this investigation. The first one is a 2node truss element in 3dimensional space called space truss. The second one is the linear 4node solid tetrahedral element in 2dimensional space, which was available in the ITK FEM module. The space truss Space trusses are frequently used to model lattice domes and aerospace structures [4]. They compose 3dimensional structures that are subjected to 3dimensional force systems. A space truss element is a long, thin and straight prismatic element. Prismatic means that its cross sections are all homogeneous by having a cross sectional area (A) and an elastic modulus (E) that remains constant through out the length of the truss. Space trusses are connected by nodes that act as frictionless ballandsocket joints [17]. In a structure composed of space trusses, any unsupported node (with no assigned boundary conditions) will result in translations in three local directions (x, y, and z). If no boundary conditions are applied to its nodes, then each truss element has six node displacements or six degrees of freedom (DOF) to establish its deformed position [4] (Figure 212). A space truss element is defined by the Young's modulus (E), cross sectional area (A), the length (L), and a Poisson's ratio (v). One important characteristic of the any truss element is that it has a preferred dimension, that is, along its longitudinal or axial direction (Figure 212). By definition, a truss element will only develop axial forces [17]. Any perpendicular forces will not change its original length (L), since the truss element does not undergo moments and its nodes act like frictionless balland socket joints. For example, Figure 213 shows a perpendicular force (fy) applied to a truss element; as a result of the force there is no difference between the original length (L) and the length succeeding the application of the force (L'). Hence, if there are not differences in lengths (L'L = 0), there were no internal forces and no displacements either (Figure 213). The space truss is considered a simple element since it deforms only along one direction (along its local xaxis) [5]. For simple elements like trusses, the direct stiffness method can be use to determine the elements stiffness matrix (ke). The direct method uses principles of equilibrium and mechanics of materials [46] to define the element stiffness matrix, which is a relation between node displacements (ui and u2) and internal forces (fi and f2) for each element e (Figure 212). The characteristic of the space truss element of only responding to axial forces makes the calculation of the element stiffness matrix relatively easy. When the element is analyzed in its local coordinate system, by definition there are no internal forces perpendicular to the truss and hence no perpendicular displacements. In other words, the displacement vector for a truss element is given by Equation 21, and the local node displacements ul and U2 are caused by the local node forces fi and f2, as seen in Figure 214. u= (21) u2 From mechanics of materials [46], the relationship between the local node displacements and the local node forces for a prismatic truss element is linear and given by Hooke's Law [46], or Equation 22, where k is the stiffness coefficient and is equal to AE/L. f = ku (22) With Equation 22 and the principles of equilibrium [28], the element stiffness matrix (ke) for a particular element e can be calculated. As Figure 214A shows, to establish the element stiffness matrix with the direct stiffness method [28], one needs to confine the element to a unit displacement at node 1 (ni) and hold node 2 (n2) fixed. The force at node 1 is equal to k* (by Equation 22). Since the element is in equilibrium, then the internal force at node 2 must been equal to k*1. A similar situation is true (Figure 214B) when node 2 is confined to a unit displacement, and the internal forces at node 2 and node 1 are equal to k and k, respectively. The element stiffness matrix is a 2 x 2 matrix defined in local coordinates by Equation 22. Sk k =AE[ 1 1 1(23) k = or ke = AE (23) k k L 1 1 For the stiffness coefficients kij the first subscript identifies the direction of the force, while the second subscript is related to the displacement. The element stiffness matrix represents the forces at the ends of the element as functions of the displacements at the nodes (Equation 24). f}=AEI ]{i1 { or }= ke L (24) f2 L 1 1 1U2 fJ2 2 The element stiffness matrix (Equation 24) defines the space truss in local coordinates. The transformation to global coordinates is necessary to represent the truss element as part of a entire structure. The transformation matrix and is represented by Equation 25. T=coso coso0 cos0 0 0 0 T = 0 0(25) 0 0 0 cos0, cos0y cos, Finally, for one element the local element stiffness matrix ke can be transformed into the global element's stiffness matrix Ke by Equation 26 (T means the transform of the matrix). K = TTkeT (26) For the space truss element, the element's stiffness matrix in global coordinates Ke results in a 6 x 6 matrix (Equation 27). cos26 cos cos ( cos cos cos2O cos c cos o cos cos O cos y o os cos co cos x cos 0 cos2 y cos y cos O Ke AE cos2 cos OcsO cos 0 cos O cos2 (27) L cos2 6 cos ( cos cos O cos Symm. cos2 2 cos 0 cos 6 Cos2 0 Following the direct stiffness method, the most important steps to define the space truss element systematically is first to define the element's stiffness matrix in local coordinates (Equation 23); then to calculate the transformation matrix T (Equation 25) necessary to swap from local to global coordinates; and finally it is important to define the element stiffness matrix in global coordinates Ke (Equation 27). In the FEM, all elements stiffness matrices in global coordinates (Ke for e = 1 to n, n is the number of elements) are assembled into a global stiffness matrix Kg that represents the force displacement relations for an entire structure composed, in this case, of space trusses. F, is a vector that that defines the external forces in global coordinates and Ug is the resulting node displacements in global coordinates (Equation 28). The solution to the system of linear equation is given by solving for the displacements in global coordinates, Ug. F, =KU, (28) The linear tetrahedron Solid elements are defined in 2dimensions and can model a solid body. The linear tetrahedron is the simplest solid element in 2dimensional space; four coordinate points define its geometry with respect to a global coordinate system. It is known to be poor for stress analysis, as results tend to be more rigid than expected [17]. The numbering of the nodes is important for generating every individual tetrahedron correctly. The numbering of the tetrahedron's nodes must be done following the righthand rule (Figure 214) and the calculation of the tetrahedron's volume must result in a positive value (Equation 29). First an initial corer it picked (node 0) and then the triangular face opposite to the initial corner must be numbered following a counterclockwise manner (nodes 1, 2 and 3). The calculation of the tetrahedron's volume gives evidence that the numbering of the nodes is done correctly. In Equation 29, the node coordinates for the tetrahedral are given by xi, yi, and zi where i is the node number. F 1 1 1 17 1 x0 x, x, x Volume = det 2 3 0 (29) 6 y, y, y, y3 6 Z0 Z1 Z2 Z3 The derivation of the tetrahedral element's stiffness matrix in complicated, in fact it rarely derived using the direct stiffness method. In turn, the principles of virtual work are often used to derive stiffness matrices for solid elements, multidimensional elements [4]. Characteristics of Stiffness Matrices The characteristics of stiffness serve to validate the setup of the ITK FEM module for the two types of elements used in this investigation. The characteristic of the element stiffness matrices for the space truss, the linear tetrahedral element and the global stiffness matrix (Kg) were analyzed for validation of the results. The element stiffness matrix can be singular; this means that the determinant is equal to zero. Also they need to be squared and symmetric and all element along the diagonal need to be greater or equal to zero [28]. In the case of the reduced form of the global stiffness matrix, it must be square, symmetric and positive definite and all its diagonal values have to be greater than zero [28]. In addition is must have and inverse so the displacement can be solved uniquely for a given force [46]. Positive definite matrices arise frequently in FEM problems. They are characterized by the condition that all their eigenvalues are positive [13]. One of the properties of a positive definite matrix is that it can be decomposed by Cholesky decomposition. If A is a symmetric positive definite matrix, then A can be factored into a product of L and LT (T means the transpose of a matrix), where L is a lower triangular matrix with positive diagonal elements. The definition of a positive definite matrix is given by Equation 210, where X is any vector in R space [13]. XTAX > 0 (210) Linear Solvers The ITK FEM library is fitted to solve image registration problems that incorporate finite elements. The ITK FEM module offers iterative solvers mainly used to solve very large matrices, common to image registration problems. The most reliable iterative solver in the module is the ItPack Solver [24], however due to its iterative nature ItPack is a very slow solver. An alternative to using ItPack is a solver based on Cholesky decomposition called Cholmod [13]. Cholmod is a set of routines for factoring sparse symmetric positive definite matrices and solving linear equations used in Matlab (The Mathworks Inc. Natick, MA). Cholmod proved to be much faster and highly reliable; as a result it was adapted to be used with the ITK FEM module. W4gi N1OW l .l rg IN I&4ll 4 kr I I f kll.. bql. Inbl7. b II Ll . gure 21. Pretarget localization and presurgical planning is done in a computer workstation. "v Off& _0'.% . Figure 22. Procedures related to framebased image guidance surgery. A) A patient with a head ring prior to scanning on a CT. B) A surgeon setting the coordinates on a stereotactic apparatus, which will be mounted onto a head ring (reference frame) to located the coordinates for targeting a specific location on the brain. Figure 23. Frameless stereotactic surgery uses fiducials that are glued to the patient's skin. The fiducials serve for registering the patient's anatomy in real space to the image datasets on a computer workstation. The tools for the surgery have similar fiducials so they can be tracked by cameras. F A Figure 24. Typical layout when frameless stereotactic surgery is employed. Both images show how crowded the OR can get by the necessity of additional computer workstations and tracking devices. Figure 25. Template designs for the spine. A) Goffin [20], B)Birnbaum [6], C) and Yoo [50]. D) The RSB template designs from the first stages on the development of customized positioning frames. Figure 26. Router guides for craniotomy. A) Selection of the customized positioning frame onto the patient's virtual 3D model. B) Selection of craniotomy contour drawn on the surface of the skull. C) Virtual model of the customized positioning frame. D) Craniotomy attachment for the customized positioning frame. E) Fabricated structure that will guide the surgeon to perform the craniotomy as planned [34] Figure 27. Biopsy guides. A and B) Graphical user interface shows the three orthogonal views from diagnostic images and the solid virtual model. The surgeon can select the target location and the graphical user interface will appropriately position the biopsy guide C) Virtual model of the customized positioning frame and the biopsy guide attachment. D) Fabricated structure that will guide the surgeon to perform the biopsy as planned [34]. Figure 28. Rapid prototyping machines. A) Rapid prototyping machine by ZCorp model 310. B) A thin layer of powder is moved to the building chamber by the gantry, which also holds the print head that distributes the binder fluid. C) The 3D model as removed from the building chamber. D) Subtractive RP machine. 4h~ 7 TT1 GUI: Neurosurgeons virtually define a target location I Step 1: GUI renders model and selects surface for frame Step 2: A) Virtual evaluation of frame prior to fabrication. B) Is the patientspecific frame design STABLE? I ..i m.a.. . N  i I, Step3: If YES, then build patientspecific reference  alI I 1 i I 1  " i . Figure 29. Chart to visualize the significance of the present work. a* ai Frontal enilen&. . Frontal Superciliary arch Glabella I / Nasion Radix Cephalic C ( audal MedialLateral Cheek (z omar;ti mcI') ,, Figure 210. Facial landmarks. Nasiofrontal angle I Nasiofacial angle / Nasiolateral angle .,' Figure 211. Facial angles [10]. The nasiolateral angle is calculated at the level of the radix. Figure 211. Facial angles [10]. The nasiolateral angle is calculated at the level of the radix. Element e 112 X 1 V2w U2 LU2 w1 X Figure 212. Space truss element has six displacements to establish its deformed position. There are ul and u2 along its xaxis, vi and v2 along its yaxis and wl and w2 along its z axis. Though the space truss is defined in 3dimensional space, it will only develop axial forces, along its local xaxis. symbol Figure 213. Application of a force perpendicular to the axis of the space truss. There is no change in length (L'L = 0) and there are no internal forces (in a real simulation the truss element will not display as in the figure). The boundary condition symbol implies, in this case, that the element is fixed in local x and z directions. f, =k*1 fi  z A z n2 x u2= I. f2 = Sn2 f2 Figure 214. Space truss in local coordinates. A truss element in its original position is denoted by the thick black line, and the same truss element in a deformed state is shown in gray. After a force fi is applied to node 1 (ni) to create a unit displacement the stiffness coefficient kn1 defined the element at n2. The reaction force at n2 is k12. Similarly, a force corresponding to a unit displacement is applied at n2 and k21 and k22 are defined. Nodes x Figure 215. The linear tetrahedron has a unique numbering scheme. k*1 L ,,I  CHAPTER 3 MATERIALS Three patients with notably dissimilar face profiles were selected as training data for the methodology designed to test uniqueness of fit. With CT datasets and the graphical user interface, RPD designer, three anatomical solid models (Figure 31) were built and name Caucasian, Asian and infant because they belonged to a Caucasian male, an Asian male and a infant respectively. Also six masks were built, two for each model (Figure 32). For each solid model, a conforming small and a large mark were built; the small masks included considerably less facial landmarks than the large masks. While the large masks included the entire orbital area and extended laterally and medially along the superciliary arches, the small mask only included the upper orbital region and extended to the level of the glabella along the nose region. All masks were created from the "painting" of the corresponding selected surface as shown in Figure 51. The six masks designed with RPD designer are the training data for the methodology to test stability and uniqueness of fit. In reality when the patient specific reference frames are applied onto a patient, the clinician will apply a force to keep the mask in their unique location. Sliding, rotation and translation of the mask on the level of contact with the patient's face signifies exaggerated displacements at the level of the attachment tips, as shown in Figure 11D. While the large mask designed for training the methodology maintain their location when exposed to real transverse forces, the small mask would not withstand normal forces at the their extremities, along the length of the superciliary arches. To provide with a real world measure of stability for the small and large masks, a center pin was place normal to the surface (or the direction on which the mask is placed onto the solid model) and at the location of the glabella on all masks. Oriented by the center pin, crosshairs were drawn on all masks along the transverse and sagittal planes as shown in Figure 3.5. Additional pins were located close to ends of the crosshairs lines as shown in Figure 35 and numbered from 1 through 4. Pins 1 and 3 were placed on the sagittal plane and pins 2 and 4 were placed on the transverse plane (Figure 36). A normal force of 5 N, applied by weight, was applied onto all pins for large and small masks (Figure 37). The large masks remain in their unique location, and sustained the 5N force; on the other hand, the small mask exhibit large displacement when the force was applied on pins 2 and 4. This resulting slippage behavior is mainly attributed to the smaller surface are and less facial irregularities covered by the small masks. All small masks withstood the 5N normal force on pins 1 and 3. In addition, the placement of all large masks is unique, this is not the case for the small masks; there is not clear unique location where the masks fit. The sliding of all small masks when experiencing a normal force on pins 2 and 4 and the lack of a clear unique fit are evidence that the small masks are unstable, on the contrary all large masks withstood the 5N force on all pins and have a unique fit, which characterize them as stable (the fitting of the masks onto the models can be seen in Chapter 6, Figures 61 through 63). In addition to the masks, two other surfaces were created to evaluate the methodology of fit. A caplike reference frame (Figure 33) was created from CT image dataset and a half sphere surface (Figure 34) was virtually created with Paraview (Paraview 1.6, Kitware Inc., Clifton Park, NY). The cap frame covers the top of the head and extends orthogonally along the coronal and sagittal planes no further than at the level of the frontal eminences. The half sphere surface is a faceted surface defined by 31 points evenly distributed (Figure 34). The half sphere needed to be defined as a faceted surface with ten angular rotation about it axis and three coplanar levels, because a smoother surface could can not be evaluated, since it results in a mechanically underconstrained system. Figure 31. The three models used to evaluate the methodology for testing unique fit. A) Caucasian, B) Asian and C) infant. Figure 32. Masks built for testing the methodology of fit. A) Caucasian small and large masks. B) Asian small and large masks. C) Infant small and large masks. 50 Figure 33. Caplike reference frame. A) Upper head model and B) caplike frame. A 1 B Figure 34. Virtually created half sphere test model. A) Top view and B) side view. Figure 35. Placement of pins for a simple test on the real stability of the large and small masks. A) Large Asian mask and B) small Asian mask. Figure 36. Coronal view of placement of the pins onto the real masks. Notice the orientation of pins 2 and 4 and difference in lateral extend between the A) large Asian mask and B) the small Asian mask. Figure 37. The method employed to demonstrate the difference in stability between the large and small masks. The weight represents a 5N force, and it is placed onto a pin to observe the resulting displacement of the small Asian mask. CHAPTER 4 METHODOLOGY FOR MEASURING STABILITY Introduction The methodology for testing stability of fit of the patientspecific frames relies on the use of the finite element method (FEM) to model the contact interactions between the patient specific frame and the patient's facial contour. The objective of the methodology is to produce a score value that correlates with the physical stability of the six masks created. The methodology quantifies the stability of fit by information from the simulated movement of the patientspecific frame as it responses to small external forces. The virtual mask is able to experience displacements because it is modeled as a rigid solid surface, which is virtually placed on top of a deformable bed of springs (Figure 41), and loaded by small external forces. The small external forces are calculated to cause a spring deflection of 1%, because the simulation does not include friction. For every small external force applied onto the rigid layer, internal forces transfer onto the bed of springs, and by the deformation of the springs, the rigid layer can either move closer to the model in some areas and move away in other areas and still sustain contact; or it can predominantly slide or rotate (Figure 42). Sliding and rotation occur when there is relative motion between the two surfaces. The simulated movement of the mask is analogous to testing the fit of a cast by placing it over its mold and applying pressure at various locations to see if it fits securely. In a similar way, the approach used in this investigation involves applying small forces in many directions and locations and monitoring the resulting motion of the mask (cast). As described in the materials section, three solid lifelike facial models that represent three different patient's anatomies were constructed with 3dimensional printing rapid prototyping. For each solid model two masks were designed with the existing user interface RPD designer, one of the masks is known to be stable and the other is known to be unstable (from information on Chapter 3). The stability information on the real masks was used as training data for the methodology designed in this investigation. The methodology designed to test stability of the patientspecific reference frames is presented in four sections; a) system components and FEM setup, b) program parameters, and c) output parameters, and d) the criteria defined to measure stability of fit of the patientspecific reference frames. As mentioned before, the task of measuring stability is critical to understanding the appropriateness of the reference surface selected and thereby appreciating the predicted accuracy and precision of the clinically designated patientspecific design. While this is only a piece of a much broader investigation that seeks to evaluate an alternative method to imageguided surgery, it is critical to this application of patientspecific tool design. System Components and FEM Setup The system components are the external forces, the rigid layer (virtual mask or frame), the bed of springs, and the boundary conditions; all are defined in the context of FEM. The STL file, output of the graphical user interface (RPDdesigner), is the only input into the system. The STL file is standard format representation of a 3dimensional surface as an assembly of planar triangles (Figure 43) [29], the point coordinate information of the vertices of each triangle are used to define the rigid layer, the bed of springs and the boundary conditions. The first component of the FEM setup is a rigid layer and it is composed of tetrahedral finite elements. This layer is virtually constructed between the selected surface and a surface extrapolated along the positive normal direction at each vertex (Figure 41 and 43). The second component of the system is a bed of springs composed of individual linear springs that are placed between the selected surface and another surface extrapolated along the normal negative direction of each vertex (Figure 41). Hence, each individual spring in the bed of springs is constructed between the point coordinates of the selected surface and the boundary condition points (or nodes). On the whole, the rigid layer represents the patientspecific reference frame (or mask), and the boundary conditions represent the patient's facial surface. With the stimulation of external forces, the bed of springs will allow for the virtual displacement of the rigid layer, simulating the stability or instability of the real mask onto the patient's facial anatomy. External Forces The external force direction is defined with respect to a main normal (Nmain). Nmain is a normal perpendicular to the coronal plane and along the sagittal plane of the skull; in other words, in the direction in which the frame will be placed onto the model (Figure 44). The angle theta (0) is the tilt angle from Nmain and p is the spin angle about Nmain. The forty nine forces are calculated by incrementing theta (0) six times, from 0 to 900 in intervals of 15; and by rotating the angle p (spin) by eight forces vectors in increments of 45 about the Nmain (six tilt angles multiplied by eight spin angles plus one normal forces is equal to forty nine forces vectors). All forty nine forces are applied onto every node of the mask with the exception of the edge nodes (Figure 45). The force magnitude is very small enough to displace a single spring by 1% and will be discussed in the program parameters section. For example, in a case of a mask similar to the one shown in Figure 45, which is composed of about 300 nodes, 49 force vectors will be placed on each node, that will result in 14700 forces applied onto the entire mask. For each of the 14700 force vectors, all 300 springs will have a resulting displacement. Rigid Layer Probably the most important characteristic of the system is the difference in strength (elastic modulus) between the linear springs, which compose the bed of springs, and the rigid layer. The rigid layer in comparison to the bed of springs must be rigid in the sense that it can translate or rotate but not deform. The rigid layer is composed entirely of four node tetrahedral elements. From the extrapolation of the point coordinates that define the planar triangles in the STL file, virtual wedges are built. In each wedge three tetrahedral elements are fitted (Figure 4 6), and created following the righthandrule as explained in Figure 215. In order to define a four node linear tetrahedron in ITKFEM module [24], four point coordinates are need, and two mechanical properties, the modulus of elasticity (E) and a Poisson ratio (u). RDR ratio. It is important to measure if the rigid characteristic of the rigid layer are maintained at every force application; for that purpose a program flag was created and called RDR ratio. The RDR ratio is an indicative of major deformations of the rigid layer. The RDR is calculated with respect to the upper surface of the rigid layer and for each external force application because every external force application has the potential of deforming the rigid layer. To calculate the RDR the maximum difference in distance between a vertex (node) i and its neighboring vertices before (La) and after (Lb) the application of a force is measured and called maximum rigidity fault (mrfi) (Figure 47). As Equation 41 shows, the RDR is the ratio of the maximum rigidity fault (mrf) at a vertex i divided by the sum of the displacement magnitude (di) and the mrfi. mrf, RDR, = m (41) mrf + d, If RDR is close to zero (or one over a small displacement), it means that the rigid layer has not deformed. If the RDR ratio is close to one, the rigid layer has taken the entire load and it has deformed at a particular node. In the situation where the rigid layer deforms, the RDR ratio will indicate so and the simulation will not be able to provide reliable results. In the simulation, if the RDR ratio exceeds 0.05 (or 5% rigidity loss) then the force was not included in the results. For all the masks and frames evaluated, the RDR ratio had a mean of 0.3% and a standard deviation of 0.2%. Bed of Springs The bed of springs is the only deformable component of the system given that the individual springs that comprise it considerably are less strong than the rigid layer. The purpose creating a bed of springs is to identify and monitor the individual springs that acquire a state of tension or compression. A linear spring has a characteristic spring stiffness referred to as the spring constant k [46]; it is related to the strength of the material of what the spring is made of. In addition, the spring has a natural length (Lo) that is the length of the spring at rest (Figure4 8a). Figure 48b shows an axial force Fy placing a spring under a state of compression and changing its length to a new length L1. Figure 48c shows the freebody diagram of the system and the spring internal force Fs, which is determined by the equation in Figure 48d. In the Finite Element Method, linear springs that reside in 3dimensional space are modeled as trusses elements in 3dimensional space. In order to define a linear spring in 3dimesions in ITKFEM module, 2 point coordinates are needed, a modulus of elasticity (E), a Poisson ratio (u) and an area (A). The area for each spring element i is the sum of one third (because each triangle has three vertices) the areas of all triangles surrounding the vertex i and denoted in Equation 41 as Ai. In this manner, even though the triangular mesh in the STL file is irregular (Figure 43); all the springs in the entire bed have the same proportional strength, because the spring stiffness coefficient for each spring element (ki) is set proportional to the area Ai, as seen in Equation 41 [46] a the elastic modulus (E) and the spring length (Lo) are kept constant. A E k = E (42) L Program Parameters Input Parameters In addition to the STL file, the program requires five input parameters that are set for convenience of coding and not to achieve a particular real world test, they are: the natural length (Lo) for all linear springs that compose the bed of springs, external force magnitude, the height of the rigid layer (h), the modulus of elasticity for the spring elements (Es) and the modulus of elasticity for the rigid surface (Er) [46]. Spring height The spring height is set to 0.1 mm, primordially because the displacements are expected to be very small, close to 0.001 or 0.002 mm. The system was design to include linear spring element which do not carry friction and when grounded in all three direction they behave as a ballandsocket joint. External force magnitude The external force magnitude is calculated by Equation 43 and selecting an expected axial displacement (AL) of 1% the spring's length. As shown by Equations 32 and 42, the force required to move a spring a certain displacement is proportional to the cross sectional area of the spring. In a similar way, in order to calculate a average external force magnitude that corresponds to a single spring's displacement of 1%, the force must be weighted by the mask's surface area per number of springs, as shown in Equation 43. Area *E AverageForce spring AL(43) spring Lo Rigid layer height The height of the rigid layer was set to 1.0 mm because increasing the high of the rigid layer increases the rigidity of the frame, and 1.0 mm was found to be the appropriate value (Table 41). The rigidity of the frame was quantified with the RDR ratio and Table 41 shows that reducing the height of the rigid layer considerably affects the rigidity of the rigid frame. For example, in Table 41, a RDR of three over one thousands signifies a 0.3% rigidity loss. In addition, the rigidity is affected by the shape (length and height) of the tetrahedral elements [41]. Since STL information is used to mesh the rigid tetrahedral layer, and the STL facets are not regular triangles, the shape of the tetrahedral elements affects the rigidity; by increasing the height of the rigid layer the tetrahedral elements behaves more rigidly [41]. Difference in elasticity The elasticity of the rigid layer (Er) is set to have a modulus of elasticity six orders of magnitude higher than that for the bed of springs (Es). Six orders of magnitude is the largest difference that can be established without compromising the numerical stability of the global stiffness matrix. The sparse matrix linear solver Cholmod [13] has an output value which signifies the conditional number of the matrix that relates to the numerical stability of the system of equations, called rcond. It was found that unreliable values resulted when rcond is equal or exceed le13. A difference of six orders of magnitude for the Er and Es maintains the rigidity and results in a reliable simulation. Table 42 shows the how the difference in elasticity between Er and Es affect the RDR ratio and Cholmod's rcond value. Output Parameters For every external force applied onto the rigid layer, there will be a resulting displacement vector (di) at every spring i that make up the bed of springs (Figure 49). The displacement vector for each individual spring is characterized by two components, the normal displacement along the axis of the spring and the tangential displacement perpendicular to the axis of the spring. The normal displacement scalar (or distance) (ni) is calculated by the dot product between the direction of the spring (Ni) and the displacement vector (di), (Figure 49). The tangential displacement scalar (ti) is calculated by the cross product between the direction of the spring (Ni) and the displacement vector (di) divided over the norm of Ni. As Equation 44 shows, the compression or tension state of a linear spring is determined by its normal distance (n). If the normal distance is less than or equal to zero, the spring is known to be in compression, and if the normal distance is greater than zero, the springs is know to be in tension. ni compression = di Ni <= 0 (44) ni tension = di Ni > 0 i is the number of springs In a similar manner, regardless of the direction of the tangential displacement vector, the tangential distance ti is the perpendicular distance from the axis of the spring (Ni) (Equation 45 and Figure 46) to the final position of the spring (point c on Figure 49). ti= Ildi x Ni/Nil i is the number of springs (45) Linear springs have an inherent characteristic of only responding to axial forces, their internal forces (which are triggered by an application of an external force) should only be measured by it normal displacement (n). This means that springs (in the simulation) that predominantly have tangential displacements are not developing internal forces and are not contributing to the stability of the frame onto the patient's model. For example, in Figure 410, the external force (F) applied onto a rigid box results in the development of internal forces in the springs aligned with the external force (F), or springs 1, 2 and 3. The tangential displacement and change in length of springs 4, 5 and 6 is due to their attachment to the rigid box, and though the spring is stretching, its normal displacements is very small (insignificant) and therefore they do are not developing noteworthy internal forces. However, the tangential displacement of an individual spring is a good indicator of the sliding, slipping or rotation of the frame on the patient's model because it gives evidence on the internal forces being developed by other springs. When there are few springs aligned with the external forces the mask experiences tangential displacements and these are expected to be high. In accordance with the inherent characteristic of linear springs, the final length of a spring (L' in Figure 49) is not a good indicator of the internal forces developed by a linear spring as a response to external forces; however, the final length does helps in determining which springs are exceeding the bounding compression region as shown in Figure 411. Criteria For this investigation two criteria were defined. The first criterion is a single spring criteria, which determines if the resulting displacement of a single spring is acceptable to aid in the definition of stability. The second criterion is a single force criterion that looks to define if the frame sustains contact and has a unique fit at the onset of a single external force. With the information from the single force criterion an overall instability score for each mask is determined. This is done by providing a percentage ratio between the forces that result in instability and the total number of forces that were applied onto the mask. The criteria was designed with the expectation that for the stable mask cases, the displacements of the springs will fall on the upper area of the sphere in Figure 411, because the external forces applied in the simulation are small. Accordingly, it is expected that the stable masks would not experience full compression or come close to the boundary condition surface. For this reason a ratio called the TNT ratio, which is an indicative of the tangential displacement to normal displacement was define to aid in differentiating between springs that are in a state of compression, but experience mainly tangential displacement. The TNT ratio is a part of the single point criterion, its definition, the determination of a threshold value and its consistency is discussed in the following sections. TNT Ratio The TNT ratio is a normalized measure of the tangential movement of an individual spring with respect to its normal displacement. The TNT ratio is an indication of the fractional slippage of an individual spring. It is calculated from each individual spring's tangential (ti) and normal (ni) displacements (Figure 46). Numerically, it is defined by Equation 46 as the tangential displacement magnitude (t) divided by the sum of the normal displacement magnitude (n) and the tangential displacement magnitude (t) (Figure 49). TNT = (46) t+n The purpose of the TNT ratio is to differentiate, for small displacements, between springs that are predominantly sliding to springs that are predominately in compression. In Figure 410 a force (F) is moving a rigid box by a unit displacement in the negative x direction and is placing three horizontal springs (springs 1, 2, and 3) in compression; since the three horizontal springs are experiencing only compression (n = 1 and t = 0), the TNT ratio for all of them is equal to zero. Assuming the box is perfectly rigid, the vertical springs (springs 4, 5 and 6) will be mostly sliding, (n = 0, t = 1) and their TNT ratio will be equal to 1 by Equation 46. After the application of the force (F), the horizontal springs will be the only ones developing internal forces to sustain the movement to the rigid box. The vertical springs in Figure 410, are not contributing to with their internal forces to the equilibrium of the rigid box, since linear springs are characterized by not responding to perpendicular forces. The displacement of these springs (springs 4, 5, and 6) result from the sliding movement of the rigid box. The TNT ratio distinguishes springs that slip (or move tangentially) like springs 4, 5 and 6 in Figure 410, from springs that compress, like springs 1, 2 and 3 in Figure 410. Furthermore, Figure 413 shows two springs in compression; Figure 413A shows a spring with a TNT ratio of 0.5 and Figure 413B shows a TNT ratio of 0.833. In the case of the latter spring, its normal displacement is in the compression direction (n<=0) and final spring length is smaller than the original length (L' state of a mask. Its behavior is mostly a result of the sliding movement of a rigid body. If a surface has irregularities that can stop the rigid frame from sliding, then the TNT ratio for a particular spring will result in small values (close to 0). In contrast, if the surfaces are smooth, the bodies can move excessively relative to the grounded (fixed) springs, resulting in springs with high TNT ratios that will approximate one (Figure 414). Single Point Criterion The single point criterion determines if an individual spring is contributing to the stability of the rigid frame onto the patient's model or not. By the single point criterion a spring that contributes to the stability of the frame by experiencing a considerable state of compression is called a "good compressor". Table 43 shows that a "good compressor" is a spring defined by a normal displacement less than or equal to zero (n<=0), a final spring's length (L') less than or equal to the original spring's length (L) (L'<=L) and a TNT ratio below or equal to 0.7. In addition, Figure 412 shows a 2dimensional interpretation of the single point criterion. In order for a spring to be considered a "good compressor", it must remain inside the red area (Figure 412). The red area is bounded by the area below the horizontal red line (including the line), which indicates that the normal displacement is less than or equal to zero (compression); the area inside the half circle (including the line), which indicates that the spring should remain inside its compression bounds and the area under the incline line, which represents a TNT ratio set to 0.70. The slope of the incline line (blue) corresponds to a normal displacement of 1 and a tangential displacement of 2.41, which corresponds to a absolute TNT ratio of 0.7. As mentioned before, it is expected that the springs have small displacements, so the TNT ratio of 0.7 (inclined line) will differentiate between springs that either have a very small normal displacement and a large tangential displacement or springs that have a very small normal and tangential displacement; in both case these springs are not considered "good compressors" because they are not contributing to the stability of the rigid frame by not developing considerable internal forces. The TNT threshold was set to 0.70 because it was defined empirically to correlate with the information about the real stability of the large masks and the instability of the small masks. For example, Figure 413B shows that there are occasions where the normal displacement shows a compression state and yet the tangential displacement is high; by defining TNT upper threshold of 0.7, a spring like in Figure 413B will be excluded from being considered a "good compressor", it has a TNT ratio of 0.833. TNT threshold The threshold for the TNT ratio was defined empirically to be 0.7 based on the training data provided by the six masks constructed with rapid prototyping and their fitting onto the solid models (Figure 71, 72, and 73). To determine this empirically the simulation was ran with TNT ratios of 0.6, 0.7 and 0.8 for large and small masks (Figure 415 and 416 respectively). In the simulation a TNT ratio of 0.70 corresponds with the physical stability of the large masks and instability of the small masks. Figure 415 shows color coded surface plot for all three large masks, Caucasian, Asian and infant at different TNT ratios (0.6, 0.7 and 0.8). The location of the color dots correspond to the location where a virtual force was applied and three or more points behaved as "good compressors" (Table 43). The color of the dots corresponds to the maximum tilt angle (for all spin angles) where at least three "good compressors" were found. For example, a red dot signifies that for a maximum tilt angle of 90 (and for all spin angles) three or more "good compressor" springs were found, that is 49 forces resulted in at least three "good compressor" springs. The first row in Figure 415 associates to a TNT of 0.6, and as the figure shows, the Asian and the infant mask are noted as false negatives, because the color coded surface plots show that the maximum tilt angle in which the forces was withheld is smaller than 900, this indicates that the application of a force resulted in an unstable configurations (this information contradicts that premise that the large masks are stable). The second and third rows, which correspond to TNT ratio of 0.7 and 0.8 respectively, concur with the physically stability of the masks; but in Figure 416, a TNT threshold of 0.8 dictates a case of false positive for the small Asian mask, which is mark as maintaining alignment for all forces at all locations. This information also contradicts the premise that the design of the three small masks results in a unstable configurations. Figure 414 shows that when a tangential displacement magnitude is 2.41 times bigger than the normal displacement magnitude the TNT ratio is equal to 0.70. Also a TNT ratio of 0.5 occurs when the tangential displacement is equal to the normal displacement, and the TNT ratio approaches one when the springs experience mostly tangential displacements (Figure 414). In summary, a TNT upper threshold was empirically defined at 0.7, this means that all springs that are in compression and exhibit TNT ratios of 0.7 or below are in a "good" compressive state. On the other hand, the springs that exceed a TNT ratio of 0.7 will be regarded as mainly slipping or sliding. Consistency of TNT threshold To test the consistency of the TNT threshold, the force magnitude was increased five and ten times, and the number of point on the mask was increase for a single mask case (Figure 4 17). When the force was increased by five and ten times there were no changes in the color coded surface plots shown in Figure 414. This means that the TNT threshold is not sensitive to force. In the case of increase number of points, the results are shown in Figure 418; which illustrates that with the exception of one single point in the upper left corer of the Asian large mask, which held stability up to a theta value of 75, the TNT threshold of 0.7 corresponded with physical results. Single Force Criterion The purpose of the single force criterion is to define on a per force basis if the frame or mask is stable or not. The criterion requires a minimum of 3 points in compression, that is 3 "good compression" points (as defined by the single point criteria) to consider a mask stable with relation to a particular force vector. Three springs in compression will overcome for translations and rotations of the frame. Translation can be recognized by evaluating the movement of a single point, but two points are needed to recognize rotation, since rotation occurs about a normal vector and at a specific point that may not move. Translations or sliding is defined as the relative motion between two surfaces and rotation is the angular motion about a common normal [45]. The criterion requires a third point because three points are necessary to identify rotation in 3dimensions. If the displacement of the rigid body results in three "good compression" points, the simulation indicates that the rigid surface is moving close to the grounded points rather than sliding or rotating, which in turn indicates that the mask and model surfaces remain in alignment. If by the onset of an external force the mask results in a configuration where there are less the three "good compression" points, then the mask is known to be unstable for that particular force. Intuitively, more than three "good compression" points and their dispersion can give evidence of a level of stability. More "good compression" points and broadly disperser leads to greater stability. Therefore, a bounding volume defined by the minimum volume of all "good compressor" points compared to the bounding volume of the entire mask, is a good indicator of the level of stability of the mask for a particular force. Box ratio: The box ratio quantifies the level of stability as defined by the single force criterion by calculating a ratio between the volume of the compression patch (cvol) and the volume of the entire frame model (fvol) (Figure 418). The box ratio is calculated on a per force basis. The compression patch volume is the minimum volume that includes all "good compression" springs (Figure 419). If the mask is unstable, with respect to a particular force (as defined by the single force criteria), then the compression patch volume is equal to 0. The small box in Figure 418 shows a compression patch volume (cvol), in comparison with the volume of the entire frame (fvol). The box ratio is numerically defined by the compression patch volume cvol over the entire frame's volume fvol (Equation 47). Box ratio = cvo (47) fvol As Equation 47 shows, a box ratio is equal to 0 when the compression patch volume is equal to zero (or the mask is unstable) (Figure 419). When the mask is found stable for a particular force, a compression patch volume takes a value greater than zero. A box ratio close to 1 indicates that the compression patch volume approximates to the volume of the entire frame (Figure 419). Box ratios closer to the value of zero correspond to smaller compression patch volumes or "good compressor" points close to each other. A score for the overall instability (instability score) of the masks can be derived from the single force criteria by providing with an percentage value for the number of forces that result in instability (box ratio equal to 0) to the total number of forces applied to the mask (Equation 48). If this value is equal to 0%, there were not external forces that made the mask loose stability, and if this value is 100% it means that all external forces disrupted the alignment of the mask onto the patient's surface. #Ins y s e of forces the result in instability 100 Instability score = x 100 (48) total number offorces applied onto the mask In summary, a total of 49 external forces will be applied to all nonedge points onto a mask (if a mask has 300 nonedge points, there will be 14700 forces applied to it). For each individual external force applied onto a mask, there will be displacements of all springs that comprise the bed of springs. The single point criteria will individually classify each spring as a "good compressor". The number of "good compression" points is recorded. If there are less than 3 "good compression" points, the mask is considered unstable, the compression patch volume is equal to zero, and the box ratio is equal to one. If there are 3 or more "good compression" springs then the mask is considered stable and the compression patch is calculated and the box ratio will have a value between 0 and 1. As an overall percentage instability score for a mask, the number of forces the result in instability is compared to the total number of forces applied onto the mask. A value of 0% means zero instability (Equation 48). Table 41. Relationship between the tetrahedral height and a measure of the rigidity of the rigid layer. Tetrahedral height 1 mm 0.1 mm 0.01 mm Mean RDR 3/1000 (0.3%) 2/100 (2%) 2/10 (20%) Table 42. Table that shows the relationship between the difference in elastic modulus between the rigid layer and the bed of springs and the rigidity measure and Cholmod's rcond value. ErEs Mean RDR Cholmod's rcond 1.0e6 3/1000 2.14e8 1.0e7 2/1000 2.14e9 1.0e8 5/1000 2.14e10 Table 43. Single point criteria. Compression Tension No Slip n <= 0 and n > 0 or L' <=L and L' >L TNT <= 0.70 "good compressor" Slip n <=0 and n > 0 or L' TNT > 0.70   Outer extrapolated Surface Selected Rigid layer Surface Bed of springs / Boundary Conditions B Figure 41. The methodology for stability of fit involves simulating a rigid layer onto of a deformable bed of springs. A) The selected surface (red) as painted on the rendered model (green). B) A cross section through the selected surfaces and anatomical model showing the idealization of the contact mechanics between the rigid frame and the solid model Figure 42. Displacement of a mask as it responds to external forces. Figure 43. Triangulation of a customized positioning device proper to STL format. Every vertex in the mask is a point coordinate in 3D space. 1 Sagittal plane I I nI I 7 I r^"^ I CIL I I I .. I I r I I . : I . I IJ I I I II I III I If I ~C oronal iplann I j _ I I I I Tilt (0 = 0 Nmain (p = 0 to 3600) Figure 44. Direction of the Nmain. Nmain is specifically chosen for each mask, and forty nine force directions are calculated by defining a tilt and a spin angle. Edge points Figure 45. Nodes point on the rigid frame. Forces are not applied on edge points. Tetrahedral 3 3 0 Tetrahedral 2 Figure 46. In a wedge three tetrahedral elements are fitted. The nodes are shown next to the black dots. The numbering of the tetrahedral elements is important and should be done as shown in the legend. Tetrahedral Nodes 1 0123 2 1234 3 2345 Maximum rigidity fault Figure 47. Schematic showing how the maximum rigidity fault is calculated. y LI L1 x // /// Fs = k(L, L,) D Fs Figure 48. The behavior of a linear spring element under the application of an axial force. A) A spring at rest, nl is the natural length of the spring. B) The same spring under compression by a force Fy. C) Free body diagram of the spring in compression. D) Equation of a linear spring. / . JId II / Figure 49. Schematic that demonstrates how each spring's displacement was evaluated by identifying the normal (ni) and tangential (ti) displacement magnitudes. The normal (ni) and tangential (ti) displacement are scalar distances. In this particular case the spring is in compression because the normal displacement magnitude ni will result in a value less than or equal to zero, also notice that the natural length of the spring has changed from L to L'. 1 2 3 5 6 F Figure 410. Schematic that shows why it is important to distinguish between the springs that are mostly in compression or tension and those that are basically slipping or have a high TNT ratio. The force (F) exerts a unit displacement on the rigid box in the x direction. As a result of the force the springs have normal (n) and tangential (t) displacements are shown in table. n<=O In n L vI 1...... ..... .......... ...L L ..... ... Figure 411. Bounding volume for a spring to be in compression. For the stable cases the springs displacements are expected to reside on the upper area of the half sphere (gray area). Springs n t TNT 1,2,3 1 0 1 4,5,6 0 1 0 ~c~5~5~5~5~5~~ TNT<=0.70 or n/t = 0.41 Figure 412. Twodimensional interpretation of the single point criteria. In order for a spring to be considered a "good compressor" it should remain inside the red area. Consider the spring in Figure 413A to correspond with spring 1 (with the tangential displacement in the opposite direction), and spring 2 corresponds to the spring in Figure 413B. n=l S0.833 TNT = 0.5 Figure 413. Schematic to demonstrate the importance of evaluating the percentage slip of TNT ratio for each spring. A) The TNT ratio in this case is close to 0.5 because the compression and tangential displacements are similar in length. B) The TNT ratio is close to 0.833, because the ratio of compression to tangential displacement is about 1 to 5. TNT = 0 t< TNT t=n I TNT =0.70 t= nx2.41 o 0 0 C: II r 7 TNT =1 t>>n Schematic to shows at the level at which the TNT threshold was defined. Large Caucasian t ( f Large Asian I * I * ** * .* S. **" C. S .< Large inrani. 7 Theta TNT = 0.6 90 750 * ,. t .."600 .. r. ... *i, 450 7 .. 30 N ( 300 TNT = 06.7 150 ' 1 #., m.,, S _. w, . TNT = 0.8 Figure 415. Empirical determination that the TNT threshold greater than 0.6 correlates with the training data, which says that the large masks are stable. Small Caucasian TNT = 0.6 ** N r .. 7....., Unstab i , TNT = 0.8 Small Asian . i Small infant r . j r  : I; f?:'+' + S Er..7 ~ 'a. ? I ''A" Theta *90 750 3 60 45 300 150 6 .: mo eni Figure 416. Empirical determination that the TNT threshold set smaller than 0.8 correlates with the training data, which says that the small masks are unstable. Decima Large Asian **. T= 6 Theta S. *90 S. 750 .. ..600 TNT= '' te 15 Decimate 10 Figure 417. A case of the large Asian mask, which shows the consistency of an upper TNT threshold of 0.7. Even when the number of springs is increase, the TNT threshold of 0.70 shows stability. 80 I N I  ( (. (U S . * *i S. t **r .i. u' I, I. Figure 418. Schematic showing the relation between the volume of the entire frame (fvol) and the compression patch volume (cvol). Instability (>3 "good compressor "points, cvo, = 0) 0.25 0.25 Stability(<= 3 "good compressor" points, cvol o0) 0.75 Box Ratio Figure 419. Box ratio interpretation. Notice that a box ratio of one is the most stable configuration where the compression patch volume (cvol) is equal to the entire volume of the frame (fvol), meaning that all points are "good compression" points. Also notice that a box ratio of zero corresponds to instability; since there are less than three "good compression" points, which result in a compression patch volume equal to zero. 81 i I CHAPTER 5 PROGRAM DESCRIPTION AND VALIDATION Program Description The program was written in C++ language and it uses open source programming libraries for image visualization: VTK (VTK 4.2, Kitware Inc., Clifton Park, NY) and, for advanced imaging processing, ITK (ITK 1.8, Kitware Inc., Clifton Park, NY). The library for advanced imaging processing has several classes dedicated to finite element method (FEM), which were used to construct the simulation of the contact problem. In addition Cholmod [13] was adapted to solve the linear system of equations particular to the FEM; also Paraview (Paraview 1.6, Kitware Inc., Clifton Park, NY) was used for the visualizations of the results. The program is organized into five main sections. The first section is dedicated to reading the input information from the STL file, introducing input parameters and preparing the data to be written into the appropriate input format ITKFEM module. The second section requires the ITKFEM module to read the input file previously created and assemble a single global stiffness matrix (Kg) in the correct format for Cholmod. The third section defines the directions of the force at each point on the rigid layer and iterates over all points while storing into arrays the output parameters. The forth section sorts over the information in the output parameters and applied the single point and global criteria; and finally the last section presents the results, in a histogram format and as a colormap in Paraview. Figure 51 is a chart flow of how the program is organized. The first section of the program starts by reading the STL file and applying an ITK filter that selects the largest region in the STL file. This is done because sometimes additional small surfaces are selected in RPD designer that are not related to the main selected surface. In general the STL surface is characterized as an irregular mesh. The STL serves the purpose well of describing a surface by dividing it in planar triangles. For this reason, another VTK filter, a decimation filter is applied to reduce the number of triangle that are close to each other and have the same dihedral angle [38]. The decimation filter does not alter the topology of the selected surface; rather it serves the purpose of reducing the number of triangles in the STL mesh [38]. However, sometimes the output of the decimation filter results in relatively large triangles, so an optional function was created that added nodes to large triangles and remeshes the surface. At this level in the program the surface does not undergo additional modification, and the nodes on the edges are saved to an array for later identification. The next step is to calculate the normal at every vertex. VTK calculates the normals at every vertex on the surface as an average of the normals of the triangles connected to the vertex; however a function was created to calculate a better average normal, which is weighted by the angles of the triangles adjoining to a particular vertex. With the appropriate normal information a top surface and bottom point cloud are created by extrapolating the point coordinate information along the normal in the positive direction and in the negative direction, respectively, as shown in Figure 52. The rigid layer is created by meshing with tetrahedral elements the space between the selected surface and the top surface (Figure 52). The bottom point cloud become boundary condition points, fixed in all 3dimensions and represents the patient's surface. The spring elements that comprise the bed of springs extend from the point cloud to the selected surface (Figure 52). Next the program calculates the surface area of the frame and with the expected displacement se to 0.001mm, determines the magnitude of the external force (Equation 52). The last step of the first section of the program involves using the point coordinate information and input parameters to write an input file for the ITKFEM module. Most of section two involves using the ITKFEM classes to read the file created, assemble a single the global stiffness matrix (Kg) that describes the entire system and the global force vector (Fg). Afterwards a function was created to read the global stiffness matrix and write it into an input format for Cholmod. The format required by Cholmod is referred as triplet format, and is a method of writing large sparse matrices in manner that saves computer space and quickens the retrieval of information [13]. The direction of forty nine force vectors is calculated in a combination of forty nine tilts (0) and spins (p) about the main normal (Nmain) as shown in Figure 412. The same small magnitude is assigned to each force vector and all forty nine are applied onto each nonedge nodes on the mask. One global stiffness matrix (Kg) is assembled for the entire system and is not altered once it is assembled; however, for each force application to a particular node (49 times the number of nonedge points) on the rigid layer, a global force vector Fg is populated. In this manner, both the global stiffness matrix Kg, and the global force vector Fg are given to the Cholmod routines to solve for the displacement vector U that characterizes the movement of the frame as a result of a single force application. From the displacement vector U, which include displacements for each nodes, the normal (n) and tangential (t) distances are calculated for each node as shown in (Figure 48). In section five, the single point criteria and the global point criteria, as described in Chapter 4, are applied to the tangential and normal scalars. In addition, the program creates histograms for the box ratio per angle rotation and creates a file for visualizing the frame with colorcoded information per angle of the maximum angle in Paraview (Paraview 1.6, Kitware Inc., Clifton Park, NY). Validation of FEM elements The validation of the program involves testing the behavior of the linear spring or space truss element and the linear tetrahedral element to assure they respond in the appropriate manner. In addition, the program is tested for a simple case of a corner frame. Linear Spring Validation The validation of linear spring element was performed by testing the displacement of the element in three orthogonal directions as shown in Figure 53. In all instances the same mechanical characteristics (E, v, A, L) where assigned and the same force vector (F) was applied (Figure 53). The results of the stiffness matrix as assembled by the ITK FEM module are shown in Equation 51, 52 and 53, which correspond to Figure 53A, 53B and 53C, respectively. The resulting displacements for each case resulted in the same value of 0.1mm in the direction of the force vector, which corresponds to the displacement as hand calculated in Equation 54, following the linear spring equation. 750 0 0 750 0 0 0 0 0 0 0 0 kA 0 00 0 00 (51) 0 0 0 0 0 0 (5 0 0 0 0 0 0 0 0 0 0 00 0 750 0 0 750 0 k'B = (52) 0 750 0 0 75 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 750 0 0 750 0 0 750 0 0 750 FL 75N2mm d= = = 0.1 mm (54) AE 0.lmm2 .1500MPa Tetrahedron Element Validation The validation of the tetrahedral element was done by comparing the ITK FEM results to those by ANSYS (Ansys Inc. Canonsburg, PA) following the example in Figure 55. Both ANSYS and ITK FEM provided the same displacement shown in Table 51. A Simple Model Validation The program was validated with a simple frame in the shape of a corner (Figure 55). When a force was applied at the location marked with a white dot and in the direction of normal direction at that point; the entire corner was placed in a state of compression as will be expected if a force physically applied onto a corner (Figure 55A). In the reverse case where the force vector is inverted and rather than pull it pushes away, the corer results in a complete tension state. Table 51. Displacement in the X, Y and Z direction from the example in Figure 54 tested with the ITKFEM module and ANSYS. Node X Y Z 1 0.0000 0.0000 0.0000 2 0.0000 0.0000 0.5200e3 3 0.0000 0.0000 0.0000 4 0.0000 0.0000 0.0000 5 0.0000 0.0000 0.5200e3 Figure 51. Flow chart of program for testing the uniqueness of fit. 88 Get main frame normal (Nmain) For all nodes (not edge nodes) For 0 = 0 to 900, p = 0 to 3600 (49 force vectors/node) Rotate force vector and assign magnitude Insert new force vector into Fg Solve Fg=KgU with Cholmod Get U and calculate: normal and tangent displacements K Yes Y No OUPUT to histograms and frame colormap 1~~   "" Section V Figure 51. Continued. Section III Section IV Apply Criterion Single Point Criteria (for each point on surface) Calculate: TNT, #T, #C If (#C > = 3 && TNT>=0.75) {GoodCompressor++} Test for Flags: Positive Definiteness and RDR Global Point Criteria (for every force application) If (GoodCompressor >=3) {Calculate Box Ratio AngleRotateCount++} I E I End  Rigid tetrahedral / Selected surface layer Bottom point cloud Figure 52. Original surface is extended in the positive normal direction to form the top surface and the points extended in the negative normal direction to become boundary condition points. E =1500 MPa v = 0.3 F=75N A= 0.1 mm2 L = 2mm Element 1 < > L Cross Sectional Area (A) /Y 4 Y Figure 53. Examples of al linear spring element or space truss tested in all 3dimensions.  30e6 MPa 0.3 1000 N Figure 54. Example of two tetrahedral elements used in the validation of the tetrahedral elements used in the program. Node Coordinates (X,Y,Z) 1 0,0,0 2 0,10,0 3 0,0, 10 4 1,0,0 5 1, 10, 0 Force application point Force direction 1350 B Figure 55. A corner frame that was used to validate the program. A) The corer frame is shown on blue and the white dot signifies the location of the force. B) The direction of the force is the normal located at 1350 from all three sides of the corner. C) The wireframe arrows show a scaled normal displacement vector for all nodes on the frame, and the small solid arrows show the displacement vector at each node. D) The force application in the opposite direction results in all normal displacements in the opposite direction to C. CHAPTER 6 RESULTS AND DISCUSSION A set of criteria were develop that showed that all large masks were found to be stable and all small masks were not stable. Also, the caplike reference frame and the half sphere surfaces were found to be not stable. In reality, all masks maintained stability when loaded in the direction and location of each mask's main normal (Nmain). In reality the conforming surfaces of all small masks do not have a unique fit, meaning that there is not one unique location where the mask fits, rather there are various locations were it could fit (Figures 61B, 62B and 63B). The small masks were designed with the expectation that they will not have a unique fit, for the purpose of validating the methodology. For all large masks there is a unique fit (Figures 61A, 62A and 63A). All large masks can withstand real lateral forces or forces 900 from their main normal (Nmain). The Asian small mask is the only small mask that can sustain contact stability when a real load is applied at large angles from its main normal (about 600 from Nmain); it is the only small mask that extends laterally along the entire superciliary arches (Figure 62B). It also extends caudally onto the nasion and sits onto the upper orbital area. The small Caucasian mask does not extend onto the orbital area and does not extend onto the nasion (Figure 61A). The small infant mask extends laterally along the superciliary arches but it does not reach over the nasion (Figure 6 3B). The facial angles of all three anatomical models were measured as defined by Figure 211; the results are shown in Table 61. One can appreciate that stability could be sustained if the nasiolateral and nasiofrontal angels are relatively small and if the nasiofacial angle is large. With smaller nasiolateral and nasiofrontal angles there is more opposing contact area to lateral forces (900 from the Nmain). The Caucasian model has the smaller nasiolateral angle, and the second smallest nasiofrontal angle and the largest nasiofacial angle (Table 61). For the simulation results, two types of graphs were produced, a) a color coded surface plot of the masks onto the solid models (Figure 65 to 69) and b) a normalized histogram for the box ratios for all values of p (p = 0 to 3600), for each tilt angle (0) (Figure 512) accumulated over all points (Figure 610 to 617). In the case of the surface models, the results are presented for each force application point where the color signifies the maximum tilt angle for which stability was found (based on the singlepoint and single stability criteria) for a full rotation (p = 0 to 3600). At each color coded point on Figure 65 and 69, a set of 49 force vectors were applied. For each force vector all springs that comprised the bed of springs were evaluated for stability. At each spring the maximum angular force (for a full rotation p = 0 to 3600), the furthest from the main normal (Nmain), at which at least 3 springs were found in compression was noted. This maximum angle is color coded in Figures 65 to 69. If the frame was found stable at every force application it will be indicated in red, demonstrating that for any force between theta equal to 0 and theta equal to 90, the mask did not slide or tilt beyond the criteria established (TNT >= 0.7). Some masks (small unstable masks) experienced some tangential movement or rotation, so they only exhibit stability at a lesser angle (i.e. theta = 45 from the main normal). For the box ratio histogram plots, the results are presented as frequency of the box ratio by individual tilt angle (theta or 0) for all p values over all points on the mask (Figures 610 to 6 17). For example, in Figure 610 for the Caucasian large mask, for a tilt angle of 0, 30% of the box ratios are equal to 0.4. For the Caucasian large mask (Figure 65A), the surface plot shows that all nonedge points sustain stability for a maximum theta of 90 and for all rotations about the Nmain (p = 0 to 3600); while the small Caucasian mask's design did not sustain contact for any full rotation for any point (Figure 65B). For the Asian large mask (Figure 66A) and for the infant large mask (Figure 67A), the surface plots show similar results to the large Caucasian mask (Figure 65A), all points sustained stability all the way to a tilt angle of 900. It is important to notice that fewer forces were applied on the large Asian and large infant masks because there are less nonedge points (Figures 66 and 67). The surface plot for the small Asian mask shows more stability per location than any other small mask. This mask extends along the superciliary arches and sits on top of the upper orbital area and extends caudally past the naison area providing a surface to judge correct alignment. The color coded surface plot for the small Asian mask agrees with the design as they show stability along the left lateral side and over the upper orbital area. On the right side less stability was found simply because the mask is not symmetric and includes less coverage over the right upper orbital area. The color coded surface plot for the small infant mask shows some point locations where the tilt angle theta sustains stability up to 600 (Figure 67). This is due to the mask extending over the upper orbital area and laterally along the superciliary arches. The single force criterion states that three points need to be in compression as defined by the single point criteria (with a TNT ratio of 0.70 or below) to have stability at a particular force. As mentioned before, if 3 or more "good compressor" springs are found, then the compression patch volume (cvo) is different than zero and the box ratio results in a value greater than zero and less than or equal to one (Figure 419). The level of stability is the highest when the box ratio is close to 1, since this indicates that the volume of the compression patch (cvo) is similar to the volume of the entire frame (fvol). The box ratio histograms for the small Caucasian mask (Figure 610) shows highest frequencies at box ratio equal to zero, this means that for most external force applications the mask was found unstable (from the single force criterion). These results correspond to the surface plot results in Figure 65B. The box ratio histogram for the small Asian mask shows a range of box ratios between 0.0 and 0.55, with most frequencies (70%) at box ratios of 0.15 (Figure 613), this means that when stability was found the compression patch volume (cvol) was small. However, for all theta values there are frequencies where the box ratio is equal to zero, this means that even though stability was found in most of the cases, there were always instances were the mask lost stability. This results are harder to relate to the surface color maps, because the as mentioned before, the surface color map show stability per location, and the box ratio histograms show stability by accumulating all forces for a particular theta value over all rotations. The infant small mask sustained stability for theta equal to 0, but had no stability for all other angles (Figure 615). At theta equal to 900, 80% of the of the box ratios are equal to zero, which signifies a high level of instability for forces applied at 900 from the main normal. The box ratio histograms for the large Caucasian mask shows no instability (zero frequency for box ratios equal to zero), and the box ratios range between 0.25 and 0.95 (Figure 610). As the theta angle (tilt angle from the Nmain) increases the box ratios decreases, showing that at higher angles the compression patch volume is smaller. Similar box ratio histogram results can be seen for the Asian large mask (Figure 612) and the infant large mask (Figure 6 14) where there is no instability or zero frequency for box ratios equal to zero. However, for the large Asian mask the box ratios are skewed to the left, they range from 0.25 to 0.65 (Figure 6 12). This means that for the large Asian mask, increasing theta values result in smaller compression patch volumes (cvol) than for the large Caucasian mask. For the large infant mask the results follow the same pattern, increasing the tilt angle of the external force results in decreasing box ratios (Figure 614). The range for the large infant mask is from 0.15 to 0.75, this is a slightly larger box ratios range than for the box ratio range for he large Asian mask (0.25 to 0.65) (Figure 612). Based on the facial angle measurements in Figure 211 and on the premise that the more pronounced facial features result in more stable masks, the ideal angle combination for a stable mask will be small nasiolateral and nasiofrontal angles and large nasiofacial angles. Out of the three large masks, the Caucasian mask corresponds to a model will the smaller nasiolateral angle and the second smallest nasiofrontal angle and the largest nasiofacial angle. In addition, the box ratio range for the Caucasian mask is the largest (0.25 to 0.95) when compared to the Asian (0.25 to 0.65) and infant (0.15 to 0.75). Figure 613 shows a profile view of all three models. The infant large mask corresponds to a model with the smallest nasiofrontal angle and the second smallest nasiolateral angle, but the smallest nasiofacial angle. The infant model has the smallest nasiofrontal angle due to a curvature on the forehead as can be seen in Figure 618C. It is likely that this pronounced curvature is responsible for the infant large mask having the second highest box ratios ranging between 0.15 and 0.75. In addition to the mask reference frames, a caplike reference frame (Figure 68B) and a half sphere surface (Figure 69) were evaluated with the methodology for testing the uniqueness of fit. In reality the caplike frame has an uncertain fit, in other words, it was difficult to find the location where it was designed to fit. The color coded surface plot results for the caplike reference frame shows that the frame sustains stability (as defined by the single point and single force criteria) all the way to a maximum theta value of 90 when forces are applied on the upper most cephalic region. However, forces no greater than 60 where sustained at the coronal and sagittal extensions of the frame (Figure 68). The box ratio histogram (Figure 616) for the cap like frame corresponds with the information on the surface plots. The box ratio histogram shows that some level of stability was found for all theta angles, but it was not sustained through out the entire surface, since about 5% of the box ratio show instability for all tilt angles. The remaining box ratios (about 70%) have values of 0.15, which signifies that the compression patch volume (cvol) was small compared to the volume of the entire mask (fvol) for the instances were stability was found. In the case of the halfsphere surface (Figure 69), the surface plot shows that there were no locations where the surface sustained any force, possibility not even in the direction of the normal force (Nmain). The box ratio histogram for the half sphere test shows that stability was indeed found for all angles (there are box ratios between 0.Oand 0.35); however, for the majority of the cases (60% to 80% frequency) no stability was found or the box ratio equaled the value of zero. In summary, Table 62 shows the instability scores for all masks and frames as defined in Equation 48. There is zero percent instability for all large masks and there is some percentage of instability for all small masks, the caplike frame and the halfsphere surface, as was expected. The least stable case is the halfsphere test, with a instability score of 94%, this high percentage instability make sense because the surface has no irregularities to oppose any external forces (besides in the main normal direction). The small Asian mask was the mask with the smallest instability score at 5%, as was mentioned previously; this is due to the design of the mask that includes the upper orbital area on one side and extends the furthest medially and laterally along the superciliary arches. This results correlate well with the color coded surface plots, were stability was found in the left upper orbital region on the mask. Table 61. Measurement of the facial angles for the three anatomical models. Model Nasiolateral Nasiofrontal Nasiofacial angle angle angle Caucasian 890 1160 400 Asian 1230 1240 290 Infant 1040 1060 230 Table 62. Instability score, presented as the percentage of forces that resulted in instability and the total forces applied on a frame. Frame Total Number of forces Total Number of Percentage applied forces applied Large Caucasian 4935 0 0% mask Large Asian mask 1575 0 0% Large infant mask 1702 0 0% Small Caucasian 874 602 69% mask Small Asian mask 1755 88 5% Small infant mask 874 334 38% Caplike frame 2695 241 9% Halfsphere test 1094 975 94% Figure 61. Placement of Caucasian masks onto solid model. A) Large Mask. B) Small Mask. Figure 62. Placement of Asian masks onto solid model. A) Large Mask. B) Small Mask. Figure 63. Placement of infant masks onto solid model. A) Large Mask. B) Small Mask. Figure 64. Placement of a caplike reference frame onto an upper head model. Theta 90 75" 600 450 30 15" I'" ;jm jj~ B4 B Figure 65. Surface plot color coded by the maximum full rotation angle for the Caucasian masks. A) Large Mask. B) Small Mask. Theta 0 E900  A 4750 600 B S150 moo A B Figure 66. Surface plot color coded by the maximum full rotation angle for the Large Mask. B) Small Mask. Asian masks. A) .6' ^,^ tt ! ItV Theta *90 750 600 450 300 150 m0o NiONv a't 7 Ji__ Figure 67. Surface plot color coded by the maximum tilt angle for the infant masks. A) Large Mask. B) Small Mask. Theta 750 600 300 150 7 A m"o Figure 68. Surface plot color coded by maximum tilt angle for the cap reference frame design B S * .... * a. . 4 Theta E90 750 d 600 ,. 450 0 150 moo * Figure 69. Surface plots color coded by maximum tilt angle for the half sphere test. ThetaO STheta15 100% O Theta30  O Theta45 90% Theta60 80% Theta75 70% Theta90 60% S50% 240% 30: 209 103 03 0.00 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 Box Ratio Figure 610. Normalized histogram of box ratios for each tilt angle (0) for the Caucasian large mask. 104  60% 50% 40% U 30 20 10 C'1 0.00 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 Box Ratio Figure 611. Normalized histogram of box ratios for each tilt angle (0) for the Caucasian small mask. 0.00 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 Box Ratio Figure 612. Normalized histogram of box ratios for each tilt angle (0) for the Asian large mask. 60%  S50% 40% U 30  20 1C0 CI 0.00 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 Box Ratio Figure 613. Normalized histogram of box ratios for each tilt angle (0) for the Asian small mask. 100% 90% 80% 70% 60% 50% 40% 30 20 100 00O 0.00 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 Box Ratio Figure 614. Normalized histogram of box ratios for each tilt angle (0) for the infant large mask. 100% 90% 80% 70% 60% 50% 40% 30  20 10 C11 A 0S Theta0 ES Thetal5 SS Theta30 O S Theta45  SS Theta60 SS Theta75 SS Theta90 0.00 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 Box Ratio Figure 615. Normalized histogram of box ratios for each tilt angle (0) for the infant small mask. 60% 6 E 50% 40% S 3C, 10 CI 0.00 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 Box Ratio Figure 616. Normalized histogram of box ratios for each tilt angle (0) for the cap reference frame design. S60% c 50% 40% U 30 20 10 C 0.00 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 Box Ratio Figure 617. Normalized histogram of box ratios for each tilt angle (0) for the half sphere test. Figure 618. Profiles of solid models. A) Caucasian. B) Asian. C) Infant. CHAPTER 7 CONCLUSION AND FUTURE WORK The results have determined a clear difference in stability of fit between the large masks and the small masks for all cases. Both the surface plots and the box ratio histograms show that the large masks are stable while the small masks are not. Furthermore, the surface plots are able to show which force location points and what angles result in an unstable fit. The box ratio histograms are able to quantify stability on a tilt angle basis; by assigning a normalized numerical value to the magnitude of the contact volume, if three compression points are present that satisfy the TNT threshold. In addition, the box ratio histograms are able to provide information about the maximum tilt angle (theta) that the entire mask can withstand over all rotations about the main mask normal. The two additional cases, the caplike reference frame and the half sphere surface have further validated the program by providing expected results when compared to predicted real space frame movements. Furthermore, these two additional surfaces have shown that the surround coverage of the head does not particularly results in stability, and hints that the facial landmarks are the best locations to uniquely fit a mask reference frame onto a patient. The training data was compared to the three models' facial landmarks and motivating relations about stability and facial angles were found, which suggests that if masks are built for patients with small nasiolateral and nasiofrontal angles and large nasiofacial angles, their masks may be highly stable. In the future more masks should be tested with the parameters defined in the methodology. Also it could be worth while to measure the facial angles of more patients and run statistical analysis to compare the relations between facial angles and stability as measured by the methodology defined in this investigation. Concerning the program, a meshing filter could be beneficial for regularly meshes the surfaces. Also the box ratio definition could be redefine for easy understanding, by simple calculating the compression patch volume (cvol) divided by the entire frame volume (fvol), this will result in values from zero to one, in which zero signifies instability and one signifies full stability. The present investigation has provided with a methodology for testing uniqueness of fit for the patientspecific reference frames that improves the new and promising technology to perform imageguided surgery. The methodology will provide the clinician with prior information of frameonpatient stability, which was not available before and is essential for allowing the clinician to effectuate his surgical plan with more confidence and accuracy. In addition the methodology was created using an open source image analysis software (ITK and VTK), and by doing so, it will be easily adapted to the already existing user interface (RPD designer also was created with ITK and VTK), which allows the surgeon to create his/her own reference frames based on diagnostic images. By including the stability measuring module in the already existing graphical user interface (RPD designer), the new technology for imageguided surgery will be more robust and a step closer towards being introduced to other medical communities. APPENDIX A PRINCIPLE OF VIRTUAL WORK AND ITK FEM MODULE This section provides additional background information related to the finite element method (FEM) for two main reasons. The methodology for testing the stability of fit of the patientspecific frames is based on FEM and the tools used for performing the FEM lack documentation. This chapter first provides background information on FEM, and then uses the direct stiffness method to explain the basic steps to derive the element stiffness matrix for a 2 dimensional truss element. It also comments on the principles of virtual work because the ITK FEM module (ITK 1.8, Kitware Inc., Clifton Park, NY) follows them to arrive to one of the element's stiffness matrices used in this investigation. Finally it discusses the ITK FEM module and its organization in an attempt to supplement for the currently available documentation. FEM Terminology for Truss Elements There is a specific technical terminology associated with FEM. The following list [4,28,46] presents the main components involved in FEM as it relates to structural analysis and will aid in the understanding of the subsequent sections. * Discretization: The process of subdividing a structures into members or finite elements with the purpose of modeling, simulating, or approximating their mechanical behavior under known forces. * Elements: In the case of structure analysis, an element is a single member of a structure. Its contribution to the structure is defined by a forcedisplacement relationship and this relationship is used to analyze the element. * Node: The purpose of the node is to connect elements, it is assumed to be infinitesimally small in size and therefore frictionless. Depending on the structure or part simulated the nodes can behave differently. For example, a node will act as a frictionless hinge in 2 dimensional truss elements, and as a frictionless balland socket in 3dimensional truss elements. * Global and Local Coordinates: In FEM, usually two coordinate systems are used. One is referred as the global coordinate system (X, Y, and Z) and it describes the overall geometry and loaddeformation relation of the entire structure. The other is referred as a local coordinate system (x, y and z) and it describes the forcedisplacement relations for each individual element. In the case of the truss element, the local coordinate system's origin is located at one end of the truss element and the local xaxis is oriented along the axis of the element. * Degrees of Freedom (DOF): Degrees of freedom represent node displacements, such as translations and rotations. When a structure is subjected to an external load, the degrees of freedom are the independent node displacements necessary to indicate the deformed shape of the structure. A 1dimensional truss element has 1 DOF per node, a 2dimensional truss element has 2 DOF per node and a space truss has 3 DOF per node. For truss elements the nodes are assumed to be frictionless joints, so the elements are not subjected to moments, and the rotations are zero. * Boundary Conditions (BC): Boundary conditions are applied at nodes to fix the model in space or to set a prescribed displacement onto the nodes. In a structural model, they represent fixed supports necessary to remove the rigid body motion of a structure. In addition, when the boundary conditions are set correctly, the total number of unknowns must be equal to, or smaller than, the number of equation necessary to describe the model under known loading conditions [4]. * External forces and internal forces: External forces are specified in global coordinates and internal forces are specified in local coordinates. Local coordinates are induced onto elements by external forces. For example, an external force is applied on a particular node on a structure, as a result the elements that compose the structure will deform and internal forces will be induced at each individual element's node. * Material Properties: Material properties define elastic characteristics of the elements. For simple truss elements they are the elastic modulus E (MPa), and the Poisson ratio v [3]. * Transformation Matrices (T): For simplification, the integration of the element is performed in local coordinates. A transformation matrix is need that allows the transition from local to global (and vice versa) coordinates. It involves the direction cosines of the local coordinate system in the global coordinate system. * Element Stiffness Matrices (ke): The stiffness matrix of an element relates the local forces to the local nodes' displacement. The element stiffness matrix is composed of stiffness coefficients. The stiffness coefficient denoted ki is explained as the force in the i direction required to cause a displacement at node j while all other node displacements are zero [4]. Principle of Virtual Work In the FEM the forcedisplacement relations can be derived from the workenergy principles; the main reason for this being that by using workenergy principles the force displacement relations can be conveniently applied to solid elements. The principle for virtual work for rigid bodies states [4]. If a rigid body, which is in equilibrium under a system of forces (and couples), is subjected to any small virtual rigidbody displacement, the virtual work done by the external forces (and couples) is zero. A virtual displacement is simply a small imaginary displacement. The principle of virtual work for rigid bodies implies that if the rigid body is in equilibrium, then the sum of forces and moments must add to zero, thereby causing the virtual work, 6We, to be equal to zero, or 6We = 0. The principle for virtual work for deformable solids state [4]. If a deformable structure, which is in equilibrium under a system of forces (and couples), is subjected to any small virtual displacement consistent with the support and continuity conditions of the structure, then the virtual external work done by the real external force (and couples) acting through the virtual external displacements (and rotations) is equal to the virtual strain energy stored in the structure. The principle mentioned above implies that virtual external work is equal to virtual internal work, and since internal work is equal to the structure strain energy then 6We = 6U. In this investigation we are using only linear elements. FEM linear analysis follows two basic assumptions [4]. * The part or model is simulated with a linear elastic material, this means that the strainstress relations behave in a linear manner and follow Hooke's Law, which states that the stress a (MPa) is proportionally related to the strain E (unitless), by the elastic modulus of the material denoted as E (MPa), or c=ES. * The model experiences deformations that are infinitesimally small; therefore the squares and higher powers of the element's deformed shapes are negligible in comparison to unity. Indeed, the equation of equilibrium can be based on the non deformed geometry of the part. ITK FEM and Principles of Virtual Work The following section is based on information from the ITK libraries [24], Aslam et al. [4] and Felippa et al. [17, 18] and has the purpose of providing background information on how the ITK FEM module sets up the element stiffness matrix, because the module lacks explanations and clear documentation. The ITK FEM module follows principles of virtual work to generate the element stiffness matrix in local coordinates. The ITK FEM module requires 1) the definition of the shape functions for each element types, 2) the derivation of the shape functions (strain displacement matrix) and 3) values related to the numerical integration of the element stiffness matrix. To get a brief idea of how these components relate to the principles of virtual work first consider a 3dimensional differential element in equilibrium and subjected to loading (Figure A 1). Following the principle of virtual work for deformable solid bodies [4], the loading on the differential element results in real stresses (o), which cause in virtual strains (6S) and internal forces that perform virtual work. So 6b = a = (Al) Yxy xy zx zx Equation Ai shows the notation for the stresses (a and z) and virtual strain that define the differential element under loading and in equilibrium. The total virtual internal work (or virtual strain energy) stored in the element can be measured by integrating over the volume [4], as shown in Equation A2. aw, = J(&xx + a + ,C + xy T y + yz, +r yzx iV (A2) V S = ISWadV (A3) Equation A2 and Equation A3 (compressed formed) defines the principle of virtual work for deformable bodies expressed in terms of stresses (a) and virtual strains (5e), they represent the virtual strain energy stored in the element. Now the focus becomes defining the local element stiffness matrix ke, in terms of stresses and strain using principles of virtual work of deformable bodies. The first step is to define the node displacements in terms of stresses and strains. Displacement Functions The node displacement for a truss element can be represented by a displacement function with the form of a linear polynomial (Equation A4) [11]. The displacement function or displacement variation i(x) needs to be defined solemnly along the local xaxis. U(x) = ao + ax (A4) The roots of the polynomial, ao and ai, are constants and must satisfy the boundary conditions, that is the displacement conditions at the nodes [11]. From Figure 21, when x is equal to 0 (at the origin of the local coordinate system), the displacement variation function i(x) is equal to ul. (O)= u = ao (A5) In the same way, when x is equal to L, i(x) is equal to u2. u(L) = u2 = ao + aL (A6) The displacement variations of the truss element in terms of node displacements, ui and u2 can be calculated by substituting Equation A5 and A6 into A4 and rearranging to yield Equation A7 [11]. iT(x) I= 1 U + kU2 (A7) L) L Shape Functions The displacement function (Equation A7) can be alternatively expressed as shape functions. There are as many shape functions for an element as there are nodes, and they have the characteristic of having a value of one at its corresponding node (N1 has a value of one at node 1) and a value of zero at the other node points [3]. For a truss element there are two shape functions denotes by N1 and N2. The roots of the polynomials in Equation A7 represent the shape functions (Equation A8). N, = ) N2 =Q (A8) For convenience the displacement variation function (i(x)) can be expressed in terms of shape functions and node displacements (Equation A9). Equation A10, is the final form or the element's displacement function (or interpolating function) and N is the element's shape function matrix, and u are the element's node displacements. [ix]= [N, N2 1 (A9) LU2 i Nu (A10) Strain Displacement Equation Getting back to the aim of representing the node displacement in term of stresses and strain, Equation A11 shows the relationship between strain (E) and displacement variation function (u(x)), for a truss element [11]. E = (A11) dx From Equation A 1, notice that a truss element only withstand strains along the local x axis. S= [d [ ] = Do (A12) dx D is called the differential operator matrix [4]. At last the strain can be related to the node displacements u by substituting Equation A10 into Equation A12. S= DNu =Bu (A13) In order to get the strain displacement matrix, B, the derivatives of the shape functions need to be calculated (Equation A13). The node displacements do not depend on the local x direction, so they are treated as a constant, while the differentiation can be applied to shape functions N and not to the node displacements u [4]. Following Hooke's law, the stresses can also be related to the end displacements (Equation A14), as required by the principle of virtual work for deformable bodies [4]. a = EBu (A14) So far both strains (Equation A13) and the stresses (Equation A14) are expressed in terms of node displacements. Local Element Stiffness Matrix All the necessary information required by the ITK FEM module is now available to determine the element's stiffness matrix, with the exception of a relation between the element end forces (fi and f2 as in Figure 31) and the node displacements u. When subjecting a truss element to the principles of virtual external work ( We), the truss element will be in equilibrium and subjected to end forces. These end forces will manifest themselves as virtual end displacements (6ul and 6u2). Equation A15 shows how the virtual external work can be defined in terms of virtual end displacements and end forces [4]. W, = f,5ug + f2 u2 (A15) W, = Su'f (A16) In symbolic form Equation A16 6uT is the transform of the virtual end displacements and f is the end force vector. Equation A16 can be substituted into Equation A3 to get Equation A 17. SuTf = IETdV (A17) v A relation between the element stiffness matrix ke and the strain displacement matrix B is finally accomplished by substituting Equations A11 and A12 into Equation A17 and simplifying by canceling 6uT (Equation A18) [4]. f = BTEBdV u= keu (A18) The element stiffness matrix ke is determined by virtual work principles as defined by Equation A19 [4]. ke = BEB dV (A19) Gauss Quadrature To evaluate the integration in Equation A18, the ITK FEM module uses a method of numerical integration called the Gauss quadrature (Equation A20) [18]. n ke = wJB EB, (A20) 1=1 In this equation, n is the number of integration points, i is the integration point index, w is the integration weight, B is the strain displacement matrix, and Ji is the jacobian determinant evaluated at the integration point. In the ITK FEM module all these parameters are necessary to define an element [18]. The determination of the element stiffness matrix is the fundamental part of the ITK FEM process; afterwards the element stiffness matrix needs to be transformed into global coordinates and assembly into the stiffness matrix that defines the complete system [28]. In the ITK FEM module after the element stiffness matrix is constructed for a particular element and transformed into the global coordinates, the remaining steps of assembly and reduction are operations that are not element dependent. The purpose of the past section was to explain the basic of parameters in referenced to the ITK FEM module to complement for the lack of documentation. ITK FEM Module ITK is open source software that provides data representation and algorithms for implementing image analysis, primarily segmentation and registration [24]. Its origins date back to 1999 when the US National Library of Medicine of the National Institute of Health granted a threeyear contract to create an opensource registration and segmentation toolkit, later it came to be know as the Insight Toolkit or ITK. The software is written in C++ and is greatly object oriented and its organization represents a pipeline where data and process are connected together. ITK has an FEM module primarily intended to be used for image registration, but complete to perform mechanical modeling. ITK provides an objectoriented FEM library to solve general FEM modeling. Classes are used to specify the geometry, behavior of the elements, apply external forces, apply boundary conditions, specify elasticity, and to solve the problems (Figure A2). The core of the library is the element class, which is divided in two parts: the definition of the element geometry, and the physics of the elements or the behavior of the element. The outcome of combining the two parts results in the creation of a new element. Figure A2 shows a section of the ITK FEM element class. Notice that from the Element class (itk::fem::Element) the space truss is defined by itk::fem::ElementStd<2,3> and the linear tetrahedron by itk::fem::ElementStd<4,3>. The first number in the class name implies the number of nodes in the element and the second number implies the dimension in which the element resides. Following Figure A2B, itk::fem::Element3DCOLinearLine is a class derived from ElementStd (ElementBaseClass), this class defines the geometry of the element by defining the following parameters. These parameters relate to the following the FEM by principles of virtual work. Integration Point and weight Number of Integration Points Shape Functions Shape Function Derivatives Get Local From Global coordinates Jacobian In the class itk::fem::Element3DStrain templates class that is used to define an element in 3D space. Element3DStrain together with Element3DCOLinearLine define the element in full. itk::fem: :Elemenet3DCOLinearLineStrain provides access to the following functions that sets the element stiffness matrix. Get Strain Displacement Matrix Get Material Matrix Get number of degrees of freedom per node The following list [24] of bullets further explains Figure A2. * itk::fem::LightObject: base class that defines the entire FEM library. * itk::fem::Element: Abstract base class created as the parent structure, like any base class no objects are created here; subclasses are derived and used to create the objects. * itk::fem::ElementStd: Templated class that automatically defines the sum of the virtual functions in the elements. Three template parameters are defined here: number of nodes, number of spatial dimension, and the template class from which this class is derived. * itk::fem::ElementNode: this class stores information required to define a node, like the node coordinates. * itk::fem::Load: Abstract base class that defines external loads that act on the FEM system. * itk::fem::LoadElement: Subclass that defines the external loads that act on the element. A force vector is defined; it has pointers that point to the element where the load acts. * itk::fem::LoadBC: Class that applies essential or dirichlet boundary conditions, specifies which degrees of freedom are fixed in the model. * itk::fem::LoadNode: Class used to specify a load on a point within an element object. A pointer to an element object, the number of the node on which the load acts, and the force vector. * itk::fem::Material: Base class that defines the element's material properties, * itk: :fem::MaterialLinearElasticity: Class that defines the material properties of the finite elements, these are: E (elastic modulus), A(area), I (moment of interia), nu (poisson ratio), h (thickness), and Rhoc (density multiplied by heat capacity). Figure A1. Differential element subjected to real stresses and in equilibrium. itk::fem::FEMLightObj ect I I I I itk::fem::Element itk::fem::Element::Node itk::fem::Load itk::fem::Material S r C i itk: :fem: :Element3DCOLinearLine i  itk: :fem: :Element3DStrain CitC::C CP CP COiCn iO etri .......................... .... ... I ... .... .. ................ .......................... ..... itk: :fem: :Element3DCOLinearLineStrain r I itk::fem: :Element3DCOLinearTetrahedron I .I ,...................................... ......................... Sitk::fem::Element3DStrain itk::fem: :Element3DCOLinearTetrahedronStrain Figure A2. The greater part of the ITK FEM library. The boxes surrounded by dashed lines represent where the geometry of the element is defined. The boxes surrounded by dotted line, represent where the physics of the problem is defined. Together they entirely define the element, in this cases the Element3DCOLinearLineStrain and the Element3DCOLinearTetrahedronStrain [24]. LIST OF REFERENCES 1. Adler JR. Surgical guidance now and in the nuture. The next generation of instrumentation. Clinical neurosurgery. Proceedings of the Congress of Neurological Surgeons. San Diego, California 2001. Ed Lippincott Williams & Wilkins 2002;105114. 2. Alker G, Kelly PJ. An overview of CT based stereotactic systems for localization of intercranial lesions. Computational Radiology 1984;8:193196. 3. Allaire PE. Basics of the finite element method. Iowa: Wm. C. Brown Publishers, 1985. 4. Aslam, Kassimali. Matrix analysis of structures. California: Brooks/Cole Publishing Company, 1999. 5. Becker AA. An introductory guide to finite element analysis. New York: ASME Press, 2004. 6. Birnbaum Klaus, S.E., Decker Nils, Prescher Andreas, Ulrich Klapper, Radermarcher Klaus. ComputerAssisted Orthopedic Surgery with individual templates and comparison to Conventional Operation Method. Spine 2001;26:365370. 7. Bova, F. Grant proposal: Imageguided surgery based on rapid prototyping models. Department of Neurological Surgery, University of Florida. Gainesville, FL. March 2003 8. Brown RA: A computer tomographycomputer graphic approach to stereotactic localization. J Neurosurg 1979; 50:715720. 9. Brown RA: .A stereotactic head from for use with CT body scanners. Invest Radiol 1979;14:300304. 10. Cheney ML. Facial surgery. Plastic and reconstructive. Baltimore: Williams & Williams, 1997. 11. Cheung YK, Yeo MF. A practical introduction to finite element analysis. London: Pitman Publishing Limited, 1979. 12. Cowin, S. C. Bone Mechanics. New York: CRC Press, 2001. 13. Davis, T A. User Guide for Cholmod: A Sparse Cholesky Factorization and Modification Package. Version 1.4. Gainesville: Department of Computer and Information Science and Engineering, 2006. 14. D'Urso PS, Atkinson RL, Lanigan MW, Earwaker WJ, Bruce IJ, Holmes A, Barker TM, Effeney DJ, Thompson RG. Stereolithographic (SL) biomodelling in craniofacial surgery. British Journal of Plastic Surgery 1998;51:522530. 15. D'Urso PS, Earwaker WJ, Barker TM, Redmond MJ, Thompson RG, Effeney DJ, Tomlinson FH. Custom cranioplasty using stereolithography and acrylic. British Journal of Plastic Surgery 2000;53:200204. 16. D'Urso US Patent # 5,752,962. 1998. 17. Felippa Carlos. Introduction to finite element methods 2006. Department of Aerospace Engineering Sciences, University of Colorado at Boulder. Last Accessed 30 Jan 2006. http://www.colorado.edu/engineering/CAS/courses.d/IFEM.d/ 18. Felippa Carlos. Nonlinear Finite Element Methods 2006. Department of Aerospace Engineering Sciences, University of Colorado at Boulder. Last Accessed 30 Jan 2006. http://www.colorado.edu/engineering/CAS/courses.d/NFEM.d/ 19. Gildenberg PL. Sterotactic Surgery: Present and Past. In Heilbrun MP, ed. Stereotactic Neurosurgery Vol 2: Concepts in Neurosurgery. Balitmore:Williams & Williams, 1988.:1 16. 20. Goffin J., V.B.K., Vander Sloten J., Van Auderkercke R., Smet M. H., Marchal G., Van Craen W., Swaelens B., Verstreken K. 3DCT based, personalized drill guide for posterior transarticular screw fixation at C1C2: Technical Note. NeuroOthopedicsl999;25:4756. 21. Goffin Jan, Van Brussel Karel, Martens Kristen, Vander Sloten Jos, Van Audekercke Remi, Smet MariaHelena. Threedimensional computer tomographybased, personalized drill guide for posterior cervical stabilization at C1C2. Spine 2001,21;12:13431347. 22. Heilbrun MP. Computed tomographyguided stereotactic system. Clinical Neurosurgery 1984;31:564581. 23. Horsley V, Clarke RH: The structure and function of the cerebellum examined by a new method. Brainl908;31:45124. 24. Ibanez L., Schroeder. W., Ng L, Cates J. The ITK Software Guide.NewYork: Kitware Inc.,2003. 25. Loftis KL, Geer CP, Daneslson KA, Slice DE, Stitze JD. Analysis of pediatric head anthropometry using computed tomography for application to head injury prediction. 2007;43:8690. 26. Mackerle J. Finite element modeling and simulations in orthopedics: a bibliography 1998 2005. Comput Methods Biomech Biomed Engin2006;9(3): 149199. 27. Maciunas RJ, Galloway RL, Latimer JW. The Application Accuracy of Stereotactic Frames. Neurosurgery 1994;35(4):682695. 28. Kim NH., Sankar B V. Lecture notes: EML 4500. Finite Elements analysis and design. Gainesville: Department of Mechanical & Aerospace Engineering, 2004. 29. Pham DT, Gault RS. A comparison of rapid prototyping technologies. International Journal of Machine tools and Manufacturing 1998;3 8:12571287. 30. Pikley WD, Pikley OH. Mechanics of Solids. New York: Schaum, 1974. 31. Popov, E P. Introduction to mechanics of solids. New Jersey: Prentice Hall,1968. 32. Raabe A, Krishnan R, Zimmermann M, Seifert V: Frameless and Framebased stereotaxy? How to choose the appropriate procedure. Zentralbl Neurochir.2003;64(1): 15. 33. Rajon DA, Bova FJ, Bhasin RR, Friedman WA. An investigation of the potential of rapid prototyping technology for imageguided surgery. Journal of Applied Clinical Medical Physics 2006;7(4):8198. 34. Rajon DA, Bova FJ, Bhasin RR, Friedman WA. Poster Presentation: Fabrication of custom surgical guides using rapidprototyping technology. University of Florida College of Medicine. Celebration of Research; March 20,2007. 35. Saw CB, Ayyangar K, Suntharalingram N: Coordinate transformation and calculation of the angular paramters for stereotactic system. Med Phy1987; 14:10421044,. 36. Schad L. Lott S, Schmitt F, Stum V, Lorenz WJ: Correction of spatial distortion in MR imaging: a prerequisite for accurate stereotaxy. Journal Computer Assisted Tomography 11(3):499505,1987. 37. Schileo E, Taddei F, Malandrino A, Cristofolini A, Viceconti M. Subjectspecific finite element models can accurately predict strain levels in long bones. Journal of Biomechanics. 2007 doi: 10.1016/j.biomech.2007.02.010 38. Schroeder W, Lorensen B. The visualization Toolkit. New Jersey: Prentice Hall, 1997. 39. SedefM, Samur E, Basdogan C. Realtime finiteelement simulation of linear viscoelastic tissue behavior based on experimental data. IEEE comput Graph Appl 2006;26(6):5868. 40. Smith KR, Frank KJ, Bucholz RD. The NeuroStationa highly accurate, minimally invasive solution to frameless stereotactic neurosurgery. Computerized Medical Imagining and Graphics 1994;18(4):247256. 41. Shewchuck JR. What is a good linear element? Eleventh International Meshing Round Table. New York: Sandia National Laboratories September 2002:115126. 42. Spiegel EA, Wycis HT. Stereoencephalotomy, Part I. New York: Grune & Shatton, 1952. 43. Spiegel EA, Wycis HT, Lee MM. A: Stereotactic apparatus for operations on the human brain. Sciencel947; 106:349350. 44. Spiegel EA, Wycis HT: Thalamotomy and related procedures. JAMA 1952;148:446451. 45. Suresh G, Elliot N, Pinson N, Sinden FW. Simulation of dynamics of interacting rigid bodies including friction I: General problem and contact model. Engineering and Computers 1994; 10:162174. 46. Timoshenka. S P, Gere J M. Mechanics of Materials. California: Thomson Brooks/Cole Engineering Division, 1984. 47. Van Cleynenbreugel J., S. F., Goffin J., Van Brussel K., Suetens P.. Imagebased planning and validation of C1C2 transarticular screw fixation using personalized drilled guides. Computer Aided Surgery. 2002;7:4148. 48. Van Staden RC, Guan H, Loo YC. Application of the finite element method in dental implant research. Comput Methods Biomech Biomed Engin. 2006;9(4):257270. 49. Yan X, Gu P. A review of rapid prototyping technologies and systems. Computer Aided Designl996;28(4):307318. 50. Yoo T. S., M.J., Chen D. T., Richardson A. C., Burgess J., Template Guided Intervention: Interactive Visualization and Design for Medical Fused Deposition Models. Proceedings of the Workshop on Interactive Medical Image Visualization and Image Analysis 2001;18:4548. 51. Zamorano LJ, NolteL, Kadi AM, Jiang Z. Interactive intraoperative localization using an infraredbased system. Stereotactic and Functional Neurosurgery 1994;63:8488 BIOGRAPHICAL SKETCH Barbara Garita Canet was born in 1977. She is the second child of Donato Garita and Noemi Canet who had six children. She grew up in San Jose, Costa Rica and attended Lincoln School from 1983 to 1995. After graduation from Lincoln School, she started her career in mechanical engineering at the Universidad de Costa Rica; however she only studied there for a semester. She was offered a scholarship for playing varsity volleyball for Florida Institute of Technology (FIT). From FIT she graduated with honors with a bachelor's degree in mechanical engineering in May 2000. In August 2000 she started a master's in biomedical engineering at the University of Florida and graduated in December 2002. The topic of her master's thesis was related to bone mechanics. From December 2002 to January 2003 she looked for a new committee chair to work on her doctorate. Finally, in January 2004 she stated her doctorate, on mechanical stability for customized instruments in imageguided surgery, and in August 2007 she completed her doctoral studies. Immediately after graduation she will rest and travel. She will then look for a job in biomedical engineering field and will write a business plan with the goal to start her own company. PAGE 1 1 GUIDED SURGERY USING RAPID PROTOT YPING PATIENTSPECIFIC GUIDES: A METHODOLOGY TO QUANTIF Y MECHANICAL STABILITY AND UNIQUENESS OF FIT By BARBARA GARITA A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2007 PAGE 2 2 2007 Barbara Garita PAGE 3 3 To Donato, Noemi, Esteban, Pablo, Ifi, Ignacio and Natalia. PAGE 4 4 ACKNOWLEDGMENTS First, I must thank my committee for taki ng the time to work with me during this endeavor. I must also express gratitude to my advisor Dr. Fra nk Bova for persevering with me throughout the time it took me to complete this research and write th is dissertation. Dr. Friedman, Dr. Gilland, Dr. Kim, and Dr. Vemuri have also been generous enough to dedicate their time and expertise to improve my work; I am grateful for their contributions and their kind support. I am grateful as well to Dr. Didier Rajon, for his time ly appraisals of my work and useful suggestions th roughout these years. I must express my profound gratitude to Dr Jonathan Earle for invaluable support and friendship, especially during the last year of my dissertation. Jeff Citty has been an immense support as well; working with him and the EFTP pr ogram for the last 4 year s has given me a true sense of belonging while at th e University of Florida. I need to express my gratitude and deep appr eciation to Serggio who always stood next to me in good times and bad times. I must acknowledge, as well the many friends and colleagues who assisted, advised and supported my research and writing effort s over the years, my gratitude goes to Ashley, Mafer, Ryanne, April, Luis, Mike, Pam, Atchar and Rolando. Finally I sincerely thank my parents and siblings for enc ouraging my decision to stay away from home for so many years in order to accomplish these academic dreams. I find the time spent separated from them irreplaceable and I hope I will get the chance to make it all up in the near future. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS...............................................................................................................4 LIST OF TABLES................................................................................................................. ..........7 LIST OF FIGURES................................................................................................................ .........8 ABSTRACT....................................................................................................................... ............12 CHAPTER 1 INTRODUCTION..................................................................................................................14 2 BACKGROUND AND SIGNIFICANCE..............................................................................19 ImageGuided Surgery........................................................................................................... .19 FrameBase and Frameless Stereotactic Approaches......................................................20 Model Based Guidance....................................................................................................22 Customized positioning frames for pedicle screw placement..................................22 Customized position frames fo r intracranial surgery...............................................24 Significance................................................................................................................... .........27 Facial Landmarks and Angles................................................................................................28 The Finite Element Method and ITK......................................................................................29 Elements Used in the Simulation....................................................................................30 The space truss.........................................................................................................30 The linear tetrahedron..............................................................................................34 Characteristics of Stiffness Matrices...............................................................................35 Linear Solvers................................................................................................................. .35 3 MATERIALS...................................................................................................................... ...47 4 METHODOLOGY FOR MEASURING STABILITY..........................................................54 Introduction................................................................................................................... ..........54 System Components and FEM Setup.....................................................................................55 External Forces................................................................................................................56 Rigid Layer.................................................................................................................... ..56 Bed of Springs.................................................................................................................58 Program Parameters............................................................................................................. ...59 Input Parameters..............................................................................................................59 Spring height............................................................................................................59 External force magnitude.........................................................................................59 Rigid layer height.....................................................................................................60 Difference in elasticity.............................................................................................60 Output Parameters...........................................................................................................60 PAGE 6 6 Criteria....................................................................................................................... .............62 TNT Ratio...................................................................................................................... ..63 Single Point Criterion......................................................................................................64 TNT threshold..........................................................................................................65 Consistency of TNT threshold.................................................................................67 Single Force Criterion.....................................................................................................67 5 PROGRAM DESCRIPTI ON AND VALIDATION..............................................................82 Program Description............................................................................................................ ...82 Validation of FEM elements...................................................................................................85 Linear Spring Validation.................................................................................................85 Tetrahedron Element Validation.....................................................................................86 A Simple Model Validation.............................................................................................86 6 RESULTS AND DISCUSSION.............................................................................................93 7 CONCLUSION AND FUTURE WORK.............................................................................109 APPENDIX A PRINCIPLE OF VIRTUAL WORK AND ITK FEM MODULE........................................111 FEM Terminology for Truss Elements.................................................................................111 Principle of Virtual Work.....................................................................................................112 ITK FEM and Principles of Virtual Work............................................................................113 Displacement Functions................................................................................................115 Shape Functions.............................................................................................................115 Strain Displacement Equation.......................................................................................116 Local Element Stiffness Matrix.....................................................................................117 Gauss Quadrature..........................................................................................................118 ITK FEM Module.................................................................................................................119 LIST OF REFERENCES.............................................................................................................123 BIOGRAPHICAL SKETCH.......................................................................................................127 PAGE 7 7 LIST OF TABLES Table page 41 Relationship between the tetrahedral height and a measure of the rigidity of the rigid layer.......................................................................................................................... ..........70 42 Table that shows the relationship between the difference in elastic modulus between the rigid layer and the bed of springs a nd the rigidity measure and Cholmods rcond value.......................................................................................................................... .........70 43 Single point criteria...................................................................................................... ......70 51 Displacement in the X, Y and Z direction from the example in Figure 54 tested with the ITKFEM module and ANSYS.....................................................................................87 61 Measurement of the facial angles for the three anatomical models...................................99 62 Instability score, presented as the percenta ge of forces that resulted in instability and the total forces applied on a frame.....................................................................................99 PAGE 8 8 LIST OF FIGURES Figure page 11 Components involved in this investigation and the prot otypes for the new technique for performing imageguided surgery with patientspecific reference guides...................18 21 Pretarget localization and presurgical planning is done in a computer workstation.......36 22 Procedures related to fram ebased image guidance surgery..............................................37 23 Frameless stereotactic surgery uses fiduc ials that are glued to the patients skin.............37 24 Typical layout when frameless stereotactic surgery is employed......................................38 25 Template designs for the spine..........................................................................................38 26 Router guides for craniotomy............................................................................................39 27 Biopsy guides.............................................................................................................. .......40 28 Rapid prototyping machines..............................................................................................41 29 Chart to visualize the significance of the present work.....................................................42 210 Facial landmarks.......................................................................................................... ......43 211 Facial angles............................................................................................................. ..........44 212 Space truss element has six displaceme nts to establish its deformed position..................45 213 Application of a force perpendicula r to the axis of the space truss...................................45 214 Space truss in local coordinates.........................................................................................46 215 The linear tetrahedron ha s a unique numbering scheme....................................................46 31 The three models used to evaluate the methodology for testing unique fit.......................49 32 Masks built for testing the methodology of fit...................................................................50 33 Caplike reference frame................................................................................................... .51 34 Virtually created half sphere test model............................................................................51 35 Placement of pins for a simple test on the real stability of the large and small masks......52 36 Coronal view of placement of the pins onto the real masks..............................................52 PAGE 9 9 37 The method employed to demonstrate the difference in stability between the large and small masks................................................................................................................ .53 41 The methodology for stability of fit i nvolves simulating a rigid layer onto of a deformable bed of springs..................................................................................................71 42 Displacement of a mask as it responds to external forces.................................................72 43 Triangulation of a customized pos itioning device proper to STL format..........................72 44 Direction of the Nmain.........................................................................................................73 45 Nodes point on the rigid frame..........................................................................................73 46 In a wedge three tetrahedral elements are fitted................................................................74 47 Schematic showing how the maximum rigidity fault is calculated...................................75 48 The behavior of a linear spring element under the application of an axial force...............75 49 Schematic that demonstrates how eac h springs displacement was evaluated by identifying the normal (ni) and tangential (ti) displacement magnitudes...........................76 410 Schematic that shows why it is important to distinguish between the springs that are mostly in compression or tension and those that are basically s lipping or have a high TNT ratio...................................................................................................................... .....77 411 Bounding volume for a spring to be in compression.........................................................77 412 Twodimensional interpretati on of the single point criteria..............................................78 413 Schematic to demonstrate the importance of evaluating the percentage slip of TNT ratio for each spring.......................................................................................................... .78 414 Schematic to shows at the level at which the TNT threshold was defined........................79 415 Empirical determination that the TNT th reshold greater than 0.6 correlates with the training data, which says that the large masks are stable...................................................79 416 Empirical determination that the TNT thre shold set smaller than 0.8 correlates with the training data, which says that the small masks are unstable........................................80 417 A case of the large Asian mask, whic h shows the consista ncy of an upper TNT threshold of 0.7............................................................................................................... ...80 418 Schematic showing the relation betw een the volume of the entire frame (fvol) and the compression patch volume (cvol)........................................................................................81 419 Box ratio interpretation.................................................................................................. ....81 PAGE 10 10 51 Flow chart of program for testing the uniqueness of fit.....................................................88 52 Original surface is extended in the positive normal direction to form the top surface and the points extended in the negative normal direction to become boundary condition points............................................................................................................... ...90 53 Examples of al linear spring element or space truss tested in all 3dimensions................90 54 Example of two tetrahedral elements us ed in the validation of the tetrahedral elements used in the program............................................................................................91 55 A corner frame that was used to validate the program......................................................92 61 Placement of Caucasian masks onto solid model............................................................100 62 Placement of Asian masks onto solid model...................................................................100 63 Placement of infant masks onto solid model...................................................................101 64 Placement of a caplike reference frame onto an upper head model...............................101 65 Surface plot color coded by the maximu m full rotation angle for the Caucasian masks.......................................................................................................................... ......102 66 Surface plot color coded by the maximum full rotation angle for the Asian masks........102 67 Surface plot color coded by the maxi mum tilt angle for the infant masks......................103 68 Surface plot color coded by maximum tilt an gle for the cap reference frame design.....103 69 Surface plots color coded by maximum tilt angle for the half sphere test.......................104 610 Normalized histogram of box ratios for each tilt angle ( ) for the Caucasian large mask........................................................................................................................... ......104 611 Normalized histogram of box ratios for each tilt angle ( ) for the Caucasian small mask........................................................................................................................... ......105 612 Normalized histogram of box ratios for each tilt angle ( ) for the Asian large mask.....105 613 Normalized histogram of box ratios for each tilt angle ( ) for the Asian small mask.....106 614 Normalized histogram of box ratios for each tilt angle ( ) for the infant large mask.....106 615 Normalized histogram of box ratios for each tilt angle ( ) for the infant small mask.....107 616 Normalized histogram of box ratios for each tilt angle ( ) for the cap reference frame design......................................................................................................................... ......107 PAGE 11 11 617 Normalized histogram of box ratios for each tilt angle ( ) for the half sphere test.........108 618 Profiles of solid models.................................................................................................. .108 A1 Differential element subjected to real stresses and in equilibrium..................................121 A2 The greater part of the ITK FEM library.........................................................................122 PAGE 12 12 Abstract of Dissertation Pres ented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy GUIDED SURGERY USING RAPID PROTOT YPING PATIENTSPECIFIC GUIDES: A METHODOLOGY TO QUANTIF Y MECHANICAL STABILITY AND UNIQUENESS OF FIT By Barbara Garita August 2007 Chair: Frank Bova Major: Biomedical Engineering Rapid prototyping technology (3dimensional pr inting) in conjunction with the virtual manipulation of computergenerated 3D anatom ical models rendered from diagnostic images (CT and MRI) is a novel techniqu e for directing imageguided surgery. By incorporating a graphical user interface, the surgeon can plan an intracranial surgical procedure by fabricating a patientspecific reference frame (masklike facial frame) through the us e of rapid prototyping technology. Early testing of these patientspecific frames revealed the necessity to include a step prior to fabrication to analyze the design of the frame for stability of fit. In response to this need, the objective in the present investig ation is to develop a method that numerically assesses the frames stability. The approach was to create a finiteelement simulation that models the patientspecific reference frame as a rigid layer com posed of linear tetrah edron elements under a deformable material simulated by a bed of linear spring elements. The program analyzes the displacement of the rigid layer by simulating extern al forces distributed around its perimeter. The stability of fit is characterized by the response of the linear springs. PAGE 13 13 Three patients with different facial anatomical features were selected as training data for testing the simulation program. For each case tw o patientspecific reference frames, one that involves more surface area and su rface variations such as co ntours of the nose bridge and another that consists of less surface area and less surface variations, were built with Rapid prototyping based on the patients 3D anatomical facial models. In all three cases the methodology correlated a stability sc ore with the stability of the real reference frames. In addition two more cases were evaluated; the firs t did not incorporate the facial structures and was design to hold onto the top of the head in a caplike manner, the second case was a virtually created half sphere. In both cases it was exp ected to find no a unique fit, and correspondingly the simulation found both cases to be unstable. PAGE 14 14 CHAPTER 1 INTRODUCTION Starting in 2003, researchers in the radiosurgery and biology lab (RSB) at the University of Floridas McKnight Brain Institute developed the capability to fabricate a patientspecific reference frame that has a potential of beco ming a new type of imageguided surgery (IGS) technique [33]. The fabrication of the patientspecific reference fr ame involves coding a series of algorithms that render a detailed 3dimensi onal lifelike anatomical model from computer tomography (CT) or magnetic resonance imaging (MRI) [6] datasets (Figure 11). Onto a patientspecific 3dimensional anatomical model a surface is selected (virtually painted) from which the patientspecific reference frames are cr eated. The ability to create, visualize, and manipulate 3dimensional virtual anatomical mo dels on a computer monitor heightened the enormous potential of utilizing medical image datase ts to create customized reference frames for medical diagnostics and treatment, specifically in the area of image guidance surgery [36]. The objective the patientspeci fic reference frame is to replace existing technology of framedbased (fixing with pins a frame to the skull) and/or frameless (using fiducials and extra tracking equipment in the operati on room) stereotactic systems [6]. Stereotactic techniques are used in neurological research and/ or surgery to direct a surgical instrument (needle or electrode) into a specific location in the brain [1]. The patientspecific reference frames (or customized guiding template) are designed to tightly fit to a specific patient s facial anatomy contour [6] and allow guidance for neurosurgical applications such as cranioto mies, biopsies, and deep brain stimulation [34]. Several clinical trials were initiated to test this new imageguided surgical approach. One of the first involved developing a customized gui ding template that fitted onto vertebral bodies and served as a guide for pedicle screw placemen t. As the templates for vertebral bodies were PAGE 15 15 tested (through a nonsignificant risk IRB approved clinical tria l) on actual patients in the operation room (OR), the research team quickly identified a problem that was defined within the context of biomedical engineering: the computer selected surface for the reference frame could not accurately predict the surface that the surgeon was able to present (because of the residual soft tissue), therefore no unique surface for the template could be identified. Eventually, the idea of using customized guiding templates for the spine was abandoned and the development of this approach for cranial procedures became th e primary focus of the design application. The idea of using customized guiding templates for st ereotactic surgery prov ed to have significant potential [33], but the problem of stability of fit was still apparent. The problem of the stability of fit for the reference frames became of great importance because of the need to add attachments to the reference frame for particular intracranial guidance surgeries (Figure 11D and 11E). As the distance between the referenc e frame and the attachments increases, the lever arm increases and small movements on the refere nce frame will results in exaggerated moments on the attachments tips (Figure 11D). For this reason, the accuracy of the patientspecific reference frames depends upon the identification of a broad and conforming patient surface. This surface must allow a reference frame mask (negativ e from the anatomical models surface) to be designed in a manner that prevents both rotations and translations. In other words, the mask must fit tightly onto the patients surface and not slide or wiggle. Hence, the initiative of the present investigation stems from the necessity to solve this stability problem by developing a numerical methodology to evaluate the patientspecific frames design by testing its stability of fit prior to fabrication with Rapid prototyping technology. The methodology developed for testing stability of fit relies on the creation of a finite element simulation program that models contact interactions between pa tientspecific frame and PAGE 16 16 the patients anatomical face model (Figure 11). The approach used is the idealization of a contact mechanics problem by mean s of the simulating the patientspecific reference frame as a rigid layer composed of linear tetrahedron elements under a defo rmable material simulated by a bed of linear springs (or truss) elements. The a pproach will simulate the movement of the rigid tetrahedral layer by the application of external fo rces onto the frames surface. As training data for the methodology, three solid lifelike anatomical face models of three patients with notably dissimilar facial profiles were created with ra pid prototyping and for each patients two frames one the is known to be stable a nd the other one known to be unsta ble were also created. The methodology for testing stability of the frames i nvolves two criteria, both based on the resulting simulated movement of the frames onto the patient s 3dimensional anatomical models. The first criterion is concern with identifying locations of excessive movement and the second criterion provides with a yes or no answer to the stability of the frame on a per force basis. It will be shown that these two criteria are able to predict the degree of in stability of the selected surfacemask construct; by identifying when the construct is able to resist transl ations and rotations to typical external force directions. The methodology to test uniqueness of fit presen ted in this investiga tion is an essential component to the new technique for imageguided surgery that involves rapid prototyping. The design of the patientspecific re ference frame has the imperative requirement to have a unique fit onto the patients anatomy, this is the majo r concern of this new imageguided surgery technology [33]. Therefore the work presented here, wh ich numerically evaluates the mechanical stability and uniqueness of fit of the patientsp ecific reference frame, will overcome this major problem and will allow a surgeon to use the new technology with more conf idence and accuracy. In this way, the new IGS technology is a step cl oser towards being dist ributed to other the PAGE 17 17 medical communities. In addition, the methodology has been written in the same language and uses the same open source libraries as the main graphical user interface, from which the patientspecific reference frames are designed, so it can be easily adapted as a module component. PAGE 18 18 Figure 11. Components i nvolved in this investigation and th e prototypes for the new technique for performing imageguided surgery with pa tientspecific reference guides. A) The graphical user interface that can create 3dimensional patient anatomical model from a series of 2dimensional medical image datasets, and the selection of a contour surface (in red) for the fabrication of a pati entspecific reference frames. B) Physical patients anatomical model created with Ra pid prototyping for evaluating the stability of fit. C) Patientspecific reference fram es that conformingly fits onto the patients model. D) Virtual model and reference frame with attachments to guide a biopsy probe. As the moment arm increases, small rotations and translation on the mask will have an exaggerated effect at the level of the attachment tip. E) Physical reference frame and attachments to guide a biopsy probe. A B C E D Moment arm PAGE 19 19 CHAPTER 2 BACKGROUND AND SIGNIFICANCE ImageGuided Surgery During the early 20th century, Sir Victor Horsle y recognized the need fo r targeting specific diseased regions in the brain without injuring neighboring healt hy structures. He hired an engineer to provide with a method to solve the problem; the engineer cam e up with the idea of using three imaginary orthogonal spatial plan es (horizontal, frontal and sagittal) [23] to specify any location in the brain. Later, during the mid 20th century, Spiegel and Wycis et al. [41] introduced the first stereotactic application. As it was originally envi sioned, the stereotactic approach relates to movement in 3dimensional space; therefore, stereotactic surgery uses a 3dimensional coordinate system to locate small ta rgets inside the brain. Sp iegel and Wycis et al. [44], began the work of locating specific regions in the brain us ing intrinsic reference points and in 1952 [43] they published the first st ereotactic map of the brain using the pineal body as the only reference point. However, this localiz ation technique proved unreliable because of anatomical differences between pa tients and as the surgeons move d further from their reference point they lost accuracy. Seeking to improve th e reliability of their technique, Spiegel and Wycis et al. developed the approach of fitting a custom cap to the patient and attaching a head ring to this cap; in this way they established a relation between the brai n atlas and a reference head ring [43]. Besides the Spiegel and Wycis custom cap a pparatus, during the following year several others investigators developed stereotactic devices for differe nt purposes and with different usage methods [18]. Some systems were based on Cartes ian coordinates, ot hers used polar coordinates for placing electrodes at target locations; while some designs moved the patients PAGE 20 20 head towards a target location; ot hers used interlocking arcs to move the apparatus to get to a specific location [18]. With the advent of computer tomography (C T), imaging of the brain became readily available. The BrownRobertsWells (BRW) system was the first to be de signed particularly for use with CT images and was the first to impl ement a computer base algorithm to calculate coordinates and a trajecto ry to a target point [22]. As a result, a specific patients anatomical information obtained through CT imaging was, for the first time, properly visualized and referenced to a head ring [36]. Later, in 1987, Schad et al. [36] introduced multimodality image registration, which increased the accuracy and tissue differentiation for stereotactic targeting. FrameBase and Frameless Stereotactic Approaches Today image guidance procedures can be divi ded in two categories, according to different referencing methods employed. These are ge nerally termed framebase and frameless stereotactic approaches. The pur pose of both is to provide the surgeon with an understanding of the anatomy beneath the visible ope rative surface. State of the art stereotactic procedures allow the surgeon to physically relate spatial informati on from a set of diagnostic images (represented by a patients 3dimensional virt ual model), onto a patients actu al, realworld anatomy in a reliable and accurate manner. In general, th is is accomplished by means of a discernible reference frame (or landmarks) that are integrat ed in both the virtual patients model and the realworld patient and correlates the informa tion amongst them. The creation of a diagnostic plan, for both framebase and frameless stereota ctic surgeries, is pe rformed on the virtual, computer generated patient, where image data sets (MR and CT) depict the anatomical information together with the rigid frame land marks for framebased or skin landmarks for frameless stereotactic surgeries (Figure 21) [6]. PAGE 21 21 Framebase stereotactic surgery is th e most accurate technique available [32], but it is not free of significant disadvantages. This represents the most accurate method because it rigidly fixes a reference frame onto the patients skull, through the use of minimally invasive pins. For example, in the case of a stereotactic frameba sed biopsy a head ring is applied to the patient prior to imaging (Figure 22A). The diagnostic scan will allow the surgeon to visualize the intracranial anatomy and create a virtual surgical plan on a co mputer workstation. With a graphical user interface the surgeon will define th e intracranial target location and the computer system will mathematically calculates the corres ponding coordinates relating the virtual plan to the real world patient anatomy (Figure 21). Th e coordinates are set on th e physical stereotactic system (Figure 22B) (which is directly fixe d onto the rigid frame), providing a guide or trajectory for the surgeon [8, 35]. The main shortcomings of framebase stereota ctic surgery include th e extend of time that the procedure takes (from to presurgical frame application to surgery completion), the specific expertise required for frame application, the hi gh cost of replacement equipment, and the discomfort that the head ring causes on the patient. Indeed, the rigid frame is a discomfort for the patient at the time of imaging and during the preparat ion of the sterile field, and it gets in the way of the surgical procedures in general. Nevert heless, framebase imageguided surgery provides the gold standard for precision and accuracy for intracranial stereotactic procedures. Frameless stereotactic surgery has been pr esented as an alternative to framebased procedures. The main difference between frameba sed and frameless stereotactic systems, is that instead of a head ring, fiducials (landmarks) are used (F igure 23). For cranial surgery these are either MR and CT identifiable markers that are temporally glued onto the patients skin [2, 51] or anatomical landmarks are defined on 3dimensi onal rendered images base upon the patients CT PAGE 22 22 or MR datasets [39]. In both cases, the procedures do not have to be completed in an extended and continuous span of time. Both framebased and frameless procedures permit the same virtual computer planning to be performed by the surgeon; however frameless stereotactic surgery requires instrumenttracking equi pment within the operating room (OR) [51]. The main drawback of the frameless approach is the necessity for registration of the fiducials with a reference frame used to track movement between the patient and the surgical tools [6]. The registration requires special cameras (electromagnetic) and additional computer workstations to be present in an already crowded operating room (F igure 24). Although electromagnetic tracking allows fo r nonline of sight track ing, the surgical t eam is still burdened with the additional computers in the OR [6], special traceable instrume nts, the need to minimize ferrous metals from with the operative field and the registration procedures. Model Based Guidance Customized positioning frames for pedicle screw placement As mentioned in the introducti on, researchers in the RSB lab initially built customized guiding templates to fit onto vertebral bodies for th e placement of pedicle screws. It was during the testing of these spine templates that the probl em of uniqueness of fit was first addressed. The following section provides background informati on on the first implementations of custom guiding templates for the spine and also provides additional information to aid in better understanding of the problem at hand. Stereotactic guidance with the aid of custom ized guiding templates initiated in 1999 for the placement of pedicle screws [20]. In 1999, Goffin et al. [20] first published a paper that explained about a guide for the posterior transart icular screw fixation of C1C2 (cervical 1 to cervical 2); and in 2001, Goffin et al. [21], Birnbaum et al. [6] and Yoo et al. [50] published papers referring again to template guides for place ment of pedicle screws in the spine (Figure 25 PAGE 23 23 AC). While Goffin et al. [20] developed templates for the fixation of the cervical spine, Birnbaum et al. [6] and Yoo et al. [50] tested particularly on the lumbar spine of cadavers and humans. All three research groups us ed rendered 3dimesional CT information to create the template guides based on the bony surface of the spine. Goffin et al. [20] printed the 3D solid models of the cervical spine with stereolithography and physically defined the correct placement of the pedicle screws, and then he translated the trajectory to a CAD (Computer Aided Design) program. Goffin et al. [21] and Birnbaum et al. [6] used computer controlled milling machines or NC (numerically cont rolled) machines to create the temp late guides; Yoo et al. used FDM (Fused Deposition modeling), a machine from Strata sys (Stratasys, Tempe, AZ), that follows the same principle of ster eolithography but is composed of a m ovable nozzles where melted plastic forms the template one layer at a time. Goffin et al. [21] manufactured the device on an acryl ate resin and the template had a horseshoe appearance that was positioned on the posterior aspect of C2 (Figure 25 A). He tested on some cadavers and two patients. On the cadavers he initially experienced rotational instability of the template and increased the co ntact area to include more vertebral support. On one of the two patient cases he was successful; while on the other the template guide did not perform well. Birnbaum et al. [6] created polycarbonate template s with NC machines (Figure 25 B), and conducted a cadaveric study compari ng computerassisted surgery with individual templates on thirteen lumbar spines. Under imageguided technique, he found two misplaced screws; while all were correctly placed when us ing the template guides. Yoos templates were constructed from ABS plastic and expanded over a large contact area comprising both transverse PAGE 24 24 processes, the laminae areas and the spinous process [50]. He tested on dry dissected vertebrae. In all cases the positions of the screws were evaluated with postope rative CT scans. A later paper by Van Cleynenbreugel et al. [47], published in 2002, focuses on the validation of the template guides for the C1C 2 transarticular screws (continuing on Goffins work). As the others had previously done, he us ed a preoperative CT to base the design of the template and the testing was done on cadavers. During the first part of the study, he ran into some serious instability issues, particularly in the direction of the application of the drilling force; the result was dramatic displacement of th e template guide on the order of 2 to 9 mm. Consequently, they altered the design of the te mplate guides to include more support and had significantly bette r results. The appearance of the templates developed by th e three researchers resembled a block that overlaps a broad surface area over the spine. Fr om the experience gained in the RSB Lab, it was found that the templates that overlap a broad area on the spine create the necessity of removing a large amount of muscle tissue to expose the ma ting bony surface of the spin e; this design is not beneficial for the patient and time consuming for the surgeon; in addition residual soft tissue causes the surface conformity between the templa te and the anatomy to be disrupted. As a consequence, the work on custom guiding te mplates for the spine was abandoned; however important knowledge was learned from the experience Figure 15 shows the template designs published by other investigators an d the initial designs of the spine template from the RSB lab. Customized position frames for intracranial surgery The first to suggest customized position fr ames for the head was DUrso et al. [16]. His approach involved implanting reference sockets on a patients skull prior to scanning. Then from the diagnostic images he create a physical 3dimensional m odel of the patient head that included the mounting sockets. The surgical trajectory was planned on the solid model by PAGE 25 25 attaching a reference arc to the mounting sockets. By correspondence, the arc was refitted onto the patient and directed the surgeon to the predetermined location in the brain. The technology developed in the RSB for performing imageguided surgery by use of patientspecific frames has been shown to have superior accuracy and precision when compared to previously documented framebased procedures [33]. Rajon et al. [33] used a simple phantom to measure the point localization accura cy of the customized guiding frames. He found that the technology could position a probe tip with an accuracy of 1.7mm. As mentioned before, the standard in terms of accur acy for imageguided surgery is the framebased approach. The overall system accuracy of a standard frame ba sed procedure was measured by Maciunas el at. [27] and determined to be 2.28mm. Their accuracy measurement is an estimated measure of the imaging accuracy and the mechanical system accur acy extrapolated for th e ideal case of zero slice thickness (CT or MR slice thickness) to remove image localization related errors. Today in the RSB lab the researchers have de veloped customized positioning frames with specific component attachments to perform cran iotomies and biopsies (Figure 26 and Figure 27). Both of these procedures are incorporated into a graphical user in terface, with the aim of presenting the technology as anothe r tool available to the neurosur geon; one where he or she can preplan a surgical procedure and fabricate his or her own patientspecific frames and surgical guiding attachments. Rapid prototyping and RPD Designer. Rapid prototyping (RP) can be accomplished using a number of approaches. One such approach is 3dimensional printing RP (Figure 28A and B), another is subtracting rapid prototyping (SRP) (Figure 28D). Both RP technologies provide a means to fabricate a physical 3dimensi onal solid model from in formation contained in a stereolithography (STL) file. The file consists of planar triangular facets that represent a 3 PAGE 26 26 dimensional surface [29]. In the case of 3dimensional RP the first step in creating the solid model is the virtual slicing of the STL file by a preprocessing program [49]. For each virtual slice a layer is printed in real space. The bottom layer is built fi rst, followed by additional layers until the last layer is built (Figure 28B). Th e RP system purchased in the RSB lab is built by ZCorp model 310 (Burlington, MA) and has a build accuracy of 0.3 mm [33]. The system creates the 3dimensional model by translating a thin layer of powder from a reservoir location onto the building chamber. Subsequently, a bind er fluid fuses the powder as specified by the slicing of the STL file. Afterw ards, the build chamber moves ve rtically downwards allowing for a new layer of powder to be placed onto the prev ious layer, and as the binder fluid is placed again it will also glue within la yers. When all layers are built th e solid model is taken out of the building chamber (Figure 28C) and infiltrated with cyanoacrylate glue to assure durability [33]. The 3dimensional solid, or in the case of this investig ation the patientspecific frame, can be gas sterilized to be used within a sterile operative field [33]. The RSB lab purchase another machine to pe rform rapid prototyping, this machine is a programmable milling machine, and instead of performing 3dimensional printing, it starts with a solid object and removes unwanted material (Figur e 28D). Because it removes material, this type of rapid prototyping is called subtractive rapid prototypi ng (SRP). There are several reasons why the researchers i nvested in the model MDX650 (Roland DG Corp., Knoxville TN); the main reasons are its speed and the capability of using a variety of materials to fabricate the reference guides. With the subtractive rapid pr ototyping machine the fabrication time was been cut in a third of the original time, from 11 hours to about 3 or 4 hours [34]. Also the reference frames created with this machin e can be sterilized at high a temperature which is much faster than cold sterilization required by the patientspecific reference frames created with the 3 PAGE 27 27 dimensional RP machine. In addition, early accura cy tests demonstrate that the precision of the SRP machine is superior to the 3D printing machine [34]. The creation of the STL file is conducted by a graphical user interf ace created in the RSB lab, called RPD designer [33]. The RPD designer has capabili ties of rendering a solid model from a set of diagnostic images. The patientsp ecific frame is created by selecting a reference surface (called selected surface) by painting th e surface of the virtual 3dimensional rendered anatomical model [33] (Figure 11). Significance There is a great potential for the new proposed idea of using diagnostic images and rapid prototyping to create a custom ized guiding frame for use in imageguided surgery. The most significant benefits of this new technology are rooted on the fact that it will eliminate the need for fixed frame application thereby reducing pati ent discomfort. It also has the ability to significantly reduce the cost of re placement parts related to framebase stereotactic surgery. In addition, the use of the patientspecific referen ce frames will eliminate registration processes required by frameless stereotactic surgery and allow for reorganization of an overwhelmingly crowded OR. The new technology will provide the su rgeon with the ability to plan and create his or her own surgical guides while reducing patient wait times during surgical planning. In addition, the surgeon will be able to remove a nd replace the patientspecific reference frame onto the patient as required, and maintaining high accur acy on a range of cranial procedures such as shunt placement, skin incision placement, bone removal during a craniotomy, biopsies, tumor resections and deep brain stimulation [34]. Besides making stereota ctic neurosurgery more minimally invasive, the patientspecific refere nce frame could maintain framebased accuracy [33]. Furthermore, the cost for creating a frame and a whole set of attach ments to perform model PAGE 28 28 based imageguided surgery is estimated to be approximately $50.00 per patient in materials cost [6]. Rajon et al. [34] has already created the software that provides the surgeon with the ability to directly build a virtual 3dimensional model from diagnostic image datasets. This program allows selecting a surface on th e virtual 3dimensional model for the customized guiding frame and appropriately includes components to effect uate his or her approach towards a patientspecific intracranial surgery. The significance of the present i nvestigation is based of the fact that some level of testing is required prior to the fabrication of the custom ized guiding frame. The results of this testing will provide the clinician with prior information of frameonpatie nt stability. In addition, it will allow the clinician to effectuate his surgic al plan with more confidence and accuracy. Facial Landmarks and Angles The stability of fit of the patientspecific fr ames (or masks) is greatly benefited by the surface characteristics of the human face. If a reference frame is to be designed based on the superficial surface of the human head, the cont ours of the human face serve the purpose of anchoring landmarks. The faci al characteristics are ideal because of various neighboring irregular and curved surfaces. Besides, the facial contour is in close proximity to the brain, the skin is smooth and in most patients, close to th e bony structure, to which the brain is naturally referenced. Therefore, it is important to obt ain a basic understanding of the surface anatomy of the face. Figure 210 illustrates the facial anatomical landmarks of importance. On the forehead one can palpate the frontal eminences, which are prominences of th e frontal bone that vary among individuals and are often unsymmetrical [10]. Caudal and medial to the frontal eminences is the glabella, a smooth somewhat triangular bony promin ence on the frontal lobe located above the nasion which is the point of in tersection between the two nasal bones and the frontal bone. PAGE 29 29 The nasion represents the root of the nose. Lateral to the glabella bilaterally, one finds the superciliary arches : two arched elevations of the frontal bone that tend to be more prominent on the medial aspect of the face. Some males have notably prominent superciliary arches that may significantly contribute to the stab ility of the mask. Superciliary aches tend to be small in females and absent in children. Caudal and lateral to the orbit, bilaterally, is the zygomatic bone, which forms the prominence, known commonly as the cheek bone. The prominence of the zygomatic bone can vary significan tly among individuals. Finally, the orbital area is also an important area of depression that cont ains the eye within the bony skull [10]. Furthermore, based on these facial surface anatomy landmarks, one can ex trapolate facial angles that are important to our discussion, as they also play an important role in determining the stability of the mask (Fig 211). The Finite Element Method and ITK The patientspecific frames are complex in geometry and shape, and the best way to predict its behavior as it response to small external forces is by applying the principles of finite element analysis (FEA). The FEM involves dividi ng (discretizing) the entire frame into smaller parts or finite elements [5]. By dividing the entire frame into smaller elements, it is possible to obtain an accurate prediction of the disp lacement of the frame by approximating the displacement of each element to a polynomial function [28]. The patientspecific frames are modeled as ri gid frames composed of linear tetrahedral elements under a deformable interface composed of truss elements (or linear springs) (Figure 52). The deformable layer, composed of truss elemen ts, is attached at one end to the rigid layer, and on the other end all the truss elements are fi xed (cannot moved in any direction x, y and z); this fixing of the elements on one end represent patients surface. PAGE 30 30 In addition, FEM has been used extensively in medically related mechanical analyses and modeling [12], such as in orthopedics [26], dentistry [48], head injury [25] and organ analysis[39]; and to evaluate stresses in human bones [37], to mention a few. For the present investigation, a linear 3dimensional truss el ement was added to the local ITK FEM (ITK 1.8, Kitware Inc., Clifton Park, NY) library because it was not available. The ITK FEM module was created for image registration and has three iterative solvers to solve large systems of linear equations; these solvers are slow and not adapted for the cases in the present investigation. A new lin ear solver, Cholmod [13] was adapted which can solve sparse systems of equations in seconds compared to hours that th e available ITKs linear solvers take. Also, the ITK library was carefully studied and manipulated to make the methodology defined in this work progress as fast as possible (Appendix A). The IT K FEM module was chosen for this assignment to provide compatibility to the already bu ilt and working graphical user interface (RPD designer). Elements Used in the Simulation Two simple elements were used for this in vestigation. The first one is a 2node truss element in 3dimensional space called space truss. The second one is the linear 4node solid tetrahedral element in 2dimensional space, which was available in the ITK FEM module. The space truss Space trusses are frequently us ed to model lattice domes and aerospace structures [4]. They compose 3dimensional structures that ar e subjected to 3dimens ional force systems. A space truss element is a long, thin and straight pris matic element. Prismatic means that its cross sections are all homogeneous by having a cross sec tional area (A) and an elastic modulus (E) that remains constant through out the length of the truss. Space trusses are connected by nodes that act as frictionless ballandsocket joints [17]. In a structure composed of space trusses, any PAGE 31 31 unsupported node (with no assigned boundary conditions) wi ll result in translations in three local directions (x, y, and z). If no boundary conditions are applied to its nodes, then each truss element has six node displacements or six degree s of freedom (DOF) to establish its deformed position [4] (Figure 212). A space truss element is defined by the Youngs modulus (E), cross sectional area (A), the length (L), and a Poissons ratio ( ). One important characteristic of the any truss element is that it has a preferred dimension, that is, along its longitudinal or axial direction (Figure 212). By definition, a truss element will only develop axial forces [17]. Any perpendicular forces will not change its original length (L), since the truss element does not undergo moment s and its nodes act lik e frictionless ballandsocket joints. For example, Figur e 213 shows a perpendicular force (fy) applied to a truss element; as a result of the force there is no di fference between the orig inal length (L) and the length succeeding the applic ation of the force (L). Hence, if there are not differences in lengths (LL = 0), there were no internal forces and no displacements either (Figure 213). The space truss is considered a simple element since it deforms only along one direction (along its local xaxis) [5]. For simple elements like trusses, the direct stiffness method can be use to determine the elements stiffness matrix ( ke). The direct method uses principles of equilibrium and mechanics of materials [46] to define the element stiffness matrix, which is a relation between node displacements (u1 and u2) and internal forces (f1 and f2) for each element e (Figure 212). The characteristic of the space truss element of only responding to axial forces makes the calculation of the element stiffness matrix relative ly easy. When the element is analyzed in its local coordinate system, by definition there are no internal forces perpendicular to the truss and hence no perpendicular displacements. In othe r words, the displacement vector for a truss PAGE 32 32 element is given by Equation 21, an d the local node displacements u1 and u2 are caused by the local node forces f1 and f2, as seen in Figure 214. 2 1u u u (21) From mechanics of materials [46], the relationship between the local node displacements and the local node forces for a prismatic tru ss element is linear and given by Hookes Law [46], or Equation 22, where k is the stiffness coefficient and is equal to AE/L. ku f (22) With Equation 22 and the principles of equilibrium [28], the element stiffness matrix ( ke) for a particular element e can be calculated. As Figure 214A shows, to establish the element stiffness matrix with the direct stiffness method [28], one needs to confin e the element to a unit displacement at node 1 (n1) and hold node 2 (n2) fixed. The force at node 1 is equal to k*1 (by Equation 22). Since the element is in equilibrium, then the internal force at node 2 must been equal to k*1. A similar situation is true (F igure 214B) when node 2 is confined to a unit displacement, and the internal forces at node 2 and node 1 are equal to k and k, respectively. The element stiffness matrix is a 2 x 2 matrix defined in local coor dinates by Equation 22. k k k kek or 1 1 1 1 L AEek (23) For the stiffness coefficients kij the first subscript identifies the direction of the force, while the second subscript is related to the displacement. The element stiffness matrix represents the forces at the ends of the element as functions of the displacements at the nodes (Equation 24). 2 1 2 11 1 1 1 u u L AE f f or 2 1 2 1u u f fek (24) PAGE 33 33 The element stiffness matrix (Equation 24) de fines the space truss in local coordinates. The transformation to global coordinates is necessar y to represent the truss element as part of a entire structure. The transformation ma trix and is represented by Equation 25. z y x z y x cos cos cos 0 0 0 0 0 0 cos cos cos T (25) Finally, for one element the local element stiffness matrix ke can be transformed into the global elements stiffness matrix Ke by Equation 26 (T means the transform of the matrix). T k T KTe e (26) For the space truss element, the elements stiffness matrix in global coordinates Ke results in a 6 x 6 matrix (Equation 27). z z y y z x y x x z z y z x z z y y y x z y y z x y x x z x y x x eSymm L AE 2 2 2 2 2 2 2 2 2cos cos cos cos cos cos cos cos cos cos cos cos cos cos cos cos cos cos cos cos cos cos cos cos cos cos cos cos cos cos cos cos cos K (27) Following the direct stiffness method, the most impor tant steps to define the space truss element systematically is first to define the elements s tiffness matrix in local c oordinates (Equation 23); then to calculate the transformation matrix T (Equation 25) necessary to swap from local to global coordinates; and finally it is important to define the element stiffness matrix in global coordinates Ke (Equation 27). In the FEM, all elements stiffness matrices in global coordinates ( Ke for e = 1 to n, n is the number of elements) are assembled into a global stiffness matrix Kg that represents the forcedisplacement relations for an en tire structure composed, in th is case, of space trusses. Fg is a vector that that defines the external forces in global coordinates and Ug is the resulting node PAGE 34 34 displacements in global coordi nates (Equation 28). The solu tion to the system of linear equation is given by solving for the di splacements in global coordinates, Ug. g g gU K F (28) The linear tetrahedron Solid elements are defined in 2dimensions and can model a solid body. The linear tetrahedron is the simplest solid element in 2dime nsional space; four coordi nate points define its geometry with respect to a global coordinate system. It is known to be poor for stress analysis, as results tend to be more rigid than expected [17]. The numbering of th e nodes is important for generating every individual tetrahedron correctly. The numbering of the tetrahedrons nodes mu st be done following the righthand rule (Figure 214) and the calculation of the tetrahed rons volume must result in a positive value (Equation 29). First an initial corner it picked (node 0) and th en the triangular face opposite to the initial corner must be numbered followi ng a counterclockwise manner (nodes 1, 2 and 3). The calculation of the tetrahedrons volume give s evidence that the numbering of the nodes is done correctly. In Equation 29, the node coor dinates for the tetrahedral are given by xi, yi, and zi where i is the node number. 0 1 1 1 1 det 6 13 2 1 0 3 2 1 0 3 2 1 0 z z z z y y y y x x x x Volume (29) The derivation of the tetrahedral elements stiffness matrix in complicated, in fact it rarely derived using the direct stiffness method. In turn, the principles of virtual work are often used to derive stiffness matrices for solid elements, multidimensional elements [4]. PAGE 35 35 Characteristics of Stiffness Matrices The characteristics of stiffness serve to vali date the setup of the ITK FEM module for the two types of elements used in this investigation. Th e characteristic of the element stiffness matrices for the space truss, the linear tetrah edral element and the global stiffness matrix ( Kg) were analyzed for validation of the results. The element stiffness matrix can be singular; th is means that the determinant is equal to zero. Also they need to be squared and symmetr ic and all element along the diagonal need to be greater or equal to zero [28]. In the case of the reduced fo rm of the global stiffness matrix, it must be square, symmetric and positive definite and all its diagonal values have to be greater than zero [28]. In addition is must have and inverse so the displacement can be solved uniquely for a given force [46]. Positive definite matrices arise frequently in FEM problems. They are characterized by the condition that all their ei genvalues are positive [13]. One of the propertie s of a positive definite matrix is that it can be decomposed by Chol esky decomposition. If A is a symmetric positive definite matrix, then A can be f actored into a product of L and LT (T means the transpose of a matrix), where L is a lower triangular matrix w ith positive diagonal elements. The definition of a positive definite matrix is given by Equa tion 210, where X is any vector in R space [13]. 0 AX XT (210) Linear Solvers The ITK FEM library is fitted to solve image re gistration problems that incorporate finite elements. The ITK FEM module offers iterative solvers mainly used to solve very large matrices, common to image registration problems. The most reliable it erative solver in the module is the ItPack Solver [24], however due to its iterative nature ItPack is a very slow solver. An alternative to using ItPack is a solver based on Cholesky decomposition called Cholmod [13]. PAGE 36 36 Cholmod is a set of routines for factoring sparse symmetric positive definite matrices and solving linear equations used in Matlab (The Mathworks Inc. Natick, MA). Cholmod proved to be much faster and highly reliable; as a result it was adapted to be used with the ITK FEM module. Figure 21. Pretarget localizat ion and presurgical planning is done in a computer workstation. PAGE 37 37 Figure 22. Procedures related to framebased image guidance surgery. A) A patient with a head ring prior to scanning on a CT. B) A surge on setting the coordina tes on a stereotactic apparatus, which will be mounted onto a h ead ring (reference frame) to located the coordinates for targeting a sp ecific location on the brain. Figure 23. Frameless stereotactic surgery uses fi ducials that are glued to the patients skin. The fiducials serve for registering the patients anatomy in real space to the image datasets on a computer workstation. The tools for the surgery have similar fiducials so they can be tracked by cameras. A B PAGE 38 38 Figure 24. Typical layout when frameless stereotactic surgery is employed. Both images show how crowded the OR can get by the necessity of additional computer workstations and tracking devices. Figure 25. Template designs for the spine. A) Goffin [20], B)Birnbaum [6], C) and Yoo [50]. D) The RSB template designs from the firs t stages on the development of customized positioning frames. A B C D PAGE 39 39 Figure 26. Router guides for craniotomy. A) Se lection of the customized positioning frame onto the patients virtual 3D model. B) Sele ction of craniotomy contour drawn on the surface of the skull. C) Virtual model of the customized positioning frame. D) Craniotomy attachment for the customized positioning frame. E) Fabricated structure that will guide the surgeon to perf orm the craniotomy as planned [34] A B C D E PAGE 40 40 Figure 27. Biopsy guides. A a nd B) Graphical user interface shows the three orthogonal views from diagnostic images and the solid virtual model. The surgeon can select the target location and the graphical user interface will appropriately position the biopsy guide C) Virtual model of the customized positioning frame and the biopsy guide attachment. D) Fabricated structure that will guide the surgeon to perform the biopsy as planned [34]. A B C D PAGE 41 41 Figure 28. Rapid prototyping machines. A) Rapid prototypi ng machine by ZCorp model 310. B) A thin layer of powder is moved to th e building chamber by the gantry, which also holds the print head that distributes the bi nder fluid. C) The 3D model as removed from the building chamber. D) Subtractive RP machine. A B C D PAGE 42 42 Figure 29. Chart to visualize the significance of the present work. GUI: Neurosurgeons virtually define a target location Step 1: GUI renders model and selects surface for frame Step 2: A) Virtual evaluation of frame prior to fabrication. B) Is the patientspecific frame design STABLE? Step3: If YES, then build patientspecific reference PAGE 43 43 Figure 210. Facial landmarks. MedialLateral Caudal Cephalic Cheek (zygomatic bone) Orbital area Nasion Glabella Radix Superciliary arch Sagittal Plane Coronal Plane Frontal eminences PAGE 44 44 Figure 211. Facial angles [10]. The nasiolateral a ngle is calculated at th e level of the radix. Nasiolateral angle Nasiofrontal angle Nasiofacial angle PAGE 45 45 Figure 212. Space truss element has six displaceme nts to establish its deformed position. There are u1 and u2 along its xaxis, v1 and v2 al ong its yaxis and w1 and w2 along its zaxis. Though the space truss is defined in 3dimensional space, it will only develop axial forces, along its local xaxis. Figure 213. Application of a fo rce perpendicular to the axis of the space truss. There is no change in length (LL = 0) and there are no internal for ces (in a real simulation the truss element will not display as in the figure). The boundary condition symbol implies, in this case, that the element is fixed in local x and z directions. n2 u1 n1 X Z Y f1 f2 x y z Element e v1 w1 u2 v2 w2 n2 n1 X Z Y fy x y z L L Boundary condition symbol PAGE 46 46 Figure 214. Space truss in local coordinates. A truss element in its original position is denoted by the thick black line, and the same truss element in a deformed state is shown in gray. After a force f1 is applied to node 1 (n1) to create a unit displacement the stiffness coefficient k11 defined the element at n2. The reaction force at n2 is k12. Similarly, a force corresponding to a unit displacement is applied at n2 and k21 and k22 are defined. Figure 215. The linear tetrahedr on has a unique numbering scheme. n2 u1=1 n1 x z y L f1 n2 u2=1 n1 x z y L f2 A B 11 k f 12 k f face 1 2 3 0 1 2 3 Nodes Z Y X PAGE 47 47 CHAPTER 3 MATERIALS Three patients with notably dissimilar face prof iles were selected as training data for the methodology designed to test uni queness of fit. With CT data sets and the graphical user interface, RPD designer, three anatomical solid models (Figure 31) were built and name Caucasian, Asian and infant because they belonged to a Caucasian male, an Asian male and a infant respectively. Also six masks were built two for each model (Figure 32). For each solid model, a conforming small and a large mark were built; the small masks included considerably less facial landmarks than the large masks. While the large masks included the entire orbital area and extended laterally and media lly along the supercilia ry arches, the small mask only included the upper orbital region and extende d to the level of the glabe lla along the nose region. All masks were created from the pai nting of the corresponding select ed surface as shown in Figure 51. The six masks designed with RPD designer are the training data for the methodology to test stability and uniqueness of fit. In reality when the patient specific reference frames are applied onto a patient, the clinician will apply a force to keep the mask in their uniqu e location. Sliding, rotation and translation of the mask on the level of contact with the patient s face signifies exaggerated displacements at the level of the attachment tips, as shown in Fi gure 11D. While the large mask designed for training the methodology maintain their location when exposed to real transverse forces, the small mask would not withstand normal forces at the their extremities, along the length of the superciliary arches. To provide with a real world measure of stability for the small and large masks, a center pin was place normal to the surfa ce (or the direction on which the mask is placed onto the solid model) and at the location of the glabella on all ma sks. Oriented by the center pin, crosshairs were drawn on all masks along the transv erse and sagittal plan es as shown in Figure PAGE 48 48 3.5. Additional pins were located close to ends of the crosshairs lines as shown in Figure 35 and numbered from 1 through 4. Pins 1 and 3 were placed on the sagittal plane and pins 2 and 4 were placed on the transverse plane (Figure 36) A normal force of 5 N, applied by weight, was applied onto all pins for large and small masks (Figure 37). The large masks remain in their unique location, and sustained the 5N force; on the other hand, the small mask exhibit large displacement when the force was applied on pins 2 and 4. This resulting slippage behavior is mainly attributed to the sma ller surface are and less facial i rregularities covered by the small masks. All small masks withstood th e 5N normal force on pins 1 and 3. In addition, the placement of all large masks is unique, this is not the case for the small masks; there is not clear unique location where the masks fit. The sliding of all small masks when experiencing a normal force on pins 2 and 4 and the lack of a clear unique fit are evidence that the small masks are unstable, on the contrary all large masks withstood the 5N force on all pins and have a unique fit, which characterize them as stable (the fitting of the masks onto the models can be seen in Chapter 6, Figures 61 through 63). In addition to the masks, two other surfaces were created to evaluate the methodology of fit. A caplike reference frame (Figure 33) was created from CT im age dataset and a halfsphere surface (Figure 34) was virtually create d with Paraview (Parav iew 1.6, Kitware Inc., Clifton Park, NY). The cap frame covers the top of the head and ex tends orthogonally along the coronal and sagittal planes no furt her than at the level of the frontal eminences. The half sphere surface is a faceted surfa ce defined by 31 points evenly distribut ed (Figure 34). The half sphere needed to be defined as a faceted surface with ten angular rotation about it axis and three coplanar levels, because a smoother surface could can not be evaluated, since it results in a mechanically underconstrained system. PAGE 49 49 Figure 31. The three models used to evalua te the methodology for test ing unique fit. A) Caucasian, B) Asian and C) infant. A B C PAGE 50 50 Figure 32. Masks built for testi ng the methodology of fit. A) Caucasian small and large masks. B) Asian small and large masks. C) Infant small and large masks. A B C PAGE 51 51 Figure 33. Caplike reference frame. A) Upper head model and B) caplike frame. Figure 34. Virtually created ha lf sphere test model. A) Top view and B) side view. A B Sagittal plane Coronal plane Upper most cephalic region A B PAGE 52 52 Figure 35. Placement of pins for a simple test on the real stability of the large and small masks. A) Large Asian mask and B) small Asian mask. Figure 36. Coronal view of placem ent of the pins onto the real ma sks. Notice the orientation of pins 2 and 4 and difference in lateral exte nd between the A) large Asian mask and B) the small Asian mask. Transverse p laneSagittal plane 1 4 2 3 A B 4 2 4 2 A B PAGE 53 53 Figure 37. The method employed to demonstrate the difference in stability between the large and small masks. The weight represents a 5N force, and it is placed onto a pin to observe the resulting displacement of the small Asian mask. PAGE 54 54 CHAPTER 4 METHODOLOGY FOR MEASURING STABILITY Introduction The methodology for testing stability of fit of the patientspecific frames relies on the use of the finite element method (FEM) to model the contact interactions between the patientspecific frame and the patients facial contour. The objective of the me thodology is to produce a score value that correlates with the physical st ability of the six masks created. The methodology quantifies the stability of fit by information from the simulated movement of the patientspecific frame as it responses to small external forces The virtual mask is able to experience displacements because it is modeled as a rigid so lid surface, which is virt ually placed on top of a deformable bed of springs (Figure 41), and loaded by small external forces. The small external forces are calculated to cause a spring deflection of 1%, becaus e the simulation does not include friction. For every small external force applied onto the rigid layer, intern al forces transfer onto the bed of springs, and by the deformation of the sp rings, the rigid layer can either move closer to the model in some areas and move away in other areas and still sust ain contact; or it can predominantly slide or rotate (F igure 42). Sliding and rotation occur when there is relative motion between the two surfaces. The simulated m ovement of the mask is analogous to testing the fit of a cast by placing it over its mold and applying pressure at various locations to see if it fits securely. In a similar way, the approach us ed in this investigati on involves applying small forces in many directions and locations and m onitoring the resulting moti on of the mask (cast). As described in the materials section, three solid lifelike facial models that represent three different patients anatomies we re constructed with 3dimensi onal printing rapid prototyping. For each solid model two masks were designed w ith the existing user interface RPD designer, one of the masks is known to be stable and the other is known to be unstable (from information PAGE 55 55 on Chapter 3). The stability information on the real masks was used as training data for the methodology designed in this investigation. The methodology designed to test stability of the patientspecific reference frames is presented in four sections; a) system component s and FEM setup, b) program parameters, and c) output parameters, and d) the criteria defined to measure stability of fit of the patientspecific reference frames. As mentioned before, the ta sk of measuring stability is critical to understanding the appropriateness of the referenc e surface selected and thereby appreciating the predicted accuracy and precision of the clinically designated patientspecific design. While this is only a piece of a much broade r investigation that seeks to ev aluate an alternative method to imageguided surgery, it is critical to this application of patients pecific tool design. System Components and FEM Setup The system components are the external forces, the rigid layer (virtual mask or frame), the bed of springs, and the boundary conditions; all are defined in the context of FEM. The STL file, output of the graphical user interface (RPDdesigner), is the onl y input into the system. The STL file is standard format representation of a 3dimensional surface as an assembly of planar triangles (Figure 43) [29], the point coordinate information of the vertices of each triangle are used to define the rigid layer, the be d of springs and the boundary conditions. The first component of the FEM setup is a ri gid layer and it is com posed of tetrahedral finite elements. This layer is virtually cons tructed between the selected surface and a surface extrapolated along the positive normal direction at each vertex (Figure 41 and 43). The second component of the system is a bed of springs com posed of individual linear springs that are placed between the selected surface and another su rface extrapolated along the normal negative direction of each vertex (Figure 41). Hence, each individual sp ring in the bed of springs is constructed between the point coordinates of the selected surface and the boundary condition PAGE 56 56 points (or nodes). On the whole, the rigid laye r represents the patientspecific reference frame (or mask), and the boundary conditions represen t the patients facial surface. With the stimulation of external forces, the bed of spri ngs will allow for the virtual displacement of the rigid layer, simulating the stability or instabil ity of the real mask onto the patients facial anatomy. External Forces The external force direction is defined with respect to a main normal (Nmain). Nmain is a normal perpendicular to the coronal plane and al ong the sagittal plane of the skull; in other words, in the direction in which the frame will be placed onto the model (Figure 44). The angle theta ( ) is the tilt angle from Nmain and is the spin angle about Nmain. The forty nine forces are calculated by incrementing theta ( ) six times, from 0 to 90 in intervals of 15; and by rotating the angle (spin) by eight forces vectors in increments of 45 about the Nmain (six tilt angles multiplied by eight spin angles plus one normal for ces is equal to forty nine forces vectors). All forty nine forces are applied onto every node of the mask with the exception of the edge nodes (Figure 45). The force magnitude is very sm all enough to displace a single spring by 1% and will be discussed in the program parameters section. For example, in a case of a mask similar to the one shown in Figure 45, which is compos ed of about 300 nodes, 49 force vectors will be placed on each node, that will result in 14700 forces applied onto the entire mask. For each of the 14700 force vectors, all 300 springs w ill have a resulting displacement. Rigid Layer Probably the most important ch aracteristic of the system is the difference in strength (elastic modulus) between the linear springs, which compose the bed of springs, and the rigid layer. The rigid layer in comparison to the bed of springs must be rigid in the sense that it can PAGE 57 57 translate or rotate but not deform. The rigid laye r is composed entirely of four node tetrahedral elements. From the extrapolation of the point coor dinates that define the planar triangles in the STL file, virtual wedges are built. In each wedge three tetrahedra l elements are fitted (Figure 46), and created following the right handrule as explained in Figur e 215. In order to define a four node linear tetrahedron in ITKFEM module [24], four point coordinates are need, and two mechanical properties, the modulus of elasticity (E) and a Poisson ratio ( ). RDR ratio. It is important to measure if the ri gid characteristic of the rigid layer are maintained at every force application; for that purpose a program flag was created and called RDR ratio. The RDR ratio is an indicative of major deformations of the rigid layer. The RDR is calculated with respect to the upper surface of the rigid layer and for each external force application because every external force applic ation has the potential of deforming the rigid layer. To calculate the RDR the maximum differ ence in distance between a vertex (node) i and its neighboring vertices before (La) and after (Lb) the application of a force is measured and called maximum rigidity fault (mrfi) (Figure 47). As Equation 41 shows, the RDR is the ratio of the maximum rigidity fault (mrf) at a ve rtex i divided by the su m of the displacement magnitude (di) and the mrfi. i i i id mrf mrf RDR (41) If RDR is close to zero (or one over a small displacement), it means that the rigid layer has not deformed. If the RDR ratio is close to one, th e rigid layer has taken the entire load and it has deformed at a particular node. In the situation where the rigid layer deforms, the RDR ratio will indicate so and the simulation will not be able to pr ovide reliable results. In the simulation, if the RDR ratio exceeds 0.05 (or 5% rigidity loss) then the force was not included in the results. For PAGE 58 58 all the masks and frames evaluated, the RDR rati o had a mean of 0.3% and a standard deviation of 0.2%. Bed of Springs The bed of springs is the only deformable component of the system given that the individual springs that comprise it considerably are less strong th an the rigid layer. The purpose creating a bed of springs is to id entify and monitor the individual springs that acquire a state of tension or compression. A linear spring has a characteristic spring stiffness referred to as the spring constant k [46]; it is related to the strength of the ma terial of what the spring is made of. In addition, the spring has a natural length (L0) that is the length of the spring at rest (Figure48a). Figure 48b shows an axial force FY placing a spring under a state of compression and changing its length to a new length L1. Figure 48c shows the fr eebody diagram of the system and the spring internal force FS, which is determined by the e quation in Figure 48d. In the Finite Element Method, linear spri ngs that reside in 3dimensiona l space are modeled as trusses elements in 3dimensional space. In order to define a linear spring in 3dimesions in ITKFEM module, 2 point coordinates are needed, a modulus of elasti city (E), a Poisson ratio ( ) and an area (A). The area for each spring element i is the sum of one third (because each triangle has three vertices) the areas of all triangles surrounding the vertex i and denoted in Equation 41 as Ai. In this manner, even though the triangular mesh in the STL file is irregul ar (Figure 43); all the springs in the entire bed have the same proportional strength, because the spring stiffness coefficient for each spring element (ki) is set proportiona l to the area Ai, as seen in Equation 41 [46] a the elastic modulus (E) and the spring length (Lo) are kept constant. o i iL E A k (42) PAGE 59 59 Program Parameters Input Parameters In addition to the STL file, the program requi res five input parameters that are set for convenience of coding and not to achieve a particular real world test, they are: the natural length (L0) for all linear springs that compose the bed of springs, external force ma gnitude, the height of the rigid layer (h), the modulus of elasticity for the spring elements (Es) and the modulus of elasticity for the rigid surface (Er) [46]. Spring height The spring height is set to 0.1 mm, primordia lly because the displacem ents are expected to be very small, close to 0.001 or 0.002 mm. Th e system was design to include linear spring element which do not carry friction and when grounde d in all three directi on they behave as a ballandsocket joint. External force magnitude The external force magnitude is calculated by Equation 43 and selecting an expected axial displacement ( L) of 1% the springs length. As s hown by Equations 32 and 42, the force required to move a spring a certain displacement is proportional to the cross sectional area of the spring. In a similar way, in order to calc ulate a average external force magnitude that corresponds to a single springs displacement of 1%, the force must be weighted by the masks surface area per number of springs, as shown in Equation 43. L L E spring area spring ce AverageFors* *0 (43) PAGE 60 60 Rigid layer height The height of the rigid layer was set to 1.0 mm because increasing the high of the rigid layer increases the rigidity of the frame, a nd 1.0 mm was found to be the appropriate value (Table 41). The rigidity of the frame was qua ntified with the RDR ratio and Table 41 shows that reducing the height of the rigi d layer considerably affects the ri gidity of the rigid frame. For example, in Table 41, a RDR of three over one thousands signifies a 0.3% rigidity loss. In addition, the rigidity is affected by the shape (l ength and height) of the tetrahedral elements [41]. Since STL information is used to mesh the rigi d tetrahedral layer, a nd the STL facets are not regular triangles, the shape of the tetrahedral elements affects the rigidity; by increasing the height of the rigid layer the tetrahed ral elements behaves more rigidly [41]. Difference in elasticity The elasticity of the rigid layer (Er) is set to have a modulus of elasticity six orders of magnitude higher than that for the bed of springs (Es). Six orders of magnitude is the largest difference that can be established without comp romising the numerical stability of the global stiffness matrix. The sparse matrix linear solver Cholmod [13] has an output value which signifies the conditional number of the matrix that relates to the num erical stability of the system of equations, called rcond. It was found that un reliable values resulted when rcond is equal or exceed 1e13. A difference of six orders of magnitude for the Er and Es maintains the rigidity and results in a reliable simulation. Table 42 shows the how the difference in elasticity between Er and Es affect the RDR ratio and Cholmods rcond value. Output Parameters For every external force applied onto the rigid layer, there will be a resulting displacement vector (di) at every spring i that make up the bed of springs (Figure 49). The displacement vector for each individual spring is characteri zed by two components, the normal displacement PAGE 61 61 along the axis of the spring and the tangential di splacement perpendicular to the axis of the spring. The normal displacement scalar (or distance) (ni) is calculated by the dot product between the direction of the spring (Ni) and the displacement vector (di), (Figure 49). The tangential displacement scalar (ti) is calculated by the cross produc t between the direction of the spring (Ni) and the displacement vector (di) divided over the norm of Ni. As Equation 44 shows, the compression or tens ion state of a linear spring is determined by its normal distance (n). If th e normal distance is less than or equal to zero, the spring is known to be in compression, and if the normal dist ance is greater than ze ro, the springs is know to be in tension. In a similar manner, regardless of the direc tion of the tangential displacement vector, the tangential distance ti is the perpendicular distance from the axis of the spring (Ni) (Equation 45 and Figure 46) to the final position of the spring (point c on Figure 49). Linear springs have an inherent characteris tic of only responding to axial forces, their internal forces (which are triggered by an a pplication of an external force) should only be measured by it normal displacement (n). This means that springs (in the simulation) that predominantly have tangential displacements ar e not developing internal forces and are not contributing to the stability of the frame onto the patients model. For example, in Figure 410, the external force (F) applied onto a rigid box re sults in the development of internal forces in the springs aligned with the external force (F), or springs 1, 2 and 3. The tangential displacement and change in length of springs 4, 5 and 6 is due to their attachment to the rigid box, and though ti = di x Ni/Ni i is the number of springs (45) ni compression = di Ni <= 0 (44) ni tension = di Ni > 0 i is the number of springs PAGE 62 62 the spring is stretching, its normal displacements is very small (insignificant) and therefore they do are not developing noteworthy internal forces. However, the tangential displacement of an individual spring is a good indica tor of the sliding, slipping or rotation of the frame on the patients model because it gives evidence on the internal forces being developed by other springs. When there are few springs aligned wi th the external forces the mask experiences tangential displacements and these are expected to be high. In accordance with the inherent characteristic of linear springs, the fi nal length of a spring (L in Figure 49) is not a good indicator of the internal forces developed by a linear spring as a response to external forces; how ever, the final length does help s in determining which springs are exceeding the bounding compression region as shown in Figure 411. Criteria For this investigation two criteria were defi ned. The first criteri on is a single spring criteria, which determines if the resulting displace ment of a single spring is acceptable to aid in the definition of stability. The s econd criterion is a single force cr iterion that looks to define if the frame sustains contact and has a unique fit at the onset of a si ngle external force. With the information from the single force criterion an overall instability score for each mask is determined. This is done by providing a percenta ge ratio between the forces that result in instability and the total number of for ces that were applied onto the mask. The criteria was designed with the expecta tion that for the stable mask cases, the displacements of the springs will fall on the upper area of the sphere in Figure 411, because the external forces applied in the simulation are smal l. Accordingly, it is expected that the stable masks would not experience full compression or co me close to the boundary condition surface. For this reason a ratio called the TNT ratio, which is an indicative of the tangential displacement to normal displacement was define to aid in differ entiating between springs that are in a state of PAGE 63 63 compression, but experience mainly tangential displacement. The TNT ratio is a part of the single point criterion, its definiti on, the determination of a thres hold value and its consistency is discussed in the following sections. TNT Ratio The TNT ratio is a normalized measure of the tangential movement of an individual spring with respect to its normal displacement. The TNT ratio is an indication of the fractional slippage of an individual spring. It is calculat ed from each individual springs tangential (ti) and normal (ni) displacements (Figure 46). Numerically, it is defined by Equation 46 as the tangential displacement magnitude (t) divided by the sum of the normal displacement magnitude (n) and the tangential displacement magnitude (t) (Figure 49). n t t TNT (46) The purpose of the TNT ratio is to differentia te, for small displacements, between springs that are predominantly sliding to springs that ar e predominately in compression. In Figure 410 a force (F) is moving a rigid box by a unit displacement in the negative x direction and is placing three horizontal springs (springs 1, 2, and 3) in compression; since the three horizontal springs are experiencing only compression (n = 1 and t = 0) the TNT ratio for all of them is equal to zero. Assuming the box is perfectly rigid, the vertical springs (spri ngs 4, 5 and 6) will be mostly sliding, (n = 0, t = 1) and their TNT ratio w ill be equal to 1 by Equation 46. After the application of the force (F), the horizontal springs will be the only ones developing internal forces to sustain the movement to the rigid box. The vertical springs in Figure 410, are not contributing to with their internal forces to th e equilibrium of the rigid box, since linear springs are characterized by not responding to perpendicular forces. The displacement of these springs (springs 4, 5, and 6) result from the sliding movement of the rigid box. The TNT ratio PAGE 64 64 distinguishes springs that slip (or move tangentially) like springs 4, 5 and 6 in Figure 410, from springs that compress, like springs 1, 2 and 3 in Figure 410. Furthermore, Figure 413 shows two springs in compression; Figure 413A shows a spring with a TNT ratio of 0.5 and Figure 413B shows a TNT ratio of 0.833. In the case of the latter spring, its normal displacement is in the compre ssion direction (n<=0) and final spring length is smaller than the original length (L PAGE 65 65 and the area under the incline line, which represents a TNT ratio set to 0.70. The slope of the incline line (blue) corresponds to a normal disp lacement of 1 and a tangen tial displacement of 2.41, which corresponds to a absolute TNT ratio of 0.7. As mentione d before, it is expected that the springs have small displacements, so the TNT ratio of 0.7 (inclined line) will differentiate between springs that either have a very sm all normal displacement and a large tangential displacement or springs that have a very sma ll normal and tangential displacement; in both case these springs are not considered good compresso rs because they are not contributing to the stability of the rigid frame by not developing considerable inte rnal forces. The TNT threshold was set to 0.70 because it was defined empirically to correlate with the in formation about the real stability of the large masks and the instability of the small masks. For example, Figure 413B shows that there are occasions where the norma l displacement shows a compression state and yet the tangential displacement is high; by defini ng TNT upper threshold of 0.7, a spring like in Figure 413B will be excluded from being considered a good compressor, it has a TNT ratio of 0.833. TNT threshold The threshold for the TNT ratio was defined em pirically to be 0.7 based on the training data provided by the six masks constructed with rapid prototyping and thei r fitting onto the solid models (Figure 71, 72, and 73). To determ ine this empirically the simulation was ran with TNT ratios of 0.6, 0.7 and 0.8 for large and sma ll masks (Figure 415 and 416 respectively). In the simulation a TNT ratio of 0.70 corresponds with the physical st ability of the large masks and instability of the small masks. Figure 415 shows color coded surface plot for all three large masks, Caucasian, Asian and infant at different TNT ratios (0.6, 0.7 and 0.8). The location of the colo r dots correspond to the location where a virtual force was applied and three or more points behaved as good PAGE 66 66 compressors (Table 43). The color of the dot s corresponds to the maximum tilt angle (for all spin angles) where at least three good compre ssors were found. For example, a red dot signifies that for a maximum tilt angle of 90 (and for all spin angles) three or more good compressor springs were found, that is 49 forces resulted in at least three good compressor springs. The first row in Figure 415 asso ciates to a TNT of 0.6, and as the figure shows, the Asian and the infant mask are noted as false negatives, because the color coded surface plots show that the maximum tilt angle in which the forces was wit hheld is smaller than 90, this indicates that the application of a force resulted in an unstable configurations (this information contradicts that premise that the large masks ar e stable). The second and thir d rows, which correspond to TNT ratio of 0.7 and 0.8 respectively, concur with the physically stability of the masks; but in Figure 416, a TNT threshold of 0.8 dictates a case of false positive for the small Asian mask, which is mark as maintaining alignment for all forces at all locations. This information also contradicts the premise that the design of the three small masks results in a unstable configurations. Figure 414 shows that when a tangential displacement magnitude is 2.41 times bigger than the normal displacement magnitude the TNT ra tio is equal to 0.70. Also a TNT ratio of 0.5 occurs when the tangential displacement is equa l to the normal displacement, and the TNT ratio approaches one when the springs experience mos tly tangential displacements (Figure 414). In summary, a TNT upper threshold was empirically define d at 0.7, this means that all springs that are in compression and exhibit TNT ratios of 0.7 or below are in a good compressive state. On the other hand, the springs that exceed a TNT ratio of 0.7 will be regarded as mainly slipping or sliding. PAGE 67 67 Consistency of TNT threshold To test the consistency of the TNT threshold, the force magnitude was increased five and ten times, and the number of point on the mask was increase for a single mask case (Figure 417). When the force was increased by five a nd ten times there were no changes in the color coded surface plots shown in Figure 414. This m eans that the TNT threshold is not sensitive to force. In the case of increase number of point s, the results are shown in Figure 418; which illustrates that with the exception of one single point in the upper left corner of the Asian large mask, which held stability up to a theta value of 75 the TNT threshold of 0.7 corresponded with physical results. Single Force Criterion The purpose of the single force criterion is to define on a per force basis if the frame or mask is stable or not. The criterion requires a minimum of 3 points in compression, that is 3 good compression points (as defined by the single poi nt criteria) to consider a mask stable with relation to a particular force vector. Three spri ngs in compression will overcome for translations and rotations of the frame. Translation can be recognized by evaluating the movement of a single point, but two points are needed to recognize rotation, since rotation occurs about a normal vector and at a specific point that may not move. Translations or sliding is defined as the relative motion between two surf aces and rotation is the angular motion about a common normal [45]. The criterion requires a third point because three points are necessary to identify rotation in 3dimensions. If the displacement of the rigid body results in three good compression points, the simulation indicates that the rigid surface is moving close to the grounded points rather than sliding or rotating, which in turn indicates that the mask and model surfaces remain in alignment. If by the onset of an external force the mask re sults in a configuration where there are less the PAGE 68 68 three good compression points, then the mask is know n to be unstable for that particular force. Intuitively, more than three good compression points and their dispersi on can give evidence of a level of stability. More good compression po ints and broadly disperser leads to greater stability. Therefore, a bounding volume de fined by the minimum volume of all good compressor points compared to the bounding volume of the entire mask, is a good indicator of the level of stability of the ma sk for a particular force. Box ratio: The box ratio quantifies the level of st ability as defined by the single force criterion by calculating a ra tio between the volume of the compression patch (cvol) and the volume of the entire frame model (fvol) (Figure 418). The box ratio is calculated on a per force basis. The compression patch volume is th e minimum volume that includes all good compression springs (Figure 419). If the mask is unstable, with respect to a particular force (as defined by the single force criteria), then the comp ression patch volume is equal to 0. The small box in Figure 418 shows a compression patch volume (cvol), in comparison with the volume of the entire frame (fvol). The box ratio is numerically defined by the compression patch volume cvol over the entire frames volume fvol (Equation 47). vol volf c ratio Box (47) As Equation 47 shows, a box ratio is equal to 0 when the compression patch volume is equal to zero (or the mask is unstable) (Figur e 419). When the mask is found stable for a particular force, a compression patch volume take s a value greater than zero. A box ratio close to 1 indicates that the compression patch volume approximates to the volume of the entire frame (Figure 419). Box ratios clos er to the value of zero corres pond to smaller compression patch volumes or good compressor points close to each other. PAGE 69 69 A score for the overall instability (instability score) of the masks can be derived from the single force criteria by providing w ith an percentage value for the number of forces that result in instability (box ratio equal to 0) to the total numbe r of forces applied to the mask (Equation 48). If this value is equal to 0%, there were not extern al forces that made the mask loose stability, and if this value is 100% it means th at all external forces disrupted the alignment of the mask onto the patients surface. 100 mask onto the es applied er of forc total numb bility t in insta the resul #of forces y score Instabilit (48) In summary, a total of 49 external forces will be applied to all none dge points onto a mask (if a mask has 300 nonedge points, ther e will be 14700 forces applied to it). For each individual external force applied onto a mask, there will be displacements of all springs that comprise the bed of springs. The single point criteria will individually classify each spring as a good compressor. The number of good compression points is recorded. If there are less than 3 good compression poin ts, the mask is considered unstable, the compression patch volume is equal to zero, and the box ratio is equal to one. If there are 3 or more good compression spring s then the mask is c onsidered stable and the compression patch is calculated and the box ra tio will have a value between 0 and 1. As an overall percentage instability score for a mask, the number of forces the result in instability is compared to the total number of forces applied onto the mask. A value of 0% means zero instability (Equation 48). PAGE 70 70 Table 41. Relationship between the tetrahedral he ight and a measure of th e rigidity of the rigid layer. Table 42. Table that shows th e relationship between the difference in elastic modulus between the rigid layer and the bed of springs a nd the rigidity measure and Cholmods rcond value. ErEs Mean RDR Cholmods rcond 1.0e6 3/1000 2.14e8 1.0e7 2/1000 2.14e9 1.0e8 5/1000 2.14e10 Table 43. Singl e point criteria. Compression Tension No Slip n <= 0 and L <= L and TNT <= 0.70 good compressor n > 0 or L > L Slip n <=0 and L PAGE 71 71 Figure 41. The methodology for stability of fit involves simula ting a rigid layer onto of a deformable bed of springs. A) The sele cted surface (red) as painted on the rendered model (green). B) A cross section thr ough the selected surfaces and anatomical model showing the idealizati on of the contact mechanics between the rigid frame and the solid model Bed of s p rin g s Boundar y Conditions Outer extrapolated Surface A Selected Surface B Ri g id la y er PAGE 72 72 Figure 42. Displacement of a mask as it responds to external forces. Figure 43. Triangulation of a customized pos itioning device proper to STL format. Every vertex in the mask is a point coordinate in 3D space. F F Springs in tension Springs in compression PAGE 73 73 Figure 44. Direction of the Nmain. Nmain is specifically chosen for each mask, and forty nine force directions are calculated by de fining a tilt and a spin angle. Figure 45. Nodes point on th e rigid frame. Forces are not applied on edge points. Edge points Sagittal plane Nmain Coronal p lane Tilt ( = 0 to90) Spin ( = 0 to 360) Nmain PAGE 74 74 Figure 46. In a wedge three tetrahedral elemen ts are fitted. The nodes are shown next to the black dots. The numbering of the tetrahed ral elements is important and should be done as shown in the legend. 4 3 0 1 2 5 Tetrahedral Nodes 1 0123 2 1234 3 2345 Tetrahedral 1 4 3 0 2 5 1 Tetrahedral 2 Tetrahedral 3 PAGE 75 75 Figure 47. Schematic showing how the maximum rigidity fault is calculated. Figure 48. The behavior of a lin ear spring element under th e application of an axial force. A) A spring at rest, nl is the natural length of the spring. B) The same spring under compression by a force Fy. C) Free body diagram of the spring in compression. D) Equation of a linear spring. Lb La Maximum rigidity fault Vertexi FY L1 k L0 ) (0 1L L k FS y x FY FS A B C D PAGE 76 76 Figure 49. Schematic that demonstrates how each springs displacement was evaluated by identifying the normal (ni) and tangential (ti) displacement magnitudes. The normal (ni) and tangential (ti) displacement are scalar distances In this particular case the spring is in compression because the normal displacement magnitude ni will result in a value less than or equal to zero, also no tice that the natural le ngth of the spring has changed from L to L. d di B ni =di Ni A C N L ti =di x Ni/Ni L PAGE 77 77 Figure 410. Schematic that shows why it is impor tant to distinguish between the springs that are mostly in compression or tension and thos e that are basically slipping or have a high TNT ratio. The force (F) exerts a unit displacement on the rigid box in the x direction. As a result of the force the springs have normal (n) and tangential (t) displacements are shown in table. Figure 411. Bounding volume for a spring to be in compression. For the stable cases the springs displacements are expected to reside on the upper area of the half sphere (gray area). F y x 1 2 3 4 56 Springs n t TNT 1, 2, 3 1 0 1 4, 5, 6 0 1 0 L<=L n<=0 L n t PAGE 78 78 Figure 412. Twodimensional interp retation of the single point criteria. In order for a spring to be considered a good compressor it should re main inside the red area. Consider the spring in Figure 413A to correspond with spring 1 (with the tangential displacement in the opposite direction), and spring 2 co rresponds to the spring in Figure 413B. Figure 413. Schematic to demonstrate the importa nce of evaluating the percentage slip of TNT ratio for each spring. A) The TNT ratio in this case is close to 0.5 because the compression and tangential displacements are similar in length. B) The TNT ratio is close to 0.833, because the ratio of comp ression to tangential displacement is about 1 to 5. n = 1 t=1 n = 1 t = 5 A B TNT = 0.5 TNT = 0.833 L<=L n<=0 TNT<=0.70 or n/t = 0.41 1 2.14 L n t 1 2 PAGE 79 79 Figure 414. Schematic to shows at the leve l at which the TNT threshold was defined. Figure 415. Empirical determina tion that the TNT threshold greater than 0.6 correlates with the training data, which says that the large masks are stable. TNT =1 t>>n TNT = 0 t< PAGE 80 80 Figure 416. Empirical determina tion that the TNT threshold set smaller than 0.8 correlates with the training data, which says that the small masks are unstable. Figure 417. A case of the large Asian mask, which shows the consis tency of an upper TNT threshold of 0.7. Even when the number of springs is increase, the TNT threshold of 0.70 shows stability. 90 75 60 45 30 15 0 Theta Decimate15 Decimate10 TNT = 0.8 TNT = 0.7 TNT = 0.6 Large Asian 90 75 60 45 30 15 0 Theta TNT = 0.7 TNT = 0.8 TNT = 0.6 Unstable Small Asian Small Caucasian Small infant PAGE 81 81 Figure 418. Schematic showing the relati on between the volume of the entire frame (fvol) and the compression patch volume (cvol). Figure 419. Box ratio interpretation. Notice that a box ratio of one is the most stable configuration where the compression patch volume (cvol) is equal to the entire volume of the frame (fvol), meaning that all points are g ood compression points. Also notice that a box ratio of zero corresponds to instab ility; since there are less than three good compression points, which result in a co mpression patch volume equal to zero. 1.0 0.0 Instability (>3 good compressor points, cvol = 0) Box Ratio 0.5 cvol = fvol cvol = 3/4 fvol cvol = 1/2 fvol cvol = fvol cvol = 0 0.75 0.25 Stability(<= 3 good compressor points, cvol 0) Most stable PAGE 82 82 CHAPTER 5 PROGRAM DESCRIPTI ON AND VALIDATION Program Description The program was written in C++ language a nd it uses open source programming libraries for image visualization: VTK (V TK 4.2, Kitware Inc., Clifton Park, NY) and, for advanced imaging processing, ITK (ITK 1.8, Kitware Inc., Clifton Park, NY). The library for advanced imaging processing has several classes dedicated to finite element method (FEM), which were used to construct the simulation of the contact problem. In addition Cholmod [13] was adapted to solve the linear system of equations particul ar to the FEM; also Paraview (Paraview 1.6, Kitware Inc., Clifton Park, NY) was used for the visualizations of the results. The program is organized into five main sectio ns. The first section is dedicated to reading the input information from the S TL file, introducing input paramete rs and preparing the data to be written into the appropriate input format ITKFEM module. The s econd section requires the ITKFEM module to read the input file previously created and a ssemble a single global stiffness matrix (Kg) in the correct format for Cholmod. The third section defines the directions of the force at each point on the rigid layer and iterates over all points while storing into arrays the output parameters. The forth s ection sorts over the information in the output parameters and applied the single point and global criteria; and finally the last section presents the results, in a histogram format and as a colormap in Paraview. Figure 51 is a chart fl ow of how the program is organized. The first section of the program starts by read ing the STL file and applying an ITK filter that selects the largest region in the STL file. This is done because sometimes additional small surfaces are selected in RPD designe r that are not related to the main selected surface. In general the STL surface is characterized as an irregula r mesh. The STL serves the purpose well of PAGE 83 83 describing a surface by dividing it in planar triang les. For this reason, another VTK filter, a decimation filter is applied to reduce the number of triangle that are close to each other and have the same dihedral angle [38]. The decimation filter does not alter the topology of the selected surface; rather it serves the purpose of reduci ng the number of triangles in the STL mesh [38]. However, sometimes the output of the decimation filte r results in relatively large triangles, so an optional function was created that added nodes to large triangles and rem eshes the surface. At this level in the program the surface does not undergo additional modification, and the nodes on the edges are saved to an array for later identificati on. The next step is to calculate the normal at every vertex. VTK calculates the normals at ev ery vertex on the surface as an average of the normals of the triangles connected to the vertex ; however a function was created to calculate a better average normal, which is weighted by the angl es of the triangles adjoining to a particular vertex. With the appropriate normal inform ation a top surface and bottom point cloud are created by extrapolating the point coordinate information along the normal in the positive direction and in the negative direction, respectively, as shown in Figure 52. The rigid layer is created by meshing with tetrahedral elem ents the space between the selected surface and the top surface (Figure 52). The bottom point cloud become boundary condition points, fixed in all 3dimensions and represents the patients surface. The spring elements that comprise the bed of springs extend from the point cloud to the selected surface (Figure 52). Next the program calculates the su rface area of the frame and with the expected displacement se to 0.001mm, determines the magn itude of the external force (Equation 52). The last step of the first sec tion of the program involves using the point coordinate information and input parameters to write an input file for the ITKFEM module. PAGE 84 84 Most of section two involves us ing the ITKFEM classes to read the file created, assemble a single the global stiffness matrix (Kg) that describes the entire syst em and the global force vector (Fg). Afterwards a function was created to read th e global stiffness matrix and write it into an input format for Cholmod. The format required by Cholmod is referred as triplet format, and is a method of writing large sparse matrices in manner that saves computer space and quickens the retrieval of information [13]. The direction of forty nine force vectors is calculated in a combination of forty nine tilts ( ) and spins ( ) about the main normal (Nmain) as shown in Figure 412. The same small magnitude is assigned to each force vector and all forty nine are applied onto each nonedge nodes on the mask. One global stiffness matrix (Kg) is assembled for the entire system and is not altered once it is assembled; however, for each force applicati on to a particular node (49 times the number of nonedge points) on the rigid la yer, a global force vector Fg is populated. In this manner, both the global stiffness matrix Kg, and the global force vector Fg are given to the Cholmod routines to solve for the displacement vector U that characterizes the movement of the frame as a result of a single force application. From the displacement vector U, which include displacements for each nodes, the normal (n) and tangential (t) dist ances are calculated for each node as shown in (Figure 48). In section five, the single point criteria a nd the global point criteria, as described in Chapter 4, are applied to the tangential and nor mal scalars. In addition, the program creates histograms for the box ratio per angle rotation and creates a file for visualizing the frame with colorcoded information per angle of the maximu m angle in Paraview (Paraview 1.6, Kitware Inc., Clifton Park, NY). PAGE 85 85 Validation of FEM elements The validation of the program involves testi ng the behavior of th e linear spring or space truss element and the linear tetrah edral element to assure they re spond in the appropriate manner. In addition, the program is tested fo r a simple case of a corner frame. Linear Spring Validation The validation of linear spring element was pe rformed by testing the displacement of the element in three orthogonal directions as show n in Figure 53. In all instances the same mechanical characteristics (E, A, L) where assigned and the same force vector (F) was applied (Figure 53). The results of the stiffness matrix as assembled by the ITK FEM module are shown in Equation 51, 52 and 53, which correspond to Figure 53A, 53B and 53C, respectively. The resulting displacements for each case resulted in the sa me value of 0.1mm in the direction of the force vector, which corres ponds to the displacement as hand calculated in Equation 54, following the linear spring equation. 0 0 0 0 0 0 0 0 0 0 0 0 0 0 750 0 0 750 0 0 0 0 0 0 0 0 0 0 0 0 0 0 750 0 0 7501Ak (51) 0 0 0 0 0 0 0 75 0 0 750 0 0 0 0 0 0 0 0 0 0 0 0 0 0 750 0 0 750 0 0 0 0 0 0 01Bk (52) PAGE 86 86 750 0 0 750 0 0 0 0 0 0 0 0 0 0 0 0 0 0 750 0 0 750 0 0 0 0 0 0 0 0 0 0 0 0 0 01Ck (53) mm MPa mm mm N AE FL d 1 0 1500 1 0 2 752 (54) Tetrahedron Element Validation The validation of the tetrahedral element was done by comparing the ITK FEM results to those by ANSYS (Ansys Inc. Canonsburg, PA) fo llowing the example in Figure 55. Both ANSYS and ITK FEM provided the same displacement shown in Table 51. A Simple Model Validation The program was validated with a simple fram e in the shape of a corner (Figure 55). When a force was applied at the location marked with a white dot and in the direction of normal direction at that point; the entire corner was placed in a state of compression as will be expected if a force physically applied onto a corner (Figur e 55A). In the revers e case where the force vector is inverted and rather th an pull it pushes away, the corner results in a complete tension state. PAGE 87 87 Table 51. Displacement in the X, Y and Z direc tion from the example in Figure 54 tested with the ITKFEM module and ANSYS. Node X Y Z 1 0.0000 0.0000 0.0000 2 0.0000 0.0000 0.5200e3 3 0.0000 0.0000 0.0000 4 0.0000 0.0000 0.0000 5 0.0000 0.0000 0.5200e3 PAGE 88 88 Figure 51. Flow chart of program for testing the uniqueness of fit. Decimate filter STL reader filter Largest region filter Assemble Kg Calculate weighted normals function Add Points and remesh function Create new top surface function Create new bottom point cloud Write FEM file Read FEM file Write Kg in triple format Section I Section II Calculates the surface area and the force magnitude Assemble Fg INPUT: STL file PAGE 89 89 Figure 51. Continued. No Yes Get U and calculate: normal and tangent displacements Apply Criterion Single Point Criteria (for each point on surface) Calculate: TNT, #T, #C If (#C > = 3 && TNT>=0.75) {GoodCompressor++} Test for Flags: Positive Definiteness and RDR Global Point Criteria (for every force application) If (GoodCompressor >=3) {Calculate Box Ratio An g leRotateCount++} For all nodes (not edge nodes) Insert new force vector into Fg Rotate force vector and assign magnitude For = 0 to 90, = 0 to 360 (49 force vectors/node) Solve Fg=KgU with Cholmod Section III Section V Get main frame normal ( Nmain) Section IV = 90 = 360 ? More nodes? Yes End OUPUT to histograms and frame colormap PAGE 90 90 Figure 52. Original surface is extended in th e positive normal direction to form the top surface and the points extended in the negative normal direction to become boundary condition points. Figure 53. Examples of al li near spring element or space truss tested in all 3dimensions. Selectedsurface Top surface Bottom point cloud Rigid tetrahedral layer X Y E =1500 MPa = 0.3 F = 75 N A= 0.1 mm2 L = 2mm Z X Y Z X Y Z (A) (B) (C) L Cross Sectional Area (A) Element 1 PAGE 91 91 Figure 54. Example of two tetrahedral elemen ts used in the validation of the tetrahedral elements used in the program. F F X Y Z 1 3 2 5 4 E = 30e6 MPa = 0.3 F = 1000 N Element 1 Element 2 NodeCoordinates (X,Y,Z) 1 0, 0, 0 2 0,10,0 3 0, 0, 10 4 1, 0, 0 5 1, 10, 0 PAGE 92 92 Figure 55. A corner frame that was used to va lidate the program. A) The corner frame is shown on blue and the white dot signi fies the location of the for ce. B) The direction of the force is the normal located at 135 from a ll three sides of th e corner. C) The wireframe arrows show a scaled normal displacement vector for all nodes on the frame, and the small solid arrows show the displacement vector at each node. D) The force application in the opposite direction re sults in all normal displacements in the opposite direction to C. A B Force application point Force direction 135 C D PAGE 93 93 CHAPTER 6 RESULTS AND DISCUSSION A set of criteria were develop that showed that all large mask s were found to be stable and all small masks were not stable. Also, the caplike reference frame and the half sphere surfaces were found to be not stable. In reality, all masks maintained stability when loaded in the direction and location of each masks main normal ( Nmain). In reality the conforming surfac es of all small masks do not have a unique fit, meaning that there is not one unique location where the mask fits, rather there are various locations were it coul d fit (Figures 61B, 62B and 63B). The small masks were designed with the expectation that they will not have a unique fit, for the purpose of validating the methodology. For all large masks there is a un ique fit (Figures 61A, 62A and 63A). All large masks can withstand real lateral forces or forces 90 from their main normal ( Nmain). The Asian small mask is the only small mask that can sustain contact stability when a real load is applied at large angles from its main normal (about 60 from Nmain); it is the only small mask that extends laterally along the en tire superciliary arches (Figure 62B). It also extends caudally onto the nasion and sits onto th e upper orbital area. The small Caucasian mask does not extend onto the orbital area and does not extend onto the na sion (Figure 61A). The small infant mask extends laterally along the superciliary arches but it does not reach over the nasion (Figure 63B). The facial angles of all thr ee anatomical models were meas ured as defined by Figure 211; the results are shown in Table 61. One can appr eciate that stability coul d be sustained if the nasiolateral and nasiofront al angels are relatively small and if the nasiofacial angle is large. With smaller nasiolateral and nasiofront al angles there is more opposing contact area to lateral forces PAGE 94 94 (90 from the Nmain). The Caucasian model has the smalle r nasiolateral angl e, and the second smallest nasiofrontal an gle and the largest nasiofacial angle (Table 61). For the simulation results, two t ypes of graphs were produced, a) a color coded surface plot of the masks onto the solid models (Figure 65 to 69) and b) a normaliz ed histogram for the box ratios for all values of ( = 0 to 360), for each tilt angle ( ) (Figure 512) accumulated over all points (Figure 610 to 617). In the case of the surface models, the result s are presented for each force application point where the color signifies the maximum tilt angl e for which stability was found (based on the singlepoint and single stability criteria) for a full rotation ( = 0 to 360). At each color coded point on Figure 65 and 69, a se t of 49 force vectors were applie d. For each force vector all springs that comprised the bed of springs were evaluated for stability. At each spring the maximum angular force (for a full rotation = 0 to 360), the furthest from the main normal ( Nmain), at which at least 3 spri ngs were found in compression wa s noted. This maximum angle is color coded in Figures 65 to 69. If the fr ame was found stable at ev ery force application it will be indicated in red, demonstrating th at for any force between theta equal to 0 and theta equal to 90 the mask did not slide or tilt beyond the criteria established (T NT >= 0.7). Some masks (small unstable masks) experienced some tangential movement or rotation, so they only exhibit stability at a less er angle (i.e. theta = 45 from the main normal). For the box ratio histogram plots, the results are presented as frequency of the box ratio by individual tilt angle (theta or ) for all values over all points on the mask (Figures 610 to 617). For example, in Figure 610 for the Caucasia n large mask, for a tilt an gle of 0, 30% of the box ratios are equal to 0.4. PAGE 95 95 For the Caucasian large mask (Figure 65A), the surface plot show s that all nonedge points sustain stability for a maximum thet a of 90 and for all rotations about the Nmain ( = 0 to 360); while the small Caucasian masks design di d not sustain contact for any full rotation for any point (Figure 65B). For the Asian large mask (Figure 66A) and fo r the infant large mask (Figure 67A), the surface plots show similar results to the large Caucasian mask (Figure 65A), all points sustained stability all the way to a tilt angle of 90. It is important to notice that fewer forces were applied on the large Asian and large infant masks because there are less nonedge points (Figures 66 and 67). The surface plot for the small Asian mask show s more stability per location than any other small mask. This mask extends along the superc iliary arches and sits on top of the upper orbital area and extends caudally past the naison area pr oviding a surface to judge correct alignment. The color coded surface plot for the small Asian mask agrees with the design as they show stability along the left lateral side and over the u pper orbital area. On the right side less stability was found simply because the mask is not symm etric and includes less coverage over the right upper orbital area. The color coded surface plot for the small infa nt mask shows some point locations where the tilt angle theta sustains stability up to 60 (Figure 67). This is due to the mask extending over the upper orbital area and late rally along the superciliary arches. The single force criterion states that three po ints need to be in compression as defined by the single point criteria (with a TN T ratio of 0.70 or below) to have stability at a particular force. As mentioned before, if 3 or more good comp ressor springs are found, then the compression patch volume (cvol) is different than zero and the box ratio results in a value greater than zero and PAGE 96 96 less than or equal to one (Figure 419). The level of stability is the highest when the box ratio is close to 1, since this indicates that the volume of the compression patch (cvol) is similar to the volume of the entire frame (fvol). The box ratio histograms for the small Cauc asian mask (Figure 610) shows highest frequencies at box ratio equal to zero, this means that for most external force applications the mask was found unstable (from the single force criterion). These resu lts correspond to the surface plot results in Figure 65B. The box ratio histogram for the small Asian ma sk shows a range of box ratios between 0.0 and 0.55, with most frequencies (70%) at box ratio s of 0.15 (Figure 613), this means that when stability was found the co mpression patch volume (cvol) was small. However, for all theta values there are frequencies where the box ratio is equal to zero, this means that even though stability was found in most of the cases, there were always instances were the mask lost stability. This results are harder to relate to the surface color maps, because the as mentioned before, the surface color map show stability per location, and the box ratio histograms show stability by accumulating all forces for a particular theta va lue over all rotations. The infant small mask sustained stability for theta equal to 0, but had no stability for all other an gles (Figure 615). At theta equal to 90, 80% of the of the box ratios are equal to zero, which signifies a high level of instability for forces applied at 90 from the main normal. The box ratio histograms for the large Cau casian mask shows no instability (zero frequency for box ratios equal to zero), and the box ratios range between 0.25 and 0.95 (Figure 610). As the theta angle (tilt angle from the Nmain) increases the box ratios decreases, showing that at higher angles the comp ression patch volume is smaller. Similar box ratio histogram results can be seen for the Asian large mask (F igure 612) and the infant large mask (Figure 6 PAGE 97 97 14) where there is no instability or zero frequency for box ratios equa l to zero. However, for the large Asian mask the box ratios are skewed to th e left, they range from 0.25 to 0.65 (Figure 612). This means that for the large Asian mask, increasing theta values result in smaller compression patch volumes (cvol) than for the large Caucasian mask. For the large infant mask the results follow the same pattern, increasing the tilt angle of the external force results in decreasing box ratios (F igure 614). The range for the large infant mask is from 0.15 to 0.75, this is a slightly larger box ratios range than for the box ratio range for he large Asian mask (0.25 to 0.65) (Figure 612). Based on the facial angle measurements in Fi gure 211 and on the premise that the more pronounced facial features result in more stable masks, the id eal angle combination for a stable mask will be small nasiolateral and nasiofrontal angles and large nasiofacial angles. Out of the three large masks, the Caucasian mask corresponds to a model will the sma ller nasiolateral angle and the second smallest nasiofront al angle and the larges t nasiofacial angle. In addition, the box ratio range for the Caucasian mask is the largest (0.25 to 0.95) when compared to the Asian (0.25 to 0.65) and infant (0.15 to 0.75). Figure 613 shows a profile view of all three models. The infant large mask corresponds to a model with the smallest na siofrontal angle and the second smallest nasiolateral angle, but the smallest nasi ofacial angle. The infant model has the smallest nasiofrontal angle due to a curvature on the forehead as can be seen in Figure 618C. It is likely that this pronounced curvature is responsible fo r the infant large mask having the second highest box ratios ranging between 0.15 and 0.75. In addition to the mask reference frames, a caplike reference frame (Figure 68B) and a half sphere surface (Figure 69) were evaluated with the met hodology for testing the uniqueness of fit. In reality the caplike frame has an uncerta in fit, in other words, it was difficult to find the PAGE 98 98 location where it was designed to fit. The co lor coded surface plot results for the caplike reference frame shows that the frame sustains st ability (as defined by th e single point and single force criteria) all the way to a maximum theta value of 90 when forces are applied on the upper most cephalic region. However, forces no greater than 60 where sustained at the coronal and sagittal extensions of the fram e (Figure 68). The box ratio hist ogram (Figure 616) for the caplike frame corresponds with the information on the surface plots. The box ratio histogram shows that some level of stability was found for all th eta angles, but it was not sustained through out the entire surface, since about 5% of the box ratio show instability for all tilt angles. The remaining box ratios (about 70%) have values of 0.15, wh ich signifies that the compression patch volume (cvol) was small compared to the volume of the entire mask (fvol) for the instances were stability was found. In the case of the halfsphere surface (Figur e 69), the surface plot shows that there were no locations where the surface sustained any force, possibility not even in the direction of the normal force ( Nmain). The box ratio histogram for the half sphere test shows that stability was indeed found for all angles (there are box ratios between 0.0and 0.35); however, for the majority of the cases (60% to 80% freque ncy) no stability was found or th e box ratio equaled the value of zero. In summary, Table 62 shows the instability scores for all masks and frames as defined in Equation 48. There is zero per cent instability for all large mask s and there is some percentage of instability for all small masks, the caplike fr ame and the halfsphere surface, as was expected. The least stable case is the halfsphere test, with a instability score of 94%, this high percentage instability make sense because the surface has no irregularities to oppose any external forces (besides in the main normal direction). The sm all Asian mask was the mask with the smallest PAGE 99 99 instability score at 5%, as was mentioned previous ly; this is due to the design of the mask that includes the upper orbital area on one side and extends the furthe st medially and laterally along the superciliary arches. This results correlate well with the color coded surface plots, were stability was found in the left upper orbital region on the mask. Table 61. Measurement of the facial a ngles for the three anatomical models. Model Nasiolateral angle Nasiofrontal angle Nasiofacial angle Caucasian 89 116 40 Asian 123 124 29 Infant 104 106 23 Table 62. Instability score, pres ented as the percentage of forces that resulted in instability and the total forces applied on a frame. Frame Total Number of forces applied Total Number of forces applied Percentage Large Caucasian mask 4935 0 0% Large Asian mask 1575 0 0% Large infant mask 1702 0 0% Small Caucasian mask 874 602 69% Small Asian mask 1755 88 5% Small infant mask 874 334 38% Caplike frame 2695 241 9% Halfsphere test 1094 975 94% PAGE 100 100 Figure 61. Placement of Caucasian masks onto so lid model. A) Large Mask. B) Small Mask. Figure 62. Placement of Asian masks onto soli d model. A) Large Mask. B) Small Mask. A B A B PAGE 101 101 Figure 63. Placement of infant masks onto soli d model. A) Large Mask. B) Small Mask. Figure 64. Placement of a caplike refe rence frame onto an upper head model. A B Coronal p lane Sagittal plane PAGE 102 102 Figure 65. Surface plot color coded by the maximum full rotation angle for the Caucasian masks. A) Large Mask. B) Small Mask. Figure 66. Surface plot color coded by the maxi mum full rotation angle for the Asian masks. A) Large Mask. B) Small Mask. 90 75 60 45 30 15 0 Theta A B 90 75 60 45 30 15 0 Theta A B PAGE 103 103 Figure 67. Surface plot color coded by the maxi mum tilt angle for the infant masks. A) Large Mask. B) Small Mask. Figure 68. Surface plot colo r coded by maximum tilt angle for the cap reference frame design 90 75 60 45 30 15 0 Theta B A 90 75 60 45 30 15 0 Theta PAGE 104 104 Figure 69. Surface plots color coded by ma ximum tilt angle for the half sphere test. 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%Frequency (%) 0.000.150.250.350.450.550.650.750.850.95 Box Ratio Theta0 Theta15 Theta30 Theta45 Theta60 Theta75 Theta90 Figure 610. Normalized histogram of box ratios for each tilt angle ( ) for the Caucasian large mask. 90 75 60 45 30 15 0 Theta PAGE 105 105 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%Frequency (%) 0.000.150.250.350.450.550.650.750.850.95 Box Ratio S Theta0 S Theta15 S Theta30 S Theta45 S Theta60 S Theta75 S Theta90 Figure 611. Normalized histogram of box ratios for each tilt angle ( ) for the Caucasian small mask. 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%Frequency (%) 0.000.150.250.350.450.550.650.750.850.95 Box Ratio Theta0 Theta15 Theta30 Theta45 Theta60 Theta75 Theta90 Figure 612. Normalized histogram of box ratios for each tilt angle ( ) for the Asian large mask. PAGE 106 106 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%Frequency (%) 0.000.150.250.350.450.550.650.750.850.95 Box Ratio S Theta0 S Theta15 S Theta30 S Theta45 S Theta60 S Theta75 S Theta90 Figure 613. Normalized histogram of box ratios for each tilt angle ( ) for the Asian small mask. 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%Frequency (%) 0.000.150.250.350.450.550.650.750.850.95 Box Ratio Theta0 Theta15 Theta30 Theta45 Theta60 Theta75 Theta90 Figure 614. Normalized histogram of box ratios for each tilt angle ( ) for the infant large mask. PAGE 107 107 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%Frequency (%) 0.000.150.250.350.450.550.650.750.850.95 Box Ratio S Theta0 S Theta15 S Theta30 S Theta45 S Theta60 S Theta75 S Theta90 Figure 615. Normalized histogram of box ratios for each tilt angle ( ) for the infant small mask. 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%Frequency (%) 0.000.150.250.350.450.550.650.750.850.95 Box Ratio S Theta0 S Theta15 S Theta30 S Theta45 S Theta60 S Theta75 S Theta90 Figure 616. Normalized histogram of box ratios for each tilt angle ( ) for the cap reference frame design. PAGE 108 108 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%Frequency (%) 0.000.150.250.350.450.550.650.750.850.95 Box Ratio S Theta0 S Theta15 S Theta30 S Theta45 S Theta60 S Theta75 S Theta90 Figure 617. Normalized histogram of box ratios for each tilt angle ( ) for the half sphere test. Figure 618. Profiles of solid models. A) Caucasian. B) Asian. C) Infant. B A C PAGE 109 109 CHAPTER 7 CONCLUSION AND FUTURE WORK The results have determined a clear difference in stability of fit between the large masks and the small masks for all cases. Both the su rface plots and the box ratio histograms show that the large masks are stable while the small masks are not. Furthermore, the surface plots are able to show which force location points and what an gles result in an unstable fit. The box ratio histograms are able to quantify stability on a tilt angle basis; by assigning a normalized numerical value to the magnitude of the contac t volume, if three compression points are present that satisfy the TNT threshold. In addition, the box ratio histograms are able to provide information about the maximum tilt angle (theta) that the entire mask can withstand over all rotations about the main mask normal. The two additional cases, the caplike reference frame and the half sphere surface have further validated the program by providing expected results when compared to predicted realspace frame movements. Furthermore, these two additional surfaces have shown that the surround coverage of the head does not particularly results in stabil ity, and hints that the facial landmarks are the best locations to uniquely f it a mask reference frame onto a patient. The training data was compared to the three models facial landmarks and motivating relations about stability and facial angles were found, which suggests that if masks are built for patients with small nasiolateral and nasiofrontal angles and larg e nasiofacial angles, their masks may be highly stable. In the future more masks should be tested w ith the parameters defined in the methodology. Also it could be worth while to measure the facial angles of more patients and run statistical analysis to compare the relations between faci al angles and stability as measured by the methodology defined in this investigation. Conc erning the program, a mesh ing filter could be PAGE 110 110 beneficial for regularly meshes the surfaces. Also the box ratio definition could be redefine for easy understanding, by simple calcula ting the compression patch volume (cvol) divided by the entire frame volume (fvol), this will result in values from zero to one, in which zero signifies instability and one signifies full stability. The present investigat ion has provided with a methodology fo r testing uniqueness of fit for the patientspecific reference frames that impr oves the new and promising technology to perform imageguided surgery. The methodology will provid e the clinician with prior information of frameonpatient stability, which was not available before and is essential for allowing the clinician to effectuate his surg ical plan with more confidence and accuracy. In addition the methodology was created using an open source imag e analysis software (ITK and VTK), and by doing so, it will be easily adapted to the alrea dy existing user interface (RPD designer also was created with ITK and VTK), which allows the su rgeon to create his/her own reference frames based on diagnostic images. By including the stab ility measuring module in the already existing graphical user interface (RPD designer), the new technology for imageguided surgery will be more robust and a step closer towards being in troduced to other medical communities. PAGE 111 111 APPENDIX A PRINCIPLE OF VIRTUAL WORK AND ITK FEM MODULE This section provides additional background in formation related to the finite element method (FEM) for two main reasons. The methodol ogy for testing the stab ility of fit of the patientspecific frames is based on FEM and the tools used for pe rforming the FEM lack documentation. This chapter first provides b ackground information on FEM, and then uses the direct stiffness method to explai n the basic steps to derive the element stiffness matrix for a 2dimensional truss element. It also comments on the principles of virtua l work because the ITK FEM module (ITK 1.8, Kitware Inc. Clifton Park, NY) follows them to arrive to one of the elements stiffness matrices used in this inves tigation. Finally it disc usses the ITK FEM module and its organization in an attempt to supplemen t for the currently available documentation. FEM Terminology for Truss Elements There is a specific technical terminology a ssociated with FEM. The following list [4,28,46] presents the main components involved in FEM as it relates to structural analysis and will aid in the understanding of the subsequent sections. Discretization: The process of subdividing a stru ctures into members or finite elements with the purpose of modeling, simulating, or a pproximating their mechan ical behavior under known forces. Elements: In the case of structure analysis, an element is a single member of a structure. Its contribution to the structure is defined by a forcedisplacement relationship and this relationship is used to analyze the element. Node: The purpose of the node is to connect el ements, it is assumed to be infinitesimally small in size and therefore frictionless. Depe nding on the structure or part simulated the nodes can behave differently. For example, a node will act as a frictionless hinge in 2dimensional truss elements, and as a frictionl ess balland socket in 3dimensional truss elements. Global and Local Coordinates: In FEM, usually two coordinate systems are used. One is referred as the global coordinate system (X, Y, and Z) and it describes the overall geometry and loaddeformation relation of the entire structure. The ot her is referred as a local coordinate system (x, y and z) and it descri bes the forcedisplacement relations for each PAGE 112 112 individual element. In the case of the truss el ement, the local coordinate systems origin is located at one end of the truss element and the lo cal xaxis is oriented along the axis of the element. Degrees of Freedom (DOF): Degrees of freedom represent node displacements, such as translations and rotations. When a structure is subjected to an external load, the degrees of freedom are the independent node displacements n ecessary to indicate the deformed shape of the structure. A 1dimensional truss elem ent has 1 DOF per node, a 2dimensional truss element has 2 DOF per node and a space truss has 3 DOF per node. For truss elements the nodes are assumed to be frictionless joints, so th e elements are not subjected to moments, and the rotations are zero. Boundary Conditions (BC): Boundary conditions are applied at nodes to fix the model in space or to set a prescribed di splacement onto the nodes. In a structural model, they represent fixed supports necessary to remove the rigid body motion of a structure. In addition, when the boundary conditions are se t correctly, the total number of unknowns must be equal to, or smaller than, the number of equation necessary to describe the model under known loading conditions [4]. External forces and internal forces: External forces are specified in global coordinates and internal forces are specified in local coor dinates. Local coordi nates are induced onto elements by external forces. For example, an external force is appl ied on a particular node on a structure, as a result the elements that co mpose the structure will deform and internal forces will be induced at each individual elements node. Material Properties: Material properties define elastic char acteristics of the elements. For simple truss elements they are the elastic modulus E (MPa), and the Poisson ratio [3]. Transformation Matrices ( T ): For simplification, the integration of the element is performed in local coordinates. A transformation matrix is need that allows the transition from local to global (and vice versa) coordinate s. It involves the direction co sines of the local coordinate system in the global coordinate system. Element Stiffness Matrices ( ke): The stiffness matrix of an element relates the local forces to the local nodes displacement. The element stiffness matrix is composed of stiffness coefficients. The stiffness coefficient denoted kij is explained as the fo rce in the i direction required to cause a displacement at node j wh ile all other node disp lacements are zero [4]. Principle of Virtual Work In the FEM the forcedisplacement relations can be derived from the workenergy principles; the main reason for this being th at by using workenergy principles the forcedisplacement relations can be conve niently applied to solid elemen ts. The principle for virtual work for rigid bodies states [4]. PAGE 113 113 If a rigid body, which is in equilibrium under a system of forces (and couples), is subjected to any small virtual rigidbody displacement, th e virtual work done by the external forces (and couples) is zero. A virtual displacement is simply a small imagin ary displacement. The principle of virtual work for rigid bodies implies that if the rigid bod y is in equilibrium, then the sum of forces and moments must add to zero, ther eby causing the virtual work, We, to be equal to zero, or We = 0. The principle for virtual work for deformable solids state [4]. If a deformable structure, which is in equili brium under a system of forces (and couples), is subjected to any small virtual displacement consistent with the support and continuity conditions of the structure, then the virtual external work done by th e real external force (and couples) acting thro ugh the virtual external displacements (and rotations) is equal to the virtual strain energy stored in the structure. The principle mentioned above implies that virtual external work is equal to virtual internal work, and since internal work is equal to the structure strain energy then We = U. In this investigation we are using only lin ear elements. FEM linear analysis follows two basic assumptions [4]. The part or model is simulated with a linear elas tic material, this means that the strainstress relations behave in a linear manner and follow Hookes Law, which states that the stress (MPa) is proportionally related to the strain (unitless), by the elastic modulus of the material denoted as E (MPa), or =E The model experiences deformations that are infinitesimally small; therefore the squares and higher powers of the elements deformed shapes are negligible in comparison to unity. Indeed, the equation of equilibr ium can be based on the non deformed geometry of the part. ITK FEM and Principles of Virtual Work The following section is based on information from the ITK libraries [24], Aslam et al. [4] and Felippa et al. [17, 18] and has the purpose of providing background information on how the ITK FEM module sets up the element stiffness matrix, because the module lacks explanations and clear documentations. PAGE 114 114 The ITK FEM module follows principles of virt ual work to generate the element stiffness matrix in local coordinates. The ITK FEM m odule requires 1) the definition of the shape functions for each element types, 2) the deriva tion of the shape functions (strain displacement matrix) and 3) values related to the numerical integration of the element stiffness matrix. To get a brief idea of how thes e components relate to the prin ciples of virtual work first consider a 3dimensional different ial element in equilibrium and subjected to loading (Figure A1). Following the principle of virtua l work for deformable solid bodies [4], the loading on the differential element results in real stresses ( ), which cause in virtual strains ( ) and internal forces that perform virtual work. zx yz xy z y x zx yz xy z y x (A1) Equation A1 shows the not ation for the stresses ( and ) and virtual strain that define the differential element under loading a nd in equilibrium. The total virt ual internal work (or virtual strain energy) stored in the element can be measured by integrating over the volume [4], as shown in Equation A2. dV WV zx zx yz yz xy xy z z y y x x e (A2) dV Wv T e (A3) Equation A2 and Equation A3 (compressed form ed) defines the principle of virtual work for deformable bodies expressed in terms of stresses ( ) and virtual strains ( ), they represent the virtual strain energy stored in the element. PAGE 115 115 Now the focus becomes defining the local element stiffness matrix ke, in terms of stresses and strain using principles of vi rtual work of deformable bodies. The first step is to define the node displacements in terms of stresses and strains. Displacement Functions The node displacement for a truss element can be represented by a displacement function with the form of a linear polynomial (Equation A4) [11]. The displacement function or displacement variation (x) needs to be defined solemnly along the local xaxis. x a a x u1 0) ( (A4) The roots of the polynomial, a0 and a1, are constants and must satisfy the boundary conditions, that is the displ acement conditions at the nodes [11]. From Figure 21, when x is equal to 0 (at the origin of the local coordina te system), the displacement variation function (x) is equal to u1. 0 1) 0 ( a u u (A5) In the same way, when x is equal to L, (x) is equal to u2. L a a u L u1 0 2) ( (A6) The displacement variations of the truss el ement in terms of node displacements, u1 and u2 can be calculated by substituting Equation A5 and A6 into A4 and rearranging to yield Equation A7 [11]. 2 11 ) ( u L x u L x x u (A7) Shape Functions The displacement function (Equation A7) can be alternatively expressed as shape functions. There are as many shape functions for an element as there are nodes, and they have PAGE 116 116 the characteristic of having a value of one at its corresponding node (N1 has a value of one at node 1) and a value of zero at the other node points [3]. For a truss element there are two shape functions denotes by N1 and N2. The roots of the polynomials in Equation A7 represent the shape functions (Equation A8). L x N L x N2 1 1 (A8) For convenience the displacement variation function ( (x)) can be expressed in terms of shape functions and node displacements (Equation A9). Equation A10, is the final form or the elements displacement function (o r interpolating function) and N is the elements shape function matrix, and u are the elements node displacements. 2 1 2 1u u N N ux (A9) = Nu (A10) Strain Displacement Equation Getting back to the aim of representing the node displacement in term of stresses and strain, Equation A11 shows the relationship between strain ( ) and displacement variation function ( (x)), for a truss element [11]. dx u dx (A11) From Equation A11, notice that a truss element only withstand strains along the local xaxis. x xu dx du = D (A12) PAGE 117 117 D is called the differential operator matrix [4]. At last the strain can be related to the node displacements u by substituting Equation A10 into Equation A12. Bu DNu (A13) In order to get the strain displacement matrix, B, the derivatives of the shape functions need to be calculated (Equation A13). Th e node displacements do not depend on the local x direction, so they are treated as a constant, while the differentiation can be applied to shape functions N and not to the node displacements u [4]. Following Hookes law, the stresses can also be related to the end displacements (Equation A14), as required by the principle of virtual work for deformable bodies [4]. Bu E (A14) So far both strains (Equation A13) and the st resses (Equation A14) are expressed in terms of node displacements. Local Element Stiffness Matrix All the necessary informati on required by the ITK FEM module is now available to determine the elements stiffness matrix, with the exception of a relation between the element end forces (f1 and f2 as in Figure 31) and the node displacements u. When subjecting a truss element to the pr inciples of virtual external work ( We), the truss element will be in equilibrium and subjected to end forces. These end forces will manifest themselves as virtual end displacements ( u1 and u2). Equation A15 shows how the virtual external work can be defined in terms of virtual end displacements and end forces [4]. 2 2 1 1u f u f eW (A15) f uT eW (A16) PAGE 118 118 In symbolic form Equation A16 uT is the transform of the virtual end displacements and f is the end force vector. Equation A16 can be substituted into Equation A3 to get Equation A17. dVv T T f u (A17) A relation between the element stiffness matrix ke and the strain displacement matrix B is finally accomplished by subs tituting Equations A11 and A12 into Equation A17 and simplifying by canceling uT (Equation A18) [4]. u k u B B fe dV Ev T (A18) The element stiffness matrix ke is determined by virtual wo rk principles as defined by Equation A19 [4]. dV Ev T B Bek (A19) Gauss Quadrature To evaluate the integration in Equation A18, the ITK FEM module uses a method of numerical integration called the Gauss quadrature (Equation A20) [18]. i T i i n i i eE w B B J1k (A20) In this equation, n is the number of integrati on points, i is the integration point index, w is the integration weight, B is the strain displacement matrix, and Ji is the jacobian determinant evaluated at the integration point. In the ITK FEM module all these parameters are necessary to define an element [18]. The determination of the element stiffness matr ix is the fundamental part of the ITK FEM process; afterwards the element stiffness matrix needs to be transformed into global coordinates PAGE 119 119 and assembly into the stiffness matrix that defines the complete system [28]. In the ITK FEM module after the element stiffness matrix is cons tructed for a particular element and transformed into the global coordinates, the remaining steps of assembly and reduction are operations that are not element dependent. The purpose of the past sec tion was to explain the ba sic of parameters in referenced to the ITK FEM module to comp lement for the lack of documentation. ITK FEM Module ITK is open source software that provides data representation and algorithms for implementing image analysis, primarily segmentation and registration [24]. Its origins date back to 1999 when the US National Library of Medicine of the National Institute of Health granted a threeyear contract to create an opensource regi stration and segmentation t oolkit, later it came to be know as the Insight Toolkit or ITK. The software is wri tten in C++ and is greatly object oriented and its organization repr esents a pipeline where data and process are connected together. ITK has an FEM module primarily intended to be used for image registration, but complete to perform mechanical modeling. ITK provides an objectoriented FEM library to solve general FEM modeling. Classes are used to specify the geometry, behavior of the elements, apply external forces, apply boundary conditions, specify elasticity, and to solve the prob lems (Figure A2). The core of the library is the element class, which is divided in two parts: the definition of the el ement geometry, and the physics of the elements or the behavior of the el ement. The outcome of combining the two parts results in the creation of a new element. Fi gure A2 shows a section of the ITK FEM element class. Notice that from the Element class (i tk::fem::Element) the space truss is defined by itk::fem::ElementStd<2,3> and the linear tetrah edron by itk::fem::ElementStd<4,3>. The first number in the class name implies the number of nodes in the element and the second number implies the dimension in which the element resides. PAGE 120 120 Following Figure A2B, itk::fem::Element3DCO LinearLine is a class derived from ElementStd (ElementBaseClass), this class defi nes the geometry of the element by defining the following parameters. These parameters relate to the following the FEM by principles of virtual work. Integration Point and weight Number of Integration Points Shape Functions Shape Function Derivatives Get Local From Global coordinates Jacobian In the class itk::fem::Element3DStrain PAGE 121 121 itk::fem::LoadBC: Class that applies essent ial or dirichlet boundary conditions, specifies which degrees of freedom are fixed in the model. itk::fem::LoadNode: Class used to specify a load on a point within an element object. A pointer to an element object, the number of th e node on which the load acts, and the force vector. itk::fem::Material: Base class that defi nes the elements material properties, itk::fem::MaterialLinearElasticity: Class that defines the material properties of the finite elements, these are: E (elastic modulus), A(area ), I (moment of interia) nu (poisson ratio), h (thickness), and Rhoc (density multiplied by heat capacity). Figure A1. Differential element subjected to real stresses and in equilibrium. y z x zx yz xy zx xy yz y x z PAGE 122 122 Figure A2. The greater part of the ITK FEM library. The boxes surrounded by dashed lines represent where the geometry of the el ement is defined. The boxes surrounded by dotted line, represent where th e physics of the problem is defined. Together they entirely define the element, in this cases the Element3DC0LinearLineStrain and the Element3DC0LinearTetrahedronStrain [24]. itk::fem::FEMLightObject itk::fem::Element itk::fem::Element::Node itk::fem::Load itk::fem::Material itk::fem::ElementStd<2,2> itk::fem::ElementStd<3,2> itk::fem::ElementStd<4,2> itk::fem::ElementStd<2,3> itk::fem::ElementStd<4,3> itk::fem::ElementStd<6,2> itk::fem::ElementStd<8,3> itk::fem::LoadBC itk::fem::LoadBCMF itk::fem::LoadElement itk::fem::LoadNode itk::fem::MaterialLineaerElasticity itk::fem::Element3DCOLinearTetrahedron itk::fem::Element3DStrain PAGE 123 123 LIST OF REFERENCES 1. Adler JR. Surgical guidance now and in the nuture. The ne xt generation of instrumentation. Clinical neurosurgery. Proceedings of the Congress of Neurological Surgeons. San Diego, California 2001. Ed Lippincott Williams & Wilkins 2002;105114. 2. Alker G, Kelly PJ. An overview of CT base d stereotactic systems for localization of intercranial lesions. Comput ational Radiology 1984;8:193196. 3. Allaire PE. Basics of the finite element me thod. Iowa: Wm. C. Brown Publishers, 1985. 4. Aslam, Kassimali. Matrix analysis of st ructures. California: Brooks/Cole Publishing Company, 1999. 5. Becker AA. An introductory guide to finite element analysis. New York: ASME Press, 2004. 6. Birnbaum Klaus, S.E., Decker Nils, Presch er Andreas, Ulrich Klapper, Radermarcher Klaus. ComputerAssisted Orthopedic Surger y with individual templates and comparison to Conventional Operation Method. Spine 2001;26:365370. 7. Bova, F. Grant proposal: Imageguided su rgery based on rapid prototyping models. Department of Neurological Surgery, Universi ty of Florida. Gainesville, FL. March 2003 8. Brown RA: A computer tomographycomput er graphic approach to stereotactic localization. J Neurosurg 1979; 50:715720. 9. Brown RA: .A stereotactic h ead from for use with CT body scanners. Invest Radiol 1979;14:300304. 10. Cheney ML. Facial surgery. Plastic and reconstructive. Baltimore: Williams & Williams, 1997. 11. Cheung YK, Yeo MF. A practical introduction to finite element analysis. London: Pitman Publishing Limited, 1979. 12. Cowin, S. C. Bone Mechanics. New York: CRC Press, 2001. 13. Davis, T A. User Guide for Cholmod: A Sp arse Cholesky Factori zation and Modification Package. Version 1.4. Gainesville: Departme nt of Computer and Information Science and Engineering, 2006. 14. DUrso PS, Atkinson RL, Lanigan MW, Earwak er WJ, Bruce IJ, Holmes A, Barker TM, Effeney DJ, Thompson RG. Stereolithographic (S L) biomodelling in craniofacial surgery. British Journal of Plastic Surgery 1998;51:522530. PAGE 124 124 15. DUrso PS, Earwaker WJ, Barker TM, Redmond MJ, Thompson RG, Effeney DJ, Tomlinson FH. Custom cranioplasty using stereolithography and acr ylic. British Journal of Plastic Surgery 2000;53:200204. 16. DUrso US Patent # 5,752,962. 1998. 17. Felippa Carlos. Introduction to finite el ement methods 2006. Department of Aerospace Engineering Sciences, University of Colo rado at Boulder. Last Accessed 30 Jan 2006. http://www.colorado.edu/engine ering/CAS/courses.d/IFEM.d/ 18. Felippa Carlos. Nonlinear Finite Elem ent Methods 2006. Department of Aerospace Engineering Sciences, University of Colo rado at Boulder. Last Accessed 30 Jan 2006. http://www.colorado.edu/engine ering/CAS/courses.d/NFEM.d/ 19. Gildenberg PL. Sterotactic Surgery: Present and Past. In Heilbrun MP, ed. Stereotactic Neurosurgery Vol 2: Concepts in Neurosur gery. Balitmore:Williams & Williams, 1988.:116. 20. Goffin J., V.B.K., Vander Sloten J., Van A uderkercke R., Smet M. H., Marchal G., Van Craen W., Swaelens B., Verstreken K. 3DCT ba sed, personalized drill guide for posterior transarticular screw fixation at C1C2: Technical Note. NeuroOthopedics1999;25:4756. 21. Goffin Jan, Van Brussel Karel, Martens Kr isten, Vander Sloten Jos, Van Audekercke Remi, Smet MariaHelena. Threedimensiona l computer tomographybased, personalized drill guide for posterior cer vical stabilization at C1C 2. Spine 2001,21;12:13431347. 22. Heilbrun MP. Computed tomographyguide d stereotactic system. Clinical Neurosurgery1984;31:564581. 23. Horsley V, Clarke RH: The structure and function of the cerebellum examined by a new method. Brain1908;31:45124. 24. Ibanez L., Schroeder. W., Ng L, Cates J. The ITK Software Guide.NewYork: Kitware Inc.,2003. 25. Loftis KL, Geer CP, Daneslson KA, Slice DE Stitze JD. Analysis of pediatric head anthropometry using computed tomography for application to head injury prediction. 2007;43:8690. 26. Mackerle J. Finite element modeling and si mulations in orthopedics : a bibliography 19982005. Comput Methods Biomech Biomed Engin2006;9(3):149199. 27. Maciunas RJ, Galloway RL, Latimer JW. The Application Accuracy of Stereotactic Frames. Neurosurgery1994;35(4):682695. 28. Kim NH., Sankar B V. Lecture notes: EML 4500. Finite Elements analysis and design. Gainesville: Department of Mechan ical & Aerospace Engineering, 2004. PAGE 125 125 29. Pham DT, Gault RS. A comparison of rapi d prototyping technologi es. International Journal of Machine tools a nd Manufacturing1998;38:12571287. 30. Pikley WD, Pikley OH. Mechanics of Solids. New York: Schaum, 1974. 31. Popov, E P. Introduction to mechanics of solids. New Jersey: Prentice Hall,1968. 32. Raabe A, Krishnan R, Zimmermann M, Seif ert V: Frameless and Framebased stereotaxy? How to choose the appropriate procedure. Zent ralbl Neurochir.2003;64(1): 15. 33. Rajon DA, Bova FJ, Bhasin RR, Friedman WA. An investigati on of the potential of rapid prototyping technology for imageguided surgery. Journal of Applied Clinical Medical Physics 2006;7(4):8198. 34. Rajon DA, Bova FJ, Bhasin RR, Friedman WA. Poster Pr esentation: Fabrication of custom surgical guides using rapidprototypi ng technology. University of Florida College of Medicine. Celebration of Research; March 20,2007. 35. Saw CB, Ayyangar K, Suntharalin gram N: Coordinate transf ormation and calculation of the angular paramters for stereot actic system. Med Phy1987;14:10421044,. 36. Schad L. Lott S, Schmitt F, Sturn V, Lorenz WJ: Correction of spatial distortion in MR imaging: a prerequisite for accurate stereota xy. Journal Computer Assisted Tomography 11(3):499505,1987. 37. Schileo E, Taddei F, Malandrino A, Cristofolin i A, Viceconti M. Subjectspecific finite element models can accurately predict stra in levels in long bo nes. Journal of Biomechanics. 2007 doi: 10.1016/j.biomech.2007.02.010 38. Schroeder W, Lorensen B. The visualizati on Toolkit. New Jersey : Prentice Hall, 1997. 39. Sedef M, Samur E, Basdogan C. Realtime fin iteelement simulation of linear viscoelastic tissue behavior based on experimental data. IEEE comput Graph Appl 2006;26(6):5868. 40. Smith KR, Frank KJ, Bucholz RD. The NeuroStationa highly accurate, minimally invasive solution to frameless stereotactic ne urosurgery. Computerized Medical Imagining and Graphics 1994;18(4):247256. 41. Shewchuck JR. What is a good linear elemen t? Eleventh International Meshing Round Table. New York: Sandia Nationa l Laboratories September 2002:115126. 42. Spiegel EA, Wycis HT. Stereoen cephalotomy, Part I. New York: Grune & Shatton, 1952. 43. Spiegel EA, Wycis HT, Lee MM. A: Stereota ctic apparatus for operations on the human brain. Science1947;106:349350. 44. Spiegel EA, Wycis HT: Thalamotomy and related procedures. JAMA 1952;148:446451. PAGE 126 126 45. Suresh G, Elliot N, Pinson N, Sinden FW. Simulation of dynamics of interacting rigid bodies including friction I: General probl em and contact model. Engineering and Computers 1994;10:162174. 46. Timoshenka. S P, Gere J M. Mechanics of Materials. California: Thomson Brooks/Cole Engineering Division, 1984. 47. Van Cleynenbreugel J., S. F., Goffin J., Van Brussel K., Suetens P.. Imagebased planning and validation of C1C2 transa rticular screw fixation using personalized drilled guides. Computer Aided Surgery. 2002;7:4148. 48. Van Staden RC, Guan H, Loo YC. Applicati on of the finite element method in dental implant research. Comput Methods Bi omech Biomed Engin. 2006;9(4):257270. 49. Yan X, Gu P. A review of rapid prototypi ng technologies and system s. Computer Aided Design1996;28(4):307318. 50. Yoo T. S., M.J., Chen D. T., Richardson A. C., Burgess J., Template Guided Intervention: Interactive Visualization and Design for Me dical Fused Deposition Models. Proceedings of the Workshop on Interactive Medical Imag e Visualization and Image Analysis 2001;18:4548. 51. Zamorano LJ, NolteL, Kadi AM, Jiang Z. Inte ractive intraoperative localization using an infraredbased system. Stereotactic and Functiona l Neurosurgery 1994;63:8488 PAGE 127 127 BIOGRAPHICAL SKETCH Barbara Garita Canet was born in 1977. She is the second child of Donato Garita and Noemi Canet who had six children. She grew up in San Jose, Costa Ri ca and attended Lincoln School from 1983 to 1995. After graduation from Lincoln School, she started her career in mechanical engineering at the Un iversidad de Costa Rica; however she only studied there for a semester. She was offered a scholarship for pl aying varsity volleyball for Florida Institute of Technology (FIT). From FIT she graduated with honors with a bachelors degree in mechanical engineering in May 2000. In August 2000 she starte d a masters in biomed ical engineering at the University of Florida and graduated in Decem ber 2002. The topic of her masters thesis was related to bone mechanics. From December 2002 to January 2003 she looked for a new committee chair to work on her doctorate. Finally in January 2004 she stated her doctorate, on mechanical stability for customized instrume nts in imageguided surgery, and in August 2007 she completed her doctoral studies. Immediately af ter graduation she will rest and travel. She will then look for a job in biomedical engineer ing field and will write a business plan with the goal to start her own company. 