A NOVEL CONSTRAINT-BASED DATA FUSION SYSTEMFOR LIMITED-ANGLE COMPUTED TOMOGRAPHYbyJEFFREY EDWIN BOYDB.Sc., The University of Calgary, 1983M.Sc., The University of Calgary, 1986A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinTHE FACULTY OF GRADUATE STUDIESDepartment of Computer ScienceWe accept this thesis as conformin to the required standardTHE UNIVERSITY OF BRITISH COLUMBIAFebruary 1994© Jeffrey Edwin Boyd, 1994In presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.(Signature)__________________Department of_________________The University of British ColumbiaVancouver, CanadaDate Z- v114iCt4-DE-6 (2/88)AbstractComputed tomography (CT) is a non-destructive evaluation technique thatreconstructs the cross section of a specimen from x-ray raysum measurements. WhereasCT reconstruction is an ill-posed inverse problem that is easily solved, limited-angle CT,where raysum data are missing for a range of angles, is more severely ill-posed and moredifficult to solve. In the limited-angle case, a priori assumptions are necessary toconstrain the problem. Specimens wider than the x-ray source to sensor spacing requirelimited-angle CT. Furthermore, if the specimen is a sandwich structure, i.e., some corematerial surrounded by load-bearing face sheets, then the face sheets must lie in the nullspace. Components in the null space do not appear in the raysum data and thus confoundCT reconstruction because there is no basis for interpolation. This thesis proposes anovel constraint-based data fusion method for limited-angle CT reconstruction ofsandwich structures. The method reduces the reliance of limited-angle CT onassumptions by using range and ultrasound measurements to constrain the solution.Fusion of the data sources results in a problem with a much smaller null space that nolonger includes the face sheets. The reduction of the null space in a manner consistentwith the specimen yields a more accurate tomographic reconstruction. Synthetic and realdata experiments show marked improvement in reconstruction accuracy achieved byusing the fusion system.IITable of ContentsAbstract iiTable of Contents iiiList of Tables viList of Figures viiAcknowledgments xiiChapter 1 Introduction 11.1 Background 51.2 Limitations of Limited-Angle CT and Proposed Solution 71.3 Numerical Methods 91.4 Experimentation and Results 121.5 Outline 16Chapter 2 Background 182.1. Constraint-Based Data Fusion 192.2. Computed Tomography 232.2.1 X-ray Absorption 242.2.2. Basic Apparatus for CT Data Acquisition 252.2.3. The Radon Transform and the Fourier Slice Theorem 262.2.4. CT Reconstruction 272.3. Limited Angle Computed Tomography 302.3.1. Ill-Posed Nature of Limited Angle Computed Tomography 312.3.2. Reconstruction from Limited-Angle Data 332.4. Ultrasound 372.4.1. Forward Models 382.4.2. Ultrasound Inversion 402.4.3. Minimal Ultrasound for Fusion with ComputedTomography 422.5 Chapter Summary 43111Chapter 3 Limited-Angle Computed Tomography for Sandwich StructuresUsing Data Fusion 453.1. Null Space of the Limited Angle Radon Transform 463.2. Data Acquisition for Limited-Angle CT Data Fusion System 493.2.1 Raysum Data Acquisition 503.2.2 Range Data 513.2.3 Thickness Data 523.3 The Fusion System 533.3.1 Data Flow 543.3.2. Mathematical Formulation 553.4 Chapter Summary 58Chapter 4 Numerical Methods for Limited-Angle Tomography System 604.1 Singular Value Decomposition 614.1.1 Properties of the Decomposition 614.1.2 Least-Squares Minimum-Norm Solution 634.1.3 Application to Limited-Angle CT 644.2 Regularization and Conjugate Gradient Method 664.2.1 Regularization 664.2.2 Parameter Selection 694.2.3 The Conjugate Gradient Method 704.2.4 Application to Limited-Angle CT 724.3 Projection onto Convex Sets 734.3.1 General Description of POCS 734.3.2 Constraint Sets for Limited-Angle CT 764.3.3 Constraint and Parameter Selection 794.3.4 Application to Limited-Angle CT 804.4 Chapter Summary 81ivChapter 5 Experimentation and Results 835.1 Description of Experiments 845.1.1 Implementation of Algorithms 85Raysum Operator 85Implementation of SVD 87Implementation of R/CG 88Implementation of POCS 905.1.2 Synthetic Data 915.1.3 RealData 945.2 Results 975.2.1 SVD Trials 985.2.2 RICG Trials 1025.2.3 POCS Trials 1095.3 Chapter Summary 117Chapter 6 Conclusions and Discussion 1196.1 Conclusions 1196.1.1 Novel Contributions 1196.1.2 Generalization of Results 1216.2 Discussion 1226.2.1 Practical Application 1226.2.2 CAD Model as a Source of Constraints 1246.2.3 More Sophisticated Constraints 1256.2.4 Compton Scatter and Ultrasound Fusion 130Bibliography 133VList of TablesTable 4.1 Summary of numerical methods considered for limited-angle computed tomography with data fusion 82Table 5.1 SVD trial results 99Table 5.2 R/CG synthetic data trial results 102Table 5.3 R/CG real data trial results 103Table 5.4 POCS synthetic data trial results 108Table 5.5 POCS trial 6 (real data) results 110Table 6.1 Summary of successive minimizations used to producereconstruction in Figure 6.1 123Table 6.2 Summary of successive POCS applications used toproduce reconstruction in Figure 6.2 124viList of FiguresFigure 1.1 Synthetic cross section image for SVD trials 13Figure 1.2 Synthetic cross section image for R/CG and POCS trials.13Figure 1.3 Cross section of plexiglass phantom simulating asandwich structure with graphite/epoxy composite facesheets and aluminum honeycomb core 14Figure 1.4 Reconstructions from synthetic data with range andultrasound constraints for SVD, RJCG, and POCS 15Figure 1.5 Reconstructions of the plexiglass phantom using RJCGand POCS 16Figure 2.1 Model of computational vision 19Figure 2.2 Clark and Yuille taxonomy of data fusion systems 20Figure 2.3 Classes of weakly-coupled data fusion systems 22Figure 2.4 Strongly-coupled data fusion systems 23Figure 2.5 Schematic of x-ray absorption 25Figure 2.6 Schematic of basic computed tomography apparatus 26Figure 2.7 Example of CT showing lines of projection on a two-by-two grid 29Figure 2.8 Example cases for which only limited-angle raysums areavailable: (a) object is very long in one direction, and (b)other objects obstruct scanning 31Figure 2.9 Effect of limited-angle Radon transform on Fourierdomain sampling: (a) range of angles sampled (shadedregion) and (b) effective sampling of Fourier domain(shaded region is sampled) 32VuFigure 2.10 Triangular region of spatial support from Sato et at. [31]method for limited-angle CT reconstruction of the letter‘A’ 36Figure 2.11 Schematic of a lossless layered ultrasonic medium: (a)top layer and (b) internal layers 39Figure 3.1 Function sr,ad(x,y); shaded region has value a and zeroeverywhere else 48Figure 3.2 Sandwich specimen with face sheets at angle 0. Thewidth of the specimen prevents raysum acquisitionoutside the range— e 0 a Note that 0 <(2r/2)— 0.so the face sheets lie almost entirely within the limited-angle Radon transform null space 49Figure 3.3 Raysum data acquisition system 50Figure 3.4 Scan method for limited-angle raysum acquisition: (a)parallel with vertical and (b) at an angle to vertical 51Figure 3.5 Scanning method modified to include laser range finderfor bounding region measurement 52Figure 3.6 Apparatus for ultrasound data acquisition 53Figure 3.7 Block diagram of strongly-coupled feed forward data fusionsystem for limited-angle CT of sandwich structures 54Figure 4.1 Example of POCS convergence 75Figure 4.2 Demonstration of the variability of POCS convergence.POCS converges to three distinct but correct solutionsdepending on the order of projection operators 75viiiFigure 5.1 Steps to a compact raysum operator: (a) hypothetical gridfor reconstruction, (b) parallel rays perpendicular to uppersurface of reconstruction, (c) parallel rays at an angle toreconstruction, and (d) set of rays, one at each scan angle 85Figure 5.2 Convolution kernel for Laplacian operator 89Figure 5.3 Convolution kernel for Laplacian operator with differentcoefficients for x and y directions 89Figure 5.4 Synthetic cross section image for SVD trials 91Figure 5.5 Synthetic cross section image for RJCG and POCS trials 92Figure 5.6 Synthetic fusion data vectors: for spatial support data only(a) known region and (b) known linear attenuation, and forspatial support and face sheet data (c) known region and (d)known linear attenuation 93Figure 5.7 Cross section of plexiglass phantom simulating a sandwichstructure with graphite/epoxy composite face sheets andaluminum honeycomb core 94Figure 5.8 Apparatus for limited-angle CT raysum data acquisitionsystem 95Figure 5.9 Reconstructed images for SVD trials: (a) trial 1reconstruction (max = 0.48, mm = -0.01), (b) trial 1 errorvector (max = 0.30, mm = 0.0), (c) trial 2 reconstruction(max = 0.55, mm = 0.0), (d) trial 2 error vector (max = 0.23,mm = 0.0), (e) trial 3 reconstruction (max = 0.42, mm =-0.04), (1) trial 3 error vector (max = 0.04, mm = 0.0), and (g)rescaled trial 3 error vector (max = 0.04, mm = 0.0) 100ixFigure 5.10 (a) Weakly-coupled reconstruction using SVD (max = 0.48,mm = -0.01), and (b) the corresponding absolute error vector(max = 0.15, mm = 0.0) 102Figure 5.11 Reconstructed images for RICO synthetic data trials: (a) trial1 reconstruction (max = 0.49, mm = -0.06), (b) trial 1 errorvector (max = 0.40, mm = 0.0), (c) trial 2 reconstruction(max = 0.55, mm = -0.05), (d) trial 2 error vector (max =0.43, mm = 0.0), (e) trial 3 reconstruction (max = 0.45, mm =-0.07), and (1) trial 3 error vector 103Figure 5.12 (a) Weakly-coupled reconstruction using RICO (max =0.49, mm = -0.05), and (b) the corresponding absoluteerror vector (max = 0.22, mm = 0.0) 105Figure 5.13 Reconstructed images for RICO real data trials: (a)a2=0.001,(b)a2=0.01,(c)a2=0.1,(d)a2=1.0,(e) a2 =10.0,(f) =100.0,and(g) a12 =0.001 anda2 = 0.1. See Table 5.3 for maxima and minima 107Figure 5.14 Reconstructed images for POCS synthetic data trials: (a) trial1 reconstruction, (b) trial 1 error vector (max = 0.4, mm =0.0), (c) trial 2 reconstruction, (d) trial 2 error vector (max =0.4, mm = 0.0), (e) trial 3 reconstruction, (f) trial 3 errorvector (max = 0.22, mm = 0.0), (g) trial 4 reconstruction, and(h) trial 4 error vector (max = 0.25, mm = 0.0) 110Figure 5.15 (a) Weakly-coupled reconstruction using POCS, and (b)the corresponding absolute error vector (max = 0.24, mm=0.0) 113Figure 5.16 Reconstructed cross section for POCS real data trial 5 114xFigure 5.17 Reconstructed images from POCS real data trial 6: (a)ER=0.1,and(b) ER—0.001 114Figure 5.18 Plots of convergence for POCS: (a) trial 6(a)- ER = 0.1and(b)trial6(b)- ER=0.001 116Figure 6.1 Reconstruction from successive RICG iterationsincorporating the assumption that pixels are either air(u = Ocm’) or plexiglass (p = 0.45cm1) 127Figure 6.2 Reconstruction from successive POCS iterationsincorporating the assumption that pixels are either air(p = Ocm’) or plexiglass (p =0. 45c&’ ) 129xiAcknowledgmentsMany people have contributed to the success of my research. I owe Jim Little, mysupervisor, a special debt of gratitude. The time he committed to me and his genuineinterest in my work helped me greatly. I am specially grateful for the opportunity hegave me to prove myself in computer science. The other members of my committee,Bob Woodham, Un Ascher, David Lowe and Jean Meloche, have also earned mygratitude both as instructors in graduate courses and in guiding my thesis research.UBC is a special place. It is not special because of the buildings or the scenery,but because of the people who live and work there. I take this opportunity to thank myfellow students for their friendship, conversation, and companionship that made my stayat UBC worthwhile.I thank the Department of National Defence for its fmancial support. At DefenceResearch Establishment Pacific, several people lent me a hand. The assistance of BillRamsbottom, of Galatea Research Inc., with data acquisition was invaluable. BillSturrock contributed his wisdom, guidance and encouragement. Ken McRae made surethat I had the best available advice with respect to ultrasound technology.During the course of my research I had the occasion to meet Heinz Bauschke andhis colleagues in the Centre for Experimental and Constructive Mathematics at SimonFraser University. Heinz’s help in understanding convex projections was not onlyextremely helpful to my project, but also inspired me to learn more about appliedmathematics.It is a long road to complete a PhD program, so I always strove to enjoy what Iwas doing. Fortunately for me there is my wife, Brenda. She not only supports myacademic endeavours, but is ever vigilant in assuring that my life is full and happy, bothat work and at home.x1Chapter 1IntroductionNon-destructive evaluation (NDE) is the field of endeavour that evaluates thestructural integrity of an object without destroying the object. Measurements of aplethora of physical phenomena, including x-ray radiation, neutron radiation, ultrasound,eddy currents, infrared radiation, and fluorescence, indicate the structural integrity ofNDE specimens. Although often obscured by detail, the ultimate goal of an NDEtechnique is to use such measurements to predict the failure of a specimen and preventexpensive or catastrophic failure.Radiographs, images produced by projecting x-rays through a specimen onto asensor, are commonplace in NDE. Each pixel in a radiograph is a raysum measurement,i.e., the integral of linear attenuation coefficients of the specimen along the ray from thex-ray source to the sensor. A single radiograph presents the three-dimensional structureof a specimen projected onto a two-dimensional image. Proper inspection requiresimages from multiple viewing angles. For example, in real-time radiography, viewingimages over continually changing angles hints at the underlying three-dimensionalstructure.2Computed tomography (CT) is an NDE technique that reconstructs a cross sectionof a specimen from its raysums. Adjacent cross sections combine to form a volumetricimage, so CT solves the problem of determining three-dimensional structure fromraysums. The proliferation of CT scanners in medicine, in spite of cost and patientexposure to radiation, attests to the diagnostic value of reconstructed cross sections.The Radon transform of a cross section is the set of all possible raysums throughthe cross section. CT reconstruction from the continuous Radon transform is an ill-posedinverse problem; the solution does not vary continuously with the data. The practicalproblem of reconstruction from a discrete sampling of the Radon transform is also ill-posed, but the ill-posedness is overcome by assuming that the reconstruction is bandlimited.Often it is not possible to measure raysums over the full range of angles. In suchcases, reconstruction is more illposed* because of the missing data. Techniques that dealwith this type of limited data are called collectively limited-angle CT and have receivedconsiderable attention. To compensate for missing data, limited-angle CT methods relyon constraints to compute reconstructions.Specimens wider than the x-ray source-to-sensor spacing prevent raysummeasurement at some angles and require limited-angle CT for reconstruction. Normallythe specimen is part of some mechanical structure and must be rigid in the presence ofbending moments. For economy of material and weight the structures use a thin load-bearing material at the outer surfaces with some core between, i.e., a sandwich structure.Common examples are I-beams, windsurfers, skis, and aircraft control surfaces. Thisthesis focuses on aircraft control surfaces that use two face sheets (typically aluminum,but some modem aircraft use graphite/epoxy composite) to carry loads on either side of* We say that a problem is more ill-posed than another if it has a larger null space, i.e., the null space has agreater number of dimensions.3an aluminum honeycomb core. The core provides a rigid separation between the facesheets.Not only do wide specimens restrict raysum acquisition, but, as this thesis shows,limited-angle CT cannot accurately reconstruct the cross section in the presence of facesheets. The face sheet structure lies almost entirely within the null space of the limited-angle Radon transform. In short:1. wide specimens require limited-angle CT,2. wide specimens are often sandwich structures with face sheets, and3. limited-angle CT does not work in the presence of face sheets.Typical constraints that make limited-angle CT possible cannot account for the facesheets but instead lead to erroneous reconstructions based on general assumptions.However, the reliance of limited-angle CT on constraints makes it an excellent candidatefor constraint-based data fusion. Constraint-based data fusion uses both a prioriassumptions and data from other sources to constrain the solution of an inverse problem.Ultrasound offers a source of additional data that complements the raysums.Whereas raysums are oblivious to discontinuities along the rays (the very property thatprevents them from dealing with the face sheets), ultrasound responds precisely to suchdiscontinuities and is capable of measuring the face sheet structure. This thesis proposesa novel method for limited-angle CT using data fusion. The method capitalizes on thecomplementary nature of x-rays and ultrasound and is specifically applicable to sandwichstructures. In the context of this thesis we consider the method as a means of producingaccurate reconstructions only. For the sake of development, we ignore the ultimateapplication, inspection to predict failure.The method fuses laser range data and ultrasound data with raysums to computean accurate reconstruction. Laser range data give the spatial support of the specimen.Ultrasound data, combined with range data, segment the reconstruction into exterior, facesheet, and interior regions. These data then restrict the null space of the forward operator4so that it does not include the face sheet structures, leading to accurate reconstructions.Data fusion alone cannot guarantee better solutions to inverse problems, but here thecomplementary nature of ultrasound and x-ray yields superior reconstructions. Thisthesis shows that the proposed limited-angle CT system computes accuratereconstructions of sandwich structures where an accurate reconstruction is not possibleotherwise.The following system of linear equations combines the raysum data with thefusion data:ER1 E 1[WJX [wxFj’where R is the forward raysum operator, x is the cross section image vector, W is adiagonal matrix, y is the raysum data and XF is a vector of linear attenuation coefficientsbased on the fusion data. An image vector of weights indicating whether or not the linearattenuation coefficient for a pixel is known from the fusion data, forms the diagonal of W.This thesis considers three methods to solve for x from the above equation:1. singular value decomposition (SVD),2. regularization and conjugate gradient method (RICG), and3. projection onto convex sets (POCS).Experimental trials with synthetic and real data apply all three methods to computereconstructions. Synthetic data trials show marked improvement in the reconstructionsobtained using the fusion data constraints. In addition, the SVD trials show that thenumber of dimensions in the null space (as indicated by the number of singular valuesequal to zero) decreases by adding the fusion constraints. Real data trials show that theproposed fusion method works with real data, and, in general, experimental resultssupport the contention that the proposed technique improves reconstruction.SVD gives accurate reconstructions from synthetic data, but, because ofcomputational limitations, it is only applicable to small problems (smaller than practicalreal data problems). RJCG gives excellent reconstructions from error free synthetic data5and is fast (allowing application to real data problems). However, in the presence oferrors in the real data, regularization forces a tradeoff between image smoothness andfidelity to data, and no degree of compromise is totally acceptable. POCS yields the bestresults of the three numerical methods considered. It produces accurate reconstructionsfrom error-free synthetic data, does not force a smoothness compromise in the presenceof errors in real data, and converges quickly. Each of the numerical methods usesdifferent constraints to arrive at a solution, which means the methods solve differentproblems. The results presented here serve only to demonstrate alternative methods forreconstruction and do not comprise a proper comparison of the methods.This thesis concludes that the complementary nature of x-ray and ultrasoundallows exploitation of constraint-based data fusion to improve limited-angle CT. Thiscomplementary nature means the constraints from the ultrasound data are not redundantbut instead reduce the null space of the problem. In particular, limited-angle CT cannotdeal with face sheets in sandwich structures. However, the ultrasound can measure theface sheets of the structure thus allowing accurate reconstruction when combined with theraysum data.1.1 BackgroundThe proposed limited-angle CT method uses constraint-based data fusion asdescribed by Clark and Yuille [9]. They describe two types of data fusion systems forsolving inverse problems. The systems are composed of modules where each moduletakes some sensory data as input and produces as output a solution or set of solutions toan inverse problem. The modules may or may not use a priori constraints.In the first type of fusion system, the modules operate independently to producesolutions. A fusion module then combines these solutions to form a single solution. Theindependent operation of the modules characterizes these weakly-coupled systems. In thesecond type of fusion system, at least one module provides constraints to another.Interaction of modules distinguishes these strongly-coupled systems. Implementation of6the interaction via constraints to modules gives the alternative appellation constraint-based.It is common practice to use a priori assumptions to constrain a module solvingan inverse problem. The constraints serve to stabilize the solution to an ill-posed or ill-conditioned inverse problem (e.g., Poggio, Torre and Koch [30] and Bertero, Poggio andTone [5]). However, the assumptions may be insufficient or invalid, thereby leading toerroneous solutions. In such cases, constraint-based data fusion is useful because it shiftsreliance from poor assumptions to measurements from other sensory sources. Althoughin general the use of other sensory data to provide constraints does not necessarilyimprove the solutions of inverse problems, if the data sources complement each other,i.e., one provides information where the other cannot, then improved results are possible.Computed tomography (CT) is an inverse problem in the field of non-destructiveevaluation (NDE). CT computes a reconstruction of the cross section of a specimen fromits raysums (line integrals measured by an x-ray system). The cross section is an imagein which the pixel values are the linear attenuation coefficients of the specimen. Thecomplete set of raysums for a cross section is its Radon transform. The Fourier slicetheorem relates the Radon transform to the Fourier transform. It states that measuringraysums is equivalent to sampling in the Fourier domain. Inversion of the continuousRadon transform is an ill-posed problem; the solution does not vary continuously with thedata. In practice, only measurement of the discrete Radon transform is possible.Reconstruction from the discrete Radon transform is also ill-posed since the null space isnon-trivial and consists of high-frequency components. Requiring that the reconstructionbe band-limited easily handles the ill-posedness.Circumstances arise that prevent raysum measurement throughout a full range ofangles leading to a special variation of the CT problem called limited-angle CT. Limitedangle CT is ill-posed because the missing raysum data leave the Fourier space undersampled. Ill-posedness in the limited-angle case is more severe than for conventional CT7(i.e., the null space is much larger for limited-angle CT) and it is not easily handled. Avariety of techniques proposed in the literature perform limited-angle CT. Withoutexception, they all rely on a priori assumptions to constrain the reconstruction and find aunique solution. In cases where a priori assumptions may be erroneous or insufficient,limited-angle CT can benefit from constraint-based data fusion by incorporating datafrom a complementary source.Ultrasound is a source of data that complements x-rays very well. X-rays measureline integrals and are insensitive to discontinuities along the line of integration. Theopposite is true for ultrasound which detects discontinuities in a layered system. Thelimited-angle CT system with data fusion proposed in this thesis exploits this property ofultrasound to measure the thickness of face sheets in sandwich structures. X-rays cannotmeasure this thickness, so ultrasound, by virtue of its complementary nature to x-rays,adds non-redundant data to the problem.In summary, CT is an ill-posed inverse problem that relies on band-limitingconstraints for solution. Limited-angle CT is further ill-posed because of incompletesampling of the Fourier domain and requires additional constraints for solution.Situations may arise where a priori assumptions used as constraints for limited-angle CTare either erroneous or insufficient. Reliance on constraints for solution suggests the useof constraint-based data fusion. Data fusion requires another source of data tocomplement the raysum data. Ultrasound provides such a source because it can measurediscontinuities that x-rays cannot.1.2 Limitations of Limited-Angle CT and Proposed SolutionThe null space of the limited-angle Radon transform is the set of all functionswith Fourier transforms equal to zero within the region of the Fourier domain sampled bythe limited-angle Radon transform. Physical structures that lie within this null spacecannot be reliably reconstructed from the limited-angle raysums. Such structures containlong edges at angles beyond those sampled by the raysums.8Sandwich structures are commonplace and lie within the null space of the limited-angle Radon transform. They economize on weight and material by using some thinmaterial at the outer surfaces of the structure with some core material sandwichedbetween. The outer surfaces carry the loads due to bending while the core servesprimarily to separate the two surfaces. Commonplace examples are I-beams, windsurfers,skis, and aircraft control surfaces. In aircraft control surfaces, two face sheets (typicallyaluminum but some modern aircraft use graphite/epoxy composite) surround analuminum honeycomb core. Face sheets carry most of the loads in the specimen whilethe honeycomb provides a rigid separation of the face sheets. An x-ray source-to-sensorspacing narrower than the face sheet width necessitates limited-angle CT. At the sametime, the face sheet that necessitates limited-angle methods lies in the null space of thelimited-angle Radon transform. Thus, the presence of face sheets in a sandwich structureconfounds accurate CT reconstruction.This thesis proposes a constraint-based data fusion system that deals with theinability of other limited-angle CT systems to reconstruct sandwich structures accurately.The system uses constraints derived from laser range data and ultrasound data toconstrain the tomographic reconstruction. Here we consider the system as a means ofproducing accurate reconstructions only. Further evaluation as a technique for inspectionand prediction of failure is ignored.There are three modules in the system. The first module takes laser range data asinput and computes a bounding box for the specimen, i.e., the spatial support of thespecimen. The output of the module is a segmentation of the reconstruction into exteriorand interior regions. Since the exterior regions must be air (which has zero linearattenuation) the segmentation provides a constraint for reconstruction.The second module takes ultrasound thickness measurements as input, and,constrained by the spatial support data, further segments the reconstruction into exterior,9interior and face sheet regions. With the known linear attenuation of the face sheetmaterial, the segmentation provides another constraint for reconstruction.The last module takes limited-angle raysum data as input and, subject to spatialsupport and face sheet constraints, reconstructs a cross-section of the image as output.This module computes its output by solving the linear system:ER] E 1[WJX [wj’where R is the discrete limited-angle raysum operator, x is the cross section image, and yis the raysum data. The segmentation of the reconstruction provides a vector of weightsthat forms the diagonal of the diagonal matrix W. A weight of one on the diagonalindicates that the corresponding pixel has a known value while a zero indicates unknown.XF is a vector of linear attenuation coefficients derived from the constraints. Entries inXF that correspond to a weight of one in W have known linear attenuation and otherentries are arbitrary.rRl r 1Let A= [ j and b [ j to give a generalized version of the above equation:W WXFAx = b.A is not square and is not full column rank. Solution for x from the above is still an ill-posed problem; if a solution exists then there is an infinity of solutions. Some a prioriassumptions are necessary to arrive at a unique solution. However, these assumptions donot lead to an erroneous reconstruction because A properly accounts for the face sheetsdue to the fusion constraints.1.3 Numerical MethodsThis thesis considers three numerical methods for solving the linear system ofequations for the proposed limited-angle CT method. These are:1. singular value decomposition (SVD),2. regularization and conjugate gradient method (RJCG), and103. projection onto convex sets (POCS).This section introduces the methods which are explained in more detail in Chapter 4.Table 4.1 summarizes the comparison of all three methods.The first method, SVD, computes the decomposition:A=UEVTwhere U, and V are orthogonal and I is diagonal. Given the decomposition, it ispossible to condition the matrix A and solve for the minimum-norm least squaressolution, , using the pseudo-inverse of I I i.e.:= VrUb.SVD has the beneficial property of determining the null space and range of A, and itssingular values (the square roots of the eigenvalues of ATA). The number of non-zerosingular values gives the rank (row and column) of ATA. Unfortunately, SVD hasO(mn2+ n3) complexity in time (where m is the number of rows in A, and n is thenumber of columns in A and the number of pixels in the reconstruction) and is too slow toapply to large problems.The second method, R/CG, combines regularization and the conjugate gradientmethod. The conjugate gradient method (CG) solves the linear system Ax = b byminimizing the objective function:E=ixTAx_xTbFc,2where c is a constant and A is positive definite. For a sparse matrix A, iterations are fastand, depending on the condition number of A, CG converges quickly to a solution. A forthe data fusion system is neither positive definite nor square. Regularization solves thisproblem by creating the new objective function:EA 1 lb 12E= [aQj[a2x0j11rAlwhere 0 makes I full column rank, a is a constant and x0 is some vector to be[aQjspecified. CG finds x that minimizes E by solving the new positive definite system ofequations:(AT+a2T0)x = ATb +a2QTx0.a2 selects a compromise between the original system of equations and the solution toOx = .2x0 With a2 set too small the system is ill-conditioned and CG fails. A largera2 gives a system that is well-conditioned but biases the solution towards Ox = 2x0.There is a variety of options in selection of 0 and x0. In experimental trialsconducted for this thesis 0 is:ED0=1 1[D2where D1 and D2 are discrete first derivative operators in the x and y directionsrespectively, and:roxo=LoThus, a2 sets the degree of compromise between a solution that conforms to the originalill-conditioned and ill-posed problem Ax = b, and a well-conditioned problem with thesolution biased to be flat.The third numerical method, POCS, uses a set of constraints where eachconstraint requires the solution to lie within a convex set. POCS also requires that theintersection of the constraint sets be non-empty. Starting from some arbitrary initialsolution, POCS iterates towards an ultimate solution by sequentially projecting thesolution onto each of the constraint sets. A non-empty intersection of the sets guaranteesconvergence. POCS subsumes the better known technique called ART (algebraicreconstruction technique). Whereas POCS applies to any convex constraint sets, ART isPOCS using only convex sets that are hyperplanes. For the limited-angle CT problem athand the convex sets are hyperslabs and spheres. In practice, with 0 as a starting point,12POCS finds a solution near the minimum norm solution. This is not guaranteed ingeneral for POCS. POCS converges quickly for the problem of interest. So POCS is afast method (comparable to CG for this problem), finds a solution like the minimum normsolution of SVD, and avoids the smoothing of RJCG. Also, the flexibility of POCSallows mixing of 12 and fits to the data, which is valuable in constraint-based datafusion.Implementation of POCS for the proposed limited-angle CT data fusion systemuses three constraint types:1. conformity to raysum data,2. conformity to fusion data, and3. conformity to amplitude limits.Previous work implements type 1 constraints, conformity to raysum data, as linearequalities, i.e., the convex set is a hyperplane within which the solution must lie. Such aconstraint restricts the solution severely by requiring it to match the data precisely. Inthis thesis, the type 1 constraints are hypersiabs surrounding a hyperplane, i.e., theintersection of two half spaces centred about the hyperplane. Solutions in the hyper planematch the data exactly while other solutions in the slab only approximately match thedata. A constant parameter determines the thickness of the hyperslabs. Type 2constraints, conformity to fusion data, are a variation on previously published work usingfull reference images for reconstruction. In this thesis, the constraint does not require afull reference image, but instead uses a partial reference image from the fusion data.Type 3 constraints, conformity to amplitude limits, are also hyperslabs requiring thatvalues in the solution lie between fixed limits.1.4 Experimentation and ResultsA series of experimental trials verifies the validity of the proposed limited-angleCT data fusion system. Each trial consists of application of one of the three numerical13algorithms to some data, to produce a reconstruction for analysis. The algorithm, itsparameters, and the data uniquely define each trial.Experiments use both synthetic and real data. Application of the forward raysumoperator to a synthetic cross section image generates synthetic raysum data. Thesynthetic cross sections imitate the plexiglass phantom used in the real data trials, whichitself imitates a graphite/epoxy sandwich structure with honeycomb core. The plexiglassphantom has thicker parts than aluminum honeycomb to avoid pushing resolution limitsof the apparatus. Because SVD can handle only small problems, SVD trials are limited tosynthetic data based on the low-resolution cross section of Figure 1.1. Synthetic data forRICG and POCS trials are based on the high-resolution image of Figure 1.2.Figure 1.2: Synthetic cross section image for R/CG and POCStrials.The plexiglass phantom used as the specimen for real data trials, shown in Figure1.3, consists of two plexiglass face sheets (3mm thick) surrounding a core of plexiglassmembers (3mm thick) perpendicular to the face sheets. Some plexiglass inserts providean anomaly in the structure similar to entrapped water. Chapter 5 describes the dataacquisition system for the real data trials. It is a modification of an x-ray inspectionFigure 1.1: Synthetic cross section image for SVD trials.14system intended to produce two-dimensional radiographs and is not ideally suited to CT.Nevertheless, the data it produces are adequate to show the validity of the method.HI 12 j...m. All dimensions in mmTTop and bottom face sheets (horizontal components)El Simulated honeycomb (vertical compontents)Simulated entrapped water (defect)Figure 1.3: Cross section of plexiglass phantom simulating asandwich structure with graphite/epoxy composite face sheetsand aluminum honeycomb core.For each numerical method, the first three trials use synthetic data with varyinglevels of fusion constraints. All methods exhibit a marked improvement in the accuracyof reconstruction when using spatial support and face sheet constraints. For SVD theerror measure improves from 50.7% to 4.8% by using all constraints. Similarly, withRICG the error measure drops from 64.0% to 6.7%, and with POCS the error measuredrops from 62.6% to 6.0%. Furthermore, SVD trials show the change in the number ofsingularities in the problem. There is a reduction from 148 singularities of 300 singularvalues for raysum data to only 13 singularities for raysum data with spatial support andface sheet constraints. The character of the problem changes significantly with theconstraints in such a way as to improve the accuracy of the reconstruction. Figure 1.4shows sample reconstructions for all three numerical methods using all the constraints.Real data trials with RJCG and POCS show the validity of the proposed systemwith real data. Although the quality of reconstruction is not on par with the synthetic15trials, this is largely due to errors in raysum measurements. Superior apparatus shouldlead to reconstructions of sufficient quality for practical application. Figure 1.5 showsreconstructions of the plexiglass phantom for RJCG and POCS.(c)Figure 1.4: Reconstructions from synthetic data with rangeand ultrasound constraints for (a) SVD, (b) RICG, and (c)POCS.Errors in the real data emphasize the effects of parameters used in reconstruction.A series of real data trials for both R/CG and POCS highlights the effects of theparameters. RICG trials show the compromise selected by the parameter a2 and POCStrials show the effect of changing constraint set size on convergence and reconstructionquality.(a)(b)16Figure 1.5: Reconstructions of the plexiglass phantom using(a) R?CG and (b) POCS.1.5 OutlineChapter 2 gives background material for the thesis starting with a review of theClark and Yuille model for constraint-based data fusion. Following that, the chaptergives an introduction to computed tomography, pointing out that CT is an ill-posedproblem that is easily managed. From CT follows the idea of limited-angle CT wherereconstruction is done from incomplete data. Limited-angle CT is ill-posed to a greaterdegree and must rely on a priori assumptions for constraints, thus making it a goodcandidate for improvement by constraint-based data fusion. Chapter 2 describesultrasound, which complements x-rays well, and, therefore, provides an excellent sourceof additional data for fusion.Chapter 3 defines the null space of the limited-angle Radon transform and showsthat face sheets in a sandwich structure lie mostly within this null space. The result isthat limited-angle CT cannot reconstruct sandwich structures properly. A novel methodfor limited-angle CT using constraint-based data fusion proposed in the chapter addresses(a)(b)17this problem and is the focus of this thesis. Chapter 3 describes the method and its datasources, ending with the mathematical formulation of the method as a system of linearequations.There are many ways to solve a system of linear equations. Chapter 4 describesthe three numerical methods used in experimentation for this thesis, including applicationof each method to the proposed limited-angle CT method. A comparison at the end ofChapter 4 highlights the strengths and weaknesses of each numerical method.Chapter 5 presents the results of experimentation that verify the validity of theproposed method. Synthetic data trials show a reduction in null space and a greatimprovement in accuracy of reconstruction with the incorporation of data fusionconstraints. Real data trials show that the method works in practice and not just witherror-free synthetic data. Experiments also show the effects of parameters onreconstruction quality.Chapter 6 summarizes the conclusions of the thesis. It also discusses possibilitiesfor further enhancements and variations of the proposed method, and its potential forapplication.18Chapter 2BackgroundAs a prelude to the proposal of a novel method for limited-angle CT in Chapter 3,this chapter presents relevant background material in four areas: constraint-based datafusion, CT, limited-angle CT, and ultrasound. Section 2.1 describes the constraint-basedapproach to data fusion. Constraint-based data fusion, so called because it uses theoutput from some inversion module as constraints for another, reduces or eliminatesreliance on a priori assumptions for constraints. The value of the fusion emerges whenthe fusion constraints replace invalid assumptions.CT is the ill-posed problem of reconstructing the cross section of an object fromits x-ray raysums. Although the problem is ill-posed, the null space is small, consistingof high-frequency oscillations. Assuming a priori that the reconstructed cross section isband-limited is a sufficient constraint for reconstruction. Section 2.2 reviewsconventional CT.The shape of an object, or its environment, may prevent acquisition of raysumdata through the full range of angles required for CT. The ensuing problem of limitedangle CT (reviewed in Section 2.3) is ill-posed too, but has a much larger null space. A19priori assumptions provide constraints beyond those of conventional CT. The reliance oflimited-angle CT on assumptions that may not always hold suggests that it can benefitfrom constraint-based data fusion.Ultrasound (reviewed in Section 2.4) is a versatile NDE technique providing datathat are complementary to x-ray raysums. Whereas x-rays are insensitive todiscontinuities along the path of radiation, ultrasound explicitly detects discontinuities inits path. It is this complementary nature that makes ultrasound potentially useful as asource of data for fusion with limited-angle CT.2.1. Constraint-Based Data FusionAs data fusion is a central part of this thesis, it is necessary to consider how toimplement the data fusion. This thesis adopts the elegant constraint-based approach ofClark and Yuille [9], summarized in this section.Representation of WorldI I I II IData Acquisition, Computational Vision,Measurement, Imaging InversionFigure 2.1: Model of computational vision.Figure 2.1 shows a model of computational vision systems. On the left-hand sideof the figure is a world that is the setting for a task. Some data acquisition processacquires measurements of the world. These measurements (or data) are images for visualsystems, but they may be raysum data or ultrasound traces for NDE systems. The goal ofthe computational vision system (right-hand side of Figure 2.1) is to invert the acquiredAcquired Data,World Measurements,Images20data to a representation of the world. Thus, computational vision is the solution ofinverse problems where the solution is a representation of the world, suitable for the taskat hand.Invariably, the inversion problem is ill-posed because it has no unique solution.One approach to dealing with this ill-posedness is to create a new problem by applyingconstraints to the original. The constraints restrict the set of solutions so that themodified problem is well-posed. Ill-posed problems occur often in computer vision. Forexamples see Poggio, Torre and Koch [30] or Bertero, Poggio and Torre [5].When there is only one source of data, a priori assumptions about the nature ofthe world, and consequently about the nature of the solution, provide constraints. Theassumptions are a convenient way of restricting solutions in the face of ambiguous dataand may not always hold. Where assumptions are not valid, data fusion offers analternative. Constraint-based data fusion systems acquire additional data fromindependent sources and use the resulting solutions to constrain the original inverseproblem. The notion of fusing data from multiple sources as constraints is the essence ofconstraint-based data fusion. It strives to improve the quality of solutions to inverseproblems by using independent data sources as constraints, without falling prey toerroneous a priori assumptions.FeedClass I Class II Class III Recurrent________Forward__________Figure 2.2: Clark and YuilIe taxonomy of data fusion systems.21Clark and Yuile give a taxonomy of data fusion systems (shown in Figure 2.2).Fusion systems consist of a set of modules where each module accepts some input dataand performs an inversion subject to some constraints. The first division in the taxonomyis between weakly-coupled and strongly-coupled systems. Weakly-coupled systems aredistinguished by their independent modules, i.e., there is no sharing of informationbetween modules. Strongly-coupled systems are the opposite; for a system to bestrongly-coupled at least one module must constrain the output of another.There are three classes of weakly-coupled systems, all shown in Figure 2.3.Modules in a class I system are stable and produce unique solutions independently. Theultimate solution produced by the system is a weighted sum of the solutions from theindividual modules, with the weights determined by the relative reliability of themodules. Solutions from more reliable modules are weighted greater than those from lessreliable ones. Class II weakly-coupled fusion modules do not produce a unique solution.It is the function of the fusion module to combine the module outputs into a uniquesolution (the fusion module may use a priori constraints to help it). An example is themeasurement of electrical resistance by separate voltage and current measurements.Neither measurement yields a unique resistance value, but the combination of the twoVdoes. In this case the fusion module performs the calculation R =—i. Class III fusioncombines classes I and II. Modules in a class III system do not produce unique solutionsso a fusion module combines the module outputs as in a class II system, while solutionsfrom the modules are not equally reliable so they are weighted as in class I.Figure 2.4 shows schematically feed-forward and recurrent strongly-coupledfusion systems. Their common feature is that at least one module constrains the output ofanother module. Feedforward is the simplest of the strongly-coupled systems. Thesystem feeds forward the output of one module to provide constraints to a second module.Output from the second module is the output of the system. In the recurrent system the22output from each module feeds back to other modules in the form of constraints. Therecurrent system is the most complex of the fusion methods.Sensor dataModule 1a1 7% w1aw1a+ w2a2Sensor dataClass IClass IIf(a, a)f(w1a,22,...,wa)a priori constraintsSensor dataa priori constraintsClass IIIFigure 2.3: Classes of weakly-coupled data fusion systems.23Sensor data Sensor dataoutput 1output output 2Feed-forward RecurrentFigure 2.4: Strongly-coupled data fusion systems.Strongly-coupled fusion is synonymous with constraint-based data fusion.Whenever an inversion module relies on a priori assumptions to constrain the solution ofan inverse problem, constraint-based data fusion can reduce reliance on assumptions. Itdoes so by replacing constraining assumptions with constraints derived from other datasources. There is, of course, no guarantee that by using fusion the inverse problem willbecome well-conditioned or even well-posed. However, with careful attention to thenature of the data and its use as constraints, it is possible to improve the quality ofsolution by avoiding reliance on possibly erroneous a priori assumptions.2.2. Computed TomographyAs x-ray radiation passes through an object the object attenuates it.Measurements of the attenuation approximate raysums of the object. Computedtomography (CT) is the inverse problem of reconstructing a cross section of an objectfrom a set of coplanar raysums. The solution, i.e., the cross section, is a distribution of xray linear attenuation coefficients.This section presents the fundamentals of CT starting with a description of x-rayabsorption and raysum measurement. The basic apparatus for CT illustrates the use of xrays to acquire raysum data for reconstruction. These data are a discrete measurement of24the Radon transform of the object, the basis for CT reconstruction. Inversion of thecontinuous Radon transform is an ill-posed problem because the solution does not varycontinuously with the data. In the case of discrete data, reconstruction is ill-posedbecause there is not a unique solution. In both cases the ill-posedness is easily handledby assuming that the solution must be band-limited.2.2.1. X-ray AbsorptionFigure 2.5 shows schematically x-ray radiation penetrating an object along a ray.The object absorbs some of the radiation, attenuating the intensity of the emergingradiation. The following equation models the attenuation process [19]:I=10e ‘ , (2-1)where I is the incident radiation intensity, I is the emerging intensity, 1 is the path of theradiation, x is the distance along that path, and u(x) is the x-ray linear attenuation whichvaries along 1. Rearranging Equation (2-1) gives:ln(-L)= —Ju(x)dx , (2-2)which shows the linear relationship between log attenuation, ln(L), and linearattenuation. The integral on the right-hand side of Equation (2-2) is called a raysum orprojection of u(x) (raysum is used for the remainder of this thesis to avoid confusionwith projections in the POCS method described in Chapter 4) and is the basis of dataacquisition for CT reconstruction.If the linear attenuation function, u(x), is piece-wise continuous, Equation (2-2)becomes:ln(L) ix, (2-3)where and zx, are the linear attenuation and the thickness of the ith component of theobject along the path of radiation. Thus x-ray raysums are linear functions of linearattenuation coefficients and a linear system adequately models x-ray absorption.252.2.2. Basic Apparatus for CT Data AcquisitionFigure 2.6 shows a schematic of the basic apparatus for a CT data acquisitionsystem [22] [34]. The acquisition procedure is:1. Collect projection data for one orientation:(a) Project a thin beam of radiation (from the source) through thespecimen and measure the output intensity (with the sensor).(b) Move the source and sensor in a direction perpendicular to thedirection of the radiation and repeat 1(a).2. Rotate the specimen and repeat step 1.3. Repeat steps 1 and 2 to collect raysums through r radians.Collimators (not shown in Figure 2.6) confine the x-rays to a thin beam so that raysumsmeasure only a thin slice of the specimen. Successive adjacent slices sample thespecimen in three dimensions if desired.The above process measures discrete samples from the Radon transform of thespecimen. This transform, described in the next section, is fundamental to CT.Figure 2.5: Schematic of x-ray absorption.26Sensor lateral movement of sensor(coordinated with source)rotational movementof specimen________lateral movement of sourceSource (coordinated with sensor)Figure 2.6: Schematic of basic computed tomographyapparatus2.2.3. The Radon Transform and the Fourier Slice TheoremThe Radon transform, first explored by Radon in 1917 [32], is fundamental to CTand forms the basis of the Fourier slice theorem. Together, the Radon transform and theFourier slice theorem provide a mathematical basis for understanding CT reconstruction.Define a Cartesian coordinate system with its origin at the centre of rotation inFigure 2.6. Based on this coordinate system, the Radon transform of a function 1u(x,y)[Rji](l,O), is [22]:[Rji](l, 0) = fii(lcos 0 + sin0, cos 0 + lsin 0)dc. (2-4)where the integral on the right-hand side is a raysum of u(x,y) along a ray at angle 0with the y axis and distance I from the origin. Thus, the Radon transform of an object isthe set of all its raysums. The apparatus of Figure 2.6 measures a discrete approximationof the Radon transform. Measurements are approximate because raysums of Equation (2-274) have infinitesimal thickness (not possible in practice), and because other practicalproblems in radiography such as scatter and beam hardening confound data acquisition.The Fourier slice theorem (also known as the projection theorem) relates theRadon transform to the Fourier transform. A short intuitive explanation of the theoremfollows; a proof is available from many sources, e.g., [22] or [34]. The separable natureof the Fourier transform allows separate computation in x and y directions. Supposecomputation begins with one-dimensional transforms in the y direction, then the zero-frequency component of the resulting transforms are raysums for rays parallel to the yaxis. Completion of the transform on the zero frequency data, i.e., raysum data, gives aslice of the Fourier transform of the specimen in the x direction. Expressedmathematically, the Fourier slice theorem is [22]:[T1Tji](co,0)= L’F,#](oJ1cos 0,—co1sin 0), (2-5)where 9 is the one-dimensional Fourier transform operator in 1, co and 0 are polarcoordinates of the two-dimensional Fourier space, and is the two-dimensionalFourier transform operator in Cartesian coordinates. The Fourier slice theorem says thatmeasuring a set of parallel raysums of an object is equivalent to sampling a slice of theFourier transform along a line through the origin, perpendicular to the direction of theraysums.2.2.4. CT ReconstructionThe Fourier slice theorem allows a short and simple explanation of why CT ispossible. According to the theorem, step one of the data acquisition procedure samples aslice of the specimen in Fourier space. Steps 2 and 3 sample slices through ir radians,thus sampling the entire Fourier space. A complete sampling of the Fourier space meansthat all the data necessary for reconstruction are available. Although tomographicreconstruction by assembling the Fourier transforms from raysums is possible, thismethod is not popular. Of the many methods available, the convolution-back projection28(or some variant) is most popular, and the work presented here focuses on other methods.Nevertheless, the Fourier slice method gives a conceptual model useful for understandinghow raysums sample a function, and the ill-posed nature of limited-angle CT.Inversion of the Radon transform was first shown by Radon [32] and his resultmay be found in more recent literature, e.g. [22] and [34]. The solution to the inverseproblem exists, but does not depend continuously on the data. Thus CT reconstruction isan ill-posed problem. In practice it is not possible to acquire the continuous Radontransform, only a finite number of samples, usually distributed uniformly throughout(1, ). The practical problem of tomographic reconstruction from discretely sampledraysum data is also ill-posed because there is no unique solution.Commonplace use of Cf suggests that its ill-posedness can be managed, as shownbelow. The following linear system models raysum acquisition:Rx=y (2-6)where y is a vector of projection data, x is a vector of linear attenuation coefficients, andR is a linear matrix operator that models x-ray raysum acquisition. (i.e. a discrete versionof the Radon transform operator, 20. Assume that R samples the Radon transformuniformly. CT reconstruction is the inverse problem of solving for x given raysum datay.The matrix R has a non-trivial null space, i.e., the set of vectors, n ,that satisfyRn 0 is not limited to the zero vector. Fortunately, although the null space is nontrivial, the non-trivial components of the null space are high-frequency oscillations.Assuming that the solution should be band-limited modifies the ill-posed problem intoone that is well-posed. This is normally a reasonable assumption because it is known apriori that the oscillations do not exist or they are not of interest even if they do exist.For reconstruction from the continuous transform, the same band-limiting assumptionyields a well-posed problem by ensuring that the solution varies continuously with thedata.29Figure 2.7: Example of CT showing lines of projection on atwo-by two grid.A small example illustrates the ill-posed behaviour of CT reconstruction from adiscrete Radon transform. Consider the two-by-two discrete image shown in Figure 2.7.Arrows show the sampled raysums which represent a uniform sampling of the Radontransform on a two-by-two grid. R for this sampling scheme is:11000011R=10100101The square roots of the eigenvalues of RTR are the singular values of R (see Chapter 4for a description of singular value decomposition, SVD) and are 2, -‘J, and 0. Thezero singular value indicates that R has a non-trivial null space as expected.Suppose we have the following 2 by 2 cross section of a specimen with pixel values:112 2=l2= 12Application of the forward operator, R, to x gives the data:332430From this value of y, the minimum-norm solution, k, (using SVD of R) is:12 12x= =1 122SVD also gives the vectors in the null-space of R. In this case there is one singular valueof zero so there is only one dimension in the null space in the direction of:-y2_y-y y2n-y2As expected, n is an oscillatory image. The reconstructed image, , is accurate becausethe image x is perpendicular to n, as the minimum-norm solution must also be. In thiscase the minimum-norm requirement selects a band-limited solution by rejecting thehigh-frequency components of the null space. If the original image, x, were to containsome component of n, then the minimum-norm solution, , would be in error. CTreconstruction from a uniform sampling is only feasible if the cross section does notcontain high-frequency components, or if those high-frequency components exist, theyare not important.2.3. Limited Angle Computed TomographyThe previous section showed that CT reconstruction is possible from a uniformlydistributed subset of the continuous Radon transform. Reconstruction is an ill-posedproblem which becomes soluble by assuming the solution is band-limited. This sectionpresents the more difficult problem of reconstruction from a set of raysums which isincomplete not simply because of discrete sampling, but also because raysums are onlyavailable for a limited range of angles. This problem arises because physical constraintsimposed on the CT apparatus, either by the specimen or its environment, prevent raysumacquisition throughout the required range of —- 0 - (see Figure 2.8 for examples31taken from [40]). The name given to this particular problem is limited-angle CT. Figure2.8 (a) is important here because it represents the situation with inspection of aircraftcontrol surfaces.X-ray source,j‘..-fiDectector(a)X-ray source,1Dectector(b)Figure 2.8: Example cases for which only limited-angleraysums are available: (a) object is very long in one direction,and (b) other objects obstruct scanning.2.3.1. Ill-Posed Nature of Limited Angle Computed TomographyThe Fourier slice theorem states that sampling a specimen with raysums isequivalent to sampling slices in Fourier space. It is obvious, therefore, that in order tosample a specimen completely one must measure raysums through a range of rradians,i.e. —- 0 -. As illustrated in Figure 2.8, circumstances arise where projectionscannot be sampled throughout a full range, i.e., —e 0 e, and 0< e <-, thus-fl32leaving a portion of the Fourier space unsampled. When the full range of angles isavailable, then the sampled Fourier space covers a disc-shaped region. In the limitedangle case, a portion of the Radon transform is missing and the sampled Fourier spacecovers only the butterfly-shaped region shown in Figure 2.9.—xFigure 2.9: Effect of limited-angle Radon transform onFourier domain sampling: (a) range of angles sampled(shaded region) and (b) effective sampling of Fourier domain(shaded region is sampled).Fourier space sampling illustrates clearly the ill-posed nature of limited-angle CT.Ill-posedness occurs because there is no unique solution, but instead of a small null spacecontaining only high-frequency oscillations, the null space also contains componentscorresponding to the unsampled Fourier space. Although both the full-range and limited-angle problems are ill-posed, we say the limited-angle problem is more ill-posed becauseit has a greater number of dimensions in its null space. Davison [101 and later Louis [23]use singular value decomposition to explore the ill-conditioned nature of limited-angleCT. Louis concludes “... the limited angle problem is much more ill-posed than the fullrange problem. ... Clearly the algorithms applied in this situation have to be carefullydesigned and it is desirable to use as much a priori knowledge as possible.” Louis alsoshows some images which are vectors in the null space of R. These null-space imagesshow what types of artifacts may appear in limited-angle CT reconstructions.(a) (b)332.3.2. Reconstruction from Limited-Angle DataThe Fourier slice theorem shows that limited-angle raysum data do not sample theFourier space completely. Therefore, any reconstruction method based on limited angledata necessarily fills in the unsampled regions of Fourier space. Possibilities for filling inare unlimited which is precisely the source of the ill-posedness cited above. A prioriassumptions can give constraints that lead to a unique solution. Current techniques forlimited-angle CT fall into these three categories distinguished by the a priori constraintsthat guide the reconstruction to a unique solution and complete the Fourier space:1. fill the unsampled Fourier space with zeros,2. impose spatial support constraints on the solution, and3. extrapolate the sampled Fourier space data.The first category, fill with zeros, is the minimum-norm solution. Of all possible valuesfor the unsampled region, all zeros contains the least amount of energy and thereforemust be the minimum-norm solution. The second category, impose spatial supportconstraints is similar. Constraining spatial support does not necessarily determine themissing data on its own. One still seeks a minimum-norm solution as in category one,but with the additional constraint added. The last category, extrapolate the Fourier space,assumes that the available data establish all important trends in the Fourier space, anassumption that is not necessarily valid. Examples from all three categories follow.Section 2.2.4 showed briefly the minimum-norm solution to a system of equationsas computed by SVD. Refer to the example in Section 2.2.4. The discrete raysumoperator for the limited angle case is:El 1 o oR=ILO 0 1 1with singular values of -‘.I, -‘.J, 0, and 0. Restricting the angle adds another singularvalue of zero to the problem and increases the number of dimensions in the null space34from one to two. Two orthogonal vectors form a basis for the null space, one identical tothat for the full-range case, i.e.:-y2y2_-y y21-y2and the other, which is entirely attributable to the limited-angle sampling, is:-y2n2 y y2 y2y2The new null space component is a step across a vertically oriented edge. If x is:112 2x= =12 12then the raysum data, y, is:E31Y= [jand the minimum-norm reconstruction, , is:_y2 y2Xy y2.y2Note that the vertical edge present in the original image is missing in the reconstructionbecause the minimum-norm solution omits the portion of x in the null space, i.e., thevertical edge. Suppose instead that x is:111 1x = =22 2235then the minimum-norm reconstruction is:11 11x= =2 222The reconstruction is accurate in this case because x is orthogonal to the null space of R.The salient point here is that minimum-norm reconstruction works so long as the functionbeing reconstructed contains no component in the null space, and the success or failure ofreconstruction is case-dependent.As one might reasonably guess, any reconstruction method intended for full-rangedata can compute the minimum-norm solution by assuming that all the missing raysumsare zero (Tuy [43]). Consequently, one does not need to use the unwieldy and slowsingular value decomposition to compute the minimum-norm reconstruction fromlimited-angle data.The second category of limited-angle CT uses the often reasonable assumptionthat an object lies within known spatial bounds. This is a common approach toimproving the quality of limited-angle CT reconstruction. GrUnbaum [16], for example,makes the a priori assumption that the specimen is contained within a rectangular regionof support and linear attenuation values outside that region must be zero. The limitedspatial support effectively interpolates the Fourier domain, and so there is an implicitassumption to justify interpolation that the Fourier space is smooth. GrUnbaum givesresults for synthetic data but there is no indication of the circumstances for which theeffective interpolation is valid.Although the GrUnbaum rectangular region improves upon minimum-normsolutions, in some circumstances it is possible to further constrain spatial support. Sato etal. [33] employ a triangular region of support to reconstruct a synthetic image of theletter ‘A’ (see Figure 2.10). They perform reconstruction using an iterative method thatalternates between the Fourier and spatial domains, applying constraints in each:36triangular support in the spatial domain and a butterfly-shaped region in the Fourierspace. They provide a proof of convergence, but not a rate of convergence. Theirmethod is, in fact, the method of projection onto convex sets (see Chapter 4) whichguarantees convergence. As is the case for Grünbaum, examples are for synthetic dataonly. Such tight restrictions on spatial support are useful but not trivial to implement inpractice.AFigure 2.10: Triangular region of spatial support from Sato etaL [33] method for limited-angle CT reconstruction of theletter ‘A’.GrUnbaum and Sato et al.. are just two examples of using spatial supportrestrictions to aid limited angle CT reconstruction. There are more examples in theliterature, such as Dusaussoy and Abdou [11] and Tam et al. [40]. It is interesting to notethat Tam’s algorithm appears to be identical to that of Sato et al.The third category of reconstruction method is interpolation of the Fourierdomain. Methods in this category extrapolate data in the Fourier domain to fill in theunsampled regions. Soumekh [36] models the Fourier domain with a periodic functionthat can be represented by a Fourier series expansion. He then makes the a prioriassumption that only a finite number of coefficients are necessary to represent the series,i.e., coefficients for high frequencies are zero or insignificant. The available Fourierspace data are sufficient to determine the finite number of coefficients. The truncatedFourier series then interpolates to fill in the missing Fourier domain data. Although theinterpolation by Soumek does not translate into simple assumptions about spatial support,it does restrict the class of objects that may be accurately reconstructed. Reconstruction37is successful with synthetic examples, but the extrapolation requires that the availabledata establish all important trends. Heiskanen et al. [20] give another example of directinterpolation of raysum data.Minimum norm solutions assume that there is no important information in thenull-space of the raysum operator. Spatial support restrictions and direct extrapolationmethods assume that the sampled data establish all trends in the specimen, i.e., everythingyou need lies in the range of the limited-angle raysum operator. The common featureamong all categories is that none can properly reconstruct a function unless the null spaceis either trivial or unimportant. In fact, the following is true for limited-angle CT:1. limited-angle CT requires filling in the unsampled Fourier space of afunction,2. there is no universally correct way to fill in the missing data, and3. a priori assumptions should guide the reconstruction as much as possible.This thesis proposes moving beyond a priori assumptions in limited-angle CT and usingdata from other sources to constrain the reconstruction, i.e., constraint-based data fusionfor limited-angle CT. Data fusion fills in the missing data correctly, even when a prioriassumptions for interpolation are not valid.2.4. UltrasoundThe previous section showed that CT reconstruction from limited-angle data is ill-posed and has a large null space. To reliably reconstruct tomographic images of aspecimen, limited-angle CT methods rely significantly on a priori assumptions toconstrain the solution. Constraint-based data fusion avoids reliance on such assumptionsby using other data sources to provide constraints. Ultrasound is a commonplace NDEtechnique that gives data complementary to x-rays and has the potential to provideconstraints to limited-angle CT. This section presents some background on ultrasound asa prelude to its use in a limited-angle CT data fusion system.382.4.1. Forward ModelsIn ultrasonic testing, a transducer transmits a high-frequency sound wave into thesurface of a specimen. The sound travels through the specimen unchanged until itencounters an interface between materials of different acoustic impedance. At theinterface some energy is reflected and some is transmitted. Reflections (echoes) at thesurface contain information about the internal structure of the specimen. Take, forexample, a single interface between two layers in a layered specimen; the upper andlower layers have acoustic impedances z1 and z2 respectively. The amplitude of areflected signal is r, the reflection coefficient for the interface, times the amplitude of theincident signal, where r, is given by [24]:r=Z21. (2-7)z2 + z1A negative value of r indicates a phase reversal. The following expression gives theechoed signal:y(t) = s(t)* w(t) + r(t), (2-8)where y(t) is the echoed signal, s(t) is the impulse response of the material, w(t) is thesound wavelet introduced at the surface, and r(t) is noise (* is the convolution operator).s(t) contains the desired structural information. It is a train of impulses, each impulsecorresponding to an interface in the specimen. The magnitude and sign of the impulsesdepends on r for the interface and the attenuation in the path of the sound to and from thatinterface. Equation (2-8) gives a simple, non-parametric (i.e. material parameters are notspecified), forward model of the ultrasonic process [25].More sophisticated parametric models define the output of an ultrasonic systembased on material properties. These models fit into two categories: (1) equal time slicesand (2) variable time slices [25]. The models assume a lossless layered medium like thatshown in Figure 2.11. m(t) is the signal injected into the specimen and y(t) is the output39as in Equation (2-8). u(t) and d(t) are the upward and downward traveling signals inthe th layer.\JroLayer iFigure 2.11: Schematic of a lossless layered ultrasonicmedium: (a) top layer and (b) internal layers.The equal time slice model divides the medium into layers of equal travel time.Each layer is a half sample period thick, which makes the two-way travel time equal toone sample period per layer. The set of reflection coefficients, r0,r1 ,r2,. .., r, where n isthe number of samples, completely defines the model. The following equations describethe ultrasonic system:d1÷(t)= (1+r)d(t—--) — ru+1(t)for internal layers, andd1(t) = (1+r0)m(t) — r0u1(t)y(t) = r0 m(t) + (1—r0)u1(t)(a)(b)u(t+---)= j=1,2,...,n(2-9)(2-10)40for the top surface.The variable time slice model divides the specimen into slices corresponding tothe actual material layers. Reflection coefficients, r0 ,r1r2,... , r, and the travel times,1’ ‘r2 ,..., ‘r,, where n is the number of layers in the specimen, completely define themodel. The following equations describe the ultrasonic system:d(t+ ‘r)=(l+r_1)d_(t)—r1_ u(t)(2-li)u(t+ r) = rd1_(t) + (1— r)u1+(t)for the internal layers, andd0(t)=m(t) (2-12)y(t) = u1(t)define the input and output (in this case for a unidirectional sensor just beneath the uppersurface). Variable time slices have the inherent advantage that the parameters corresponddirectly to the structure of the specimen.2.4.2. Ultrasound InversionBoth the non-parametric (Equation (2-8)) and the parametric (Equations (2-9) to(2-12)) models have corresponding inverse problems reviewed here.Non-parametric modeling leads to the inverse problem of ultrasonicdeconvolution. Equation (2-8) expressed with linear matrix operators is:y=Ws+r, (2-13)where y is the vector of output data, s is the vector of impulses, and r is an error vector.W is a matrix operator which performs a convolution with the ultrasound wavelet. Theleast-squares approach to solving Equation (2-13) gives the optimization problem ofminimizing the quadratic norm of the residual, Ws—y, i.e., minimizing the objectivefunction E(s) where:E(s) = Iws—(2-14)41Unfortunately W is not full column rank and the optimization is ill-posed; the solution isnot unique. One way to deal with the ill-posedness is to modify the problem byregularization. Rather than solve for s from Equation (2-13), solve for s from:Ewl Fyi[D]S[OJ. (2-15)Solution of Equation (2-15) by least-squares minimizes the objective function:E(s)= liws— y112 + A7NDsII2. (2-16)rwiD is chosen so that the matrix [Dj is full column rank. Minimizing E(s) in (2-16)gives the following symmetric positive definite system:(WTW+A2D)s= WTy, (2-17)which is essentially the classic Wiener deconvolution [24]. In Wiener deconvolution theoperator D corresponds to multiplication by the noise spectrum in the frequency domain.A more recent method called 4-norm deconvolution uses 4 norms rather than 12and assumes that s(t) contains impulses only [45]. Based on experimentation,developers of the 4-norm method state that the 4 norm gives superior results and is fasterwhen compared to the12-norm, and that the assumption that s(t) contains only impulsesleads to improved temporal resolution.Non-parametric inversion only recovers the impulse response, s(t), and no directstructural information. Inversion with a parametric model allows recovery of structuralinformation in the form of an acoustic impedance profile. There are two basic methodsfor parametric inversion: direct inversion and optimization.The first, direct inversion, solves for the uppermost interface from the first datapoint (Equation (2-10)). It then solves for subsequent layers from subsequent data pointsand solutions to upper layers in top-to-bottom order. There are two major problems withdirect solution. First, it assumes that the data are deconvolved, i.e., the data are the seriesof impulses s(t) rather than the convolved data y(t). Second, the method is ill-42conditioned. Small errors in a layer propagate to lower layers and grow rapidly inamplitude.In the second type of parametric inversion, optimization, one creates an objectivefunction:E(z)= IIU(z)— II2 (2-18)where z is the impedance profile and U(z) is the forward ultrasonic process, and thenminimizes it among the possible z. U(z) is non-linear so E(z) is non-convex anddifficult to minimize. The minimum may not be unique and so the optimization is illposed. In spite of these problems, optimization is superior to direct inversion.Examples of parametric inversion are given by Goupillaud [151, who uses anequal time slice model for ultrasonic inversion and more recently Habibi-Ashrafi andMendel [18], Mendel and Goutsias [26], and Zala and McRae [45] who use a variabletime slice model. Most current literature on ultrasonic inversion treats the material aslossless although in reality this is rarely the case. Usually the specimen absorbs somesound energy in a manner analogous to x-ray absorption. The absorption coefficient isdifficult to measure precisely and ultrasonic inversion is extremely sensitive to it. Soultrasonic inversion for an absorbing medium is even more complex than for animaginary lossless medium. Zala and Churchill [46] give an excellent review ofultrasonic inversion methods as well as the practical problems encountered when themedium is not lossless.2.4.3. Minimal Ultrasound for Fusion with Computed TomographyUltrasound data complement x-ray data because they measure features that x-rayscannot. Whereas raysums are insensitive to discontinuities along the path of radiation,ultrasound explicitly detects discontinuities in its path. This complementary nature of x -ray and ultrasound portends their effective fusion in limited-angle CT.43Ultrasound inversion for a layered lossless specimen is difficult but possible. Inthe more realistic situation of an absorptive specimen ultrasonic inversion is tenuous; thesolution is extremely sensitive to the absorption coefficient which is generally unknownand highly variable.Fortunately, CT reconstruction of sandwich specimens by fusion with other datasources does not require the full potential of ultrasound. Instead, as shown in Chapter 3,only measurements of the face sheets are necessary. As long as the external layer iscomposed of a homogenous slab, thickness measurement only requires identification oftwo reflections in the ultrasound signal: one from the top surface and one from thebottom. The two reflections are distinct making their identification easy and thicknessmeasurement simple.There are several reasons why thickness measurement might not be simple.Delaminations in a composite face sheet can prevent detection of bottom surfacereflections. Face sheets in some regions of a specimen may have multiple layers and sothe reflection from the bottom must be separated from other reflections. In spite of theseconfounding factors, it is reasonable to assume that ultrasound can provide sufficientthickness measurements of the face sheet of a sandwich specimen to allow accurateinterpolation of thickness data over the entire face sheet.2.5 Chapter SummarySection 2.1 of this chapter introduced the constraint-based data fusion model ofClark and Yuille [9]. Their elegant approach implements fusion by using the output ofone data inversion module as constraints for another. Constraint-based data fusionprovides a coherent framework within which to build data fusion systems.The basic inverse problem of x-ray raysum inversion is known as CT (computedtomography). Conventional CT is ill-posed, but the ill-posedness is easily handled bysuppressing high-frequency components in the reconstructed solution. Practical44considerations often dictate that raysum data cannot be collected over the full range ofangles in the Radon transform. In these limited-angle cases the inverse problem is stillill-posed, but the null space is much larger and reconstruction requires additionalconstraints. Many methods exist to solve limited-angle CT. They all rely on a prioriassumptions (which may not be valid) to constrain the solution. This reliance portends amore-accurate solution from constraint-based data fusion when a priori assumptions arenot valid.In the field of non-destructive evaluation ultrasound is commonplace. It providesa source of sensory data that complements x-rays and may provide constraints to limited-angle CT. Many types of inspection are possible with ultrasound, but the limited-angleCT problem of this thesis requires only the ability to measure the thickness of a face sheetin a sandwich structure.Chapter 3 shows that limited-angle CT fails to sample adequately the face sheetstructure of a sandwich specimen causing limited-angle CT reconstruction methodsconstrained by generalized a priori assumptions to fail. The chapter then proposes anovel system for limited-angle CT that fuses ultrasound data with x-ray data to accuratelyreconstruct sandwich structures.45Chapter 3Limited-Angle Computed Tomographyfor Sandwich Structures Using DataFusionSection 3.1 of this chapter defines the null space of the limited-angle Radontransform. The nature of the null space suggests that wide structures withcorrespondingly thin Fourier transforms will lie in the null space and be invisible inlimited-angle raysum data. Face sheets in sandwich structures have this property.Therefore, accurate limited-angle CT reconstruction of sandwich structures is notpossible because the face sheets are not properly measured and interpolation of the data isinvalid.The limitations of limited-angle CT present an opportunity to use constraint-baseddata fusion to improve reconstruction. Sections 3.2 and 3.3 introduce a novel stronglycoupled feed-forward fusion system for limited-angle CT. The system exploits thecomplementary nature of x-ray and ultrasound by using the latter to measure face sheetfeatures that are otherwise invisible in raysum data. Fusing face sheet measurements intoreconstruction removes the face sheets from the null space. Thus, the fusion system46produces an accurate reconstruction of a sandwich structure where it would otherwise beimpossible.Section 3.2 describes the data acquisition process for the three sources of data:range finder, ultrasound, and x-ray. Section 3.3 shows how the fusion system combinesthe multiple data sources, including the mathematical formulation of the applicableinverse problem.3.1. Null Space of the Limited Angle Radon TransformDefine the limited-angle Radon transform of the cross section Jt(x, y),[e/_U1(l, 0), as:[eIh1(l,O)[!RJM(t,O): —6 o+e, 0< 9< ,r/2. (3-1)The limited-angle transform is identical to the Radon transform except that it is definedover a limited range of angles only. By definition, a function p, is in the null space ofthe Re if [RLi J(l,O) = 0. Any structure in a specimen that lies in the null space doesnot contribute to the data measured by the transform. In this sense, such structures areinvisible to data acquisition. It is important, therefore, to know the extent of the nullspace and which structures are contained within it. The remainder of this section definesthe null space of the limited-angle Radon transform and an important class of structuresthat lie almost entirely within it.As stated above, by definition ,u,, e if and only if:1(1,0) = 0. (3-2)Taking the Fourier transform of Equation (3-2) gives:LF1Q9u1((o,O)=0. (3-3)By the Fourier slice theorem (see Chapter 2, Section 2.2.3):[.Tt9eImn](co1,0) = ](a cos 0, —o sin 0), (3-4)and[.Til(coo)= 0, — e 0 e, tanO -s-. (35)47The case where cot, coy, = 0 requires special consideration. The limited-angle transformdoes sample this case, and, as for Equations (3-2) through (3-5):[9uJ(O,0) = 0. (3-6)Therefore the null-space is:O e,(3-7)tanO=—, [u}(O,0)=0(oxIn words, the null space is the set of all functions whose Fourier transform is zero withinthe sampled regions of the Fourier domain.Equation (3-7) illustrates the effect of limiting the range of angles. As edecreases the portion of a function (in Fourier space) that must be zero for the function tolie in the null space also decreases and the null space becomes larger. A practicalconsequence of a larger null space is that more of the structure of a specimen lies withinit and is not sampled by the limited-angle Radon transform. Although CT reconstructionis ill-posed for both large and small null spaces, reconstruction is more difficult with thelarge null space because more structural information is missing.Equation (3-7) also hints at what types of structures lie in the null space. As arule, narrow structures have wide Fourier transforms and wide structures have narrowFourier transforms. One expects, therefore, that a wide structure will lie in the null space,provided that it is oriented so that its narrow Fourier transform lies in the unsampledregion. An example is the face sheets of sandwich structure. It is well known that theouter extremities of a structure carry stresses due to bending. Therefore, to achieveeconomy of weight and material, designers often create structures to carry loads in a thinouter shell separated by some core material, i.e., a sandwich structure. One example isthe control surfaces of an airplane. The outer shell, i.e., the face sheets, carries the loadswhile the purpose of the core is to provide a rigid connection between the face sheets. Ifthe structure is wide, as is the case for aircraft components, it does not allow acquisition48of full range raysum data, and furthermore, the face sheet structure lies almost entirely inthe limited-angle Radon transform null space. We now show that this is the case.The function Srad (x, y), shown in Figure 3.1, represents a hypothetical facesheet and is given by:d . d,_Ja, ----(xsin+ycosØ)—r--Srt,adXY)— 2 2 , ( - )0, otherwisei.e., a rectangular function in the y-axis direction with amplitude a and width d, shifted rfrom the origin in the y-axis direction and rotated by angle •. The Fourier transform ofSrad (x, y), [Yl(,ySrI,ad 1( cot, wa,), is:[9,ySrpad ](o), Oi),) = c(O )ad sinc(rdo )e’°’, (3-9)whereri — r05 _sin1Fax1[o;f[sin cos J[oj’i.e., the Fourier transform is a sinc function along the y axis, phase shifted (for r) androtated by angle .Figure 3.1: Function Sr,ad(X,Y); shaded region has value a and zeroeverywhere else.Figure 3.2 shows a sandwich structure with two face sheets surrounding somecore material. The width of the structure restricts the range of angles to — e 0 ex49where e <- — . Thus, with the exception of the zero-frequency components, the facesheets lie within the null space of the limited-angle Radon transform. The presence of theface sheets confounds reconstruction. Their position in the null space eliminates thepossibility of accurate reconstruction based on general a priori assumptions, so specificinformation about the face sheets is essential. Data fusion solves this problem byincorporating face sheet measurements gleaned from other data sources. The remainderof this chapter describes a novel data fusion system for limited-angle CT that acquiresface sheet data and fuses them with raysum data in reconstruction.3.2. Data Acquisition for Limited-Angle CT Data Fusion SystemClearly, limited-angle CT for sandwich structures based solely on raysum data isnot feasible. The face sheets of the sandwich lie in the limited-angle Radon transformnull space so that the assumptions necessary for interpolation are insufficient. Thissection proposes a novel system for limited-angle CT, starting with a description of thedata acquisition. The system is based on raysum acquisition with modifications to allowyFigure 3.2: Sandwich specimen with face sheets at angle . Thewidth of the specimen prevents raysum acquisition outside the range— 9 0 e. Note that e <(,r/2) — , so the face sheets lie almostentirely within the limited-angle Radon transform null space.50face sheet measurements from range finders and ultrasound. A description of a strongly-coupled feed-forward data fusion system that combines the data follows in the nextsection.3.2.1 Raysum Data AcquisitionFigure 3.3 shows schematically the raysum data acquisition apparatus for thinobjects. The apparatus sits on a large C-shaped manipulator (the C-arm) that allows thex-ray source (mounted on one end of the C) and the sensor (mounted on the other end ofthe C) to maneuver about a thin specimen.Linear arrayx-ray tubeFigure 3.3: Raysum data acquisition system.To acquire an image, the C- arm sweeps the array across an area. A frame bufferassembles the image from the lines of data acquired by the array. Image acquisition isnot important for CT, but a carefully arranged set of raysums is. Figure 3.4 shows,looking at the linear array end-on, the C-arm motion during data acquisition. The source-sensor array scans the specimen along a line as if acquiring an image, but only pixelsfrom the centre of the array are used. The projection angle (the angle between verticaland the source-sensor plane) is initially zero. Subsequent scans acquire projection datafor angles in the range —0 0 9.I II Ithin specimenIIIIII I“I raysum data51x-ra sourcex-ray source/ \ direction of scan/k\)\\\ \\.\\\Nspecimen specimen• p direction of scanlinear array linear array(a) (b)Figure 3.4: Scan method for limited-angle raysum acquisition:(a) parallel with vertical and (b) at an angle to vertical.3.2.2 Range DataAn important assumption often used in limited-angle CT is that of limited spatialsupport for the specimen, i.e., the specimen fits into a bounding region. Constraint-baseddata fusion improves solutions of inverse problems by substituting reliance onconstraining assumptions with reliance on measured data. It is sensible, therefore, tomeasure the bounding shape of the specimen rather than rely on assumptions if possible.Fortunately, devices exist that measure precisely the bounding shape of an object.An example is a laser range finder, an instrument that measures the distance from theinstrument to an opaque surface in front of it. It projects a small dot of laser light ontothe surface of the specimen. The optical system of the range finder detects the dot andtriangulates its position. By sweeping the dot along a line and triangulating at eachposition, the range finder acquires range data along a line. Whereas the linear array x-raydetector collects raysum data along a line, a laser range finder collects range data along aline. To collect range data over an area, the range finder sweeps across the area in amanner analogous to the linear array x-ray detector when acquiring an image. Oneexample of a commercially available laser range finder is the Saturn-2000 built by ServoRobot [35J. Two range finders mounted on the C-arm apparatus, one beside the x-ray52source and the other beside the detector, can measures the bounding region of thespecimen with a single sweep of the C-arm. Figure 3.5 shows the modified apparatus foracquisition of raysum and range data.laserrange finderx-ray sourcedirection of scanspecimenlinear array laserx-ray detector • range finderFigure 3.5: Scanning method modified to include laser rangefinder for bounding region measurement.Laser range finders are not the only method of determining the bounding region,but they are available and are accurate. An alternative is stereo vision or a modifiedstereo system. Two video cameras mounted in place of the range finders can acquire asequence of images as the apparatus sweeps across the specimen. From these images itshould be possible to compute range data, i.e., compute depth from motion where themotion is known. This idea is similar to stereo vision but whereas stereo uses only twocameras to compute depth, the sweeping motion of the C-arm provides many images.3.2.3 Thickness DataMeasurement of the bounding region of the specimen by a range finder gives theposition of the outer side of each face sheet. To completely measure the face sheet it isalso necessary to know the thickness of the face sheet. Ultrasound is well suited for thismeasurement.To acquire thickness measurements use a set of ultrasound pulse-echo traces takenalong the face sheet in the scan direction. Ideally the traces have two reflections, onefrom the outer surface and one from the inner surface. From the position of the two53reflections calculate the time for the sound to travel between the two surfaces. Assumingthat the speed of sound in the face sheet material is known compute the thickness of theface sheet from the speed and travel time.Figure 3.6 shows the apparatus for ultrasound data acquisition. Conceptually, theultrasound system sits on the C-arm with the range finders and x-ray system. Practicalrequirements for ultrasound (e.g. water squirters) and x-ray (high-voltage circuits)prohibit their coexistence. Therefore, the ultrasound apparatus is separate from the x-rayand range apparatus.ultrasonictransducerultrasoundenergy in_________________couplant direction of scanspecimenultrasoundenergy incouplantultrasonictransducerFigure 3.6: Apparatus for ultrasound data acquisition.It is possible that thickness measurement may fail occasionally. This can occurfor many reasons, e.g., scattering obstructions, defects within the face sheet, or internalstructures which confound the reflections. Fortunately, the thickness of the face sheets isnot rapidly varying and it is a straightforward process to interpolate missing thicknessvalues.3.3 The Fusion SystemThis section introduces a novel data fusion system for limited-angle CT ofsandwich structures based on the apparatus shown in Section 3.2. Section 3.3.1 shows54the flow of data which is essentially a feed-forward strongly-coupled data fusion system,and Section 3.3.2 shows the mathematical formulation of the relevant inverse problem.3.3.1 Data FlowA strongly-coupled feed-forward data fusion system processes the data acquiredby the apparatus of Section 3.2. Figure 3.7 shows a block diagram of the system.input data__________________map of spatial supportrange data— Bounding RegionComputationa prioriassumptions Imap of spatial supportand face sheet valuesultrasound data _eShee_}__a prioriassumptionsreconstructedcross-section of specimenraysum data IIFigure 3.7: Block diagram of strongly-coupled feed-forward datafusion system for limited-angle CT of sandwich structures.The first module of the system, bounding region computation, takes range data asinput and produces a map of the spatial support for the specimen. The map segments thecross section into exterior regions and regions within the specimen.The second module, face sheet thickness computation, takes ultrasound data asinput and computes a map segmenting the specimen into exterior regions, face sheetCT Reconstruction55regions, and interior regions. The ultrasound data, constrained by the known speed ofsound, give face sheet thickness, and spatial support constraints define the positions ofthe face sheets. Combined, the face sheet positions and face sheet thickness give a mapof the face sheet regions. Regions that are not exterior nor in a face sheet must be, bydefault, in the specimen interior.The map of the exterior and face sheet regions constrains the third module, CTreconstruction, whose input and output are the raysum data and the reconstructed cross-section respectively. As shown in Section 3.1, the face sheet structures are invisible tolimited-angle raysums, and thus confound accurate reconstruction. The complementarynature of x-ray and ultrasound data allow the segmentation based on range and ultrasoundto constrain the reconstruction to overcome this limitation. Whereas with the raysumdata alone the reconstruction fills the unsampled Fourier space based on invalid a prioriassumptions, the new constraints restrict the solutions to conform to the face sheets.These constraints do not necessarily make the problem well-posed so a prioriassumptions are still necessary for reconstruction, but these assumptions no longer lead toerroneous results. Wide internal structures with edges parallel or near parallel to the facesheets are still a problem. In practice though, internal structures are much narrower thanthe face sheets so more of their Fourier transform lies in the sampled Fourier space,allowing better interpolation from a priori assumptions.3.3.2. Mathematical FormulationAcquisition of raysums is represented by the discrete linear system:Rx=y, (3-10)where x is a vector of linear attenuation coefficients in a specimen and y is a vector ofraysums. R is a matrix representing the discrete Radon transform. Each raysummeasurement is given by the inner product:y1=(r,x), (3-11)56where y is the ith raysum and r is the th row of R. The th element of r., r,j , is thelength of the jth ray that passes through the th element of x. For limited-angle raysumacquisition, only rows corresponding to the raysums that are possible to measure are in R.The available raysums determine R and the number of singular values equal to zero.From Section 3.1 we know that R has a non-trivial null space containing the face sheetstructures. Therefore, accurate CT reconstruction based on Equation (3-10) alone is ill-posed.Range data and ultrasound data provide additional constraints so thatreconstruction recovers the face sheet structure. The range finder gives a boundingregion for the specimen. This leads to a map of regions inside and outside the specimenand the following system of equations:WX=WXF (3-12)where:W—_diag(w1,w2..., J,...,w)11 x. outside specimenw.H. .. ,and‘ 0 xj inside specimen— fl1air xj outside specimen— 10 x1 inside specimenThe subscript F denotes fusion. is the linear attenuation coefficient for air and iszero. XF is a partial reference image based on the fusion data. The range finder datainversion gives the elements on the diagonal of W. Intermediate values between zero andone on the diagonal indicate that a particular element lies across a boundary betweenregions.The ultrasound module gives the thickness of the face sheets. Face sheetthickness and the location of the outer side of each face sheet determine which elementsof x are exterior, which are in a face sheet, and which are interior. This gives a newmatrix W and partial reference image XF where:5711 x. in face sheet or exteriorw.=. . .. and[0 x in specimen interiorIair xj outside specimen= #face sheet xj in face sheet0 in specimen interior#face sheet is the linear attenuation coefficient of the face sheet material.Equations (3-10) and (3-12) combine to form the system:[Z]x=[vJ] (3-13)Solving for x from Equation (3-13) yields a solution that is constrained by Equation (3-12) and therefore constrained by the range and ultrasound data. Although Equation (3-rR12) further constrains the problem, the reconstruction is still ill-posed. The matrix [wis not necessarily square and is not full column rank. Equation (3-13) eliminates the facerRlsheet structures from the null space of , allowing a reconstruction that properlyaccounts for the face sheets. Interior structures with edges parallel or nearly parallel tothe face sheets are still at least partly in the null space.It is possible to weight the fusion data in Equation (3-13) to a greater or lesserdegree by using other values on the diagonal of W, where a large weight indicates agreater degree of confidence in the constraint. For example, one can choose the weightsto be the reciprocals of the standard deviations of constraint measurements [39]. Theunitary weights are arbitrary here, but they happen to produce good results using singularvalue decomposition and regularization (see next chapter). Weights in W are irrelevant tothe method of projection onto convex sets (also in the next chapter) where the parametersthat control constraint set size determine the relative weighting of constraints.In the context of the data fusion system proposed here, strongly-coupled feedforward data fusion consists of:1. forming the system of equations, and582. solving for x from Equation (3-13).It remains to specify a method of solution. Chapter 4 presents and compares threepossible methods.3.4 Chapter SummaryThe null space of the limited-angle Radon transform is the set of all functionswhose Fourier transform is zero over the region sampled by the transform. This nullspace specifically includes the face sheets of sandwich structures, with the exception ofthe zero-frequency component. Therefore, an important part of a sandwich structure isinvisible to limited-angle raysums and accurate reconstruction is not possible fromraysum data alone.A novel limited-angle CT system, proposed in this chapter, uses constraint-baseddata fusion to solve this problem. The system specifies data acquisition methods forrange data, ultrasound data, and x-ray raysum data. A strongly-coupled feed-forwarddata fusion system combines the data from the different sources.The range data define a bounding region for the specimen. Outside that regionlinear attenuation must be zero. Range data also locate one side of each face sheet.Ultrasound data give thickness measurements over the face sheet. Thicknessmeasurements constrained by the face sheet locations determine which regions of thespecimen are face sheet. Face sheet regions have a known linear attenuation. A part ofthe specimen that is neither outside the bounding region nor in the face sheets is, bydefault, in the interior of the specimen. No assumptions are made about linearattenuation in the interior.The segmentation of the specimen into exterior, face sheet, and interior regionsgives a set of constraints expressed by the linear system:Wx = WXF59where W is a diagonal matrix whose elements indicate whether or not an element of x isin a region where linear attenuation is known. XF is a partial reference image based onthe fusion data.Raysum acquisition is modeled by the linear system:Rx=y.Range data and ultrasound data constrain the reconstruction from raysums in thecombined linear system:FRi FIx=I[WJ [WxFData fusion in this context consists of building the above system of equations and solvingrRlfor x. The matrix I I is not square, nor is it full column rank, so the inversion is still ill[WJposed. However, the new system reduces the size of the null space so that it no longercontains the face sheets. Only internal structures with edges parallel or nearly parallel tothe face sheets are in the null space. In practice, internal structures are narrow enoughthat limited-angle data should recover them adequately.The next chapter shows three numerical methods suitable for solving for x fromEquation (3-13).60Chapter 4Numerical Methods for Limited-AngleTomography SystemChapter 3 introduced a novel method for limited-angle computed tomography thatexploits data fusion to properly reconstruct a cross section of a specimen with face sheets.The method ultimately requires solution of the linear system of equations:Ax = bwhere:ER1 FyA=I I and b=I[WJ LWXFThe nature of A dictates that solving for x is an ill-posed problem; if a solution exists, it isnot unique.A plethora of methods exists for solving linear systems. This chapter reviews thefollowing three methods suitable for the limited-angle CT fusion problem:1. singular value decomposition (SVD),2. regularization and the conjugate gradient method (RJCG), and3. projection onto convex sets (POCS).Note that these are not the only methods possible.61A brief description of each method follows, including details of application to theproblem of interest, and its advantages and disadvantages. Section 4.4 summarizes themethods in tabular form. Results of application of each method to synthetic data, andapplication of RICG and POCS to real data are in Chapter 5.4.1 Singular Value DecompositionStrang [38] suggests that singular value decomposition (SVD) “... is not nearly asfamous as it should be.” Whether or not this is true, SVD is an excellent tool for analysisof linear systems. This presentation of SVD, based on Golub and Vanloan [13], Stoerand Bulirshch [37], Forsythe, Malcolm and Moler [12], Strang [38], and Press et al. [31],focuses on the decomposition itself rather than the algorithm that computes it.4.1.1 Properties of the DecompositionStart with the generic linear system:Ax=b. (4-1)We know nothing about A except that it is m x n and m n. For m < n fill A with rowsof zeros to get m — n. SVD decomposes the matrix A into the product of three matrices:A = UEVT. (4-2)The following tableau elucidates the decomposition:u vTA= . (4-3)mXm mxn nxnU and V are both orthogonal, but the significance of SVD lies in the structure of :01 0O•20 . (4-4)062The values cr, relate to the eigenvalues of ATA; o, = , where is the ftheigenvalue of ATA. Since ATA is positive semi-definite, A. 0 and 0. Becausethe last m— n rows of I are zeros, the last m — n columns of U are unnecessary.Ignoring these rows and columns gives the abbreviated decomposition:U I VTA= , (4-5)mXn nXn nXnwhere:01= . . (4-6)0Rearranging Equation (4-2) gives:AV=UI, (4-7)from VT = I. Individual columns in Equation (4-7) are:Av1 = o,u. (4-8)Consider Equation (4-8) for each of two categories for1) =OAv =0, soE ilV(A)2) o,>O1A—v1 = u, soe Range(A).The following items summarize important aspects of the decomposition:1. If =0 for any j then the matrix ATA is singular and the linear systemof Equation (4-1) is under determined,632. The null space of A, W(A), and its orthogonal complement, 9V.1-(A), aregiven by:= span(v : = 0), and= span(v : 0),and3. The range of A, Range(A), and its orthogonal complement, Range’(A)are given by:Range(A) = span(u : 0), andRange’ (A) = span(u1 : cr, = 0).SVD indicates whether or not a system of equations is under-determined, and it computesorthogonal bases for the null space and range of a matrix and their orthogonalcomplements.4.1.2 Least-Squares Minimum-Norm SolutionIn addition to finding the null space and range of a matrix, SVD leads directly tothe least-squares minimum-norm solution of an under-determined linear system. Theleast-squares minimum-norm solution for x, is:= vEub, (4-9)i.e., is the shortest x that minimizes IlAx — b112 and is unique. Equation (4-9) usesVTV= I and UTU = I, and the pseudo-inverse I given by:0= ., (4-10)owhere:0, cY=O-,.+— 1 4-ll—, otherwiseJjOften in practice o 0, but is very small. The condition number of A, cond(A),where:64max(ojcond(A)=. (4-12)min(01)cY Odescribes this situation. If min(T) is small then cond(A) is large and the matrix is illcrOconditioned. Either the machine precision or the accuracy of the data determine themaximum tolerable condition number. Let e be the limiting precision, then the conditionnumber should not exceed YE’ i.e.:max(oj /cond(A) = (4-13)mln(0.)cTjOyields a well-conditioned problem.SVD provides the opportunity to condition the problem by throwing out any cthat is too small and replacing it with a zero. Consequently, the null space increases toinclude precisely the part of the problem that was ill-conditioned at the cost of a loss ofrange. Redefinition of the pseudo inverse to condition the matrix gives:0, <Emax(cT)÷— 1——, otherwise -Jj4.1.3 Application to Limited-Angle CTComputation of the tomographic reconstructions proposed in Chapter 3 usingSVD follows these steps:rR1. build the matrix A= [2. use a canned SVD algorithm, e.g. svdcmp () from [31], to compute U, V,and the diagonal of .3. set1, cY<Emax(J)—otherwisefor j=l,2,...,n.654. use a canned SVD back-substitution algorithm, e.g. svbksb () from [31],to compute the least-squares minimum-norm solution, k, from U, V and ,anddatab=[WXFRepeat step 4 as needed for different data provided A does not change.SVD has O(mn2+ n3) complexity in time. Current computer technology is tooslow to tackle CT problems of a practical size with SVD and third-order complexitystaves off the day when it will. Storage requirements are also exorbitant. In theexperiments reported in Chapter 5 the image size is 72 x 200 giving n = 14,400. Thematrix A must contain at least n2 = 207,360,000 elements. Using a 32-bit floating pointrepresentation for each element, A occupies about 790 Mbytes of memory. SVD requirestwo arrays of this size. Arrays of this size are not practical using current technology andthis is only a moderately sized problem. Images of n = 200 x 200=40,000 are routine inCT.It may be possible to capitalize on sparsity in A and reduce both computation timeand storage. Unfortunately, sparsity in A does not translate directly to sparsity in U andV. Currently no efficient methods exist for SVD of sparse systems. Should suchtechniques come to fruition, SVD has the advantage that once a decomposition iscomplete, solutions are easily computed for different data vectors b, i.e., do thedecomposition once to compute U, V and I, taking a long time if necessary, andrepeatedly solve for any number of different b’s quickly.To summarize, SVD has the following advantages:1. definition of range and null space,2. controlled conditioning of the problem, and3. well-posedness due to unique least-squares minimum-norm,and the following disadvantages:1. third-order computational complexity, and2. exorbitant storage requirements.66If one is willing to sacrifice resolution of reconstruction SVD can perform CTreconstruction. The resolution will be too low for practical applications, but thedecomposition gives a useful analysis of the problem.4.2 Regularization and Conjugate Gradient Method4.2.1 RegularizationAs shown in the previous section, SVD gives a least-squares solution toAx = b,i.e., it minimizes the quadratic objective function E = IlAx— b112. The minimum is notunique because A has a non-trivial null space, so SVD fmds the shortest x that minimizesE = lAx— b112 (least-squares minimum-norm solution).An alternative is to modify the system so that the least-squares solution is uniquewithout the minimum-norm constraint. Tikhonov and Arsenin [41] describe one suchtechnique, called regularization, to solve singular and ill-conditioned systems of linearalgebraic equations. This section examines regularization in the context of the proposedlimited-angle CT system.Begin with the linear system:Ax=b. (4-15)To have a unique least squares solution it is necessary for A to be full column rank.When A is not full column rank, regularization modifies A by adding more rows to give:rAl rb 1[aQj[aQx0j’rAiwhere a > 0, and Q is some matrix. If Q is full column rank, then the matrix I I is[a2Jalso full column rank and the least squares problem of Equation (4-16) is well-posed. IfrA1Q is nearly full column rank, then it may be sufficient to make [Qj full column rank.To see that this is indeed the regularization proposed by Tikhonov and Arsenin, look at67the quadratic objective function minimized in the least-squares solution of Equation (4-16):EA 1 Eb 12E= [aQj’La2x0j‘ (4-17)which expands toE IlAx— b112 +a2llQ(x — xo)112. (4-18)Compare Equation (4-18) to Tikhonov and Arsenin’s functionals:= lIAz—ü112+ allz—z’2,a >0, (4-19)andMa [z, ü, A] = NAz — jill2 + aQ[z], a > 0. (4-20)The addition of the regularizing functional (the last terms in Equations (4-19) and (4-20))is equivalent to adding more rows to Equation (415)*.Now consider the elements of the regularizing functional, a, £2 and x0. If a =0Equation (4-16) reduces to the original singular system of equations, Equation (4-15). Onthe other hand, if a —> 00 Equation (4-16) reduces to:Qx = £2x0,which is not the problem of interest. So a varies through a continuum. At one end asmall a gives a problem that is close to the original but is ill-conditioned, while at theother end, a large a gives a different problem but one that is well-conditionedVarying a makes a tradeoff between conditioning and fidelity of solution to the data.We now examine two options for selecting £2. The first option is £2 = I, theEA1identity matrix, which is full column rank so I I is also full column rank. This leads to[al]the following variations of the problem:* Tikhonov and Arsenin use a where I use a2.** In this case ill-posedness and ill-conditioning are different degrees of the same property. a = 0 gives anill-posed problem, whereas a> 0 but very small gives a well-posed problem that is ill-conditioned.681) if x0 = 0 then IIQ(x — x0 )112 = IIxN2. Regularizing in this case gives a wellposed problem by biasing the solution towards the origin.2) if x0 = , where x is a vector of i’s and i is the average value of x, thenIIQ(x — x0 )112 = Bx — JI2. Regularizing in this case gives a well posedproblem by biasing toward the mean. This approach appears in the CTliterature [1] [21] [22].ED1The second option is Q= ‘where D1 is a discrete partial derivative in the xLD2Jaxis direction in the reconstruction and is a discrete partial derivative in the y -axisrol rD1ldirection. In this case we consider only x0= [ j. The matrix [ j is not full column0rank but is nearly full column rank. To see this consider the corresponding system ofpartial differential equations:df(x,y)—odxdf(x,y)—dyfor which the solution is:f(x,y) = c,where c is a constant. There is only one degree of freedom in the solution, that is theED1selection of the constant. If A can resolve the single degree of freedom in thenLD2]AI3 is full column rank and least squares solution of Equation (4-16) is well-posed.D2rD1l rolNote that in this case, if a is too large then Equation (4-16) reduces to [ jx = [ j and0Ax = b diminishes from the problem. The single degree of freedom in the solution is notresolved and the least-squares problem is ill-conditioned.The first option, Q = I, biases the solution towards x0. Therefore, when x0 =0the solution is drawn towards a vector of zeros. We know a priori that such a bias is not69warranted. Likewise, when x0 = the solution is biased towards the mean. This bias,although more sensible than the previous one, is still not warranted by knowledge of thespecimen. Values of x are known to be of certain values only and the bias drawssolutions to some value known not to exist in the specimen (except by chance). Althoughit is possible to use the fusion data to avoid some of this problem by making £2 I— W,the problem still exists for the interior regions of the reconstruction.EDi roiThe second option, £2=and x0=, biases the solution to smooth images.LD2J L°JAlthough we know that the cross section is not perfectly smooth, it is smooth over mostlocal regions, i.e., the specimen is piece-wise constant. Although there is a penaltyincurred by smoothing local regions containing discontinuities, this option is morepalatable than the first. Experiments reported in Chapter 5 use the second option for thelimited-angle CT problem.4.2.2 Parameter SelectionSelecting a value for the parameter a is an important consideration. If a is toosmall then the problem is ill-conditioned, but if a is too large then the solution deviatesgreatly from the data. In general, we want to have a as small as possible while stillhaving a well-conditioned problem. Experiments presented in Chapter 5 use a trial-and-error approach to select a. Since the conjugate gradient method solves the problemquickly (see Section 4.2.3), it is reasonable to compute solutions for a range of a. Thensimply use the smallest value of a for which CG converges to a reasonable solution. Areasonable solution is not so smooth that important detail is invisible, but does not haveconfounding artifacts due to ill-conditioning. Admittedly, the process is very subjective.Cross validation [13] [14] is an alternative technique for selection of anotexplored in this thesis. If, for example, £2 = I and x0 = 0 then cross validation selectsa2 by minimizing the cross-validation weighted square error [13]:C(a2)=wk[akxk(a2)_bk]2.m k=170Paraphrasing Golub and Van Loan [13], minimization of C(a2) is tantamount tochoosing a2 such that the final solution does is not overly dependent on any singlemeasurement. Intuitively, cross validation minimizes the sensitivity to the data. Crossvalidation is possible for other Qs and x0s; see Bates and Wahba [3] for example.Although the trial-and-error method is satisfactory for the work presented in this thesis,cross validation may warrant future consideration should selection of the regularizationparameter become critical, bearing in mind that cross validation has relatively highcomputational costs [3].4.2.3 The Conjugate Gradient MethodHaving regularized the problem, the resulting symmetric positive definite systemrequires solution. The conjugate gradient method (CG), recommended by Artzy, Elfving,and Herman [1] for fast convergence in CT, is one of many possible methods for solvingthe system. This section presents a summary of CG. For more details see Axelsson andBarker [2], Press et at.. [31], or Pierre [29].The following steps constitute the CG algorithm for finding the least-squaressolution of Ax = b, where A is positive-definite [11 [31]:(1) g0=hb—Ax(2) s1 = Ah= (g,,h1)(h,s)(3) x÷1=x+Ah(4) (direct method) g11 = b — Ax1÷(indirect method) g11=g. —5 — (g1÷,+) r — ((g,÷1 — g),g+1()— (g,g)°— (gg)(6) h11 = g1 + y1h(7) go to step (2) for next iteration.71Note that in the context of CG, A is not[R],instead A = RTR + WTW +a2QTQ (seeSection 4.2.4). Steps (1) and (4) initialize and compute the sequence g0,g1 ,g2,...,g_1,where g. is the negated gradient of E at x, i.e., g. —VE(x1). With initial conditionsdefined in step (1), steps (5) and (6) produce a sequence of mutually A-conjugate vectorsh0,h12...,h i.e.:hAh1=O, ij,where the h’s are non-trivial. Because the h’s are mutually A-conjugate they arenecessarily linearly independent. Therefore, to minimize E it is sufficient to do nseparate minimizations along the independent directionsh0,h12...,h_ (steps (2) and(3)).Each iteration requires multiplication by the matrix A, which has 0(n2)computational complexity. The complete minimization requires n iterations for a netcomplexity of 0(n3). Two factors make 0(n3) pessimistic though. First, often A issparse and multiplication by A can be faster than 0(n2). Second, requiring n iterations isa worst case scenario. Axelsson and Barker [2] give the more precise bound:p(e) .-gc(A) 1n)+ 1,where p(e) is the number of iterations required to converge to precision 8. ic(A) is thespectral condition number of A defined as:Fc(A) = maxIininwhere A., and are the maximum and minimum eigenvalues of A. Clusteredeigenvalues lead to further reductions in the number of iterations. We have no measureof the eigenvalues for A, but Davison [10] and Louis [23] show that singular values for Rare clustered, which leads to clustered eigenvalues for RTR. Results in Chapter 5indicate that convergence is much better than the worst case, so it is likely that the72limited-angle CT problem here benefits from both the condition number and theclustering of eigenvalues.4.2.4 Application to Limited-Angle CTApplication of regularization to the limited-angle CT problem (with data fusion)results in the modified system of equations:A bD1 x= 0 , (4-21)D2 0or more precisely:R yW WXFx = . (4-22)D1 00Minimizing the objective function:E= IIRx—+ IIW(x — xF)II + + a2IIDx, (4-23)or equivalently solving the following symmetric positive-definite linear system:(RTR+WT+a2D1T+a2DT)x= RTy+ WTXF. (4-24)finds the least squares solution to Equation (4-22).To solve the limited-angle data fusion problem by regularization:1. Implement algorithms to compute the forward raysum operator, R, and theback projection operator, RT.2. Implement algorithms to compute W = WT and D1TDI + D2TD (a discreteLaplacian operator).3. Compute the vector RTy + WTXF.4. Use CG to solve Equation (4-24).In summary, the advantages of regularization and CG are that it is fast for limitedangle CT, and it gives a well-posed and well-conditioned problem. The disadvantage isthat regularization makes a tradeoff between fidelity of solution to the data and73conditioning of the problem. If one is willing to concede fidelity to achieve a well-posedand well-conditioned problem, regularization and CG offer a practical approach tolimited-angle CT.4.3 Projection onto Convex SetsProjection onto convex sets (POCS) is a flexible method, able to incorporate avariety of constraints. As such it is particularly useful for constraint-based data fusion.The following discussion of POCS is based on Brègman [7], Youla and Webb [44],Gubin et al. [17], and Bauschke and Borwein [4]. POCS is used extensively for CT andsubsumes the method known as ART (algebraic reconstruction technique) [8] [22] [42].See Censor and Herman [8] for a review of projection methods, or Bauschke andBorwein [4] for a detailed survey.4.3.1 General Description of POCSPOCS assumes a set of m constraints on the domain of A, all of the form:xeC,, (4-25)for 1 i m, where C, is a convex set. Combined, the constraints restrict the solution tolie in the convex set C0, i.e.:xeC0=(ThC. (4-26)POCS makes the important assumption that C0 is non-empty. Each convex set C, has acorresponding projection operator P, of the form:P.z= [point in c. closest to z, z C, (4-27)otherwiseApplication of P, to some point z finds the point in C, closest to z. Points within C, arefixed points of F,. Let 1’, = 1 + A1(P —1) be the relaxed operator corresponding to F,,with relaxation coefficient, A, in the range 0 A 2. The relaxed projection operators(and the projection operators for 2 = 1) are non-expansive mappings, i.e.:BT,x—T,yII lix—y.74The non-expansive property leads to the convergence properties of POCS.Consider the parallel combination of a set of non-expansive mappings:=yTwhere the y are weights, 0 y 1 and y =1. The sequence {x} = {T;x0}converges weakly to a point in C0. Furthermore, if x is in Euclidean space then theconvergence is strong. Note that if the y1 vary between iterations, then becomes thesequential combination of the T,:Tseq_TmTm_i••Tiby appropriate selection of the y. This property gives flexibility, guaranteeing that themethod will converge for both parallel and sequential variations, as well as for any otherallowed values of y1. For arbitrary convex constraint sets there is no guarantee of rate ofconvergence. However, for constraints such as hyperslabs, or when the sets overlap (seeYoula and Webb [44], Theorem 3) then convergence is at a linear rate*, although there isno way to know the value for /3. For the purposes of this research we set 2 = 1 always,and select the y to give sequential application of the projection operators. Therefore, wecompute the sequence {x} = {Px0}where P PmPmi••• ‘iFigure 4.1 illustrates POCS convergence. Two constraints require that thesolution lie in the intersection of two half planes. Arrows in Figure 4.1 trace the progressof {Px0} towards x. Convergence is strong, because x is in two-dimensional Euclideanspace, and linear, because of the constraint set overlap.* By linear convergence we mean II’c — xli a/3, where 0 /3 < 1 and a 0. In some referencesthis is called geometric convergence.75C.’xo_V __XIPX0 2C2Figure 4.1: Example of POCS convergence.\//Figure 4.2: Demonstration of the variability of POCSconvergence. POCS converges to three distinct but correctsolutions depending on the order of projection operators.It is tempting to relate the POCS solution to the minimum-norm solution of SVD.if x0 = 0 then P0x9,where I’ is the projection onto the intersection C0, is the minimumnorm solution. Despite this temptation, POCS (as described here) does not computeP0x but only guarantees that it finds a solution that lies within C0. Figure 4.2 illustratesC•276this. Each of three paths converges to a point common to three convex sets. Only theorder of application of the projection operators distinguishes between the paths. Eachpath leads to a different but correct solution, but only one solution is P0x,and this is bycoincidence only. The ultimate solution varies with the order of projection operators. Ifthe minimum-norm property is essential, then one may use Hildreth’s method [8] orBrègman’s method [7] [8].4.3.2 Constraint Sets for Limited-Angle CTApplication of POCS to the limited-angle CT problem with data fusion requiresspecification of convex constraints and corresponding projection operators. Theflexibility of POCS leads to diversity of constraint sets for CT. Three different constrainttypes are sufficient for the problem at hand. These are:1. fit to raysum data,2. fit to range and ultrasound data, and3. amplitude constraints.Equation (4-28), below, forms the basis of the first constraint type:Rx=y. (4-28)Consider each raysum measurement, y, independently. Each measurement constrainsthe solution to lie within the hyperplane defined by:(r,x)—y1=0, (4-29)where r is the th row of R, and (.,.) denotes an inner product. For m raysummeasurements there are m hyperplanes each of which is a convex set given by:CR ={x:(r1,x)_y =o}, (4-30)where the subscript R refers to raysum. ART is the POCS method that uses only this typeof constraint set. The corresponding projection operator is found in many sources datingback to Brègman [7] [8] [27] [28]. The projection operator for CR. is:77PRz= r,, z CR• (4-3 1)otherwiseThe constraints of Equation (4-30) require a precise fit to data. Because it is essential toPOCS that C0 be non-empty, such precise constraints are too restrictive. Considerinstead this modified linear constraint, a hypersiab:CR {x: (r,x) — yj ER}, (4-32)where the allowed margin for error in the fit to raysum data, ER, determines the thicknessof the slab. Equation (4-32) breaks down into the two intersecting half planes:CR ={x: (rl,x)—yIER}(4-33)CR = {x: y — (r,x) ER}The modified projection operators are:y. +ER —(r.,z)= + (r1 , r) r1, (r , z) — y > ER (434)z, otherwiseand— ER (r1,z)PR.bZZ + (r1,) r, y — (r,z) > ER (4-35)z, otherwiseEquations (4-34) and (4-35) are simple variations of Equation (4-3 1). The modifiedconstraint sets allow control over the margin of error, ER, to vary the size of the sets andensure that C0 is non-empty.There are two ways to incorporate the second constraint type, fit to fusion data.The first is analogous to the fit to raysum data. Fusion constraints originate from thelinear equation:WX=WXF, (4-36)where W is a diagonal matrix with elements on the diagonal equal to one or zero. Ignoreall the trivial rows of W because they provide no useful constraints. The remaining rowsdefine the convex sets:78CF = {x: w1(x,— xF ) E} (4-37)where the subscript F refers to fusion. The corresponding projection operator is:EF £Z:zj=X+, ZjXF>Wi WiPFZ= z:z, =XF——,—z >—. (4-38)Wi Wiz otherwiseThe projection operator above clips pixel values to lie within a region about the partialreference image XF, with the size of the region determined by 8F and w.A second approach to fusion constraints originates with Oskoui-Fard and Stark[27]. They use a full reference image to constrain limited-angle CT. The concept issimilar to the fusion problem here. Oskoui-Fard and Stark assume that they have acomplete reference image, i.e., W = I, and define the convex set:CF = {x : lix — XFII EF}. (4-39)The corresponding projection operator is:XF + EFz— XF, liz — XFII> EFPFZ= liz — XFII . (4-40)z otherwiseUnfortunately, the assumption of a full reference image is not practical and, in thelimited-angle CT problem at hand, the reference image from the range and ultrasounddata can only be partial. Redefinition of the constraint set for a partial reference gives:CF = {x : llW(x — xF)ll eF}, (4-41)with the corresponding projection operator:(I—W)z+WxF+EF W(z—xF) llW(z—XF)il>CFPFz = llwz — xF)ll . (4-42)z otherwiseFor reasons explained in Section 4.3.3 the latter form of constraint is best for the limitedangle CT fusion system.Amplitude constraints are the third and last constraint type. It is physicallyimpossible to have a negative linear attenuation coefficient. Also, one can usually79determine a maximum allowable linear attenuation based on knowledge of what materialsare in the object. Therefore, the amplitude of a pixel in a reconstruction must lie betweena lower and an upper bound, giving the following convex constraint sets:CA={x:axIb, Oa<b}, (4-43)where the subscript A refers to amplitude. The CA are also hyperslabs. Theircorresponding projection operator is:z <az > b (4-44)otherwiseAmplitude constraints amount to simply clipping the data so that pixel values lie betweenthe prescribed upper and lower bounds.4.3.3 Constraint and Parameter SelectionThe previous section describes two versions of data fusion constraints. Thedifference between the two is not trivial. To illustrate this, consider the problem:• Solve for x from Ax = b by minimizing the objective functionE = lAx— bIIwhere ll•L denotes an 1 norm. This problem is equivalent to the linear programmingproblem:• Minimize e subject to the constraints(a1,x)—bJE, i=l,2,...,m.The constraints above are identical to the raysum constraints from the previous section.This leads to the following observations:1. with the raysum constraints, POCS is solving something like the abovelinear programming problem for Rx=y, and2. the optimum value for 8R (i.e., CR as small as possible with C0 nonempty) gives the l solution to Rx=y subject, of course, to the otherconstraints in the problem.80It is possible to solve the linear progranmiing problem and find an optimal value for ER(e.g., Brègman [7]) but it is simpler instead to estimate what 8R should be based on theaccuracy of the data. If the estimate of 8R is too large, then POCS may converge to apoor solution. On the other hand, if 8R is too small, POCS converges to a region of goodsolutions and then oscillates within it. The following summarizes selection of ER:1. find the optimal value of 8R by solving the linear programming problem,2. instead estimate 8R based on known accuracy of the data, and3. if estimating 8R’ err on the side of ER too small to keep the ultimatesolution within a region of good solutions.Now consider the fusion data constraints. While the raysum constraints findsomething like an l solution to Rx=y, the fusion data constraints of the form:CF ={x:IIvv(x—xF)II LF},find an 12 solution to Wx = WXF. As with ER, it is easiest to estimate a value for EFbased on the accuracy of the data, while keeping in mind that it is best to err on the sideof EF too small.The difference between the l and 12 constraints is significant. While the l. normis sensitive to any deviation in the data, the 12 norm can effectively ignore a large localdeviation by global averaging. The significance emerges when the different data sourcesare not consistent. An fit to the raysum data forces consideration of all raysummeasurements without averaging out local deviations. Should the raysum data not agreewith the partial reference image, an 12 fit allows the reconstruction to deviate locally fromthe reference so long as the global fit remains. This fusion task demonstrates theflexibility of POCS in facilitating the combination of l, and 12 constraints in the sameproblem.4.3.4 Application to Limited-Angle CTLimited-angle CT reconstruction using POCS follows these steps:811. setx0=O.2. perform an iteration by computing:= ‘A1 ‘R,,.1 “‘RX0.3. repeat step 2 to compute the series {x } = {P’x0} until convergence.Advantages of the POCS method are:1. it is fast for the problem at hand,2. it incorporates varied constraints easily and blends well with constraint-based data fusion, and3. it has linear convergence under many conditions.Disadvantages of the POCS method are:1. its solution depends arbitrarily on the order of application of projectionoperators, and2. it requires care in the selection of constraint sets.In practice the POCS proves to be fast and versatile making it an excellent choice for CTwith data fusion.4.4 Chapter SummaryThis chapter presents three numerical algorithms to solve the proposed limited -angle CT problem with data fusion. Table 4.1 below summarizes each of the threemethods.The next chapter shows results from each of the three methods explored here.SVD, because of its computational complexity, is restricted to small sample problems. Itdoes, however, show exactly what happens to the limited-angle CT problem with theproposed data fusion. Regularization with CG and POCS are both adequate for practicalsized problem. Results in Chapter 5 show that sacrifices in accuracy made byregularization to achieve a well-conditioned problem are too great. The ill-posedness ofthe limited-angle problem is not easily handled by regularization when the data are less82than ideal. POCS proves to be the winner in spite of its ambiguity about the solution. Itis important to remember that because each of the numerical methods incorporatesdifferent constraints, they are, in essence, solving different problems.MethodProperty SVD CG POCSUniqueness minimum norm regularization not unique - solutiondepends on orderof operators andinitial conditionsConditioning eliminate small use regularizing enlarge constraint setssingular values functional to get C0 non-emptyAdvantages • computes null • fast (for the task at • fast (for the task atspace and range hand) hand)• can condition the • control over • incorporates a wideproblem conditioning variety of convexconstraints• gives well-posedproblem • converges at a linearrateDisadvantages • slow • makes a tradeoff • solution depends onbetween fidelity to order of operators• exorbitant storage data andconditioning• must ensure C0 isnon-emptyConstraints• minimize lAx— b112 ‘ = ‘O • amplitude limits2 ax.band lxN . .• minimizerA 1 [b 12I Ix—I[aQjOther • projects data onto • projects data onto • projects ontorange of A EA 1 constraint sets inrange of [2j the domain of ATable 4.1: Summary of numerical methods considered forlimited-angle computed tomography with data fusion.83Chapter 5Experimentation and ResultsThis chapter presents experimental results that establish the validity of theproposed limited-angle CT data fusion system. In addition, the experiments explore thecapacity of each of the three numerical methods described in Chapter 4 to deal withconfounding errors in the data.The experiments consist of a series of trials. Each trial applies one of the threereconstruction algorithms to some data to produce a reconstructed cross section foranalysis. The algorithm, its parameters, and the data distinguish trials. Section 5.1describes these distinguishing features, including salient features of the algorithmimplementations and the applicable parameters. The section then describes the sources ofdata, both synthetic and real. Synthetic data trials verify the validity of the proposedlimited-angle CT system, while the real data trials show practical application of thesystem. Errors in the real data highlight the capacity of each algorithm to deal withconfounding errors.84Section 5.2 presents the trial results organized by algorithm. The first three trialsfor each algorithm compute reconstructions based on varying fusion constraints. Thesetrials show a significant improvement in accuracy due to the fusion constraints. SVDshows, in addition to improved accuracy, the changing nature of the problem in terms ofthe number of singular values equal to zero.Whereas SVD is limited to small problems, regularization combined with theconjugate gradient method (RICG) and POCS can handle large problems and so areapplied to real data. Trials based on RJCG and POCS, and real data, show the success ofthe proposed limited-angle CT system with real data. They also show the capacity of thereconstruction algorithms to cope with erroneous data, and the response to varyingreconstruction parameters.The salient conclusion of the experimentation is that, when fused with limited-angle CT, constraints based on ultrasound face sheet data improve the accuracy ofreconstruction. Although measurement errors confound reconstruction of accurate crosssections from real data, the limited-angle CT system does, nevertheless, work in practice.Of the algorithms used, POCS proves to be superior because it produces accuratereconstructions, is less sensitive to its parameters, and does not force a compromisebetween image smoothness and conformity to data.5.1 Description of ExperimentsExperimentation presented in this thesis consists of a set of trials where each trialinvolves:1. acquisition of some data,2. application of a reconstruction algorithm (SYD, RICG, or POCS) to thedata, and3. analysis of the algorithm output.These items uniquely describe any trial:1. the algorithm,siuwpq3rqutJY3AtSioiudoujnsAiijooiqoi3j‘jopudotUflSJqijoaunbaiuon,nj1suoaIojspoqwpc.iurnuaiq1livJO11LdO11ltU3IUBUSipipuo‘sCi.ijojs(p)puc‘uopn.gsuoiouipsA.ipiiind()‘uopnhJsuoJjowj.ins.iddno11Ifl!PUd1d sAiu(q)‘uO!pflisuo.IoJqjodAq(i’):oidownsAiupudiuoiosdjg:j’ainj(p)(0)4,4 S•g4444.4.I•.——aa*a*——S4%4s*4%4444444..%SS-‘‘4%’’’’’(R)‘I,---..a--I----FI4I---—-%b.—-.(-sJI,?%11144(q)£.IJjJ.’Jliii...I•.ia.•aaIIIIIIII•UAIospwquopqoJojsuonipuoouiddospui‘suonipuoolu!‘sJIu1mJpssnuoniuuijdwiqoiqoiqMUTuuu.iqisoquospuqpu‘jowjodownsAioqJouoiuosaidarioidwooiuiquospAqsisuonossqJuiquouonoruisuooaiqojouoiiiuujdunoqoj.uusioidownsAii.ioqSUI1IJIJOIVJourspuionipu/soq‘wpjosoinospuiuupuopqojouonnuouioidwijosainwjuTioqsAqswiTsoqsquospuonoossTqinippui‘siuriud86the distance an individual ray passes through an element of the reconstruction. Thematrix required to represent all rows of the raysum operator completely is huge(207,360,000 elements for a 200 pixel by 72 pixel reconstruction). However, mostelements of each row are zero, i.e., R is sparse, and there are patterns in the operator thatlead to a compact representation.Figure 5.1(a) shows a hypothetical grid for reconstruction. Elements of theraysum operator rows are the distances a particular ray travels through each box in thegrid. A clipping algorithm (from computer graphics) computes these distances. Figure5.1(b) shows the set of rays corresponding to a scan across the specimen parallel tovertical. A scan at an angle to vertical measures raysums for the set of rays shown inFigure 5.1(c). Rays for a scan at a constant angle are parallel, and parallel raysums areidentical to each other except that the columns (of the reconstruction grid) are offsetcorresponding to the motion of the scan. With the exception of the column offsets, theset of rays in Figure 5.1(d), where there is one ray for each scan angle, completelydescribes all rays used in sampling. Therefore, to get a compact representation of theraysum operator:1. record only one row of the operator for each scan angle, and2. record only non-zero elements.The compact raysum operator is a list of vectors with one element in the list for each scanangle. Each vector is a list of non-zero distances as well as the row and columncoordinates for the corresponding grid elements. To reconstruct a row of the raysumoperator quickly from the compact representation:1. assume the row consists of all zeros,2. select the vector for the appropriate angle from the compactrepresentation,3. apply the appropriate column offset to the column coordinates of thevector, and874. put the offset elements of the compact representation into thereconstructed row.R has one row for each combination of scan angle and column offset, but excludes thosecombinations that produce rays passing beyond the width of the grid. Only SVD requiresthe full raysum operator representation. For other algorithms, the compact representationallows quick recovery of the elements needed for computation.Implementation of SVDrRlSVD requires an explicit representation of the matrix [wj• Implementation ofthe proposed limited-angle CT reconstruction follows these steps:1. declare matrices for A, V, and I (note that since I is diagonal only theelements on the diagonal are represented),2. put rows of R into rows of A from the compact raysum operator asdescribed above,3. put rows of W into rows of A,4. if A has fewer rows than columns add sufficient trivial rows to make Asquare,5. call SVD subroutine, svdcmp () [31],6. set non-zero singular values that are too small to zero, i.e., condition thematrix,7. create the vector b= V 1 and[Wx F]8. compute the reconstruction from the decomposition and b withsvbksb() [31].Because svdcmp () specifically deals with singular matrices it does not require anyspecial precautions.88Implementation of R/CGSection 4.2.3 gives the conjugate gradient method for finding the least squaressolution to the positive definite system:Ax = b.Two things are necessary to implement the method:1. an algorithm to compute Ax, and2. an algorithm to compute b.From Equation (4-24) we have:Ax = (RTR+ WTW+a2DrD1+a2DTD)x, andb = RTy + WTXF.Break the computation of Ax down into three separate computations: RTRx, WTWx,and (a2D1TD +a2DD)x.To compute RTRx, first compute Rx using the compact raysum operator. Rx isa column vector, each element of which is the inner product of a row of R and x.Computation of the inner products in Rx is straight forward from the compactrepresentation of R.RT is the back summation operator. The ii” element of the vector RTb is:[RTb].i.e., [RTbJ, is the weighted sum of all the raysums for rays that pass through thecorresponding grid element. Computation on this basis is inefficient, but a more efficientapproach arises from observing that each raysum contributes to a set of grid elements inthe back summation. That set of elements is the same as the set that contributed to theraysum. The weight applied to a raysum for its contribution to an element of the backsummation is the same as the weight for that element in the original raysum. Thus theback summation can be computed in a manner analogous to the raysum. For theexperimentation here, the vector RTRx is computed by a two step process: first Rx andthen RT(Rx) . The compact raysum operator allows fast computation of both steps.89The matrix W is diagonal so WT = W and WTW = diag(w, wy,. .., w,). From thisit is easy compute the vector WTWx.+a2DD is a discrete Laplacian operator. Convolution with the kernelshown in Figure 5.2 best summarizes its implementation. Should different a2 for the xand y directions be necessary, the more general kernel of Figure 5.3 applies.a2a2 4a2 a2Figure 5.2: Convolution kernel for Laplacian operator.-a2 La 2—a1 2—a1-aFigure 5.3: Convolution kernel for Laplacian operator withdifferent coefficients for x and y directions.Using RTRx, WTWx and (a2DrD1+a2DD)x, compute Ax from:Ax = (RTR+ WTW+a2DrD1+a2DD)x= RTx + WWx +(a2DD1+a2DD)xComputation of b is similar to Ax; the back summation operator gives RTy andthe identity WT = W gives WTx. Add these two vectors to get b = RTy + WTXF.Initial and stopping conditions are important to an iterative procedure such as CG.Although any arbitrary starting condition, x0 , is acceptable, these experiments use= WXF (or x0 =0 for trials omitting fusion data) to make maximum use of the fusiondata to avoid unnecessary iterations. Iteration stops when the norm of the residual,IlAx — bI, falls below a specified threshold.90Implementation of POCSPOCS computes the sequence {; } = {Tx0} where T = Tm Tmi ,..., T and=1 + — 1). There is no known method for determining optimum values for therelaxation coefficients A, so experiments here use = 1 for all i. This reduces to thesequence {x} = {Px0} where P= and P is the projection operator ontothe th constraint set.Three types of projection operators are necessary for the three different types ofconstraints. The first is the raysum constraints. Equations (4-34) and (4-35) give the twoprojection operators required for each row of the raysum operator R. The salient featureof these equations is the inner products (r1 , z) and (r1 , r,). The compact raysum operatorallows quick computation of these inner products.The second constraint type is the data fusion constraint. Experiments presentedhere use Equation (4-42) to implement the data fusion constraints. The computation ofthe projection is straightforward; it involves only computing the inner products of vectorsthat are directly available.The third, and last, constraint is the amplitude constraint. Its implementation istrivial. Simply clip all the elements of the solution to lie within the specified upper andlower bounds.As with conjugate gradient, it is important to specify starting and stoppingconditions. Unlike the conjugate gradient method, the starting point for POCS changesthe ultimate solution. Although POCS does not generally compute the minimum-normsolution, these experiments use x0 = 0 in order to arrive at a solution near the minimum-norm solution. This starting condition does not guarantee such an outcome, but workswell in practice. POCS iteration does not introduce spurious signals into the solution sox0 = 0 is a sensible starting point. Iteration terminates when the norm ix, — x-_1 fallsbelow a specified level, i.e., the solution stops moving (or slows sufficiently). Note that91unlike the conjugate gradient, this stopping condition does not guarantee that the solutionis good, just that the solution is not going anywhere.5.1.2 Synthetic DataThis section describes data synthesis for synthetic data trials. Data synthesis startswith an arbitrary cross section image as the subject of an experiment. Application of theraysum operator to the image yields the raysum data for the trial. The cross section alsodefines the spatial support data and the face sheet data directly.Application to aircraft parts provides the initial motivation for this research.Accordingly, aircraft parts, in particular sandwich structures with graphite/epoxycomposite face sheets around aluminum honeycomb core, are the basis for synthetic data.Thin aluminum in the honeycomb core pushes the resolution limits of the real dataacquisition apparatus (described in the next section). To avoid problems with resolution,the real-data experiments use a plexiglass phantom. In the phantom, two plexiglass facesheets surround a set of vertical plexiglass members simulating honeycomb. Allplexiglass parts are 3 mm thick which, with a reconstruction resolution of 0.5 mm by0.5 mm, translates to 6 pixels thick in the reconstruction. Real honeycomb is thinner than1 pixel. Synthetic data in these experiments simulate the conditions of the plexiglassphantom.Figure 5.4: Synthetic cross section image for SVD trials.Both time and numerical precision limit the SVD problem size to 10 pixels by 30pixels. Figure 5.4 shows the synthetic image for SVD trials. Two parallel face sheets and92honeycomb are 2 pixels thick. Within the honeycomb cells, two blocks mimic theplexiglass inserts of the real phantom.RJCG and POCS trials use 72 pixels by 200 pixels. At the specified resolution of0.5 mm by 0.5 mm this corresponds to a cross section 36 mm thick by 100 mm wide. Lowresolution restricts SVD trials to parallel face sheets, but RJCG and POCS trials have nosuch restriction. Non-parallel face sheets make the trials more interesting and they alsoquell possible doubts that results may rely on parallel face sheets. Figure 5.5 shows thesynthetic cross section for RICG and POCS synthetic data trials. Face sheets and verticalmembers are all 6 pixels thick.Figure 5.5: Synthetic cross section image for RICG and POCStrials.To conduct synthetic data trials it is necessary to synthesize data corresponding tothe three data sources. Application of the forward raysum operator gives the syntheticraysum data. All synthetic data trials compute raysum data this way using singleprecision floating point arithmetic. CT reconstruction is tolerant of unbiased noise on thedata, so trials with unbiased noise added are uninformative. In the absence of a bettermodel for measurement errors, synthetic data trials use error-free data.Range and ultrasound data give a map of regions that have known linearattenuation coefficients and the known values. Trials using raysum and spatial supportdata only assume regions external to the specimen are known and have linear attenuationcoefficients of zero. With the addition of face sheet data, trials assume external and facesheet regions have known linear attenuation values: zero for external and linearattenuation of plexiglass for face sheets. Figure 5.6 shows the synthetic fusion data as93images. Values from the image of known regions form the diagonai of W while thevalues of known linear attenuation coefficients form XF. For synthetic trials, theassumed linear attenuation of plexiglass is 0.40 cm’.Figure 5.6: Synthetic fusion data vectors: for spatial supportdata only (a) known region and (b) known linear attenuation,and for spatial support and face sheet data (c) known regionand (d) known linear attenuation.(a)(b)(d)945.1.3 RealDataReal data come from measurements of a plexiglass phantom simulating asandwich structure with graphite/epoxy composite face sheets and aluminum honeycombcore. As mentioned in the previous section, the thin honeycomb core pushes beyond theresolution of the experimental apparatus. So for the sake of experimentation, componentsof the phantom are thicker than what is actually found in real aluminum structures.Figure 5.7 shows schematically the cross section of the plexiglass phantom.All dimensions in mm3 12ITop and bottom face sheets (horizontal components)[El Simulated honeycomb (vertical components)Simulated entrapped water (defect)Figure 5.7: Cross section of plexiglass phantom simulating asandwich structure with graphite/epoxy composite face sheetsand aluminum honeycomb core.Figure 5.8 shows the apparatus for raysum data acquisition. The apparatusoriginates from a real-time radiography system custom-built by Philips for inspection ofaircraft control surfaces. The x-ray tube and linear array are mounted on a C-armmanipulator system (see section 3.2.1). The x-ray tube on one end of the C-arm projectsa cone of radiation through the specimen to a sensor on the other end of the arm. The Cshape allows the source and sensor to maneuver about and inspect various aircraft parts.95Linear arrayx-ray detector________12 bits 8 bits_/LXXW14Iwith frame grabberrLl_____I I 8-bit datax-ray tubeFigure 5.8: Apparatus for limited-angle CT raysum dataacquisition system.In the original configuration, the sensor was a real-time x-ray image intensifierand a video monitor displayed the radiograph in a shielded control room. The imageintensifier has since been replaced by a Thompson linear array detector. Whereas theintensifier acquires images of an area, the linear array acquires a single line of data. Toacquire an image the C-arm sweeps the array across an area. A frame buffer assemblesthe acquired lines of data into an image. The array is 1024 elements across with an inter-element spacing of 0.2 mm and a width of 0.5 mm. Analog to digital converters digitizethe intensity data to 12 bits but at present the frame buffer uses only the eight mostsignificant bits. The frame buffer produces an RS-170 video signal to display images ona video monitor. A PC-based video frame grabber acquires digital data by sampling theframe buffer RS-170 output.The experimental trials use the following procedure to acquire raysum data:1. sweep the C-arm across the specimen at a constant angle to the specimenand acquire data at 0.5 mm interval spacing,2. with the frame buffer, assemble the data into an image,3. digitize the image with the PC-based frame grabber,964. compute the average of the centre column of the digitized image and thetwo columns adjacent to it (to compensate for anomalies created by thelinear array), and5. use the data from the averaged column as a set of parallel raysums.In summary, the x-ray tube projects x-rays onto the linear array. The array createsan electrical signal from the x-rays and digitizes the signal to twelve bits. Eight bits ofthe twelve go to frame buffer which in turn produces an RS-170 video image. A PC-based frame grabber digitizes the image. Three adjacent columns of the digitized imageform a set of parallel raysums.There is a variety of error sources inherent in the apparatus. There is no filtrationof the x-ray radiation other than that inherent in the manufacture of the tube, so beamhardening is a problem. The linear array requires only a fan of radiation, so any excessradiation contributes to scatter.Raysum measurements are not uniformly sensitive to x-ray intensitymeasurements. Equation (2-3) describes the raysum as a function of intensity. For ahomogenous specimen Equation (2-3) becomes:Iut = —ln(—),10where I is the measured intensity, I is the initial intensity, u is the linear attenuationcoefficient, and t is the material thickness. The sensitivity of jit to I is:d(#t) — 1dl— I’which is a function of I. At I 225 the sensitivity is - 1/225 = -0.0044. A one quantumchange in I gives a change in ut of 0.0044. For plexiglass (u = 0.45 cm1 at 60 kV) thiscorresponds to a change in material thickness of 0.01 cm. However, at an intensity of 50,the sensitivity is -0.02, corresponding to a 0.044cm change in plexiglass thickness, aboutthe resolution of the reconstructions. When I is low, i.e., the specimen is thick, the dataare more sensitive to intensity errors. The sensitivity of the measurements to a quantum97change in intensity can be of the order of the reconstruction resolution. Certainly thisdata acquisition system needs improvement, but it is suitable for verifying the proposedlimited-angle CT system.The region reconstructed in experiments is 72 rows by 200 columns which, at aresolution of 0.5 mm by 0.5 mm, corresponds to a specimen 36 mm thick. The plexiglassphantom is 36 mm thick so no pixels in the reconstruction are external to the specimen.The region of reconstruction matches precisely the spatial support of the specimen and,for this experiment only, it is not necessary to acquire range data.The design of the phantom allows face sheet thickness measurement directly witha micrometer. Therefore, ultrasound measurement is not necessary. * The 3 mm facesheet thickness (as measured by the micrometer) translates to 6 pixels. So for real datatrials, the region of known linear attenuation consists of the top and bottom six rows ofthe reconstruction. The linear attenuation for these rows is that of plexiglass. At the xray tube setting of 60 kV used for these experiments, the linear attenuation coefficient ofplexiglass is 0.45 cm1.5.2 ResultsThis section presents experimental results. The experiments consist of a series oftrials wherein each trial involves application of one numerical method to one data set toproduce a reconstruction. The following three items distinguish trials:1. the algorithm,2. the parameters for the algorithm, and3. the data.For each numerical method, there are synthetic data trials followed by real datatrials. The first three synthetic data trials for each algorithm produce reconstructionsfrom:* In fact, if one were to use ultrasound to measure face sheet thickness, the micrometer measurementswould be used to calibrate the ultrasound measurements.981. raysum data only,2. raysum and spatial support data, and3. raysum, spatial support and face sheet data.Results of these three trials show the improvement in reconstruction accuracy with theincorporation of constraints.Real data trials show (for RJCG and POCS):1. the success of the proposed limited-angle CT with real data,2. the ability of the algorithms to cope with errors in the data, and3. the effect of parameters.Real data are more effective than synthetic data for items two and three above. Syntheticdata lead too easily to good reconstructions for a wide range of parameter values, evenwith Gaussian noise added. However, errors in the real data make the effects of theparameters obvious.An error measure is essential for comparing the accuracy of the different methods.Experimental results presented here use the measure, e, defined by [28]:e=X_x1OO%, (5-1)NxIIwhere x is the true cross section and is the reconstructed cross section.5.2.1 SVD TrialsLimitations of the SVD algorithm restrict SVD experimentation to trials based onthe low resolution cross section of Figure 5.4. Raysum data for these trials cover sevenscan angles from 6O0 to 600 in steps of 20°. Spatial support constraints assume the topand bottom rows of the reconstruction are known to be zero. Face sheet constraintsassume the second and third rows from the top and bottom are known to be plexiglass.As for all the synthetic data trials, SVD synthetic data trials give results forvarying levels of constraints. SVD gives a compelling illustration of the change in thereconstruction problem with added constraints by recording:991. the number of raysum rows in A,2. the number of spatial support rows in A,3. the number of face sheet rows in A,4. the number of singular values equal to zero, and5. the error measure.Figure 5.9 shows the reconstructions for each of the three trials and the correspondingabsolute error vectors. Grey scales in the images are set to allow direct comparison of thereconstructions and the errors, but the maximum and minimum values (listed in thecaption) differ for each trial. Figure 5.9(g) shows the absolute error vector for trial 3rescaled to show more detail. Table 5.1 summarizes the SVD trial results.Trial Description Rows of A Number of ErrorNumber Singularities Measuree(%)1 raysums only • 152 raysum 148 50.7• 148 trivial2 raysums and • 152 raysum 88 31.5spatial support • 60 spatial support• 88 trivial3 raysums, • 152 raysum 13 4.8spatial support, and • 60 spatial supportface sheets • 120 face sheet• 0 trivialTable 5.1: SVD trial results.Results in Table 5.1 show how spatial-support and face-sheet constraints add rowsto the matrix A. The number of singularities (singular values equal to zero) in A reflectsthe degree to which the problem is ill-posed. Reliance on raysum data alone is thepoorest scenario with 148 singularities out of 300 singular values. Spatial supportconstraints improve upon the raysum data, reducing the number of singularities to 88, butstill leaving the problem ill-posed. Incorporation of face-sheet constraints further reduces100the number of singularities to 13. Although the use of data fusion still results in an ifiposed problem, it markedly reduces the number of singularities and size of the null space.Figure 5.9: Reconstructed images for SVD trials: (a) trial 1reconstruction (max = 0.48, mm = -0.01), (b) trial 1 error vector(max = 0.30, mm = 0.0), (c) trial 2 reconstruction (max = 0.55,mm = 0.0), and (d) trial 2 error vector (max = 0.23, mm = 0.0).Continued on next page.(a)(c)(d)101(g)Figure 5.9 (continued): Reconstructed images for SVD trials:(e) trial 3 reconstruction (max = 0.42, mm = -0.04), (1) trial 3error vector (max = 0.04, mm = 0.0), and (g) rescaled trial 3error vector (max = 0.04, mm = 0.0).Not only does the data fusion reduce the size of the null space, it improvesaccuracy because the reduced null space does not contain the face sheet structure.However, even after addition of the fusion constraints the null space is not trivial and it ispossible to have a specimen with components in the remaining null space. Theminimum-norm solution would not be accurate because it omits the null-spacecomponents. Such components would be wide horizontal edges in the interior of thespecimen.(e)102From data that are readily available in trials 1 through 3 it is possible to compute aClass I wealdy-coupled fusion reconstruction, from:= WXF + (I — W)kwhere x is the solution to Rx=y, i.e., the reconstruction in Figure 5.9(a). Intuitively,the weekly-coupled reconstruction is the partial reference image pasted onto of theunconstrained reconstruction. Figure 5.10 shows this weakly-coupled reconstruction andits absolute error vector. The error measure for the reconstruction is 13.9% (versus 4.8%for the strongly-coupled reconstruction). Clearly, the strongly-coupled constrained-basedsystem produces superior results.5.2.2 RICG TrialsAll RICG trials compute 72 pixel by 200 pixel reconstructions. The first threetrials use synthetic data based on Figure 5.5, while the remaining trials use real data frommeasurements of the plexiglass phantom of Figure 5.7. Both synthetic and real data trialsuse raysums for 13 scan angles from -60° to 60° in steps of 10°.(a)(b)Figure 5.10: (a) Weakly-coupled reconstruction using SVD(max = 0.48, mm = -0.01), and (b) the corresponding absoluteerror vector (max = 0.15, mm = 0.0).103Figure 5.11 shows the reconstructions for the first three trials and thecorresponding absolute error vectors. As for all synthetic data trials, the first three usevarying levels of constraints. The parameter a2 is set to a constant value, a2 = 0.001.The stopping condition is lAx— bil < 0.1. Table 5.2 summarizes the RJCG synthetic datatrials.]Iiii ii_(c)—0.10Figure 5.11: Reconstructed images for RICG synthetic data trials: (a) trial 1reconstruction (max = 0.49, mm = -0.06), (b) trial 1 error vector (max = 0.40,mm = 0.0), (c) trial 2 reconstruction (max = 0.55, mm = -0.05), and (d) trial 2error vector (max = 0.43, mm = 0.0). Continued on next page.0.55104(f)Figure 5.11 (continued): Reconstructed images for RICG syntheticdata trials: (e) trial 3 reconstruction (max = 0.45, mm = -0.07), and (1)trial 3 error vector.Trial Description Regularization Iterations ErrorNumber Parameter Measurea2 e(%)1 raysums only 0.00 1 10 64.02 raysums and 0.001 12 55.3spatial support3 raysums, 0.001 13 6.7spatial support, andface sheetsTable 5.2: RICG synthetic data trial results.The images of Figure 5.11 show subjectively the improvement in accuracyattained by incorporating fusion constraints. Error measures in Table 5.2 show thisimprovement quantitatively. As was the case for SVD results, with RICG too there is agreat improvement in reconstruction accuracy for trials using data fusion.(e)I105In trial three, the error measure for RICG is comparable to that for SVD (6.7%versus 4.8%). The small difference is not significant and may be due to the difference inthe data and problem size. Figure 5.11 exhibits an oscillatory pattern in the interior of thereconstruction. The reason for its existence is well known. Finite precision means thatthe synthetic data do not precisely match the forward model. Consequently, RICO usesER1part of the null space of [wj (an oscillatory pattern) to match the data. A larger valuefor cx2 prevents this, but at the expense of smoothing the solution. The reconstructions ofFigure 5.11 are about optimum, based on trial-and-error selection of a2.Figure 5.12 shows the reconstruction of the synthetic cross section using weakly -coupled fusion. The error measure for the weakly-coupled reconstruction is 19.6%(versus 6.7% for the strongly coupled reconstruction), which indicates again thesuperiority of the strongly-coupled method.(b)Figure 5.12: (a) Weakly-coupled reconstruction using RICG (max =0.49, mm = -0.05), and (b) the corresponding absolute error vector(max = 0.22, mm = 0.0).The remaining RICO trials use real data to show:1. the success of the method with real data,(a)1062. the capacity of RICG to handle errors in the data, and3. effects of a2 on the quality of reconstruction.There are seven real data trials, designated 4(a) through 4(g). Trials 4(a) through 4(f) areidentical except that the parameter a2 varies from 0.001 to 100.0 in decade steps. Trial4(g) reconstructs the image using a different a2 for the x and y directions. All trials use astopping condition of IAx — bil < 0.1. A smaller tolerance does not significantly affect thereconstructions. Figure 5.13 shows the reconstructed images for these trials on the samegrey scale to allow for easy comparison. Table 5.3 summarizes the results.Trial Description Regularization Iterations Maximum MinimumNumber Parameter ( cmj ( cmja24 (a) • real data 0.001 20 0.71 -0.30• all constraints4 (b) • real data 0.01 20 0.71 -0.27• all constraints4 (c) • real data 0.1 25 0.68 -0.20• all constraints4(d) • real data 1.0 40 0.62 -0.17• all constraints4 (e) • real data 10.0 55 0.50 -0.05• all constraints4(f) • realdata 100.0 102 0.38 0.12• all constraints4 (g) • real data a12 0.001 24 0.68 -0.22• all constraints a2 =0.1Table 5.3: RJCG real data trial results.Reconstructions in Figure 5.13, in particular 5.13(c) and 5.13(g), resemble thecross section of the phantom in Figure 5.7. It is fair to conclude that the proposed fusionmethod for limited-angle computed tomography is successful with real data. One is107forced to admit, however, that at least for RJCG, the quality of reconstruction from realdata is poorer than from synthetic data.Figure 5.13: Reconstructed images for R!CG real data trials:(a) a2 = 0.001, (‘o) a2 = 0.01, (c) a2 = 0.1, and (d) a2 = 1.0.See Table 5.3 for maxima and minima. Continued on next(a)(b)(c)(d)page.108(g)Figure 5.13 (continued): Reconstructed images for RJCG realdata trials: (e) a2 =10.0, (f) a2 = 100.0, and (g) a12 = 0.001and a2 = 0.1. See Table 5.3 for maxima and minima.RICG does not appear to handle errors in the data well. In the synthetic case, thery 1 rR1vector b = I I is, within the limits of numerical precision in the range of I I. Such[WxF] [WJis not the case with the real data. Consequently, the errors exaggerate the effects of illconditioning and larger null space components appear in the reconstruction.The parameter a2 does not help deal with errors in this case. It merely allows oneto select the degree of compromise in the solution. Figure 5.13(a), for example, showsone extreme with a2 = 0.001. The image is sharp, but oscillations are severe asindicated by the maximum and minimum image values (they should be 0.45 cm’ and(e)(f)1090.00 cm’ respectively). At the other extreme, Figure 5.13(f) shows the reconstruction fora2 = 100.0. In this case the solution is so heavily biased towards a flat image that almostall structural detail is lost. Certainly this would not work with real honeycomb which ismuch thinner.The ill-posed nature of this limited-angle CT problem requires more smoothing inthe y direction than for x, so one can use a different a2 for each direction. Figure 5.13(g)shows the reconstruction for a12 = 0.001 and a2 = 0.1. It supports the conclusion thatthe proposed method works, but none of the images of Figure 5.13 is particularlysatisfying.Aside from selecting the degree of compromise in the solution, a2 affects thenumber of iterations required. Smoother reconstmctions require more iterations.Reasonable solutions like those for trials 4(c) and 4(g) take about 25 iterations and do notimpose an outrageous computational burden.5.2.3 POCS TrialsThere is a total of seven POCS trials. Three synthetic data trials show theimprovement in reconstruction accuracy due to fusion constraints, and a fourth shows theeffect of changing the stopping condition. A real data trial shows the success of theproposed limited angle CT method with real data. Two more real data trials illustrate theeffect of altering the parameter 8RFigure 5.14 shows the reconstructions for the four POCS synthetic data trials andthe corresponding absolute error vectors. Because POCS employs amplitude constraints,the reconstructions all have a maximum value of 0.40 and a minimum value of 0.0. Table5.4 summarizes the results. The first three synthetic data trials do reconstructions basedon varying degrees of fusion constraints (as done for SVD and RJCG). The fourth trialshows the effect of changing the stopping condition for the algorithm.110Figure 5.14: Reconstructed images for POCS synthetic data trials:(a) trial 1 reconstruction, (b) trial 1 error vector (max = 0.4, mm =0.0), (c) trial 2 reconstruction, and (d) trial 2 error vector (max = 0.4,mm = 0.0). Continued on next page.(a)(b)(c)(d)111tEtliri_(g)I0’)Figure 5.14 (continued): Reconstructed images for POCS syntheticdata trials: (e) trial 3 reconstruction, (f) trial 3 error vector (max =0.22, mm = 0.0), (g) trial 4 reconstruction, and (h) trial 4 error vector(max = 0.25, mm = 0.0).0.00(e)‘**t ‘.0)112Trial Description Parameters Iterations ErrorNumber Measuree(%)1 raysums only LR—O.OO1 13 62.6£F°1stop at11x — xi_i II < 0.12 raysums and 8R=0.001 33 38.9spatial support 8FO.lstop at— x1_ < 0.13 raysums, LR—.0.001 9 6.0spatial support, and eF=O.lface sheets stop at— x_1 < 0.14 raysums, ER—0.001 146 5.5spatial support, and EF=Olface sheets stop at11x — xiiII < 0.001Table 5.4: POCS synthetic data trial results.Trials 1 through 3 exhibit the same marked improvement in accuracy ofreconstruction observed for SVD and RJCG. Trial 3 shows that POCS computes areconstruction comparable in accuracy to SVD and RICG (6.0% for POCS versus 4.8%for SVD and 6.7% for RICG). POCS, however, does not make a compromise towards aflat solution, while it is on par with RICG computationally (i.e., it converges in about thesame number of iterations with about the same amount of work per iteration). Of thethree methods, only POCS manages to provide an accurate solution quickly withoutsmoothing.Only nine iterations are necessary for trial 3. This is not a lot of iterations, whichcompels one to ask whether or not stricter stopping conditions will improve the accuracyfurther. Trial 4 computes a reconstruction based on the stricter stopping condition113— x,_1 <0.001. The result is a large increase in the number of iterations (from 9 to146) as expected, but only a slight improvement in accuracy (6.0% down to 5.5%).Obviously the stricter stopping condition has not paid off in a significant improvement inaccuracy.Figure 5.15 shows the reconstruction of the synthetic cross section using weakly -coupled fusion. As is the case for both SVD and RICG, the weakly-coupledreconstruction is poorer than the strongly-coupled one (error measures of 18.0% versus6.0%), showing the superiority of the strongly-coupled method.(b)Figure 5.15: (a) Weakly-coupled reconstruction using POCS, and (b)the corresponding absolute error vector (max = 0.24, mm = 0.0).Figure 5.16 shows the reconstruction for trial 5, based on real data. Parametersfor the trial are = 0.01 and eF = 0.1. With a stopping condition of — x11 1 < 0.1 thetrial requires 30 iterations. The reconstruction is qualitatively superior to those for RJCG.Structural detail is clearly visible and the image is not smoothed. However, there aresome artifacts visible in the interior of the reconstruction. These artifacts are aninevitable consequence of measurement errors. Although practical application requires(a)114better reconstructions, the POCS method shows that if the errors are reduced, a sharpaccurate reconstruction is possible.Figure 5.16 Reconstructed cross section for POCS real datatrial 5.Two further trials with POCS and real data show the effect of varying 8R’ theextent to which POCS forces the solution to fit the raysum data. Large ER relaxes the fitto data (appropriate for unreliable raysums) while a small ER enforces a tight fit to theraysum data. Figure 5.17 shows the reconstructed cross sections from the real data withER = 0.1 and ER = 0.001. Table 5.5 summarizes the results.(a)(b)Figure 5.17: Reconstructed images from POCS real data trial6: (a) ER = 0.1, and (b) ER = 0.001.115Trial Description Parameters IterationsNumber6 (a) raysums, 8R =0.1 7spatial support, and EF=O.lface sheets stop at11x — xi_i II < 0.16 (b) raysums, ER0.001 32spatial support, and EFO.lface sheets stop at— x1_ < 0.1Table 5.5: POCS trial 6 (real data) results.A large ER enlarges the constraint sets to ensure a non-empty intersection, andleads to faster convergence. Table 5.5 indicates that, as expected, POCS converges fastertowards a solution with ER large. Figure 5.17 show that the reconstruction is of a higherquality with the small ER.Plots of convergence for trial 6 in Figure 5.18 are instructive in understanding theeffects of ER. There is one plot each for trials 6(a) and 6(b) with two lines in each plot.One line in each plot has a value for each integer on the x axis, 1, 2, 3, ..., withcorresponding y-axis values showing:1. change in x from beginning of iteration 1 to after projection onto raysumconstraints,2. change in x from after projection onto raysum constraints to afterprojection onto fusion constraints in iteration 1,3. change in x from after projection onto fusion constraints to after projectiononto amplitude constraints in iteration 1,4. same as (1) but for iteration 2,5. same as (2) but for iteration 2,6. same as (3) but for iteration 2,etc.cDI.0Cl)Cci)CC-)C0D0Cl)U)D)CCt5-cC)iteration116The other line in each plot has data points for x-axis values 3, 6, 9, 12, ..., withcorresponding y-axis values indicating the change in x between full iterations. So oneline shows the movement in x within iterations (sub-iterations) while the other showschanges in x for full iterations.101100l010210110010-15 10 15 20iteration(a)sub iterationfull iteration —a-———I -20 40 60 800(b)Figure 5.18: Plots of convergence for POCS: (a) trial 6(a) -SR = 0.1 and (b) trial 6(b)- R = 0.001.The plot for trial 6(a) shows oscillations in sub-iteration movement that quicklysubside along with movement between full iterations. Trial 6(b), on the other hand,117shows sub-iteration movement that continues while the solution slows between iterations.Plot 6(a) suggest a non-empty intersection of constraints while 6(b) appears to have anempty intersection. In 6(b) the solution jumps around in a regular pattern within aniteration while the solution does not move much between iterations. These plots do notprovide conclusive evidence though that the intersections are non-empty and emptyrespectively.In spite of the apparent empty constraint intersection in trial 6(b), thereconstruction is good. Other than wasted iteration time, there is no penalty for making6R too small; the solution just bounces around in a region of good solutions. Iterationsare fast so it seems wise to select £R to be small. Only when the solution does notconverge does it make sense to set £R larger.5.3 Chapter SummaryThe salient conclusion of this chapter is that data fusion constraints markedlyimprove the accuracy of reconstruction produced by the limited-angle CT system. Allthree reconstruction methods exhibit improved accuracy by including spatial supportconstraints and further improvement by including face sheet constraints. In particular, theSVD trials show a decrease in the number of singularities due to the added constraints,and a corresponding reduction in the size of the null space. The fusion constraints do notarbitrarily reduce the null space, but they reduce the null space in a manner consistentwith measurements of the specimen. This consistency in the reduction of the null spaceresults in the improved accuracy of reconstruction.Although the source of real data used here is far from ideal, real data trials serveto show that the method is not restricted to synthetic data, but works with real data. Inparticular, POCS gives good reconstructions with the structure and anomalies clearlyvisible.RICG real data trials illustrate the compromises associated with the regularizationparameter a2. A large value of a2 leads to a smooth reconstruction. A small a2 does118not smooth as much, but conditioning becomes a problem. Selecting a2 proves to be acompromise between smoothing errors and conditioning errors.POCS gives better reconstructions than RJCG. POCS reconstructions are sharpand free of the oscillatory patterns exhibited by RICO. There is, however, someanomalous clutter in the interior of the image due to errors in the data. POCS is alsorelatively insensitive to its parameters. One can make ER too small without degrading thereconstruction and suffering only an increase in the number of iterations required.119Chapter 6Conclusions and DiscussionSection 6.1 of this chapter summarizes the important conclusions of this thesis,i.e.: identification of the inability of limited-angle CT to correctly reconstruct sandwichstructures, and the efficacy of the novel system using data fusion with ultrasound toovercome the inability. The section then generalizes the conclusions to cover a broaderfield of application. Section 6.2 extrapolates from the conclusions by discussingpossibilities for further development with respect to practical application of the limited-angle CT system, use of CAD models for constraints, other types of constraints, andapplication of fusion to Compton back scatter imaging.6.1 Conclusions6.1.1 Novel ContributionsCurrent limited-angle CT techniques rely on general a priori assumptions toconstrain reconstruction. This thesis contributes the novel observation that the success orfailure of reconstruction depends on the validity of these assumptions for the specificspecimen. It is not sufficient to claim that some limited-angle CT method is successfulwithout defining the types of specimens that yield accurate reconstruction. In particular,120current methods cannot accurately reconstruct a sandwich structure, as shown in Section3.1. The methods fail because the face sheets of sandwich structures lie almost entirely inthe limited-angle Radon transform null space. Consequently, there is no valid basis forinterpolation of the data to fill in missing data. This failure of limited-angle CT issignificant because sandwich structures are commonplace; they carry loads due tobending moments economically with respect to weight and material.Ultrasound measurements complement x-ray data well. Whereas x-rays areinsensitive to discontinuities along the direction of the raysum, ultrasound specificallylocates discontinuities along its path. This thesis describes a novel method for limited-angle CT that exploits the complementary nature of x-ray and ultrasound by fusingultrasound measurements into the reconstruction process. The method yields accuratereconstructions of sandwich structures where conventional limited-angle CT cannot.Ultimately, the fusion system requires formation and solution of the system oflinear equations:ER1 FI Ix=I[WJ [WxFRaysum data give the equations Rx=y and range and ultrasound data give Wx = WXF,where XF is a partial reference image and W is a diagonal matrix. The partial referenceimage contains the linear attenuation coefficients for the face sheet and exterior regions.Weights on the diagonal of W indicate whether or not a pixel in XF is known (i.e., in aface sheet or exterior) or unknown (i.e., in the interior). The range and ultrasoundequations remove the face sheet structures from the null space so that face sheets do notconfound reconstruction, thus allowing accurate reconstruction where it would otherwisebe impossible.Previous limited-angle CT work has used a full reference image [28], but the useof a partial reference image as a constraint set for POCS is novel. Therefore, the mappingof the fusion problem here to the constraint set of Equation (4-41) is also novel.1216.1.2 Generalization of ResultsThis thesis focuses on a specific problem and a specific type of specimen.Despite the narrow focus, the conclusions have a broader implication. The novel CTmethod succeeds because the data fusion focuses directly on the limitations of theoriginal problem. One cannot expect that the ad hoc addition of data to an inverseproblem will improve the solution. Additional constraints that do not focus onuncertainty in the original problem are not beneficial.For example, suppose that rather than use ultrasound data, the proposed methodused neutron radiation. Although neutrons have certain properties that complement xrays, they are not more sensitive to discontinuities. Since it is the failure of the x-rays tolocate discontinuities along the rays that confounds reconstruction, one should expect noimprovement due to fusion with neutron data (other than the cancellation of unbiasederrors by redundant measurements). Another example is CT, which is itself a fusionproblem that fuses a set of two-dimensional radiographs to produce a three-dimensionalreconstruction. However, as Chapter 2 points out, just fusing views from a few angles isnot sufficient. A full range of angles is necessary for reconstruction. Data acquired forone angle complement well data acquired at a perpendicular angle. This observation hasimplications for computational vision methods that derive scene descriptions frommultiple images, the opaque version of CT • It is important to select camera views tocomplement each other, i.e., an image should contain information that the others do nothave.Therefore, a key to success with data fusion is to ensure that additional datasources address uncertainty in the original problem. Ultrasound addresses well theweaknesses of x-rays, so the fusion method succeeds.* The opaque problem differs from CT because opacity provides additional constraints, but leads to acorrespondence problem instead.1226.2 DiscussionThis section presents four areas for further development of work in this thesis.First, it considers topics vital to practical application of the proposed limited-angle CTmethod, including improved apparatus and ability to detect anomalies. Second, Section6.2.2 discusses the use of CAD models (available for many NDE specimens) asconstraints. Third is an examination of other possibilities for constraints, specificallynon-convex ones. Fourth, and last, is a discussion speculating about fusion of Comptonback scatter data and ultrasound data.6.2.1 Practical ApplicationThe apparatus used in experiments is a patchwork of devices designed for digitalradiograph acquisition and is not well suited to the more rigorous requirements of CT.Although it is sufficient to perform the experiments presented in this thesis, practicalapplication of the proposed limited-angle CT method requires improved apparatus.Improvements to the apparatus must address sources of error in the data. CTtolerates small amounts of unbiased noise by canceling it out. Therefore, the smallamount of photon noise in the data is not an issue in the experimentation. Should photonnoise be a problem with a modified apparatus, longer integration times can compensate.Scattered radiation is a source of errors related to the specimen and is difficult to model.The best approach to elimination of scatter is collimation of the x-ray beam. Collimationprevents unnecessary radiation from reaching the specimen and reduces the subsequentscattered radiation. Whereas medical systems rely on collimation to reduce patientexposure to radiation, conventional radiography does not normally use collimation. Theapparatus used for experiments in this thesis does not have the collimation that isessential for practical application. Chapter 5, Section 5.1.3, shows the sensitivity ofraysum measurements to radiation intensity. The quantization interval in eight-bitintensity data corresponds to a large change in the raysum values. Twelve-bit data canimprove this sensitivity, provided that the additional bits are significant. A saner123acquisition process that avoids sending intensity signals through multiple amplifiers,digital-to-analogue and analogue-to-digital converters will help to reduce errors inintensity data and increase the number of significant bits. The linear array scanner isexcellent at its intended job, acquisition of digital radiographs, but it is not the best devicefor CT. For practical application of the proposed method, a more appropriate x-raysensor is essential.This thesis uses the 12 norm of the error, e, as a measure of accuracy where:e= lix — X X 100%.lixilx is the true cross section while * is the reconstruction. is an 12 norm. Because the 12norm averages the square of the errors over the entire solution, it gives a global measureof accuracy. For example, if * is identical to x everywhere except for a large error at onesample, the 12 norm will average the error at that sample over the whole image.Consequently, the 12 norm indicates a good global match between x and * while ignoringthe large local error. A greater number of samples averages out large local errors to agreater degree.A global error measure is well suited to some tasks, e.g., pattern recognition.However, for inspection tasks, where small local anomalies are important, an errormeasure sensitive to local errors is better. An norm of the error has such localproperties. For the previous example, the large local error determines the maximum erroryielding a large l, norm. The l norm of the error is important in identification ofanomalies because it determines the size of local deviation that is significant. Forexample, if one expects a pixel in a reconstruction to be 0.5 cm1 but observes instead0.51 cm’, the difference can only be interpreted to be significant if the reconstruction hasbetter than 0.01 cm1 accuracy in the 1,,, sense. In the 12 sense there is no way to know if alocal deviation is significant.Section 4.3.3 of Chapter 4 points out that the POCS constraints lead to a pseudol norm solution. However, minimizing the residual, as POCS does, is not necessarily124the same as minimizing the error. The question of how to minimize the local error meritsconsideration.The proposed limited-angle CT method has potential for application beyondinspection of sandwich structures. Often an outer shell obscures an internal structure justas face sheets obscure the interior of a sandwich structure. If such an obscured structurealso requires limited-angle techniques, fusion with ultrasound will be helpful. Theultrasound allows one to mathematically peal away the outer shell so that it does notconfound reconstruction of the interior. Thus, the method is not limited to controlsurfaces, but also allows inspection of some internal structures in situ that wouldotherwise require disassembly.6.2.2 CAD Model as a Source of ConstraintsWith the advent of computer-aided design (CAD) it has become commonplace todesign structures with computers. Consequently, many structures targeted for inspectionhave CAD models available. The question then arises: can CAD models provide asubstitute for ultrasound data in the proposed limited-angle CT method? In short, theanswer is yes, depending on two factors.First, it is essential that face sheet data from the CAD model register preciselywith the raysum data. This is not a trivial concern, but such registration is within therealm of current technology.Second, one must consider the effects of errors in the CAD model, i.e.,discrepancies between the CAD model and the true structure. Such errors may arise frommanufacture or from changes in the structure during service. The success of the methodthen depends upon the ability of the range and raysum data to correct the model. The nullspace of the limited-angle Radon transform provides a guide to determining what errorsare correctable. Raysum data cannot correct errors that are wide and thin, like the facesheets, while they can correct narrow ones. Therefore, small local variations in the face125sheet data should not be a problem. A global error in face sheet thickness or a wideanomaly in the face sheet will confound reconstruction.POCS provides an advantage here. The mixture of 12 and fits allowed byPOCS means that raysum data can make local corrections to the CAD data. While the lfit forces a degree of local conformity to the raysum data, some local errors in the CADmodel are tolerable, provided that they are not so large that they affect the global fit.Substitution of ultrasound data with a CAD model offers important benefits.Elimination of ultrasound allows simultaneous acquisition of all the data, and thuscontributes a major time savings. Also, some structures may be difficult to test withultrasound. In such cases the CAD model makes it easier to carry out a limited-angle CTinspection.6.2.3 More Sophisticated ConstraintsThe proposed limited-angle CT method relies primarily on the data for solutionbut also uses some a priori assumptions. Solution by RJCG assumes a smooth solution(with dubious success) and POCS assumes fixed amplitude limits for the solution. Thenumerical methods examined require either a convex objective function (for CG) orconvex constraint sets (for POCS). However, there are other legitimate assumptions thatmerit consideration and do not fit directly into CG or POCS. Two assumptions aboutgraphite/epoxy and aluminum honeycomb sandwiches not considered in the previouschapters are:1. a finite set of materials in the specimen, and2. vertical continuity in the honeycomb.Assumption one is reasonable because one usually knows exactly what materials arepossible in the specimen. In the case of the plexiglass phantom, only plexiglass and airare possible while in a graphite/epoxy and aluminum honeycomb sandwich, onlygraphite/epoxy, aluminum, and air are possible. One can also allow for the existence ofmaterials in typical anomalies such as water or oxidization products. A finite set of126possible materials leads to non-convex constraints in both CG and POCS. Assumptiontwo is reasonable because, with the exception of anomalies, a material in the interior hasidentical material above and below it. Vertical continuity leads to convex constraints andis not special, but it has the potential to improve reconstruction for thin honeycomb byreflecting known honeycomb properties in the constraints. Here we focus on the moredifficult problem of implementation of the non-convex constraints.In the context of RJCG, constraining the solution to consist of a finite set ofmaterials requires that the objective function contain local minima at positionscorresponding to allowed materials. With many local minima, the objective function isnon-convex and CG is not guaranteed to find the global minimum, or even a closeapproximation. One method for minimizing such non-convex objective functions is thegraduated non-convexity algorithm (GNC) [6] described as follows:(1) create a new objective function with a parameter, say a, such that asa — oo the function is convex, and as a —> o the function becomes theoriginal non-convex function,(2) set a large and solve the effectively convex minimization,(3) reduce a and solve the locally convex minimization at the previousminimum,(4) repeat step (3) until a is sufficiently small.GNC performs a non-convex minimization by successively minimizing the locallyconvex objective function for decreasing values of a.Realizing that GNC only ever needs to minimize a locally convex objectivefunction suggests a simplification. Restrict the problem to minimizing the followingquadratic objective function:127R 2-WxE= a x—a FaD1 0aD2 0where all variables and parameters are as described in Chapter 4. The followingprocedure performs a minimization akin to GNC but only solves linear systems:(1) minimize E with a large, i.e., biased heavily towards a smooth solutionand ignoring XF,(2) if a pixel in the solution is near one of the allowed material values, assumethat the pixel is that material and modify W and XF accordingly,(3) reduce a and repeat the minimization,(4) repeat steps (2) and (3) until a is sufficiently small.The procedure starts by finding a very smooth solution. As the smoothness constraintrelaxes, the solution migrates towards one that contains only allowed values.Consideration of only one material value for each pixel avoids non-convexity, leads to afast linear algorithm, and arrives at essentially the same solution as GNC. Figure 6.1shows a reconstruction of the plexiglass phantom using the above method assuming thatthe solution contains only air (i’ = Ocm’) and plexiglass (u = O.45cm’). Thereconstruction is the result of the eleven quadratic optimizations summarized in Table6.1. No attempt is made to find an optimal schedule for reducing a.Figure 6.1: Reconstruction from successive RJCG iterationsincorporating the assumption that pixels are either air(4u = Ocm’) or plexiglass (u = O.45cm’).128Optimization a2 - 1 Number ofNumber Iterations1 8 1/8 572 4 1/4 423 2 1/2 274 1 1 125 1/2 2 106 1/4 4 137 1/8 8 298 1/16 16 659 1/32 32 11110 1/64 64 7411 1/128 128 5Table 6.1: Summary of successive minimizations used toproduce reconstruction in Figure 6.1.POCS can also solve the non-convex problem. Define the constraint set CM. as:CM = {x : lxi— Xml 6M’ Xmwhere 9( is the set of allowed values in the solution. CM is a potentially non-convex setconsisting of regions about solutions containing allowed values, with the parameter EMdetermining the size of the regions. If EM is large the regions overlap and CM is convex,but if EM is small the regions do not overlap and CM is non-convex. The correspondingprojection operator is:z:z =Xm+EM Z >Xm+EM XmflnfliflMZeSlzjXl9PMz= z:z, =XmEM Z <XmEM xmm1ninuzelzjl(.z otherwisePOCS finds a solution in a manner analogous to GNC by varying EM as follows:(1) start with EM large and perform the reconstruction starting from x =0,(2) reduce EM and repeat the reconstruction starting from the previoussolution,(3) repeat step (2) until EM is sufficiently small.The method works because the constraint set CM is always locally convex, i.e., once theprojection operator has determined which material a pixel matches best, it considers onlythat material. As EM diminishes it may be necessary to increase the parameter ER toensure that the intersection of constraints is non-empty. Figure 6.2 shows a129reconstruction based on the above method using 8R = 0.1, CF = 0.1, andEM = 1, , i,.. . , stopping when 1x — x11 11<0.1. Note that with the assumption of afinite set of materials the amplitude constraint is redundant and therefore it is omitted.Table 6.2 summarizes the iterations in the reconstruction.POCS Solution EM Number ofNumber Iterations1 1 12 1/2 13 1/4 14 1/8 45 1/16 96 1/32 97 1/64 48 1/128 29 1/256 210 1/512 2Table 6.2: Summary of successive POCS applications used toproduce reconstruction in Figure 6.2.In summary, it is possible to incorporate non-convex constraints in thereconstruction process. Figures 6.1 and 6.2 show the results of using such constraintswith RJCG and POCS respectively. It is not clear from this evidence that assuming afinite set of possible materials improves the accuracy of reconstruction. More accuratemeasurements using superior apparatus may be necessary to gain a real benefit from suchconstraints.Figure 6.2: Reconstruction from successive POCS iterationsincorporating the assumption that pixels are either air(u = Ocm’) or plexiglass (u = 0.45cm’).1306.2.4 Compton Scatter and Ultrasound FusionThis section speculates about the potential for application of data fusion toCompton back scatter imaging and ultrasound. Figure 6.3 shows schematically aCompton back scatter imaging system. An x-ray source projects a collimated beam intothe specimen. As the specimen absorbs the x-rays, it emits scattered radiation along thebeam. Compton effect scattering radiates in all directions, including back towards thesource. An x-ray pin-hole camera directs radiation from varying depths in the specimento an array of sensors producing intensity signals related to the amount of scatteroriginating at each depth. Philips COMSCAN, a commercial version of the apparatus,scans the specimen over two dimensions to produce a volumetric image. The methodpotentially offers high contrast images of anomalies not visible to conventionalradiography*, and requires access to only one side of the specimen. See Chapter 2 for adescription of ultrasound.At present, there are no algorithms for solving the inverse problem posed by backscatter imaging. Instead, commercial devices, such as COMSCAN, simply scale the rawintensity data to produce an image. The inverse problem is difficult for two reasons:1. linear attenuation varies greatly depending on position in the path fromsource to sensor because:a. the spectral content of the incident radiation is different from thescattered radiation, andb. beam hardening changes the spectral content of both incident andscattered radiation,and2. radiation intensity diminishes with depth into the specimen, so noise andsensitivity increase.* Theoretically, the contrast for a void in a homogenous specimen is infinite, but in practice spatialresolution limits contrast.131A possible approach to solving the inverse problem is to reconstruct the specimen onelayer at a time, starting from the top and working down. The top layer is easilyreconstructed because there is only air above it. Subsequent layers depend only onattenuation values for layers above so a downward step-wise procedure can work.Compton scatter emittedalong the path of radiationFigure 5.3: Apparatus for Compton back scatter imaging.For ideal lossless media, ultrasound inversion is tractable, but in reality, materialsare not lossless and ultrasound inversion is difficult. Losses confound ultrasoundinversion in two ways:1. losses are inconsistent and unknown, and2. losses cause signal strength to diminish with depth and reduce accuracy.For example, the effects of losses can make it difficult to distinguish between the normalreflection from an interface deep in the specimen and the reflection off of an air gap inthe interface (an anomaly).Both back scatter and ultrasound inversion can operate independently. However,they are both prone to errors with increasing depth in the specimen. The two methods do132not complement each other as in the proposed CT method, but their fusion may still bebeneficial. For example, a strongly-coupled recurrent fusion system can combine thelayer-by-layer back scatter reconstruction with a similar ultrasound algorithm. As thereconstruction proceeds downwards, the algorithm must decide what material (or whatinterface in the ultrasound model) to put in the reconstruction. When both back scatterand ultrasound data concur the algorithm proceeds in a straight-forward manner. Whenone has ambiguous data the other can suggest the correct interpretation. For example,when ultrasound cannot determine if an interface has an air gap in it, it can assume thereis no air gap if the back scatter data do not indicate a drop in scattered radiation.Likewise, a drop in scattered radiation deep in the specimen may be due to a small voidor just a random variation in the data. If it is a void, however, there should be acorresponding reflection in the ultrasound data.It is not certain that such a fusion system will lead to a superior inspectiontechnique. Nevertheless, if the system works, it will constitute an important contributionto NDE.133Bibliography[1] Artzy, E, Elfving, T, and Herman, G.T., “Quadratic Optimization for ImageReconstruction, II,” Computer Graphics and Image Processing, vol. 11, pp 242-261, 1979.[2] Axeisson, 0. and Barker, V.A., Finite Element Solution of Boundary ValueProblems Theory and Computation, Academic Press, Orlando, 1984.[3] Bates, D.M. and Wahba, G., “Computational Methods for Generalized Cross-Validation with Large Data Sets,” Treatment of Integral Equations byNumerical Methods, C.T.H. Baker and G.F. Miller editors, Academic Press,London, 1982.[4] Bauschke, H.H. and Borwein, J.M., “On Projection Algorithms for SolvingConvex Feasibility Problems,” Simon Fraser University, Department ofMathematics and Statics Research Report No. 93-12, June, 1993.[5] Bertero, M., Poggio, T.A. and Torre, V., “Ill-Posed Problems in Early Vision,”Proceedings of the IEEE, Vol. 76, No. 8, pp 869-889, August 1988.[6] Blake, A., and Zisserman, A., Visual Reconstruction, MIT Press, Cambridge,MA, 1987.[7] Brègman, L.M., “The Method of Successive Projection for Finding a CommonPoint of Convex Sets,” translated by J.G. Ceder, Dokiady, Vol. 162, No. 3, pp688-692, 1965.[8] Censor, Y. and Herman, G.T., “Row-Generation Methods for Feasibility andOptimization Problems Involving Sparse Matrices and Their Applications,”Sparse Matrix Proceedings 1978, SIAM, Philadelphia, pp 197-2 19, 1978.[9] Clark, J.J. and Yuille, A.L., Data Fusion for Sensory Information ProcessingSystems, Kiuwer Academic Publishers, Boston, 1990.134[10] Davison, M.E., “The Ill-Conditioned Nature of the Limited Angle TomographyProblem,” SIAM Journal ofApplied Mathematics, Vol. 43, No. 2, pp 428-448,1983.[11] Dusaussoy, N.J. and Abdou, I.E., “The Extended MENT Algorithm: A MaximumEntropy Type Algorithm Using Prior Knowledge for ComputerizedTomography,” IEEE Transactions on Signal Processing, Vol. 39, No. 5, pp 1164-1180, 1991.[12] Forsythe, G.E., Malcolm, M.A. and Moler, C.B., Computer Methods forMathematical Computations, Prentice-Hall, Englewood Cliffs, N.J., 1977.[13] Golub, G.H. and Van Loan, C.F., Matrix Computations, The Johns HopkinsUniversity Press, Baltimore, Maryland, 1983.[14] Golub, G.H., Heath, M., and Wahba, G., “Generalized Cross Validation as aMethod of Choosing a Good Ridge Parameter,” Technometrics, Vol. 21, pp 215-213, 1979.[15] Goupillaud, P.L., “An Approach to Inverse Filtering of Near-Surface LayerEffects from Seismic Records,” Geophysics, Vol. 26, No. 6, pp 754-760, 1961.[16] GrUnbaum, F.A., “A Study of Fourier Space Methods for Limited Angle ImageReconstruction,” Numer. Funct. Anal. and Optimiz., Vol. 2, No. 1, pp 3 1-42,1980.[17] Gubin, L.G., Polyak, B.T., and Raik, E.V., “The Method of Projections forFinding the Common Point of Convex Sets.” U.S.S.R Computational Mathematicsand Mathematical Physics, Vol. 7, No. 6, pp 1-24, 1967.[18] Habibi-Ashrafi, F. and Mendel, J.M., “Estimation of Parameters in LosslessLayered Media Systems,” IEEE Transactions on Automatic Control, Vol. AC-27,No. 1, pp 31-48, 1982.[19] Halmshaw, R., Industrial Radiology: Theory and Practice, Applied SciencePublishers, London, 1982.[20] Heiskanen, K.A., Rhim, H.C. and Monteiro, J.M., “Computer Simulations ofLimited Angle Tomography of Reinforced Concrete,” Cement and ConcreteResearch, Vol. 21, pp 625-634, 1991.[21] Herman, G.T. and Lent, A., “Quadratic Optimization for Image Reconstruction I,”Computer Graphics and Image Processing, Vol 5, pp 3 19-332, 1976.[22] Herman, G.T., Image Reconstruction from Projections: The Fundamentals ofComputerized Tomography, Academic Press, New York, 1980.[23] Louis, A.K., “Incomplete Data Problems in X-Ray Computerized Tomography,”Numerische Mathematik, Vol. 48, pp 25 1-262, 1986.[24] McRae, K.I., “Deconvolution Techniques for Ultrasonic Imaging of AdhesiveJoints”, Materials Evaluation, pp 1380-1384, 1990.135[25] Mendel, J.M. and Habibi-Ashrafi, F., “A Survey of Approaches to SolvingInverse Problems for Lossless Layered Media Systems”, IEEE Transaction onGeoscience and Remote Sensing, Vol. GE-18, No. 4, pp 320-330, 1980.[261 Mendel, J.M. and Goutsias, J., “One-Dimensional Normal-Incidence Inversion: ASolution Procedure for Band-Limited and Noisy Data,” Proceedings of the IEEE,Vol. 74, No. 3, pp 401-4 14, 1986.[27] Oskoui-Fard, P, and Stark, H., “Tomographic Image Reconstruction Using theTheory of Convex Projections,” IEEE Transactions on Medical Imaging, Vol. 7,No. 1, March 1988.[28] Oskoui-Fard, P, and Stark, H., “A Comparative Study of Three ReconstructionMethods for a Limited-View Computer Tomography Problem,” IEEETransactions on Medical Imaging, Vol. 8, No. 1, March 1989.[29] Pierre, D.A., Optimization Theory with Applications, Dover Publications, NewYork, 1986.[30] Poggio, T., Torre, V. and Koch, C., “Computational Vision and RegularizationTheory,” Nature, Vol. 26, pp 3 14-319, September 1985.[31] Press, W.H., Plannery, B.P., Teukolsky, S.A. and Vetterling, W.T., NumericalRecipes in C the Art of Scientific Computing, Cambridge University Press,Cambridge, 1990.[32] Radon, J., “Uber die Bestimmung von Funktionen durch ihre Integralwerte langsgewisser Mamiigfaltigkeiten,” Ber. Verb. Saechs. Akak. Wiss., Leipzig, Math.Phys. Kl., Vol. 69, pp 262-277, 1917.[33] Sato, T., Norton, S.J., Linzer, M., Ikeda, 0. and Hirama, M., “Tomographic ImageReconstruction from Limited Projections using Iterative Revisions in Image andTransform Spaces,” Applied Optics, Vol. 20, No. 3, pp 395-399, 1981.[34] Scudder, H. J., “Introduction to Computer Aided Tomography”, Proceedings ofthe IEEE, Vol. 66, No. 6, pp 628-637, 1978.[35] Servo-Robot, Saturn-2000 Integrated 3-D Laser Vision System ProductDescription, Montreal, Canada, circa 1989.[36] Soumekh, M., “Image Reconstruction from Limited Projections using the AngularPeriodicity of the Fourier Transform of the Radon Transforms of an Object,”AcousticalImaging, Vol. 13, Proceedings of the Thirteenth InternationalSymposium, Minneapolis, MN, USA, pp 31-42, October 1983.[37] Stoer, J.and Bulirsch, R., Introduction to Numerical Analysis, Springer-Verlag,New York, 1980.[38] Strang, G., Linear Algegra and Its Applications, 2nd ed., Accademic Press,1980.[39] Strang, G., Introduction to Applied Mathematics, Wellesley-Cambridge Press,Wellesley, MA, 1986.136[40] Tam, K.C., Eberhard, J.W. and Mitchell, K.W., “Incomplete-Data CT ImageReconstructions in Industrial Applications,” IEEE Transactions on NuclearScience, Vol. 37, No. 3, pp 1490-1499, 1990.[41] Tikhonov, A.N. and Arsenin, V.Y., Solutions of Ill-Posed Problems, translationeditor F. John, V.H. Winston and Sons, Wahsington, D.C., 1977.[42] Trummer, M.R., “A Note on the ART of Relaxation,” Computing, Vol. 33, pp349-352, 1984.[43] Tuy, H., “Reconstruction of a Three-Dimensional Object from a Limited Range ofViews,” Journal ofMathematical Analysis and Applications, Vol. 80, pp 598-616,1981.[44] Youla, D.C. and Webb, H., “Image Restoration by the Method of ConvexProjections: Part 1 - Theory,” IEEE Transactions on Medical Imaging, Vol. MI1, No. 2, October 1982.[45] Zala, C.A. and McRae, K.I., “An Optimization Method for Acoustic ImpedanceEstimation of Layered Structures using Prior Knowledge”, Acoustical Imaging,Vol. 18, pp 363-372, 1991.[46] Zala, C.A. and Churchill, L., “Inversion of Ultrasonic Traces,” DREPContractor’s Report Series 90-18, Victoria, B.C, Canada, 1990.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A novel constraint-based data fusion system for limited-angle...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A novel constraint-based data fusion system for limited-angle computed tomography Boyd, Jeffrey Edwin 1994
pdf
Page Metadata
Item Metadata
Title | A novel constraint-based data fusion system for limited-angle computed tomography |
Creator |
Boyd, Jeffrey Edwin |
Date Issued | 1994 |
Description | Computed tomography (CT) is a non-destructive evaluation technique that reconstructs the cross section of a specimen from x-ray raysum measurements. Whereas CT reconstruction is an ill-posed inverse problem that is easily solved, limited-angle CT, where raysum data are missing for a range of angles, is more severely ill-posed and more difficult to solve. In the limited-angle case, a priori assumptions are necessary to constrain the problem. Specimens wider than the x-ray source to sensor spacing require limited-angle CT. Furthermore, if the specimen is a sandwich structure, i.e., some core material surrounded by load-bearing face sheets, then the face sheets must lie in the null space. Components in the null space do not appear in the raysum data and thus confound CT reconstruction because there is no basis for interpolation. This thesis proposes a novel constraint-based data fusion method for limited-angle CT reconstruction of sandwich structures. The method reduces the reliance of limited-angle CT on assumptions by using range and ultrasound measurements to constrain the solution. Fusion of the data sources results in a problem with a much smaller null space that no longer includes the face sheets. The reduction of the null space in a manner consistent with the specimen yields a more accurate tomographic reconstruction. Synthetic and real data experiments show marked improvement in reconstruction accuracy achieved by using the fusion system. |
Extent | 3525442 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-04-07 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0051601 |
URI | http://hdl.handle.net/2429/6877 |
Degree |
Doctor of Philosophy - PhD |
Program |
Computer Science |
Affiliation |
Science, Faculty of Computer Science, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1994-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1994-893281.pdf [ 3.36MB ]
- Metadata
- JSON: 831-1.0051601.json
- JSON-LD: 831-1.0051601-ld.json
- RDF/XML (Pretty): 831-1.0051601-rdf.xml
- RDF/JSON: 831-1.0051601-rdf.json
- Turtle: 831-1.0051601-turtle.txt
- N-Triples: 831-1.0051601-rdf-ntriples.txt
- Original Record: 831-1.0051601-source.json
- Full Text
- 831-1.0051601-fulltext.txt
- Citation
- 831-1.0051601.ris
Full Text
Cite
Citation Scheme:
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>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0051601/manifest