UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations  Vendor Digitized Files   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
HIGH PRECISION TREATMENT FOR FRACTIONATED STEREOTACTIC RADIOTHERAPY By CHINGCHONG YANG A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1994 ACKNOWLEDGEMENTS I would like to express my sincere appreciation to the members of my supervisory committee for their kind help and guidance throughout this work. I would also like to express special thanks to my committee chairman, Dr. Francis Bova, for his invaluable suggestions and sincere instruction through this research. I would also like to give special thanks to Dr. Carl Crane for invaluable initiation about the solving matrix problem. Also I would like to thank Dr. William Friedman for his kind support for the whole project and Dr. William Mendenhall for his comments on the radiobiological aspects of this project. Special thanks should be given to Dr. Lawrence Fitzgerald for his constant encouragement and comments throughout the research period. I also would like to thank Professor James Tulenko and Dr. William Properzio for discussions, corrections, and encouragement. I would like to give special thanks to Ms. Pam Mydock for her excellent job in proofreading the draft of this manuscript. I also could not forget my fellow graduate students and the physicists, engineers, and doctors at the Shands Cancer Center who encouraged me to go on for my graduate career: Kevin Kalbaugh, Linda Ewald, Rene Craig, Ruth Ann Good, Fred Harmon, Sanford Meeks, Soon Huh, Dale Moss, Don Dubois, Revlon Williams, Duane Sandene and Drs. Nancy Mendenhall, Jim Parsons, Michael Sombeck, Jatinder Palta, and Chihray Liu. Finally I would like to thank my sweet Ph.D. wife to has accompanied me through the years without any complaint and my parents for their full support, both spiritually and financially. TABLE OF CONTENTS P ACKNOWLEDGEMENTS ................................. TABLE OF CONTENTS .................................. ABSTRA CT .......................................... CHAPTERS 1 INTRODUCTION .................................. High Precision Radiation Therapy ........................ Radiobiology for Fractionated Treatment . Patient Repositioning ................ Fractionated High Precision Radiotherapy . 2 REVIEW OF LITERATURE .......... Stereotactic Radiosurgery Systems . Radiobiology for Fractionated Treatment . Localization Errors ............ Immobilization Methods ......... Fractionated Stereotactic Radiotherapy Systems Refixation Techniques .......... Relocalization Methods .......... 3 THE UF STEREOTACTIC SYSTEMS .... The UF Stereotactic Radiosurgery System .. Hardware Characteristics ......... Treatment Planning Patient Treatment Procedures ..... . . . The UF Fractionated Stereotactic Radiotherapy System . Hardware Characteristics . . Localization and Refixation Procedures . iv ages . ii .iv .vii 1 1 3 4 5 S7 S7 S9 12 13 15 16 18 21 21 21 22 24 24 24 27 I . Patient Treatment ................... ........... 29 4 TARGET RELOCALIZATION BY 3D CONTOUR POINTS ....... 30 Statement of Problems .................. ... .......... 30 Methodology .................... .................. 32 5 3D COORDINATE MEASURING SYSTEMS ................ 39 Laser Tracking Interferometer System ....................... 39 Camera System .......................... ........... 41 Articulated Arm System .............................. 49 Frequency Matching System .................. .......... 50 Ultrasonic System .................................... 51 Conclusion ................... ................... .52 6 THE FITTING OF 3D REFERENCE DATA SETS ............. 53 Basic Principles ...................... ............ 53 The Iterative M ethod ................................ 57 The Singular Value Decomposition (SVD) Algorithm ............. 58 The Vector Analysis of Robotic Method ................. .... 66 The Optimization Methods ............................... 68 Hooke and Jeeves Pattern Search Method ............... 69 The Downhill Simplex Method in Three Dimensions ......... 72 The Simulated Annealing Optimization Method in Three Dimensions ............................ 77 Conclusion ..................... ..................... .82 7 ANALYSIS OF RELOCALIZATION ERRORS ................ 84 Comparison of Algorithms ...............................84 Phantom Image Tests .................................. 98 3D Camera System Data Tests .......................... 124 Correlation of CT Image and 3D Camera Data Sets ............. 131 8 THE HARDWARE SYSTEMS AND THE VERIFICATION PROCEDURES .. ........ ......... ... ........... ....... .. ...... ...136 System Characteristics ................................ 136 High Precision Camera System ..................... 136 The RealTime Digitizing Probe ................ .... 142 The Patient Immobilization System ................... 145 The Patient Bite Plate Verification System ............... 145 The Verification Procedure for Fractionated Treatment . ... 146 Conclusion ....................................... 156 9 CONCLUSION ................... ................. 157 Conclusion ......................................... 157 Future Work ...................... ................ 159 APPENDICES A COORDINATE TRANSFORMATION ..................... 163 B COMPUTER PROGRAMS ............. ................ 171 REFERENCES ............. ........................... 269 BIOGRAPHICAL SKETCH ............ .. ................ 278 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy HIGH PRECISION TREATMENT FOR FRACTIONATED STEREOTACTIC RADIOTHERAPY By ChingChong Yang April, 1994 Chairperson: Frank J. Bova Major Department: Nuclear Engineering Sciences Fractionated Stereotactic Radiotherapy (FSR) refers to high precision, multifraction, externalbeam irradiation that utilizes accurate refixation techniques. Three optimization algorithms were developed to calculate sixdegreesoffreedom motion parameters for patient repositioning in FSR. These algorithms provide the parameters for patient rotation and translation. Simulation results have shown a maximum target displacement error of 0.4 mm and an orientation error of 0.01 degree. Phantom image and realtime 3D data studies have been performed to verify the accuracy of the process. A new repositioning system is proposed that provides for high rotational and translational precision in fractionated treatment. An independent verification procedure is also developed and verified to assure the high treatment precision in the intracranial region. CHAPTER 1 INTRODUCTION High Precision Radiation Therapy Historically, a large single dose divided into many small fractions has been employed in cancer treatment. This is to maximize the tumor dose and to take advantage of the differential radiobiological responses between tumor and normal tissues. In order to achieve the constant treatment accuracy for each treatment course, various fixation methods have been developed to maintain reproducible patient treatment position. Depending upon the different treatment area, treatment errors from several millimeters up to a centimeter have been observed [Rab85]. In order to achieve high precision in radiation therapy positioning, a treatment system with high target localization and irradiation accuracies has been developed. The newly developed irradiation technique is termed stereotactic radiosurgery since it performs noncoplanar radiation treatment without invasive open surgery for intracranial diseases. With the use of a reference frame attached to the patient's skull, stereotactic radiosurgery allows high precision for both target localization and treatment. The fixation technique of radiosurgery is invasive and it is thus not easy to perform dose fractionation. Therefore, a large, single fraction of exposure is currently administered in a stereotactic radiosurgery procedure, thus the radiobiological advantages of fractionation have not yet been realized. 2 Stereotactic radiosurgery is used to treat patients who are not good candidates for conventional neurosurgery. A single high dose delivery to a vascular lesion has been acknowledged as an effective method to cause lesion thromboses. Through the invasive fixation of a haloshaped of metal immobilization device, the patient's anatomy relative to the BrownRobertsWells (BRW) neurosurgical stereotactic reference ring and frame system is fixed, and a CT image set with multiple slices is performed without removing the head ring. The images provide geometrical relationships of the head relative to the head ring and thus any point in the head can be assigned an X, Y, Z coordinate relative to the known geometry on the ring. Once the target coordinate related to the head ring is known, the treatment planning process can begin. After interactive optimization of treatment planning, the dose prescription is selected and patient treatment proceeds. Radiation is delivered to a point whose reference coordinates are preserved through the planning and treatment processes since the head ring remains fixed on the patient until treatment is complete. With a single fraction stereotactic treatment, the prescribed dose can be delivered exactly to the target area. However, once the ring is removed from the patient, the relative coordinates disappear and then fractionation becomes impossible. The second fraction of radiation could be provided by repeating the same invasive localization procedure. Clinically, this is impractical. Therefore, stereotactic radiosurgery achieves high precision but sacrifices the radiobiological differentiation between tumor and normal tissues. Stereotactic radiosurgery manipulates the treatment planning design by direct implementation of dosimetric information on the image basis. This virtual simulation 3 procedure involves no patient contact while obtaining the best noncoplnar radiation beam combinations on the computer. Virtual simulation differs from most of the radiation oncology routine practices by replacing conventional treatment simulation with 3D image data basis and computer software. The same geometry used in imaging the patient will be used in treatment, and the data sets provide the basic information for beam delivery. Application of virtual simulation requires the capabilities of transferring the treatment planning geometry from the computer to the treatment position in a way which is accurate, reproducible, and efficient [She90]. The mathematical projection of radiation beams to the target can be transferred to the exact position under the treatment machine since the mathematical simulation process on computer conveys the same information in the real treatment condition. Radiobiology for Fractionated Treatment The single fraction treatment for radiosurgery of Arterial Venous Malformation (AVM) has been proved to be clinically very effective [Fri92a]. The AVM suffers sufficient cellular damage from a single dose of radiation to the vascular structure that over a period of many months the abnormal vessels thrombose. However, the use of a single large dose for the treatment of malignancies is contrary to all the radiotherapy experience that has been accumulated over the past decade. In the treatment of a malignancy, the therapeutic ratio can be greatly improved by fractionation. Fractionation increases the cellular depopulation of the tumor because of the reoxygenation while reducing the damage to critical late responding normal tissues. The biological principle of fractionated treatment of tumors is based on the experimental evidence that there is 4 a difference in shape between the doseresponse characteristics of early responding tissues and tumors and late responding tissues [Hal93]. Therefore, treating malignant tumors with a single fraction will result in a suboptimal therapeutic gain between local tumor control and the normal tissue late effects even for small tumors. An improved therapeutic ratio is expected if the treatments are fractionated. Patient Repositioning Radiation therapy in cancer management is a complex process that involves many procedures to achieve the best treatment results. However, during the process, inaccuracies will occur due to a variety of procedures, such as localization of tumor volume, treatment setups, and immobilization and repositioning techniques. Routine radiation therapy is usually administered under multifractionated treatments. During the treatment process, special attention is required in four different areas. First, the prescribed doses should be defined based upon clinical research and continuous evaluation; this is beyond the scope of our discussion. Second, the machine output and dose calculation should be correctly specified. In this area, errors are already minimized due to the improvement of dosimetry systems, NIST1 traceable ionization chambers, an improved dosimetry protocol such as TG21 [Tas83], and redundant calibration procedures by medical physicists. The third and the fourth requirements involve the technical setup of the daily treatment and accurate dose delivery to the tumor for the duration of each treatment. Among these steps, tumor localization and radiation field misalignment are the most errorprone factors for radiation therapy. Errors in 'NIST: National Institute of Standards and Technology, Gaithersburg, MD. 5 localization can lead to reduced doses at the sites of tumors or to excessive doses to normal tissues. If the target is missed, doses will be delivered to normal tissues, which increases the complication rate of normal tissues and decreases the dose required at the target. Because the target receives a lower than prescribed dose, tumor control can be adversely affected. Incorrect placement of the radiation field relative to patient anatomy results in the same effects described above. This misalignment primarily comes from patient movement and repositioning errors. Rabinowitz et al. [Rab85] have analyzed the simulation and port films to determine the field misalignment errors. By using the existing immobilization techniques, the average radiation field discrepancies for consecutive films vary from 3.5 mm to 9.2 mm for the head and neck region and the thorax region. Since the bony structures of the head and neck region could be strongly immobilized, an average error range of 3.5 mm is the minimum field alignment error of the whole body area according to study results. The localization method of stereotactic radiosurgery are potentially useful to assure patient positioning reproducibility during each treatment fraction since the reported accuracy is higher than that of any treatment techniques applied in the generalized routine radiation therapy. Fractionated High Precision Radiotherapy Fractionated high precision radiotherapy refers to high precision, multifraction, externalbeam irradiation using stereotactic localization methods. In stereotactic procedures, with a single fraction radiation treatment, the prescribed dose can be accurately delivered to the target area. In order to achieve high precision fractionated 6 radiotherapy, it is necessary to develop a technique that allows repeated fixation of patient position for multifractionated radiation delivery schemes. Radiosurgery has been successfully implemented by hardware and software to determine the target size and location, to perform treatment planning, and to deliver proper dose by the plan. However, for certain tumors, the single fraction treatment may provide inferior clinical results when compared to conventional radiation therapy. On the other hand, conventional fractionated radiation therapy may introduce setup errors not present in the stereotactic treatments with only single irradiation delivery. The high precision multifractioned radiation treatment could greatly increase the possibility of local control for various benign and malignant tumors of the brain and skull such as acoustic neuroma, meningioma, pituitary adenoma, malignant glioma, locally recurrent nasopharyngeal carcinoma, and chordoma and chondrosarcoma of the skull base and cervical spine. Patient repositioning is an important factor to determine success of the tumor control for fractionated treatment. Therefore, the improvement of patient repositioning for fractionated stereotactic radiotherapy is a critical step in achieving better local control. CHAPTER 2 REVIEW OF LITERATURE Stereotactic Radiosurgery Systems Stereotactic radiosurgery has four qualities that distinguish it from routine radiation therapy. They are 1) a target volume that does not lend itself to routine therapeutic simulation techniques, 2) the necessity for true three dimensional treatment planning, 3) stricter isocenter criteria for gantry and patient motions, and 4) small beam size [Bov90a]. Because of the complex structures and critical organs in the brain, high precision and accuracy are critical for treating intracranial diseases. In the 1950s, the Swedish neurosurgeon Lars Leksell built a stereotaxic instrument for open brain operations. He then modified this unit and irradiated the target through a large number of small beams by fixing a semicircular stereotaxic frame to the patient's skull. He named this method radiosurgery. The original radiation source used was a 280 kV Xray tube [Lek51]. This source of radiation was later replaced by multiple converging cobalt sources, initially 179 and later 201 beams. This device, called the gamma knife, has been in clinical use since 1968 [Lek71]. In 1958 and 1968, other teams in Berkeley [Tob64] and Boston [Kje77] gained experience in the use of charged particles for treatment. These beams of protons or helium ions were produced by synchrocyclotrons, taking advantage of the rapid beam 8 falloff characteristics (i.e., Bragg peaks) to spare adjacent healthy tissues. While these methods are effective, they need very expensive, delicate equipment and are only available in a few large medical or research centers. In 1974 the use of modern linear accelerators (LINACs) as a radiosurgical alternative was proposed [Lar74]. For the past decade clinical studies involving radiosurgery with LINACs have been extensively pursued. Betti [Bet84] reported the first usage of a LINAC as a standard treatment machine for radiosurgery. Colombo [Col85], Hartman [Har85], Sturm [Stu87], and Lutz [Lut88] subsequently have reported modified LINAC radiosurgery techniques to deliver highly focused radiation to a target. However, most of the groups used the stereotactic reference ring by attaching it to the patient's skull and only delivered a single dose fraction for high precision treatment of intracranial targets. In LINACbased dose delivery systems, the accuracy and precision of the LINAC gantry and the patient support rotations are generally contained within 2 mm diameter sphere. Since they are seldom concentric, the combined inaccuracy can exceed a +2 mm radius sphere [Bov90a]. These errors are significantly greater than those of target localization, and potentially could limit the precision of the procedure. Friedman and Bova [Fri89a] at the University of Florida developed a three dimensional sliding bearing system to couple an auxiliary collimator to the LINAC head. An overall system accuracy of 0.2 mm with a standard deviation of 0.1 mm was achieved. This level of accuracy and precision will permit significant increases in the spacial accuracy of image localization before the dose delivery system becomes the limiting factor in radiosurgery. To date, this is the most precise and 9 accurate system of those described in the literature. It is commercially available as the Philips SRS 200 Stereotactic Radiosurgery System. With this radiosurgery system, a high precision treatment of radiation beams has been successfully applied clinically to irradiate inoperable intracranial targets. Using a routine simulator, the treatment beam cannot visualize the patient. Routine stereotactic radiosurgery is a true application of virtual simulation in that the mathematical description of the patient is utilized to plan the treatment process. Then the treatment is based on the virtual planning process, in which the radiation beams are manipulated mathematically to irradiate proper patient anatomy. This virtual process with 3D patient information and Beam's EyeView [Goi83] function could improve the precision of radiation delivery relative to the correct patient anatomy. Radiobiology for Fractionated Treatment In 1914, Schwarz pointed out that a single dose radiation treatment might not be efficient because some of the tumor cells were in different stages of radiosensitivity, and there was a better chance that multiple exposures could hit the cell in a radiosensitive phase [Per87]. However, the basic parameters of dose versus time are complex functions. Moss noted the following advantages of dose fractionation [Mos79]: (1) There is a reduction in the number of hypoxic cells through cell killing and reoxygenation, (2) There is a reduction in the absolute number of tumor cells by the initial fractions with the initial killing of the betteroxygenated cells, (3) Blood vessels compressed by a growing cancer are decompressed as the cancer shrinks, permitting better oxygenation, (4) Fractionation exploits the difference in recovery rate between normal tissues and neoplastic tumors, and (5) The acute effects of single doses of radiation can be decreased with fractionation. The basic parameters in dosetime considerations constitute a complex function, which interprets the interdependence of total dose, time, and number of fractions to produce a biological effect in the tumor. The fractionation schemes to generate expected biological responses of both tumor and normal tissues are closely related to the four "Rs" of radiation [Hol91, Per87]: (1) Repair of sublethal damage of potentially lethal damage, (2) Reproductions of cells between each fraction, (3) Redistribution of cells throughout the cell cycle, and (4) Reoxygenation after radiation exposures. The repair of sublethal damage produces the clinical observation that larger total doses can be tolerated when the treatment is divided into multiple small fractions [Hol91]. This biological factor indicates that the fractionation of radiation therapy could introduce clinical tumor control benefits because the repair mechanisms of normal tissues and tumor are differentiated [Gei87]. A small fraction dose allows normal tissues to recover to a larger fraction of the cell population for the next irradiation, while the less recovered tumor cells continue to be damaged by the following irradiation. Fractionation treatment attempts to deliver radiation therapy beams to achieve the total destruction of 11 the proliferative capability of the tumor cells, while leaving a certain level of function of the normal tissues or organs involved. A linearquadratic survival relationship (a/f3 ratio) of the cell responses to radiation is used to assess the fractionated doses (Fig. 21). This characteristic curve is represented as follows: Log, S = ad + 3d2 (2.1) S is the survival fraction of cells; a represents the linear, nonseparable component of log cell kill; and 3 represents more effective killing after some repair process has been eliminated with increasing dose [Per87]. cell ac cell kill survival j ? cell kill a/a Dose per fraction Figure 21: LinearQuadratic model for multifraction dose assessment [Per87] AVM, for instance, is a lateresponse tissue with small a/O3 ratio. Fractionation will not generate benefits for destruction of this lateresponding tissues. However, for tumor control, the a/f3 ratio is larger, which means that the dose fractionation spares late responding tissues more than early responding tissue. The biological benefits are then observed. In order to perform a fractionated procedure in the stereotactic environment, it 12 is important to identify the target location repeatedly and precisely. Fractionated treatment of radiotherapy requires reposition of patients for each treatment course relative to a specified coordinate system. Clinically, the error ranges have been extensively examined and evaluated. Many methods have been created to minimize the errors. Localization Errors Accurate daily positioning of the patient in the treatment position and reduction of patient movement during treatment are essential to deliver the prescribed dose and achieve the planned dose distribution. In the study of anatomic and geometric precision achieved in the treatment of mantle and pelvis port films, one third of the localization errors were caused by patient movement [Hau72]. It then shows that the difficulty of positioning the customized shielding blocks is caused by the movement of the external landmarks on the surface of the patient, which are used to indicate the correct positioning of the shields. When the landmarks move, the shields move in relation to the internal anatomy, and localization errors occur [Mar74]. Hodgkin's disease provides a difficult metastasis site with resulting frequent field placement errors. Marks et al. reviewed the verification films to detect geometric miss or localization error and found that the maximum number of 44 errors resulted from a total of 77 verification films during one study period. A maximum of 55% of the localization error has been observed [Mar74]. Those localization errors result from the external landmark changes because of movement of skin and subcutaneous tissues over underlying skeleton and musculature [Mar76]. Rabinowitz et al. studied the accuracy in all radiation fields alignment in clinical practice. They found that the mean worstcase 13 discrepancy averaged over all treatment sites was 7.7 mm, the lowest mean value was 3.5 mm in the head and neck region with proper fixation, and the highest mean value was 9.2 mm in the thorax region [Rab85]. Errors of over 1 cm were found in 10% to 23% of the radiation fields and an increase of infield relapse rate from 7% to 33% was found, when the radiation beams did not encompass the tumor region adequately [Hu189]. Byhardt et al. [Byh78] separated field misplacement into four categories: field malposition, field malrotation, patient malposition, and block malposition. An overall field placement error was 15% (10% if the error limit was set at 1.0 cm or greater). Field malposition, which means cephalocaudad and/or transverse field shift, was the most common field misplacement and composed 74% of all errors. Goitein [Goi86] summarized the reasons for nonuniform dose delivery and concluded that some degree of field misplacement is virtually inevitable. He then suggested a margin of approximately 1.5x o2+p 2 where a. is the standard deviation of patient motion and oa is the standard deviation of the penumbra of the dose profile. According to clinical surveys in routine radiation therapy, several millimeters up to centimeter range inaccuracy in relocalization are observed for fractionated treatment. Those inaccuracies increase normal tissue complications as well as decrease tumor control. Immobilization Methods The precision required in radiation treatment is most critical in management of tumors in the head and neck area where the proximity of the optical nerve structure, brain, and cervical spinal cord necessitates a high degree of precision in both patient 14 positioning and immobilization. However, to improve treatment accuracy, many repositioning and immobilization techniques have been created. Khan [Kha84] mentioned some general guidelines for maintaining the patient treatment position, such as a bite block system, a head clamp and mask for head and neck treatment, and a partial body mold or cast for irradiation of other body parts. Choice of the technique depends on the location of treatment fields. Williamson [Wil79] developed a styrofoam block for defining posterior landmarks in position on a simulator to improve the lateral treatment accuracy. Verhey et al. [Ver82] developed a number of immobilization schemes that permit precise daily positioning of patients for radiation therapy. They studied different patient positions that indicate the movement from a mean of 1.4 mm for different treatment areas. A daily reproducibility of patient position using laser alignment and pretreatment radiographs for final verification of positions indicates that the initial laser alignment can be used to position the patient within 2.2 mm + 1.4 mm. By using a combination of pretreatment radiographs and rather simple contoured plastic masks or casts, holding intracranial treatment movement to less than 2 mm is possible in over 80% of the treatment studied. This result indicates that rigid immobilization devices can improve the radiation therapy, especially in the head and neck regions [Ver82]. Thermal plastic (Aquaplast") has been employed for head and neck treatment because of its low cost and easy utilization [Ger82]. In clinical practice, physicians tend to use line marks and thermoplastic together to avoid possible treatment field localization errors. Custom molds constructed with the use of polyurethane formed to patient 15 contours are also becoming major aids to immobilization and repositioning. Some typical sites, including breast and upper body for Hodgkin's disease, need molds for strong support of body structures. Another immobilization device is the vacuumform body immobilizer, which is often used for rigid body positioning and easy adjustment during the treatment schemes. Some institutions developed an elaborate casting technique to immobilize and reposition patients during radiation treatment. This technique is not popular because there is a potential for movement inside the cast due to patient weight loss or the regression of tumor. All those techniques try to reposition the patient properly; however, they rely on individual skill and differ from institution to institution. As discussed in the previous section, utilization of these immobilization techniques results in an average minimum inaccuracy of 3.5 mm, which means there is room for improving the treatment precision. Fractionated Stereotactic Radiotherapy Systems As mentioned previously, stereotactic radiosurgery is a high precision radiation therapy technique to treat intracranial disorders. Invasive and single fraction radiosurgery has been proved an effective methodology to treat inoperable brain diseases. However, in order to differentiate tumor dose response with normal tissues, fractionated stereotactic radiotherapy may be the preferable choice. Therefore, patient repositioning with a high degree of accuracy to spare critical structures inside the brain becomes the ultimate goal. Stereotactic radiotherapy requires a high degree of precision because of the very steep dose gradients used, typically decreasing the dose from 90% to 50% in 16 2 mm. The localization errors from several millimeters to a centimeter in routine clinical applications are inappropriate for this type of treatment, especially due to the close proximity of critical structures. For example, routine radiation treatment with a field size of 15 cm x 15 cm and a +2 mm error margin is mostly possible and will not greatly increase the complication rate in most body sites. However, with a small treatment field of 2 cm x 2 cm and a +2 mm error margin near the critical structure such as the brain stem, this error might cause complications simply because the penumbra region of the beam profile might cover part of the critical organ. Fractionated stereotactic radiotherapy is initiated to achieve two goals: delivering the highest accumulative tumor dose to the target area throughout the treatment course and minimizing the beam portal errors by accurate localization of the tumor with respect to the treatment beams. Refixation Techniques To perform fractionated stereotactic radiotherapy, several fixation techniques had to be developed. The easiest method is to use a special head frame or mask for immobilizing the patient's daily treatment position. Lyman et al. [Lym89] developed a mask and stereotactic frame combination technique for charged particle radiation treatment. Usually in clinical environments facial or bony landmarks are matched to correct possible relocalization error. The mask has to be made giving particular attention to the bony structures such as the mandible, zygomatic arches, frontal ridges, and nose. The mask is divided in two halves coronally and requires careful molding procedures. As Lyman discovered, the mask is fit over the skin and hair of the head and relies on friction for immobilization; it is found that some slippage can occur for large rotations 17 around the localization position because the proton beam can only be delivered at a fixed angle. Thus, the port angles of beams will be limited in order to prevent the large rotational errors of the patient. The immobilization of the patient's head in a reproducible manner within the stereotactic frame is based on the visualization of landmarks, for instance, nose, chin, and frontal processes. The proper fit of the mask will be under continuous evaluation from procedure, to procedure which might require extra time for patient immobilization [Lym89]. With the system, they reported an error of approximately 1 mm in repeated sessions in each of three orthogonal directions could be achieved. This is only the repositioning error, not the total treatment error. Another factor influencing patient position is the potential for patient movement during approximately 1 hour of diagnostic procedures. This effect will distort the coordinates that were registered in the images. For stereotactic radiosurgery, which needs high treatment accuracy, using only a mask without other supporting devices seems to have some potential problems. Hariz and his coworkers [Har90] designed the Laitinen stereoadapter for fractionated radiotherapy. The noninvasive stereoadapter shows a high degree of reproducibility between repeated mountings. The total maximum error between isocentrical positioning of the target does not exceed 3 mm. Delannes et al. [Del91] uses Laitinen's stereoadapter to immobilize the patient's head. This noninvasive head frame depends on the patient's facial and head structures. A 1 mm spatial coordinate precision of reproducibility is reported. Gill et al. in UK use a noninvasive relocatable frame to localize patient position by mechanically and 18 physically constraining the patient skull for stereotactic external beam therapy [Gil91]. This frame is also employed for fractionated radiotherapy, with overall mean displacements of 1.2 mm for AP and 1.8 mm for lateral directions [Gra91]. These measurements are repositioning inaccuracies; they do not include any localization errors. However, the reproducibility of the positioning of the head for each of these procedures depends on the closeness of the fit of the stereoadapter, the uniqueness of each patient's cranial and facial features, and the cooperation of the patient. The amount of soft tissue and the uniqueness of cranial features are related and greatly affect the reproducibility of positioning [Lym89]. Jones [Jon90] applied three implanted gold markers to define the reference coordinate system, which promoted fractionated therapy. The total potential misalignment of this method in all directions could be up to 2 mm. Relocalization Methods Pelizzari and Chen [Pel87,89,91] developed a method to allow accurate registration of different image modalities (CT, MR, PET). Utilizing a nonlinear least squares search, the "hat" image (image to be transformed) is translated and rotated to minimize the differences from the "head" image (reference image). This is the socalled "Hat" and "Head" method. From 3D reconstructions of the patient's face and skull, data points that are far apart can be used to match these two rigid bodies. An error range of 1 mm is generated by careful manipulation of the digitizer system. This method could be applied to reposition the patient's head for stereotactic radiotherapy. Using a three dimensional digitizer, the relative coordinates of the new patient position can be obtained with respect to the reference position. Thus the localizer marks on a patient 19 head can be realigned to match treatment room laser planes for gantry and treatment couch setups. Kato et al. [Kat91] have developed a frameless, armless neurological system to localize the target inside the brain. They used a 3D digitizer with lowfrequency magnetic field to detect the 3D coordinates and the three orientation angles. Mean positioning errors from 1.7 mm up to 4.0 mm were reported when the digitized probe was tested. However, when the magnetic digitizer system interferes with metals, coordinate shifts as much as several centimeters were observed. This system can locate the intracranial target with careful manipulation. In comparison with the conventional head frame systems, the accuracy is still under investigation. Barnett et al. [Bar93] utilized a frameless, armless stereotactic localization system in brain tumor surgery. The localization system consists of a localizing wand that emits ultrasonic pulses that can be detected by a tablemounted array of microphones. With the triangulation of the emitter signals, the 3D localization of the wand tip can be determined with respect to the microphone detector arrays. In the operating room environment, using four spark pairs, the reproducibility of a point in space is +0.6 mm standard deviation. The mean error for measuring distance within a 1000 cm cube is 1.1 mm + 1.0 mm, which represents 1.0% + 0.7% in space. However, for the overall determination of accuracy, the mean linear error in localizing 19 targets in 48 patients is 4.8 mm + 2.1 mm standard deviation. This error term refers to the image data of 3 mm thickness computed tomograph or 2 mm voxel magnetic resonance imaging volume acquisition. 20 As previously described, researchers in other clinics used landmarks to properly reposition the patient's treatment position. Those landmark positions were retrieved from visual localization or by using plain radiographs. Errors can be analyzed by using radiographs or a portal image system [Mee90]. All the fractionated treatment courses are based upon the error estimation of the treatment position. Therefore, the relocalization procedure becomes a critical issue to realign the tumor for fractionated treatment. All the fractionated radiotherapy methods have potential patient repositioning errors. Errors from millimeters to centimeters have been discussed previously and are based on the observed results. Those methods need to pay particular attention to bony landmarkers; experienced manipulation is required to correct possible errors. The drawbacks of those repositioning systems are prolonged procedures which introduce potential patient movements. For high precision fractionated stereotactic radiotherapy, a fast, accurate, and easy technique would be the ideal alternative. CHAPTER 3 THE UF STEREOTACTIC SYSTEMS The UF Stereotactic Radiosurgery System Hardware Characteristics A modified radiosurgery system of the BRWm stereotactic ring (a conventional neurosurgery device) has been designed at the University of Florida and commercialized as the Philips SRS 200 radiosurgery system [Bov90a]. It is a three gimbal bearing system; the head stand and bearingholder system are incorporated into a single portable piece of equipment (Fig 3.1). Instead of attaching the BRW ring to the treatment couch, it is attached to the head stand, which delivers better accuracy in the patient's treatment position. This head stand provides high accuracy micrometer adjustment in positioning the localized tumor target center to coincide with the system isocenter. The swing arm system provides rigid rotation and support of the auxiliary collimators. The auxiliary collimators are to define the treatment beam area and minimize the penumbra effects. The set of circular collimators ranges from 5 mm to 40 mm in diameter and with proper divergence delivers superior treatment volume selection and optimization. When the heavy LINAC gantry head rotates towards the horizontal treatment position, the gantry sags and tilts, displacing the central beam below the isocenter. A gimbal bearing system is designed to decouple the torque of the gantry from the auxiliary 21 Figure 31: University of Florida Stereotactic Radiosurgery System ([Fri92a], page 834, used with author's permission) collimator and ensures that the sag in the LINAC does not result in angular deviation of the collimated beam from the target point [Fig 32, Fri92a]. Treatment Planning Procedures The stereotactic radiosurgery procedures are currently performed on an outpatient basis. Before starting the diagnostic procedures (Computed Tomography, CT; angiography; or Magnetic Resonance, MR), the BRW ring has to be fixed to the patient's head under local anaesthesia. CT is the standard imaging modality for target localization; if uncertainty exists, angiograph or MR are utilized to visualize the target shape. All the localization procedures refer to the BRW coordinates, which are fixed with respect to the patient's internal and external anatomy. After finishing the imaging sequence (CT or MR), these images are transferred to the treatment planning system (Sun" Workstation) and the localization process starts. If an angiograph has been performed to locate the target, the biplane films with sixteen fiducial markers (eight per I I .... ____, Figure 32: LINAC gantry sag was taken care of by the gimbal bearing (A) system, therefore the mechanical integrity of 0.2 mm + 0.1 mm is maintained. ([Fri92], page 146, used under the author's permission) AP/lateral projection) are used to define the reference coordinates [Sid87]. Then the radiopaque markers are digitized into the treatment planning computer for the calculation of a geometrical center and a center of mass that closely matches the nidus [Bov91]. The best matching pair should be used as the center point of target. A stateoftheart treatment planning software was developed to perform rapid dose calculation and to display isodose distributions. A standard nine arc plan, shown in Table 31, is the starting plan for the final optimized plan. Varying such factors as 24 table angles, collimator sizes, or beam weighting parameters changes the isodose distribution on three orthogonal views. Multiple isocenters are specified if the target is extended or irregularly shaped. The visual optimization of the treatment plan is considered a key factor in obtaining the best plan, since it has to take sophisticated planning interaction into consideration. Patient Treatment After the best treatment plan has been obtained and the adequate treatment information is collected, the stereotactic radiosurgery accessories are then attached to the LINAC and the target isocenter is set on the head stand. Setting up the UF radiosurgery system takes less than 15 minutes. An independent phantom is then set up to match the isocenter position, and xray images of the phantom target locations are taken at various table and gantry combinations [Win88]. If the images show target position in the center of collimator field within +0.2 mm, the head stand is considered at the right treatment position. The patient is brought to the room, and the BRW halo is attached to the head stand. Treatment then proceeds as defined in the treatment planning procedures. At the end of treatment, the BRW ring is removed and patient is free to leave. The UF Fractionated Stereotactic Radiotherapy System Hardware Characteristics The University of Florida stereotactic radiotherapy system has the following goals: Table 31: Standard nine arc treatment plan Arc Collimator Table Gantry Gantry Arc Number Size Angle Start Angle Stop Angle Weight 1 10 mm 10 30 130 1 2 10 mm 300 30 130 1 3 10 mm 50 30 130 1 4 10 mm 700 300 130 1 5 10 mm 3500 230 330 1 6 10 mm 3300 230 3300 1 7 10 mm 3100 2300 3300 1 8 10 mm 2900 2300 330 1 9 10 mm 2700 230 3300 1 26 (1) It should be frameless. The traditional BRW ring system in the LINAC radiosurgery provides rigid coordinate transformations that bring the tumor center to the machine isocenter. The BRW coordinate system is not valid for fractionated treatment if the ring is removed. Instead of the use of invasive pins on the patient's skull, the patient will be held in position with the aid of four soft rubber head posts and soft pressure pad systems. To minimize the rotational error corrections, the diagnostic procedures (CT, MR, and angiography) will be carried out with the patient positioned in the same immobilization device that will be used for treatment. The head holder for CT and MR has a convergent set of rods embedded into the base and provides a method whereby the axial increment of the CT or MR coordinate can be verified by means other than the scanner readouts. The head holder also allows the attachment of two orthogonal planes with fiducial markers, permitting the reconstruction of the target position as described by Sidden et al. [Sid87]. (2) It should achieve the highest treatment accuracy. In order to obtain the highest accuracy of the UF refixation technique, this device will be incorporated into the existing UF radiosurgery system. Since the UF radiosurgery system provides the best mechanical treatment accuracy among the LINACbased radiosurgery systems, this add on refixation device will deliver the highest treatment accuracy. (3) It should correlate to the tumor shape as closely as possible. Experience shows that many tumors are irregularly shaped. By simply translating the patient with respect to the laser or couch, readings for fractionated therapy may cause target misalignment either in a linear direction or in orientation. To overcome the problem and 27 to correlate the tumor shape as close as possible for different treatment courses, a six degreeoffreedom parameter estimation process is necessary to realign the patient for correct treatment position and orientation. Algorithms for patient realignment with translation and rotation operations have been developed as part of this work. This software provides keys to the treatment of patients with intracranial abnormalities. A mechanical device will be created to provide the necessary mechanism for repositioning purpose. Therefore, not only the isocenter treatment position can be properly aligned but also the target shape can be matched with the minimum normal tissue exposed to the high precision, high isodose gradient treatment. (4) It should be easy to use. Since a lengthy treatment preparation procedure may cause potential patient discomfort and movement, this system should be easy to adapt to the UF radiosurgery system and should not be complicated to operate. Setup time of this device should be short so that the patient can tolerate the treatment without difficulty. Daily setup should also be easy for the technologist to achieve optimum treatment and to do so economically. Localization and Refixation Procedures The procedure of frameless repeat fixation and treatment starts with a dental bite plate system. This custom bite plate will be fabricated by a prosthodontist with three fiducial markers in place for later verification purpose. This procedure is performed prior to the diagnostic scans. After the bite plate system is placed in the correct position, the patient will go through the CT or MR scanning process. Each image not only contains the 3D information of the patient's anatomy but also the coordinates of the fiducial markers relative to the 3D image data sets. After the patient scanning procedures are done, a similar treatment planning session as described in the UF radiosurgery papers [Fri89a, Bov90] is performed to calculate dose distribution of multiple noncoplanar arcs. The dose to the target is then maximized while the dose to the normal tissues and critical structures is minimized. The position of isocenter(s) could be obtained as well as those of fiducial markers. The patient is then positioned at the LINAC with the UF radiosurgery system for strict isocenter control. The patient is immobilized with the soft pad headclamping system which allows little or no head movement during the treatment procedure. This immobilization device should not interfere with any fiducial markers and anatomical reference points. Then the fiducial markers are digitized with a 3D camera digitizing system which is well calibrated with respect to the machine isocenter. Software programs are executed to obtain translation vectors and rotational angles required to realign the patient's treatment position. A 3D translation and threeaxis rotation system is designed and adapted into the existing University of Florida radiosurgery system as a subsystem to perform the patient realignment function. Once the patient is properly repositioned according to the calculation results, the modified 3D image correlation program which was originally initiated by Pelizzari/Chen [Pel91] at the University of Chicago or the new programs developed at the University of Florida can be used to compute and ensure the final treatment coordinates by using a portable personal computer. After this verification procedure is finished, treatment with a regular stereotactic treatment plan may proceed. Patient Treatment After the patient treatment position is "fine tuned" to the correct coordinates, a verification system consisting of three spheres will be employed. The current location of these spheres, which are incorporated into the basic CT/MR data sets are then verified prior to treatment. The bite plate system differs from others used in that it is not used for repositioning and does not have any realignment pressures applied, thereby maintaining its true geometry relative to the patient's anatomy. The fiducial markers are then digitized again to ensure that the final coordinates are valid. The patient treatment position actually can be verified at any time when there is suspected patient movement. If minor errors exist, a decision can be made as to the adjustment or relocalization of the treatment target position. The daily treatment will start and the radiation is then delivered. Since the final position of the fiducial markers on the bite plate system has the greatest sensitivity on the treatment position errors, this verification system will certainly provide the greatest treatment accuracy. A plain film procedure can also be utilized to verify the final treatment accuracy if necessary. This procedure would provide yet another independent technique for the fractionated treatment schemes. CHAPTER 4 TARGET RELOCALIZATION BY 3D CONTOUR POINTS Statements of Problems Radiosurgery can achieve extreme accuracy by immobilizing patient's position with strong fixation (i.e. stainless steel pins into the skull). Single fraction treatments are designed to generate maximum damage to the target with high precision. However, they cannot take advantage of the maximum benefits of fractionated treatments. Some tumors, such as a large pituitary tumors or craniopharyngiomas, are more sensitive to radiation than are normal brain tissues when multifractions of radiation are delivered. Therefore, fractionated, stereotactic radiosurgery would be the ideal treatment choice. While treating intracranial tumors, traditional radiosurgical treatment might place the optical nerve at risk, as tumors may directly abut that structure. Utilization of routine radiation therapy techniques would involve large fields, encompassing much more normal brain tissue than would utilization of radiosurgery techniques. An ideal compromise would be the fractionated, stereotactic technique [Spi91]. The traditional BRW stereotactic frame and other fixation frames fix directly to the skull by means of pins and anchor posts. Several new stereotactic frames have been developed which locate to the bony anatomy and are acceptable for repeat fixation for fractionated therapy [Gil91, Har90, Hou91, Jon90, Lym89]. To achieve better tumordose response and to 31 reduce dose to the normal brain tissues, the fractionated technique is very important for tumor treatments other than AVM. The basic principle for fractionated treatment is to reposition the target isocenter back to the machine isocenter for each treatment fraction. Registration of the patient's contours or the fiducial markers is an important step in realigning for repositioning. Most systems perform repositioning only by repeating the treatment settings from laser or treatment couch readings [Gil91, Hei89]. Due to limited accuracy of the devices, this might put some of the critical structures at risk. In order to fully realign the treatment position, not only is the translation back to the isocenter necessary but also the 3D orientation of the tumor should be matched. Therefore, with six degreeoffreedom fitting procedures (translation and rotation) instead of using three degreeoffreedom parameters (translation only), the tumor shape can be matched as closely as possible. Correlation of two sets of data points is a manytoone mapping problem mathematically. Calculation of translational and rotational parameters is done by solving nonlinear equations. Theoretically, no analytical solution is available for these equations. A leastsquared approach is the favored choice for solving this kind of problem. Since all the 3D data points contain certain types of digitizing and systematic errors, finding the best fit of these two different 3D data sets is also equivalent to minimizing the error terms. Verification of the treatment position is also an important factor in identifying the correct treatment geometry. A separate system should be supplied to perform independent treatment verification. Error analysis of the 3D data sets and the target position and orientation should also be performed. Methodology The methodology for patient repositioning at the University of Florida is based upon the calculation of six degreeoffreedom parameters of two registration data sets which come from different fractionated schedules. Mapping of those anatomical points and contour lines becomes the key issue in calculating the patient's realigned treatment location. Most of the researchers only perform the translation function to realign the patients. According to the references [Del91, Gil91, Har90, Har85, Hei89, Jon90, Pel 91, Tho91], the treatment precision ranges from 1 mm to 5 mm which depends upon the accuracies of laser and treatment couch, the patient's motion, and the operator's skill. In order to correlate the "true" tumor position and orientation, rotational parameters are necessary to fully realign the target, especially in the 3D orientation of the irregular target shape. Kessler et al. [Kes91] has described four different methodologies to implement the calculation of transformation matrix for two different 3D data sets. They are (a) point matching, (b) line matching, (c) surface matching, and (d) interactive matching. In the point matching method, the transformation parameters are determined by establishing a onetoone response between a set of point landmarkers visible in both fractionated treatments. These points can be either internal anatomical landmarkers or externally placed fiducials. Internal landmarkers include structures such as optical nerve foramen, the center of orbits, and the odontoid process, etc. External markers depend on the image modalities and selected fiducials. For the study, the external markers are 33 arranged in a fixed pattern for repeat localization. In situations where one of the image studies covers only a limited extent of the cranium, it is difficult to place fiducial markers to ensure the visibility. Line pattern markers are suitable for the continuous scan which can provide the location on each image data set. The line markers can also be marked at a specified position which correlates the image coordinates. Surface matching method uses the edgethresholding algorithm to extract the outer contour lines from the geometrical feature of patient anatomy. Since the patient head is a rigid structure, the contour information should be invariant over time. The head model is created by forming a pseudosolid model through the triangulation of the contour lines. The transformation is solved by minimizing the mismatch between two surfaces, which is so called the head/hat method by Pelizzari/Chen [Pel89, 91]. The transformation matrix is achieved by manipulating the hat (by varying the transformation parameters) to properly fit (by minimizing the mismatch) onto the head model. This method provides six degreeoffreedom fitting of two data sets. However, the precision for minimizing the mismatch is affected by the inherent tolerance of the hardware system and the computing work load is heavy. The complexity of determining the transformation parameters is proportional to the multiplication of the number of triangles formed on the head model and the data points on the hat model. In the interactive matching technique, the transformation parameters relating two 3D data sets are determined manually using a graphics utility. The method also needs a relatively high computing power to assure that the information can be retrieved and displayed during a short period of time. 34 We intend to employ both the anatomical points and the patient surface contour line to perform and estimate the mapping function. A vertical facial contour line of patient midplane is selected as well as the anatomical points are recorded to form a complete data set for patient texture information. This data set will be used to correlate the corresponding points or line with proper identification on CT/MR images. In order to correctly set the treatment coordinates for a patient's fractionated treatment, several algorithms have been developed to identify the six degreeoffreedom fitting parameters. These algorithms will be compared with the head/hat concept and identify the validity of the calculations. The following approach explains the procedure for completing the calculation task in finding transformation parameters. First the patient's midline is digitized as one reference data set; secondly several obvious anatomical points are also digitized as part of the mapping data set. Those data points are integrated as one 3D reference set of points, which represent the 3D structure with respect to the patient's own coordinates (this coordinate system will be the same as that of the LINAC treatment coordinates). Fig. 41 shows the conceptual drawing of registration of the 3D contour points. Since the CT/MR data sets have finite resolution of the pixel size and image scanning slice thickness, the data sets always contain a discontinuous step function for each pixel value. Therefore, a spline fitting function is necessary and is incorporated to resample the midline contour for data mapping purpose. A piecewise cubic spline fitting function is applied to obtain the smooth and continuous curve for the patient's midline. The spline fitting procedure will be considered as follows [DeB78, Wol92]: the 2D PATIENT IMNOBILIZATION DURING TREATMENT A RIGHT EAR TRAGUS B RIGHT EYE LATERAL COMMISSURE C RIGHT EYE MEDIAL COIMISSURE D CONTOUR T.INE SOFT PRESSURE PAD Figure 41: The conceptual drawing of patient contour and anatomical points registration process (courtesy of Dr. Bova) spline fitting procedure is first selected because of the inherent simplicity in data fitting process. Given a 2D point sets of (xk, y) for 0 spline consists of n1 cubic polynomials. They pass through the supplied points, which are known as the controlpoints. We now derive the piecewise interpolating polynomials. The ku polynomial piece, fk, is defined to pass through two consecutive input points on 9 'e % 36 the ranges [xk, xk+,]. Furthermore, f, is joined at Xk (k=l,...,n2) such that fk (first derivative of the function) and f"k (second derivative of the function) are continuous. The interpolating polynomial fk is then given as fk (x) =a3 (xxk) 3+a2(xxk) 2+a (xxk) +a (41) The four coefficients of f, can be defined in terms of data points and their first and second derivatives. If the end point derivatives are not known, the "notaknot" condition is assumed to imply that the first and the last interior point knots are not active for calculation [DeB78]. These two knots have the same derivatives as those of their neighboring knots. By applying the spline fitting function a nice and smooth contour curve can be generated to truly represent a patient's midline structure. The 2D concept of a spline can be extended into 3D coordinates. The contour data points in the CT/MR images can be resampled again as the comparative standard for fractionated treatments. The logical steps to form the data set of anatomical points and contour lines are as Fig. 42. A set of computer programs were written to perform data sorting, resampling and reweighting of the initial 3D data set. After the 3D reference data set is saved as well as the treatment planning information, a separate 3D data set which follows the same contour lines and anatomical points is acquired from the 3D digitizer for mapping. Six degreeoffreedom fitting parameters are then calculated and compared with different algorithms in order to find out the best correlation of these two 3D data sets. A real time 3D digitizer is incorporated into the data acquiring system for the fractionated treatment schedules. This system is well calibrated with respect the LINAC Figure 42: Sequence for generating the reference CT/MR data sets coordinate system. Because of the high resolution and precision, the data sets from this camera system provide reliable information which has the minimum perturbation of the system accuracy. A flow diagram sketched in Fig. 43 is presented to integrate various steps for the 3D data correlation process. An immobilization device is also created to maintain patient's treatment position while the radiation is being delivered. This add on device has the function to tune in the correct coordinates for the fractionated treatments, remembering that this device cannot interfere with the treatment fields. Therefore, the dosimetry considerations will be the same as that in single fraction radiosurgery. After the patient has been properly Sort the CT/MR image data sets into increasing axial coordinates Spline fitting of the data sets Weigh the data sets for both anatomical points and contour lines Figure 42: Procedures for the calculation of 6D fitting parameters of the Fractionated Stereotactic Radiotherapy realigned and immobilized, a separate verification bite plate system consisting of three spheres will be employed. The correct location of these spheres were previously incorporated into the basic CT/MR imaging data sets prior to treatment. This system differs from the other systems used is that it is not utilized for repositioning and does not have any realignment pressures applied. Therefore, the true geometry relative to the patient's anatomy (i.e., tumor location) can be maintained. Once the patient's treatment position and orientation are verified as desired, fractionated treatment can start by using the calculated treatment planning information. CHAPTER 5 3D COORDINATE MEASURING SYSTEMS Laser Tracking Interferometer System The laser interferometer uses light interference principles to precisely measure the linear displacement or velocity of a body. A laser produces a beam of coherent, monochromatic light that is passed through a beam splitter. Part of the beam is reflected towards a fixed mirror and the other part of the beam is transmitted through the beam splitter toward a moving mirror. The light from both the fixed and moving mirrors is reflected back toward the beam splitter, which is designed to recombine the beam as illustrated in Fig. 51. The photodetector of the detection system will sense an alternating intensity. One cycle of intensity change represents a mirror displacement of onehalf a wavelength of the light from the laser. The accuracy of the interferometer is affected by any influences that alter the wavelength of light. Therefore, any errors caused by the interference should be properly corrected or compensated. If properly constructed, the laser interferometer could be used to measure 3D coordinates in space. Systems that combine laser interferometers, optics and electromechanical devices to perform coordinate measurements have been developed by several research centers [Bro86, Lau85, Zhu92]. They use both length and angle measurements to determine the target location. Relative distance measurements provided Fixed Mirror  Figure 51: Basic layout of a laser interferometer by laser interferometers have extremely high resolution. However, several factors should be considered for coordinate measurement. First, the interferometer readings only provide relative displacement information; there should exist one reference system in which relative displacements are measured. Second, the direction of the incident beam must be determined. Third, since it is difficult to position the mirror center, the mirror center offset, which is the distance from the mirror center to the intersection point of the incident beam with the mirror surface, should be properly compensated. However, the accuracy of a coordinate measuring machine based on a laser tracking system is primarily determined by the geometry of the tracking mirror system and its control. The major error sources are gimbal axis misalignment and mirror center offset. There are two primary reasons for not using the laser interferometer as the first choice to measure the object coordinates. The first is that it is difficult to establish a reference position. Since the interferometer collects the information of the reflected laser Beam Splitter ,. Photodetector Output Signal 7 I h Laser e L yt 2CLCU Moving Mirror 41 beam, the measurements have to start again if the light beam is interrupted at any time. The second reason is that the interferometer is expensive. A typical commercial system with necessary optics costs more than $30,000. This cost does not include the calibration subsystem. Camera System Camera systems have been widely utilized to capture scenes and perform image processing with various applications. In order to apply the systems in 3D coordinate measurements, the concept of stereo vision is applied. Mapping a 3D scene onto an image plane is a manytoone transformation. That is, an image point does not uniquely determine the location of a corresponding world point. The stereo vision technique could retrieve the missing depth information which provides 3D coordinate information. Stereo vision involves obtaining two separate imaging views of an object of interest. In Fig. 52, the distance between the centers of two lenses is the baseline (B). Two image points (xl,yl) and (x2,y2) are given to calculate the world coordinate (X,Y,Z) of the corresponding point. If the cameras are identical and perfectly aligned, the 3D information could be retrieved from the stereo vision system. A typical vision system consists of a lens and a light sensitive array. Light from the image is then focused on the array by the lens. The light sensitive array consists of a large number of discrete cells that sense the frequency of light which strikes them. These cells are referred to as pixels. The frequency information from each pixel is digitized and a number is assigned to the sensed frequency. Using the information from each cell, the image could be reconstructed and analyzed. For experimental purposes, image 1 Figure 52: Model of the stereo imaging process a program has been written to verify the accuracy of the stereo vision concept of the camera system from the film dosimetry setup. This program utilizes the functions imbedded in the film dosimetry project [Dub93] to execute the basic image grabbing operations. From those acquired images, a set of dots are analyzed to calculate the target location from the theory which was developed by Siddon [Sid87] to identify the BRW' coordinates in plain radiographs. This camera system uses a Data Translation DT2851 high resolution frame grabber board and a DT2858 frame processor board to capture and process images. A lucite plate is made as an calibration device in order to calculate the accuracy and resolution of the imaging system. It has four dots which are separated at 6 cm intervals and one dot on the cross hair of the diagonal lines. In order to examine (x1,yl) Image 2 (x2,y2) world point optical axis 43 the concept of stereo vision, a calibration experiment of the image resolution was performed as shown in Fig. 53. First, the lucite plate was positioned on top of the light box and those dots were easily captured from the frame grabber into the image plane. Second, in order to calibrate this camera system with respect to a specified physical dimension, a calibration factor which represents the actual dimension of each pixel (millimeter per pixel) is entered. This process helps to define the image resolution according to the user's selection. Third, a Sobel image processing filtering technique [Gon87] was performed to accurately identify the geometrical centers of the dots. The external video input signal is supplied by a CIDTEC CID2710 Charge Injected Device (CID) camera. The CID2710 sensor chip consists of a high resolution matrix with 776 horizontal and 512 vertical pixel elements. The active area is approximately 60 mm2. A 512x480 pixel image is stored in the buffer memory on the image processing board, and displayed on the monitor. The center point of the calibration plate is calculated with respect to the four corner points. The error magnitude of the X and Y coordinates of the center point is an integer of the pixel resolution. Therefore, the accuracy of coordinates from those fiducial points depends on the pixel resolution. The resolution in the angio localization procedure of the UF stereotactic radiosurgery system is about 0.0254 mm. If the camera system had higher resolution than that of the existing angio localization procedure, the stereo vision concept will be available for the determination of coordinates of the intracranial target. Calibration Plate Figure 53: Film dosimetry system setup (camera, light box, and stand). A calibration plate is designed and analyzed for the system resolution. Table 51: Pixel accuracy estimation of different calibration factors Calibration factor Pixel accuracy of x Pixel accuracy of y (mm/pixel) direction direction x cal' = 0.005 mm 0.297776 mm 0.233384 mm y_cal' = 0.005 mm x cal = 0.05 mm 0.297019 mm 0.233392 mm y_cal = 0.05 mm x_cal = 0.5 mm 0.296985 mm 0.232635 mm y_cal = 0.5 mm *: where x_cal and y_cal represent the user defined pixel resolution. 46 After analyzing the images, the Table 51 shows the result of pixel accuracy of the camera system. The resolution determines the positions of the fiducial markers. Therefore, they will affect the solution of the determination of target coordinates. From Table 51 a noticeable characteristic is observed that the accuracy of each pixel actually is determined by the CID sensor chip where the image was formed. User defined pixel resolution has no effect on the camera system. For a different calibration distance of the lucite plate, the pixel accuracy of the image is fixed with any changes of calibration factors. Therefore, the stereo vision method for the experimental purpose is limited by the resolution of the image chip unless we purchase high accuracy photosensor chip, for example, 1024x1024 array instead of 512x512 array. The resolution of the array plays an important role in measuring the 3D coordinates. For example, assume that an array consists of 512 rows and 512 columns. If the lens system is focused on a 152.4 cm (5 ft) by 152.4 cm (5 ft) workspace, each pixel would be averaging the light from 0.3 cm by 0.3 cm square. If a resolution of 0.05 mm is required for a given identification procedure and the vision system to be used has an 512 by 512 CID array, the measurement volume will be limited to no more than a few centimeters. Therefore, it would be impractical using this passive imaging system to obtain the 3D coordinates from the fiducial markers on the images. Only the active imaging method (i.e. with an active light or laser source at the measuring position) can accurately locate the target position if the camera system were selected. The reason is that for the passive camera method the image quality is responsible for the determination of the 3D coordinate of the selected points. However, extra image correlation and processing efforts are always 47 necessary to extract the coordinate information. Uncertainty of the definition for matching points on the image and localization errors might occur while performing the image correlation procedure. The stereo vision concept has also been verified by Fitzgerald et al. [Fit75] by using two stereo shifted radiographs instead of using cameras. Xray radiographs theoretically carry infinite resolution. The shifted radiograph technique is using a fiducial jig allows the shift parameters to be derived from information contained in the films. By translating the Xray linearly, the depth information of the object can be calculated. However, there are potential errors when small shifts are used. This method explains another possibility of calculation error if the two cameras are not separated far enough. Besides the limited resolution of camera sensors, the geometrical structure of camera setups also play an important role in 3D coordinate calculations. Several commercial companies such as Pixsys", Inc. and Northern Digital", Inc. have successfully applied the stereo vision concept with an active light source to manufacture 3D optical digitizer products. In this commercial system several charge coupled device (CCD) cameras are positioned at different directions to capture images from any angle since only two cameras might miss some of the views. Northern Digital", Inc. manufactures a camera type of measuring system which is called OPTOTRAK/3000 [Nor93]. This system has the capability of using a photosensitive detector that is sensitive to light in the infrared range. The sensor is capable of accurately determining the location of the center of a spot of light that has been projected onto its surface. The detector is an analog rather than digital device so the resolution is 48 theoretically infinite. The device is also highly linear, which enhances the overall accuracy of the total system. These systems operate by the following procedures. Several infrared Light Emitting Diodes (LEDs) are attached to the points for measurement. Then the system controller switches each LED on in sequence and the light from the LED is projected through the lens onto the photodiode. The XY coordinates on the image plane of the detector are proportional to the XY location of the LEDs in the field of view. Each LED could be located in approximately 50 /sec so the system can record the trajectory of those LEDs. Therefore, those cameras can track and report the 3D coordinate location of the LED. Those data sets can be transferred into the workstation or PC for 3D graphics and for image correlation. According to the manufacturer's specification, they have 0.01 mm resolution at 2.5 m, with an accuracy of 0.1 mm in the XY direction, and 0.15 mm in depth which is perpendicular to the XY plane. Pixsys" also utilizes a array of CCD cameras which are placed above the patient table. These cameras track and report the 3D coordinates of tiny LEDs which are built into the instruments such as stylus, forceps, biopsy needles or ultrasound transducers. The 3D coordinates are reported up to 100 times per second. If the CT or MR data sets are preloaded to the computer, then the corresponding images from these data sets can be displayed from the locations of LEDs which shows the instruments' location relative to anatomical landmarks. However, camera resolution is 0.3 mm in one cubic meter measuring volume and 0.6 mm in two cubic meter measuring volume, which is worse than the Northern Digital camera system. Articulated Arm System Articulated arm systems for 3D coordinate measurements have been widely used in aerodynamics or machinery. These devices are manufactured to ensure precise motion along the desired axes and they are instrumented to determine the joint displacement and orientation to a high degree of accuracy. The basic structure of this type of mechanical system consists of joints and arms; a set of potentiometers or encoders are also necessary to read out the spatial position and orientation. Such devices typically have a resolution of 0.025 mm and a total accuracy on the order of 0.125 mm [Moo91]. Articulated type of coordinate measurement systems are typically designed for inspection of parts and assemblies. In the clinical uses, such coordinate measurement systems are quite cumbersome because the heavy weight and geometrical structures often interfere with either the diagnostic or therapy procedures. Watanabe et al. [Wat87] utilized a 3D digitizer for CT stereotaxic surgery. This digitizer has multiple aluminum arms and joints which consist of high resolution potentiometers. Analog signals from the potentiometer are linear with respect to the angle of the joint. The 3D coordinates of the tip of the arm are calculated from the joint angles and each segment length. Another group in Germany developed a robot type of system in house to perform computer assisted surgery. Adams et al. [Ada90] developed an electromechanical robot for 3D coordinate measurements. This digitizer can interact with CT image data sets and perform a realtime surgery. The system has an accuracy of +1 mm which is typical in robotic standards. All the articulated type of systems have the following characteristics: (1) Heavy and fixed mechanical systems, cumbersome in clinical settings; (2) Lengthy time for measuring a large number of contour points; A company named Shooting Star Technology" had developed a inexpensive 6D (translation plus rotation) tracking localization system ADL1". However, the resolution is at the order of 0.635 mm for displacement and 0.3 degrees for orientation, which depends on the position in work volume. Apparently it is not suitable for high precision clinical uses. Frequency Matching System Instead of using mechanical arms and joints, PolhemusT, Inc. developed a real time, 3D data generation and manipulation system for motion tracking. This system utilizes a radiofrequency transmitter and receiver for emitting and sensing the electromagnetic signals. From the geometrical relationship, the controller automatically sends out and records the 3D information. Because the system can track objects in real time, if the sensor is attached to a stylus, the 3D coordinate information can express the object outlines which are being traced. The accuracy of this digitizing system ranges from 3.3 mm (0.13 inches) RMS (Root Mean Square) from 101.6 mm (4 inches) to 381 mm (15 inches) sourcetosensor separation. As the separation increases, the positional accuracy degrades linearly. Although this system has reasonable accuracy for measuring 3D coordinates, the drawback is that large metallic objects, such as desks or cabinets, located near the source or the sensor adversely affect the performance of the system. If the patient is positioned under the treatment machine, the metallic parts and the electromagnetic pulses will 51 interfere with the proper reading of the 3D coordinates. Therefore, this system is excluded for our application. Ultrasonic System The ultrasonic system for determination of 3D coordinates of defined points is also clinically available in neurosurgical operations. The operational principles of the ultrasonic navigation system are similar to those of the frequency matching system. There are several commercial products which perform the rangefinding function as well as calculate the spatial coordinates of the selected points. Science Accessories Corp. (Southport, CT) developed the Model GP83DT ultrasonic rangefinding device as a 3D navigation system. This range finder is computer controlled and calculates distances between sound emitters and microphones by measuring the transit time of an acoustic pulse traveling from a sound emitter to the microphones. The ultrasonic system was examined with respect to all the possible errors, either systematically and randomly, to quantify the possibility of clinical use. Friets et al. [Fri89b] studied the different errors in this ultrasonic system. The range finder error is 0.5 mm + 0.2 mm with a 0.2 mm standard deviation. The average determination error of actual microphone separation is 0.2 mm + 0.1 mm with 0.1 mm standard deviation. Transformation of the emitter signals to the image plane introduces up to 4.3 mm + 2.3 mm with a standard deviation for the orientation point (the furthest point from the focal plane of the system). The minimum error of the system while locating a point from the real operating environment in CT space is 2.2 mm + 1.0 mm with a standard deviation of 0.5 mm. They claimed that the phantom registration can achieve as low as 2 mm and 3 mm display error of the 52 contours (matrix transformation error). Environmental parameter such as temperature can also influence the accuracy of the ultrasonic system because the medium for carrying the signals is air. Change of the medium characteristics influences the response curve of the ultrasonic system. Clinically, the ultrasonic system has been proven as one of the candidates to perform real time image registration during operations. Conclusion After a careful selection procedure, the Northern DigitalT camera has been selected as the 3D coordinate measuring device for the following reasons: (1) High resolution at certain long distances. For example, it has a resolution of 0.01 mm at 2.5 m distance. (2) Precalibrated 3D coordinate system which eliminates user calibration process. Once the camera is properly calibrated, the results are stable indefinitely. (3) Sampling rate of 3,500 markers per second and 3D data sets available in real time for viewing during data collection. From using the real time digitized probe, contour line data and the selected anatomical points can be obtained by an experienced user within a short time, perhaps a few minutes. Therefore, the calculation of the transformation matrix between the patient scanning image data and the data set from the patient treatment setup location is also possible in the fast and accurate manner. (4) Interface with PC, which makes an easy operating front end in the clinical environment. The required information can be retrieved and computed by simple command line language which needs less human interference. CHAPTER 6 THE FITTING OF 3D REFERENCE DATA SETS Basic Principles While performing motion analysis of a rigid body, the estimation of the 3D motion parameters which are described as rotation and translation is important. This analysis can be very useful in many applications such as scene analysis, motion prediction and trajectory planning in the defense industry. For fractionated stereotactic radiotherapy the patient needs to be set up at the correct treatment position during the multiple treatment courses. Solving the refixation problem requires the matching of 3D surface and anatomical data points on the rigid object at two different times. Therefore, the patient relocalization procedures are similar to that of the estimation of motion parameters. Determining the position and orientation of the head for each treatment fraction is an important aspect of patient relocalization. The mathematical background for solving the refixation parameters of position and orientation is indeed the process of fitting two sets of 3D data points to the patient's contour. Since the patient's head can be considered as a rigid object, the 3D points on the patient's head could provide necessary information for the calculation of repositioning parameters. When the patient experiences the initial stereotactic radiotherapy treatment, the CT/MR scanning data sets and the corresponding coordinates of contour points on the 54 patient's face provide the geometrical relation of tumor location. However, during fractionated radiotherapy, the well calibrated 3D localizer provides coordinate information for those fiducial points on the patient's surface. How to correlate these two sets of data in order to reposition the tumor location to match the machine isocenter without using any invasive relocalization method becomes a key issue to our research. The methodology of the contour match solution is to find the best fit via translation and rotation of an optical 3D data set with a corresponding 3D CT/MR data set. In order to determine the patient relocalization parameters, the following mathematical problem is encountered. The initial diagnostic scans of the patient consist of 3D geometrical relationship between the tumor and patient surface structures. This is a complete standard data set for treatment planning. If patient comes in for the fractionated treatment in another time frame, the interior anatomy can not be observed and only the surface information can be retrieved. Rigid body transformation indicated that if the patient surface information of two treatment time schemes have the best fit and then the tumor should be at the correct treatment coordinates. It is equivalent to find the best correlation parameters for two sets of 3D data. Therefore, assume that we are given two 3D point sets p,=[xi,Yi,zi]t, p'i=[x'i,y'i,z'ji], i= 1,2,3,...,N, (here pi and p'i are considered as 3x1 column matrices.) The two sets of points are related by: p1 = Rpi + T + Ni (6.1) where R is 3x3 rotation matrix, T is a translation vector (3x1 column matrix), and N, is a noise factor. The rotation matrix and translation vectors are answers for best correlation of the data sets but with highly nonlinear functions, which means difficult to solve analytically (no exact answers are available). Huang et al. [Hua86] has shown that rotation and translation can be decoupled by a centroid coincidence theorem. This theorem uses the differentiating technique of Lagrangian multipliers and provides groundwork for the determination of rotation and translation parameters (three from rotation and three from translation) separately. Assume that R and T generate the leastsquare best fit of the following equation: 2= IIp (Rp+T) 112 (6.2) i=1 Therefore by minimizing the above leastsquare term the centroids of {p',} and {Rp,+T} can be proved to be coincident, i.e., the centroid after rotation and translation is the same as that of the leastsquare solution of proper rotation and translation. Mathematically, the centroid coincidence theorem can be shown as follows: p' = p (6.3) n p p (6.4) l= + wphr n ips t p + o (6.5)i n i=i where n is the total number of contour data points. This centroid coincidence theorem implies that we can simplify the original least square problem by translating the point sets first, and the rotation parameters can be solved accordingly. From the above theorem, we clearly understand that if we want to find the solutions of rotation and translation parameters, the following logical procedures should be performed: (1) Simplify the original leastsquare problem by translating the point set {p,} by the following equation: T Ap'p (6.6) where T, is the translation vector. Notice that the T, is not the translation vector which gives the best leastsquared fit answer. It is only the estimation of the centroid difference between two 3D data sets. The accurate translation vector should be calculated after the rotational manipulation, which indicates the rotational sequence is before the translation. The order is important because that matrix multiplication is not commutative. (2) From the centroid coincidence theorem, the coordinates of two centroid points are calculated. Then the 3D points of the two are translated by subtracting the corresponding values from the new origin. (3) Performing the fitting method to acquire the rotation parameters, i.e., find the R vector to minimize the following equation: n E= \ IIQ'RQi\\2 (6.7) i=1 where Qi PiP (6.8) Qip pi (6.9) (4) The translation vector is calculated by utilizing the following relation: T=p'Rp (6.10) (5) The six variables which relate rotation and translation vectors are the best fitting of the leastsquare function of two 3D point data sets. The rationale for using the centroid coincidence theorem is to decouple the translation and rotation operations into two separate parts. Using this technique can reduce the parameter searching space from 6D (rotation and translation) to 3D (rotation only). The calculation time and complexity of the problem is also reduced. In other words, rotational angles are actually the primary unknowns in our mathematical approach since the translation vectors can be calculated in a fast and accurate way after least squared rotational manipulation. In order to solve the leastsquare minimization problem several methods are developed. The selection of different methods totally depends upon personal choice. The detailed description of different approaches are introduced. The six fitting parameters can be calculated to realign patient treatment positions by following the centroid decoupling methodology and the results are compared. The Iterative Method Huang et al. [Hua86] utilized an exhaustive iteration method to search for the converged answer. He decomposed the rotation matrix into three successive rotations around zaxis, xaxis, and yaxis, respectively. The rotation matrix can be described as: R=R, (1C) Rx ((<) Rz (0) (6.11) where is the rotation angle on yaxis 4 is the rotation angle on xaxis 0 is the rotation angle on zaxis The least square term can be rewritten as: N S2(,4)q Q)= R,() R, ()R,() 12 (6.12) i=1 This approach minimizes leastsquare fitting errors with respect to one variable at a time while restraining the other two axes with fixed angles. The method transforms a 3D problem to a 2D problem and the error will eventually converge to the minimum. The iterative method needs initial guesses of three rotation angles for the function value. Using the projection technique onto a 2D plane where the angles are reassumed, a new leastsquare value can be calculated with a new rotation angle. The diagrammatic concept is shown in Fig. 61. By revolving through three different orthogonal axes, the new leastsquare value can be calculated and minimized at each plane. When the function value approaches or is less than some preset threshold, the process is terminated. The three final rotation angle for each axis is the solution. According to the simulation work which was performed by Huang et al., it was discovered that the function value of the leastsquare term does have several saddle points. They claimed that even though the saddle points exist, the algorithm will find the minimum value, since the minimum value was calculated in 2D for each step. The Singular Value Decomposition (SVD) Algorithm The SVD algorithm for matrix operation is a powerful technique for dealing with sets of equations or matrices that are either singular or else numerically very close to Y (A1,B1,C1) (A2,B2,C2) I (0,0,0) nul 2 Z ' \ Figure 61: Using projection plane to calculate the iteration angle for each rotation singular. It is also the method of choice for solving most linear leastsquares problems where the number of equations is larger than the number of unknowns. This method starts with a set of simultaneous equations that can be reformatted into a simple matrix equation: Ax=b (6.13) where A is a matrix, b and x are vectors. Equation (6.13) defines A as a linear mapping 60 matrix from the vector space x to the vector space b. For example, consider the problem of fitting a polynomial p(x)=al+a2x+a x2+... +aNxN of degree Nl through M data points (xi,y), i=1,...,M. The coefficients a1,...,aN must satisfy the M linear equations to deliver a minimized leastsquared solution. Equation (6.14) shows a matrix multiplication form for solving the linear equations. Inverting a matrix and getting parameters for fitting a set of linear equations are 2 N1 a1 1x x, x ... x al 2 N1 XM XM... M a straightforward by using standard numerical methods such as continuous iteration, Gaussian elimination or QR decomposition [Ken89]. However, if these linear equations are close to linear dependence or have been classified as rank deficit, inverting a matrix becomes a very tedious and unstable procedure. The above algorithms sometimes fail to calculate the right matrix and the answers are easily diverged. However, the SVD numerical analysis method explicitly constructs orthonormal bases for the matrix and decomposes the matrix into proper mathematical format. Therefore, inverting the decomposed mathematical matrix form becomes much simpler because these terms are well defined and no divergence will occur even if the matrix function is close to singular. A simplified representation for the linear transformation A defined by the m x n matrix A by introducing two orthonormal bases can be easily computed. A decomposed matrix multiplication of A can be found at the cost of finding a great deal of information 61 about the eigenvalue system related to A. The arbitrary orthonormal set of vectors [u,,...,uJ which is used for the domain space n x n are generated while the other arbitrary orthonormal set of vectors [v1,...,vj is used in the domain space m x m, i.e., U== [Ul, ., u] (6.15) v= [v, ..., vmI (6.16) Then any m x n matrix A whose number of rows m is larger than or equal to the number of columns n (points are larger than unknowns), can be written as a product: A=UT.EV (6.17) where A is an (m x n) matrix U and V are (n x n) and (m x m) unitary orthonormal matrices E = E (wj) satisfies wj = 0 if i j The decomposition of matrix A can always be done, no matter how singular the matrix is. A proof that any matrix could have SVD matrices is beyond the scope of this paper, interested readers should refer to the reference on matrix multiplication [Gol89]. The SVD algorithm organizes the following representation for the inverse matrix of A. The inverse matrix A' is immediately solved as follows: AT = VaE1 UT (6.18) Thus x and b can be related as: x= vE1 (UTb) (6.19) 62 If b and x have some systematic errors involved or more equations than unknowns, the problem turns out to be that there is no linear mapping between b and x. The closest answer we encounter finds the leastsquares solution. The calculation of leastsquare minimization, which relates the two sets of data points with a best estimated rotation matrix, is based on the SVD algorithm. To apply the SVD algorithm in solving the rotation matrix in the subspace of these data points, the following method approaches the goal in generating a unique and minimum matrix for rotation estimation. The use of the SVD algorithm to estimate the leastsquared parameters is based upon the following discussion. Suppose that A EIR" is a data matrix obtained by performing a certain set of experiments. If the same set of experiments is performed again, then a different data matrix, BEIR', is obtained. The possibility that B can be rotated into A with the rotation matrix R is explored by solving the following problem: minimize II ARB IIF subject to RTR = Ip This is a typical application for solving the leastsquares problem. If the leastsquared answer is achieved, the rotation matrix provides the best match between two sets of data points. Note that if RE IR" is orthogonal, then expanding the above equation we have: IARBfII. N = (ARB) T(ARB) (.20) /1 N =J (ATA+BTB2ATRB) i=1 =Trace(ATA) +Trace(BTB) 2 Trace (ATRB) where Trace is defined as and a, is the diagonal element of the matrix A. N Trace(A) =x aii (6.21) i=1 Therefore, the above equation for solving leastsquares minimization is equivalent to the problem of maximizing the Trace(ATRB). Maximizing R can be found by calculating the SVD of ATB. Indeed, if UT(ATB)V = E = diag(a, ....,ap) is the SVD of the matrix and we can define the orthogonal matrix Z by Z = VTRU, then: Trace(RA B) = Trace(RUEVT) = Trace(ZE) P P (6.22) = E z1iai ai i=1 i=1 In order to prove that equation (6.22) is valid, we have to prove the following lemma first: Lemma: For any positive definite matrix AAT, and any orthonormal matrix B, the following equation exists. Trace(BAAT) Proof of Lemma: Let a, be the i' column of A. Then: Trace(BAAT) =Trace(ATBA) = af (Bai) (6.24) 1 From the Schwarz inequality, the following equation is valid, ar (Bai) Hence, Trace (BAA T) <.E aTai=Trace (AA T) (6.26) Clearly, the upper bound is attained by setting the R = VUT for Z = I,. Therefore, the procedures to maximize the trace of equation (6.23) can be explained as the following algorithm. This algorithm to solve the overestimated system can be summarized in several logical steps. The following explanation about the SVD algorithm can be applied to solve the rotation matrix for leastsquared matching in two 3D data sets. Given A and B in IR", the following algorithm finds an orthogonal REIR" such that II ARB II is minimum. (1) Calculate C = ATB (2) Compute the SVD UEVT = C, then UTCV = E. Save U and V. (3) Obtain R = VUT (4) Then we have: RATB=VUTUEVT=VEVT (6.27) We know that the trace of RATB has the largest value. Then the VUT is the rotation matrix we want and is the orthogonal polar factor of ATB. From the above algorithm, the matrix R which would minimize II ARB  p is calculated. This matrix delivers the closest solution for matching two 3D data sets A and B. The matrix R could be further transformed into three orthogonal Euler transformations which are as follows: where c represents cosine function and s represents sine function. Respectively, a, f3 and cacp sacy+caspsy sasycaspcy Q=sacp cacysaspsy casy+saspcy (6.28) sp cpsy cpcy 7 are three Euler angles for Yaw, Pitch and Roll angles in the orthogonal coordinate system. With proper matrix manipulation, the three rotation Euler angles can be calculated from the relationship contained inside the matrix function. If we assume the resulting matrix as another form such as the equation (6.29), the angles are easily calculated inside the matrix. rll r12 r13 Q=r21 r22 r23 (6.29) r31 r32 r33 where a = Atan2(r21,rll) 0 = Atan2(r31,sqrt(rll112+r212) y = Atan2(r32,r33) when 0 = 90 => a = 0, y = Atan2(r12,r22) 3 = 90 = > a = 0, y = Atan2(r12,r22) If the experimental data have some types of systematic errors (such as digitization errors), this rotation matrix will apparently give the best estimation of the rotation angles. Arun et al. [Aru87] developed a noniterative algorithm which involves the SVD of a 3x3 matrix to solve the leastsquare fitting of two 3D point sets. The purpose was to compare different algorithms in matching two 3D data points and the conclusion was that the SVD algorithm is the most stable and the fastest among iterative methods. 66 Pelizzari et al. [Pel91] has applied the SVD algorithm to successfully relocalize the patient's treatment position. A system of aligning two sets of 3D patient contour data was reported. The surface fitting method determines the transformation between two coordinate systems by finding the transformation which best matches models of the same surface as defined in each of the coordinate systems. Fragments of the source codes obtained from Pellizzari and Chen which were written in C on VAX/780" have been modified and adapted into the University of Florida PC system. This program provides the leastsquares solution of the 3D point matching problem. This source code has been extensively modified because of system differences and the specification parameters of six degrees fitting. This modified software package will be considered as the standard protocol for comparing the calculated results which are obtained from the methods of optimization techniques. The Vector Analysis of Robotic Method A potentially useful method for 3D point fitting is robotic vector analysis. Since sixdegree motion parameters are generally an important topic in robotic analysis, a geometrical approach is performed to solve the matching problem. Similar to robotic manipulators, the mechanisms for solving translation vectors and rotation angles are dominated by highly nonlinear trigonometric equations. The approach here is to accomplish the solutions by relating the translational and rotational mechanisms directly to the trigonometric formulation. This procedure provides deeper understanding of geometrical representation of implicit correspondence between two data sets using homogeneous matrices. 67 The vector analysis starts with initiation of a certain type of mechanism (please refer to the following Fig. 62). From the designed mechanism, we can analyze the geometrical relationship of different joints. The mechanical setups of each arm and joint are confined in advance to limit the movements. Thus the translation and rotation parameters can be calculated with the kinematic equations by both forward and inverse analyses. The manipulator study is to detail the algorithm which models the six degree offreedom serial manipulator into a single closed loop spatial mechanism of mobility one. The algorithm is explained and organized as follows: (1) First determine the location of three fiducial marks with respect to a fixed coordinate system. Determine relationship between the "head" coordinate system and the fixed coordinate system, (FTH). (2) Perform the forward analysis to determine the relationship between the fixed coordinate system and the head ring system. Let us define the fixed coordinate system as (FT6) and invert the coordinate as (6TF). (3) Determine the relationship between the fixed coordinate system and the head coordinate system by (6TH) = (6TF)(H). (4) Determine the required position and orientation of the head ring system so that the reference points are in the necessary position. (5) Then perform the reverse analysis to determine the joint values to move reference points to the required position. If properly defined for all the joint values, such as limitation of each slide length and each rotation angle, it is possible only one solution will fall within allowable joint limits. Table Figure 62: Kinematic model of the apparatus in solving the translational and rotational parameters However, the above approach of robotic vector analysis requires exact mapping of two data sets. If the noise vector in equation (6.1) becomes larger than a certain level, this algorithm will not solve the correct coordinates of the 3D data points. Therefore, error analysis is important while performing coordinate transformation. In other words, if this algorithm is applied, the optimization procedure should be incorporated into the process to obtain the leastsquared errors of the two data sets. For demonstration purposes, the program runs on a Silicon Graphics" workstation with the geometric mechanism displayed. The Optimization Methods Utilizing optimization techniques to solve the highly nonlinear leastsquared problem is a promising approach. The purpose of this approach is to apply optimization ~L 8/ 2 69 algorithms and to find the optimum parameters for automatic rotation and translation for patient relocalization. Numerical optimization techniques have been widely used in industry to improve production efficiency [Van84]. The optimization method provides improved accuracy, less computation time, and logical design procedures. Since we are dealing with a multidimensional optimization problem (i.e., three rotation and three translation parameters), the simplicity of onedimensional optimization has been lost. The optimization criteria is to minimize the object function, which is the leastsquared difference between two 3D patient contour and anatomical point data sets. The numerical optimization techniques are heavily dependent on personal taste. A method which does not require any derivatives of the design function was chosen, which could save a lot of computation time. This direct approach depends upon the converging criterion of the leastsquare function for two 3D data points. Three different optimization techniques are performed to calculate and compare the rotation and translation parameters for repositioning the patient's treatment position. Therefore, the irregular shape of a tumor for fractionated treatment is closely correlated during different treatment courses. Hooke and Jeeves Pattern Search Method The Hooke and Jeeves pattern search technique is a climbing method that does not require the use of derivatives. The algorithm has ridge following properties and is based on the premise that any set of design moves that have been successful during early experiments will be worth trying again. The method is based on the assumption of a simple unimodality of the object function value and is used to find the minimum of a multivariable, unconstrained function of the form: Object function=F(x1,x2, ... .x) (6.30) The simple logic diagram for this method can be seen in Fig. 63 and the detailed mathematical chart is shown as Fig. 64. The algorithm proceeds as follows. First, a base point in the feasible design space is chosen together with exploration step sizes. For example, three Euler angles are assumed zero as the initial guessing values. Next, an exploration is performed at a given increment along each of the independentvariable directions following the logic shown in Fig. 65 and Fig. 66. Figure 63: Logic diagram of Hooke and Jeeves optimization algorithm Whenever a functional improvement (improvement of the object function value) Hookes and Jeeves Pattern Search Mlntmization Algorithm 71 is obtained, a new temporary base point is established. Once this exploration is complete, a new base point is established and a new "pattern move" takes place. This pattern move consists of an exploration along a line between the new base point and the previous base point. Mathematically, this exploration is defined as: (k+l) = (k+ +a (Bk+1 B (6.31) B!,0 Bi i a B)( where Bi,o0'1 becomes a new temporary base point or "head." In this expression i is the variable index, k is the linear stage index, and a is an acceleration factor that is greater than or equal to 1.0. Once the new temporary base point has been found, an exploration about this point is instituted to see if a better base point can be found. If the answer is "yes", the pattern proceeds with this improved base. When the step size has been decreased below a predetermined value and still no substantial change in the merit value of object function can be achieved, the procedure terminates. Because of its simplicity to program, this method is one of the most popular techniques for computer implementation. The merit function is set to calculate the minimum errors between two data sets. In calculating the rotation parameters, the leastsquared minimum value of the 3D coordinates between two translated data points has been chosen as the object function. The object function is expressed as the following format: Object Function = F(xi,yi,zi) n (6.32) i= iJ(xix2,)2(y1 Y2,i)2 (1,i 2, 2 i=1 where n is the total number of contour points. 72 If the Hooke and Jeeves optimization method can search for the minimum value of the above object function, the spatial planes which are formed by two sets of 3D data points should be exactly aligned. That means the two sets of 3D points are matched by proper translational and rotational operations. The calculated parameters (ax, Ay, Az, a, 0, y) can perform the six degreeoffreedom coordinate transformation in the working space. The Downhill Simplex Method in Three Dimensions The Downhill Simplex method is due to Nelder and Mead [Nel65]. A simplex is an ndimensional, closed geometric figure in space that has linear line edges which intersects at (n+1) vertices. For example, in two dimensions this figure would be a triangle. In three dimensions it is a tetrahedron. Search schemes are based on the simplex utilizing the function values of each vertex. The basic move in this method is to find a reflection point of the Simplex and generate a new vertex point and a new simplex. The choice of reflection direction and the new vertex points relies on the location of the worst point in the simplex. The new reflected point is called the complement of the worst point [Pre92]. If the algorithm oscillates back and forth instead of going toward to the extremum, the second worst point should be used to locate the complement point. The Downhill Simplex method only requires function evaluation rather than derivative evaluation. The simplex points are manipulated by reflection, contraction, and expansion. This simplex method adapts itself to the local landscape, elongating down long inclined planes, changing direction on encountering a valley at an angle, and contracting in the HOOKE AND JEEVES ZEROORDER OPTIMIZATION METHOD Figure 64: Hooke and Jeeves optimization process (X is input parameters, N is dimension of searching space, P is factor of step size decrement, and A is the increment on each dimension which can be justified) PATTERN SEARCH Figure 65: Pattern search in the Hooke and Jeeves optimization process EXPLORATORY MOVES (MOVE IN '+' DIRECTION) Figure 66: Exploratory move in the pattern search approach of the Hooke and Jeeves optimization process 76 neighborhood of a minimum [Nel65]. By repeating the reflection, contraction or expansion procedures, this action contracts the simplex towards to lowest point, and will eventually bring all points to the minimum value. The simple diagram which describes the optimization method is shown in Fig. 67. Figure 67: The diagram of Downhill Simplex optimization algorithm There are two difficult situations in this type of optimization technique. First, the answers might be the local minimum points instead of the global minimum points. This difficulty has been found in using the Simplex method on a fourdimensional surface having a long, curved, valley with extremely steep sides; along the valley bottom the function varies considerably compared with the required accuracy for the minimum DownhillSmplex Geometrical Adaptation Minimization Algorithm 77 function value [Nel65]. In other words, the calculated results may be trapped by the saddle points. Secondly, computer execution time might not be economical if too many function variables are involved. To solve the rotation matrix and translation vector for patient relocalization in fractionated stereotactic radiotherapy, the patient position will be closely repositioned by the laser before fine adjustment. It means that the motion parameters we seek are fairly close to the global minimum points. Chance of being trapped at a local minimum point will be decreased to an acceptable level. Besides, we only have six design variables for repositioning the patient's head, therefore the PC is a suitable tool for the optimization calculation. The detailed mathematical description of Downhill Simplex method in a flow chart format can be seen in Fig. 68. The Simulated Annealing Optimization Method in Three Dimensions The object function in the determination of the rotation parameters (three Euler angles) consists of nonlinear variables and may have multiple minima. The condition of nonlinear variables excludes certain minimization methods such as linear or quadratic programming, which need a specific object function format. The condition of multiple minima is the limitation of some quick optimization algorithms because they might yield a solution at the nearest local minimum value. A parameter estimation method that is resistant to becoming trapped by local extrema is named Simulated Annealing. The usefulness of the Simulated Annealing method comes from the possibility of accepting a "uphill" move in the parameter space. That is, instead of moving towards to the minimum value, the searching process sometimes jumps up to the higher function value, which is against the search trend. This method prevents the function value from being DOWNHILL SIMPLEX METHOD MOVE ALL POINTS 1/2 THE DISTANCE TOWARD THE BEST POINT I Figure 68: The Downhill Simplex optimization algorithm for the calculation of 6D fitting parameters 79 trapped in the local minimum and provides the necessary search to the global minimum. The optimization process for this geometrical minimization problem can be observed in Figure 69: The diagram of Simulated Annealing optimization algorithm the Fig. 69. The basic idea of the Simulated Annealing method is also applicable in solving the threedimensional leastsquared minimization problems for the calculation of rotational parameters. With continuous three dimensional spaces, this algorithm will ideally find the global minimum in the presence of many local minima. The four elements required by the Metropolis procedure in Simulated Annealing algorithm are as follows [Met53, Boh86]: (1) The value of the object function, in our case, which is the leastsquared value of the Simulated Annealing Geometrical Annealing Minimization Algorithm difference from two 3D data sets. (2) The system state at the point x, which is the vector at that specified state. (3) The control parameter T, which is similar to the temperature, with an annealing schedule by which it is gradually reduced. T is also an important factor to specify the "cooling schedule" or "cooling curve". (4) A generator of random changes, taking a random step from x to x+ax. Simulated Annealing consists of three basic steps sequentially to constitute one iteration of the procedure [Met53, Kir83]: (1) A random generator function, which starts with the current configuration or initial value of object function, generates a new function value. The magnitude of all the changes to the variables of each iteration is called step size, which can be characterized inside the program according to the user's preferences. (2) An acceptance function which accepts or rejects the new configuration based on the changes AV of the object function. If the object function value decreases, the new configuration is always accepted; otherwise it is randomly accepted with a probability analogous to the Boltzmann distribution exp(AV/T), where T, the temperature, controls the degree of "uphill climbing" of the object function. (3) An updated function which determines the step size and temperature is used in the next iteration procedure. There is an easy way to perform the Simulated Annealing procedure. It is possible to combine simulated annealing with a version of Downhill Simplex searching that allows upper and lower bounds on all parameters [Pre92]. The Simulated Annealing 81 procedure can come close to the true global minimum in parameter space, while the simplex search, now begun near the global minimum, can search for a true global minimum without fear of falling into the local minimum value. The implementation procedure of the Metropolis scheme (cooling curve) is by adding a positive, logarithmically distributed random number (it will always give out the positive term), proportional to the annealing temperature T, to the stored function value associated with every vertex of the simplex geometry. Then subtract a similar random variable from the function value of every new point that is tried to replace the old vertex. Like the standard Metropolis process, this method always accepts a true downhill step, but sometimes accepts the uphill one. The selection of downhill or uphill steps is to provide the random process and the local minimum trap can be avoided. At a finite value of T, the simplex expands to a scale that approximates the size of the region that can be reached at this temperature, and then executes a stochastic, tumbling Brownian motion within that region, sampling new, approximately random, points as it does so. If the temperature is reduced sufficiently slowly, it becomes highly likely that the simplex will shrink into the region containing the lowest relative minimum encountered. The Simulated Annealing algorithm is a fairly new method for industrial or even medical uses. However, this method provides several extremely attractive features, which supercede other optimization algorithms. First, it is not easy fooled by the "quenching" procedure which is used for almost all the other optimization methods. Provided that sufficient general reconfigurations are 82 given, it wanders around freely among those local minima of depth less than T. Therefore, as T decreases, the number of qualifying for frequent visits of the minimum values are gradually reduced. Second, changes that cause the greatest energy differences are shifted over when the control parameters are large. The potential advantage of this method is that the step size adjusts itself automatically, expanding or contracting to a scale that approximates the size of the configuration space that can be reached at a given value T. This decision becomes more permanent as T gets lower. Smaller refinements are then achieved to the real solution. Simulated Annealing of the Downhill Simplex algorithm has been applied to determine the rotation angles in the 3D data fitting procedures. In essence, Simulated Annealing is used to provide the starting vertices for a Downhill Simplex search, thus minimizing the risk of simplex searching being trapped in a local minimum. The converging criteria is that optimization stops when the difference in the object function values of the highest and lowest configurations in the simplex is less than 0.001 [Pre92]. The Simulated Annealing algorithm provides a comparative basis for calculating the rotation parameters of the three Euler angles. Conclusion SVD algorithm is a standard technique to solve the leastsquared overestimated system in the numerical analysis. Solving rotation matrix is a good representation of the mathematical method for 3D data correlation. Different optimization methods have been utilized to design the engineering products or analyze various system performance. 83 However, they have not been examined to solve the 3D data correlation problems. The relationship between two patient treatment position can be evaluated and crosscalibrated through different optimization methods. Along with the SVD algorithm, systematic studies and performance evaluation are pursued. The following chapters will compare the effectiveness and convenience to correlate patient treatment position. In a clinical oriented environment, these optimization algorithms provide a very useful tool to perform patient repositioning and refixation for high precision fractionated radiotherapy treatments. CHAPTER 7 ANALYSIS OF RELOCALIZATION ERRORS Comparison of Algorithms There are four different algorithms being compared in this section, namely SVD algorithm, Hooke and Jeeves algorithm, Downhill Simplex algorithm and Simulated annealing algorithm. These methods provide a solid basis for calculation of the six degreesoffreedom fitting parameters. In order to observe the results of these four algorithms and compare their potential clinical utilization, two studies were performed. One study used the same clinical data sets from Pelizzari [Pel89, Pe192] and compared the error terms of those algorithms. The leastsquared object function values for those algorithms during optimization are also considered and plotted to demonstrate the manipulation process. The other study randomized the input data set to generate sets of stochastic outputs and performed inverse calculations for finding the bestfit parameters. A total of one hundred random outputs from the same input data were produced to perform the statistical analysis. Reverse mapping of the input and output data sets reveals the stability and feasibility of the optimization algorithms. The analysis of relocalization errors are summarized as follows : (1) Analysis of clinical digitized data set of patient: 85 Two sets of clinical patient digitization data provided by Pelizzari [Pel92] were utilized in this study. The two data sets consist of independent 3D points from a PolhemusT digitizer for different time periods. The first data set carried 18 points of data, and the second carried 51 points by randomly selecting anatomical points of the patient surfaces. The calculation results of the error estimation are shown in Table 71 and Table 72. From the results, the optimization algorithms potentially reveal smaller error terms. Since each digitized point includes different sources of errors, then the mean error, mean square error and the root mean square error are the most appropriate forms to analyze the merit function of the 6D fitting process. The optimization procedures are basically the searching steps for the minimum function values. The leastsquared function values of the two data sets can be approximated and minimized through the minimization algorithms. Continuous evaluation of the function values shows the 6D fitting parameters are improved during the calculation process and the closest answers can be obtained once the preset tolerance value is reached. Fig. 71 and Fig. 72 show that the Hooke and Jeeves algorithm is always searching for the better function value through each iteration. The pattern search algorithm looks for the sharpest decrease of function value and finds the best object minimization value within specific tolerance, which is the difference between two consecutive object function values during the searching process (defined as 0.0001 in the calculation algorithms). The function value is finally stabilized and reaches the minimum. This method is effective and fast if the object function is not very 86 complicated. The 6D fitting parameters are easily extracted from the program and show the best match for the two 3D data sets. Observation of the graph shows that a sharp ridge occurs in the first few iterations and then the object function stabilizes for the further evaluation. This is the typical ridge characteristic of the Hooke and Jeeves algorithm. This quick converging method can bring the object function value to its minimum within a very short time and introduces no singular problem (the invert matrix does not exist) which appears often in matrix manipulation. In order to observe object function changes in the Downhill Simplex method, the same methodology used in the Hooke and Jeeves algorithm is applied to analyze the trend in each iteration step. Fig. 73 and Fig. 74 show the function value changes during the iteration process. Fluctuation of the function values during each iteration step is observed simply because of the reflection, extrapolation and contraction procedures which occur on the simplex vertices. Geometrically, the function values will finally search for the minimum state and calculate for the best 6D fitting parameters. A simulated annealing process is also incorporated into the optimization algorithm to verify the calculated results, which are the best fitting parameters for the two data sets. The inverse calculation might deliver ambiguous answers. For example, instead of moving a patient to the bestfit values, the calculation might give a false local minimum that does not necessarily relocate the patient's treatment position and orientation. With the possible global searching characteristics of the annealing algorithm, the final calculation can be compared as a bench mark for the optimization algorithm. For each iteration loop a total number of 25 random iteration steps are performed with 87 three random number generating processes. Inside the first loop of iteration, a function value is selected (not necessarily the smallest) as a starting base. Then the random process continues and compares with the previous object function value. If there were no improvement, maintain the original function value. Otherwise, the function value is replaced with the new randomized value. This process will continue until a preset threshold or tolerance is reached and the whole procedure is terminated. The random operation provides the perturbation which brings out the global minimum object function value. Fig. 75 to Fig. 78 demonstrate the annealing process which each iteration loop consists of twenty five iterations and the function value settles down after a long updown random approach. (2) Random output data set analysis: Another approach is that one set of the digitized data is considered as the seed. To compute and identify the merits of each algorithm, a random rotation and translation program was developed and incorporated to simulate the random translation and rotation characteristics. From the utilization of a random number generator [Pre92], the seed points are stochastically translated and rotated to all the possible combinations. For simulation purposes, the following parameters for the transformation are assumed to manipulate the input data points. They are (A) 50 mm < Random Translation (three orthogonal axes) < 50 mm (B) 20 degrees < Random Rotation (three Euler angles) < 20 degrees (C) Gaussian Noise: Mean = 0; Variance = 1 The selection criteria is based on the reasonable assumption of the patient location Object Function Evaluation (1 8 (Hooke and Jeeves algorithm) points) 20 40 60 80 100 Iteration Number Figure 7.1: Function value changes algorithm) through each iteration (Hooke and Jeeves Object Function Evaluation (5 1 (Hooke and Jeeves algorithm) points) 0 20 40 60 80 100 Iteration Number Figure 72: Function value changes algorithm) through each iteration (Hooke and Jeeves Object Function Evaluation (1 8 (Downhill Simplex algorithm) 10 9 8 7 6 5 3 2 0 0 10 20 30 40 50 60 Iteration Number points) Figure 73: Function values changes through each iteration (Downhill Simplex algorithm) Object Function Evaluation (5 1 (Downhill Simplex algorithm) 0 10 20 30 40 50 Iteration Number points) 60 70 80 Figure 74: Function value changes through each iteration (Downhill Simplex algorithm) Object Function Evaluation (18 points) (Simplex Annealing algorithm) 7 6 04 .2 4  L 2 1 0 0 20 40 60 80 100 120 140 Iteration Loop Figure 75: Function value changes through each iteration loop (Simplex Annealing algorithm) 93 Object Function Distribution (18 pts) (Simplex Annealing algorithm) 7 o 4 (1 i~3 0 5 10 15 20 25 Iteration Number Between Each Loop Figure 76: The function value changes within each iteration loop of Simplex Annealing algorithm 