If you have trouble viewing this, try the pdf of this post. You can download the code used to produce the figures in this post.

edit: Sep 9, 2011 I extended the code to include convex polygons and renamed the function as *CTProjSim* .

I previously discussed the rationale, the C++ implmentation, and the the Matlab interface for a computed tomography projection simulator. In this post, I discuss a Matlab-only implementation of a simulator. The simulator is includes ellipses and convex polygons.The projection lines are assumed parallel but it is simple and can be (fairly) easily extended to other object types and geometries.

The geometry for the ellipses is shown in Fig. 1↓. With this geometry, the length of the intersection of the projection line with an ellipse is
*T* = ⎧⎨⎩
2*ab*√(*s*^{2}_{m} − *s*^{2}) ⁄ *s*^{2}_{m}
|*s*| ≤ *s*_{m}
0
|*s*| > *s*_{m}
where *s* is the distance of the line *L* from the center of the ellipse and *θ* is the angle from the perpendicular of the line to the principal axis of the ellipse. Also, the maximum offset of the ellipse perpendicular to the projection line is *s*^{2}_{m} = *a*^{2}cos^{2}(*θ*) + *b*^{2}sin^{2}(*θ*).The code of the function *CTProjSim* is a straightforward implementation of these formulas.

I used the function to compute the projections of the Shepp-Logan phantom and reconstructed with my CTrecon function. The results are shown in Fig. 2↓.

We can use concepts from my previous posts on intersection of line segments and polylines to compute the intersection of a line with a polygon. I assume the polygon is convex so the line can intersect with at most two segments of the polygon. The geometry is shown in Fig. 3↓. In two dimensions, we can specify a line by a vector **n** through the origin and perpendicular to the line. We can write **n** as *n***n**̂, where **n**̂ is a unit length vector. If the line passes through the origin, *n* = 0, but we still know its direction from **n**̂.

A point on the line is
**r**_{L} = *n***n**̂ + *t***s**̂
while a point on the line segment is
**r**_{S} = P_{1} + *u***d**.
At the intersection,
**r**_{L} = *n***n**̂ + *t***s**̂ = **r**_{S} = P_{1} + *u***d**.
Taking the dot product of this equation with **n**̂
*n* = **n**̂⋅**P**_{1} + *u***n**̂⋅d
Solving for *u*
*u* = (*n* − **n**̂⋅**P**_{1})/(**n**̂⋅d)

To find the intersection of a line with a polygon, we can test 0 ≤ *u* < 1 for all the segments of the polygon sides. If there are two intersections, then the line overlaps the polygon. We can compute the intersect length *T* as
*T* = |**r**_{intersect, 1} − **r**_{intersect, 2}|
Multiplying *T* by the (possibly vector) density gives the line integral.

These formulas are implemented in the *CTProjSim* along with the code for the intersection with ellipses. Line integrals are linear so we can add the line integrals for different shapes for each projection line.

I created the projections of a hexagon embedded in an ellipse and reconstructed them to produce the image in Fig. 4↓. The sharp vertexes of the polygon cause large aliasing artifacts in the reconstruction. I reduced them somewhat by filtering them with the ’*filter*_*type*’, ’*hamming*’, ’*freq*_*cutoff*’, 0.9 parameters in *CTrecon*. You can download the code to reproduce the figures in this post.

Last edited Sep. 9, 2011

© 2011 by Aprend Technology and Robert E. Alvarez

Linking is allowed but reposting or mirroring is expressly forbidden.