High precision treatment for fractionated stereotactic radiotherapy

MISSING IMAGE

Material Information

Title:
High precision treatment for fractionated stereotactic radiotherapy
Physical Description:
vii, 278 leaves : ill. , photos ; 29 cm.
Language:
English
Creator:
Yang, Ching-Chong, 1962-
Publication Date:

Subjects

Genre:
bibliography   ( marcgt )
theses   ( marcgt )
non-fiction   ( marcgt )

Notes

Thesis:
Thesis (Ph. D.)--University of Florida, 1994.
Bibliography:
Includes bibliographical references (leaves 269-277).
Statement of Responsibility:
by Ching-Chong Yang.
General Note:
Typescript.
General Note:
Vita.

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 001981551
notis - AKF8482
oclc - 31920336
System ID:
AA00003232:00001

Full Text








HIGH PRECISION TREATMENT FOR
FRACTIONATED STEREOTACTIC RADIOTHERAPY















By

CHING-CHONG 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 3-D CONTOUR POINTS ....... 30

Statement of Problems .................. ... .......... 30
Methodology .................... .................. 32

5 3-D 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 3-D 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
3-D Camera System Data Tests .......................... 124
Correlation of CT Image and 3-D Camera Data Sets ............. 131

8 THE HARDWARE SYSTEMS AND THE VERIFICATION PROCEDURES
.. ........ ......... ... ........... ....... .. ...... ...136

System Characteristics ................................ 136
High Precision Camera System ..................... 136
The Real-Time 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

Ching-Chong Yang

April, 1994


Chairperson: Frank J. Bova
Major Department: Nuclear Engineering Sciences

Fractionated Stereotactic Radiotherapy (FSR) refers to high precision,

multifraction, external-beam irradiation that utilizes accurate refixation techniques. Three

optimization algorithms were developed to calculate six-degrees-of-freedom 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 real-time 3-D 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 halo-shaped of metal immobilization device, the patient's anatomy relative

to the Brown-Roberts-Wells (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 non-coplnar radiation beam

combinations on the computer. Virtual simulation differs from most of the radiation

oncology routine practices by replacing conventional treatment simulation with 3-D 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 dose-response 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 TG-21 [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 error-prone 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,

external-beam 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 X-ray

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 LINAC-based 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 3-D patient information and Beam's Eye-View [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 better-oxygenated 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 dose-time 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 linear-quadratic survival relationship (a/f3

ratio) of the cell responses to radiation is used to assess the fractionated doses (Fig. 2-1).

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 2-1: Linear-Quadratic model for multifraction dose assessment [Per87]


AVM, for instance, is a late-response tissue with small a/O3 ratio. Fractionation

will not generate benefits for destruction of this late-responding 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 worst-case






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 in-field 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 cephalo-caudad 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 vacuum-form 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 non-invasive 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 so-called

"Hat" and "Head" method. From 3-D 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 3-D digitizer with low-frequency

magnetic field to detect the 3-D 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 table-mounted array of microphones. With

the triangulation of the emitter signals, the 3-D 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 bearing-holder 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 3-1: 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 3-2, 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 3-2: 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 state-of-the-art treatment planning software was developed to perform rapid

dose calculation and to display isodose distributions. A standard nine arc plan, shown

in Table 3-1, 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 x-ray 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 3-1: 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 LINAC-based 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

degree-of-freedom 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 3-D information of the patient's anatomy but also the coordinates of the









fiducial markers relative to the 3-D 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 non-coplanar 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 head-clamping

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 3-D 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 3-D translation and three-axis 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 3-D 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 3-D 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 tumor-dose 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 3-D

orientation of the tumor should be matched. Therefore, with six degree-of-freedom

fitting procedures (translation and rotation) instead of using three degree-of-freedom

parameters (translation only), the tumor shape can be matched as closely as possible.

Correlation of two sets of data points is a many-to-one 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 least-squared approach is the favored choice for solving this kind of problem. Since

all the 3-D data points contain certain types of digitizing and systematic errors, finding

the best fit of these two different 3-D 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 3-D 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 degree-of-freedom 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 3-D 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 3-D 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 one-to-one 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 edge-thresholding 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 pseudo-solid 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 degree-of-freedom 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

3-D 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 degree-of-freedom 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 3-D reference set of points, which represent the 3-D structure

with respect to the patient's own coordinates (this coordinate system will be the same as

that of the LINAC treatment coordinates). Fig. 4-1 shows the conceptual drawing of

registration of the 3-D 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 2-D





















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 4-1: 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 2-D point sets of (xk, y) for 0
spline consists of n-1 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,...,n-2) 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 (x-xk) 3+a2(x-xk) 2+a (x-xk) +a (4-1)

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 "not-a-knot"

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 2-D concept of a spline

can be extended into 3-D 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. 4-2.

A set of computer programs were written to perform data sorting, resampling and

reweighting of the initial 3-D data set. After the 3-D reference data set is saved as well

as the treatment planning information, a separate 3-D data set which follows the same

contour lines and anatomical points is acquired from the 3-D digitizer for mapping. Six

degree-of-freedom fitting parameters are then calculated and compared with different

algorithms in order to find out the best correlation of these two 3-D data sets.

A real time 3-D digitizer is incorporated into the data acquiring system for the

fractionated treatment schedules. This system is well calibrated with respect the LINAC

































Figure 4-2: 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. 4-3 is presented to integrate various

steps for the 3-D 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 4-2: Procedures for the calculation of 6-D 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
3-D 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. 5-1. The photodetector of the detection system will sense an

alternating intensity. One cycle of intensity change represents a mirror displacement of

one-half 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 3-D 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 5-1: 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

sub-system.

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 3-D coordinate

measurements, the concept of stereo vision is applied. Mapping a 3-D scene onto an

image plane is a many-to-one 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 3-D coordinate information.

Stereo vision involves obtaining two separate imaging views of an object of

interest. In Fig. 5-2, 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 3-D 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 5-2: 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. 5-3. 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 5-3: Film dosimetry system setup (camera, light box, and stand). A calibration
plate is designed and analyzed for the system resolution.








Table 5-1: 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 5-1 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 5-1 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 3-D 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 3-D 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 3-D 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. X-ray 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 X-ray 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 set-ups also play an important role in 3-D 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 3-D 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 X-Y coordinates on the image plane of the

detector are proportional to the X-Y 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 3-D coordinate

location of the LED. Those data sets can be transferred into the workstation or PC for

3-D 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 X-Y

direction, and 0.15 mm in depth which is perpendicular to the X-Y plane.

Pixsys" also utilizes a array of CCD cameras which are placed above the patient

table. These cameras track and report the 3-D coordinates of tiny LEDs which are built

into the instruments such as stylus, forceps, biopsy needles or ultrasound transducers.

The 3-D 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 3-D 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 3-D 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 3-D 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 3-D

coordinate measurements. This digitizer can interact with CT image data sets and

perform a real-time 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 6-D

(translation plus rotation) tracking localization system ADL-1". 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, 3-D 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 3-D information. Because the system can track objects in real

time, if the sensor is attached to a stylus, the 3-D 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) source-to-sensor

separation. As the separation increases, the positional accuracy degrades linearly.

Although this system has reasonable accuracy for measuring 3-D 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 3-D coordinates. Therefore, this system is

excluded for our application.

Ultrasonic System

The ultrasonic system for determination of 3-D 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 range-finding function as well

as calculate the spatial coordinates of the selected points. Science Accessories Corp.

(Southport, CT) developed the Model GP-8-3DT ultrasonic range-finding device as a 3-D

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 3-D 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) Pre-calibrated 3-D 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 3-D 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 3-D REFERENCE DATA SETS


Basic Principles


While performing motion analysis of a rigid body, the estimation of the 3-D

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 3-D

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 3-D data points to the patient's contour. Since the patient's head can

be considered as a rigid object, the 3-D 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 3-D 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 3-D data set with a corresponding 3-D 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 3-D 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 3-D data. Therefore, assume that we are

given two 3-D 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 least-square best fit of the following equation:


2= IIp- (Rp+T) 112 (6.2)
i=1

Therefore by minimizing the above least-square 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 least-square 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 least-square 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 least-squared fit answer. It is only the estimation of the centroid difference

between two 3-D 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 3-D 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 Pi-P (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 least-square function of two 3-D 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 6-D (rotation and translation) to 3-D (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 least-square 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 z-axis, x-axis, and y-axis, respectively. The rotation matrix can be described as:

R=R, (1C) Rx ((<) Rz (0) (6.11)


where is the rotation angle on y-axis








4 is the rotation angle on x-axis

0 is the rotation angle on z-axis

The least square term can be rewritten as:

N
S2(,4)q Q)= -R,() R, ()R,() 12 (6.12)
i=1

This approach minimizes least-square fitting errors with respect to one variable

at a time while restraining the other two axes with fixed angles. The method transforms

a 3-D problem to a 2-D 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 2-D plane where the angles are reassumed,

a new least-square value can be calculated with a new rotation angle. The diagrammatic

concept is shown in Fig. 6-1. By revolving through three different orthogonal axes, the

new least-square 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 least-square 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 2-D 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 6-1: Using projection plane to calculate the iteration angle for each rotation

singular. It is also the method of choice for solving most linear least-squares 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:

A-x=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 N-l through M data

points (xi,y), i=1,...,M. The coefficients a1,...,aN must satisfy the M linear equations

to deliver a minimized least-squared 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 N-1 a1
1x x, x ... x- al



2 N-1
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.E-V (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:


A-T = VaE-1 UT (6.18)

Thus x and b can be related as:


x= v-E-1- (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 least-squares solution. The calculation of

least-square 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 least-squared 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 A-RB IIF subject to RTR = Ip

This is a typical application for solving the least-squares problem. If the least-squared

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:


IA-RBfII.
N
= (A-RB) T(A-RB) (.20)
/--1
N
=J (ATA+BTB-2ATRB)
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 least-squares 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 over-estimated system can be summarized in several

logical steps. The following explanation about the SVD algorithm can be applied to

solve the rotation matrix for least-squared matching in two 3-D data sets.

Given A and B in IR", the following algorithm finds an orthogonal REIR"

such that II A-RB 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 A-RB | p is

calculated. This matrix delivers the closest solution for matching two 3-D 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 sasy-caspcy
Q=-sacp cacy-saspsy 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 non-iterative algorithm which involves the SVD

of a 3x3 matrix to solve the least-square fitting of two 3-D point sets. The purpose was

to compare different algorithms in matching two 3-D 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 3-D 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 least-squares solution of the 3-D 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 3-D point fitting is robotic vector analysis. Since

six-degree 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 non-linear 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. 6-2). 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-

of-freedom 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 6-2: 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 3-D 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 least-squared 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 non-linear least-squared

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 one-dimensional optimization

has been lost. The optimization criteria is to minimize the object function, which is the

least-squared difference between two 3-D 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 least-square function for two 3-D 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. 6-3 and the detailed

mathematical chart is shown as Fig. 6-4. 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 independent-variable

directions following the logic shown in Fig. 6-5 and Fig. 6-6.


Figure 6-3: 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 least-squared minimum value of the 3-D

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(xi-x2,)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 3-D data

points should be exactly aligned. That means the two sets of 3-D points are matched by

proper translational and rotational operations. The calculated parameters (ax, Ay, Az,

a, 0, y) can perform the six degree-of-freedom 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 n-dimensional, 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


ZERO-ORDER OPTIMIZATION METHOD


Figure 6-4: 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 6-5: Pattern search in the Hooke and Jeeves optimization process









EXPLORATORY MOVES


(MOVE IN '+' DIRECTION)


Figure 6-6: 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. 6-7.


Figure 6-7: 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 four-dimensional 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


Downhill-Smplex 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. 6-8.

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 6-8: The Downhill Simplex optimization algorithm for the calculation of 6-D
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 6-9: The diagram of Simulated Annealing optimization algorithm


the Fig. 6-9. The basic idea of the Simulated Annealing method is also applicable in

solving the three-dimensional least-squared 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 least-squared value of the


Simulated Annealing Geometrical Annealing Minimization Algorithm








difference from two 3-D 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 3-D 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 least-squared over-estimated

system in the numerical analysis. Solving rotation matrix is a good representation of the

mathematical method for 3-D 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 3-D data correlation problems. The

relationship between two patient treatment position can be evaluated and cross-calibrated

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

degrees-of-freedom 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 least-squared 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 best-fit 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 3-D 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 7-1

and Table 7-2. 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 6-D fitting process.

The optimization procedures are basically the searching steps for the minimum

function values. The least-squared function values of the two data sets can be

approximated and minimized through the minimization algorithms. Continuous

evaluation of the function values shows the 6-D fitting parameters are improved during

the calculation process and the closest answers can be obtained once the pre-set tolerance

value is reached.

Fig. 7-1 and Fig. 7-2 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 6-D fitting parameters are easily extracted from the program and show

the best match for the two 3-D 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. 7-3 and Fig. 7-4 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 6-D 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 best-fit 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 pre-set

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. 7-5 to Fig. 7-8 demonstrate the annealing process which each iteration loop

consists of twenty five iterations and the function value settles down after a long up-down

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 7-2: 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 7-3: 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 7-4: 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 7-5: 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 7-6: The function value changes within each iteration loop of Simplex
Annealing algorithm