feat: compress marche

main.py

@@ -0,0 +1,429 @@
+import io
+from types import NoneType
+import obja.obja as obja
+import numpy as np
+import argparse
+
+
+def sliding_window(l: list):
+    head = l[:-1]
+    tail = l[1:]
+    return zip(head, tail)
+
+
+class MAPS(obja.Model):
+    """_summary_
+
+    Args:
+        obja (_type_): _description_
+    """
+
+    def __init__(self):
+        super().__init__()
+        self.deleted_faces = set()
+
+    def one_ring(self, index: int) -> tuple[list[int], list[int]]:
+        """ Return the corresponding 1-ring
+
+        Args:
+            index (int): index of the 1-ring's main vertex
+
+        Returns:
+            list[int]: ordered list of the 1-ring vertices
+        """
+
+        # Find the 1-ring faces
+        ring_faces, ring_face_indices = [], []
+        for face_index, face in enumerate(self.faces):
+            if face == None:
+                continue
+
+            if index in (face.a, face.b, face.c):
+                ring_faces.append(face)
+                ring_face_indices.append(face_index)
+
+        # Initialize the ring
+        start_index = ring_faces[0].a if ring_faces[0].a != index else ring_faces[0].b
+        ring = [start_index]
+        ring_faces.pop(0)
+
+        # Select the indexes of the ring in the right order
+        while len(ring_faces) > 0:
+            broke = False
+            prev_index = ring[-1]
+            for i, face in enumerate(ring_faces):
+                if prev_index in (face.a, face.b, face.c):
+                    # Found the face that correspond to the next vertex
+                    current_index = (  # select the next vertex from the face
+                        face.a if face.a != index and face.a != prev_index else
+                        face.b if face.b != index and face.b != prev_index else
+                        face.c
+                    )
+                    ring.append(current_index)
+                    ring_faces.pop(i)
+                    broke = True
+                    break
+
+            if not broke:
+                raise ValueError(
+                    f"Vertex {prev_index} is not in the remaining faces {ring_faces}. Origin {ring} on {index}")
+
+        return ring, ring_face_indices
+
+    def compute_area_curvature(self, index: int) -> tuple[float, float, int]:
+        """ Compute area and curvature the corresponding 1-ring
+
+        Args:
+            index (int): index of the 1-ring's main vertex
+
+        Returns:
+            tuple[float, float]: area and curvature
+        """
+        area_sum = 0
+        curvature_sum = 0
+        one_ring_vertices, _ = self.one_ring(index)
+
+        p1 = self.vertices[index]  # the center of the one-ring
+        for index2, index3 in sliding_window(one_ring_vertices):
+            p2 = self.vertices[index2]  # the second vertice of the triangle
+            p3 = self.vertices[index3]  # the third vertice of the triangle
+            M = np.array([  # build the matrix, used to compute the area
+                [p1[0], p2[0], p3[0]],
+                [p1[1], p2[1], p3[1]],
+                [p1[2], p2[2], p3[2]],
+            ])
+            area = abs(np.linalg.det(M) / 2)  # compute the area
+            area_sum += area
+
+            curvature = 4 * area / (  # compute the curvature
+                np.linalg.norm(p1-p2) *
+                np.linalg.norm(p2-p3) *
+                np.linalg.norm(p3-p1)
+            )
+            curvature_sum += curvature
+
+        return area_sum, curvature_sum, len(one_ring_vertices)
+
+    def compute_priority(self, lamb: float = 0.5, max_length: int = 12) -> list[float]:
+        """ Compute selection priority of vertices (0 -> hight priority ; 1 -> low priority)
+
+        Args:
+            lamb (float, optional): convex combination factor. Defaults to 0.5.
+            max_length (int, optional): 1-ring maximum length to be prioritary. Defaults to 12.
+
+        Returns:
+            list[float]: priority values
+        """
+        # Compute area and curvature for each vertex
+        areas, curvatures, ring_lengths = [], [], []
+        for i in range(len(self.vertices)):
+            if type(self.vertices[i]) != NoneType:
+                area, curvature, ring_length = self.compute_area_curvature(i)
+            else:
+                area, curvature, ring_length = -1.0, -1.0, None
+            areas.append(area)
+            curvatures.append(curvature)
+            ring_lengths.append(ring_length)
+        # Get maxes to normalize
+        max_area = max(areas)
+        max_curvature = max(curvatures)
+
+        # Compute priorities
+        priorities = []
+        for a, k, l in zip(areas, curvatures, ring_lengths):
+            if l != None and l <= max_length:
+                # Compute priority
+                priority = lamb * a / max_area + \
+                    (1.0 - lamb) * k / max_curvature
+            else:
+                # Vertex with low priority
+                priority = 2.0
+            priorities.append(priority)
+
+        return priorities
+
+    def select_vertices(self) -> list[int]:
+        """ Select vertices for the current level reduction
+
+        Returns:
+            list[int]: selected vertices
+        """
+        print("Selecting vertices...")
+
+        # Order vertices by priority
+        priorities = self.compute_priority()
+        vertices = [i[0]
+                    for i in sorted(enumerate(priorities), key=lambda p: p[1])]
+
+        selected_vertices = []
+        while len(vertices) > 0:
+            # Select prefered vertex
+            vertex = vertices.pop(0)  # remove it from remaining vertices
+            if priorities[vertex] == 2.0:
+                continue
+            selected_vertices.append(vertex)
+
+            # Remove neighbors
+            # for face in remaining_faces:
+            for face in self.faces:
+                if face == None:
+                    continue
+
+                face_vertices = (face.a, face.b, face.c)
+                if vertex in face_vertices:
+
+                    # Remove face and face's vertices form remainings
+                    # remaining_faces.remove(face)
+                    for face_vertex in face_vertices:
+                        if face_vertex in vertices:
+                            vertices.remove(face_vertex)
+
+        for ver in selected_vertices:
+            ring, _ = self.one_ring(ver)
+            for v in selected_vertices:
+                for r in ring:
+                    if r == v:
+                        print(f"fuck:{ver},{ring}")
+
+        print("Vertices selected.")
+        return selected_vertices
+
+    def project_polar(self, index: int) -> list[np.ndarray]:
+        """ Flatten the 1-ring to retriangulate
+
+        Args:
+            index (int): main vertex of the 1-ring
+
+        Returns:
+            list[np.ndarray]: list the cartesian coordinates of the flattened 1-ring projected in the plane
+        """
+        ring, _ = self.one_ring(index)
+        radius, angles = [], []
+        teta = 0.0  # cumulated angles
+        for index1, index2 in sliding_window(ring + [ring[0]]):
+            r = np.linalg.norm(self.vertices[index] - self.vertices[index1])
+            teta += self.compute_angle(index1, index, index2)  # add new angle
+            radius.append(r)
+            angles.append(teta)
+        angles = [2 * np.pi * a / teta for a in angles]  # normalize angles
+        coordinates = [np.array([r * np.cos(a), r * np.sin(a)])
+                       for r, a in zip(radius, angles)]  # parse polar to cartesian
+
+        return coordinates, ring
+
+    def compute_angle(self, i: int, j: int, k: int) -> float:
+        """ Calculate the angle defined by three points
+
+        Args:
+            i (int): previous index
+            j (int): central index
+            k (int): next index
+
+        Returns:
+            float: angle defined by the three points
+        """
+        a = self.vertices[i]
+        b = self.vertices[j]
+        c = self.vertices[k]
+        u = a - b
+        v = c - b
+        u /= np.linalg.norm(u)
+        v /= np.linalg.norm(v)
+        res = np.dot(u, v)
+
+        return np.arccos(np.clip(res, -1, 1))
+
+    def clip_ear(self, index: int) -> tuple[list[obja.Face], int]:
+        """ Retriangulate a polygon using the ear clipping algorithm
+
+        Args:
+            index (int): index of 1-ring
+
+        Returns:
+            tuple[list[obja.Face], int]: list the triangles
+        """
+        polygon, ring = self.project_polar(index)
+
+        faces = []  # the final list of faces
+
+        indices = [(local_i, global_i)
+                   for local_i, global_i in enumerate(ring)]  # remainging vertices
+
+        node_index = 0
+        while len(indices) > 2:
+            # Extract indices
+            local_i, global_i = indices[node_index - 1]
+            local_j, global_j = indices[node_index]
+            local_k, global_k = indices[node_index + 1]
+
+            # Extract verticies
+            prev_vert = polygon[local_i]
+            curr_vert = polygon[local_j]
+            next_vert = polygon[local_k]
+
+            is_convex = MAPS.is_convex(prev_vert, curr_vert, next_vert)
+            is_ear = True
+            if is_convex:  # the triangle needs to be convext to be an ear
+                # Begin with the point next to the triangle
+                test_node_index = (node_index + 2) % len(indices)
+                while indices[test_node_index][0] != local_i and is_ear:
+                    test_vert = polygon[indices[test_node_index][0]]
+                    is_ear = not MAPS.is_inside(prev_vert,
+                                                curr_vert,
+                                                next_vert,
+                                                test_vert)
+                    test_node_index = (test_node_index + 1) % len(indices)
+            else:
+                is_ear = False
+
+            if is_ear:
+                faces.append(obja.Face(global_i, global_j, global_k))
+                indices.pop(node_index)  # remove the point from the ring
+
+            node_index = (node_index + 2) % len(indices) - 1
+
+        return faces
+
+    def is_convex(prev_vert: np.ndarray, curr_vert: np.ndarray, next_vert: np.ndarray) -> bool:
+        """ Check if the angle less than pi
+
+        Args:
+            prev_vert (np.ndarray): first point
+            curr_vert (np.ndarray): middle point
+            next_vert (np.ndarray): last point
+
+        Returns:
+            bool: angle smaller than pi
+        """
+        a = prev_vert - curr_vert
+        b = next_vert - curr_vert
+        dot = a[0] * b[0] + a[1] * b[1]
+        det = a[0] * b[1] - a[1] * b[0]
+        angle = np.arctan2(det, dot)
+        if angle < 0.0:
+            angle = 2.0 * np.pi + angle
+        internal_angle = angle
+        return internal_angle >= np.pi
+
+    def is_inside(a: np.ndarray, b: np.ndarray, c: np.ndarray, p: np.ndarray) -> bool:
+        """ Check if p is in the triangle a b c
+
+        Args:
+            a (np.ndarray): point one
+            b (np.ndarray): point two
+            c (np.ndarray): point three
+            p (np.ndarray): point to check
+
+        Returns:
+            bool: if the point to check is in a b c
+        """
+        # Compute vectors
+        v0 = c - a
+        v1 = b - a
+        v2 = p - a
+
+        # Compute dot products
+        dot00 = np.dot(v0, v0)
+        dot01 = np.dot(v0, v1)
+        dot02 = np.dot(v0, v2)
+        dot11 = np.dot(v1, v1)
+        dot12 = np.dot(v1, v2)
+
+        # Compute barycentric coordinates
+        denom = dot00 * dot11 - dot01 * dot01
+        if abs(denom) < 1e-20:
+            return True
+        invDenom = 1.0 / denom
+        u = (dot11 * dot02 - dot01 * dot12) * invDenom
+        v = (dot00 * dot12 - dot01 * dot02) * invDenom
+
+        # Check if point is in triangle
+        return (u >= 0) and (v >= 0) and (u + v < 1)
+
+    def contract(self, output: io.TextIOWrapper) -> None:
+        """ Compress the 3d model
+
+        Args:
+            output (io.TextIOWrapper): Output file descriptor
+        """
+        operations = []
+
+        # while len(self.vertices) > 64:
+        for i in range(16):
+            selected_vertices = self.select_vertices()  # find the set of vertices to remove
+
+            faces_to_add, faces_to_remove = [], []
+            while len(selected_vertices) > 0:
+                vertex = selected_vertices.pop(0)  # get next vertex
+                print(f"     {len(selected_vertices)}      ", end='\r')
+
+                # Extract ring faces
+                _, ring_faces = self.one_ring(vertex)
+                faces_to_remove += ring_faces
+
+                # Apply retriangulation algorithm
+                faces = self.clip_ear(vertex)
+                faces_to_add += faces
+
+                # Edit the first faces
+                for i in range(len(faces)):
+                    operations.append(
+                        ('ef', ring_faces[i], self.faces[ring_faces[i]]))
+                    self.faces[ring_faces[i]] = faces[i]
+
+                # Remove the last faces
+                for i in range(len(faces), len(ring_faces)):
+                    operations.append(
+                        ('af', ring_faces[i], self.faces[ring_faces[i]]))
+                    self.faces[ring_faces[i]] = None
+
+                # Remove the vertex
+                operations.append(('av', vertex, self.vertices[vertex]))
+                self.vertices[vertex] = None
+
+        # REgister remaining vertices and faces
+        for i, face in enumerate(self.faces):
+            if face != None:
+                operations.append(('af', i, face))
+        for i, vertex in enumerate(self.vertices):
+            if type(vertex) != NoneType:
+                operations.append(('av', i, vertex))
+
+        # To rebuild the model, run operations in reverse order
+        operations.reverse()
+
+        # Write the result in output file
+        output_model = obja.Output(output)
+
+        for (op, index, value) in operations:
+            if op == 'av':
+                output_model.add_vertex(index, value)
+            elif op == 'af':
+                output_model.add_face(index, value)
+            elif op == 'ev':
+                output_model.edit_vertex(index, value)
+            elif op == 'ef':
+                output_model.edit_face(index, value)
+
+
+def main(args):
+    """ Run MAPS model compression
+
+    Args:
+        args (Namespace): arguments (input and output path)
+    """
+    model = MAPS()
+    model.parse_file(args.input)
+
+    with open(args.output, 'w') as output:
+        model.contract(output)
+
+
+if __name__ == '__main__':
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument('-i', '--input', type=str, required=True)
+    parser.add_argument('-o', '--output', type=str, required=True)
+    args = parser.parse_args()
+
+    main(args)

src/main.py

@@ -1,265 +0,0 @@
-from sympy import RisingFactorial
-from obja import obja
-import numpy as np
-import itertools
-
-
-def sliding_window(iterable, n=2):
-    iterables = itertools.tee(iterable, n)
-
-    for iterable, num_skipped in zip(iterables, itertools.count()):
-        for _ in range(num_skipped):
-            next(iterable, None)
-
-    return zip(*iterables)
-
-
-class MAPS(obja.Model):
-    """_summary_
-
-    Args:
-        obja (_type_): _description_
-    """
-
-    def __init__(self):
-        super().__init__()
-        self.deleted_faces = set()
-
-    def one_ring(self, index: int) -> list[int]:
-        """ Return the corresponding 1-ring
-
-        Args:
-            index (int): index of the 1-ring's main vertex
-
-        Returns:
-            list[int]: ordered list of the 1-ring vertices
-        """
-
-        # Find the 1-ring faces
-        ring_faces = []
-        for face in enumerate(self.faces):
-            if index in (face.a, face.b, face.c):
-                ring_faces.append(face)
-
-        # Initialize the ring
-        start_index = ring_faces[0][0] if ring_faces[0][0] != index else ring_faces[0][1]
-        ring = [start_index]
-        ring_faces.pop(0)
-
-        # Select the indexes of the ring in the right order
-        while len(ring_faces) > 0:
-            prev_index = ring[-1]
-            for i, face in enumerate(ring_faces):
-                if prev_index in (face.a, face.b, face.c):
-                    # Found the face that correspond to the next vertex
-                    current_index = (  # select the next vertex from the face
-                        face.a if face.a != index and face.a != prev_index else
-                        face.b if face.b != index and face.b != prev_index else
-                        face.c
-                    )
-                    ring.append(current_index)
-                    ring_faces.pop(i)
-                    break
-
-        return ring
-
-    def compute_area_curvature(self, index: int) -> tuple[float, float]:
-        """ Compute area and curvature the corresponding 1-ring
-
-        Args:
-            index (int): index of the 1-ring's main vertex
-
-        Returns:
-            tuple[float, float]: area and curvature
-        """
-        area_sum = 0
-        curvature_sum = 0
-        one_ring_vertices = self.one_ring(index)
-
-        p1 = self.vertices[index]  # the center of the one-ring
-        for index2, index3 in sliding_window(one_ring_vertices, 2):
-            p2 = self.vertices[index2]  # the second vertice of the triangle
-            p3 = self.vertices[index3]  # the third vertice of the triangle
-            M = np.array([  # build the matrix, used to compute the area
-                [p1[0], p2[0], p3[0]],
-                [p1[1], p2[1], p3[1]],
-                [p1[2], p2[2], p3[2]],
-            ])
-            area = np.linalg.det(M) / 2  # compute the area
-            area_sum += area
-
-            curvature = 4 * area / (  # compute the curvature
-                np.linalg.norm(p1-p2) *
-                np.linalg.norm(p2-p3) *
-                np.linalg.norm(p3-p1)
-            )
-            curvature_sum += curvature
-
-        return area_sum, curvature_sum, len(one_ring_vertices)
-
-    def compute_priority(self, lamb: float = 0.5, max_length: int = 12) -> list[float]:
-        """ Compute selection priority of vertices (0 -> hight priority ; 1 -> low priority)
-
-        Args:
-            lamb (float, optional): convex combination factor. Defaults to 0.5.
-            max_length (int, optional): 1-ring maximum length to be prioritary. Defaults to 12.
-
-        Returns:
-            list[float]: priority values
-        """
-        # Compute area and curvature for each vertex
-        areas, curvatures, ring_lengths = [], [], []
-        for i in range(len(self.vertices)):
-            area, curvature, ring_length = self.compute_area_curvature(i)
-            areas.append(area)
-            curvatures.append(curvature)
-            ring_lengths.append(ring_length)
-
-        # Get maxes to normalize
-        max_area = max(areas)
-        max_curvature = max(curvatures)
-
-        # Compute priorities
-        priorities = []
-        for a, k, l in zip(areas, curvatures, ring_lengths):
-            if l <= max_length:
-                # Compute priority
-                priority = lamb * a / max_area + (1 - lamb) * k / max_curvature
-            else:
-                # Vertex with low priority
-                priority = 1.0
-            priorities.append(priority)
-
-        return priorities
-
-    def select_vertices(self) -> list[int]:
-        """ Select vertices for the current level reduction
-
-        Returns:
-            list[int]: selected vertices
-        """
-        remaining_faces = self.faces.copy()
-
-        # Order vertices by priority
-        priorities = self.compute_priority()
-        vertices = [i[0]
-                    for i in sorted(enumerate(priorities), key=lambda p: p[1])]
-
-        selected_vertices = []
-        while len(vertices) > 0:
-            # Select prefered vertex
-            vertex = vertices.pop(0)  # remove it from remaining vertices
-            selected_vertices.append(vertex)
-
-            # Remove neighbors
-            for face in remaining_faces:
-                face_vertices = (face.a, face.b, face.c)
-                if vertex in face_vertices:
-
-                    # Remove face and face's vertices form remainings
-                    remaining_faces.remove(face)
-                    for face_vertex in face_vertices:
-                        if face_vertex in vertices:
-                            vertices.remove(face_vertex)
-
-        return selected_vertices
-
-    def project_polar(self, index: int) -> list[np.ndarray]:
-        """ Flatten the 1-ring to retriangulate
-
-        Args:
-            index (int): main vertex of the 1-ring
-
-        Returns:
-            list[np.ndarray]: list the cartesian coordinates of the flattened 1-ring projected in the plane
-        """
-        ring = self.one_ring(index)
-        radius, angles = [], []
-        teta = 0.0  # cumulated angles
-        for index1, index2 in sliding_window(ring, 2):
-            r = np.linalg.norm(self.vertices[index], self.vertices[index1])
-            teta += self.compute_angle(index1, index, index2)  # add new angle
-            radius.append(r)
-            angles.append(teta)
-        angles = [2 * np.pi * a / teta for a in angles]  # normalize angles
-        coordinates = [np.array([r * np.cos(a), r * np.sin(a)])
-                       for r, a in zip(radius, angles)]  # parse polar to cartesian
-
-        return coordinates
-
-    def compute_angle(self, i: int, j: int, k: int) -> float:
-        """ Calculate the angle defined by three points
-
-        Args:
-            i (int): previous index
-            j (int): central index
-            k (int): next index
-
-        Returns:
-            float: angle defined by the three points
-        """
-        a = self.vertices[i]
-        b = self.vertices[j]
-        c = self.vertices[k]
-        u = a - b
-        v = c - b
-        u /= np.linalg.norm(u)
-        v /= np.linalg.norm(v)
-        res = np.dot(u, v)
-        
-        return np.arccos(res)
-
-    # def contract(self, output):
-    #     """
-    #     Decimates the model stupidly, and write the resulting obja in output.
-    #     """
-    #     operations = []
-
-    #     for (vertex_index, vertex) in enumerate(self.vertices):
-    #         operations.append(('ev', vertex_index, vertex + 0.25))
-
-    #     # Iterate through the vertex
-    #     for (vertex_index, vertex) in enumerate(self.vertices):
-
-    #         # Iterate through the faces
-    #         for (face_index, face) in enumerate(self.faces):
-
-    #             # Delete any face related to this vertex
-    #             if face_index not in self.deleted_faces:
-    #                 if vertex_index in [face.a,face.b,face.c]:
-    #                     self.deleted_faces.add(face_index)
-    #                     # Add the instruction to operations stack
-    #                     operations.append(('face', face_index, face))
-
-    #         # Delete the vertex
-    #         operations.append(('vertex', vertex_index, vertex))
-
-    #     # To rebuild the model, run operations in reverse order
-    #     operations.reverse()
-
-    #     # Write the result in output file
-    #     output_model = obja.Output(output, random_color=True)
-
-    #     for (ty, index, value) in operations:
-    #         if ty == "vertex":
-    #             output_model.add_vertex(index, value)
-    #         elif ty == "face":
-    #             output_model.add_face(index, value)
-    #         else:
-    #             output_model.edit_vertex(index, value)
-
-
-def main():
-    """
-    Runs the program on the model given as parameter.
-    """
-    model = MAPS()
-    model.parse_file('example/suzanne.obj')
-
-    # with open("example/suzanne.obja", "w") as output:
-    #     model.contract(output)
-
-
-if __name__ == '__main__':
-    np.seterr(invalid='raise')
-    main()