diff --git a/src/aspire/abinitio/J_sync.py b/src/aspire/abinitio/J_sync.py
new file mode 100644
index 0000000000..b8458d9807
--- /dev/null
+++ b/src/aspire/abinitio/J_sync.py
@@ -0,0 +1,247 @@
+import logging
+
+import numpy as np
+from numpy.linalg import norm
+
+from aspire.utils import J_conjugate, all_pairs, all_triplets, tqdm
+from aspire.utils.random import randn
+
+logger = logging.getLogger(__name__)
+
+
+class JSync:
+    """
+    Class for handling J-synchronization methods.
+    """
+
+    def __init__(
+        self,
+        n,
+        epsilon=1e-2,
+        max_iters=1000,
+        seed=None,
+    ):
+        """
+        Initialize JSync object for estimating global handedness synchronization for a
+        set of relative rotations, Rij = Ri @ Rj.T, where i <= j = 0, 1, ..., n.
+
+        :param n: Number of images/rotations.
+        :param epsilon: Tolerance for the power method.
+        :param max_iters: Maximum iterations for the power method.
+        :param seed: Optional seed for power method initial random vector.
+        """
+        self.n_img = n
+        self.epsilon = epsilon
+        self.max_iters = max_iters
+        self.seed = seed
+
+    def global_J_sync(self, vijs):
+        """
+        Global J-synchronization of all third row outer products. Given 3x3 matrices vijs, each
+        of which might contain a spurious J (ie. vij = J*vi*vj^T*J instead of vij = vi*vj^T),
+        we return vijs that all have either a spurious J or not.
+
+        :param vijs: An (n-choose-2)x3x3 array where each 3x3 slice holds an estimate for the corresponding
+        outer-product vi*vj^T between the third rows of the rotation matrices Ri and Rj. Each estimate
+        might have a spurious J independently of other estimates.
+
+        :return: vijs, all of which have a spurious J or not.
+        """
+
+        # Determine relative handedness of vijs.
+        sign_ij_J = self.power_method(vijs)
+
+        # Synchronize vijs
+        vijs_sync = vijs.copy()
+        for i, sign in enumerate(sign_ij_J):
+            if sign == -1:
+                vijs_sync[i] = J_conjugate(vijs[i])
+
+        return vijs_sync
+
+    def power_method(self, vijs):
+        """
+        Calculate the leading eigenvector of the J-synchronization matrix
+        using the power method.
+
+        As the J-synchronization matrix is of size (n-choose-2)x(n-choose-2), we
+        use the power method to compute the eigenvalues and eigenvectors,
+        while constructing the matrix on-the-fly.
+
+        :param vijs: (n-choose-2)x3x3 array of estimates of relative orientation matrices.
+
+        :return: An array of length n-choose-2 consisting of 1 or -1, where the sign of the
+            i'th entry indicates whether the i'th relative orientation matrix will be J-conjugated.
+        """
+
+        # Set power method tolerance and maximum iterations.
+        epsilon = self.epsilon
+        max_iters = self.max_iters
+
+        # Initialize candidate eigenvectors
+        n_vijs = vijs.shape[0]
+        vec = randn(n_vijs, seed=self.seed)
+        vec = vec / norm(vec)
+        residual = 1
+        itr = 0
+
+        # Power method iterations
+        logger.info(
+            "Initiating power method to estimate J-synchronization matrix eigenvector."
+        )
+        while itr < max_iters and residual > epsilon:
+            itr += 1
+            # Note, this appears to need double precision for accuracy in the following division.
+            vec_new = self._signs_times_v(vijs, vec).astype(np.float64, copy=False)
+            vec_new = vec_new / norm(vec_new)
+            residual = norm(vec_new - vec)
+            vec = vec_new
+            logger.info(
+                f"Iteration {itr}, residual {round(residual, 5)} (target {epsilon})"
+            )
+
+        # We need only the signs of the eigenvector
+        J_sync = np.sign(vec, dtype=vijs.dtype)
+
+        return J_sync
+
+    def sync_viis(self, vijs, viis):
+        """
+        Given a set of synchronized pairwise outer products vijs, J-synchronize the set of
+        outer products viis.
+
+        :param vijs: An (n-choose-2)x3x3 array where each 3x3 slice holds an estimate for the corresponding
+        outer-product vi*vj^T between the third rows of the rotation matrices Ri and Rj. Each estimate
+        might have a spurious J independently of other estimates.
+
+        :param viis: An n_imgx3x3 array where the i'th slice holds an estimate for the outer product vi*vi^T
+        between the third row of matrix Ri and itself. Each estimate might have a spurious J independently
+        of other estimates.
+
+        :return: J-synchronized viis.
+        """
+
+        # Synchronize viis
+        # We use the fact that if v_ii and v_ij are of the same handedness, then v_ii @ v_ij = v_ij.
+        # If they are opposite handed then Jv_iiJ @ v_ij = v_ij. We compare each v_ii against all
+        # previously synchronized v_ij to get a consensus on the handedness of v_ii.
+        _, pairs_to_linear = all_pairs(self.n_img, return_map=True)
+        for i in range(self.n_img):
+            vii = viis[i]
+            vii_J = J_conjugate(vii)
+            J_consensus = 0
+            for j in range(self.n_img):
+                if j < i:
+                    idx = pairs_to_linear[j, i]
+                    vji = vijs[idx]
+
+                    err1 = norm(vji @ vii - vji)
+                    err2 = norm(vji @ vii_J - vji)
+
+                elif j > i:
+                    idx = pairs_to_linear[i, j]
+                    vij = vijs[idx]
+
+                    err1 = norm(vii @ vij - vij)
+                    err2 = norm(vii_J @ vij - vij)
+
+                else:
+                    continue
+
+                # Accumulate J consensus
+                if err1 < err2:
+                    J_consensus -= 1
+                else:
+                    J_consensus += 1
+
+            if J_consensus > 0:
+                viis[i] = vii_J
+        return viis
+
+    def _signs_times_v(self, vijs, vec):
+        """
+        Multiplication of the J-synchronization matrix by a candidate eigenvector.
+
+        The J-synchronization matrix is a matrix representation of the handedness graph, Gamma, whose set of
+        nodes consists of the estimates vijs and whose set of edges consists of the undirected edges between
+        all triplets of estimates vij, vjk, and vik, where i<j<k. The weight of an edge is set to +1 if its
+        incident nodes agree in handednes and -1 if not.
+
+        The J-synchronization matrix is of size (n-choose-2)x(n-choose-2), where each entry corresponds to
+        the relative handedness of vij and vjk. The entry (ij, jk), where ij and jk are retrieved from the
+        all_pairs indexing, is 1 if vij and vjk are of the same handedness and -1 if not. All other entries
+        (ij, kl) hold a zero.
+
+        Due to the large size of the J-synchronization matrix we construct it on the fly as follows.
+        For each triplet of outer products vij, vjk, and vik, the associated elements of the J-synchronization
+        matrix are populated with +1 or -1 and multiplied by the corresponding elements of
+        the current candidate eigenvector supplied by the power method. The new candidate eigenvector
+        is updated for each triplet.
+
+        :param vijs: (n-choose-2)x3x3 array, where each 3x3 slice holds the outer product of vi and vj.
+
+        :param vec: The current candidate eigenvector of length n-choose-2 from the power method.
+
+        :return: New candidate eigenvector of length n-choose-2. The product of the J-sync matrix and vec.
+        """
+
+        # All pairs (i,j) and triplets (i,j,k) where i<j<k
+        n_img = self.n_img
+        triplets = all_triplets(n_img)
+        pairs, pairs_to_linear = all_pairs(n_img, return_map=True)
+
+        # There are 4 possible configurations of relative handedness for each triplet (vij, vjk, vik).
+        # 'conjugate' expresses which node of the triplet must be conjugated (True) to achieve synchronization.
+        conjugate = np.empty((4, 3), bool)
+        conjugate[0] = [False, False, False]
+        conjugate[1] = [True, False, False]
+        conjugate[2] = [False, True, False]
+        conjugate[3] = [False, False, True]
+
+        # 'edges' corresponds to whether conjugation agrees between the pairs (vij, vjk), (vjk, vik),
+        # and (vik, vij). True if the pairs are in agreement, False otherwise.
+        edges = np.empty((4, 3), bool)
+        edges[:, 0] = conjugate[:, 0] == conjugate[:, 1]
+        edges[:, 1] = conjugate[:, 1] == conjugate[:, 2]
+        edges[:, 2] = conjugate[:, 2] == conjugate[:, 0]
+
+        # The corresponding entries in the J-synchronization matrix are +1 if the pair of nodes agree, -1 if not.
+        edge_signs = np.where(edges, 1, -1)
+
+        # For each triplet of nodes we apply the 4 configurations of conjugation and determine the
+        # relative handedness based on the condition that vij @ vjk - vik = 0 for synchronized nodes.
+        # We then construct the corresponding entries of the J-synchronization matrix with 'edge_signs'
+        # corresponding to the conjugation configuration producing the smallest residual for the above
+        # condition. Finally, we the multiply the 'edge_signs' by the cooresponding entries of 'vec'.
+        v = vijs
+        new_vec = np.zeros_like(vec)
+        pbar = tqdm(desc="Computing signs_times_v", total=len(triplets))
+        for i, j, k in triplets:
+            ij = pairs_to_linear[i, j]
+            jk = pairs_to_linear[j, k]
+            ik = pairs_to_linear[i, k]
+            vij, vjk, vik = v[ij], v[jk], v[ik]
+            vij_J = J_conjugate(vij)
+            vjk_J = J_conjugate(vjk)
+            vik_J = J_conjugate(vik)
+
+            conjugated_pairs = np.where(
+                conjugate[..., np.newaxis, np.newaxis],
+                [vij_J, vjk_J, vik_J],
+                [vij, vjk, vik],
+            )
+            residual = np.stack([norm(x @ y - z) for x, y, z in conjugated_pairs])
+
+            min_residual = np.argmin(residual)
+
+            # Assign edge weights
+            s_ij_jk, s_ik_jk, s_ij_ik = edge_signs[min_residual]
+
+            # Update multiplication of signs times vec
+            new_vec[ij] += s_ij_jk * vec[jk] + s_ij_ik * vec[ik]
+            new_vec[jk] += s_ij_jk * vec[ij] + s_ik_jk * vec[ik]
+            new_vec[ik] += s_ij_ik * vec[ij] + s_ik_jk * vec[jk]
+            pbar.update()
+        pbar.close()
+
+        return new_vec
diff --git a/src/aspire/abinitio/__init__.py b/src/aspire/abinitio/__init__.py
index 3685479732..8ec67c0adc 100644
--- a/src/aspire/abinitio/__init__.py
+++ b/src/aspire/abinitio/__init__.py
@@ -1,16 +1,13 @@
-from .commonline_base import CLOrient3D
-
 # isort: off
+from .J_sync import JSync
 from .commonline_utils import (
-    cl_angles_to_ind,
-    estimate_third_rows,
-    complete_third_row_to_rot,
-    estimate_inplane_rotations,
+    build_outer_products,
+    g_sync,
 )
+from .commonline_base import CLOrient3D
 from .commonline_sdp import CommonlineSDP
 from .commonline_lud import CommonlineLUD
 from .commonline_irls import CommonlineIRLS
-from .sync_voting import SyncVotingMixin
 from .commonline_sync import CLSyncVoting
 from .commonline_sync3n import CLSync3N
 from .commonline_c3_c4 import CLSymmetryC3C4
diff --git a/src/aspire/abinitio/commonline_base.py b/src/aspire/abinitio/commonline_base.py
index b179a1d7df..2947195769 100644
--- a/src/aspire/abinitio/commonline_base.py
+++ b/src/aspire/abinitio/commonline_base.py
@@ -10,6 +10,8 @@
 from aspire.utils import Rotation, complex_type, fuzzy_mask, tqdm
 from aspire.utils.random import choice
 
+from .commonline_utils import _generate_shift_phase_and_filter
+
 logger = logging.getLogger(__name__)
 
 
@@ -310,8 +312,8 @@ def build_clmatrix_host(self):
 
         # Prepare the shift phases to try and generate filter for common-line detection
         r_max = pf.shape[2]
-        shifts, shift_phases, h = self._generate_shift_phase_and_filter(
-            r_max, max_shift, shift_step
+        shifts, shift_phases, h = _generate_shift_phase_and_filter(
+            r_max, max_shift, shift_step, self.dtype
         )
 
         # Apply bandpass filter, normalize each ray of each image
@@ -409,8 +411,8 @@ def build_clmatrix_cu(self):
         #
         # Note the CUDA implementation has been optimized to not
         # compute or return diagnostic 1d shifts.
-        _, shift_phases, h = self._generate_shift_phase_and_filter(
-            r, self.max_shift, self.shift_step
+        _, shift_phases, h = _generate_shift_phase_and_filter(
+            r, self.max_shift, self.shift_step, self.dtype
         )
         # Transfer to device, dtypes must match kernel header.
         shift_phases = cp.asarray(shift_phases, dtype=complex_type(self.dtype))
@@ -558,8 +560,8 @@ def _get_shift_equations_approx(self):
         # applied to maximize the common line calculation. The common-line filter
         # is also applied to the radial direction for easier detection.
         r_max = pf.shape[2]
-        _, shift_phases, h = self._generate_shift_phase_and_filter(
-            r_max, self.offsets_max_shift, self.offsets_shift_step
+        _, shift_phases, h = _generate_shift_phase_and_filter(
+            r_max, self.offsets_max_shift, self.offsets_shift_step, self.dtype
         )
 
         d_theta = np.pi / n_theta_half
@@ -696,34 +698,6 @@ def _estimate_num_shift_equations(self, n_img):
 
         return n_equations
 
-    def _generate_shift_phase_and_filter(self, r_max, max_shift, shift_step):
-        """
-        Prepare the shift phases and generate filter for common-line detection
-
-        The shift phases are pre-defined in a range of max_shift that can be
-        applied to maximize the common line calculation. The common-line filter
-        is also applied to the radial direction for easier detection.
-
-        :param r_max: Maximum index for common line detection
-        :param max_shift: Maximum value of 1D shift (in pixels) to search
-        :param shift_step: Resolution of shift estimation in pixels
-        :return: shift phases matrix and common lines filter
-        """
-
-        # Number of shifts to try
-        n_shifts = int(np.ceil(2 * max_shift / shift_step + 1))
-
-        # only half of ray, excluding the DC component.
-        rk = np.arange(1, r_max + 1, dtype=self.dtype)
-
-        # Generate all shift phases
-        shifts = -max_shift + shift_step * np.arange(n_shifts, dtype=self.dtype)
-        shift_phases = np.exp(np.outer(shifts, -2 * np.pi * 1j * rk / (2 * r_max + 1)))
-        # Set filter for common-line detection
-        h = np.sqrt(np.abs(rk)) * np.exp(-np.square(rk) / (2 * (r_max / 4) ** 2))
-
-        return shifts, shift_phases, h
-
     def _generate_index_pairs(self, n_equations):
         """
         Generate two index lists for [i, j] pairs of images
diff --git a/src/aspire/abinitio/commonline_c2.py b/src/aspire/abinitio/commonline_c2.py
index 7bd3495b84..febd63951a 100644
--- a/src/aspire/abinitio/commonline_c2.py
+++ b/src/aspire/abinitio/commonline_c2.py
@@ -3,17 +3,20 @@
 import numpy as np
 from scipy.linalg import eigh
 
-from aspire.abinitio import (
-    CLSymmetryC3C4,
-    complete_third_row_to_rot,
-    estimate_third_rows,
-)
+from aspire.abinitio import CLOrient3D, JSync
+from aspire.abinitio.sync_voting import _syncmatrix_ij_vote_3n
 from aspire.utils import J_conjugate, Rotation, all_pairs
 
+from .commonline_utils import (
+    _complete_third_row_to_rot,
+    _estimate_third_rows,
+    _generate_shift_phase_and_filter,
+)
+
 logger = logging.getLogger(__name__)
 
 
-class CLSymmetryC2(CLSymmetryC3C4):
+class CLSymmetryC2(CLOrient3D):
     """
     Define a class to estimate 3D orientations using common lines methods for molecules with C2 cyclic symmetry.
 
@@ -44,7 +47,6 @@ def __init__(
         shift_step=1,
         epsilon=1e-3,
         max_iters=1000,
-        degree_res=1,
         min_dist_cls=25,
         seed=None,
         mask=True,
@@ -59,24 +61,19 @@ def __init__(
         :param max_shift: Maximum range for shifts as a proportion of resolution. Default = 0.15.
         :param shift_step: Resolution of shift estimation in pixels. Default = 1 pixel.
         :param epsilon: Tolerance for the power method.
-        :param max_iter: Maximum iterations for the power method.
-        :param degree_res: Degree resolution for estimating in-plane rotations.
+        :param max_iters: Maximum iterations for the power method.
         :param min_dist_cls: Minimum distance between mutual common-lines. Default = 25 degrees.
         :param seed: Optional seed for RNG.
         :param mask: Option to mask `src.images` with a fuzzy mask (boolean).
             Default, `True`, applies a mask.
         """
+
         super().__init__(
             src,
-            symmetry="C2",
             n_rad=n_rad,
             n_theta=n_theta,
             max_shift=max_shift,
             shift_step=shift_step,
-            epsilon=epsilon,
-            max_iters=max_iters,
-            degree_res=degree_res,
-            seed=seed,
             mask=mask,
             **kwargs,
         )
@@ -84,16 +81,10 @@ def __init__(
         self.min_dist_cls = min_dist_cls
         self.epsilon = epsilon
         self.max_iters = max_iters
-        self.degree_res = degree_res
         self.seed = seed
         self.order = 2
 
-    def _check_symmetry(self, symmetry):
-        symmetry = symmetry.upper()
-        if symmetry != "C2":
-            raise NotImplementedError(
-                f"Only C2 symmetry supported. {symmetry} was supplied."
-            )
+        self.J_sync = JSync(src.n, self.epsilon, self.max_iters, self.seed)
 
     def build_clmatrix(self):
         """
@@ -120,8 +111,8 @@ def build_clmatrix(self):
 
         # Prepare the shift phases and generate filter for common-line detection.
         r_max = pf.shape[2]
-        shifts, shift_phases, h = self._generate_shift_phase_and_filter(
-            r_max, self.max_shift, self.shift_step
+        shifts, shift_phases, h = _generate_shift_phase_and_filter(
+            r_max, self.max_shift, self.shift_step, self.dtype
         )
         n_shifts = len(shifts)
 
@@ -230,7 +221,7 @@ def estimate_rotations(self):
         viis = np.vstack((np.eye(3, dtype=self.dtype),) * self.n_img).reshape(
             self.n_img, 3, 3
         )
-        vis = estimate_third_rows(vijs, viis)
+        vis = _estimate_third_rows(vijs, viis)
 
         logger.info("Estimating in-plane rotations and rotations matrices.")
         Ris = self._estimate_inplane_rotations(vis, Rijs, Rijgs)
@@ -253,7 +244,7 @@ def _estimate_relative_viewing_directions(self):
         self.build_clmatrix()
 
         # Step 2: Calculate relative rotations associated with both mutual common lines.
-        Rijs, Rijgs = self._estimate_all_Rijs_c2(self.clmatrix)
+        Rijs, Rijgs = self._estimate_all_Rijs_c2()
 
         # Step 3: Inner J-synchronization
         Rijs, Rijgs = self._local_J_sync_c2(Rijs, Rijgs)
@@ -282,7 +273,7 @@ def _global_J_sync(self, Rijs, Rijgs):
         vijs = (Rijs + Rijgs) / 2
 
         # Determine relative handedness of vijs.
-        sign_ij_J = self._J_sync_power_method(vijs)
+        sign_ij_J = self.J_sync.power_method(vijs)
 
         # Synchronize relative rotations
         for i, sign in enumerate(sign_ij_J):
@@ -307,7 +298,7 @@ def _estimate_inplane_rotations(self, vis, Rijs, Rijgs):
         H = np.zeros((self.n_img, self.n_img), dtype=complex)
         # Step 1: Construct all rotation matrices Ris_tilde whose third rows are equal to
         # the corresponding third rows vis.
-        Ris_tilde = complete_third_row_to_rot(vis)
+        Ris_tilde = _complete_third_row_to_rot(vis)
 
         pairs = all_pairs(self.n_img)
         for idx, (i, j) in enumerate(pairs):
@@ -350,23 +341,36 @@ def _estimate_inplane_rotations(self, vis, Rijs, Rijgs):
     # Secondary Methods for computing outer product #
     #################################################
 
-    def _estimate_all_Rijs_c2(self, clmatrix):
+    def _estimate_all_Rijs_c2(self):
         """
         Estimate the two sets of relative rotations, Rijs and Rijgs, between pairs
         of images using the voting method.
 
-        :param clmatrix: 2 x n_img x n_img array holding two sets of mutual common-lines
-            between pairs of images.
         :return: Relative rotations, Rijs and Rijgs.
         """
         k_list = np.arange(self.n_img)
-        n_theta = self.n_theta
         pairs = all_pairs(self.n_img)
         Rijs = np.zeros((len(pairs), 3, 3), dtype=self.dtype)
         Rijgs = np.zeros((len(pairs), 3, 3), dtype=self.dtype)
         for idx, (i, j) in enumerate(pairs):
-            Rijs[idx] = self._syncmatrix_ij_vote_3n(clmatrix[0], i, j, k_list, n_theta)
-            Rijgs[idx] = self._syncmatrix_ij_vote_3n(clmatrix[1], i, j, k_list, n_theta)
+            Rijs[idx] = _syncmatrix_ij_vote_3n(
+                self.clmatrix[0],
+                i,
+                j,
+                k_list,
+                self.n_theta,
+                self.hist_bin_width,
+                self.full_width,
+            )
+            Rijgs[idx] = _syncmatrix_ij_vote_3n(
+                self.clmatrix[1],
+                i,
+                j,
+                k_list,
+                self.n_theta,
+                self.hist_bin_width,
+                self.full_width,
+            )
 
         return Rijs, Rijgs
 
diff --git a/src/aspire/abinitio/commonline_c3_c4.py b/src/aspire/abinitio/commonline_c3_c4.py
index c581005c38..34d0ee4763 100644
--- a/src/aspire/abinitio/commonline_c3_c4.py
+++ b/src/aspire/abinitio/commonline_c3_c4.py
@@ -1,30 +1,23 @@
 import logging
 
 import numpy as np
-from numpy.linalg import eigh, norm, svd
+from numpy.linalg import norm, svd
 
-from aspire.abinitio import (
-    CLOrient3D,
-    SyncVotingMixin,
-    estimate_inplane_rotations,
-    estimate_third_rows,
-)
+from aspire.abinitio import CLOrient3D, JSync
+from aspire.abinitio.sync_voting import _syncmatrix_ij_vote_3n
 from aspire.operators import PolarFT
-from aspire.utils import (
-    J_conjugate,
-    Rotation,
-    all_pairs,
-    all_triplets,
-    anorm,
-    cyclic_rotations,
-    trange,
+from aspire.utils import J_conjugate, Rotation, all_pairs, anorm, trange
+
+from .commonline_utils import (
+    _estimate_inplane_rotations,
+    _estimate_third_rows,
+    _generate_shift_phase_and_filter,
 )
-from aspire.utils.random import randn
 
 logger = logging.getLogger(__name__)
 
 
-class CLSymmetryC3C4(CLOrient3D, SyncVotingMixin):
+class CLSymmetryC3C4(CLOrient3D):
     """
     Define a class to estimate 3D orientations using common lines methods for molecules with
     C3 and C4 cyclic symmetry.
@@ -91,6 +84,8 @@ def __init__(
         self.degree_res = degree_res
         self.seed = seed
 
+        self.J_sync = JSync(src.n, self.epsilon, self.max_iters, self.seed)
+
     def _check_symmetry(self, symmetry):
         if symmetry is None:
             raise NotImplementedError(
@@ -116,10 +111,17 @@ def estimate_rotations(self):
         vijs, viis = self._global_J_sync(vijs, viis)
 
         logger.info("Estimating third rows of rotation matrices.")
-        vis = estimate_third_rows(vijs, viis)
+        vis = _estimate_third_rows(vijs, viis)
 
         logger.info("Estimating in-plane rotations and rotations matrices.")
-        Ris = estimate_inplane_rotations(self, vis)
+        Ris = _estimate_inplane_rotations(
+            vis,
+            self.pf,
+            self.max_shift,
+            self.shift_step,
+            self.order,
+            self.degree_res,
+        )
 
         self.rotations = Ris
 
@@ -145,7 +147,7 @@ def _estimate_relative_viewing_directions(self):
         Riis = self._estimate_all_Riis_c3_c4(sclmatrix)
 
         # Step 4: Calculate relative rotations
-        Rijs = self._estimate_all_Rijs_c3_c4(self.clmatrix)
+        Rijs = self._estimate_all_Rijs_c3_c4()
 
         # Step 5: Inner J-synchronization
         vijs, viis = self._local_J_sync_c3_c4(Rijs, Riis)
@@ -168,51 +170,13 @@ def _global_J_sync(self, vijs, viis):
 
         :return: vijs, viis all of which have a spurious J or not.
         """
-        n_img = self.n_img
 
         # Determine relative handedness of vijs.
-        sign_ij_J = self._J_sync_power_method(vijs)
-
-        # Synchronize vijs
-        for i, sign in enumerate(sign_ij_J):
-            if sign == -1:
-                vijs[i] = J_conjugate(vijs[i])
-
-        # Synchronize viis
-        # We use the fact that if v_ii and v_ij are of the same handedness, then v_ii @ v_ij = v_ij.
-        # If they are opposite handed then Jv_iiJ @ v_ij = v_ij. We compare each v_ii against all
-        # previously synchronized v_ij to get a consensus on the handedness of v_ii.
-        _, pairs_to_linear = all_pairs(n_img, return_map=True)
-        for i in range(n_img):
-            vii = viis[i]
-            vii_J = J_conjugate(vii)
-            J_consensus = 0
-            for j in range(n_img):
-                if j < i:
-                    idx = pairs_to_linear[j, i]
-                    vji = vijs[idx]
-
-                    err1 = norm(vji @ vii - vji)
-                    err2 = norm(vji @ vii_J - vji)
-
-                elif j > i:
-                    idx = pairs_to_linear[i, j]
-                    vij = vijs[idx]
-
-                    err1 = norm(vii @ vij - vij)
-                    err2 = norm(vii_J @ vij - vij)
-
-                else:
-                    continue
-
-                # Accumulate J consensus
-                if err1 < err2:
-                    J_consensus -= 1
-                else:
-                    J_consensus += 1
-
-            if J_consensus > 0:
-                viis[i] = vii_J
+        vijs = self.J_sync.global_J_sync(vijs)
+
+        # Determine relative handedness of viis, given synchronized vijs.
+        viis = self.J_sync.sync_viis(vijs, viis)
+
         return vijs, viis
 
     #################################################
@@ -266,8 +230,8 @@ def _self_clmatrix_c3_c4(self):
         # Compute the correlation over all shifts.
         # Generate Shifts.
         r_max = pf.shape[-1]
-        shifts, shift_phases, _ = self._generate_shift_phase_and_filter(
-            r_max, max_shift_1d, shift_step
+        shifts, shift_phases, _ = _generate_shift_phase_and_filter(
+            r_max, max_shift_1d, shift_step, self.dtype
         )
         n_shifts = len(shifts)
 
@@ -356,17 +320,21 @@ def _estimate_all_Riis_c3_c4(self, sclmatrix):
 
         return Riis
 
-    def _estimate_all_Rijs_c3_c4(self, clmatrix):
+    def _estimate_all_Rijs_c3_c4(self):
         """
         Estimate Rijs using the voting method.
         """
-        n_img = self.n_img
-        n_theta = self.n_theta
-        pairs = all_pairs(n_img)
+        pairs = all_pairs(self.n_img)
         Rijs = np.zeros((len(pairs), 3, 3))
         for idx, (i, j) in enumerate(pairs):
-            Rijs[idx] = self._syncmatrix_ij_vote_3n(
-                clmatrix, i, j, np.arange(n_img), n_theta
+            Rijs[idx] = _syncmatrix_ij_vote_3n(
+                self.clmatrix,
+                i,
+                j,
+                np.arange(self.n_img),
+                self.n_theta,
+                self.hist_bin_width,
+                self.full_width,
             )
 
         return Rijs
@@ -449,201 +417,3 @@ def _local_J_sync_c3_c4(self, Rijs, Riis):
             vijs[idx] = opts[min_idx]
 
         return vijs, viis
-
-    #######################################
-    # Secondary Methods for Global J Sync #
-    #######################################
-
-    def _J_sync_power_method(self, vijs):
-        """
-        Calculate the leading eigenvector of the J-synchronization matrix
-        using the power method.
-
-        As the J-synchronization matrix is of size (n-choose-2)x(n-choose-2), we
-        use the power method to compute the eigenvalues and eigenvectors,
-        while constructing the matrix on-the-fly.
-
-        :param vijs: (n-choose-2)x3x3 array of estimates of relative orientation matrices.
-
-        :return: An array of length n-choose-2 consisting of 1 or -1, where the sign of the
-        i'th entry indicates whether the i'th relative orientation matrix will be J-conjugated.
-        """
-
-        # Set power method tolerance and maximum iterations.
-        epsilon = self.epsilon
-        max_iters = self.max_iters
-
-        # Initialize candidate eigenvectors
-        n_vijs = vijs.shape[0]
-        vec = randn(n_vijs, seed=self.seed)
-        vec = vec / norm(vec)
-        residual = 1
-        itr = 0
-
-        # Power method iterations
-        logger.info(
-            "Initiating power method to estimate J-synchronization matrix eigenvector."
-        )
-        while itr < max_iters and residual > epsilon:
-            itr += 1
-            # Note, this appears to need double precision for accuracy in the following division.
-            vec_new = self._signs_times_v(vijs, vec).astype(np.float64, copy=False)
-            vec_new = vec_new / norm(vec_new)
-            residual = norm(vec_new - vec)
-            vec = vec_new
-            logger.info(
-                f"Iteration {itr}, residual {round(residual, 5)} (target {epsilon})"
-            )
-
-        # We need only the signs of the eigenvector
-        J_sync = np.sign(vec)
-
-        return J_sync
-
-    def _signs_times_v(self, vijs, vec):
-        """
-        Multiplication of the J-synchronization matrix by a candidate eigenvector.
-
-        The J-synchronization matrix is a matrix representation of the handedness graph, Gamma, whose set of
-        nodes consists of the estimates vijs and whose set of edges consists of the undirected edges between
-        all triplets of estimates vij, vjk, and vik, where i<j<k. The weight of an edge is set to +1 if its
-        incident nodes agree in handednes and -1 if not.
-
-        The J-synchronization matrix is of size (n-choose-2)x(n-choose-2), where each entry corresponds to
-        the relative handedness of vij and vjk. The entry (ij, jk), where ij and jk are retrieved from the
-        all_pairs indexing, is 1 if vij and vjk are of the same handedness and -1 if not. All other entries
-        (ij, kl) hold a zero.
-
-        Due to the large size of the J-synchronization matrix we construct it on the fly as follows.
-        For each triplet of outer products vij, vjk, and vik, the associated elements of the J-synchronization
-        matrix are populated with +1 or -1 and multiplied by the corresponding elements of
-        the current candidate eigenvector supplied by the power method. The new candidate eigenvector
-        is updated for each triplet.
-
-        :param vijs: (n-choose-2)x3x3 array, where each 3x3 slice holds the outer product of vi and vj.
-
-        :param vec: The current candidate eigenvector of length n-choose-2 from the power method.
-
-        :return: New candidate eigenvector of length n-choose-2. The product of the J-sync matrix and vec.
-        """
-
-        # All pairs (i,j) and triplets (i,j,k) where i<j<k
-        n_img = self.n_img
-        triplets = all_triplets(n_img)
-        pairs, pairs_to_linear = all_pairs(n_img, return_map=True)
-
-        # There are 4 possible configurations of relative handedness for each triplet (vij, vjk, vik).
-        # 'conjugate' expresses which node of the triplet must be conjugated (True) to achieve synchronization.
-        conjugate = np.empty((4, 3), bool)
-        conjugate[0] = [False, False, False]
-        conjugate[1] = [True, False, False]
-        conjugate[2] = [False, True, False]
-        conjugate[3] = [False, False, True]
-
-        # 'edges' corresponds to whether conjugation agrees between the pairs (vij, vjk), (vjk, vik),
-        # and (vik, vij). True if the pairs are in agreement, False otherwise.
-        edges = np.empty((4, 3), bool)
-        edges[:, 0] = conjugate[:, 0] == conjugate[:, 1]
-        edges[:, 1] = conjugate[:, 1] == conjugate[:, 2]
-        edges[:, 2] = conjugate[:, 2] == conjugate[:, 0]
-
-        # The corresponding entries in the J-synchronization matrix are +1 if the pair of nodes agree, -1 if not.
-        edge_signs = np.where(edges, 1, -1)
-
-        # For each triplet of nodes we apply the 4 configurations of conjugation and determine the
-        # relative handedness based on the condition that vij @ vjk - vik = 0 for synchronized nodes.
-        # We then construct the corresponding entries of the J-synchronization matrix with 'edge_signs'
-        # corresponding to the conjugation configuration producing the smallest residual for the above
-        # condition. Finally, we the multiply the 'edge_signs' by the cooresponding entries of 'vec'.
-        v = vijs
-        new_vec = np.zeros_like(vec)
-        for i, j, k in triplets:
-            ij = pairs_to_linear[i, j]
-            jk = pairs_to_linear[j, k]
-            ik = pairs_to_linear[i, k]
-            vij, vjk, vik = v[ij], v[jk], v[ik]
-            vij_J = J_conjugate(vij)
-            vjk_J = J_conjugate(vjk)
-            vik_J = J_conjugate(vik)
-
-            conjugated_pairs = np.where(
-                conjugate[..., np.newaxis, np.newaxis],
-                [vij_J, vjk_J, vik_J],
-                [vij, vjk, vik],
-            )
-            residual = np.stack([norm(x @ y - z) for x, y, z in conjugated_pairs])
-
-            min_residual = np.argmin(residual)
-
-            # Assign edge weights
-            s_ij_jk, s_ik_jk, s_ij_ik = edge_signs[min_residual]
-
-            # Update multiplication of signs times vec
-            new_vec[ij] += s_ij_jk * vec[jk] + s_ij_ik * vec[ik]
-            new_vec[jk] += s_ij_jk * vec[ij] + s_ik_jk * vec[ik]
-            new_vec[ik] += s_ij_ik * vec[ij] + s_ik_jk * vec[jk]
-
-        return new_vec
-
-    @staticmethod
-    def g_sync(rots, order, rots_gt):
-        """
-        Every estimated rotation might be a version of the ground truth rotation
-        rotated by g^{s_i}, where s_i = 0, 1, ..., order. This method synchronizes the
-        ground truth rotations so that only a single global rotation need be applied
-        to all estimates for error analysis.
-
-        :param rots: Estimated rotation matrices
-        :param order: The cyclic order asssociated with the symmetry of the underlying molecule.
-        :param rots_gt: Ground truth rotation matrices.
-
-        :return: g-synchronized ground truth rotations.
-        """
-        assert len(rots) == len(
-            rots_gt
-        ), "Number of estimates not equal to number of references."
-        n_img = len(rots)
-        dtype = rots.dtype
-
-        rots_symm = cyclic_rotations(order, dtype).matrices
-
-        A_g = np.zeros((n_img, n_img), dtype=complex)
-
-        pairs = all_pairs(n_img)
-
-        for i, j in pairs:
-            Ri = rots[i]
-            Rj = rots[j]
-            Rij = Ri.T @ Rj
-
-            Ri_gt = rots_gt[i]
-            Rj_gt = rots_gt[j]
-
-            diffs = np.zeros(order)
-            for s, g_s in enumerate(rots_symm):
-                Rij_gt = Ri_gt.T @ g_s @ Rj_gt
-                diffs[s] = min([norm(Rij - Rij_gt), norm(Rij - J_conjugate(Rij_gt))])
-
-            idx = np.argmin(diffs)
-
-            A_g[i, j] = np.exp(-1j * 2 * np.pi / order * idx)
-
-        # A_g(k,l) is exp(-j(-theta_k+theta_l))
-        # Diagonal elements correspond to exp(-i*0) so put 1.
-        # This is important only for verification purposes that spectrum is (K,0,0,0...,0).
-        A_g += np.conj(A_g).T + np.eye(n_img)
-
-        _, eig_vecs = eigh(A_g)
-        leading_eig_vec = eig_vecs[:, -1]
-
-        angles = np.exp(1j * 2 * np.pi / order * np.arange(order))
-        rots_gt_sync = np.zeros((n_img, 3, 3), dtype=dtype)
-
-        for i, rot_gt in enumerate(rots_gt):
-            # Since the closest ccw or cw rotation are just as good,
-            # we take the absolute value of the angle differences.
-            angle_dists = np.abs(np.angle(leading_eig_vec[i] / angles))
-            power_g_Ri = np.argmin(angle_dists)
-            rots_gt_sync[i] = rots_symm[power_g_Ri] @ rot_gt
-
-        return rots_gt_sync
diff --git a/src/aspire/abinitio/commonline_cn.py b/src/aspire/abinitio/commonline_cn.py
index 297bf04a6d..9adddb37ac 100644
--- a/src/aspire/abinitio/commonline_cn.py
+++ b/src/aspire/abinitio/commonline_cn.py
@@ -3,7 +3,7 @@
 import numpy as np
 from numpy.linalg import norm
 
-from aspire.abinitio import CLSymmetryC3C4, cl_angles_to_ind, complete_third_row_to_rot
+from aspire.abinitio import CLOrient3D, JSync
 from aspire.operators import PolarFT
 from aspire.utils import (
     J_conjugate,
@@ -17,10 +17,18 @@
 )
 from aspire.utils.random import Random, randn
 
+from .commonline_utils import (
+    _cl_angles_to_ind,
+    _complete_third_row_to_rot,
+    _estimate_inplane_rotations,
+    _estimate_third_rows,
+    _generate_shift_phase_and_filter,
+)
+
 logger = logging.getLogger(__name__)
 
 
-class CLSymmetryCn(CLSymmetryC3C4):
+class CLSymmetryCn(CLOrient3D):
     """
     Define a class to estimate 3D orientations using common lines methods for molecules with
     Cn cyclic symmetry, with n>4.
@@ -65,22 +73,24 @@ def __init__(
 
         super().__init__(
             src,
-            symmetry=symmetry,
             n_rad=n_rad,
             n_theta=n_theta,
             max_shift=max_shift,
             shift_step=shift_step,
-            epsilon=epsilon,
-            max_iters=max_iters,
-            degree_res=degree_res,
-            seed=seed,
             mask=mask,
             **kwargs,
         )
 
+        self._check_symmetry(symmetry)
+        self.epsilon = epsilon
+        self.max_iters = max_iters
+        self.degree_res = degree_res
+        self.seed = seed
         self.n_points_sphere = n_points_sphere
         self.equator_threshold = equator_threshold
 
+        self.J_sync = JSync(src.n, self.epsilon, self.max_iters, self.seed)
+
     def _check_symmetry(self, symmetry):
         if symmetry is None:
             raise NotImplementedError(
@@ -102,7 +112,27 @@ def estimate_rotations(self):
 
         :return: Array of rotation matrices, size n_imgx3x3.
         """
-        super().estimate_rotations()
+        vijs, viis = self._estimate_relative_viewing_directions()
+
+        logger.info("Performing global handedness synchronization.")
+        vijs, viis = self._global_J_sync(vijs, viis)
+
+        logger.info("Estimating third rows of rotation matrices.")
+        vis = _estimate_third_rows(vijs, viis)
+
+        logger.info("Estimating in-plane rotations and rotations matrices.")
+        Ris = _estimate_inplane_rotations(
+            vis,
+            self.pf,
+            self.max_shift,
+            self.shift_step,
+            self.order,
+            self.degree_res,
+        )
+
+        self.rotations = Ris
+
+        return self.rotations
 
     def _estimate_relative_viewing_directions(self):
         logger.info(f"Estimating relative viewing directions for {self.n_img} images.")
@@ -123,8 +153,8 @@ def _estimate_relative_viewing_directions(self):
 
         # Generate shift phases.
         r_max = pf.shape[-1]
-        shifts, shift_phases, _ = self._generate_shift_phase_and_filter(
-            r_max, self.max_shift, self.shift_step
+        shifts, shift_phases, _ = _generate_shift_phase_and_filter(
+            r_max, self.max_shift, self.shift_step, self.dtype
         )
         n_shifts = len(shifts)
 
@@ -286,6 +316,31 @@ def _compute_cls_inds(self, Ris_tilde, R_theta_ijs):
                 cij_inds[i, j, :, 1] = c2s
         return cij_inds
 
+    def _global_J_sync(self, vijs, viis):
+        """
+        Global J-synchronization of all third row outer products. Given 3x3 matrices vijs and viis, each
+        of which might contain a spurious J (ie. vij = J*vi*vj^T*J instead of vij = vi*vj^T),
+        we return vijs and viis that all have either a spurious J or not.
+
+        :param vijs: An (n-choose-2)x3x3 array where each 3x3 slice holds an estimate for the corresponding
+        outer-product vi*vj^T between the third rows of the rotation matrices Ri and Rj. Each estimate
+        might have a spurious J independently of other estimates.
+
+        :param viis: An n_imgx3x3 array where the i'th slice holds an estimate for the outer product vi*vi^T
+        between the third row of matrix Ri and itself. Each estimate might have a spurious J independently
+        of other estimates.
+
+        :return: vijs, viis all of which have a spurious J or not.
+        """
+
+        # Determine relative handedness of vijs.
+        vijs = self.J_sync.global_J_sync(vijs)
+
+        # Determine relative handedness of viis, given synchronized vijs.
+        viis = self.J_sync.sync_viis(vijs, viis)
+
+        return vijs, viis
+
     @staticmethod
     def relative_rots_to_cl_indices(relative_rots, n_theta):
         """
@@ -300,8 +355,8 @@ def relative_rots_to_cl_indices(relative_rots, n_theta):
         c1s = np.array((-relative_rots[:, 1, 2], relative_rots[:, 0, 2])).T
         c2s = np.array((relative_rots[:, 2, 1], -relative_rots[:, 2, 0])).T
 
-        c1s = cl_angles_to_ind(c1s, n_theta)
-        c2s = cl_angles_to_ind(c2s, n_theta)
+        c1s = _cl_angles_to_ind(c1s, n_theta)
+        c2s = _cl_angles_to_ind(c2s, n_theta)
 
         inds = np.where(c1s >= n_theta // 2)
         c1s[inds] -= n_theta // 2
@@ -333,7 +388,7 @@ def generate_candidate_rots(n, equator_threshold, order, degree_res, seed):
             while counter < n:
                 third_row = randn(3)
                 third_row /= anorm(third_row, axes=(-1,))
-                Ri_tilde = complete_third_row_to_rot(third_row)
+                Ri_tilde = _complete_third_row_to_rot(third_row)
 
                 # Exclude candidates that represent equator images. Equator candidates
                 # induce collinear self-common-lines, which always have perfect correlation.
diff --git a/src/aspire/abinitio/commonline_d2.py b/src/aspire/abinitio/commonline_d2.py
index f8022d3db9..e1730bf35e 100644
--- a/src/aspire/abinitio/commonline_d2.py
+++ b/src/aspire/abinitio/commonline_d2.py
@@ -10,6 +10,8 @@
 from aspire.utils.random import randn
 from aspire.volume import DnSymmetryGroup
 
+from .commonline_utils import _generate_shift_phase_and_filter
+
 logger = logging.getLogger(__name__)
 
 
@@ -131,8 +133,8 @@ def _compute_shifted_pf(self):
         # Generate shift phases.
         r_max = pf.shape[-1]
         max_shift_1d = np.ceil(2 * np.sqrt(2) * self.max_shift)
-        shifts, shift_phases, _ = self._generate_shift_phase_and_filter(
-            r_max, max_shift_1d, self.shift_step
+        shifts, shift_phases, _ = _generate_shift_phase_and_filter(
+            r_max, max_shift_1d, self.shift_step, self.dtype
         )
         self.n_shifts = len(shifts)
 
diff --git a/src/aspire/abinitio/commonline_sync.py b/src/aspire/abinitio/commonline_sync.py
index 6fcc8daddc..6ed774958b 100644
--- a/src/aspire/abinitio/commonline_sync.py
+++ b/src/aspire/abinitio/commonline_sync.py
@@ -2,14 +2,15 @@
 
 import numpy as np
 
-from aspire.abinitio import CLOrient3D, SyncVotingMixin
+from aspire.abinitio import CLOrient3D
+from aspire.abinitio.sync_voting import _rotratio_eulerangle_vec, _vote_ij
 from aspire.utils import nearest_rotations
 from aspire.utils.matlab_compat import stable_eigsh
 
 logger = logging.getLogger(__name__)
 
 
-class CLSyncVoting(CLOrient3D, SyncVotingMixin):
+class CLSyncVoting(CLOrient3D):
     """
     Define a class to estimate 3D orientations using synchronization matrix and voting method.
 
@@ -201,9 +202,11 @@ def _syncmatrix_ij_vote(self, clmatrix, i, j, k_list, n_theta):
         :return: The (i,j) rotation block of the synchronization matrix
         """
 
-        _, good_k = self._vote_ij(clmatrix, n_theta, i, j, k_list)
+        _, good_k = _vote_ij(
+            clmatrix, n_theta, i, j, k_list, self.hist_bin_width, self.full_width
+        )
 
-        rots = self._rotratio_eulerangle_vec(clmatrix, i, j, good_k, n_theta)
+        rots = _rotratio_eulerangle_vec(clmatrix, i, j, good_k, n_theta)
 
         if rots is not None:
             rot_mean = np.mean(rots, 0)
diff --git a/src/aspire/abinitio/commonline_sync3n.py b/src/aspire/abinitio/commonline_sync3n.py
index ed7ca94048..7be4b37c2c 100644
--- a/src/aspire/abinitio/commonline_sync3n.py
+++ b/src/aspire/abinitio/commonline_sync3n.py
@@ -6,14 +6,15 @@
 from numpy.linalg import norm
 from scipy.optimize import curve_fit
 
-from aspire.abinitio import CLOrient3D, SyncVotingMixin
+from aspire.abinitio import CLOrient3D
+from aspire.abinitio.sync_voting import _syncmatrix_ij_vote_3n
 from aspire.utils import J_conjugate, all_pairs, nearest_rotations, random, tqdm, trange
 from aspire.utils.matlab_compat import stable_eigsh
 
 logger = logging.getLogger(__name__)
 
 
-class CLSync3N(CLOrient3D, SyncVotingMixin):
+class CLSync3N(CLOrient3D):
     """
     Define a class to estimate 3D orientations using common lines Sync3N methods (2017).
 
@@ -957,8 +958,14 @@ def _estimate_all_Rijs_host(self, clmatrix):
         Rijs = np.zeros((len(self._pairs), 3, 3))
 
         for idx, (i, j) in enumerate(tqdm(self._pairs, desc="Estimate Rijs")):
-            Rijs[idx] = self._syncmatrix_ij_vote_3n(
-                clmatrix, i, j, np.arange(n_img), n_theta
+            Rijs[idx] = _syncmatrix_ij_vote_3n(
+                clmatrix,
+                i,
+                j,
+                np.arange(n_img),
+                n_theta,
+                self.hist_bin_width,
+                self.full_width,
             )
 
         return Rijs
diff --git a/src/aspire/abinitio/commonline_utils.py b/src/aspire/abinitio/commonline_utils.py
index bfd1e5ce5f..45352bf026 100644
--- a/src/aspire/abinitio/commonline_utils.py
+++ b/src/aspire/abinitio/commonline_utils.py
@@ -4,12 +4,12 @@
 from numpy.linalg import eigh, norm
 
 from aspire.operators import PolarFT
-from aspire.utils import Rotation, all_pairs, anorm, tqdm
+from aspire.utils import J_conjugate, Rotation, all_pairs, anorm, cyclic_rotations, tqdm
 
 logger = logging.getLogger(__name__)
 
 
-def estimate_third_rows(vijs, viis):
+def _estimate_third_rows(vijs, viis):
     """
     Find the third row of each rotation matrix given a collection of matrices
     representing the outer products of the third rows from each rotation matrix.
@@ -59,39 +59,69 @@ def estimate_third_rows(vijs, viis):
     return vis
 
 
-def estimate_inplane_rotations(cl_class, vis):
+def _generate_shift_phase_and_filter(r_max, max_shift, shift_step, dtype):
     """
-    Estimate the rotation matrices for each image by constructing arbitrary rotation matrices
-    populated with the given third rows, vis, and then rotating by an appropriate in-plane rotation.
+    Prepare the shift phases and generate filter for common-line detection
 
-    :cl_class: A commonlines class instance.
-    :param vis: An n_imgx3 array where the i'th row holds the estimate for the third row of
-    the i'th rotation matrix.
+    The shift phases are pre-defined in a range of max_shift that can be
+    applied to maximize the common line calculation. The common-line filter
+    is also applied to the radial direction for easier detection.
+
+    :param r_max: Maximum index for common line detection.
+    :param max_shift: Maximum value of 1D shift (in pixels) to search.
+    :param shift_step: Resolution of shift estimation in pixels.
+    :param dtype: dtype for shift phases and filter.
+    :return: shift phases matrix and common lines filter.
+    """
+
+    # Number of shifts to try
+    n_shifts = int(np.ceil(2 * max_shift / shift_step + 1))
+
+    # only half of ray, excluding the DC component.
+    rk = np.arange(1, r_max + 1, dtype=dtype)
+
+    # Generate all shift phases
+    shifts = -max_shift + shift_step * np.arange(n_shifts, dtype=dtype)
+    shift_phases = np.exp(np.outer(shifts, -2 * np.pi * 1j * rk / (2 * r_max + 1)))
+    # Set filter for common-line detection
+    h = np.sqrt(np.abs(rk)) * np.exp(-np.square(rk) / (2 * (r_max / 4) ** 2))
+
+    return shifts, shift_phases, h
 
+
+def _estimate_inplane_rotations(vis, pf, max_shift, shift_step, order, degree_res):
+    """
+    Estimate the rotation matrices for each image of a cyclically symmetric molecule by
+    constructing arbitrary rotation matrices populated with the given third rows, vis, and
+    then rotating by an appropriate in-plane rotation.
+
+    :param vis: An n_imgx3 array where the i'th row holds the estimate for the third row of
+        the i'th rotation matrix.
+    :param pf: The polar Fourier transform of the source images, shape (n_img, n_theta/2, n_rad).
+    :param max_shift: Maximum range for shifts (in pixels) for estimating in-plane rotations.
+    :param shift_step: Shift step (in pixels) for estimating in-plane rotations.
+    :param order: Cyclic order.
+    :param degree_res: Resolution (in degrees) of in-plane rotation to search over.
     :return: Rotation matrices Ris and in-plane rotation matrices R_thetas, both size n_imgx3x3.
     """
-    pf = cl_class.pf
-    n_img = cl_class.n_img
-    n_theta = cl_class.n_theta
-    max_shift_1d = cl_class.max_shift
-    shift_step = cl_class.shift_step
-    order = cl_class.order
-    degree_res = cl_class.degree_res
+    n_img = vis.shape[0]
+    dtype = vis.dtype
+    n_theta = pf.shape[1] * 2
 
     # Step 1: Construct all rotation matrices Ri_tildes whose third rows are equal to
     # the corresponding third rows vis.
-    Ri_tildes = complete_third_row_to_rot(vis)
+    Ri_tildes = _complete_third_row_to_rot(vis)
 
     # Step 2: Construct all in-plane rotation matrices, R_theta_ijs.
     max_angle = (360 // order) * order
     theta_ijs = np.arange(0, max_angle, degree_res) * np.pi / 180
-    R_theta_ijs = Rotation.about_axis("z", theta_ijs, dtype=cl_class.dtype).matrices
+    R_theta_ijs = Rotation.about_axis("z", theta_ijs, dtype=dtype).matrices
 
     # Step 3: Compute the correlation over all shifts.
     # Generate shifts.
     r_max = pf.shape[-1]
-    shifts, shift_phases, _ = cl_class._generate_shift_phase_and_filter(
-        r_max, max_shift_1d, shift_step
+    shifts, shift_phases, _ = _generate_shift_phase_and_filter(
+        r_max, max_shift, shift_step, dtype
     )
     n_shifts = len(shifts)
 
@@ -100,7 +130,7 @@ def estimate_inplane_rotations(cl_class, vis):
     # and theta_i in [0, 2pi/order) is the in-plane rotation angle for the i'th image.
     Q = np.zeros((n_img, n_img), dtype=complex)
 
-    # Reconstruct the full polar Fourier for use in correlation. cl_class.pf only consists of
+    # Reconstruct the full polar Fourier for use in correlation. pf only consists of
     # rays in the range [180, 360), with shape (n_img, n_theta//2, n_rad-1).
     pf = PolarFT.half_to_full(pf)
 
@@ -108,87 +138,83 @@ def estimate_inplane_rotations(cl_class, vis):
     pf /= norm(pf, axis=-1)[..., np.newaxis]
 
     n_pairs = n_img * (n_img - 1) // 2
-    with tqdm(total=n_pairs) as pbar:
-        idx = 0
-        # Note: the ordering of i and j in these loops should not be changed as
-        # they correspond to the ordered tuples (i, j), for i<j.
-        for i in range(n_img):
-            pf_i = pf[i]
-
-            # Generate shifted versions of images.
-            pf_i_shifted = np.array(
-                [pf_i * shift_phase for shift_phase in shift_phases]
-            )
-
-            Ri_tilde = Ri_tildes[i]
-
-            for j in range(i + 1, n_img):
-                pf_j = pf[j]
+    pbar = tqdm(total=n_pairs)
+    # Note: the ordering of i and j in these loops should not be changed as
+    # they correspond to the ordered tuples (i, j), for i<j.
+    for i in range(n_img):
+        pf_i = pf[i]
 
-                Rj_tilde = Ri_tildes[j]
+        # Generate shifted versions of images.
+        pf_i_shifted = np.array([pf_i * shift_phase for shift_phase in shift_phases])
 
-                # Compute all possible rotations between the i'th and j'th images.
-                Us = np.array(
-                    [Ri_tilde.T @ R_theta_ij @ Rj_tilde for R_theta_ij in R_theta_ijs]
-                )
+        Ri_tilde = Ri_tildes[i]
 
-                # Find the angle between common lines induced by the rotations.
-                c1s = np.array([[-U[1, 2], U[0, 2]] for U in Us])
-                c2s = np.array([[U[2, 1], -U[2, 0]] for U in Us])
+        for j in range(i + 1, n_img):
+            pf_j = pf[j]
 
-                # Convert from angles to indices.
-                c1s = cl_angles_to_ind(c1s, n_theta)
-                c2s = cl_angles_to_ind(c2s, n_theta)
+            Rj_tilde = Ri_tildes[j]
 
-                # Perform correlation, corrs is shape n_shifts x len(theta_ijs).
-                corrs = np.array(
-                    [
-                        np.dot(pf_i_shift[c1], np.conj(pf_j[c2]))
-                        for pf_i_shift in pf_i_shifted
-                        for c1, c2 in zip(c1s, c2s)
-                    ]
-                )
+            # Compute all possible rotations between the i'th and j'th images.
+            Us = np.array(
+                [Ri_tilde.T @ R_theta_ij @ Rj_tilde for R_theta_ij in R_theta_ijs]
+            )
 
-                # Reshape to group by shift and symmetric order.
-                corrs = corrs.reshape((n_shifts, order, len(theta_ijs) // order))
+            # Find the angle between common lines induced by the rotations.
+            c1s = np.array([[-U[1, 2], U[0, 2]] for U in Us])
+            c2s = np.array([[U[2, 1], -U[2, 0]] for U in Us])
+
+            # Convert from angles to indices.
+            c1s = _cl_angles_to_ind(c1s, n_theta)
+            c2s = _cl_angles_to_ind(c2s, n_theta)
+
+            # Perform correlation, corrs is shape n_shifts x len(theta_ijs).
+            corrs = np.array(
+                [
+                    np.dot(pf_i_shift[c1], np.conj(pf_j[c2]))
+                    for pf_i_shift in pf_i_shifted
+                    for c1, c2 in zip(c1s, c2s)
+                ]
+            )
 
-                # For each pair of lines we take the maximum correlation over all shifts.
-                corrs = np.max(np.real(corrs), axis=0)
+            # Reshape to group by shift and symmetric order.
+            corrs = corrs.reshape((n_shifts, order, len(theta_ijs) // order))
 
-                # corrs[i] is the set of correlations for theta_ij in [2pi * i / order, 2pi * (i + 1) / order).
-                # Due to symmetry, each corrs[i] represents correlations over identical pairs of lines.
-                # With that in mind, we average over corrs[i] and find the max correlation.
-                # This produces an index corresponding to theta_ij in the range [0, 2pi/order).
-                corrs = np.mean(np.real(corrs), axis=0)
-                max_idx_corr = np.argmax(corrs)
+            # For each pair of lines we take the maximum correlation over all shifts.
+            corrs = np.max(np.real(corrs), axis=0)
 
-                theta_ij = degree_res * max_idx_corr * np.pi / 180
+            # corrs[i] is the set of correlations for theta_ij in [2pi * i / order, 2pi * (i + 1) / order).
+            # Due to symmetry, each corrs[i] represents correlations over identical pairs of lines.
+            # With that in mind, we average over corrs[i] and find the max correlation.
+            # This produces an index corresponding to theta_ij in the range [0, 2pi/order).
+            corrs = np.mean(np.real(corrs), axis=0)
+            max_idx_corr = np.argmax(corrs)
 
-                Q[i, j] = np.exp(-1j * order * theta_ij)
+            theta_ij = degree_res * max_idx_corr * np.pi / 180
 
-                pbar.update()
+            Q[i, j] = np.exp(-1j * order * theta_ij)
 
-                idx += 1
+            pbar.update()
+    pbar.close()
 
-        # Populate the lower triangle and diagonal of Q.
-        # Diagonals are 1 since e^{i*0}=1.
-        Q += np.conj(Q).T + np.eye(n_img)
+    # Populate the lower triangle and diagonal of Q.
+    # Diagonals are 1 since e^{i*0}=1.
+    Q += np.conj(Q).T + np.eye(n_img)
 
-        # Q is a rank-1 Hermitian matrix.
-        eig_vals, eig_vecs = eigh(Q)
-        leading_eig_vec = eig_vecs[:, -1]
-        logger.info(f"Top 3 eigenvalues of Q (rank-1) are {str(eig_vals[-3:][::-1])}.")
+    # Q is a rank-1 Hermitian matrix.
+    eig_vals, eig_vecs = eigh(Q)
+    leading_eig_vec = eig_vecs[:, -1]
+    logger.info(f"Top 3 eigenvalues of Q (rank-1) are {str(eig_vals[-3:][::-1])}.")
 
-        # Calculate R_thetas.
-        R_thetas = Rotation.about_axis("z", np.angle(leading_eig_vec ** (1 / order)))
+    # Calculate R_thetas.
+    R_thetas = Rotation.about_axis("z", np.angle(leading_eig_vec ** (1 / order)))
 
-        # Form rotation matrices Ris.
-        Ris = R_thetas @ Ri_tildes
+    # Form rotation matrices Ris.
+    Ris = R_thetas @ Ri_tildes
 
-        return Ris
+    return Ris
 
 
-def complete_third_row_to_rot(r3):
+def _complete_third_row_to_rot(r3):
     """
     Construct rotation matrices whose third rows are equal to the given row vectors.
     For vector r3 = [a, b, c], where [a, b, c] != [0, 0, 1], we return the matrix
@@ -237,7 +263,19 @@ def complete_third_row_to_rot(r3):
     return rots
 
 
-def cl_angles_to_ind(cl_angles, n_theta):
+def _cl_angles_to_ind(cl_angles, n_theta):
+    """
+    Map 2D direction vectors to discretized angular indices.
+
+    For each 2D vector [x, y] in `cl_angles`, compute its polar angle
+    and find the nearest of `n_theta` polar ray indices.
+
+    :param cl_angles: Array of shape (n, 2) of [x, y] values corresponding
+        to the commonline induced by a pair of rotations.
+    :param n_theta: Resolution of polar rays.
+
+    :return: int or array of length n of commonline indice in the range [0, n_theta - 1].
+    """
     thetas = np.arctan2(cl_angles[:, 1], cl_angles[:, 0])
 
     # Shift from [-pi,pi] to [0,2*pi).
@@ -251,3 +289,104 @@ def cl_angles_to_ind(cl_angles, n_theta):
         ind = ind.flat[0]
 
     return ind
+
+
+def g_sync(rots, order, rots_gt):
+    """
+    Given ground truth rotations, synchronize estimated rotations over
+    symmetry group elements.
+
+    Every estimated rotation might be a version of the ground truth rotation
+    rotated by g^{s_i}, where s_i = 0, 1, ..., order. This method synchronizes the
+    ground truth rotations so that only a single global rotation need be applied
+    to all estimates for error analysis.
+
+    :param rots: Estimated rotation matrices
+    :param order: The cyclic order asssociated with the symmetry of the underlying molecule.
+    :param rots_gt: Ground truth rotation matrices.
+
+    :return: g-synchronized ground truth rotations.
+    """
+    assert len(rots) == len(
+        rots_gt
+    ), "Number of estimates not equal to number of references."
+    n_img = len(rots)
+    dtype = rots.dtype
+
+    rots_symm = cyclic_rotations(order, dtype).matrices
+
+    A_g = np.zeros((n_img, n_img), dtype=complex)
+
+    pairs = all_pairs(n_img)
+
+    for i, j in pairs:
+        Ri = rots[i]
+        Rj = rots[j]
+        Rij = Ri.T @ Rj
+
+        Ri_gt = rots_gt[i]
+        Rj_gt = rots_gt[j]
+
+        diffs = np.zeros(order)
+        for s, g_s in enumerate(rots_symm):
+            Rij_gt = Ri_gt.T @ g_s @ Rj_gt
+            diffs[s] = min([norm(Rij - Rij_gt), norm(Rij - J_conjugate(Rij_gt))])
+
+        idx = np.argmin(diffs)
+
+        A_g[i, j] = np.exp(-1j * 2 * np.pi / order * idx)
+
+    # A_g(k,l) is exp(-j(-theta_k+theta_l))
+    # Diagonal elements correspond to exp(-i*0) so put 1.
+    # This is important only for verification purposes that spectrum is (K,0,0,0...,0).
+    A_g += np.conj(A_g).T + np.eye(n_img)
+
+    _, eig_vecs = eigh(A_g)
+    leading_eig_vec = eig_vecs[:, -1]
+
+    angles = np.exp(1j * 2 * np.pi / order * np.arange(order))
+    rots_gt_sync = np.zeros((n_img, 3, 3), dtype=dtype)
+
+    for i, rot_gt in enumerate(rots_gt):
+        # Since the closest ccw or cw rotation are just as good,
+        # we take the absolute value of the angle differences.
+        angle_dists = np.abs(np.angle(leading_eig_vec[i] / angles))
+        power_g_Ri = np.argmin(angle_dists)
+        rots_gt_sync[i] = rots_symm[power_g_Ri] @ rot_gt
+
+    return rots_gt_sync
+
+
+def build_outer_products(n, dtype):
+    """
+    Builds sets of outer products of 3rd rows of rotation matrices.
+    This is a helper function used in commonline testing.
+
+    :param n: Number of 3rd rows to construct outer product from.
+    :param dtype: dtype of outputs
+
+    :return: tuple of (vijs, viis, gt_vis), where vijs are the pairwise
+        outer products of gt_vis and viis are self outer products of gt_vis.
+    """
+    # Build random third rows, ground truth vis (unit vectors)
+    gt_vis = np.zeros((n, 3), dtype=dtype)
+    for i in range(n):
+        np.random.seed(i)
+        v = np.random.randn(3)
+        gt_vis[i] = v / norm(v)
+
+    # Find outer products viis and vijs for i<j
+    nchoose2 = int(n * (n - 1) / 2)
+    vijs = np.zeros((nchoose2, 3, 3), dtype=dtype)
+    viis = np.zeros((n, 3, 3), dtype=dtype)
+
+    # All pairs (i,j) where i<j
+    pairs = all_pairs(n)
+
+    for k, (i, j) in enumerate(pairs):
+        vijs[k] = np.outer(gt_vis[i], gt_vis[j])
+
+    for i in range(n):
+        viis[i] = np.outer(gt_vis[i], gt_vis[i])
+
+    return vijs, viis, gt_vis
diff --git a/src/aspire/abinitio/sync_voting.py b/src/aspire/abinitio/sync_voting.py
index b866f1e8ff..f7dad8ff11 100644
--- a/src/aspire/abinitio/sync_voting.py
+++ b/src/aspire/abinitio/sync_voting.py
@@ -7,310 +7,326 @@
 logger = logging.getLogger(__name__)
 
 
-class SyncVotingMixin(object):
+def _syncmatrix_ij_vote_3n(
+    clmatrix, i, j, k_list, n_theta, hist_bin_width, full_width, sigma=3.0
+):
     """
-    SyncVotingMixin is a mixin implementing methods for the synchronization voting algorithm
-    which are shared by CLSynVoting and CLSymmetryC3C4
+    Compute the (i,j) rotation block of the synchronization matrix using voting method
+
+    Given the common lines matrix `clmatrix`, a list of images specified in k_list
+    and the number of common lines n_theta, find the (i, j) rotation block Rij.
+
+    :param clmatrix: The common lines matrix
+    :param i: The i image
+    :param j: The j image
+    :param k_list: The list of images for the third image for voting algorithm
+    :param n_theta: The number of points in the theta direction (common lines)
+    :param hist_bin_width: Bin width in smoothing histogram (degrees).
+    :param full_width: Selection width around smoothed histogram peak (degrees).
+        `adaptive` will attempt to automatically find the smallest number of
+        `hist_bin_width`s required to find at least one valid image index.
+    :param sigma: Voting contribution smoothing factor. Default is 3.0.
+
+    :return: The (i,j) rotation block of the synchronization matrix
     """
+    alphas, good_k = _vote_ij(
+        clmatrix,
+        n_theta,
+        i,
+        j,
+        k_list,
+        hist_bin_width,
+        full_width,
+        sigma=sigma,
+        sync=True,
+    )
+
+    angles = np.zeros(3)
+
+    if alphas is not None:
+        angles[0] = clmatrix[i, j] * 2 * np.pi / n_theta + np.pi / 2
+        angles[1] = np.mean(alphas)
+        angles[2] = -np.pi / 2 - clmatrix[j, i] * 2 * np.pi / n_theta
+        rot = Rotation.from_euler(angles).matrices
+
+    else:
+        # This is for the case that images i and j correspond to the same
+        # viewing direction and differ only by in-plane rotation.
+        # We set to zero as in the Matlab code.
+        rot = np.zeros((3, 3))
+
+    return rot
+
+
+def _rotratio_eulerangle_vec(clmatrix, i, j, good_k, n_theta):
+    """
+    Compute the rotation that takes image i to image j
 
-    def _syncmatrix_ij_vote_3n(self, clmatrix, i, j, k_list, n_theta):
-        """
-        Compute the (i,j) rotation block of the synchronization matrix using voting method
-
-        Given the common lines matrix `clmatrix`, a list of images specified in k_list
-        and the number of common lines n_theta, find the (i, j) rotation block Rij.
-
-        :param clmatrix: The common lines matrix
-        :param i: The i image
-        :param j: The j image
-        :param k_list: The list of images for the third image for voting algorithm
-        :param n_theta: The number of points in the theta direction (common lines)
-        :return: The (i,j) rotation block of the synchronization matrix
-        """
-        alphas, good_k = self._vote_ij(clmatrix, n_theta, i, j, k_list, sync=True)
-
-        angles = np.zeros(3)
-
-        if alphas is not None:
-            angles[0] = clmatrix[i, j] * 2 * np.pi / n_theta + np.pi / 2
-            angles[1] = np.mean(alphas)
-            angles[2] = -np.pi / 2 - clmatrix[j, i] * 2 * np.pi / n_theta
-            rot = Rotation.from_euler(angles).matrices
-
-        else:
-            # This is for the case that images i and j correspond to the same
-            # viewing direction and differ only by in-plane rotation.
-            # We set to zero as in the Matlab code.
-            rot = np.zeros((3, 3))
-
-        return rot
-
-    def _rotratio_eulerangle_vec(self, clmatrix, i, j, good_k, n_theta):
-        """
-        Compute the rotation that takes image i to image j
-
-        Given a common lines matrix, where the index of each common line
-        is in the range of n_theta and a list of good image k from voting results.
-
-        :param clmatrix: The common lines matrix
-        :param i: The i image
-        :param j: The j image
-        :param good_k: The list of good images k from voting algorithm
-        :param n_theta: The number of points in the theta direction (common lines)
-        :return: The rotation matrix that takes image i to image j for good index of k.
-        """
-
-        if i == j:
-            return []
-
-        # Prepare the theta values from the differences of common line indices
-        # C1, C2, and C3 are unit circles of image i, j, and k
-        # cl_diff1 is for the angle on C1 created by its intersection with C3 and C2.
-        # 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_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, sync=False
-        )
+    Given a common lines matrix, where the index of each common line
+    is in the range of n_theta and a list of good image k from voting results.
 
-        if len(c_alpha) == 0:
-            return None
-        alpha = np.arccos(c_alpha)
-
-        # Convert the Euler angles with ZYZ conversion to rotation matrices
-        angles = np.zeros((alpha.shape[0], 3))
-        angles[:, 0] = clmatrix[i, j] * 2 * np.pi / n_theta + np.pi / 2
-        angles[:, 1] = alpha
-        angles[:, 2] = -np.pi / 2 - clmatrix[j, i] * 2 * np.pi / n_theta
-        r = Rotation.from_euler(angles).matrices
-
-        return r[good_idx, :, :]
-
-    def _vote_ij(self, clmatrix, n_theta, i, j, k_list, sync=False):
-        """
-        Apply the voting algorithm for images i and j.
-
-        clmatrix is the common lines matrix, constructed using angular resolution,
-        n_theta. k_list are the images to be used for voting of the pair of images
-        (i ,j).
-
-        :param clmatrix: The common lines matrix
-        :param n_theta: The number of points in the theta direction (common lines)
-        :param i: The i image
-        :param j: The j image
-        :param k_list: The list of images for the third image for voting algorithm
-        :param sync: Perform 180 degree ambiguity synchronization.
-        :return: (alpha, good_k), angles and list of all third images
-            in the peak of the histogram corresponding to the pair of
-            images (i,j)
-        """
-
-        if i == j or clmatrix[i, j] == -1:
-            return None, []
-
-        # Some of the entries in clmatrix may be zero if we cleared
-        # them due to small correlation, or if for each image
-        # we compute intersections with only some of the other images.
-        #
-        # Note that as long as the diagonal of the common lines matrix is
-        # -1, the conditions (i != j) && (j != k) are not needed, since
-        # if i == j then clmatrix[i, k] == -1 and similarly for i == k or
-        # j == k. Thus, the previous voting code (from the JSB paper) is
-        # correct even though it seems that we should test also that
-        # (i != j) && (i != k) && (j != k), and only (i != j) && (i != k)
-        #  as tested there.
-        cl_idx12 = clmatrix[i, j]
-        cl_idx21 = clmatrix[j, i]
-        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]
-        cl_idx32 = clmatrix[k_list, j]
-
-        # Prepare the theta values from the differences of common line indices
-        # C1, C2, and C3 are unit circles of image i, j, and k
-        # cl_diff1 is for the angle on C1 created by its intersection with C3 and C2.
-        # 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 = cl_idx13 - cl_idx12
-        cl_diff2 = cl_idx23 - cl_idx21
-        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, sync=sync
-        )
+    :param clmatrix: The common lines matrix
+    :param i: The i image
+    :param j: The j image
+    :param good_k: The list of good images k from voting algorithm
+    :param n_theta: The number of points in the theta direction (common lines)
 
-        if np.any(np.abs(cos_phi2) - 1 > 1e-12):
-            logger.warning(
-                f"Globally Consistent Angular Reconstruction (GCAR) exists"
-                f" numerical problem: abs(cos_phi2) > 1, with the"
-                f" difference of {np.abs(cos_phi2)-1}."
-            )
-        cos_phi2 = np.clip(cos_phi2, -1, 1)
-
-        # Store angles between i and j induced by each third image k.
-        phis = cos_phi2
-        # Sore good indices of l in k_list of the image that creates that angle.
-        inds = k_list[good_idx]
-
-        if phis.shape[0] == 0:
-            return None, []
-
-        # Parameters used to compute the smoothed angle histogram.
-        ntics = int(180 / self.hist_bin_width)
-        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 = getattr(self, "sigma", 3.0)  # get from class if avail
-
-        # Compute the histogram of the angles between images i and j
-        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
-        # angle between i and j that is at most 10 off the true angle.
-        # Even for debugging, don't put a value that is smaller than two
-        # tics, since the peak might move a little bit due to wrong k images
-        # that accidentally fall near the peak.
-        peak_idx = angles_hist.argmax()
-
-        if self.full_width == -1:
-            # Adaptive width  (MATLAB)
-            # Look for the estimations in the peak of the histogram
-            w_theta_needed = 0
-            idx = []
-            while sum(idx) == 0:
-                w_theta_needed += self.hist_bin_width  # widen peak as needed
-                idx = np.abs(angles - angles_grid[peak_idx]) < w_theta_needed
-            if w_theta_needed > self.hist_bin_width:
-                logger.info(
-                    f"Adaptive width {w_theta_needed} required for ({i},{j}), found {sum(idx)} indices."
-                )
-        else:
-            # Fixed width
-            idx = np.abs(angles - angles_grid[peak_idx]) < self.full_width
-
-        good_k = inds[idx]
-        alpha = np.arccos(phis[idx])
-
-        return alpha, good_k.astype("int")
-
-    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
-
-        Given C1, C2, and C3 are unit circles of image i, j, and k, compute
-        resulting cos values of rotation angles between i an j images when both
-        of them are intersecting with k.
-
-        To ensure that the smallest singular value is big enough, controlled by
-        the determinant of the matrix,
-           C=[  1  c1  c2 ;
-               c1   1  c3 ;
-               c2  c3   1 ],
-        we therefore use the condition below
-               1+2*c1*c2*c3-(c1^2+c2^2+c3^2) > 1.0e-5,
-        so the matrix is far from singular.
-
-        :param cl_diff1: Difference of common line indices on C1 created by
-            its intersection with C3 and C2
-        :param cl_diff2: Difference of common line indices on C2 created by
-            its intersection with C1 and C3
-        :param cl_diff3: Difference of common line indices on C3 created by
-            its intersection with C2 and C1
-        :param n_theta: The number of points in the theta direction (common lines)
-        :param sync: Perform 180 degree ambiguity synchronization.
-        :return: cos values of rotation angles between i and j images
-            and indices for good k
-        """
-
-        # Calculate the theta values from the differences of common line indices
-        # C1, C2, and C3 are unit circles of image i, j, and k
-        # theta1 is the angle on C1 created by its intersection with C3 and C2.
-        # theta2 is the angle on C2 created by its intersection with C1 and C3.
-        # theta3 is the angle on C3 created by its intersection with C2 and C1.
-        theta1 = cl_diff1 * 2 * np.pi / n_theta
-        theta2 = cl_diff2 * 2 * np.pi / n_theta
-        theta3 = cl_diff3 * 2 * np.pi / n_theta
-
-        c1 = np.cos(theta1)
-        c2 = np.cos(theta2)
-        c3 = np.cos(theta3)
-
-        # Each common-line corresponds to a point on the unit sphere. Denote the
-        # coordinates of these points by (Pix, Piy Piz), and put them in the matrix
-        #   M=[ P1x  P2x  P3x ;
-        #       P1y  P2y  P3y ;
-        #       P1z  P2z  P3z ].
-        #
-        # Then the matrix
-        #   C=[  1  c1  c2 ;
-        #       c1   1  c3 ;
-        #       c2  c3   1 ],
-        # where c1, c2, c3 are given above, is given by C = M.T @ M.
-        # For the points P1, P2, and P3 to form a triangle on the unit sphere, a
-        # necessary and sufficient condition is for C to be positive definite. This
-        # is equivalent to
-        #       1+2*c1*c2*c3-(c1^2+c2^2+c3^2) > 0.
-        # However, this may result in a triangle that is too flat, that is, the
-        # angle between the projections is very close to zero. We therefore use the
-        # condition below
-        #       1+2*c1*c2*c3-(c1^2+c2^2+c3^2) > 1.0e-5.
-        # This ensures that the smallest singular value (which is actually
-        # controlled by the determinant of C) is big enough, so the matrix is far
-        # from singular. This condition is equivalent to computing the singular
-        # values of C, followed by checking that the smallest one is big enough.
-
-        cond = 1 + 2 * c1 * c2 * c3 - (np.square(c1) + np.square(c2) + np.square(c3))
-        good_idx = np.nonzero(cond > 1e-5)[0]
-
-        # Calculated cos values of angle between i and j images
-        if sync:
-            # MATLAB
-            cos_phi2 = (c3[good_idx] - c1[good_idx] * c2[good_idx]) / (
-                np.sqrt(1 - c1[good_idx] ** 2) * np.sqrt(1 - c2[good_idx] ** 2)
-            )
+    :return: The rotation matrix that takes image i to image j for good index of k.
+    """
 
-            #  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.
-            TOL_idx = 1e-12
+    if i == j:
+        return []
 
-            # Select only good_idx
-            theta1 = theta1[good_idx]
-            theta2 = theta2[good_idx]
-            theta3 = theta3[good_idx]
+    # Prepare the theta values from the differences of common line indices
+    # C1, C2, and C3 are unit circles of image i, j, and k
+    # cl_diff1 is for the angle on C1 created by its intersection with C3 and C2.
+    # 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_diff3 = clmatrix[good_k, j] - clmatrix[good_k, i]  # for theta3
 
-            # Check sync conditions
-            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)
-
-            # Apply sync
-            cos_phi2[align180] = -cos_phi2[align180]
-        else:
-            # Python
-            cos_phi2 = (c3[good_idx] - c1[good_idx] * c2[good_idx]) / (
-                np.sin(theta1[good_idx]) * np.sin(theta2[good_idx])
+    # Calculate the cos values of rotation angles between i an j images for good k images
+    c_alpha, good_idx = _get_cos_phis(cl_diff1, cl_diff2, cl_diff3, n_theta, sync=False)
+
+    if len(c_alpha) == 0:
+        return None
+    alpha = np.arccos(c_alpha)
+
+    # Convert the Euler angles with ZYZ conversion to rotation matrices
+    angles = np.zeros((alpha.shape[0], 3))
+    angles[:, 0] = clmatrix[i, j] * 2 * np.pi / n_theta + np.pi / 2
+    angles[:, 1] = alpha
+    angles[:, 2] = -np.pi / 2 - clmatrix[j, i] * 2 * np.pi / n_theta
+    r = Rotation.from_euler(angles).matrices
+
+    return r[good_idx, :, :]
+
+
+def _vote_ij(
+    clmatrix, n_theta, i, j, k_list, hist_bin_width, full_width, sigma=3.0, sync=False
+):
+    """
+    Apply the voting algorithm for images i and j.
+
+    clmatrix is the common lines matrix, constructed using angular resolution,
+    n_theta. k_list are the images to be used for voting of the pair of images
+    (i ,j).
+
+    :param clmatrix: The common lines matrix
+    :param n_theta: The number of points in the theta direction (common lines)
+    :param i: The i image
+    :param j: The j image
+    :param k_list: The list of images for the third image for voting algorithm
+    :param hist_bin_width: Bin width in smoothing histogram (degrees).
+    :param full_width: Selection width around smoothed histogram peak (degrees).
+        `adaptive` will attempt to automatically find the smallest number of
+        `hist_bin_width`s required to find at least one valid image index.
+    :param sigma: Voting contribution smoothing factor. Default is 3.0.
+    :param sync: Perform 180 degree ambiguity synchronization.
+
+    :return: (alpha, good_k), angles and list of all third images
+        in the peak of the histogram corresponding to the pair of
+        images (i,j)
+    """
+
+    if i == j or clmatrix[i, j] == -1:
+        return None, []
+
+    # Some of the entries in clmatrix may be zero if we cleared
+    # them due to small correlation, or if for each image
+    # we compute intersections with only some of the other images.
+    #
+    # Note that as long as the diagonal of the common lines matrix is
+    # -1, the conditions (i != j) && (j != k) are not needed, since
+    # if i == j then clmatrix[i, k] == -1 and similarly for i == k or
+    # j == k. Thus, the previous voting code (from the JSB paper) is
+    # correct even though it seems that we should test also that
+    # (i != j) && (i != k) && (j != k), and only (i != j) && (i != k)
+    #  as tested there.
+    cl_idx12 = clmatrix[i, j]
+    cl_idx21 = clmatrix[j, i]
+    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]
+    cl_idx32 = clmatrix[k_list, j]
+
+    # Prepare the theta values from the differences of common line indices
+    # C1, C2, and C3 are unit circles of image i, j, and k
+    # cl_diff1 is for the angle on C1 created by its intersection with C3 and C2.
+    # 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 = cl_idx13 - cl_idx12
+    cl_diff2 = cl_idx23 - cl_idx21
+    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 = _get_cos_phis(cl_diff1, cl_diff2, cl_diff3, n_theta, sync=sync)
+
+    if np.any(np.abs(cos_phi2) - 1 > 1e-12):
+        logger.warning(
+            f"Globally Consistent Angular Reconstruction (GCAR) exists"
+            f" numerical problem: abs(cos_phi2) > 1, with the"
+            f" difference of {np.abs(cos_phi2)-1}."
+        )
+    cos_phi2 = np.clip(cos_phi2, -1, 1)
+
+    # Store angles between i and j induced by each third image k.
+    phis = cos_phi2
+    # Sore good indices of l in k_list of the image that creates that angle.
+    inds = k_list[good_idx]
+
+    if phis.shape[0] == 0:
+        return None, []
+
+    # Parameters used to compute the smoothed angle histogram.
+    ntics = int(180 / hist_bin_width)
+    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. `sigma` ensures that the width of the density
+    # estimation kernel is roughly 10 degrees. For 15 degrees, the
+    # value of the kernel is negligible.
+
+    # Compute the histogram of the angles between images i and j
+    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
+    # angle between i and j that is at most 10 off the true angle.
+    # Even for debugging, don't put a value that is smaller than two
+    # tics, since the peak might move a little bit due to wrong k images
+    # that accidentally fall near the peak.
+    peak_idx = angles_hist.argmax()
+
+    if full_width == -1:
+        # Adaptive width  (MATLAB)
+        # Look for the estimations in the peak of the histogram
+        w_theta_needed = 0
+        idx = []
+        while sum(idx) == 0:
+            w_theta_needed += hist_bin_width  # widen peak as needed
+            idx = np.abs(angles - angles_grid[peak_idx]) < w_theta_needed
+        if w_theta_needed > hist_bin_width:
+            logger.info(
+                f"Adaptive width {w_theta_needed} required for ({i},{j}), found {sum(idx)} indices."
             )
+    else:
+        # Fixed width
+        idx = np.abs(angles - angles_grid[peak_idx]) < full_width
+
+    good_k = inds[idx]
+    alpha = np.arccos(phis[idx])
+
+    return alpha, good_k.astype("int")
+
+
+def _get_cos_phis(cl_diff1, cl_diff2, cl_diff3, n_theta, sync=False):
+    """
+    Calculate cos values of rotation angles between i and j images
+
+    Given C1, C2, and C3 are unit circles of image i, j, and k, compute
+    resulting cos values of rotation angles between i an j images when both
+    of them are intersecting with k.
+
+    To ensure that the smallest singular value is big enough, controlled by
+    the determinant of the matrix,
+       C=[  1  c1  c2 ;
+           c1   1  c3 ;
+           c2  c3   1 ],
+    we therefore use the condition below
+           1+2*c1*c2*c3-(c1^2+c2^2+c3^2) > 1.0e-5,
+    so the matrix is far from singular.
+
+    :param cl_diff1: Difference of common line indices on C1 created by
+        its intersection with C3 and C2
+    :param cl_diff2: Difference of common line indices on C2 created by
+        its intersection with C1 and C3
+    :param cl_diff3: Difference of common line indices on C3 created by
+        its intersection with C2 and C1
+    :param n_theta: The number of points in the theta direction (common lines)
+    :param sync: Perform 180 degree ambiguity synchronization.
+
+    :return: cos values of rotation angles between i and j images
+        and indices for good k
+    """
+
+    # Calculate the theta values from the differences of common line indices
+    # C1, C2, and C3 are unit circles of image i, j, and k
+    # theta1 is the angle on C1 created by its intersection with C3 and C2.
+    # theta2 is the angle on C2 created by its intersection with C1 and C3.
+    # theta3 is the angle on C3 created by its intersection with C2 and C1.
+    theta1 = cl_diff1 * 2 * np.pi / n_theta
+    theta2 = cl_diff2 * 2 * np.pi / n_theta
+    theta3 = cl_diff3 * 2 * np.pi / n_theta
+
+    c1 = np.cos(theta1)
+    c2 = np.cos(theta2)
+    c3 = np.cos(theta3)
+
+    # Each common-line corresponds to a point on the unit sphere. Denote the
+    # coordinates of these points by (Pix, Piy Piz), and put them in the matrix
+    #   M=[ P1x  P2x  P3x ;
+    #       P1y  P2y  P3y ;
+    #       P1z  P2z  P3z ].
+    #
+    # Then the matrix
+    #   C=[  1  c1  c2 ;
+    #       c1   1  c3 ;
+    #       c2  c3   1 ],
+    # where c1, c2, c3 are given above, is given by C = M.T @ M.
+    # For the points P1, P2, and P3 to form a triangle on the unit sphere, a
+    # necessary and sufficient condition is for C to be positive definite. This
+    # is equivalent to
+    #       1+2*c1*c2*c3-(c1^2+c2^2+c3^2) > 0.
+    # However, this may result in a triangle that is too flat, that is, the
+    # angle between the projections is very close to zero. We therefore use the
+    # condition below
+    #       1+2*c1*c2*c3-(c1^2+c2^2+c3^2) > 1.0e-5.
+    # This ensures that the smallest singular value (which is actually
+    # controlled by the determinant of C) is big enough, so the matrix is far
+    # from singular. This condition is equivalent to computing the singular
+    # values of C, followed by checking that the smallest one is big enough.
+
+    cond = 1 + 2 * c1 * c2 * c3 - (np.square(c1) + np.square(c2) + np.square(c3))
+    good_idx = np.nonzero(cond > 1e-5)[0]
+
+    # Calculated cos values of angle between i and j images
+    if sync:
+        # MATLAB
+        cos_phi2 = (c3[good_idx] - c1[good_idx] * c2[good_idx]) / (
+            np.sqrt(1 - c1[good_idx] ** 2) * np.sqrt(1 - c2[good_idx] ** 2)
+        )
+
+        #  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.
+        TOL_idx = 1e-12
+
+        # Select only good_idx
+        theta1 = theta1[good_idx]
+        theta2 = theta2[good_idx]
+        theta3 = theta3[good_idx]
+
+        # Check sync conditions
+        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)
+
+        # Apply sync
+        cos_phi2[align180] = -cos_phi2[align180]
+    else:
+        # Python
+        cos_phi2 = (c3[good_idx] - c1[good_idx] * c2[good_idx]) / (
+            np.sin(theta1[good_idx]) * np.sin(theta2[good_idx])
+        )
 
-        return cos_phi2, good_idx
+    return cos_phi2, good_idx
diff --git a/src/aspire/utils/__init__.py b/src/aspire/utils/__init__.py
index ae781d823f..ae896ebeb8 100644
--- a/src/aspire/utils/__init__.py
+++ b/src/aspire/utils/__init__.py
@@ -1,5 +1,4 @@
 from .types import complex_type, real_type, utest_tolerance  # isort:skip
-
 from .coor_trans import (  # isort:skip
     mean_aligned_angular_distance,
     cart2pol,
diff --git a/tests/test_commonline_utils.py b/tests/test_commonline_utils.py
new file mode 100644
index 0000000000..64d015d1dd
--- /dev/null
+++ b/tests/test_commonline_utils.py
@@ -0,0 +1,118 @@
+import numpy as np
+import pytest
+
+from aspire.abinitio import JSync
+from aspire.abinitio.commonline_utils import (
+    _complete_third_row_to_rot,
+    _estimate_third_rows,
+    build_outer_products,
+)
+from aspire.utils import J_conjugate, Rotation, randn, utest_tolerance
+
+DTYPES = [np.float32, np.float64]
+
+
+@pytest.fixture(params=DTYPES, ids=lambda x: f"dtype={x}", scope="module")
+def dtype(request):
+    return request.param
+
+
+def test_estimate_third_rows(dtype):
+    """
+    Test we accurately estimate a set of 3rd rows of rotation matrices
+    given the 3rd row outer products vijs =  vi @ vj.T and viis = vi @ vi.T.
+    """
+    n_img = 20
+
+    # `build_outer_products` generates a set of ground truth 3rd rows
+    # of rotation matrices, then forms the outer products vijs =  vi @ vj.T
+    # and viis = vi @ vi.T.
+    vijs, viis, gt_vis = build_outer_products(n_img, dtype)
+
+    # Estimate third rows from outer products.
+    # Due to factorization of V, these might be negated third rows.
+    vis = _estimate_third_rows(vijs, viis)
+
+    # Check if all-close up to difference of sign
+    ground_truth = np.sign(gt_vis[0, 0]) * gt_vis
+    estimate = np.sign(vis[0, 0]) * vis
+    np.testing.assert_allclose(ground_truth, estimate, rtol=1e-05, atol=1e-08)
+
+    # Check dtype passthrough
+    assert vis.dtype == dtype
+
+
+def test_complete_third_row(dtype):
+    """
+    Test that `complete_third_row_to_rot` produces a proper rotations
+    given a set of 3rd rows.
+    """
+    # Build random third rows.
+    r3 = randn(10, 3, seed=123).astype(dtype)
+    r3 /= np.linalg.norm(r3, axis=1)[..., np.newaxis]
+
+    # Set first row to be identical with z-axis.
+    r3[0] = np.array([0, 0, 1], dtype=dtype)
+
+    # Generate rotations.
+    R = _complete_third_row_to_rot(r3)
+
+    # Check dtype passthrough
+    assert R.dtype == dtype
+
+    # Assert that first rotation is the identity matrix.
+    np.testing.assert_allclose(R[0], np.eye(3, dtype=dtype))
+
+    # Assert that each rotation is orthogonal with determinant 1.
+    assert np.allclose(
+        R @ R.transpose((0, 2, 1)), np.eye(3, dtype=dtype), atol=utest_tolerance(dtype)
+    )
+    assert np.allclose(np.linalg.det(R), 1)
+
+
+def test_J_sync(dtype):
+    """
+    Test that the J_sync `power_method` returns a set of signs indicating
+    the set of relative rotations that need to be J-conjugated to attain
+    global handedness consistency, and that `global_J_sync` returns the
+    ground truth rotations up to a spurious J-conjugation.
+    """
+    n = 25
+    rots = Rotation.generate_random_rotations(n, dtype=dtype).matrices
+
+    # Generate ground truth and randomly J-conjugate relative rotations,
+    # keeping track of the signs associated with J-conjugated rotations.
+    n_choose_2 = (n * (n - 1)) // 2
+    signs = np.random.randint(0, 2, n_choose_2) * 2 - 1
+    Rijs_gt = np.zeros((n_choose_2, 3, 3), dtype=dtype)
+    Rijs_conjugated = np.zeros((n_choose_2, 3, 3), dtype=dtype)
+    ij = 0
+    for i in range(n - 1):
+        Ri = rots[i]
+        for j in range(i + 1, n):
+            Rj = rots[j]
+            Rijs_gt[ij] = Rij = Ri.T @ Rj
+            if signs[ij] == -1:
+                Rij = J_conjugate(Rij)
+            Rijs_conjugated[ij] = Rij
+            ij += 1
+
+    # Initialize JSync instance with default params.
+    J_sync = JSync(n)
+
+    # Perform power method and check that signs are correct up to
+    # multilication by -1. Also check dtype pass-through.
+    signs_est = J_sync.power_method(Rijs_conjugated)
+    np.testing.assert_allclose(signs[0] * signs, signs_est[0] * signs_est)
+    assert signs_est.dtype == dtype
+
+    # Perform global J sync and check that rotations are correct up to
+    # a spurious J conjugation. Also check dtype pass-through.
+    Rijs_sync = J_sync.global_J_sync(Rijs_conjugated)
+
+    # If the first is off by a J, J-conjugate the whole set.
+    if np.allclose(Rijs_gt[0], J_conjugate(Rijs_sync[0])):
+        Rijs_sync = J_conjugate(Rijs_sync)
+
+    np.testing.assert_allclose(Rijs_sync, Rijs_gt)
+    assert Rijs_sync.dtype == dtype
diff --git a/tests/test_orient_symmetric.py b/tests/test_orient_symmetric.py
index d7c4c4716f..ed9c5a6904 100644
--- a/tests/test_orient_symmetric.py
+++ b/tests/test_orient_symmetric.py
@@ -1,17 +1,15 @@
 import numpy as np
 import pytest
-from numpy import pi, random
-from numpy.linalg import det, norm
 
 from aspire.abinitio import (
     CLSymmetryC2,
     CLSymmetryC3C4,
     CLSymmetryCn,
-    cl_angles_to_ind,
-    complete_third_row_to_rot,
-    estimate_third_rows,
+    build_outer_products,
+    g_sync,
 )
 from aspire.abinitio.commonline_cn import MeanOuterProductEstimator
+from aspire.abinitio.commonline_utils import _cl_angles_to_ind
 from aspire.source import Simulation
 from aspire.utils import (
     J_conjugate,
@@ -19,8 +17,6 @@
     all_pairs,
     cyclic_rotations,
     mean_aligned_angular_distance,
-    randn,
-    utest_tolerance,
 )
 from aspire.volume import CnSymmetricVolume
 
@@ -127,7 +123,7 @@ def test_estimate_rotations(n_img, L, order, dtype):
     rots_gt = src.rotations
 
     # g-synchronize ground truth rotations.
-    rots_gt_sync = cl_symm.g_sync(rots_est, order, rots_gt)
+    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.
@@ -141,9 +137,7 @@ def test_relative_rotations(n_img, L, order, dtype):
     src, cl_symm = source_orientation_objs(n_img, L, order, dtype)
 
     # Estimate relative viewing directions.
-    cl_symm.build_clmatrix()
-    cl = cl_symm.clmatrix
-    Rijs = cl_symm._estimate_all_Rijs_c3_c4(cl)
+    Rijs = cl_symm._estimate_all_Rijs_c3_c4()
 
     # Each Rij belongs to the set {Ri.Tg_n^sRj, JRi.Tg_n^sRjJ},
     # s = 1, 2, ..., order. We find the mean squared error over
@@ -326,8 +320,8 @@ def test_self_commonlines(n_img, L, order, dtype):
     # Get angle difference between scl_gt and scl.
     scl_diff1 = scl_gt - scl
     scl_diff2 = scl_gt - np.flip(scl, 1)  # Order of indices might be switched.
-    scl_diff1_angle = scl_diff1 * 2 * pi / n_theta
-    scl_diff2_angle = scl_diff2 * 2 * pi / n_theta
+    scl_diff1_angle = scl_diff1 * 2 * np.pi / n_theta
+    scl_diff2_angle = scl_diff2 * 2 * np.pi / n_theta
 
     # cosine is invariant to 2pi, and abs is invariant to +-pi due to J-conjugation.
     # We take the mean deviation wrt to the two lines in each image.
@@ -339,7 +333,7 @@ def test_self_commonlines(n_img, L, order, dtype):
     min_mean_angle_diff = scl_idx.choose(scl_diff_angle_mean)
 
     # Assert scl detection rate is 100% for 5 degree angle tolerance
-    angle_tol_err = 5 * pi / 180
+    angle_tol_err = 5 * np.pi / 180
     detection_rate = np.count_nonzero(min_mean_angle_diff < angle_tol_err) / len(scl)
     assert np.allclose(detection_rate, 1.0)
 
@@ -484,45 +478,6 @@ def test_global_J_sync(n_img, dtype):
         assert np.allclose(viis, viis_sync)
 
 
-@pytest.mark.parametrize("dtype", [np.float32, np.float64])
-def test_estimate_third_rows(dtype):
-    n_img = 20
-
-    # Build outer products vijs, viis, and get ground truth third rows.
-    vijs, viis, gt_vis = build_outer_products(n_img, dtype)
-
-    # Estimate third rows from outer products.
-    # Due to factorization of V, these might be negated third rows.
-    vis = estimate_third_rows(vijs, viis)
-
-    # Check if all-close up to difference of sign
-    ground_truth = np.sign(gt_vis[0, 0]) * gt_vis
-    estimate = np.sign(vis[0, 0]) * vis
-    assert np.allclose(ground_truth, estimate)
-
-
-@pytest.mark.parametrize("dtype", [np.float32, np.float64])
-def test_complete_third_row(dtype):
-    # Build random third rows.
-    r3 = randn(10, 3, seed=123).astype(dtype)
-    r3 /= norm(r3, axis=1)[..., np.newaxis]
-
-    # Set first row to be identical with z-axis.
-    r3[0] = np.array([0, 0, 1], dtype=dtype)
-
-    # Generate rotations.
-    R = complete_third_row_to_rot(r3)
-
-    # Assert that first rotation is the identity matrix.
-    assert np.allclose(R[0], np.eye(3, dtype=dtype))
-
-    # Assert that each rotation is orthogonal with determinant 1.
-    assert np.allclose(
-        R @ R.transpose((0, 2, 1)), np.eye(3, dtype=dtype), atol=utest_tolerance(dtype)
-    )
-    assert np.allclose(det(R), 1)
-
-
 @pytest.mark.parametrize("dtype", [np.float32, np.float64])
 def test_dtype_pass_through(dtype):
     L = 16
@@ -558,31 +513,6 @@ def build_self_commonlines_matrix(n_theta, rots, order):
     return scl_gt
 
 
-def build_outer_products(n_img, dtype):
-    # Build random third rows, ground truth vis (unit vectors)
-    gt_vis = np.zeros((n_img, 3), dtype=dtype)
-    for i in range(n_img):
-        random.seed(i)
-        v = random.randn(3)
-        gt_vis[i] = v / norm(v)
-
-    # Find outer products viis and vijs for i<j
-    nchoose2 = int(n_img * (n_img - 1) / 2)
-    vijs = np.zeros((nchoose2, 3, 3), dtype=dtype)
-    viis = np.zeros((n_img, 3, 3), dtype=dtype)
-
-    # All pairs (i,j) where i<j
-    pairs = all_pairs(n_img)
-
-    for k, (i, j) in enumerate(pairs):
-        vijs[k] = np.outer(gt_vis[i], gt_vis[j])
-
-    for i in range(n_img):
-        viis[i] = np.outer(gt_vis[i], gt_vis[i])
-
-    return vijs, viis, gt_vis
-
-
 def test_mean_outer_product_estimator():
     """
     Manully run MeanOuterProductEstimator for prebaked inputs.
@@ -642,6 +572,6 @@ def _gt_cl_c2(n_theta, rots_gt):
                 U = Ri.T @ g @ Rj
                 c1 = np.array([-U[1, 2], U[0, 2]])
                 c2 = np.array([U[2, 1], -U[2, 0]])
-                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)
+                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