In situ time series alignment class

[ceobrero]

What this PR does...

Note: updated from PR #119. TimeSeries class is now object-oriented and serialized. Pytests included.

Takes an unaligned time series image stack and applies 2D cross-correlation to correct for drift. The image stack is preprocessed (padded and edge-blended) upon instantiation. The subpixel shifts are determined on an image stack that is filtered via 2D convolution to accentuate features (edges) to align with.

Example of applying 2D cross-correlation on a preprocessed time series image stack:

binned_pad.align_stack(running_average_frames = 10, sigma_edge = 0.7, sf_val = 0.5)

Visualizing the image stack mean before and after alignment:

Example Notebook:
dev_time_series_example.ipynb

[ceobrero]

@bobleesj updated this time series alignment PR, ready for review!

[bobleesj]

Thanks @ceobrero

General comment on performance since for times-series, performance is the primary concern (just like in DriftCorrection).

• We want to do FFT at once across all images with vectorized operations instead of applying FFT across each loop.
  sth like:

 # FFT entire stack at once, multiply, IFFT back
      F = np.fft.fft2(images, axes=(-2, -1))
      F_filtered = F * edge_kernel  # Broadcasting: (N,H,W) * (H,W)
      filtered = np.abs(np.fft.ifft2(F_filtered, axes=(-2, -1)))

• For filtering, etc. we want to utilize fourier if possible. Right now I see there are 6 convolutions are being called but we can do this in a single pass.

Sth like:

      kx = np.fft.fftfreq(H)[:, None]
      ky = np.fft.fftfreq(W)[None, :]
      gaussian = np.exp(-2 * (np.pi * sigma_edge)**2 * (kx**2 + ky**2))
      gradient_mag = np.sqrt((2 * np.pi * kx)**2 + (2 * np.pi * ky)**2)
      edge_kernel = gaussian * gradient_mag

• Since it's time-series and it must be done sequentially w/ moving average, we probalby want to have some sort of progress:

from tqdm import tqdm
  for a0 in tqdm(range(nframes - 1), desc="Aligning"):

• Parabolic fitting vs.. DFT upsampling (already in cross_correlation_shift), which one do we choose?

Also, if you could share, a performance chart would be good. So that we know how long it would take process a batch of images across size dimensions.

[bobleesj]

sf_val is a bit hard to guess what this parameter means without looking at the API. the intenral variable names can be also improved.

[bobleesj]

also, if we end-up using a simple variable name, we want to write as compcat as possible like writing a series of equation, etc without uncessasry blank lines (please use the below as example only)

def _gaussian_1d_kernel(sigma: float, device, dtype):
    # Match your r range: [-ceil(2*sigma), ..., ceil(2*sigma)]
    radius = int(torch.ceil(torch.tensor(2.0 * sigma)).item())
    r = torch.arange(-radius, radius + 1, device=device, dtype=dtype)
    k = torch.exp(-(r * r) / (2.0 * sigma * sigma))
    k = k / k.sum()  # normalize (optional but usually desirable)
    return k  # shape (K,)