DEVELOPMENT OF A LOW COST X-RAY COMPUTED TOMOGRAPHY SYSTEM By Ahmed Zaki Alsinan A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Electrical Engineering—Master of Science 2013 ABSTRACT DEVELOPMENT OF A LOW COST X-RAY COMPUTED TOMOGRAPHY SYSTEM By Ahmed Zaki Alsinan Computed Tomography (CT) has been one of the most exciting imaging modalities developed in the biomedical field. Although the problem of imaging and image reconstruction arises in many scientific fields, its most prominent applications are found in the field of diagnostic medicine. Imaging from x-ray projections is the process of reconstructing an image of the two-dimensional distribution of x-ray attenuation coefficient from measurements of projection data. Based on the well-tested Filtered Back Projection approach, a new CT imaging system is developed to enable the reconstruction of objects from multiple projection data obtained from a GE vertical-stationary x-ray imaging system. Although this system is not designed for CT purposes, in this thesis we present a low cost scheme for obtaining multiple view projection, from which to reconstruct the object. There are many artifacts introduced by the instrumentation of the system for CT imaging. These issues are investigated and algorithms are developed to compensate for these artifacts caused mainly by mechanical errors during rotation. Modifications in adaptive filtering algorithms are also incorporated for further improving the reconstruction results. Thereby, an accurate and low-cost CT imaging system is realized. ACKNOWLEDGMENTS I would like to extend my sincere gratitude to my advisor Dr. Lalita Udpa for her tremendous guidance and support throughout my studies here at Michigan State University. I thank her and the committee members Dr. Satish Udpa and Dr. Robert McGough for their encouragement, motivation and constructive criticism. The assistance and innovation of Mr. Brian Wright, Research Equipment Technologist of ECE department, is greatly acknowledged. Many thanks go to my colleagues in the Nondestructive Evaluation Laboratory. In particular, I am grateful to Seyed Morteza Safdarnejad for his suggestions and discussions in signal processing. I wish to express my heartfelt gratitude to my mother Iftekhar AbuAlsaud and my father Zaki Alsinan, who both have been the greatest inspiration in my life. I owe special thanks to my siblings for their encouragement. I cannot articulate how thankful I am to my wife Zahraa for her love and care. iii TABLE OF CONTENTS LIST OF FIGURES .......................................................................................................................................... vi Chapter 1 Introduction .............................................................................................................................. 1 1.1 Overview .......................................................................................................................................................... 1 1.2 Scope and Objectives .................................................................................................................................. 3 Chapter 2 Background ............................................................................................................................... 5 2.1 CT Physical Principles ................................................................................................................................ 6 2.2 CT Mathematical Principles .................................................................................................................... 8 2.2.1 Mathematical Background ........................................................................................................... 8 2.2.2 Radon Transform ......................................................................................................................... 11 Chapter 3 Theory ........................................................................................................................................ 16 3.1 Image Reconstruction Definition ....................................................................................................... 16 3.2 The Fourier Slice Theorem ................................................................................................................... 18 3.2.1 Proof of the Fourier Slice Theorem ....................................................................................... 19 3.3 Filtered Back Projection Algorithm .................................................................................................. 23 3.3.1 FBP Filtering................................................................................................................................... 26 3.3.2 FBP Numerical Implementation ............................................................................................. 28 3.4 Sampling Geometry ................................................................................................................................. 30 3.4.1 Parallel Beam ................................................................................................................................. 30 3.4.2 Fan Beam ......................................................................................................................................... 32 Chapter 4 Experimental Setup .............................................................................................................. 36 4.1 Imaging System ......................................................................................................................................... 36 4.2 Rotation Control System........................................................................................................................ 38 4.3 CT Reconstruction Algorithm .............................................................................................................. 40 4.4 Experimental Data ................................................................................................................................... 44 Chapter 5 Results and Discussion ....................................................................................................... 48 5.1 Overview ....................................................................................................................................................... 48 5.2 Experimental Data Initial Results ..................................................................................................... 49 5.2.1 First Object...................................................................................................................................... 49 5.2.2 Second Object ................................................................................................................................ 50 5.2.3 Third Object .................................................................................................................................... 51 5.3 Problems & Issues ..................................................................................................................................... 53 5.3.1 Mechanical Issues......................................................................................................................... 53 5.4 Experimental Data Final Results ....................................................................................................... 68 5.4.1 First Object...................................................................................................................................... 68 5.4.2 Second Object ................................................................................................................................ 69 5.4.3 Third Object .................................................................................................................................... 70 Chapter 6 Conclusions and Future Work ......................................................................................... 72 iv 6.1 6.2 Conclusions .................................................................................................................................................. 72 Future Work................................................................................................................................................ 73 APPENDIX ...................................................................................................................................................... 74 BIBLIOGRAPHY ............................................................................................................................................ 81 v LIST OF FIGURES Figure 2.1 Physical Model ............................................................................................................................. 7 Figure 2.2 Two-dimensional Rectangular Function (left) and its CSFT (right) .................... 10 Figure 2.3 Line Representation in Polar Coordinates ..................................................................... 12 Figure 2.4 Line Representation in Polar Coordinates 2 ................................................................. 13 Figure 2.5 Line Representation in Polar Coordinates 3 ................................................................. 13 Figure 2.6 Coordinate Rotation Interpretation ................................................................................. 14 Figure 3.1 Projection Interpretation...................................................................................................... 17 Figure 3.2 Projection in Frequency Domain ....................................................................................... 19 Figure 3.3 Fourier Slice Theorem Diagram ......................................................................................... 22 Figure 3.4 Rectangular vs. Polar grids................................................................................................... 23 Figure 3.5 Band limited High-Pass Filter ............................................................................................. 26 Figure 3.6 Convolution of Two Functions............................................................................................ 27 Figure 3.7 Numerical Example ................................................................................................................. 28 Figure 3.8 Filtered Back Projection Performed on Numerical Example .................................. 29 Figure 3.9 Parallel Beam Geometry........................................................................................................ 31 Figure 3.10 Fan Beam Geometry................................................................................................................ 31 Figure 3.11 Fan Beam Projection I ............................................................................................................ 33 Figure 3.12 Fan Beam Projection II .......................................................................................................... 35 Figure 4.1 General Electric Vertical X-Ray Inspection System .................................................... 37 Figure 4.2 Stepper Motor Controller ..................................................................................................... 38 Figure 4.3 X-Ray Rotation Stage Control System .............................................................................. 39 vi Figure 4.4 Synthetic Data – Phantom .................................................................................................... 40 Figure 4.5 High-Pass Filter Applied to Projections .......................................................................... 41 Figure 4.6 Phantom Sinogram (top) and Filtered Sinogram (bottom)..................................... 42 Figure 4.7 Phantom Reconstruction – Projection Number is indicated................................... 43 Figure 4.8 Objects Scanned: (a) side view and (b) top view......................................................... 44 Figure 4.9 First Object ................................................................................................................................. 45 Figure 4.10 Second Object ............................................................................................................................ 46 Figure 4.11 Third Object ............................................................................................................................... 47 Figure 5.1 First Example (a) Projection at 00 (b) Sinogram ......................................................... 49 Figure 5.2 Second Example (a) Projection at 00 (b) Out-of-Phase Sinogram......................... 50 Figure 5.3 Third Example (a) Projection at 00 (b) Out-of-Phase Sinogram ............................ 52 Figure 5.4 Phantom ...................................................................................................................................... 54 Figure 5.5 Phantom Sinogram .................................................................................................................. 55 Figure 5.6 Phantom Sinogram – Shifted by 10 and 6 pixels ......................................................... 56 Figure 5.7 Phantom Sinogram – Shifted by 5 and 3 pixels ............................................................ 56 Figure 5.8 Reconstruction Result of Shifted Sinogram by 10 and 6 pixels ............................. 57 Figure 5.9 Reconstruction Result with +3% Rotation Error ........................................................ 58 Figure 5.10 Reconstruction Result with -5% Rotation Error ......................................................... 59 Figure 5.11 Sinogram From Rotation Error with σ=0.05 ................................................................. 60 Figure 5.12 Reconstructing Result of Corrupted Projections ......................................................... 60 Figure 5.13 Rotation Error Schematic ..................................................................................................... 61 Figure 5.14 (a) Reconstruction Result (b) Filtered Result............................................................... 63 Figure 5.15 Over-Lapping Regions of Interest ..................................................................................... 64 vii Figure 5.16 SMC Grip and Object Layout at Different Angles ......................................................... 65 Figure 5.17 Corresponding Cross-Section Due to Tilted Insertion Angle .................................. 65 Figure 5.18 (a) Sinogram (b) Threshold Sinogram............................................................................. 66 Figure 5.19 Sinograms of (a) Wave-Like (b) Straight-Edge ............................................................ 67 Figure 5.20 (a) First Object (b) Sinogram (c) Reconstruction (d) Filtered Reconstruction 68 Figure 5.21 (a) Second Object (b) Corrected Sinogram .................................................................... 69 Figure 5.22 (a) Third Object (b) Sinogram (c) Corrected Reconstruction ................................ 71 viii Chapter 1 Chapter 1 Introduction Introduction 1.1 Overview Over its more than 40-year history, computed tomography (CT) has been one of the most exciting imaging modalities developed in the medical field. In recent years, x-ray CT technology has grown tremendously. There have been significant advancements achieved in CT scanners with regard to both physical components as well as new and efficient reconstruction algorithms. Overall, there are five generations of commercial CT scanners [1]. The term tomography refers to the process of imaging a cross-section of an object from either transmitted or reflected data. It is derived from two Greek terms; tomos, meaning slice, and graphein, meaning to write. By illuminating an object from many 1 different directions, reflected or transmitted data are obtained. The ultimate aim is to employ these data to reconstruct an image of the object. Although the problem of imaging arises in many scientific fields, e.g. medicine, astronomy and non-destructive testing, its most prominent applications are found in the field of diagnostic medicine. CT is one of these applications, which has revolutionized diagnostic radiology in the twentieth century. In particular, CT has enabled medical doctors to view internal organs with unprecedented precision while maintaining patient’s safety. Sir Godfrey N. Hounsfield constructed first CT scanner in 1972 at EMI Central Research Laboratories. On the other hand, Allan M. Cormack simultaneously developed and designed a similar process. Atkinson Morley’s Hospital in Wimbledon, England, held the first CT xray machine. Hounsfield and Cormack were awarded the Nobel Prize in medicine in 1979 for their independent development in CT [2]. The solution to the problem of how to reconstruct a function from its projections was presented by Radon in 1917. However, Hounsfield’s first-invented CT scanner demonstrated that it was possible to compute high-quality cross-sectional images with a reasonable accuracy. Furthermore, Hounsfield showed that it is possible to process a large number of measurements with fairly complex mathematical operations and generate an accurate image. Since then, many advances have been made in CT scanner technology as well as in the algorithms used for reconstruction. A well-known example is the Filtered Back Projection (FBP) algorithms, developed by Ramachandran and Lakshminarayanan but popularized by Shepp and Logan. Many CT scanners were built based on these algorithms as they considerably reduced processing time for reconstruction and produced images that were numerically accurate [3]. 2 1.2 Scope and Objectives In this thesis, computed tomography algorithm is developed to enable the reconstruction of objects scanned by an imaging device that is not manufactured for CT scanning purposes. After brief review of CT history and background, its basic physical principles are discussed in Chapter 2. CT mathematical principles are also discussed with further details. In particular, Radon Transform is presented and illustrated. Chapter 3 defines the engineering problem of image reconstruction. The Fourier Slice Theorem, which relates Radon and Fourier transforms, is discussed and its proof is presented. In addition, the well-tested Filtered Back Projection method is discussed in detail with focus on its filtering approach, which is proved numerically. Furthermore, two CT sampling geometries, namely parallel-beam and fan-beam, are introduced in this chapter. However, the mathematics of the fan-beam geometry is presented with further details. A discussion of the experimental setup developed in this work is provided in Chapter 4. Both the imaging system used as well as the rotation control system developed, are explained. The proposed CT reconstruction algorithm is explained and implemented on a synthetic set of data. Objects selected to be projected and analyzed are described at the end of Chapter 4. 3 In Chapter 5, results of experimental data are displayed. Problems mechanical issues encountered while acquiring projection data as well as developing the algorithm are discussed with visual examples. Proposed solutions are also presented in details and the improvements in reconstruction results are shown in this chapter. Concluding remarks and suggested directions for future work are discussed in Chapter 6. END 4 Chapter 2 Chapter 2 Background Background German physicist Wilhelm Conrad Röntgen discovered x-rays in 1895. Röntgen was investigating the external effects from different types of vacuum tube equipment when an electrical discharge is passed through them. While experimenting with an electric discharge through a highly evacuated tube, he observed that a phosphor screen located at some distance from the tube fluoresced. He placed his hand to stop this fluorescence and was able to see a shadowed image of his bones [4]. Röntgen’s experiment may be considered the first x-ray radiography attempt. Superposition and conspicuity due to overlapping structures are among fundamental limitations of the conventional radiograph. The problem is due to the fact that a three-dimensional object is compressed to a two-dimensional image. Thus underlying structures are superimposed resulting in visibility reduction. Recognition of these limitations led to the development of conventional tomography and consequently the development of computed tomography [1]. 5 In this chapter, we briefly discuss the basic physics of computed tomography. In section 2.1, a discussion is provided to relate the physical measured projections to a mathematical model. The mathematical principles of CT are presented in section 2.2. Subsection 2.2.1 reviews image representation in image processing perspective. Radon Transform is discussed in subsection 2.2.2. Two approaches to understand Radon Transform are presented. 2.1 CT Physical Principles X-rays are electromagnetic radiations whose wavelengths range from . These wavelengths correspond to frequencies to to . X-rays are described by their intensities, i.e. number of photons. X-ray photons that leave a source reach the detector towards which they are directed in vacuum. However, if an object is placed between the source and the detector, some photons are removed from the beam; either absorbed or scattered [2]. In CT, it is essential to measure detected photons, i.e. it is required to physically measure the projections. An x-ray beam focused by using a pinhole collimator is indicated in Figure 2.1 below. The beam travels along straight lines towards the detector. In general, it is assumed that there is neither diffraction nor refraction. In addition, it is assumed that the electrons in the beam are all of the same energy, i.e. they are monochromatic or have the same wavelength [5]. 6 Figure 2.1 Physical Model (“For interpretation of the references to color in this and all other figures, the reader is referred to the electronic version of this thesis.”) Figure 0.1 Physical Model The number of photons is a Poisson random variable: { where is the number of photons at depth Poisson random variable at position distance the object } into the material and is the mean of . We can relate the photons travelling through a , the average number of photons at location (2.1) entering, and the absorption coefficient of by the following differential equations: (2.2) 7 This can be solved as: ∫ (2.3) (2.4) ∫ ( ) (2.5) ∫ Equation (2.3) is one representation of Beer-Lambert’s law, which describes intensity change as a function of x-ray energy. In this particular example, Equation (2.3) states that the number of photons at a depth x equals to the number of photons entering multiplied by some attenuation factor. While the left side of equation (2.5) is measured, the right side represents the line integrals discussed in subsection 2.2.2 [6]. 2.2 CT Mathematical Principles 2.2.1 Mathematical Background The one-dimensional continuous-time function can be represented in Frequency domain by using the Continuous Time Fourier Transform (CTFT); 8 ∫ where (2.6) is the angular frequency measured in radians per second. To represent two-dimensional images in Frequency domain, let be a function of two independent variables is defined by using and ; then its Fourier transform Continuous Space Fourier Transform (CSFT); ∫ (2.7) This two-dimensional transform can be considered as two one-dimensional transforms by splitting the exponential into two parts, with respect to ∫ In general, and respectively; ∫ (2.8) is a complex function. For example, the rectangle function can be transformed into Frequency domain as follows [3]: ∫ ∫ (2.9) 9 ∫ ∫ Figure 0.2 Two-dimensional Rectangular Function (left) and its CSFT (right) 10 2.2.2 Radon Transform The Radon transform is an integral transform that integrates a function over straight lines. According to Radon Theory, given an unknown function of an object, the line integral of the function along a straight line at some angle is called the projection or Radon Transform. It was first introduced by Austrian mathematician Johann Radon in 1917. There are many ways to explain the Radon transform. Figures 2.3 through 2.5 demonstrate a simplified explanation. In rectangular x-y plane, the standard equation for the dashed-line is: (2.10) where is the line’s slope and is the y-intercept. However, the line equation might be written in terms of parameters and { } 11 as: (2.11) Figure 0.3 Line Representation in Polar Coordinates Let the two-dimensional function Let of be defined over the region shown in Figure 2.4. be the forward projection at angle over the line , which is equivalent to the integral defined as: ∫ where (2.12) is a differential element of the line. One may express equation (2.12) in terms of and . With the help of the delta function, equation (2.12) may be written as: ∫ ∫ (2.13) 12 Hence, for any and , equation (2.13) will integrate Equation (2.13) is called the Radon transform of Figure 0.4 [7]. Line Representation in Polar Coordinates 2 Figure 2.5 shows the line integrals along Figure 0.5 over over the line forming a projection along . Line Representation in Polar Coordinates 3 13 . Another way to show Radon Transform is displayed in Figure 2.6. The geometric interpretation is a rotation by angle applied to the x-y rectangular coordinate system to produce a new r-z coordinate system, as shown in Figure 2.6. Figure 0.6 Coordinate Rotation Interpretation If we assume that the rotation is performed in counter-clockwise fashion, rotation matrix can be written as: [ ] (2.14) Therefore, the rotation to obtain the new coordinate system can be represented mathematically as: 14 * + Note that * + (2.15) is orthonormal. Equation (2.12) can be rewritten as: ∫ ( * +) ∫ ([ ] * +) ∫ ([ ]) (2.16) ∫ (2.17) From Equations (2.13) or (2.17), measured projections can be obtained. The CT reconstruction problem can now be stated as follows. If an object’s measured line integrals are provided, how do we estimate the spatial distribution of the object’s attenuation? In other words, now that we have measured the projections, how do we reconstruct the object? In the next chapter, we discuss one solution to this problem. END2 15 Chapter 3 Chapter 3 Theory Theory 3.1 Image Reconstruction Definition Image reconstruction from projections is defined as “the process of producing an image of a two-dimensional distribution (usually of some physical property) from estimates of its line integrals along a finite number of lines of known locations” [2]. Computed Tomography’s aim is to obtain information regarding the nature of the material inside an object. In general, CT reconstruction algorithms can be classified into two main categories: direct and iterative. Filtered Back Projection (FBP) is an example of a direct algorithm. Algebraic Reconstruction Techniques (ART) as well as Statistical Image Reconstruction Techniques (SIRT) are two examples of iterative algorithms [5]. 16 Many imaging systems can measure projections displacement through an object with density at angle and . This is demonstrated in Figure 3.1 below. Figure 0.1 Projection Interpretation Theoretically, projections are collected at every angle to note that where and displacement . It is essential corresponds to the center of rotation of , where . CT scanners are used to produce a sinogram, the two-dimensional array of data containing the projections of an object. From this sinogram, a two-dimensional image of the x-ray attenuation coefficient distribution in a cross-section of the object, is produced. Furthermore, the reconstruction of a series of parallel cross-sections would allow for displaying the inside shape and structure of an object [2]. Thus, the inverse 17 problem of CT, reconstructs the object, i.e. , from the measured projections, i.e. . The key for solving this inverse problem is the Fourier Slice Theorem which is discussed in section 3.2. 3.2 The Fourier Slice Theorem An important relation between Radon and Fourier transforms are presented by the Fourier Slice Theorem. The theorem relates the Fourier transform of a projection to the Fourier transform of the object along a radial. The theorem is stated by the following equations: { { } (3.1) } (3.2) then (3.3) By taking the CTFT of the projection at a particular angle frequency domain. Furthermore, if we take CSFT of a 2D object The theorem states that projections in frequency domain 18 , we obtain we obtain equals to in . since is in polar coordinates as displayed in Figure 3.2. Figure 0.2 3.2.1 Projection in Frequency Domain Proof of the Fourier Slice Theorem One straightforward proof to the Fourier Slice Theorem can be shown for the parallel projection case, i.e. . In this case, while and Equation (2.12) becomes: ∫ then, using Fourier transform, 19 (3.4) ∫ (3.5) ∫[ ∫ ] ∫ ∫ Therefore, by CSFT rotation property, this result must hold for any A more formal proof [6] can be presented by first finding the CTFT of ∫ as follows: (3.6) ∫ [∫ ∫ ∫ . ( ( By change of variables; 20 * +) * +) ] * + * + This results in: ∫ ∫ (3.7) ∫ ∫ The concepts of the Fourier Slice Theorem are summarized in a diagram displayed in Figure 3.3. 21 AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAaaaaa AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAaaaaa AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAaaaaa AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAaaaaa AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAaaaaa AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAaaaaa AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAaaaaa AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAaaaaa AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAaaaaa AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAaaaaa AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAaaaaa Figure 0.3 Fourier Slice Theorem Diagram From Figure 3.3, it is concluded that while the Fourier Slice Theorem provides a significant theoretical result, it does not provide a practical procedure for the image reconstruction problem. Practical implementations require a different approach. On the other hand, in order to implement the theorem into a practical solution, one must interpolate from these radial points to the points on a square grid as shown in Figure 3.4. This interpolation is computationally expensive because it is a nonhomogeneous interpolation. 22 Figure 0.4 Rectangular vs. Polar grids Furthermore, this calculation involves solving a large set of simultaneous equations that leads to unstable solutions. Moreover, since the density of the radial points becomes sparser as one reaches further away from the center, the interpolation error also becomes larger. Therefore, there is greater error of the high frequency components in the image, which results in degradation [3]. In section 3.3, a well-known practical procedure that overcomes these problems, Filtered Back Projection algorithm, is discussed in details. 3.3 Filtered Back Projection Algorithm Filtered Back Projection algorithm is being used in many straight ray tomography applications. It has been shown to be both accurate and amenable to fast implementation. 23 The Fourier Slice Theorem is the basis of this algorithm. Mathematically, it can be shown as follows. In order to compute the inverse CSFT of in polar coordinates, we must use the Jacobian of the polar coordinate transformation; | | (3.8) Therefore, substituting in Equation (2.7): ∫ ∫ | | ∫∫ ∫[ ∫ ⏟ ⏟ ⏞ | | ] (3.9) We can also explain Equation (3.9) in terms of convolution. We first define a new function where , 24 ∫ [ ∫| | ⏟ ] ∫| | {| | } (3.10) where {| |} ∫ (3.11) (3.12) From Equation (3.12), it is concluded that the back projection might be achieved by integrating back projections from angles to . Because of the convolution in Equation (3.10), FBP is sometimes referred to as Convolution Back Projection algorithms. The algorithms may be summarized in three steps as expressed in Equation (3.9); i. Physically measuring projections ii. Filtering the projections iii. Back projecting filtered-projections. 25 3.3.1 FBP Filtering The frequency response of filter in Equation (3.10) is given by: | | (3.13) However, real filters are band limited for some cut-off frequency as shown in Figure 3.5. Figure 0.5 Band limited High-Pass Filter Therefore, Equation (3.13) becomes: [ ( 26 ) ( )] (3.14) (3.15) Conceptually, one issue arises in the Filtered Back Projection procedure. The filtering performed here demonstrates that the zero-frequency term, i.e. the DC term, would be multiplied by zero. This implies that the reconstructed image has a zero average value. However, as shown in Chapter 5, images were reconstructed correctly. This apparent paradox might be explained by recalling an important aspect of the convolution theory. When two functions are convolved, the resulting support is equal to the sum of the supports of the individual functions as seen in Figure 3.6. This is independent of the operation domain; spatial or frequency. It is true that the entire reconstructed image, which has double the support of the original image, indeed has an average value of zero. However, the part of the image that is of interest contains exactly the correct positive values. The part of no interest that has negative values results in zero average value [8]. Figure 0.6 Convolution of Two Functions 27 3.3.2 FBP Numerical Implementation The Filtered Back Projection algorithm may also be shown numerically. Let’s consider the following object shown in Figure 3.7: Figure 0.7 Projections are performed at Numerical Example different angles and the results are found accordingly. It can be shown that this object can be reconstructed from these projections by implementing the Filtered Back Projection algorithm. 28 Figure 0.8 Filtered Back Projection Performed on Numerical Example We start from the first projection at and sum the projection at , , and respectively. This is representative of back-projecting. Furthermore, we subtract the initial projection values, i.e. , and then divide by 29 . The subtraction in this case is equivalent to high-pass filtering. This results in reconstructing the original 2X2 image. Please refer to Figure 3.8. 3.4 Sampling Geometry As mentioned in Chapter 1, there are five generations of CT scanners. Sampling the Radon Transform of the attenuation coefficients as well as the filtered back projection implementations vary in these scanners. Data collected in the first-generation and secondgeneration scanners sample a single projection from a set of parallel rays. The thirdgeneration and fourth-generation scanners sample a single projection from one focused point. The latter method of data collection is referred to as fan beam projection. Most recent scanners employ sampling geometry known as cone beam. 3.4.1 Parallel Beam As its name suggests, parallel-beam scanners have a source-detector combination that linearly scan over the projection length at each angle. After rotating by a certain angular interval, another linear scan is performed over the corresponding projection length. This type of sampling geometry is the simplest. Hence, our discussions in Chapters 2 and 3 were based on parallel beam geometry displayed in Figure 3.9. 30 Figure 0.9 Figure 0.10 Parallel Beam Geometry Fan Beam Geometry 31 3.4.2 Fan Beam Due to its importance in our experiment presented in Chapter 4, fan beam geometry is discussed in further detail. This type of geometry is possible if a single source is placed in a fixed position relative to a linear array of detectors. The fan-beam focal point is the x-ray source. This is shown in Figure 3.10. Thus, the projection process in this case is faster compared to the parallel case. Mathematical equations for this type of geometry may be derived by extending the discussion in section 3.3. In other words, the parallel-beam scenario may be converted to suit the fan-beam case. This mathematical derivation is based on [9]. First we introduce a fan beam shape in Figure 3.11 below. Recall that any object can be reconstructed mathematically using Equations (3.9) and (3.10) and substituting by : ∫ ∫ ∫ (3.16) We can write Equation (3.16) over and convert rectangular components into their corresponding polar form as follows: 32 (3.17) ∫ ∫ From basic trigonometry, In addition, the focal length Therefore, ∫ ∫ Figure 0.11 Fan Beam Projection I 33 is related to by: Equation (3.17) can be expressed in the convolution form with respect to . For a given reconstruting point , we define and . Equation (3.17) becomes: ∫ ∫ (3.18) A special property of the ramp filter is used as follows: ( ) Using the defination of the ramp filter kernel: ∫ | | ∫ | | =( ) ∫ | =( ) ∫ | | =( | ) If we donate, 34 (3.19) ( ) Then, the fan-beam convolution back-projection algorithm is obtained as: ∫ (3.20) ∫ Figure 0.12 Fan Beam Projection II END3 35 Chapter 4 Chapter 4 Experimental Setup Experimental Setup 4.1 Imaging System There are various imaging systems that can measure projections and displacement through an object of density at angle . However, it can be shown that complete and valid CT reconstruction may be performed using a vertical x-ray machine. The machine employed is General Electric Vertical X-Ray Inspection System. This inspection system is a self-contained turnkey system used to perform non-destructive testing. The system is computer controlled for repeatability and accuracy. The system has state-of the-art digital image enhancement and analysis functions, VISTAPLUS V. which offers real-time filtering of the projection data [10]. The vertical-beam X-ray inspection system includes a 150kV micro-focus x-ray source, a 3-axis manipulator that allows zooming for higher magnification and manual 36 manipulation as needed. All of these subcomponents are supplied in a compact X-ray protection cabinet with pneumatic sliding door and a hinged service door as shown in Figure 4.1. The system is designed to meet current international radiation safety standards and complies with 21 CER 1020.40 [10]. Figure 0.1 General Electric Vertical X-Ray Inspection System 37 4.2 Rotation Control System In standard CT scanning, detectors and sensors rotate while the object of interest remains stationary. However, since the imaging device used is a vertical scan, neither component is rotatable. Instead, the object itself was rotated mechanically while both the detector and sensor remained stationary. To achieve precise rotation, a Stepper Motor Controller (SMC), shown in Figure 4.2, was designed and placed inside the inspection system. Figure 0.2 Stepper Motor Controller 38 Furthermore, the SMC can be programmed to rotate by user interface software that was developed in LabView. As shown in Figure 4.3, the software provides options to specify the step size in degrees, as well as the direction of rotation, clock-wise or counter-clockwise. The software sends the correct serial string to the SMC that runs a micro-stepping motor. The motor drives a rotary table that has a three jaw chuck mounted on the face. The three jaw chuck allows for different sample diameters. Figure 0.3 X-Ray Rotation Stage Control System 39 4.3 CT Reconstruction Algorithm Based on the discussion presented in section 3.3 of this thesis, direct reconstruction algorithm was developed in MATLAB. The algorithm is based on the well-tested Filtered Back Projection method and is included in the Appendix. In the algorithm development phase, a synthetic set of data, shown in Figure 4.4, was used. -300 -200 Y-axis -100 0 100 200 300 -300 Figure 0.4 -200 -100 0 X-axis 100 200 300 Synthetic Data – Phantom The first step was to measure the projections of this Phantom. Based on these projections, a sinogram was formed. According to the FBP derivation presented in section 3.3, a High-Pass Filter, displayed in Figure 4.5, is achieved. This filter, a one-dimensional 40 function, is applied to the each row of the sinogram, a two-dimensional function, and the result before and after is shown in Figure 4.6. Clearly the projections’ edges are sharpened after filtering. Filtered projections can then be added according to the FBP algorithm. It can be observed from Figure 4.7 that back-projecting with just 46 projections of synthetic set of data produces a reliable reconstruction. 0.25 0.2 0.15 Magnitude 0.1 0.05 0 -0.05 -0.1 -0.15 -40 -30 -20 Figure 0.5 -10 0 Sample Number 10 20 High-Pass Filter Applied to Projections 41 30 40 Angle (degrees) 0 50 100 150 Angle (degrees) 0 -300 -200 -100 0 Displacement 100 200 300 -300 -200 -100 0 Displacement 100 200 300 50 100 150 Figure 0.6 Phantom Sinogram (top) and Filtered Sinogram (bottom) 42 Back Projection of Angle 0 Reconstruction using 2 Backprojections -300 -200 y coordinate 0 -200 -100 -100 -300 -100 y coordinate -300 -200 y coordinate Reconstruction using 4 Backprojections 0 100 100 100 200 200 200 300 300 300 -300 -200 -100 0 100 x coordinate Reconstruction using 8 Backprojections 200 300 0 -300 -200 -100 0 100 Reconstruction using 16 Backprojections x coordinate 200 -300 300 -100 -100 300 -300 -200 -100 200 300 -100 0 0 0 100 100 100 200 200 200 300 300 300 -300 -200 -100 0 100 coordinate Reconstructionxusing 64 Backprojections 200 300 -300 -200 -100 0 100 x using 128 Reconstructioncoordinate Backprojections 200 300 -200 -100 -200 -100 0 x coordinate 100 -300 y coordinate -300 y coordinate 200 -200 y coordinate -200 -100 0 100 x coordinate Reconstruction using 32 Backprojections -300 y coordinate -300 -200 y coordinate -300 -200 0 1 2 4 8 16 32 64 128 0 100 100 200 200 300 300 -300 -200 -100 0 x coordinate 100 Figure 0.7 200 300 -300 -200 -100 0 x coordinate 100 200 300 Phantom Reconstruction – Projection Number is indicated 43 4.4 Experimental Data In order to validate the overall system, three objects, shown in Figure 4.8, were selected to be examined and reconstructed. Both symmetric and asymmetric shapes with different material intensities were considered. Furthermore, the objective was to accomplish reconstruction using minimum number of projections. The GE Vertical X-Ray Inspection System was used to obtain x-ray scans of these objects. An example of these projection scans of each object is displayed in Figures 4.9 through 4.11. (a) Figure 0.8 (b) Objects Scanned: (a) side view and (b) top view 44 Figure 0.9 First Object 45 Figure 0.10 Second Object 46 Figure 0.11 Third Object 47 Chapter 5 Chapter 5 Results and Discussion Results and Discussion 5.1 Overview In this chapter, the algorithm was applied to experimental sets of data and analyzed. Different shapes and internal materials were considered in selecting the objects to be scanned. Sets of projection data were obtained using the proposed low cost X-ray CT system. The experimental procedure is as explained in Chapter 4. A number of issues were observed while acquiring data with this system. Observations and problems encountered are discussed in detail below. In addition, proposed solutions are presented and validated by examples. Figures 5.1 through 5.3 show the object of interest, its sinogram and reconstructed images of objects 1, 2 and 3 respectively. In each image, the sonogram representing experimental data exhibits anomalies due to mechanical problems associated with the CT system. These problems and compensation algorithms for addressing these problems are described in section 5.3. 48 5.2 Experimental Data Initial Results 5.2.1 First Object (a) (b) (c) (d) Figure 0.1 First Example (a) Projection at 00 (b) Sinogram (c) Reconstruction (d) Filtered Reconstruction 49 5.2.2 Second Object 50 100 150 200 250 300 (a) (b) 350 100 50 50200 100 250 300 300 350 350 400 400 450 600 200 250 500 150 200 400 100 150 300 450 100 Figure 0.2 200 (c) 300 400 100 200 (d) 300 400 Second Example (a) Projection at 00 (b) Out-of-Phase Sinogram (c) Reconstruction (d) Filtered Reconstruction 50 700 800 900 5.2.3 Third Object (a) 50 100 150 200 250 300 350 100 200 300 400 (b) 500 51 600 700 800 900 100 200 300 400 500 600 700 800 900 200 400 (c) 600 800 200 400 (d) 600 800 100 200 300 400 500 600 700 800 900 Figure 0.3 Third Example (a) Projection at 00 (b) Out-of-Phase Sinogram Reconstruction with (c) Rotation Error and (d) Off-Center Error 52 5.3 Problems & Issues 5.3.1 Mechanical Issues By observing sinograms in Figures 5.1 through 5.3, it can be clearly seen that the sinograms formed by the vertical system projections have some visual artifacts. In particular, the sinogram edges in Figure 5.1 (b) are not straight and are curved or wavy. In a conventional CT system, the sinogram of a cylindrical object should have straight edges. In addition, sinograms in Figures 5.2 (b) and 5.3 (b) are not constructed appropriately. Reconstruction based on these sinograms results in prone to error. As can be observed in Figures 5.1 through 5.3, the regions of interest were not reconstructed properly and exhibit a smearing effect. These issues resulting in imperfect sinograms are largely introduced by errors in the mechanical rotation stage. Mechanical rotation errors caused by the Stepper Motor Controller (SMC) can be described as follows: 1. It is observed that the SMC does not rotate a complete 360-degree in one rotation. 2. Eccentric rotation of the object because of the imperfection of SMC grip. In order to understand the effect of these problems on the sinogram, a virtual phantom object was considered as shown in Figure 5.4 and the corresponding sinogram was 53 generated using a model. Figure 5.5 depicts the ideal sinogram resulting from projections of the virtual phantom. Figure 0.4 Phantom 54 Figure 0.5 Phantom Sinogram The effects of placing the phantom off center and also rotating it with error were studied. The corresponding sinograms with error were generated by shifting the phantom down and right by 10 and 6 pixels respectively. Furthermore, it is then rotated around the center of the image, rather than the center of the cylinder, to obtain the projections. Figure 5.6 shows the resultant sinogram. Figure 5.7 shows the sinogram for the same phantom shifted down and right by 5 and 3 pixels respectively. 55 Figure 0.6 Phantom Sinogram – Shifted by 10 and 6 pixels Figure 0.7 Phantom Sinogram – Shifted by 5 and 3 pixels 56 Reconstruction result from the sinogram of shifted cylinder, shown in Figure 5.6, is computed and displayed in Figure 5.8. As can be seen, the reconstruction is accurate and shifted according to the shift of the target. Thus, the following might be concluded based on this study; when the sinogram is obtained with off-centered projections, the reconstruction is accurate in the sense that the reconstruction is also shifted accordingly. Figure 0.8 Reconstruction Result of Shifted Sinogram by 10 and 6 pixels The effect of errors in mechanical rotation on reconstruction is studied next. For this purpose, it is first assumed that the SMC rotation has a constant error. Therefore if it is desired to rotate by degrees, the actual rotation is is rotated by a full 360 degrees in degrees. Thus, when the object steps, it actually rotates more or less than 360 degrees. To explain this in an example, the reconstruction of the phantom is shown in Figure 5.9 for 57 a error. As may be observed in this figure, the region of interest, i.e. the reconstruction of objects appears smeared. Figure 5.10 shows the reconstruction result for a rotation error where the effect is more severe. Figure 0.9 Reconstruction Result with +3% Rotation Error 58 Figure 0.10 Reconstruction Result with -5% Rotation Error It is also possible that the rotational error might not be a constant and cumulative error, but a random error for each step of the SMC rotation. In the next test, for each rotation of the phantom, an error with Gaussian distribution of mean of zero and standard deviation equal to is added to the desired rotation. Figure 5.11 shows the sinogram resulting from such an error with . As the sinogram suggests, this level of error is unrealistic, but shows how the result will look like in the extreme case. Figure 5.12 shows the result of reconstructing such corrupted projections. 59 Figure 0.11 Figure 0.12 Sinogram From Rotation Error with σ=0.05 Reconstructing Result of Corrupted Projections 60 Based on these studies, it is concluded that the imperfect reconstructions shown in Figures 5.1 through 5.3 are due to a cumulative rotation error, which is not random in each step of the SMC mechanical grip. In fact, by careful examination of 360 projections with 1 degree step size, it can be seen that the SMC does not rotate back to its initial position. In particular, the first (00) and last (3600) projections are expected to be the same theoretically. However, the first project is very different compared to the last projection. This verifies that the SMC rotation mechanism has some mechanical back-lash error. While the physical improvement of the SMC components is expensive, compensation algorithms to overcome these problems was the logical alternative. To compensate for the cumulative rotation error, a correction factor was obtained and tested in MATLAB. Let us suppose that the SMC rotates by N degrees in a full rotation instead of 360-degree, as shown in Figure 5.13. Figure 0.13 Rotation Error Schematic 61 Therefore, one way to correct this error in rotation is as follows. If the projection is performed with angle , then back-projection should be performed at an angle . After extensive testing, N was found to be around 350-degree for the SMC used in NDEL. The improved reconstruction results and their original counterparts are shown in Figure 5.14 on the next page. 62 (a) (b) (c) (d) Figure 0.14 (a) Reconstruction Result (b) Filtered Result (c) Result after Mechanical Error Fixed (d) Filtered Result As may be seen in Figure 5.14 above, for each of the holes there are two circles representing the shape. Figure 5.15 offers a zoomed-in version of the reconstruction of the holes. The true hole and smeared hole are circled in red for clarity. 63 Figure 0.15 Over-Lapping Regions of Interest This phenomenon was investigated. The problem may be due to the angle of the object with the SMC grip. When an object is inserted in the SMC, the angle it forms with the grip is different from 90 degrees, as shown in Figure 5.16. Thus, when the object rotates, indications of the holes on the projections have size and shape variations as indicated in Figure 5.17. 64 Figure 0.16 Figure 0.17 SMC Grip and Object Layout at Different Angles Corresponding Cross-Section Due to Tilted Insertion Angle 65 To compensate for this problem, the algorithms were trained to correct the formed sinograms with the following steps: 1. Threshold sinogram 2. Find where { } { } 3. 4. 5. 𝑛𝑖 50 100 (a) 150 Figure 0.18 200 250 (b) 300 (a) Sinogram (b) Threshold Sinogram 66 350 400 (a) Figure 0.19 (b) Sinograms of (a) Wave-Like (b) Straight-Edge All these corrections along with corrections for the mechanical rotation error result in much more accurate reconstructions as shown in subsection 5.4. 67 5.4 Experimental Data Final Results 5.4.1 First Object (a) (c) Figure 0.20 (b) (d) (a) First Object (b) Sinogram (c) Reconstruction (d) Filtered Reconstruction 68 5.4.2 Second Object 50 100 150 200 250 300 (a) 100 (b) 350 100 200 300 400 500 600 700 50 200 100 300 150 200 400 250 300 500 350 600 400 450 700 (c) 100 800 Figure 0.21 100 200 (a) Second Object (b) Corrected Sinogram 300 400 500 600 700 800 (c) Reconstruction (d) Filtered Reconstruction 69 (d) 200 300 400 800 900 5.4.3 Third Object (a) 50 100 150 200 250 300 350 100 200 300 (b) 400 70 500 600 700 100 200 300 400 500 600 700 Figure 0.22 100 200 300 400 (c) 500 600 700 (a) Third Object (b) Sinogram (c) Corrected Reconstruction 71 Chapter 6 Chapter 6 Conclusions and Future Work Conclusions and Future Work 6.1 Conclusions  This thesis demonstrates the feasibility of converting a GE vertical-stationary x-ray imaging system into a dependable Computed Tomography imaging system.  The GE vertical-stationary x-ray imaging system was instrumented with a rotating table that holds a test sample. Multiple projections of the sample were obtained by turning the table by controlled angles, driven by external software.  Rotational and Mechanical Errors encountered resulted in artifacts in the sonogram. These issues were studied and compensation algorithms were developed for addressing different types of artifacts  The reconstruction of 2D slices using the calibrated projection data were seen to be accurate. 72  In summary, a low cost CT system was built to provide accurate reconstructions of test object in 3D. 6.2 Future Work  More extensive testing and evaluation of the algorithms needs to be done on realistic projection data. We are currently working with the Civil Engineering department to apply the reconstruction algorithm to imaging of concrete samples  Extension of the algorithm to full 3D reconstructions – This can be accomplished using two approaches, namely direct 3D reconstruction (FDK algorithm) and by slice by slice reconstruction. The implementation of FDK algorithm in the presence of issues described in this thesis will be undertaken as future work.  Fully automated system that integrates the rotation of target, projection data acquisition and reconstruction, will be developed next.  Comparison of results with that obtained using a real X-ray CT system needs to be done and evaluated quantitatively. 73 APPENDIX 74 % Name: Image Reconstruction Algorithms % Author: Ahmed Alsinan % Last modified: 11/30/2013 % clc;;clear all;close all %% Obtaining DATA ProjctionNum=360; %how many projections out of 180 for i=1:ProjctionNum rowNum=800; %Select the row number on the 2D image for which reconstruction is desired ReadImage=imread(strcat('E:\TwoNails\',num2str(i),'.bmp ')); ReadImage = ReadImage(:,1:end-13); ConvertedReadImage = im2double(ReadImage); %Convert uint8 to double [ImageRows,ImageColumns]=size(ReadImage); % Forming Sinogram SinogramMatrix(i,:)=ConvertedReadImage(rowNum,:); % Compute filter for CBP, in time domain L=size(ReadImage,1); [h,n] = AZCBPFilter(L); % To be able to apply the cut-off frequency, the filter is transfered % to freq domain, some parts is cut off (by setting to zero) and then % back to the time domain h_f = fft(h); h_f(end*1/10:end*9/10) = 0; h = real(ifft(h_f)); % Apply filter to each projection angle 75 FilteredProjectionsTmp=conv(SinogramMatrix(i,:),h); FilteredProjections(i,:)=FilteredProjectionsTmp(L+1:L+I mageColumns); end %% To compensate for the off-centric projection. Offcentric projection does not cause problem temp = SinogramMatrix>mean(SinogramMatrix(:)); temp = temp(:,1:end/2); figure;imagesc(temp) clear SinogramMatrix2 for i=1:size(temp,1) startInd(i) = find(temp(i,:)==0,1,'first'); SinogramMatrix2(i,:) = SinogramMatrix(i,(startInd(i)-80):(startInd(i)80+800)); end SinogramMatrix2 = SinogramMatrix2(:,1:end-20); SinogramMatrix = SinogramMatrix2; %Uncomment to correct for tilt figure;imagesc(SinogramMatrix2) for i=1:ProjctionNum FilteredProjectionsTmp=conv(SinogramMatrix(i,:),h); FilteredProjections(i,:)=FilteredProjectionsTmp(L+1:L+I mageColumns); end %% % % While filtering each row of the sinogram using colvolution, pixels at the % % edge will have invalid value, due to not having valid neighbors to be % % used in convolution calculations. Set these pixel's value to zero (to have better display). % FilteredProjections(:,1:20) = 0;FilteredProjections(:,end-19:end) = 0; % figure;imagesc(FilteredProjections) % %% % % Implement CBP 76 % [ImgRow,ImgColumns] = size(FilteredProjections); % tmp = zeros(ImgColumns,ImgColumns); % CBP=tmp; % result2 = zeros(size(CBP)); % for j=1:ProjctionNum % temp = repmat(FilteredProjections(j,:),ImgColumns,1); % result2 = result2 + imrotate(temp,(j1)*(360*.95/ProjctionNum),'crop'); % end % % CBP=pi*CBP/ImageRows; % figure;imagesc(result2),axis square, % figure;imagesc(result2(end/4:end*3/4,end/4:end*3/4)),ax is square % figure;imshow(result2(end/4:end*3/4,end/4:end*3/4),[.05,.3]) 77 % Name: Effect of Eccentric % Author: Ahmed Alsinan % Last modified: 11/28/2013 of1 = 0; of2 = 0; %offset of rotation center c11 = 201+of1; c12 = 201+of2; r1 = 40; c21 = 301+of1; c22 = 201+of2; r2 = 30; c31 = 201+of1; c32 = 301+of2; r3 = 20; c41 = 201+of1; c42 = 201+of2; r4 = 150; rotationErrPercentage = 5; %Accumaltive error of rotation sigma = 0;%.05; %random error of rotation/ I = zeros(401,401); for i=1:size(I,1) for j=1:size(I,1) if (((i-c11)^2+(j-c12)^2)