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