1#!/usr/bin/env python
2#
3# This file is part of libigl, a simple c++ geometry processing library.
4#
5# Copyright (C) 2017 Sebastian Koch <s.koch@tu-berlin.de> and Daniele Panozzo <daniele.panozzo@gmail.com>
6#
7# This Source Code Form is subject to the terms of the Mozilla Public License
8# v. 2.0. If a copy of the MPL was not distributed with this file, You can
9# obtain one at http://mozilla.org/MPL/2.0/.
10import sys, os
11
12# Add the igl library to the modules search path
13sys.path.insert(0, os.getcwd() + "/../")
14import pyigl as igl
15
16from shared import TUTORIAL_SHARED_PATH, check_dependencies, print_usage
17
18dependencies = ["glfw"]
19check_dependencies(dependencies)
20
21
22def append_mesh(C_vis, F_vis, V_vis, V, F, color):
23    F_vis.conservativeResize(F_vis.rows() + F.rows(), 3)
24    F_vis.setBottomRows(F.rows(), F + V_vis.rows())
25    V_vis.conservativeResize(V_vis.rows() + V.rows(), 3)
26    V_vis.setBottomRows(V.rows(), V)
27    C_vis.conservativeResize(C_vis.rows() + F.rows(), 3)
28    colorM = igl.eigen.MatrixXd(F.rows(), C_vis.cols())
29    colorM.rowwiseSet(color)
30    C_vis.setBottomRows(F.rows(), colorM)
31
32
33def update(viewer):
34    global V, F, T, W, slice_z, overlay
35    plane = igl.eigen.MatrixXd([0, 0, 1, -((1 - slice_z) * V.col(2).minCoeff() + slice_z * V.col(2).maxCoeff())])
36    V_vis = igl.eigen.MatrixXd()
37    F_vis = igl.eigen.MatrixXi()
38    J = igl.eigen.MatrixXi()
39    bary = igl.eigen.SparseMatrixd()
40    igl.marching_tets(V, T, plane, V_vis, F_vis, J, bary)
41    W_vis = igl.eigen.MatrixXd()
42    igl.slice(W, J, W_vis)
43    C_vis = igl.eigen.MatrixXd()
44    igl.parula(W_vis, False, C_vis)
45
46    if overlay == 1:  # OVERLAY_INPUT
47        append_mesh(C_vis, F_vis, V_vis, V, F, igl.eigen.MatrixXd([[1., 0.894, 0.227]]))
48    elif overlay == 2:  # OVERLAY_OUTPUT
49        append_mesh(C_vis, F_vis, V_vis, V, F, igl.eigen.MatrixXd([[0.8, 0.8, 0.8]]))
50
51    viewer.data().clear()
52    viewer.data().set_mesh(V_vis, F_vis)
53    viewer.data().set_colors(C_vis)
54    viewer.data().set_face_based(True)
55
56
57def key_down(viewer, key, modifier):
58    global overlay, slice_z
59
60    if key == ord(' '):
61        overlay = (overlay + 1) % 3
62    elif key == ord('.'):
63        slice_z = min(slice_z + 0.01, 0.99)
64    elif key == ord(','):
65        slice_z = max(slice_z - 0.01, 0.01)
66
67    update(viewer)
68
69    return False
70
71
72if __name__ == "__main__":
73    keys = {"space": "toggle showing input mesh, output mesh or slice through tet-mesh of convex hull",
74            ". / ,": "push back/pull forward slicing plane"}
75
76    print_usage(keys)
77
78    V = igl.eigen.MatrixXd()
79    BC = igl.eigen.MatrixXd()
80    W = igl.eigen.MatrixXd()
81    T = igl.eigen.MatrixXi()
82    F = igl.eigen.MatrixXi()
83    G = igl.eigen.MatrixXi()
84
85    slice_z = 0.5
86    overlay = 0
87
88    # Load mesh: (V,T) tet-mesh of convex hull, F contains facets of input
89    # surface mesh _after_ self-intersection resolution
90    igl.readMESH(TUTORIAL_SHARED_PATH + "big-sigcat.mesh", V, T, F)
91
92    # Compute barycenters of all tets
93    igl.barycenter(V, T, BC)
94
95    # Compute generalized winding number at all barycenters
96    print("Computing winding number over all %i tets..." % T.rows())
97    igl.winding_number(V, F, BC, W)
98
99    # Extract interior tets
100    Wt = sum(W > 0.5)
101    CT = igl.eigen.MatrixXi(Wt, 4)
102    k = 0
103    for t in range(T.rows()):
104        if W[t] > 0.5:
105            CT.setRow(k, T.row(t))
106            k += 1
107
108    # find bounary facets of interior tets
109    igl.boundary_facets(CT, G)
110
111    # boundary_facets seem to be reversed...
112    G = G.rowwiseReverse()
113
114    # normalize
115    W = (W - W.minCoeff()) / (W.maxCoeff() - W.minCoeff())
116
117    # Plot the generated mesh
118    viewer = igl.glfw.Viewer()
119    update(viewer)
120    viewer.callback_key_down = key_down
121    viewer.launch()
122