In [6]:

# -*- coding: utf-8 -*-
#  Copyright 2021 - 2022 United Kingdom Research and Innovation
#  Copyright 2021 - 2022 The University of Manchester
#  Copyright 2021 - 2022 Technical University of Denmark
#
#  Licensed under the Apache License, Version 2.0 (the "License");
#  you may not use this file except in compliance with the License.
#  You may obtain a copy of the License at
#
#      http://www.apache.org/licenses/LICENSE-2.0
#
#  Unless required by applicable law or agreed to in writing, software
#  distributed under the License is distributed on an "AS IS" BASIS,
#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#  See the License for the specific language governing permissions and
#  limitations under the License.
#
#   Authored by:    Jakob S. Jørgensen (DTU)
#                   Edoardo Pasca (UKRI-STFC)
#                   Laura Murgatroyd (UKRI-STFC)
#                   Gemma Fardell (UKRI-STFC)

Introductory demo of Core Imaging Library (CIL)¶

3D laboratory micro-CT, cone-beam data of a walnut¶

This exercise walks through the steps needed to load in and reconstruct by FDK a 3D cone-beam dataset of a walnut, acquired by laboratory micro-CT.

Learning objectives are:

Load and investigate a TXRM data set.
Apply CIL's TransmissionAbsorptionConverter.
Compute FDK reconstruction using CIL and compare with a reconstruction made from fewer projections.
Write out the reconstructed image to a TIFF stack.

This example requires the dataset walnut.zip from https://zenodo.org/record/4822516 :

https://zenodo.org/record/4822516/files/walnut.zip

If running locally please download the data and update the 'path' variable below.

In [7]:

path = '/mnt/materials/SIRF/Fully3D/CIL/Walnut'

In [8]:

# remove some annoying warnings
import logging
logger = logging.getLogger('dxchange')
logger.setLevel(logging.ERROR)

First import all of the modules we will need:

In [9]:

import os
from cil.io import ZEISSDataReader, TIFFWriter
from cil.processors import TransmissionAbsorptionConverter, Slicer
from cil.recon import FDK
from cil.utilities.display import show2D, show_geometry
from cil.utilities.jupyter import islicer, link_islicer

Load the 3D cone-beam projection data of a walnut:

In [10]:

filename = os.path.join(path, "valnut_2014-03-21_643_28/tomo-A/valnut_tomo-A.txrm")

data = ZEISSDataReader(file_name=filename).read()

The data is loaded in as a CIL AcquisitionData object:

In [11]:

type(data)

We can call print for the data to get some basic information:

In [12]:

print(data)

Note how labels refer to the different dimensions. We infer that this data set contains 1601 projections each size 1024x1024 pixels.

In addition to the data itself, AcquisitionData contains geometric metadata in an AcquisitionGeometry object in the geometry field, which can be printed for more detailed information:

In [13]:

print(data.geometry)

CIL can illustrate the scan setup visually from the AcquisitionData geometry:

In [14]:

show_geometry(data.geometry);

We can use the dimension labels to extract and display 2D slices of data, such as a single projection:

In [15]:

show2D(data, slice_list=('angle',800));

From the background value of 1.0 we infer that the data is transmission data (it is known to be already centered and flat field corrected) so we just need to convert to absorption/apply the negative logarithm, which can be done using a CIL processor, which will handle small/large outliers:

In [16]:

data = TransmissionAbsorptionConverter()(data)

We can also use a basic interactive viewer which will allow us to scroll through the various projections. islicer has a very similar interface to show2D, but not quite the same.

In [17]:

islicer(data)

We again take a look at a slice of the data, now a vertical one to see the central slice sinogram after negative logarithm:

In [18]:

show2D(data, slice_list=('vertical', 512));

CIL supports different back-ends for which data order conventions may differ. Here we use the FDK algorithm from CIL's recon module. FDK is filtered back-projection for cone beam data. By default, the recon module uses TIGRE as a back-end, and requires us to permute the data array into the right order:

In [19]:

data.reorder(order='tigre')

The data is now ready for reconstruction. To set up the FDK algorithm we must specify the size/geometry of the reconstruction volume. Here we use the default one:

In [20]:

ig = data.geometry.get_ImageGeometry()

We can then create the FDK algorithm and reconstruct the data:

In [21]:

fdk =  FDK(data, ig)
recon = fdk.run()

In [22]:

show2D(recon, slice_list=[('vertical',512), ('horizontal_x', 512)], fix_range=(-0.01, 0.06));

We can also interact with the data with islicer

In [23]:

islicer(recon, direction='vertical', minmax=(-0.01, 0.06))

In [24]:

islicer(recon, direction='horizontal_x', minmax=(-0.01, 0.06))

We can save the reconstructed volume to disk for example as a stack of TIFFs:

In [26]:

# save_base_path = os.getcwd()
# save_path = os.path.join(save_base_path, 'walnut')
# os.makedirs(save_path)

# TIFFWriter(data=recon, file_name=os.path.join(save_path, "out"), compression='uint16').write()

We now demonstrate the effect of reducing the number of projections on the FDK reconstruction.

In [ ]:

from cil.processors import Slicer

reduce_factor = 10

data_reduced = Slicer(roi={'angle': (0,-1,reduce_factor)})(data)

ig = data_reduced.geometry.get_ImageGeometry()
fdk =  FDK(data_reduced, ig)
recon_reduced = fdk.run()

We show the same slices as before:

In [ ]:

show2D(recon_reduced, slice_list=[('vertical',512), ('horizontal_x', 512)], fix_range=(-0.01,0.06));

We could also visually check the difference between the 2 reconstructions, interactively with islicer or compare side-by-side the reconstructions with show2D.

In [ ]:

sl1 = islicer(recon, minmax=(-0.01, 0.06))
sl2 = islicer(recon_reduced, minmax=(-0.01, 0.06))
link_islicer(sl1, sl2)

In [ ]:

show2D([recon, recon_reduced, recon, recon_reduced], 
       title=['Full data', 'Reduced data', 'Full data', 'Reduced data'], 
       slice_list=[('vertical',512), ('vertical',512), ('horizontal_x', 512), ('horizontal_x', 512)], 
       fix_range=(-0.01,0.06));

In the vertical slice of the data, a ring can be seen. This is an edge artifact common in FBP reconstruction. We can remove it from our image using a circular mask.

In [ ]:

recon.apply_circular_mask(radius=0.8, in_place=True)

show2D(recon);

Alternatively, the ring can be prevented by padding the AcquisitionData before reconstruction, to ensure that all data is within the field of view. We'll learn about padding in 02_intro_sandstone_parallel_roi