The Mesh

1class Mesh:
2    def __init__(self):
3        self.Vertices = [] # list of Point3d
4        self.Faces = []    # list of MeshFace (vertex index tuples)
5        self.Normals = []  # list of Vector3d for shading

1class Mesh:
2    def __init__(self):
3        self._vertices = []
4        self._faces = []
5
6    def get_vertices(self):
7        return self._vertices
8
9    def get_faces(self):
10        return self._faces

1def Translate(self, vector: 'Vector3d'):
2    """Translates the entire mesh by modifying the vertex pool."""
3    for v in self.Vertices:
4        v.X += vector.X
5        v.Y += vector.Y
6        v.Z += vector.Z

1def Scale(self, factor: float, origin: 'Point3d'):
2    """Scales the entire vertex pool relative to an origin."""
3    for v in self.Vertices:
4        v.X = origin.X + (v.X - origin.X) * factor
5        v.Y = origin.Y + (v.Y - origin.Y) * factor
6        v.Z = origin.Z + (v.Z - origin.Z) * factor

1import math
2
3def Rotate(vertices: list[Point3d], angle_deg: float, axis: Vector3d, origin: Point3d) -> list[Point3d]:
4    """Rotates all vertices around an axis through an origin point."""
5    rad = math.radians(angle_deg)
6    cos_a = math.cos(rad)
7    sin_a = math.sin(rad)
8    
9    # Normalize axis vector
10    axis_len = (axis.X**2 + axis.Y**2 + axis.Z**2)**0.5
11    ux, uy, uz = axis.X / axis_len, axis.Y / axis_len, axis.Z / axis_len
12    
13    rotated = []
14    for v in vertices:
15        # Translate to origin
16        x, y, z = v.X - origin.X, v.Y - origin.Y, v.Z - origin.Z
17        
18        # Rodrigues' rotation formula
19        dot = x*ux + y*uy + z*uz
20        rx = x*cos_a + (uy*z - uz*y)*sin_a + ux*dot*(1 - cos_a)
21        ry = y*cos_a + (uz*x - ux*z)*sin_a + uy*dot*(1 - cos_a)
22        rz = z*cos_a + (ux*y - uy*x)*sin_a + uz*dot*(1 - cos_a)
23        
24        # Translate back
25        rotated.append(Point3d(rx + origin.X, ry + origin.Y, rz + origin.Z))
26        
27    return rotated

1import math
2
3def Deform(self, amplitude: float):
4    """Deforms the mesh vertices using a wave function."""
5    for v in self.Vertices:
6        # Displace Z based on X and Y
7        v.Z += math.sin(v.X * 5) * math.cos(v.Y * 5) * amplitude

1def RayTriangleIntersect(orig: Point3d, dir: Vector3d, v0: Point3d, v1: Point3d, v2: Point3d) -> Point3d | None:
2    """Möller–Trumbore ray-triangle intersection math."""
3    edge1 = Vector3d(v1.X - v0.X, v1.Y - v0.Y, v1.Z - v0.Z)
4    edge2 = Vector3d(v2.X - v0.X, v2.Y - v0.Y, v2.Z - v0.Z)
5    h = CrossProduct(dir, edge2)
6    a = DotProduct(edge1, h)
7    if -1e-6 < a < 1e-6:
8        return None
9    f = 1.0 / a
10    s = Vector3d(orig.X - v0.X, orig.Y - v0.Y, orig.Z - v0.Z)
11    u = f * DotProduct(s, h)
12    if u < 0.0 or u > 1.0:
13        return None
14    q = CrossProduct(s, edge1)
15    v = f * DotProduct(dir, q)
16    if v < 0.0 or u + v > 1.0:
17        return None
18    t = f * DotProduct(edge2, q)
19    if t > 1e-6:
20        return Point3d(orig.X + dir.X * t, orig.Y + dir.Y * t, orig.Z + dir.Z * t)
21    return None
22
23def MeshRaycast(mesh: Mesh, ray_orig: Point3d, ray_dir: Vector3d) -> list[Point3d]:
24    """Intersects a ray against all faces in a mesh."""
25    hits = []
26    for face in mesh.Faces:
27        v0 = mesh.Vertices[face.A]
28        v1 = mesh.Vertices[face.B]
29        v2 = mesh.Vertices[face.C]
30        pt = RayTriangleIntersect(ray_orig, ray_dir, v0, v1, v2)
31        if pt:
32            hits.append(pt)
33    return hits

1def Weld(vertices: list[Point3d], faces: list[tuple[int, ...]], tolerance: float = 0.001) -> tuple[list[Point3d], list[tuple[int, ...]]]:
2    """Merges coincident vertices within tolerance and updates face indices.
3    Enables smooth shading across face boundaries.
4    """
5    unique_vertices = []
6    index_map = {}
7    
8    for i, v in enumerate(vertices):
9        # Find if a matching vertex already exists
10        match_idx = -1
11        for j, uv in enumerate(unique_vertices):
12            dist = ((v.X - uv.X)**2 + (v.Y - uv.Y)**2 + (v.Z - uv.Z)**2)**0.5
13            if dist < tolerance:
14                match_idx = j
15                break
16        if match_idx == -1:
17            unique_vertices.append(v)
18            index_map[i] = len(unique_vertices) - 1
19        else:
20            index_map[i] = match_idx
21            
22    # Rebuild faces with updated indices
23    new_faces = []
24    for f in faces:
25        new_faces.append(tuple(index_map[idx] for idx in f))
26        
27    return unique_vertices, new_faces

1def SubdivideFaceMidpoint(vertices: list[Point3d], faces: list[tuple[int, int, int]]) -> tuple[list[Point3d], list[tuple[int, int, int]]]:
2    """Splits each triangular face into 4 smaller triangles by finding edge midpoints."""
3    new_vertices = list(vertices)
4    new_faces = []
5    midpoint_cache = {}
6    
7    def get_midpoint(i1, i2):
8        key = tuple(sorted((i1, i2)))
9        if key not in midpoint_cache:
10            p1 = new_vertices[i1]
11            p2 = new_vertices[i2]
12            mid_pt = Point3d((p1.X + p2.X)/2.0, (p1.Y + p2.Y)/2.0, (p1.Z + p2.Z)/2.0)
13            new_vertices.append(mid_pt)
14            midpoint_cache[key] = len(new_vertices) - 1
15        return midpoint_cache[key]
16        
17    for f in faces:
18        a, b, c = f[0], f[1], f[2]
19        ab = get_midpoint(a, b)
20        bc = get_midpoint(b, c)
21        ca = get_midpoint(c, a)
22        
23        new_faces.append((a, ab, ca))
24        new_faces.append((b, bc, ab))
25        new_faces.append((c, ca, bc))
26        new_faces.append((ab, bc, ca))
27        
28    return new_vertices, new_faces

1def LaplacianSmooth(vertices: list[Point3d], faces: list[tuple[int, int, int]], iterations: int, lambda_factor: float = 0.5) -> list[Point3d]:
2    """Applies Laplacian smoothing to the mesh vertices by averaging neighbor positions."""
3    curr_vertices = list(vertices)
4    
5    # 1. Build adjacency list of vertex indices
6    neighbors = {i: set() for i in range(len(vertices))}
7    for f in faces:
8        for idx in f:
9            for other_idx in f:
10                if idx != other_idx:
11                    neighbors[idx].add(other_idx)
12                    
13    # 2. Iterate smoothing steps
14    for _ in range(iterations):
15        next_vertices = []
16        for i, v in enumerate(curr_vertices):
17            n_indices = neighbors[i]
18            if not n_indices:
19                next_vertices.append(v)
20                continue
21            
22            # Compute centroid of neighbor coordinates
23            sum_x = sum(curr_vertices[n].X for n in n_indices)
24            sum_y = sum(curr_vertices[n].Y for n in n_indices)
25            sum_z = sum(curr_vertices[n].Z for n in n_indices)
26            count = len(n_indices)
27            
28            centroid_x = sum_x / count
29            centroid_y = sum_y / count
30            centroid_z = sum_z / count
31            
32            # Interpolate towards centroid
33            nx = v.X + (centroid_x - v.X) * lambda_factor
34            ny = v.Y + (centroid_y - v.Y) * lambda_factor
35            nz = v.Z + (centroid_z - v.Z) * lambda_factor
36            next_vertices.append(Point3d(nx, ny, nz))
37            
38        curr_vertices = next_vertices
39        
40    return curr_vertices

1def EdgeCollapse(vertices: list[Point3d], faces: list[tuple[int, int, int]], v_keep: int, v_remove: int) -> tuple[list[Point3d], list[tuple[int, int, int]]]:
2"""Naive implementation"    
3"""Collapses the edge from v_remove to v_keep, reducing mesh complexity."""
4    p_keep = vertices[v_keep]
5    p_remove = vertices[v_remove]
6    
7    # Place kept vertex at midpoint of collapsed edge
8    vertices[v_keep] = Point3d(
9        (p_keep.X + p_remove.X) / 2.0,
10        (p_keep.Y + p_remove.Y) / 2.0,
11        (p_keep.Z + p_remove.Z) / 2.0
12    )
13    
14    new_faces = []
15    for f in faces:
16        # Re-map references to collapsed vertex
17        a = v_keep if f[0] == v_remove else f[0]
18        b = v_keep if f[1] == v_remove else f[1]
19        c = v_keep if f[2] == v_remove else f[2]
20        
21        # Keep face only if it is not degenerate (has 3 distinct vertices)
22        if a != b and b != c and c != a:
23            new_faces.append((a, b, c))
24            
25    return vertices, new_faces

1def ComputeFaceSignedVolume(v0: Point3d, v1: Point3d, v2: Point3d) -> float:
2    """Computes signed volume of the tetrahedron formed by the face and origin."""
3    cross_x = v1.Y * v2.Z - v1.Z * v2.Y
4    cross_y = v1.Z * v2.X - v1.X * v2.Z
5    cross_z = v1.X * v2.Y - v1.Y * v2.X
6    
7    dot = v0.X * cross_x + v0.Y * cross_y + v0.Z * cross_z
8    return dot / 6.0
9
10def ComputeMeshVolume(vertices: list[Point3d], faces: list[tuple[int, int, int]]) -> float:
11    """Computes total volume of a closed watertight mesh using divergence theorem."""
12    total_volume = 0.0
13    for f in faces:
14        v0 = vertices[f[0]]
15        v1 = vertices[f[1]]
16        v2 = vertices[f[2]]
17        total_volume += ComputeFaceSignedVolume(v0, v1, v2)
18    return abs(total_volume)

1def FindNakedEdges(faces: list[tuple[int, int, int]]) -> list[tuple[int, int]]:
2    """Finds all boundary (naked) edges that belong to exactly one face."""
3    edge_counts = {}
4    
5    for f in faces:
6        # Sort indices to treat edge (u, v) same as (v, u)
7        edges = [
8            tuple(sorted((f[0], f[1]))),
9            tuple(sorted((f[1], f[2]))),
10            tuple(sorted((f[2], f[0])))
11        ]
12        for e in edges:
13            edge_counts[e] = edge_counts.get(e, 0) + 1
14            
15    # Boundary edges have exactly one face reference
16    naked_edges = [e for e, count in edge_counts.items() if count == 1]
17    return naked_edges

1def FlipMesh(faces: list[tuple[int, int, int]], normals: list[Vector3d]) -> tuple[list[tuple[int, int, int]], list[Vector3d]]:
2    """Inverts the mesh so that inside becomes outside by reversing winding and negating normals."""
3    # 1. Reverse face winding order (swap index B and C)
4    flipped_faces = [(f[0], f[2], f[1]) for f in faces]
5    
6    # 2. Negate normal vectors
7    flipped_normals = [Vector3d(-n.X, -n.Y, -n.Z) for n in normals]
8    
9    return flipped_faces, flipped_normals

1def IsWatertight(faces: list[tuple[int, int, int]]) -> bool:
2    """A mesh is watertight (closed) if every edge is shared by exactly 2 faces."""
3    edge_counts = {}
4    for f in faces:
5        edges = [
6            tuple(sorted((f[0], f[1]))),
7            tuple(sorted((f[1], f[2]))),
8            tuple(sorted((f[2], f[0])))
9        ]
10        for e in edges:
11            edge_counts[e] = edge_counts.get(e, 0) + 1
12            
13    # Check if there are any boundary (naked) edges
14    for edge, count in edge_counts.items():
15        if count != 2:
16            return False # Non-manifold or open boundary edge found
17            
18    return True

1def UnifyNormals(faces: list[tuple[int, int, int]]) -> list[tuple[int, int, int]]:
2    """Traverses mesh topology to synchronize winding order (consistent orientation)."""
3    # 1. Build edge adjacency map to track face connections
4    edge_to_faces = {}
5    for i, f in enumerate(faces):
6        edges = [(f[0], f[1]), (f[1], f[2]), (f[2], f[0])]
7        for e in edges:
8            key = tuple(sorted(e))
9            if key not in edge_to_faces:
10                edge_to_faces[key] = []
11            edge_to_faces[key].append((i, e))
12            
13    # 2. Breadth-First Search (BFS) to orient winding consistently
14    visited = [False] * len(faces)
15    queue = [0] # start with first face
16    visited[0] = True
17    
18    while queue:
19        curr_idx = queue.pop(0)
20        curr_face = faces[curr_idx]
21        curr_edges = [(curr_face[0], curr_face[1]), (curr_face[1], curr_face[2]), (curr_face[2], curr_face[0])]
22        
23        for e in curr_edges:
24            key = tuple(sorted(e))
25            neighbors = edge_to_faces.get(key, [])
26            for n_idx, n_edge in neighbors:
27                if visited[n_idx]:
28                    continue
29                
30                # If neighbor shares the edge, they must traverse it in opposite directions
31                # If they traverse in same direction, neighbor winding is inverted and must flip
32                n_face = faces[n_idx]
33                if e[0] == n_edge[0]: # same direction traversal
34                    faces[n_idx] = (n_face[0], n_face[2], n_face[1]) # swap B and C index
35                    
36                visited[n_idx] = True
37                queue.append(n_idx)
38                
39    return faces

1def ClosestVertex(test_pt: Point3d, vertices: list[Point3d]) -> Point3d:
2    """Finds the closest vertex in a mesh to a given test point."""
3    min_dist = float('inf')
4    best_pt = None
5    
6    for pt in vertices:
7        dist = ((test_pt.X - pt.X)**2 + (test_pt.Y - pt.Y)**2 + (test_pt.Z - pt.Z)**2)**0.5
8        if dist < min_dist:
9            min_dist = dist
10            best_pt = pt
11            
12    return best_pt

1def ClosestFaceCenter(test_pt: Point3d, mesh: Mesh) -> Point3d:
2    """Finds the closest face center in a mesh to a given test point."""
3    min_dist = float('inf')
4    best_pt = None
5    
6    for face in mesh.Faces:
7        a = mesh.Vertices[face.A]
8        b = mesh.Vertices[face.B]
9        c = mesh.Vertices[face.C]
10        
11        center = Point3d(
12            (a.X + b.X + c.X) / 3.0,
13            (a.Y + b.Y + c.Y) / 3.0,
14            (a.Z + b.Z + c.Z) / 3.0
15        )
16        
17        dist = ((test_pt.X - center.X)**2 + (test_pt.Y - center.Y)**2 + (test_pt.Z - center.Z)**2)**0.5
18        if dist < min_dist:
19            min_dist = dist
20            best_pt = center
21            
22    return best_pt

1def ClosestPointOnEdge(test_pt: Point3d, edges: list[tuple[Point3d, Point3d]]) -> Point3d:
2    """Finds the closest point on the mesh edge network to a test point."""
3    min_dist = float('inf')
4    best_pt = None
5    
6    for a, b in edges:
7        ab = Vector3d(b.X - a.X, b.Y - a.Y, b.Z - a.Z)
8        ap = Vector3d(test_pt.X - a.X, test_pt.Y - a.Y, test_pt.Z - a.Z)
9        
10        ab_len_sq = ab.X**2 + ab.Y**2 + ab.Z**2
11        if ab_len_sq == 0:
12            continue
13            
14        t = DotProduct(ap, ab) / ab_len_sq
15        t = max(0.0, min(1.0, t))
16        
17        pt = Point3d(a.X + t*ab.X, a.Y + t*ab.Y, a.Z + t*ab.Z)
18        dist = ((test_pt.X - pt.X)**2 + (test_pt.Y - pt.Y)**2 + (test_pt.Z - pt.Z)**2)**0.5
19        
20        if dist < min_dist:
21            min_dist = dist
22            best_pt = pt
23            
24    return best_pt

1def ClosestPointOnTriangle(p: Point3d, a: Point3d, b: Point3d, c: Point3d) -> Point3d:
2    """Finds the closest point on a triangle ABC to point P using vector projection."""
3    ab = Vector3d(b.X - a.X, b.Y - a.Y, b.Z - a.Z)
4    ac = Vector3d(c.X - a.X, c.Y - a.Y, c.Z - a.Z)
5    ap = Vector3d(p.X - a.X, p.Y - a.Y, p.Z - a.Z)
6    
7    d1 = DotProduct(ab, ap)
8    d2 = DotProduct(ac, ap)
9    if d1 <= 0.0 and d2 <= 0.0:
10        return a
11        
12    bp = Vector3d(p.X - b.X, p.Y - b.Y, p.Z - b.Z)
13    d3 = DotProduct(ab, bp)
14    d4 = DotProduct(ac, bp)
15    if d3 >= 0.0 and d4 <= d3:
16        return b
17        
18    vc = d1 * d4 - d3 * d2
19    if vc <= 0.0 and d1 >= 0.0 and d3 <= 0.0:
20        v = d1 / (d1 - d3)
21        return Point3d(a.X + ab.X * v, a.Y + ab.Y * v, a.Z + ab.Z * v)
22        
23    cp = Vector3d(p.X - c.X, p.Y - c.Y, p.Z - c.Z)
24    d5 = DotProduct(ab, cp)
25    d6 = DotProduct(ac, cp)
26    if d6 >= 0.0 and d5 <= d6:
27        return c
28        
29    vb = d5 * d2 - d1 * d6
30    if vb <= 0.0 and d2 >= 0.0 and d6 <= 0.0:
31        w = d2 / (d2 - d6)
32        return Point3d(a.X + ac.X * w, a.Y + ac.Y * w, a.Z + ac.Z * w)
33        
34    va = d3 * d6 - d5 * d4
35    if va <= 0.0 and (d4 - d3) >= 0.0 and (d5 - d6) >= 0.0:
36        w = (d4 - d3) / ((d4 - d3) + (d5 - d6))
37        return Point3d(b.X + (c.X - b.X) * w, b.Y + (c.Y - b.Y) * w, b.Z + (c.Z - b.Z) * w)
38        
39    denom = 1.0 / (va + vb + vc)
40    v = vb * denom
41    w = vc * denom
42    return Point3d(a.X + ab.X * v + ac.X * w, a.Y + ab.Y * v + ac.X * w, a.Z + ab.Z * v + ac.Z * w)
43
44def ClosestPoint(mesh: Mesh, test_pt: Point3d) -> Point3d:
45    """Finds the closest location on the entire mesh surface."""
46    min_dist = float('inf')
47    best_pt = None
48    
49    for face in mesh.Faces:
50        a = mesh.Vertices[face.A]
51        b = mesh.Vertices[face.B]
52        c = mesh.Vertices[face.C]
53        pt = ClosestPointOnTriangle(test_pt, a, b, c)
54        dist = ((test_pt.X - pt.X)**2 + (test_pt.Y - pt.Y)**2 + (test_pt.Z - pt.Z)**2)**0.5
55        if dist < min_dist:
56            min_dist = dist
57            best_pt = pt
58            
59    return best_pt