Advent of tomosipo: s006_3d_data
The subject of today is 3D data. As for the s002_data2d post, this is going to be short, because data objects are rarely useful in tomosipo.
First, we import the volume and projection geometry from the ASTRA code example:
import astra
import numpy as np
import tomosipo as ts
# Create a 3D volume geometry.
# Parameter order: rows, colums, slices (y, x, z)
vol_geom = astra.create_vol_geom(64, 48, 32)
# 2 projection directions, along x and y axis resp.
V = np.array([[ 1,0,0, 0,0,0, 0,1,0, 0,0,1],
[0,1,0, 0,0,0, -1,0,0, 0,0,1]],dtype=np.float)
# 32 rows (v), 64 columns (u)
proj_geom = astra.create_proj_geom('parallel3d_vec', 32, 64, V)
# Import geometries into tomosipo
vg = ts.from_astra(vol_geom)
pg = ts.from_astra(proj_geom)
# We create an operator, which we use later for its domain_shape
# and range_shape properties.
A = ts.operator(vg, pg)
Volume data
In ASTRA, volume data can be created filled with
- zeros (this is the default)
- a fixed value (3.0 in the example below)
- a numpy array (which is silently converted to float32 if it contains float64 values)
# initialized to zero
v0 = astra.data3d.create('-vol', vol_geom)
# initialized to 3.0
v1 = astra.data3d.create('-vol', vol_geom, 3.0)
# initialized to a matrix. A may be a single or double array.
# Coordinate order: slice, row, column (z, y, x)
array_astra = np.zeros((32, 64, 48))
v2 = astra.data3d.create('-vol', vol_geom, array_astra)
In tomosipo, the same is possible:
# 1. Create 0-filled data from volume geometry
vd0 = ts.data(vg)
# 2. Create 2..0-filled data from volume geometry
vd1 = ts.data(vg, 2.0)
# 3. Create data filled with specific array
# - with warning when float64 is automatically converted:
array = np.random.normal(size=A.domain_shape)
vd2 = ts.data(vg, array)
# - without warning with float32:
array32 = array.astype(np.float32)
vd3_32 = ts.data(vg, array32)
/lib/python3.6/site-packages/tomosipo/links/numpy.py:25: UserWarning:
The parameter initial_value is of type float64; expected `np.float32`.
The type has been Automatically converted.
Use `ts.link(x.astype(np.float32))' to inhibit this warning.
f"The parameter initial_value is of type {initial_value.dtype};
expected `np.float32`. "
In these examples, and throughout the source code of tomosipo, a volume geometry
is often abbreviated as vg and a projection geometry as pg. For data, we use
vd and pd.
Projection data
In ASTRA, there is a difference between projection and volume data. Tomosipo can automatically infer this from the provided geometry.
s0 = astra.data3d.create('-proj3d', proj_geom)
# Initialization to a scalar or zero works exactly as with a volume.
# Initialized to a matrix:
# Coordinate order: row (v), angle, column (u)
array_sino = np.zeros((32, 2, 64))
s1 = astra.data3d.create('-proj3d', proj_geom, array_sino)
In tomosipo, not much is different from the volume geometry example:
pd0 = ts.data(pg)
print(array_sino.dtype)
# When you repeatedly create a 64-bit data object, the automatic
# conversion warning may go away at some point, as it does here
pd1 = ts.data(pg, array_sino)
float64
Retrieving data
In ASTRA, you may retrieve a mutable instance of the data using:
# Retrieve data:
R = astra.data3d.get(v1)
The same is possible in tomosipo:
R = vd1.data
# We can mutate the data like so:
R[:] = 2.0
# More natural would probably be:
vd1.data[:] = 2.0
Note the [:] in vd1.data[:]: this is needed and tomosipo will raise an error
if you forget.
Cleaning up
This is how you clean up your data in ASTRA:
# Delete all created data objects
astra.data3d.delete(v0)
astra.data3d.delete(v1)
astra.data3d.delete(v2)
astra.data3d.delete(s0)
astra.data3d.delete(s1)
And this is how you do it in tomosipo:
with v0, v1, v2, s0, s1:
pass
As you can see, cleaning up data in tomosipo is not very straightforward. In fact, it is not advised to use data objects for most tasks. For compatibility with ASTRA, and some historical reasons, however, they have been included in tomosipo. For more info, see When to use data objects in tomosipo.
Thanks to Francien Bossema for spotting a typo in this post!
