Notebook

PaganinProcessor examples¶

The code in this notebook contains some examples and comparisons of using the PaganinProcessor phase retrieval methods in CIL

This notebook requires CIL v24.1.0 or greater, check the version below

In [ ]:

import cil
print(cil.__version__)

Load some dependencies from CIL

In [ ]:

from cil.utilities import dataexample
from cil.processors import PaganinProcessor, Slicer, Binner, TransmissionAbsorptionConverter, Padder, RingRemover
from cil.utilities.display import show2D
from cil.recon import FDK, FBP
from cil.utilities.jupyter import islicer
from cil.io.utilities import HDF5_utilities
from cil.framework import AcquisitionGeometry, AcquisitionData

This notebook also requires TomoPy, numpy and matplotlib

In [ ]:

from tomopy.prep.phase import retrieve_phase
import numpy as np
import matplotlib.pyplot as plt 

Parallel beam data¶

First we test the PaganinProcessor with parallel beam data. In the following cell

Get a test parallel beam dataset
Perform a filtered back projection reconstruction using FBP()

In [ ]:

data = dataexample.SIMULATED_PARALLEL_BEAM_DATA.get()
data.reorder(order='tigre')
data.geometry.config.units = 'um'
data_abs = -1*data.log()
ig = data.geometry.get_ImageGeometry()
fbp =  FBP(data_abs, ig)
recon = fbp.run(verbose=0)

Next we repeat the above steps with the PaganinProcessor

Choose the parameters to be used in the phase retrieval

Delta and beta are the real and complex part of the material refractive index. These can be found for x-ray wavelengths at https://henke.lbl.gov/optical_constants/getdb2.html. Here we just use some exagerated values to demonstrate the effect.
The experiment peak energy in default units eV

In [ ]:

delta = 1
beta = 0.002
energy = 40000

Run the phase retrieval using the PaganinProcessor which is implemented based on Paganin 2002. The processor returns the material retrieved thickness, removing the effect of phase in the image

$$ T(x,y) = - \frac{1}{\mu}\ln\left (\mathcal{F}^{-1}\left (\frac{\mathcal{F}\left ( M^2I_{norm}(x, y,z = \Delta) \right )}{1 + \alpha\left ( k_x^2 + k_y^2 \right )} \right )\right ) $$

where

$T$, is the sample thickness,
$\mu = \frac{4\pi\beta}{\lambda}$ is the material linear

attenuation coefficient where $\beta$ is the complex part of the material refractive index and $\lambda=\frac{hc}{E}$ is the probe wavelength,

$M$ is the magnification at the detector,
$I_{norm}$ is the input image which is expected to be the

normalised transmission data,

$\Delta$ is the propagation distance,
$\alpha = \frac{\Delta\delta}{\mu}$ is a parameter determining

the strength of the filter to be applied in Fourier space where $\delta` is the real part of the deviation of the material refractive index from 1

$k_x, k_y = \left ( \frac{2\pi p}{N_xW}, \frac{2\pi q}{N_yW}

\right )$ where $p$ and $q$ are co-ordinates in a Fourier mesh in the range $-N_x/2$ to $N_x/2$ for an image with size $N_x, N_y$ and pixel size $W$.

In the following cell:

Run the PaganinProcessor to retrieve $T(x,y)$ from the test dataset
Reconstruct using FBP

In [ ]:

processor = PaganinProcessor(delta=delta, beta=beta, energy=energy)
processor.set_input(data)
thickness = processor.get_output(override_geometry={'propagation_distance':10})
fbp =  FBP(thickness, ig)
recon_thickness = fbp.run(verbose=0)

# calculate mu to get recon_attenuation with the same scaling as the original image
attenuation = thickness*processor.mu
fbp =  FBP(attenuation, ig)
recon_attenuation = fbp.run(verbose=0)

Next we test the phase retrieval using PaganinProcessor(full_retrieval=False). In this implementation, the same filter is applied in Fourier space but the $-log()$ is not applied. $$ I_{filt} = \mathcal{F}^{-1}\left (\frac{\mathcal{F}\left ( I(x, y,z = \Delta) \right )} {1 - \alpha\left ( k_x^2 + k_y^2 \right )} \right ) $$ This gives flexibility to apply a Paganin-like filter but doesn't require data that has already been converted from transmission to absorption.

In [ ]:

# Run PaganinProcessor as a filter using full_retrieval=False on the absorption data and reconstruct
processor = PaganinProcessor(delta=delta, beta=beta, energy=energy, full_retrieval=False)
processor.set_input(data_abs)
filtered_image = processor.get_output(override_geometry={'propagation_distance':10})
fbp =  FBP(filtered_image, ig)
recon_filter = fbp.run(verbose=0)

For comparison run Tomopy phase retreival with raw data, then convert to absorption and reconstruct

In [ ]:

tomopy_alpha = (1/(delta/beta))/(4*np.pi**2)
data_tomopy = data.copy()
data_tmp = retrieve_phase(data.array, pixel_size=processor.pixel_size, dist=processor.propagation_distance, energy=energy/1000, alpha=tomopy_alpha)
data_tomopy.fill(data_tmp)
data_tomopy = -1*data_tomopy.log()
ig = data_tomopy.geometry.get_ImageGeometry()
fbp =  FBP(data_tomopy, ig)
recon_tomopy = fbp.run(verbose=0)

Also run Tomopy phase retreival with absorption data and reconstruct

In [ ]:

tomopy_alpha = (1/(delta/beta))/(4*np.pi**2)
data_tomopy_abs = data_abs.copy()
data_tmp = retrieve_phase(data_abs.array, pixel_size=processor.pixel_size, dist=processor.propagation_distance, energy=energy/1000, alpha=tomopy_alpha)
data_tomopy_abs.fill(data_tmp)
ig = data_tomopy_abs.geometry.get_ImageGeometry()
fbp =  FBP(data_tomopy_abs, ig)
recon_tomopy_abs = fbp.run(verbose=0)

Compare the reconstructions

In [ ]:

vertical_slice = 67
show2D([recon, recon_thickness, recon_attenuation, recon_filter, recon_tomopy, recon_tomopy_abs],
        title=['Original image', 'Phase retrieval - thickness', 'Phase retrieval - scaled by mu', 
               'Phase retrieval - full_retrieval=False', 'Tomopy phase retrieval transmission', 'Tomopy phase retrieval absorption'], 
        axis_labels = ('horizontal_y', 'horizontal_x'), num_cols=3, slice_list=('vertical',vertical_slice))

Zoom in on the reconstructions

In [ ]:

vertical_slice = 67
x_range = slice(50,90)
y_range = slice(50,90)

show2D([recon.array[vertical_slice,x_range,y_range], recon_thickness.array[vertical_slice,x_range,y_range], recon_attenuation.array[vertical_slice,x_range,y_range], recon_filter.array[vertical_slice,x_range,y_range], recon_tomopy.array[vertical_slice,x_range,y_range], recon_tomopy_abs.array[vertical_slice,x_range,y_range]],
title=['Original image', 'Phase retrieval - thickness', 'Phase retrieval - scaled by mu', 'Phase retrieval - full_retrieval=False', 'Tomopy phase retrieval transmission', 'Tomopy phase retrieval absorption'], 
axis_labels = ('horizontal_y', 'horizontal_x'), num_cols=3)

Compare a cross-section through the reconstruction

In [ ]:

fig, axs = plt.subplots(1,2, figsize=(12,5))
ax = axs[0]
vertical_slice = 67
y_slice = 70
x_range = range(50,90)
ax.plot(x_range, recon.array[vertical_slice,y_slice,x_range])
ax.plot(x_range, recon_attenuation.array[vertical_slice,y_slice,x_range])
ax.plot(x_range, recon_filter.array[vertical_slice,y_slice,x_range])
ax.plot(x_range, recon_tomopy.array[vertical_slice,y_slice,x_range],'--')
ax.plot(x_range, recon_tomopy_abs.array[vertical_slice,y_slice,x_range],'--')

ax.set_xlabel('Horizontal x')
ax.set_ylabel('Intensity')
ax.set_title('Line Profile at horizontal_y=' + str(y_slice) + ', vertical slice=' + str(vertical_slice))

ax = axs[1]
x_slice = 70
y_range = range(50,90)
ax.plot(y_range, recon.array[vertical_slice,y_range,x_slice])
ax.plot(y_range, recon_attenuation.array[vertical_slice,y_range,x_slice])
ax.plot(y_range, recon_filter.array[vertical_slice,y_range,x_slice])
ax.plot(y_range, recon_tomopy.array[vertical_slice,y_range,x_slice],'--')
ax.plot(y_range, recon_tomopy_abs.array[vertical_slice,y_range,x_slice],'--')

ax.set_xlabel('Horizontal y')
ax.set_ylabel('Intensity')
ax.set_title('Line Profile at horizontal_x=' + str(x_slice) + ', vertical slice=' + str(vertical_slice))
ax.legend(['Original', 'Phase retrieved - scaled by mu', 'Phase retrieved - full_retrieval=False', 'TomoPy on transmission data', 'TomoPy on absorption data'])

We see that all methods blur the result in comparison to the original reconstruction. The scaled phase retrieval in CIL matches the Tomopy method performed on transmission data and the filter in CIL matches the Tomopy method performed on absorption data

We can approximate the signal to noise of each reconstruction as the mean divided by the standard deviation

In [ ]:

print("Original reconstruction SNR = " + str(recon.mean()/recon.array.std()))
print("Phase retrieved reconstruction SNR = " + str(recon_attenuation.mean()/recon_attenuation.array.std()))
print("Phase retrieved (full_retrieval=False) reconstruction SNR = " + str(recon_filter.mean()/recon_filter.array.std()))
print("TomoPy on transmission data reconstruction SNR = " + str(recon_tomopy.mean()/recon_tomopy.array.std()))
print("TomoPy on absorption data reconstruction SNR = " + str(recon_tomopy_abs.mean()/recon_tomopy_abs.array.std()))

In all cases, the phase retrieval improves the SNR

Cone beam data¶

In [ ]:

delta = 1
beta = 0.0001
energy = 40000

With cone beam data, the magnification $M$ has an effect on the phase retrieval
$ T = -\frac{1}{\mu}\ln(F^{-1}\frac{F(M^2 I_{norm}(x,y,z=\Delta))}{1+\frac{\Delta\lambda\delta}{4\pi\beta}(k_x^2+k_y^2)/M})$
The $M^2$ on top means sometimes we get a number larger than 1 inside the $\ln$

Get some cone beam data and perform reconstruction without phase retrieval

In [ ]:

data = dataexample.SIMULATED_CONE_BEAM_DATA.get()
data.geometry.config.units = 'um'
print('Magnification = ' + str(data.geometry.magnification))
data.reorder(order='tigre')
data_abs = -1*data.log()
ig = data.geometry.get_ImageGeometry()
fdk =  FDK(data_abs, ig)
recon = fdk.run(verbose=0)

Run phase retrieval on raw data and reconstruct

In [ ]:

processor = PaganinProcessor(delta=delta, beta=beta, energy=energy)
processor.set_input(data)
thickness = processor.get_output()
recon_thickness = fdk.run(verbose=0)
# calculate mu to get recon_attenuation with the same scaling as the original image
attenuation = thickness*processor.mu
fdk =  FDK(attenuation, ig)
recon_attenuation = fdk.run(verbose=0)

Run PaganinProcessor as a filter using get_output(full_retrieval=False) on the absorption data and reconstruct

In [ ]:

processor = PaganinProcessor(delta=delta, beta=beta, energy=energy, full_retrieval=False)
processor.set_input(data_abs)
filtered_image = processor.get_output()

fdk =  FDK(filtered_image, ig)
recon_filter = fdk.run(verbose=0)

TomoPy can only be used with parallel beam data so we do not use it for comparison here

Compare the reconstructions

In [ ]:

vertical_slice = 67
x_range = slice(50,90)
y_range = slice(50,90)
show2D([recon.array[vertical_slice,x_range,y_range], recon_thickness.array[vertical_slice,x_range,y_range], recon_attenuation.array[vertical_slice,x_range,y_range], recon_filter.array[vertical_slice,x_range,y_range]],
title=['Original image', 'Phase retrieval - thickness', 'Phase retrieval - scaled by mu', 'Phase retrieval - full_retrieval=False'],
axis_labels = ('horizontal_y', 'horizontal_x'))

Compare the cross-sections through the reconstructions

In [ ]:

fig, axs = plt.subplots(1,2, figsize=(12,5))
ax = axs[0]
vertical_slice = 67
y_slice = 70
x_range = range(50,90)
ax.plot(x_range, recon.array[vertical_slice, y_slice, x_range])
ax.plot(x_range, recon_attenuation.array[vertical_slice, y_slice, x_range])
ax.plot(x_range, recon_filter.array[vertical_slice, y_slice, x_range])

ax.set_xlabel('Horizontal x')
ax.set_ylabel('Intensity')
ax.set_title('Line Profile at horizontal_y=' + str(y_slice) + ', vertical slice=' + str(vertical_slice))

ax = axs[1]
x_slice = 70
y_range = range(50,90)
ax.plot(y_range, recon.array[vertical_slice,y_range,x_slice])
ax.plot(y_range, recon_attenuation.array[vertical_slice,y_range,x_slice])
ax.plot(y_range, recon_filter.array[vertical_slice,y_range,x_slice])

ax.set_xlabel('Horizontal y')
ax.set_ylabel('Intensity')
ax.set_title('Line Profile at horizontal_x=' + str(x_slice) + ', vertical slice=' + str(vertical_slice))
ax.legend(['Original', 'Phase retrieval - scaled by mu', 'Phase retrieval - full_retrieval=False'])

plt.tight_layout()

We can see that both methods blur the original image. The phase retrieval method becomes negative because of the values > 1 are passed into the negative log. This may be an indication that the full phase retrieval is not valid for this experimental setup, in which case using full_retrieval = False may be more useful as it just applies a Paganin-like filter to the data

Compare the signal to noise ratio

In [ ]:

print("Original reconstruction SNR = " + str(recon.mean()/recon.array.std()))
print("Phase retrieved reconstruction SNR = " + str(recon_attenuation.mean()/recon_attenuation.array.std()))
print("Phase retrieved (full_retrieval=False) reconstruction SNR = " + str(recon_filter.mean()/recon_filter.array.std()))

The negative values in the phase retrival skew the result but with full_retrieval=False the SNR is improved

Generalised Paganin method¶

The generalised Paganin method is implemented in CIL following the description in https://iopscience.iop.org/article/10.1088/2040-8986/abbab9
When features in the image are close to the Nyquist frequency of the system, a more generalised form of the Pagnin filter can be used which preserves these high frequency features while still boosting SNR. This may have a similar effect to applying an unsharp mask after the normal Paganin phase retrieval.

$$ T(x,y) = -\frac{1}{\mu}\ln\left (\mathcal{F}^{-1}\left (\frac{ \mathcal{F}\left ( M^2I_{norm}(x, y,z = \Delta) \right )}{1 - \frac{2 \alpha}{W^2}\left ( \cos(Wk_x) + \cos(Wk_y) -2 \right )} \right ) \right ) $$

Choose delta and beta

In [ ]:

delta = 1
beta = 0.001

Get the simulated parallel data and perform the reconstruction without phase retrieval

In [ ]:

data = dataexample.SIMULATED_PARALLEL_BEAM_DATA.get()
data.geometry.config.units = 'um'
data.reorder(order='tigre')
data_abs = -1*data.log()
ig = data.geometry.get_ImageGeometry()
fbp =  FBP(data_abs, ig)
recon = fbp.run(verbose=0)

Run phase retrival using the original Paganin method and reconstruct

In [ ]:

processor = PaganinProcessor(delta=delta, beta=beta, energy=energy)
processor.set_input(data)
thickness = processor.get_output(override_geometry={'propagation_distance':10})
fbp =  FBP(thickness, ig)
recon_thickness = fbp.run(verbose=0)

# calculate mu to get recon_attenuation with the same scaling as the original image
attenuation = thickness*processor.mu
fbp =  FBP(attenuation, ig)
recon_attenuation = fbp.run(verbose=0)

Run phase retrieval on the data using the generalised Paganin method and reconstruct

In [ ]:

processor = PaganinProcessor(delta=delta, beta=beta, energy=energy, filter_type='generalised_paganin_method')
processor.set_input(data)
thickness_GPM = processor.get_output(override_geometry={'propagation_distance':10})
fbp =  FBP(thickness_GPM, ig)
recon_thickness_GPM = fbp.run(verbose=0)

# calculate mu to get recon_attenuation with the same scaling as the original image
attenuation_GPM = thickness_GPM*processor.mu
fbp =  FBP(attenuation_GPM, ig)
recon_attenuation_GPM = fbp.run(verbose=0)

Plot cross-sections through the results

In [ ]:

fig, axs = plt.subplots(1,2, figsize=(12,5))
ax = axs[0]
vertical_slice = 67
y_slice = 70
x_range = range(50,90)
ax.plot(recon.array[vertical_slice, y_slice, x_range])
ax.plot(recon_attenuation.array[vertical_slice, y_slice, x_range])
ax.plot(recon_attenuation_GPM.array[vertical_slice, y_slice, x_range])

ax.set_xlabel('Horizontal x')
ax.set_ylabel('Intensity')
ax.set_title('Line Profile at horizontal_x=' + str(x_slice) + ', vertical slice=' + str(vertical_slice))
ax.legend(['Original', 'Phase retrieval', 'Phase retrieval - GPM'])

ax = axs[1]
x_slice = 70
y_range = range(50,90)
ax.plot(recon.array[vertical_slice,y_range,x_slice])
ax.plot(recon_attenuation.array[vertical_slice,y_range,x_slice])
ax.plot(recon_attenuation_GPM.array[vertical_slice,y_range,x_slice])

ax.set_xlabel('Horizontal y')
ax.set_ylabel('Intensity')
ax.set_title('Line Profile at horizontal_x=' + str(x_slice) + ', vertical slice=' + str(vertical_slice))

Check the SNR

In [ ]:

print("Original reconstruction SNR = " + str(recon.mean()/recon.array.std()))
print("Phase retrieved reconstruction SNR = " + str(recon_attenuation.mean()/recon_attenuation.array.std()))
print("Phase retrieved with GPM reconstruction SNR = " + str(recon_attenuation_GPM.mean()/recon_attenuation_GPM.array.std()))

The GPM has slightly improved resolution of the sample features while maintaining the SNR boost

TomoBank example¶

This example uses dataset tomo_00068 from TomoBank https://tomobank.readthedocs.io/en/latest/source/data/docs.data.phasecontrast.html#wet-sample which can be retrieved using: wget https://g-a0400.fd635.8443.data.globus.org/tomo_00068/tomo_00068.h5

The data were collected at Syrmep beamline of the Elettra synchotron. A description of the experiment is given in https://link.springer.com/chapter/10.1007/978-3-319-19387-8_70

Load the file, you may need to change the filename to the path where you downloaded it

In [ ]:

filename = 'tomo_00068.h5' 
data = HDF5_utilities.read(filename=filename, dset_path='/exchange/data')

Construct a CIL AcquisitionData object using the parameters in https://tomobank.readthedocs.io/en/latest/source/data/docs.data.phasecontrast.html#wet-sample

In [ ]:

pixel_size = 0.0041 #mm
propagation_distance = 150 #mm
angles = HDF5_utilities.read(filename=filename, dset_path='/exchange/theta')
ag = AcquisitionGeometry.create_Parallel3D(detector_position=[0, propagation_distance, 0], units='mm').set_panel([np.shape(data)[2],np.shape(data)[1]], pixel_size=pixel_size).set_angles(angles)
data = AcquisitionData(data, deep_copy=False, geometry = ag)
data.reorder(order='tigre')

Apply a centre of rotation correction

In [ ]:

data.geometry.set_centre_of_rotation(1463.5-(data.shape[2]/2), distance_units='pixels')

Use islicer to view all slices so we can choose a region to crop

In [ ]:

islicer(data)

Crop the data

In [ ]:

processor = Slicer(roi={'horizontal':(600,2500,1)})
processor.set_input(data)
data = processor.get_output()

Check the cropped data looks sensible

In [ ]:

show2D(data)

Next we bin a few angles

In [ ]:

print('Data shape before: ' + str(data.shape))
processor = Binner(roi={'angle':(None, None, 3)})
processor.set_input(data)
data = processor.get_output()
print('Data shape after: ' + str(data.shape))

Run phase retrieval on the raw data

Parameters from https://tomobank.readthedocs.io/en/latest/source/data/docs.data.phasecontrast.html#wet-sample

In [ ]:

delta = 1
beta = 1e-1
energy = 14000

processor = PaganinProcessor(delta=delta, beta=beta, energy=energy)
processor.set_input(data)
thickness = processor.get_output()

# calculate mu to get recon_attenuation with the same scaling as the original image
data_phase = thickness*processor.mu

Run phase retrieval using the generalised Paganin method on the raw data

In [ ]:

processor = PaganinProcessor(delta=delta, beta=beta, energy=energy, 
                             filter_type='generalised_paganin_method')
processor.set_input(data)
thickness = processor.get_output()

# calculate mu to get recon_attenuation with the same scaling as the original image
data_phase_generalised = thickness*processor.mu

Get a slice of each data set

In [ ]:

vertical_slice = 900
data_phase = data_phase.get_slice(vertical=vertical_slice)

In [ ]:

vertical_slice = 900
data_phase_generalised = data_phase_generalised.get_slice(vertical=vertical_slice)

For comparison just run TransmissionAbsorptionConverter on the same slice of the original data

In [ ]:

data_slice = data.get_slice(vertical=vertical_slice)

processor = TransmissionAbsorptionConverter()
processor.set_input(data_slice)
processor.get_output(out=data_slice)

Perform some extra processing steps on both datasets

In [ ]:

# Pad the data
ig = data_slice.geometry.get_ImageGeometry()
padsize = 1000
data_slice = Padder.linear_ramp(pad_width={'horizontal': padsize}, end_values=0)(data_slice)
data_phase = Padder.linear_ramp(pad_width={'horizontal': padsize}, end_values=0)(data_phase)
data_phase_generalised = Padder.linear_ramp(pad_width={'horizontal': padsize}, end_values=0)(data_phase_generalised)

# Ring remover
N_decompositions = 20
wavelet_filter_name = 'db20'
sigma = 5.5

processor = RingRemover(N_decompositions, wavelet_filter_name, sigma, info=False) 
processor.set_input(data_slice)
data_slice = processor.get_output()

processor = RingRemover(N_decompositions, wavelet_filter_name, sigma, info=False) 
processor.set_input(data_phase)
data_phase = processor.get_output()

processor = RingRemover(N_decompositions, wavelet_filter_name, sigma, info=False) 
processor.set_input(data_phase_generalised)
data_phase_generalised = processor.get_output()

# Reconstruct
fbp =  FBP(data_slice, ig)
recon = fbp.run(verbose=0)

fbp =  FBP(data_phase, ig)
recon_phase = fbp.run(verbose=0)

fbp =  FBP(data_phase_generalised, ig)
recon_phase_g = fbp.run(verbose=0)

Compare the reconstructions

In [ ]:

show2D([recon, recon_phase, recon_phase_g], ['Original','Paganin Method', 'Generalised Paganin Method'], num_cols=3)

Plot a difference map between the Paganin Method and Generalised Paganin Method reconstructions

In [ ]:

show2D(recon_phase.array[700:950,700:950]-recon_phase_g.array[700:950,700:950], title='Paganin Method - Generalised Paganin Method', cmap='seismic')

We see some of the high resolution details (e.g. the edges) are preserved in the GPM

Compare SNR

In [ ]:

print("Original reconstruction SNR = " + str(recon.mean()/recon.array.std()))
print("Phase retrieved reconstruction SNR = " + str(recon_phase.mean()/recon_phase.array.std()))
print("Phase retrieved with GPM reconstruction SNR = " + str(recon_phase_g.mean()/recon_phase_g.array.std()))