[WIP] match pixel size in sliding window on GPU

[LorenzLamm]

Rescale tomograms using Fourier Cropping / Padding in a sliding window approach. This can be done efficiently on GPU.

Joint work with @JoelVO

[codecov-commenter]

⚠️ Please install the to ensure uploads and comments are reliably processed by Codecov.

Codecov Report

Attention: Patch coverage is 0% with 58 lines in your changes missing coverage. Please review.

Project coverage is 0.00%. Comparing base (bea5ae0) to head (65e2e09).
Report is 10 commits behind head on main.

Files	Patch %	Lines
...processing/pixel_size_matching/match_pixel_size.py	0.00%	58 Missing ⚠️

❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files

@@          Coverage Diff           @@
##            main     #76    +/-   ##
======================================
  Coverage   0.00%   0.00%            
======================================
  Files         40      46     +6     
  Lines       1411    1670   +259     
======================================
- Misses      1411    1670   +259     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[alisterburt]

Hey guys - just FYI no pressure to use but torch-fourier-rescale is on PyPI for 2D/3D Fourier cropping/padding

https://github.com/alisterburt/torch-fourier-rescale/tree/main/src/torch_fourier_rescale

[alisterburt]

The plan is to move that and other image processing utilities

• torch-image-lerp
• torch-grid-utils
• torch-fourier-slice
• torch-subpixel-crop
• torch-fourier-shell-correlation
• torch-fourier-shift
  into teamtomo once I've written docs, the API's are all simple/consistent and packages can already be depended on - any feedback very welcome

[LorenzLamm]

The plan is to move that and other image processing utilities

• torch-image-lerp
• torch-grid-utils
• torch-fourier-slice
• torch-subpixel-crop
• torch-fourier-shell-correlation
• torch-fourier-shift
  into teamtomo once I've written docs, the API's are all simple/consistent and packages can already be depended on - any feedback very welcome

Ooh that's cool!! Thanks a lot @alisterburt .
That sounds super useful to have. I'll look into these, especially the rescaling and see whether we can replace ours with it :)

[LorenzLamm]

Hey @alisterburt , I just tested your torch-fourier-rescale package and it's working very nicely!

I am wondering what's your thought on smooth decays of the Fourier coefficients, like the cosine decay we have here:

membrain-seg/src/membrain_seg/tomo_preprocessing/matching_utils/px_matching_utils.py

Line 244 in 77f8b03

exponential_filter = cosine_ellipsoid_filter(

@rdrighetto pointed me to the problem that by cropping / padding the Fourier coefficients, some artifacts may arise, which could be avoided by some smoothing.

I see you don't have them in your package. I have never encountered any of these, though, so I'm wondering whether they are actually necessary in practice.

[alisterburt]

Yeah I agree these things can be useful but the right choice of filter is very situation dependent, often we pad and smooth circular mask images in real space instead - I would probably recommend dropping to a lower level of torch-fourier-rescale and use _rescale_rfft_2d/3d on your reweighted rfft if you want to do something like this

https://github.com/alisterburt/torch-fourier-rescale/blob/ab9cc41d3b4fa12f377a9d3f7f27fa469a63c35e/src/torch_fourier_rescale/fourier_rescale_3d.py#L46-L51

I'll expose those two functions at the base of the package so they can be imported/used without looking private too, sound good? 🙂

import torch
import torch.fft
from torch_fourier_rescale import fourier_rescale_rfft_3d
from my_package import reweight

volume = torch.rand((32, 32, 32))

# fftshift so that origin of phases in DFT is in the center of the image
# this maintains the center of the image at the position of DC component pre/post rescaling
volume = torch.fft.fftshift(volume, dim=(-3, -2, -1))
dft = torch.fft.rfftn(image, dim=(-3, -2, -1))
dft = torch.fft.fftshift(dft, dim=(-3, -2))

# do the reweighting
dft = reweight(dft)
dft = rescale_rfft_3d(dft)

# transform back to real space and recenter
dft = torch.fft.ifftshift(dft, dim=(-3, -2))
rescaled_volume = torch.fft.irfftn(dft, dim=(-3, -2, -1))
rescaled_volume = torch.fft.ifftshift(rescaled_image, dim=(-3, -2, -1))

[rdrighetto]

Whenever you apply a hard windowing operation in Fourier space (think of it as a top-hat filter), the hard edges of the window will cause ringing artifacts in real-space. Fourier cropping/padding is essentially a hard window. In the case of padding, the hard edge might be where your actual observed data meets the padding zeros. In practice the effect may be very small or even negligible. In normal tomograms you probably won't see the ringing because it's hidden in the actual density. But you might see it if you resample for example a binary segmentation and don't apply any threshold to re-binarize it afterwards.

I suggest doing a few tests in the context where this function would be used (if you haven't already) and if no noticeable artifact is observed, don't bother about it.
But if you wanna be absolutely on the safe side the cosine ellipsoid filter might be a good choice.

[alisterburt]

@rdrighetto sorry I didn't look carefully and thought you guys were talking about a fourier filter.

I totally agree that smooth real space masking to ensure there are no wraparound artifacts is important. API wise it feels like this should be up to the downstream user as it can be an important decision - a torch-shapes package would be great though for the various shapes people might want for real/fourier space masks... we should do that 🙂

there might be some useful things in https://github.com/teamtomo/libtilt/tree/main/src/libtilt/shapes too

edit: not thinking clearly - @rdrighetto I'm rereading and realising it is a fourier filter

[rdrighetto]

Now I'm confused because I was indeed talking about a Fourier filter 😅

[alisterburt]

hah! I realised as you posted - sorry, this is what I get for looking at things in the morning rush.

Are there cases we don't want this filtering or where users should think carefully about the form of their fourier filter? I'm thinking in SPA/STA where we don't want to accidentally attenuate a bunch of high frequency info

[LorenzLamm]

Cool, thanks for exposing the functions, Alister. We'll use those then.

But I agree that users should be careful with this depending on the application. For segmentation purposes, I think it's fine if high frequencies are attenuated, but of course for STA that's dangerous.

So I guess for general purpose, it should not be default?

[alisterburt]

@LorenzLamm sounds good! Have fun :-)