feat: ca marche 1 PE

main.py

@@ -1,10 +1,12 @@
 import io
-from types import NoneType
+from math import floor
 import obja.obja as obja
 import numpy as np
 import argparse
+from time import time
+
+from rich.progress import track, Progress
 
-from rich.progress import track
 
 def cot(x: float):
     sin_x = np.sin(x)
@@ -20,6 +22,55 @@ def sliding_window(l: list, n: int = 2):
     return res
 
 
+class Edge:
+    def __init__(self, a, b):
+        self.a = min(a, b)
+        self.b = max(a, b)
+
+        self.face1 = None
+        self.face2 = None
+
+        self.fold = 0.0
+        self.curvature = 0.0
+
+    def __eq__(self, __o: object) -> bool:
+        return self.a == __o.a and self.b == __o.b
+
+
+class Face:
+    def __init__(self, a, b, c):
+        self.a = a
+        self.b = b
+        self.c = c
+
+        self.normal = np.zeros(3)
+
+    def to_obja(self):
+        return obja.Face(self.a, self.b, self.c)
+
+    def __eq__(self, __o: object) -> bool:
+        if __o is None:
+            return False
+
+        return self.a == __o.a and self.b == __o.b and self.c == __o.c
+
+
+class Vertex:
+    def __init__(self, pos):
+        self.pos = pos
+
+        self.vertex_ring = []
+        self.face_ring = []
+
+        self.normal = np.zeros(3)
+
+        self.area = 0.0
+        self.curvature = 0.0
+
+    def to_obja(self):
+        return self.pos
+
+
 class MAPS(obja.Model):
     """_summary_
 
@@ -29,7 +80,100 @@ class MAPS(obja.Model):
 
     def __init__(self):
         super().__init__()
-        self.deleted_faces = set()
+
+    def parse_file(self, path):
+        super().parse_file(path)
+        for i, vertex in enumerate(self.vertices):
+            self.vertices[i] = Vertex(vertex)
+        for i, face in enumerate(self.faces):
+            self.faces[i] = Face(face.a, face.b, face.c)
+
+    def update(self):
+        self.update_edges()
+        self.update_rings()
+        self.update_normals()
+        self.update_area_curvature()
+
+    def update_edges(self):
+        self.edges = {}
+
+        remaining_faces = self.faces.copy()
+        while None in remaining_faces:
+            remaining_faces.remove(None)
+
+        for face in track(self.faces, description='Update edges'):
+            if face is None:
+                continue
+
+            for a, b in sliding_window([face.a, face.b, face.c], n=2):
+                new_edge = Edge(a, b)
+                if self.edges.get(f"{new_edge.a}:{new_edge.b}") is None:
+                    new_edge.face1 = face
+
+                    if face in remaining_faces:
+                        remaining_faces.remove(face)
+
+                    for face2 in remaining_faces:
+                        face2_vertices = (face2.a, face2.b, face2.c)
+                        if not (a in face2_vertices and b in face2_vertices):
+                            continue
+                        new_edge.face2 = face2
+                        break
+                    if new_edge.face2 is None:
+                        print('ooooooooooooooooooooooo')
+                    self.edges[f"{new_edge.a}:{new_edge.b}"] = new_edge
+
+    def update_rings(self):
+        for i, vertex in enumerate(self.vertices):
+            if vertex is None:
+                continue
+
+            vertex_ring, face_ring = self.one_ring(i)
+            vertex.vertex_ring = vertex_ring
+            vertex.face_ring = face_ring
+
+    def update_area_curvature(self):
+        for i, vertex in enumerate(self.vertices):
+            if vertex is None:
+                continue
+
+            area, curvature = self.compute_area_curvature(i)
+            vertex.area = area
+            vertex.curvature = curvature
+
+        self.feature_edges = []
+        for edge in self.edges.values():
+            edge.fold = np.dot(edge.face1.normal, edge.face2.normal)
+
+            if edge.fold < 0.5:
+                self.feature_edges.append(edge)
+
+    def update_normals(self):
+        for face in self.faces:
+            if face is None:
+                continue
+
+            p1 = self.vertices[face.a].pos
+            p2 = self.vertices[face.b].pos
+            p3 = self.vertices[face.c].pos
+            u = p2 - p1
+            v = p3 - p1
+            n = np.cross(u, v)
+            n /= np.linalg.norm(n)
+
+            face.normal = n
+
+            self.vertices[face.a].normal += n
+            self.vertices[face.b].normal += n
+            self.vertices[face.c].normal += n
+
+        for vertex in self.vertices:
+            if vertex is None:
+                continue
+
+            norm = np.linalg.norm(vertex.normal)
+            if norm != 0:
+                vertex.normal /= norm
 
     def one_ring(self, index: int) -> tuple[list[int], list[int]]:
         """ Return the corresponding 1-ring
@@ -40,11 +184,13 @@ class MAPS(obja.Model):
         Returns:
             list[int]: ordered list of the 1-ring vertices
         """
+        if self.vertices[index] is None:
+            return None, None
 
         # Find the 1-ring faces
         ring_faces, ring_face_indices = [], []
         for face_index, face in enumerate(self.faces):
-            if face == None:
+            if face is None:
                 continue
 
             if index in (face.a, face.b, face.c):
@@ -52,7 +198,9 @@ class MAPS(obja.Model):
                 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
+        start_index = (ring_faces[0].a if ring_faces[0].a != index and ring_faces[0].c != index else
+                       ring_faces[0].b if ring_faces[0].a != index and ring_faces[0].b != index else
+                       ring_faces[0].c)
         ring = [start_index]
         ring_faces.pop(0)
 
@@ -79,7 +227,7 @@ class MAPS(obja.Model):
 
         return ring, ring_face_indices
 
-    def compute_area_curvature(self, index: int) -> tuple[float, float, int]:
+    def compute_area_curvature(self, index: int) -> tuple[float, float]:
         """ Compute area and curvature the corresponding 1-ring
 
         Args:
@@ -88,15 +236,20 @@ class MAPS(obja.Model):
         Returns:
             tuple[float, float]: area and curvature
         """
-        area_sum = 0
-        laplace_sum = 0
-        one_ring_vertices, _ = self.one_ring(index)
+        if self.vertices[index] is None:
+            return None, None
 
-        teta = 0.0
-        p1 = self.vertices[index]  # the center of the one-ring
-        for index1, index2, index3 in sliding_window(one_ring_vertices, n=3):
-            p2 = self.vertices[index1]  # the second vertice of the triangle
-            p3 = self.vertices[index2]  # the third vertice of the triangle
+        ring = self.vertices[index].vertex_ring
+        p1 = self.vertices[index].pos
+        n1 = self.vertices[index].normal
+
+        area_sum = 0
+        curvature = 0
+        for index1, index2 in sliding_window(ring, n=2):
+            # the second vertice of the triangle
+            p2 = self.vertices[index1].pos
+            p3 = self.vertices[index2].pos  # the third vertice of the triangle
+            n2 = self.vertices[index1].normal
             M = np.array([  # build the matrix, used to compute the area
                 [p1[0], p2[0], p3[0]],
                 [p1[1], p2[1], p3[1]],
@@ -105,24 +258,17 @@ class MAPS(obja.Model):
             area = abs(np.linalg.det(M) / 2)  # compute the area
             area_sum += area
 
-            teta += self.compute_angle(index1, index, index2)
+            edge_curvature = np.dot(n2 - n1, p2 - p1) / \
+                np.linalg.norm(p2 - p1)**2
+            edge_curvature = abs(edge_curvature)
+            edge_key = f"{min(index, index1)}:{max(index, index1)}"
+            self.edges[edge_key].curvature = edge_curvature
 
-            laplace = self.compute_laplace(index, index1, index2, index3)
-            laplace_sum += laplace
+            curvature += edge_curvature
 
-        K = (2 * np.pi - teta) / area * 3
-        H = np.linalg.norm(laplace_sum) / 4 / area * 3
-        curvature = abs(H - np.sqrt(H*H - K)) + abs(H + np.sqrt(H*H - K))
-        curvature = K
+        curvature /= len(ring)
 
-        return area_sum, curvature, len(one_ring_vertices)
-
-    def compute_laplace(self, i: int, j: int, a: int, b: int) -> np.ndarray:
-        alpha = self.compute_angle(i, a, j)
-        beta = self.compute_angle(i, b, j)
-        cot_sum = cot(alpha) + cot(beta)
-        vec = self.vertices[j] - self.vertices[i]
-        return cot_sum * vec
+        return area_sum, curvature
 
     def compute_priority(self, lamb: float = 0.0, max_length: int = 12) -> list[float]:
         """ Compute selection priority of vertices (0.0 -> hight priority ; 1.0 -> low priority)
@@ -134,33 +280,25 @@ class MAPS(obja.Model):
         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)
+        max_area = max(
+            [vertex.area for vertex in self.vertices if vertex is not None])
+        max_curvature = max(
+            [vertex.curvature for vertex in self.vertices if vertex is not None])
 
         # Compute priorities
         priorities = []
-        for a, k, l in zip(areas, curvatures, ring_lengths):
-            if l != None and l <= max_length:
+        for vertex in self.vertices:
+            if vertex is not None and len(vertex.vertex_ring) < max_length:
                 # Compute priority
-                priority = lamb * a / max_area + \
-                    (1.0 - lamb) * k / max_curvature
+                priority = (
+                    lamb * vertex.area / max_area +
+                    (1.0 - lamb) * vertex.curvature / max_curvature
+                )
             else:
                 # Vertex with low priority
                 priority = 2.0
             priorities.append(priority)
 
-        # return np.random.rand(len(priorities))
         return priorities
 
     def select_vertices(self) -> list[int]:
@@ -169,38 +307,50 @@ class MAPS(obja.Model):
         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:
+        with Progress() as progress:
+            task = progress.add_task('Select vertices', total=len(vertices))
+            while not progress.finished:
+                # Select prefered vertex
+                vertex = vertices.pop(0)  # remove it from remaining vertices
+                progress.advance(task)
+
+                if priorities[vertex] == 2.0:
                     continue
 
-                face_vertices = (face.a, face.b, face.c)
-                if vertex in face_vertices:
+                incident_count = 0
+                for feature_edge in self.feature_edges:
+                    if vertex in (feature_edge.a, feature_edge.b):
+                        incident_count += 1
 
-                    # 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)
+                if incident_count > 2:
+                    continue
 
-        print("Vertices selected.")
-        return selected_vertices
+                selected_vertices.append(vertex)
+
+                # Remove neighbors
+                # for face in remaining_faces:
+                for face in self.faces:
+                    if face is 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)
+                                progress.advance(task)
+
+        return selected_vertices[:floor(1.0 * len(selected_vertices))]
 
     def project_polar(self, index: int) -> list[np.ndarray]:
         """ Flatten the 1-ring to retriangulate
@@ -211,11 +361,12 @@ class MAPS(obja.Model):
         Returns:
             list[np.ndarray]: list the cartesian coordinates of the flattened 1-ring projected in the plane
         """
-        ring, _ = self.one_ring(index)
+        ring = self.vertices[index].vertex_ring
         radius, angles = [], []
         teta = 0.0  # cumulated angles
         for index1, index2 in sliding_window(ring):
-            r = np.linalg.norm(self.vertices[index] - self.vertices[index1])
+            r = np.linalg.norm(
+                self.vertices[index].pos - self.vertices[index1].pos)
             teta += self.compute_angle(index1, index, index2)  # add new angle
             radius.append(r)
             angles.append(teta)
@@ -236,9 +387,9 @@ class MAPS(obja.Model):
         Returns:
             float: angle defined by the three points
         """
-        a = self.vertices[i]
-        b = self.vertices[j]
-        c = self.vertices[k]
+        a = self.vertices[i].pos
+        b = self.vertices[j].pos
+        c = self.vertices[k].pos
         u = a - b
         v = c - b
         u /= np.linalg.norm(u)
@@ -256,45 +407,88 @@ class MAPS(obja.Model):
         Returns:
             tuple[list[obja.Face], int]: list the triangles
         """
-        polygon, ring = self.project_polar(index)
+        polygon_, ring_ = self.project_polar(index)
+        
+        main_v = []
+        for i, r in enumerate(ring_):
+            for feature_edge in self.feature_edges:
+                feat_edge_vertices = (feature_edge.a, feature_edge.b)
+                if r in feat_edge_vertices and index in feat_edge_vertices:
+                    main_v.append(i)
 
+        if len(main_v) < 2:
+
+            polygons_rings = [(polygon_, ring_)]
+
+        else:
+
+            v1 = ring_[main_v[0]]
+            v2 = ring_[main_v[1]]
+            ring1, ring2 = [], []
+            polygon1, polygon2, = [], []
+
+            start = ring_.index(v1)
+            while ring_[start] != v2:
+                ring1.append(ring_[start])
+                polygon1.append(polygon_[start])
+                start += 1
+                start %= len(ring_)
+            ring1.append(ring_[start])
+            polygon1.append(polygon_[start])
+
+            start = ring_.index(v2)
+            while ring_[start] != v1:
+                ring2.append(ring_[start])
+                polygon2.append(polygon_[start])
+                start += 1
+                start %= len(ring_)
+            ring2.append(ring_[start])
+            polygon2.append(polygon_[start])
+
+            polygons_rings = [(polygon1, ring1), (polygon2, ring2)]
+        
         faces = []  # the final list of faces
 
-        indices = [(local_i, global_i)
-                   for local_i, global_i in enumerate(ring)]  # remainging vertices
+        for polygon, ring in polygons_rings:
 
-        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]
+            indices = [(local_i, global_i)
+                       for local_i, global_i in enumerate(ring)]  # remainging vertices
 
-            # Extract verticies
-            prev_vert = polygon[local_i]
-            curr_vert = polygon[local_j]
-            next_vert = polygon[local_k]
+            node_index = 0
+            cycle_counter = 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]
 
-            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
+                # Extract verticies
+                prev_vert = polygon[local_i]
+                curr_vert = polygon[local_j]
+                next_vert = polygon[local_k]
 
-            if is_ear:
-                faces.append(obja.Face(global_i, global_j, global_k))
-                indices.pop(node_index)  # remove the point from the ring
+                is_convex = MAPS.is_convex(prev_vert, curr_vert, next_vert)
+                is_ear = True
+                if is_convex or cycle_counter > len(indices):  # 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
+                    cycle_counter += 1
 
-            node_index = (node_index + 2) % len(indices) - 1
+                if is_ear:
+                    faces.append(Face(global_i, global_j, global_k))
+                    indices.pop(node_index)  # remove the point from the ring
+                    cycle_counter = 0
+
+                node_index = (node_index + 2) % len(indices) - 1
 
         return faces
 
@@ -353,26 +547,36 @@ class MAPS(obja.Model):
 
         # Check if point is in triangle
         return (u >= 0) and (v >= 0) and (u + v < 1)
-    
+
     def truc(self, output):
+        self.update()
         priorities = self.compute_priority()
-        min_p = min(priorities)
-        priorities = [p - min_p for p in priorities]
-        max_p = max(priorities)
-        # colors = [priorities[face.a] + priorities[face.b] + priorities[face.c] for face in self.faces]
-        # min_c = min(colors)
-        # colors = [c - min_c for c in colors]
-        # max_c = max(colors)
+        # min_p = min(priorities)
+        # priorities = [p - min_p for p in priorities]
+        # max_p = max(priorities)
+        colors = [priorities[face.a] + priorities[face.b] +
+                  priorities[face.c] if face is not None else 0.0 for face in self.faces]
+        min_c = min(colors)
+        colors = [c - min_c for c in colors]
+        max_c = max(colors)
         operations = []
         for i, face in enumerate(self.faces):
             if face != None:
-                # r, g, b = colors[i] / max_c, 1.0, 1.0
-                # operations.append(('fc', i, (r, g, b)))
+                r, g, b = colors[i] / max_c, 1.0, 1.0
+                c = 0
+                for x in (face.a, face.b, face.c):
+                    for feature_edge in self.feature_edges:
+                        if x in feature_edge:
+                            c += 1
+                            break
+                if c > 1:
+                    r, g, b = 1.0, 0.0, 0.0
+                operations.append(('fc', i, (r, g, b), 0, 0, 0))
                 operations.append(('af', i, face, 0, 0, 0))
         for i, vertex in enumerate(self.vertices):
-            if type(vertex) != NoneType:
-                r, g, b = priorities[i] / max_p , 1.0, 1.0
-                operations.append(('av', i, vertex, r, g, b))
+            if vertex is None:
+                # r, g, b = priorities[i] / max_p , 1.0, 1.0
+                operations.append(('av', i, vertex, 1.0, 1.0, 1.0))
         operations.reverse()
 
         # Write the result in output file
@@ -396,8 +600,7 @@ class MAPS(obja.Model):
                     file=output
                 )
 
-
-    def contract(self, output: io.TextIOWrapper) -> None:
+    def compress(self, output: io.TextIOWrapper, final_only: bool) -> None:
         """ Compress the 3d model
 
         Args:
@@ -406,41 +609,45 @@ class MAPS(obja.Model):
         operations = []
 
         # while len(self.vertices) > 64:
-        for i in range(1):
+        for _ in range(2):
+            self.update()
             selected_vertices = self.select_vertices()  # find the set of vertices to remove
 
-            for vertex in track(selected_vertices, description="compression"):
-                # print(f"     {len(selected_vertices)}      ", end='\r')
+            for v_index in track(selected_vertices, description="Compression"):
 
                 # Extract ring faces
-                _, ring_faces = self.one_ring(vertex)
+                ring_faces = self.vertices[v_index].face_ring
 
                 # Apply retriangulation algorithm
-                faces = self.clip_ear(vertex)
+                faces = self.clip_ear(v_index)
 
                 # Edit the first faces
                 for i in range(len(faces)):
-                    # operations.append(
-                    #     ('ef', ring_faces[i], self.faces[ring_faces[i]]))
+                    if not final_only:
+                        operations.append(
+                            ('ef', ring_faces[i], self.faces[ring_faces[i]].to_obja()))
                     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]]))
+                    if not final_only:
+                        operations.append(
+                            ('af', ring_faces[i], self.faces[ring_faces[i]].to_obja()))
                     self.faces[ring_faces[i]] = None
 
                 # Remove the vertex
-                # operations.append(('av', vertex, self.vertices[vertex]))
-                self.vertices[vertex] = None
+                if not final_only:
+                    operations.append(
+                        ('av', v_index, self.vertices[v_index].to_obja()))
+                self.vertices[v_index] = 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))
+            if face is not None:
+                operations.append(('af', i, face.to_obja()))
+        for i, v_index in enumerate(self.vertices):
+            if v_index is not None:
+                operations.append(('av', i, v_index.to_obja()))
 
         # To rebuild the model, run operations in reverse order
         operations.reverse()
@@ -477,7 +684,9 @@ def main(args):
     model.parse_file(args.input)
 
     with open(args.output, 'w') as output:
-        model.truc(output)
+        model.compress(output, args.final)
+    # with open(args.output, 'w') as output:
+    #     model.truc(output)
 
 
 if __name__ == '__main__':
@@ -485,6 +694,7 @@ if __name__ == '__main__':
     parser = argparse.ArgumentParser()
     parser.add_argument('-i', '--input', type=str, required=True)
     parser.add_argument('-o', '--output', type=str, required=True)
+    parser.add_argument('-f', '--final', type=bool, default=False)
     args = parser.parse_args()
 
     main(args)