Disagrement between relion_project and from libtilt.projection.project_fourier

[rsanchezgarc]

Hi,

I have compared what you get from libtilt.projection.project_fourier against what you get from relion_project and I am seeing some disagreements that are difficult to spot by the naked eye, but that you can tell easily by computing the difference. For instance,

I used this code

import os
import mrcfile
import torch
from scipy.spatial.transform import Rotation as R
from libtilt.projection import project_fourier
from starstack import ParticlesStarSet
from matplotlib import pyplot as plt


dirname = os.path.expanduser("~/cryo/data/preAlignedParticles/EMPIAR-10166/data")
particleIdx = 0

vol = mrcfile.read(os.path.join(dirname,"allparticles_reconstruct.mrc"))
gt_projections = ParticlesStarSet(os.path.join(dirname, "projections/proj_noCTF.star"))

anglesDegs, xyShiftAngs = gt_projections.getPose(particleIdx)
rotation_matrices = R.from_euler("ZYZ", anglesDegs, degrees=True).as_matrix()

vol = torch.FloatTensor(vol)
rotation_matrices = torch.FloatTensor(rotation_matrices).unsqueeze(0)
proj = project_fourier(
    volume = vol,
    rotation_matrices = rotation_matrices,
    rotation_matrix_zyx = False,
    pad = True,
)

img =  gt_projections[particleIdx][0]

proj = (proj-proj.mean())/proj.std()
img = (img-img.mean())/img.std()


diff = torch.abs(proj -img)
print(diff.sum())

f, axes = plt.subplots(1, 3)
axes[0].imshow(proj.cpu().squeeze(0).numpy(), cmap="gray", label="libtilt")
axes[0].set_title("libtilt")
axes[1].imshow(img, cmap="gray", label="relion")
axes[1].set_title("relion")
axes[2].imshow(diff.cpu().squeeze(0).numpy(), cmap="gray", label="abs diff")
axes[2].set_title("abs diff")

plt.show()

and this Relion command

scipion run relion_project --i allparticles_reconstruct.mrc  --ang 1000particles.star --o projections/proj_noCTF 

Do you have any ideas of why the results are not the same? Could it be perhaps the definition of the centre of coordinates?

[alisterburt]

huh - thanks for doing the comparison, I hadn't done it!

It does look visually like there is some small shift - maybe by [1, 1], are they ~identical when [-1, -1] is applied to the RELION image?

I don't expect any differences for centering as I stick to the conventions of RELION for the volume center (where the DC component of the DFT would be) but that is indeed what it looks like...

Stupid question, no shifts in the star file for RELION?

[alisterburt]

having a think, there is a rotation center for the 3D volume and one for the 2D grid of coordinates - I use the same convention for both but possibly RELION doesn't?

Will have to look at their code carefully to figure it out

[alisterburt]

question worth asking, does a projection by eulers 0, 0, 0 in RELION give the same as summing the volume down the Z dimension? This is true for my code afaik

[rsanchezgarc]

The shifts were not 0 in my first experiment (little bug on my side). Sorry for that. But still. They don't get the same results. I will post things in a second

[rsanchezgarc]

Still not 100% the same

[alisterburt]

cool! I'm glad the centering seems right - this is super interesting

What is the magnitude of these errors?

Some potential sources of difference:

• does RELION pre-correct for linear interpolation of the volume by dividing by a sinc^2 in relion_project? I know they do during refinement and I do here
• I think RELION doesn't take info from the corners of the Fourier transform, I do. Can you check the difference in Fourier space to check?

[alisterburt]

Something worth checking too, if you do a reconstruction from 1000 projections in each program how do the FSCs compare? The visual aspect of the libtilt projections looks better to me although it is subtle, what do you think?

[rsanchezgarc]

cool! I'm glad the centering seems right - this is super interesting

What is the magnitude of these errors?

Some potential sources of difference:

* does RELION pre-correct for linear interpolation of the volume by dividing by a sinc^2 in relion_project? I know they do during refinement and I do here

* I think RELION doesn't take info from the corners of the Fourier transform, I do. Can you check the difference in Fourier space to check?

The mean difference is 0.3. Some pixels are quite different. The images are normalized (mean=0, std=1)

Here you have the fourier differences. There is a nasty grid in them

f_img = torch.fft.fftshift(torch.fft.fft2(torch.FloatTensor(img)))

f, axes = plt.subplots(1, 3)
f.suptitle(f"anglesDegs {anglesDegs}, xyShiftAngs {xyShiftAngs}")
axes[0].imshow(f_proj.cpu().abs().log().numpy(), cmap="gray", label="libtilt")
axes[0].set_title("libtilt")
axes[1].imshow(f_img.cpu().abs().log().numpy(), cmap="gray", label="relion")
axes[1].set_title("relion")
axes[2].imshow((f_proj-f_img).cpu().abs().log().numpy(), cmap="gray", label="abs diff")
axes[2].set_title("abs diff")
plt.show()```

[rsanchezgarc]

Something worth checking too, if you do a reconstruction from 1000 projections in each program how do the FSCs compare? The visual aspect of the libtilt projections looks better to me although it is subtle, what do you think?

The FSC is nearly 1 in for all frequencies

[alisterburt]

really interesting - the nasty grid present in the difference looks to be present in the RELION image and not the libtilt image... but rotated 😆

So what do we have so far:

• the projections are slightly different
• the libtilt projections look more normal than the RELION projections
• both give good FSC's after a proj/backproj cycle
• the centering appears identical from each

Something that could measure which are more 'correct' is the rate of increase of signal at nyquist for each projection routine?

[rsanchezgarc]

Something that could measure which are more 'correct' is the rate of increase of signal at nyquist for each projection routine?

I am not very sure if I follow you here.

There is something that bugs me. I implemented a simple template matching protocol and I am not getting a perfect cross-correlation match between the Relion projections and libtilt projections when the angular difference is small. This is the code:

import os
from functools import lru_cache

import numpy as np
import mrcfile
import scipy.stats
import torch
from scipy.spatial.transform import Rotation as R
from libtilt.projection import project_fourier
from starstack import ParticlesStarSet
from matplotlib import pyplot as plt

dirname = os.path.expanduser("~/cryo/data/preAlignedParticles/EMPIAR-10166/data")
vol = mrcfile.read(os.path.join(dirname,"allparticles_reconstruct.mrc")
gt_projections = ParticlesStarSet(os.path.join(dirname, "projections/one_single_orientation_and_local_env.star")) 
anglesDegs, xyShiftAngs = gt_projections.getPose(range(len(gt_projections)))

rotation_matrices = R.from_euler("ZYZ", anglesDegs, degrees=True).as_matrix()
vol = torch.FloatTensor(vol)
rotation_matrices = torch.FloatTensor(rotation_matrices)
projs = project_fourier(
    volume = vol,
    rotation_matrices = rotation_matrices,
    rotation_matrix_zyx = False,
    pad = True,
)

@lru_cache()
def compare(i,j):
    out = scipy.signal.correlate2d(projs[i] - projs[i].mean(),  gt_projections[j] -  gt_projections[j].mean(), boundary='symm', mode="same").max()
    return out

out = np.eye(projs.shape[0])
for i in range(projs.shape[0]):
    for j in range(projs.shape[0]):
        out[i,j] = compare(i,j)

>> print(np.concatenate([anglesDegs, xyShiftAngs], -1))
\# Rot Tilt Psi, DeltaXAngst DeltaYAngs
[[45.  45.  45.  -0.   0. ]
 [45.5 45.5 45.5 -0.   0. ]
 [44.5 44.5 44.5 -0.   0. ]
 [45.5 46.  45.  -0.   0. ]
 [46.  46.  46.  -0.   0. ]
 [46.  44.  46.  -0.   0. ]]

>> print(out)
[[28660.291 28359.844 28471.395 28304.654 28069.785 28207.973]
 [28443.434 28328.039 28232.643 28188.611 28066.979 28158.293]
 [28544.984 28224.215 28544.799 28194.643 27897.143 28076.162]
 [28409.473 28208.508 28224.457 28222.594 27942.877 27952.312]
 [28226.586 28141.637 27979.811 27998.314 28061.146 28070.945]
 [28314.439 28188.473 28112.061 27966.465 28029.811 28485.125]]


>>  out.argmax(-1)

array([0, 0, 0, 0, 0, 5])

projs are libtilt projections and gt_projections are relion ones

The differences may be small, but still, something that we may need to worry about, since they can affect alignment accuracy, don't you think so?

[alisterburt]

Something that could measure which are more 'correct' is the rate of increase of signal at nyquist for each projection routine?

I am not very sure if I follow you here.

I mean that you could do reconstructions using libtilt from 100, 200, 300, 400, 500 particles and for each calculate the FSC with the ground truth volume - which reconstruction becomes 'better' as measured by the FSC more quickly is the one with the more accurate projections

I'm sorry I don't understand what I would expect to see in that matrix to show me which projections were better - I still only see that they are slightly different which we already agree on

[rsanchezgarc]

The matrix us the cross correlation between projections with relion and libtilt at different (but similar angles) In the ideal case, the diagonal should contain the largest values, meaning that the most similar relion abd libtilt images are the one with the same angle But this only happens for two rows, hence my concerns

[alisterburt]

right but it doesn't actually help us to answer which are 'better' - given that there is weirdness in the FTs of the RELION projections and no obvious centering issues I currently think it's more likely that the problem is with relion_project than with the libtilt projections... the experiment I proposed above would answer which projections are more accurate 🙂