From 23064caca9f3e3146819e0e1ab38e089a6463fc1 Mon Sep 17 00:00:00 2001
From: Anthony Mercer <anton@praviel.com>
Date: Wed, 21 Oct 2026 04:31:03 -0700
Subject: [PATCH] Compute cotangent weights correctly for integer mesh vertices

---
 src/potpourri3d/mesh.py  |  8 ++++++++
 test/potpourri3d_test.py | 38 ++++++++++++++++++++++++++++++++++++++
 2 files changed, 46 insertions(+)

diff --git a/src/potpourri3d/mesh.py b/src/potpourri3d/mesh.py
index e03e9e3..003438b 100644
--- a/src/potpourri3d/mesh.py
+++ b/src/potpourri3d/mesh.py
@@ -184,6 +184,11 @@ def cotan_laplacian(V, F, denom_eps=0.):
     validate_mesh(V, F, force_triangular=True)
     nV = V.shape[0]
 
+    integer_coordinates = np.issubdtype(V.dtype, np.integer)
+    if integer_coordinates:
+        # Keep coordinate differences exact, including unsigned and wide integers.
+        V = V.astype(np.int64 if V.dtype.itemsize < 8 else object)
+
     mat_i = []
     mat_j = []
     mat_data = []
@@ -195,6 +200,9 @@ def cotan_laplacian(V, F, denom_eps=0.):
         inds_k = F[:,(i+2)%3]
         vec_ki = V[inds_i,:] - V[inds_k,:]
         vec_kj = V[inds_j,:] - V[inds_k,:]
+        if integer_coordinates:
+            vec_ki = vec_ki.astype(float)
+            vec_kj = vec_kj.astype(float)
         dots = np.sum(vec_ki * vec_kj, axis=1)
         cross_mags =  np.linalg.norm(np.cross(vec_ki, vec_kj), axis=1)
         cotans = 0.5 * dots / (cross_mags + denom_eps)
diff --git a/test/potpourri3d_test.py b/test/potpourri3d_test.py
index 33d9450..5e74c30 100644
--- a/test/potpourri3d_test.py
+++ b/test/potpourri3d_test.py
@@ -153,6 +153,44 @@ def test_mesh_vector_heat(self):
                 self.assertEqual(logmap.shape[0], V.shape[0])
                 self.assertEqual(logmap.shape[1], 2)
     
+    def test_integer_cotan_laplacian(self):
+        cases = [(np.int8, 10), (np.uint8, 10),
+                 (np.int16, 100), (np.uint16, 100),
+                 (np.int32, 100000), (np.uint32, 100000),
+                 (np.int64, 100000), (np.uint64, 100000),
+                 (np.float32, 100), (np.float64, 100)]
+        # Integrated linear basis gradients on a right triangle of sides 3 and 4.
+        expected = np.array([[25. / 24., -2. / 3., -3. / 8.],
+                             [-2. / 3., 2. / 3., 0.],
+                             [-3. / 8., 0., 3. / 8.]])
+        for dtype, scale in cases:
+            V = np.array([[0, 0, 0], [3 * scale, 0, 0],
+                          [0, 4 * scale, 0]], dtype=dtype)
+            original = V.copy()
+            for winding in [[0, 1, 2], [0, 2, 1]]:
+                F = np.array([winding])
+                for eps in [0., 1.]:
+                    with self.subTest(dtype=dtype, winding=winding, eps=eps):
+                        L = pp3d.cotan_laplacian(V, F, denom_eps=eps)
+                        factor = 12. * scale * scale / (12. * scale * scale + eps)
+                        np.testing.assert_allclose(L.toarray(), expected * factor,
+                                                   rtol=1e-6, atol=1e-8)
+                        np.testing.assert_array_equal(V, original)
+
+    def test_integer_cotan_laplacian_large_offset(self):
+        expected = np.array([[25. / 24., -2. / 3., -3. / 8.],
+                             [-2. / 3., 2. / 3., 0.],
+                             [-3. / 8., 0., 3. / 8.]])
+        for dtype in [np.int64, np.uint64]:
+            V = np.array([[0, 0, 0], [3, 0, 0], [0, 4, 0]], dtype=dtype)
+            V += 2**60
+            original = V.copy()
+            for winding in [[0, 1, 2], [0, 2, 1]]:
+                with self.subTest(dtype=dtype, winding=winding):
+                    L = pp3d.cotan_laplacian(V, np.array([winding]))
+                    np.testing.assert_allclose(L.toarray(), expected)
+                    np.testing.assert_array_equal(V, original)
+
     def test_mesh_cotan_laplace(self):
 
         V, F = pp3d.read_mesh(os.path.join(asset_path, "bunny_small.ply"))