search

ComputationalCryoEM / ASPIRE-Python

Algorithms for Single Particle Reconstruction
http://spr.math.princeton.edu
GNU General Public License v3.0
46 stars 21 forks source link

Problem in projecting 3D volume to 2D images #196

Closed junchaoxia closed 10 months ago

junchaoxia commented 4 years ago

We decided that the rotation matrices to generate the 2D images from a 3D map need to be transposed so that the x and y axes of images will not be reversed, https://github.com/ComputationalCryoEM/ASPIRE-Python/blob/master/src/aspire/volume/__init__.py#L87. This will result in rebuilding many related unit tests but will make consistent image representation between Cov3D and other submodules such as Cov2D and orientation estimation. This is a related issue to #122.

janden commented 4 years ago

We decided that the rotation matrices to generate the 2D images from a 3D map need to be transposed so that the x and y axes of images will not be reversed,

Those are two separate things. Transposing the rotation matrices means that we're considering R^T instead of R. Reversing the x and y axes corresponds to just permuting the two axes after we've projected the image. These are two separate conventions (you can transpose R without reversing x and y and vice versa).

As you say, we should adopt the convention throughout that projecting a volume is given by the integral along z of the volume after applying R^T to the coordinates, like this:

image

We should also adopt the convention of not transposing the axes of the images (which I think we do not in most places, but some of the legacy code assumes that we do).

junchaoxia commented 4 years ago

When I switched to R^T, the following tests failed from the develop branch.

tests/test_PolarBasis2D_new.py::PolarBasis2DTestCase::testPolarBasis2DEvaluate_t ------------------------------------------------------------------ live log call ------------------------------------------------------------------ 2020/08/20 17:12:53 Represent 2D image in a polar Fourier grid FAILED [ 27%] tests/test_PolarBasis2D_new.py::PolarBasis2DTestCase::testPolarBasis2DUnit ------------------------------------------------------------------ live log call ------------------------------------------------------------------ 2020/08/20 17:12:54 Represent 2D image in a polar Fourier grid FAILED [ 27%] tests/test_anisotropic_noise.py::SimTestCase::testAnisotropicNoisePSD ------------------------------------------------------------------ live log call ------------------------------------------------------------------ 2020/08/20 17:12:59 Applying forward transformations in pipeline 2020/08/20 17:12:59 All forward transformations applied 2020/08/20 17:12:59 Loaded 512 images 2020/08/20 17:13:03 Applying forward transformations in pipeline 2020/08/20 17:13:03 All forward transformations applied 2020/08/20 17:13:03 Loaded 512 images 2020/08/20 17:13:08 Applying forward transformations in pipeline 2020/08/20 17:13:08 All forward transformations applied 2020/08/20 17:13:08 Loaded 512 images 2020/08/20 17:13:12 Applying forward transformations in pipeline 2020/08/20 17:13:12 All forward transformations applied 2020/08/20 17:13:12 Loaded 512 images FAILED [ 28%] tests/test_anisotropic_noise.py::SimTestCase::testAnisotropicNoiseVariance ------------------------------------------------------------------ live log call ------------------------------------------------------------------ 2020/08/20 17:13:17 Applying forward transformations in pipeline 2020/08/20 17:13:17 All forward transformations applied 2020/08/20 17:13:17 Loaded 512 images 2020/08/20 17:13:22 Applying forward transformations in pipeline 2020/08/20 17:13:22 All forward transformations applied 2020/08/20 17:13:22 Loaded 512 images FAILED
tests/test_coor_trans_0528.py::UtilsTestCase::testRegistRots FAILED [ 45%]

tests/test_covar3d.py::Covar3DTestCase::testCovar3D1 ------------------------------------------------------------------ live log call ------------------------------------------------------------------ 2020/08/20 17:45:40 Running Covariance Estimator 2020/08/20 17:45:41 Applying forward transformations in pipeline 2020/08/20 17:45:41 All forward transformations applied 2020/08/20 17:45:41 Loaded 512 images 2020/08/20 17:47:03 Applying forward transformations in pipeline 2020/08/20 17:47:03 All forward transformations applied 2020/08/20 17:47:03 Loaded 512 images 2020/08/20 17:48:14 Computing kernel 2020/08/20 17:50:33 Computing non-centered Fourier Transform FAILED

tests/test_mean_estimator.py::MeanEstimatorTestCase::testAdjoint ------------------------------------------------------------------ live log call ------------------------------------------------------------------ 2020/08/20 17:55:59 Expanding 3D map in a spatial-domain Fourier–Bessel basis using the direct method. 2020/08/20 17:56:03 Applying forward transformations in pipeline 2020/08/20 17:56:03 All forward transformations applied 2020/08/20 17:56:03 Loaded 512 images 2020/08/20 17:56:09 Applying forward transformations in pipeline 2020/08/20 17:56:09 All forward transformations applied 2020/08/20 17:56:09 Loaded 512 images 2020/08/20 17:56:11 Determined adjoint mappings. Shape = (98,) FAILED

tests/test_mean_estimator.py::MeanEstimatorTestCase::testEstimate ------------------------------------------------------------------ live log call ------------------------------------------------------------------ 2020/08/20 17:56:11 Expanding 3D map in a spatial-domain Fourier–Bessel basis using the direct method. 2020/08/20 17:56:15 Applying forward transformations in pipeline 2020/08/20 17:56:15 All forward transformations applied 2020/08/20 17:56:15 Loaded 512 images 2020/08/20 17:56:20 Applying forward transformations in pipeline 2020/08/20 17:56:20 All forward transformations applied 2020/08/20 17:56:20 Loaded 512 images 2020/08/20 17:56:23 Determined adjoint mappings. Shape = (98,) 2020/08/20 17:56:23 Computing kernel 2020/08/20 17:56:27 Computing non-centered Fourier Transform 2020/08/20 17:56:27 Delta 0.0797800937599991 (target 1.9990477404123765e-06) 2020/08/20 17:56:27 Delta 0.024037038259842568 (target 1.9990477404123765e-06) 2020/08/20 17:56:27 Delta 0.012623049238812842 (target 1.9990477404123765e-06) 2020/08/20 17:56:27 Delta 0.0025506936924427157 (target 1.9990477404123765e-06) 2020/08/20 17:56:27 Delta 0.0009657581819514878 (target 1.9990477404123765e-06) 2020/08/20 17:56:27 Delta 0.00041791009634349597 (target 1.9990477404123765e-06) 2020/08/20 17:56:28 Delta 0.0002269496934411429 (target 1.9990477404123765e-06) 2020/08/20 17:56:28 Delta 9.307806938833651e-05 (target 1.9990477404123765e-06) 2020/08/20 17:56:28 Delta 5.4286792500520105e-05 (target 1.9990477404123765e-06) 2020/08/20 17:56:28 Delta 2.7287394577282653e-05 (target 1.9990477404123765e-06) 2020/08/20 17:56:28 Delta 7.299243059307118e-06 (target 1.9990477404123765e-06) 2020/08/20 17:56:28 Delta 3.4908303772196455e-06 (target 1.9990477404123765e-06) 2020/08/20 17:56:28 Delta 9.687487285353967e-07 (target 1.9990477404123765e-06) FAILED

tests/test_mean_estimator.py::MeanEstimatorTestCase::testOptimize1 ------------------------------------------------------------------ live log call ------------------------------------------------------------------ 2020/08/20 17:56:28 Expanding 3D map in a spatial-domain Fourier–Bessel basis using the direct method. 2020/08/20 17:56:28 Computing kernel 2020/08/20 17:56:32 Computing non-centered Fourier Transform 2020/08/20 17:56:32 Delta 0.08045686382616518 (target 2.019646365141193e-06) 2020/08/20 17:56:32 Delta 0.024370262028064155 (target 2.019646365141193e-06) 2020/08/20 17:56:32 Delta 0.012836099220420493 (target 2.019646365141193e-06) 2020/08/20 17:56:32 Delta 0.0026842082176493933 (target 2.019646365141193e-06) 2020/08/20 17:56:32 Delta 0.0011181476346979794 (target 2.019646365141193e-06) 2020/08/20 17:56:32 Delta 0.0005297912932382796 (target 2.019646365141193e-06) 2020/08/20 17:56:32 Delta 0.00030197758269178034 (target 2.019646365141193e-06) 2020/08/20 17:56:32 Delta 0.00011546263461169434 (target 2.019646365141193e-06) 2020/08/20 17:56:32 Delta 6.455763152319231e-05 (target 2.019646365141193e-06) 2020/08/20 17:56:32 Delta 2.7378844855547492e-05 (target 2.019646365141193e-06) 2020/08/20 17:56:32 Delta 7.488630708193439e-06 (target 2.019646365141193e-06) 2020/08/20 17:56:32 Delta 3.4345387022287738e-06 (target 2.019646365141193e-06) 2020/08/20 17:56:32 Delta 1.0171736247984685e-06 (target 2.019646365141193e-06) FAILED
[ 72%] tests/test_mean_estimator.py::MeanEstimatorTestCase::testOptimize2 ------------------------------------------------------------------ live log call ------------------------------------------------------------------ 2020/08/20 17:56:32 Expanding 3D map in a spatial-domain Fourier–Bessel basis using the direct method. 2020/08/20 17:56:32 Computing kernel 2020/08/20 17:56:36 Computing non-centered Fourier Transform 2020/08/20 17:56:36 Computing Preconditioner kernel 2020/08/20 17:56:36 Circularizing kernel 2020/08/20 17:56:36 Circularizing dimension 0 2020/08/20 17:56:36 Circularizing dimension 1 2020/08/20 17:56:36 Circularizing dimension 2 2020/08/20 17:56:36 Delta 0.05739158666598119 (target 2.019646365141193e-06) 2020/08/20 17:56:36 Delta 0.015647322709427272 (target 2.019646365141193e-06) 2020/08/20 17:56:36 Delta 0.005219114855058752 (target 2.019646365141193e-06) 2020/08/20 17:56:36 Delta 0.001291617823187844 (target 2.019646365141193e-06) 2020/08/20 17:56:36 Delta 0.00036078116932163024 (target 2.019646365141193e-06) 2020/08/20 17:56:36 Delta 8.134238727098225e-05 (target 2.019646365141193e-06) 2020/08/20 17:56:36 Delta 2.3238160937359644e-05 (target 2.019646365141193e-06) 2020/08/20 17:56:36 Delta 4.298898358441506e-06 (target 2.019646365141193e-06) 2020/08/20 17:56:36 Delta 1.3244415990542887e-06 (target 2.019646365141193e-06) FAILED

tests/test_simulation.py::SimTestCase::testSimulationImages ------------------------------------------------------------------ live log call ------------------------------------------------------------------ 2020/08/20 18:02:01 Appending a NoiseAdder to generation pipeline FAILED [ 86%] tests/test_simulation.py::SimTestCase::testSimulationImagesDownsample ------------------------------------------------------------------ live log call ------------------------------------------------------------------ 2020/08/20 18:02:06 Appending a NoiseAdder to generation pipeline 2020/08/20 18:02:06 Setting max. resolution of source = 8 FAILED [ 86%] tests/test_simulation.py::SimTestCase::testSimulationImagesNoisy ------------------------------------------------------------------ live log call ------------------------------------------------------------------ 2020/08/20 18:02:10 Appending a NoiseAdder to generation pipeline 2020/08/20 18:02:13 Applying forward transformations in pipeline 2020/08/20 18:02:13 All forward transformations applied 2020/08/20 18:02:13 Loaded 512 images FAILED

tests/test_white_noise.py::SimTestCase::testWhiteNoise ------------------------------------------------------------------ live log call ------------------------------------------------------------------ 2020/08/20 18:02:17 Determining Noise variance in batches of 512 2020/08/20 18:02:22 Applying forward transformations in pipeline 2020/08/20 18:02:22 All forward transformations applied 2020/08/20 18:02:22 Loaded 512 images 2020/08/20 18:02:26 Applying forward transformations in pipeline 2020/08/20 18:02:26 All forward transformations applied 2020/08/20 18:02:26 Loaded 512 images 2020/08/20 18:02:26 Noise variance = 0.0031053894192200108 FAILED

I expect that test_orient_sync.py will fail too in the orientation branch. We can start to fix them after the PR #163 was merged.

janden commented 4 years ago

Ok so why are they failing exactly? Is it because we have the rotations hardcoded in the tests? It might be necessary to transpose them in the tests then to make things pass.

junchaoxia commented 4 years ago

When we change from R to R^T. The resulted 2D images will switch the x, y. We need to make the corrections to related unit tests.

janden commented 4 years ago

The resulted 2D images will switch the x, y.

Why are the images getting transposed? I thought this was dealing with transposing the rotation matrices, not the images.

We need to make the corrections to related unit tests.

Yes if we have hardcoded rotations, those may need to be replaced with their transposes for certain tests.

junchaoxia commented 4 years ago

In line 87 as below: https://github.com/ComputationalCryoEM/ASPIRE-Python/blob/master/src/aspire/volume/__init__.py#L87

We discussed that rot_matrices[i, :, :] @ pts should be rot_matrices[i, :, :].T @ pts for the change from R to R^T. This resulted projections will switch x, y. So if the hardcoded unit tests depend on this, we need to make change. I believe some of them just need switching x,y in results of unit tests. I will try after the big PR is merged.

janden commented 4 years ago

Changing from R to R^T shouldn't invert x and y. Are you sure that's what's happening?

garrettwrong commented 2 years ago

@j-c-c , this relates to some discussions we've had and also your current discussions with Joakim so I am assigning to you.

janden commented 2 years ago

First thing here would be to clarify what convention the method is actually doing (i.e., is it using R or R^T? in our z-y-x convention, does it transpose x and y when projecting?). Then we can decide whether we want to change it and look at the consequences.

garrettwrong commented 1 year ago

Adding a note from our chat. This diff keeps the tests passing, but as you can see it is a little strange. Maybe that helps when combined with the stuff you've been looking at already.

diff --git a/src/aspire/volume/volume.py b/src/aspire/volume/volume.py
index c7738db9..fcf5729d 100644
--- a/src/aspire/volume/volume.py
+++ b/src/aspire/volume/volume.py
@@ -658,19 +658,19 @@ def rotated_grids(L, rot_matrices):
         Frequencies are in the range [-pi, pi].
     """

-    grid2d = grid_2d(L, indexing="xy", dtype=rot_matrices.dtype)
+    grid2d = grid_2d(L, indexing="yx", dtype=rot_matrices.dtype)
     num_pts = L**2
     num_rots = rot_matrices.shape[0]
     pts = np.pi * np.vstack(
         [
+            np.zeros(num_pts, dtype=rot_matrices.dtype),
             grid2d["x"].flatten(),
             grid2d["y"].flatten(),
-            np.zeros(num_pts, dtype=rot_matrices.dtype),
         ]
     )
     pts_rot = np.zeros((3, num_rots, num_pts), dtype=rot_matrices.dtype)
     for i in range(num_rots):
-        pts_rot[:, i, :] = rot_matrices[i, :, :] @ pts
+        pts_rot[:, i, :] = rot_matrices[i, :, :] @ pts[::-1]

     pts_rot = pts_rot.reshape((3, num_rots, L, L))
j-c-c commented 11 months ago

Adding some findings here. This is the diff that produces the same projection images as Relion:

diff --git a/src/aspire/volume/volume.py b/src/aspire/volume/volume.py
index c7738db9..365cca80 100644
--- a/src/aspire/volume/volume.py
+++ b/src/aspire/volume/volume.py
@@ -658,7 +658,7 @@ def rotated_grids(L, rot_matrices):
         Frequencies are in the range [-pi, pi].
     """

-    grid2d = grid_2d(L, indexing="xy", dtype=rot_matrices.dtype)
+    grid2d = grid_2d(L, indexing="yx", dtype=rot_matrices.dtype)
     num_pts = L**2
     num_rots = rot_matrices.shape[0]
     pts = np.pi * np.vstack(
@@ -672,7 +672,7 @@ def rotated_grids(L, rot_matrices):
     for i in range(num_rots):
         pts_rot[:, i, :] = rot_matrices[i, :, :] @ pts

-    pts_rot = pts_rot.reshape((3, num_rots, L, L))
+    pts_rot = pts_rot.reshape((3, num_rots, L, L))[::-1]

     return pts_rot

I'll attach a minimal script that compares Relion generated projections with ASPIRE Simulation projections. For details on how the Relion projections were generated see issue #1026.

Here is a link to the starfile for the RelionSource: https://drive.google.com/drive/folders/1fATzVomr5Q2oe-eNAlGbnQULDH5UXRe8?usp=sharing

It appears that the contrast is flipped between the two sets of images, but the orientations of the images are aligned. rln_asp_proj_comp_script.py.txt euler_angles.npy.txt

j-c-c commented 10 months ago

This issue was resolved in #1030. To summarize, the rotations provided to the Volume.project() method rotate the underlying grid to produce projection images. In other words, they correspond to R.T from the formula in the second comment.

It is also noteworthy to mention that our projections now correspond to projections generated by relion given the same set of rotations.