Reconstruct Liver with TomoPy

Overview

This tutorial demonstrates 3d scalar field reconstruction using data derived from penetrating waves. 

Approximate time: ??? (I never completed this, though the instructions are complete and should work. I invite you to try it from beginning to end and update how long it took you!)

Dependencies

Input

Importing

Preprocessing

Reconstruction

Visualizing

Dependencies

FFmpeg, Anaconda, TomoPy, matplotlib, NumPy, cv2, and glob

Install FFmpeg, a command line tool that can be used for video processing and conversion.

Install Anaconda to manage package installations for TomoPy. It also necessitates running your script from a virtual environment.

In terminal, run the following to install the TomoPy packages into a new Conda environment :

conda create --name tomopy --channel conda-forge tomopy

Activating this virtual environment will make the installed packages accessible from your path. In terminal:

conda activate tomopy

Synchrotron facilities will often produce data in particular formats. Installing DXchange will read in these kind of files. In terminal:

conda install -c conda-forge dxchange

NumPy is used for optimized array processing, matplotlib is used for visualizations, glob is used to load in the directory of frames, and cv2 is used to read in the images.

Create a new file for your Python script. The first line should import all the dependencies:

import tomopy, dxchange, numpy, matplotlib, cv2, glob

Input

The dataset for this tutorial is an x-ray tilt series of a liver. Download liver.avi from Google Drive. 

The file you downloaded can be viewed in Google Drive or with VLC.

This data was produced by a sensor reading in the x-ray waves that successfully penetrated through the liver. These projections were taken at a variety of angles, equally spaced. The projections are bundled into a video, but we need individual frames, so the first step will be to compile the projections into a directory of images. FFmpeg can do this easily for a variety of formats. We must use a lossless image format such as PNG, GIF, or TIFF. In the same directory as the video.file, in terminal:

ffmpeg -i liver.avi frame_%03d.png

The %03d dictates that the ordinal number of each output image will be formatted using 3 digits.

Organize the produced frames into a folder named projections and store it in the same directory as your Python script.

Importing 

Now that we have a directory of projections, the next step is bringing them into Python in a format suitable for TomoPy. 

First, lets compile a list containing the paths to each frame:

image_files = glob.glob('projections/*')

Make a function to read in each path. Take a path as input, read it in as a cv2 image, then convert it into a NumPy array

def read_to_numpy(path):

    img = cv2.imread(path)

    return numpy.array(img)

Apply the function to each path in image_files. This can be done with map but requires appropriate retyping. Map applies a function to each obect in an iterable. It does this lazily by returning a map object which doesn't actually compute until iterated through. We can use list to iterate, but then make it a NumPy array again

proj = numpy.array(list(map(read_to_numpy, image_files)))[:,:,:,0]

We then make a list called theta containing the angles corresponding to each projection, each 1º apart

theta = tomopy.angles(proj.shape[0])

Preprocessing

Next, we linearize the data by taking the minus log of each pixel. This helps ensure stable processing

proj = tomopy.minus_log(proj)

Though the projections are the only required parameters for TomoPy’s reconstruction method, supplying the rotation center greatly increases reconstruction accuracy. To find a good approximation of the center use find center

rot_center = tomopy.find_center(proj, theta, init=290, ind=0, tol=0.5)

The output of this function represents the center as a proportion of the frame (e.g. if the image is 400 wide and the center is at 100, the center would be .25)

Reconstruction

The TomoPy recon algorithm returns the reconstructed object as a 3d NumPy array

recon = tomopy.recon(proj, theta, center=rot_center, algorithm='gridrec', sinogram_order=False)

Visualization

Let's see the results! Circle masking helps display contrast with matplotlib because coloring is based on relative differences

recon = tomopy.circ_mask(recon, axis=0, ratio=0.95)

matplotlib.plt.imshow(recon[0, :, :])

matplotlib.plt.show()