Jul 21, 2020

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

### Creating data

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.

### Changing data

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 data.

vd0.data[:] = 0.0
vd0.data[:] = array32


### Retrieving data

Use ASTRA’s get method to get data:

R = astra.data2d.get(v2)


Or tomosipo’s data property:

vd0.data


### Freeing data

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 block:

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 domain_shape and 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 domain and 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

del x


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.

## Summary

You rarely need data objects in tomosipo, but when you need them

• use ts.data
• use a with block for cleanup
• display data with tomosipo.qt.display