From d1ccfc40be6f1077c1a6ceccd562644d2f795958 Mon Sep 17 00:00:00 2001 From: gdamms Date: Mon, 3 Oct 2022 13:31:38 +0200 Subject: [PATCH] feat: compress marche --- main.py | 429 ++++++++++++++++++++++++++++++++++++++++++++++++++++ src/main.py | 265 -------------------------------- 2 files changed, 429 insertions(+), 265 deletions(-) create mode 100644 main.py delete mode 100644 src/main.py diff --git a/main.py b/main.py new file mode 100644 index 0000000..3b7cde8 --- /dev/null +++ b/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) diff --git a/src/main.py b/src/main.py deleted file mode 100644 index 7439efe..0000000 --- a/src/main.py +++ /dev/null @@ -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()