From 5f44c7f8f3eac4d14f274ea478f0e25b9b225b78 Mon Sep 17 00:00:00 2001
From: Alex Kaszynski <akascap@gmail.com>
Date: Tue, 7 May 2024 12:09:50 -0600
Subject: [PATCH] remove scipy and clean up create mesh

---
 pyacvd/clustering.py | 97 +++++++++++++++++++++++++++++---------------
 setup.py             |  2 +-
 2 files changed, 66 insertions(+), 33 deletions(-)

diff --git a/pyacvd/clustering.py b/pyacvd/clustering.py
index b5fa967..1e5750a 100644
--- a/pyacvd/clustering.py
+++ b/pyacvd/clustering.py
@@ -4,7 +4,8 @@
 
 import numpy as np
 import pyvista as pv
-from scipy import sparse
+from pyvista import ID_TYPE
+from vtkmodules.vtkCommonDataModel import vtkCellArray
 
 from pyacvd import _clustering
 
@@ -196,41 +197,77 @@ def cluster_centroid(cent, area, clusters):
     return cval.T
 
 
-def create_mesh(mesh, area, clusters, cnorm, flipnorm=True):
-    """Generates a new mesh given cluster data
+def polydata_from_faces(points: np.ndarray, faces: np.ndarray) -> pv.PolyData:
+    """
+    Generate a polydata from a faces array containing no padding and all triangles.
 
-    moveclus is a boolean flag to move cluster centers to the surface of their
-    corresponding cluster
+    Parameters
+    ----------
+    points : np.ndarray
+        Points array.
+    faces : np.ndarray
+        ``(n, 3)`` faces array.
+
+    Returns
+    -------
+    PolyData
+        New mesh.
 
     """
-    faces = mesh.faces.reshape(-1, 4)
-    points = mesh.points
-    if points.dtype != np.double:
-        points = points.astype(np.double)
+    if faces.ndim != 2:
+        raise ValueError("Expected a two dimensional face array.")
+
+    pdata = pv.PolyData()
+    pdata.points = points
+
+    carr = vtkCellArray()
+    offset = np.arange(0, faces.size + 1, faces.shape[1], dtype=ID_TYPE)
+    carr.SetData(pv.numpy_to_idarr(offset, deep=True), pv.numpy_to_idarr(faces, deep=True))
+    pdata.SetPolys(carr)
+    return pdata
+
+
+def create_mesh(
+    mesh: pv.PolyData,
+    area: np.ndarray,
+    clusters: np.ndarray,
+    cnorm: np.ndarray,
+    flipnorm: bool = True,
+) -> pv.PolyData:
+    """
+    Generate a new mesh given cluster data and other parameters.
 
-    # Compute centroids
-    ccent = np.ascontiguousarray(cluster_centroid(points, area, clusters))
+    Parameters
+    ----------
+    mesh : pyvista.PolyData
+        Input mesh object that needs to be reshaped.
+    area : numpy.ndarray
+        An array representing the area of each face of the mesh.
+    clusters : numpy.ndarray
+        An array representing clusters in the mesh.
+    cnorm : numpy.ndarray
+        An array representing the centroids of the clusters.
+    flipnorm : bool, default: True
+        If ``True``, flip the normals of the faces.
 
-    # Create sparse matrix storing the number of adjacent clusters a point has
-    rng = np.arange(faces.shape[0]).reshape((-1, 1))
-    a = np.hstack((rng, rng, rng)).ravel()
-    b = clusters[faces[:, 1:]].ravel()  # take?
-    c = np.ones(len(a), dtype="bool")
+    Returns
+    -------
+    pyvista.PolyData
+        Returns a PolyData object representing the new generated mesh.
 
-    boolmatrix = sparse.csr_matrix((c, (a, b)), dtype="bool")
+    """
+    faces = mesh._connectivity_array.reshape(-1, 3)
+    points = mesh.points.astype(np.float64, copy=False)
 
-    # Find all points with three neighboring clusters.  Each of the three
-    # cluster neighbors becomes a point on a triangle
-    nadjclus = boolmatrix.sum(1)
-    adj = np.array(nadjclus == 3).nonzero()[0]
-    idx = boolmatrix[adj].nonzero()[1]
+    # Compute centroids
+    ccent = np.ascontiguousarray(cluster_centroid(points, area, clusters))
 
-    # Append these points and faces
-    points = ccent
-    f = idx.reshape((-1, 3))
+    # Ignore faces that do not connect to different clusters
+    f_clus = np.sort(clusters[faces], axis=1)
+    mask = np.all(np.diff(f_clus, axis=1) != 0, axis=1)
+    f = f_clus[mask]
 
-    # Remove duplicate faces
-    f = f[unique_row_indices(np.sort(f, 1))]
+    f = f[unique_row_indices(f)]  # Remove duplicate faces
 
     # Mean normals of clusters each face is build from
     if flipnorm:
@@ -251,11 +288,7 @@ def create_mesh(mesh, area, clusters, cnorm, flipnorm=True):
         mask = agg < 0.0
         f[mask] = f[mask, ::-1]
 
-    # Create vtk surface
-    triangles = np.empty((f.shape[0], 4), dtype=f.dtype)
-    triangles[:, -3:] = f
-    triangles[:, 0] = 3
-    return pv.PolyData(points, triangles.ravel())
+    return polydata_from_faces(ccent, f)
 
 
 def unique_row_indices(a):
diff --git a/setup.py b/setup.py
index 4fb284b..a099c4d 100644
--- a/setup.py
+++ b/setup.py
@@ -55,6 +55,6 @@ def read(*paths):
         "Programming Language :: Python :: 3.12",
     ],
     python_requires=">=3.7",
-    install_requires=["pyvista>=0.37.0", "numpy", "scipy"],
+    install_requires=["pyvista>=0.37.0", "numpy"],
     keywords="vtk uniform meshing remeshing, acvd",
 )