import igl
import numpy as np
from iglhelpers import *

V = igl.eigen.MatrixXd();
F = igl.eigen.MatrixXi();
igl.read_triangle_mesh("../tutorial/shared/fertility.off", V, F);

# Alternative discrete mean curvature
HN = igl.eigen.MatrixXd()
L = igl.eigen.SparseMatrixd()
M = igl.eigen.SparseMatrixd()
Minv = igl.eigen.SparseMatrixd()


igl.cotmatrix(V,F,L)
igl.massmatrix(V,F,igl.MASSMATRIX_TYPE_VORONOI,M)

igl.invert_diag(M,Minv)

# Laplace-Beltrami of position
HN = -Minv*(L*V)

# Extract magnitude as mean curvature
H = HN.rowwiseNorm()

# Compute curvature directions via quadric fitting
PD1 = igl.eigen.MatrixXd()
PD2 = igl.eigen.MatrixXd()

PV1 = igl.eigen.MatrixXd()
PV2 = igl.eigen.MatrixXd()

igl.principal_curvature(V,F,PD1,PD2,PV1,PV2)

# Mean curvature
H = 0.5*(PV1+PV2)

viewer = igl.viewer.Viewer()
viewer.data.set_mesh(V, F)

# Compute pseudocolor
C = igl.eigen.MatrixXd()
igl.parula(H,True,C)

viewer.data.set_colors(C)

# Average edge length for sizing
avg = igl.avg_edge_length(V,F)

# Draw a blue segment parallel to the minimal curvature direction
red = p2e(np.array([[0.8,0.2,0.2]]))
blue = p2e(np.array([[0.2,0.2,0.8]]))

viewer.data.add_edges(V + PD1*avg, V - PD1*avg, blue)

# Draw a red segment parallel to the maximal curvature direction
viewer.data.add_edges(V + PD2*avg, V - PD2*avg, red)

# Hide wireframe
viewer.core.show_lines = False

viewer.launch();