### Install libigl Python Bindings

Source: https://github.com/libigl/libigl-python-bindings/blob/main/README.md

Install the libigl package using pip. This is the primary method for users to get started with the library.

```bash
python -m pip install libigl
```

--------------------------------

### Setup Python Virtual Environment and Install cibuildwheel

Source: https://github.com/libigl/libigl-python-bindings/blob/main/README.md

Set up a Python virtual environment and install the cibuildwheel tool for testing cibuildwheel locally.

```bash
/Library/Frameworks/Python.framework/Versions/3.11/bin/python3.11 -m venv venv-official-3.11
source venv-official-3.11/bin/activate
python -m pip install cibuildwheel
```

--------------------------------

### Install Meshplot for Visualization

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tut-chapter0.ipynb

Install the meshplot library from a specific GitHub archive for visualization purposes. This is a dependency for plotting examples.

```bash
pip install https://github.com/skoch9/meshplot/archive/0.4.0.tar.gz
```

--------------------------------

### Read Mesh for Quadratic Programming Example

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tutorials.ipynb

Reads a mesh file for a quadratic programming example. This snippet is a precursor to more complex QP problems.

```python
#TODO: Check why results differ, add interactivity

v, f, _ = igl.read_off(os.path.join(root_folder, "data", "cheburashka.off"))
```

--------------------------------

### Example Binding Skeleton

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/contributing.md

This is an example of a binding skeleton file. You can find existing skeletons in subfolders of `attic` or copy an existing binding if a new file is needed.

```cpp
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>

#include <igl/decimate.h>

namespace py = pybind11;

// This is a placeholder for the actual binding code.
// The real binding code would involve py::module_ and py::class_ definitions.
// For example, to bind a function like igl::decimate:
// 
// PYBIND11_MODULE(libigl_python_bindings, m) {
//     m.doc() = "pybind11 plugin for libigl";
//     m.def("decimate", &igl::decimate<Eigen::MatrixXd, Eigen::MatrixXf, Eigen::MatrixXi>,
//           py::arg("vertices"), py::arg("faces"), py::arg("target_faces"));
// }

```

--------------------------------

### Install libigl Python Bindings

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/index.md

Install the libigl Python bindings using pip. This command is used to add the library to your Python environment.

```bash
python -m pip install libigl
```

--------------------------------

### Install meshplot

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tutorials.ipynb

Install the meshplot library from a specific tarball URL using pip.

```python
python -m pip install https://github.com/skoch9/meshplot/archive/0.4.0.tar.gz
```

--------------------------------

### Shape Deformation Setup

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tutorials.ipynb

Loads a mesh and prepares it for shape deformation by swapping axes and defining selection constraints. This is a common setup for various deformation techniques.

```python
v, f = igl.read_triangle_mesh(os.path.join(root_folder, "data", "decimated-max.obj"))
v[:,[0, 2]] = v[:,[2, 0]] # Swap X and Z axes
u = v.copy()
s = igl.read_dmat(os.path.join(root_folder, "data", "decimated-max-selection.dmat"))
b = np.array([[t[0] for t in [(i, s[i]) for i in range(0, v.shape[0])] if t[1] >= 0]]).T
```

--------------------------------

### LSCM Parametrization Example

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tutorials.ipynb

Demonstrates LSCM parametrization by fixing two boundary points and then computing the UV coordinates. The resulting parametrization can be visualized in 2D or 3D.

```python
b = np.array([2, 1])
bnd = igl.boundary_loop(f)
b[0] = bnd[0]
b[1] = bnd[int(bnd.size / 2)]
bc = np.array([[0.0, 0.0], [1.0, 0.0]])
_, uv = igl.lscm(v, f, b, bc)
p = plot(v, f, uv=uv, shading={"wireframe": False, "flat": False}, return_plot=True)

@interact(mode=['3D','2D'])
def switch(mode):
    if mode == "3D":
        plot(v, f, uv=uv, shading={"wireframe": False, "flat": False}, plot=p)
    if mode == "2D":
        plot(uv, f, uv=uv, shading={"wireframe": True, "flat": False}, plot=p)
```

--------------------------------

### faces_first (Reindexing Example)

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/igl_docs.md

Example of reindexing faces using the IM output from faces_first.

```APIDOC
## faces_first (Reindexing Example)

### Description
Example of reindexing faces using the IM output from the `faces_first` function.

### Method
Not specified (likely a function call)

### Endpoint
Not applicable (Python function)

### Parameters
#### Path Parameters
None

#### Query Parameters
None

#### Request Body
None

### Request Example
```python
T = T.unaryExpr(bind1st(mem_fun( static_cast<VectorXi::Scalar& (VectorXi::*)(VectorXi::Index)>(&VectorXi::operator())), &IM)).eval();
```

### Response
#### Success Response (200)
- **T** (array) - Reindexed faces.

#### Response Example
```json
[[0, 1, 2]]
```
```

--------------------------------

### LSCM Parametrization Setup

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tutorials.ipynb

Initializes mesh loading for Least Squares Conformal Maps (LSCM) parametrization. This snippet demonstrates reading a mesh file, which is the first step before computing the LSCM energy matrices.

```python
v, f = igl.read_triangle_mesh(os.path.join(root_folder, "data", "camelhead.off"))
```

--------------------------------

### Editable Install of Bindings

Source: https://github.com/libigl/libigl-python-bindings/blob/main/README.md

Perform an editable (incremental) build of the libigl Python bindings. This command is used for development and modification of the bindings.

```bash
CMAKE_BUILD_PARALLEL_LEVEL=10  python -m pip install --no-build-isolation --config-settings=editable.rebuild=true -Cbuild-dir=build -ve.
```

--------------------------------

### Binding a Class Example

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/contributing.md

This is a placeholder for binding a class. Binding classes is discouraged; instead, consider hiding the class internally and exposing only the final computation. Look at `classes.cpp` for an example if necessary.

```cpp
// Example from classes.cpp (conceptual)
// PYBIND11_MODULE(libigl_python_bindings, m) {
//     py::class_<MyClass>(m, "MyClass")
//         .def(py::init<int>())
//         .def("compute", &MyClass::compute);
// }

```

--------------------------------

### Import Libraries for Skinning

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tut-chapter7.ipynb

Imports necessary libraries for mesh processing, numerical operations, and visualization. Ensure these are installed before running.

```python
import igl
import scipy as sp
import numpy as np
import meshplot as mp

import os

root_folder = os.getcwd()
```

--------------------------------

### Generate and Visualize Multiple Iso-Surfaces

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tut-chapter6.ipynb

This example demonstrates generating multiple iso-surfaces by iterating through a range of isovalues and then visualizing selected surfaces. It loads tet mesh data, computes scalar fields, extracts meshes for different isovalues, and plots them.

```python
tv = np.load(os.path.join(root_folder, "data", "marching_cube_tv.npy"))
tt = np.load(os.path.join(root_folder, "data", "marching_cube_tt.npy"))
s = np.linalg.norm(tv, axis=1)

svs = []
sfs = []
for i in np.linspace(0.05, 0.75, 15):
    sv, sf, _, _ = igl.marching_tets(tv, tt, s, i)
    svs.append(sv)
    sfs.append(sf)


i = 0
for t in [3, 8, 11]:
    if i == 0:
        p = subplot(svs[t], sfs[t], s = [1, 3, i])
    else:
        subplot(svs[t], sfs[t], s = [1, 3, i], data=p)
    i += 1

p
```

--------------------------------

### Common Assertions Example

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/contributing.md

This snippet demonstrates common input assertions for size and type checking, as found in `common.h`. Ensure your bindings include similar checks.

```cpp
#include <Eigen/Core>

// Example assertion for matrix dimensions
// assert_valid_dimensions(matrix, expected_rows, expected_cols);
// Example assertion for data type
// assert_double_or_float(matrix);
// assert_int(matrix);

```

--------------------------------

### Active Set Solver Example

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tutorials.ipynb

Solves a constrained optimization problem using the active set method. Requires defining matrices for equality and inequality constraints, as well as lower and upper bounds.

```python
b = np.array([[2556]])
bc = np.array([[1.0]])
l = igl.cotmatrix(v, f)
m = igl.massmatrix(v, f, igl.MASSMATRIX_TYPE_VORONOI)
minv = sp.sparse.diags(1 / m.diagonal())
q = l @ (minv @ l)
bz = np.zeros((v.shape[0], 1))
lx = np.zeros((v.shape[0], 1))
ux = np.ones((v.shape[0], 1))
beq = np.array([[0.08]])
aeq = sp.sparse.csc_matrix(m.diagonal())
aieq = sp.sparse.csc_matrix((0, 0))
bieq = np.array([])
z = igl.active_set(q, bz, b, bc, aeq, beq, aieq, bieq, lx, ux, max_iter=8)
plot(v, f, z[1])
```

--------------------------------

### project_to_line

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/igl_docs.md

Projects points onto a line defined by a start point and a direction vector, and computes the squared distance to the line.

```APIDOC
## project_to_line

### Description
PROJECT_TO_LINE project points onto vectors, that is find the parameter t for a point p such that proj_p = (y-x).*t, additionally compute the squared distance from p to the line of the vector, such that p - proj_p² = sqr_d

### Method
`project_to_line(p: array, s: array, d: array)`

### Parameters
#### Path Parameters
- None

#### Query Parameters
- None

#### Request Body
- **P** (array) - #P by dim list of points to be projected
- **S** (array) - size dim start position of line vector
- **D** (array) - size dim destination position of line vector

### Request Example
```python
[T,sqrD] = project_to_line(P,S,D)
```

### Response
#### Success Response (200)
- **T** (array) - #P by 1 list of parameters
- **sqrD** (array) - #P by 1 list of squared distances

#### Response Example
```python
# Example response structure (actual values depend on input)
[T,sqrD] = ...
```
```

--------------------------------

### ARAP Parametrization Initialization

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tut-chapter4.ipynb

This snippet demonstrates initializing ARAP parametrization using harmonic parametrization as a starting point. It maps the boundary to a circle and then computes the harmonic mapping for internal vertices before solving with ARAP. Requires numpy and os imports.

```python
v, f  = igl.read_triangle_mesh(os.path.join(root_folder, "data", "camelhead.off"))

## Find the open boundary
bnd = igl.boundary_loop(f)

## Map the boundary to a circle, preserving edge proportions
bnd_uv = igl.map_vertices_to_circle(v, bnd)

## Harmonic parametrization for the internal vertices
uv = igl.harmonic(v, f, bnd, bnd_uv, 1)

arap = igl.ARAP(v, f, 2, np.zeros(0))
uva = arap.solve(np.zeros((0, 0)), uv)

p = subplot(v, f, uv=uva, shading={"wireframe": False, "flat": False}, s=[1, 2, 0])
p = subplot(uva, f, uv=uva, shading={"wireframe": False, "flat": False}, s=[1, 2, 1], data=p)

# @interact(mode=['3D','2D'])
# def switch(mode):
#     if mode == "3D":
#         plot(v, f, uv=uva, shading={"wireframe": False, "flat": False}, plot=p)
#     if mode == "2D":
#         plot(uva, f, uv=uva, shading={"wireframe": True, "flat": False}, plot=p)
```

--------------------------------

### Calculate Geodesic Distances and Visualize

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tutorials.ipynb

This example calculates geodesic distances from a single source vertex to all other vertices in the mesh. It then visualizes these distances by coloring the mesh based on the sine of the normalized distances. Ensure `numpy` is imported.

```python
v, f = igl.read_triangle_mesh(os.path.join(root_folder, "data", "armadillo.obj"))

## Select a vertex from which the distances should be calculated
vs = np.array([0])
##All vertices are the targets
vt = np.arange(v.shape[0])

d = igl.exact_geodesic(v, f, vs, vt)#, fs, ft)

strip_size = 0.1
##The function should be 1 on each integer coordinate
c = np.abs(np.sin((d / strip_size * np.pi)))
plot(v, f, c, shading={"wireframe": False})
```

--------------------------------

### Harmonic Parametrization Example

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tutorials.ipynb

Computes harmonic parametrization for a mesh. It first finds the boundary loop, maps it to a circle, and then computes the harmonic functions for the internal vertices. The result can be used to apply procedural textures.

```python
v, f  = igl.read_triangle_mesh(os.path.join(root_folder, "data", "camelhead.off"))
## Find the open boundary
bnd = igl.boundary_loop(f)

## Map the boundary to a circle, preserving edge proportions
bnd_uv = igl.map_vertices_to_circle(v, bnd)

## Harmonic parametrization for the internal vertices
uv = igl.harmonic(v, f, bnd, bnd_uv, 1)
v_p = np.hstack([uv, np.zeros((uv.shape[0],1))])

p = plot(v, f, uv=uv, shading={"wireframe": False, "flat": False}, return_plot=True)

@interact(mode=['3D','2D'])
def switch(mode):
    if mode == "3D":
        plot(v, f, uv=uv, shading={"wireframe": False, "flat": False}, plot=p)
    if mode == "2D":
        plot(v_p, f, uv=uv, shading={"wireframe": True, "flat": False}, plot=p)
```

--------------------------------

### Polyharmonic Deformation Example

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tutorials.ipynb

Deforms a flat domain into a bump by solving various k-harmonic PDEs. Adjusts the maximum height and the harmonic order (k).

```python
v, f = igl.read_triangle_mesh(os.path.join(root_folder, "data", "bump-domain.obj"))
u = v.copy()

# Find boundary vertices outside annulus
vrn = np.linalg.norm(v, axis = 1)
is_outer = [vrn[i] - 1.00 > -1e-15 for i in range(v.shape[0])]
is_inner = [vrn[i] - 0.15 < 1e-15 for i in range(v.shape[0])]
in_b = [is_outer[i] or is_inner[i] for i in range(len(is_outer))]

b = np.array([i for i in range(v.shape[0]) if (in_b[i])]).T
bc = np.zeros(b.size)

for bi in range(b.size):
    bc[bi] = 0.0 if is_outer[b[bi]] else 1.0

c = np.array(is_outer)

p = plot(u, f, c, shading={"wire_width": 0.01, "colormap": "tab10"}, return_plot=True)

@interact(z_max=(0.0, 1.0), k=(1, 4))
def update(z_max, k):
    z = igl.harmonic(v, f, b, bc, int(k))
    u[:, 2] = z_max * z
    p.update_object(vertices=u)
```

--------------------------------

### Define Mesh Vertices and Faces

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tut-chapter0.ipynb

Example of defining mesh geometry using numpy arrays for vertices (V) and triangle connectivity (F).

```python
V = np.array([
    [0., 0, 0],
    [1, 0, 0],
    [1, 1, 1],
    [2, 1, 0]
])

F = np.array([
    [0, 1, 2],
    [1, 3, 2]
])
```

--------------------------------

### Project Points onto Line

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/igl_docs.md

Projects a set of points onto a line defined by a start point and a direction vector. It also computes the squared distance from each point to the line.

```python
# P  #P by dim list of points to be projected
# S  size dim start position of line vector
# D  size dim destination position of line vector
[T,sqrD] = project_to_line(P,S,D)
```

--------------------------------

### Get libigl Package Version

Source: https://github.com/libigl/libigl-python-bindings/blob/main/README.md

Retrieve the installed version of the libigl package within your Python code. This method is available since version 2.5.4.dev0.

```python
import importlib.metadata
libigl_version = importlib.metadata.version('libigl')
```

--------------------------------

### Import libigl and dependencies

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tutorials.ipynb

Import necessary libraries including igl, scipy, numpy, meshplot, and git for data downloading. It also includes logic to clone a data repository if it doesn't exist.

```python
import igl
import scipy as sp
import numpy as np
from meshplot import plot, subplot, interact

import os
root_folder = os.getcwd()
#root_folder = os.path.join(os.getcwd(), "tutorial")
data_folder = os.path.join(root_folder,"/data")
# pip install gitpython
from git import Repo
if not os.path.isdir(data_folder):
    Repo.clone_from("https://github.com/libigl/libigl-tutorial-data.git", data_folder)
```

--------------------------------

### Fetch and Configure nanobind Dependency

Source: https://github.com/libigl/libigl-python-bindings/blob/main/CMakeLists.txt

Downloads and integrates the nanobind library using FetchContent. Ensure the specified Git repository and tag are accessible.

```cmake
# Enable FetchContent to download dependencies at configure time
include(FetchContent)

# Download and set up nanobind
FetchContent_Declare(
  nanobind
  GIT_REPOSITORY https://github.com/wjakob/nanobind.git
  GIT_TAG        v2.7.0
)
FetchContent_MakeAvailable(nanobind)
```

--------------------------------

### Fetch and Configure libigl Dependency

Source: https://github.com/libigl/libigl-python-bindings/blob/main/CMakeLists.txt

Downloads and integrates the libigl library using FetchContent. The specific Git tag ensures a consistent version is used.

```cmake
FetchContent_Declare(
  libigl
  GIT_REPOSITORY https://github.com/libigl/libigl.git
  GIT_TAG 678e1fff76815e0c4c5d1f025ee2129181cc7d86
)
FetchContent_MakeAvailable(libigl)
```

--------------------------------

### project_to_line_segment

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/igl_docs.md

Projects points onto a line segment defined by start and end points, and computes the squared distance to the segment.

```APIDOC
## project_to_line_segment

### Description
PROJECT_TO_LINE_SEGMENT project points onto vectors, that is find the parameter t for a point p such that proj_p = (y-x).*t, additionally compute the squared distance from p to the line of the vector, such that p - proj_p² = sqr_d

### Method
`project_to_line_segment(p: array, s: array, d: array)`

### Parameters
#### Path Parameters
- None

#### Query Parameters
- None

#### Request Body
- **P** (array) - #P by dim list of points to be projected
- **S** (array) - size dim start position of line vector
- **D** (array) - size dim destination position of line vector

### Request Example
```python
[T,sqrD] = project_to_line_segment(P,S,D)
```

### Response
#### Success Response (200)
- **T** (array) - #P by 1 list of parameters
- **sqrD** (array) - #P by 1 list of squared distances

#### Response Example
```python
# Example response structure (actual values depend on input)
[T,sqrD] = ...
```
```

--------------------------------

### Comparison of Skinning Techniques Visualization

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tut-chapter7.ipynb

Visualizes and compares Rigid Skinning, Linear Blend Skinning (LBS), and Direct Delta Mush (DDM) skinning techniques side-by-side. An interactive slider allows animating through different frames to observe the differences in deformation.

```python
viewer_comp = mp.Viewer({})
# Rigid Blue
offset_rigid = np.array([-90.0, 0.0, 0.0])
rigid = viewer_comp.add_mesh(v + offset_rigid, f, np.array([0.0, 0.5, 0.0]))

# LBS Red
offset_lbs = np.array([0.0, 0.0, 0.0])
lbs = viewer_comp.add_mesh(v + offset_lbs, f, np.array([0.5, 0.0, 0.0]))

# DDM Green
offset_ddm = np.array([90.0, 0.0, 0.0])
ddm = viewer_comp.add_mesh(v + offset_ddm, f, np.array([0.0, 0.5, 0.5]))


@mp.interact(frame=(1, num_frames))
def update_frame(frame):
    frame = frame - 1
    pose = anim[:, frame].reshape(parents.shape[0] * 4, 3, order='F')
    rigid_deformed_mesh = rigid_deform_mesh(pose)
    lbs_deformed_mesh = lbs_deform_mesh(pose)
    ddm_deformed_mesh = ddm_deform_mesh(pose)
    
    viewer_comp.update_object(oid=rigid, vertices=rigid_deformed_mesh + offset_rigid)
    viewer_comp.update_object(oid=lbs, vertices=lbs_deformed_mesh + offset_lbs)
    viewer_comp.update_object(oid=ddm, vertices=ddm_deformed_mesh + offset_ddm)
    
viewer_comp._renderer
```

--------------------------------

### Plot Mesh with N-RoSy Field and Singularities

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tutorials.ipynb

Plots a mesh with an N-RoSy field and its singularities. Constrained faces are colored red. Requires interactive plotting setup.

```python
v, f = igl.read_triangle_mesh(os.path.join(root_folder, "data", "bumpy.off"))
#print(points_n)

# Constrained faces id
b = np.array([[0]])

# Highlight in red the constrained faces
c = np.ones((f.shape[0], 3))
for i in range(b.size):
    c[b[i]] = np.array([1, 0, 0])

p = plot(v, f, c)
pn_id = pp_id = e_id = None


# Constrained faces representative vector
bcv = np.array([[1.0, 1.0, 1.0]])
avg = igl.avg_edge_length(v, f)
bc = igl.barycenter(v, f)

# Plots the mesh with an N-RoSy field and its singularities on top
# The constrained faces (b) are colored in red.
@interact(n=(1, 10))
def plot_mesh_nrosy(n=1):
    global pn_id, pp_id, e_id, bcv, avg, b, bc

    r, s = igl.nrosy(v, f, b, bcv, np.array([[]], dtype=np.int64), np.array([[]]), np.array([[]]), n, 0.5)
    
    # Expand the representative vectors in the full vector set and plot them as lines
    y = representative_to_nrosy(v, f, r, n)
    be = np.zeros((bc.shape[0] * n, 3))
    for i in range(bc.shape[0]):
        for j in range(n):
            be[i * n + j] = bc[i]

    if e_id:
        p.remove_object(e_id)
    e_id = p.add_lines(be, be + y * (avg / 2))

    # Plot the singularities as colored dots (red for negative, blue for positive)
    points_n = []
    points_p = []
    for i in range(0, s.size):
        if s[i] < -0.001:
            points_n.append(v[i])
        elif s[i] > 0.001:
            points_p.append(v[i])
    
    if pp_id and pn_id:
        p.remove_object(pn_id)
        p.remove_object(pp_id)
    if len(points_n) > 0:
        pn_id = p.add_points(np.array(points_n), c="red", shading={"point_size": 2.0})
    if len(points_p) > 0:
        pp_id = p.add_points(np.array(points_p), c="blue", shading={"point_size": 2.0})
```

--------------------------------

### Run cibuildwheel for macOS

Source: https://github.com/libigl/libigl-python-bindings/blob/main/README.md

Execute cibuildwheel to build wheels for macOS, specifying the output directory and platform. This is part of testing cibuildwheel locally.

```bash
CIBW_BUILD="cp311-*" python -m cibuildwheel --output-dir wheelhouse --platform macos
```

--------------------------------

### Quadratic energy minimization with fixed boundary conditions

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tut-chapter2.ipynb

Minimizes a quadratic energy function subject to fixed boundary values and optional equality constraints. Uses libigl's `min_quad_with_fixed` function.

```python
v, f = igl.read_triangle_mesh(os.path.join(root_folder, "data", "cheburashka.off"))

## Two fixed points: Left hand, left foot should have values 1 and -1
b = np.array([4331, 5957])
bc = np.array([1., -1.])
B = np.zeros((v.shape[0], 1))

## Construct Laplacian and mass matrix
L = igl.cotmatrix(v, f)
M = igl.massmatrix(v, f, igl.MASSMATRIX_TYPE_VORONOI)
Minv = sp.sparse.diags(1 / M.diagonal())

## Bi-Laplacian
Q = L @ (Minv @ L)

## Solve with only equality constraints
Aeq = sp.sparse.csc_matrix((0, 0))
Beq = np.array([])
_, z1 = igl.min_quad_with_fixed(Q, B, b, bc, Aeq, Beq, True)

## Solve with equality and linear constraints
Aeq = sp.sparse.csc_matrix((1, v.shape[0]))
Aeq[0,6074] = 1
Aeq[0, 6523] = -1
Beq = np.array([0.])
_, z2 = igl.min_quad_with_fixed(Q, B, b, bc, Aeq, Beq, True)

## Normalize colors to same range
min_z = min(np.min(z1), np.min(z2))
max_z = max(np.max(z1), np.max(z2))
z = [(z1 - min_z) / (max_z - min_z), (z2 - min_z) / (max_z - min_z)]

## Plot the functions
p = subplot(v, f, z[0], shading={"wireframe":False}, s=[1, 2, 0])
subplot(v, f, z[1], shading={"wireframe":False}, s=[1, 2, 1], data=p)
```

--------------------------------

### Compute and plot normalized Gaussian curvature

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tut-chapter1.ipynb

Normalizes Gaussian curvature by the area of each vertex's Voronoi cell to get an integral average. This requires computing the mass matrix.

```python
m = igl.massmatrix(v, f, igl.MASSMATRIX_TYPE_VORONOI)
minv = sp.sparse.diags(1 / m.diagonal())
kn = minv.dot(k)
plot(v, f, kn)
```

--------------------------------

### project

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/igl_docs.md

Projects 3D object points to 2D screen space using model, projection, and viewport matrices.

```APIDOC
## project

### Description
Project 3D object points to 2D screen space.

### Method
`project(v: array, model: array, proj: array, viewport: array)`

### Parameters
#### Path Parameters
- None

#### Query Parameters
- None

#### Request Body
- **v** (array) - #V by 3 list of object points
- **model** (array) - model matrix
- **proj** (array) - projection matrix
- **viewport** (array) - viewport vector

### Request Example
```python
# Assuming v, model, proj, viewport are defined
P = project(v, model, proj, viewport)
```

### Response
#### Success Response (200)
- **P** (array) - #V by 3 list of screen space points

#### Response Example
```python
# Example response structure (actual values depend on input)
P = ...
```
```

--------------------------------

### Visualize Mesh Subdivision Methods

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tut-chapter6.ipynb

Compares the 'Coarse', 'Upsample', and 'Loop' subdivision methods by plotting them side-by-side. Requires initial mesh loading.

```python
ov, of = igl.read_triangle_mesh(os.path.join(root_folder, "data", "decimated-knight.off"))
uv, uf = igl.upsample(ov, of)
lv, lf = igl.loop(ov, of)

p = subplot(ov, of, shading={"wireframe": True}, s=[1, 3, 0])
subplot(uv, uf, shading={"wireframe": True}, s=[1, 3, 1], data=p)
subplot(lv, lf, shading={"wireframe": True}, s=[1, 3, 2], data=p)
p
```

--------------------------------

### Configure libigl Build Options

Source: https://github.com/libigl/libigl-python-bindings/blob/main/CMakeLists.txt

Sets various options to control which features of libigl are built. Enable or disable these based on project requirements.

```cmake
# Download and set up libigl
option(LIBIGL_COPYLEFT_CORE        "Build target igl_copyleft::core"       ON)
option(LIBIGL_COPYLEFT_CGAL        "Build target igl_copyleft::cgal"       ON)
option(LIBIGL_EMBREE               "Build target igl::embree"              ON)
option(LIBIGL_COPYLEFT_TETGEN      "Build target igl_copyleft::tetgen"     ON)
option(LIBIGL_RESTRICTED_TRIANGLE  "Build target igl_restricted::triangle" ON)
```

--------------------------------

### Project Points onto Line Segment

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/igl_docs.md

Projects a set of points onto a line segment defined by a start and end point. It computes the parameter t along the segment and the squared distance to the line containing the segment.

```python
# P  #P by dim list of points to be projected
# S  size dim start position of line vector
# D  size dim destination position of line vector
[T,sqrD] = project_to_line_segment(P,S,D)
```

--------------------------------

### Get Unique Edges from Mesh Faces

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/igl_docs.md

Constructs a list of unique edges from a mesh's face list. Supports triangle or tet faces. The output is a list of edges represented as pairs of vertex indices.

```python
V, F, _ = igl.read_off("test.off")
E = igl.edges(F)
```

--------------------------------

### As-Rigid-As-Possible (ARAP) Parametrization

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tutorials.ipynb

Implements ARAP parametrization by first finding the boundary loop, mapping it to a circle, computing harmonic parametrization for internal vertices, and then solving for ARAP UV coordinates. The result can be plotted in 2D or 3D.

```python
v, f  = igl.read_triangle_mesh(os.path.join(root_folder, "data", "camelhead.off"))
bnd = igl.boundary_loop(f)
bnd_uv = igl.map_vertices_to_circle(v, bnd)
uv = igl.harmonic(v, f, bnd, bnd_uv, 1)
arap = igl.ARAP(v, f, 2, np.zeros((0, 0)))
uva = arap.solve(np.zeros((0, 0)), uv)
p = plot(v, f, uv=uva, shading={"wireframe": False, "flat": False}, return_plot=True)

@interact(mode=['3D','2D'])
def switch(mode):
    if mode == "3D":
        plot(v, f, uv=uva, shading={"wireframe": False, "flat": False}, plot=p)
    if mode == "2D":
        plot(uva, f, uv=uva, shading={"wireframe": True, "flat": False}, plot=p)
```

--------------------------------

### Compute Gradient Operator and Apply to Surface Function

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tut-chapter1.ipynb

Loads a mesh and scalar data, computes the gradient operator G, applies it to the data to get triangle gradients, and visualizes the results by drawing lines representing the gradient vectors.

```python
v, f = igl.read_triangle_mesh(os.path.join(root_folder, "data", "cheburashka.off"))
u = igl.read_dmat(os.path.join(root_folder, "data", "cheburashka-scalar.dmat"))

g = igl.grad(v, f)
gu = g.dot(u).reshape(f.shape, order="F")

gu_mag = np.linalg.norm(gu, axis=1)
p = plot(v, f, u, shading={"wireframe":False}, return_plot=True)

max_size = igl.avg_edge_length(v, f) / np.mean(gu_mag)
bc = igl.barycenter(v, f)
bcn = bc + max_size * gu
p.add_lines(bc, bcn, shading={"line_color": "black"})
```

--------------------------------

### Precompute ARAP Solver

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tut-chapter3.ipynb

Initializes the As-Rigid-As-Possible (ARAP) solver with the mesh geometry and boundary constraints. This precomputation step optimizes subsequent solve operations.

```python
# Precomputation
arap = igl.ARAP(v, f, 3, b)
```

--------------------------------

### Mesh Data Smoothing with Laplacian and Hessian Energies

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tut-chapter6.ipynb

Applies Laplacian and Hessian energy smoothing to noisy mesh data. Requires libigl and SciPy. The example visualizes the original, noisy, and smoothed functions, highlighting differences in boundary behavior.

```python
v, f = igl.read_triangle_mesh(os.path.join(root_folder, "data", "beetle.off"))
e = igl.edges(f)

# Constructing an exact function to smooth
z_exact = v[0:, 2] + 0.5 * v[0:, 1] + v[0:, 1] * v[0:, 1] + v[0:, 2] * v[0:, 2] * v[0:, 2]
    
# Make the exact function noisy
s = 0.2 * (np.max(z_exact) - np.min(z_exact))
np.random.seed(5)
z_noisy = z_exact + s * np.random.rand(*z_exact.shape)

# Constructing the squared Laplacian and squared Hessian energy
l = igl.cotmatrix(v, f)
m = igl.massmatrix(v, f, igl.MASSMATRIX_TYPE_BARYCENTRIC)

m_inv_l = sp.sparse.linalg.spsolve(m, l)
ql = l.T @ m_inv_l
qh = igl.hessian_energy(v, f)

# Solve to find Laplacian-smoothed and Hessian-smoothed solutions
al = 8e-4;
zl = sp.sparse.linalg.spsolve(al * ql + (1 - al) * m, al * m.dot(z_noisy))
ah = 5e-6;
zh = sp.sparse.linalg.spsolve(ah * qh + (1 - ah) * m, ah * m.dot(z_noisy))

# Calculate isolines
ilx_v, ilx_e = igl.isolines(v, f, z_exact, 30)
iln_v, iln_e = igl.isolines(v, f, z_noisy, 30)
ill_v, ill_e = igl.isolines(v, f, zl, 30)
ilh_v, ilh_e = igl.isolines(v, f, zh, 30)


#TODO add edges to subplot
p = subplot(v, f, z_exact, s=[2, 2, 0])
# p.view.add_edges(ilx_v, ilx_e)

subplot(v, f, z_noisy, s=[2, 2, 1], data=p)
# p.view.add_edges(iln_v, iln_e)

subplot(v, f, zl, s=[2, 2, 2], data=p)
# p.view.add_edges(ill_v, ill_e)

subplot(v, f, zh, s=[2, 2, 3], data=p)
# p.view.add_edges(ilh_v, ilh_e)
p

# e_id = p.add_edges(ilx_v, ilx_e)
# @interact(mode=['Original', 'Noisy', 'Biharmonic smoothing (0-Neumann)', 'Biharmonic smoothing (Natural Hessian)'])
# def switch(mode):
#     global e_id
#     p.remove_object(e_id)
#     if mode == "Original":
#         p.update_object(colors=z_exact)
#         e_id = p.add_edges(ilx_v, ilx_e)
#     if mode == "Noisy":
#         p.update_object(colors=z_noisy)
#         e_id = p.add_edges(iln_v, iln_e)
#     if mode == "Biharmonic smoothing (0-Neumann)":
#         p.update_object(colors=zl)
#         e_id = p.add_edges(ill_v, ill_e)
#     if mode == "Biharmonic smoothing (Natural Hessian)":
#         p.update_object(colors=zh)
#         e_id = p.add_edges(ilh_v, ilh_e)
```

--------------------------------

### Planarize Quad Mesh and Visualize Planarity

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tutorials.ipynb

Planarizes a quad mesh using a local/global approach and visualizes the planarity of quads with color. Requires interactive plotting.

```python
# Load a quad mesh generated by a conjugate field
vqc, fqc, _ = igl.read_off(os.path.join(root_folder, "data", "inspired_mesh_quads_Conjugate.off"))

# Convert it to a triangle mesh
fqc_tri = np.zeros((fqc.shape[0] * 2, 3), dtype="int64")
 הראשונה[:fqc.shape[0]] = fqc[:, :3]
fqc_tri[fqc.shape[0]:, 0] = fqc[:, 2]
fqc_tri[fqc.shape[0]:, 1] = fqc[:, 3]
fqc_tri[fqc.shape[0]:, 2] = fqc[:, 0]

# Planarize it
vqc_p = igl.planarize_quad_mesh(vqc, fqc, 100, 0.005)

# Calculate a color to each quad that corresponds to its planarity
planarity = igl.quad_planarity(vqc, fqc)
planarity_p = igl.quad_planarity(vqc_p, fqc)

c = np.concatenate([planarity, planarity])
c_p = np.concatenate([planarity_p, planarity_p])

p = plot(vqc, fqc_tri, c, shading={"normalize":[min(np.min(c), np.min(c_p)), max(np.max(c), np.max(c_p))]})

@interact(mode=['Curved','Planar'])
def switch(mode):
    if mode == "Curved":
        p.update_object(colors=c)
    if mode == "Planar":
        p.update_object(vertices=vqc_p, colors=c_p)
```

--------------------------------

### Solve QP with Active Set Solver

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tut-chapter2.ipynb

Use the active set solver to optimize discrete biharmonic kernels. Requires setting up the quadratic term (q), linear term (bz), equality constraints (b, bc, aeq, beq), inequality constraints (aieq, bieq), and bounds (lx, ux).

```python
#TODO: Check why results differ, add interactivity

v, f, _ = igl.read_off(os.path.join(root_folder, "data", "cheburashka.off"))

# One fixed point on belly
b = np.array([[2556]])
bc = np.array([[1.0]])

# Construct Laplacian and mass matrix
l = igl.cotmatrix(v, f)
m = igl.massmatrix(v, f, igl.MASSMATRIX_TYPE_VORONOI)
minv = sp.sparse.diags(1 / m.diagonal())

# Bi-Laplacian
q = l @ (minv @ l)

# Zero linear term
bz = np.zeros((v.shape[0], 1))

# Lower and upper bound
lx = np.zeros((v.shape[0], 1))
ux = np.ones((v.shape[0], 1))

# Equality constraint constrains solution to sum to 1
beq = np.array([[0.08]])
aeq = sp.sparse.csc_matrix(m.diagonal())

# Empty inequality constraints
aieq = sp.sparse.csc_matrix((0, 0))
bieq = np.array([])

z = igl.active_set(q, bz, b, bc, aeq, beq, aieq, bieq, lx, ux, max_iter=8)
plot(v, f, z[1])
```

--------------------------------

### Project 3D Points to Screen Space

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/igl_docs.md

Projects a list of 3D object points to 2D screen space using model, projection, and viewport matrices.

```python
# V  #V by 3 list of object points
# model  model matrix
# proj  projection matrix
# viewport  viewport vector
P = project(v, model, proj, viewport)
```

--------------------------------

### Solver Classes

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/igl_docs.md

APIs for solver classes like ARAP, BBW, SLIM, and shapeup.

```APIDOC
## class ARAP

### Description
ARAP solver class.

### Method
N/A (Python class method)

### Endpoint
N/A (Python class method)

### Parameters
#### Path Parameters
N/A

#### Query Parameters
N/A

#### Request Body
N/A

### Request Example
N/A

### Response
#### Success Response (200)
N/A

#### Response Example
N/A

## class BBW

### Description
BBW solver class.

### Method
N/A (Python class method)

### Endpoint
N/A (Python class method)

### Parameters
#### Path Parameters
N/A

#### Query Parameters
N/A

#### Request Body
N/A

### Request Example
N/A

### Response
#### Success Response (200)
N/A

#### Response Example
N/A

## class SLIM

### Description
SLIM solver class.

### Method
N/A (Python class method)

### Endpoint
N/A (Python class method)

### Parameters
#### Path Parameters
N/A

#### Query Parameters
N/A

#### Request Body
N/A

### Request Example
N/A

### Response
#### Success Response (200)
N/A

#### Response Example
N/A

## class shapeup

### Description
shapeup solver class.

### Method
N/A (Python class method)

### Endpoint
N/A (Python class method)

### Parameters
#### Path Parameters
N/A

#### Query Parameters
N/A

#### Request Body
N/A

### Request Example
N/A

### Response
#### Success Response (200)
N/A

#### Response Example
N/A
```

--------------------------------

### Initialize Plot and Perform ARAP Solves

Source: https://github.com/libigl/libigl-python-bindings/blob/main/tutorial/tut-chapter3.ipynb

Initializes a plot with the mesh and its pseudocolors, then iteratively solves for new vertex positions using the ARAP method with modified boundary conditions and updates the plot.

```python
# Plot the mesh with pseudocolors
p = subplot(v, f, c, s=[1, 4, 0])
for k in range(3):
    t= 1 + k*3
    bc = np.zeros((b.size, v.shape[1]))
    for i in range(0, b.size):
        bc[i] = v[b[i]]
        if s[b[i]] == 0:
            r = mid[0] * 0.25
            bc[i, 0] += r * np.sin(0.5 * t * 2 * np.pi)
            bc[i, 1] = bc[i, 1] - r + r * np.cos(np.pi + 0.5 * t * 2 * np.pi)
        elif s[b[i]] == 1:
            r = mid[1] * 0.15
            bc[i, 1] = bc[i, 1] + r + r * np.cos(np.pi + 0.15 * t * 2 * np.pi)
            bc[i, 2] -= r * np.sin(0.15 * t * 2 * np.pi)
        elif s[b[i]] == 2:
            r = mid[1] * 0.15
            bc[i, 2] = bc[i, 2] + r + r * np.cos(np.pi + 0.35 * t * 2 * np.pi)
            bc[i, 0] += r * np.sin(0.35 * t * 2 * np.pi)

    vn = arap.solve(bc, v)
    subplot(vn, f, c, s=[1, 4, k+1], data=p)
```