1b-Alignment_template_scattering¶
In the following template the scattering signal is used for the alignment of the data. For this, the normalization of the transmission signal needs to be done individually for each projection which is different compared to the scattering data.
Import packages & libraries
[1]:
%matplotlib ipympl
from tqdm import tqdm
import numpy as np
import os
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
import h5py
import pickle
from ipywidgets import interact
from data_processing.dataset import Dataset
from data_processing.filtered_back_projection import load_FBP_data_scattering, load_FBP_data
Step 1 - Load args from previous step.¶
In addition, we will adjust the normalization behaviour to account for e.g. fluctuating beam intensities, to be done projection by projection
[3]:
#--- User Input to load the exported args file --------------------#
#filepath_args = '/das/work/p18/p18877/Christian/work_dir_Christian/Au_rod_short/args_Au_rod_short.pickle'
#filepath_args = '/data/visitors/formax/20220566/2023031408/process/reconstructions/9425L_4w_XHP_NT/args_9425L_4w_XHP_NT.pickle' # Irene implant
filepath_args = '/data/visitors/formax/20220566/2023031408/process/work_dir_Christian/big_brain/args_big_brain.pickle'
#--- End Input to load the exported args file --------------------#
## ---------- Don't change this if not needed ----------------- ##
with open(filepath_args,'rb') as fid:
args = pickle.load(fid)
frame_id_range = args['frame_id_range']
for key in args.keys():
print(key)
base_path
work_directory
integration_folder
frameID
sample_name
air_id
proposal
visit
bkg_scan
detector_name
normkey
valid_pixels
fast_scan_direction
slow_scan_direction
projection_direction
detector_angle_origin
detector_angle_pos_dir
norm_transmission
flat_field_level
correct_background
tomodata
beamline
inner_rotation_axis
inner_rotation_key
first_rotation_indexes
tilt_axis
tilt_key
data_sorting
data_index_origin
principal_rotation_right_handed
secondary_rotation_right_handed
detector_angle_0
detector_angle_right_handed
offset_positive
air_transmission
snake
phi_det
q
norm_sum
air_scattering
multiprocessing
frame_id_range
q_rois
Load beamline specific loader functions based on args input
[4]:
#do not modify
print(f"Loading beamline specific functions from {args['beamline']}_utils")
# Import the cSAXS loader functions
if args['beamline'] == 'csaxs':
from data_processing.cSAXS_utils import metadata_reader_complete as metadata_reader
from data_processing.cSAXS_utils import transmission_loader_mcs as transmission_loader
from data_processing.cSAXS_utils import scattering_loader as scattering_loader_eiger
from data_processing.cSAXS_utils import create_args
# Import the PX loader functions
elif args['beamline'] == 'px':
from data_processing.PX_utils import metadata_reader_json as metadata_reader
from data_processing.PX_utils import transmission_loader_eiger as transmission_loader
from data_processing.PX_utils import scattering_loader_eiger
from data_processing.PX_utils import create_args
# Import the PX loader functions
elif args['beamline'] == 'formax':
from data_processing.ForMAX_utils import metadata_reader
from data_processing.ForMAX_utils import transmission_loader
from data_processing.ForMAX_utils import scattering_loader_eiger
from data_processing.ForMAX_utils import create_args
Loading beamline specific functions from formax_utils
Load the dataset
Loading of data, might take a couple of minutes depending on the allocation (recommended 8 cores or more).
GPUs are a bit slower here. With 8-12 cores it takes around 1min for full dataset
[5]:
%%time
## ---------- Don't change this if not needed ----------------- ##
ds = Dataset(frame_id_range, metadata_reader, transmission_loader, scattering_loader_eiger, **args)
CPU times: user 335 ms, sys: 365 ms, total: 700 ms
Wall time: 7.01 s
Step 2 - Visualize absorption data¶
[6]:
#------------------------------------------------- User Input for Plot ------------------------------------------------#
#Histogram bins and range
bins = 500
histrange = (-0.01, 0.1)
vmin = None # Default None
vmax = None # Default None
#------------------------------------------------- End Input for Plot ------------------------------------------------#
#------------------------------------------------- Plotting --------------------------------------------------#
fig, axs = plt.subplots(1,2, figsize = (10,5))
#update the colorbar for plot
cax = None
i = 0
axs[1].hist(ds.stack[i].absorption_coeff.flatten(), bins = bins)#, range = histrange)
@interact(i=(0, len(ds.stack) -1))
def plot(i):
axs[0].clear()
axs[1].clear()
# Modify vmin, vmax and set them to fixed values to avoid autoscaling
#if vmin is None and vmax is None:
#vmin, vmax = np.percentile(ds.stack[i].transmission_factor ,[10,90])
img = axs[0].imshow(ds.stack[i].absorption_coeff, cmap = 'Greys_r', vmin = vmin, vmax = vmax)
axs[0].set_title(f'absorption for projection #{i}')
#-----------------update colorbar-----------------#
global cax
if cax is None:
cax = fig.colorbar(img, ax=axs[0]).ax
else:
fig.colorbar(img, cax=cax)
#----------------------------------#
# histogram
axs[1].hist(ds.stack[i].absorption_coeff.flatten(), bins = bins, range=histrange)
axs[1].set_title(f'histogram absorbance for projection #{i}')
# Change ranges
axs[1].set_xlim(histrange[:])
axs[1].grid()
Step 3 - Make tomogram for 0 tilt (FBP projections)¶
The very first step for the tomography requires you to find the center of rotation for the sinograms.
In the future, this can be automated for sure. There are various concepts or algorithms for this. cSAXS codes for SASTT alignment and ptycho tomo alignment has a range of scripts there from which ideas could be taken and implemented here!
Note, the second plot shows the angle used for the sinogram. There is a field in
args['rotation_zero_tilt']from which these angles are taken now. In case this is off, this might has to be adapted manually. Note: This is already created in 0-Plotting_template, but can easily changed here by accessing the dictionary args and changing the input
Note: It laods the scattering data with multiple cores if available from allocation. It might take a few minutes, but will print out time statistics due to the line %%time
[7]:
%%time
#------------------------------------------------- User Input for qranges ------------------------------------------------#
q_range = None #[0.1, 0.15] # Choose a q range, or use None for full q range
# To allow multi processing of data
args['multiprocessing'] = 1
#------------------------------------------------- End Input for qranges ------------------------------------------------#
# if q_range is None:
# q_range = [args['q'][0], args['q'][-1]]
## ---------- Don't change this if not needed ----------------- ##
sinograms, tomo_angles = load_FBP_data_scattering(ds, q_range)
CPU times: user 8.91 s, sys: 19.5 s, total: 28.4 s
Wall time: 3min 23s
Step 4.0 - Plot sinograms, find center of rotation in sinograms for 0 tilt data¶
Adjust the position of the red line to match the center of the sinograms
[16]:
#--- User Input, find center of rotation --------------------#
# correct value until red line is in center of sinogram, this can for sure be improved with some automization!
off_center_value = 1
#--- End Input, find center of rotation --------------------#
# Plot sinograms from 0 tilt
image_shape = sinograms.shape[:2]
args['off_center_tomo_axis'] = off_center_value
fig, axs = plt.subplots(1,2,figsize = (10,5))
if args['inner_rotation_axis'][-1] == args['fast_scan_direction'][-1]:
print(f"Sample is rotating around the fast-scanning direction, '{args['fast_scan_direction'][-1]}'")
axs[0].set_title('Example sinogram')
axs[0].imshow(sinograms[:, image_shape[1]//2,:])
axs[0].plot([0 ,sinograms.shape[2]-1], [image_shape[0]//2 + off_center_value,
image_shape[0]//2 + off_center_value], color="red", linewidth=2)
axs[0].set_ylabel(f" Scan direction {args['fast_scan_direction']}")
axs[0].set_xlabel(args['inner_rotation_key'])
else:
print(f"Sample is rotating around the slow-scanning direction, {args['slow_scan_direction'][-1]}")
fig = plt.figure(figsize = (5,5))
axs[0].set_title('Example sinogram')
axs[0].imshow(sinograms[image_shape[0]//2,:,:])
axs[0].plot([0 ,sinograms.shape[2]-1], [image_shape[1]//2 + off_center_value,
image_shape[1]//2 + off_center_value], color="red", linewidth=2)
axs[0].set_ylabel(f" Scan direction {args['fast_scan_direction']}")
axs[0].set_xlabel(args['inner_rotation_key'])
# get angles from current selection
rot_angles, tilt_angles = ds.get_angles()
rot_0_range = slice(args['first_rotation_indexes'][0], args['first_rotation_indexes'][-1]+1)
# Plot angles for sinogram
axs[1].plot(tomo_angles, '-o', color='black', label = 'rotation')
ax2 = axs[1].twinx()
ax2.plot(tilt_angles[rot_0_range], '-o', color='red', label = 'tilt')
ax2.set_ylim([0,45])
ax2.set_ylabel('tilt angle')
axs[1].set_title('Rotation angles from list of first rotation indexes')
axs[1].set_xlabel('projection number')
axs[1].set_ylabel('tomo rotation angle')
axs[1].legend()
ax2.legend(loc='lower right')
plt.show()
Sample is rotating around the slow-scanning direction, y
FBP reconstruction at 0 tilt
No additional alignment
[17]:
#--- User Input, adjust color scale --------------------#
off_center_value = args['off_center_tomo_axis'] #It picks it up from the section above!
# Modify vmin, vmax to adjust the colorbar in the images below, otherwise it is autoscaling
vmin = None
vmax = None
#--- End Input, adjust color scale --------------------#
from data_processing.filtered_back_projection import FBP_recon
# Do FBP reconstruction
tomogram_FBP = FBP_recon(sinograms, tomo_angles, off_center_value = off_center_value, **args)
volume_shape = tomogram_FBP.shape
# Plot orthogonal central slices
fig, axs = plt.subplots(1,3,figsize = (10,6))
@interact(slice_xy=(0, volume_shape[2] -1))
def plot(slice_xy):
axs[0].imshow(tomogram_FBP[:,:,slice_xy], vmin = vmin, vmax = vmax)
axs[0].set_xlabel('y')
axs[0].set_ylabel('x')
axs[0].set_title(f"slice_xy, z={slice_xy}")
@interact(slice_xz=(0, volume_shape[1] -1))
def plot(slice_xz):
axs[1].imshow(tomogram_FBP[:,slice_xz,:], vmin = vmin, vmax = vmax)
axs[1].set_xlabel('x')
axs[1].set_ylabel('z')
axs[1].set_title(f"slice_xz, y={slice_xz}")
@interact(slice_yz=(0, volume_shape[0] -1))
def plot(slice_yz):
axs[2].imshow(tomogram_FBP[slice_yz,:,:], vmin = vmin, vmax = vmax)
axs[2].set_xlabel('z')
axs[2].set_ylabel('y')
axs[2].set_title(f"slice_yz, x={slice_yz}")
Create a valid mask
In this step you can create a mask for the reconstruction. This is useful to for instance mask out needles or some high absorption regions from within the sample
Mads wrote a function to automatically detects high absorption regions which can be used
Alternatively, one can manually exclude regions. In the example below, we had some issues at one tilt angle (crashing into a pipe while measuring), which is visible in the data. The few lines of code mask out the top 19 rows for this specific tilt angle (tile > 35)
[29]:
from data_processing.filtered_back_projection import make_mask_array
#------------------------------------------------- User Input for automatic masking ------------------------------------------------#
if 'mask_options' in args:
mask_options = args['mask_options']
else:
mask_options = {'threshold': np.inf, 'extendwidth': 0} # np.inf means no masking. the default value is factor 1000 absorption
args['mask_options'] = mask_options
#------------------------------------------------- End Input for automatic masking ------------------------------------------------#
# Mask creation based on mask_options
masks = make_mask_array(ds, index_list = range(len(ds.stack)), mask_options = mask_options)
#------------------------------------------------- User Input for manual masking ------------------------------------------------#
#Here you can manually adjust the mask for each projection
# Example 1
masks[:10,:,:] = 0 # This for instance masks the first 5 rows, uncomment for giving it a try
masks[:,0,:] = 0
masks[:,-1,:] = 0
masks[-1,:,:] = 0
# Example 2
#In this case for instance, I would like to exclude the top rows by adjusting the mask from the highest tilt
#Find relevant angles
#rot_angles, tilt_angles = ds.get_angles()
#index = [i for i, tilt in enumerate(tilt_angles) if tilt > 35]
# Exclude them by setting the mask to 0
#masks[:19,:,index] = 0
#------------------------------------------------- End Input for manual masking ------------------------------------------------#
from data_processing.filtered_back_projection import make_mask_array
fig, axs = plt.subplots(1,3, figsize = (10,4))
cax = None
@interact(proj_nr=(0, len(ds.stack) -1))
def plot(proj_nr):
axs[0].clear()
axs[1].clear()
img = axs[0].imshow(ds.stack[proj_nr].absorption_coeff, cmap = 'Greys_r')#,vmin =0 , vmax = 1)#sample dependent
axs[0].set_title(f'Absorption coeff projection#{proj_nr}')
global cax
if cax is None:
cax = fig.colorbar(img, ax=axs[0]).ax
else:
fig.colorbar(img, cax=cax)
img = axs[1].imshow(ds.stack[proj_nr].absorption_coeff*masks[:,:,proj_nr], cmap = 'Greys_r')
axs[1].set_title(f'Absorption coeff projection masked')
img = axs[2].imshow(masks[:,:,proj_nr], interpolation = None, vmin = 0, vmax=1)
axs[2].set_title(f'Mask')
419it [00:00, 4829.23it/s]
Step 4 - full dataset¶
Geometry parameters are being defined based on the args that were created. This may differ from beamline to beamline. It will be important to set those once properly for different measurements
This step might already be the first that only runs on a GPU
Note: This step may take some (fully integrated data needs to be loaded) - it might take time, a few minutes with multiple cores allocated (12minutes for 419 projections brain)
[30]:
%%time
## ---------- Don't change this section if not needed ----------------- ##
# Load all the absorption data (also non zero tilt)
sinograms_full, _ = load_FBP_data_scattering(ds, q_range, index_list = range(len(ds.stack)), mask_needle = False, mask_options = args['mask_options'])
# Make all the rotation matrices
from data_processing.filtered_back_projection import str_to_vec
from data_processing.geometry_parameters import _Rx, _Ry, _Rz
# Define a function that Generate the rotation matrix for each projection
def r(metadata):
# Read out rotation angles and convert to radians
alpha = metadata['roty']*np.pi/180
beta = metadata['rotx']*np.pi/180
# Construct and return the rotation matrix
return _Rx(beta) @ _Ry(alpha)
rot_mat = np.zeros((len(ds.stack), 3, 3))
for ii, proj in enumerate(ds.stack):
rot_mat[ii, :,:] = r(proj.metadata)
from data_processing.geometry_parameters import GeometryParameters
geom = GeometryParameters(rot_mat,
args['phi_det'],
volume_shape,
image_shape,
p_direction_0 = str_to_vec(args['projection_direction']),
j_direction_0 = str_to_vec(args['slow_scan_direction']),
k_direction_0 = str_to_vec(args['fast_scan_direction']),
detector_direction_origin = str_to_vec(args['detector_angle_origin']),
detector_direction_positive90 = str_to_vec(args['detector_angle_pos_dir']),
)
volume_shape = tomogram_FBP.shape
image_shape = sinograms_full.shape[:2]
from data_processing.projectors import ProjectorNumba
projector = ProjectorNumba(geom,
volume_shape,
image_shape,
mask = masks,
)
No offsets given
CPU times: user 1min 10s, sys: 2min 20s, total: 3min 31s
Wall time: 12min 23s
Back project tilt angles
Note: Up until now, there has still not been any alignment applied to individual projections
[31]:
## ---------- Don't change this section if not needed ----------------- ##
print(tomogram_FBP.shape)
FBP_forward_proj = projector.project_stack(np.ascontiguousarray(tomogram_FBP))
#TODO Fixing an error with the current implem. of the numba projector. Should be fixed in the projector iself.
if len(FBP_forward_proj.shape) == 4:
FBP_forward_proj = FBP_forward_proj.reshape(FBP_forward_proj.shape[:-1])
(87, 131, 87)
100%|██████████| 419/419 [00:06<00:00, 66.69it/s]
Check the output
In this step it is import to check and compare that back projected projections and raw data rotate and tilt in the same direction. If not, this might hint towards tomography related paramters in args that are note correctly set
The visualization below should help with this
[42]:
# -------- This section is mostly for visualization, one could for instance modify vmin, vmax to avoid autoscaling --------- #
# Initialize plot
vmin = None
vmax = None
fig, axs = plt.subplots(1,2, figsize = (12,6))
cax = None
cax2 = None
@interact(proj_nr=(0, len(ds.stack) -1))
def plot(proj_nr):
axs[0].clear()
print(f"roty = {ds.stack[proj_nr].metadata['roty']:.1f}, rotx = {ds.stack[proj_nr].metadata['rotx']:.1f}")
axs[0].set_title(f'Forward projected\nFBP recon. Proj #{proj_nr}')
img = axs[0].imshow(FBP_forward_proj[:,:,proj_nr], cmap = 'Greys_r',vmin=vmin, vmax=vmax)
axs[0].set_xlabel(args['fast_scan_direction'])
axs[0].set_ylabel(args['slow_scan_direction'])
global cax
if cax is None:
cax = fig.colorbar(img, ax=axs[0]).ax
else:
fig.colorbar(img, cax=cax)
axs[1].clear()
axs[1].set_title(f'Measured proj #{proj_nr}')
img2 = axs[1].imshow(sinograms_full[:,:,proj_nr], cmap = 'Greys_r', vmin=vmin, vmax=vmax)
axs[1].set_xlabel(args['fast_scan_direction'])
axs[1].set_ylabel(args['slow_scan_direction'])
global cax2
if cax2 is None:
cax2 = fig.colorbar(img2, ax=axs[1]).ax
else:
fig.colorbar(img2, cax=cax2)
Step 4.1 - Alignment of projections¶
At this moment, there are two possible ways implemented.
First one: You can choose the FOV by cropping the data in a first step, also to avoid for instance regions where the needle comes into the FOV and then align projection with subpixel precision.
Secondly, an option that uses the mask but can not do subpixel shifts. The cropped FOV could also be applied for the second option.
They produce two different outputs: (1) shift_1 and (2) shift_2. Afterwards you may choose with which to continue
Some more background on the procedure: Both alignments use the skimage.registration package and cross correlates the image by comparing forward projected data from FBP to the raw data. It computes the shifts to optimize the position for each projection (Mads will be able to explain this in much more detail..)
4.1a Subpixel alignment
The red rectangle shows you the region that you picked for the alignment
OUTPUT: shift_1
[44]:
sinograms_full.shape
[44]:
(131, 85, 419)
[48]:
crop_region = (slice(10, -10), slice(None), slice(None))
sinograms_full[crop_region].shape
print(crop_region[1].start, crop_region[0].start)
None 10
[49]:
#------------------------------------------------- User Input for chosing region for alignment ------------------------------------------------#
#Choose here values to crop the FOV for alignment based pixels in y direction (y being the scanning direction parallel to the tomography axis)
crop_region = (slice(10, -10), slice(None), slice(None)) # For no cropping slice(None, None) otherwise positive and negative the distance from top and bottom
#------------------------------------------------- End Input for chosing region for alignment ------------------------------------------------#
fig, axs = plt.subplots(1,1, figsize = (6,6))
height = sinograms_full[crop_region].shape[0]
width = sinograms_full[crop_region].shape[1]
@interact(proj_nr=(0, len(ds.stack) -1))
def plot(proj_nr):
axs.imshow(sinograms_full[:,:,proj_nr])
if crop_region[0].start ==None:
corner_right = 0
else:
corner_right = crop_region[0].start
if crop_region[1].start ==None:
corner_left = 0
else:
corner_left = crop_region[1].start
rect = Rectangle((corner_left, corner_right), width, height, fill = False, color = "red", linewidth = 2)
axs.add_patch(rect)
axs.set_xlabel('-x')
axs.set_ylabel('y')
Compute subpixel shift corrections
[50]:
from skimage.registration import phase_cross_correlation as ss_phase_cross_correlation
from data_processing.utils import normal_cross_corr
#init vector
shifts_1 = []
# Loop through images
for proj_nr in tqdm(range(len(ds.stack))):
#measured_image = np.pad(full_abs_data[:,:,proj_nr],pad)
measured_image = sinograms_full[crop_region][:,:,proj_nr]
#shifts_1.append(normal_cross_corr(FBP_forward_proj[:,:,proj_nr], measured_image))
shifts_1.append(ss_phase_cross_correlation(FBP_forward_proj[crop_region][:,:,proj_nr],
measured_image, upsample_factor = 100,)[0])
shifts_1 = np.array(shifts_1)
#Plot result
plt.figure()
plt.plot(shifts_1[:,0])
plt.plot(shifts_1[:,1])
plt.grid()
plt.legend(['shift index 1; dy', 'shift index 2; dx'])
plt.xlabel('Projection number')
plt.ylabel('Shift (pixels)')
plt.show()
100%|██████████| 419/419 [00:02<00:00, 175.90it/s]
4.1b Subpixel alignment
The red rectangle shows you the region that you picked for the alignment
OUTPUT: shift_2
[51]:
# sub-pixel shift correction
from skimage.registration import phase_cross_correlation as ss_phase_cross_correlation
from data_processing.utils import normal_cross_corr
from tqdm import tqdm
#init vector
shifts_2 = []
# Loop through images
for proj_nr in tqdm(range(len(ds.stack))):
#measured_image = np.pad(full_abs_data[:,:,proj_nr],pad)
measured_image = sinograms_full[:,:,proj_nr]
#shifts_1.append(normal_cross_corr(FBP_forward_proj[:,:,proj_nr], measured_image))
shifts_2.append(ss_phase_cross_correlation(FBP_forward_proj[:,:,proj_nr], measured_image, reference_mask = masks[:,:,proj_nr], moving_mask = masks[:,:,proj_nr] ))
if np.mean(masks[ :, :, proj_nr]) < 0.5:
shifts_2[-1] = shifts_2[-2]
if np.any(np.abs(shifts_2[-1])>15):
shifts_2[-1] = shifts_2[-2]
shifts_2 = np.array(shifts_2)
#Plot result
plt.figure()
plt.plot(shifts_2[:,0])
plt.plot(shifts_2[:,1])
plt.grid()
plt.legend(['shift index 1; dy', 'shift index 2; dx'])
plt.xlabel('Projection number')
plt.ylabel('Shift (pixels)')
plt.show()
100%|██████████| 419/419 [00:03<00:00, 114.59it/s]
Step 5 - SIRT reconstructions¶
This block for sure only runs with GPUs. Be sure that you have the right allocation. It reconstructs the full tomo using the SIRT algorithm. You can choose the number of iterations below.
Pick from one of the alignment option - 1.) shifts_1 - subpixel shifts (potentially cropped FOV to avoid needle in alignment) - 2.) shifts_2 - no subpixel shifts but with mask - Note: running this section may take some time
[52]:
#------------------------------------------------- User Input for applying shifts ------------------------------------------------#
# Here, one can choose between picking either shifts_1 or shifts_2 based on the 2 different alignemnt procedures
shift = shifts_1 #shifts_2
# This parameter defines the number of iterations for the SIRT alignment algorithm, typically 50 is a good value
n_iter_SIRT = 100 # 50 is a good default value
#------------------------------------------------- User Input for applying shifts ------------------------------------------------#
# Do SIRT tomographic reconstructions using the full dataset
from data_processing.projectors import ProjectorAstra
# recon with shifts
geom_2 = GeometryParameters(rot_mat,
args['phi_det'],
volume_shape,
image_shape,
p_direction_0 = str_to_vec(args['projection_direction']),
j_direction_0 = str_to_vec(args['slow_scan_direction']),
k_direction_0 = str_to_vec(args['fast_scan_direction']),
detector_direction_origin = str_to_vec(args['detector_angle_origin']),
detector_direction_positive90 = str_to_vec(args['detector_angle_pos_dir']),
offsets_j = shift[:,0],
offsets_k = shift[:,1],
)
projector_2 = ProjectorAstra(geom_2,
volume_shape,
image_shape,
mask = masks,
)
sirt_recon = projector_2.SIRT_recon(sinograms_full.transpose((1,0,2)),n_iter_SIRT,vmin = 0)
# recon without shifts
projector_1 = ProjectorAstra(geom,
volume_shape,
image_shape,
mask = masks,
)
sirt_recon_no_align = projector_1.SIRT_recon(sinograms_full.transpose((1,0,2)),n_iter_SIRT,vmin = 0)
Step 5.1 Check quality of SIRT reconstruction
The next step is important to check the quality of the reconstructions for the alignment. - The plot shows you the central slice through the tomogram, namely xy, yz and xz, for 3 different reconstructions: - 1. FBP reconstruction for the first rotation - 2. SIRT reconstruction of the full dataset with aligned projections - 3. SIRT reconstruction of the full dataset without alignment applied
The SIRT reconstructions with alignment should give the best results!
[53]:
#------------------------------------------------- User Input for plotting ------------------------------------------------#
# Adjust here which slice is shown, initial one here is the central slide
slice_shown = (volume_shape[0]//2, volume_shape[1]//2, volume_shape[2]//2) #(volume_shape[0]//2, volume_shape[1]//2, volume_shape[2]//2)
## Change [1,99] to adjust
vmin, vmax = np.percentile(tomogram_FBP[:,:,:], [1,99])# change [1,99] to adjust the colorscheme
#------------------------------------------------- End Input for plotting ------------------------------------------------#
print(volume_shape)
fig, axs = plt.subplots(3,3, figsize = (9,10))
@interact(proj_nr_zy=(0, tomogram_FBP.shape[0] -1))
def plot(proj_nr_zy):
axs[0,0].clear()
axs[1,0].clear()
axs[2,0].clear()
# FBP data
axs[0,0].imshow(tomogram_FBP[proj_nr_zy,:,:], vmin = vmin, vmax = vmax)
#axs[0,0].set_title('First rotation FBP')
axs[0,0].set_xlabel('y')
axs[0,0].set_ylabel('z')
# SIRT data
axs[1,0].imshow(sirt_recon[proj_nr_zy,:,:], vmin = vmin, vmax = vmax)
#axs[1,0].set_title('Full dataset SIRT aligned projections')
axs[1,0].set_xlabel('y')
axs[1,0].set_ylabel('z')
# SIRT data, no alignment of projections
axs[2,0].imshow(sirt_recon_no_align[proj_nr_zy,:,:], vmin = vmin, vmax = vmax)
#axs[2,0].set_title('Full dataset SIRT no alignment')
axs[2,0].set_xlabel('y')
axs[2,0].set_ylabel('z')
@interact(proj_nr_xy=(0, tomogram_FBP.shape[1] -1))
def plot(proj_nr_xy):
axs[0,1].clear()
axs[1,1].clear()
axs[2,1].clear()
# FBP data
axs[0,1].imshow(tomogram_FBP[:, proj_nr_xy, :], vmin = vmin, vmax = vmax)
#axs[0,1].set_title('First rotation FBP')
axs[0,1].set_xlabel('y')
axs[0,1].set_ylabel('x')
# SIRT data
axs[1,1].imshow(sirt_recon[:, proj_nr_xy, :], vmin = vmin, vmax = vmax)
#axs[1,1].set_title('Full dataset SIRT aligned projections')
axs[1,1].set_xlabel('y')
axs[1,1].set_ylabel('x')
# SIRT data, no alignment of projections
axs[2,1].imshow(sirt_recon_no_align[:, proj_nr_xy, :], vmin = vmin, vmax = vmax)
#axs[2,1].set_title('Full dataset SIRT no alignment')
axs[2,1].set_xlabel('y')
axs[2,1].set_ylabel('x')
@interact(proj_nr_xz=(0, tomogram_FBP.shape[2] -1))
def plot(proj_nr_xz):
axs[0,2].clear()
axs[1,2].clear()
axs[2,2].clear()
# FBP data
axs[0,2].imshow(tomogram_FBP[:, :, proj_nr_xz], vmin = vmin, vmax = vmax)
axs[0,2].set_title('First rotation FBP')
axs[0,2].set_xlabel('z')
axs[0,2].set_ylabel('x')
# SIRT data
axs[1,2].imshow(sirt_recon[:, :, proj_nr_xz], vmin = vmin, vmax = vmax)
axs[1,2].set_title('Full dataset SIRT aligned projections')
axs[1,2].set_xlabel('z')
axs[1,2].set_ylabel('x')
# SIRT data, no alignment of projections
axs[2,2].imshow(sirt_recon_no_align[:, :, proj_nr_xz], vmin = vmin, vmax = vmax)
axs[2,2].set_title('Full dataset SIRT no alignment')
axs[2,2].set_xlabel('z')
axs[2,2].set_ylabel('x')
plt.show()
# Add some comments explanations about what to look for here!
(87, 131, 87)
Step 6 - Save alignment parameters¶
[54]:
# Ouput the useful files
import json
import pickle
import os
# Make data directory
dir_name = args['work_directory'] + f"/{args['sample_name']:s}"
try:
os.mkdir(dir_name)
except FileExistsError:
print(f'Directory {dir_name:s} already exists.')
pass
# Save shift values
shift_directory = {'offsets_j':[s for s in shifts_1[:,0]],
'offsets_k':[s for s in shifts_1[:,1]]}
with open(dir_name + '/shift_values.json','w') as fid:
json.dump(shift_directory, fid)
# Save masks array
with open(dir_name + '/mask_arrays.npy','wb') as fid:
pickle.dump(masks, fid)
# Save absorption tomogram
with open(dir_name + '/tomogram.npy','wb') as fid:
pickle.dump(sirt_recon, fid)
Directory /data/visitors/formax/20220566/2023031408/process/work_dir_Christian/big_brain already exists.
[ ]: