Comments (14)
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?
from libtilt.
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
from libtilt.
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
from libtilt.
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
from libtilt.
Still not 100% the same
from libtilt.
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?
from libtilt.
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?
from libtilt.
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()```
from libtilt.
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
from libtilt.
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?
from libtilt.
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?
from libtilt.
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
from libtilt.
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
from libtilt.
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 🙂
from libtilt.
Related Issues (19)
- github issue filesystem
- FSC implementation is no longer correct HOT 1
- zyx handling in `central_slice_coordinates` doesn't work correctly HOT 1
- calculate_ctf support Tensors as input HOT 4
- backproject_fourier.py backproject_fourier(...) do_gridding_correction argument is not used HOT 1
- Add resolution threshold to correlate_dft_2d HOT 1
- grid flip is slow in central slice grid generation HOT 2
- use vmap instead of list comprehension in real space projection code
- add soft edge code breaks if device is cuda
- add 3d rescaling
- extract_central_slices_rfft performance HOT 4
- 2D +3D unified affine transforms API HOT 7
- 2D to 1D projection for common lines
- alignment API HOT 5
- sample dft 2d
- Defining generic wrapper decorators HOT 3
- unittesting device agnostic code HOT 3
- add 2D/3D rotations
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from libtilt.