UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Analytic three dimensional image reconstruction from projections Kinahan, Paul Eugene 1988

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-UBC_1988_A7 K56.pdf [ 5.48MB ]
Metadata
JSON: 831-1.0085037.json
JSON-LD: 831-1.0085037-ld.json
RDF/XML (Pretty): 831-1.0085037-rdf.xml
RDF/JSON: 831-1.0085037-rdf.json
Turtle: 831-1.0085037-turtle.txt
N-Triples: 831-1.0085037-rdf-ntriples.txt
Original Record: 831-1.0085037-source.json
Full Text
831-1.0085037-fulltext.txt
Citation
831-1.0085037.ris

Full Text

ANALYTIC THREE DIMENSIONAL IMAGE RECONSTRUCTION FROM PROJECTIONS B y P a u l Eugene K i n a h a n B . A . S c , The Universi ty of B r i t i s h Co lumbia , 1985 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF T H E REQUIREMENTS FOR T H E D E G R E E OF MASTER OF APPLIED SCIENCE in T H E FACULTY OF GRADUATE STUDIES DEPARTMENT OF PHYSICS We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA A p r i l 1988 © Pau l Eugene Kinahan , 1988 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department The University of British Columbia Vancouver, Canada DE-6 (2/88) Abstract Thi s work presents an analytic three dimensional image reconstruction a lgor i thm that was developed for a proposed volume-imaging P E T scanner. The development of the a lgor i thm was motivated by the scanner's abil i ty to collect an order of magnitude more data than current P E T systems and the lack of an efficient a lgori thm that could use the extra data. The a lgor i thm is based on an extension of the Recovery Operator of Orlov[68] and operates by convolution i n object space. This method of operation sets it apart from other analytic direct image reconstruction algorithms that rely on Fourier transforms. The algori thm is tested wi th ideal data and parameters that are appropriate to the new P E T scanner. The results of the test show that the a lgor i thm behaves as expected except for a 17% overshoot i n the reconstructed value i n one area. A n explanat ion of this artifact is suggested, al though not verified. Final ly , the efficacy of the algori thm is demonstrated by proving that it is func-t ional ly equivalent to Fourier transform methods. Contents Abstract ii List of Figures v i Acknowledgements vii i 1 Basics of P E T 1 1.1 Introduction 1 1.2 Phys ica l Processes of PET 2 1.3 Introduction to Image Reconstruction on P E T scanners 7 1.4 Limi ta t ions of P E T 9 2 Three Dimensional P E T 13 2.1 The Case for Volume Imaging Systems 13 2.2 Previous Work i n Three Dimensional Image Reconstruct ion 19 3 Reconstruction A l g o r i t h m 23 3.1 A n a l y t i c F o r m 23 3.1.1 Mathemat ica l Descript ion of the Prob lem 23 3.1.2 T w o Dimensional Image Reconstruction 24 3.1.3 Three Dimensional Image Reconstruction 29 3.1.4 Theory of Orlov 36 in 3.1.5 A New Life for Orlov 's Recovery Operator 37 3.1.6 Ana ly t i c F i l t e r Form 40 3.1.7 L imi t s of Integration 41 3.2 Discrete F o r m 43 3.2.1 Discrete Reconstruction 43 3.2.2 Three Dimensional Backprojection 46 3.2.3 Calcula t ion of Dig i t a l F i l ter 49 4 Implementation 64 4.1 Reconstruction Details 64 4.2 Results 66 5 Equivalence to Fourier Methods 72 5.1 Fourier Transform of the Reconstruction Fi l te r 72 5.2 Compar ison to Frequency Space Fil ters 75 Summary 79 Bibl iography 82 A P S F Derivat ion 96 A . l Introduction 96 A . 2 Detector Geometry 97 A . 3 Derivat ion of P S F 98 B Derivation of L(T X XT) 100 C Derivatives of L(6f,8T) 106 C . l F i r s t Derivative 106 C.2 T h e Second Derivative 107 iv D Fi l ter P r o g r a m List ing 109 E Reconstruction P r o g r a m Listing 129 v List of Figures 1.1 Basic physical processes of P E T (not to scale) 4 1.2 Por t ion of r ing of detectors of typical P E T system 6 2.1 Cross section showing shielding for a planar imaging detector 15 2.2 I l lustrat ion of paral lax error 16 2.3 Detector geometry of TRIUMF PET scanner 18 3.1 T w o dimensional object and one of its one dimensional projections . . . 25 3.2 F i l t e r ing and backprojection of projection data 27 3.3 Coordinates of the direction vector r 30 3.4 T y p i c a l volume imaging P E T scanner geometry 32 3.5 Var i a t ion of detector P S F wi th point source posi t ion 34 3.6 M i n i m u m acceptance angle for an invariant detector P S F 35 3.7 Truncated spherical detector and projection plane geometry 38 3.8 Bi l inear interpolation 48 3.9 Possible intersections of the line = 9 w i th a p ixe l 53 3.10 Regions of integration in a pixel for I2ij-(9r) 54 3.11 Sketch of the second derivative of 1/L{9[\ 9T) 56 3.12 Integration scheme over a pixel for i"3 t )j ;(0 T) 58 3.13 D i g i t a l filter for 9T = 90° and $ = 10° 61 3.14 Negative of the filter w i th the center 9 x 9 pixels set to zero 62 v i 4.1 The x-y plane of p(x) 67 4.2 T h e z-x plane of p(x) 68 4.3 Values of p(x) along the x-axis 70 4.4 Values of p(x) along the z-axis 71 5.1 T w o dimensional and three dimensional coordinates i n frequency space . 77 B . l Geometry of the L(? x xT) function 101 B.2 Intersection of the (a, 6, c) plane and the region G 103 v i i Acknowledgements I acknowledge the following people whose contributions have made this thesis better than what i t would otherwise be. D r . Joel G . Rogers, my research supervisor, who provided me w i t h everthing that a graduate student could need. It is his knowledge and patient guidance that have made this work possible. D r . R i cha rd R . Johnson, m y faculty supervisor, who has supplied me w i t h valuable encouragement and advice on many matters. D r . R o n a l d Harrop , who has followed my work carefully, and made important suggestions and corrections. Drs . T o m R u t h and Wayne M a r t i n , who cheerfully answered my questions about the rest of P E T . M y Friends of Mickey Mouse; Chr is Backhouse, Connie Baxter , Eraser Duncan , Shell i and W i l l i e Fajber, Marce l Lefrancois, and Dave MacQui s t an , who have helped me wi th many hours of often painless diversions. M y co-workers i n the U B C / T R I U M F P E T group; E d Grochowski , K e n Pederson, G o r d Riley, and Ne i l W i l k i n s o n . These guys helped make work intresting. The Science C o u n c i l of B r i t i s h Columbia , which has provided me wi th financial support for two years. M y Mothe r and Father, whose love, support, and encouragement make everything possible. v i i i Chapter 1 Basics of P E T 1.1 Introduction Posi t ron Emiss ion Tomography (PET)[100] is an imaging methodology pr imar i ly used i n nuclear medicine. It is a system that can determine the posi t ion and amount of a radioactive tracer present at each location wi th in a solid object. A P E T system can determine the locat ion and amount of tracer mater ial an order of magnitude more accurately than other nuclear medicine imaging systems. T h i s is done by taking advantage of the properties of the part icular radioactive tracers used. Tomography (from the Greek: tomo - slice) is a branch of the field of image reconstruction. Image reconstruction is common to many disciplines including radio and solar astronomy, electron microscropy, industr ial nondestructive testing, radar and sonar systems, and several types of medical scanners. The common basis these fields have is the use of tomography to deduce information about the internal nature of an object solely from external views, called projections, of the object. In nuclear medicine, image reconstruction as used by P E T allows the study of i n vivo biochemical reaction processes i n objects as complicated as the human brain. P E T systems are currently only used as research tools to s tudy (mainly human) bio-1 chemistry, but there is a large amount of investigation into the possible roles of a P E T scanner as a diagnostic a id for diseases such as Alzheimer 's [62], Huntington 's [41], and Parkinson's [12]. T h e potential of P E T stems from its abil i ty to determine internal function as op-posed to internal structure. Structure can be measured w i t h an X - r a y Computer ized A x i a l Tomography (CAT) scanner or a Nuclear Magnet ic Resonance (NMR) scanner, but neither of these methodologies can determine what internal processes are occuring. A l though the next three sections give some of the necessary background for the rest of this thesis, the best general introduction to P E T is the Scientific American article by Ter-Pogossian et a l . [100]. M o r e detailed introductions are the works by Sorenson and Phelps [93] and Phelps and Mazz io t t a [74]. 1.2 Physical Processes of PET The basic process of P E T is the decay of a radioactive isotope by the emmision of a positron, and the subsequent mutual annihilat ion of the posi tron and an electron. The positron, also called a / 3 + particle, can be considered to be a posit ively charged electron, a j3~ particle. W h e n the positron and the electron are annihi lated, they are converted into energy i n the form of two 511 ki lo electron-Volt (keV) gamma rays that travel away at 180° to each other. This reaction is symbolized by, P+ + P~ —> 2 7 The positron comes from the decay of specific radioisotopes such as 18F, 1 5 0 , 13N, and 1 1 C , wi th half-lives ranging from minutes to hours. These radioactive isotopes are used to label compounds by replacing either naturally occuring isotopes or other atoms. The purpose of the P E T system is to determine the amount and locat ion of /3+ + /3~ annihilations, and this information is then interpreted as being the same for the labelled 2 compound. For example, the compound 18FDG is a labelled analogue of glucose that undergoes a metabolism i n the bra in similar to that of glucose. Since glucose is a metabolic fuel, by determining the distr ibution of 18FDG w i th in a bra in , researchers can investigate metabolic act ivi ty i n different regions on the b ra in [93,92,91]. Th i s method, called tracer kinetic studies, has been used w i t h other nuclear medicine systems, but to a much lesser degree of accuracy. The key process that sets P E T apart is that the two gamma rays from the annihilat ion travel away i n opposite directions along a ' l ine of interaction' . These photons can be detected i n coincidence by de-vices around the annihi lat ion point i n a PET scanner detector as shown i n figure 1.1. Note that i n figure 1.1 the size of the nuclei and the positron, and the distances they travel, are greatly exagerated relative to the size of the detector circle. B y using its two detected endpoints, the line of interaction can be measured more accurately and efficiently than wi th other nuclear medicine scanners. Since the gamma ray photons are travelling at the speed of light, they w i l l arrive at the detector wi th in ~ 1 0 - 9 sec. of each other. This coincidence detection allows the separation of different annihi lat ion events in time. B y detecting enough events it is possible to measure the location and quantity of positron emitters. A l though the fundamental process used by PET is straightforward, the P E T systems themselves are not. Current ly there are approximately 60 P E T scanners w i th no two being exactly alike, w i th many quite different i n method of operation and appearance. It is possible, though, to give a description of a typical system, and some common variations. A standard system would consist of one or more rings of crystals capable of detect-ing the annihi lat ion photons. The rings are stacked to form a cy l indr ica l volume to enable detection of gamma rays from a three dimensional object inside the rings. The crystals, such as sodium iodide (Nal ) or bismuth germinate ( B G O ) , work by absorb-3 4 ing the single 511 k e V photon and then releasing the energy as a shower of thousands of secondary) photons [64,35,73]. E a c h detector r ing can consist of as many as 1200 crystals [11] or as few as 6 [65]. A s shown i n figure 1.2 the secondary photons pass through some form of optical coupling, such as lucite, to the surface of a series of photo mult ipl ier tubes ( P M T s ) . The entire assembly is shielded from outside visible or U V light because of the sensitivity of the P M T s . Usually, the surfaces of the crystal not opt ical ly coupled to anything are covered w i th a reflective mater ial i n order to redirect more secondary photons to the P M T s . The P M T s can detect ind iv idua l secondary photons and convert the energy of the photons to a series of electrical pulses that identify which crystal absorbed the original gamma ray. Some P E T systems have been bui l t around two rotating parallel plates that sweep out a cyl indr ica l volume [5,2]. The plates are M u l t i W i r e Propor t ional Chambers ( M W P C s ) i n which the gamma ray causes an avalanche breakdown i n a gas, and the avalanche charges are picked up by a series of closely spaced wires. The electrical signals from the detector subsystem are analyzed to determine co-incidence i n t ime of two gamma rays, and the resulting information is stored i n a computer as a single coincidence event. In addi t ion some machines attempt to mea-sure the difference i n arr ival times between the the two gamma rays i n order to localize the annihi la t ion point on the line of interaction [64,35,92]. Such Time-Of-Fl igh t ( T O F ) systems have had mixed results because of the added complexity bo th of measuring and incorporat ing the T O F information. W i t h the stored information of coincidence events ( typical ly 100,000 to 10,000,000 events) it is possible to recover the internal dis t r ibut ion of positron emitters, and thus the radioactive tracer dis t r ibut ion. Th i s w i l l be shown i n chapter 3. The above describes the role that a P E T system fulfills w i th in the overall framework of P E T studies. The two other major areas of a P E T study are the preparation of 5 > Figure 1.2: Por t ion of r ing of detectors of typical P E T system the labelled compounds, which requires a cyclotron [106,39], and the mathematical model l ing of the compound's reaction sequence i n order to analyze the data from the P E T scanner [47,61]. A l l three of these components are necessary to perform a nuclear medicine P E T study [74]. Th i s thesis is concerned wi th the last step i n the P E T scan phase of the P E T study. Th i s is the step where a l l the coincidence events that are stored i n the computer are transformed into the reconstructed internal dis t r ibut ion of positron emitters - this is the process of image reconstruction. 1.3 Introduction to Image Reconstruction on PET scanners Thi s section gives a brief introduction to image reconstruction as it applies to PET. The theory of three dimensional image reconstruction w i l l be covered i n excruciating detail i n chapters 3 to 5. Image reconstruction is a branch of applied mathematics that deals w i th recovering the internal details of an object from projections through the object. Projections are the integrations of internal distributions of some physical quanti ty on a straight line through an object. W h e n the reconstruction is performed by a computer, the general process is sometimes called Computer ized Tomography ( C T ) . In P E T scanners, projections are formed from the stored lists of coincidence events, and a computer performs an image reconstruction to recover the internal tracer dis t r ibut ion. Th i s is akin to determining the structure of a three dimensional chinese puzzle inside a locked cupboard by taking a series of X - r a y photographs of the cupboard at different angles. F r o m this point on, the lists of coincidence events w i l l be referred to as projection data or, simply, projections. The problem of analytically determining a function from its projections was first 7 solved i n 1917 by Johan Radon [78]. Unfortunately his work was overlooked un-t i l after the same problem had been re-solved by others almost 60 years later. B y the mid-1960's several electron microscopists had realized that they could determine the internal structure of smal l objects by using Fourier transforms of the projections [30,79]. Since this was before Fast Fourier Transforms ( F F T s ) were used, the trans-forms were time-consuming to compute and were only useful for objects w i th high degrees of rotat ional symmetry, such as bacteriophages. In 1967 Bracewell , a radio astronomer, showed how to perform an exact image reconstruction from projections without resorting to Fourier transforms [8]. Th i s work was soon duplicated by others 1 such as Ber ry and Gibbs [6]. Since then the study of image reconstruction has been driven by several pract ical applications, resulting in its establishment as a separate field of study. The or iginal tomographs developed for medicine employed X-rays to image a two dimensional plane, or slice, perpendicular to the long axis of the patient [9]. The image reconstructed from the collected projection data gave a picture of the internal dis t r ibut ion of tissue and bone wi th in the slice. A series of slices next to each other were then used to bu i ld a picture of the three dimensional tissue and bone distr ibut ion inside the patient. For P E T systems the same method of successive slices is used to b u i l d a three dimensional picture of tracer distr ibution. The major difference between P E T and CAT scanners is that CAT scanners transmit X-rays through an object and look at the amount absorbed, whereas P E T scanners look at the emission of gamma rays from sources inside the object. P E T system designers are constantly t ry ing to improve resolution and sensitivity [64,35,99]. A broad definition of resolution is the abil i ty to distinguish, or resolve, two 1 T h i s re-inventing of the wheel was a common occurence in the early days of C T , possibly due in part to the interdisciplinary nature of the field. 8 or more objects from a single larger object; so a system wi th a better resolution can dist inguish separate points placed more closely together. Sensitivity can be thought of as the systems abil i ty to receive data; thus a system wi th a higher sensitivity can detect more events i n a fixed time period. The first few generations of P E T scanners used two dimensional image reconstruction for each slice, but as the increasing ex-perience of scanner designers has led to improved resolution and sensitivity, it has been realized that it is useful to aquire data that cannot be used by two dimensional reconstruction theory. A P E T system that shows the need for direct three dimensional image reconstruction is being developed at T R I U M F 2 . Th i s w i l l be discussed i n more detail i n the next chapter. 1.4 Limitations of P E T Since i t is important not to lose sight of the overall picture, several of the most impor-tant l imi t i ng factors w i l l be listed i n this section. Th i s w i l l be done i n the order that information passes through a typical P E T system as described i n section 1.1, start-ing w i t h the tracer decay by positron emission and ending wi th the display of the reconstructed image. A s the labelled compound decays, the positron is ejected from the isotope nucleus w i th a non-zero velocity [93]. Th i s means that the positron w i l l combine wi th an electron at a varying distance from the actual tracer. Th i s positron 'wander' has a gaussian dis t r ibut ion about its emission point from the tracer that has a F u l l W i d t h Ha l f M a x i m u m ( F W H M ) of 1 m m to 5 m m depending on the labelled compound. A n added effect of the positron's velocity is that the two annihi lat ion gamma rays vary i n the angle that they travel apart. The separation angle has a gaussian variat ion about 2 T h e TRI-University Meson Facility located in Vancouver B . C . Canada. 9 180° w i t h a F W H M of ~ 0.5°. B o t h of these effects contribute to a fundamental l imi t of how accurately the line of interaction can be measured [64,35,33]. After the emission and before the detection of the two 511 k e V gamma rays, the object itself can degrade the accuracy of measurement by its interaction w i th the gamma rays. Th i s influence is main ly through self absorption and Compton scatter of the gamma rays by the object. Self absorption occurs when one or both of the gamma rays is absorbed by the object. A manifestation of this is seen i n the centres of strong self-absorbers where the measured act ivi ty is lower than the true activity. If one or bo th of the gamma rays is scattered wi th in the object but s t i l l detected, then an incorrect line of interaction is recorded. B o t h scatter and self-absorption can be corrected for w i th varying degrees of success [93,64,48,38,16,102]. A t the stage where the gamma rays are detected by the P E T scanner, the l imi t ing factors become entirely machine dependent, w i th the most important factors being de-tection efficiency, intrinsic resolution, and detector dead time [35]. Detector efficiency refers to the probabi l i ty that a gamma ray passing through a detector w i l l interact w i th the detector and thus be observed. The probabil i ty of interaction can vary from about 10% for M W P C based detectors to 95% for systems wi th large B G O crystals [64,35,99,34]. The intrinsic detector resolution determines how accurately the point of interaction between the gamma ray and the detector can be known, and varies from about 1 m m to 15 m m . Usual ly one scinti l lat ion crystal is coupled to one P M T so that by looking at the outputs of the P M T s it is possible to find out which crystal absorbed the gamma ray, so i n this method the intrinsic resolution is the crystal size. Th i s ac-curacy can also be a function of several other factors, such as the amount of light from the scintal lat ion and how much of it makes it to the P M T , and the efficiency of the P M T i n converting the secondary U V photons into electrical signals [43]. These signals are processed by the detector electronics, and the processing of each event causes the 10 detector to be shut down for a period known as the dead time, during which no other events can be detected [60,59]. The two dominant factors i n system dead t ime are the t ime taken to reject s ingle 3 events, and the t ime taken to process a coincidence event to the point where the information can be transfered to the computer. Addi t ional ly , the detector electronics cannot always differentiate between true coincidence events and random coincidences, where a random coincidence is a detection of two single events wi th in the coincidence time window [64,14,73]. Th i s is an error, since the two single events are the results of annihilations that are close i n t ime but not necessarily i n space. Mechanical structure affects performance i n many, often subtle, ways. The more obvious ways are the physical posit ioning of detector units, geometrical efficiency, and shielding. Th i s shielding, which is usually lead or i ron to efficiently absorb gamma rays, is used to reduce the detection of unwanted events [64,93,73,58]. These events are either gamma rays originat ing outside the scanners F i e l d O f V i e w ( F O V ) , internally scattered gamma rays, or events that supply information that cannot be used i n the image reconstruction process. Th i s last case, and geometrical efficiency, w i l l be discussed further i n the next chapter. In addit ion, shielding can scatter gamma rays s imilar to the prompt scatter caused by the object. Once the detected events are stored as projections, the factors affecting perfor-mance are quantization and statistics [93,18]. In the mathematics of analytic image reconstruction the projections are approximated as true integrals of act ivi ty along lines through the object. The data is actually a set of sums of ind iv idua l events and the accuracy of the approximation depends on the number of events collected. Since the fundamental process is nuclear decay, which has Poisson statistics, the error of approx-3Single events are those where a second gamma ray was not detected within the allowed coincidence time window. The ratio of single to coincidence events ranges from 50:1 to 100:1. 11 imat ion is proport ional to 1/y/n, where n is the number of events. The error can be reduced by scanning for longer periods, but this may not be possible because of patient considerations or isotope decay. Increasing sensitivity always improves statistics, but may also increase the amount of scattered gamma rays that are detected, as w i l l be shown i n section 2.1. Quant iza t ion is the subdivision of a continous variable into a series of discrete subunits, and always degrades accuracy [93,46,57]. It is necessary i n a P E T system because the image reconstruction is performed by a computer that can only work wi th discrete quantities. Quant izat ion shows itself i n several places; first i n the discretiza-t ion of the projection data. If the smallest measurable region contains an even smaller strong point source, then the act ivi ty of the point source w i l l be averaged over the entire region, thus giving a false lower value. Th i s is known as the par t ia l volume effect. The second instance of quantization error is i n the manner that the data is manipulated dur ing the image reconstruction process. Final ly , the reconstructed i m -age can be presented i n different degrees of quantization. A l though quantizat ion can usually be increased to the point where it no longer l imits performance, this is often not done because of the consequent increases i n needed computer memory and compu-tat ion time. The effects of quantization on the accuracy of image reconstruction w i l l be discussed i n section 4.2, but analyzing the effects of the other l imi t ing factors on the reconstructed image is beyond the scope of this work. 12 Chapter 2 Three Dimensional P E T Some of the concepts introduced i n the previous chapter w i l l be used i n this chapter to show the advantages of scanning an entire volume at once instead of scanning an object as a series of separate planes. Some of the problems of volume imaging and their solutions w i l l also be discussed. Since three dimensional imaging leads to three dimensional image reconstruction theory, the chapter ends w i t h a survey of current three dimensional reconstruction algorithms. 2.1 The Case for Volume Imaging Systems Standard two dimensional image reconstruction algorithms only use projection data that a l l lie i n a single plane [9,84,3,52,87,105,45,44]. Pa r t ly for this reason the cur-rent generation of P E T scanners form three dimensional images by reconstructing a series of adjacent slices and then stacking the planes together to form a volume image [10,20,19]. The separate planes, or slices, are called transaxial slices since they are perpendicular to the long axis of the patient. Cross-plane events are those i n which the two gamma rays are detected i n different detector planes. In other words, cross-plane data cannot be assigned to a single detector plane and thus cannot be used by standard two dimensional image reconstruction algorithms. 13 To keep detectors from picking up cross-plane events, annular shielding rings called inter-slice septa are placed between adjacent detector rings as shown i n figure 2.1 [73]. Note that the detector surface has cyl indr ical symmetry about the centre line. Septa are also used wi th in a slice to absorb gamma rays that could introduce paral lax error [64]. Th i s error arises when scinti l lat ion crystals are deep i n comparison to the desired intr insic resolution as shown i n figure 2.2. Since the interaction point wi th in the crystal is not known, there can be a large error i n the calculated line of interaction for a gamma ray that arrives at an oblique angle to the crystal . Th i s problem can be solved by determining the depth of interaction of the gamma ray i n the crystal or by making the crystal shallow i n comparison to the intrinsic resolution. The latter method also greatly reduces the probabil i ty of interaction between the crystal and the gamma ray and thus reduces overall system sensitivity. There is often an exception made to the rule of not using cross-plane data. Th i s is the practice of using data that is detected on two adjacent planes to reconstruct a plane ly ing midway between the two detector planes [73]. Th i s type of image is called an inter-slice plane as opposed to a true plane 1 . Inter-slice planes do increase the number of planes, and thus the axial intrinsic resolution, but the sensitivity, resolution, and F O V are usually dramatical ly different from those of true planes. Th i s can lead to substantial inaccuracies i n a three dimensional image formed by stacking an alternating series of true and inter-slice planes. In a mul t i - r ing scanner the only other method of improving the axia l resolution is to make the crystals smaller i n the axial direction. Referring to figure 2.1, using smaller crytals is equivalent to br inging the inter-slice septa closer together. Th i s de-creases the total volume i n the scanner's F O V unless more detector rings are added, consequently increasing complexity and cost. Addi t ional ly , experience w i t h P E T scan-1 Formed from events detected by a single detector ring. 14 m / y J s~ s J s\ object dxis o end absorbed 1 - 7 7 7 7 y Y ' J J j / s J J. V./ / / / / / Yff/SYf/ V////////////, 0U~ 4 Figure 2.1: Cross section showing shielding for a planar imaging detector. 15 i O ^ i o A' / if /_ \ / / / Figure 2.2: I l lustration of parallax error 16 ners has led to the general rule that improving intrinsic resolution by using smaller discrete c r y s t a l / P M T detector units leads to reduced system sensitivity. A n i l lustrat ion of this trade-off can be seen by comparing two prototype single r ing P E T scanners. The Donner 600 scanner [36] has a best case resolution of 2.6mm and a corresponding sensitivity of 7100 events/sec per / x C i / m l of act ivi ty i n a standard cyl indr ica l flood phantom. Compare this to the P C R - I scanner [96] that has a best case resolution of 4 .5mm and a sensitivity of 45,000 events/sec per ^ C i / m l of activity. T h e way to escape this bottleneck and increase both resolution and sensitivity is to remove the inter-slice septa. Th i s increases the solid angle of coverage inside the detector and leads to an increase i n the scanner's geometric efficiency [64,73,2]. The geometric efficiency measures how well a scanner can receive data from a point source moved everywhere inside the detectors F O V . B y removing septa, system designers can keep resolutions of ~ 5 m m and increase sensitivity 5 to 20 times [24], but there are three problems that immediately arise 2: parallax error, increased scatter, and how to incorporate the addi t ional data into the image reconstruction algori thm. The first problem has been solved by a new P E T detector unit developed at T R I U M F [83,80,82]. Th i s unit uses a large (15 x 15 x 2.5 cm) N a l crystal that is optical ly coupled to a posit ion sensitive P M T . B y analyzing the secondary photon pattern the three dimensional interaction point of a gamma ray inside the crystal can be located. Th i s obviates the need for using separate rings and septa, and a complete scanner is being developed by arranging 9 of these units i n a cyl indr ica l pattern as shown i n figure 2.3 3 . A n addi t ional important benefit is that intrinsic resolution is now the same i n a l l directions, since the interaction point is now measured as a continuous variable as opposed to a discrete variable. 2 Which is why the septa were there in the first place 3 Figure 2.3 is courtesy of R. Moore, TRIUMF. 17 > 18 Thi s abi l i ty to measure data that has not been directionally selected is an important one. The human bra in and other organs are complicated three dimensional structures and are not designed to to be conveniently looked at i n a series of planes. A n y data collection system should be able to produce reconstructed images that have equal sensitivity and resolution i n a l l directions. It is known that the amount of scattered events that are detected increases faster than the rate of detection of true coincidences as the geometrical efficiency increases [73]. Th i s has been verified by recent Monte Car lo simulations [82,94,101] of the effects of increasing geometric efficiency. Since scatter is largely featureless for extended sources, gains i n statitstics are more important than reducing scatter. So the net effect of increasing the scanner's total solid angle is s t i l l an improvement i n the signal to noise ratio. The last problem of being able to use a l l the addi t ional cross-plane data is the subject of this thesis. 2.2 Previous Work in Three Dimensional Image Reconstruction Direct three dimensional image reconstruction differs from indirect reconstruction i n that the three dimensional volume image is formed al l at once, instead of as a stacked series of separate two dimensional image reconstructions. A n attempt to use cross-plane data i n an indirect three dimensional image recon-struct ion has been made by Daube-Witherspoon and Muehllenher [28]. The basis of the method is to assign the cross-plane event to the plane midway between the two planes i n which the gamma rays were detected, or to average the event information between a l l the intervening image planes. After this a standard image reconstruction 19 is performed on the image data. Tests of this method have shown that even though statist ical quality is improved, artifacts can be present i n the final three dimensional image [28]. Several of the most common methods of direct three dimensional image reconstruc-t ion fall under the heading of iterative reconstruction techniques. Algor i thms such as the Image Space Reconstruction A l g o r i t h m [27,103,29], also developed by Daube-Wi therspoon and Muehllenher, and other series expansion methods [55,15,69,49,4] essentially a l l rely on many iterations of sucessively more accurate reconstructions. Th i s is different from analytic reconstruction methods that perform a complete image reconstruction i n one operation where no guess work is required. In general, these iterative algorithms work by t ry ing to estimate an object that w i l l give rise to the ob-served projections. The differences between the observed projections and projections formed from the estimated object are calculated and are used to modify the estimated object, at which point the process is repeated. The iterations of modifying the esti-mated object are stopped when the differences between the sets of projections are less than some pre-determined amount. These algorithms differ main ly i n their method of using the differences between the projections as a guide for modifying the estimated object. The advantages of iterative reconstruction methods are that they can use cross-plane data and a pr ior i knowledge of the object, and i n some cases they have low noise generation characteristics. T w o relatively minor disadvantages are that edge artifacts and random noise form i n the reconstructed image at h igh i teration counts [99], and that it is often difficult to know when to stop the i teration process. The major disadvantage is the computation time required for an image reconstruction. There are two dimensional iterative image reconstruction algorithms that do not require inordinate amounts of computat ion time [99,98,53,88,89,50,56], but three dimensional 20 implementations of these algorithms require many hours of computat ion for a single image [27]. Th i s amount of t ime for an image reconstruction is not acceptable, as a P E T system can be expected to scan up to four times a day [85]. The slowness of the iterative algorithms combined wi th the advantages of us-ing cross-plane data has resulted i n efforts to develop an analyt ic three dimen-sional image reconstruction algorithm. Some of these efforts have been directed to-wards the cone-beam data collection method tr ied on C A T scanners [32,66,72,90], but these have been based on algorithms developed for the parallel-beamdata collection method characteristic of P E T . These methods rely on the use of Fourier transforms [97,70,21,104,95,26,77,71,22,86,76]. These algorithms w i l l be out l ined i n the next chap-ter, but the point to be made here is that i n general any analytic reconstruction method is much faster than an iterative one. Fourier transform techniques have the drawback that since a computer implemen-ta t ion requires the use of Discrete Fourier Transforms ( D F T s ) , they also suffer from the side effects of D F T s . These effects are poor noise characteristics and the possi-bi l i ty of aliasing artifacts. Addi t ional ly , while Fourier methods can use some of the cross-plane data, they cannot use a l l of the available data. The amount of useable data depends on the relative geometries of the object and the P E T scanner system. T w o new methods of three dimensional image reconstruction that are promising have been devised w i t h specific detector geometries i n m i n d . B o t h are intended for use w i th two rotat ing parallel M W P C detectors and both can use most of the detected data. The method of Defrise et a l . [31] is essentially a multi-pass Fourier transform technique similar to that of C h o et al . [21] that requires more t ime than a single analyt ic image reconstruction, but should be considerably faster than a true iterative method. The only algori thm to date that uses a l l projection data i n a single pass is that of Clack et a l . [23]. This last method is based on a series of two dimensional 21 image reconstructions at different angles, and is the most promising of the alternative methods surveyed. The final method we w i l l discuss, and the one on which this work is based, is the generalized analytic Recovery Operator of Or lov [67,68]. The recovery operator, which w i l l be explained i n the next chapter, is an analytic reconstruction technique that has the same requirement of an invariant detector P S F , but does not rely on D F T s . Tha t is, it works entirely i n object, or image, space as opposed to frequency space. It w i l l be shown i n chapter 3 how the recovery operator can be adapted to a P E T scanner system. 22 Chapter 3 Reconstruction Algorithm 3.1 Analytic Form Before developing a solution to the problem of direct three dimensional image recon-struct ion, it is necessary to describe the problem mathematically. T h e notat ion used is based on that of Orlov [67,68] and Denton et al . [32]. 3.1.1 Mathematical Description of the Problem A n object of interest has an unknown internal dis t r ibut ion of posi t ron emit t ing ra-dioactive tracer. Th i s distr ibution is completely described by the function p(x), where x is a three dimensional vector. A s mentioned above, the projection data, p(xT; f ) 1 , are approximated as the integrals of tracer act ivi ty along a straight line through the ob-ject. Th i s approximation w i l l be treated as exact from this point on. The relationship between the original dis tr ibut ion and the projections is described by the equation, /oo p(x + tr)dt (3.1) -oo Semicolons are used throughout to delimit variables from parameters. For example in p(xT;r), xT is the variable and r is treated as a parameter. 23 ) where xT — X — (X-T)T is a two dimensional projection of the three dimensional vector x onto the projection plane S T and t i s a dummy integration variable along the direction f . The projection plane is the plane passing through the origin and perpendicular to the unit direction vector f . F rom this, the projection data is a two dimensional data set, which lies on the projection plane, that has xT as a variable and f as a parameter. In other words, a different value of f picks a new two dimensional data set. The last few sentences describe a three dimensional object w i th two dimensional projections, but it is easier to il lustrate w i th a two dimensional object and a corresponding one dimensional set of projections, as is done i n figure 3.1. In figure 3.1 it is assumed that p(x) = 0 outside the object of interest. In order to extend the geometry of figure 3.1 to the three dimensional image reconstruction problem, it is only necessary to let x be a three dimensional vector and let xT and f be two dimensional vectors, w i t h corresponding changes i n the dimensionality of p(x) and p(xT;f). In this case the projection data, p(xT; f ) , is considered to be a four dimensional data set wi th each two dimensional projection plane located by the two dimensional direction vector f . The problem to be solved can now be described as follows; given the measured four dimensional projection data set, p(xr;f), invert equation 3.1 to recover the original three dimensional tracer distr ibution, p(x). Th i s problem arises out of several other fields such as solar astronomy, seismology, and sonar. Conversely, any solution can be used by the same fields. Before leaping i n to tackle the problem, it is instructive to consider the simpler two dimensional image reconstruction problem. 3.1.2 Two Dimensional Image Reconstruction In two dimensional image reconstruction the direction vector, r , is completely de-scribed by the polar angle 6 as shown i n figure 3.1. In this case the solution for the 24 Figure 3.1: T w o dimensional object and one of its one dimensional projections 25 original source dis t r ibut ion is [84], dd \ dt'p(t-t';6)f(t') (3.2) where x = t cos 9, y = t sin 9, and f(t) = 6lD(t)-l/t2 where <$i.o(i) is the one dimensional Dirac-del ta function. The variable t above is the one dimensional equivalent to xT. The operation inside the square brackets is called convolution and is represented by an asterisk, *, as i n , Equa t ion 3.2 can be interpeted as a one dimensional convolution filtering operation followed by a backprojection that is repeated for a l l values of 9 between —7r/2 and n/2. This is i l lustrated i n figure 3.2 for one value of 9 and a part icular value of t = t0. The filtering operation wi th the filter function f(t) removes the distort ion caused by the backprojection operation. If the projections are directly backprojected without filtering, the two dimensional result, called S(iT) is where ** represents two dimensional convolution. The two dimensional image formed by direct backprojection is the original object convolved w i th the function In the discussion so far i n this section, the mathematical analysis assumes that the plane being imaged is infinitely th in . In reality this cannot be or else no data would be detected, so i n practice the data is collected from a plane of finite wid th . Th i s means /oo dx'f(x - x')g(x') •oo —or 26 II on Figure 3.2: F i l te r ing and backprojection of projection data 27 that the resulting reconstructed two dimensional image is assumed to have a constant value along a line perpendicular to the object plane. A t this point the scanner can be described using linear systems theory, that is, given a known input the output can always be determined. The input is related to the output by the transfer function or Point Spread Funct ion ( P S F ) . Th i s definition comes from setting the input to the system, p(x), to be the Dirac-del ta function, 6(x), since for any function / ( x ) , f(x) * *<5(x) = /(#), so S ( x ) = 6(x) * = --jr \x\ \x\ In other words this is the resulting image if a point source is placed i n the detector. The definition of P S F w i l l also be used i n the next section. It is possible to reverse the order of filtering and backprojection i n equation 3.2 by start ing w i t h the backprojected image i n the above equation and using a two dimensional convolution filter as i n the BackProject and F i l t e r ( B P F ) a lgor i thm [42]. ^(x) = S ( x ) * * | < 5 2 ^ ) - | 4 p } (3.3) where ( ^ ( x ) is the two dimensional Dirac-del ta function. The value of p(x) i n equa-tions 3.2 and 3.3 are identical, since the filtering and backprojection operations are linear and so their order is interchangeable. The Fi l te r and BackProject ( F B P ) method of equation 3.2 is used more often because it requires less computer memory. It also allows certain other operations such as those used i n tracer kinetic modeling to be performed directly on the projection data [13]. F r o m Fourier transform theory it is known that the convolution of two functions is equivalent to the mul t ip l ica t ion ( in frequency space) of their Fourier transforms [7]. For example, i f F{) is the Fourier transform operator, and F(f)=F{f(x)}, G{f) = T{g{x)} 28 then F{f{x)*g{x)}=F{f)G{f) so the convolution operation is equivalent to f(x)*g(x) = F-1{F(f)G(f)} where ^ r - 1 { } is the inverse Fourier transform operator. In this manner convolution filtering can be replaced wi th Fourier transform filtering, at least i n the ideal case of no quantiat ion error and perfect data. In pract ical implementations, convolution methods are used because of advantages of speed and low noise [63]. 3.1.3 Three Dimensional Image Reconstruction In three dimensional image reconstruction the unit direction vector, f, ranges over the two spherical coordinate angles 6T and <f>T as shown i n figure 3.3. If 0 < 6T < ir and 0 < <f>T < 7T, then the three dimensional source dis t r ibut ion, p(x), can be recovered from the set of projections, given by equation 3.1, by an extension of equation 3.2 to [32], p(x) = fV d<j>T f d6T[p(xT; 0T, 4>T) * * / ( £ . ) ] (3.4) Jo Jo where the reconstruction filter is / ( x T ) = S2D(\xT\) - — A n d s imilar i ly to equation 3.3 the result of backprojecting a l l the projections, £ ( £ ) , is given by, s ( x ) = l" d<t>T r derp{xT-eT^T) Jo Jo then p(x) is also given by [32], p(x) = E ( x ) * * * ( < 5 3 D ( x ) - 1 J i I l (3.5) x 29 ( y Figure 3.3: Coordinates of the direction vector r 30 where * * * is the three dimensional convolution operator and 8ZD(X) is the three dimen-sional Dirac-del ta function. Analogous to the two dimensional image reconstruction case, the convolutions shown above can be replaced wi th Fourier transforms [26]. Since most P E T scanners are cyl indr ica l i n nature, the projection coordinate 0T cannot range over 0 to 7r, but must lie wi th in some smaller range. For a typical scanner this range is TT/2 — ip < 8T < TT/2 + tp where ip is a parameter of the scanner geometry as is i l lustrated i n figure 3.4. Note that the origin of the coordinate system is at the centre of the cyl indr ical region. A s xp —> 0, the detector becomes a two dimensional planar detector and i f xp = 7r/2 the detector becomes an infinitely long cylinder that can detect a l l possible gamma rays. F r o m this it is easy to see how the geometrical efficiency is t ied to the parameter xp. The crux of the three dimensional image reconstruction problem is that for 0 < xp < 7r/2, none of the preceeding equations i n this section can be used to directly recover p(x). T h i s is because the reconstruction filters discussed cannot work w i th a detector system that does not uniformly sample an object i n a l l directions, nor expl ic i t ly take into account the non-planar nature of the images. In other words, the filters have the same form wherever they are applied i n equations 3.2 to 3.5. F r o m equation 3.1 and figure 3.4 it can be seen that the projection data gathered for a typica l volume imaging scanner w i l l have different qualities depending on the direction vector f . Appropr ia te reconstruction filters based on the Fourier transform method of three dimensional image reconstruction have been derived for the case where 0 < xp < IT/2 [70,26,86]. However these algorithms s t i l l require a spatially invariant P S F for the detector. To i l lustrate spatial invariance, imagine a detector that forms a complete sphere. A point source moving around inside the sphere appears equally bright wherever it is located. Th i s is because the amount of data that the detector can receive from 31 Figure 3.4: T y p i c a l volume imaging P E T scanner geometry 32 the object does not depend on the posit ion of the point source. For this imprac t ica l 2 scanner it can easily be shown that the P S F of this detector is l / | x ] 2 [68]. The case of a spatially variant P S F is i l lustrated by imagining a point source inside the detector shown i n figure 3.4. In this case i t is shown i n appendix A that the P S F , called h(x) , depends on posit ion and i n spherical coordinates is given by, h(r,9,<j>) = h(r,e) = ^rect(^j^j where f 1 i f |* | < 1/2 rect(x) = < 0 otherwise and r , 6, and <f> are standard spherical coordinates. Th i s is slightly different from the result quoted by Colsher [26], since he also included an arbitrary scaling factor of F r o m the above result it can be seen that as the point source moves around, its apparent brightness w i l l seem to vary because of the changes i n detector solid angle (called the acceptance angle) subtended. Th i s is demonstrated i n figure 3.5 for two locations of a point source. Note that figure 3.5 shows a side view of the detector surface shown i n figure 3.4. It is possible to obtain an invariant P S F by l imi t ing the acceptance angle at a l l points i n the object to the m i n i m u m angle that is determined at the edge of the object. For a spherically symmetric object, w i th radius R, centred at the or igin as i n figure 3.6, the m i n i m u m acceptance solid angle is 27n/>m m, where V'min = V7 + a rcs in ( i l cos i / ' / i2 ( j ) (3-6) and R^ is the detector radius as shown i n figure 3.6. W i t h V'min a s the acceptance angle, any point source located wi th in the bounds of the object w i l l appear to have the same strength regardless of its locat ion. B y using 2 A spherical PET scanner would be too large to be practical 33 Z AX'S Figure 3.5: Var ia t ion of detector P S F wi th point source posit ion. 34 ( { Figure 3.6: Minimum acceptance angle for an invariant detector PSF. 35 V'min a s a fixed acceptance angle, data that comes from the centre of the detector wi th an angle 6T < TT/2 — V>niin or 6T > TT/2 + V'min cannot be used. A l s o from figure 3.6 and equation 3.6 it can be seen that as R —• R^, V'min ~* 0- Th i s means that for any object of finite size, requiring a spatially invariant P S F means rejecting detected data. In some cases up to 75% of the available data is lost [24,23]. To summarize the material presented i n this section, it is possible to perform an analyt ic three dimensional image reconstruction of p(x) using Fourier transform methods by discarding enough data to obtain a spatially invariant P S F for the detector. Convolut ion methods had only been developed for the case where ip = 0 or ip = TT/2. Thi s leads to the next section that starts the investigation into using convolution methods for direct analytic three dimensional image reconstructions. 3.1.4 Theory of Orlov In 1975, the Russ ian electron microscopist S. S. Orlov showed that p(x), of dimension n , can be recovered from a complete set of projections, p(xT; f ) of dimension n — 1, by the following relationship [67,68], J ^ J * w s w ( 3 ' 7 ) G S T where, for three dimensional image reconstruction, S T is the plane through the origin perpendicular to r ; G is the area on the unit sphere that decribes the set of a l l f vectors, or projection directions, and L(f x x T ) is a geometric term, "which depends only on the direction of its argument, is equal to the length of the arc of the great circle ly ing inside the region G and passing through the endpoints of the vectors f and x • | x | _ 1 . " [68]. Since the set of projections is defined by the region G , Orlov also showed that for G to contain a complete set of projections, it must contain the domain of at least one semi-great circle. 36 To il lustrate the definitions of the terms i n equation 3.10 it is convenient to model the detector geometry as a truncated sphere of unit radius as shown i n figure 3.7. The results derived for a spherical-type detector can easily be converted for use wi th a cy l indr ica l detector. Note that i n figure 3.7 two views of the projection plane and its coordinates are shown, the correct one w i t h only the ly axis visible, and the same plane rotated about the ly axis by 90° i n order to show the lx vector and the relationship between the x, f, and / vectors. The region G i n the figure is s imply the por t ion of a hemisphere ly ing wi th in the bounds of TT/2 — xp < 6T < -rc/2. The coordinate system used i n the projection plane is defined by the following, / = XT — lxlx "1" lyly ly = T X lx F r o m this, lx, ly, and f form a right hand coordinate system since x = xT + tr. O n the projection plane the rectangular and polar coordinates are given by (lx, ly) and (/, 9{). Also shown i n figure 3.7 is an arc the length of which is the value of the function L(f x xT). The explicit form of L(9 x xT) is derived i n appendix B . The validi ty of equation 3.7 was demonstrated for the case where tp = IT/2 for a complete spherical detector. In this case L(j x xT) = IT and equation 3.7 reduces to equation 3.4. If xp < 7r/2, from figure 3.7 it is obvious that L(r x xT) w i l l not be a constant. 3.1.5 A New Life for Orlov's Recovery Operator In this section equation 3.7 is reformulated so that it can be used i n the case where xp < TT/2. The first step is to rewrite the inner integral of equation 3.7 as a two 37 Figure 3.7: Truncated spherical detector and projection plane geometry 38 dimensional convolution, that is / dS' ''*\ = p(xT;f) * *—- _ J \xT \L{r x xT) \XT\L(T x xT) (3.9) The outer integral of equation 3.7 is a backprojection into a volume of the two dimensional result of the inner integral. The entire equation can be then rewrit ten as the Laplac ian of a three dimensional backprojection of a two dimensional convolution, 1 G p(xT;r) * *—f xT \L(T X XT) After bringing the laplacian operator inside the integral, the inner term becomes 1 p(xT-r) * *—77 \xT \L(T x x T ) . Now since for any two functions / ( x ) and g(x) [7], V 2 [/(x) * *g(x)] = [V2f(x)]**g(x) = / (z )**[V 2 <7(*) ] equation 3.7 can be rewrit ten as, p(x) = J dGT [p(xT;r)**f(xT;r)] where / ( x T ; f ) --(3.10) (3.11) 4-7T2 \\XT'\L(T X xT)J ' Since the projections are first convolved wi th f(xT; f ) before being backprojected, f(xT; f ) is referred to as the reconstruction filter. Equations 3.10 and 3.11 form the basis of the new algori thm as follows; the desired source density dis tr ibut ion, p(x), is recovered by first filtering the measured two dimensional projections and then back-projecting the result into the three dimensional reconstruction volume. Th i s algorithm is analogous to the standard two dimensional image reconstruction method based on 39 equation 3.2. In two dimensional reconstruction, filtering is done by convolution over one linear dimension and the backprojection is done over one angular dimension. In three dimensional reconstruction, filtering is by convolution over two linear dimen-sions and backprojection is over two angular dimensions. Us ing this method of three dimensional image reconstruction, p(x) can be recovered without the use of Fourier transforms, and the advantages of convolution filtering over Fourier transform filtering i n the two dimensional case are expected to exist i n the three dimensional case. 3.1.6 Analytic Filter Form In order to use equation 3.10 as a reconstruction algori thm, it is necessary to find the explicit form of the reconstruction filter given by equation 3.11. In appendix B it is shown that the function L(T X XT) can be expressed i n the coordinates of figure 3.4 and figure 3.7 as, L(T X xr) = L(d,;$r) = < 2 arcsin [ , s ' n ^ ) , otherwise 7r, i f | cos 9i I > cos ip11 sin 8r sint/> \Jl — cos 2 9\ s i n 2 9T (3.12) In the same coordinate system, the filter function, f(xT; f ) , of equation 3.11 becomes, {mm) <3-13) and using the the cyl indr ica l coordinate system (l,6i,t) w i th the Laplac ian operator given by, idi Ydi) + Pde? + at equation 3.13 now becomes / ( W r ) = 4 T T 2 / 3 1 d2 + L(9,;0T) d9t2 \L(6r,9T) (3.14) 40 The explicit form of f(l,9i;9T) is obtained by substi tut ing equation 3.12 i n for L(6f, 9r) i n equation 3.14. The expansion of the second derivative is done i n appendix C and is too lengthy to present here. However several features of the two dimensional reconstruction filter are already apparent upon inspection of equation 3.14. The filter has an overall — l / / 3 dependence that is modified by a term that varies only as a function of the polar angle 9\. There are also two different types of singularities present, one at the origin where / = 0, and the other along the lines where | cos#;| = cost/^/l s in# T | . These singularities w i l l be discussed i n section 3.2. The filter also has no variat ion wi th respect to <f>T, which is expected because of the symmetry of the detector about the 2-axis as shown i n figure 3.4. The filter is also symmetr ical about 9T = 7r/2 and 9\ = 0 ,±7r /2 ,7 r , and both of these symmetries w i l l be used later to simplify the calculat ion of the digi ta l reconstruction filter. 3.1.7 Limits of Integration The last theoretical hurdle remaining before equation 3.10 can be used is to determine the appropriate l imits of integration. W r i t i n g out the two dimensional convolution expl ic i t ly i n rectangular coordinates gives, where f(lx, ly; 9T) is obtained from equation 3.14 by converting from polar to rectangu-lar coordinates i n the projection plane. In the above equation the l imits of integration are ± o o but the — l / / 3 term i n equation 3.14 means that i n a pract ical implementation it is possible to neglect contributions for lx' and / y ' larger than some pre-determined The integration over the region G is a double integral over the two spherical angles p(lx, ly; 9T, <j>T) * *f(lx, ly-8T) = max-41 that describe f . ' J ddr d(f>T sin 8T G It is not immediately obvious what the l imits of 9T and (j)T are. Since, however, the inner integral of equation 3.10 is a convolution of the projections w i t h a spatially invariant filter, and from the discussion i n section 3.1.3, it follows that the P S F of the detector i n figure 3.7 must be made spatially invariant. For an object centered at the origin that does not extend further than a distance R, and from the discussion preceding equation 3.6, the effective acceptance angle is given by V'min = V> + arcsin(i2cos ij>/R^). Th i s means that the l imi ts of 6T are TT/2 ± */>min. T h e l imi ts of <f>T can be determined from figure 3.7 by not ing that where lx and ly are the reflections of the coordinates of lx and ly about the origin. In other words, a l l the projection data is contained, without dupl icat ion, i n a set of projections described b y 3 , Equa t ion 3.15 is the basis for an algori thm that uses the same cross-plane data as the three dimensional reconstruction methods based on Fourier transforms that were discussed i n section 3.1.3. However the derivation of the two dimensional filter and the reconstruction of the three dimensional image occur entirely i n object space. Th i s means that the use of D F T s can be avoided. 3 A n equivalent choice is n/2 — ip < 9T < TT/2 and 0 < <j>T < 2TT. p(lx, ly', #t, 4>r) = P(lx, ly', K — 0T,7T + (j)T) T T / 2 - 0 < 0T < TT/2 + l/f, 0 < (f>T < 7T In its full glory, equation 3.10 now becomes, (3.15) 42 3.2 Discrete Form In any pract ical si tuation, variables that are continuous i n theory can only be measured and manipulated i n a discrete manner. The projections p(lx, ly;6r,</>T) are measured as averages over bins (or pixels), so the integrals of equation 3.15 must be approximated by a set of discrete summations, and the reconstructed estimate of p(x) must be presented as a three dimensional array of values (or voxels). Since the algori thm w i l l manipulate discrete data, equation 3.15 must be rewrit ten to reflect this. In the three remaining subsections of this chapter, 3.2.1 lays the groundwork for the discrete image reconstruction, 3.2.2 develops the backprojection algori thm, and 3.2.3 contains the calculat ion of the digi ta l reconstruction filter. 3.2.1 Discrete Reconstruction The first step i n the discrete image reconstruction is to express the integrals i n equa-t ion 3.15 as a set of sums of integrals over each b in . = Ufa A O l 9 ddrSmdrY^J dcf>T (3.16) • E E , dlj I , dly'p{lx-lx',ly-ly'-dTAr)f{UAy';eT) , i, ji Jli -d/2 Jlj -d/2 where W and 1/ are the centers of the bins i n the projection plane and 9m and (f>n are the centers of the projection direction bins. E a c h b i n is referenced by four indices, the first two, %' and j', describe the location of the b in i n a projection plane, and the last two indices, m and n , describe the orientation of the projection plane. If each b in is of square size d x d i n the projection plane and each projection plane is separated by 43 A9T i n 9T and A<j>T i n <j>T, then the region of each integration i n equation 3.16 is li - d/2 < lx < U + d/2, where U = id lj - d/2 < ly < lj + d/2, where lj = jd 0m-A9T/2 < 9T < 9m + A9T/2, where 9m = mA8r <f>n — A(f>T/2 < <f>T < (f>n + Acf>T/2, where <f>n = nA<f)T The m a x i m u m and m i n i m u m values of each index are determined as follows, ^max = Jmax^ = jmax.d ?min = J min = - 2 r a a x V'min = r n m a x A 0 T  m m i n = _ m m a x T = "maxA?i>T n m i n = 0 If 'max is finite and A9T, A(f>T, and d are not infinitesimally smal l , which means that the projections are discrete, equation 3.16 represents an estimate of p(x). If p(lx, ly; 9T, <f>T) is frequency bandl imited such that it is adequately sampled by bins of size A8T x A<f>T x dxd then it can be represented by an average value that is constant over each of the bins. In order for p(lx, ly; 9T, <f>T) to be bandl imited, it must have no spatial frequencies greater than some frequency fmax. The condit ion that must be satisfied for p(lx, ly; 9T, cj>r) to be adequately sampled is d < l / ( 2 / m a x ) , which is known as the Nyquis t sampling criterion [7]. The equivalent cri terion for sampling i n 9T and 4>T is that A9T, A<j>T < l / ( 2 7 2 / m a x ) where R is the radius of the object [8]. It is not possible for any object of finite size to be bandl imited, but p(lx,ly;9T,(f>T) could be bandl imited i f i t were possible to apply a low pass filter to p(x) before the projection samples were made. 44 If the above cri teria are met, then it is possible to define an exact band l imi ted reconstruction of p(x), called p(x), N o w since p(lx — lx\ ly — ly; 9r, <j>r) is constant over each of the dlx' and dly integrals and only depends on the indices i — i', j — j', m, and n , then <f>n+A<t>T/2 A(j>r/2 d<f>7 (3.17) rli+d/2 rlj'+d/2 E E P i - i ' - H w l J / 0 dl* h> J t dly f(lx\h''^r) •, j, Jli —d/2 Jlj —d/2 where p , j ; m > n are the bandl imited measured projections for each value of m,n). Since f(lx, ly', 9T) can be determined apriori , as w i l l be shown i n section 3.2.3, equation 3.17 can be rewriten as, y ' ^ J9m-A0T/2 TZ?Un-Aj>T/2 12 Pi-i'J-j';m,n fi,j;(&r) (3.18) where fi,j;(9T) is a two dimensional digi ta l filter given by, h i m = Ju,_ ' di h'+d/2 flj'+d/2 (3.19) The evaluation of the above equation is performed i n section 3.2.3. The double sum inside the square brackets of equation 3.18 is a two dimensional discrete convolution filtering and can be wri t ten as, Pi,j;m,n = Pi,r,m,n * *fij.(0T), w h e r e p * J ; m n is the filtered projection and it is understood that the convolution is over the two indices i and j , and f is now a discrete function of m and n. 45 3.2.2 Three Dimensional Backprojection The sums over m and n i n equation 3.18 are a three dimensional backprojection of the discrete filtered projections, p*j.mtn, where i and j are functions of x, m , and n . The projections are constant over each d9T and d<f>T integral and / ; j ; ( 0 T ) has no <pT dependence, so that each integral over d<f>T reduces to r<f>n+A<f>T/2 J<t>n-A<pT/2 Similar ly, i n practice the interval A9r is smal l enough so that sin# T and fi,j-,(9T) are slowly varying over each interval and so, f$m+AeT/2 L A a dBr s i n °r = A9T sin 0m. J9m-A9T/2 It is possible to now write equation 3.18 as P(*) = E E [P«-J;m,n * */.,i;m] , m n where /,-,j;m = A<f>TA9Tsm9mfitj.(9T). To complete the description of the backprojection process it is necessary to convert p(x) to a discrete form. The reconstructed image is manipulated by a computer so it is stored as a three dimensional array of discrete values, ratb>c, where r is the value of the image at the point (or voxel) indexed by the indices a, b, and c. The idea behind backprojection is that for each point x the value of p(x) is given by the sum of a l l the contributions to that point from al l the filtered projections. This is i l lustrated for two dimensional backprojection i n figure 3.2 There are several different methods used for backprojection, each w i t h its own mer-its and disadvantages, and the trade-off is the usual one of accuracy vs. computat ion t ime [9,42]. The method that was chosen for the current image reconstruction algo-r i t h m is as follows; the coordinate of the center of the voxel i n the three dimensional 46 reconstruction volume, r a )(, ) C, istransformed to the coordinates lx and ly on the (m, n) projection plane and the value of the filtered projection at that point is added to r 0 )(, ) C. T h i s is then repeated for a l l filtered projection planes. The coordinate transformation is i l lustrated i n figure 3.7, where x is the point to be reconstructed, $ T and <j)r describe the orientation of the filtered projection plane, and I = (lx,ly) is the transformation of x onto the plane. The value of the filtered projection at the point / is to be added to r a ( ) j C (pointed to by if). G iven the point (x,y,z) that is indexed by (a ,6,c) and the direction vector (0T, (f>T\ the needed coordinates (lx, ly) are found from, lx = —x sin <J>t + y cos <f>T ly = — x cosdT cos <f>T — y cos6T sin<^T + z s'mdT Since the filtered projection plane, p * j . m i n , is a discrete two dimensional histogram and it is unlikely that the transformed coordinate (lx, ly) w i l l be at the precise center of a histogram b in , some means of interpolation is required. Three of the most common methods are, i n order of increasing accuracy and computation time, nearest neighboor, bil inear interpolat ion, and bicubic spline. The bilinear interpolat ion method was cho-sen as a good compromise [40] between speed and accuracy, and is the method used i n the reconstruction program listed i n appendix E . The bil inear interpolation method [75] is shown i n figure 3.8, where the values Pf, -P2*, P 3*, and PX are located at the centers of four adjacent bins (separated by unit distance) of p * j ; m n for a given m and n. The four points shown are those that surround the point (lx,ly). Us ing bilinear interpolation, the value of the filtered projection at the point (lx,ly) is given by, P\lx,ly)m,n = (1 " dx)(l - dy)Pt^n + dx(l - dy)P2*>mtn +(i-dx)dyp;.mtn + dxdyp;.mtn A summary of the backprojection algori thm used i n the reconstruction program of + Figure 3.8: Bilinear interpolation 48 appendix E is as follows; h For each plane p * j ; m i „ (indexed by m and n) loop h For each voxel r 0 ) i ) C (indexed by a, b, and c) loop 1 Given (m,n) and (a, 6, c), transform a? onto the projection plane to get 2 Given (lx,ly) find the interpolated value p*(lx, ly)m,n H End loop H End loop 3.2.3 Calculation of Digital Filter W h a t now remains, i n order to complete the discrete version of the reconstruction algo-r i t hm, is the calculat ion of the digi tal filter i n equation 3.19. The transformation from the polar coordinates i n equation 3.14 to the rectangular coordinates i n equation 3.19 is given by, lx = / cos Q\ ly = / sin 9i The difficulties encountered i n calculating fi,j-,{QT) from equation 3.19 are due to the last term i n equation 3.14, because of the discontinuities i n the first derivative along the lines 0/ = ±a rccos (cosV>rn in / s in# T ) . To see this i n the region 0 < 9\ < TT/2, rewrite l/L(0i;dT) as TO = ^ " < * - * t a k r ; ) ( 3 - 2 0 ) 49 where 0 = arccos(cosV'min/ sin# T) is the angle where L(0i;0T) changes form, u(x) is the Heaviside step function definied by ' 1, i f x > 0 0, i f x < 0, u(x) = < and L(9f,9T) = 2arcs in sin tp-mm t ^ / l — cos 2 0/ s i n 2 #T which is defined to be 7r for 8, <8. N o w using the relationship [54] — {u(x - a)f(x)} = u(x - a)faf(x) + *>(x ~ «)/(«) the first derivative of l/L(8i;9T) can be calculated as, and since L(8; 8T) = 7r, the last term is identically zero everyhere. For the second derivative, d2 1 08,2 \L(8,;8r) = 6(9, - 8) d8,L(8r,8T) Using the above result i n equation 3.15 yields 1 (I f(l,0r,Or) = 4 7 T 2 / 3 \7T +6(0t - 0) + u(0i - 0) d 1 ldO,L(0r,Or)\0lJ§ 1 \ \ Th i s means that the digi ta l filter fi,j-,(8T) can be given as the integral of / ( / , 0j; 8T) over a square (d x d) p ixe l . rli+d/2 rlj+d/2 / y s ( * r ) = A „ d/, A dlyf(lA;0T) Jli-d/2 Jlj-d/2 (3.21) 50 where I2 = lx2 + ly2 and 0/ = arc tan({ y / / x ) . To separate the integral i n equation 3.21 into manageable parts, /,-,j ;(0T) is treated as the sum of three separate integrals, = IliJt{0r) + I2ith(dT) + / 3 , J ; ( 0 T ) , (3.22) where 1 fU+d/i ftj+d/2 i Aj /  47T 3 JU-d/2 Jlj-d/2 ^ V I3 anc 1 rli+d/2 rlj-U+d/2 rlj+d/2 dlv I. dl lj-d/2 V I3 + 6(6, - 0) d 0i=0; d8i L(8r,0T). The two terms I l t i j . (# T ) and J2, i j ;(e9T) w i l l be evaluated expl ici t ly below, but / 3 , ) j ; ( 0 r ) is too complicated for analytic integration, and must evaluated by numerical approxi-mat ion. E a c h term w i l l now be considered i n turn . For each term, the pixel of consideration is demarked by the four corners ( x i , y i ) , faiiSfe), ( x 2 , y i ) , and ( x 2 , y 2 ) , where xx = l{- d/2, x 2 = U + d/2 Vi = h-d/2, j/2 = lj + d/2. The direct evaluation of i l , j ; ( 0 T ) gives. Jl.\i;(*r) = y=V2 4?r3 xy (3.23) x—xl J y=y\ 51 In the case where x x — 0 then 4 T T 3 x2y V=V2 and i f t/i = 0 instead, then K M = 4 T T 3 y / * 2 + y2 2 XT/2 ar=xi For Z2j ij ;(0 T) and Z3i,j;(0 T), i t is necessary to consider the four different ways that the line 9[ = 9 can intersect the pixel being integrated over, as shown i n figure 3.9. In evaluat ing72, J ; (# T ) , it suffices to use one of the cases shown i n figure 3.9, which is enlarged upon i n figure 3.10, where the pixel is divided into three regions of integration. In region Ri the integrand i n the expression for / 2 , j ; ( 0 T ) is zero since 9\ < 9. For region i?2, I2ij.(9T) has the form, a r rr. ~ _n y=vt 1 rVb ryI tan U n dlv l7T 3 Jya V Jx, dlx ( x 2 + y 2 ) 3 / 2 / 4?r3 cos 9 x2y y where ya = x x tan 9 and J/& = x 2 tan 9. The functional form of I2ij.(9T) over the region i?3 is similar to i"l t i j ; (# T ), so the complete form of I2,j;(# T) for the case shown i n figure 3.11 is, f l X=X2 -i V=J/2^ 4 T T 3 x 2 y COS0 y y=yb y=y* Vx2 + y2 xy X—X\ (3.24) y=yb J The forms of J 2 t j . ( 0 T ) i n the other three cases shown i n figure 3.9 are similar to that of equation 3.24 w i t h minor variations to take into account how the line 9i = 9 intersects the pixel . The last and most difficult of the terms to evaluate is / 3 , J ; ( 0 T ) . F r o m the form of the second derivative of l/L(9i\9T) shown i n appendix C , it is obvious the the integrations needed to calculate i"3,j ;(#T) must be done numerically. Star t ing wi th 52 53 Figure 3.10: Regions of integration i n a pixel for I2 , j . (0 T ) . 54 fb I dx6(x — x0)f(x) = < Ja [25], / ( x 0 ) , i f £o G [a, 6] 0, i f x 0 £ [a, b] the term i n J3j,j ; (0 T ) involving 6(8i — 9) is straightforward to evaluate. It turns out, however, not to be necessary to do because of the following. B y analyzing the form of the second derivative of \/L{9i\ 9T), shown i n appendex C , it can be seen that, and l i m 1 - c o . [d9i2 \L(9 + t;9T)J This is i l lustrated i n figure 3.11. Because of the behaviour of the second derivative as 9\ —• 9+, any straightforward numerical integration technique w i l l underestimate the contr ibut ion of the second derivative term to 73,j ;(# T). In order to circumvent this problem, a different approximation approach that takes into account the behaviour of the second derivative and incorporates the above Dirac-del ta function was adopted. Since a f I \ o d9l\L{9l-9T)j at 9\ = 0 or 7r / 2 , it can be seen that d2 0 2 \L(9r,9T)Jj and since the second derivative term is zero for 9, < 9 and 6, = 7r/2, /•8—t+ot de,2 so that [9—e+ct d9i KL(9t d2 );0T))} /*-«+« {d9<2 {m;9T) 39S \L(9l-9T)i [d9t L(9t;9T) 9i=9—e+ct (3.25) 55 / / / Figure 3.11: Sketch of the second derivative of l/L(0f, 0T). 56 Thi s means that i f there is an angular integration range of a that includes the value 9, = 9, it is possible to replace the integral wi th a non-singular term that avoids dealing directly w i th the singularities of the second derivative. The difficulty now arises that the integration i n equation 3.21 is over a square pixel - but it is necessary to perform the integration on polar coordinates i n order to avoid the singularity at 9i = 9. Th i s difference i n coordinate systems can be reconciled by approximating a th in strip about the line singularity to be an annular segment as shown i n figure 3.12. Since I3ij.(9T) is to be calculated by a two dimensional numerical integration, such as Simpson's rule, the integrand would normally be evaluated at the gr id points shown, and, N N nutfr) * E E W s,t 1 / ,„ ^ ( 1 & 2 3=0 t=0 + 6(9, - 9) (3.26) 4 T T 2 / 3 V V"' ~J\L(9,;9T) ' d9,2L(9,;9T) d 1 [d9, L(9r>9T)\ where h = d/N is the spacing between gr id points i n the pixel , d = x-i — x% = y2 — j / i , ws>t are the appropriate quadrature weights, and / and 9, come from lx = x x + sh and ly = yi + th (so that I2 = I2 + ly2 and 9, = axctan(ly/lx)). B u t because of the nature of the second derivative of \/L(9,\9T) as shown by figure 3.11, the accuracy of the integration i n equation 3.26 would depend i n an arbitrary way on the spacing of the gr id points near 9, = 9. Th i s is the underestimation of J3 X ) j ; (0 T ) mentioned above that can now be corrected for by using equation 3.25 to calculate the contribution to the integral from those points i n equation 3.26 that fall inside the strip shown i n figure 3.12. The narrower the strip, the smaller the deviations from a true annular segment are. For the annular segment h < I < l2, 9 — e < 9, < 9 — e + a, where / i = Xi cos 9 and l2 = x2 cos 9, the part of the sum i n equation 3.26 for gr id points 57 Figure 3.12: Integration scheme over a pixel for I3 t i j ; (0 T ) . 58 1 W ^ T 1 47T2P \ v ' /d9,2 L(9,-6T) + 8(9, - 9) * _ _ L L ° ~ £ + a <» I a 2 _ 1 ~ 4 ^ 4 P J9-t 1 \d9,2 L(9,;9T 1 ' ' follows, ' d 1 y d9iL(9r,9T)_ L 4 7 r 2 / J i = f l [ ^ ( M O N o w by let t ing e —• 0, the final form for J 3 , J ; ( 0 T ) is, 1 (u(9,-9) 9,=9—e+a N N s=0 t=0 4 T T 2 ; W2 4 T T 2 / 3 V ^ ( ^ ; ^ T ) a I + u(8, -(9 + a)) d2 d9,2L(9,-9T)t (3.27) d9,L(9,;9T)\d[=d+a where it is understood that w'st — 0 for s and t inside the segment. The accuracy of the numerical approximation depends asymptotical ly on the value of N w i t h the usual trade-off between speed and accuracy. F r o m equation 3.25 it is seen that equation 3.27 should be relatively insensitive to the value of the parameter a as long as a is smal l enough that the segment is close to being annular. A l so , from the geometry of figure 3.12, a must be larger than d/l\ so that it is guaranteed that some gr id points w i l l lie i n the segment. In trials, by evaluating equation 3.27 for different values of a, it was indeed found that for N = 50, J3 t i j ; (0 T ) d id not change significantly for d/lx < a < 5d/lv Beyond TV = 30 the value of ^ J2j / « ' , J ; ( ^ T ) changed l i t t le w i t h increasing N, but even N = 30 results i n evaluating first term of equation 3.27 at 900 gr id points for each pixel . To reduce the computat ional burden it is possible to use much smaller values of N for pixels i n regions where f(l,9,;9T) is slowly varying. The gridding algori thm adopted was to use JV = 8 unless the p ixe l was w i th in two pixels of the or igin or the line 9, = 8, i n which case N = 30 was used. It was found that this gave the same values for IZitj.(9T) as when N = 30 was used for a l l pixels. For a 41 x 41 digi ta l filter this 59 method reduced the computation t ime by ~ 80%. The weights w3tt i n equation 3.27 were the ones from the extended Simpson's rule [75]. N o w by using equations 3.22, 3.23, 3.24, and 3.27 it is possible to compute fi,j-,(8T) for i,j, > 0. To calculate fij;(8T) for i and /or j < 0 it is only necessary to observe the two-fold symmetry of L(8i;8T), i n equation 3.12, about 8, = 0 , ± 7 r / 2 , a n d TT. The values of fi,j-,(8T) for such i and j are then found by reflection of the values for i and/or j > 0 about the lx and/or ly axes. If one of i or j is zero then / ( / , 8,; 8T) is symmetric about the axis through the middle of the pixel and fitj.(9T) only has to be calculated for the port ion where 0 < 8[ < TT/2, and then doubled. The last part of fi,j;(8r) to be calculated is the origin, or /O,O;(0T)- Because of the — 1/P term i n f(l,8i;6T), it is not possible to directly calculate /o,o;(0i-)- The center pixel can be obtained by imposing the condit ion, /oo too dlx / dlyf(lx,ly;8T) = 0 -oo J—oo Thi s follows from the property that the integral over a l l space of the filter represents the value of the Fourier transform of f(lx,ly;8T) evaluated at zero spatial frequency [7], and by comparison w i th the frequency space filters i n chapter 5, the filter must be zero at zero frequency [45]. F r o m this comes the relationship, £ £ A ; ; ( 0 r ) = o or fo,o-,(8T) = - E £ A i ; ( ' r ) (3-28) («J)^(o,o) Figure 3.13 shows a histogram plot of the two dimensional digi ta l filter, fi,j-(8T), for 8T = 90° and tp = 1 0 ° 4 . Because of the magnitude of /o,o ;(0T) no detail is visible, so i n figure 3.14 the center 9 x 9 pixels are set to zero and the histogram is flipped upside 4 A value of xp appropriate to the new PET machine 60 Figure 3.13: Digital filter for 0T = 90° and ip = 10' 61 Figure 3.14: Negative of the filter w i th the center 9 x 9 pixels set to zero 62 down (by mut ip ly ing it by —1) i n order to show more detail . The ridges at 9, = ± 1 0 ° are due to the discontinuity at the same angle that is shown i n figure 3.11 that is i n t u rn due to the discontinuity i n slope of the function L(T x xT) i n equation 3.12. B y looking at fij.(9T) vs. xp and 9T it was noticed that the overall shape of fi,j;(9T) depends more strongly on xp than on 9T. Unfortunately it is difficult to develop an intuit ive understanding of how changes i n the shape of the filter affect the reconstructed image. It is apparent that the calculation of the digi ta l filter is non- t r iv ia l and time con-suming. However once the detector parameter xp and the projection angles 9m and <pn are known, fi,j;(9T) need only be calculated once and then stored i n computer memory. The computer program used to calculate fij-(9T) is l isted i n appendix D . 63 > Chapter 4 Implementation This chapter reports on the results of using the three dimensional reconstruction pro-gram listed i n appendix E . M i n o r modifications were made to the program i n order to provide a simple test case for the reconstruction algori thm and to speed up the test cycle. 4.1 Reconstruction Details The object chosen to provide a test for the algori thm was a sphere of radius R w i th a uniform density, ' 1, i f |x | < R 0, i f \x\ > R T h e projections of a uniform sphere at any angle are given by, p(x) = p(xT;r) = p(l,di;8r,<f>T) = 2y/R2 - |/~12, if \1\<R 0, i f |?| > R so the discrete projections are, ' 2y/R2-(i2+j2)d2, i f i2 + j2 < (R/d)2 Pi,j;m,n — * o, i f i 2 + p > (R/dy (4.1) 64 where i and j are the pixel index numbers on the projection plane. The parameter V'min w a s chosen to be 10° and projections of the sphere where taken at, where A6T = V'min/S and A<f>T = TT/60. The reconstruction volume was 40 x 40 x 40 voxels, each of unit size, and the sphere had a radius of 10 units. The projection planes were 40 x 40 pixels of unit size, and because of the symmetry of the sphere a l l the projections were identical. The chosen values for V'min #m were used by the filter generation program listed i n appendix D to generate a 41 x 41 digi ta l reconstruction filter. In order to reduce aliasing artifacts, it is necessary to bandl imit an object before sampling it . T h i s is not possible wi th a P E T scanner because the data are acquired discretely. Th i s means that the projections are sampled before any bandl imi t ing can be done, al though a low-pass Hamming type filter w i l l help to reduce any resultant aliasing. Instead of first low-pass filtering the projections, since the filtering operations are linear it is equivalent to low-pass filter the reconstruction filter,/,j;m, and then use this composite filter on the projections. Th i s shifts some of the computat ional burden from the reconstruction program to the filter generation program; thus saving time during an image reconstruction. To construct the two dimensional low-pass filter, the most straightforward method is to calculate the equivalent one dimensional filter and use it to form a rotat ionally symmetr ic two dimensional filter [37] that has a frequency space form of [75], where fn is the Nyquist frequency l imi t given by / „ = 1 / (2A) , A is the sampling interval , and the rect() function is defined on page 33. The first term of the R H S of 6m = 71-/2 + m A 0 T , m = - 3 , . . . , 3 (7 angles) <f>n = nA(j>r, n = 0 , . . . , 59 (60 angles) (4.2) 65 equation 4.2 is the low-pass filter term that abruptly cuts off any frequency greater than / „ to avoid aliasing. The second term on the R H S is the Hamming window function that causes W(f) to ro l l off smoothly to zero as / — > • / „ , i n order to reduce r ippl ing effects i n object space [17], and may also help to reduce aliasing that is al-ready present. A l though the choice of the window function is arbitrary, there is l i t t le difference between the effects of windowing functions wi th similar functional forms. To find the object space version of equation 4.2, it is straightforward to inverse Fourier transform W(f) to get, w(x) = / „ ( s inc (2 / n x) + l / 2 s i n c ( 2 / n x - 1) + l / 2 s i n c ( 2 / n x + 1)) and since the pixels have unit size, A = 1 and / „ =. 1/2. To make the filter two —# dimensional v i a rotat ional symmetry, x is replaced wi th / = |/ | , w(l) = ^sinc(Z) -(- j s i n c ( / - 1) + js inc(7 + 1) (4.3) & ft fr Th i s low-pass filter is used to filter each of the two dimensional reconstruction filters. 4.2 Results A n ideal three dimensional reconstruction of the sphere would return an object that could be described by equation 4.1. Since a band-l imited reconstruction is being performed, sharp boundaries cannot be recovered. W i t h this i n m i n d the expected ideal reconstruction is an object that has uniform density for |x | < R — 8/2, is zero for |x | < R + 6/2, and has a smooth transi t ion between the two regions over the range R — 6/2 < \x\ < R + 8/2. The wid th , 6, of the transi t ion region depends on the nature of the low-pass filter used (equation 4.2). A histogram plot of the x-y plane of the band-l imited reconstruction, p(x), of the sphere is shown i n figure 4.1. The z-x plane histogram is shown i n figure 4.2. Values 66 Figure 4.1: The x-y plane of p(x). 67 Figure 4.2: The z-x plane of p(x). 68 of p(x) along the x-axis are shown i n figure 4.3 and the values along the z-axis are shown i n figure 4.4. F r o m figure 4.3 it can be seen that the transi t ion wid th is 6 ~ 4. The values on the x-y plane i n figure 4.1 are seen to be rotat ionally symmetric so that the values of p{x) along any line i n the x-y plane passing though the or igin are identical to those shown i n figure 4.3. The value of p(x) does not go completely to zero for |x | > 10 and this is a consequence of equation 3.28. Since the reconstruction filter is not infinite i n extent, the filter w i l l have a slightly non-zero value at zero frequency that leads to an overall offset i n the reconstructed image. In a l l other respects the results shown in figures 4.1 and 4.3 are as expected. In figures 4.2 and 4.4 there is a deviation from circular symmetry. The largest manifestation of this lack of symmetry is a 17% undershoot of the reconstructed values on the 2-axis, just outside the radius of the sphere. Th i s undershoot is an unexpected artifact that is also believed due to the manner of implementation of equation 3.28, because the contributions to /O,O ;(#T) from / , J ; ( # T ) outside the 41 x 41 range of the filter were neglected. Because of this, the undershoot artifact is believed to be an implementat ion error rather than a theoretical one. Th i s is borne out by the results of chapter 5, where it is shown that the algori thm is functionally equivalent to analytic Fourier transform methods that do not show such an artifact. 69 ) 1.5 1.0-0.5-o . o -.5 Figure 4.3: Values of p(x) along the x-axis 70 < { 1.5 1.0-0.5-0.0H - . 5 - 2 0 - 1 0 Figure 4.4: Values of p(x) along the z-axis 71 Chapter 5 Equivalence to Fourier Methods T h e three dimensional reconstruction algori thm i n chapter 3 is derived i n , and only op-erates i n , object space. Previously developed analytic three dimensional image recon-struction methods rely on Fourier transform methods [97,70,21,104,95,26,77,71,22,86,76] In order to see how the object space reconstruction filter compares to the earlier fre-quency space filters, the two dimensional object space filter is Fourier transformed into its frequency space form. 5.1 Fourier Transform of the Reconstruction Filter To begin, F(R, 9T) is defined as the two dimensional Fourier transform of the filter function, f(l,9i;9T), that is defined by equation 3.14 where R and $ are the polar coordinates on the transformed plane. F(R,^,8T) = ^2D{f(l,9i;9T)}. N o w since for any two dimensional function g(r,9) [7], T2D { V M r , 0 ) } = - 4 7 r V 2 ^ 2 D {g(r,9)} 72 the Fourier transform of equation 3.13 is given by F(R,$;$T) = R2F2D {lL(e,;9T)} and wr i t ing out the transform explici t ly i n polar coordinates [7], ,2* roo -2-KiRlcOs(6l - <5) F(RWr) = R°l M,Jo 141 . (5.!) Because of the form of the function L(9,\9T) i n equation 3.12, it does not seem possible to directly perform the integrations of equation 5.1. Instead, the propeties of the Dirac-del ta function are used to simplify the evaluation of F(R,$;9T). To start w i th , equation 5.1 is rewritten as r /"k f2ir i /-oo P—2mRlcos(9i — $) F(R, $ ; 9T) = R2 \J d9i + j d9,\ J dl — so tir re F(R,$-9T) = R2 / d9, / Jo Jo dl L(9,;9T) -2mRl cos(9i - $) e-2iviRl cos(9, + T T - $ ) + L(9,;9T) L(9, + 7r;9T) F r o m equation 3.12 it can be seen that L(9, + TT; 9t) = L(6i\ 9T), and since cos(0/ + 7r — $ ) = — cos(0; — $ ) , the second integrand above can be rewrit ten so that, r-K roo F(R, $; 9T) = R2 / dd, / dl Jo Jo e-2mRlcos(9i - $) e2-KiRl cos(9, - $) + L(9,-9r) L(9,-9T) Now by setting / —• — I m the second integrand, so that dl —• — dl • di e-2iriRicos(9' - - d i e - 2 7 v i R U o s ( e ' -$) U p o n reversing the l imits i n the last integral and combining the two integrations wi th respect to /, the result is, F(R,$;0T) = R2 [d9, ^Q ) J°° dl e~2niR1 co<9' " (5.2) 73 Recal l ing the integral form of the Dirac-del ta function defined by [25] 6(x) = f ° d t e ~ 2 v i x t . J—oo Thi s means, H dl e - 2 7 c i R l co<9> - * ) = T^S(cos(e, - $ ) ) . (5.3) J-oo \R\ where R is the polar radia l coordinate so \R\ = R. The term <5(cos(#; — $) ) can be expanded to a simpler form by the relationship [25], *(/(*)) = E ^ ( * - * < ) , where f'(x) is the derivative of f(x) and X{ are the simple zeros of f(x). Note that if any of the zeros of f(x) are greater than first order, then the expression 6(f(x)) has no meaning. For 9, 6 [0,7r], the zeros of cos(0/ — $ ) are 8, — $ ± 7r/2, and since | - s i n ( ± 7 r / 2 ) | = 1, <5(cos(0, - $ ) ) = 6(6, - $ - T T / 2 ) + 8(6, - $ + T T / 2 ) , and equation 5.3 becomes f00 dl e ~ 2 7 r i R l cos(d' - * ) = 1 (<5(0, - $ - T T / 2 ) + £(0, - $ + T T / 2 ) ) . Us ing the above result i n equation 5.2 yields Recal l ing that [54], t dx f(x)8(x - x n ) - l / ( X 0 ) ' i f X° E [°' b]  J a I 0, i f x0 £ [a, b] the first integral i n equation 5.4 is Joddl L(9,-9T) - \ L($+1/2;9T)> * - T / 2 < * < * / 2 0, otherwise. 74 Simi la r ly for the second integral i n equation 5.4, f MW ~ * + * / 2 ) _ L{*-\/2-eTy Jo 1 L($r,eT) I 0 ) otherwise i. N o w once again using the relationship L(9, + n; 9T) = L(9i; 9T), the above two results can be combined to get, Subst i tu t ing i n equation 3.12 for L(9i;9T) yields the frequency space form of the two dimensional reconstruction filter, F(R,$;9T) = { ' R Y, i f sin $ > cos i / > m i n / s m R 2 arcsin sinV> m in -, otherwise (5-5) , \ A - s i n 2 $ s i n 2 0 T / 5.2 Comparison to Frequency Space Filters T o see how the frequency space form of the two dimensional reconstruction filter compares to previously derived frequncy space filters, it is necessary to briefly review Fourier reconstruction methods. These fall into two classes, two dimensional filtering and three dimensional filtering. In two dimensional filtering the a lgori thm is as follows; the projections are Fourier transformed to frequency space and are then filtered by mul t ip ly ing them by the appropriate two dimensional filter. A t this point the three dimensional image can be reconstructed i n two ways, either by inverse transforming the filtered projections to object space and then backprojecting them, or by first •placing the filtered projections into the reconstruction volume (in frequency space), and then inverse transforming the three dimensional result [95]. The equivalence of these last two steps is a consequence of the Fourier-slice theoem that states; i f E(u,v,w) is the three dimensional Fourier transform of p(x,y,z) such 75 that Z(u,v,w) = F3D{p(x,y,z)} and i f P(lu, lv; 9T, <f>T) = T2D {p( / x , ly; 9T, (j>T)} is the two dimensional filtered projection, where the relationships between (lu,lv) and (u, v, w) are analogous to those between (lx, ly) and (x, y, z) i n equation 3.8, then P ( / u , lv; $T, (f>r) = H(—lu sin<j>T — lv cos <f>T cos 8r, lu cos <^>T — /„ sin <f>T cos # r, /„ s in0 T ) . Stated i n words, the Fourier transform of a projection of an object is identical to a slice, at the same orientation as the projection, through the three dimensional Fourier transform of the object [71]. For three dimensional filtering methods the order of filtering and backprojection is reversed. Fi rs t the projection data is backprojected into the three dimensional reconstruction volume, where it is Fourier transformed into frequency space and then filtered by mul t ip l ica t ion w i th a three dimensional filter. The result is then inverse transformed into object space to become the three dimensional reconstructed image. For a l l (except possibly one) of the published frequency space methods the two dimensional frequency space filter can be related to a three dimensional version by the Fourier-slice theorem. In other words the appropriate two dimensional frequency-space filter is a slice at the appropriate angle through the origin of the three dimensional filter. The possible exception to this is the composite filter derived by R a et a l . [77], where the analytic filter form is different from al l other two dimensional filters [70,26,86]. To construct a three dimensional filter from equation 5.5 it is only necessary to incorporate the spherical projection angle, 9T, and the planar polar coordinates R and $ into a three dimensional frequency space coordinate system. F r o m figure 5.1 it is seen that c o s © = s i n $ s i n # T [1], where 0 is the spherical polar coordinate. Now since \Jl — s i n 2 $ s i n 2 9T = | s i n © | and R is the three dimensional radia l coordinate, 76 ( Figure 5.1: T w o dimensional and three dimensional coordinates i n frequency space 77 the three dimensional version of the frequency space filter equation 5.5 is i R i f | 0 | > $ , i f | 0 | < $ 7T ' F(R, 0 ) = (5.6) T h i s is identical to the frequency space filter derived independently by Pelc [70], C o l -sher [26], and Schorr et a l . [86]. The filter derived by Pelc requires some minor trigono-metric substitutions i n order to show that it is the same as equation 5.6. Since equation 5.6 shows that the Recovery Operator of Orlov [68] is equivalent to frequency space reconstruction methods, it also independently confirms the val idi ty of the three dimensional object space image reconstruction algori thm. 78 Summary There is one major advantage to bui ld ing a volume imaging P E T system; that of in-creasing sensitivity while keeping resolution constant. The two pr incipal problems that such a system has to overcome are parallax error and the difficulty of incorpo-rat ing a l l the collected data into the reconstructed image. The first problem has been solved by a new position-sensitive detector developed at T R I U M F that has the abil i ty to measure depth-of-interaction information. A proposed P E T system incorporates 15 of these 15 x 15 x 2.5 cm detector units into a cyl indrical r ing to form a volume-imaging scanner that has the abil i ty to detect cross-plane gamma rays. Previous work i n three dimensional image reconstruction that can use a l l of the cross-plane information mostly falls into two classes; iterative reconstruction and ana-lyt ic reconstruction using Fourier transforms. Iterative algorithms have the two advan-tages of being able to use a l l cross-plane data and the abil i ty to use apr ior i knowledge of the object being reconstructed. The pr imary disadvantage of iterative algorithms is that they require too much computat ion time. The Fourier transform algorithms require much less computational t ime than itera-tive algorithms. B y using Fourier transforms, the requirement of a spatially-invariant scanner P S F means that the amount of collected data that can actually be used is l imi ted by the relative geometries of the P E T scanner and the object i n its F O V . In order to make Fourier transform algorithms computationally efficient, they employ F F T filtering methods that add more computational noise than convolution filtering 79 methods. The new algori thm presented i n this thesis was based on an extension of Orlov's Recovery Operator. In its original form the Recovery Operator could only reconstruct objects for which there were projections i n a l l directions. To put it into a form suitable for use by the new P E T scanner, it was first reformulated as a F B P method, analogous to standard two dimensional image reconstruction, and then extended to work w i th incomplete three dimensional coverage. B y analyzing the angular l imits of integration for the backprojection part of the algori thm, it was shown that the new algori thm can use the same amount of cross-plane data as the Fourier transform algorithms. A n extension to the a lgori thm has been proposed by my colleagues and me [51,81] that allows the use of a l l the cross-plane data. Th i s extension is suited for three dimensional image reconstruction algorithms that operate i n object space. In the discrete form of the algorithm, a band-l imited form was created that dealt w i th the extensive singularities present i n the analytic form of the reconstruction filter. The calculat ion of the digi ta l filter is complicated and time consuming, but it only has to be calculated once and then stored i n computer memory for retrieval during the reconstruction of an image. The new algor i thm was tested by reconstructing a sphere from perfect projection data that was low-pass filtered wi th a Hamming window,. The sphere was faithfully reconstruced except for a 17% undershoot i n the reconstructed value on the 2-axis just past the edge of the sphere. This artifact is believed due to the method of normalizing the filter at its centre. Th i s algori thm has not yet been tested w i th data that has noise added. To demonstrate that the overshoot artifact is not an error i n analytic form of the algori thm, a result was derived showing for the first t ime the functional equivalence of 80 Orlov 's Recovery Operator and the Fourier transform algorithms. 81 ) Bibliography [1] CRC Standard Mathematical Tables. The Chemica l Rubber Company, seven-teeth edit ion, 1969. [2] Med ica l imaging for the masses. Staff Ar t i c le i n New Scientist, page 33, 7 M a y 1987. [3] M D Altschuler , Y Censor, G T Herman, A Lent , R M Lewi t t , S N Srihar, H Tuy, and J K Udupa . Mathemat ica l apects of image reconstructions from projec-tions. In L N K a n a l and A Rosenfeld, editors, Progress in Pattern Recognition, pages 323-375, Nor th -Hol land Publ ishing Company, 1981. [4] M D Altschuler and G T Herman. Fully-three-dimensional image reconstruction using series expansion methods. In A B B r i l l , R R Pr ice , W J M c l a i n , and M W Landy, editors, Fifth International Conference on Information Processing in Medical Imaging, pages 125-140, January 1978. He ld at Vanderbuil t Univer-sity, Nashvil le , Tennessee. June 27 - Ju ly 1 1977. Publ ished by the Biomedical Compu t ing Technology Information Centre at Oak Ridge Nat ional Laboratory, Oak Ridge, Tennessee 37830. [5] J E Bateman, J F Connolly, G J Stephenson, and A C Flesher. The Rutherford App le ton labori tory 's mark I mult iwire proport ional chamber. Nuclear Instru-ments and Methods in Physics Research, 225:209-231, 1984. 82 [6] M V Ber ry and D F Gibbs . The interpretation of optical projections. Proceedings of the Royal Society of London, Series A, 314:143-152, 1970. [7] R N Bracewell . The Fourier Transform and its Applications. M c G r a w - H i l l , second edit ion, 1978. [8] R N Bracewell and A C Riddle . Inversion of fan-beam scans i n radio astronomy. Astophysics Journal, 150:427-434, 1967. [9] R A Brooks and G D i C h i r o . Principles of computer assisted tomography ( C A T ) i n radiographic and radioisotopic imaging. Physics in Medicine and Biology, 32(5):689-732, 1976. [10] T F Budinger and G T Gulberg. Three-dimensional reconstruction i n nuclear medicine emission imaging. IEEE Transactions on Nuclear Science, NS-21(3):2-20, 1974. [11] C A B u r n h a m , D E Kaufman , D A Chesler, C W Stearns, D R Wolfson, and G L Brownel l . Cy l ind r i ca l P E T detector design. To appear in : IEEE Transactions on Nuclear Science, February 1988. [12] D B Calne, J W Langston, W R W M a r t i n , A J Stoessl, T J R u t h , M J A d a m , B D Pate, and M Schulzer. Posi t ron emission tomography after M P T P : Ob-servations relating to the cause of Parkinson's disease. Nature, 317:246-248, 1985. [13] R E Carson. A max imum likelyhood method for region-of-intrest evaluation i n emission tomography. Journal of Computer Assisted Tomography, 10(4):654-663, 1986. 83 [14] M E Casey and E J Hoffman. Quant i ta t ion i n positron emission tomography:7. A technique to redice noise i n accidental coincidence measurements and coincidence efficiency claibrat ion. Journal of Computer Assisted Tomography, 10(5):845-850, 1986. [15] Y Censor. F in i te series-expansion reconstruction methods. Proceedings of the IEEE, 71(3):409-419, 1983. [16] L T Chang . A method for attenuation correction i n radionuclide computed tomography. IEEE Transactions on Nuclear Science, NS-25(l) :638-643, 1978. [17] D A Chesler and S J Riederer. R ipp le suppression during reconstruction in transverse tomography. Physics in Medicine and Biology, 20(4):632-636, 1975. [18] D A Chesler, S J Riederer, and N J Pelc. Noise due to photon counting statistics i n computed X - r a y tomography. Journal of Computer Assisted Tomography, l ( l ) : 6 4 - 7 4 , 1977. [19] Z H C h o and J R Burger. Construct ion restoration and enhancement of 2 and 3-dimensional images. IEEE Transactions on Nuclear Science, NS-24(2):886-899, 1977. [20] Z H Cho , J K Chan , E L H a l l , R P Kruger , and D G McCaughy . A comparitive study of 3-d image reconstruction algorithms wi th reference to number of projec-tions and noise filtering. IEEE Transactions on Nuclear Science, NS-22:344-358, February 1975. [21] Z H Cho , J B R a , and S K H i l a l . True three-dimensional reconstruction ( T T R ) - App l i ca t ion of a lgori thm toward full u t i l iza t ion of oblique rays. IEEE Trans-actions on Medical Imaging, MI-2(1):6-18, 1983. 84 [22] G C h u and K - C T a m . Three-dimensional imaging i n the positron camera using fourier techniques. Physics in Medicine and Biology, 22(2):245-265, 1977. [23] R Clack , D Townsend, and M Defrise. A n algori thm for three-dimensional re-construction incorporating cross-plane rays. Submit ted to: IEEE Transactions on Medical Imaging, December 1987. [24] R Clack, D Townsend, and A Jeavons. Increased sensitivity and field of view for a rotat ing positron camera. Physics in Medicine and Biology, 29(11):1421-1431, 1984. [25] C Cohen-Tannoudji , B D i u , and F Laloe. Quantum Mechanics, pages 1468-1479. Volume two, Wiley-Interscince, 1978. [26] J G Colsher. Fu l ly three-dimensional positron emission tomography. Physics in Medicine and Biology, 25(1): 103-105, 1980. [27] M E Daube-Witherspoon and G Muehllehner. A n iterative image space recon-struction algori thm suitable for volume E C T . IEEE Transactions on Medical Imaging, MI-5(2):61-66, 1986. [28] M E Daube-Witherspoon and G Muehllehner. Treatment of axia l data i n three-dimensional P E T . Journal of Nuclear Medicine, 28(11):1717-1724, 1987. [29] A R de Pier ro . O n the convergence of the iterative image space reconstruc-t ion algori thm for volume E C T . IEEE Transactions on Medical Imaging, M I -6(2):174-175, 1987. [30] D J de Rosier and A K l u g . Reconstruction of three dimensional structures from electron micrographs. Nature, 217:130-134, 1968. 85 [31] M Defrise, S K u i j k , and F Deconinick. A new three dimensional reconstruction method for positron cameras using plane detectors. Physics in Medicine and Biology, 33(1):43-51, 1988. [32] R V Denton, B Friedlander, and A J Rockmore. Direct three-dimensional image reconstruction from divergent rays. IEEE Transactions on Nuclear Science, N S -26(5):4695-4703, 1979. [33] S E Derenzo. Mathemat ica l removal of positron range blurr ing i n a h igh reso-lu t ion tomography. IEEE Transactions on Nuclear Science, NS-33(l) :565-569, 1986. [34] S E Derenzo. Monte carlo calculations of the detection efHciecy of arrays of N a l ( T l ) , B G O , C s F , Ge, and plastic detectors for 511 kev photons. IEEE Trans-actions on Nuclear Science, NS-28(1):131-136, 1981. [35] S E Derenzo. Recent developments i n positron emission tomograph (pet) instru-mentation. In SPIE vol.671: Physics and Engineering of Computerized Multidi-mensional Imaging and Procesing, pages 232-243, 1986. [36] S E Derenzo, R H Huesman, J L Cahoon, A B Geyer, W W Moses, D C Uber , T Vule t ich , and T F Budinger . A positron tomograph w i th 600 B G O crystals and 2.6 m m resolution. To appear i n : IEEE Transactions on Nuclear Science, February 1988. [37] D E Dudgeon and R M Mersereau. Multidimensional Digital Signal Processing. Prent ice-Hal l , 1984. [38] L Er iksson and Z H Cho. A simple absorption correction i n positron (annihilat ion gamma coincidence detection) transverse axial tomography. Physics in Medicine and Biology, 21(3):429-433, 1976. [39] J S Flower and A P Wolf. Posi t ron emitter-labelled compounds: Priorit ies and problems. In M E Phelps, J C Mazz io t ta , and H R Schelbert, editors, Positron Emission Tomography and Autoradiography, pages 391-492, Raven Press, 1986. [40] G Harauz and F P Ottensmeyer. Interpolation i n computing forward projections i n direct three-dimensional reconstruction. Physics in Medicine and Biology, 28(12):1419-1429, 1983. [41] M R Hayden, W R W M a r t i n , A J Stoessl, C Cla rk , S Hollenberg, M J A d a m , W A m m a n n , R Harrop, J Rogers, T R u t h , C Sayre, and B D Pate. Posi t ron emission tomography i n the early diagnosis of Huntington's disease. Neurology, 36(7):888-894, 1986. [42] G T Herman. Image Reconstruction from Projections. Academic Press, 1980. [43] E J Hoffman, M Dah lbom, A R R i c c i , and I N Weinberg. Examina t ion of the role of detection systems i n quanti tat ion and image quali ty i n P E T . IEEE Transactions on Nuclear Science, NS-33(l) :420-424, 1986. [44] E J Hoffman and M E Phelps. Posi t ron emission tomography: Principles and quanti tat ion. In M E Phelps, J C Mazz io t t a , and H R Schelbert, editors, Positron Emission Tomography and Autoradiography, pages 237-286, Raven Press, 1986. [45] B K P Horn . Densi ty reconstruction using arbitrary ray-sampling schemes. Pro-. ceedings of the IEEE, 66(5):551-562, 1978. [46] S-C Huang, E J Hoffman, M E Phelps, and D E K u h l . Quant i ta t ion i n positron emission tomography:3. Effect of sampling. Journal of Computer Assisted To-mography, 4(6):819-826, 1980. 87 ) [47] S-C Huang and M E Phelps. Principles of tracer kinetic modeling i n positron emission tomography and autoradiography. In M E Phelps, J C Mazz io t t a , and H R Schelbert, editors, Positron Emission Tomography and Autoradiography, pages 287-390, Raven Press, 1986. [48] R H Huesman, S E Derenzo, J L Cahoon, A B Geyer, W W Moses, D C Uber , T Vule t ich , and T F Budinger. Orb i t ing transmission source for positron tomog-raphy. To appear in : IEEE Transactions on Nuclear Science, February 1988. [49] A Imiya and H Ogawa. Direct reconstruction of three dimensional image from its line integrals. In SPIE vol.671: Physics and Engineering of Computerized Multidimensional Imaging and Procesing, pages 42-49, 1986. [50] L Kaufman. Implementing and accelerating the E M algori thm for posi tron emis-sion tomography. IEEE Transactions on Medical Imaging, MI-6(1):37-51, 1987. [51] P E K i n a h a n , J G Rogers, R Harrop, and R R Johnson. Three-dimensional image reconstruction i n object space. To appear in : IEEE Transactions on Nuclear Science, February 1988. [52] R M Lewi t t . Reconstruction algorithms: transform methods. Proceedings of the IEEE, 71(3):390-408, 1983. [53] R M Lewit t and G Muehllehner. Accelerated iterative reconstruction for positron emission tomography based on the E M algori thm for max imum likelyhood esti-mat ion. IEEE Transactions on Medical Imaging, MI-5( l ) :16-22 , 1989. [54] M J L igh th i l l . Introduction to Fourier Analysis and Generalized Functions. Cambridge Universi ty Press, 1958. 88 [55] C B L i m , A Cheng, D P B o y d , and R S Hatner. A 3-d iterative reconstruction method for stationary planar positron cameras. IEEE Transactions on Nuclear Science, NS-25(1):196-201, 1978. [56] J Llacer , A Veklerov, and E Veklerov. Towards a pract ical implementat ion of the M L E algor i thm for positron emission tomography. IEEE Transactions on Nuclear Science, NS-33(l) :468-477, 1986. [57] A Macovsk i . Physica l problems of computerized tomography. Proceedings of the IEEE, 71(3):373-378, 1983. [58] D A Mankoff and G Muehllehner. Performance of positron imaging systems as a function of energy threshold and shielding depth. IEEE Transactions on Medical Imaging, MI-3( l ) :18-24 , 1984. [59] D A Mankoff, G Muehllehner, and J S K a r p . The effect of detector performance on high countrate P E T imaging wi th a tomograph based on positron-sensitive detectors. To appear in : IEEE Transactions on Nuclear Science, February 1988. [60] B M Mazoyer , M S Roos, and R H Huesman. Dead t ime correction and counting statistics for positron tomograph. Physics in Medicine and Biology, 30(5):385-399, 1985. [61] J C Mazz io t t a and M E Phelps. Posi t ron emission tomography sudies of the bra in . In M E Phelps, J C Mazz io t ta , and H R Schelbert, editors, Positron Emission Tomography and Autoradiography, pages 493-579, Raven Press, 1986. [62] P L McGeer , H K a m o , R Harrop, D K B L i , H Tuokko, E G McGeer , M J A d a m , W A m m a n n , B L Beattie, D B Calne, W R W M a r t i n , B D Pate, J G Rogers, T J R u t h , C I Sayre, and A J Stoessl. Posi t ron emission tomography i n patients wi th 89 cl inical ly diagnosed Alzheimer 's disease. Canadian Medical Association Journal, 134:597-607, 1986. [63] R M Mersereau. Direct fourier transform techniques i n 3-d image reconstruction. Computers in Medicine and Biology, 6:247-258, 1976. [64] G Muehllehner and J S K a r p . Posi t ron emission tomography imaging - technical considerations. Seminars in Nuclear Medicine, X V I ( l ) : 3 5 - 5 0 , 1980. [65] G Muehllehner, J S K a r p , D A Mankoff, D Beerbohm, and C E Ordonez. Design and performance of a new positron emission tomograph. To appear i n : IEEE Transactions on Nuclear Science, February 1988. [66] O Nalc ioglu and Z H Cho. Reconstruction of 3-d objects from cone bean pro-jections. Proceedings of the IEEE, 66(11):1584-1585, 1978. [67] S S Orlov. Theory of three dimensional reconstruction. I conditions for a com-plete set of projections. Soviet Physics Crystallography, 20(2):312-314, 1976. [68] S S Orlov. Theory of three dimensional reconstruction. II the recovery operator. Soviet Physics Crystallography, 20(4):429-433, 1976. [69] L M Pecora. 3d tomographic reconstruction from 2d data using spherical har-monics. IEEE Transactions on Nuclear Science, NS-34(2):642-650, 1987. [70] N J Pelc . A Generalized Filtered Backprojection Algorithm for Three Dimen-sional Reconstruction. P h D thesis, Harvard School of Pub l i c Heal th , 1979. [71] N J Pelc and D A Chesler. Ut i l i za t ion of cross-plane rays for three-dimensional reconstruction by filtered back-projection. Journal of Computer Assisted To-mography, 3(3):385-395, 1979. 90 [72] F C Peyr in . The generalized back projection theorm for cone beam reconstruc-t ion . IEEE Transactions on Nuclear Science, NS-32(4):1512-1519, 1985. [73] M E Phelps, E J Hoffman, S-C Huang, and D E K u h l . Posi t ron emission tomog-raphy: present and future design alternatives. In IAEA-SM-247/92, 1980. [74] M E Phelps and J C Maz io t t a . Posi t ron emission tomography: human bra in function and biochemistry. Science, 228(4701):799-809, 1985. [75] W H Press, B P Flannery, S A Teukolsky, and W T Vetter l ing. Numerical Recipes. Cambridge Universi ty Press, 1986. [76] J B R a and Z H Cho . Generalized true three-dimensional reconstruction algo-r i t hm. Proceedings of the IEEE, 69(5):668-670, 1981. [77] J B R a , C B L i m , Z H Cho, S K H i l a l , and J Core l l . A true three-dimensional reconstruction algori thm for the spherical positron emission tomograph. Physics in Medicine and Biology, 27(l) :37-50, 1982. [78] J Radon . O n the determination of functions from their integral values along certain manifolds. IEEE Transactions on Medical Imaging, MI-5(5):170-176, 1986. Th i s is a translation by P C Parks of the original German text by Johann R a d o n published i n Berichte der Sachsischen Akadamie der Wissenschaft, V o l . 69, pp. 262-277, 1917. (Session of A p r i l 30,1917). [79] G N Ramachandran and A V Lakshminarayanan. Three-dimensional reconstruc-t ion from radiographs and electron micrographs: Appl ica t ions of convolutions instead of fourier transforms. Proceedings of the National Academy of Science of USA, 68(9):2236-2240, 1971. 91 [80] J G Rogers. Testing an improved scinti l lat ion camera for P E T and S P E C T . IEEE Transactions on Nuclear Science, NS-33(l) :519-522, 1986. [81] J G Rogers, R Harrop, and P E K inahan . The theory of three-dimensional image reconstruction for P E T . Physics in Medicine and Biology, MI-6(3):239-243, September 1987. [82] J G Rogers, R Harrop, P E Kinahan , N A W i l k i n s o n , P W Doherty, and D P Saylor. Conceptual design of a whole body P E T machine. To appear i n : IEEE Transactions on Nuclear Science, February 1988. [83] J G Rogers, D P Saylor, R Harrop, X G Yao, C V M Leitao, and B D Pate. Design of an efficient position sensitive gamma ray detector for nuclear medicine. Physics in Medicine and Biology, 31(10):1061-1090, 1986. [84] A Rosenfeld and A C K a k . Digital Picture Processing, chapter 8. Volume I, Academic Press, second edition, 1982. [85] T R u t h . Private communication. [86] B Schorr, D Townsend, and R Clack. A general method for three-dimensional filter computat ion. Physics in Medicine and Biology, 28(9):1009-1019, 1983. [87] L A Shepp and B F Logan. The fourier reconstruction of a head section. IEEE Transactions on Nuclear Science, NS-21:21-43, 1974. [88] L A Shepp and Y V a r d i . M a x i m u m l ikelyhood reconstruction for emission to-mography. IEEE Transactions on Medical Imaging, MI-1(2):113-122, 1982. [89] L A Shepp, Y V a r d i , J B R a , S K H i l a l , and Z H Cho . M a x i m u m likelyhood P E T wi th real data. IEEE Transactions on Nuclear Science, NS-31(2):910-912, 1984. 92 [90] B D Smi th . Image reconstruction from cone-beam projections: necessary and sufficient conditions and reconstruction methods. IEEE Transactions on Medical Imaging, MI-4( l ) :14-25 , 1985. [91] L Sokoloff. Cerebral circulat ion, energy metabolism, and protein synthesis: gen-eral characteristics and priciples of measurement. In M E Phelps, J C Mazz io t t a , and H R Schelbert, editors, Positron Emission Tomography and Autoradiogra-phy, pages 1-71, Raven Press, 1986. [92] L Sokoloff. Local iza t ion of functional act ivi ty i n the central nervous system by measurement of glucose ut i l iza t ion w i th radioactive deoxyglucose. Journal of Cerebral Blood Flow and Metabolism, 1:7-36, 1981. [93] J A Sorenson and M E Phelps. Physics in Nuclear Medicine, chapter 20. Grune and Straton, second edition, 1987. [94] C W Stearns, C A Burnham, D A Chesler, and G L Brownel l . Simulat ion studies for cyl indr ica l posi tron tomography. To appear i n : IEEE Transactions on Nuclear Science, February 1988. [95] C W Stearns, D A Chesler, and G L Brownel l . Three dimensional image re-construction i n the fourier domain. IEEE Transactions on Nuclear Science, NS-34(l) :374-378, 1987. [96] C W Stearns, D A Chesler, J E Ki r sch , and G L Brownel l . Quanti tat ive imaging w i t h the M G H analog r ing positron tomograph. To appear i n : IEEE Transac-tions on Nuclear Science, February 1988. [97] K C T a m , G C h u , V Perez-Mendez, and C B L i m . Three-dimensional recon-struction i n planar positron cameras using fourier deconvolution of generalized tomograms. IEEE Transactions on Nuclear Science, NS-25(1):152-159, 1978. 93 ) [98] E Tanaka. A fast reconstruction algori thm for stationary positron emission to-mography based on a modified E M algorithm. IEEE Transactions on Medical Imaging, MI-6(2):98-105, 1987. [99] E Tanaka. Recent progress on single photon and positron emission tomography - from detectors to algorithms. IEEE Transactions on Nuclear Science, N S -34(l):313-320, 1987. [100] M M Ter-Pogossian, M E Raichle, and B E Sobel. Positron-emission tomography. Scientific American, 243(4):170+, October 1980. [101] C J Thompson. The effect of col l imat ion on scatter fraction i n multi-slice P E T . To appear i n : IEEE Transactions on Nuclear Science, February 1988. [102] C J Thompson, A Dagher, D N Lunney, S C Strother, and A C Evans. A technique to reject scattered radiat ion i n P E T transmission scans. In SPIE vol.671: Physics and Engineering of Computerized Multidimensional Imaging and Procesing, 1986. [103] D M Ti t te r ington. O n the iterative image space reconstruction algori thm for E C T . IEEE Transactions on Medical Imaging, MI-6( l ) :52-56 , 1987. [104] D Townsend, R Clack, R Maganin i , P Frey, and A Donath . Image reconstruction for a rotat ing positron tomograph. IEEE Transactions on Nuclear Science, N S -30(l):594-600, 1983. [105] B K Vainshtein and S S Orlov. Theory of the recovery of functions from their projections. Soviet Physics - Crystallography, 17(2):213-216, 1972. [106] A P Wolf. Special characteristics and potential for radiopharmaceuticals for posi tron emission tomograph. Seminars in Nuclear Medicine, XVI (1 ) :2 -12 , 94 1981. 95 Appendix A PSF Derivation T h i s appendix contains the derivation of the Point Spread Funt ion ( P S F ) for a rota-t ional ly symmetric solid angle detector, such as a rotat ing parallel-plate P E T scanner. Th i s derivation follows that of Orlov [68]. A . l Introduction F r o m linear systems theory [7], for the system shown below, h(x) i f / ( x ) 1 is the input to a linear and invariant system that has a P S F given by h(x), then the output, g{x), is given by, g(x) = f(x) * h(x). (A.l) 1 In this appendix vectors are denoted by an arrow, and unit vectors are denoted by a caret. 96 where * represents convolution of the appropriate dimensionality. For the case of interest, f(x) is the three-dimensional source density distr ibution which is being imaged, g(x) is the backprojection of the measured projections, and h(x) is the P S F of the solid angle detector which we want to calculate. The backprojection, g(x), is formed by first taking the two-dimensional projections of f(x) at different angles and then projecting back the same values (constant along each ray) through a three-dimensional volume. A t each point i n the volume the values of a l l the projection rays passing through each point are summed up, which results i n the backprojected image g(x). O u r purpose is to find h(x) such that equation A . l describes the same operation mathematically. A.2 Detector Geometry To model the geometry of a rotating-plate detector, i t is convenient to use a unit sphere w i th the top and bo t tom end caps removed as shown i n figure 3.6 The object to be imaged, / ( i f ) , is located entirely wi th in the detector surface, where ip is the cutoff angle of the detectors and f is a unit vector perpendicular to the projection or f-plane. Thus the f-plane is determined by the two spherical coordinate angles (0, <f) and the truncated sphere of figure 3.6 describes a l l possible orientations of the unit vector f . The two-dimensional projections, p(l; r), are the integrals of f(x) along the rays parallel to f so that for each plane f can be considered to be a parameter. The projection values for a given f are thus given by (A.2) 97 A.3 Derivation of PSF To construct the backprojection, g(x), the projection data are integrated over a l l possible values of f for TT/2 — if> < 9 < n. r2v fir/2 q(x) = / / p(I: r) sin 9d9dd> *v J J<i>=o Je=*/2-4, v ' N o w since p(\\ f) — p(l; — f), /•2TT rir/2+Tp g(x) = l/2 / p( I; f ) sin 9d9d<f> 7«*=o h-Tr/2-d> /•  i-k/ 4> /0=O Je=ir/2-iJ ~ Recal l ing the definition of the function rect(x), rect(x) = < 1 i f | x | < l / 2 0 i f |x | > 1/2 so that rect 9-TT/2 2tp 1 i f |^| < -0 dr T T / 2 0 i f \9\ > ip ± T T / 2 equation A . 4 is now rewriten as and using equation A . 2 , r2ir riz g(x) = 1/2 / / rect J<h=Q J9=0 '6 - T T / 2 \ y°° or <7* = / / / J rect / 1 ' 2 J<t>=0 J9=0 Jt=-oo t1 \ 2lp (A.3) (A.4) 0(x) = 1/2 T rect ( ^ ) p{T; r) sin9d8d<f> (A.5) f(x + tr)dtsin9d9d<j> (A.6) t2dt sin 9d9d<f> (A.7) where the change i n l imits of integration is justified since i n equation A . 6 the integra-t ion is over a l l space twice. 98 N o w by considering (t, 9, (f>) as a spherical coordinate system, it can be converted to a rectangular coordinate system, 5? — (x',y',z'), by, ( t , M ) - ( * ' , y ' , 0 tr —> x1 t2 -• |^|2 6 -• 6' t2dt sin 9 d9d<f> - • dx'dy'dz' where 9' — arccos j , z I . Equat ion A . 3 now becomes \yjx'2 + y'2) m - L L L * ^ ^ ) " * " (A-8) Recognizing the right hand side of equation A . 8 as the three-dimensional convolu-t ion of f(x) w i t h the function 1 (9 - T T / 2 ' —rect 2t/> J where r is understood to be the spherical radia l coordinate, that is r 2 = \x\2. Now rewri t ing equation A . 8 using convolution notation, g(x) = f(x) * I r e c t (^^) • (A.9) F ina l ly , by comparing equations A . l and A . 9 it can be seen that n = n(r,&, <p) — —rect r2 \ 2ip Thi s is the P S F of the solid angle detector for a given cutoff angle ip that is described on page 33. 99 Appendix B Derivation of L(f x xr) The definition of the L{9 x x T ) function supplied by Orlov [68] is , L(? x x T ) = - j - j dGT, 8{r' • (x x r ) ) . ( B . l ) G The verbal definition of L(r x x T ) by Orlov is quoted i n section 3.1.4. In the devel-opment of the three dimensional image reconstruction algori thm an explicit form for L{T X X t ) is needed for the case where the region G is the surface of a truncated hemisphere, as shown i n figure B . l Since S(T'• (X x ?)) = 0 unless ?', x and f are coplanar, the integral i n equation B . l sweeps out the length of the port ion of the semi-great circle on the plane that is i n the region G . Th i s length is shown as the solid line part of the semi-great circle. The dashed line part cannot contribute to the value of L(r x x T ) since it lies outside the region G . A l so note i n figure B . l that x = x / | x | and the parameter ip is used to indicate where the unit hemisphere is truncated along the z-axis. There is an addit ional constraint added that x • f = 0 because the vector of intrest is xT = x — (x • f)?, which is the projection of the vector x onto the plane perpendicular to f . The statement of the problem is now, given x , f, ift, and the constraint x • f = 0, find L{r x x T ) which is the total length of the solid line shown i n figure B . l . Th i s can be simplified to finding the length of the dashed line, £ , between points P i and P2, 100 XL z Figure B . l : Geometry of the L(r x xT) function 101 since L(f x xT) = TT — £ The curved line segments lie on the plane where f ' • ( i x f ) = 0. Th i s plane has a normal unit vector given by (a, b, c) = x x r, so the plane is called the (a, b, c) plane. The points P i and P 2 are the intersection of the plane z = sin ip and the semi-great circle ly ing on the (a, 6, c) plane. Looking now at the (a,b,c) plane i n figure B .2 , the straight line segment from P i to P 2 has a length of d and since the circle has unit radius the angle between the line segments OP\ and OP2 is £ . Th i s means that £ = 2arcsin((f/2) and i f P i = (x\,y\,z\) and P 2 = ( x 2 , y 2 , z2), d = y/(xx - x 2 ) 2 + (yx - y2)2 + (zx - z2)2 but Zi = z2 = smV>, so I ( f x xT) = TV - £ = T T - 2arcs in Q\/(xi - x2)2 + (yi - y 2 ) 2 ) (B.2) and the problem now reduces to finding the points ( x i , y i ) and ( x 2 , y 2 ) . The points Pi and P 2 satisfy three conditions: they lie on the unit radius semi-great circle, they are on the z = sin ip plane, and they are also on the (a, 6, c) plane. L i s t ing these constraints i n the above order, x2 + y 2 + z2 = 1 z — simp (B.3) ax + by + cz = 0 The last equation above comes from the definition of a plane passing through the origin that has a unit normal of (a, b, c) [1]. Now i n order to solve equations B.3 for (x, y, z) it is necessary to first find the values a, b, and c. Us ing standard spherical coordinates as shown i n figure 3.6, x = ( s i n ^ c o s ^ s i n & r s i n ^ c o s f l ; , ; ) f = (sin# T cos <f>T, sin# T sin</>T, cos 8T) 102 o 1 Figure B . 2 : Intersection of the (a, 6, c) plane and the region G . 103 Another usefull constraint for a plane through the origin and coplanar w i t h both x and f is [1], x y z sin 8X cos (j>x sin 9X sin <f>x cos 8X = 0 sin 9T cos <f>T sin 6T sin >^T cos 8T Thi s can be rewrit ten as, ax + by + cz = 0 where a = b = c = sin 8X sin cos 8T — sin 0 T sin <f>T cos 0 r sin 0 T cos <f>T cos 0^ — sin8 X cos </>x cos 8T (B-4) sin 0^ cos <j>x sin 0 T sin <^ T — sin 8X sin ^ sin 8T cos ^>T = s i n 0 r s i n 0 T sin(^>x — <^ )T) F r o m equations B.3 and B.4 it is now possible to solve for P\ and P2 i n terms of the coordinates (8X,</>X) of x and (8T,<f>T) of f . In the solution it useful to note that since (a, 6, c) = x x f and £ • f = 0, a2 + b2 + c2 = 1. Us ing this, there are two solutions for (x, y, z) and, \\l{xx - x2)2 + (y! - y 2 ) 2 = C = ^ 1 -Th i s means that, s i n 2 V' cos 2 0 X + cos 2 4>2 I L(T x x T ) = L(9X; 9T) = 7r — 2arcs in and since, i f 0 < a: < 1 [1], V 1 -s i n 2 ip \ cos 2 8X + cos 2 4>x ir — 2arcs in (Vl — z 2 ) == 2arcsin(a;) 104 the explicit form of L(f x x T ) i s L(9X;8T) = 2arcs in f , 2 • 2 , ) • \ VCOS J 0 X + COS^  / In the above equation however, the arcsin function becomes undefined i f s i n 2 tp > cos 2 6X + cos 2 6r. In figure B . l this is equivalent to the semi-great circle ly ing entirely w i th in the region G so that the points P x and P2 do not exist. B u t i n this case C = 0 and L(9X; 9T) = TT, so 7r, i f s i n 2 ip > cos 2 8X + cos 2 9T L(9x;9T) = l ( • I \ (B.5) J I 2 arcsin , s i n ^ , i f s i n 2 tp < cos 2 9X + cos 2 9T J \ y cos 2 9X + cos 2 9T J Note that for sin2t/> = cos 2 9X + cos 2 9T, the points Pi and P2 are coincident and L(9X;9T) = 7r, thus L(9X;9T) is continuous but has a discontinuity i n slope. It is necessary i n chapter 3 to have the function L(9X;9T) expressed i n the coor-dinates of the projection plane as shown i n figure 3.6. Us ing spherical trigonome-try similar to that of figure 5.1, the angle 9X is related to 9\ and 8T by the relation cos 9X = — s i n 0 j s i n 0 T . Th i s means that cos 2 9X + cos 2 9T = 1 — cos 2 9, s i n 2 9T and T, i f | cosi/>| < I COS0j| • I cos# T | L(9,;9T)={ ( • , \ I 2arcs in — 2 i 2 J 1 i f | cosV>| > | cos0, | • | c o s 0 r | \ V 1 - c o s 2 0 , s i n 2 0 T / Th i s is the result quoted i n section 3.1.6. 105 Appendix C Derivatives of L($i',6T) Thi s appendix shows the calculation of the derivatives of L(9i;9T) w i th respect to F r o m appendix B , L(8i; 9r) is given by, L(9,;9T) = < 7T, i f | cos ip I < I cos Oi I • I cos 9TI 2 arcsin j , s i n ^ \ ^ | C Q g ^, | > | c o s 9i\ • \ cos 9T\ cos281 s i n 2 0 T / In the case where L(8i; 9r) — 7r, a l l derivatives are zero and it is only necessary to to find the derivatives where L(9f,9T) ^ ir. Note that this means there is a discontinuity of the derivatives at | cost/>| = | cos#/| • | cos# T | . For consistency w i th chapter 3, a new function 9T) is defined as, L(9r,9T) = 2arcsin | S m ^ ) . \yfl- cos 2 8i s i n 2 9Tj C . l First Derivative To start wi th ; 6 r ) = S i n ^ • A ( 1 _ cos 2 9, s i n 2 0 T ) " 1 / 2 ^ / , _ sin2V> dei V cos 2 9i s i n 2 0 T 106 and so F ina l ly , d 2 2 1 / 2 s ing; c o s 0 / s i n 2 0 T —-(1 — cos 61 sin 6T) ' = - ——. 9 „ .... der ' 2 ( 1 - c o s 2 0, s i n 2 0 T ) 3 / 2 1^ s in if> cos 2 6i s i n 2 0 T sin 0/ cos 0/ s i n 2 0 T "2(1 - cos 2 0/ s i n 2 0 T ) 3 / 2 sin 20/ s i n 2 0 T sin ip —L(9t; 9T) = , d(*i (1 - cos 2 0/ s i n 2 0T) V cos 2 V> - cos 2 0/ s i n 2 0 T C.2 The Second Derivative Using the above result and continuing on, WL{0i'%Or) = lk d9, [A(9l]9T)B(9l]9T)\ or 50/ where by definition, and Q2 €S O —1(6,; 9T) = QfAtfr, 9r)B(9,; 0 T) + A(6,; 0 T) ^ ( 8 , ; 9T) ddi A ( 0 , ; 0 T ) = sin 20/ s i n 2 0 T sin i/> "(1 - cos 2 91 s i n 2 0 T ) 5(0/ ; 0 T) = (cos 2 </> - cos 2 0/ s i n 2 0 T ) " 1 / 2 -The derivatives of A(0/ ; 0 T) and B(8\; 0 T) are given by, — A(a a\— 2 c o s 2 ^ / ( 1 - cos 2 0 / s in 2 0 T ) s in 2 0 T s inz /> - s i n 2 20 / s in 4 0 T sin ^ — A ( 0 / ; 0 T ) _ * ( l - c o s 2 0 / s i n 2 0 T ) 2 and JLB(9-9 \ = sin20, s i n 2 0 T d8t y h t ) 2(cos 2 tp - cos 2 0/ s i n 2 0 T ) 3 / 2 ' 107 The expression for the second derivative of L{9[\ #T)is too long to write out properly, but it is expl ici t ly expressed by the above five equations. To usehe results of this section i n equation 3.27, it is only necessary to note that, d ( i \ i a d9, \m-9T)j L(9r,6Tf d6, and L(9,-9T) 108 Appendix D Filter Program Listing 109 program g e n e r a t e _ f l i t e r C • C GENERATE_F1LTER C Generates a set of 2D reconstruction f i l t e r s . C The f i l t e r s are n x n and there are t of them, so each value of t C represents a d i f f e r e n t 2D f i l t e r . C 1. The f i l t e r parameters ( n , t , p s i , p i x e l _ s i z e ) are read by get_parameters() C from the f i l e par f i l e . C 2. An array f o r s t o r i n g the f i l t e r s i s created and passed with the parameters C to c a l c u l a t e _ f i l t e r ( ) , which does a l l the actual c a l c u l a t i o n s . C 3. A unique filename f o r the f i l t e r set i s created by make_filename() C and passed to o u t p u t _ f i l t e r ( ) , which copies the unformatted 4-byte C f i l t e r values to the f i l e . C 4. An ASCII formated output of the f i l t e r (for humans) along with the C f i l t e r parameters and filenames i s produced by hardcopy(). C 5. F i n a l l y the array storage space i s released. C C History: 2nd. version written l9-feb-1988 Paul Kinahan c c i m p l i c i t none integer n,t,filter_name_len,max_name_len,par_name_len,wordsize, + l i b $ g e t _ v m , l i b $ f r e e _ v m , s t a t u s , f i l t e r _ a d d r , f i l t e r _ s i z e parameter(max_name_len=50,wordsize=4) r e a l p s i , p i x e l _ s i z e character*(max_name_len) f i l t e r _ f i l e , p a r _ f i l e c a l l get_parameters(n,t,psi,pixel_size,par_flle,par_name_len) f i l t e r _ s i z e = (n**2)*t*wordsize ~ s t a t u s = l i b $ g e t _ v m ( f i l t e r _ s i z e , f i l t e r _ a d d r ) i f ( . n o t . s t a t u s ) c a l l l i b $ s i g n a l ( % v a l ( s t a t u s ) ) c a l l c a l c u l a t e _ f i l t e r ( n , t , p s i , p i x e l _ s i z e , % v a l ( f i l t e r _ a d d r ) ) c a l l m a k e _ f i l e n a m e ( n , t , p s i , p i x e l _ s i z e , f i l t e r f i l e , f i l t e r _ n a m e _ l e n ) c a l l o u t p u t _ f i l t e r ( n , t , p s i , p i x e l _ s i z e , % v a l ( f I l t e r _ a d d r ) , + f i l t e r _ f i l e , f i l t e r _ n a m e _ l e n ) c a l l h a r d c o p y ( n , t , p s i , p i x e l _ s i z e , % v a l ( f i l t e r addr), + p a r _ f i l e , p a r _ n a m e _ l e n , f i l t e r _ f T i e , f i l t e r _ n a m e _ l e n ) s t a t u s = l i b $ f r e e _ v m ( f i l t e r _ s i z e , f i l t e r _ a d d r ) i f ( . n o t . s t a t u s ) c a l l l i b $ s i g n a l ( % v a l ( s t a t u s ) ) end 110 subroutine c a l c u l a t e _ f i l t e r ( n > t , p s i , p i x e l _ s i z e , f i l t e r ) C C CALCULATE_FILTER C Does the actual c a l c u l a t i o n of the set of 2D reconstruction f i l t e r s C as explained i n my t h e s i s . Only the f i r s t quadrant i s c a l c u l a t e d C since the f i l t e r s are symmetrical about the x- and y-axes (and about C theta_tau = 90deg). The f i l t e r values are c a l c u l a t e d f o r each p i x e l C by c a l c u l a t e _ p i x e l ( ) , which also uses the weighting arrays f o r numerical C i n t e g r a t i o n . I f the p i x e l i s close to the o r i g i n or the l i n e of C s i n g u l a r i t i e s , t h e f i n e - g r i d d i n g array i s passed, otherwise the course-C gridding array i s passed. C C History: 2nd version written 4-mar-1988 PK C C i m p l i c i t none l o g i c a l use_fine_gridding integer n , t , t i , i , j , f i n e divs,course divs,numpix,totalpix, + timer_addr,lib$Init_timer,lib$show_timer,lib$free_timer parameter(fine_divs=50,course_divs=8) r e a l t h e t a _ t a u , p s i , p i x e l _ s i z e , t h e t a _ b a r , x i , y i , p i , + fine_weights(0:fine_divs),course_weights(0:course_divs), + f i l t e r ( 0 : n - l , 0 : n - l , 0 : t - l ) , c a l c u l a t e _ p i x e l , q u a d l , p o s _ a x e s parameter(pi=3.1415927) ! Set up ID arrays of weights for numerical i n t e g r a t i o n (Simpsons) ! weights f o r f i n e gridding of p i x e l fine_weights(0) =1.0 ! endpoints fine_weights(fine_divs) =1.0 do i=l,fine_divs-3,2 fine_weights(i) =4.0 ! even points fine_weights(i+l) =2.0 I odd points enddo fine_we i g h t s ( f i n e _ d i v s - l ) =4.0 ! penultimate point I weights for course gridding of p i x e l course_weights(0) =1.0 ! endpoints course_weights(course_divs) =1.0 do i=l,course_divs-3,2 course_weights(i) =4.0 ! even points course_weights(i+l) =2.0 ! odd points enddo course_weights(course_divs-l) =4.0 ! penultimate point t o t a l p i x = (n**2)*t type ' (/,a,i5,a,/) *, ' There are a t o t a l o f , + t o t a l p i x , ' p i x e l s to c a l c u l a t e ' c a l l l i b $ init_timer(timer_addr) 1 to s t a r t c o l l e c t i n g program s t a t i s t i c s do t i = 0 , t - l ! Calculate f i l t e r values for each theta_tau. theta_tau = pi/2.0 - p s i * f l o a t ( t i ) / f l o a t ( t - l ) theta_bar=acos(cos(psi)/sln(theta_tau)) ! c r i t i c a l angle do j=0,n-l ! Calculate f i l t e r values f o r 1st. Quadrant do i=0,n-1 numpix = i + 1 + j*n + ti*n**2 type '(a,i5)','+ Doing p i x e l number: ',numpix x i = f l o a t ( i ) * p i x e l _ s i z e y i = f l o a t ( j ) * p i x e l _ s i z e use_fine_gridding = ( ( i . l e . 2 ) . a n d . ( j . l e . 2 ) ) ! i f near o r i g i n or 111 + . o r . ( ( y i * c o s ( t h e t a _ b a r ) - x i * s i n ( t h e t a _ b a r ) ) . l e . p i x e l _ s i z e ) I c r i t i c a l i f ( u s e _ f i n e _ g r i d d i n g ) then 1 angle f i l t e r ( i , j , t i ) = c a l c u l a t e _ p i x e l ( x i , y i , t h e t a ^ t a u , p s i , + t h e t a _ b a r , p i x e l _ s i z e , f i n e d T v s , f i n e weights) e l s e ! use course gridding to speed up computation time f i l t e r ( i , j , t i ) = c a l c u l a t e _ p i x e l ( x i , y i , t h e t a _ t a u , p s i , + theta_bar,pixel_size,course_divs,course_weights) endif enddo enddo 1 set centre of f i l t e r t o negative of the r e s t of f i l t e r , c a l u l a t i n g axes ! seperately guadl =0.0 pos axes = 0.0 do 5=l»n~l do i = l , n - l quadl = quadl + f i l t e r ( i , j , t i ) enddo enddo do i = l , n - l pos_axes = pos_axes + f i l t e r ( 0 , i , t i ) + f i l t e r ( i , 0 , t i ) enddo f i l t e r ( 0 , 0 , t i ) = -(4*quadl + 2*pos_axes) c a l l window(filter,n,ti) ! to low-pass f i l t e r the reconstruction f i l t e r enddo ! end of loop f o r each value of theta_tau type * (•a)',,+ Finished • type ' ( / , a ) ' , ' F i l t e r generation program s t a t i s t i c s ' type ' ( a , / ) ' , ' • c a l l lib$show_timer(timer_addr) c a l l lib$free_timer(timer_addr) type »(a)',' ' end 112 r e a l function c a l c u l a t e _ p i x e l ( x i , y i , t h e t a _ t a u , p s i , t h e t a _ b a r , + pixel_size,divs,weights) C C CALCULATE_PIXEL C Calculates the value of a s i n g l e p i x e l of width p i x e l _ s i z e centred C at ( x i , y i ) as explained i n my t h e s i s . C The c a l c u l a t i o n i s broken into 3 parts depending on where the p i x e l C i s i n r e l a t i o n to the c r i t i c a l angle theta_bar. The 4 corners of the C p i x e l are given by x l , x 2 , y l , and y2. C C History: 2nd version written 19-feb-1988 PK C-i m p l i c i t none l o g i c a l i n t e r s e c t i o n integer divs r e a l x i , y i , t h e t a _ t a u , p s i , t h e t a _ b a r , p i x e l _ s i z e , + xl,x2,yl,y2,theta_max,thetajmin,weights(0:divs), + integral_l,integral_2,integral_3,term_l,term_2,tenn_3 c a l c u l a t e _ p i x e l =0.0 1 quickly return 0.0 i f we're at the o r i g i n i f ((xi.eq.0.0).and.(yi.eq.0.0)) goto 999 ! set up i n i t i a l values x l = amaxl(0.0,xi - pixel_size/2.0) ! i n case we are on y-axis x2 = x i + p i x e l _ s i z e / 2 . 0 y l = amaxl(0.0,yi - pixel_size/2.0) Jin case we are on x-axis y2 = y i + p i x e l _ s i z e / 2 . 0 theta_min = atan2(yl,x2) 1 min value of t h e t a _ l i n p i x e l theta_max = atan2(y2,xl) ! wax value of t h e t a _ l i n p i x e l ! subdivide c a l c u l a t i o n depending on where c r i t i c a l ray i s ! r e l a t i v e to p i x e l i f (theta_bar.gt.theta_max) then I i n t e r s e c t i o n =.false. term_l = i n t e g r a l _ l ( x l , x 2 , y l , y 2 ) tenn_2 = 0.0 term_3 = 0.0 elseif((theta_min.le.theta_bar).and.(theta_bar.le.theta_inax))then i n t e r s e c t i o n =.true. term_l = i n t e g r a l _ l ( x l , x 2 , y l , y 2 ) term_2 = integral_2(xl,x2,yl,y2,theta_bar) tem>_3 = integral_3(xl,x2,yl,y2,theta tau,psi,theta_bar,divs, + weights,intersection) e l s e i f (theta_min.gt.theta_bar) then i n t e r s e c t i o n =.false. term_l = i n t e g r a l _ l ( x l , x 2 , y l , y 2 ) term_2 = term_l term_3 = integral_3(xl,x2,yl,y2,theta tau,psi,theta_bar,divs, + weights,intersection) endif c a l c u l a t e _ p i x e l = term_l+term_2+term_3 999 continue 1 goto (gasp) point f o r p i x e l at o r i g i n end 113 r e a l function i n t e g r a l _ l ( x l , x 2 , y l , y 2 ) C C INTEGRAL_1 C Calculates the II term (see thesis) f o r a p i x e l with corners at x l , x 2 , y l , C and y2. C Tests r e a l values f o r zero since i t can be assigned by the c a l l i n g routine. C C History: written 19-feb-1988 PK C-i m p l i c i t none r e a l xl,x2,yl,y2,pi parameter(pi = 3.1415927) i f ( ( x l . l t . 0 . 0 ) . o r . ( y l . l t . 0 . 0 ) ) then p r i n t *,*%USR-F-KaughtyVal Negative p i x e l coordinates passed 1 + 'to routine i n t e g r a l _ l ( ) ' c a l l e x i t e l s e i f (xl.eq.0.0) then ! p i x e l i s on y-axis i n t e g r a l _ l = (1.0/(2.0*pi**3))*( sqrt(x2**2 4 y2**2)/(X2*y2) + - sqrt(x2**2 + yl**2)/(x2*yl) ) e l s e i f (yl.eq.0.0) then ! p i x e l i s on x-axis i n t e g r a l _ l = (1.0/(2.0*pi**3))*( sqrt(x2**2 + y2**2)/(x2*y2) + - s g r t ( x l * * 2 + y2**2)/(xl*y2) ) else ! p i x e l i s somewhere i n 1st quadrant i n t e g r a l _ l = (1.0/(4.0*pi**3))*( sqrt(x2**2 + y2**2)/(x2*y2) + - s q r t ( x l * * 2 + y2**2)/(xl*y2) + - sqrt(x2**2 + yl**2)/(x2*yl) + + s q r t ( x l * * 2 + y l * * 2 ) / ( x l * y l ) ) endif end 114 r e a l function integral_2(xl,x2,yl,y2,theta_bar) C C INTEGRAL_2 C Calculates the 12 term (see thesis) f o r a p i x e l with corners at x l , x 2 , y l , c and y2. C Tests r e a l values f o r zero since i t can be assigned by the c a l l i n g routine. C C History: written 19-feb-1988 PK C i m p l i c i t none r e a l xl,x2,yl,y2,theta bar,pi,integral_l,y_start,y_end,x_end, + upper_y_intersection,lower_y_intersection parameter(pi~= 3.1415927) t decide upon region of evaluation upper_y_intersection = x2*tan(theta_bar) lower_y_intersection = xl*tan(theta_bar) x_end = aminl(x2,y2/tan(theta_bar)) y _ s t a r t = amaxl(yl,lower_y_intersection) y_end = aminl(y2,upper_y_intersection) 1 evaluate part of p i x e l defined by i n t e r s e c t i o n with ray i f ( ( x l . l t . 0 . 0 ) . o r . ( y l . l t . 0 . 0 ) ) then p r i n t *,'%USR-F-NaughtyVal Negative p i x e l coordinates passed ', + • 'to routine i n t e g r a l _ 2 ( ) ' c a l l e x i t e l s e i f (xl.eq.0.0) then ! p i x e l i s on y-axis integral_2 = -(1.0/(2.0*pi**3))*cos(theta_bar) + *(1.0/y_end - 1.0/y_start) e l s e i f (yl.eq.0.0) then ! p i x e l i s on x-axis integral_2 = (1.0/(2.0*pi**3)) + *( (sqrt(y_end**2 + xl**2)/(xl*y_end)-cos(theta_bar)/y_end) + - ( s g r t ( y _ s t a r t * * 2 + x l * * 2 ) / ( x l * y _ s t a r t ) - c o s ( t h e t a _ b a r ) / y _ s t a r t ) ) e l s e ! i t must be that we are somewhere el s e i n the 1st quadrant integral_2 = (1.0/(4.0*pi**3)) + *( (sqrt(y_end**2 + xl**2)/(xl*y_end)-cos(theta_bar)/y_end) + - ( s q r t ( y _ s t a r t * * 2 + x l * * 2 ) / ( x l * y _ s t a r t ) - c o s ( t h e t a _ b a r ) / y _ s t a r t ) ) endif ! i f there i s a rectangular region l e f t , pass i t to i n t e g r a l _ l ( ) 1 to evaluate and multiply by 2 i f the p i x e l i s on an axis i f (y_end.lt.y2) + integral_2 ~ integral_2 - integral_l(xl,x2,y_end,y2) i f ((xl.eq.0.0).or.(yl.eq.0.0)) + integral_2 = 2*integral_2 end 115 r e a l function integral_3(xl,x2,yl,y2,theta_tau,psi,theta_bar, + ~ divs,weights,intersection) C C INTEGRAL_3 C Calculates the 13 term (see thesis) f o r a p i x e l with corners at x l , x 2 , y l , C and y2. C Uses a 2D Simpson's r u l e f o r numerical i n t e g r a t i o n . The 2D quadrature C weights are c a l c u l a t e d by multiplying the appropriate ID weights passed C i n the ID array weights(0:divs). The l o g i c a l v a r i a b l e ' i n t e r s e c t i o n ' i s C true i f the c r i t i c a l ray (theta_bar) passes through the p i x e l , c C History: written 19-feb-1988 PK C i m p l i c i t none l o g i c a l i n t e r s e c t i o n integer i , j , d i v s r e a l lx,ly,xl,x2,yl,y2,theta_tau,psi,theta bar,alpha_scale, + h,k,rl,r2,r,alpha,theta_l,weights(0:divs),grid_weight, + l _ f c n , l _ f c n _ d l , l _ f c n _ d 2 , p i parameter(alpha_scale=l.0,pi=3.1415927) h •= (x2 - x l ) / f l o a t ( d i v s ) ! d e l t a x g r i d width k = (y2 - y l ) / f l o a t ( d i v s ) ! d e l t a y g r i d width i f (intersection) then I c a l c u l a t e alpha r l = s q r t ( amaxl(xl,yl/tan(theta_bar))**2 + + amaxl(yl,xl*tan(theta_bar))**2 ) r2 = sqrt( aminl(x2,y2/tan(theta_bar))**2 + + aminl(y2,x2*tan(theta_bar))**2 ) alpha = alpha_scale*2*asin(amaxl(h,k)/(rl + r 2 ) ) ! to be sure to get el s e J the r i g h t distance alpha = 0.0 endif int e g r a l _ 3 =0.0 do j=0,divs ! evaluate using 2d Simpsons r u l e do i=0,divs l x = x l + ( f l o a t ( i ) ) * h ! g r i d points l y = y l + ( f l o a t ( j ) ) * k r = s q r t ( l x * * 2 + ly**2) ! polar coordinate of g r i d point t h e t a _ l = atan2(ly,lx) grid_weight = (h*k/9.0)*weights(i)*weights(j) integral_3 = integral_3 + grid_weight*-l.0/(4.0*(pi**2)*(r**3)) + *(l_fcn(theta_l,theta_tau,psi,theta_bar) + + l_fcn_d2(theta_l,theta_tau,psi,theta_bar+alpha)) enddo enddo i f (intersection) then ! add any contribution from the s i n g u l a r i t y integral_3 = integral_3 + (1.0/(4.0*pi**2))*(1.0/r2 - 1.0/rl) + *l_fcn_dl(theta_bar+alpha,theta_tau,psi,theta_bar) endif i f ((xl.eq.0.0).or.(yl.eq.0.0)) integral_3=2*integral_3 ! since on axis end 116 r e a l function l_fcn(theta_l,theta_tau,psi,theta_bar) C • C L_FCN C Calculates the Orlov L() function term as given i n t h e s i s . C C History: written 19-feb-1988 PK c i m p l i c i t none r e a l theta_l,theta_tau,psi,theta_bar,pi parameter(pi=311415927) i f (theta_l.le.theta_bar) then l f c n = p i el s e l _ f c n = 2*asin(sin(psi) + /sqrt(1.0-(cos{theta_l)*s i n ( t h e t a t a u ) ) * *2)) endif end 117 { r e a l function l_fcn_dl(theta_l,theta_tau,psi,theta_bar) c . C L_FCN_D1 C Calculates the value of the f i r s t d e r i v a t i v e of 1/L(), where L() i s C the Orlov function term as given i n t h e s i s . Statement functions are C used to abbreviate f i n a l term C C History: written 19-feb-1988 PK c i m p l i c i t none r e a l theta_l,theta_tau,psi,theta_bar, + d l f d 2 , d 3 , f l , f 2 , f 3 , a , a d l , b , b d l , l , l d l , l d 2 ! dummy arguements 1 define statement functions ! note that f o r dummy arguments: ! d l = t h e t a _ l ! d2 = theta_tau ! d3 = p s i (when used) 1 also note that a d l , b d l , l d l are the f i r s t d e r i v a t i v e s of a,b,Sl w.r.t d l I and " " ld2 i s " second " 1 " " fl ( d l , d 2 ) = 1.0 - (cos(dl)**2)*(sin(d2)**2) f2(dl,d2,d3) «= cos(d3)**2 - (cos(dl)**2)*(sin(d2)**2) f3(dl,d2) = sin(2*dl)*(sin(d2)**2) I(dl,d2,d3) = 2 * a s i n ( s i n ( d 3 ) / s q r t ( f l ( d l , d 2 ) ) ) ! the L() function again a(dl,d2,d3) = - sin ( d 3 ) * f 3 ( d l , d 2 ) / f l ( d l , d 2 ) b(dl,d2,d3) = 1.0/sqrt(f2(dl,d2,d3)) Idl(dl,d2,d3) = a(dl,d2,d3)*b(dl,d2,d3) ! 1st d e r i v a t i v e of L() 1 - End of statement function d e f i n i t i o n s i f (theta_l.le.theta_bar) then l _ f c n _ d l =0.0 else l _ f c n _ d l = - l d l ( t h e t a _ l , t h e t a _ t a u , p s i ) + / ( l ( t h e t a l , t h e t a _ t a u . p s i ) * * 2 ) endif end 118 ( r e a l function l_fcn_d2(theta_l,theta_tau,psi,theta_bar) C : C L_FCN_D2 C Calculates the value of the second d e r i v a t i v e of 1/L(), where L() i s C the Orlov function term as given i n t h e s i s , statement functions are C used t o abbreviate f i n a l term C C History: written 19-feb-1988 PK C i m p l i c i t none r e a l theta_l,theta_tau,psi,theta_bar, + d l , d 2 , d 3 , f I , f 2 , f 3 , a , a d l , b , b d l , l , l d l , l d 2 1 dummy arguements 1 define statement functions 1 note that f o r dummy arguments: I d l = t h e t a _ l ! d2 = theta_tau 1 d3 = p s i (when used) ! also note that a d l , b d l , l d l are the f i r s t d e r i v a t i v e s of a,b,&l w.r.t d l ! and " " ld2 i s " second •• 1 " " fl ( d l , d 2 ) = 1.0 - (cos(dl)**2)*(sin(d2)**2) f2(dl,d2,d3) = cos(d3)**2 - (cos(dl)**2)*(sin(d2)**2) f3(dl,d2) = sin(2*dl)*(sin(d2)**2) I(dl,d2,d3) = 2 * a s i n ( s i n ( d 3 ) / s q r t ( f l ( d l , d 2 ) ) ) 1 the L() function again a(dl,d2,d3) = - sin ( d 3 ) * f 3 ( d l , d 2 ) / f l ( d l , d 2 ) adl(dl,d2,d3) - - 2*sin(d3)*(sin(d2)**2)*cos(2*dl)/f1(dl,d2) + + s i n ( d 3 ) * ( f 3 ( d l , d 2 ) / f l ( d l , d 2 ) ) * * 2 b(dl,d2,d3) = 1.0/sqrt(f2(dl,d2,d3)) bdl(dl,d2,d3) = - f3(dl,d2)/(2*(f2(dl,d2,d3)**(1.5))) Idl(dl,d2,d3) = a(dl,d2,d3)*b(dl,d2,d3) ! 1st d e r i v a t i v e of L() ld2(dl,d2,d3) = adl(dl,d2,d3)*b(dl,d2,d3) ! 2nd d e r i v a t i v e of L() + + a(dl,d2,d3)*bdl(dl,d2,d3) I - End of statement function d e f i n i t i o n s i f (theta_l.le.theta_bar) then l_fcn_d2 = 0.0 els e l_fcn_d2 = (1.0/(l(theta_l,theta_tau,psi)**2)) + *( 2.0*ldl(theta_l,theta_tau,psi)**2 + / l ( t h e t a _ l , t h e t a _ t a u , p s i ) + - ld2(theta_l,theta_tau,psi) ) endif end 119 subroutine window(array,n,t) C C WINDOW C Low-pass f i l t e r s the 2D array indexed by 'a r r a y ( , , t ) ' by using 2D C convolution. The low-pass f i l t e r i s the object space version of a C Hamming window. The o r i g i n a l array only contains the f i r s t quadrant C values so values i n the other 3 quadrants are obtained by r e f l e c t i o n C of the i n d i c i e s 13 and j3 about the axes. C C History: written 9-oct-1987 c i m p l i c i t none integer n,i,j,i2,j2,i3,j3,nmax,tmax,t parameter(nmax=51,tmax=10) r e a l array(0:n-l,0:n-l),copy(0:nmax-l,0:nmax-l),r,w,sinc,pi parameter(pi=3.141592654) if(n.gt.nmax) then type *,' Dummy array i n routine WINDOW too small' c a l l e x i t endif do j=0,n-l ! f o r f i r s t quadrant do i=0,n-1 copy(i,j)=0 do j2 = - ( n - l ) , n - l ! f o r e n t i r e plane do i 2 = - ( n - l ) , n - l j3 = j2 i3 = i2 i f (j2.1t.O) j3 = -j2 i f (i2.1t.O) i3 = -i2 r •= s q r t ( f l o a t ( i - i 2 ) **2 + f lo a t (j - j 2) **2) w = 0.5*sinc(r) + 0.25*(sine(r-1.0) + sinc(r+1.0)) ! Hamming term copy( i . j ) = copy(i,j) + array(i3,j3)*w enddo enddo enddo enddo do j=0,n-l 1 put f i l t e r e d values i n o r i g i n a l array do i=0,n-l a r r a y ( i , j ) = copy(i,j) enddo enddo end C-r e a l function sinc(x) i m p l i c i t none r e a l x ,pi parameter (pi=3.141592654) i f (x.eg.0.0) then sine = 1.0 else sine = s i n ( p i * x ) / ( p i * x ) endif end 120 subroutine o u t p u t _ f i l t e r ( n , t , p i x e l s i z e , p s i , f i l t e r , f i l e n a m e , + filename_len) C C OUTPUT_FILTER C Writes out the unformatted 3D f i l t e r array to the f i l e 'filename'. C The f i l t e r parameters are written on the f i r s t l i n e of the f i l e , with C the f i l t e r values s t a r t i n g i n the second l i n e . C This f i l e i s the main goal of the f i l t e r generation. C C History: written 21-feb-1988 PK C i m p l i c i t none l o g i c a l e x i s t s integer n,t,i,j,k,index,filename_len,lun,1ib$get_lun,status, + max_name_len,length parameter(max_name_len=50) r e a l f i l t e r ( 0 : n - l , 0 : n - l , 0 : t - l ) , p i , p s i , p i x e l _ s i z e parameter(pi=3.1415927) character*(max_name_len) filename,newname,fullname character answer*20,char1*1,blank*1 parameter(blank=' ') 1 check with user i f f i l e name i s n ' t unique 10 inquire(file=filename(l:filename_len),exist=exists) i f (exists) then 20 type '(3a,/,a,$)',• The f i l e ', + filename(1:filename_len),' already e x i s t s ' , + ' Do you want to use a d i f f e r e n t f i l e name? [N]: ' read '(a20)',answer c h a r l *= answer (1:1) i f ((charl.ne.'Y').and.(charl.ne.' •).and.(charl.ne.'N')) then type '( a ) ' , ' Type YES, NO, or return f o r d e f a u l t [NO]' go to 20 e l s e i f ((charl.ne.'N').and.(charl.ne.' ')) then type 1 ( a , $ ) ' , ' Enter new f i l e name: ' read '(q,a)', length,newname i f (length.gt.max_name_len) then t y p e ' ( a , i 2 , a ) 1 , ' Maximum name length i s ', + max_name_len,' characters' go to 20 els e filename=newname(1:filename_len) go to 10 endif e l s e type •(3a) ', ' A new version of ', + filename(1:filename_len), + ' w i l l be created' endif endif status=LIB$GET_LUN(lun) i f ( . n o t . s t a t u s ) c a l l LIB$SlGNAL(%val(status)) open(unit=lun,file=filename(1:filename_len),status='NEW', + form=•UNFORMATTED') write(lun) 2 * n - l , 2 * t - l , p s i , p i x e l _ s i z e ! convert to f u l l values write(lun) ( ( ( f i l t e r f i , j , k ) , i = 0 , n - l ) , j = 0 , n - l ) , k = 0 , t - l ) close(lun) 121 ! Get the f u l l f i l e name inquire(file=filename(1:filename_len),name=fullname) filename_len = index(fullname,blank)-1 filename(1:filename_len) = fullname(1:filename_len) type '(2a)',' F i l t e r data i s i n f i l e «, + filename(1:filename_len) end 122 subroutine get_parameters(n,t,psi,pixel_size,filename,name_len) c C GET_PARAMETERS C Reads f i l t e r generation parameters i n from the f i l e 'filename'. C The parameters are checked and then converted f o r use by the routine C c a l c u l a t e _ f i l t e r ( ) . C C History: written 21-feb-1988 PK C - f o r now the parameter f i l e i s set to 'FILTER.PAR' C i m p l i c i t none l o g i c a l e x i s t s integer n,n_max,t,t_max,naroe_len,max_name_len,lun,status, + two,one,zero parameter(max_name_len=50,n_max=101,t_max=ll,zero=0,one=l,two=2) r e a l p s i , p i x e l _ s i z e parameter(pi=3.1415927) character*(max_name_len) filename,answer,newname,blank*l parameter(blank=• ') 1 include condition codes and declare external routines include '($FORDEF)' integer LIB$GET_LUN 1 get the f i l e containg f i l t e r parameters filename = 'FILTER.PAR' ! default f i l e name len = index(filename,blank)-1 10 continue inquire(file=filename(1:name_len),exist=exists,name=newname) i f ( . n o t . e x i s t s ) then type '(3a)',' Parameter f i l e ',filename(l:name_len), + ' does not e x i s t or i s i n use' 20 continue type ,(2a,S)*,' Input parameter filename ', + '(or return to quit) : ' read '(q,a)', name_len,answer i f ( a n s w e r ( l : l ) . e q . ' ') go to 1010 if(name_len.gt.max_name len) then type '(a , i 2 , a ) ' , ' F i l e name must be l e s s than ' + ,max_name_len,' characters' goto 20 endif filename(l:name_len) = answer(l:name_len) goto 10 else name_len = index(newname,blank)-1 filename(1:name_len)=newname(1:name_len) endif type '(/,2a)',' F i l t e r parameters from f i l e ', + filename(l:name_len) ! Open f i l e and read parameter values, l i n e s with a non-blank 1st character 1 are comment l i n e s and are skipped over status=LIB$GET_LUN(lun) i f ( . n o t . s t a t u s ) c a l l LIB$SIGNAL(%val(status)) open(unit=lun,file=filename(l:name_len),carriagecontrol='LIST', + status='OLD',err=1000) 123 c a l l skip_comments(lun,filename,name_len) read(lun,*,err=1000) n if(n.gt.n_max) then type l ( 2 a , i 3 ) , f ' Parameter n i s greater than max allowed', + * value of *,n_max goto 1000 elseif(mod(n,two).ne.one) then type '( a ) ' , ' Parameter n must be odd' goto 1000 endif c a l l skip_comments(lun,filename,name_len) readflun,*,err=1000) t if(t.gt.t_max) then type '(2a,i2)',' Parameter t i s greater-than max allowed 1, + • value of ',tjmax goto 1000 elseif(mod(t,two).ne.one) then type ' ( a ) ' , ' Parameter t must be odd' goto 1000 endif c a l l skip_comments(lun,filename,name_lenj read(unit=lun,fmt='(f5.2) 1,err=1000) p s i i f ( ( p s i . I t . 1 . 0 ) . o r . ( p s i . g t . 8 5 . 0 ) ) then type '(2a)',' Parameter p s i i s out of the allowed range', + ' of 1 to 85 degrees' goto 1000 endif c a l l skip_comments(lun,filename,name_len) read(unit=lun,fmt=*(f7.5)',err=1000) p i x e l _ s i z e i f ( p i x e l _ s i z e . l e . z e r o ) then type '(2a)',' Parameter p i x e l s i z e cannot be l e s s than', + • or equal to zero' goto 1000 endif 1 Convert values f o r use by routine c a l c u l a t e _ f i l t e r ( ) since values read ! i n are for e n t i r e reconstruction and we only meed t o c a l c u l a t e 1st. 1 quadrant and f o r theta_tau . l e . 90deg n = (n+l)/2 ! since we only need to do 1st- quadrant t = (t+l)/2 1 because of symmetry about theta=90deg p s i = psi*(pi/180.0) ! convert to radians goto 1020 1 skip e r r o r handling J E r ror handling 1000 continue 1 e r r o r e x i t point type '( a ) ' , ' F a t a l error opening or reading parameter f i l e ' 1010 continue ! non-error e x i t point type 1 ( a ) ' , 1 E x i t i n g f i l t e r generation program' c a l l e x i t 1020 continue ! normal end of routine close(lun) end 124 subroutine skip_comments (lun, filename, name_len) i 1 Lines with a non-blank f i r s t character are comment l i n e s . This p o s i t i o n s f i l e 1 so that the next non-comment l i n e i s the next one to be read i n . i m p l i c i t none l o g i c a l comment_line integer lun,name_len character f i1ename * 3 0 , i n _ l i n e * 2 6 5 comment_line =.true. do while (comment_line) read(unit=lun,fmt='(a)•,err=l000,end=1000) i n l i n e i f ( i n _ l i n e ( l : l ) . e q . • ') comment_line =.false, enddo 1 Go back 1 l i n e since we j u s t read a data l i n e backspace(unit=lun,err=1000) goto 1010 1 skip e r r o r e x i t 1000 continue ! error e x i t type •(3a)',' F a t a l error t r y i n g to skip over comment l i n e s i n ' , + ' parameter file',filename(1:name^len) type ' ( a ) ' , ' E x i t i n g f i l t e r generation program' c a l l e x i t 1010 continue ! non-error e x i t point end 125 subroutine make_filename(n,t,psi,pixel_size,filename,name_len) c C MAKE_FILENAME C Makes a filename using the parameters n, t, p s i , and p i x e l _ s i z e - i n that C order. p i x e l _ s i z e i s scaled so that p i x e l _ s i z e = l -> "D010" C and p s i i s i n degrees, n and t are set to the f u l l reconstruction values C C History: written 7-mar-1988 PK C i m p l i c i t none integer int_psi,n,t,int_px,name_len,max_name_len,index parameter(max name_len=50) r e a l p s i , p i , p T x e l _ s i z e , p i x e l _ s c a l e _ f a c t o r parameter(pixel_scale_factor=10.0,pi=3.1415927) character*10 char_psi,char px,char_n,char_t 1 integers (*4) have 10 d i g i t s character*(max_name_len) filename, b l a n k * l parameter(blank= * ') ! Convert values t o ASCII characters write(unit=char_n,fmt='(ilO)') 2*n-l i f (char_n(9:9).eq.• ') char_n(9:9) = '0' i f (char_n(8:8).eq.• ') char_n(8:8) = '0' write(unit=char_t,fmt='(110)') 2 * t - l i f (char_t(9:9).eq.• ') char_t(9:9) = '0' i n t _ p s i = nint(psi*(180.0/pi)) write(unit=char_psi,fmt='(ilO)') i n t _ p s i i f (char_psi(9:9).eq.' ') char_psi(9:9) = '0' int_px = n i n t ( p i x e l _ s i z e * p i x e l _ s c a l e _ f a c t o r ) ! t h i s w i l l cause values write(unit=char_px,fmt='(ilO) ' J int_px ! too large or small i f (char_px(8:8).eq.' ') char_px(8:8) = '0' ! to be written as "000" i f (char_px(9:9).eq.' •) char_px(9:9) = '0' i f (char_px(10:10).eq.• ') char_px(10:10) = '0' filename=•N•//cha r_n(8:10)//'T'//char_t(9:10)// + ,P'//char_psi(9:10)// ,D'//char_px(8:10)//'.SPE' name_len = index(filename,blank)-1 end 126 subroutine h a r d c o p y ( n , t , p s i , p i x e l _ s i z e , f i l t e r , p a r _ f i l e , + par_name_len,out_file,out_name_len) C C HARDCOPY C Outputs a formatted copy of the f i l t e r to a s p e c i f i e d f i l e . Each theta C index (angle) i s p r i n t e d on a new page along with the f i l t e r parameters C and f i l e names. C C History: written 4-mar-1988 c i m p l i c i t none integer n,t,ti,status,lun,i,j,max_width,name_len,index, + par_name_len,out_name_len,max_naroe_len parameter(max_name_len=50) r e a l f i l t e r ( 0 : n - 1 , 0 : n - 1 , 0 : t - 1 ) , p i x e l _ s i z e , p s i , d e l t a _ t h e t a , + pi,theta_tau,theta_bar,scale_factor,edgesum,penedge parameter(pi=3.1415927,scale_factor=l.0) character*(max_name_len) filename,answer,char1*1,par_flie, + out_file,blank*l,fullname parameter(blank=' ') ! Declare external routines and condition codes integer LIB$GET_LUN include '($F0RDEF) 1 filename = 'FILTER.LOG' ! Set to t h i s f o r now name_len = index(filename,blank)-1 ! Open f i l e and write f i l t e r data to i t , a new page f o r each angle status=LIB$GET_LUN(lun) i f ( . n o t . s t a t u s ) c a l l LIB$SIGNAL(%val(status)) open(unit=lun,file=filename(l:name_len),status= ,NEW) delta_theta = ( p s i / f l o a t ( t - l ) ) * ( 1 8 0 . 0 / p i ) ! convert to degrees ! Start b i g loop. do t i = 0 , t - l ! Calculate f i l t e r values f o r each theta_tau. theta_tau = 90.0 - ( p s i * f l o a t ( t i ) / f l o a t ( t - 1 ) ) * ( 1 8 0 . 0 / p i ) ! i n degrees theta_bar = acosd(cos(psi)/sind(theta_tau)) 1 c r i t i c a l angle " " ! Page header write(unit=lun,fmt='(a,//,a,f6.3,a)') + '1',' F i l t e r data f o r p s i = ',psi*(180.0/pi),' degrees' write(unit=lun,fmt='(/,2a,/,2a,/,a,i3,a,i3,a,i2,a,f6.3,2a,f7.4)') + • Parameters are from f i l e : •,par_file(l:par_name_len), + ' Data copied to section f i l e : ', + out_file(1:out_name_len), + • F i l t e r s i z e i s ',2*n-l,' x , , 2 * n - l , ' x \ 2 * t - l , + ' Delta theta_tau = 1,delta_theta,'deg.', + * P i x e l s i z e = ' , p i x e l _ s i z e write(unit=lun,fmt='(/,a,f5.2,a,f7.1)') + ' F i r s t quadrant of f i l t e r f o r theta_tau = ',theta_tau, + ' degrees. Data has been m u l t i p l i e d b y*,scale_factor write(unit=lun,fmt='(a,f5.2,//)') + ' C r i t i c a l angle theta_bar = *,theta_bar 127 1 P r i n t out f i l t e r data max_width «= min(n,12) do j = n - l , 0 , - l write(unit=lun,fmt='(/,lx,12,12(lx,g9.2))•) + j , ( f i l t e r ( i , j , t i ) * s c a l e _ f a c t o r , i = 0 , m a x _ w i d t h - l ) enddo write(unit=lun,fmt='(/,12(8x,i2))') + (i,i=0,max_width-l) enddo 1 End of Big Loop close(lun) ! Get the f u l l f i l e name in q u i r e ( f ile=filename(1:name_len),name=fullname) name_len «• index(fullname,blank)-1 filename(l:name_len) •= fullname(l:name_len) type , ( 2 a ) * , ' F i l t e r generation i n f o i s i n log f i l e *, + filename(l:name_len) end 128 Appendix E Reconstruction Program Listing 129 program recon C R E C O N C Performs a 3D image reconstruction by the method described i n my t h e s i s . C The p r o j e c t i o n data i s i n a 4D array i n the f i l e 'data.dat'. The C reconstruction f i l t e r i s i n a 3D array i n the f i l e ' f i l t e r . d a t ' . The C reconstructed image i s stored i n a 3D array i n the f i l e 'image.dat'. c c The f i r s t l i n e of each f i l e contains the necessary parameters f o r that C data set. The data i s a l l stored as unormatted 4-byte r e a l values. C The arrays are a l l handled with subroutines. Variables that begin with C d (for data) r e f e r to the data array. S i m i l a r i l y f o r f (for f i l t e r ) and D i (for image). c C 1. The parameters from the data and f i l t e r f i l e s are checked to be sure C that the arrays are of the r i g h t s i z e C 2. Space f o r s t o r i n g the 3 arrays i s created and the data and f i l t e r C are read i n . C 3. The arrays and the parameters are passed to reconloop(), which does C a l l the work. C 4. The image array returned from reconloop() i s copied to a f i l e and the C array storage space i s a l l released. C C History: 2nd version written 10-sep-1987 Paul Kinahan C C i m p l i c i t none l o g i c a l e x i s t s integer n,dn,dt,dp,fn,ft,dsize,fsize,isize,daddr,faddr,iaddr, + lunl/l/,lun2/2/,lun3/3/,status,wordsize/4/, + lib$get_vm,lib$free_vm,lib$init_timer,lib$show_timer, + lib$free_timer,timer_addr/0/ r e a l p s i , d p s i , f p s i character*(*) d a t a f i l e , f i l t e r f i l e , i m a g e f i l e p a r a m e t e r ( d a t a f i l e = • d a t a . d a t ' , f i l t e r f i l e = ' f i l t e r . d a t • , + imagefile='image.dat') status=lib$init_timer(timer_addr) ! s t a r t c o l l e c t i n g program s t a t i s t i c s , i f ( . n o t . s t a t u s ) c a l l l i b $ s i g n a l ( % v a l ( s t a t u s ) ) i n q u i r e ( f i l e = d a t a f i l e , e x i s t = e x i s t s ) i f ( . n o t . e x i s t s ) c a l l quit('No data f i l e 1 ) i n q u i r e ( f i l e = f i l t e r f i l e , e x i s t = e x i s t s ) i f ( . n o t . e x i s t s ) c a l l quit('No f i l t e r f i l e ' ) open(lunl,file=datafile,status='old',form='unformatted') read(lunl) dn,dp,dt,dpsi open (lun2 , f i l e = f i l t e r f i l e , status= • old', form= 'unformatted.') read(lun2) f n , f t , f p s i i f ( ( d n . n e . f n ) . o r . ( d t . n e . f t ) . o r . ( d p s i . n e . f p s i ) ) + c a l l q u i t ( ' F i l t e r and data f i l e parameters do not match') n = dn ! use these values - since they are the same as fn and f p s i p s i = dpsi d s i z e = (dn**2)*dp*dt*wordsize J s i z e of data array status = lib$get_vm(dsize,daddr) ! get storage space i f ( . n o t . s t a t u s ) c a l l l i b $ s i g n a l ( % v a l ( s t a t u s ) ) c a l l readdat(lunl,dn,dp,dt,%val(daddr)) c l o s e ( l u n l ) 130 I f s i z e <= ((2*fn)**2)*((ft+l)/2)*wordsize 1 s i z e of f i l t e r array status = lib$get_vm(fsize,faddr) ! get storage space i f ( . n o t . s t a t u s ) c a l l l i b S s i g n a l ( % v a l ( s t a t u s ) ) c a l l r e a d f i l ( l u n 2 , f n , f t , % v a l ( f a d d r ) ) close(lun2) i s i z e = (n**3)*wordsize ! s i z e of image array status = lib$get_vm(isize,iaddr) ! get storage space i f ( . n o t . s t a t u s ) c a l l l i b S s i g n a l ( % v a l ( s t a t u s ) ) c a l l zeroimage(n,%val(iaddr)) I zero image array c a l l reconloop(n,fn,dp,dt,psi,%val(daddr),%val(faddr), + %val(iaddr)) open(lun3,file=imagefile,status='new',form='unformatted *) write(lun3) n c a l l copyimage(n,%val(iaddr),lun3) close(lun3) status = lib$free_vm(dsize,daddr) ! free up a l l temporary storage i f ( . n o t . s t a t u s ) c a l l l i b S s i g n a l ( % v a l ( s t a t u s ) ) status = lib$free_vm(fsize,faddr) i f ( . n o t . s t a t u s ) c a l l l i b $ s i g n a l ( % v a l ( s t a t u s ) ) status = lib$free_vm(isize,iaddr) i f ( . n o t . s t a t u s ) c a l l l i b S s i g n a l ( % v a l ( s t a t u s ) ) type ' ( a ) ' , ' Program s t a t i s t i c s ' status=lib$show_timer(timer_addr) ! p r i n t out program s t a t i s t i c s i f ( . n o t . s t a t u s ) c a l l l i b S s i g n a l ( I v a l ( s t a t u s ) ) status=lib$free_timer(timer_addr) i f ( . n o t . s t a t u s ) c a l l l i b S s i g n a l ( % v a l ( s t a t u s ) ) end ! of main routine C A u x i l l i a r y routines used by main routine c subroutine q u i t ( o u t s t r i n g ) i m p l i c i t none character*(*) o u t s t r i n g type '(2x,a)', o u t s t r i n g c a l l e x i t end subroutine readdat(lun,nl,n2,n3,array) i m p l i c i t none integer lun,nl,n2,n3,i,j,k,1 r e a l array(nl,nl,n2,n3) read(lun) ( ( ( ( a r r a y ( i , j , k , 1 ) , i = l , n l ) , j = l , n l ) , k = l , n 2 ) , l = l , n 3 ) end subroutine r e a d f i l ( l u n , n , t , a r r a y ) i m p l i c i t none integer l u n , n , t , i , j , k r e a l array(2*n,2*n,0:(t-l)/2) read(lun) (((array(i,j,k),i=l,2*n),j=l,2*n),k=0,(t-l)/2) end C subroutine zeroimage(n,array) i m p l i c i t none 131 integer n , i , j , k r e a l array(n,n,n) do k=l,n do j=l,n do i=l,n a r r a y ( i , j , k ) = 0.0 enddo enddo enddo end subroutine copyimage(n,array,lun) i m p l i c i t none integer l u n , n , i , j , k r e a l array(n,n,n) write(lun) ( ( ( a r r a y ( i , j , k ) , i = l , n ) , j = l , n ) , k = l , n ) end 132 subroutine reconloop(n,fn,p,t,psi,data,filter,image) C C RECONLOOP C This f i l t e r s the data with the f i l t e r and returns the reconstructed C 3D image. C C The f i l t e r i n g and backprojection takes place one p r o j e c t i o n at a time. C Phi_tau ranges between 0 and 180 deg. and theta_tau ranges between 90-psi C and 90+psi. The sets of 2D f i l t e r data are indep. of phi_tau and are C symmetrical about 90 deg. i n theta_tau. C C Two temporary arrays are created f o r manipulating the data and are C passsed to the c a l l e d routines by t h e i r addresses. C C History : written 20-jun-1987 PK C i m p l i c i t none integer n, fn,p,t,pn,psize,paddr,rn,ressize,resaddr,status, + " done,tot,pindex,tindex,tmax,tmin, + lib$get_vm,lib$free_vm,wordsize/4/ r e a l d a t a ( n , n , 0 : ( p - 1 ) , - ( t - l ) / 2 : ( t - l ) / 2 ) , + f i l t e r ( 0 : ( f n - l ) / 2 , 0 : ( f n - l ) / 2 , 0 : ( t - l ) / 2 ) , + image(n,n,n), + phi_tau,delta_phi_tau,theta_tau,delta_theta_tau,pi,psi parameter(pi=3.141592653) pn = n + 2*fn - 2 ! dimension of padding array (pn x pn) psize = (pn**2)*wordsize ! s i z e of padding array status = lib$get_vm(psize,paddr) ! get storage space i f ( . n o t . s t a t u s ) c a l l l i b S s i g n a l ( % v a l ( s t a t u s ) ) c a l l zero2darray(pn,%val(paddr)) J i n i t i a l i z e to zero rn = n + fn - 1 ! dimension of r e s u l t array (rn x rn) r e s s i z e = (rn**2)*wordsize ! s i z e of r e s u l t array status = lib$get_vm(ressize,resaddr) ! get storage space i f ( . n o t . s t a t u s ) c a l l l i b $ s i g n a l ( % v a l ( s t a t u s ) ) c a l l zero2darray(rn,%val(resaddr)) ! i n i t i a l i z e to zero tmin = - ( t - l ) / 2 ! max value of theta_tau tmax = ( t - l ) / 2 ! min value of theta_tau delta_theta_tau = psi/float(tmax) delta_phi_tau = p i / f l o a t ( p ) t o t = t*p type , ( a , i 4 , a ) , , , There are ',t o t , ' projections to do' done = 1 type ' ( a , i 4 ) ' , ' doing projection '.done ! once now f o r overwriting. C Reconstruction loop: f i l t e r and backproject one p r o j e c t i o n at a time. C Because f i l t e r s have r o t a t i o n a l symmetry about z-axis, we a l t e r theta l a s t C i n the loop. The theta index (i) ranges from max. neg. value of theta (-psi) C to max. pos. value of theta (+psi) do tindex=tmin,tmax ! repeat f o r each d i f f e r e n t f i l t e r theta_tau = pi/2.0 - f l o a t ( t i n d e x ) * d e l t a _ t h e t a _ t a u do pindex=0,p-l phi_tau = float(pindex)*delta_phi_tau type '(a,i4)','+ doing p r o j e c t i o n '.done done «= done + 1 c a l l f i l t e r d a t a ( p i n d e x , t i n d e x , n , f n , p , t , d a t a , f i l t e r , 133 + pn,%val(paddr),rn,%val(resaddr)) c a l l backproject(rn,%val(resaddr),phi_tau,theta_tau,n,image) enddo enddo type '(a)',*+ Finished ' end ! of f i l t e r i n g routine C • C A u x i l l i a r y routines C subroutine zero2darray(n,array) i m p l i c i t none integer n , i , j r e a l array(n,n) do j=l,n do i=l,n array(i,j)=0.0 enddo enddo end 134 i subroutine filter_data(p_index,t_index,data_ndim,filter_ndim, + data_pdim,data_tdim,data,fliter,pad_ndim, + pad,result_ndim,result) c C FILTER_DATA C F i l t e r s the 2D r e a l data i n data(,,,) by f i r s t putting the data i n C the over-sized zero-padded array pad(,), and then performing a 2D C convolution with the array f i l t e r ( , , ) with the f i l t e r e d data C going in t o the array r e s u l t ( , ) . This e n t i r e operation i s done f o r C only one value of the higher dimension i n d i c i e s (t_index,p_index). C C Arrays:-data(,,,) i s dimensioned s . t . data(,,0,0) r e f e r s t o phi=0 and C theta=pi/2 C - f i l t e r ( , , ) i s dimensioned s . t . f i l t e r ( 0 , 0 , 0 ) r e f e r s t o theta=pi/2 C and x=y=0 C Note that f i l t e r _ t _ i n d e x = abs(t_index) C C Because of 4-fold f i l t e r symmetry, we only use 1st. quadrant of f i l t e r C values (thus the dimensioning of the array f i l t e r ( , , , ) ) . C C The r.atching of the array dimensions i s c r i t i c a l , but since t h i s routine C can be c a l l e d very often, we'll leave the dimensional checking to the C c a l l i n g routine C C History: written 20-jun-1987 PK C i m p l i c i t none integer data_ndim,data pdim,data_tdim,filter_ndim,pad_ndim, + r e s u l t ndim,t_Index,p_index,fti,pad_offset, + l i , l j , T , j , f i l t e r _ o f f s e t , f o C Arrays r e a l data(data_ndim,data ndim, 0:data_pdim-l, + -(data_tdTm-l)/2:(data_tdim-l)/2) r e a l pad(pad_ndim,pad_ndim) r e a l f i l t e r ( 0 : f i l t e r _ n d i m - l , 0 : f i l t e r _ n d i m - l , 0:(data_tdim-l)/2) r e a l result(result_ndim,result_ndim) C Abbreviations f t i = abs(t_index) ! since f i l t e r i s symmetric about 90 degrees pad_offset - 2 * f i l t e r _ n d i m - 2 ! using f u l l f i l t e r s i z e here f i l t e r o f f s e t = f i l t e r _ n d i m - 1 ! using f u l l f i l t e r s i z e here fo = f ! l t e r _ o f f s e t C Copy data i n t o centre of zero-padded array do j=l,data_ndim do i=l,data_ndim pad(i+pad_offset,j+pad_offset) = data(i,j,p_index,t_index) enddo enddo C F i l t e r do l j = l , r e s u l t _ n d i m do l i = l , r e s u l t _ n d i m ! f i r s t add c o n t r i b u t i o n at o r i g i n of f i l t e r r e s u l t ( l i , l j ) - f i l t e r ( 0 , 0 , f t i ) * p a d ( l i + f o , l j + f o ) 135 1 then add contributions from axes, using 2-fold symmetry along axes do i = l , f i l t e r _ n d i m - l ! x-axis f i r s t r e s u l t ( l i , l j ) = r e s u l t ( l i , l j ) + + f i l t e r ( i , 0 , f t i ) * ( p a d ( l i + f o - i , l j + f o ) + + pad(li+fo+i,lj+fo)) enddo do j = l , f i l t e r _ n d i m - l ! y-axis r e s u l t ( l i , l j ) = r e s u l t ( l i , l j ) + + f i l t e r ( 0 , j , f t i ) * ( pad(li+fo,1j+fo-j) + + pad(li+fo,lj+fo+j)) enddo ! f i n a l l y do the r e s t of the f i l t e r , only using the 1st. quadrant 1 because of 4-fold f i l t e r symmetry do j = l , f i l t e r _ n d i m - l do i = l , f i l t e r _ n d i m - l r e s u l t ( l i , l j ) = r e s u l t ( l i , l j ) + + f i l t e r f i , j , f t i ) * ( p a d ( l i + f o - i , l j + f o - j ) + + pad(li+fo-i,lj+fo+j) + + pad(li+fo+i,lj+fo-j) + + pad(li+fo+i,lj+fo+j)) enddo enddo enddo enddo end 136 ( subroutine backproject(m,plane,phi,theta,n,volume) C :--C BACKPROJECT C Backproject a 2D array (at a given theta and phi) int o a 3D array by C forward-projecting the center of each voxel onto the 2D array and C using b i l i n e a r i n t e r p o l a t i o n . See t h e s i s f o r explanation of the algorithm C C History: written 23-jun-1987 PK C-i m p l i c i t none integer m , n , l i , l j , i , j , k r e a l plane(m,m),volume(n,n,n),theta,phi,lx,ly,dx,dy,p,v p = 1.0 ! p i x e l s i z e . For now there i s no reason to d i f f e r from t h i s v = 1.0 I voxel s i z e . do k=l,n do j=l,n do i=l,n 1 f i n d forward projected point on plane l x = (v/p)*( -(i-(float(n)+1.0)/2.0)*sin(phi) + +(j-(float(n)+1.0)/2.0)*cos(phi) ) l y = (v/p)*( -(i-(float(n)+1.0)/2.0)*cos(theta)*cos(phi) + -(j-(float(n)+1.0)/2.0)*cos(theta)*sin(phi) + +(k-(float(n)+1.0)/2.0)*sin(theta) ) ! f i n d nearest gridpoints l i = i n t ( l x + (float(m)+1.0)/2.0) l j <= i n t ( l y + (float(m)+1.0)/2.0) ! i n t e r p o l a t i o n parameters dx = (lx - f l o a t ( l i ) ) + (float(m)+1.0)/2.0 dy = ( l y - f l o a t ( l j ) ) + (float(m)+1.0)/2.0 ! add inte r p o l a t e d value to voxel volume(i,j,k) = volume(i,j,k) + + (1.0-dx)*(1.0-dy)*plane(li,lj) + + (dx)*(1.0-dy)*plane(li+l,lj) + + (1.0-dx)*(dy)*plane(li,lj+l) + + (dx ) * ( d y ) * p l a n e ( l i + l , l j + l ) enddo enddo enddo end 137 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0085037/manifest

Comment

Related Items