diff --git a/main.py b/main.py
index 9ec2584..ac05733 100644
--- a/main.py
+++ b/main.py
@@ -4,11 +4,20 @@ import obja.obja as obja
 import numpy as np
 import argparse
 
+from rich.progress import track
 
-def sliding_window(l: list):
-    head = l[:-1]
-    tail = l[1:]
-    return zip(head, tail)
+def cot(x: float):
+    sin_x = np.sin(x)
+    if sin_x == 0:
+        return 1e16
+    return np.cos(x) / sin_x
+
+
+def sliding_window(l: list, n: int = 2):
+    k = n - 1
+    l2 = l + [l[i] for i in range(k)]
+    res = [(x for x in l2[i:i+n]) for i in range(len(l2)-k)]
+    return res
 
 
 class MAPS(obja.Model):
@@ -80,13 +89,14 @@ class MAPS(obja.Model):
             tuple[float, float]: area and curvature
         """
         area_sum = 0
-        curvature_sum = 0
+        laplace_sum = 0
         one_ring_vertices, _ = self.one_ring(index)
 
+        teta = 0.0
         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
+        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
             M = np.array([  # build the matrix, used to compute the area
                 [p1[0], p2[0], p3[0]],
                 [p1[1], p2[1], p3[1]],
@@ -95,17 +105,27 @@ class MAPS(obja.Model):
             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
+            teta += self.compute_angle(index1, index, index2)
 
-        return area_sum, curvature_sum, len(one_ring_vertices)
+            laplace = self.compute_laplace(index, index1, index2, index3)
+            laplace_sum += laplace
 
-    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)
+        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
+
+        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
+
+    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)
 
         Args:
             lamb (float, optional): convex combination factor. Defaults to 0.5.
@@ -140,6 +160,7 @@ class MAPS(obja.Model):
                 priority = 2.0
             priorities.append(priority)
 
+        # return np.random.rand(len(priorities))
         return priorities
 
     def select_vertices(self) -> list[int]:
@@ -178,13 +199,6 @@ class MAPS(obja.Model):
                         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
 
@@ -200,7 +214,7 @@ class MAPS(obja.Model):
         ring, _ = self.one_ring(index)
         radius, angles = [], []
         teta = 0.0  # cumulated angles
-        for index1, index2 in sliding_window(ring + [ring[0]]):
+        for index1, index2 in sliding_window(ring):
             r = np.linalg.norm(self.vertices[index] - self.vertices[index1])
             teta += self.compute_angle(index1, index, index2)  # add new angle
             radius.append(r)
@@ -339,6 +353,49 @@ class MAPS(obja.Model):
 
         # Check if point is in triangle
         return (u >= 0) and (v >= 0) and (u + v < 1)
+    
+    def truc(self, output):
+        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)
+        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)))
+                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))
+        operations.reverse()
+
+        # Write the result in output file
+        output_model = obja.Output(output)
+
+        for (op, index, value, r, g, b) in operations:
+            if op == 'av':
+                output_model.add_vertex_rgb(index, value, r, g, b)
+            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)
+            elif op == 'fc':
+                print('fc {} {} {} {}'.format(
+                    len(output_model.face_mapping),
+                    value[0],
+                    value[1],
+                    value[2]),
+                    file=output
+                )
+
 
     def contract(self, output: io.TextIOWrapper) -> None:
         """ Compress the 3d model
@@ -349,12 +406,11 @@ class MAPS(obja.Model):
         operations = []
 
         # while len(self.vertices) > 64:
-        for i in range(16):
+        for i in range(1):
             selected_vertices = self.select_vertices()  # find the set of vertices to remove
 
-            while len(selected_vertices) > 0:
-                vertex = selected_vertices.pop(0)  # get next vertex
-                print(f"     {len(selected_vertices)}      ", end='\r')
+            for vertex in track(selected_vertices, description="compression"):
+                # print(f"     {len(selected_vertices)}      ", end='\r')
 
                 # Extract ring faces
                 _, ring_faces = self.one_ring(vertex)
@@ -364,18 +420,18 @@ class MAPS(obja.Model):
 
                 # Edit the first faces
                 for i in range(len(faces)):
-                    operations.append(
-                        ('ef', ring_faces[i], self.faces[ring_faces[i]]))
+                    # 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]]))
+                    # 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]))
+                # operations.append(('av', vertex, self.vertices[vertex]))
                 self.vertices[vertex] = None
 
         # Register remaining vertices and faces
@@ -401,6 +457,14 @@ class MAPS(obja.Model):
                 output_model.edit_vertex(index, value)
             elif op == 'ef':
                 output_model.edit_face(index, value)
+            elif op == 'fc':
+                print('fc {} {} {} {}'.format(
+                    index,
+                    value[0],
+                    value[1],
+                    value[2]),
+                    file=output
+                )
 
 
 def main(args):
@@ -413,7 +477,7 @@ def main(args):
     model.parse_file(args.input)
 
     with open(args.output, 'w') as output:
-        model.contract(output)
+        model.truc(output)
 
 
 if __name__ == '__main__':