Jul 29, 2020

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:

1. import from astra
2. convert the parallel beam geometry to a vector geometry in tomosipo
3. 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:

1. import from astra
2. convert the parallel beam geometry to a vector geometry in tomosipo
3. 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!