libigl_Martin/python/tutorial/702_WindingNumber.py

import sys, os

# Add the igl library to the modules search path
sys.path.insert(0, os.getcwd() + "/../")
import pyigl as igl

from shared import TUTORIAL_SHARED_PATH, check_dependencies, print_usage

dependencies = ["viewer"]
check_dependencies(dependencies)


def append_mesh(C_vis, F_vis, V_vis, V, F, color):
    F_vis.conservativeResize(F_vis.rows() + F.rows(), 3)
    F_vis.setBottomRows(F.rows(), F + V_vis.rows())
    V_vis.conservativeResize(V_vis.rows() + V.rows(), 3)
    V_vis.setBottomRows(V.rows(), V)
    C_vis.conservativeResize(C_vis.rows() + F.rows(), 3)
    colorM = igl.eigen.MatrixXd(F.rows(), C_vis.cols())
    colorM.rowwiseSet(color)
    C_vis.setBottomRows(F.rows(), colorM)


def update(viewer):
    global V, F, T, W, slice_z, overlay
    plane = igl.eigen.MatrixXd([0, 0, 1, -((1 - slice_z) * V.col(2).minCoeff() + slice_z * V.col(2).maxCoeff())])
    V_vis = igl.eigen.MatrixXd()
    F_vis = igl.eigen.MatrixXi()
    J = igl.eigen.MatrixXi()
    bary = igl.eigen.SparseMatrixd()
    igl.slice_tets(V, T, plane, V_vis, F_vis, J, bary)
    W_vis = igl.eigen.MatrixXd()
    igl.slice(W, J, W_vis)
    C_vis = igl.eigen.MatrixXd()
    igl.parula(W_vis, False, C_vis)

    if overlay == 1:  # OVERLAY_INPUT
        append_mesh(C_vis, F_vis, V_vis, V, F, igl.eigen.MatrixXd([[1., 0.894, 0.227]]))
    elif overlay == 2:  # OVERLAY_OUTPUT
        append_mesh(C_vis, F_vis, V_vis, V, F, igl.eigen.MatrixXd([[0.8, 0.8, 0.8]]))

    viewer.data.clear()
    viewer.data.set_mesh(V_vis, F_vis)
    viewer.data.set_colors(C_vis)
    viewer.data.set_face_based(True)


def key_down(viewer, key, modifier):
    global overlay, slice_z

    if key == ord(' '):
        overlay = (overlay + 1) % 3
    elif key == ord('.'):
        slice_z = min(slice_z + 0.01, 0.99)
    elif key == ord(','):
        slice_z = max(slice_z - 0.01, 0.01)

    update(viewer)

    return False


if __name__ == "__main__":
    keys = {"space": "toggle showing input mesh, output mesh or slice through tet-mesh of convex hull",
            ". / ,": "push back/pull forward slicing plane"}

    print_usage(keys)

    V = igl.eigen.MatrixXd()
    BC = igl.eigen.MatrixXd()
    W = igl.eigen.MatrixXd()
    T = igl.eigen.MatrixXi()
    F = igl.eigen.MatrixXi()
    G = igl.eigen.MatrixXi()

    slice_z = 0.5
    overlay = 0

    # Load mesh: (V,T) tet-mesh of convex hull, F contains facets of input
    # surface mesh _after_ self-intersection resolution
    igl.readMESH(TUTORIAL_SHARED_PATH + "big-sigcat.mesh", V, T, F)

    # Compute barycenters of all tets
    igl.barycenter(V, T, BC)

    # Compute generalized winding number at all barycenters
    print("Computing winding number over all %i tets..." % T.rows())
    igl.winding_number(V, F, BC, W)

    # Extract interior tets
    Wt = sum(W > 0.5)
    CT = igl.eigen.MatrixXi(Wt, 4)
    k = 0
    for t in range(T.rows()):
        if W[t] > 0.5:
            CT.setRow(k, T.row(t))
            k += 1

    # find bounary facets of interior tets
    igl.boundary_facets(CT, G)

    # boundary_facets seem to be reversed...
    G = G.rowwiseReverse()

    # normalize
    W = (W - W.minCoeff()) / (W.maxCoeff() - W.minCoeff())

    # Plot the generated mesh
    viewer = igl.viewer.Viewer()
    update(viewer)
    viewer.callback_key_down = key_down
    viewer.launch()