PACT The Salzburg NTN-Method for the Radon Transform Deliverable D5Z-2 (Contract GZ 308.930/1-IV/3/93) Responsible Partner: Research Institute for Softwaretechnology, University of Salzburg Contributing Partners: Authors: Erika Hausenblas Version: 1.0 Date: October 31 1994 Status: release Con dentiality: con dential PACT Abstract The radon transform arise in natural way, if X-rays traverses a body, like in a computer tomograph. Thus, to nd out the internal structur of a body in computer tomography corresponds to reconstructions a function from their radon transform. This report deals with a parallized numerical reconstruction algorithm, in particular the ltered backprojection. The ltered backprojection uses high dimensional integration. We test numbertheoretical methods of Hlawka and Koksma applied to the algorithm. Further, to decreasing the executing time, we implemented a parallelized the algorithm. D4Z-2/Rel 1.0/October 31 1994 1 CONTENTS PACT Contents 1 Introduction 2 The Radon Transform 2.1 2.2 2.3 2.4 2.5 The Radon Transform : : : : : : : : : : : : : : : : : : The Filtered Backprojection Algorithm : : : : : : : : Numbertheoretical Methods in Numerical Integration : The In uence of the Parameters : : : : : : : : : : : : Applications of the Radon Transform : : : : : : : : : : 3 4 : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 3 The Software Tool: RADON 3.1 A short Overview of the Objects 10 : : : : : : : : : : : : : : : : : : : : : : : : : 4 Parallel Virtual Machine (PVM) 4.1 Implementation with PVM and load balancing : : : : : : 4.2 PVM and load balancing applied to the program package 4.3 The Improve of the Executing Time : : : : : : : : : : : : 5 Results 6 Conclusions 7 References D4Z-2/Rel 1.0/October 31 1994 4 5 7 8 8 10 12 : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 12 13 13 14 17 18 2 Introduction 1 PACT Introduction The computerized tomography is a technique to reconstruct internal structurs of a body from data collected outside the body. This means, computer tomography provides a very good method to investigate the internal structur of an object without destroying it. In particular this method is very important in medicine. Theoretically it is possible to reconstruct the internal structur without errors. There is only one problem. If the values of the radon transform have been obtained by physical measuremant, one cannot assume, that the data are complete. Therefore it is very important to nd algorithm which yields to good approximation of the original object in spite of discrete or incomplete datas. This paper deals with a certain reconstruction algorithm, the ltered backprojection. To reconstruct a point, a lter function is convoluted by the radon transform. Thus, the result is determined by numerical integration. The executing time of the numerical implementation can be improved by parallelizing. In the program this is implemented by using PVM. Further, the common algorithm uses equidistant points for the integration. But, instead of equidistant points, one can use good lattice points explored by Hlawka and Koksma. In this paper we test a new lter function combined with numerical integration based on good lattice points and compare the result with other algorithms. The paper can be divided in the following sections: First we explore the radon transform and the algorithm we deal with. Since it is not possible, we only give a motivation and no complete introduction. In the next section we present the software tool written for research and testing the new algorithm. Because of the executing time, some parts of the programm are written in PVM. The two last sections deal with the result. We present the picture, compare the algorithm and give at least some conclusions. D4Z-2/Rel 1.0/October 31 1994 3 The Radon Transform PACT 2 The Radon Transform 2.1 The Radon Transform The radon transform of a three-dimensional body is the parallel projection of the denseness function onto a plane. For simplicity in this paper the radon transform will be only treated for two-dimensional objects. To generalize the radon transform, you have to take into account the fact that a two-dimensional object is a crosscut of a three dimensional one. So in order to obtain a three dimensional image one scans the body layer by layer. A two dimensional object can be described by a function f (x; y ) : IR  IR ! IR+ , where f (x; y ) is the denseness of the object in the point (x; y ). The Radon transform R (s) of a point (; s) is given by the parallel projection onto the corresponding line L at an angle  with the x-axis, evaluated at a distance s from the origin. The projection onto the line L can be described by the following integral: p (s) = Z (L (s))? f (x(u); y (u)) du where L? denotes a line orthogonal to L and through the point L (s), s the distance to the origin. Therefore each point of the state space [0;  ]  IR corresponds to a value p (t). The Radon transform of the function f (x; y ) is the two-dimensional function (Rf ) (; t) : [0;  ]  IR + ! IR , with (Rf ) (; t) = p (t) The problem is to reconstruct the object, respective its denseness function f (x; y ), given the radon transform. The solution is given by the inversion formula of the radon transform. Johann Radon developed an explicit inversion formula in 1917 [Ra17]. If there are available all values of R (t), the inversion formula works quite well. But this formula leads to numerical instabilities, if there are only a nite number of lines L. The problem now is, that in physics the condition on a continuous sample space is violated. That means that one can measure the Radon transform only on a nite number of points, for example on a lattice Mn = f(si ; i ) j i = 1


ng. Further the datas are often incomplete or it is expensive to get a great number of datas. Therefore it is of practical importance to nd algorithms which yield a good approximation of the denseness function, utilizing only few number of points. D4Z-2/Rel 1.0/October 31 1994 4 The Radon Transform PACT 2.2 The Filtered Backprojection Algorithm The radon transform (Rf )(;
) is a parallel projection onto the line L through the origin with the angle . Along the line L the values of the denseness function are integrated. To test the algorithm, we rst create a phantom picture of superposed ellipses and then we determine the radon transform. The next gure shows an example of the radon transform of an ellipse. The angle, drawn in the ellipse is 6 . 1 0.75 0.5 0.25 0 -0.25 -0.5 -0.75 -1 -0.75 -0.5 -0.25 0 0.25 0.5 0.75 1 The ellipse 400 300 1 200 100 0 0 s -2 0 angel The corresponding radon transform -1 2 Note, that the value of the denseness function in the point (x; y ) appears at the point on the line L of a xed angle . The function value f (x0 ; y0) at a point (x0 ; y0) can be represented by a convolution of f with a delta function1 . If the function is available only on s() 1 The delta function is roughly spoken zero everywhere but of one point x . In this point, the function is 0 not de ned resp. in nite. The most important property of a delta function x0 is g   xo D4Z-2/Rel 1.0/October 31 1994 = Z 0 1 g(x)x0 (x) dx = g(x0 ) 5 The Radon Transform PACT a nite set of points fxi ; i = 1


ng, one can approximate the function g by convolution of the sample function with an approximated delta function. ! n X 1 g ~ (x) =   g (x ) (x) appr appr n i i Most interpolation formulas can be written as a convolution with an approximation of the delta function. It can be compared with smearing the exactly value in the neighborhood of each point. Such an approximative delta function is called low pass lter. Now, the idea behind the ltered backprojection is to convolute the reconstructed function with an appropriate lter. The trick is to convolute rst the radon transform with the lter, and then by the second step to reconstruct the function by the exact formula of Johann Radon. This can be done, because both operations, calculating the radon transform and convolution, commute. This signi cates, that it doesn't matter to go in the left diagram rst right and then down or rst down and then right. Because the inversion radon transform is unique, the second diagram commutes also. The advantage is that the convolution of two radon transform is just a multiplication and much more quicker to determine than the convolution of the functions themselfs. f Rf f*d Rf f R(f*d)= (Rf)*(Rd) (Rf)*(Rd) f*d Thus, to reconstruct the value, rst the radon transform is multiplicated by a lter and then the inversion formula of J. Radon is applied: Z  Z p2 (1) f~(x; y ) = p g(x;y) (; s) (Rf )(; s) ds d 0 2 where g(x;y) is the kernel of the lter. The kind of lter function depends on the algorithm and vice versa. We have implemented two lter functions: 1. Filter b = 15 0.5 0.4 0.3 g1 0.2 = cos b(2bzz)2? 1 0.1 0 -1 -0.5 0 0.5 1 where g : IR ! IR is an arbitrary continuous function. D4Z-2/Rel 1.0/October 31 1994 6 The Radon Transform PACT and 2. Filter b = 15 60 40 20 = g2 s 0 b2 Z 1 0 cos (zbu) u(1 ? u2 ) du -20 -1 -0.5 0 0.5 1 where z = zx;y (; s) = x sin() + y cos() ? s is a function of (; s) and denotes the distance of L (s) to the intersection point of P. The rst smears the value around the point much more than the second. Thus, the reconstruction picture by the rst lter is much more blury than the of the second lter. Chapter 2.4 deals with the in uence of the bandparameter. 2.3 Numbertheoretical Methods in Numerical Integration The integration in 1 can only be determined by numerical h p p imethods, because it is impossible to measure the values of (; s) on all points in ? 2; 2  [0;  ]. Usually, this numerical integration is done on a grid with equidistant points: i si = ni  i = 1


n ; p = (2 ni ? 1) 2 fappr (x; y ) = 1 1 n ns n X ns X i=1 j =1 i= 1


ns ; R g(x;y)(i ; si )( f )(i; si) In the theory of numerical integration the method of good lattice points is well known. Good lattice points are a sequence of points with small discrepancy. It exceeds the scope of this paper to point out numbertheoretical methods in numerical integration. For a good introduction, please see: [Hl61], [Hl79].In this paper we used good lattice points created by Fibonacci numbers2 .   i Fn?1 i =  i = 1


Fn ; Fn 2 The n'th Fibonacci number is given by the following recursion: F0 = 1; F1 = 1; Fn = Fn?2 + Fn?1 D4Z-2/Rel 1.0/October 31 1994 7 The Radon Transform PACT p = (2 Fi ? 1) 2 i = 1


F ; If we use good lattice points instead of equidistant points in the numerical integration, the error will decrease. si n n 2.4 The In uence of the Parameters To illustrate the in uence of the band parameter, in appendix 1 a sequence of reconstructions is presented, varying only the band parameter. Varying the band parameter corresponds to the approximation the  function better or worse. The lter is stretched by increasing the band parameter and is contracted by decreasing the parameter. 1 0.5 10 g2 0 8 -0.5 -1 6 b -1 4 0 s 2 1 Thus, if there are only a few radon transforms available, the band parameter must be choosen low. If it can be assumed, that there are many points of the radon transform, the band parameter can be choosen higher. Thus the reconstructed picture is bleary. If the band parameter increases, also the accurancy of the pictures increases. But the artefact caused by the numerical method increases as well. Further depends the band parameter on the lter kernel. For example, the rst lter needs a higher band parameter than the second. The in uence is easy seen in the pictures, which are presented in appendix 1. 2.5 Applications of the Radon Transform The most common application of the radon transform is computerized tomography. The radon transform arises in a natural way, if X-rays traverse a body. A cross section of the D4Z-2/Rel 1.0/October 31 1994 8 The Radon Transform PACT human body is scanned by a thin X-ray beam whose intensity loss is recorded by a detector. Let f (s) be the denseness of the tissue along the line L? (s); s 2 IR. The amount of absorbed intensity
I in a small distance
s at the point L? (s) is proportional to
I =I = f (s)
s . Therefore the intensity of the beam after having passed the body along L? is X-Rays I0 Body = exp ? Z L?  ( ) ds ! f s Detector Thus the reconstruction problem of computerized tomography is just the inversion problem of the radon transform. Further over the years the physical method of computerized tomography became to be used in a wider sense, namly for any technique of reconstruction internal structures in a body or tissue from data collected by detectors outside the body, for example in electron microscopy and electromagnetic scattering. Therefore computerized tomography is a method which is used in many elds like physics, engineers, bioscientists, medicin and many others. D4Z-2/Rel 1.0/October 31 1994 9 The Software Tool: RADON PACT 3 The Software Tool: RADON The program is built up like a library. The aim is to give a tool for research in the radon transform, especially to investigate in numbertheoretical approaches. To test an algorithm it is necessary to know the Radon transform exactly. Further it is important, that the radon transform can be quickly calculated. This conditions are given by a picture of superposed ellipses, called the phantom. The third advantage of superimposed ellipses is that they look like a crosscut of a brain. This picture is called phantom picture, because it is not real. Thus the test algorithm is divided in the following steps: 1. Create the Phantom Picture De ne the appropiate ellipses 2. Draw the Phantom Picture 3. Transform the Phantom Picture Evaluate the radon transform 4. Reconstruct the Object 5. Determine the Error The error of the reconstructed picture is the 2-norm of the di erence of the original. Note that after reconstructing the values of the denseness function, the picture will be prepared for visualization. This will be done by stretching the range of the function on the intervall [0; 255]. Therefore, rst the range of the phantom picture will be prepared like the reconstructed picture and afterwards the error will be determined. vuu 1 1 X X ( ( =t nx error nx ny ny i=1 j =1 f x i ; yj ) ? frec (xi ; yj ))2 (2) The phantom picture, the picture and the picture of the radon transform are objects of the program package. The procedures to execute these steps are implemented as functions of the corresponding object. For example the procedure phantom:: radon transform( radon picture ) determines the radon transform and saves the values in the picture radon picture. attention: We restrict the range of the phantom picture on [?1; 1] or [0; 1]. 3.1 A short Overview of the Objects This section contains a short description of the objects and corresponding functions. 1. Object Phantom: The object Phantom is roughly spoken an array of ellipses. Further, the object possesses a number of functions, e.g. with the function add() you can add one ellipse to the phantom. Thus you create your phantom picture by adding ellipses one by one. The D4Z-2/Rel 1.0/October 31 1994 10 The Software Tool: RADON PACT next step is to draw the picture. This can be done by the function draw( picture ). There are also functions to determine the radon transform and to assign the radon transform to an appropiate object. (a) draw( picture): draws the phantom in the given pictures (b) pvm draw( picture): the same procedure like above. the only di erence is, it uses PVM (c) radon transform( radon picture): determines the radon picture (d) pvm radon transform( radon picture): determines the radon picture using PVM 2. Object ellipse: an ellipse is described by its middle point, axis and angle. Initializing an ellipse can be done by: ellipse ell1(m1,m2,a,b,phi,color) or if the ellipse exists already: ell1.initialize(m1,m2,a,b,phi,color). The other functions, like radon transform( radon picture), of the object are private and are used by functions of the object phantom. b m2 b Phi m1 3. Object picture: It is just an array of integer values. Further, a subclass belongs to the object. (a) error(): determines the error described in equation 2 (b) print(int,int): creates a le, which can be visualized by the program tool xv the subclass lter: The lter contains also a second array, but of double, with the same size like the other array. The functions are: D4Z-2/Rel 1.0/October 31 1994 11 PACT Parallel Virtual Machine (PVM) (a) reconstruct( radon picture ): reconstructs a picture of the given radon tramsform and entries the values to the double array. (b) entry(): carries over the array of double into the array of integer and prepares the values to visualize. This means, the range is streched on the range [0; 255]. There are some possibilities implemented for visualization. attention: the error function is calculated by the visualized picture. 4. Object radon: The object radon stores the evaluated values of the radon transform. A subclass belongs also to the object: the object radon rand. This subclass possesses the same functions as the superclass but of two initialization procedure. It treats the case of good lattice points. The most functions are overloaded. The functions are: (a) Fibonacci(): belongs to radon rand. Determines the good lattice points based on Fibonacci number and stores it. (b) randomize(): belongs to radon rand. Throws random points. (c) set zero(int,int): sets the given point in the array to zero. If the rst entry is minus one, the procedure sets all the values of the given value by the second argument to zero and vice versa. By random rand you can only choose the number of points which should be set to zero. 4 Parallel Virtual Machine (PVM) The Parallel Virtual Machine is a software system that enables a collection of heterogeneous computers to be used as a coherent and exible parallel computational resource. Especially it can be used for the development and execution of large concurrent or parallel applications that consist of many interacting, but relatively independent components [Ge93]. PVM was developed initially at the Emory University and Oak Ridge National Laboratory. The individual computers may be shared- or local-memory multiprocessors, vector supercomputers or scalar workstations, that can be interconnected by a variety of networks (ethernet, FDDI, etc.). User programs written in C, C++ or Fortran are provided access to PVM through the use of calls to PVM library routines for functions such as process initiation, message transmission and reception. The PVM system handles message routing, data conversion for incompatible architectures and other tasks that are necessary for operation in a heterogeneous network environment. For more details see [Ge92]. 4.1 Implementation with PVM and load balancing Since the radon transform or reconstruction for di erent points (s ;  ) of the lattice for a point (x ; y ) of the picture are evaluated independently, the outer loop may be distributed among i i i i D4Z-2/Rel 1.0/October 31 1994 12 Parallel Virtual Machine (PVM) PACT the processors. On a homogeneous network, where the user disposes 100% of the machines' capacity on each node, this can be done by partition the lattice into equal parts for each node. If there are other processes running or if a heterogeneous network is used, a processor with high workload or a slower one would cause a delay in the whole computation. With a ner partitioning of the lattice it is possible to react on changes in the systems performance. Therefore the algorithm is implemented using the master{slave paradigm and the asychronous pool of tasks methodology. 4.2 PVM and load balancing applied to the program package The program of determining the radon transform of a given phantom and reconstructing the picture contains three time consuming procedures: phantom::draw( ), phantom::transform( ) and radon::reconstruct( ). Therefore for each procedure there exists a version using PVM with load balancing. The rough structure of the functions are given below: 1. phantom::pvm draw: First the master send the size of the picture to each slave. The next step is to distribute the ellipses among the slaves. The task of each slave is to draw the received ellipse. After he nished, he will send the picture back to the master. If there are ellipses left, he gets the next one. The master collects all pictures and sums up all together. 2. phantom::pvm radon transform: The task is divided in exactly the same manner like in phantom::pvm draw. Each slave calculates the radon transform of the given ellipse. The master adds one by one to his own radon picture and orders the slave to calculate another ellipse if it is necessary. 3. lter::pvm reconstruct: By reconstructing a picture, the denseness function of each point in the unit square can be reconstructed independently. Therefore each slave can determine the value of the denseness function in one point. But the amount of communication time increases if the tasks are to small. Because of this each slave gets the task to reconstruct a column of the picture. First the master distributes the radon transform and the size of the picture to reconstruct to all. After this each slave gets the number of the column he has to evaluate. The next step, after a slave nished his task, he will send back the result and get a new number of column if there is some left. 4.3 The Improve of the Executing Time Because all time consuming functions are also implemented in PVM, you can neglect the amount of time needed by the other procedures. Further, the time a program using PVM depends on the number of hosts. Thus, the executing time of a program using PVM is in proportion to the time without PVM like one to the number of hosts. Evaluationg a picture with 150x150 pixels, 82x82 points of radon transform needs 48 minutes, using PVM and four hosts, only 22 minutes. D4Z-2/Rel 1.0/October 31 1994 13 Results 5 PACT Results To test the algorithm we reconstuct two sequences of pictures. The pictures of the rst serie are approximated by an algorithm described in the PhD-thesis of Michael Revers (algorithm 1). The pictures of the second serie are evaluated by an new method which is based on an article of Michael Revers (algorithm 2). Further, one serie can be divided into three kinds of reconstruction. At each reconstruction we used di erent grids for the numerical integration.  Equidistance Points  = ni  i = 1


n; i p s = (2 ni ? 1) 2 i = 1


n; i  Randomly Choosen Points  = ni  i = [rand()=RAND MAX ]; p s = (2 ni ? 1) 2 i = [rand()=RAND MAX ]; i i  Good Lattice Points   i F ? 1 = F  i = 1


F ; p s = (2 Fi ? 1) 2 i = 1


F ; n i i n n n n Note, that the angle is only running from 0 to  . The points of [; 2 ] are given by the symmetric relation (Rf )(; s) = (Rf )( + ; ?s). Thus, in the program we identi ed the two points and reconstruct the pictures where  runs over [0; 2 ]. For testing the algorithm we create two phantom pictures, phantom 1 and phantom 2. You can see the pictures in appendix 2. Further, there are presented pictures, reconstructed by equidistant points, random choosen points and Fibonacci points. We took the choice of the best band parameter. Thereas the band parameter was assumed to be eleven. Between 11 and 15, the error is minimized. D4Z-2/Rel 1.0/October 31 1994 14 Results PACT fibonacci points equidistance points 5000 5000 4000 4000 3000 3000 2000 2000 1000 1000 0 1 2 4 3 5 6 7 random points 0 8000 8000 7000 7000 6000 6000 5000 5000 4000 4000 3000 3000 2000 2000 2 4 3 5 7 6 all 1000 1000 0 1 1 2 3 4 5 6 7 0 0 1 2 3 4 5 6 7 Attention: Please, note that the band parameter 1 in the picture corresponds truely the band parameter 10, further 2 corresponds 11 and so on. If you look at the corresponding pictures you will notice, that to make out the object reconstructed by randomly choosen points is much more dicult than the others. By physical damage, it can happen that some measurement points fail. To simulate this we have set some points of the radon transform to zero. 1. First, we set all points of one angle to zero. This could be regarded as the failure of all beams at one xed angle. We took two di erent angles: 10 and 30 degrees (pixels 150x150, 82x82 points and b=13). The picture is presented in appendix 3. The same D4Z-2/Rel 1.0/October 31 1994 15 Results PACT event cannot happen by Fibonacci numbers. Thus, to simulate an equivalent event, we erase 82 points of the radon transform (pixels 150x150, 8765 points and b=13). 2. Here only one point is set to zero. You can see the picture also in appendix 3. We took the same choice of the number of pixels, points and band parameter. To compare the robustness of the two algorithms, we have done the same with a picture, reconstructed by Fibonacci numbers. D4Z-2/Rel 1.0/October 31 1994 16 Conclusions 6 PACT Conclusions To research the di erent between algorithm using bonacci numbers and equidistant points is given in the next paper. the two advantages which arise by good lattic points points are given below. 1. To improve the error rate, it is necassary to add some points. Using good lattice, you can add points one by one. This means, if you have already n measure points, you can add m other ones, where the m can be choosen arbitrary. If the points are equidistant and you want n + m points, you have to choose between two possibilities: to cast the old points and to measure new n + m points or, to have a grid, which consists no more of equidistant points. To choose the rst alternative is to waste a lot of time. If you choose the second, there will appear artefacts because the domain of integration is no more treated equal. 2. If some points fail because of a physical damage or something else, the algorithm using Fibonacci numbers yields to a lower error than the algorithm using equidistant points. Thus, the Fibonacci numbers lead to an algorithm which is much more robust than the other. The cause is the absence of a regular structure. Looking at a grid of equidistance points, you will note, that you can rotate it only in angles of 4 ; 2 ; 34 ;  . If you single out a circle of a lattice, generated by Fibonacci numbers, you will not observe a signi cant di erence, rotationg the circle in any angle you want. The circles are all symmetric. Thus, one point fails, there is no regularity or high structure destroyed. Acknowledgment: I am gratefull to Andreas Uhl for his assistance and advices in UNIX, printing images and PVM. Further, I have to thank Michael Revers for his help. D4Z-2/Rel 1.0/October 31 1994 17 References 7 PACT References Na85 Natterer, Frank: The Mathematics of Computerized Tomographie, John Wiley and Sons (1985) Ge93 Geist, A.: PVM3 beyond network computation In J. Volkert, editor Parallel Computing, Lectures Notes in Computer Sience 73, p. 194-203, Springer (1993) Ge92 Geist, A. and V.S. Sunderam: Network based on concurrent computing on the PVM system: Concurrency: Practice and Experience 4(4), p.293-3115 (1992) He80 Helgason, S.:The Radon Transform(1980) Birkhaeuser Berlin Hl61 Hlawka, E.: Funktionen von beschraenkter Variation in der Theorie der Gleichverteilung, Annali di Mathematica,(1961) Hl79 Hlawka, E.: Theorie der Gleichverteilung, Annali di Mathematica,(1979) Ra17 Radon, Johann: Ueber die Bestimmung von Funktionen durch ihre Intergralwerte laengs gewisser Mannigfaltigkeiten (1917)Berichte der Saechsische Akademie der Wissenschaften, Leipzig, Math.-Phys., p. 262-267,69 Re90b Revers, Michael: Anwendung zahlentheoretischer Methoden fuer die Numerik der Computertomographie, PhThesis,(1990), Salzburg Re94 Revers, Michael: Numeric Integration of the Radon Transform on Classes Es in multiple ( nite) dimensions (1994), preprint Sa88 Sanz, Jorge L.C. and Eric B. Hinkle and Anil K.Jain: Radon and Projection TransformBased Computer Vision, Springer Series in Information Sciences (1988) D4Z-2/Rel 1.0/October 31 1994 18