SIRT reconstruction with GPUΒΆ

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
# -------------------------------------------------------------------
# This script performs a complete reconstruction workflow.
# The reconstruction algorithm used is the SIRT performed on a GPU.
# -------------------------------------------------------------------
import numpy as np
import neutompy as ntp

# set pixel size in cm
pixel_size  = 0.0029

# set the last angle value of the CT scan: np.pi or 2*np.pi
last_angle = 2*np.pi

# read dataset containg projection, dark-field, flat-field images and the projection at 180 degree
proj, dark, flat, proj_180 = ntp.read_dataset()

# normalize the projections to dark-field, flat-field images and neutron dose
norm, norm_180 = ntp.normalize_proj(proj, dark, flat, proj_180=proj_180, dose_draw=True, crop_draw=True)

# rotation axis tilt correction
norm = ntp.correction_COR(norm, norm[0], norm_180)

# clean up memory
del dark; del flat; del proj; del proj_180

# remove outliers, set the optimal radius and threshold
norm = ntp.remove_outliers_stack(norm, radius=1, threshold=0.018, outliers='dark', out=norm)
norm = ntp.remove_outliers_stack(norm, radius=3, threshold=0.018, outliers='bright', out=norm)

# perform minus-log transform
norm =  ntp.log_transform(norm, out=norm)

# remove stripes in sinograms
norm = ntp.remove_stripe_stack(norm, level=4, wname='db30', sigma=1.5, out=norm)

# define the array of the angle views in radians
angles = np.linspace(0, last_angle, norm.shape[0], endpoint=False)

# SIRT reconstruction with 100 iterations using GPU
print('> Reconstruction...')
rec    = ntp.reconstruct(norm, angles, 'SIRT_CUDA', parameters={"iterations":100}, pixel_size=pixel_size)

# select the directory and the prefix file name of the reconstructed images to save.
recon_dir = ntp.save_filename_gui('', message = 'Select the folder and the prefix name for the reconstructed images...')

# write the reconstructed images to disk
ntp.write_tiff_stack(recon_dir, rec)