[FEATURE] Align different reconstructions with cross-correlation

[ewbellec]

Recently, I had to align different scans of Pd particles during the alpha-beta transition. The particle was getting eaten little by little so the center of mass (COM) is always different for each scans and you can't realign the particles between each other using only the COM. I used a cross-correlation instead.
Mouad had the same issue of aligning different Bragg reconstructions of a Pt particle with several dislocations. For each Bragg, the dislocation holes are at different position and thus the COM is always slightly different.
I made a function to realign them based on the support. It looks a bit complicated but it's very simple :

• make a list of objects to align with each other (make sure the objects have all the same voxel size)
• force the same array shape on all objects (if needed)
• find the shift regarding the reference object using skimage.registration.phase_cross_correlation
• shift the objects using scipy.ndimage.shift (no sub-pixel shifts here)
  my codes :

def create_support(obj, threshold_module,
                  fill_support=False) :
    
    module = np.abs(obj)
    support = (module > threshold_module * np.max(module))
    
    if fill_support :
        support_convex = np.zeros(support.shape)
        for axis in range(support.ndim):
            support_cum = np.cumsum(support, axis=axis)
            support_cum_inv = np.flip(np.cumsum(np.flip(support,axis=axis), axis=axis), axis=axis)
            support_combine = support_cum * support_cum_inv
            support_convex[support_combine != 0] = 1
        return support_convex
    else:
        return support

def force_same_shape(obj_list,
                     verbose = False) :
    
    '''
    In case objects in obj_list don't have the same shape, 
    this function forces an identical shape by padding objects with 0s
    :obj_list: a list of objects of shape (number of objects, individual object shape)
    '''
    
    shape_list = [obj.shape for obj in obj_list]
    
    if np.all([shape == shape_list[0] for shape in shape_list]) :
        if verbose :
            print('All objects already have the same shape')
        return obj_list
    
    forced_shape = np.max(shape_list, axis=0)
    
    for n, obj in enumerate(obj_list):
        pad = (np.array(forced_shape) - np.array(obj.shape))
        padding = [(p//2, p//2 + p%2) for p in pad]
        obj_list[n] = np.pad(obj, padding,mode='constant',constant_values=(0,))
        
    if verbose :
        print(f'All objects have now the shape : {forced_shape}')
        
    return np.array(obj_list)

from skimage.registration import phase_cross_correlation
import scipy
def realign_object_list(obj_list,
                        ref_index=0,
                        threshold_module=.15, fill_support=False) :
    
    '''
    Align all objects in obj_list using the supports and a phase_cross_correlation.
    Limited to integer pixels shift. No sub-pixel shifts.
    :obj_list: a list of objects of shape (number of objects, individual object shape)
    :ref_index: index of the reference object 
                (for example with ref_index=0, the first object is the reference position)
    :threshold_module: threshold (between 0 and 1) used to create the support.
    :fill_support: If True, fill holes in the support (mostly for particles with dislocations)
    '''
    
    obj_ref = np.copy(obj_list[ref_index])
    support_ref = create_support(obj_ref, threshold_module, fill_support=fill_support)
        
    obj_list_shift = np.zeros(obj_list.shape)
    for n, obj in enumerate(obj_list):
        support = create_support(obj, threshold_module, fill_support=fill_support)
        shift, error, diffphase = phase_cross_correlation(support_ref, support)
        obj_list_shift[n] += scipy.ndimage.shift(np.abs(obj), shift)
        
    return obj_list_shift

Then simply use the code

obj_list = [obj1, obj2] # or even more objects
obj_list = force_same_shape(obj_list, verbose=True)
obj_list_shift = realign_object_list(obj_list)

[ewbellec]

And here's the results on Mouad's particle.
With only center of mass alignement (first Bragg reconstruction modulus in Blue, the second in Red) :

With the cross-correlation :

Clearly much better.
For my particles in alpha/beta transition it's even more important since part of the particle is eaten in beta.

[clatlan]

Hello Ewen,
Thank you for this. My question is: is this a feature that should be accessible when running a standard analysis ? Or is this something that should be available in the cdiutils API? The way I see it is that once you have processed your scans, you make a new notebook to analyse and interpret the processed data and eventually make your publication figures. At this stage you need to align your reconstructions.

If it is only a matter of aligning reconstructions after the traditional process, it's easy!

[ewbellec]

It's only useful for later analysis, not really during the reconstruction process.
Maybe somewhere if you include multi-Bragg in cdi-utils but I wouldn't put this into the reconstruction for single Bragg. It's not a necessity and might lead to problems.
I also wanted to give the function for you since it's quite useful for particles transitioning partially to beta.
Maybe don't put it in cdi-utils for now and you can take back this code whenever you think it will be useful.

[clatlan]

It can be included in cdiutils.utils. You can even do it yourself :)