Stacking EELS spectrum images to create a time/temperature series

[mseifner]

Hi,

I recently tried to stack my EELS spectrum images acquired at different temperatures to create a temperature series. My goal is to create a 4D dataset where the 2D images of my sample can be aligned and the regions cropped to a common region of interest. I have not found a solution yet. Do you have a suggestion of how to approach this goal?

Best,
Michael

[ericpre]

I would try to do the following:

1. Make a stack of your spectrum images using hs.stack or hs.load (use a globe pattern to load several files and they should be stacked automatically)
2. Calculate the shifts using estimate_shift - you may need to sum on the signal axis and use as_signal2D to set the axes to use to estimate the shifts as function of temperature.
3. Transpose the original 4 dimensional signal (similar as previous step but without taking the sum) and apply the shifts from previous step. You may need to add another dimension to the shifts array to match the additional navigation axis
4. Transpose back to the correct signal dimension.

[CSSFrancis]

Hi Micheal,

Let me clean up and post some of the code we worked on!

In the meantime https://github.com/pyxem/pyxem/blob/main/pyxem/signals/insitu_diffraction2d.py#L212-L236

This is the relevant lines in pyxem (for doing this in 5D STEM)

@ericpre I can make this as an example maybe as I think this is something that might become more of a common desire

[ericpre]

An example would be good indeed.

I am wondering if it would be useful to add some argument to estimate_shift and/or align2D to facilitate this by selecting or ignoring some axes. It this is generic enough, do you think it would possible to simplify / replace https://github.com/pyxem/pyxem/blob/main/pyxem/signals/insitu_diffraction2d.py#L212-L236 with estimate_shift and/or align2D. Are they caveats that wouldn't make this possible?

[CSSFrancis]

I just used a simulated dataset of a cross in real space with 4 Dimensions. This should be equivilent to an EELS dataset with 4 Dimensions.

The alignment will bring everything into line and seems to work quite well. The only problem is that it requires that you load each individual time series into memory which (can be) a bit of an issue if you are limited by RAM. In this case there are more complicated ways to do this.

import numpy as np
import hyperspy.api as hs
from scipy import ndimage as ndi

signals = []
for i in range(10):
    s = np.zeros(( 60,60,128))
    s[30+i,:]=1
    s[:, 30+i]=1
    signals.append(hs.signals.Signal1D(s))

time_series = hs.stack(signals) 
virtual_images = time_series.sum(axis=(-1)).transpose(signal_axes=(0,1))   

shifts = virtual_images.estimate_shift2D()
shifts = hs.signals.Signal1D(shifts)
transposed = time_series.transpose(3)

def align_time_series(data, shift, order=1):
    data_t = ndi.affine_transform(
        data, np.identity(3), offset=(shift[0], shift[1], 0), order=order
    )
    return data_t


shifted = transposed.map(align_time_series, shift =shifts , inplace=False)
aligned = shifted.transpose(1)

[mseifner]

Hi,

Thank you very much to the both of you for your input! I tested the code on my dataset and it worked well for a limited amount of spectrum images. As addressed before, the RAM limitation became a problem after adding more and more spectrum images to the stack. You mentioned a more complicated way to perform the stacking in such a scenario. Could you please elaborate on the basic idea behind the complicated approach?

Best,
Michael

[CSSFrancis]

Thank you very much to the both of you for your input! I tested the code on my dataset and it worked well for a limited amount of spectrum images. As addressed before, the RAM limitation became a problem after adding more and more spectrum images to the stack.

So this method should work assuming you can fit one full time step of the eels data into memory.

I think you want to load all of the data lazily. Then stack the data which should return a lazy signal. Then save everything probably as a .zspy. That's just a nice step to make sure you save it efficiently. One thing to note is to make sure that the chunking is done along the time axis with a chunk size of 1. I think after stacking this should be the default. The rechunk function can also be used.

Then reload the data with lazy=True. And run the functions we discussed above.

You mentioned a more complicated way to perform the stacking in such a scenario. Could you please elaborate on the basic idea behind the complicated approach?

There are even less memory intensive ways of doing this but that requires that you spend some time saving and resaving the data using the rechunker package. In that case you have to chunk in the energy, and time axes and then resample the spanned x,y.

Let me know if you run into issues. I can make a more detailed explanation with code examples.

[mseifner]

Thanks a lot! I will try the suggested approach and let you know in case anything is unclear.