Čech complex.

cech.py

import numpy as np
from scipy.spatial import Delaunay

def compute_cech_complex(points, radius):
    """
    Compute the Čech complex from a set of points.

    Args:
        points (numpy.ndarray): An array of shape (n, m) where n is the number of points and m is the dimensionality.
        radius (float): The radius parameter for Čech complex.

    Returns:
        cech_complex (list of sets): A list of sets representing the Čech complex.
    """
    n = len(points)
    cech_complex = [set() for _ in range(n)]

    # Create a Delaunay triangulation
    tri = Delaunay(points)

    for i in range(n):
        for j in range(i + 1, n):
            # Check if the distance between points i and j is less than the radius
            if np.linalg.norm(points[i] - points[j]) <= radius:
                # For each vertex in the Delaunay simplex containing points i and j, add it to the Čech complex
                simplex_indices = tri.find_simplex([points[i], points[j]])
                for simplex_index in simplex_indices:
                    cech_complex[i].add(simplex_index)
                    cech_complex[j].add(simplex_index)

    return cech_complex


def plot_cech_complex(points, radius):
    # Compute the Čech complex
    cech_complex = compute_cech_complex(points, radius)

    # Create a 3D plot
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    # Plot points
    ax.scatter(points[:, 0], points[:, 1], points[:, 2], c='b', marker='o', label='Points')

    # Project and plot simplices onto the xy-plane
    for i, neighbors in enumerate(cech_complex):
        for j in neighbors:
            simplex = points[Delaunay(points).simplices[j]]
            projected_simplex = simplex[:, :2]  # Projection onto the xy-plane
            simplex_edges = list(zip(projected_simplex, np.roll(projected_simplex, -1, axis=0)))
            for edge in simplex_edges:
                ax.plot(*zip(*edge), 'r-')

    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    plt.title('Projection of Čech Complex in 3D (Projected onto the XY-plane)')
    plt.legend()

    plt.show()

 # Generate random 3D points
np.random.seed(0)
points = np.random.rand(20, 3)
radius = 0.2  # Adjust the radius as needed

plot_cech_complex(points, radius)