ctf model ramp weighting correction and updated interpolation center (#…

…141)

* fix interpolation center by using voltools for tilt rotation

* add exact weighting

* fix for non-ctf

* fix overlap weighting with a exact filter along the y-axis (tilt-axis)

* overlap frequency should be further out

* set weighting to numbers

* switch to use slice width

* fix use of weighting after putting along y axis

* try with D as slicewidth

* continue cutoff outside of edge

* switch to ramp weighting up to the overlap frequency

* add 0 frequency at 1 / len(tilt_angles)

* leave 0 frequency as data is anyway normalized by subtracting mean, equal to setting 0 freq in fourier transform to 0

* correctly decrease weights

* make options for contructing weights based on reconstruction method

* update version; remove optional for reconstruction method - weighting up to overlap freq should be most accurate

* add comment

* Update src/pytom_tm/weights.py

Co-authored-by: Sander Roet <[email protected]>

---------

Co-authored-by: Sander Roet <[email protected]>

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "pytom-match-pick"
-version = "0.4.3"
+version = "0.5.0"
 description = "PyTOM's GPU template matching module as an independent package"
 readme = "README.md"
 license = {file = "LICENSE"}

src/pytom_tm/weights.py

@@ -2,8 +2,10 @@
 import numpy.typing as npt
 import logging
 import scipy.ndimage as ndimage
+import voltools as vt
 from typing import Optional, Union
 from pytom_tm.io import UnequalSpacingError
+from itertools import pairwise
 
 
 constants = {
@@ -488,6 +490,20 @@ def _create_tilt_weighted_wedge(
     q_grid = radial_reduced_grid(shape)
     tilt_weighted_wedge = np.zeros((image_size, image_size, image_size // 2 + 1))
 
+    # create ramp weights to correct tilt summation for overlap
+    tilt_increment = min([abs(x - y) for x, y in pairwise(tilt_angles)])
+    # Crowther freq. determines till what point adjacent tilts overlap in Fourier space
+    overlap_frequency = 1 / (tilt_increment * image_size)
+    freq_1d = np.abs(np.arange(
+        -image_size // 2 + image_size % 2,
+        image_size // 2 + image_size % 2, 1.
+    )) / (image_size // 2) * .5  # multiply with .5 for nyquist frequency
+    ramp_filter = freq_1d / overlap_frequency
+    ramp_filter[ramp_filter > 1] = 1  # linear increase up to overlap frequency
+
+    # generate 2d weights along the tilt axis
+    ramp_weighting = np.tile(ramp_filter[:, np.newaxis], (1, image_size))
+
     for i, alpha in enumerate(tilt_angles):
         if ctf_params_per_tilt is not None:
             ctf = np.fft.fftshift(
@@ -507,19 +523,22 @@ def _create_tilt_weighted_wedge(
                     ctf
                 ),
                 axis=1
-            )
+            ) * ramp_weighting
         else:
-            tilt[:, :, image_size // 2] = 1
+            tilt[:, :, image_size // 2] = ramp_weighting
 
         # rotate the image weights to the tilt angle
         rotated = np.flip(
-            ndimage.rotate(
+            vt.transform(
                 tilt,
-                np.rad2deg(alpha),
-                axes=(0, 2),
-                reshape=False,
-                order=3
-            )[:, :, : image_size // 2 + 1]
+                rotation=(0, alpha, 0),
+                rotation_units='rad',
+                rotation_order='rxyz',
+                center=(image_size // 2, ) * 3,
+                interpolation='filt_bspline',
+                device='cpu'
+            )[:, :, :image_size // 2 + 1],  # crop back z-axis to reduced Fourier form
+            axis=2
         )
 
         # weight with exposure and tilt dampening
@@ -536,7 +555,8 @@ def _create_tilt_weighted_wedge(
                     rotated *
                     np.cos(alpha)  # apply tilt-dependent weighting
             )
-        tilt_weighted_wedge = np.maximum(tilt_weighted_wedge, weighted_tilt)
+
+        tilt_weighted_wedge += weighted_tilt
 
     tilt_weighted_wedge[q_grid > cut_off_radius] = 0