Advent of tomosipo: s005_3d_geometry

Jul 29, 2020

Allard Hendriksen

7 minute read

We have finally moved on to 3D geometries! ASTRA supports parallel-beam and cone-beam geometries. For both geometries, it supports a regular geometry that is useful in most cases, and a “vector” geometry that can be used to model arbitrary geometries.

Regular parallel beam

In ASTRA, you can define a regular 3D parallel beam geometry like so. This one rotates in a full circle around the z-axis (rather than 180°) stopping at 48 locations. The pixels are square of size 1.0, and the detector is 32 pixels high and 64 pixels wide.

import astra
import numpy as np
# There are two main 3d projection geometry types: cone beam and parallel beam.
# Each has a regular variant, and a 'vec' variant.
# The 'vec' variants are completely free in the placement of source/detector,
# while the regular variants assume circular trajectories around the z-axis.

# -------------
# Circular Parallel beam
# -------------
# Parameters: width of detector column, height of detector row, #rows, #columns
angles = np.linspace(0, 2*np.pi, 48, False)
proj_geom = astra.create_proj_geom('parallel3d', 1.0, 1.0, 32, 64, angles)

We can import this geometry in tomosipo, and create the geometry using ts.parallel:

import tomosipo as ts
# 1. import from astra
pg_from_astra = ts.from_astra(proj_geom)
# 2. create it ourselves
pg = ts.parallel(angles=angles, shape=(32, 64))
# pixel size is 1.0 by default in tomosipo, but we can
# specify the size of the detector as well:
pg = ts.parallel(angles=angles, shape=(32, 64), size=(32, 64))

assert pg == pg_from_astra      # ensure geometries are equal
pg

ts.parallel(
    angles=array([0.        , 0.13089969, 0.26179939, 0.39269908, 0.52359878,
       0.65449847, 0.78539816, 0.91629786, 1.04719755, 1.17809725,
       1.30899694, 1.43989663, 1.57079633, 1.70169602, 1.83259571,
       1.96349541, 2.0943951 , 2.2252948 , 2.35619449, 2.48709418,
       2.61799388, 2.74889357, 2.87979327, 3.01069296, 3.14159265,
       3.27249235, 3.40339204, 3.53429174, 3.66519143, 3.79609112,
       3.92699082, 4.05789051, 4.1887902 , 4.3196899 , 4.45058959,
       4.58148929, 4.71238898, 4.84328867, 4.97418837, 5.10508806,
       5.23598776, 5.36688745, 5.49778714, 5.62868684, 5.75958653,
       5.89048623, 6.02138592, 6.15228561]),
    shape=(32, 64),
    size=(32, 64),
)

And this is how the result looks like:

from tomosipo.qt import animate

animate(ts.scale(1/32) * pg)

Arbitrary parallel beam geometry

We can define the same geometry using an ASTRA vector geometry:

# We generate the same geometry as the circular one above.
vectors = np.zeros((len(angles), 12))
for i in range(len(angles)):
  # ray direction
  vectors[i,0] = np.sin(angles[i])
  vectors[i,1] = -np.cos(angles[i])
  vectors[i,2] = 0

  # center of detector
  vectors[i,3:6] = 0

  # vector from detector pixel (0,0) to (0,1)
  vectors[i,6] = np.cos(angles[i])
  vectors[i,7] = np.sin(angles[i])
  vectors[i,8] = 0;

  # vector from detector pixel (0,0) to (1,0)
  vectors[i,9] = 0
  vectors[i,10] = 0
  vectors[i,11] = 1

# Parameters: #rows, #columns, vectors
proj_geom = astra.create_proj_geom('parallel3d_vec', 32, 64, vectors)

Now, there are three ways to obtain the same geometry in tomosipo:

import from astra
convert the parallel beam geometry to a vector geometry in tomosipo
create a parallel_vec geometry from scratch.

# 1. import:
pg_vec_from_astra = ts.from_astra(proj_geom)
# 2. convert tomosipo parallel geometry to tomosipo parallel vector geometry:
pg_to_vec = pg.to_vec()
# 3. create from scratch using ts.parallel_vec
ray_dir = np.zeros((48, 3))
ray_dir[:, 0] = 0.0             # z-direction
ray_dir[:, 1] = -np.cos(angles) # y-direction
ray_dir[:, 2] = np.sin(angles)  # x-direction

# vector from detector pixel (0,0) to (1,0)
det_v = np.zeros((48, 3))
det_v[:, 0] = 1.0               # z-direction (rest is zero)

# vector from detector pixel (0,0) to (0,1)
det_u = np.zeros((48, 3))
det_u[:, 0] = 0.0               # z-direction
det_u[:, 1] = np.sin(angles)    # y-direction
det_u[:, 2] = np.cos(angles)    # x-direction

pg_vec = ts.parallel_vec(
    shape=(32, 64),
    ray_dir=ray_dir,
    det_pos=np.zeros((48, 3)),
    det_v=det_v,
    det_u=det_u,
)

# Ensure that our geometries are the same
assert pg_vec_from_astra == pg_to_vec  # 1. == 2.
assert pg_vec == pg_to_vec             # 3. == 2.

Note that tomosipo always orders the axes Z, Y, X.

Regular circular cone beam

This is how to create a circular cone beam geometry in ASTRA:

# Parameters: width of detector column, height of detector row, #rows, #columns,
#             angles, distance source-origin, distance origin-detector
angles = np.linspace(0, 2*np.pi, 48, False)
proj_geom = astra.create_proj_geom('cone', 1.0, 1.0, 32, 64, angles, 1000, 0)

And this is how to do it in tomosipo

# 1. import from astra
pg_from_astra = ts.from_astra(proj_geom)
# 2. using ts.cone
pg = ts.cone(
    angles=angles,
    shape=(32, 64),
    src_orig_dist=1000.0,
    src_det_dist=1000.0,        #
)
assert pg == pg_from_astra
print(pg)

ts.cone(
    angles=array([0.        , 0.13089969, 0.26179939, 0.39269908, 0.52359878,
       0.65449847, 0.78539816, 0.91629786, 1.04719755, 1.17809725,
       1.30899694, 1.43989663, 1.57079633, 1.70169602, 1.83259571,
       1.96349541, 2.0943951 , 2.2252948 , 2.35619449, 2.48709418,
       2.61799388, 2.74889357, 2.87979327, 3.01069296, 3.14159265,
       3.27249235, 3.40339204, 3.53429174, 3.66519143, 3.79609112,
       3.92699082, 4.05789051, 4.1887902 , 4.3196899 , 4.45058959,
       4.58148929, 4.71238898, 4.84328867, 4.97418837, 5.10508806,
       5.23598776, 5.36688745, 5.49778714, 5.62868684, 5.75958653,
       5.89048623, 6.02138592, 6.15228561]),
    shape=(32, 64),
    size=(32, 64),
    src_orig_dist=1000.0,
    src_det_dist=1000.0,
)

And this is how the result looks like:

from tomosipo.qt import animate
animate(ts.scale(1 / 256) * pg)

Arbitrary cone beam geometry

This is how to create a cone vector geometry in ASTRA:

vectors = np.zeros((len(angles), 12))
for i in range(len(angles)):
	# source
	vectors[i,0] = np.sin(angles[i]) * 1000
	vectors[i,1] = -np.cos(angles[i]) * 1000
	vectors[i,2] = 0

	# center of detector
	vectors[i,3:6] = 0

	# vector from detector pixel (0,0) to (0,1)
	vectors[i,6] = np.cos(angles[i])
	vectors[i,7] = np.sin(angles[i])
	vectors[i,8] = 0

	# vector from detector pixel (0,0) to (1,0)
	vectors[i,9] = 0
	vectors[i,10] = 0
	vectors[i,11] = 1

# Parameters: #rows, #columns, vectors
proj_geom = astra.create_proj_geom('cone_vec', 32, 64, vectors)

Now, there are three ways to obtain the same geometry in tomosipo:

import from astra
convert the parallel beam geometry to a vector geometry in tomosipo
create a cone_vec geometry from scratch.

# 1. import:
pg_vec_from_astra = ts.from_astra(proj_geom)
# 2. convert tomosipo cone geometry to tomosipo cone vector geometry:
pg_to_vec = pg.to_vec()
# 3. create from scratch using ts.cone_vec
src_pos = np.zeros((48, 3))
src_pos[:, 0] = 0.0                     # z-direction
src_pos[:, 1] = -1000 * np.cos(angles)  # y-direction
src_pos[:, 2] = 1000 * np.sin(angles)   # x-direction

# vector from detector pixel (0,0) to (1,0)
det_v = np.zeros((48, 3))
det_v[:, 0] = 1.0               # z-direction (rest is zero)

# vector from detector pixel (0,0) to (0,1)
det_u = np.zeros((48, 3))
det_u[:, 0] = 0.0               # z-direction
det_u[:, 1] = np.sin(angles)    # y-direction
det_u[:, 2] = np.cos(angles)    # x-direction

pg_vec = ts.cone_vec(
    shape=(32, 64),
    src_pos=src_pos,
    det_pos=np.zeros((48, 3)),
    det_v=det_v,
    det_u=det_u,
)

# Ensure that our geometries are the same
assert pg_vec_from_astra == pg_to_vec  # 1. == 2.
assert pg_vec == pg_to_vec             # 3. == 2.

Helical cone beam

# First create a normal circular cone beam geometry
pg_ccb = ts.cone(angles=96, src_orig_dist=2.0, src_det_dist=4.0)

# Create a translation transformation consisting of 96 steps that
# starts by moving an object downward in the z-direction by 2.0
# units, and ends by moving it upward by 2.0 units.
translation = np.zeros((96, 3))
translation[:, 0] = np.linspace(-2, 2, 96)  # z-direction
T = ts.translate(translation)

# You can apply a geometric transformation by multiplying:
pg_helical = T * pg_ccb.to_vec()
# Note that we have converted to cone beam geometry to a vector geometry before
# translating. We do this because a regular cone beam geometry cannot represent
# the irregular source and detector positions of a helical geometry. The
# translation does not fail when you forget this: tomosipo converts the geometry
# under the hood, and gives a warning that the geometry has been converted to a
# vector geometry.

# The two geometries can be visualized as follows:
animate(pg_ccb, pg_helical)

The easiest way to obtain helical geometry in ASTRA, is by exporting from tomosipo as follows:

pg_astra = ts.to_astra(pg_helical)  # option 1
pg_astra = pg_helical.to_astra()    # option 2

Note that that the helical geometry as defined above is not that useful for tomography. First of all, it contains only a single detector pixel. In addition, the translation in the Z-direction might be too large to obtain a good reconstruction. These choices did improve the visualization of the geometries, which was the goal of this exercise.

Summary

All ASTRA geometries can be created and respresented in tomosipo as well. The geometries have a nice textual representation and can be shown as an animation. In addition, tomosipo has facilities for manipulating geometries, which makes it quite easy to create a helical cone beam geometry, for instance.

Thanks to Francien Bossema for proof-reading!

Allard Hendriksen | CI @ CWI