From 0c79065d1142df13f08db364144ef5049ee56c74 Mon Sep 17 00:00:00 2001 From: gdamms Date: Mon, 17 Oct 2022 15:51:37 +0200 Subject: [PATCH] feat: ca marche 1 PE --- main.py | 482 ++++++++++++++++++++++++++++++++++++++++---------------- 1 file changed, 346 insertions(+), 136 deletions(-) diff --git a/main.py b/main.py index ac05733..7699bf4 100644 --- a/main.py +++ b/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)