@prefix vivo: . @prefix edm: . @prefix ns0: . @prefix dcterms: . @prefix skos: . vivo:departmentOrSchool "Applied Science, Faculty of"@en, "Mechanical Engineering, Department of"@en ; edm:dataProvider "DSpace"@en ; ns0:degreeCampus "UBCV"@en ; dcterms:creator "Angus, Edward J."@en ; dcterms:issued "2015-11-05T03:35:20"@en, "2015"@en ; vivo:relatedDegree "Master of Applied Science - MASc"@en ; ns0:degreeGrantor "University of British Columbia"@en ; dcterms:description """The greatest cost in sawmills is contained in uncut logs. Significant increases in profit stand to be made if logs are properly processed so that the amount of defect and feature-free products is maximized. This has driven research into internal sensing such as computed tomography (CT) for use in sawmills. Traditional CT systems from the medical industry are ill-suited for industrial use and provide distorted reconstructions when used with less than perfect or partially incomplete data. Industrial CT systems in areas such quality control have embraced so-called Algebraic Reconstruction Techniques (ART) for robustness to errant data, but these systems are difficult to provide reconstructions rapidly enough for sawmill use. A system is proposed here for an ART scanning arrangement specific for log-scanning use. The image space's voxel pattern is defined to reflect the geometry of a log's internal features. This greatly reduces the number of unknowns without a loss of information and allows for faster reconstruction. The geometry of the voxels makes exact calculation and storage of the series solution basis fast and practical for multi-slice scans. Data scaling and normalization eliminate unnecessary voxels in the reconstruction. It also makes the reconstruction scheme tolerant of rigid body motion radially and circumferentially about the cone beam source. The voxel numbering scheme means that knot features are contained in voxels that are numerically close for quick registration. A camera-based detector system was implemented to collect radiographs of logs. Logs were translated and rotated in an x-ray cone beam and their position was monitored by an optical encoder. The detector was activated at the proper intervals to yield radiographic data. The series solution basis for this helical movement was constructed. An iterative solver and selective low-pass filtering was found to provide good reconstruction results. Segmentation was implemented to demonstrate use of the reconstructions for lumber processing. Eccentric spiral motion tests demonstrated the effectiveness of the data scaling and normalization process. Results are provided that point towards efficient automated registration of features."""@en ; edm:aggregatedCHO "https://circle.library.ubc.ca/rest/handle/2429/55136?expand=metadata"@en ; skos:note "Robust CT Scanning of Logs with Feature-Tailored VoxelsbyEdward J. AngusB.S. Mechanical Engineering, Iowa State University, 2013a thesis submitted in partial fulfillmentof the requirements for the degree ofMaster of Applied Scienceinthe faculty of graduate and postdoctoral studies(Mechanical Engineerng)The University of British Columbia(Vancouver)October 2015© Edward J. Angus, 2015AbstractThe greatest cost in sawmills is contained in uncut logs. Significant increases in profit stand tobe made if logs are properly processed so that the amount of defect and feature-free productsis maximized. This has driven research into internal sensing such as computed tomography(CT) for use in sawmills. Traditional CT systems from the medical industry are ill-suited forindustrial use and provide distorted reconstructions when used with less than perfect or partiallyincomplete data. Industrial CT systems in areas such quality control have embraced so-calledAlgebraic Reconstruction Techniques (ART) for robustness to errant data, but these systemsare difficult to provide reconstructions rapidly enough for sawmill use.A system is proposed here for an ART scanning arrangement specific for log-scanning use.The image space's voxel pattern is defined to reflect the geometry of a log's internal features.This greatly reduces the number of unknowns without a loss of information and allows forfaster reconstruction. The geometry of the voxels makes exact calculation and storage of theseries solution basis fast and practical for multi-slice scans. Data scaling and normalizationeliminate unnecessary voxels in the reconstruction. It also makes the reconstruction schemetolerant of rigid body motion radially and circumferentially about the cone beam source. Thevoxel numbering scheme means that knot features are contained in voxels that are numericallyclose for quick registration.A camera-based detector system was implemented to collect radiographs of logs. Logs weretranslated and rotated in an x-ray cone beam and their position was monitored by an opticalencoder. The detector was activated at the proper intervals to yield radiographic data. Theseries solution basis for this helical movement was constructed. An iterative solver and selectivelow-pass filtering was found to provide good reconstruction results. Segmentation was imple-mented to demonstrate use of the reconstructions for lumber processing. Eccentric spiral motiontests demonstrated the effectiveness of the data scaling and normalization process. Results areprovided that point towards efficient automated registration of features.iiPrefaceAll the work presented here was carried out at the Renewable Resources Laboratory at theUniversity of British Columbia's Point Grey campus.I was the principal investigator and was responsible for all major areas of data collection andanalysis. Dr. Gary Schajer, my supervisor, was responsible for the concept of a feature-specifictimber algebraic CT algorithm. Dr. Anthony An was responsible for the seminal constructionof the experimental apparatus and control system, as well as the rotational algorithms. Iterativesolvers from AIR Tools were developed by Christian Hansen of Denmark Technical University.Production of the helical reconstruction algorithm and its synthesis with modified smoothingiterative solvers was my own.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiNomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .xvii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Overview of log grading and scaling . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Sensing technologies in sawmills . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.3 X-ray sensing in sawmills . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.4 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 Computed tomography overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.1 Basics of computed tomography . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.1.1 Computed tomography setup . . . . . . . . . . . . . . . . . . . . . . . . . 72.1.2 Beer's law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.2 Overview of CT terminology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.3 The formal CT reconstruction problem and solution . . . . . . . . . . . . . . . . 11iv2.3.1 The Radon transform and problem statement . . . . . . . . . . . . . . . . 112.3.2 Filtered backprojection overview . . . . . . . . . . . . . . . . . . . . . . . 152.3.3 Series expansion methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.3.4 Comparing methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.4 Three dimensional reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213 The feature-tailored voxel approach . . . . . . . . . . . . . . . . . . . . . . . . . . 223.1 Overview of the method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.1.1 Method description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.1.2 Feature-tailored voxel geometries . . . . . . . . . . . . . . . . . . . . . . . 223.1.3 Data binning, centering, and normalization . . . . . . . . . . . . . . . . . 243.2 Calculation of basis for a single slice . . . . . . . . . . . . . . . . . . . . . . . . . 313.2.1 Path length basis function geometry by voxel arrangement view . . . . . . 313.2.2 Path length permutations by view . . . . . . . . . . . . . . . . . . . . . . 343.3 Calculation of the basis functions for multiple slices . . . . . . . . . . . . . . . . . 373.3.1 Multi-slice basis function geometry . . . . . . . . . . . . . . . . . . . . . 373.3.2 Permutations by view for multi-slice rotational reconstruction . . . . . . . 403.3.3 Permutations by view for multi-slice helical reconstruction . . . . . . . . 433.4 Advantages of the method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 474 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 504.1 Detector assembly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 524.2 X-ray source specifications and shielding . . . . . . . . . . . . . . . . . . . . . . . 564.3 Log carriage and locomotion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 584.3.1 Log carriage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 584.3.2 Belt drive rotation and helical movement . . . . . . . . . . . . . . . . . . 594.3.3 Control scheme, hardware, and sensing . . . . . . . . . . . . . . . . . . . . 604.4 Test specimens . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 624.4.1 Physical phantom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 624.4.2 Tree specimens . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 645 Single slice reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 665.1 Direct solution via the normal equations . . . . . . . . . . . . . . . . . . . . . . . 665.2 Synthetic phantom reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . 67v5.2.1 Modified centered finding . . . . . . . . . . . . . . . . . . . . . . . . . . . 695.2.2 Tikhonov regularization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 715.3 Real data reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 725.3.1 Calibration of the detector . . . . . . . . . . . . . . . . . . . . . . . . . . . 725.3.2 Real log data sample and reconstruction . . . . . . . . . . . . . . . . . . . 746 Multiple slice reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 776.1 Computational considerations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 776.1.1 Storage of the basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 776.1.2 Ill-posedness of the basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 776.2 Iterative methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 786.2.1 Randomized Kaczmarz's method . . . . . . . . . . . . . . . . . . . . . . . 786.2.2 Simultaneous algebraic reconstruction technique . . . . . . . . . . . . . . 806.2.3 Selective regularization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 816.2.4 Algorithm comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 836.3 Rotational multiple slice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 846.3.1 Phantom reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 846.4 Helical scanning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 906.4.1 Effect of limited views . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 906.4.2 Phantom reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 906.4.3 Real log reconstruction of dry log . . . . . . . . . . . . . . . . . . . . . . . 956.4.4 Real log reconstruction of green log . . . . . . . . . . . . . . . . . . . . . . 987 Segmentation of log internal features . . . . . . . . . . . . . . . . . . . . . . . . . 1027.1 Reconstruction of lumber surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . 1027.2 Feature identification via data plots . . . . . . . . . . . . . . . . . . . . . . . . . . 1068 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1098.1 Limitations of the method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1098.1.1 Assumptions made . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1098.1.2 Difficulty of integration into mill . . . . . . . . . . . . . . . . . . . . . . . 1098.2 Future work and considerations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1108.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113viList of Tables4.1 Angular θ coordinates of phantom knots . . . . . . . . . . . . . . . . . . . . . . . 636.1 Slice B knot locations - measured and reconstructed . . . . . . . . . . . . . . . . 876.2 Centered spiral reconstructed phantom knots . . . . . . . . . . . . . . . . . . . . 92viiList of Figures1.1 Linear Radiograph . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.2 Multi-source detector system and output from Seger et al. . . . . . . . . . . . . . 41.3 CT detector configuration and output from Seger et al. . . . . . . . . . . . . . . 52.1 CT configuration showing detector signal . . . . . . . . . . . . . . . . . . . . . . . 82.2 Object rotated to discern bodies B, C . . . . . . . . . . . . . . . . . . . . . . . . 82.3 Beer-Lambert Law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.4 Typical CT measurement geometries . . . . . . . . . . . . . . . . . . . . . . . . . 102.5 Outline of CT configurations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.6 Distributed density function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.7 The Radon Transform coordinate space . . . . . . . . . . . . . . . . . . . . . . . 132.8 Graphical representation of Radon transform . . . . . . . . . . . . . . . . . . . . 142.9 Filtered Backprojection Steps Visualization . . . . . . . . . . . . . . . . . . . . . 162.10 Series expansion representation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183.1 Medical grade CT reconstruction of softwood specimen from Szathmary . . . . . 233.2 Voxel reconstruction cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 243.3 Log setup in 2D detector fan-beam case . . . . . . . . . . . . . . . . . . . . . . . 253.4 Cylindrical binning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.5 Uniform density in a bounded circle and its Radon Transform . . . . . . . . . . . 273.6 Relating elliptical Radon Transform width to height . . . . . . . . . . . . . . . . 283.7 Vertical log shift . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 293.8 Horizontal log shift . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303.9 Data normalization scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303.10 Annular Case (a) Pixel intersection (b) Gview for annular case . . . . . . . . . . . 323.11 Sector Case (a) Pixel intersection (b) Gview for sector case . . . . . . . . . . . . . 33viii3.12 Combined pixel ray intersection . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343.13 Annular single slice total basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 353.14 Sector rotation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 363.15 Sector single slice total basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 363.16 Combined single slice total basis . . . . . . . . . . . . . . . . . . . . . . . . . . . 373.17 3D CT arrangement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 383.18 Voxels intersected in multi-slice reconstruction . . . . . . . . . . . . . . . . . . . . 393.19 Multiple slice ray entrance and exit . . . . . . . . . . . . . . . . . . . . . . . . . . 403.20 Path length matrix for 3D view Gview_MS . . . . . . . . . . . . . . . . . . . . . . 413.21 Multi-slice rotational G construction . . . . . . . . . . . . . . . . . . . . . . . . . 423.22 Motion of the log in view registers . . . . . . . . . . . . . . . . . . . . . . . . . . 443.23 Representation of Gview_spiral . . . . . . . . . . . . . . . . . . . . . . . . . . . . 453.24 Helical basis matrix G construction . . . . . . . . . . . . . . . . . . . . . . . . . . 463.25 Cartesian segmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 483.26 Combined voxel numbering breakdown . . . . . . . . . . . . . . . . . . . . . . . . 494.1 Experimental Setup Layout . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 514.2 Log scanning setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 524.3 Andor iXon 897 layout from An . . . . . . . . . . . . . . . . . . . . . . . . . . . 534.4 Detector box . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 544.5 Detector box diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 544.6 Barrel correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 554.7 X-ray tube . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 564.8 Rough alignment via laser . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 574.9 Sufficient fine alignment via perspective image . . . . . . . . . . . . . . . . . . . . 584.10 Log carriage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 594.11 Timing belt conveyors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 604.12 Log scanning control scheme from An . . . . . . . . . . . . . . . . . . . . . . . . 614.13 Encoder on shaft . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 614.14 Phantom dimensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 634.15 Log phantom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 644.16 Douglas fir samples (left) dry log (right) green log . . . . . . . . . . . . . . . . . . 655.1 Synthetic phantoms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68ix5.2 Center artifact in reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . 695.3 Center finding example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 705.4 Tikhonov regularization matrix for curvature . . . . . . . . . . . . . . . . . . . . 715.5 Phantom Artifact Correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 725.6 Fiberboard calibration blocks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 735.7 Calibration data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 745.8 Single Slice data sample from An . . . . . . . . . . . . . . . . . . . . . . . . . . . 755.9 Single Slice reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 766.1 Kaczmarz algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 796.2 Overdetermined Kaczmarz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 806.3 Selective regularization voxel indexing . . . . . . . . . . . . . . . . . . . . . . . . 826.4 Rotational multi-slice data capture . . . . . . . . . . . . . . . . . . . . . . . . . . 856.5 Centered multi-slice reconstruction slices . . . . . . . . . . . . . . . . . . . . . . . 866.6 Centered multi-slice knot directions . . . . . . . . . . . . . . . . . . . . . . . . . . 876.7 Off-center counterweights . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 886.8 Rigid body motion multi-slice data . . . . . . . . . . . . . . . . . . . . . . . . . . 886.9 Off-center multi-slice reconstruction slices . . . . . . . . . . . . . . . . . . . . . . 896.10 Centered spiral reconstruction data . . . . . . . . . . . . . . . . . . . . . . . . . . 916.11 Centered spiral reconstruction of physical phantom . . . . . . . . . . . . . . . . . 926.12 Off-center spiral reconstruction of physical phantom . . . . . . . . . . . . . . . . 936.13 Planar cross sections of centered spiral reconstruction . . . . . . . . . . . . . . . . 946.14 Segmentation rendering of physical phantom . . . . . . . . . . . . . . . . . . . . . 946.15 Knot clusters in dried douglas fir . . . . . . . . . . . . . . . . . . . . . . . . . . . 956.16 Combined pattern dry log slices - centered spiral . . . . . . . . . . . . . . . . . . 966.17 Sector pattern dry log slices- centered spiral . . . . . . . . . . . . . . . . . . . . . 966.18 Dry log major planes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 976.19 Dry Log segmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 986.20 Wet log centered spiral data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 996.21 Wet log major planes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1006.22 Wet log cross sectional slice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1016.23 Fresh log segmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1017.1 Plane sawn boards . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103x7.2 Dry log plane sawn boards . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1037.3 Fresh log plane sawn boards . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1047.4 Quarter sawn boards . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1047.5 Dry log quarter sawn boards . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1057.6 Fresh log quarter sawn boards . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1057.7 Surface corresponding to internal knot and crack . . . . . . . . . . . . . . . . . . 1077.8 Surface corresponding to five knot cluster and crack . . . . . . . . . . . . . . . . 1077.9 Fourth annulus voxel ring of five knot cluster . . . . . . . . . . . . . . . . . . . . 1088.1 Green log reconstruction with superimposed conical slices . . . . . . . . . . . . . 110xiNomenclatureα pixel-ray column inclination to xy planeβ material x-ray mass attenuation coefficient∆θ angular sampling spacing of Radon Transform∆d ray spacing of Radon samplingλ relaxation parameter of Kaczmarz iterateC n matrix column rotation operator by n placesDn Derivative operator with respect to the nth variableHn Hilbert Transform operator with respect to the nth variableR Radon transform operatorR−1 Inverse Radon Transform operatorΨ x-ray cone beam angleρ (x) one dimensionsal density functionρ j density function value of voxel jθ ′ coordinate of polar input to Radon Transform~ρ density vecotr for reconstruction~ytotal Radon Transform samples for all views~yview Radon Transform samples for one θ′xiic fractional scaler of innermost annular voxeldΨ angle between each equal-angular pixel rayDc diagonal matrix made up of column sums of GDL distance from x-ray source to log centerDr diagonal matrix made up of row sums of GDs distance from x-ray source to detector screenE radius of domain of the Radon Transform in the image spacef (r, φ) polar Radon transform operandf (x, y) Cartesian Radon transform operandf i threshold function for selective smoothingG path length basisgi j path length intersection of ray i with voxel jGsub basis sub component from pixel-column intersection with a given sliceGview_MS basis component from one view of multi-slice reconstructionGview_spiral the component of the path length basis from one view in spiral reconstructionGview basis matrix component from one viewI attenuated x-ray signal intensityI0 unattenuated x-ray signal intensityK Radon space ordinateL Radon space abscissam number of Radon Transform viewsn number of samples per Radon Transform viewxiiiNannuli number of annular voxel elementsNelem number of vertical detector elements on EMCCDNsector s number of sector voxel elementsNsl ices_view number of slices within the cone beam in one data captureNtotal total number of slices reconstructedNviews number of views around logO Radon space originP1, P2 furthest z coordinate extent of cone beam on logpi hyperplane described by the ith row of GPangnorm normalized log pixel-ray radiusPang number of pixel-rays in fan beamPcol number of pixel-ray columnsRn Hilbert Space of degree nri coefficients of the ith row of GRlog physical diameter of the logRnorm normalized log radiusRpix number of pixel rays subtended by the logt ′ length along Radon space abscissavi voxel value of ith index in smoothing clusterw smoothing weight for selective regularizationWD height of the detector screenWz new active detector widthxivyi sample of Radon Transformz′ coordinate along Kt slice thickness valuexvAcknowledgementsI would like to thank first and foremost my supervisor Dr. Gary Schajer, whose perspective andcountenance were always sure to steer my efforts in the right direction and with meaning. AlsoDr. Anthony An, without whom none of the work here would be possible. Finally, thanks tomy lab-mates for their support and suggestions, insights and skill which made the task at handthat much easier.xviDedicationFor my parents, whose wisdom is appreciated more by the year.xviiChapter 1Introduction1.1 Overview of log grading and scalingThe sawmill industry currently operates on an economic model in which significant funds arespent on the raw material for both softwood and hardwood operations. In any sawmill themost valuable asset is contained in the form of unprocessed material or felled logs, which oftenaccount for more than three-quarters of all expenses [1]. As such, there is a significant economicmotivation to maximize the amount of profit extracted from each individual log by optimizingthe type and quality of product made from it. In order to make the most of the raw material,sawmills today employ persons known as graders to both grade and scale the logs [2]. Gradingrefers to the process by which features of the wood such as size, shape, knots, cracks, etc.are taken into account in order to assign a specimen with a certain grade reflecting its overallquality. The process of scaling is then performed in which the quantity of said product able tobe extracted from the log is determined, and another label applied to the sample. Oftentimes,the grade of the wood may reflect certain minimum standards for scaling as well i.e. a log meantfor siding panels of a certain dimension must have some minimum diameter.Clearly, the downside of visual inspection is the variability between graders and the inabilityto detect features and defects that are not visible or indicated on the surface before the first cut.Rot, voids, and internal knots (limbs lost earlier in the tree's life and grown over) are all difficultto detect and may cause both the grade and scale of the log to be overestimated. This in turnwill lead to improper processing of the sample and maximum value will not be realized. Incorrectscaling and grading necessitates increased investment in quality assurance and monitoring sothat inferior product does not make its way out of the operation. The manipulation of logs1via roll out does not fit with the industrial process flow of the modern sawmill in which timelyprocessing is of the utmost importance [1, 3].1.2 Sensing technologies in sawmillsIt is estimated that the accurate identification of features within logs and subsequent smartprocessing may increase lumber yield value upwards of 15% [4]. The application of variousnondestructive sensing technologies in log grading and scaling has accelerated over the pastthirty years, with many technologies being developed and applied. These run the gamut fromultrasonic sensors to a wide variety of electromagnetic wave based systems, with the mostsuccessful implementations those that are both mechanically and computationally robust [2].The need for reliable and economic operation is paramount in the function of an industrial logscanner, and several technologies have developed. Furthermore, the technologies must operatequickly and determine the ideal product and/or orientation of a log in almost real time, as theflow of material on line may be as high as 3m/s and can not be impeded or profit is lost [5].The sensing technologies currently in use in modern sawmills for cut and uncut wood maybroadly be classified into external technologies (benefiting mainly the scaling process of logs)as well as internal technologies by which further information regarding the interior of pieces isobtained. Technologies that can tell something about the piece before the first cut maximizethe profit potential from the log and have the most promise. Of the external technologies,laser based topological scanners have been widely studied and are implemented in industry todimension logs and scale them from their overall form. These scanners are able to identifyexternally protruding knots from which inner information may be inferred and cut or gradingdecisions made, and are the most common sensors currently in sawmills [6].The internal sensing technologies that have been implemented include ultrasound, interfer-ometry vibrational sensing, microwave, nuclear magnetic resonance, x-ray and more. Of these,x-ray, and particularly x-ray computed tomography, has shown the most promise in being ef-fectively implemented for feature detection [2, 7, 8]. This is due to several factors, not leastof which is fresh lumber's x-ray mass attenuation qualities being close to that of the humanbody, allowing for much of the work done in the medical realm to be implemented towards anindustrial sawmill scanner. X-rays are also relatively inexpensive and easy to generate, and thusare a logical choice for a sawmill technology that is not only effective but affordable as well [2].Lastly, in the extreme industrial environment of a sawmill, x-ray is attractive as no physicalcontact need be made (unlike ultrasound in which various gels are implemented), it need not2be implemented near the log surface, and it is comparatively noise robust when compared totechniques such as interferometry. X-ray does have downside in the form of power consumptionas well has real safety concerns requiring shielding in an industrial setting.1.3 X-ray sensing in sawmillsMuch work has been put into x-ray technologies for application in sawmills, with a generalemphasis on application to hardwood operations in North America and softwood products inEurope. The area of x-ray log defect scanning may be broadly partitioned into single directionx-ray scanning, multiple source-detector arrangements (in effect partial computed tomography),as well as full-field computed tomography (CT) technologies. Single direction x-ray scanninginvolves 2D radiography in which a single x-ray source and detector are used to obtain informa-tion regarding the interior of the log, and is used in some sawmills today in conjunction withlaser surface scanning to assemble a representation of the piece's overall form [7]. Figure 1.1shows the scanning arrangement for this case as well as what a representative data capture maylook like.LogSourceDetector(a) (b)Figure 1.1: Linear RadiographIt's clear that this arrangement, while being rapid and low cost, also has some drawbacks.The inability to differentiate depth and the superposition of features make accurate measurementdifficult. This lack of clear demarcation of feature orientation has been the impetus for both themulti-source and CT systems. In larger systems, such as shown in Figure 1.2, multiple sourcesand detectors are used yielding additional information about the locations of knots and thelike [5, 9]. These systems have implemented up to eight source-detector pairings, together withthe consequent increased cost, energy consumption, and safety considerations of many x-ray3sources. Furthermore, they are often incapable of resolving all the needed details about theinternal log geometry. Seger et al. produced a two source-detector pair system comprising ofmultiple 1D detectors arrayed together and oriented at 90º with respect to each other aboutthe log's longitudinal axis. The reconstruction results of that system are displayed in Figure1.2(b) and are represented by a transverse cross section of the log. It shows that althoughknots may be displayed clearly, the sapwood heartwood details are lost as would be any othercircumferentially dominant feature.(a) (b)Multiple detectorsMultiple sourcesFigure 1.2: Multi-source detector system and output from Seger et al.The remedy for the limitations of the previous two systems has been available for some timein the form of computed tomography, as introduced by Hounsfield and Cormack in the late1960's and early 70's. The technology has matured profoundly since its first implementations in1971, almost entirely due to research in the medical sector. An x-ray source is revolved about aregion in space so that a detector is always in view, a minimum number of views are acquired asdictated by features such as the x-ray beam profile and reconstruction algorithm being used, andan inverse problem is solved so that a full reconstruction of the density variations is obtained.It is CT that has generated the most interest in the forestry research community for its 3-Dcapabilities. Figure 1.3(a) demonstrates this relative motion of the detector and source, and aCT reconstruction of the same transverse slice of 1.2 [5]. Here the heartwood sapwood boundaryis clearly visible.4(a) (b)Revolving SourceFigure 1.3: CT detector configuration and output from Seger et al.Still, there are many barriers-to-entry of CT technology to the conditions present and op-erational demands within a sawmill. Medical inspired log CT scanners require a great deal ofprecision and accuracy in reconstruction that causes them to be too slow for implementation ina real industrial line setting, too delicate for continued service (near continuous operation is areality for most major mills), and far too costly for many small-medium sawmills to purchase[1, 6]. These considerations are not serious impediments in the medical setting, but the lumberindustry is still waiting on a product that can fulfill these demands.1.4 MotivationRotating source-detector gantry CT sawmill scanners developed from medical-style CT scannersare challenging for sawmill use. In these machines, the rotating x-ray source or rotating source-detector assembly maintains power and data transfer via complex transmission rings, which willrequire attention for continuous operation. They are subject to corruption of the reconstructiondue to movement of the log within the scan area that may ensue from the relatively uncontrolledenvironment of an industrial sawmill. Although capable of reconstructions with ~1mm resolu-tion, they struggle to meet the time requirement for implementation in a sawmill (estimated tobe one minute to reconstruct a four meter long log buck) [7]. Furthermore, current offeringsare large, standalone systems that are capital intensive and do not present much flexibility forincorporation into current lines without a large amount of production line reorganization.With these factors in mind, it is the intention of this work to lead toward a computedtomography system that is well suited for use in sawmill-specific applications. It will accomplishthis by being mechanically and computationally simple, robust, as well as devised with cost inmind. The ideal product will feature some immunity to perturbations of the log during scan time5and deliver data packaged to make fast decisions on processing and cutting. Lastly, a premiumwill be put on the system's adaptability to different sawmill environments and variations in rawmaterial geometry, as it will require retro-fitting into existing mills.There are several key differences in logs as a reconstruction object from other conventionallyscanned bodies. First, because a log is a non-living rigid body, it may be manipulated inways not possible for medical CT applications. The idea to rotate the log itself and not theCT equipment is a possibility. The log will not be harmed by higher ionizing radiation levels,which makes algorithms that sacrifice reconstruction speed for low dosage not necessary. Lastly,much is known about the way in which logs of a given tree species grow. This means logs'overall shapes and the forms of their features don't change much overall. That fact may beexploited to make for faster CT reconstruction. The work here will present a CT algorithm thatis specifically constructed for scanning softwood logs based off of individual log's commonalitythat will use this idea. It is first necessary, however, to introduce the fundamentals of howcomputed tomography works.6Chapter 2Computed tomography overview2.1 Basics of computed tomography2.1.1 Computed tomography setupBefore formally introducing the computed tomographic reconstruction problem and its solutions,it is helpful to give a basic introduction. The purpose of CT is to determine the spatial locationsof features inside a body by using the penetrative power of x-rays. A typical measurementsystem comprises an x-ray source, an object to be reconstructed, and an x-ray detector. Thex-ray detector measures the x-ray intensity. The x-rays are partially absorbed when they passthrough denser objects within the body, so in effect a 'shadow' is cast on the detector at thoseregions. The spatially distributed nature of the detector means that the relative positions ofthese shadows within the object are known. Figure 2.1 shows this by plotting the spatial signalof the detector looking at body A from the side. As can be seen, the detector picks up on thepresence of the two objects B and C, but is unable to determine that there are two distinctbodies. It is necessary for another view, shown in Figure 2.2 to be taken to distinguish theirpositions. The rotation of body A means both objects are seen and block the rays. If the processis repeated, and the object rotated in the detector field for a number of intervals, the interiorof object A may be completely reconstructed.The example shows that the arrangement is capable of finding the location of the interiorobjects, but does not explain how the densities of these objects and the 'shadow' are correlated.This connection is found in the Beer-Lambert law or Beer's law from classical optics, and is afundamental building block of CT.7Figure 2.1: CT configuration showing detector signalFigure 2.2: Object rotated to discern bodies B, C2.1.2 Beer's lawWhen an x-ray beam encounters an object, it is partially absorbed. The Beer-Lambert lawstates that the ratio of the attenuated power intensity of a ray to that of the un-attenuatedray through the same space is equal to the exponent of the line integral of density along the8path divided by an attenuation coefficient β intrinsic to the attenuating substance. This isrepresented in equation form by 2.1.II0= exp(−ˆρ (x)βdx)(2.1)x(x),I I0Figure 2.3: Beer-Lambert LawFigure 2.3 illustrates this relationship for a single slab of attenuating material with theproperties ρ, β. The line integral of the density through the slab is what is of interest here, sothat the density of the material may be found. In order to find this, Equation 2.1 is solved forthis quantity by taking the natural logarithm of both sides and rearranging to yield Equation2.2. Note that in this case, the coordinate x is referring to the distance along the ray from itspoint of intersection in the direction of propagation.ˆρ (x) dx = −βln(II0)(2.2)It is the energy intensity of the rays that is measured by the detector. The constant β maybe determined by calibrating the system with a substance of known density and dimensions.Thus by taking a radiograph of the object in question to get I and an un-attenuated 'reference'image to obtain I0, an approximation of area density of the object is readily determined. Thisinformation can then be used to reconstruct the form of the object, as will be shown.2.2 Overview of CT terminologyComputed tomographic reconstruction machines may be broadly classified based on the geome-try of the x-ray source and construction of the detector assembly, as well as the relative motion9of the assembly to the object and the dimensional aspect of the detector used. When the termsthird and fourth generation CT are used, it refers to the physical arrangement of the x-raysource with respect to the detector. A fourth generation assembly features a rotating sourcewithin a static detector ring, while a third generation scanner has both source and detectorrotating together (Figure 2.4). For a third generation machine, the detector may be either flat(termed a 'flat panel' detector) or curved as shown.(a) (b)Figure 2.4: Typical CT measurement geometries(a) third generation (b) fourth generationThe earliest x-ray scanners used non-divergent x-rays and were called parallel beam. Inboth scanner arrangements of Fig. 2.4, the x-rays emanating from the source are divergent andform a fan beam. If the detector is multi-slice, a cone beam is made and the machine iscalled 7th generation. The movement of the object being scanned along the transverse axis ofthe source trajectory indicates whether the reconstruction is rotational or helical (also knownas spiral) CT. In rotational CT, there is no translation of the object, while helical featurestranslation. For cases in which the speed of the scan is a factor, helical cone beam CT is theclear choice for quick data acquisition. The algorithms change for rotational vs. helical, singlevs. multiple slice reconstruction, but broadly fall into either analytical techniques related tofiltered backprojection (FBP) or series expansion techniques commonly called algebraic recon-struction techniques (ART) in the literature. Figure 2.5 contains three illustrations of the CTarrangements used in this research. Figure 2.5 (a) shows a single slice fan beam flat panel 3rdgeneration measurement system, the fan beam results from only one line of detectors beingactive. With all detectors active, the system becomes a 7th generation multi-slice cone beamscanner as shown in 2.5 (b). Finally, introducing log translation relative to the detector yieldsthe 7th generation helical cone beam scanner arrangement of 2.5 (c).10(a) (b) (c)Figure 2.5: Outline of CT configurations(a) 3rd generation fan beam single slice (b) 7th generation multi-slice cone beam (c) 7thgeneration cone beam helical2.3 The formal CT reconstruction problem and solution2.3.1 The Radon transform and problem statementThis discussion will consider the case of parallel rays, as it is the most straightforward and theconcepts can be applied to more advanced cases found in other references. The fundamentalcomputed tomography inversion problem for x-ray attenuation has its foundation in what isknown as the Radon Transform, introduced by Johann Radon at the beginning of the lastcentury. Consider a function f (x, y) that represents some attribute about which more specificinformation is desired e.g. density in the case of standard x-ray computed tomography (Figure2.6). It is the case that the spatial (x and y correspondence) form of the function is whollyunknown, and must be determined from information contained within its in-plane 'views'. Inorder to define these views, the function is first represented in polar coordinates as f (r, φ), inwhich r =√x2 + y2and φ = arctan(yx), x = r cos φ and y = r sin φ. The density at a given point(r ′, φ′) is then defined as f (r ′, φ′).11yxf(r', ')'r'f(r, )Figure 2.6: Distributed density functionConsider now the line L = θ ′ in the plane. The angle θ ′ may be thought to correspond toa single view of the function f . Looking at Figure 2.7, it can be seen that a new coordinatesystem (t, z) may be defined corresponding to this view, with t denoting the length coordinatealong L from the origin to some orthogonal line K (the intersection labeled as point P) and zthe length coordinate along K .12yxLt'KPOr'tz(r, )z'Figure 2.7: The Radon Transform coordinate spaceThese new arguments (t, θ) of the Radon transform are related to the coordinates r and φ bythe relation of Equation 2.3. It maps each point in R2 into the appropriate Radon coordinate tfor 0 ≤ θ ≤ 2pi. Figure 2.7 shows how for a given point (r, φ), the coordinate t ′ corresponds toK for selected view θ = θ ′, this in turn will determine the value z′.t = r cos (θ − φ) (2.3)Thus f (r, φ) may be expressed f(√t2 + z2, φ + arctan(zt))and the Radon Transform de-noted by the operator R acting on f , R f (t, θ) is given by Eqs. 2.4.Rf (t, θ) =ˆ ∞−∞f(√t2 + z 2 , φ + arctan( zt))dz (2.4)Rf (0 , θ) =ˆ ∞−∞f(z, φ +pi2)dzThe physical interpretation of this formula is that for a single view, the Radon Transformyields the continuum of line integrals orthogonal to L (that is in the same direction as K) overthe distributed function f (r, φ) . Figure 2.8 demonstrates this by plotting these line integralsfrom for the arbitrary density function in the plane for line L.13Figure 2.8: Graphical representation of Radon transformThe domain of the Radon Transform is limited by the physical detector size so that itis confined to a circle about the origin of radius E. The picture reconstruction area is thustraditionally held as the square of dimensions√2E × √2E inscribed within the circle. Thedomain may then be expressed as −E ≤ t ≤ E and 0 ≤ θ ≤ 2pi. It should be noted that thetransform is periodic for one half-rotation for the case of parallel forward projected rays, andso R f (t, θ) = R f (t, θ + pi). The analytic form of the Radon Transform is continuous, but aswith any physical system of a finite number of detectors the sampling is discrete. A total ofup to n rays at spacing ∆d = E/n may therefore be sampled per view for a total of m viewsthat are spaced ∆θ apart. Figure 2.8 represents this by the discrete dots yi along the transform.The data from a single view ~yview ∈ Rn may then be concatenated with the samplings from allnon-redundant views to form ~ytotal ∈ Rm×n. This leads to the formal problem statement:Given the sampling ~ytotal, approximate f (r, φ)The actual samplings of the Radon Transform or the line integrals of Rf (t, θ) may befound using Beer's law. Referring back to eq. 2.2, it is easily seen that the equation is simplystating the value of the Radon Transform for a given t, θ and that in this case the coordinate14x is supplanted by the coordinate z from Figure 2.7. As was presented in 2.2, the two mostcommon methods of solving the problem are the backprojection family of methods and the seriesexpansion methods.2.3.2 Filtered backprojection overviewUnderstanding of the Filtered Back-Projection (FBP) method is important for differentiatingand recognizing the advantages of the method of reconstruction described in this thesis comparedwith the industry standard, so a brief introduction is presented here. Exhaustive sources onanalytical methods may be found by Natterer and Epstein [10, 11]. It has historically been thecase that the medical community has relied on the filtered backprojection techniques for themajority of implementations, because it is computationally less burdensome than ART for thesame reconstruction resolution.In brief, the filtered backprojection techniques may be thought of as three separate steps.The first is the forward projection of rays through the image reconstruction circle as detailedpreviously. The Radon signal for each view is then convolved with a filter, often a ramp filterin the frequency domain. This filtered signal is then backprojected through the image space forthe view, and the density function estimated at a point in R2 by taking the summation of thebackprojections from all pertinent views. Filtering is necessary in order to offset the 'blurring'effect resulting from the smearing that takes place during the backprojection. A graphicalrepresentation is presented here for visualization of the method.15Figure 2.9: Filtered Backprojection Steps Visualization(a) Forward projection produces the radon transform (b) Radon transform is filtered (c)Filtered transform is back projected thru image space (d) & (e) Process is repeated for mviews within pi to provide reconstructionIn the case of 2D beam geometries, this backprojection provides an analytically exact methodfor producting the inverse Radon Transform R−1 . The individual steps, presented withoutderivation, may be summarily described as first taking the derivative of the Radon transform(2.9 (a)) applied to f (r, φ)) with respect to the first variable (D), then the Hilbert transformHof this with respect to the first variable. This yields the filtered form of the forward projection,and it may be proved that this is equivalent to carrying out the aforementioned convolution (b).The backprojection operator B of Equation 2.5 (shown for arbitrary function q(r, φ)) may thenbe applied (c), and the entire operation multiplied by −(12pi)which is a step commonly referedto as normalization, which ensures that no single backprojected view is overrepresented. Theoperation is then repeated as θ ′ varies for m views in steps (d) and (e).Bq (r, φ) =ˆ pi0q(r cos(θ − φ), θ)dθ (2.5)Altogether, this yields Equation 2.6. This states that for a function in the xy plane that16is expressed in terms of its polar coordinates , the inverse Radon transform R−= −(12pi)BHDmay be determined such thatR−1 (Rf (r, φ)) = f (r, φ) and produce the desired estimation of thedensity function. These steps here are described in the spatial domain, but are usually executedin the Fourier domain in practice for computational convenience. Once the density distributionfunction f (r, φ) is determined for many (r, φ), it may be converted into the Cartesian coordinates(x, y) of the pixel picture space C ∈ R2.BHDq(r, φ)= −1piˆ pi ˆ ∞−∞D(t, θ)r cos(θ − φ) − t dtdθ (2.6)The filtered backprojection method has historically been favored by the medical communityfor its ease of implementation and computational efficiency. As the reconstruction is essentiallya summation of all backprojected views, easily calculable thanks in part to the simple naturesof the Hilbert transform and backprojection formula, it is very computationally efficient. Dis-cretization of the image space takes place 'at the end' of the method in the sense that the voxels(pixels of the image space with an associated thickness) are segmented and assigned value onlyafter summation. It is therefore an attractive option for CT. It does have some limitations offlexibility in application, to be discussed shortly.2.3.3 Series expansion methodsThe series expansion or algebraic reconstruction technique is a well established method, indeedit was the method used by Hounsfield in the first generation scanner. It is simpler to conceptu-alize than the highly analytic FBP method. This method may be considered to 'start' with thediscretization as opposed to ending with it as in the case of the filtered backprojection method.For the image space C with an arbitrarily oriented ray through it, the Beer-Lambert law ofequation 2.2 may be discretized so that the integrand on the left hand side may be representedas a summation of the segmented intersection path lengths gi j of ray i with voxels j = 1, 2, ...n2multiplied by the densities ρ j of each respective voxel. This discretized version of 2.2 may thenbe written as equation 2.7. Here, yi is again the sampled Radon Transform of Figure 2.8. Figure2.10 shows the ray passing through the reconstruction space, with its intersection path lengthwith voxel j emphasized.n2∑j=1gi j ρ j = yi (2.7)17nnray iFigure 2.10: Series expansion representationThis can be done for all rays detected in order to setup the system of linear equationspresented in 2.8. The matrix G contains the path lengths of all intersections of each ray with allthe voxels in the image reconstruction space. The right hand side is now the Radon transformsampling for all rays, ~ytotal . This large system of linear equations correspondingly has manyways to go about solving or optimizing the system in order to estimate a form of the densityvector ~ρ that best approximates the real form of the solution. Another way to look at it is thatthe column space of the matrix G represents a linearly independent basis for the image spacethat ~ρ acts upon to best approximate the reconstructed cross section.G~ρ = ~ytotal (2.8)In practice, the calculation of G for all rays at the prescribed directions θ is no small feat, anda significant amount of cost in computing power is incurred. The way that this has been handledin the past is by simply approximating the basis functions in G using different schemes to getas close to the actual path lengths as possible. The most rudimentary of these is simply settingthe value of all intersected voxels in the path length matrix as the length of one side of a voxelin the reconstruction space. More sophisticated techniques have been proposed, that take intoaccount the changing aspect of the voxels based off of the change in viewing angle and the ray'sinclination and spacing in the plane. It is from these non-exact bases that the technique getsit's given series expansion name, as the density function is being approximated via the series18expression of these non-exact bases. Exact methods have been found to be computationallyexpensive and slow to implement [12]. In all cases, it is the nature of G that it is highly sparse,overdetermined, and diagonally dominant. Equation 2.8 has been solved in practice applying awide variety of solution schemes, including direct solvers such as least squares, QR factorization,or singular value decomposition, as well as iterative methods ranging from classical techniquessuch as Gauss-Seidel to more recent approaches like GMRES or the conjugate gradient method.To ensure an optimization close to the true form of the density vector, the column space mustbe sufficiently well-conditioned to allow for practicable inversion. This consideration often playsa large part in the determination of the solving system used.2.3.4 Comparing methodsWhat the series expansion method benefits from in conceptual simplicity, it loses out on incomputational cost. In modern medical CT scanners where sub millimeter resolution is necces-sary for sucessful prognoses, the size of G can grow considerably, and the inversion can takea prohibitvely long time to implement in an industrial scanner. This is in large part why theFBP derived methods have proven so resilient as the industry standard. Still, the backprojec-tion methods also have concerns. It is necessary for the Radon transform to be taken regularlywithin an interval 0 to pi, ie the algorithm is highly sensitive to only a limited number of viewsbeing available. If there are not sufficient views, the reconstructed image will exhibit heavyartifacts. This is not as much of an issue with ART, it is much more resilient in working withless than 180◦coverage [13]. This is not a great concern for single slice reconstruction, but willbe significant for fast helical scanning. Also, FBP is more sensitive to noise in sampling datathan Algebraic Reconstruction.The series expansion methods also maintain an advantage in being easy to implement forany generalized source-detector geometry and configuration. This makes for a great deal offlexibility in defining the path length geometry as well as reconstruction geometry within thesource-detector plane as will be seen later. In the case for filtered backprojection ands its deriva-tive analytic methods, new algorithms must be developed as the case is changed from parallelbeam to divergent beam geometry, or if the detector shape changes (flat vs. curved). This isbecause the analytic form of the Radon Transform must be acquired to enable formulation of thebackprojection operator. In ART, these changes are trivial and simple geometric adjustmentsmay be done to enable Eq. 2.8 to still be used without major modification, so that the form ofthe Radon Transform need never be obtained. As the implementation of a log scanning system19into an existing production line necessities a retrofitting that may take many forms in geometryand componentry, the robustness of ART's adaptability is attractive [14]. Additionally, ARThas shown promise in several studies at yielding not only accurate spatial feature detection,but also a superior ability to accurately find voxel density values when compared to analyticmethods [14].Both algorithms in conventional use are similar in that they are implemented with the imagespace C being referenced to an absolute physical space. That is to say, the algorithms reconstructwhatever passes through the reconstruction space in the radius E circle of Figure 2.8. As thematerial to be reconstructed moves through space, all the recorded Radon transform arguments(t ′, θ ′)correspond to the stationary world coordinates(r, φ)with the origin O in Figure 2.7 ofboth systems corresponding to the center of rotation of the source-detector pair. This allowsthe CT arrangement to be applicable to a great variety of body shapes and sizes, as anythingmay be constructed so long as it fits in the source-detector circle and does not move betweenviews. Both algorithms have also been traditionally implemented with reconstruction spacebeing defined to comprise a Cartesian grid. With a fine enough voxel pattern, this allows thealgorithm to reconstruct irregular shapes to a high degree of accuracy, important in a medicalsetting.2.4 Three dimensional reconstructionThree dimensional CT, termed 7th generation for both the pure rotational multi-slice andhelical cases, may take one of two forms for both types of algorithm families. The first approachinvolves taking the two dimensional algorithms of the last section and applying them repeatedlyin the transverse direction. This may be done by an actual translation step before rotation (ef-fectively decoupling the translation and rotation from happening concurrently) or by translatingand rotating simultaneously and interpolating the pitch of the source detector plane helix intosingle reconstruction planes periodically. This last technique, commonly referred to as 180LI or360LI (LI here to designate linear interpolation) necessitates that pitch be tight enough forthe interpolation along the helix to be usable. These algorithms are typically used with eithersingle line detectors or two dimensional detectors with only a few lines arrayed together, as thetwo dimensional algorithms of Section 2.3 are only effective for small lateral angles. Multipledetector lines do help to relax the source-detector pitch required.Clearly, the extended 2D algorithms are suboptimal for applications in which speed of re-construction is a consideration, and the proliferation of wide-area x-ray detection technology20makes their use obsolete as they cannot function on these wide data sets. Also, the requirementfor a pitch close to the size of a reconstructed slice puts significant mechanical limitations onthe system in order to acquire data quickly. Fully volumetric formulations of the forward pro-jection of lines for cone beam geometries (analogous to the Radon Transform of the previoussection) and its inverse have been developed for both circular and helical source path cases. Asis the situation in 2D reconstruction case, the 3D analytical projection-backprojection methodsare computationally elegant, fast, and accurate when implemented correctly. Expectedly, thesame problems are present as were in planar reconstruction. Limited number of views nowcauses artifacts in both the transverse and longitudinal directions, and changes in geometryor scanning scheme now require fundamental changes to the underlying backprojection math-ematics [15]. This makes the algorithms highly suited to reconstructions in which the areabeing scanned and aspects of the reconstruction will be unaltered throughout service life and indifferent implementation.The series expansion method for fully three dimensional reconstruction is fundamentallyunchanged. The size of the path length matrix of eq. 2.8 is now much larger than before,but as long as a careful account of the geometry of the path lengths is kept, the inversionproblem may be solved as was done previously. Thus a system implementing ART in volumetricreconstruction may be quickly adapted to a variety of detector shapes and ray-pixel binningarrangements without careful consideration of the effects on the quality of reconstruction. Itis for this reason that the development of this work's novel industrial scanner the simplicity ofART is best suited.2.5 SummaryThis chapter has briefly outlined the major ideas that have advanced CT in the past forty years.It is often the case that computational speed and robustness are at odds with accurate, cleanreconstruction in an industrial setting. The work ahead aims to blend the best features of whathas been discussed to achieve the present goals. The mechanical simplicity of using one detectorand one source arrangement will be made feasible through manipulation of the log as opposedto the scanning hardware. The quick data acquisition of a multi-slice, helical cone beam setupwill allow for fast reconstruction. Lastly, the noise resistant and simple nature of ART can beimplemented rapidly with special considerations towards the shape of the features in the logsbeing scanned.21Chapter 3The feature-tailored voxel approach3.1 Overview of the method3.1.1 Method descriptionThe method proposed here is to implement a 7th generation helical cone-beam flat panelCT scanner utilizing series expansion reconstruction to detect knots and features in logs. Itwill do this utilizing several key features that differentiate it from systems that have beenwidely implemented in the medical industry and make it specifically for use with logs. First, itimplements an unconventional voxel pattern that closely follows the geometry of an actual log'sfeatures of interest. It also uses a scheme to detect, normalize, and center detector data to allowfor rigid-body motion in the reconstruction plane, as well as enabling fast basis construction andsolving times. In order to work up to this helical scanner scheme, it is useful to first consider asingle slice fan beam scanner and then graduate to a multi-slice rotational cone beam scanner.With these two cases established, the more complex helical scan can be developed.3.1.2 Feature-tailored voxel geometriesLooking at the log cross section of Fig. 3.1, it is clear that the features affecting wood valuein coniferous trees have some basic geometries. First, the logs are nominally circular in crosssection and have very mild eccentricities compared to hardwood specimens [16]. This circulargeometry is reflected in the form of the heartwood-sapwood boundary that separates the denseouter living part of the tree from the less dense mature heartwood. The knots are branches thatgrow outwards from the center of the tree over its lifetime. All knots start from the center of22the tree as the branches form successively at the upper part of the tree when it is still small.If a branch has been cut or broken off early during the tree's life, the knot in that area will begrown over and become internalized. Thus it is never the case that a knot begins in the interiorof the heartwood or sapwood and extends outwards without intersecting the center, but a knotmay start at the center and be truncated in the interior.Figure 3.1: Medical grade CT reconstruction of softwood specimen from SzathmaryThe first and most obvious alteration to the standard reconstruction methods it to utilizethe entirety of the reconstruction space as the image space. Instead of discarding the areasoutside image space square, the image is simply defined to have a circular boundary of radius E.This means that there is no wasted space in the reconstruction as the circular image containsthe log cross section entirely. The voxel geometries of the derivations of Sections 2.3.2 and 2.3.3are Cartesian in both 2D and 3D with each voxel being square and cubic respectively. Thismeans that it would take a large number of voxels to be able to reconstruct the shape of theheartwood sapwood boundary and knots within a log. If instead the voxels are concentric circlesor a 'bulls-eye' arrangement, it takes only a few voxels to be able to demarcate the heartwoodsapwood transition. A segmented circle pattern with sectors on the scale of the angular size ofa knot means that relatively few sectors may accurately reconstruct the position of all knots.These arrangements are shown in Figure 3.2. The logical step is to use a superposition of both23patterns in order to get the benefits of both, yielding the combined geometry.Figure 3.2: Voxel reconstruction casesThis polar voxel arrangement has been implemented and studied before in CT, by Jian et al.in the single slice case and Thibedau et al. for volumetric purely rotational reconstruction, butnot for helical lumber CT which presents unique opportunities due to the way in which woodforms [17][18]. While retaining the same amount of feature resolution for a given log section,the spatial voxel resolution may be reduced by multiple orders of magnitude to make seriesexpansion reconstruction practical, bringing all of its comparative advantages to bear to thesystem under consideration. Another big advantage is that the reconstructed image basis maybe formed by calculating the path lengths of each ray's intersection with the voxels. Becausethe voxel pattern is symmetric with respect to the origin, the basis calculation may be solvedquickly for one view and then reused, as will be shown later.3.1.3 Data binning, centering, and normalizationIn order for the entire reconstruction space to be equal in size to the image space, the datacollected must be manipulated so that no empty space is included in the reconstruction problem.This has the added benefit of making the system immune to small movement of the log in thedetector view since the reconstruction space follows the log.24Figure 3.3: Log setup in 2D detector fan-beam caseConsider the system shown in Figure 3.3, that of a flat panel 3rd generation fan beamscanner with a log in view. The distance from the x-ray source to the detector is taken to be Dsand the distance from the source to the center of the log is taken to be DL . The real, physicalradius of the log in the fan beam arrangement is Rlog and the detector is taken to be of heightWD divided into Nelem detector elements. This detector subtends the fan beam at an angle of2Ψ. Each detector element may be thought of as corresponding to a pixel-ray that extendsfrom the source straight to the detector. In the context of this research the detector elementsare linearly spaced and present in a number in excess of what is needed by the voxel densitypattern. It is therefore possible to re-bin the pixel rays that lie linearly along the detector intoPang pixel rays that intersect equal-angularly at an interval dΨ = 2Ψ/Pang at the x-ray source.The equal-angularly spaced pixel rays are the arrangement that would be measured with acurved detector of radius Ds . The relation found in eq. 3.1 details the linear spacing of theseequal-angular pixels along the flat detector, and the relation is shown in Figure 3.4. Binningthe rays like this will allow for the correction of small displacements parallel to the detector aswill be shown. This is important because the harsh nature of a sawmill makes it very difficultto manipulate and orient a log on the production line.25Di =WD2− Ds sin (Ψ − idΨ) (3.1)dWDDiyxPangFigure 3.4: Cylindrical binningBefore the log can be centered and the extraneous detectors eliminated, it is necessary tofind the location of the center of the log in the plane and approximate the log's size. The generalscheme for locating the log center within the detector view is to first approximate the geometriccenter of the log by the centroid of it's Radon Transform signal, then estimating the radius ofthe log by approximating it's Radon signal as a semi-ellipse. This method uses all of the dataand has been observed to be more robust in practice than edge detection methods, which involvejust a few measurements near the outside surface. The method is best understood for the caseof a uniform density log within the detector view. Consider the circle shown in figure 3.5 andits Radon Transform shown on the ordinate. Using the coordinates of 2.7, the density functionbounded by the circle with radius Rlog and constant density ρ′ is eq. 3.1.3.ρ (r, φ) =ρ′ r ≤ Rpix0 r > Rpix(3.2)For the view θ = 0, it is clear that the parallel beam Radon Transform in this case is equal in26form to the equation of an ellipse as the line integrals parallel to the y axis have the arc lengthvalues z = 2√R2pix − x2. In this setup's case this is an approximation as the non-parallel natureof the divergent fan-beam distorts the signal to be not exactly elliptical. Still, for small enoughvalues of Ψ the approximation holds well enough to give good results.x = ty = zRf(t,0) =(r,φ)Figure 3.5: Uniform density in a bounded circle and its Radon TransformThe relation between the semi-ellipse's centroid height and it's major axis dimension is givenby Cy = 4h3pi , and its area by A =piRpixh2 from geometry. The signal may be thought of as seriesof pixel bins each containing a sample of the Radon Transform. Using the additive nature ofsecond moment of area, the centroid may be determined as Cy =∑b2i2∑bi; the area of course beingA =∑bi . Manipulation of the equations leads to cancelation of h and eq. 3.3, readily providingan estimate of the radius Rpix . of the log. This may then be related via simple geometry to anestimate of the log's actual radius, Rlog , by eq. 3.4.27Figure 3.6: Relating elliptical Radon Transform width to heightRpix =163pi2(∑bi )2∑b2i(3.3)Rlog = Dl sin(RpixdΨ)(3.4)Figure 3.7 shows the log arrangement from before but now with a small vertical displacementof the cross section in the reconstruction plane. Accordingly, the measured Radon Transformis shifted on the detector as well by a small value e. If the displacement is small, it canbe approximated to be along the circumference of the circle of radius DL with the source atits center. Because the pixel-rays are now equal-angular in nature, the Radon Transform ofthe log may be found in the detector field and re-centered simply by shifting the center bye = Pang/2−Cx . In addition, the data outside of Cx ±Rpix may be truncated, causing the signalto take up the entire reconstruction space.28Figure 3.7: Vertical log shiftThis now leads to the consideration of horizontal motion as shown in Fig. 3.8. As the logmoves closer to the source, its Radon projection will take up a greater amount of the detectorarea and vice versa for movement towards the detector. To accommodate this situation it isassumed that for small angles of the diverging x-ray beams, the path length basis column spacewill not change considerably. It is then possible to scale the truncated data to a normalized casefor which the path lengths have already been computed for a user defined number of normalizedangular pixels Pang_norm. This normalized case in turn has its own associated 'real' log, witha radius of Rnorm. Since the basis weight (i.e. the Radon Transform sampling) has a directrelation to the linear size of the radius of the density function, expression 3.5 may be used toscale the data to the normalized case. The steps are shown in Figure 3.9. These steps areconvenient in that they allow for a predetermined path length matrix to be used, so that nopath lengths need to be calculated at runtime of the reconstruction inversion.29Figure 3.8: Horizontal log shiftFigure 3.9: Data normalization schemeRf (t, θ)scaled = Rf (t, θ)RlogRnorm(3.5)With the data conditioned so that the predetermined path length matrix may be used,the system of equations of eq. 2.8 may be constructed. The system of equations is highlyoverdetermined, with most of the work going into assembling the path length matrix G. The30steps used to calculate G change whether the reconstruction is single vs multiple slice, purerotational vs helical, and depending on the voxel arrangement used.3.2 Calculation of basis for a single sliceThe calculation of the basis for each of the voxel arrangements follows a similar process. First,the path lengths for the intersection of each ray and each voxel must be determined and assem-bled into a path length matrix Gview . The matrix arrangement is determined by the numberingof the voxels and rays, with the voxel numbering patterns shown in Figs. 3.10-3.12 and the raysin the fanbeam starting from the top, least to greatest as in Figure 3.4. The next step is toadvance the voxel arrangement by a rotation dθ = 2pi/Nsector s . The rays in these additionalviews add rows to the path length matrix, and their intersections with the voxel pattern areconcatenated onto the bottom of the previous view. This will be explored in more depth witheach case.3.2.1 Path length basis function geometry by voxel arrangement viewThe simplest case is the annular model. First, the annular voxel sizes must be determined in astandard way for a user chosen number of annuli, Nannuli . This is done by specifying a startinginnermost voxel radius as a percentage cRnorm of the normalized radius Rnorm of the standardlog view. From there equation 3.6 may be used to ensure that the remaining annular voxels areequal in area, which is a convenient way to standardize the method for determining their sizes fora desired number. Intersection path lengths are determined by considering a ray i intersectingwith the annular voxel j. These intersections are shown in Fig. 3.10 for a given ray by thedots along the pixel-ray line. The path length of a given ray in a voxel may be calculated bycalculating the path lengths within the outer circles and then subtracting the path lengths fromwithin the inner circles. Doing this for all rays Pang in the fan beam for one view is necessaryto establish the path length matrix for the first view's worth of equations.31Figure 3.10: Annular Case (a) Pixel intersection (b) Gview for annular caseRj =√( j − 1) R2norm(1 − c2)(Nannuli − 1) + c2R2norm (3.6)For the sector model, the path length intersections of each ray in a view with the voxelboundaries must again be determined. The total number of sectors is constrained here to beeven, as well as the number of ray pixels in the fan beam. This ensures that no pixel willfall along a voxel boundary and means that only half the path lengths must be calculated,and then their symmetric path lengths can be simply inserted into Gview . Figure 3.11 showsthe intersection of an arbitrary ray with this voxel pattern. Note that the voxel numberinggoes counterclockwise, with the first sector taken to be that one in the near vertical position.Calculation of the intersection path lengths is accomplished by defining a point O on the raythat intersects the origin orthogonally. For each sector, the line comprising the leading edge ofeach sector is recorded, giving the intersection point pj . The differences between the lengths li jgives the entries into the path length matrix.32Figure 3.11: Sector Case (a) Pixel intersection (b) Gview for sector caseThe combined model algorithm is expectedly a combination of the path length combinationsof the previous two, with each ray now having many more intersections than was seen previously.The voxel indexing now starts from the center and goes outwards from the innermost near-topsector voxel and continues outwards in the counterclockwise direction. Now path lengths ininner voxels are determined by subtracting the path lengths of the smaller inner sectored circlesfrom the outer larger ones. The voxel numbering means that the view matrix is a horizontalconcatenation of grouped sector matrices, with each submatrix having Nannuli columns, or onecolumn for each annulus.33ray i(a) (b)Sector 1 Sector 2 NsectorsPangNannuliFigure 3.12: Combined pixel ray intersection3.2.2 Path length permutations by viewIn the case of the annular voxel pattern, there is no shifting of voxels as the log rotates. Thismeans that the path length matrix G for Nviews is simply Gview stacked on top of itself Nviewsnumber of times (Figure 3.13).34NviewsGviewGviewGviewFigure 3.13: Annular single slice total basisThe case of the sector model is more complicated. The nature of the rotation is shown in3.14. After each data sample, the log is rotated away from the detector by one sector. Thismeans that the path lengths that were determined to construct Gview can still be used, but theindices of the values must change before the matrix can be concatenated on. The rotation inthe physical log corresponds to a rotation of the columns of Gview by one column to the left.This may be done for as many views as desired to fully constrain the system. This is shownin Figure 3.15, in which C n represents this matrix operation to the nth rotation. It has beenshown by Turbell et al. that for a fanbeam pixel arrangement, the CT data becomes redundantafter γ = pi + 2Ψ angular coverage [15]. More views than this do not introduce new information,but will serve to simply constrain the system more.35Figure 3.14: Sector rotationNviewsGviewC1[Gviews]CNviews[Gviews]Figure 3.15: Sector single slice total basisFor the combined log model, it is the submatrices that are rotated, that is each additionalconcatenation has a rotation to the left by Nannuli columns. This is shown in Figure 3.16.Again, only γ views/sector rotations are needed to get a fully redundant data set. Still, moreviews do not hurt the reconstruction and are easily accommodated into the path length matrix36for the single slice view as size is less of an issue and speed of scan is not yet a factor.NviewsGviewCNannuli[Gview]CNviews*Nannuli[Gview]Figure 3.16: Combined single slice total basis3.3 Calculation of the basis functions for multiple slices3.3.1 Multi-slice basis function geometryThe case of multi-slice reconstruction in which path length calculations carried out on a 3Dmodel instead of on a single plane is more complicated. The largest difference here is that nowthe pixel-rays may pass through more than one slice at a time, meaning that the slices areinterconnected in the subsequent matrix inversion. This makes the inversion much larger andsusceptible to instabilities.37Figure 3.17: 3D CT arrangementThere are now two dimensions that the pixel rays extend into, as shown in Figure 3.17. Thecoordinate space is arranged with the origin at the x-ray source and the x, y, and z axes asshown. The detector is now taken to be a 2D plane at a distance Ds along the y axis. Thetotal number of pixels Pview in a given image capture is now equal to Pcol pixel columns in thez direction of the log and Pang pixels that are in the x direction. While the angular pixels areagain cylindrically distributed, the longitudinal pictures are linearly spaced in the z direction,so that the pixel columns may be thought of as extending outwards like pages of an open bookand form a pyramidal shape with the detector screen. The total number of log slices that areencountered by the cone beam in a view is taken to be Nsl ice_view . The height of the detectoris still WD, but the z dimensional width is now Wz , which is less than or equal to WD. Thisdimension is dependent on the reconstructed slice thickness t and Nsl ice_view . The dependencyis shown in Figure 3.18. As can be seen, the cone-beam intersects the log (taken here to bea regular cylinder of radius Rlog , estimated from a view as detailed in Section 3.1.3) with thefurthest z coordinates being located at points P1, P2 on the detector side of the log. Using similartriangles, the value of Wz may be determined. Note that for the extreme reconstructed slices on38either side of the log there are areas in which the cone beam will never touch, shown in red. Theslices that contain untouched voxels are shown in light orange. These conical regions will notbe included in the inversion, and must be kept track of in order to eliminate the correspondingvoxels from eq. 2.8.Figure 3.18: Voxels intersected in multi-slice reconstructionThe algorithm for the intersection of each pixel column with each slice makes use of the sameprocess as the single slice case, but with some key differences in its operation. Symmetry is usedby necessitating that that Pang , Plong , and Nsl ice_view are all even so that only one 'quadrant'of path lengths must be calculated, and then mirrored across the x and z axes for the other paths.First, a logical matrix Glog is constructed that records whether each pixel column intersects agiven slice. If a slice is intersected, the parametric coefficients of the lines that correspond tothe pixel-rays in that column are solved for so that the x and y coordinates of both entrancesand exits of the rays with the slice are found (Fig. 3.19). These point's coordinates are thenstored in an array. The block in Gview that corresponds to the pixel column's intersectionwith the slice is then determined by calling the corresponding path length algorithm from the39previous section, now with an input that keeps track of the ray's new entrance and exit into theslice and alters the path lengths accordingly. The ray's inclination α with regards to the x-yplane is then found and the path lengths found previously are scaled by Gi j_new = Gi j/ cos αto account for their oblique nature.Figure 3.19: Multiple slice ray entrance and exitIt is the case that the size of G becomes prohibitively large in multi-slice reconstructionwhen stored in full. It is necessary then to utilize sparse matrix storage in order to minimizecomputational waste. This scheme is implemented with MATLAB 2013b sparse tools in thiswork. Here, all zero elements in the path length matrix are not stored, but only non-zeroelements. The non-zero's have three associated arrays, one with a registry of row indices, onewith column indices, and another containing the path length in double precision. This cutsdown greatly on the amount of storage required and keeps the program from unnecessary zeromultiplications.3.3.2 Permutations by view for multi-slice rotational reconstructionThe path length matrix permutations for purely rotational multi-slice reconstruction are similarto that of single slice. The main difference is that now the path length matrix Gview_MS for40one data capture is made up of numerous submatrices Gsub , one for every pixel column-sliceintersection. Each Gsub for the annular, sector, and combined cases is identical in size to therespective case of Gview from the single slice cases in Figs. 3.13, 3.15, and 3.16 and is organizedsimilarly. Figure 3.20 shows how the Gsub 's make up the total view matrix Gview_MS . Thevalue Nvoxel s_sl iceis Nannuli for the annular voxel case and Nsector , Nsector ×Nannuli for sectorand combined reconstruction respectively.Slice 1 Slice Nslices_viewPixel Col 1Pixel Col PcolGsubFigure 3.20: Path length matrix for 3D view Gview_MSAfter Gview_MS is constructed, it may again be reused in such a way as to make assembly ofG quick and memory friendly. G for the annular case is just a vertical repetition/concatenationof the Gview_MS matrix. The rotational operations of the single-slice sector and combined voxelcases are again used, but are now enacted on the sub-matrices Gsub that make up Gview_MS .41For all voxel cases, the entire path length matrix is constructed as though every voxel in theNsl ice_view set were intersected, paying no attention to those voxels that are within the red zoneof 3.18. After it is assembled, the zero columns of G are truncated, as well as the correspondingvoxels in the vector~ρ. Those intersected voxels in the incomplete slices (the orange slices of Fig.3.18) are included in the system inversion but are discarded afterwards.Gview_MSC rot_k[Gsub]Concatenate for k=N_viewsk=0k=N_viewsGview_MSFigure 3.21: Multi-slice rotational G construction423.3.3 Permutations by view for multi-slice helical reconstructionIn helical multi-slice reconstruction, the total number of slices reconstructed Ntotal is greaterthan the Nsl ices_view contained in one image. The value of Ntotal is determined by the physicallength of the log L and the slice thickness t, so that Ntotal = Lt . Instead of thinking of each ofthe Nsl ices_view positions of 3.18 as a real, physical part of the log in space, it is easier nowto consider it as a register into which a given slice of the log can enter and exit between datacaptures. In order to keep track of the log ends as it moves into the cone beam, it is takenthat only the first register is filled in the first image, leaving the rest empty. From there, thelog advances in the -z direction one slice thickness at a time and one sector's rotation at atime. The last image is taken when slice Ntotal is in the last register, meaning that all sliceshave received Nsl ices_view number of views. Figure 3.22 illustrates the manner in which the logmoves into and out of the cone beam. The helical reconstruction complicates the computationof the path length matrix as there are pixel columns on both ends that will never intersect areconstructed voxel, shown in red, as well pixel columns that intersect the log as well as excessspace in the view area. This buffer space is shown in orange in Figure 3.22, and is handled aswas done in the rotational multi-slice where it is computed but then discarded. This must betaken into consideration in constructing G, but as was done with the path length matrix of purerotational multi-slice, it is first assembled as though every pixel-ray captured by the detector isto be included in the inversion.43Figure 3.22: Motion of the log in view registers(a),(b), & (c) motion of the log into the viewspace (d) the log in view (e) & (f) log exiting viewThe memory saving technique of reusing path lengths from Gview_MS can be extended44into another dimension now that the log translates one slice per image. The first step is toconstruct the first data capture's matrix Gview_spiral . This is essentially a matrix of dimensionsPang × Pcol rows by Ntotal × Nvoxel s_sl ice columns that is identical in the first Nsl ices_viewcolumns to the values of Gview_MS for the same geometry. The rest of the matrix is taken aszero because none of these slices are in view; the form of the matrix is presented in Figure 3.23.All zerosFigure 3.23: Representation of Gview_spiralNow for each additional view of the log, two nested rotations are enacted, one that acts onthe entirety of the column space of Gview_MS with form C −nk∗Nvoxels_slice for n slice registerstranslation between data captures on the kth concatenation, and one that acts only on theGsub 's as before with column rotation rot_k corresponding to the sector and combined casesas for multi-slice rotational reconstruction. The nature of this matrix manipulation is shownin Figure 3.24. After the full path length matrix G is formed, it is necessary to truncate theempty pixel columns. A script is run taking Glogic as input and from it determines which pixelindices are not contacting a reconstructed slice for the entrance and exit of the log into the conebeam. The rows corresponding to these columns are eliminated, making for a large part of thecolumn space of G to be null. This column space is trimmed away to bring the correct numberof voxels into the inversion problem. After the inversion is completed, voxels that contain data45outside the actual log (the orange slices) may finally be discarded.All zerosAll zerosAll zerosConcatenate for k=NviewsC rot_k[Gsub]C -nk*Nvoxels_slice[Gview_spiral]k=0k=1k=NviewsFigure 3.24: Helical basis matrix G construction463.4 Advantages of the methodThe traditional method of ART takes the complete path matrix as a logical matrix containingsome average length dimension in those voxels encountered by rays. It lacks accuracy and isundesirable in a system already subject to reduced resolution. Exact path length ray tracingschemes developed by Siddon and improved upon by Jacob have been available for some time,but are costly for large cubic voxel projections [12, 15]. Fortunately, the symmetric nature ofthe voxel pattern helps to make accurate basis calculation fast and relatively simple. Jian andThibedau have shown that utilizing polar series expansion methods has led to speedups by afactor ranging from five to twenty in basis construction times [17, 18]. Since the basis may bestored prior to runtime, this is not a critical factor but is most welcome.One of the major advantages of the system is how the data manipulation enables an accuratereconstruction even with log motion in the cone beam. Traditional CT systems have no suchmechanism and thus will present a blurring artifact in any region that was displaced even slightly.This in turn will lead the reconstruction greatly distort the reconstructed shape and location ofany features of interest or lose them. In the uncontrolled environment of a sawmill, vibrationswill abound and mechanical manipulation of the log will not always be optimal. Differencesbetween log sizes will make the data location scheme valuable as no physical adjustment of thesource or detector is needed to compensate for a bigger or smaller log. The centering methodpresented here is robust when compared to standard edge detection methods (Sobel, Canny,etc.) as the entirety of the CT data is used to determine the geometric center quickly and withless processing. This is important in a sawmill, as loose bark/detritus, bulging knots, etc willlead to ill-defined edges of the actual log in frame.The other half of CT reconstruction is the problem of segmentation, or the identification offeatures within a three dimensional data set in order to perform analysis or act. This in manycases may be just as large of a problem as the CT reconstruction itself. In a traditional 3D dataset, features may occupy voxels in a haphazard way so that the indexing of voxels containingparts of the feature are not numerically convenient. As an example, consider the case of Figure3.25 in which a reconstruction is contained in Cartesian 3D with the voxels numbered by columnand then by slice extending into the page. The log contained in the reconstruction space spansmultiple rows, columns, and slices of voxels in a way that is unknown. Thus any segmentationsoftware must look through data that may be large indices of separation apart from each otherand recognize that this is the same feature of interest. This requires special algorithms that atthe very least must rearrange any data vector of density voxel values into slices as a starting47point, locate the feature from there, and then use the directionality of the detected feature toplace it in reference to the reconstructed body's overall orientation.512 by 512 sliceVoxel jVoxel j + ~500Figure 3.25: Cartesian segmentationFor knot detection and segmentation from this work's process, the problem is less daunting.The voxel indexing of Figure 3.12 means that a knot will occupy a space of closely numberedvoxels for a thoughtfully scaled pattern (meaning the sectors are close to the physical size ofan average knot). This will prove useful in allowing for quicker identification and orientationof the knot in the log, making decisions of how to rotate the log before cutting faster. Theadvantageous ordering of voxels closely within features will be explored later in full for severalreal reconstruction cases. The cropped knot image from before is shown in Figure 3.26, nowwith an overlaid polar voxel pattern into which the same density data would be reconstructed.Clearly, it will be easier to identify areas of higher density in these larger numerically adjacentbins.48j+1j+2j+3Voxel jFigure 3.26: Combined voxel numbering breakdownLastly, the helical scan system is very flexible in its construction when compared to thatof commercially available machines. The nature of the 'dual' rotation of the matrices inGview_spiral along with the robustness of ART means that the total path length matrix maybe modified for any amount data capture intervals, so long as the pitch length of the log's move-ment remains unchanged and the amount of translation and rotation is accounted. This aspectwill be helpful in future work taking on the challenges of log manipulation.49Chapter 4Experimental setupThe setup used to implement these algorithms is an upgraded version to the one used by An, andconsists of a sample of log approximately one meter long mounted to a railed carriage assemblywith shafts, which allows for both rotation and translation of the sample within an x-ray conebeam. The cone beam projects thru the log onto a detector box that is parallel to the track andlocated at its center. This means that the source and detector can remain stationary and onlythe log carriage assembly needs to move. The detector interfaces with a PC that triggers datacapture by monitoring log rotational position via an optical encoder. The PC is also able to drivethe log's rotation in the cone beam via a stepper motor and timing belt connected to a pulleyon the log shaft. A different configuration using two stepper motors allows for translation, ortranslation along with rotation. These may be made to stop on data captures for the exposuretime. This PC as well as the x-ray generator, cooler, and controller are located behind leadshielding to protect the operator during tests. Figure 4.1 presents a plan view of the assembly,showing the overall layout of the equipment in the room for helical scanning. Figure 4.2 showsphotos of the setup, with the various components labeled.50Figure 4.1: Experimental Setup Layout51(a) (b)DetectorCarriageLogTrackShieldingX-ray tubeFigure 4.2: Log scanning setup(a) view of detector (b) viewing towards source4.1 Detector assemblyThe x-ray detector for this research is the same used by An. et al. in their experimental trials[8, 19]. The central part of the x-ray detector used in this research is an Andor iXon 897 electronmultiplying charge coupled device (EMCCD) camera. This camera uses a quantum cascadingof electrons in its pixel registers that make it much more sensitive than a standard CCD camera.It is capable of detecting single photon counts on the sensor surface while operating at very lowtemperatures. The sensor itself comprises 512 by 512 pixels. The camera is contained in adetector box that is carefully fitted to block entry of any ambient light into the detector box.The system is novel in comparison to conventional CT systems because it makes use of a staticdetector and accordingly a static source. This will necessitate rotation of the log as will bediscussed, and presents some advantages. An image of the camera setup is shown in Figure 4.3[19].52Figure 4.3: Andor iXon 897 layout from AnClearly visible are the camera connectors and apertureOn one side of the box there is a 60cm by 60cm gadolinium oxysulfide scintillator screenfrom Kasei Optonix that is in the path of the X-ray cone beam. This screen was selected froma number of options in the research of An and was found to be the brightest for a given x-rayenergy level, while minimizing blur [19]. The camera is oriented at 90 degrees to the scintillatorscreen and inside a lead-shielded box so that no x-rays may directly enter the aperture andcreate optical noise. A mirror is used so that the camera gets a central view of the scintillatorscreen. Figure 4.4 shows the detector box with its door opened as well as the scintillator screen,Fig. 4.5 shows a diagram of the box interior.53(a) (b)Figure 4.4: Detector box(a) door open (b) scintillator screenlead shieldcameramirrorscintillatorFigure 4.5: Detector box diagramDue to spatial constraints, the camera is very close to the mirror and a wide angle lensis used in order to capture the entirety of the scintillator. This leads to barrel distortion ofthe captured image, which is easily corrected after data capture. A large black and white54checkerboard was placed centrally and normal to the screen's lateral and vertical dimensions.The highly contrasting corner points of these squares are easily located via software, and theirundistorted locations solved for based off of the screen's size. A standard barrel correctionalgorithm was then be implemented to correct the image.(a) (b)Figure 4.6: Barrel correction(a) before (b) after barrel correctionThe detector assembly is crude when compared to medical-grade flat panel systems usedtoday, but presents some advantages. Apart from the obvious advantage of cost, a detectorof this type presents interesting possibilities due to its ability to scale. By using a differentscintillator and setup geometry, the system can be made larger to accommodate logs that areof a greater diameter. This in turn makes the detector more cost advantageous as scaling a flatpanel system requires the purchasing of multiple detectors. This is mainly due to the lack oflarge area flat panel detectors currently available on the market. Finally the fact that the lateraldimension of the detector is 'wide' (capable of capturing many slices in a single image' lightensthe mechanical burden on the system by making very fast rotation of the log (or a very smallpitch in the case of helical reconstruction) less necessary for rapid scans. This is important asmanipulation of a full sized log is a challenge unto itself. The static detector is desirable thoughbecause it greatly cuts down on the mechanical complexity of the CT hardware. Traditional 3rdand 4th generation derived systems are often in need of maintenance in order to upkeep theirpower and data transmission components. This is not necessary if it is the reconstructed body55that is rotating.4.2 X-ray source specifications and shieldingThe x-ray system used consists of a Phillips MG-160 x-ray generator/control panel combination,a Comet MXR 160 x-ray tube, and a Bernard water based cooling system. The generator iscapable of producing x-rays up to 160kVp while x-ray tube current is adjustable up to 10mA.The x-ray tube is specified to work in this range ast well. The focal spot of the x-ray tube isspecified as 0.4 - 1.5mm, with emanating cone beam subtending an angle of twenty degrees.The tube is shown in Figure 4.7.Figure 4.7: X-ray tubeAperature and shutter visibleFor the basis determined in Section 3.3 to be valid, it is necessary to get the x-ray tubealigned as closely as possible to normal to the plane of the scintillator and centered on itssurface. It is also critical that the x-ray cone beam be normal to the longitudinal axis ofthe log's translation, and so there are two alignments that are necessary. Alignment with thedetector was done roughly by first making a small jig to hold a laser pointer that would interfacewith the tube and could then be placed flat on the front of the detector box in its middle (Fig4.8). When the laser coincided with the same source points on the tube and box, the tube was56close to being in alignment. Further centering of the tube on the detector as well as adjustingthe main axis of the track was accomplished via trial and error. A 'hooped phantom' centeredon the detector and held on the log carriage's shafts made exaggerated perspective errors inimages apparent. Between images, the source could be manipulated slightly until satisfactoryalignment was obtained (Fig 4.9).(a) (b)Figure 4.8: Rough alignment via laser(a) source (b) detector57Figure 4.9: Sufficient fine alignment via perspective imageThe samples used in the experiments are of a density that necessitates voltages in excessof 100kVp and currents over 5mA. In order to protect from radiation back scatter, the setupwas assembled in a basement room in which measurements were made and a shielding planwas designed. Calculations outlined by National Council on Radiation Protection document147 on the safe construction of x-ray systems led to the setup of shown in part (b) of Figure4.2. Three of the four walls of the room were solid concrete of sufficient thickness to require noshielding. The fourth wall was composed of non-structural gypsum board, and so seven panelswere constructed of lead shielding measuring 1.6mm in thickness. University Radiation RiskManagement safety personnel were consulted and implemented appropriate instrumentation toconfirm that the safety shielding was sufficient.4.3 Log carriage and locomotion4.3.1 Log carriageThe rotating log is mounted on shafts attached to two carts that ride along a wheeled track,shown in Figure 4.10. The log sample must be a circular cylindrical object approximately onemeter in length and around 30cm in diameter. The track is made up of an MDF base on which58two 5m aluminum u-channels run. The x-ray detector and source are placed transversely to thistrack arrangement so that the shafts are roughly in the center of the scintillator view area. Thelog is mounted to the shaft via two MDF plates that have four 5/16 mounting holes that maybe used to attach the log sample via lag screws or bolts. To account for misalignment betweenthe center of the log and the two shafts, low speed universal joints are used between the MDFplates and the shaft. These joints are capable of handling up to 15 degrees of misalignment andensure that the log does not bind. Figure 4.10 shows the arrangement of the log carriage systemand track as well as the mounting blocks and u-joints.(a) (b)Figure 4.10: Log carriage(a) carts visible (b) u-joint and mounting plate4.3.2 Belt drive rotation and helical movementThe setup uses stepper motors, pulleys, and timing belts to rotate or rotate plus translatethe log sample in rotational and helical movement respectively. In the case of purely rotationalmovement, one stepper motor directly drives the pulley shown in 4.10(b). The carriage assemblyis moved so that the features of interest are within the detector view for purely rotational multi-slice and centered if single slice reconstruction is desired.For the case of helical motion, two 1.8 degree resolution size 34 stepper motors are employedto translate the log as well as provide rotational movement. This is done by linking the cart totwo separate timing belts, one that provides simple translation and the other imparting rotation59to the shaft pulley. These components are shown in Figure 4.11. The locomotion timing belt ishighlighted in the dashed green line, and the rotational belt in the solid red line. The rotationaltiming belt is not connected to any motor directly, but simply rides along the pulley systemhighlighted in the figure. This mechanically links the amount that the log rotates to how far ittranslates, eliminating the necessity for a control system to monitor this. Stepper motor controlis handled via an ATmega 103L 8-bit microcontroller that provides constant speed options fromω =.1rev/s to 1.6 rev/s on 18 tooth 3/8 pitch pulleys [19]. The diameter of the shaft pulleyDp determines the helical pitch of the log's motion to be it's circumference p = piDp and whencoupled with the resolution of the sector voxel arrangement the reconstruction slice thicknesscan be determined by t = 2piNsector s Dp . In this research, two pulleys were used to study theeffects changing the log pitch and thus voxel resolution - one 63 tooth pulley and another 36tooth pulley. Both pulleys had standard 3/8 pitches to interface with the 1 HTD timing beltproperly.(a) (b) (c)Figure 4.11: Timing belt conveyors(a) far side (b) helical belt (c) near side4.3.3 Control scheme, hardware, and sensingFigure 4.12 shows the control diagram implemented in linking the locomotion, sensing, and datacapture systems [19]. A US Digital S1 shaft encoder with 720 count/rev resolution was used tomonitor how much the log rotates (Fig 4.13). The user is able to specify an angular intervalupon which the camera is triggered for a data capture, the camera exposure time and electronmultiplying gain of the EMCCD, and whether or not the system should delay log movementduring capture and for what length of time. This angular value corresponds with the desiredangular dimension of the sector voxels that the log will be reconstructed into. This rotational60based control scheme eliminates the necessity to monitor the log's movement in time as it isentirely position dependent, a large upside for simplifying the control scheme. Data acquisitionis via a US Digital USB4 DAQ that reads the encoder and provides control of the ATmega103L. All tests were run via a Desktop PC featuring an Intel Core 2 CPU at 2.89 GHz with8GB RAM.Figure 4.12: Log scanning control scheme from AnFigure 4.13: Encoder on shaft614.4 Test specimens4.4.1 Physical phantomInitial scanning was done utilizing a physical phantom or a body of known density distributionto provide a known entity against which to compare reconstructions. This physical phantomwas constructed to have features of interest nominally resembling those that would be found ina log. The phantom consists of a doubled 12 (32.5cm) diameter concrete forming tube thatis reinforced by fiberglass strapping tape and approximately 90cm long. Inside this tube is aconcentric 6 (16.2cm) diameter concrete forming tube that has runs the entire length. Thecardboard tube material has a density of approximately 700kg/m3. The center of the phantomconsists of a 1.5 square wooden shaft that is imbedded into two endcaps comprised of doubled3/4 ply, providing the phantom with structural strength resistance to torsional shear. Thisleaves 84cm of length in which to insert features. Dried and compacted wheat kernel was usedin the innermost cylinder as a filler, having a nominal density of 820kg/m3. In the annularcavity between the outer and inner tubes, damp cedar sawdust was incrementally poured andtamped until it filled the entire area with a density averaging about 322kg/m3. These fillermaterials were packed around a wedge shaped cylinder of blue housing insulation foam, meantto emulate a radial 'crack' running the length of the log and with a density of 25.5kg/m3. Alongwith the crack feature, three separate planar 'knot cluster' features were inserted into the tubesat seperate locations. These clusters are made of sector shaped pieces of 1 medium densityfiberboard (density 760kg/m3) subtending on average 15 degrees and are arranged as is shownin Figure 4.14. The radial coordinates are given in Table 4.1 The images in Fig. 4.15 showsone such knot cluster as it was laid into the phantom as well as the phantom tube. Thesewere measured by marking the outside of the tube as the knots were laid into it, and thenmeasuring the markings circumferentially from a reference on the outside (the positive x axis inthe diagram).62A B CFigure 4.14: Phantom dimensionsKnot Slice A Slice B Slice C1 99 0 2122 200 78 3183 240 198 3534 332 229 //5 // 300 //Table 4.1: Angular θ coordinates of phantom knotsAll values in degrees from x axis63Figure 4.15: Log phantom4.4.2 Tree specimensTwo real tree specimens were scanned and reconstructed as well. The samples included one dryspecimen as well as a green specimen which is representative of what would be encounteredin a sawmill. Both specimens are of the Douglas fir species The dry specimen measured 90cmlong and ~27cm in diameter, while the fresh specimen measured 80cm long and ~28cm indiameter. The dry specimen featured several knot clusters as well as a longitudinal crack fromrapid drying, unlikely to be present in an actual sawmill log except at the very ends. The greenlog was sealed completely upon cutting using PVA wood glue to ensure that no such splittingoccurred. The dry log was much less dense on average and was used to simply confirm theexperiment's effectiveness at identifying log features, as good data was easy to collect. The wetlog was used to show the setup is still valid for highly dense fresh logs as would be seen insawmills. Figure 4.16 contains a picture of both log samples, clearly showing the crack in thedry log and the externally visible knots on both samples.64Figure 4.16: Douglas fir samples (left) dry log (right) green log65Chapter 5Single slice reconstruction5.1 Direct solution via the normal equationsIn the case of single slice reconstruction, the path length basis matrix is relatively small (onthe order of 103equations to solve). It is the case that it is still well conditioned enough thatdirect methods may be utilized to solve equation 2.8. Here, the choice has been made to utilizethe traditional least squares approach. It may be shown that the vector−→ρ ′ that minimizes theerror of the residuals squared of a model J~ytotal − G~ρ′K2, that is the least squares solution, isthe solution of the normal equation for the system:GTG−→ρ ′ = GT−→y total (5.1)As is clear, the normal equation is simply constructed by taking both sides of the systemequation and left-multiplying by the transpose of the basis matrix, GT . It is the case thatGTG is a diagonally dominant positive-definite sparse matrix of size Nvoxel s_sl ice. To solvefor the vector−→ρ ′, Cholesky factorization is employed as it is applicable to this type of matrix.Equation 5.1 may be written as LL∗−→ρ ′ = GT~ytotal , in which L is a lower triangular matrix andL∗is its conjugate transpose. The system may then be readily solved in O(n2) steps and thedensity vector found. For rapid solution, the factorization L may be saved for the normalizedparameters and used for further CT data at a later time.665.2 Synthetic phantom reconstructionTests of the algorithm were conducted using synthetic data to evaluate the center finding schemeof Section 3.1.3 and the subsequent inversion. Density vectors for several phantom log crosssections were built so that they had distributions similar to that of what would be found in reallogs, similar in concept to the ubiquitous Shepp-Logan used for brain imaging. Some examplesare shown in Figure 5.1. Different cases were tested including purely annular/circular features(a and b), origin-symmetric sector cases (c), and asymmetric cases with combinations of features(c and d). These grayscale images correspond to density distributions with Nsector s = 36 andNannuli = 18. Voxels with high intensity are shown in lighter shades, and less dense regions aredarker. The images shown here correspond to realistic log densities ranging from 350kg/m3to550kg/m3and are scaled on a 255 grayscale as a fraction of 650kg/m3. The synthetic log physicalradius was taken to be .17m, to reflect a realistic value close to the physical sepcimens.67(a) (b)(c) (d)(e)Figure 5.1: Synthetic phantomsThese density vectors were multiplied by an appropriate path length matrix G correspondingto the voxel case, Pang =92, Rnorm =.1667m, Ds= 2.1m and Dl = 1.67m to form a realisticsynthetic Radon transform vector ~ytotal . This vector was segmented into Nviews = Nsector s =36vectors of length Pang (values that were set by the size of G). These vectors were then insertedinto zero arrays of size 128×1, simulating a 4×1 binning of the 512 element sensor to be used inreal trials. the values Rnorm = Rlog and Pang= 92 were chosen as they made it unnecessary toreverse the steps of Figure 3.9 (it was determined that a log of radius .1667m would subtend 92pixel-rays in a 128 equal-angular fanbeam). These synthetic data arrays were then fed back into68the center-finding/normalization algorithm and assembled into a vector ~ysyn. The inversion wasthen completed using the method as discussed in Section 5.1.It was found that for the phantoms not displaying symmetry with respect to the origin,the center finding algorithm would incorrectly displace the center of many views by a smallamount (on the order of 1-2% in the worst cases). Although small, this center displacementwas enough to cause problems in reconstruction. For the direct solution of the combined modelof an asymmetric phantom, this error presented itself as alternating intensities of light anddark sectors strongest at the origin. This 'starburst' artifact propagates through the sectorsto the outer annuli of the reconstruction. This is shown for phantom (d) of Figure 5.1 in Fig.5.2. The reconstruction of (a) exhibits the artifact due to center miscalculations, while (b) is areconstruction with the center explicitly defined at the correct position.(a) (b)Figure 5.2: Center artifact in reconstruction(a) artifact (b) forced center5.2.1 Modified centered findingThe error in finding the center of the log arises when the centroid of the Radon Transform doesnot coincide exactly with the geometric center of the log cylinder. This error occurs when thecross-section contains a substantially non-symmetric density distribution. The non-symmetrycan be reduced (but not entirely eliminated) by penalizing those parts of the radon signal thatare larger in magnitude while increasing the relative size of the lower parts of the signal. A69simple approach is to take the square root of the Radon Transform, which in effect biases thecentroid away from the areas of high area density. In practice, it was found that results wereimproved even further by taking the cube root, subsequent gains were not readily observablefor higher order roots. Figure 5.3 shows an example of this effect. The initial Radon Transformsignal for the synthetic phantom (e) of figure 5.1 taken at θ = θ ′ as shown in part (a). Thelarger, high density 'knot' region causes the ordinate centroid to lie at 63.7, slightly less thangeometric center which is represented by the vertical dashed green line at pixel 64. Takingthe cubed root of the transform as shown in (c) effectively makes the predicted centroid lie atprecisely pixel 64. Symmetric views are unaffected by the process, and on the whole it servesto provide a more accurate center prediction.(a)(b) (c)Figure 5.3: Center finding example(a) Ground truth f with θ = θ ′ marked (b) R f (t, θ ′) (c) cubed root of R f (t, θ ′)705.2.2 Tikhonov regularizationThe modified center-finding algorithm is effective at biasing the center in the correct direction formost views, but it is not perfect. The reconstructions still feature faint traces of alternating light-dark sector elements about the center due to very small misalignment. In an effort to alleviatethese areas of high contrast, Tikhonov smooth regularization in the circumferential directionwas applied to the center voxels in the path length matrix before inversion. This techniqueworks to make the form of the intensity's curvature among adjacent sectors adjustable in orderto allow for varying degrees of smoothness in the reconstruction [20]. Smoothing the voxels'curvature allows for mild regularization in order to not lose the features of interest.The smoothness of the voxels in the central circle is analogous to the second derivativeof the density function in the circumferential direction. This second derivative may be rep-resented in matrix form for a numerically sequential arrangement of elements by the matrixin Figure 5.4. In the case of taking the density function in sequence about the central annulifor the arrangement in Figure 3.16, this matrix must be adjusted for the non-sequential num-bering in this direction. This effectively yields a very sparse tri-diagonal Nvoxel s × Nvoxel smatrix with 2's located at entries (p × Nannuli, p × Nannuli ), p = 0, 1, 2, 3... and -1's at(p × Nannuli, (p × Nannuli + 1)modNsector s ) and (p × Nannuli, (p × Nannuli + (Nsector s − 1))modNsector s ).-1 2 -1-1 2 -1-1 2 -1-1 2 -1c =Figure 5.4: Tikhonov regularization matrix for curvatureFollowing the method, this derivative matrix may be pre-multiplied by it's transpose to bedimensionally consistent with GTG and multiplied by a coefficient α in order to modulate its71effects on the inversion. This yields the altered system of normal equations contained in eq. 5.2.GTG ~ρ′ + αcT c = GT~ytotal (5.2)There are many schemes for optimization of the parameter α, to achieve adequate noisereduction while avoiding distortion of the solution by excess smoothing. For this case, a roughneighborhood of acceptable values was determined by starting with a very low coefficient (10−6)and raised by orders of magnitude while observing changes in the mean-2 norm relative error of~ρ′. When excessive blurring caused the reconstruction to degenerate, the coefficient was reduced.A coefficient of α ∼ 10−3 was ultimately found to be sufficient to introduce adequate smoothingwithout causing significant distortion. Figure 5.5 displays the synthetic ground truth phantom(e) of Fig. 5.1 with the initial artifact in (a) and corrected using the modified center findingand regularization in (b). The mean-2 norm error in this case was reduced from approximately11% to 4%. This example may be extrapolated to the multi-slice case, and shows the positiveeffect that regularization may have on the ART method of reconstruction.(a) (b) (c)Figure 5.5: Phantom Artifact Correction(a) artifact present (b) Corrected center (c) Regularization and corrected center5.3 Real data reconstruction5.3.1 Calibration of the detectorIn order for the inversion to give accurate density values for the reconstructed f (r, φ), the x-ray detector must be calibrated to correlate the CT data, − ln(II0), with the sampled RadonTransform,∑gi j ρ j . Looking back at equation 2.2, it is apparent that this amounts to findingthe mass attenuation coefficient β that corresponds the reconstructed log's wood. For this work,72the method of An performed at 100kVp, 5mA x-ray power was closely followed [19]. Low densityfiberboard, a wood-based material of known density and dimension, was used to perform thiscalibration. The fiberboard samples were square blocks measuring approximately 15cm per sideand 2cm thick, with an approximate density of 500kg/m3. A table was set up in front of thedetector screen so that the blocks could be centrally located in the most intense part of thecone beam. Starting with one block a sequence of images was taken, incrementing the numberof blocks by one between images. This was done from one to twenty-five blocks in sequentialorder. Figure 5.6 demonstrates the manner in which these blocks were stacked for calibration.Figure 5.6: Fiberboard calibration blocksA reference image was taken with no objects in view to provide values for I0, and thecorresponding mean value of pixels within the block stack were used to find I. In order forthe calibration to provide the correct curve for the relationship between the CT data and theRadon Transform value of the attenuated rays, the same power of x-rays must be used for thecalibration and the reconstruction data images. It is critical though that the detector not besaturated for calibration, and so for a given power of x-ray (determined by the power needed topenetrate the log being scanned) the EM gain of the camera was adjusted until saturation wasnot present. Figure 5.7 demonstrates the un-attenuated signal in (a) along with the attenuatedsignal for 5 blocks in (b) and 14 blocks in (c).73(a) (b) (c)Figure 5.7: Calibration data(a) unattenuated (b) 5 blocks (c) 14 blocksThe calibrations were performed at 140kVp, 5mA and 160kVp, 5mA, which were the twoenergy levels used for this work. It is the case that the x-ray intensity of the source usedfluctuates by as much as 5% over the course of several seconds, and so for each attenuation stepten images were taken and their values averaged. The corners of the image outside of the centralcone beam were then used to scale the data sets to be 1:1 with the reference image in intensity.It was found that past fifteen blocks at both energy levels the resolution of the detector limitedany further calibration. A polynomial curve was fit in this case to the first fifteen data pointsto yield a correspondence between CT data and the Radon sample.5.3.2 Real log data sample and reconstructionIn order further to assess the center finding scheme and solution regularization approach, asingle slice reconstruction of the data taken by An for the same log (the dried Douglas fir)was performed, using An's detector calibration. The data corresponded to a capture time of1s at 100 kVp and with a camera gain of 1. The Radon transform vector was taken as thepixel column corresponding to a four knot cluster as is shown in Figure 5.8 by a dashed line[19]. The data set consisted of 36 images (equal to Nsector s), each taken at a rotation of 10degrees from one another and the log axis of rotation remaining unchanged in the detectorview. The reconstructions are presented in Figure 5.9. It was found that the annular and sectorcases were consistent with the reconstructions of An, but the combined case presented muchless prominent artifacts and a more stable reconstruction, pointing to the advantages of thetechniques of Section 5.2.74Figure 5.8: Single Slice data sample from An75(a) (b)(c) (d)Figure 5.9: Single Slice reconstruction(a) annular reconstruction Nannuli = 18 (b) sector reconstruction Nsector s = 36 (c) combinedreconstruction (d) combined reconstruction from An76Chapter 6Multiple slice reconstruction6.1 Computational considerationsCareful management of the computational requirements of storing the basis and the robustnessof the solution becomes a necessity when the reconstruction extends into multiple slices. Thesingle slice method of storing the basis in its complete form as well as implementing a directsolver would prove to be unfeasible due to the much greater size of the mathematical solution.6.1.1 Storage of the basisAlthough it is the case that G is highly sparse, this does not necessarily ensure that GTG issimilarly sparse. This makes the implementing the Cholesky solver of Chapter 5 very computa-tionally burdensome and slow due to memory writing and rewriting on a standard system (8GBRAM) as the number of slices gets larger. Full storage of G is also impossible, as the size of thepath length matrix considering null values is much larger than GTG. For this work, Matlab'ssparse matrix format was used in order to store G in a more efficient manner, ensuring goodruntime performance and easy manipulation of the basis. In this scheme, each nonzero elementof G is stored together with its matrix coordinate values (i, j).6.1.2 Ill-posedness of the basisCalculation of the normal equations also increases the condition number of the basis, effectivelysquaring it [21]. It is therefore undesirable to use the method in the rotational and helicalreconstruction cases as the condition number increases very quickly with additional slices. It was77found that the normal equation basis became effectively singular when the modest resolution 36sector, 18 annuli reconstruction method of An was extended into multiple slice reconstruction fora full-field view dataset (45 slices of .71cm thickness each to fill the detector). This necessitatedutilization of more conventional iterative solvers that have been used for some time in ART.6.2 Iterative methodsTwo iterative methods were explored in this work for the solution of the inversion problem,Randomized Kaczmarz and the Simultaneous Algebraic Reconstruction Technique (SART).6.2.1 Randomized Kaczmarz's methodKaczmarz's method is a solution belonging to the Projection onto Convex Space (POCS) familyof solvers, and has an established history as the preeminent solver for conventional iterativeART CT. It is preferred for its generality and adaptability, as the basis need not be square orpositive definite, as is the case for other iterative methods (Gauss-Seidel, Conjugate Gradient,etc.) The algorithm works as follows: the solution of the system of equations G ~ρ′ = ~ytotal , inwhich again A is a matrix of dimension m × n, may be thought of as the intersection of quantitym hyperlanes pi = ri · ~ρ in n-dimensional space, in which ri is the ith row of G. A starting 'guess'for the solution vector~ρ′0 |0 is made in Rn. The first row of the basis matrix, r1 is taken into thealgorithm and an orthogonal projection of the solution vector is made onto p1 = r1~ρ′0 |0. Thisprojection process is then repeated for the entirety of the m constraints, at which point it isconsidered that the algorithm has gone through one iteration and the process is repeated froma new starting point,~ρ′0 |1. Thus the first subscript therefore corresponds to the 'sub-iterative'step of a projection and the second subscript to the outer algorithm iterative step. It may beshown that the iterative step, the updated solution vector ~ρ′i |k+1 for the (k +1)th iteration, maybe given by~ρ′i |k+1 = ~ρi |k +yi − pi‖ri ‖22ri (6.1)Where i = 1, 2, ...m denotes the sub-iterative projections. The method may be graphicallyrepresented in the case of n = 2 as the intersection of multiple lines in a plane (Figure 6.1). Thehighly overdetermined system of equations is represented by quantity m lines {l1, l2, ...lm } thatintersect at the solution point ~ρ. As the algorithm operates, the initial guess projects along thepath shown by the dotted lines until it coincides with ~ρ.78ρρ'0|0l1l2l1l2l3ρ'1|0ρ'2|0ρ'3|0 = ρ'0|1l3xyxy(a) (b)Figure 6.1: Kaczmarz algorithm(a) Starting guess~ρ′0 |0 and solution ~ρ (b) First iterationIn the case of this research, in which the basis is ill-conditioned and overdetermined, manysolutions may be equivalently correct according to the standard algorithm. This can be visu-alized in the plane as the case in which multiple constraints of the system do not coincide ata single solution point but rather have a zone of intersection. Intuitively, it seems reasonablethat in this case the best choice is the geometric center of such a set of points, and convergencetowards this best choice is made with the use of a relaxation parameter, λ. Equation 6.2 showsthe relaxation parameter incorporated into equation 6.1. Its effect is graphically illustrated inFigure 6.2. The solution vector in this case will make its way towards ~ρ, the best solution forthe system.~ρ′i |k+1 = ~ρi |k + λyi − pi‖ri ‖22ri (6.2)79xy'i|j'i|j+1'i|j'i|j+1(a) (b)Figure 6.2: Overdetermined Kaczmarz(a) Standard algorithm (b) With relaxation parameterVershynin and Strohmer have improved upon the standard Kaczmarz algorithm by using aversion that selects the next hyperplane for projection based on a weighting probability distribu-tion over all rows proportional to the Euclidean norm squared of each row, that is ‖ri ‖22. It hasbeen shown to converge at an exponential rate to the equivalent least squares solution indepen-dent of the number of constraining equations based on the work by Strohmer and Vershynin andconfirmed by Needell [22, 23]. They conclude that the weighted selection of projections has thebenefit of convergence to the true solution in fewer iterations than possible with conventionalKaczmarz.6.2.2 Simultaneous algebraic reconstruction techniqueAnother variant of the standard Kaczmarz algorithm is the Simultaneous Algebraic Reconstruc-tion Technique or SART. In this implementation introduced by Andersen and Kak, the problemof an inconsistent overdetermined system is accommodated by taking into consideration thepotential projections onto all hyperplanes of a given view of the system, and then taking theaverage of these projection vectors as the sub-iterate step [13, 24]. Looking back at the basisconstruction of Section 3.3, this means that the rows of G are iterated down one at a time andthe projected solution vectors ~ρ′i |k are stored. When the sub-iterate index meets the condition80i = n(Pang × Pcol), n = 1, 2, ...Nviews the step of equation 6.3 is taken to update ~ρ′k (shown herewith indexing corresponding to the first view). Once all views are considered, one iteration istaken to have been performed, that is when ~ρ′k+Nviews (Pang×Pcol ) is calculated.~ρ′k+1 = ~ρ′k +∑Pang×Pcoli=1(~ρ′i |k − ~ρ′k)Pang × Pcol (6.3)By performing the updates in this manner a smoothing effect is obtained that greatly cutsdown on high frequency noise in the reconstruction. This smoothing is not enforced regulariza-tion, but rather an integral aspect of the algorithm. It does lead to slower feature distinctionthan the randomized Kaczmarz algorithm, however.Hansen et al. has shown that the algorithm for this scheme may be vectorized and im-plemented very simply, eliminating the need for looping [25]. It is presented here withoutderivation:~ρ′k+1 = ~ρk + D−1r GTD−1c (~ytotal − G~ρk ) (6.4)Here, Dr = diag(‖gi ‖1) and Dc = diag (g j1) are diagonal square matrices of dimensionm × m and n × n made up of the row and columnar sums of the coefficients in G respectively[25].6.2.3 Selective regularizationThe effectiveness of regularization on solution stability was demonstrated in the case of singleslice reconstruction, and so an effort was made to incorporate a similar scheme into the multipleslice and helical reconstruction algorithms. It is stated by Herman that for ART reconstructionsin which the form of the image function is made up of large distinct regions of different densities,good results may be accomplished through the implementation of so-called selective smoothing.In this process, which is a 3-D expansion of the planar case presented by Herman, each voxel isexamined upon iteration of the solver and compared to its nearest neighbors. The primary voxelv1 is surrounded in three dimensions by six voxels containing a common side v2 − v5, v10, v11and twelve voxels containing a common edge v6 − v9, v12 − v19. This geometry and numberingis shown in Figure 6.3, with the primary voxel in yellow, its shared face voxels in blue, and itscommon-edge voxels in red.81v1v2v3v4v5v6v7v8v9v10 v11v12v13v14v15v16v17v18v19Figure 6.3: Selective regularization voxel indexingThis voxel indexing is shown for a generic interior voxel that is not on the perimeter orcenter of the reconstructed log; it is surrounded on all sides by other voxels. In those caseswhere a voxel is not in contact with other voxels, a careful record must be kept to ensure thatextra voxels are not included in the regularization.Establishing a record of each voxel and its neighbors, a threshold value t is defined for thereconstruction which represents the difference in relative voxel values above which smoothingwill take place. This threshold function is 'activated' via the relationf i =1, if |vi − vt | ≤ t0, otherwise(6.5)The activated threshold function may then be inserted into equation 6.6 to provide a regu-larized value for v′1. The values w1, w2 and w3 are the smoothing weights of the regularizationthat affect the directionality of the regularization and the degree to which small incongruitiesare blended. The first weight corresponds to the voxel under consideration, the second to thecommon-face voxels, and the third to the common-edge voxels.82v′1 =w1v1 + w2∑5,10,11i=2 f ivi + w3∑9,12−19i=6 f iviw1 + w2∑5,10,11i=2 f i + w3∑9,12−19i=6 f i(6.6)Herman notes that determining the relative values of the smoothing weights that yield thebest results is a heuristic procedure and that it should be used carefully so as not to distortimportant features of the true solution. The threshold parameter enables smoothing to takeplace in areas of small differences, assuring that large important features are not lost in theprocess.6.2.4 Algorithm comparisonConvergence optimization of projection solvers and the tuning of ending conditions is an inten-sive process, and has oftentimes been performed heuristically for specific applications [13, 26, 27].For this work, the random Kaczmarz solver and vectorized SART solver developed by ChristianHansen et al. in the AIR Tools Matlab package were combined with a script implementing theregularization scheme of Section 6.2.3 to reconstruct the multi-slice rotational and spiral cases.Execution was carried out in Matlab 2013b on an Intel Core i7 2.8GHz CPU with 16 GB RAM.A data set of 'perfect data' was constructed by constructing a synthetic density vector of a100 slice, 10 annuli, 48 sector geometry of a desired 'log-like' distribution ~ρsyn interacting witha 92 Pang by 42 Pcol pixel rays to get G, and then taking ~ysyn=G~ρsyn. This data vector ~ysynis perfect in the sense that all the data are perfectly centered and scaled exactly so that onlyincongruities due to the algorithms used will be detected. Per recommendation by Herman, theiterated relaxation parameter in the random Kaczmarz was set to the under-relaxed value of.05 [28]. Herman heuristically found that good results for a noisy reconstruction was obtainedusing smoothing weights of w1 = 9, w2 = 4 and w3 = 2, and so these values were used in thescript attached to both solvers. It was found that the solvers converged on the ground truthsolution with negligible error for both solvers with the number of iterations for the randomized,under-relaxed Kaczmarz at two and SART at twenty iterations.Due to standard Matlab m-code being a non-compiled, looping is extremely inefficient in theenvironment. This means that even though fewer iterations were required for convergence fromthe randomized Kaczmarz solver, execution of it took dramatically longer than implementationof the vectorized expression in eq. 6.4. As such, almost all reconstructions carried out inthe analysis of the multi-slice and helical reconstructions to be presented were carried out withSART as reconstruction quality of both was seen to be equivalent. If implemented in a compiled83script, it may well be the case that randomized Kaczmarz would prove to be faster. Optimizationand tuning of the solvers described, including development of smart stopping criterion can becounted amongst further work that must be performed as the system matures.6.3 Rotational multiple slice6.3.1 Phantom reconstructionTests were constructed on the physical phantom tube of Section 4.4.1 in order to assess theeffectiveness of the detector and reconstruction scheme. It was found that the dense grain-filledinner portion of the phantom necessitated an x-ray power of 140kVp and 5mA tube current. AnEM camera gain of 30, 100ms exposure with a 200ms pause on capture was found to give theclearest picture. The reconstruction was set have 48 sectors, giving a sector angle of 7.5 degrees(half that of the average knot feature dimension), and 10 annuli starting from an innermostannulus equal to 16% of the outer dimension of the log. The normalized diameter of thereconstruction was set to have a value of .333m. The values of DL and of DS were measured tobe 1.625m and 2.125m respectively after the source was aligned.As with the single slice reconstruction, the number of views in the reconstruction was set toequal the number of reconstructed sectors. Although in the case of purely rotational reconstruc-tion the value for the slice thickness t is decoupled from the number of sectors and pulley size,.71cm was used to correspond to the slice thickness that would be obtained from the spiral casefor the reconstruction geometry. Several data captures are shown in Figure 6.4. The centralknot cluster in the image is the Slice B of Figure 4.14.84Figure 6.4: Rotational multi-slice data captureLooking at the data images, it is clear that at the lateral edges the intensity of the conebeam significantly decreases. It was found that reconstructions were affected when the multipleslice reconstruction contained data from these areas, and so the number of slices Nsl ices waslimited to be 42 in order to capture only the best quality region of the cone beam (roughly themiddle 2/3rd of the detector). The number of angular pixels Pang in this case was heuristicallychosen to be 100 and pixel columns Pcol to be 74 in order to provide an adequate sampling ofall voxels in the inversion.Using the altered SART-smoothed algorithm described previously yielded a reconstructionthat could then be evaluated based on the dimensional positions of reconstructed features vs.the known phantom specifications. The reconstructed voxel values could also be observed vsthe known density of the physical phantom to assess the calibration and scaling of data. Figure6.5 presents two cross sections of the multi-slice reconstructed phantom.85(a) (b)Figure 6.5: Centered multi-slice reconstruction slices(a) clear cross sectional region (b) knot cluster BAs can be seen, the algorithm effectively reconstructs the form of the physical phantom,showing the dense grain inner tube as well as the outer packed sawdust annulus. The grayscaleindicates that the density values are reasonable based on the measured densities of the insertedfeatures and is expectedly uniform. Unfortunately, the system is unable to differentiate clearlythe inserted MDF knot features within the grain filled cylinder (the measured density differencebetween these two materials was 80kg/m3). The innermost annulus is less dense than its neigh-bors due to the structural dried timber core of the phantom. The transition annulus betweenthe grain and sawdust as well as the outermost annulus display a higher density than might beexpected, but this is simply due to the denser cardboard tube being present in these locations.Still - the reconstruction is not perfect, and interestingly there are some artifacts present aroundthe crack feature as well as the knots.A script was written to measure the angular locations of the five knots of Slice B in thereconstructed log as they were presented in Table 4.1. Taking the median direction of each knotfrom its form in the outer sawdust annulus, the locations of the knots could be located relativeto the coordinates of Fig 4.14. Detection of the knots was performed manually in this case, thedirection vectors for the knots were determined by selecting their edges in the axes of Figure6.6 (the dashed red lines) and then averaging each pair to define the knot directions (the solidred lines). Table 6.1 presents the knot locations of the reconstruction compared to measuredvalues. They show to be accurate for several knots compared to the measured phantom, with two86showing discrepancies from the measured phantom on the order of one sector angle, 2pi/Nsector s .Given the relatively imprecise assembly of the phantom (the tamping was performed so as to notdisturb the knot positions, but some displacement is likely) and the crudeness of measurementof their positions, the results are encouraging.Figure 6.6: Centered multi-slice knot directionsKnot Measured Reconstructed1 0 72 78 783 198 1924 229 2305 300 300Table 6.1: Slice B knot locations - measured and reconstructedAll values in degrees from x axis.The purely rotational multi-slice algorithm was also tested in the case of rigid body displace-ment of the log within the cone beam. This was accomplished by attaching the log-off center onthe carriage system, using two counterweights shown in Figure 6.7 to ensure that no backlashof the timing belt or stalling of the motors would occur. The geometry of the system-detectorarrangement allowed for a maximum eccentricity of one inch of the log axis from the shaft. The87series of images in Figure 6.8 demonstrates the image set; the phantom can clearly be seen todrift up and down over the series of images as it rotates.Figure 6.7: Off-center counterweights(a) (b)Figure 6.8: Rigid body motion multi-slice data(a) log at max displacement upwards (b) log rotated at max downward displacement88The cross sections presented previously are again shown in Figure 6.9. The reconstructionhas experienced very slight degradation in some areas due to the discrepancies introduced intothe small angle approximations of the cone-beam made in the normalization scheme. Still thislevel of distortion is not prohibitive and the knots are still visible, as well as their locationsbeing comparable to the reconstruction for the centered multi-slice reconstruction. The knotsof Slice B in Figure 6.9(b) again fall within one sector angle of their measured values.(a) (b)Slice B EccentricSlice B CenteredPhantom CenteredPhantom Eccentric(c) (d)Figure 6.9: Off-center multi-slice reconstruction slices(a) clear cross sectional region eccentric (b) knot cluster B eccentric (c) clear centered (d)B centeredThe trials prove the efficacy of the multi-slice rotational algorithm for providing an accurate89density reconstruction of the phantom for the combined voxel pattern. As the detector is notwide enough to accommodate multiple knot clusters, lateral measurement from the phantom isnot possible. It should be stated that minimum optimization was performed on the values ofPang and Pcol , it was simply taken sufficient that these values provided in excess of two pixel-rays for even the smallest central voxels, and that each slice received at least two pixel columnsof intersection. This ensured a sufficiently overdetermined system to give the best results giventhe relatively noisy data (as compared to commercial systems).6.4 Helical scanningUsing the methods of Section 3.3.3, the helical algorithm may be readily implemented for re-constructing the physical phantom.6.4.1 Effect of limited viewsA key difference between purely rotational reconstruction in which Nviews = Nsector s and helicalreconstruction is that the detector is not wide enough for this condition to be met (i.e. t×Nsector sis greater than the greatest possible Wz). This means that each voxel will not be able to receivethe same level of overdetermination that was possible in the purely rotational case. UsingTurbell's criterion of γ = pi + 2Ψ being the minimum angle of rotation needed for complete dataacquisition of a slice, it is found for the geometry under consideration that this amounts to ~200degrees of rotation, or in the case of 48 sectors Nsl ices_view = 26. With the pulley setting theslice thickness to be .71cm, this means that the log will rotate 26 × (360/48) = 195 degrees in adistance of 26 × .0071 meters, giving a helical wavelength of .34m.6.4.2 Phantom reconstructionThe data set now consists of many more images, for the .71cm slice thickness set by the pulleydiameter and sector geometry means that approximately 84/.71 = 118 slices will be neededto reconstruct the entirety of the phantom. With Nsl ices_view = 26, this means that in total26 + 118 = 144 images are required. Figure 6.10 demonstrates the nature of the data set byshowing the far ends of the log as it enters and exits the active cone beam (marked by thevertical red lines). By definition only the first register is 'filled' in (a) and only the last in(c). It is possible for Nsl ices_view to be greater without data entering the dim region of thecone beam (giving a greater degree of coverage), but preliminary trials showed this additional90overdetermination to not help the solution's convergence. As such, all helical tests in this thesisuse the minimum set by Turbell's rule to ensure the active detector width remains in the bestquality region.Figure 6.10: Centered spiral reconstruction dataPhantom entering and exiting the cone beamThe SART-smoothed reconstruction gives good performance on the inversion, performing itin ~8s. The three knot cluster slices obtained from the inversion are shown in Figure 6.11, againshowing clearly demarcated and uniform areas of materials inside the phantom, along with themean knot directions. The same script was run as before in order to determine the positions ofthe knots within the phantom, the values for these are given in Table 6.2. Again, the resultsseem reasonable given the likely error in measuring the phantom, and are close to the measuredvalues. The adjacent slices in which the respective knot clusters A, B, and C appeared could beaveraged in order to give a slice value for their position in the phantom. The number of slicesfrom the end cap of the phantom to these could then be counted in order to assess the lateralpositioning of these knot slices compared to the measured values. It was found that their meanpositions occurred in slices 17, 55, and 92 respectively. Taking into account the slice thickness,this corresponds to locations for knot clusters A, B, and C at 12, 39, and 66cm respectively.Looking back at the measured values of Figure 4.14, this confirms that the algorithm is puttingreconstructed objects back into their proper locations.91Slice A Slice B Slice CFigure 6.11: Centered spiral reconstruction of physical phantomKnot Slice A Slice B Slice C1 93 6 2032 210 70 3003 243 191 3494 341 229 //5 // 300 //Table 6.2: Centered spiral reconstructed phantom knotsAll values in degrees from x axis.The same test was conducted on a set of spiral data in which the log was displaced off-centeras was done in the rotational case. Again the maximum possible value of 1 eccentricity wasused. The forms of the reconstructed slices for the off center reconstruction are shown in Figure6.12. It is the case that the center finding and normalization algorithm again prove effectiveat reconstructing the knot cluster positions, and analysis confirmed that the knot locationvalues to be as accurate as was the case for the centered spiral reconstruction, as well as theirlateral positions. It is clear that within the allowable eccentric limits of the setup, the smallangle approximation for the path lengths holds and the reconstruction of the phantom is veryminimally impacted.92Slice B Slice CSlice AFigure 6.12: Off-center spiral reconstruction of physical phantomThe series of slice images for the reconstructions were converted into NIfTI file-format medi-cal images, paying respect to their scaling and spacing to preserve real-dimension reconstructiongeometry [29]. The centered and eccentric spiral reconstructed phantoms were visualized in theopen-source segmentation software ITK-Snap, allowing for complete visualization parallel to theyz and xz planes of Figure 4.14 [30]. This showed the uniformity of the inner grain core alongwith outer sawdust annulus. Screen captures from the software for the central are presented inFig 6.13 for the centered spiral reconstruction. It is clear from these images that the boundarybetween the dense center cylinder and the outer annulus is the noisiest location of the recon-struction. Using ITK-Snap, it is possible to produce a manual segmentation of the reconstructedphantom (the ambiguity of the grain cylinder prohibits use of automatic segmentation). A ren-dering of this segmentation is presented in Figure 6.14. With the algorithm appearing to supplyaccurate spatial reconstruction, real log reconstructions could be undertaken with confidence.93(a) (b)Figure 6.13: Planar cross sections of centered spiral reconstruction(a) yz plane (b) xz plane of Fig. 4.14Figure 6.14: Segmentation rendering of physical phantomKnots are shown in green, the crack object in pink. Grain and sawdust in translucent brown.946.4.3 Real log reconstruction of dry logThe dried Douglas fir log reconstructed in the single slice case by An was the first naturallog reconstructed. The specimen is a good case study for the algorithm because it containsmany knots, and is easily penetrated by the cone beam at 140kVp and 5mA due to it beingfully dried. There are several large knot clusters visible externally, as well as multiple smallerinternal knots. The knot clusters in this particular tree are also nearly planar, that is they havean almost normal branch angle to the longitudinal dimension of the log. Examples of these arevisible in the radiograph data sample of Figure 6.15.Figure 6.15: Knot clusters in dried douglas firExternal knot cluster in red ellipse, internal knot cluster in yellow ellipseIn order to capture these small internal knots, the spiral script was run with a 72 sector, 15annular reconstruction with Nsl ices_view = 40. The geometry now dictates that the slices be.48cm thick, making 180 slices necessary to capture the entire log. This high resolution makesfor the reconstruction to be ~3x as large as that of the physical phantom, and as expected the 20iteration SART-smoothed script finishes in more time, approximately 30s. Several reconstructedknot clusters are presented in Figure 6.16. The scaling here is the same as that for the phantomreconstruction, and the color-bars indicate the density being within the range expected for dry95softwood (~500kg/m3).(a) (b) (c)Figure 6.16: Combined pattern dry log slices - centered spiral(a) five large knot cluster (b) internal knot (c) four large knot clusterThe same reconstruction may be undertaken on the dry log using a lower resolution for thesake of speed. If it is the case that only large knots are of interest and that the heartwood-sapwood boundary is not critical, the same data set may be broken down into a sector-only casein which each slice is comprised of 72 sectors. This reduces the system of unknowns by a factorof 15, and as expected the convergence in this case takes much less time, <8s. In reality, abalance between speed and quality of detail is always the tradeoff of any CT system. Note from6.17 that while larger externally visible knots are clearly recorded, smaller internal knots are lostin the clear wood (slice (b)). The optimization of this balance of resolution vs computationaltime must be varied according to the geometry of the reconstructed tree.(a) (b) (c)Figure 6.17: Sector pattern dry log slices- centered spiral(a) five large knot cluster (b) lost internal knot (c) four large knot cluster96The Dried douglas fir log may be cut along its longitudinal major planes in the same way aswas the phantom in order to examine the inter-slice boundaries for artifacts. The reconstructionspresented in Figure 6.18 illustrate clearly two longitudinal orthogonal planes along the lengthof the sample. The heartwood sapwood boundary as well as elongated spike knots and captiveface knots are visible.(a) (b)Figure 6.18: Dry log major planes(a) xz plane (b) yz planeIn this case, the discrete areas of high density in the knots make it is possible to use theautomatic level-set segmentation algorithm of ITK-Snap in order to extract the knots from therest of the log without user input. This was performed, and the reconstructed knot structureis presented in Figure 6.19. The clear wood (both heartwood and sapwood) is presented intranslucent yellow, the knots again in green, and the large twisting crack in pink. The planarnature of the knots is clearly evident.97Figure 6.19: Dry Log segmentationClear wood is yellow, knots are green. The radial crack is in pink.6.4.4 Real log reconstruction of green logA more challenging task is reconstruction of the green log, as it is of a much higher density.This makes x-ray penetration more difficult. It was found that the x-ray source need be run atits highest power for this, 160kVp and 5mA. A sample of the data is presented in Figure 6.20,shown to the same scaling as that of the dry log data (Fig 6.15).98(a) (b) (c)Figure 6.20: Wet log centered spiral data(a) log entering cone beam (b) log in view - non circular form evident (c) log exitingThe green Douglas fir sample differs from the dried one in that it has large distinct knotsthat extend to the surface. They are also less grouped than in the dry sample. Finally, theknots in this log have a large branch angle, that is they have grown at an inclination to thelongitudinal (z) axis of the log. Also challenging with this log is an area that is not truly circular,but exhibits a locally flattened recess in which one of the knots is growing (visible in (b) of Fig6.20). The system was put to the test by reconstructing the wet log in 60 sectors and 15 annuli.The geometry necessitated a slice thickness of .57cm and therefore ~140 slices for a completereconstruction. Reconstruction takes ~15s. To handle the non circular region of the log, theestimated log radius Rpix is scaled by 96% in order to truncate the outer most irregularities. Themajor planes of this fresh log are shown below in the captures from ITK-Snap. The intensityscale is the same as was used previously, which makes it evident that this log is more dense (itcontains a mean density of 750kg/m3, with sapwood averaging 900kg/m3).99(a) (b)Figure 6.21: Wet log major planes(a) xz plane (b) yz planeNote as well from Figure 6.21 the large branch angle, this makes showing the knots by sliceless useful, as only a small part of each knot is observable in a planar transverse slice. As mostconiferous trees grow with a consistent branch angle according to the geographical latitudeat which the tree was harvested, this presents some opportunities that will be discussed. Atransverse slice displaying a partial knot intersection is given in Figure 6.22. A rendering of acomplete segmentation is given in Figure 6.23. The ambiguity of the knots from the sapwood inthe outer rings of the tree make a distinct segmentation of this region via automatic methodsimpossible, and so the knots are truncated in this region in order to clearly show the sapwoodboundary.100Figure 6.22: Wet log cross sectional sliceBranch angle causes knot to pass through transvese slice.Figure 6.23: Fresh log segmentation101Chapter 7Segmentation of log internal features7.1 Reconstruction of lumber surfacesThe inherent value of CT reconstructions of logs is accurate visualization and decision makingabout the processing that will optimize yield of knot and defect free lumber. The results for boththe green and dried logs may be manipulated to demonstrate this. One of the most straight-forward processing techniques in an industrial softwood mill is the plane-sawing of dimensionallumber and paneling. This is when a series of parallel cuts are made in a given direction at aninterval to produce boards of varying widths. Take for example the production of siding of 1thickness. Accounting for planing (i.e. adding half an inch to each thickness dimension yieldinga rough cut of 1.5) the dried and fresh logs yield seven boards each. They consist of two 1×6's,two 1×8's, and three 1×10's. A diagram detailing the rough cut lumber's arrangement withinthe log is presented in Figure 7.1. The corresponding boards produced from the dry log arepresented in Figure 7.2 and those that come from the fresh log are shown in Figure 7.3. Theimages show clearly sapwood-heartwood transitions as well as the locations of spike and faceknots. Defects are clearly shown in the case of the dry log.1021234567Figure 7.1: Plane sawn boards1 2 3 4 5 6 7Figure 7.2: Dry log plane sawn boards1031 2 3 4 5 6 7Figure 7.3: Fresh log plane sawn boardsVisualizations of other lumber products are possible as well. In the case of quarter-sawnboards, boards are cut at alternating right angles to one another in each quadrant of the log.This makes it so the grain runs transversely to the thickness dimension of the board and reduceswarping as the piece dries out. A diagram of five boards of one and a half inch rough cut thicknessis presented below for the reconstructed real logs. The corresponding pieces from the dry logand green log are shown in Figs. 7.5 and 7.6. As expected, all boards present the transitionfrom heartwood to sapwood characteristic of lumber produced in this manner.12345Figure 7.4: Quarter sawn boards1041 2 3 4 5Figure 7.5: Dry log quarter sawn boards1 2 3 4 5Figure 7.6: Fresh log quarter sawn boards1057.2 Feature identification via data plotsIn the implementation of a log scanning package into an industrial softwood mill, there will notbe time for a human operator to make decisions on what log orientation and combinations ofcuts will yield the highest profit. Thus scanning and accurate reconstruction is only the firstpart of the problem. The various voxel values in the solution vector ~ρ must be assembled andthe features contained recorded in space in order for machine decision making. Segmentationof a large number of disparate voxels based on thresholding alone is problematic. It is clearfrom Fig 7.3 that individual voxels contained in the sapwood and knotty regions of the logmay be of a high enough intensity so that knot features identified by simply their magnitudewill be extremely distorted. This points to the necessity of more discerning smart algorithmsto carry out the actual segmentation of the knots. Conventional segmentation algorithms arepowerful tools for delineating the extents of features, but have traditionally required an operatorto specify a smaller volume of interest and 'seed' the body in a step known as registration. Itis known that for medical applications, no machine has thus far been able to match humans inthe process of registration, in large part due to the complexity of human organ shapes and theirabnormalities [31].Fortunately the a priori knowledge of the way in which knots are present in logs combinedwith the expedient indexing of the voxel pattern serves to make knots readily identifiable forsegmentation seeding. Additionally, the data truncation and normalization ensures a volume ofinterest with no empty space. If the completed reconstruction vector is iterated through slice byslice and the annuli examined starting at the center and working outwards, it can be shown thatthe knots are present as local maxima in the density plots. This may be shown by looking at thedry log reconstruction slices of Figures 7.7 and 7.8, and then their corresponding density plots,presented here as surface plots over consecutive annuli within the slice. This is convenientlyaccomplished by simply taking the indexing of the reconstructed voxel density values as theyare outputed by the algorithm and plotting them side by side by sector in ascending order. Across section of the five knot cluster surface is presented in Figure 7.9, corresponding to thefourth annulus of voxels in that slice.106(a) (b)Figure 7.7: Surface corresponding to internal knot and crack(a) (b)Figure 7.8: Surface corresponding to five knot cluster and crack107Figure 7.9: Fourth annulus voxel ring of five knot clusterNote the five distinct peaks corresponding to knot locations.The clear local maxima that make up the knots may be recorded and their locations reg-istered as the coordinates from which segmentations should be seeded. Delineation may thenoccur based in the reconstruction space from which convergence is sure to be accelerated due tonumerically close neighbor voxels. The knot features may then be projected onto a single planeand a decision matrix implemented in order to orient the cut pattern so as to maximize valueyield based off of predetermined factors. This implementation is of course a project unto itself.108Chapter 8Conclusion8.1 Limitations of the method8.1.1 Assumptions madeThe assumption of small displacements having a negligible effect on the series solution basis sothat the log may simply be centered and scaled will clearly deteriorate for very large movements.An interesting case study could be made utilizing a smaller diameter physical phantom log sothat the limits of this assumption could be evaluated.The cone beam in this implementation was kept to a distinct central portion of the detector.This had two effects, one of keeping the CT data within the area of highest signal, and alsominimizing the error due to small misalignment between the source and log and log and detector.The perspective alignment method of Figure 4.9 is sufficient for using only this smaller centralregion. Undoubtedly, the increased noise in reconstruction encountered when using the fulldetector width is in part to increasing perspective error.8.1.2 Difficulty of integration into millThe greatest challenge in this scheme's implementation into a sawmill are the constraints ofspeed on the acquisition of data, reconstruction, and the preprocessing-registration-segmentationprocess. Given the constant operation nature of modern sawmills, the challenges of hardwareto meet the constant service requirements of the CT scanner are formidable. Mechanisms mustbe designed in order to obtain data captures rotated around half of the log, and these in turnmust be capable of long service life with minimal maintenance. Detector width, reconstruction109geometry, and physical detector selected must be coordinated to put a premium on speed toallow time for segmentation and log manipulation.8.2 Future work and considerationsA major breakthrough in the method will be the development of non-planar slices that reflectthe geometry of the branch angles in a given specimen. The essence of this is the extension of thefeature tailored voxel concept into another dimension to provide the commensurate advantagesas are available in the planar reconstruction case. The branch angle of a tree is determinedby a great number of factors, including the environmental circumstances of the stand the treewas harvested from. Much research has been put into the environmental growth factors as theyrelate to branch angle, and these may be combined with sensing in order to optimize the conicalplanar geometry and spacing. Consider the image of Figure 8.1. It is clear that with alteredslice geometries, fewer slices and thus voxels in total will be needed to capture an equivalentnumber of features in the green log.Figure 8.1: Green log reconstruction with superimposed conical slicesFor the hardware required for implementation, promising work has been done in the area ofsecurity CT scanning which has usage demands similar to those of sawmills. It is possible that110Electron Beam CT may be an attractive alternative to more conventional and fragile contactring rotating gantry systems. Another possibility is a fixed cone beam orientation as used inthis work's experimental setup. This would require development of a system that may rapidlytranslate and rotate a log without any line down-time.Finally, development of machine decision making schemes must be implemented in orderto take the density vector output from the presented algorithm and make cut decisions. Thismust take into account not only the features of interest present in the log but also the desiredproduct to be made and its quantity. The output must have rapid effect on log manipulationand orientation they are sent to different areas in the plant to be processed.8.3 ConclusionRaw material in the form of logs is the principal operating cost in modern sawmills. Researchinto the use of penetrative sensing in sawmills to maximize the profit from individual logs hasbeen driven by this. Computed tomography has been explored as a possibility in realizinga system that can fully map out features of interest within logs and optimize orientations forcutting lumber. Implementations thus far have been costly, delicate, and overly precise for use insawmills. The robust Algebraic Reconstruction Technique (ART) methods of CT reconstructionhave gained some traction in industrial use, but are slower due to large numbers of genericCartesian voxels.In order to speed up ART, a customized polar voxel pattern geometry in which the shapeand arrangement of the voxel structure reflects the features of interest within logs was developed.This greatly reduces the necessary resolution required for determination of features such as knots,the heartwood/sapwood boundary, and physical defects such as cracks. Careful manipulationof the log by requiring one sector rotation and one slice thickness translation in between datacaptures allows for efficient construction of the series solution basis by reusing the path lengthsof just one view, saving computation time and allowing for rapid storage. A data scaling andnormalization scheme allowed for precomputation of the basis before runtime and had the effectof making the inversion robust to the effects of log movement circumferentially and radiallywithin the cone beam. Scaling also made it so no voxels within the reconstruction space wouldbe wasted on empty space.Synthetic data arranged so as to mimic real collected data pointed towards alteration of thecenter finding scheme as well as the need for solution filtering. Two iterative solving techniques,Randomized Kaczmarz and Simultaneous Algebraic Reconstruction Technique (SART), were111combined with a selective regularization scheme to provide fast reconstruction. Radiographicdata taken with an experimental setup was used in order to demonstrate the effectiveness of thebasis and solution scheme on a physical log phantom, a dried log, and a green log. The physicalphantom confirmed that no large scale distortion was introduced into the reconstructions. Thereal log reconstructions and subsequent segmentations show the promise of the setup for lumberquality sorting. The numerical indexing output of the reconstructions is shown to be convenientfor quick registration and segmentation of knot features.112Bibliography[1] Carl Gustav Lundahl. Optimized Processes in Sawmills Optimized Processes in Sawmills.PhD thesis, Lulea University of Technology, 2007.[2] Alfred Rinnhofer, Alexander Petutschnigg, and Jean-Philippe Andreu. Internal log scan-ning for optimizing breakdown. Computers and Electronics in Agriculture, 41(1-3):721,December 2003.[3] Mauricio Varas, Sergio Maturana, Rodrigo Pascual, Ignacio Vargas, and Jorge Vera.Scheduling production for a sawmill: A robust optimization approach. International Jour-nal of Production Economics, 150:3751, April 2014.[4] Suchendra M. Bhandarkar, Xingzhi Luo, Richard Daniels, and E. William Tollner. Anovel feature-based tracking approach to the detection, localization, and 3-D reconstructionof internal defects in hardwood logs using computer tomography. Pattern Analysis andApplications, 9(2-3):155175, 2006.[5] Maria Magnusson Seger and Per Erik Danielsson. Scanning of logs with linear cone-beamtomography. Computers and Electronics in Agriculture, 41(1-3):4562, 2003.[6] Wenshu Lin and Jingxin Wang. An integrated 3D log processing optimization system forhardwood sawmills in central Appalachia, USA. Computers and Electronics in Agriculture,82:6174, 2012.[7] Stefan M. Stängle, Franka Brüchert, Antti Heikkila, Timo Usenius, Arto Usenius, andUdo H. Sauter. Potentially increased sawmill yield from hardwoods using X-ray computedtomography for knot detection. Annals of Forest Science, in press:5765, 2014.[8] Yuntao An and Gary S. Schajer. Geometry-based CT scanner for measuring logs in sawmills.Computers and Electronics in Agriculture, 105:6673, 2014.113[9] Gary S Schajer, Peter Lister, Sencer Alkan, and Gabor Szathmary. Low Resolution CTScanning for Aspen Log Sorting. 2006.[10] F. Natterer. The Mathematics of Computerized Tomography. John Wiley & Sons Ltd,Chichester, 1986.[11] Charles L. Epstein. Introduction to the Mathematics of Medical Imaging. Pearson Educa-tion, Inc., Upper Saddle River, 2003.[12] R L Siddon. Fast calculation of the exact radiological path for a three-dimensional CTarray. Medical physics, 12(2):252255, 1985.[13] Avinash C. Kak. Algebraic Reconstruction Algorithms. Computerized Tomographic Imag-ing, pages 275296, 2001.[14] G T Herman and L B Meyer. Algebraic reconstruction techniques can be made computa-tionally efficient [positron emission tomography application]. IEEE transactions on medicalimaging, 12(3):600609, 1993.[15] Henrik Turbell. Cone-Beam Reconstruction Using Filtered Backprojection. Number 672.2001.[16] Gabor Szathmary. CT Cross Section of Fir, Image, 2013.[17] Christian Thibaudeau, Jean-Daniel Leroux, Rejean Fontaine, and Roger Lecomte. Fully3D iterative CT reconstruction using polar coordinates. Medical Physics, 40(11):111904,2013.[18] Li Jian, Li Litao, Cong Peng, Song Qi, and Wu Zhifang. Rotating polar-coordinate ARTapplied in industrial CT image reconstruction. NDT and E International, 40(4):333336,2007.[19] Yuntao An. Coarse-Resolution CT Scannin for Sawmill Logs Sorting and Grading. PhDthesis, University of British Columbia, 2013.[20] A.N. Tikhonov. Numerical Methods for the Solution of Ill-Posed Problems. Springer, 1995.[21] Cleve Moler. Numerical Computing with Matlab. The Mathworks Inc., Natick, MA, elec-tronic edition, 2013.114[22] Deanna Needell and Joel A. Tropp. Paved with good intentions: Analysis of a randomizedblock Kaczmarz method. Linear Algebra and Its Applications, 441(August):199221, 2014.[23] Thomas Strohmer and Roman Vershynin. A randomized Kaczmarz algorithm with expo-nential convergence. J. Fourier Anal. Appl., 15(1):262278, 2009.[24] AH Andersen and AC Kak. Simultaneous algebraic reconstruction technique (SART): asuperior implementation of the art algorithm. Ultrasonic imaging, 6(1):8194, 1984.[25] Per Christian Hansen and Maria Saxild-Hansen. AIR tools - A MATLAB package of alge-braic iterative reconstruction methods. Journal of Computational and Applied Mathematics,236(8):21672178, 2012.[26] Ming Jiang and Ge Wang. Convergence of the Simultaneous Algebraic ReconstructionTechnique (SART). IEEE Transactions on Image Processing, 12(8):957961, 2003.[27] Anders H. Andersen. Algebraic reconstruction in CT from limited views. IEEE Transactionson Medical Imaging, 8(1):5055, 1989.[28] Gabor T. Herman. Fundamentals of Computerized Tomography - Image Reconstructionfrom Projections. Springer-Verlag London Limited, New York, 2nd editio edition.[29] Jimmy Shen. Tools for NIfTI and ANALYZE image. Accessed7/04/15. URL http://www.mathworks.com/matlabcentral/fileexchange/8797-tools-for-nifti-and-analyze-image.[30] Paul A. Yushkevich, Joseph Piven, Heather Cody Hazlett, Rachel Gimpel Smith, SeanHo, James C. Gee, and Guido Gerig. User-guided 3D active contour segmentation ofanatomical structures: Significantly imprved efficiency and reliability. Neuroimage, 31(3):111628, 2006.[31] Jayaram K. Udupa. 3D Imaging: Principles and Approaches. In Jayaram K Udupa andGabor T. Herman, editors, 3D Imaging in Medicine, chapter 1, pages 3637. CRC PressLLC, Washington D.C., 2nd edition, 2000.115"@en ; edm:hasType "Thesis/Dissertation"@en ; vivo:dateIssued "2015-11"@en ; edm:isShownAt "10.14288/1.0216012"@en ; dcterms:language "eng"@en ; ns0:degreeDiscipline "Mechanical Engineering"@en ; edm:provider "Vancouver : University of British Columbia Library"@en ; dcterms:publisher "University of British Columbia"@en ; dcterms:rights "Attribution-NonCommercial-NoDerivs 2.5 Canada"@* ; ns0:rightsURI "http://creativecommons.org/licenses/by-nc-nd/2.5/ca/"@* ; ns0:scholarLevel "Graduate"@en ; dcterms:title "Robust CT scanning of logs with feature-tailored voxels"@en ; dcterms:type "Text"@en ; ns0:identifierURI "http://hdl.handle.net/2429/55136"@en .