4-Investigate_recon_template | mumott v 0.2.1

• This template is currently still in a quite fragile state. It includes several approaches to look at different data set

• 1. Checking single reconstructed qbin from mumott, compare with raw dat. to does: - compute DOO, Orientation, and Mean Int, make export for paraview

• 2. Checking q resolved reconstructed data from mumott, compare with raw data: same as above, but for q resolved data

Packages and libraries

#---- Do not modify ----#
%matplotlib ipympl
import h5py
import os
import pickle
import json
from tqdm import tqdm

import numpy as np
from IPython.display import display
from ipywidgets import interact

import matplotlib.pyplot as plt
from matplotlib import cm, colors
from matplotlib.patches import Rectangle
from mpl_toolkits.axes_grid1 import make_axes_locatable

# Mumott related inports
from mumott.data_handling import DataContainer, TransformParameters
from mumott.output_handling import OutputHandler, LiveViewHandler, ProjectionViewer
from mumott.optimization import OptimizationParameters, Optimizer, Regularizer, RegularizationParameters

OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.

Case 1 - Only one 1 bin

Step 1 - Loading data

Load first args from file

#---- Modify input for your sample! ----#
filepath_args =  '/das/work/p18/p18877/Christian/work_dir_Christian/Au_rod_short/args_Au_rod_short.pickle'
#--- End Input to load the exported args file ---------------------#

import pickle
## ---------- Don't change this if not needed ----------------- ##
with open(filepath_args,'rb') as fid:
    args = pickle.load(fid)

#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 csaxs_utils

Fix path for data loading of Mumott

path = os.path.join(args['work_directory'] , args['sample_name'])

input_type = 'h5'
# ---- modify for your sample ---- #
name = f'dataset_q_0.047_0.066.{input_type}'
ouput_path = os.path.join(args['work_directory'] , args['sample_name'], 'output')

Step 2 - Input on shape of the data required

#Adjust shape of the data here
volume_shape = [53, 53, 116]
projection_shape = (len(args['scan_id_range']), 53, 116 ,8) # tuple with entries (number of projections, axis , axis, n_directions)

Load reconstructed and raw data

#### ---- Do not modify ---- ####
# Load all files in output_path
output_fname =[]
for file in os.listdir(ouput_path):
    if file.endswith('h5') and file.startswith('dataset_q_'):
        output_fname.append(file)


# Just looks for data in the output folder with the same dataset_q_ start.
# In case a specific scan should be loaded, it is recommended to give it a manual filename here
import re
output_data = [fname for fname in output_fname if re.match(name[:-3], fname)][0]

#### ---- If automatic reading fails or wants to be skipped ...
#### ---- User input if wanted ---- ####

# ouput_path = ''
# output_data = ''

#### ---- Do not modify ---- ####

# Create output dictionary
data = {}
recons = {}

# Init transform parameters
transform_parameters = TransformParameters(
    data_sorting = args['data_sorting'], #(1, 0, 2)
    data_index_origin = args['data_index_origin'], #(0, 0),
    principal_rotation_right_handed = args['principal_rotation_right_handed'], #True,
    secondary_rotation_right_handed = args['secondary_rotation_right_handed'], #True,
    detector_angle_0 = args['detector_angle_0'], #(1, 0),
    detector_angle_right_handed = args['detector_angle_right_handed'], #=False,
    offset_positive = args['offset_positive']) #(True, True))



with h5py.File(os.path.join(ouput_path, output_data),'r') as gh:
    data[f"recon"] = np.copy(gh['output_data/synthetic_data']).reshape(projection_shape)
    recons[f"recon_mean"] = np.copy(gh['output_data/mean_of_amplitude']).reshape(volume_shape)
    recons[f"recon_DOO"] = np.copy(gh['output_data/relative_std_of_amplitude']).reshape(volume_shape)
    recons[f"recon_EV"] = np.copy(gh['output_data/eigenvectors']).reshape((*volume_shape,3,3))
    try:
        data_container = DataContainer(data_path=path,
                                       data_filename=name,
                                       data_type=input_type)
    except FileNotFoundError:
        print('No data file found!')
    data_container.transform(transform_parameters)
    # Create correct volume_shape
    data_container.stack.volume_shape = np.array([data_container.stack[0].diode.shape[0],
                                              *data_container.stack[0].diode.shape])
    # Create reconstruction parameters
    reconstruction_parameters = data_container.reconstruction_parameters
    data[f"raw"] = reconstruction_parameters.data.reshape(projection_shape)

Step 3 Compare synthetic vs. raw data

fig, axs = plt.subplots(1,3,figsize=(10,4))

print(data['raw'].shape, data['recon'].shape)


@interact(proj_nr=(0,data['recon'].shape[0]-1))
def plot(proj_nr):

    vmin, vmax = np.nanpercentile(np.sum(data['raw'][proj_nr,...], axis=-1), [20,90])
    axs[0].imshow(np.sum(data['raw'][proj_nr,...], axis=-1), norm=colors.LogNorm(vmin = vmin, vmax = vmax))
    axs[0].set_title('raw projection')

    axs[1].imshow(np.sum(data['recon'][proj_nr,...], axis=-1), norm=colors.LogNorm(vmin = vmin, vmax = vmax))
    axs[1].set_title('synthetic projection, recon')

    img=axs[2].imshow(np.abs((np.sum(data['recon'][proj_nr,...],axis=-1)-
                            np.sum(data['raw'][proj_nr,...],axis=-1))/
                           np.sum(data['raw'][proj_nr,...],axis=-1)),vmin = 0,vmax = 0.2)
    axs[2].set_title('normalized difference abs((recon-raw)/raw)')

(289, 53, 116, 8) (289, 53, 116, 8)

Step 4 - Take a look at the data, mean intensity and DOO

# Adjust colorscal if wanted, None is default and will autoscale
vmin = None
vmax = None

threshold = 1.2 # in percentage

mask = recons[f"recon_mean"] > threshold*np.mean(recons[f"recon_mean"])
print(threshold*np.mean(recons[f"recon_mean"]))

print(recons[f"recon_mean"].shape, recons[f"recon_DOO"].shape)

cax = None
cax2 = None

# Plot orthogonal central slices
fig, axs = plt.subplots(2,3,figsize = (12,6))
@interact(slice_xy=(0, recons[f"recon_mean"].shape[2] -1))
def plot(slice_xy):
    img = axs[0,0].imshow(recons[f"recon_mean"][:,:,slice_xy], vmin = vmin, vmax = vmax)

    global cax
    if cax is None:
        divider = make_axes_locatable(axs[0,0])
        cax = divider.append_axes("right", size="5%", pad=0.05)
        cax = fig.colorbar(img, cax = cax, ax=axs[0,0]).ax
    else:
        fig.colorbar(img, cax=cax)

    axs[0,0].set_xlabel('y')
    axs[0,0].set_ylabel('x')
    axs[0,0].set_title(f"mean intensity slice z = {slice_xy}")
    vmin2, vmax2 = np.nanpercentile(recons[f"recon_DOO"][:,:,slice_xy][mask[:,:,slice_xy]],[0,95])
    img2 = axs[1,0].imshow(recons[f"recon_DOO"][:,:,slice_xy]*mask[:,:,slice_xy],
                           cmap=cm.inferno,
                           vmin = vmin2, vmax = vmax2)

    global cax2
    if cax2 is None:
        divider = make_axes_locatable(axs[1,0])
        cax2 = divider.append_axes("right", size="5%", pad=0.05)
        cax2 = fig.colorbar(img2, cax = cax2, ax=axs[1,0]).ax
    else:
        fig.colorbar(img2, cax=cax2)

    axs[1,0].set_title(f"'Masked Degree of orientation'")

@interact(slice_xz=(0, recons[f"recon_mean"].shape[1] -1))
def plot(slice_xz):
    axs[0,1].imshow(recons[f"recon_mean"][:,slice_xz,:], vmin = vmin, vmax = vmax)
    axs[0,1].set_xlabel('z')
    axs[0,1].set_ylabel('x')
    axs[0,1].set_title(f"slice y= {slice_xz}")
    vmin2, vmax2 = np.nanpercentile(recons[f"recon_DOO"][:,slice_xz,:][mask[:,slice_xz,:]],[0,95])
    axs[1,1].imshow(recons[f"recon_DOO"][:,slice_xz,:]*mask[:,slice_xz,:],
                    cmap=cm.inferno,
                    vmin = vmin2, vmax = vmax2)

@interact(slice_yz=(0, recons[f"recon_mean"].shape[0] -1))
def plot(slice_yz):
    axs[0,2].imshow(recons[f"recon_mean"][slice_yz,:,:], vmin = vmin, vmax = vmax)

    axs[0,2].set_xlabel('z')
    axs[0,2].set_ylabel('y')
    axs[0,2].set_title(f"slice x= {slice_yz}")
    vmin2, vmax2 = np.nanpercentile(recons[f"recon_DOO"][slice_yz,:,:][mask[slice_yz,:,:]],[0,95])
    print(vmin2,vmax2)
    axs[1,2].imshow(recons[f"recon_DOO"][slice_yz,:,:]*mask[slice_yz,:,:],
                    cmap=cm.inferno,
                    vmin = vmin2, vmax = vmax2)#, norm = colors.LogNorm(vmin = vmin2, vmax = vmax2))

88.4266977724216
(53, 53, 116) (53, 53, 116)

Case 2 - Q resolved data

Step 1 - Loading data

Load first args from file

#---- Modify input for your sample! ----#
filepath_args = '/data/visitors/formax/20220566/2023031408/process/reconstructions/implant_9425L_4w_XHP_NT/args_implant_9425L_4w_XHP_NT.pickle'
#--- End Input to load the exported args file ---------------------#

import pickle
## ---------- Don't change this if not needed ----------------- ##
with open(filepath_args,'rb') as fid:
    args = pickle.load(fid)

Fix path for data loading of Mumott

path = os.path.join(args['work_directory'] , args['sample_name'])

input_type = 'h5'
# ---- modify for your sample ---- #
name = f'dataset_q_0.047_0.066.{input_type}'
ouput_path = os.path.join(args['base_path'] , 'process/Software_repo_final/', 'output_q_res_recon')
#ouput_path = os.path.join(args['work_directory'] , args['sample_name'], 'output_q_res_recon')

Step 2 - Input on shape of data

#Adjust shape of the data here
volume_shape = [55, 55, 124] # Could be automated too, otherwise check in hdf5 output file shape of vectors easily
projection_shape = (len(args['scan_id_range']), 55, 124 ,8) # tuple with entries (number of projections, axis , axis, n_directions)

Load reconstructed and raw data

Note: The sorting of data into the dictionary here is not yet smart. It loads data as it finds it in the directory, such that the indizes of the loaded data need to be found later within the next cells

#### ---- Do not modify ---- ####
import re
# Load all files in output_path
output_files = os.listdir(ouput_path)
raw_files = os.listdir(path)

data = {}

input_type = 'h5'

# Init transform parameters
transform_parameters = TransformParameters(
    data_sorting = args['data_sorting'], #(1, 0, 2)
    data_index_origin = args['data_index_origin'], #(0, 0),
    principal_rotation_right_handed = args['principal_rotation_right_handed'], #True,
    secondary_rotation_right_handed = args['secondary_rotation_right_handed'], #True,
    detector_angle_0 = args['detector_angle_0'], #(1, 0),
    detector_angle_right_handed = args['detector_angle_right_handed'], #=False,
    offset_positive = args['offset_positive']) #(True, True))

print(f"Start loading of data, total number of files in folder is {len(output_files)}")
for ii, file in tqdm(enumerate(output_files)):
    if file.endswith('h5'):
        indices = [i for i, x in enumerate(raw_files) if x == file.replace('_output.h5','')]
        #print(ii, file.replace('_output.h5',''),files[indices[0]])
        with h5py.File(os.path.join(ouput_path,output_files[ii])) as gh:
            #print(list(gh['output_data']))
            data[f"recon_{ii}"] = np.copy(gh['output_data/synthetic_data']).reshape(projection_shape)
            data_container = DataContainer(data_path=path,
                                       data_filename=raw_files[indices[0]],
                                       data_type=input_type)
            data_container.transform(transform_parameters)
            # Create correct volume_shape
            data_container.stack.volume_shape = np.array([data_container.stack[0].diode.shape[0],
                                                      *data_container.stack[0].diode.shape])
            # Create reconstruction parameters
            reconstruction_parameters = data_container.reconstruction_parameters

            data[f"raw_{ii}"] = reconstruction_parameters.data.reshape(projection_shape)

Start loading of data, total number of files in folder is 65

65it [09:27,  8.73s/it]

Step 3 - Compare synthetic vs. raw data

fig, axs = plt.subplots(1,3,figsize=(15,4))

print(data['raw_1'].shape, data['recon_1'].shape)

# cax=None

q_range_list = []
q_range_list = [float(value.split('_')[2]) for value in output_files]
q_range_list.sort()
q_ranges = np.array(q_range_list)

#@interact(proj_nr=(0,data['recon_0'].shape[0]-1))
@interact(qbin=(0,len(q_ranges)-1), proj_nr=(0,data['recon_0'].shape[0]-1))
def plot(proj_nr, qbin):
    filename = f'dataset_q_{q_ranges[qbin]:.3f}'
    index = [ii for ii, x in enumerate(output_files) if x.__contains__(filename)][0]
    vmin, vmax = np.nanpercentile(np.sum(data[f"raw_{index}"][proj_nr,...], axis=-1), [20,90])
    axs[0].imshow(np.sum(data[f"raw_{index}"][proj_nr,...], axis=-1),vmin = vmin, vmax = vmax)
    axs[0].set_title(f"raw projection, qbin {qbin}")

    axs[1].imshow(np.sum(data[f"recon_{index}"][proj_nr,...], axis=-1), vmin = vmin, vmax = vmax)
    axs[1].set_title(f"synthetic projection (recon), qbin {qbin}")

    img=axs[2].imshow(np.abs((np.sum(data[f"recon_{index}"][proj_nr,...],axis=-1)-
                            np.sum(data[f"raw_{index}"][proj_nr,...],axis=-1))/
                           np.sum(data[f"raw_{index}"][proj_nr,...],axis=-1)),vmin = 0,vmax = 0.2)
    # global cax
    # if cax is None:
    #     cax = fig.colorbar(img, ax=axs[2], aspect=2).ax
    # else:
    #     fig.colorbar(img, cax=cax, aspect=2)
    axs[2].set_title(f"normalized difference abs((recon-raw)/raw)")

(281, 55, 124, 8) (281, 55, 124, 8)

Step 4 - Compare 1D scattering from recons and raw

def compute_q_scattering(data, proj, xcor, ycor, q_range_list):
    scattering_raw = np.zeros(q_range_list.shape)
    scattering_recon = np.zeros(q_range_list.shape)
    qvals = np.zeros(q_range_list.shape)
    for ii, q in enumerate(q_range_list):
        q_range = q_range_list[ii:ii+2]
        #filename = f'dataset_q_{q_range[0]:.3f}_{q_range[1]:.3f}'
        filename = f'dataset_q_{q_range[0]:.3f}'
        index = [jj for jj, x in enumerate(output_files) if x.__contains__(filename)]
        index = index[0]
        scattering_recon[ii] = np.nansum(data[f"recon_{index}"][proj, ycor, xcor,:])
        scattering_raw[ii] = np.nansum(data[f"raw_{index}"][proj, ycor, xcor,:])
        qvals[ii] = q
    return qvals, scattering_recon, scattering_raw

# Create list of q based on output
q_range_list = []
q_range_list = np.array([float(value.split('_')[2]) for value in output_files].sort())
#q_range_list.sort()

qbin = 38
proj = 0

q_range_list = []
q_range_list = [float(value.split('_')[2]) for value in output_files]
q_range_list.sort()
q_ranges = np.array(q_range_list)

q_range = q_ranges[qbin]

filename = f'dataset_q_{q_range:.3f}'
index = [ii for ii, x in enumerate(output_files) if x.__contains__(filename)][0]

################################ plotting
fig, axs = plt.subplots(2,2,figsize=(10,8))

img = np.sum(data[f"recon_{index}"][proj,:,:,:],axis=2)
vmin, vmax = np.percentile(img,[45,95])

img_ref = np.sum(data[f"raw_{index}"][proj,:,:,:],axis=2)

axs[0,0].imshow(img, picker=True, vmin  = vmin, vmax = vmax)
axs[0,0].set_title(f"Reconstruction sim. proj {proj:3d}: q range {q_range:.2f} Ang-1")
axs[1,0].imshow(img_ref, vmin  = vmin, vmax = vmax)
axs[1,0].set_title(f"Raw projection {proj:3d}: q range {q_range:.2f} Ang-1")

ycor = 39
xcor = 28
def onpick(event):
    mouseevent = event.mouseevent
    xcor = int(mouseevent.xdata)
    ycor = int(mouseevent.ydata)
    axs[0,1].clear()
    axs[1,0].clear()
    axs[1,1].clear()
    axs[0,0].clear()
    #axs[1,1].plot(data[f"recon_{index}"][proj,ycor,xcor,:],'o-',color='black',label='recon')
    #axs[1,1].plot(data[f"raw_{index}"][proj,ycor,xcor,:],'o-',color='cyan',label='raw')
    #axs[1,1].legend()


    axs[0,0].imshow(img, picker=True, vmin  = vmin, vmax = vmax)
    axs[1,0].imshow(img_ref, vmin  = vmin, vmax = vmax)
    axs[0,0].plot(xcor, ycor, '+', ms = 8, color='r')
    axs[1,0].plot(xcor, ycor, '+', ms = 8, color='r')
    qvals, scat_1D_recon, scat_1D_raw = compute_q_scattering(data,proj,xcor,ycor,q_ranges)
    axs[0,1].loglog(qvals,scat_1D_recon,'-o',color='black',label='recon')
    axs[0,1].loglog(qvals,scat_1D_raw,'-o',color='cyan',label='raw')
    axs[0,1].legend()
    axs[1,1].plot(qvals,scat_1D_raw*qvals**2,'-o',color='cyan')
    axs[1,1].plot(qvals,scat_1D_recon*qvals**2,'-o',color='black')
    axs[1,1].set_xlim(0,0.4)
    axs[1,1].set_ylabel('Iq2')

fig.canvas.mpl_connect('pick_event',onpick)

9

to do: 3D view clicking tool for 1D scattering, export to paraview, 3D reciprocal map for voxel