Chips and 3D Reconstruction Algorithms in X-ray Tomography Ryan A. - - PowerPoint PPT Presentation

AIP2009 Vienna

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) : x2 + y2 < 1}.

AIP2009 Vienna

Outline

1

Calculating π-lines with Chips

2

Chip Applications

3

Filtered Backprojection and Backprojection Filtration Methods

4

Outlook

AIP2009 Vienna

Helical CT does not reconstruct a point from all available source positions. x w IPI(x) y(st) y(sb) y(s) Instead our reconstruction at a point x depends on the projection data available from s ∈ [sb, st] = Iπ(x) the π-interval

• f x.

AIP2009 Vienna

Here sb and st are separated by no more than 2π and x lies on the line through y(sb) and y(st). We call this line the π-line. x w IPI(x) y(st) y(sb) 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

• Iπ(x)

G(x, s) ds ≈

• j

τ(Iπ(x), j∆s)G(x, j∆s). 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) = [sb, st]. If τ(Iπ(x), s) = 1Iπ(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], τ(Iπ(x), s) =    1/2 s ≈ sb, st 1 s ∈ [sb + ∆s, st − ∆s] s / ∈ (sb − ∆s, st + ∆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

x y z

• 0.5

0.5

• 0.5

0.5 5 10 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

• r 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, s3) we have d(C(s), x) = z − ρ sin(s − γ)((s + α)3 − (s − α)3) 2

• 1 − ρ2 cos2(s − γ)

+ (s + α)3 + (s − α)3 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

2

Chip Applications

3

Filtered Backprojection and Backprojection Filtration Methods

4

Outlook

AIP2009 Vienna

Here is a reconstruction, by Katsevich’s helical formula, of a piecewise constant function in the plane x3 = 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 x3 = 0. Katsevich calls the error a comet tail artifact.

AIP2009 Vienna

Comet Tails

Smooth function The artifact appears with smooth functions. We have found the artifact appears strongly with discontinuous phantoms with slightly misaligned data.

AIP2009 Vienna

Comet Tails

Shift of 1/2 detector bin The artifact appears with smooth functions. We have found the artifact appears strongly with discontinuous phantoms with slightly misaligned data.

AIP2009 Vienna

Let RBP(s) = {x : s ∈ Iπ(x)} be the region of backprojection

• f 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 x3 = 0.

AIP2009 Vienna

Consider the chip anchored at s = 0 projected onto the plane x3 = 0. The π-lines are parallel to the x2 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). fk RBP(−1) fk − fk−1 Here fk is the intermediate reconstruction up to source position sk = −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). fk RBP(2.2) fk − fk−1 Here fk is the intermediate reconstruction up to source position sk = 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 x3 = 0

• case. The artifact appears strongest in the direction of the

function’s π-lines.

AIP2009 Vienna

Outline

1

Calculating π-lines with Chips

2

Chip Applications

3

Filtered Backprojection and Backprojection Filtration Methods

4

Outlook

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

• Iπ(x)

2π ρ(x, α, s) ∂g ∂s + ∂g ∂α

• 1

sin(α∗ − α) dα ds. (2) 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

Hf(z) =

• Iπ(x)

ν(z, s) ∂g ∂s + ∂g ∂α

• ds

z ∈ π-line of x. (3) 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 sb or st 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

Reconstruction Formula Comments

Both formulas are exact and use data that is truncated along the helical axis. Each formula has been extended to more general scanning curves. The two methods are not related by interchanging the

• rder of integration.

Both rely heavily on Hilbert transforms.

AIP2009 Vienna

Reconstruction Formula Comments

The BPF formula cannot reconstruct a function at a point without calculating the function along the entire π-line of the point. Katsevich states that his FBP method is faster.

AIP2009 Vienna

Organizing the BPF

There is a dilemma in utilizing the BPF method. We typically reconstruct a volume by decomposing it into planes, Pi, parallel to the helical axis. Using a π-line for each point in Pi ignores the fact that the π-line will pass through other planes Pj. Want to avoid redundant calculations of f on π-lines that intersect more than one plane Pi. The chip provides a natural coordinate system for the helix by identifying the anchor point, π-interval endpoint, and the distance along π-line.

AIP2009 Vienna

Lets address how we can use chips to simplify the implementation of the BPF.

1 Decompose the volume of the region of interest into chips

• f size M by N such that each row in the chip corresponds

to points on the same π-line.

2 For each source position backproject the derivative data for

all points in RBP(s).

3 Apply the Hilbert inversion formula to every row of every

chip.

4 Interpolate from the chip coordinate system to the

cartesian system of the region of interest.

AIP2009 Vienna

Is the FBP faster than the BPF? Fix our spatial step size and reconstruct the entire support of a f. Suppose we vary the number of source positions and detector element size. The BPF cannot reconstruct at a point without reconstructing over the entire π-line of the point. On a π-line, the BPF will require the same number of backprojections as the FBP but will require only a single filtering step. The BPF will be faster on a chip than the FBP method. In 3D tomography the goal is to reconstruct a volume, not a slice. If we decompose the helical cylinder into chips, the BPF method will be faster than the FBP .

AIP2009 Vienna

Time Trials

768 Source Positions per turn, 192 by 384 detector elements. FBP has 138 filtering curves. Reconstruction on a 512 by 512 pixels on a chip. FBP runtime BPF runtime 131 sec 78 sec

AIP2009 Vienna

Time Trials

bpf fbp The accuracy for reconstructing discontinuous functions is identical.

AIP2009 Vienna

Chip Shortcomings

FBP would be the preferred choice for reconstructing a single slice of a volume or a ROI within the support of f. x y z

• 0.5

0.5

• 0.5

0.5 5 10 A difficulty with the chips is the grid we are reconstructing over does not have uniform spacing in z.

AIP2009 Vienna

But we can put restrictions on the pitch of the helix, the scanning radius and the spatial step size ∆x to ensure that for each plane Pi we can find a chip such that the distance between the chip and Pi is bounded by ∆x. The idea is that a chip can forced to behave as Pi around the helical axis for some pitch. The smaller ∆x gets the smaller the pitch must be.

AIP2009 Vienna

Outline

1

Calculating π-lines with Chips

2

Chip Applications

3

Filtered Backprojection and Backprojection Filtration Methods

4

Outlook

AIP2009 Vienna

We can extend the variable pitch case to the general source curve where we can define a generalized helical axis a(s) and local coordinate system u(s), v(s), and w(s). x a(s) Iπ(x) y(st) y(sb) y(s) u(s) v(s) w(s)

AIP2009 Vienna

We demonstrate our technique on the following curve as a proof of concept. y(s) = (R cos s + cos(2πs) cos s, sin(2πs), R sin s + cos(2πs) sin s) a(s) = (R cos s, 0, R sin s) The source curves spins around a circle in the x − z plane. z y x Chip with anchor s = 5.69 and R = 2.

AIP2009 Vienna

Summary

Calculate the π-lines for general scanning curves. Chips are a natural coordinate system for helical tomography.

Identified all the points in the region of backprojection. Describe comet tail artifacts on chips. BPF can be faster than FBP .

Thank You