Fourier Ring/Shell Correlation module #634

daurer · 2025-09-16T07:58:39Z

This is a clean starting point based on the initial draft #606 for FRC calculation from @jcesardasilva and additional code added for FSC calculation by @Letizia97 (see also #605)

TODO

Homogenize metrics.py with two main FRC/FSC functions at the top and all utiliy functions below (@jcesardasilva )
Make sure it works for both 2d and 3d case (@jcesardasilva )
Try to merge utilitiy functions that are used for both 2D and 3D (@jcesardasilva )
Upload custom 2d alignment code to replace scikit-image (@gfev )
Automated intersection determination based on 1-bit and 0.5 bit thresholds (@gdev)
Create a tutorial notebook that uses the new module and generates nice plots (@Letizia97 @jcesardasilva )
Looks into performance optimsation (@pierrethibault )

jcesardasilva · 2025-09-16T11:48:48Z

I noticed now that there are two types of apodization in 3D:

A square-like Hanning window for the transverse cross-section of the tomogram and one circle-like for the axial cross-sections, given that the information outside that circle in tomography is basically noise or artifacts.
A square-like Hanning window in all three dimensions (a cubic-like apodization window function if you prefer). So, this ignores that, in the tomography, the information that matters in each axial cross-section is only inside a circle.

Do we need to keep both? If so, we need to implement three different kinds of apodization, the two I mentioned above and the one for the 2D images only.

pierrethibault · 2025-09-16T13:41:17Z

ptypy/utils/metrics.py

+        window1D2 = hanning_like_window(nr, apod_width)
+        window1D3 = hanning_like_window(nc, apod_width)
+
+        windowaxial = np.outer(window1D2, window1D3)


Why not use two nested outer products here?

Maybe I did it wrong, but in my case, it did not work. Could you make a quick draft of what you meant, please?

pierrethibault · 2025-09-16T13:45:20Z

ptypy/utils/metrics.py

+
+    # Apodization of the borders
+    window = apodization(input1, apod_width)
+    img1_apod = input1 * window


The input images probably aren't "whitened". Here, multiplying with a window sends the edge values to 0. I think that artifacts would be reduced if we use instead img1_apod = window * (input1 - input1.mean()) + input1.mean().

Supposing I understood your remark well, the idea of the implemented apodization here was to transition slowly to 0 at the edges, thereby avoiding a sharp transition. If the edge values are not truly zero, given that the discrete FFTs produce numerous replicas around, you may end up with the classical (input) x (rect function), leading to artifacts when implementing the discrete FFTs.

My point is that sending the signal to 0 actually can introduce artifacts if the average signal is everywhere far from 0, especially if the apodization width is small. I tried to make a drawing to explain that...

Great, thanks. Now I see what you meant. Well, the apodization steps are now coded in independent auxiliary functions. Therefore, we can easily change them to play with different approaches without having to rewrite the main code. One could try to do what you said.

Just one question: you mentioned earlier the idea of doing img1_apod = window * (input1 - input1.mean()) + input1.mean(). But, I noticed you still use the function window, which will bring the values to zero at the borders. If I understand correctly, we are supposed to remove the window from the multiplication, and do only img1_apod = (input1 - input1.mean()) + input1.mean(), right? Or would you happen to have another window function in mind? Or, by window, you mean the natural windowing (selecting the region within the Nyquist frequency range) we do when performing the discrete FFTs ?

pierrethibault · 2025-09-16T13:50:20Z

ptypy/utils/metrics.py

+
+    print("Calculating the correlation...")
+    index = ringthickness(F1) # indexes for the ring thickness
+    for ii in range(len(f)):


There is a way to make this faster using numpy's digitize and bincount. I have an example code I can show you.

I'd be happy to review your example code, which could save me time if I were to do it myself. Does the bincount work in both 2D and 3D arrays?

jcesardasilva · 2025-09-16T16:12:55Z

I made a couple of changes in the apodization functions, given their importance. You can see the commit 22e8129
Tomorrow, I will keep working on the functions.

I think we could merge fourieringcorrelation and fouriershellcorrelation, giving the capabilities to the function to identify if the images are 2D or 3D. Except if you want to keep them separated. Please let me know.

jcesardasilva · 2025-09-17T07:12:45Z

Just in case this can be of interest, I've just found out that scikit-image has a very nice and efficient window function that works in both 2D and 3D without problem: https://github.com/scikit-image/scikit-image/blob/v0.25.2/skimage/filters/_window.py#L10-L131
But of course, this would add an extra dependency.

jcesardasilva · 2025-09-17T08:01:55Z

There are now two functions, frc for the 2D FRC, and frc for the 3D FSC. They are shorter and a bit less repetitive thanks to auxiliary functions.
I have not removed the alignment routine yet, leaving the task for someone else who has access to @pierrethibault 's custom code. Commit dfcd35d
I also included changes in the Jupyter notebook, showing that both frc and fsc seem to be working locally. I still need a tomographic dataset to test the frc, but I need to wait for the unarchiving of my data. If you have some available, please go ahead with the test on your data, and I can help if needed.

daurer · 2025-09-17T12:23:55Z

Just a few comments for cleaning up at the end. We should move any example code/notebooks into the templates/notebooks folder and remove the extra stuff that has been added to .gitignore

daurer · 2025-09-17T13:39:39Z

ptypy/utils/metrics.py

+def fouriercorr_plot(FRC, T, fn, plot_name):
+    """
+    Routine to plot the FRC curves
+
+    Parameters
+    ----------
+    FRC : array-like
+        1D array containing the FRC values
+
+    T : array-like
+        1D array containing the 1-bit threshold
+
+    fn : array-like
+        1D array containing the normalized frequencies    
+    """
+
+    plt.figure()
+    plt.clf()
+    plt.plot(fn, FRC.real, "-b", label="FRC")
+    plt.plot(fn, T, "--r", label="1 bit threshold")
+    plt.legend()
+    plt.xlim(0, 1)
+    plt.ylim(0, 1.1)
+    plt.xlabel("Spatial frequency/Nyquist [normalized units]")
+    plt.ylabel("Magnitude [normalized units]")
+    plt.show()
+    plt.savefig(plot_name)
+    return None


I think we should remove this from the module, can instead show how to plot in the example notebooks. Should also remove the matplotlib import consequently

I've just removed the plotting part and the matplotlib import.

jcesardasilva · 2025-09-17T16:14:59Z

I removed the dependence on matplotlib and the plotting section. I also adjusted the docstring style in compute_intersection to match the other docstrings.

One crucial point is that the user needs the pixel value so that dividing it by the intersection value allows for an estimate of the spatial resolution. At any point in the FSC/FRC calculations, we are inputting the pixel size. Will PtyPy provide the pixel size for the FSC/FRC calculation at any stage, or will it be up to the user to determine it? In any case, it is worth showing in the tutorial that one needs the pixel size to estimate the resolution based on the intersection value.

…n compute_intersection to match the style typically used in ptypy

jcesardasilva and others added 8 commits September 16, 2025 08:55

first draft of the FRC

bd9b6b4

fix a few mistakes and adding a jupyter notebook with an example

b08d800

Add fsc

5c634d1

temp

6b1ec75

Add FSC computation

57bcf76

Refactoring

4df50f0

Crop volumes

1915474

Show the curves on same plot

c1c6db5

daurer added the Hackathon2025 label Sep 16, 2025

daurer changed the base branch from master to dev September 16, 2025 07:58

cleaning up

e3a1ba7

pierrethibault mentioned this pull request Sep 16, 2025

Add a routine for Fourier Ring correlation to the new metrics.py in utils #606

Closed

pierrethibault reviewed Sep 16, 2025

View reviewed changes

first iteration loop on the files to improve the apodization functions

22e8129

second iteration on the FSC/FRC with separated functions and shorter

dfcd35d

Letizia97 and others added 3 commits September 17, 2025 12:52

Remove not needed

34084ba

Added function to find interesction between two curves

5d3ee43

Add notebook with FSC on random volumes

e374c96

daurer commented Sep 17, 2025

View reviewed changes

remove plot part and matplotlib dependency; also modified docstring i…

ac4f6f5

…n compute_intersection to match the style typically used in ptypy

jcesardasilva force-pushed the fsc branch from 57a8c1c to ac4f6f5 Compare September 17, 2025 16:26

Fourier Ring/Shell Correlation module #634

Are you sure you want to change the base?

Fourier Ring/Shell Correlation module #634

Uh oh!

Conversation

daurer commented Sep 16, 2025 • edited by jcesardasilva Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

jcesardasilva commented Sep 16, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jcesardasilva Sep 17, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jcesardasilva commented Sep 16, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

jcesardasilva commented Sep 17, 2025

Uh oh!

jcesardasilva commented Sep 17, 2025

Uh oh!

daurer commented Sep 17, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jcesardasilva commented Sep 17, 2025

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

6 participants

daurer commented Sep 16, 2025 •

edited by jcesardasilva

Loading

jcesardasilva commented Sep 16, 2025 •

edited

Loading

jcesardasilva Sep 17, 2025 •

edited

Loading

jcesardasilva commented Sep 16, 2025 •

edited

Loading