stashing, estimate rotations now matches MATLAB

Author: garrettwrong

src/aspire/abinitio/commonline_sync3n.py

@@ -810,12 +810,11 @@ def _syncmatrix_ij_vote_3n(self, clmatrix, i, j, k_list, n_theta):
         :param n_theta: The number of points in the theta direction (common lines)
         :return: The (i,j) rotation block of the synchronization matrix
         """
-        good_k = self._vote_ij(clmatrix, n_theta, i, j, k_list)
+        alpha, good_k = self._vote_ij(clmatrix, n_theta, i, j, k_list)
 
-        angle = self._rotratio_eulerangle_vec(clmatrix, i, j, good_k, n_theta)
         angles = np.zeros(3)
 
-        if angle is not None:
+        if alpha is not None:
             # # # BAD
             # # Convert the Euler angles with ZYZ conversion to rotation matrices
             # angles[0] = clmatrix[i, j] * 2 * np.pi / n_theta + np.pi / 2
@@ -824,7 +823,7 @@ def _syncmatrix_ij_vote_3n(self, clmatrix, i, j, k_list, n_theta):
             # rot = Rotation.from_euler(angles).matrices
 
             angles[0] = clmatrix[i, j] * 2 * np.pi / n_theta - np.pi
-            angles[1] = angle
+            angles[1] = np.mean(alpha)
             angles[2] = np.pi - clmatrix[j, i] * 2 * np.pi / n_theta
             rot = sprot.from_euler("ZXZ", angles).as_matrix()
 

src/aspire/abinitio/sync_voting.py

@@ -35,11 +35,13 @@ def _rotratio_eulerangle_vec(self, clmatrix, i, j, good_k, n_theta):
         # cl_diff2 is for the angle on C2 created by its intersection with C1 and C3.
         # cl_diff3 is for the angle on C3 created by its intersection with C2 and C1.
         cl_diff1 = clmatrix[i, good_k] - clmatrix[i, j]  # for theta1
-        cl_diff2 = clmatrix[j, good_k] - clmatrix[j, i]  # for - theta2
+        cl_diff2 = clmatrix[j, good_k] - clmatrix[j, i]  # for theta2
         cl_diff3 = clmatrix[good_k, j] - clmatrix[good_k, i]  # for theta3
 
         # Calculate the cos values of rotation angles between i an j images for good k images
-        c_alpha, good_idx = self._get_cos_phis(cl_diff1, cl_diff2, cl_diff3, n_theta)
+        c_alpha, good_idx = self._get_cos_phis(
+            cl_diff1, cl_diff2, cl_diff3, n_theta, sync=True
+        )
         if len(c_alpha) == 0:
             return None
         alpha = np.arccos(c_alpha)
@@ -82,6 +84,7 @@ def _vote_ij(self, clmatrix, n_theta, i, j, k_list):
         k_list = k_list[
             (k_list != i) & (clmatrix[i, k_list] != -1) & (clmatrix[j, k_list] != -1)
         ]
+
         cl_idx13 = clmatrix[i, k_list]
         cl_idx31 = clmatrix[k_list, i]
         cl_idx23 = clmatrix[j, k_list]
@@ -94,12 +97,16 @@ def _vote_ij(self, clmatrix, n_theta, i, j, k_list):
         # cl_diff3 is for the angle on C3 created by its intersection with C2 and C1.
         cl_diff1 = cl_idx13 - cl_idx12
         # bad or just a trig identity?
-        cl_diff2 = cl_idx21 - cl_idx23
-        # cl_diff2 = cl_idx23 - cl_idx21  # theta2 = (clmatrix(j,K)-clmatrix(j,i)) * 2*pi/L;
+        # cl_diff2 = cl_idx21 - cl_idx23
+        cl_diff2 = (
+            cl_idx23 - cl_idx21
+        )  # theta2 = (clmatrix(j,K)-clmatrix(j,i)) * 2*pi/L;
         cl_diff3 = cl_idx32 - cl_idx31
 
         # Calculate the cos values of rotation angles between i an j images for good k images
-        cos_phi2, good_idx = self._get_cos_phis(cl_diff1, cl_diff2, cl_diff3, n_theta)
+        cos_phi2, good_idx = self._get_cos_phis(
+            cl_diff1, cl_diff2, cl_diff3, n_theta, sync=True
+        )
 
         if np.any(np.abs(cos_phi2) - 1 > 1e-12):
             logger.warning(
@@ -115,28 +122,23 @@ def _vote_ij(self, clmatrix, n_theta, i, j, k_list):
         inds = k_list[good_idx]
 
         if phis.shape[0] == 0:
-            return []
+            return None, []
 
         # Parameters used to compute the smoothed angle histogram.
         ntics = int(180 / self.hist_bin_width)
-        angles_grid = np.linspace(0, 180, ntics, True)
+        angles_grid = np.linspace(0, 180, ntics + 1, True)
         # Get angles between images i and j for computing the histogram
         angles = np.arccos(phis[:]) * 180 / np.pi
+
         # Angles that are up to 10 degrees apart are considered
         # similar. This sigma ensures that the width of the density
         # estimation kernel is roughly 10 degrees. For 15 degrees, the
         # value of the kernel is negligible.
         sigma = 3.0
 
         # Compute the histogram of the angles between images i and j
-        squared_values = np.add.outer(np.square(angles), np.square(angles_grid))
-        angles_hist = np.sum(
-            np.exp(
-                (2 * np.multiply.outer(angles, angles_grid) - squared_values)
-                / (2 * sigma**2)
-            ),
-            0,
-        )
+        angles_distances = angles_grid[None, :] - angles[:, None]
+        angles_hist = np.sum(np.exp(-(angles_distances**2) / (2 * sigma**2)), axis=0)
 
         # We assume that at the location of the peak we get the true angle
         # between images i and j. Find all third images k, that induce an
@@ -169,10 +171,11 @@ def _vote_ij(self, clmatrix, n_theta, i, j, k_list):
             idx = np.abs(angles - angles_grid[peak_idx]) < self.full_width
 
         good_k = inds[idx]
+        alpha = np.arccos(phis[idx])
 
-        return good_k.astype("int")
+        return alpha, good_k.astype("int")
 
-    def _get_cos_phis(self, cl_diff1, cl_diff2, cl_diff3, n_theta):
+    def _get_cos_phis(self, cl_diff1, cl_diff2, cl_diff3, n_theta, sync=False):
         """
         Calculate cos values of rotation angles between i and j images
 
@@ -242,7 +245,29 @@ def _get_cos_phis(self, cl_diff1, cl_diff2, cl_diff3, n_theta):
 
         # Calculated cos values of angle between i and j images
         cos_phi2 = (c3[good_idx] - c1[good_idx] * c2[good_idx]) / (
-            np.sin(theta1[good_idx]) * np.sin(theta2[good_idx])
+            np.sqrt(1 - c1[good_idx] ** 2) * np.sqrt(1 - c2[good_idx] ** 2)
         )
 
+        if sync:
+            # % Some synchronization must be applied when common line is
+            # % out by 180 degrees.
+            # % Here fix the angles between c_ij(c_ji) and c_ik(c_jk) to be smaller than pi/2,
+            # % otherwise there will be an ambiguity between alpha and pi-alpha.
+            # ind1 = ( (theta1 > pi+TOL_idx) | (theta1<-TOL_idx & theta1>-pi) );
+            # ind2 = ( (theta2 > pi+TOL_idx) | (theta2<-TOL_idx & theta2>-pi) );
+            # align180 = (ind1 & ~ind2 ) | (~ind1 & ind2 );
+            # cos_phi2(align180) = -cos_phi2(align180);
+            TOL_idx = 1e-12
+            theta1 = theta1[good_idx]
+            theta2 = theta2[good_idx]
+            theta3 = theta3[good_idx]
+            ind1 = (theta1 > (np.pi + TOL_idx)) | (
+                (theta1 < -TOL_idx) & (theta1 > -np.pi)
+            )
+            ind2 = (theta2 > (np.pi + TOL_idx)) | (
+                (theta2 < -TOL_idx) & (theta2 > -np.pi)
+            )
+            align180 = (ind1 & ~ind2) | (~ind1 & ind2)
+            cos_phi2[align180] = -cos_phi2[align180]
+
         return cos_phi2, good_idx