Citation |

- Permanent Link:
- https://ufdc.ufl.edu/UFE0012401/00001
## Material Information- Title:
- Guided Surgery Using Rapid Prototyping Patient-Specific Guides A Methodology to Quantify Mechanical Stability and Uniqueness of Fit
- Creator:
- Garita, Barbara
- Place of Publication:
- [Gainesville, Fla.]
Florida - Publisher:
- University of Florida
- Publication Date:
- 2007
- Language:
- english
- Physical Description:
- 1 online resource (127 p.)
## Thesis/Dissertation Information- Degree:
- Doctorate ( Ph.D.)
- Degree Grantor:
- University of Florida
- Degree Disciplines:
- Biomedical Engineering
- Committee Chair:
- Bova, Frank J.
- Committee Members:
- Friedman, William A.
Gilland, David R. Vemuri, Baba C. Kim, Nam Ho - Graduation Date:
- 8/11/2007
## Subjects- Subjects / Keywords:
- Asians ( jstor )
Boxes ( jstor ) Coordinate systems ( jstor ) Histograms ( jstor ) Infants ( jstor ) Linear programming ( jstor ) Stiffness ( jstor ) Stiffness matrix ( jstor ) Three dimensional modeling ( jstor ) Trusses ( jstor ) Biomedical Engineering -- Dissertations, Academic -- UF contact, fem, guided, image, itk, mechanics, prototyping, rapid, surgery, vtk City of Gainesville ( local ) - Genre:
- bibliography ( marcgt )
theses ( marcgt ) government publication (state, provincial, terriorial, dependent) ( marcgt ) born-digital ( sobekcm ) Electronic Thesis or Dissertation Biomedical Engineering thesis, Ph.D.
## Notes- Abstract:
- Rapid prototyping technology (3-dimensional printing) in conjunction with the virtual manipulation of computer-generated 3-D anatomical models rendered from diagnostic images (CT and MRI) is a novel technique for directing image-guided surgery. By incorporating a graphical user interface, the surgeon can plan an intracranial surgical procedure by fabricating a patient-specific reference frame (mask-like facial frame) through the use of rapid prototyping technology. Early testing of these patient-specific 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 finite-element 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 patient-specific 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 3-D 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 cap-like 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. ( en )
- General Note:
- In the series University of Florida Digital Collections.
- General Note:
- Includes vita.
- Bibliography:
- Includes bibliographical references.
- Source of Description:
- Description based on online resource; title from PDF title page.
- Source of Description:
- This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
- Thesis:
- Thesis (Ph.D.)--University of Florida, 2007.
- Local:
- Adviser: Bova, Frank J.
- Statement of Responsibility:
- by Barbara Garita.
## Record Information- Source Institution:
- UFRGP
- Rights Management:
- Copyright Garita, Barbara. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
- Resource Identifier:
- 660162369 ( OCLC )
- Classification:
- LD1780 2007 ( lcc )
## UFDC Membership |

Downloads |

## This item has the following downloads:
garita_b.pdf
garita_b_Page_112.txt garita_b_Page_033.txt garita_b_Page_089.txt garita_b_Page_069.txt garita_b_Page_065.txt garita_b_Page_086.txt garita_b_Page_012.txt garita_b_Page_096.txt garita_b_Page_068.txt garita_b_Page_017.txt garita_b_Page_035.txt garita_b_Page_052.txt garita_b_Page_107.txt garita_b_Page_084.txt garita_b_Page_055.txt garita_b_Page_081.txt garita_b_Page_123.txt garita_b_Page_095.txt garita_b_Page_085.txt garita_b_Page_099.txt garita_b_Page_057.txt garita_b_Page_058.txt garita_b_Page_018.txt garita_b_Page_125.txt garita_b_Page_059.txt garita_b_Page_115.txt garita_b_Page_121.txt garita_b_Page_031.txt garita_b_Page_048.txt garita_b_Page_070.txt garita_b_Page_117.txt garita_b_Page_083.txt garita_b_Page_053.txt garita_b_Page_036.txt garita_b_Page_009.txt garita_b_Page_067.txt garita_b_Page_022.txt garita_b_Page_002.txt garita_b_Page_015.txt garita_b_Page_029.txt garita_b_Page_016.txt garita_b_Page_098.txt garita_b_Page_043.txt garita_b_Page_118.txt garita_b_Page_040.txt garita_b_Page_103.txt garita_b_Page_082.txt garita_b_Page_122.txt garita_b_Page_003.txt garita_b_Page_030.txt garita_b_Page_032.txt garita_b_Page_087.txt garita_b_Page_075.txt garita_b_Page_097.txt garita_b_Page_051.txt garita_b_Page_037.txt garita_b_Page_073.txt garita_b_Page_025.txt garita_b_Page_126.txt garita_b_Page_013.txt garita_b_Page_108.txt garita_b_Page_041.txt garita_b_Page_094.txt garita_b_Page_056.txt garita_b_Page_110.txt garita_b_Page_063.txt garita_b_Page_062.txt garita_b_Page_080.txt garita_b_Page_079.txt garita_b_Page_042.txt garita_b_Page_113.txt garita_b_Page_114.txt garita_b_Page_038.txt garita_b_Page_007.txt garita_b_Page_120.txt garita_b_Page_014.txt garita_b_Page_054.txt garita_b_Page_024.txt garita_b_Page_020.txt garita_b_Page_008.txt garita_b_Page_034.txt garita_b_Page_109.txt garita_b_Page_091.txt garita_b_Page_039.txt garita_b_Page_061.txt garita_b_Page_021.txt garita_b_Page_045.txt garita_b_Page_004.txt garita_b_Page_127.txt garita_b_Page_093.txt garita_b_Page_026.txt garita_b_Page_077.txt garita_b_Page_006.txt garita_b_Page_005.txt garita_b_Page_076.txt garita_b_Page_044.txt garita_b_Page_028.txt garita_b_Page_060.txt garita_b_Page_019.txt garita_b_Page_111.txt garita_b_Page_116.txt garita_b_Page_047.txt garita_b_Page_001.txt garita_b_Page_104.txt garita_b_pdf.txt garita_b_Page_072.txt garita_b_Page_090.txt garita_b_Page_050.txt garita_b_Page_023.txt garita_b_Page_101.txt garita_b_Page_100.txt garita_b_Page_049.txt garita_b_Page_027.txt garita_b_Page_078.txt garita_b_Page_011.txt garita_b_Page_105.txt garita_b_Page_092.txt garita_b_Page_064.txt garita_b_Page_102.txt garita_b_Page_106.txt garita_b_Page_071.txt garita_b_Page_119.txt garita_b_Page_010.txt garita_b_Page_074.txt garita_b_Page_124.txt garita_b_Page_066.txt garita_b_Page_046.txt garita_b_Page_088.txt |

Full Text |

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 5-1. 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 | 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 4-16. 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 4-17. 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' TNT<=0.70 or n/t = -0.41 Figure 4-12. Two-dimensional 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 4-13A to correspond with spring 1 (with the tangential displacement in the opposite direction), and spring 2 corresponds to the spring in Figure 4-13B. n=l S0.833 TNT = 0.5 Figure 4-13. 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. based image-guided 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 3-dimensional model from diagnostic image datasets. This program allows selecting a surface on the virtual 3-dimensional 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 frame-on-patient 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 patient-specific 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 2-10 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. 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 ball-and-socket joint. External force magnitude The external force magnitude is calculated by Equation 4-3 and selecting an expected axial displacement (AL) of 1% the spring's length. As shown by Equations 3-2 and 4-2, 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 4-3. Area *E AverageForce spring AL(43) spring Lo Frontal enilen&. .---- Frontal Superciliary arch Glabella I / Nasion Radix Cephalic C ( audal Medial-Lateral Cheek (z omar;ti mcI') ,, Figure 2-10. Facial landmarks. 6-17 Normalized histogram of box ratios for each tilt angle (0) for the half sphere test........ 108 6-18 P profiles of solid m models ......................................................................... ....................108 A-1 Differential element subjected to real stresses and in equilibrium. ...............................121 A-2 The greater part of the ITK FEM library ................................................................... 122 .6'- ^,^ tt- ! ItV Theta *90 750 600 450 300 150 m0o NiONv a't 7 Ji__ Figure 6-7. 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 6-8. Surface plot color coded by maximum tilt angle for the cap reference frame design B 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 2-1. Pre-target localization and pre-surgical planning is done in a computer workstation. 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 2-12). 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 2-12). 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 ball-and- socket joints. For example, Figure 2-13 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 2-13). The space truss is considered a simple element since it deforms only along one direction (along its local x-axis) [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 2-12). 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 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 6-1. 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 6-2. 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% Cap-like frame 2695 241 9% Half-sphere test 1094 975 94% tested (through a non-significant 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 1-1D and 1-1E). 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 1-1D). For this reason, the accuracy of the patient-specific 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 patient-specific 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 patient-specific frame and To Donato, Noemi, Esteban, Pablo, Ifi, Ignacio and Natalia. 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 6-13. 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 6-14. Normalized histogram of box ratios for each tilt angle (0) for the infant large mask. 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 patient-specific reference frame that has a potential of becoming a new type of image-guided surgery (IGS) technique [33]. The fabrication of the patient-specific reference frame involves coding a series of algorithms that render a detailed 3-dimensional life-like anatomical model from computer tomography (CT) or magnetic resonance imaging (MRI) [6] datasets (Figure 1-1). Onto a patient-specific 3-dimensional anatomical model a surface is selected (virtually "painted") from which the patient-specific reference frames are created. The ability to create, visualize, and manipulate 3-dimensional 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 patient-specific reference frame is to replace existing technology of framed-based (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 patient-specific 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 image-guided 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 Figure 3-1. The three models used to evaluate the methodology for testing unique fit. A) Caucasian, B) Asian and C) infant. 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 patient-specific reference frames that improves the new and promising technology to perform image-guided surgery. The methodology will provide the clinician with prior information of frame-on-patient 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 image-guided surgery will be more robust and a step closer towards being introduced to other medical communities. The element stiffness matrix (Equation 2-4) 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 2-5. T=coso coso0 cos0 0 0 0 T = 0 0(2-5) 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 2-6 (T means the transform of the matrix). K = TTkeT (2-6) For the space truss element, the element's stiffness matrix in global coordinates Ke results in a 6 x 6 matrix (Equation 2-7). 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 (2-7) 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 2-3); then to calculate the transformation matrix T (Equation 2-5) necessary to swap from local to global coordinates; and finally it is important to define the element stiffness matrix in global coordinates Ke (Equation 2-7). 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 CHAPTER 4 METHODOLOGY FOR MEASURING STABILITY Introduction The methodology for testing stability of fit of the patient-specific 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 patient-specific 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 4-1), 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 4-2). 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 life-like facial models that represent three different patient's anatomies were constructed with 3-dimensional 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 / --. JId II / Figure 4-9. 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'. Force application point Force direction 1350 B Figure 5-5. 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. element is given by Equation 2-1, and the local node displacements ul and U2 are caused by the local node forces fi and f2, as seen in Figure 2-14. u= (2-1) 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 2-2, where k is the stiffness coefficient and is equal to AE/L. f = ku (2-2) With Equation 2-2 and the principles of equilibrium [28], the element stiffness matrix (ke) for a particular element e can be calculated. As Figure 2-14A 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 2-2). 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 2-14B) 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 2-2. S-k -k =AE[ 1 1 -1(2-3) k = or ke = AE (2-3) -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 2-4). f}=AEI -]{i1 { or }= ke L (2-4) f2 L -1 1 1U2 fJ2 2 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 6-15. 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 6-16. Normalized histogram of box ratios for each tilt angle (0) for the cap reference frame design. points (or nodes). On the whole, the rigid layer represents the patient-specific 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 4-4). 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 4-5). 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 4-5, 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 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 4-18). 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 4-19). 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 4-18 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 4-7). Box ratio = cvo (4-7) fvol As Equation 4-7 shows, a box ratio is equal to 0 when the compression patch volume is equal to zero (or the mask is unstable) (Figure 4-19). 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 4-19). Box ratios closer to the value of zero correspond to smaller compression patch volumes or "good compressor" points close to each other. * 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 A-1. Differential element subjected to real stresses and in equilibrium. 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 --------------------------------------- ---- Outer extrapolated Surface Selected Rigid layer Surface Bed of springs / Boundary Conditions B Figure 4-1. 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 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 3-dimensional 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 2-node truss element in 3-dimensional space called space truss. The second one is the linear 4-node solid tetrahedral element in 2-dimensional 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 3-dimensional structures that are subjected to 3-dimensional 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 ball-and-socket joints [17]. In a structure composed of space trusses, any Figure 2-4. 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 2-5. 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. 1 2 3 5 6 --F Figure 4-10. 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 4-11. 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~~ 3-7 The method employed to demonstrate the difference in stability between the large and sm all m asks. ...........................................................................53 4-1 The methodology for stability of fit involves simulating a rigid layer onto of a deform able bed of springs .......................................................... .. ............... 71 4-2 Displacement of a mask as it responds to external forces. ..............................................72 4-3 Triangulation of a customized positioning device proper to STL format........................72 4-4 D direction of the N main ...................................... ......... ....................... ......... 73 4-5 N odes point on the rigid fram e ................................................ .............................. 73 4-6 In a wedge three tetrahedral elem ents are fitted ..................................... .................74 4-7 Schematic showing how the maximum rigidity fault is calculated .................................75 4-8 The behavior of a linear spring element under the application of an axial force ..............75 4-9 Schematic that demonstrates how each spring's displacement was evaluated by identifying the normal (ni) and tangential (ti) displacement magnitudes...........................76 4-10 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 4-11 Bounding volume for a spring to be in compression ......................................................77 4-12 Two-dimensional interpretation of the single point criteria ............................................78 4-13 Schematic to demonstrate the importance of evaluating the percentage slip of TNT ratio for each spring ......................... ......... .. ........... .. ............78 4-14 Schematic to shows at the level at which the TNT threshold was defined......................79 4-15 Empirical determination that the TNT threshold greater than 0.6 correlates with the training data, which says that the large masks are stable .............................................79 4-16 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 4-17 A case of the large Asian mask, which shows the consistency of an upper TNT thresh old of 0 .7 ............................................................................ 80 4-18 Schematic showing the relation between the volume of the entire frame (fvol) and the com pression patch volume e (cvol) ......................................................... ............... 81 4-19 B ox ratio interpretation. ......................... ....................... .......... .... .... 1 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 Brown-Roberts-Wells (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. Frame-Base and Frameless Stereotactic Approaches Today image guidance procedures can be divided in two categories, according to different referencing methods employed. These are generally termed frame-base 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 3-dimensional virtual model), onto a patient's actual, real-world 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 real-world patient and correlates the information amongst them. The creation of a diagnostic plan, for both frame-base 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 frame-based or skin landmarks for frameless stereotactic surgeries (Figure 2-1) [6]. less than or equal to one (Figure 4-19). 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 6-10) 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 6-5B. 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 6-13), 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 6-15). 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 6-10). 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 6-12) and the infant large mask (Figure 6- 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 x-axis 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 1-dimensional truss element has 1 DOF per node, a 2-dimensional 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 force-displacement relations can be derived from the work-energy principles; the main reason for this being that by using work-energy principles the force- displacement relations can be conveniently applied to solid elements. The principle for virtual work for rigid bodies states [4]. location where it was designed to fit. The color coded surface plot results for the cap-like 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 6-8). The box ratio histogram (Figure 6-16) 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 half-sphere surface (Figure 6-9), 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 6-2 shows the instability scores for all masks and frames as defined in Equation 4-8. There is zero percent instability for all large masks and there is some percentage of instability for all small masks, the cap-like frame and the half-sphere surface, as was expected. The least stable case is the half-sphere 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 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. f, =k*1 fi -- z A z n2 x u2= I. f2 = Sn2 f2 Figure 2-14. 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 2-15. The linear tetrahedron has a unique numbering scheme. k*1 L ,,I -- Tetrahedral 3 3 0 Tetrahedral 2 Figure 4-6. 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 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 6-11. 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 6-12. Normalized histogram of box ratios for each tilt angle (0) for the Asian large mask. CHAPTER 2 BACKGROUND AND SIGNIFICANCE Image-Guided 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 3-dimensional 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 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 4-1). The rigidity of the frame was quantified with the RDR ratio and Table 4-1 shows that reducing the height of the rigid layer considerably affects the rigidity of the rigid frame. For example, in Table 4-1, 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 le-13. A difference of six orders of magnitude for the Er and Es maintains the rigidity and results in a reliable simulation. Table 4-2 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 4-9). The displacement vector for each individual spring is characterized by two components, the normal displacement 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 patient-specific 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 patient-specific 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. 5-1 Flow chart of program for testing the uniqueness of fit.............................................88 5-2 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 5-3 Examples of al linear spring element or space truss tested in all 3-dimensions ...............90 5-4 Example of two tetrahedral elements used in the validation of the tetrahedral elem ents used in the program ................................................. ............................... 91 5-5 A corer frame that was used to validate the program ............................................... 92 6-1 Placem ent of Caucasian m asks onto solid m odel ........................................ .................100 6-2 Placement of Asian masks onto solid model .............. ..........................................100 6-3 Placement of infant masks onto solid model .............. ..........................................101 6-4 Placement of a cap-like reference frame onto an upper head model............................101 6-5 Surface plot color coded by the maximum full rotation angle for the Caucasian m a sk s ................... ........................................................... ................ 10 2 6-6 Surface plot color coded by the maximum full rotation angle for the Asian masks........102 6-7 Surface plot color coded by the maximum tilt angle for the infant masks ....................103 6-8 Surface plot color coded by maximum tilt angle for the cap reference frame design .....103 6-9 Surface plots color coded by maximum tilt angle for the half sphere test.......................104 6-10 Normalized histogram of box ratios for each tilt angle (0) for the Caucasian large mask..................... ...................................... 104 6-11 Normalized histogram of box ratios for each tilt angle (0) for the Caucasian small mask..................... ...................................... 105 6-12 Normalized histogram of box ratios for each tilt angle (0) for the Asian large mask. ....105 6-13 Normalized histogram of box ratios for each tilt angle (0) for the Asian small mask.....106 6-14 Normalized histogram of box ratios for each tilt angle (0) for the infant large mask. ....106 6-15 Normalized histogram of box ratios for each tilt angle (0) for the infant small mask.....107 6-16 Normalized histogram of box ratios for each tilt angle (0) for the cap reference frame d e sig n ........ ........ .......................................................................... 10 7 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 A-7 represent the shape functions (Equation A-8). N, = ) N2 =Q (A-8) For convenience the displacement variation function (i(x)) can be expressed in terms of shape functions and node displacements (Equation A-9). Equation A-10, 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. [i-x]= [N, N2 1 (A-9) LU2 i Nu (A-10) Strain Displacement Equation Getting back to the aim of representing the node displacement in term of stresses and strain, Equation A-11 shows the relationship between strain (E) and displacement variation function (u(x)), for a truss element [11]. E = (A-11) dx From Equation A- 1, notice that a truss element only withstand strains along the local x- axis. S= [d [ ] = Do (A-12) dx or MR datasets [39]. In both cases, the procedures do not have to be completed in an extended and continuous span of time. Both frame-based and frameless procedures permit the same virtual computer planning to be performed by the surgeon; however frameless stereotactic surgery requires instrument-tracking 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 2-4). Although electromagnetic tracking allows for non-line 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 C1-C2 (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 2-5 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 75N-2mm d= =- = 0.1 mm (5-4) 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 5-5. Both ANSYS and ITK FEM provided the same displacement shown in Table 5-1. A Simple Model Validation The program was validated with a simple frame in the shape of a corner (Figure 5-5). 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 5-5A). 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. displacements in global coordinates (Equation 2-8). The solution to the system of linear equation is given by solving for the displacements in global coordinates, Ug. F, =KU, (2-8) The linear tetrahedron Solid elements are defined in 2-dimensions and can model a solid body. The linear tetrahedron is the simplest solid element in 2-dimensional 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 right-hand rule (Figure 2-14) and the calculation of the tetrahedron's volume must result in a positive value (Equation 2-9). 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 2-9, 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 (2-9) 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, multi-dimensional elements [4]. 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 4-8b shows an axial force Fy placing a spring under a state of compression and changing its length to a new length L1. Figure 4-8c shows the free-body diagram of the system and the spring internal force Fs, which is determined by the equation in Figure 4-8d. In the Finite Element Method, linear springs that reside in 3-dimensional space are modeled as trusses elements in 3-dimensional space. In order to define a linear spring in 3-dimesions 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 4-1 as Ai. In this manner, even though the triangular mesh in the STL file is irregular (Figure 4-3); 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 4-1 [46] a the elastic modulus (E) and the spring length (Lo) are kept constant. A E k = E (4-2) L Maximum rigidity fault Figure 4-7. Schematic showing how the maximum rigidity fault is calculated. y LI L1 x // /// Fs = -k(L, L,) D Fs Figure 4-8. 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. Element e -112 -X 1 V2w U2 LU2 w1 X Figure 2-12. Space truss element has six displacements to establish its deformed position. There are ul and u2 along its x-axis, vi and v2 along its y-axis and wl and w2 along its z- axis. Though the space truss is defined in 3-dimensional space, it will only develop axial forces, along its local x-axis. symbol Figure 2-13. 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. D is called the differential operator matrix [4]. At last the strain can be related to the node displacements u by substituting Equation A-10 into Equation A-12. S= DNu =Bu (A-13) In order to get the strain displacement matrix, B, the derivatives of the shape functions need to be calculated (Equation A-13). 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 A-14), as required by the principle of virtual work for deformable bodies [4]. a = EBu (A-14) So far both strains (Equation A-13) and the stresses (Equation A-14) 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 3-1) 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 A-15 shows how the virtual external work can be defined in terms of virtual end displacements and end forces [4]. W, = f,5ug + f2 u2 (A-15) W, = Su'f (A-16) -"v Off& _0'.% . Figure 2-2. Procedures related to frame-based 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 2-3. 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 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 4-6). Numerically, it is defined by Equation 4-6 as the tangential displacement magnitude (t) divided by the sum of the normal displacement magnitude (n) and the tangential displacement magnitude (t) (Figure 4-9). TNT = (4-6) 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 4-10 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 4-6. 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 4-10, 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 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 4-13B 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 4-13B 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 7-1, 7-2, and 7-3). 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 4-15 and 4-16 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 4-15 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 Figure 2-7. 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]. 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 right-hand-rule as explained in Figure 2-15. 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 4-7). As Equation 4-1 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 (4-1) 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 Figure 6-1. Placement of Caucasian masks onto solid model. A) Large Mask. B) Small Mask. Figure 6-2. Placement of Asian masks onto solid model. A) Large Mask. B) Small Mask. Table 4-1. 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 4-2. 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. Er-Es Mean RDR Cholmod's rcond 1.0e6 3/1000 2.14e-8 1.0e7 2/1000 2.14e-9 1.0e8 5/1000 2.14e-10 Table 4-3. 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 Figure 4-2. Displacement of a mask as it responds to external forces. Figure 4-3. Triangulation of a customized positioning device proper to STL format. Every vertex in the mask is a point coordinate in 3D space. 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:200-204. 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. 3D-CT based, personalized drill guide for posterior transarticular screw fixation at C1-C2: Technical Note. Neuro-Othopedicsl999;25:47-56. 21. Goffin Jan, Van Brussel Karel, Martens Kristen, Vander Sloten Jos, Van Audekercke Remi, Smet Maria-Helena. Three-dimensional computer tomography-based, personalized drill guide for posterior cervical stabilization at C1-C2. Spine 2001,21;12:1343-1347. 22. Heilbrun MP. Computed tomography-guided stereotactic system. Clinical Neurosurgery 1984;31:564-581. 23. Horsley V, Clarke RH: The structure and function of the cerebellum examined by a new method. Brainl908;31:45-124. 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:86-90. 26. Mackerle J. Finite element modeling and simulations in orthopedics: a bibliography 1998- 2005. Comput Methods Biomech Biomed Engin2006;9(3): 149-199. 27. Maciunas RJ, Galloway RL, Latimer JW. The Application Accuracy of Stereotactic Frames. Neurosurgery 1994;35(4):682-695. 28. Kim NH., Sankar B V. Lecture notes: EML 4500. Finite Elements analysis and design. Gainesville: Department of Mechanical & Aerospace Engineering, 2004. A rr C D E Figure 1-1. Components involved in this investigation and the prototypes for the new technique for performing image-guided surgery with patient-specific reference guides. A) The graphical user interface that can create 3-dimensional patient anatomical model from a series of 2-dimensional medical image datasets, and the selection of a contour surface (in red) for the fabrication of a patient-specific reference frames. B) Physical patient's anatomical model created with Rapid prototyping for evaluating the stability of fit. C) Patient-specific 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. - 30e6 MPa 0.3 1000 N Figure 5-4. 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 LIST OF TABLES Table page 4-1 Relationship between the tetrahedral height and a measure of the rigidity of the rigid lay er ................... ............................................................. ................ 7 0 4-2 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 4-3 Single point criteria .................. .................. .................. ........... ......... .... 70 5-1 Displacement in the X, Y and Z direction from the example in Figure 5-4 tested with the ITKFEM module and ANSYS................................. ................................... 87 6-1 Measurement of the facial angles for the three anatomical models..............................99 6-2 Instability score, presented as the percentage of forces that resulted in instability and the total forces applied on a fram e .................................. ............... ............... 99 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 3-1) 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 3-2). 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 5-1. 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 1-1D. 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 A-C). 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 3-dimesional CT information to create the template guides based on the bony surface of the spine. Goffin et al. [20] printed the 3-D 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 2-5 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 computer-assisted surgery with individual templates on thirteen lumbar spines. Under image-guided 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 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 image-guided 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. 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 4-9) 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 4-11. 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 4-11, 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 Frame-base 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 frame-based biopsy a head ring is applied to the patient prior to imaging (Figure 2-2A). 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 2-1). The coordinates are set on the physical stereotactic system (Figure 2-2B) (which is directly fixed onto the rigid frame), providing a guide or trajectory for the surgeon [8, 35]. The main shortcomings of frame-base stereotactic surgery include the extend of time that the procedure takes (from to pre-surgical 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, frame-base image-guided surgery provides the gold standard for precision and accuracy for intracranial stereotactic procedures. Frameless stereotactic surgery has been presented as an alternative to frame-based procedures. The main difference between frame-based and frameless stereotactic systems, is that instead of a head ring, fiducials (landmarks) are used (Figure 2-3). 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 3-dimensional rendered images base upon the patient's CT 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;105-114. 2. Alker G, Kelly PJ. An overview of CT based stereotactic systems for localization of intercranial lesions. Computational Radiology 1984;8:193-196. 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. Computer-Assisted Orthopedic Surgery with individual templates and comparison to Conventional Operation Method. Spine 2001;26:365-370. 7. Bova, F. Grant proposal: Image-guided surgery based on rapid prototyping models. Department of Neurological Surgery, University of Florida. Gainesville, FL. March 2003 8. Brown RA: A computer tomography-computer graphic approach to stereotactic localization. J Neurosurg 1979; 50:715-720. 9. Brown RA: .A stereotactic head from for use with CT body scanners. Invest Radiol 1979;14:300-304. 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:522-530. 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 three-year contract to create an open-source 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 object-oriented 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 A-2). 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 A-2 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. I, I. Figure 4-18. 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 4-19. 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 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 cap-like 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 6-1B, 6-2B and 6-3B). 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 6-1A, 6-2A and 6-3A). 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 6-2B). 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 6-1A). 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 2-11; the results are shown in Table 6-1. 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 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 5-3. In all instances the same mechanical characteristics (E, v, A, L) where assigned and the same force vector (F) was applied (Figure 5-3). The results of the stiffness matrix as assembled by the ITK FEM module are shown in Equation 5-1, 5-2 and 5-3, which correspond to Figure 5-3A, 5-3B and 5-3C, 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 5-4, following the linear spring equation. 750 0 0 -750 0 0 0 0 0 0 0 0 kA 0 00 0 00 (5-1) 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 = (5-2) 0 -750 0 0 75 0 0 0 00 0 0 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 re-meshes 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 5-2. The rigid layer is created by meshing with tetrahedral elements the space between the selected surface and the top surface (Figure 5-2). The bottom point cloud become boundary condition points, fixed in all 3-dimensions 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 5-2). 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 5-2). 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. compressors" (Table 4-3). 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 4-15 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 4-16, 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 4-14 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 4-14). 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. Nasiofrontal angle I Nasiofacial angle / Nasiolateral angle .,' Figure 2-11. Facial angles [10]. The nasiolateral angle is calculated at the level of the radix. Figure 2-11. Facial angles [10]. The nasiolateral angle is calculated at the level of the radix. 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 A-4) [11]. The displacement function or displacement variation i(x) needs to be defined solemnly along the local x-axis. U(x) = ao + ax (A-4) 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 2-1, 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 (A-5) In the same way, when x is equal to L, i(x) is equal to u2. u(L) = u2 = ao + aL (A-6) The displacement variations of the truss element in terms of node displacements, ui and u2 can be calculated by substituting Equation A-5 and A-6 into A-4 and rearranging to yield Equation A-7 [11]. iT(x) I= 1- U + -kU2 (A-7) L) L Shape Functions The displacement function (Equation A-7) can be alternatively expressed as shape functions. There are as many shape functions for an element as there are nodes, and they have 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 age-Guided Surgery .................. ................ ............... ........ .. .......... .............. 19 Frame-Base 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 dimensional surface [29]. In the case of 3-dimensional 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 2-8B). 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 3-dimensional 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 2-8C) and infiltrated with cyanoacrylate glue to assure durability [33]. The 3-dimensional solid, or in the case of this investigation the patient-specific 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 3-dimensional printing, it starts with a solid object and removes unwanted material (Figure 2-8D). 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 MDX-650 (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 patient-specific reference frames created with the 3- Three patients with different facial anatomical features were selected as training data for testing the simulation program. For each case two patient-specific 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 3-D 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 cap-like 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. 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 patient-specific frame design STABLE? I ..i m.a.. -. N-- - i I, Step3: If YES, then build patient-specific reference - alI I 1 i- I 1 - " -i . Figure 2-9. Chart to visualize the significance of the present work. a* ai 1---------------------------------- Sagittal plane I I nI I 7 I r^"^ I CIL I I I .. I I r I I .- : I --. I I-J 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 4-4. 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 4-5. Nodes point on the rigid frame. Forces are not applied on edge points. LIST OF FIGURES Figure pe 1-1 Components involved in this investigation and the prototypes for the new technique for performing image-guided surgery with patient-specific reference guides................... 18 2-1 Pre-target localization and pre-surgical planning is done in a computer workstation.......36 2-2 Procedures related to frame-based image guidance surgery.................... ..................37 2-3 Frameless stereotactic surgery uses fiducials that are glued to the patient's skin.............37 2-4 Typical layout when frameless stereotactic surgery is employed.............................. 38 2-5 Tem plate designs for the spine ................................................ .............................. 38 2-6 R outer guides for craniotom y ................................................. ............................... 39 2-7 Biopsy guides ............................................................... ..... ...... ........ 40 2-8 R apid prototyping m machines ........................................... .................. ............... 41 2-9 Chart to visualize the significance of the present work ................. ............... ..............42 2-10 Facial landmarks ........................ ........ .. .... .... ........ .... .... 43 2-11 Facial angles ................................................... 44 2-12 Space truss element has six displacements to establish its deformed position ................45 2-13 Application of a force perpendicular to the axis of the space truss ..............................45 2-14 Space truss in local coordinates. ............................................................. .....................46 2-15 The linear tetrahedron has a unique numbering scheme .............................................46 3-1 The three models used to evaluate the methodology for testing unique fit .....................49 3-2 Masks built for testing the methodology of fit............ .................. ..................50 3-3 C ap -lik e referen ce fram e......................................................................... .....................5 1 3-4 Virtually created half sphere test m odel ........................................ ........................ 51 3-5 Placement of pins for a simple test on the real stability of the large and small masks......52 3-6 Coronal view of placement of the pins onto the real masks ............................................52 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. 3.5. Additional pins were located close to ends of the crosshairs lines as shown in Figure 3-5 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 3-6). A normal force of 5 N, applied by weight, was applied onto all pins for large and small masks (Figure 3-7). 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 6-1 through 6-3). In addition to the masks, two other surfaces were created to evaluate the methodology of fit. A cap-like reference frame (Figure 3-3) was created from CT image dataset and a half- sphere surface (Figure 3-4) 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 3-4). 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. In symbolic form Equation A-16 6uT is the transform of the virtual end displacements and f is the end force vector. Equation A-16 can be substituted into Equation A-3 to get Equation A- 17. SuTf = IETdV (A-17) v A relation between the element stiffness matrix ke and the strain displacement matrix B is finally accomplished by substituting Equations A-11 and A-12 into Equation A-17 and simplifying by canceling 6uT (Equation A-18) [4]. f = BTEBdV u= keu (A-18) The element stiffness matrix ke is determined by virtual work principles as defined by Equation A-19 [4]. ke = BEB dV (A-19) Gauss Quadrature To evaluate the integration in Equation A-18, the ITK FEM module uses a method of numerical integration called the Gauss quadrature (Equation A-20) [18]. n ke = wJB EB, (A-20) 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 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 5-1 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 Figure 6-3. Placement of infant masks onto solid model. A) Large Mask. B) Small Mask. Figure 6-4. Placement of a cap-like reference frame onto an upper head model. distinguishes springs that slip (or move tangentially) like springs 4, 5 and 6 in Figure 4-10, from springs that compress, like springs 1, 2 and 3 in Figure 4-10. Furthermore, Figure 4-13 shows two springs in compression; Figure 4-13A shows a spring with a TNT ratio of 0.5 and Figure 4-13B 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 4-14). 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 4-3 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 4-12 shows a 2-dimensional 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 4-12). 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 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 4-14. 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 4-18; 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 3-dimensions. 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 Following Figure A-2B, 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 A-2. * 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. Figure 3-3. Cap-like reference frame. A) Upper head model and B) cap-like frame. A 1 B Figure 3-4. Virtually created half sphere test model. A) Top view and B) side view. Figure 3-2. 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 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 3-dimensional 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 = (A-l) Yxy xy zx zx Equation A-i 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 A-2. aw, = J(&xx + a + ,C + xy T y + yz, +r yzx iV (A-2) V S = ISWadV (A-3) Equation A-2 and Equation A-3 (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. 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 4-9). 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 4-4 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 (4-4) 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 4-5 and Figure 4-6) to the final position of the spring (point c on Figure 4-9). ti= Ildi x Ni||/||Nil i is the number of springs (4-5) 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 4-10, 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 patient's anatomical face model (Figure 1-1). The approach used is the idealization of a contact mechanics problem by means of the simulating the patient-specific 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 life-like 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' 3-dimensional 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 image-guided surgery that involves rapid prototyping. The design of the patient-specific reference frame has the imperative requirement to have a unique fit onto the patient's anatomy, this is the major concern of this new image-guided surgery technology [33]. Therefore the work presented here, which numerically evaluates the mechanical stability and uniqueness of fit of the patient-specific 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 Figure 2-6. 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] 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:162-174. 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.. Image-based planning and validation of C1-C2 transarticular screw fixation using personalized drilled guides. Computer Aided Surgery. 2002;7:41-48. 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):257-270. 49. Yan X, Gu P. A review of rapid prototyping technologies and systems. Computer Aided Designl996;28(4):307-318. 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:45-48. 51. Zamorano LJ, NolteL, Kadi AM, Jiang Z. Interactive intraoperative localization using an infrared-based system. Stereotactic and Functional Neurosurgery 1994;63:84-88 |

Full Text |

PAGE 1 1 GUIDED SURGERY USING RAPID PROTOT YPING PATIENT-SPECIFIC 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 Image-Guided Surgery........................................................................................................... .19 Frame-Base 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 4-1 Relationship between the tetrahedral height and a measure of the rigidity of the rigid layer.......................................................................................................................... ..........70 4-2 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 4-3 Single point criteria...................................................................................................... ......70 5-1 Displacement in the X, Y and Z direction from the example in Figure 5-4 tested with the ITKFEM module and ANSYS.....................................................................................87 6-1 Measurement of the facial angles for the three anatomical models...................................99 6-2 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 1-1 Components involved in this investigation and the prot otypes for the new technique for performing image-guided surgery with patient-specific reference guides...................18 2-1 Pre-target localization and pre-surgical planning is done in a computer workstation.......36 2-2 Procedures related to fram e-based image guidance surgery..............................................37 2-3 Frameless stereotactic surgery uses fiduc ials that are glued to the patients skin.............37 2-4 Typical layout when frameless stereotactic surgery is employed......................................38 2-5 Template designs for the spine..........................................................................................38 2-6 Router guides for craniotomy............................................................................................39 2-7 Biopsy guides.............................................................................................................. .......40 2-8 Rapid prototyping machines..............................................................................................41 2-9 Chart to visualize the significance of the present work.....................................................42 2-10 Facial landmarks.......................................................................................................... ......43 2-11 Facial angles............................................................................................................. ..........44 2-12 Space truss element has six displaceme nts to establish its deformed position..................45 2-13 Application of a force perpendicula r to the axis of the space truss...................................45 2-14 Space truss in local coordinates.........................................................................................46 2-15 The linear tetrahedron ha s a unique numbering scheme....................................................46 3-1 The three models used to evaluate the methodology for testing unique fit.......................49 3-2 Masks built for testing the methodology of fit...................................................................50 3-3 Cap-like reference frame................................................................................................... .51 3-4 Virtually created half sphere test model............................................................................51 3-5 Placement of pins for a simple test on the real stability of the large and small masks......52 3-6 Coronal view of placement of the pins onto the real masks..............................................52 PAGE 9 9 3-7 The method employed to demonstrate the difference in stability between the large and small masks................................................................................................................ .53 4-1 The methodology for stability of fit i nvolves simulating a rigid layer onto of a deformable bed of springs..................................................................................................71 4-2 Displacement of a mask as it responds to external forces.................................................72 4-3 Triangulation of a customized pos itioning device proper to STL format..........................72 4-4 Direction of the Nmain.........................................................................................................73 4-5 Nodes point on the rigid frame..........................................................................................73 4-6 In a wedge three tetrahedral elements are fitted................................................................74 4-7 Schematic showing how the maximum rigidity fault is calculated...................................75 4-8 The behavior of a linear spring element under the application of an axial force...............75 4-9 Schematic that demonstrates how eac h springs displacement was evaluated by identifying the normal (ni) and tangential (ti) displacement magnitudes...........................76 4-10 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 4-11 Bounding volume for a spring to be in compression.........................................................77 4-12 Two-dimensional interpretati on of the single point criteria..............................................78 4-13 Schematic to demonstrate the importance of evaluating the percentage slip of TNT ratio for each spring.......................................................................................................... .78 4-14 Schematic to shows at the level at which the TNT threshold was defined........................79 4-15 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 4-16 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 4-17 A case of the large Asian mask, whic h shows the consista ncy of an upper TNT threshold of 0.7............................................................................................................... ...80 4-18 Schematic showing the relation betw een the volume of the entire frame (fvol) and the compression patch volume (cvol)........................................................................................81 4-19 Box ratio interpretation.................................................................................................. ....81 PAGE 10 10 5-1 Flow chart of program for testing the uniqueness of fit.....................................................88 5-2 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 5-3 Examples of al linear spring element or space truss tested in all 3-dimensions................90 5-4 Example of two tetrahedral elements us ed in the validation of the tetrahedral elements used in the program............................................................................................91 5-5 A corner frame that was used to validate the program......................................................92 6-1 Placement of Caucasian masks onto solid model............................................................100 6-2 Placement of Asian masks onto solid model...................................................................100 6-3 Placement of infant masks onto solid model...................................................................101 6-4 Placement of a cap-like reference frame onto an upper head model...............................101 6-5 Surface plot color coded by the maximu m full rotation angle for the Caucasian masks.......................................................................................................................... ......102 6-6 Surface plot color coded by the maximum full rotation angle for the Asian masks........102 6-7 Surface plot color coded by the maxi mum tilt angle for the infant masks......................103 6-8 Surface plot color coded by maximum tilt an gle for the cap reference frame design.....103 6-9 Surface plots color coded by maximum tilt angle for the half sphere test.......................104 6-10 Normalized histogram of box ratios for each tilt angle ( ) for the Caucasian large mask........................................................................................................................... ......104 6-11 Normalized histogram of box ratios for each tilt angle ( ) for the Caucasian small mask........................................................................................................................... ......105 6-12 Normalized histogram of box ratios for each tilt angle ( ) for the Asian large mask.....105 6-13 Normalized histogram of box ratios for each tilt angle ( ) for the Asian small mask.....106 6-14 Normalized histogram of box ratios for each tilt angle ( ) for the infant large mask.....106 6-15 Normalized histogram of box ratios for each tilt angle ( ) for the infant small mask.....107 6-16 Normalized histogram of box ratios for each tilt angle ( ) for the cap reference frame design......................................................................................................................... ......107 PAGE 11 11 6-17 Normalized histogram of box ratios for each tilt angle ( ) for the half sphere test.........108 6-18 Profiles of solid models.................................................................................................. .108 A-1 Differential element subjected to real stresses and in equilibrium..................................121 A-2 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 PATIENT-SPECIFIC 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 (3-dimensional pr inting) in conjunction with the virtual manipulation of computer-generated 3-D anatom ical models rendered from diagnostic images (CT and MRI) is a novel techniqu e for directing image-guided surgery. By incorporating a graphical user interface, the surgeon can plan an intracranial surgical procedure by fabricating a patient-specific reference frame (mask-like facial frame) through the us e of rapid prototyping technology. Early testing of these patient-specific 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 finite-element 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 patient-specific 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 3-D 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 cap-like 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 patient-specific reference frame that has a potential of beco ming a new type of image-guided surgery (IGS) technique [33]. The fabrication of the patient-specific reference fr ame involves coding a series of algorithms that render a detailed 3-dimensi onal life-like anatomical model from computer tomography (CT) or magnetic resonance imaging (MRI) [6] datasets (Figure 1-1). Onto a patient-specific 3-dimensional anatomical model a surface is selected (virtually painted) from which the patient-specific reference frames are cr eated. The ability to create, visualize, and manipulate 3-dimensional 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 patient-speci fic reference frame is to replace existing technology of framed-based (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 patient-specific 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 image-guided 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 non-significant 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 1-1D and 1-1E). 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 1-1D). For this reason, the accuracy of the patient-specific 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 patient-specific 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 tient-specific frame and PAGE 16 16 the patients anatomical face model (Figure 1-1). The approach used is the idealization of a contact mechanics problem by mean s of the simulating the patient-specific 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 life-like 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 3-dimensional 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 image-guided surgery that involves rapid prototyping. The design of the patient-specific 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 image-guided surgery technology [33]. Therefore the work presented here, wh ich numerically evaluates the mechanical stability and uniqueness of fit of the patient-sp 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 1-1. Components i nvolved in this investigation and th e prototypes for the new technique for performing image-guided surgery with pa tient-specific reference guides. A) The graphical user interface that can create 3dimensional patient anatomical model from a series of 2-dimensional medical image datasets, and the selection of a contour surface (in red) for the fabrication of a pati ent-specific reference frames. B) Physical patients anatomical model created with Ra pid prototyping for evaluating the stability of fit. C) Patient-specific 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 Image-Guided 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 3-dimensional 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 Brown-Roberts-Wells (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. Frame-Base 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 frame-base 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 3-dimensional virt ual model), onto a patients actu al, real-world 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 real-world patient and correlates the informa tion amongst them. The creation of a diagnostic plan, for both frame-base 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 frame-based or skin landmarks for frameless stereotactic surgeries (Figure 2-1) [6]. PAGE 21 21 Frame-base 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 frame-ba sed biopsy a head ring is applied to the patient prior to imaging (Figure 2-2A). 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 2-1). Th e coordinates are set on th e physical stereotactic system (Figure 2-2B) (which is directly fixe d onto the rigid frame), providing a guide or trajectory for the surgeon [8, 35]. The main shortcomings of frame-base stereota ctic surgery include th e extend of time that the procedure takes (from to pre-surgical 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, frame-base image-guided surgery provides the gold standard for precision and accuracy for intracranial stereotactic procedures. Frameless stereotactic surgery has been pr esented as an alternative to frame-based procedures. The main difference between frame-ba sed and frameless stereotactic systems, is that instead of a head ring, fiducials (landmarks) are used (F igure 2-3). 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 3-dimensi 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 frame-based and frameless procedures permit the same virtual computer planning to be performed by the surgeon; however frameless stereotactic surgery requires instrument-tracking 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 2-4). Although electromagnetic tracking allows fo r non-line 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 C1-C2 (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 2-5 PAGE 23 23 A-C). 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 3-dimesional CT information to create the template guides based on the bony surface of the spine. Goffin et al. [20] printed the 3-D 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 2-5 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 computer-assisted surgery with individual templates on thirteen lumbar spines. Under image-guided 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 C1-C 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 1-5 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 3-dimensional 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 image-guided surgery by use of patient-specific frames has been shown to have superior accuracy and precision when compared to previously documented frame-based 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 image-guided surgery is the frame-based 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 2-6 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 3-dimensional printing RP (Figure 2-8A and B), another is subtracting rapid prototyping (SRP) (Figure 2-8D). Both RP technologies provide a means to fabricate a physical 3-dimensi 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 3-dimensional 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 2-8B). 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 3-dimensional 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 2-8C) and infiltrated with cyanoacrylate glue to assure durability [33]. The 3-dimensional solid, or in the case of this investig ation the patient-specific 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 3-dimensional printing, it starts with a solid object and removes unwanted material (Figur e 2-8D). 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 MDX-650 (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 3-D 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 patient-sp ecific frame is created by selecting a reference surface (called selected surface) by painting th e surface of the virtual 3-dimensional rendered anatomical model [33] (Figure 1-1). 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 image-guided 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 patient-specific 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 patient-specific 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 patient-specific refere nce frame could maintain frame-based accuracy [33]. Furthermore, the cost for creating a frame and a whole set of attach ments to perform model PAGE 28 28 based image-guided 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 3-dimensional model from diagnostic image datasets. This program allows selecting a surface on th e virtual 3-dimensional 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 frame-on-patie 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 patient-specific 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 2-10 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 patient-specific 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 patient-specific 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 3-dimensional 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 2-node truss element in 3-dimensional space called space truss. The second one is the linear 4-node solid tetrahedral element in 2-dimensional 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 3-dimensional structures that ar e subjected to 3-dimens 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 ball-and-socket 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 2-12). 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 2-12). 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 ball-andsocket joints. For example, Figur e 2-13 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 (L-L = 0), there were no internal forces and no displacements either (Figure 2-13). The space truss is considered a simple element since it deforms only along one direction (along its local x-axis) [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 2-12). 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 2-1, an d the local node displacements u1 and u2 are caused by the local node forces f1 and f2, as seen in Figure 2-14. 2 1u u u (2-1) 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 2-2, where k is the stiffness coefficient and is equal to AE/L. ku f (2-2) With Equation 2-2 and the principles of equilibrium [28], the element stiffness matrix ( ke) for a particular element e can be calculated. As Figure 2-14A 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 2-2). 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 2-14B) 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 2-2. k k k kek or 1 1 1 1 L AEek (2-3) 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 2-4). 2 1 2 11 1 1 1 u u L AE f f or 2 1 2 1u u f fek (2-4) PAGE 33 33 The element stiffness matrix (Equation 2-4) 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 2-5. z y x z y x cos cos cos 0 0 0 0 0 0 cos cos cos T (2-5) Finally, for one element the local element stiffness matrix ke can be transformed into the global elements stiffness matrix Ke by Equation 2-6 (T means the transform of the matrix). T k T KTe e (2-6) For the space truss element, the elements stiffness matrix in global coordinates Ke results in a 6 x 6 matrix (Equation 2-7). 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 (2-7) 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 2-3); then to calculate the transformation matrix T (Equation 2-5) necessary to swap from local to global coordinates; and finally it is important to define the element stiffness matrix in global coordinates Ke (Equation 2-7). 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 2-8). 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 (2-8) The linear tetrahedron Solid elements are defined in 2-dimensions and can model a solid body. The linear tetrahedron is the simplest solid element in 2-dime 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 right-hand rule (Figure 2-14) and the calculation of the tetrahed rons volume must result in a positive value (Equation 2-9). 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 2-9, 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 (2-9) 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, multi-dimensional 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 2-10, where X is any vector in R space [13]. 0 AX XT (2-10) 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 2-1. Pre-target localizat ion and pre-surgical planning is done in a computer workstation. PAGE 37 37 Figure 2-2. Procedures related to frame-based 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 2-3. 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 2-4. 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 2-5. 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 2-6. 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 2-7. 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 2-8. 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 2-9. 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 patient-specific frame design STABLE? Step3: If YES, then build patient-specific reference PAGE 43 43 Figure 2-10. Facial landmarks. Medial-Lateral Caudal Cephalic Cheek (zygomatic bone) Orbital area Nasion Glabella Radix Superciliary arch Sagittal Plane Coronal Plane Frontal eminences PAGE 44 44 Figure 2-11. 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 2-12. Space truss element has six displaceme nts to establish its deformed position. There are u1 and u2 along its x-axis, v1 and v2 al ong its y-axis and w1 and w2 along its zaxis. Though the space truss is defined in 3-dimensional space, it will only develop axial forces, along its local x-axis. Figure 2-13. Application of a fo rce perpendicular to the axis of the space truss. There is no change in length (L-L = 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 2-14. 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 2-15. 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 3-1) 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 3-2). 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 5-1. 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 1-1D. 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 3-5 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 3-6) A normal force of 5 N, applied by weight, was applied onto all pins for large and small masks (Figure 3-7). 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 6-1 through 6-3). In addition to the masks, two other surfaces were created to evaluate the methodology of fit. A cap-like reference frame (Figure 3-3) was created from CT im age dataset and a halfsphere surface (Figure 3-4) 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 3-4). 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 3-1. 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 3-2. 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 3-3. Cap-like reference frame. A) Upper head model and B) cap-like frame. Figure 3-4. 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 3-5. 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 3-6. 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 3-7. 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 patient-specific 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 patient-specific 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 4-1), 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 4-2). 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 life-like facial models that represent three different patients anatomies we re constructed with 3-dimensi 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 patient-specific 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 patient-specific 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 patient-specific design. While this is only a piece of a much broade r investigation that seeks to ev aluate an alternative method to image-guided surgery, it is critical to this application of patient-s 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 3-dimensional surface as an assembly of planar triangles (Figure 4-3) [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 4-1 and 4-3). 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 4-1). 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 patient-specific 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 4-4). 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 4-5). 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 4-5, 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 -hand-rule as explained in Figur e 2-15. 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 4-7). As Equation 4-1 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 (4-1) 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 4-8b shows an axial force FY placing a spring under a state of compression and changing its length to a new length L1. Figure 4-8c shows the fr ee-body diagram of the system and the spring internal force FS, which is determined by the e quation in Figure 4-8d. In the Finite Element Method, linear spri ngs that reside in 3-dimensiona l space are modeled as trusses elements in 3-dimensional space. In order to define a linear spring in 3-dimesions 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 4-1 as Ai. In this manner, even though the triangular mesh in the STL file is irregul ar (Figure 4-3); 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 4-1 [46] a the elastic modulus (E) and the spring length (Lo) are kept constant. o i iL E A k (4-2) 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 ball-and-socket joint. External force magnitude The external force magnitude is calculated by Equation 4-3 and selecting an expected axial displacement ( L) of 1% the springs length. As s hown by Equations 3-2 and 4-2, 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 4-3. L L E spring area spring ce AverageFors* *0 (4-3) 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 4-1). The rigidity of the frame was qua ntified with the RDR ratio and Table 4-1 shows that reducing the height of the rigi d layer considerably affects the ri gidity of the rigid frame. For example, in Table 4-1, 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 1e-13. A difference of six orders of magnitude for the Er and Es maintains the rigidity and results in a reliable simulation. Table 4-2 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 4-9). 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 4-9). 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 4-4 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 4-5 and Figure 4-6) to the final position of the spring (point c on Figure 4-9). 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 4-10, 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 (4-5) ni compression = di Ni <= 0 (4-4) 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 4-9) 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 4-11. 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 4-11, 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 4-6). Numerically, it is defined by Equation 4-6 as the tangential displacement magnitude (t) divided by the sum of the normal displacement magnitude (n) and the tangential displacement magnitude (t) (Figure 4-9). n t t TNT (4-6) 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 4-10 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 4-6. 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 4-10, 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 4-10, from springs that compress, like springs 1, 2 and 3 in Figure 4-10. Furthermore, Figure 4-13 shows two springs in compression; Figure 4-13A shows a spring with a TNT ratio of 0.5 and Figure 4-13B 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 4-13B 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 4-13B 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 7-1, 7-2, and 7-3). 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 4-15 and 4-16 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 4-15 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 4-3). 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 4-15 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 4-16, 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 4-14 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 4-14). 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 4-14. 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 4-18; 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 3-dimensions. 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 4-18). 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 4-19). 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 4-18 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 4-7). vol volf c ratio Box (4-7) As Equation 4-7 shows, a box ratio is equal to 0 when the compression patch volume is equal to zero (or the mask is unstable) (Figur e 4-19). 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 4-19). 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 4-8). 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 (4-8) In summary, a total of 49 external forces will be applied to all non-e dge points onto a mask (if a mask has 300 non-edge 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 4-8). PAGE 70 70 Table 4-1. Relationship between the tetrahedral he ight and a measure of th e rigidity of the rigid layer. Table 4-2. 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. Er-Es Mean RDR Cholmods rcond 1.0e6 3/1000 2.14e-8 1.0e7 2/1000 2.14e-9 1.0e8 5/1000 2.14e-10 Table 4-3. 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 4-1. 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 4-2. Displacement of a mask as it responds to external forces. Figure 4-3. 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 4-4. 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 4-5. 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 4-6. 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 4-7. Schematic showing how the maximum rigidity fault is calculated. Figure 4-8. 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 4-9. 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 4-10. 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 4-11. 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 4-12. Two-dimensional 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 4-13A to correspond with spring 1 (with the tangential displacement in the opposite direction), and spring 2 co rresponds to the spring in Figure 4-13B. Figure 4-13. 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 4-14. Schematic to shows at the leve l at which the TNT threshold was defined. Figure 4-15. 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 4-16. 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 4-17. 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 4-18. Schematic showing the relati on between the volume of the entire frame (fvol) and the compression patch volume (cvol). Figure 4-19. 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 5-1 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 re-m 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 5-2. 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 3-dimensions and represents the patients surface. The spring elements that comprise the bed of springs extend from the point cloud to the selected surface (Figure 5-2). 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 5-2). 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 4-12. The same small magnitude is assigned to each force vector and all forty nine are applied onto each non-edge 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 non-edge 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 4-8). 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 color-coded 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 5-3. In all instances the same mechanical characteristics (E, A, L) where assigned and the same force vector (F) was applied (Figure 5-3). The results of the stiffness matrix as assembled by the ITK FEM module are shown in Equation 5-1, 5-2 and 5-3, which correspond to Figure 5-3A, 5-3B and 5-3C, 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 5-4, 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 (5-1) 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 (5-2) 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 (5-3) mm MPa mm mm N AE FL d 1 0 1500 1 0 2 752 (5-4) 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 5-5. Both ANSYS and ITK FEM provided the same displacement shown in Table 5-1. A Simple Model Validation The program was validated with a simple fram e in the shape of a corner (Figure 5-5). 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 5-5A). 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 5-1. Displacement in the X, Y and Z direc tion from the example in Figure 5-4 tested with the ITKFEM module and ANSYS. Node X Y Z 1 0.0000 0.0000 0.0000 2 0.0000 0.0000 -0.5200e-3 3 0.0000 0.0000 0.0000 4 0.0000 0.0000 0.0000 5 0.0000 0.0000 -0.5200e-3 PAGE 88 88 Figure 5-1. 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 re-mesh 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 5-1. 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 5-2. 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 5-3. Examples of al li near spring element or space truss tested in all 3-dimensions. 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 5-4. 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 5-5. 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 cap-like 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 6-1B, 6-2B and 6-3B). 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 6-1A, 6-2A and 6-3A). 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 6-2B). 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 6-1A). 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 2-11; the results are shown in Table 6-1. 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 6-1). 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 6-5 to 6-9) and b) a normaliz ed histogram for the box ratios for all values of ( = 0 to 360), for each tilt angle ( ) (Figure 5-12) accumulated over all points (Figure 6-10 to 6-17). 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 single-point and single stability criteria) for a full rotation ( = 0 to 360). At each color coded point on Figure 6-5 and 6-9, 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 6-5 to 6-9. 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 6-10 to 617). For example, in Figure 6-10 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 6-5A), the surface plot show s that all non-edge 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 6-5B). For the Asian large mask (Figure 6-6A) and fo r the infant large mask (Figure 6-7A), the surface plots show similar results to the large Caucasian mask (Figure 6-5A), 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 non-edge points (Figures 6-6 and 6-7). 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 6-7). 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 4-19). 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 6-10) 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 6-5B. 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 6-13), 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 6-15). 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 6-10). 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 6-12) 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 6-14). 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 6-12). Based on the facial angle measurements in Fi gure 2-11 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 6-13 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 6-18C. 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 cap-like reference frame (Figure 6-8B) and a half sphere surface (Figure 6-9) were evaluated with the met hodology for testing the uniqueness of fit. In reality the cap-like 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 cap-like 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 6-8). The box ratio hist ogram (Figure 6-16) 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 half-sphere surface (Figur e 6-9), 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 6-2 shows the instability scores for all masks and frames as defined in Equation 4-8. There is zero per cent instability for all large mask s and there is some percentage of instability for all small masks, the cap-like fr ame and the half-sphere surface, as was expected. The least stable case is the half-sphere 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 6-1. 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 6-2. 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% Cap-like frame 2695 241 9% Half-sphere test 1094 975 94% PAGE 100 100 Figure 6-1. Placement of Caucasian masks onto so lid model. A) Large Mask. B) Small Mask. Figure 6-2. Placement of Asian masks onto soli d model. A) Large Mask. B) Small Mask. A B A B PAGE 101 101 Figure 6-3. Placement of infant masks onto soli d model. A) Large Mask. B) Small Mask. Figure 6-4. Placement of a cap-like refe rence frame onto an upper head model. A B Coronal p lane Sagittal plane PAGE 102 102 Figure 6-5. Surface plot color coded by the maximum full rotation angle for the Caucasian masks. A) Large Mask. B) Small Mask. Figure 6-6. 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 6-7. Surface plot color coded by the maxi mum tilt angle for the infant masks. A) Large Mask. B) Small Mask. Figure 6-8. 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 6-9. 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 6-10. 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 6-11. 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 6-12. 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 6-13. 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 6-14. 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 6-15. 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 6-16. 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 6-17. Normalized histogram of box ratios for each tilt angle ( ) for the half sphere test. Figure 6-18. 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 cap-like 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 patient-specific reference frames that impr oves the new and promising technology to perform image-guided surgery. The methodology will provid e the clinician with prior information of frame-on-patient 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 image-guided 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 patient-specific 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 force-displacement 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 ball-and socket in 3-dimensional 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 load-deformation relation of the entire structure. The ot her is referred as a local coordinate system (x, y and z) and it descri bes the force-displacement 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 x-axis 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 1-dimensional truss elem ent has 1 DOF per node, a 2-dimensional 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 force-displacement relations can be derived from the work-energy principles; the main reason for this being th at by using work-energy 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 rigid-body 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 strain-stress 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 3-dimensional 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 (A-1) Equation A-1 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 A-2. dV WV zx zx yz yz xy xy z z y y x x e (A-2) dV Wv T e (A-3) Equation A-2 and Equation A-3 (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 A-4) [11]. The displacement function or displacement variation (x) needs to be defined solemnly along the local x-axis. x a a x u1 0) ( (A-4) 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 2-1, 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 (A-5) In the same way, when x is equal to L, (x) is equal to u2. L a a u L u1 0 2) ( (A-6) The displacement variations of the truss el ement in terms of node displacements, u1 and u2 can be calculated by substituting Equation A-5 and A-6 into A-4 and rearranging to yield Equation A-7 [11]. 2 11 ) ( u L x u L x x u (A-7) Shape Functions The displacement function (Equation A-7) 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 A-7 represent the shape functions (Equation A-8). L x N L x N2 1 1 (A-8) For convenience the displacement variation function ( (x)) can be expressed in terms of shape functions and node displacements (Equation A9). Equation A-10, 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 (A-9) = Nu (A-10) Strain Displacement Equation Getting back to the aim of representing the node displacement in term of stresses and strain, Equation A-11 shows the relationship between strain ( ) and displacement variation function ( (x)), for a truss element [11]. dx u dx (A-11) From Equation A-11, notice that a truss element only withstand strains along the local xaxis. x xu dx du = D (A-12) 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 A-10 into Equation A-12. Bu DNu (A-13) In order to get the strain displacement matrix, B, the derivatives of the shape functions need to be calculated (Equation A-13). 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 A-14), as required by the principle of virtual work for deformable bodies [4]. Bu E (A-14) So far both strains (Equation A-13) and the st resses (Equation A-14) 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 3-1) 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 A-15 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 (A-15) f uT eW (A-16) PAGE 118 118 In symbolic form Equation A-16 uT is the transform of the virtual end displacements and f is the end force vector. Equation A-16 can be substituted into Equation A-3 to get Equation A17. dVv T T f u (A-17) A relation between the element stiffness matrix ke and the strain displacement matrix B is finally accomplished by subs tituting Equations A-11 and A-12 into Equation A-17 and simplifying by canceling uT (Equation A-18) [4]. u k u B B fe dV Ev T (A-18) The element stiffness matrix ke is determined by virtual wo rk principles as defined by Equation A-19 [4]. dV Ev T B Bek (A-19) Gauss Quadrature To evaluate the integration in Equation A-18, the ITK FEM module uses a method of numerical integration called the Gauss quadrature (Equation A-20) [18]. i T i i n i i eE w B B J1k (A-20) 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 three-year contract to create an open-source 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 object-oriented 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 A-2). 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 A-2 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 A-2B, 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 A-1. 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 A-2. 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;105-114. 2. Alker G, Kelly PJ. An overview of CT base d stereotactic systems for localization of intercranial lesions. Comput ational Radiology 1984;8:193-196. 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. Computer-Assisted Orthopedic Surger y with individual templates and comparison to Conventional Operation Method. Spine 2001;26:365-370. 7. Bova, F. Grant proposal: Image-guided su rgery based on rapid prototyping models. Department of Neurological Surgery, Universi ty of Florida. Gainesville, FL. March 2003 8. Brown RA: A computer tomography-comput er graphic approach to stereotactic localization. J Neurosurg 1979; 50:715-720. 9. Brown RA: .A stereotactic h ead from for use with CT body scanners. Invest Radiol 1979;14:300-304. 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:522-530. 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:200-204. 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. 3D-CT ba sed, personalized drill guide for posterior transarticular screw fixation at C1-C2: Technical Note. NeuroOthopedics1999;25:47-56. 21. Goffin Jan, Van Brussel Karel, Martens Kr isten, Vander Sloten Jos, Van Audekercke Remi, Smet Maria-Helena. Three-dimensiona l computer tomography-based, personalized drill guide for posterior cer vical stabilization at C1-C 2. Spine 2001,21;12:1343-1347. 22. Heilbrun MP. Computed tomography-guide d stereotactic system. Clinical Neurosurgery1984;31:564-581. 23. Horsley V, Clarke RH: The structure and function of the cerebellum examined by a new method. Brain1908;31:45-124. 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:86-90. 26. Mackerle J. Finite element modeling and si mulations in orthopedics : a bibliography 19982005. Comput Methods Biomech Biomed Engin2006;9(3):149-199. 27. Maciunas RJ, Galloway RL, Latimer JW. The Application Accuracy of Stereotactic Frames. Neurosurgery1994;35(4):682-695. 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:1257-1287. 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: Frame-less and Frame-based stereotaxy? How to choose the appropriate procedure. Zent ralbl Neurochir.2003;64(1): 1-5. 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):81-98. 34. Rajon DA, Bova FJ, Bhasin RR, Friedman WA. Poster Pr esentation: Fabrication of custom surgical guides using rapid-prototypi 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:1042-1044,. 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):499-505,1987. 37. Schileo E, Taddei F, Malandrino A, Cristofolin i A, Viceconti M. Subject-specific 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. Real-time fin ite-element simulation of linear viscoelastic tissue behavior based on experimental data. IEEE comput Graph Appl 2006;26(6):58-68. 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):247-256. 41. Shewchuck JR. What is a good linear elemen t? Eleventh International Meshing Round Table. New York: Sandia Nationa l Laboratories September 2002:115-126. 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:349-350. 44. Spiegel EA, Wycis HT: Thalamotomy and related procedures. JAMA 1952;148:446-451. 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:162-174. 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.. Image-based planning and validation of C1-C2 transa rticular screw fixation using personalized drilled guides. Computer Aided Surgery. 2002;7:41-48. 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):257-270. 49. Yan X, Gu P. A review of rapid prototypi ng technologies and system s. Computer Aided Design1996;28(4):307-318. 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:45-48. 51. Zamorano LJ, NolteL, Kadi AM, Jiang Z. Inte ractive intraoperative localization using an infrared-based system. Stereotactic and Functiona l Neurosurgery 1994;63:84-88 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 image-guided 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. |