SIAM News Blog
SIAM News

# A Reconstruction Algorithm is a Key Enabling Technology for a New Ultrafast CT Scanner

A sawmill’s goal is to produce as much valuable wood as possible from forest logs. Due to natural variability, every log is different, though most distinguishing characteristics are visible only after sawing. As a result, the process of deciding how to best cut a log is mainly driven by the log’s external appearance, which significantly reduces the ability to exploit the raw material’s true value. For this reason, knowing the real internal characteristics of a log before deciding its use has remained but a dream for sawmills.

Figure 1. Plane ΠFLT(s,x) through x,y(s), y(s1), and y(s2) contains vectors β(s,x), e(s,x) and is perpendicular to u(s,x). Image credit: Alexander Katsevich.
The possible creation of an image depicting the internal details of each log during production would allow for the simulation of various production strategies from which to choose. Typically, a large number of decisions must be made about each stem that enters a sawmill: how to buck it into shorter logs, which boards to produce, and how to best rotate and position the log in front of the saws. Advanced simulation of different solutions could improve these decisions. Many trials have demonstrated that the introduction of this technology could improve the productivity of a sawmill by between 5% and 20%. If we consider a large sawmill producing wood with a value of 70 million Euro/year, the uplift of 5% in production value would mean 3.5 million Euro/year.

In 2008, Microtec developed a first prototype of a computed tomography (CT) scanner for log imaging. It used a 180 kilovolt (kV), 10 milliamp (mA) X-ray source, rotated at 2.8 rev/sec, and implemented spiral tomography with a small cone angle of 0.5°. An approximate algorithm by Feldkamp, Davis, and Kress (also known as the FDK algorithm, see [3]) was implemented for the tomographic inversion, but the maximal scan speed was 5 m/min, much less than desired. Reaching the necessary speed with the same cone angle would have required an impossible gantry rotation speed. While wide cone angle scanning would allow one to reach the desired speed, there were no algorithms that would reconstruct efficiently and accurately from such data. On the one hand, approximate algorithms, like the FDK algorithm, are efficient, but produce artifacts at large cone angles. An accurate wide cone angle reconstruction with iterative algorithms, on the other hand, is prohibitively computationally expensive. However, the theoretically exact and efficient Katsevich inversion algorithm [6-8], the first enabling solution for implementing wide cone angle scanning developed in 2001, provides the key step to demonstrating the feasibility of a CT scanner for wood [4].

The Katsevich algorithm is theoretically exact and provides accurate reconstruction at any cone angle. It also has a very efficient filtered back projection (FBP) structure; it is easily parallelizable and can be implemented on a graphics processing unit (GPU). Every cone beam projection is used only once, an advantage implying that an image slice can be sent to the user as soon as all projection data for that slice has been collected and processed, without having to wait for collection of all projection data for a log. This property is very useful for a scanner, which operates in a continuous fashion. The original version of the algorithm, which was developed for an ideal constant pitch helix, was later generalized to include variable pitch helices and even more general curves [5, 9-10]. The resulting algorithm is still theoretically exact and maintains the FBP structure. Accounting for a variable pitch is important in saw mills because the belt speed can change quite frequently.

Figure 2. Image of the CT Log scanner. © Microtec.

Let us briefly describe the dynamic pitch algorithm of [9]. Let

$C := \{y \in \mathbb{R}^3 : y_1 = R\cos(s), y_2 = R\sin(s), y_3 = \psi(s), s \in \mathbb{R}\} \tag1$

denote a dynamic pitch helix. Here $$\psi(s)$$ is a three times continuously-differentiable function that maps the parameter $$s$$ to the $$y_3$$ (or $$z$$) position of the X-ray source on the helix. In particular, for a traditional, constant pitch helix, we have $$\psi(s) = s h/(2\pi)$$, where $$h$$ is the distance advanced by the source along the $$y_3$$ axis in one rotation. Here we assume that $$\psi$$ satisfies the condition $$\psi'(s)+\psi'''(s)>0.$$

Let $$U$$ be an open set strictly inside the helix:

$\bar{U} \subset \{ x \in \mathbb{R}^3 : x_1^2 + x_2^2 < R^2 \}\tag2$

The cone beam transform of a function $$f$$ supported on $$U$$ is $$D_f(y,\Theta) := \int_{0}^{\infty} f(y + \Theta t) \textrm {d}t$$, where $$\Theta \in S^2$$ and $$S^2$$ is the unit sphere in $$\mathbb{R}^3$$. Furthermore, let $$\beta(s,x)$$ denote the unit vector that “points” to $$x$$ from the source $$y(s)$$, i.e., $$\beta(s,x) = \frac{x - y(s)}{|x - y(s)|}$$ for $$x \in U$$ and $$y(s) \in C$$. PI lines are an important concept in constant pitch case. A PI line is a line segment connecting two points on the helix that are separated by less than $$2 \pi$$ [1]. As in the constant pitch case [2], [9] shows that for any $$x \in U$$ there exists a unique PI line of a dynamic pitch helix containing the point $$x$$. We let $$s = s_b(x)$$ and $$s = s_t(x)$$ denote the two source positions that define the PI line containing $$x$$, and $$I_{PI}(x) := [s_b(x), s_t(x)]$$ is the corresponding PI parametric interval (see Figure 1).

Consider filtering planes $$\Pi_{FLT}(s,s_2)$$ defined by three points along the helical trajectory $$s < s_1 < s_2$$ or $$s_2 < s_1 < s$$, with normal vector

$u(s,s_2) = \frac{(y(s_1) - y(s)) \times (y(s_2) - y(s))}{|(y(s_1) - y(s)) \times (y(s_2) - y(s))|} \textrm{sgn}(s_2 - s), \ 0 < | s_2 - s | < 2\pi. \tag3$

Here we assume that $$s_1 = (s + s_2)/2$$. For each $$x$$ and $$s\in I_{PI}(x)$$, find $$s_2$$ such that $$x\in \Pi_{FLT}(s,s_2)$$ (see Figure 1). The dynamic pitch reconstruction formula for $$f \in C_0^\infty(U)$$ is

$f(x) = -\frac{1}{2\pi^2} \int_{I_{PI}(x)} \frac{1}{|x-y(s)|} \int_{0}^{2\pi} \frac{\partial}{\partial q} D_f(y(q), \Theta(s,x,\gamma))|_{q=s} \frac{\textrm {d}{\gamma}}{\sin \gamma} \textrm {d}{s}, \tag4$

where

$\Theta(s,x,\gamma) := \cos \gamma \beta(s,x) + \sin \gamma e(s,x),\ e(s,x) := \beta(s,x) \times u(s,x). \tag5$

Figure 3. Reconstructed cross-section of a log showing identified defects and a virtual cut. © Microtec.
The condition $$\psi'+\psi'''>0$$ does not always hold in practical situations, and a suitably-adapted approximate version is used when that condition is violated.

On the hardware side, a wide X-ray detector with a surface of 1.25 m2 enables a wide cone angle scanner, which is capable of reconstructing at 180 m/min with a rotational speed of 4 rev/sec. This implies a helical pitch of 66 cm and a cone angle of about 30°. The resolution required for the reconstructed images is pretty low relative to other CT applications: 1 × 1 mm in the plane orthogonal to the axis of the helix and 10 mm along the axis. For this reason, Microtec has developed custom X-ray sensors with a sensitive area of 7 × 2 mm, assembled in an array of 47 × 758 pixels. The sensors were read at 1,500 frames per second (fps), to be back-projected in a volume of 900 × 900 × 60 voxels. As long as the back-projection of each voxel requires 20 floating point operations, the required computational power is 1.5 teraflops, with 1.2 terabytes/sec of memory access. The implementation of the algorithm in CUDA [11] using three GPUs GTX690 in parallel delivered the first computational solution  to produce a system capable of generating a continuous stream of 900 × 900-pixel slice images: one slice per 1 cm of log, equivalent to a slice every 3.3 msec.

Figure 2 is an image of a CT Log scanner, Figure 3 shows a typical reconstructed cross-section of a log with identified defects and a superimposed virtual cut, and Figure 4 provides an example of an optimal virtual cut based on the CT image.

A number of challenges regarding mechanical constraints, data transfer, safety regulations, image processing, and optimization algorithms have been addressed. Easily installable in any sawmill, the CT Log helps loggers get the most out of each sawn tree. Currently five sawmills in Europe and the Americas have a CT Log able to scan and optimize the processing of all logs at each mill.

Figure 4. Example of an optimized virtual cut obtained from CT reconstruction of a log. © Microtec.

Alexander Katsevich and Federico Giudiceandrea received the 2016 Marcus Wallenberg Prize for their development of a CT scanner for whole tree logs. They received their diploma from His Majesty the King of Sweden at a ceremony in Stockholm, Sweden this past October. The Marcus Wallenberg Prize recognizes path-breaking scientific achievements, which contribute significantly to broadening knowledge and to technical development within the fields that benefit forestry and forest industries.

Acknowledgments: The work of Alexander Katsevich on the algorithms reported in this paper was supported by NSF under the grants DMS-0104033 and DMS-0505494, and by grants from General Electric Healthcare and Microtec. Katsevich is also grateful for the University of Central Florida’s Office of Research, and iTomography’s help and support in facilitating the collaboration with Microtec.

References
[1] Danielsson, P.E., Edholm, P., Eriksson, J., & Magnusson, S.M. (1997). Towards exact reconstruction for helical cone-beam scanning of long objects. A new detector arrangement and a new completeness condition. In D.W. Townsend and P.E. Kinahan (Eds.), Proceedings from the 1997 Meeting on Fully 3D Image Reconstruction in Radiology and Nuclear Medicine (pp. 141-144). Pittsburgh, PA.

[2] Defrise, M., Noo, F., & Kudo, H. (2000).  A solution to the long-object problem in helical cone-beam tomography. Physics in Medicine and Biology, 45, 623-643.

[3] Feldkamp, L.A., Davis, L.C., & Kress, J.W. (1984). Practical cone-beam algorithm. JOSA A, 1(6), 612-619.

[4] Giudiceandrea, F., Ursella, E., & Vicario, E. (2011). A high speed CT scanner for the sawmill industry. Proceedings of the 17th International Non-Destructive Testing and Evaluation of Wood Symposium (pp. 105-112). Sopron, Hungary.

[5] Kapralov, M., & Katsevich, A. (2006). A one-PI algorithm for helical trajectories that violate the convexity condition. Inverse Problems, 22, 2123-2143.

[6] Katsevich, A. (2002). Analysis of an exact inversion algorithm for spiral cone-beam CT. Physics in Medicine and Biology, 47, 2583-2598.

[7] Katsevich, A. (2002). Theoretically exact filtered backprojection-type inversion algorithm for Spiral CT. SIAM Journal on Applied Mathematics, 62, 2012-2026.

[8] Katsevich, A. (2004). An improved exact filtered backprojection algorithm for spiral computed tomography. Advances in Applied Mathematics, 32, 681-697.

[9] Katsevich, A., Basu, S., & Hsieh, J. (2004). Exact filtered backprojection reconstruction for dynamic pitch helical cone beam computed tomography. Physics in Medicine and Biology, 49, 3089-3103.

[10] Katsevich, A., & Kapralov, M. (2007). Filtered backprojection inversion of the cone beam transform for a general class of curves. SIAM Journal on Applied Mathematics, 687, 334-353.

[11] Nickolls, J., Buck, I., Garland, M., & Skadron, K. (2008). Scalable parallel programming with CUDA. Queue – GPU Computing, 6(2), 40-53.

 Federico Giudiceandrea co-founded the company Microtec in 1980. Since then, he has been President, CEO, and Director of Research and Development at Microtec, leading the development of a number of products applying computer vision and X-rays to wood and fruit processing. He has received many awards, including the Schweighofer Prize in 2013. Alexander Katsevich joined the Department of Mathematics at the University of Central Florida (UCF) in 1996, and has been a full professor there since 2008. Since 2011 he has also served as Chief Technology Officer of iTomography Corporation, a spin-off from UCF that he co-founded to develop and commercialize new and existing CT technologies. Enrico Ursella co-founded SeeLab, a spin-off applying computer vision to wood scanning, in 1994. He has been with the Microtec group since 1997, where he has worked mainly in computer vision industrial applications. He led the team developing the CT Log tomographic scanner.