Modal Moment fitting and Heat Clustering

examples/notebooks/example_heat_clustering.ipynb

@@ -0,0 +1,157 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import simkit as sk\n",
+    "import igl\n",
+    "from simkit.filesystem import get_data_directory"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# -- 2D Meshes -- #\n",
+    "filename = '2d/beam/beam.obj'\n",
+    "filename = '2d/bingby/bingby.obj'\n",
+    "filename = '2d/cthulu/cthulu.obj'\n",
+    "\n",
+    "[V, _, _, T, _, _]= igl.readOBJ(get_data_directory() + filename)\n",
+    "V = V[:, :2]  # Keep only the first two dimensions for 2D meshes\n",
+    "print('Vertices:', V.shape)\n",
+    "print('Triangles:', T.shape)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "k= 5\n",
+    "mode = 'elem'\n",
+    "mu = np.ones((T.shape[0],1))\n",
+    "data = sk.heat_clustering_precomp(V, T, mu=mu, k=k, mode=mode)\n",
+    "parts, C, dist, clusters, centroids = sk.heat_clustering_solve(data, max_iters=100)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import polyscope as ps\n",
+    "import polyscope.imgui as psim\n",
+    "\n",
+    "if T.shape[1] == 3:\n",
+    "    elem_type = \"faces\"\n",
+    "else:\n",
+    "    elem_type = \"cells\"\n",
+    "if mode == 'elem':\n",
+    "    part_type = elem_type\n",
+    "else:\n",
+    "    part_type = \"vertices\"\n",
+    "\n",
+    "def visualize_heat_clustering(V, T, C, clusters, centroids, dist, parts, mu=None):\n",
+    "    ps.init()\n",
+    "\n",
+    "    # Register mesh\n",
+    "    if T.shape[1] == 3:\n",
+    "        msh = ps.register_surface_mesh(\"mesh\", V, T)\n",
+    "    else:\n",
+    "        msh = ps.register_volume_mesh(\"mesh\", V, T)\n",
+    "\n",
+    "    # Register centroids and initial cluster labeling\n",
+    "    pnts = ps.register_point_cloud(\"centroids\", C)\n",
+    "    msh.add_scalar_quantity(\"clusters\", parts, enabled=True, cmap=\"rainbow\", defined_on=part_type)\n",
+    "\n",
+    "    # Optional stiffness visualization\n",
+    "    if mu is not None:\n",
+    "        msh.add_scalar_quantity(\"mu\", mu.flatten(), defined_on=elem_type, cmap=\"rainbow\")\n",
+    "\n",
+    "    # Set up interactive sliders for cluster index and iteration\n",
+    "    current_k = 0\n",
+    "    current_iter = 0\n",
+    "\n",
+    "    def callback_func():\n",
+    "        nonlocal current_k, current_iter\n",
+    "\n",
+    "        # Distance field for cluster k\n",
+    "        changed, k = psim.SliderInt(\"k\", current_k, v_min=0, v_max=dist.shape[1]-1)\n",
+    "        if changed:\n",
+    "            current_k = k\n",
+    "            msh.add_scalar_quantity(\"d\", dist[:, k], defined_on=part_type, enabled=True, cmap=\"rainbow\")\n",
+    "\n",
+    "        # Cluster labeling at iteration\n",
+    "        changed, iter = psim.SliderInt(\"iter\", current_iter, v_min=0, v_max=len(clusters)-1)\n",
+    "        if changed:\n",
+    "            current_iter = iter\n",
+    "            pnts = ps.register_point_cloud(\"centroids\", centroids[iter])\n",
+    "            pnts.update_point_positions(centroids[iter])\n",
+    "            msh.add_scalar_quantity(\"clusters\", clusters[iter], enabled=True, cmap=\"rainbow\", defined_on=part_type)\n",
+    "\n",
+    "    ps.set_user_callback(callback_func)\n",
+    "    ps.show()\n",
+    "    ps.remove_all_structures()\n",
+    "\n",
+    "visualize_heat_clustering(V, T, C, clusters, centroids, dist, parts, mu=mu)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Next version try modifying the stiffness so that there that the top part is stiffer.\n",
+    "BC = igl.barycenter(V, T)\n",
+    "\n",
+    "min_y = np.min(BC[:, 1])\n",
+    "max_y = np.max(BC[:, 1])\n",
+    "min_x = np.min(BC[:, 0])\n",
+    "max_x = np.max(BC[:, 0])\n",
+    "y_range = max_y - min_y\n",
+    "x_range = max_x - min_x\n",
+    "\n",
+    "stiffness = np.ones((T.shape[0], 1))\n",
+    "top_mask = (BC[:, 1] > 0.5 * y_range + min_y)\n",
+    "middle_mask = (BC[:, 0] > 0.4 * x_range + min_x) & (BC[:, 0] < 0.6 * x_range + min_x)\n",
+    "mask = top_mask & middle_mask\n",
+    "stiffness[mask] = 100.0\n",
+    "k = 8\n",
+    "data = sk.heat_clustering_precomp(V, T, mu=stiffness, k=k, mode=mode)\n",
+    "parts, C, dist, clusters, centroids = sk.heat_clustering_solve(data)\n",
+    "visualize_heat_clustering(V, T, C, clusters, centroids, dist, parts, mu=stiffness)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "simkit",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.16"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

examples/run_elastics_reduced.py

@@ -2,6 +2,7 @@
 import igl
 import polyscope as ps
 import scipy as sp
+# import gpytoolbox as gpy
 
 from simkit import (
     deformation_jacobian,
@@ -20,28 +21,38 @@
     view_cubature,
     view_displacement_modes,
     view_scalar_modes)
+
 from simkit.filesystem import get_data_directory
 from simkit.sims.elastic import ElasticROMMFEMSim, ElasticROMMFEMSimParams, ElasticROMFEMSim, ElasticROMFEMSimParams
+from simkit.spectral_moment_fitting_cubature import spectral_moment_fitting_cubature
 
 [X, _, _, T, _, _] = igl.readOBJ(get_data_directory() + "2d/cthulu/cthulu.obj")
+# [X, T] = gpy.regular_square_mesh(30, 30)
 X = X[:, 0:2]
 
 X = X / max(X.max(axis=0) - X.min(axis=0))
 
 dim = X.shape[1]
 
 m = 10
-k = 100
 [W, E,  B] = skinning_eigenmodes(X, T, m)
+
+k = 100
 [cI, cW, labels] = spectral_cubature(X, T, W, k, return_labels=True)
 
+# p = 1 will give you m * k total cubature points (roughly)
+# p = 2 will give you about m choose 2 per cluster.
+# For m=k=10, you end up with up to 650 samples.
+k = 10
+p = 2
+[cI, cW, labels] = spectral_moment_fitting_cubature(X, T, W, k, return_labels=True, p=p)
 
 # B =  np.identity(X.shape[0]*dim)
 # cI = np.arange(T.shape[0])
 # cW = volume(X, T) 
 # view_scalar_modes(X, T, W)
 # view_displacement_modes(X, T, B, a=1)
-# view_cubature(X, T, cI, cW, labels)
+view_cubature(X, T, cI, cW, labels)
 
 A = np.zeros(((W.shape[1], dim, dim+1)))
 A[:, :dim, :dim] = np.identity(dim)

simkit/__init__.py

@@ -12,18 +12,24 @@
 
 from .massmatrix import massmatrix
 from .dirichlet_penalty import dirichlet_penalty
+from .dirichlet_laplacian import dirichlet_laplacian
 
 from .grad import grad
 from .volume import volume
 
 from .project_into_subspace import project_into_subspace
 from .pairwise_displacement import pairwise_displacement
+from .pairwise_distance import pairwise_distance
 
 from .volume import volume
 from .deformation_jacobian import deformation_jacobian
 from .polar_svd import polar_svd
 from .rotation_gradient import rotation_gradient_F
 from .skinning_eigenmodes import skinning_eigenmodes
 from .spectral_cubature import spectral_cubature
+from .spectral_clustering import spectral_clustering
 from .gravity_force import gravity_force
-from .cluster_grouping_matrices import cluster_grouping_matrices
+from .cluster_grouping_matrices import cluster_grouping_matrices
+from .average_onto_simplex import average_onto_simplex
+from .heat_distance import heat_distance_solve, heat_distance_precompute
+from .heat_clustering import heat_clustering_precomp, heat_clustering_solve, HeatClusteringData