BEATS CT reconstruction pipeline - Paganin phase retrieval

Minimal TomoPy reconstruction pipeline for phase contrast data collected at the BEATS beamline of SESAME.

Created on: 23.05.2021
Last update: 09.10.2023

• By: Gianluca Iori, 2023

• Code license: MIT

• Narrative license: CC-BY-NC-SA

Type Ctrl + Enter on a single cell to run it.

Load experiment data

Enter the sample_name and the correct output_dir

sample_name = "bee_yazeed-20231001T170032"
work_dir = "/mnt/PETRA/SED/BEATS/IH/"+sample_name
h5file = work_dir+"/"+sample_name+".h5"

output_dir = "/mnt/PETRA/SED/BEATS/IH/scratch/tmp/"
recon_dir = output_dir+sample_name+"/recon/"
cor_dir = output_dir+sample_name+"/cor/"

Load the complete dataset (or)

# projs, flats, darks, theta = dxchange.read_aps_32id(h5file, exchange_rank=0)

Read a portion of the dataset

• sino controls the vertical detector lines to read - sino=(1360, 2160, 1)

• proj defines the range of projections - proj=(1, 1001, 1)

sino=(1360, 2160, 1)
proj=(1, 1001, 1)

projs, flats, darks, _ = dxchange.read_aps_32id(h5file, exchange_rank=0, sino=sino, proj=proj)
theta = np.radians(dxchange.read_hdf5(h5file, 'exchange/theta', slc=(proj,)))

print("Dataset size: ", projs[:, :, :].shape[:], " - dtype: ", projs.dtype)
print("Flat fields size: ", flats[:, :, :].shape[:])
print("Dark fields size: ", darks[:, :, :].shape[:])
print("Theta array size: ", theta.shape[:])

INFO:dxchange.reader:Data successfully imported: /mnt/PETRA/SED/BEATS/IH/bee_yazeed-20231001T170032/bee_yazeed-20231001T170032.h5
INFO:dxchange.reader:Data successfully imported: /mnt/PETRA/SED/BEATS/IH/bee_yazeed-20231001T170032/bee_yazeed-20231001T170032.h5
INFO:dxchange.reader:Data successfully imported: /mnt/PETRA/SED/BEATS/IH/bee_yazeed-20231001T170032/bee_yazeed-20231001T170032.h5
INFO:dxchange.reader:Data successfully imported: /mnt/PETRA/SED/BEATS/IH/bee_yazeed-20231001T170032/bee_yazeed-20231001T170032.h5
INFO:dxchange.reader:Data successfully imported: /mnt/PETRA/SED/BEATS/IH/bee_yazeed-20231001T170032/bee_yazeed-20231001T170032.h5

Dataset size:  (1000, 800, 2560)  - dtype:  uint16
Flat fields size:  (102, 800, 2560)
Dark fields size:  (102, 800, 2560)
Theta array size:  (1000,)

At any time you can take a look at your 3D data with ru.plotmidplanes(data)

ru.plot_midplanes(projs)

Flat field correction

Normalize the image background.

projs = tomopy.normalize(projs, flats, darks, ncore=ncore, averaging='median')
ru.plot_midplanes(projs)

Phase retrieval

delta_beta = 250 # ratio between real and imaginary part of the refractive index
alpha=1./(4*3.141592**2 * delta_beta)
print("alpha: ", alpha)

alpha:  0.00010132122580089835

projs_phase = tomopy.retrieve_phase(projs[:, :, :],
                                    pixel_size=1e-4*0.65,
                                    dist=2,
                                    energy=10,
                                    alpha=alpha,
                                    ncore=ncore,
                                    pad=True)
# nchunk=None,

ru.plot_midplanes(projs_phase)

Center Of Rotation (COR)

Save images reconstructed with COR range

cor_range = [1260, 1300, 1]

tomopy.write_center(projs, theta, cor_dir, cor_range)

INFO:tomopy.recon.algorithm:Reconstructing 40 slice groups with 36 master threads...

View them in Fiji

# os.system(Fiji_exe_stack + cor_dir+'{:04.2f}'.format(COR[0])+'.tiff &')

Manually insert the best COR

COR = 1301

Automatic detect COR

# COR = tomopy.find_center(projs, theta, init=projs.shape[2]/2, tol=0.5) # ind=200, 
# print("Automatic detected COR: ", COR, " - tomopy.find_center")

COR = tomopy.find_center_vo(projs, ncore=ncore)
print("Automatic detected COR: ", COR, " - tomopy.find_center_vo")

Automatic detected COR:  1301.5  - tomopy.find_center_vo

Reconstruction

GPU reconstruction with the ASTRA toolbox

Algorithm	fbp CUDA ASTRA

options = {'proj_type': 'cuda', 'method': 'FBP_CUDA'}

recon = tomopy.recon(projs_phase,
                     theta,
                     center=COR,
                     algorithm=tomopy.astra,
                     options=options,
                     ncore=1)

INFO:tomopy.recon.algorithm:Reconstructing 1 slice groups with 1 master threads...

Post-processing

Apply circular mask

recon = tomopy.circ_mask(recon,
                         axis=0,
                         ratio=0.8,
                         ncore=ncore)

ru.plot_midplanes(recon)

Write reconstructed dataset

Write output tiff stack as float32

fileout = recon_dir+'slice.tiff'
dxchange.writer.write_tiff_stack(recon, fname=fileout, axis=0, digit=4, start=0, overwrite=True)

Open virtual stack in ImageJ

# os.system(Fiji_exe_stack + fileout + ' &')