Matlab CTF

[garrettwrong]

Relates to #1309

This python class will allclose reproduce MATLAB cryo_CTF_Relion filter values up to a transpose. Can flip the transpose with the indexing grid param, but my guess is we will want this ordering for the convolution application to the FT of our C order image data... TBD.

In the process of porting I noticed that our (radial) grids appears to be off by a factor of exactly 2. This is likely the scaling factor I saw comparing them before.

It may be possible to reconcile our CTFFilter, but for the time being we will at least have this one.

I will try to integrate with phase_flip etc next so I can compare the application of the filter.

from scipy.io import loadmat
import numpy as np
import aspire

n=256
V=300
DF1=2.334469921900000e+03
DF2=2.344594921900000e+03
theta= 36.700001
Cs=2.0
pixA=1.4
A=0.1

flt = aspire.operators.filters.m_CTFFilter(
    pixel_size=pixA,
    voltage=V,
    defocus_u=DF1*10,  # Angstrom
    defocus_v=DF2*10,  # Angstrom
    defocus_ang=theta, # radians
    Cs=Cs,
    alpha=A,
    B=0,
)
h = flt.evaluate_grid(n)
#print(h)

# matlab ref
# >>>  h=cryo_CTF_Relion(n,V,DF1,DF2,theta,Cs,pixA,A);

m_h = loadmat('h.mat')['h']
np.testing.assert_allclose(h, m_h)

[garrettwrong]

Removed the transpose mentioned above. Left the prior remark as strikethrough.

At a0769 achieving ~0.02% between the MATLAB and Python phase_flip application for full resolution (360x360) 10028. I suspect a chunk of that diff is actually something like a pixel_size rounding as I don't expect the FFT implementations would be off that much between numpy as matlab.

[garrettwrong]

The changes in this branch matches (diff 0 at single precision) MATLAB workflow's phase_flip for the first ten 10028 images at full 360 pixels.

I am comparing for the entire stack and will report back. Planning to discuss how we want to formally incorporate the changes.

[codecov]

Codecov Report

❌ Patch coverage is 94.73684% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 90.60%. Comparing base (87c5976) to head (d6e2c12).
⚠️ Report is 14 commits behind head on develop.

Files with missing lines	Patch %	Lines
src/aspire/operators/filters.py	93.75%	1 Missing ⚠️

Additional details and impacted files

@@           Coverage Diff            @@
##           develop    #1315   +/-   ##
========================================
  Coverage    90.60%   90.60%           
========================================
  Files          133      133           
  Lines        14394    14394           
========================================
  Hits         13042    13042           
  Misses        1352     1352           

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

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[garrettwrong]

For the entire EMPIAR 10028 star file, the max abs error across all values wrt MATLAB phase_flip implementation came to 2.3841858e-07.

[garrettwrong]

I think I've managed to adopt the MATLAB CTFFilter using our omega grids and keeping most of the radial "optimizations" by implementing some re-normalization.

I had tried moving from passing omega grids to L as the input to evaluate, but this caused all sorts of problems with radial codes (ie cov2d etc) across the repo. If we want to address that architecturally it can come as a separate PR from replicating MATLAB.

[j-c-c]

Looks good! Just have one question.

[janden]

Just one question.

[garrettwrong]

Okay, found one small dtype regression in final testing. For an example of how sensitive these codes are to singles vs doubles... Simply computing grid_2d associated arctan and sqrts in singles vs doubles takes the reference correspondence from a max abs diff of 2.3841858e-07 to 0.06. Yikes.

Going from unshifted signal code (discuessed in the PR) to shifted (centered wrapper) is just a tad slower, but maintained single precision, so I went ahead and changed that as well in case it matters for some other code I'm not aware of.

I will merge this once the tests pass and move onto a combined preprocessing workflow example to send to Yoel for his blessing.