Chips and 3D Reconstruction Algorithms in X-ray Tomography Ryan A. Hass and Adel Faridani Department of Mathematics Oregon State University Supported by NSF grant DMS-0709495 July 23rd, 2009 AIP2009 Vienna
Computed Tomography Our goal in computed tomography (CT) is to recover a function f from its line integrals. This is the framework of a CAT scan. Here the line integrals are measured by a rotating detector array and an x-ray source. AIP2009 Vienna
Our path around the object is a helix y ( s ) = [ cos ( s ) , sin ( s ) , hs ] T S y ( s ) The object to be imaged is supported inside the helical cylinder S = { ( x , y , z ) : x 2 + y 2 < 1 } . AIP2009 Vienna
Outline 1 Calculating π -lines with Chips Chip Applications 2 Filtered Backprojection and Backprojection Filtration Methods 3 Outlook 4 AIP2009 Vienna
Helical CT does not reconstruct a point from all available source positions. w y ( s t ) I PI ( x ) x y ( s b ) y ( s ) Instead our reconstruction at a point x depends on the projection data available from s ∈ [ s b , s t ] = I π ( x ) the π -interval of x . AIP2009 Vienna
Here s b and s t are separated by no more than 2 π and x lies on the line through y ( s b ) and y ( s t ) . We call this line the π -line . w y ( s t ) I PI ( x ) x y ( s b ) y ( s ) The projection data from within a π -interval gives us a 180 ◦ view of the point x . It is known that we need at least 180 ◦ view to recover x [Quinto]. AIP2009 Vienna
Lemma π -intervals exist and are unique [Defrise, Noo, and Kudo 2000]. The 3D formulas hold with more general curves that act as a helix locally. Lemma π -intervals exist and are unique for generalized helical source curves [Katsevich and Kapralov 2007]. AIP2009 Vienna
In helical CT we compute the backprojection term � � G ( x , s ) ds ≈ τ ( I π ( x ) , j ∆ s ) G ( x , j ∆ s ) . I π ( x ) j We do not have m , n such that I π ( x ) = [ m ∆ s , n ∆ s ] for all x . AIP2009 Vienna
Our choice of τ depends on I π ( x ) = [ s b , s t ] . If τ ( I π ( x ) , s ) = 1 I π ( x ) ( s ) then we are don’t need to calculate the π -intervals. Instead we can use a geometric characterization of the detector coordinates of x [Katsevich 2004]. However this characterization, or choice of τ , will lead to higher errors, artifacts, and a lower rate of convergence. AIP2009 Vienna
A proper choice of τ is, [Noo Pack Heuscher 2003], 1 / 2 s ≈ s b , s t τ ( I π ( x ) , s ) = 1 s ∈ [ s b + ∆ s , s t − ∆ s ] . 0 s / ∈ ( s b − ∆ s , s t + ∆ s ) This interpolating function allows us to use the trapezoidal rule for the backprojection. We must calculate I π ( x ) prior to the backprojection. AIP2009 Vienna
We follow the work of [Izen 2007] to calculate the π -interval for a point x when the source curve is a helix. z y x A chip C ( s ) is a set of all π -lines with a common midpoint in their π -interval. The common midpoint is the chip’s anchor s . AIP2009 Vienna
We follow the work of [Izen 2007] to calculate the π -interval for a point x when the source curve is a helix. z y x A chip C ( s ) is a set of all π -lines with a common midpoint in their π -interval. The common midpoint is the chip’s anchor s . AIP2009 Vienna
Helical chips We can decompose the helical scanning cylinder into disjoint chips and determine which chip contains x . This yields I π ( x ) . For the standard helix all the chips are rotated and translated copies of the chip anchored at s = 0. We lose this relation if the pitch or radius of the helix varies with s . We have easily extended the idea of a chip to helical curves with a variable pitch. In each case the π -lines on a chip lie in parallel planes. AIP2009 Vienna
Stack of chips 10 5 z 0 0.5 0.5 0 0 -0.5 -0.5 y x A stack of chips for a variable pitch helix. But how did we identify a chip for any x ? AIP2009 Vienna
Calculating π -lines Suppose y ( s ) is given and admits unique π -lines for all x ∈ S . Let d ( C ( s ) , x ) be the distance from C ( s ) to x . Find C ( s ) such that x ∈ C ( s ) 1 Given s determine d ( C ( s ) , x ) by either an explicit formula or by minimizing the distance over π -lines in C ( s ) to x . 2 Numerically find s such that d ( C ( s ) , x ) = 0 3 Identify α , I π ( x ) = [ s − α, s + α ] . Our method is different than Izen’s. AIP2009 Vienna
For x = ( ρ, γ, z ) and y ( s ) = ( cos s , sin s , s 3 ) we have d ( C ( s ) , x ) = z − ρ sin ( s − γ )(( s + α ) 3 − ( s − α ) 3 ) + ( s + α ) 3 + ( s − α ) 3 1 − ρ 2 cos 2 ( s − γ ) � 2 2 (1) where α = cos − 1 ( ρ cos ( s − γ )) . We have yet to show that d − 1 ( · , x ) exists and is smooth. AIP2009 Vienna
Outline 1 Calculating π -lines with Chips Chip Applications 2 Filtered Backprojection and Backprojection Filtration Methods 3 Outlook 4 AIP2009 Vienna
Here is a reconstruction, by Katsevich’s helical formula, of a piecewise constant function in the plane x 3 = 0. Katsevich calls the error a comet tail artifact. AIP2009 Vienna
Here is a reconstruction, by Katsevich’s helical formula, of a piecewise constant function in the plane x 3 = 0. Katsevich calls the error a comet tail artifact. AIP2009 Vienna
Comet Tails The artifact appears with smooth functions. We have found the artifact appears strongly with discontinuous phantoms with slightly misaligned data. Smooth function AIP2009 Vienna
Comet Tails The artifact appears with smooth functions. We have found the artifact appears strongly with discontinuous phantoms with slightly misaligned data. Shift of 1/2 detector bin AIP2009 Vienna
Let RBP ( s ) = { x : s ∈ I π ( x ) } be the region of backprojection of s . The RBP ( s ) is the set of the points that we reconstruct from the current source position. The boundary of RBP ( s ) is important. Here the source position is at the endpoint of the π -interval. In other words we have projection data along the π -line. AIP2009 Vienna
RBP ( − 1 ) RBP ( 2 . 2 ) A description of this set is useful in describing the location of comet tail artifacts. It is difficult to get a description of RBP ( s ) in the plane x 3 = 0. AIP2009 Vienna
Consider the chip anchored at s = 0 projected onto the plane x 3 = 0. The π -lines are parallel to the x 2 axis. RBP ( − 2 ) RBP ( . 75 ) On each chip, RBP ( s ) expands towards the anchor point and then retracts. AIP2009 Vienna
A description of RBP ( s ) on any chip is now possible. Lemma Let C ( t ) be a chip anchored at t and suppose y ( s ) = ( cos s , sin s , hs ) . Then RBP ( s ) ∩ C ( t ) is S ∩ { x ∈ C ( t ) | I π ( x ) = [ t − φ, t + φ )] , | t − s | < φ < π } AIP2009 Vienna
We know that the helical cylinder S can be decomposed into chips and we can describe the region of backprojection on each chip. Thus we have a description of the region of backprojection for the entire helical cylinder S . This theory extends directly to any curve where we have a notation of chips, includes circle plus line. But what about the comet tail? AIP2009 Vienna
The idea is the artifact occurs at the boundary of RBP ( s ) . This corresponds to the cut off of the integration over I π ( x ) . f k RBP ( − 1 ) f k − f k − 1 Here f k is the intermediate reconstruction up to source position s k = − 1. AIP2009 Vienna
The idea is the artifact occurs at the boundary of RBP ( s ) . This corresponds to the cut off of the integration over I π ( x ) . f k RBP ( 2 . 2 ) f k − f k − 1 Here f k is the intermediate reconstruction up to source position s k = 2 . 2. AIP2009 Vienna
Same object but reconstruction on the chip anchored at s = 0. The artifact is now much easier to describe than the x 3 = 0 case. The artifact appears strongest in the direction of the function’s π -lines. AIP2009 Vienna
Outline 1 Calculating π -lines with Chips Chip Applications 2 Filtered Backprojection and Backprojection Filtration Methods 3 Outlook 4 AIP2009 Vienna
CT geometry Our X-ray data is given by g ( s , α, w ) , where s is the position, α is the fan angle component, w is the detector height component. w ( α, w ) α y ( s ) AIP2009 Vienna
Reconstruction formula Katsevich’s formula in curved detector coordinates is � 2 π � � ∂ g ∂ s + ∂ g � 1 ρ ( x , α, s ) sin ( α ∗ − α ) d α ds . (2) ∂α I π ( x ) 0 Inner integral represents filtering the measured data along fixed curves on the detector surface for each source position. Outer integral backprojects the filtered data over source positions within the π -interval, I π ( x ) . We call such a formula a filtered backprojection method (FBP). AIP2009 Vienna
X. Pan’s formula in curved detector coordinates is � ∂ g � � ∂ s + ∂ g H f ( z ) = ν ( z , s ) ds ∂α (3) I π ( x ) z ∈ π -line of x . First backproject the differentiated data but apply no filtering. Invert the Hilbert data along the π -line of x with a weighted Hilbert transform. A constant that depends on the projection data at s b or s t is needed. The formula for the inversion comes from classical integral equations. We call such a formula a backprojection filtration (BPF), π -line filtration or a differentiated backprojection method. AIP2009 Vienna
Recommend
More recommend
