Advent of tomosipo: s002_data2d

Jul 21, 2020

Allard Hendriksen

3 minute read

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

Allard Hendriksen | CI @ CWI