normals¶

`occulus.normals` ¶

Surface normal estimation and orientation for point clouds.

Normal vectors are required for several downstream operations including Poisson surface reconstruction and point-to-plane ICP registration.

Available functions

:func:estimate_normals — PCA-based normal estimation from local neighbourhoods
:func:orient_normals_to_viewpoint — flip normals to face a known viewpoint

All implementations use pure NumPy and SciPy. No optional dependencies required.

`estimate_normals(cloud, *, radius=None, max_nn=30)` ¶

Estimate surface normals via PCA on local neighbourhoods.

For each point, the max_nn nearest neighbours within radius are collected and a PCA is performed. The eigenvector corresponding to the smallest eigenvalue of the local covariance matrix is taken as the surface normal. Normals are not oriented — use :func:orient_normals_to_viewpoint to flip them toward a known viewpoint.

Parameters:

Name	Type	Description	Default
`cloud`	`PointCloud`	Input point cloud. Returned cloud will have the `normals` attribute populated with unit vectors.	required
`radius`	`float \| None`	Neighbourhood search radius. If `None`, a radius equal to twice the mean nearest-neighbour distance is used.	`None`
`max_nn`	`int`	Maximum number of neighbours to use per point, by default 30. More neighbours → smoother normals.	`30`

Returns:

Type	Description
`PointCloud`	A new cloud of the same concrete type with `normals` set.

Raises:

Type	Description
`OcculusValidationError`	If the cloud has fewer than 3 points (PCA is not meaningful).

Source code in src/occulus/normals/__init__.py

def estimate_normals(
    cloud: PointCloud,
    *,
    radius: float | None = None,
    max_nn: int = 30,
) -> PointCloud:
    """Estimate surface normals via PCA on local neighbourhoods.

    For each point, the ``max_nn`` nearest neighbours within ``radius`` are
    collected and a PCA is performed. The eigenvector corresponding to the
    smallest eigenvalue of the local covariance matrix is taken as the surface
    normal. Normals are not oriented — use :func:`orient_normals_to_viewpoint`
    to flip them toward a known viewpoint.

    Parameters
    ----------
    cloud : PointCloud
        Input point cloud. Returned cloud will have the ``normals`` attribute
        populated with unit vectors.
    radius : float | None, optional
        Neighbourhood search radius. If ``None``, a radius equal to twice the
        mean nearest-neighbour distance is used.
    max_nn : int, optional
        Maximum number of neighbours to use per point, by default 30.
        More neighbours → smoother normals.

    Returns
    -------
    PointCloud
        A new cloud of the same concrete type with ``normals`` set.

    Raises
    ------
    OcculusValidationError
        If the cloud has fewer than 3 points (PCA is not meaningful).
    """
    if cloud.n_points < 3:
        raise OcculusValidationError(
            f"estimate_normals requires at least 3 points, got {cloud.n_points}"
        )
    if max_nn < 3:
        raise OcculusValidationError(f"max_nn must be at least 3, got {max_nn}")

    xyz = cloud.xyz
    tree = KDTree(xyz)

    # Auto-radius: 2× the mean of the nearest-neighbour distances
    if radius is None:
        nn_distances, _ = tree.query(xyz, k=2, workers=-1)
        radius = float(nn_distances[:, 1].mean() * 2.0)
        logger.debug("estimate_normals: auto radius=%.6f", radius)

    normals = np.zeros_like(xyz)

    # Query neighbours for all points at once then process individually
    neighbour_lists = tree.query_ball_point(xyz, r=radius, workers=-1)

    for i, neighbours in enumerate(neighbour_lists):
        if len(neighbours) < 3:
            # Degenerate neighbourhood — use z-up as fallback
            normals[i] = [0.0, 0.0, 1.0]
            continue

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

        pts = xyz[neighbours]
        centred = pts - pts.mean(axis=0)
        cov = centred.T @ centred / len(pts)
        _eigenvalues, eigenvectors = np.linalg.eigh(cov)
        # Smallest eigenvalue → most perpendicular direction = normal
        normals[i] = eigenvectors[:, 0]

    # Normalise to unit length (guard against zero-length normals)
    norms = np.linalg.norm(normals, axis=1, keepdims=True)
    norms = np.where(norms < 1e-12, 1.0, norms)
    normals = normals / norms

    logger.info("estimate_normals: computed normals for %d points", cloud.n_points)

    return _copy_with_normals(cloud, normals)

`orient_normals_to_viewpoint(cloud, viewpoint)` ¶

Flip normals so that they face a known viewpoint (scanner position).

A normal is flipped if the dot product with the vector from the point to the viewpoint is negative (i.e., the normal points away).