This is going to be a short post, because data objects are rarely useful in tomosipo.
First, we setup the geometries from the ASTRA code example:
import astra import tomosipo as ts import numpy as np # Create ASTRA geometries vol_geom = astra.create_vol_geom(256, 256) angles = np.linspace(0,np.pi,180,False) proj_geom = astra.create_proj_geom('parallel', 1.0, 384, angles) # Create tomosipo geometries vg = ts.volume(shape=(1, 256, 256)) pg = ts.parallel(angles=180, shape=(1, 384))
Transcribing the example
ASTRA works with data objects that must be created before use and freed after use. This script creates some data objects:
# Create volumes # initialized to zero v0 = astra.data2d.create('-vol', vol_geom) # initialized to 3.0 v1 = astra.data2d.create('-vol', vol_geom, 3.0) # initialized to a matrix. A may be a single, double or logical (0/1) array. array = np.random.normal(size=(256,256)) v2 = astra.data2d.create('-vol', vol_geom, array) # Projection data s0 = astra.data2d.create('-sino', proj_geom) # Initialization to a scalar or a matrix also works, exactly as with a volume.
The same can be accomplished in tomosipo:
# Create a 3D data object initialized to zero vd0 = ts.data(vg) # Initialized to 3.0: vd1 = ts.data(vg, 3.0) # Initialized to a 3D array: array64 = np.random.normal(size=(1, 256, 256)) vd2 = ts.data(vg, array64) # The previous code issues a warning, because array contains float64 values. # ASTRA operates only on float32 values. The array has been copied and converted # to 32 bits. # This code does not issue warnings: array32 = array64.astype(np.float32) vd2b = ts.data(vg, array32) # For projection geometries, the code is exactly the same. There is no need to # specify '-vol' or '-sino'. This is automatically inferred from the type of the # geometry.
Both tomosipo and ASTRA convert 64-bit floating point values to 32-bit floating point values. Tomosipo warns when this happens.
ASTRA has a
store method to update its data:
# set to zero astra.data2d.store(v0, 0) # set to a matrix astra.data2d.store(v2, array)
In tomosipo, the linked data can be accessed through the
data property. The
property itself cannot be overwritten, so you need to use
[:] to overwrite the
vd0.data[:] = 0.0 vd0.data[:] = array32
get method to get data:
R = astra.data2d.get(v2)
In ASTRA do not forget to free data:
# Free memory astra.data2d.delete(v0) astra.data2d.delete(v1) astra.data2d.delete(v2) astra.data2d.delete(s0)
In tomosipo, you can make sure that a data object is freed using a
with vd0: # start using vd0 vd0.data[:] += 1.0 # Do something with data array_clone = np.clone(vd0.data) # At this point, the data is deleted. # use array_clone from here on array_clone -= 1.0
When to use data objects in tomosipo
In tomosipo, it is rarely useful to create data objects, since forward and backprojection can occur on arrays directly — without the need to create an intermediate data object.
Operators have a
range_shape property that specify the
shape of their domain and range. This can be used as follows:
A = ts.operator(vg, pg) x = np.ones(A.domain_shape, dtype=np.float32) y = np.zeros(A.range_shape, dtype=np.float32)
The reason that the properties use the
range terms, is because
the backprojection operator has these same properties:
W = A.T y = np.zeros(A.domain_shape, dtype=np.float32) x = np.zeros(A.range_shape, dtype=np.float32)
Numpy arrays are automatically freed when they fall out of use. To force their removal, use
There is one use case for data objects that is useful for interactive use: displaying data in a window. You can use the following:
from tomosipo.qt import display # Create hollow box phantom vd = ts.phantom.hollow_box(ts.data(vg)) # Display it in an interactive window display(vd)
The slicer in the window allows you to slice through the volume.
You rarely need data objects in tomosipo, but when you need them
- use a
withblock for cleanup
- display data with