registration

occulus.registration

Point cloud registration — alignment of multiple scans.

Primary interface::

from occulus.registration import icp, ransac_registration, RegistrationResult
from occulus.registration import align_scans, compute_fpfh

Algorithms

• :func:icp — point-to-point or point-to-plane ICP (auto-selects based on normals)
• :func:icp_point_to_point — explicit point-to-point ICP
• :func:icp_point_to_plane — explicit point-to-plane ICP (requires target normals)
• :func:compute_fpfh — 33-dim FPFH feature descriptors for global registration
• :func:ransac_registration — feature-matching + RANSAC global alignment
• :func:align_scans — sequential multi-scan alignment

AlignmentResult

Result of multi-scan alignment.

Attributes:

Name	Type	Description
transformations	list[NDArray[float64]]	List of 4×4 transformations, one per input cloud. The first cloud's transformation is always the identity (it is the reference).
pairwise_results	list[RegistrationResult]	Pairwise registration results between consecutive scans.
global_rmse	float	Mean inlier RMSE across all pairwise registrations.

Source code in src/occulus/registration/global_registration.py

class AlignmentResult:
    """Result of multi-scan alignment.

    Attributes
    ----------
    transformations : list[NDArray[np.float64]]
        List of 4×4 transformations, one per input cloud. The first cloud's
        transformation is always the identity (it is the reference).
    pairwise_results : list[RegistrationResult]
        Pairwise registration results between consecutive scans.
    global_rmse : float
        Mean inlier RMSE across all pairwise registrations.
    """

    def __init__(
        self,
        transformations: list[NDArray[np.float64]],
        pairwise_results: list[RegistrationResult],
    ) -> None:
        self.transformations = transformations
        self.pairwise_results = pairwise_results
        self.global_rmse = (
            float(np.mean([r.inlier_rmse for r in pairwise_results])) if pairwise_results else 0.0
        )

RegistrationResult dataclass

Result of an ICP or global registration.

Attributes:

Name	Type	Description
transformation	NDArray[float64]	4×4 rigid transformation matrix (source → target frame).
fitness	float	Fraction of source points with a correspondence within max_correspondence_distance. Range [0, 1]; higher is better.
inlier_rmse	float	Root-mean-square error of inlier correspondences (in coordinate units).
converged	bool	Whether the algorithm met its convergence criterion.
n_iterations	int	Number of iterations actually performed.

Source code in src/occulus/registration/icp.py

@dataclass
class RegistrationResult:
    """Result of an ICP or global registration.

    Attributes
    ----------
    transformation : NDArray[np.float64]
        4×4 rigid transformation matrix (source → target frame).
    fitness : float
        Fraction of source points with a correspondence within
        ``max_correspondence_distance``. Range [0, 1]; higher is better.
    inlier_rmse : float
        Root-mean-square error of inlier correspondences (in coordinate units).
    converged : bool
        Whether the algorithm met its convergence criterion.
    n_iterations : int
        Number of iterations actually performed.
    """

    transformation: NDArray[np.float64]
    fitness: float
    inlier_rmse: float
    converged: bool
    n_iterations: int = 0

align_scans(clouds, *, voxel_size=0.5, refine_with_icp=True, max_correspondence_distance=1.0)

Align multiple scans into a common coordinate system.

Performs sequential pairwise ICP registration — each scan is aligned to the previous one. The first cloud is the reference (identity transform).

For best results, pre-process clouds with :func:~occulus.filters.voxel_downsample and :func:~occulus.normals.estimate_normals before calling this function.

Parameters:

Name	Type	Description	Default
clouds	list[PointCloud]	At least 2 PointCloud objects to align. The first cloud is the fixed reference.	required
voxel_size	float	Voxel size used for temporary downsampling before feature-based coarse alignment, by default 0.5.	0.5
refine_with_icp	bool	Whether to run ICP after coarse alignment, by default True.	True
max_correspondence_distance	float	Maximum correspondence distance for ICP refinement, by default 1.0.	1.0

Returns:

Type	Description
AlignmentResult	Per-cloud transformations and pairwise quality metrics.

Raises:

Type	Description
OcculusRegistrationError	If fewer than 2 clouds are provided.

Source code in src/occulus/registration/global_registration.py

def align_scans(
    clouds: list[PointCloud],
    *,
    voxel_size: float = 0.5,
    refine_with_icp: bool = True,
    max_correspondence_distance: float = 1.0,
) -> AlignmentResult:
    """Align multiple scans into a common coordinate system.

    Performs sequential pairwise ICP registration — each scan is aligned
    to the previous one. The first cloud is the reference (identity transform).

    For best results, pre-process clouds with
    :func:`~occulus.filters.voxel_downsample` and
    :func:`~occulus.normals.estimate_normals` before calling this function.

    Parameters
    ----------
    clouds : list[PointCloud]
        At least 2 PointCloud objects to align. The first cloud is the
        fixed reference.
    voxel_size : float, optional
        Voxel size used for temporary downsampling before feature-based
        coarse alignment, by default 0.5.
    refine_with_icp : bool, optional
        Whether to run ICP after coarse alignment, by default ``True``.
    max_correspondence_distance : float, optional
        Maximum correspondence distance for ICP refinement, by default 1.0.

    Returns
    -------
    AlignmentResult
        Per-cloud transformations and pairwise quality metrics.

    Raises
    ------
    OcculusRegistrationError
        If fewer than 2 clouds are provided.
    """
    if len(clouds) < 2:
        raise OcculusRegistrationError("align_scans requires at least 2 clouds")

    from occulus.registration.icp import icp as _icp  # avoid circular at module level

    transforms: list[NDArray[np.float64]] = [np.eye(4, dtype=np.float64)]
    results: list[RegistrationResult] = []

    # Cumulative transformation — chain pairwise results
    cumulative = np.eye(4, dtype=np.float64)

    for i in range(1, len(clouds)):
        target = clouds[i - 1]
        source = clouds[i]

        logger.info("align_scans: aligning cloud %d → %d", i, i - 1)
        result = _icp(
            source,
            target,
            max_correspondence_distance=max_correspondence_distance,
            method="auto",
        )
        results.append(result)
        cumulative = result.transformation @ cumulative
        transforms.append(cumulative.copy())

    return AlignmentResult(transforms, results)

compute_fpfh(cloud, radius, *, max_nn=100)

Compute Fast Point Feature Histogram (FPFH) descriptors.

FPFH is a 33-dimensional descriptor that encodes the local geometry around each point using angular features derived from surface normals. Requires normals on the input cloud.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	Input cloud. Must have normals (cloud.has_normals must be True).	required
radius	float	Neighbourhood search radius for feature computation.	required
max_nn	int	Maximum number of neighbours per point, by default 100.	100

Returns:

Type	Description
NDArray[float64]	FPFH descriptor array of shape (N, 33).

Raises:

Type	Description
OcculusValidationError	If the cloud has no normals.

Source code in src/occulus/registration/global_registration.py

def compute_fpfh(
    cloud: PointCloud,
    radius: float,
    *,
    max_nn: int = 100,
) -> NDArray[np.float64]:
    """Compute Fast Point Feature Histogram (FPFH) descriptors.

    FPFH is a 33-dimensional descriptor that encodes the local geometry
    around each point using angular features derived from surface normals.
    Requires normals on the input cloud.

    Parameters
    ----------
    cloud : PointCloud
        Input cloud. Must have normals (``cloud.has_normals`` must be ``True``).
    radius : float
        Neighbourhood search radius for feature computation.
    max_nn : int, optional
        Maximum number of neighbours per point, by default 100.

    Returns
    -------
    NDArray[np.float64]
        FPFH descriptor array of shape (N, 33).

    Raises
    ------
    OcculusValidationError
        If the cloud has no normals.
    """
    if not cloud.has_normals:
        raise OcculusValidationError(
            "compute_fpfh requires normals. Run occulus.normals.estimate_normals() first."
        )
    assert cloud.normals is not None

    xyz = cloud.xyz
    normals = cloud.normals
    n = len(xyz)
    tree = KDTree(xyz)

    # Simplified Point Feature Histogram (SPFH) for each point
    spfh = np.zeros((n, 33), dtype=np.float64)
    neighbour_lists = tree.query_ball_point(xyz, r=radius, workers=-1)

    for i, neighbours in enumerate(neighbour_lists):
        if len(neighbours) < 2:
            continue

        idx = np.array(neighbours)
        if len(idx) > max_nn:
            dists = np.linalg.norm(xyz[idx] - xyz[i], axis=1)
            idx = idx[np.argpartition(dists, max_nn)[:max_nn]]

        # Exclude self
        idx = idx[idx != i]
        if len(idx) == 0:
            continue

        # Darboux frame for each neighbour pair (i, j)
        u = normals[i]  # (3,)
        diff = xyz[idx] - xyz[i]  # (k, 3)
        dist = np.linalg.norm(diff, axis=1, keepdims=True)
        dist = np.where(dist < 1e-12, 1.0, dist)
        v_vecs = diff / dist  # (k, 3)
        w_vecs = np.cross(u, v_vecs)  # (k, 3)

        # Three angles: alpha, phi, theta
        nj = normals[idx]  # (k, 3)
        alpha = np.einsum("ij,ij->i", w_vecs, nj)  # (k,)
        phi = np.einsum("ij,j->i", v_vecs, u)  # (k,) — u is (3,)
        theta = np.arctan2(
            np.einsum("ij,ij->i", np.cross(u, v_vecs), nj),
            np.einsum("ij,ij->i", v_vecs, nj),
        )  # (k,)

        # Bin into 11 bins each (→ 33-dim descriptor)
        bins = 11
        spfh[i, :bins], _ = np.histogram(alpha, bins=bins, range=(-1.0, 1.0))
        spfh[i, bins : 2 * bins], _ = np.histogram(phi, bins=bins, range=(-1.0, 1.0))
        spfh[i, 2 * bins :], _ = np.histogram(theta, bins=bins, range=(-np.pi, np.pi))

        # Normalise
        row_sum = spfh[i].sum()
        if row_sum > 0:
            spfh[i] /= row_sum

    # FPFH = SPFH + weighted sum of SPFH of neighbours
    fpfh = spfh.copy()
    for i, neighbours in enumerate(neighbour_lists):
        idx = np.array(neighbours)
        idx = idx[idx != i]
        if len(idx) == 0:
            continue
        dists = np.linalg.norm(xyz[idx] - xyz[i], axis=1)
        dists = np.where(dists < 1e-12, 1e-12, dists)
        weights = 1.0 / dists
        fpfh[i] += (weights[:, None] * spfh[idx]).sum(axis=0) / weights.sum()

        row_sum = fpfh[i].sum()
        if row_sum > 0:
            fpfh[i] /= row_sum

    logger.debug("compute_fpfh: computed 33-dim descriptors for %d points", n)
    return fpfh

icp_point_to_plane(source, target, *, max_correspondence_distance=1.0, max_iterations=50, tolerance=1e-06, init_transform=None)

Point-to-plane ICP registration.

Minimises the sum of squared distances projected onto target surface normals. Requires normals on the target cloud. Typically converges faster and with better accuracy than point-to-point near the solution.

Parameters:

Name	Type	Description	Default
source	PointCloud	Source cloud to be transformed.	required
target	PointCloud	Target cloud with normals (fixed reference).	required
max_correspondence_distance	float	Maximum distance for a pair to be a valid correspondence.	1.0
max_iterations	int	Maximum ICP iterations.	50
tolerance	float	Convergence tolerance on transformation change.	1e-06
init_transform	NDArray[float64] | None	Initial 4×4 transformation guess.	None

Returns:

Type	Description
RegistrationResult	Registration result.

Raises:

Type	Description
OcculusValidationError	If target cloud does not have normals.
OcculusRegistrationError	If too few correspondences are found.

Source code in src/occulus/registration/icp.py

def icp_point_to_plane(
    source: PointCloud,
    target: PointCloud,
    *,
    max_correspondence_distance: float = 1.0,
    max_iterations: int = 50,
    tolerance: float = 1e-6,
    init_transform: NDArray[np.float64] | None = None,
) -> RegistrationResult:
    """Point-to-plane ICP registration.

    Minimises the sum of squared distances projected onto target surface
    normals. Requires normals on the target cloud. Typically converges
    faster and with better accuracy than point-to-point near the solution.

    Parameters
    ----------
    source : PointCloud
        Source cloud to be transformed.
    target : PointCloud
        Target cloud with normals (fixed reference).
    max_correspondence_distance : float, optional
        Maximum distance for a pair to be a valid correspondence.
    max_iterations : int, optional
        Maximum ICP iterations.
    tolerance : float, optional
        Convergence tolerance on transformation change.
    init_transform : NDArray[np.float64] | None, optional
        Initial 4×4 transformation guess.

    Returns
    -------
    RegistrationResult
        Registration result.

    Raises
    ------
    OcculusValidationError
        If target cloud does not have normals.
    OcculusRegistrationError
        If too few correspondences are found.
    """
    if not target.has_normals:
        raise OcculusValidationError(
            "Point-to-plane ICP requires normals on the target cloud. "
            "Use occulus.normals.estimate_normals(target) first."
        )
    assert target.normals is not None  # guaranteed above

    transform = _init_transform(init_transform)
    src_xyz = _apply_transform(source.xyz, transform)
    tree = KDTree(target.xyz)

    prev_transform = transform.copy()
    converged = False

    for iteration in range(max_iterations):
        distances, indices = tree.query(src_xyz, k=1, workers=-1)
        indices = indices.ravel()
        distances = distances.ravel()

        inlier_mask = distances <= max_correspondence_distance
        if inlier_mask.sum() < 6:
            raise OcculusRegistrationError(
                f"Too few correspondences ({inlier_mask.sum()}) at iteration {iteration}."
            )

        src_pts = src_xyz[inlier_mask]
        tgt_pts = target.xyz[indices[inlier_mask]]
        tgt_nrm = target.normals[indices[inlier_mask]]

        # Point-to-plane linear system (Segal et al.)
        # Build 6-unknown system [r1 r2 r3 t1 t2 t3] via least squares
        A, b = _build_p2plane_system(src_pts, tgt_pts, tgt_nrm)
        try:
            x, _, _, _ = np.linalg.lstsq(A, b, rcond=None)
        except np.linalg.LinAlgError as exc:
            raise OcculusRegistrationError(
                f"Point-to-plane linear system is singular at iteration {iteration}: {exc}"
            ) from exc

        # Extract rotation (small-angle approximation → exact rotation via Rodrigues)
        rx, ry, rz, tx, ty, tz = x
        delta = np.eye(4)
        delta[:3, :3] = _rodrigues_to_matrix(np.array([rx, ry, rz]))
        delta[:3, 3] = [tx, ty, tz]
        transform = delta @ transform
        src_xyz = _apply_transform(source.xyz, transform)

        change = np.linalg.norm(transform - prev_transform, ord="fro")
        if change < tolerance:
            converged = True
            logger.debug("ICP p2plane converged at iteration %d", iteration)
            break
        prev_transform = transform.copy()

    fitness, rmse = _compute_metrics(src_xyz, target.xyz, max_correspondence_distance)
    logger.info(
        "ICP p2plane: fitness=%.3f rmse=%.4f converged=%s iter=%d",
        fitness,
        rmse,
        converged,
        iteration + 1,
    )
    return RegistrationResult(
        transformation=transform,
        fitness=fitness,
        inlier_rmse=rmse,
        converged=converged,
        n_iterations=iteration + 1,
    )

icp_point_to_point(source, target, *, max_correspondence_distance=1.0, max_iterations=50, tolerance=1e-06, init_transform=None)

Point-to-point ICP registration using SVD.

Minimises the sum of squared Euclidean distances between each source point and its closest target point. The optimal rigid transformation is solved analytically via Singular Value Decomposition.

Parameters:

Name	Type	Description	Default
source	PointCloud	Source cloud to be transformed.	required
target	PointCloud	Target cloud (fixed reference).	required
max_correspondence_distance	float	Maximum distance for a pair to be a valid correspondence.	1.0
max_iterations	int	Maximum ICP iterations.	50
tolerance	float	Convergence tolerance on transformation change.	1e-06
init_transform	NDArray[float64] | None	Initial 4×4 transformation guess.	None

Returns:

Type	Description
RegistrationResult	Registration result.

Raises:

Type	Description
OcculusRegistrationError	If too few correspondences are found to compute a transformation.

Source code in src/occulus/registration/icp.py

def icp_point_to_point(
    source: PointCloud,
    target: PointCloud,
    *,
    max_correspondence_distance: float = 1.0,
    max_iterations: int = 50,
    tolerance: float = 1e-6,
    init_transform: NDArray[np.float64] | None = None,
) -> RegistrationResult:
    """Point-to-point ICP registration using SVD.

    Minimises the sum of squared Euclidean distances between each source
    point and its closest target point. The optimal rigid transformation is
    solved analytically via Singular Value Decomposition.

    Parameters
    ----------
    source : PointCloud
        Source cloud to be transformed.
    target : PointCloud
        Target cloud (fixed reference).
    max_correspondence_distance : float, optional
        Maximum distance for a pair to be a valid correspondence.
    max_iterations : int, optional
        Maximum ICP iterations.
    tolerance : float, optional
        Convergence tolerance on transformation change.
    init_transform : NDArray[np.float64] | None, optional
        Initial 4×4 transformation guess.

    Returns
    -------
    RegistrationResult
        Registration result.

    Raises
    ------
    OcculusRegistrationError
        If too few correspondences are found to compute a transformation.
    """
    transform = _init_transform(init_transform)
    src_xyz = _apply_transform(source.xyz, transform)
    tree = KDTree(target.xyz)

    prev_transform = transform.copy()
    converged = False

    for iteration in range(max_iterations):
        distances, indices = tree.query(src_xyz, k=1, workers=-1)
        indices = indices.ravel()
        distances = distances.ravel()

        inlier_mask = distances <= max_correspondence_distance
        if inlier_mask.sum() < 3:
            raise OcculusRegistrationError(
                f"Too few correspondences ({inlier_mask.sum()}) at iteration {iteration}. "
                "Increase max_correspondence_distance or improve initial alignment."
            )

        src_inliers = src_xyz[inlier_mask]
        tgt_inliers = target.xyz[indices[inlier_mask]]

        delta_R, delta_t = _svd_rigid(src_inliers, tgt_inliers)

        # Apply incremental transformation
        delta = np.eye(4)
        delta[:3, :3] = delta_R
        delta[:3, 3] = delta_t
        transform = delta @ transform
        src_xyz = _apply_transform(source.xyz, transform)

        # Convergence check
        change = np.linalg.norm(transform - prev_transform, ord="fro")
        if change < tolerance:
            converged = True
            logger.debug("ICP p2p converged at iteration %d (change=%.2e)", iteration, change)
            break
        prev_transform = transform.copy()

    fitness, rmse = _compute_metrics(src_xyz, target.xyz, max_correspondence_distance)
    logger.info(
        "ICP p2p: fitness=%.3f rmse=%.4f converged=%s iter=%d",
        fitness,
        rmse,
        converged,
        iteration + 1,
    )
    return RegistrationResult(
        transformation=transform,
        fitness=fitness,
        inlier_rmse=rmse,
        converged=converged,
        n_iterations=iteration + 1,
    )

ransac_registration(source, target, source_features, target_features, *, max_correspondence_distance=1.5, ransac_n=3, max_iterations=100000, confidence=0.999)

Global registration via feature-matching and RANSAC.

Matches FPFH descriptors between source and target to find feature correspondences, then uses RANSAC to robustly estimate the rigid transformation.

Parameters:

Name	Type	Description	Default
source	PointCloud	Source cloud.	required
target	PointCloud	Target cloud.	required
source_features	NDArray[float64]	FPFH descriptors for source cloud (N, 33).	required
target_features	NDArray[float64]	FPFH descriptors for target cloud (M, 33).	required
max_correspondence_distance	float	Maximum Euclidean distance for an inlier correspondence, by default 1.5.	1.5
ransac_n	int	Number of random correspondences per RANSAC hypothesis, by default 3.	3
max_iterations	int	Maximum RANSAC iterations, by default 100_000.	100000
confidence	float	Desired success probability for RANSAC, used for adaptive early exit, by default 0.999.	0.999

Returns:

Type	Description
RegistrationResult	Best transformation found. converged=True if the inlier ratio exceeds the RANSAC confidence criterion.

Raises:

Type	Description
OcculusValidationError	If feature array shapes do not match their clouds.
OcculusRegistrationError	If no valid transformation is found.

Source code in src/occulus/registration/global_registration.py

def ransac_registration(
    source: PointCloud,
    target: PointCloud,
    source_features: NDArray[np.float64],
    target_features: NDArray[np.float64],
    *,
    max_correspondence_distance: float = 1.5,
    ransac_n: int = 3,
    max_iterations: int = 100_000,
    confidence: float = 0.999,
) -> RegistrationResult:
    """Global registration via feature-matching and RANSAC.

    Matches FPFH descriptors between source and target to find feature
    correspondences, then uses RANSAC to robustly estimate the rigid
    transformation.

    Parameters
    ----------
    source : PointCloud
        Source cloud.
    target : PointCloud
        Target cloud.
    source_features : NDArray[np.float64]
        FPFH descriptors for source cloud (N, 33).
    target_features : NDArray[np.float64]
        FPFH descriptors for target cloud (M, 33).
    max_correspondence_distance : float, optional
        Maximum Euclidean distance for an inlier correspondence,
        by default 1.5.
    ransac_n : int, optional
        Number of random correspondences per RANSAC hypothesis, by default 3.
    max_iterations : int, optional
        Maximum RANSAC iterations, by default 100_000.
    confidence : float, optional
        Desired success probability for RANSAC, used for adaptive early exit,
        by default 0.999.

    Returns
    -------
    RegistrationResult
        Best transformation found. ``converged=True`` if the inlier ratio
        exceeds the RANSAC confidence criterion.

    Raises
    ------
    OcculusValidationError
        If feature array shapes do not match their clouds.
    OcculusRegistrationError
        If no valid transformation is found.
    """
    if source_features.shape[0] != source.n_points:
        raise OcculusValidationError(
            f"source_features rows ({source_features.shape[0]}) != "
            f"source n_points ({source.n_points})"
        )
    if target_features.shape[0] != target.n_points:
        raise OcculusValidationError(
            f"target_features rows ({target_features.shape[0]}) != "
            f"target n_points ({target.n_points})"
        )

    # Build a KD-tree on target feature space to find descriptor matches
    feat_tree = KDTree(target_features)
    _distances, target_idx = feat_tree.query(source_features, k=1, workers=-1)
    target_idx = target_idx.ravel()

    # All candidate correspondences (source_i → target_idx[i])
    corr_src = np.arange(source.n_points)
    corr_tgt = target_idx

    rng = np.random.default_rng()
    best_transform = np.eye(4, dtype=np.float64)
    best_inliers = 0
    n_corr = len(corr_src)

    if n_corr < ransac_n:
        raise OcculusRegistrationError(
            f"Too few correspondences ({n_corr}) for RANSAC (need {ransac_n})"
        )

    # Adaptive RANSAC iteration count
    inlier_ratio = 1.0
    max_adaptive_iter = max_iterations

    for iteration in range(max_adaptive_iter):  # noqa: B007
        # Sample ransac_n correspondences
        sample = rng.choice(n_corr, size=ransac_n, replace=False)
        src_pts = source.xyz[corr_src[sample]]
        tgt_pts = target.xyz[corr_tgt[sample]]

        # Estimate transform from minimal sample
        try:
            R, t = _svd_rigid_3pt(src_pts, tgt_pts)
        except Exception:
            continue

        T = np.eye(4)
        T[:3, :3] = R
        T[:3, 3] = t

        # Count inliers over all correspondences
        src_transformed = _apply_transform(source.xyz[corr_src], T)
        tgt_pts_all = target.xyz[corr_tgt]
        dists = np.linalg.norm(src_transformed - tgt_pts_all, axis=1)
        inlier_count = int((dists <= max_correspondence_distance).sum())

        if inlier_count > best_inliers:
            best_inliers = inlier_count
            best_transform = T
            inlier_ratio = best_inliers / n_corr

            # Adaptive early exit (Hartley & Zisserman)
            max(1e-9, 1.0 - inlier_ratio)
            k_adaptive = np.log(1 - confidence) / np.log(1 - (inlier_ratio**ransac_n) + 1e-12)
            max_adaptive_iter = min(max_iterations, int(k_adaptive) + 1)

    if best_inliers == 0:
        raise OcculusRegistrationError("RANSAC failed to find any valid transformation.")

    # Refine on all inliers
    src_t = _apply_transform(source.xyz[corr_src], best_transform)
    tgt_all = target.xyz[corr_tgt]
    dists = np.linalg.norm(src_t - tgt_all, axis=1)
    inlier_mask = dists <= max_correspondence_distance

    src_inliers = source.xyz[corr_src[inlier_mask]]
    tgt_inliers = target.xyz[corr_tgt[inlier_mask]]
    R_ref, t_ref = _svd_rigid_3pt(src_inliers, tgt_inliers)
    best_transform[:3, :3] = R_ref
    best_transform[:3, 3] = t_ref

    fitness, rmse = _compute_metrics(
        _apply_transform(source.xyz, best_transform), target.xyz, max_correspondence_distance
    )
    converged = inlier_ratio >= (1.0 - confidence) * 10  # heuristic

    logger.info(
        "RANSAC: best_inliers=%d/%d fitness=%.3f rmse=%.4f",
        best_inliers,
        n_corr,
        fitness,
        rmse,
    )
    return RegistrationResult(
        transformation=best_transform,
        fitness=fitness,
        inlier_rmse=rmse,
        converged=converged,
        n_iterations=iteration + 1,
    )