This page was generated from /doc/source/ipynb/tomopy.ipynb. Interactive online version: Binder badge

Reconstruction with TomoPy

Here is an example of how to use the gridrec [B5] reconstruction algorithm with TomoPy [A1].

First install TomoPy.

[1]:
import tomopy

Tomographic data input in TomoPy is supported by DXchange.

[2]:
import dxchange

Matplotlib provides plotting of the result in this notebook. Paraview or other tools are available for more sophisticated 3D rendering.

[3]:
import matplotlib.pyplot as plt

Import and activate Python’s built in logging module if desired. It may print something helpful.

[4]:
import logging
logging.basicConfig(level=logging.INFO)

This data set file format follows the APS beamline 2-BM and 32-ID data-exchange definition. Other file format readers for other synchrotrons are also available with DXchange.

[5]:
proj, flat, dark, theta = dxchange.read_aps_32id(
    fname='../../../source/tomopy/data/tooth.h5',
    sino=(0, 2),  # Select the sinogram range to reconstruct.
)
INFO:dxchange.reader:Data successfully imported: /home/dching/Documents/tomopy/source/tomopy/data/tooth.h5
INFO:dxchange.reader:Data successfully imported: /home/dching/Documents/tomopy/source/tomopy/data/tooth.h5
INFO:dxchange.reader:Data successfully imported: /home/dching/Documents/tomopy/source/tomopy/data/tooth.h5
INFO:dxchange.reader:Data successfully imported: /home/dching/Documents/tomopy/source/tomopy/data/tooth.h5

Plot the sinogram

[6]:
plt.imshow(proj[:, 0, :])
plt.show()
../_images/ipynb_tomopy_11_0.png

If the angular information is not avaialable from the raw data you need to set the data collection angles. In this case, theta is set as equally spaced between 0-180 degrees.

[7]:
if theta is None:
    theta = tomopy.angles(proj.shape[0])

Perform the flat-field correction of raw data:

\[\frac{proj - dark} {flat - dark}\]
[8]:
proj = tomopy.normalize(proj, flat, dark)

Calculate $ -log(proj) $ to linearize transmission tomography data.

[9]:
proj = tomopy.minus_log(proj)

Tomopy provides various methods ([B11], [B24], [B15]) to find the rotation center.

[10]:
rot_center = tomopy.find_center(proj, theta, init=290, ind=0, tol=0.5)
INFO:tomopy.recon.rotation:Trying rotation center: [290.]
INFO:tomopy.recon.rotation:Function value = 2.014651
INFO:tomopy.recon.rotation:Trying rotation center: [304.5]
Reconstructing 1 slice groups with 1 master threads...
Reconstructing 1 slice groups with 1 master threads...
Reconstructing 1 slice groups with 1 master threads...
INFO:tomopy.recon.rotation:Function value = 2.076837
INFO:tomopy.recon.rotation:Trying rotation center: [275.5]
INFO:tomopy.recon.rotation:Function value = 2.259117
INFO:tomopy.recon.rotation:Trying rotation center: [297.25]
INFO:tomopy.recon.rotation:Function value = 1.920647
INFO:tomopy.recon.rotation:Trying rotation center: [304.5]
Reconstructing 1 slice groups with 1 master threads...
Reconstructing 1 slice groups with 1 master threads...
Reconstructing 1 slice groups with 1 master threads...
INFO:tomopy.recon.rotation:Function value = 2.076837
INFO:tomopy.recon.rotation:Trying rotation center: [293.625]
INFO:tomopy.recon.rotation:Function value = 1.939667
INFO:tomopy.recon.rotation:Trying rotation center: [300.875]
INFO:tomopy.recon.rotation:Function value = 1.997986
INFO:tomopy.recon.rotation:Trying rotation center: [295.4375]
Reconstructing 1 slice groups with 1 master threads...
Reconstructing 1 slice groups with 1 master threads...
Reconstructing 1 slice groups with 1 master threads...
INFO:tomopy.recon.rotation:Function value = 1.908336
INFO:tomopy.recon.rotation:Trying rotation center: [293.625]
INFO:tomopy.recon.rotation:Function value = 1.939667
INFO:tomopy.recon.rotation:Trying rotation center: [296.34375]
INFO:tomopy.recon.rotation:Function value = 1.906685
INFO:tomopy.recon.rotation:Trying rotation center: [297.25]
Reconstructing 1 slice groups with 1 master threads...
Reconstructing 1 slice groups with 1 master threads...
Reconstructing 1 slice groups with 1 master threads...
INFO:tomopy.recon.rotation:Function value = 1.920647
INFO:tomopy.recon.rotation:Trying rotation center: [295.890625]
INFO:tomopy.recon.rotation:Function value = 1.906942
Reconstructing 1 slice groups with 1 master threads...

Reconstruct using the gridrec algorithm. Tomopy provides various reconstruction and provides wrappers for other libraries such as the ASTRA toolbox.

[11]:
recon = tomopy.recon(proj, theta, center=rot_center, algorithm='gridrec', sinogram_order=False)
Reconstructing 2 slice groups with 2 master threads...

Mask each reconstructed slice with a circle.

[12]:
recon = tomopy.circ_mask(recon, axis=0, ratio=0.95)
[13]:
plt.imshow(recon[0, :, :])
plt.show()
../_images/ipynb_tomopy_24_0.png
[ ]: