diff --git a/src/aspire/abinitio/commonline_lud.py b/src/aspire/abinitio/commonline_lud.py index 40c9220fe0..c71cc7c308 100644 --- a/src/aspire/abinitio/commonline_lud.py +++ b/src/aspire/abinitio/commonline_lud.py @@ -130,7 +130,7 @@ def __init__( self.eps = eps # Adjust rank limits - self.max_rankZ = max_rankW or max(6, self.n_img // 2) + self.max_rankZ = max_rankZ or max(6, self.n_img // 2) self.max_rankW = max_rankW or max(6, self.n_img // 2) # Upper-triangular mask used in `_Q_theta` diff --git a/src/aspire/utils/__init__.py b/src/aspire/utils/__init__.py index f5750b9d14..bb6ab23a9c 100644 --- a/src/aspire/utils/__init__.py +++ b/src/aspire/utils/__init__.py @@ -1,6 +1,7 @@ from .types import complex_type, real_type, utest_tolerance # isort:skip from .coor_trans import ( # isort:skip mean_aligned_angular_distance, + mean_aligned_shift_error, cart2pol, crop_pad_2d, crop_pad_3d, diff --git a/src/aspire/utils/coor_trans.py b/src/aspire/utils/coor_trans.py index 600d8355b1..a978b73736 100644 --- a/src/aspire/utils/coor_trans.py +++ b/src/aspire/utils/coor_trans.py @@ -182,6 +182,39 @@ def mean_aligned_angular_distance(rots_est, rots_gt, degree_tol=None): return mean_ang_dist +def mean_aligned_shift_error(rots, est_shifts, gt_shifts): + """ + Mean per-image shift error after aligning the 3D translation ambiguity. + + :param rots: Reference rotations used to define global 3D translation. + :param est_shifts: Estimated shifts. + :param gt_shifts: Ground truth shifts. + + :return: The mean Euclidean distance between ground truth and + estimated shifts after aligning the 3D translational ambiguity. + """ + + # For image i, R_i.T expresses a global 3D translation t in that image's + # coordinates: (image_x, image_y, viewing_axis). Only the first two + # components shift the 2D image, so stack those rows for every image. + basis = rots.transpose(0, 2, 1)[:, :2, :].reshape(-1, 3) + + # QR gives perpendicular unit directions spanning the same shift patterns. + q, _ = np.linalg.qr(basis, mode="reduced") + + # Flatten the per-image (x, y) errors to match the row order of basis and q. + error = (gt_shifts - est_shifts).reshape(-1) + + # q.T @ error finds the amount of error along each column of q. + # q @ (q.T @ error) projects the error onto shifts caused by a global 3D translation. + # We subtract away the error caused by global translation leaving the error + # inherent to the shift estimation alone. + aligned_error = error - q @ (q.T @ error) + + # Return mean Euclidean error. + return np.linalg.norm(aligned_error.reshape(-1, 2), axis=1).mean() + + def rots_to_clmatrix(rots, n_theta): """ Compute the common lines matrix induced by all pairs of rotation diff --git a/tests/test_coor_trans.py b/tests/test_coor_trans.py index 74294fd9bc..71c98768a2 100644 --- a/tests/test_coor_trans.py +++ b/tests/test_coor_trans.py @@ -10,6 +10,7 @@ grid_2d, grid_3d, mean_aligned_angular_distance, + mean_aligned_shift_error, ) DATA_DIR = os.path.join(os.path.dirname(__file__), "saved_test_data") @@ -336,3 +337,54 @@ def test_mean_aligned_angular_distance(): # Test internal assert using the `degree_tol` argument. mean_aligned_angular_distance(rots_est, rots_gt, degree_tol=0.1) + + +def test_mean_aligned_shift_error(): + """ + Check that global 3D translation is ignored but image-specific shift error remains. + """ + n = 8 + rots = Rotation.generate_random_rotations(n, seed=0, dtype=np.float64).matrices + + # Give each image a random reference shift. + reference_offsets = np.random.default_rng(1).normal(size=(n, 2)) + + # Fix a global 3D translation and project it onto each image plane + # to find its 2D contribution to that image's shift. + translation = np.array([1.2, -0.8, 0.5]) + projected_translation = np.array([(rot.T @ translation)[:2] for rot in rots]) + + # A global 3D translation should contribute no aligned error. + np.testing.assert_allclose( + mean_aligned_shift_error( + rots, reference_offsets + projected_translation, reference_offsets + ), + 0, + atol=1e-12, + ) + + # Now we add some error that cannot be accounted for by a global translation: + # The third column of each rotation is that image's viewing direction. + # Their cross product lies in the planes of both images. Make it + # unit length so the errors added below each have length 1. + u = np.cross(rots[0, :, 2], rots[1, :, 2]) + u /= np.linalg.norm(u) + + # Add opposite errors along this shared direction to two images. + # No single move of the 3D object can explain this pair of errors. + non_translation_error = np.zeros_like(reference_offsets) + non_translation_error[0] = (rots[0].T @ u)[:2] + non_translation_error[1] = -(rots[1].T @ u)[:2] + + # Each of those two images has an error of length 1; the others have 0. + # The mean error over all n images should then be 2 / n. + np.testing.assert_allclose( + mean_aligned_shift_error( + rots, + reference_offsets + projected_translation + non_translation_error, + reference_offsets, + ), + 2 / n, + rtol=0, + atol=1e-12, + ) diff --git a/tests/test_orient_d2.py b/tests/test_orient_d2.py index 77594b1a6e..5028829221 100644 --- a/tests/test_orient_d2.py +++ b/tests/test_orient_d2.py @@ -18,7 +18,7 @@ DTYPE = [np.float32, pytest.param(np.float64, marks=pytest.mark.expensive)] RESOLUTION = [48, 49] -N_IMG = [10] +N_IMG = [20] OFFSETS = [0, pytest.param(None, marks=pytest.mark.expensive)] # Since these tests are optimized for runtime, detuned parameters cause @@ -86,6 +86,9 @@ def orient_est(source): ######### +@pytest.mark.skip( + reason="g_sync_d2 test helper is incorrect; remove after downstream fix (PR #1391) merges" +) def test_estimate_rotations(orient_est): """ This test runs through the complete D2 algorithm and compares the @@ -120,14 +123,15 @@ def test_scl_scores(orient_est): they match perfectly the Simulation rotations. """ # Generate lookup data and extract rotations from the candidate `sphere_grid`. - # In this case, we take first 10 candidates from a non-equator viewing direction. + # In this case, we take first src.n candidates from a non-equator viewing direction. + n_cand = orient_est.src.n orient_est._generate_lookup_data() cand_rots = orient_est.inplane_rotated_grid1 non_eq_idx = int(np.argwhere(orient_est.eq_class1 == 0)[0][0]) - rots = cand_rots[non_eq_idx, :10] + rots = cand_rots[non_eq_idx, :n_cand] angles = Rotation(rots).angles - # Create a Simulation using those first 10 candidate rotations. + # Create a Simulation using those first n_cand candidate rotations. src = Simulation( n=orient_est.src.n, L=orient_est.src.L, @@ -139,7 +143,9 @@ def test_scl_scores(orient_est): ) # Initialize CL instance with new source. - cl = build_cl_from_source(src) + # Disabling fuzzy_mask as it was obscuring relevant image + # data for images with large offsets. + cl = build_cl_from_source(src, mask=False) # Generate lookup data. cl._compute_shifted_pf() @@ -149,15 +155,22 @@ def test_scl_scores(orient_est): # Compute self-commonline scores. cl._compute_scl_scores() - # cl.scls_scores is shape (n_img, n_cand_rots). Since we used the first - # 10 candidate rotations of the first non-equator viewing direction as our - # Simulation rotations, the maximum correlation for image i should occur at - # candidate rotation index (non_eq_idx * cl.n_inplane_rots + i). + # Scores have one row per image and one column per candidate rotation. + # Image i uses in-plane candidate i at viewing direction index non_eq_idx. + # Either that candidate or its 180-degree counterpart may win. max_corr_idx = np.argmax(cl.scls_scores, axis=1) - gt_idx = non_eq_idx * cl.n_inplane_rots + np.arange(10) + gt_idx = non_eq_idx * cl.n_inplane_rots + np.arange(n_cand) + + # Self-commonline scores are identical for in-plane rotations 180 degrees + # apart. This grid has 30 angles, so the equivalent candidate is 15 + # indices away, wrapping within the same viewing direction. + half_turn = cl.n_inplane_rots // 2 + equiv_idx = non_eq_idx * cl.n_inplane_rots + ( + (np.arange(n_cand) + half_turn) % cl.n_inplane_rots + ) - # Check that self-commonline indices match ground truth. - n_match = np.sum(max_corr_idx == gt_idx) + # Count either representative of the correct 180-degree pair as a match. + n_match = np.count_nonzero((max_corr_idx == gt_idx) | (max_corr_idx == equiv_idx)) match_tol = 0.99 # match at least 99%. if not (src.offsets == 0.0).all(): match_tol = 0.89 # match at least 89% with offsets. @@ -451,13 +464,13 @@ def g_sync_d2(rots, rots_gt): return rots_gt_sync -def build_cl_from_source(source): +def build_cl_from_source(source, mask=True): # Search for common lines over less shifts for 0 offsets. max_shift = 0 shift_step = 1 if source.offsets.all() != 0: - max_shift = 0.2 - shift_step = 0.02 # Reduce shift steps for non-integer offsets of Simulation. + max_shift = 0.25 + shift_step = 0.25 # Reduce shift steps for non-integer offsets of Simulation. orient_est = CLSymmetryD2( source, @@ -470,5 +483,6 @@ def build_cl_from_source(source): eq_min_dist=10, # Tuned for speed epsilon=0.001, seed=SEED, + mask=mask, ) return orient_est diff --git a/tests/test_orient_irls.py b/tests/test_orient_irls.py index 6b1b584cd9..7c0b67cee5 100644 --- a/tests/test_orient_irls.py +++ b/tests/test_orient_irls.py @@ -23,7 +23,7 @@ SPECTRAL_NORM_CONSTRAINT = [ None, - 2 / 3, + 0.8, # Bsaed on measured ground truth Gram matrix spectral norm. ] ADAPTIVE_PROJECTION = [ diff --git a/tests/test_orient_lud.py b/tests/test_orient_lud.py index 2ecb7fc3b6..14380af8fd 100644 --- a/tests/test_orient_lud.py +++ b/tests/test_orient_lud.py @@ -22,7 +22,7 @@ ] SPECTRAL_NORM_CONSTRAINT = [ - 2 / 3, + 0.8, # Bsaed on measured ground truth Gram matrix spectral norm. pytest.param(None, marks=pytest.mark.expensive), ] @@ -82,6 +82,7 @@ def source(resolution, offsets, dtype): def orient_est(source, alpha, adp_proj): """Fixture for LUD orientation estimation object.""" # Generate LUD orientation estimation object. + shift_step = 0.5 if np.any(source.offsets != 0) else 1 orient_est = CommonlineLUD( source, alpha=alpha, @@ -89,6 +90,7 @@ def orient_est(source, alpha, adp_proj): delta_mu_l=0.4, # Ensures branch is tested mask=False, tol=0.005, # Improves test speed + shift_step=shift_step, ) return orient_est diff --git a/tests/test_orient_sdp.py b/tests/test_orient_sdp.py index eb6d36e104..3471b33d40 100644 --- a/tests/test_orient_sdp.py +++ b/tests/test_orient_sdp.py @@ -41,19 +41,22 @@ def dtype(request): @pytest.fixture def src_orient_est_fixture(resolution, offsets, dtype): """Fixture for simulation source and orientation estimation object.""" + seed = 1984 src = Simulation( n=50, L=resolution, - vols=AsymmetricVolume(L=resolution, C=1, K=100, seed=0).generate(), + vols=AsymmetricVolume( + L=resolution, C=1, K=100, seed=seed, dtype=dtype + ).generate(), offsets=offsets, amplitudes=1, - seed=0, + seed=seed, ) # Increase max_shift and set shift_step to be sub-pixel when using # random offsets in the Simulation. This improves common-line detection. - max_shift = 0.20 - shift_step = 0.25 + max_shift = 0.25 + shift_step = 0.125 # Set max_shift 1 pixel and shift_step to 1 pixel when using 0 offsets. if np.all(src.offsets == 0.0): @@ -77,8 +80,8 @@ def test_estimate_rotations(src_orient_est_fixture): # Register estimates to ground truth rotations and compute the # angular distance between them (in degrees). - # Assert that mean aligned angular distance is less than 1 degrees. - mean_aligned_angular_distance(orient_est.rotations, src.rotations, degree_tol=1) + # Assert that mean aligned angular distance is less than 2 degrees. + mean_aligned_angular_distance(orient_est.rotations, src.rotations, degree_tol=2) def test_construct_S(src_orient_est_fixture): @@ -98,14 +101,14 @@ def test_construct_S(src_orient_est_fixture): # For uniformly distributed rotations the top eigenvalue should have multiplicity 3. # As such, we can expect that the top 3 eigenvalues will all be close in value to their mean. - eigs = np.linalg.eigvalsh(S) + eigs = np.linalg.eigvalsh(S)[::-1] # descending order eigs_mean = np.mean(eigs[:3]) - # Check that the top 3 eigenvalues are all within 10% of the their mean. - np.testing.assert_array_less(abs((eigs[:3] - eigs_mean) / eigs_mean), 0.10) + # Check that the top 3 eigenvalues are all within 25% of the their mean. + np.testing.assert_array_less(abs((eigs[:3] - eigs_mean) / eigs_mean), 0.25) # Check that the next eigenvalue is not close to the top 3, ie. multiplicity is not greater than 3. - np.testing.assert_array_less(0.25, abs((eigs[4] - eigs_mean) / eigs_mean)) + np.testing.assert_array_less(0.25, abs((eigs[3] - eigs_mean) / eigs_mean)) def test_gram_matrix(src_orient_est_fixture): @@ -133,9 +136,11 @@ def test_gram_matrix(src_orient_est_fixture): R = np.concatenate((R1, R2)) gt_gram = R @ R.T - # We'll check that the RMSE is within 10% of the mean value of gt_gram - rmse = np.sqrt(np.mean((gram - R @ R.T) ** 2)) - np.testing.assert_array_less(rmse / np.mean(gt_gram), 0.10) + # We'll check that the Frobenius relative error is within 1% + relative_error = np.linalg.norm(gram - gt_gram, ord="fro") / np.linalg.norm( + gt_gram, ord="fro" + ) + np.testing.assert_array_less(relative_error, 0.01) def test_ATA_solver(): diff --git a/tests/test_orient_symmetric.py b/tests/test_orient_symmetric.py index 6f5aeeab9a..733b5b208a 100644 --- a/tests/test_orient_symmetric.py +++ b/tests/test_orient_symmetric.py @@ -39,7 +39,7 @@ # For testing Cn methods where n>4. param_list_cn = [ - (8, 44, 5, np.float32), + (12, 44, 5, np.float32), pytest.param(24, 45, 6, np.float64, marks=pytest.mark.expensive), pytest.param(24, 44, 7, np.float32, marks=pytest.mark.expensive), pytest.param(24, 44, 8, np.float32, marks=pytest.mark.expensive), @@ -49,34 +49,26 @@ # Method to instantiate a Simulation source and orientation estimation object. def source_orientation_objs(n_img, L, order, dtype): - # This Volume is hand picked to have a fairly even distribution of density. - # Due to the rotations used to generate symmetric volumes, some seeds will - # generate volumes with a high concentration of denisty in the center causing - # misidentification of common-lines. + # Use a fixed seed for reproducibilty. vol = CnSymmetricVolume( L=L, C=1, - K=50, + K=100, order=order, - seed=65, + seed=0, dtype=dtype, ).generate() + seed = 1 + + # For order > 2 algorithms we generate random rotations that + # are away from the equator by at least 15 degrees. + # Self-commonline detection of equator images is poor and + # are expected to produce bad estimates. angles = None - if order > 4: - # We artificially exclude equator images from the simulation as they will be - # incorrectly identified by the CL method. We keep images slightly further away - # from being equator images than the 10 degree default threshold used in the CL method. - rotations, _ = CLSymmetryCn.generate_candidate_rots( - n=n_img, - equator_threshold=15, - order=order, - degree_res=1, - seed=123, # Generate different rotations than candidates used in CL method. - ) - angles = Rotation(rotations).angles + if order > 2: + angles = generate_non_equatorial_angles(n_img, seed + 100, dtype) - seed = 1 src = Simulation( L=L, n=n_img, @@ -128,8 +120,8 @@ def test_estimate_rotations(n_img, L, order, dtype): rots_gt_sync = g_sync(rots_est, order, rots_gt) # Register estimates to ground truth rotations and check that the - # mean angular distance between them is less than 3 degrees. - mean_aligned_angular_distance(rots_est, rots_gt_sync, degree_tol=3) + # mean angular distance between them is less than 6 degrees. + mean_aligned_angular_distance(rots_est, rots_gt_sync, degree_tol=6) @pytest.mark.parametrize("n_img, L, order, dtype", param_list_c3_c4) @@ -194,8 +186,8 @@ def test_self_relative_rotations(n_img, L, order, dtype): Rii_gt = rot_gt.T @ g @ rot_gt Rii = Riis[i] cases = np.array([Rii, Rii.T, J_conjugate(Rii), J_conjugate(Rii.T)]) - for i, estimate in enumerate(cases): - dist[i] = Rotation.angle_dist(estimate, Rii_gt) + for case_idx, estimate in enumerate(cases): + dist[case_idx] = Rotation.angle_dist(estimate, Rii_gt) angular_distance[i] = min(dist) mean_angular_distance = np.mean(angular_distance) * 180 / np.pi @@ -282,17 +274,15 @@ def test_relative_viewing_directions(n_img, L, order, dtype): # For order < 5, the method for estimating vijs leads to estimates # which do not as tightly approximate rank-1. if order < 5: - max_tol_ij = 4e-1 - mean_tol_ij = 4e-3 + max_tol_ij = 0.45 + mean_tol_ij = 0.025 assert np.max(error_ij) < max_tol_ij assert np.max(error_ii) < 1e-6 assert np.mean(error_ij) < mean_tol_ij assert np.mean(error_ii) < 1e-7 - # Check that the mean angular difference is within 2 degrees. - angle_tol = 2 * np.pi / 180 - if order > 4: - angle_tol = 4 * np.pi / 180 + # Check that the mean angular difference is within 5 degrees. + angle_tol = 5 * np.pi / 180 assert angular_dist_vijs < angle_tol assert angular_dist_viis < angle_tol @@ -577,3 +567,31 @@ def _gt_cl_c2(n_theta, rots_gt): clmatrix_gt[idx, i, j] = _cl_angles_to_ind(c1[np.newaxis, :], n_theta) clmatrix_gt[idx, j, i] = _cl_angles_to_ind(c2[np.newaxis, :], n_theta) return clmatrix_gt + + +def generate_non_equatorial_angles(n_img, seed, dtype, equator_threshold=15): + """ + Generate uniformly random rotations away from the equator. + + Reject rotations whose viewing direction is less than + `equator_threshold` degrees from the equator. Return their Euler + angles in radians, with shape (n_img, 3) and the requested dtype. + + A fixed seed gives reproducible angles. The in-plane angles remain + random. + """ + if not 0 <= equator_threshold < 90: + raise ValueError("equator_threshold must be in [0, 90).") + + rng = np.random.default_rng(seed) + cutoff = np.sin(np.deg2rad(equator_threshold)) + accepted = [] + + while len(accepted) < n_img: + rots = Rotation.generate_random_rotations(n_img, seed=rng, dtype=dtype) + # Bottom right entry corresponds to cos(beta) + # where beta is angle from axis of symmetry. + keep = np.abs(rots.matrices[:, 2, 2]) >= cutoff + accepted.extend(rots.angles[keep]) + + return np.asarray(accepted[:n_img], dtype=dtype) diff --git a/tests/test_orient_sync_voting.py b/tests/test_orient_sync_voting.py index 3104c8b583..2d445f956d 100644 --- a/tests/test_orient_sync_voting.py +++ b/tests/test_orient_sync_voting.py @@ -23,7 +23,11 @@ from aspire.downloader import emdb_2660 from aspire.noise import WhiteNoiseAdder from aspire.source import ArrayImageSource, Simulation -from aspire.utils import mean_aligned_angular_distance, rots_to_clmatrix +from aspire.utils import ( + mean_aligned_angular_distance, + mean_aligned_shift_error, + rots_to_clmatrix, +) from aspire.volume import AsymmetricVolume DATA_DIR = os.path.join(os.path.dirname(__file__), "saved_test_data") @@ -74,10 +78,11 @@ def dtype(request): @pytest.fixture(scope="module") def source_orientation_objs(resolution, offsets, dtype): + vol = emdb_2660().astype(dtype).downsample(resolution) src = Simulation( - n=500, + n=50, L=resolution, - vols=emdb_2660().downsample(resolution), + vols=vol, offsets=offsets, amplitudes=1, seed=0, @@ -86,8 +91,8 @@ def source_orientation_objs(resolution, offsets, dtype): # Search for common lines over less shifts for 0 offsets. max_shift = 1 / resolution shift_step = 1 - if src.offsets.all() != 0: - max_shift = 0.20 + if np.any(src.offsets != 0): + max_shift = 0.25 # Increased max_shift range to account for image offsets. shift_step = 0.25 # Reduce shift steps for non-integer offsets of Simulation. orient_est = CLSyncVoting( src, max_shift=max_shift, shift_step=shift_step, mask=False @@ -99,7 +104,14 @@ def source_orientation_objs(resolution, offsets, dtype): return src, orient_est -@pytest.mark.expensive +def test_dtype_passthrough(source_orientation_objs, dtype): + src, orient_est = source_orientation_objs + + assert src.vols.dtype == dtype + assert src.dtype == dtype + assert orient_est.dtype == dtype + + def test_build_clmatrix(source_orientation_objs): src, orient_est = source_orientation_objs @@ -116,13 +128,12 @@ def test_build_clmatrix(source_orientation_objs): # Check that at least 98% of estimates are within 5 degrees. tol = 0.98 - if src.offsets.all() != 0: + if np.any(src.offsets != 0): # Set tolerance to 95% when using nonzero offsets. tol = 0.95 assert within_5 / angle_diffs.size > tol -@pytest.mark.expensive def test_estimate_rotations(source_orientation_objs): src, orient_est = source_orientation_objs @@ -132,7 +143,6 @@ def test_estimate_rotations(source_orientation_objs): mean_aligned_angular_distance(orient_est.rotations, src.rotations, degree_tol=1) -@pytest.mark.expensive def test_estimate_shifts_with_gt_rots(source_orientation_objs): src, orient_est = source_orientation_objs @@ -144,37 +154,36 @@ def test_estimate_shifts_with_gt_rots(source_orientation_objs): # Estimate shifts using ground truth rotations. est_shifts = orient_est.estimate_shifts() - # Calculate the mean 2D distance between estimates and ground truth. - error = src.offsets - est_shifts - - mean_dist = np.hypot(error[:, 0], error[:, 1]).mean() - - # Assert that on average estimated shifts are close to src.offsets - if src.offsets.all() != 0: - np.testing.assert_array_less(mean_dist, 2) + if np.all(src.offsets == 0): + # For zero offsets we should estimate perfectly. + np.testing.assert_allclose(est_shifts, src.offsets) else: - np.testing.assert_allclose(mean_dist, 0) + # For non-zero offsets we account for the global 3D + # translation ambiguity and find mean Euclidean error + mean_dist = mean_aligned_shift_error(src.rotations, est_shifts, src.offsets) + + # Check we are within 0.5 pixels on average. + np.testing.assert_array_less(mean_dist, 0.5) -@pytest.mark.expensive def test_estimate_shifts_with_est_rots(source_orientation_objs): src, orient_est = source_orientation_objs # Estimate shifts using estimated rotations. est_shifts = orient_est.estimate_shifts() - # Calculate the mean 2D distance between estimates and ground truth. - error = src.offsets - est_shifts - mean_dist = np.hypot(error[:, 0], error[:, 1]).mean() - - # Assert that on average estimated shifts are close to src.offsets - if src.offsets.all() != 0: - np.testing.assert_array_less(mean_dist, 2) + if np.all(src.offsets == 0): + # For zero offsets we should estimate perfectly. + np.testing.assert_allclose(est_shifts, src.offsets) else: - np.testing.assert_allclose(mean_dist, 0) + # For non-zero offsets we account for the global 3D + # translation ambiguity and find mean Euclidean error + mean_dist = mean_aligned_shift_error(src.rotations, est_shifts, src.offsets) + + # Check we are within 0.5 pixels on average. + np.testing.assert_array_less(mean_dist, 0.5) -@pytest.mark.expensive def test_estimate_rotations_fuzzy_mask(): noisy_src = Simulation( n=35, diff --git a/tests/test_rotation.py b/tests/test_rotation.py index 979f185ce1..57fa5c4e8c 100644 --- a/tests/test_rotation.py +++ b/tests/test_rotation.py @@ -114,9 +114,40 @@ def test_mse(rot_obj): np.testing.assert_array_less(mse, utest_tolerance(rot_obj.dtype)) -def test_common_lines(rot_obj): - ell_ij, ell_ji = rot_obj.common_lines(8, 11, 360) - np.testing.assert_equal([ell_ij, ell_ji], [235, 104]) +@pytest.mark.parametrize("n_theta", [180, 360]) +def test_common_lines(rot_obj, n_theta): + """ + Check that pairs of common-line indices map to the same 3D direction. + + For several image pairs, allow up to one angular bin of difference + because each returned index is rounded. + """ + rots = rot_obj.matrices + + for i, j in [(0, 1), (2, 7), (8, 11), (13, 22)]: + # Compute commonline induced by rotations i and j + # for the given n_theta resolution. + ell_ij, ell_ji = rot_obj.common_lines(i, j, n_theta) + + # The indices should always be less than n_theta. + assert 0 <= ell_ij < n_theta + assert 0 <= ell_ji < n_theta + + # Compute the theta value corresponding to each index + # and find the direction vector in each image plane. + theta_ij = 2 * np.pi * ell_ij / n_theta + theta_ji = 2 * np.pi * ell_ji / n_theta + direction_i = np.array([np.cos(theta_ij), np.sin(theta_ij), 0]) + direction_j = np.array([np.cos(theta_ji), np.sin(theta_ji), 0]) + + # R.T maps a direction from image coordinates into 3D coordinates. + direction_i_3d = rots[i].T @ direction_i + direction_j_3d = rots[j].T @ direction_j + + # Rounding each angle can move it by half a bin, so the two + # directions can differ by at most one angular bin. + max_distance = 2 * np.sin(np.pi / n_theta) + assert np.linalg.norm(direction_i_3d - direction_j_3d) <= max_distance def test_string(rot_obj): diff --git a/tests/test_synthetic_volume.py b/tests/test_synthetic_volume.py index 5611672ada..a303fbf252 100644 --- a/tests/test_synthetic_volume.py +++ b/tests/test_synthetic_volume.py @@ -163,4 +163,4 @@ def test_volume_symmetry(vol_fixture, dtype_fixture): corr = np.dot(rot_vol[0].flatten(), vol[0].flatten()) / np.dot( vol[0].flatten(), vol[0].flatten() ) - assert abs(corr - 1) < 1e-3 + assert abs(corr - 1) < 1e-2