| 18 | using namespace PyMesh; |
| 19 | |
| 20 | void EdgeDihedralAngleAttribute::compute_from_mesh(Mesh& mesh) { |
| 21 | using EdgeMap = DupletMap<std::array<size_t, 2>>; |
| 22 | const size_t dim = mesh.get_dim(); |
| 23 | const size_t num_faces = mesh.get_num_faces(); |
| 24 | const size_t vertex_per_face = mesh.get_vertex_per_face(); |
| 25 | |
| 26 | if (dim != 3) { |
| 27 | throw NotImplementedError( |
| 28 | "Dihedral angle computation supports 3D mesh only."); |
| 29 | } |
| 30 | |
| 31 | if (!mesh.has_attribute("face_normal")) |
| 32 | mesh.add_attribute("face_normal"); |
| 33 | Matrix3Fr face_normals = |
| 34 | MatrixUtils::reshape<Matrix3Fr>(mesh.get_attribute("face_normal"), |
| 35 | num_faces, dim); |
| 36 | if (!face_normals.allFinite()) { |
| 37 | std::cerr << "Warning: Face normal contains nan or inf!" << std::endl; |
| 38 | std::cerr << " Dihedral angle computation would be unreliable." |
| 39 | << std::endl; |
| 40 | } |
| 41 | |
| 42 | EdgeMap edge_map; |
| 43 | for (size_t i=0; i<num_faces; i++) { |
| 44 | const auto f = mesh.get_face(i); |
| 45 | for (size_t j=0; j<vertex_per_face; j++) { |
| 46 | const Duplet key(f[j], f[(j+1)%vertex_per_face]); |
| 47 | edge_map.insert(key, {i, j}); |
| 48 | } |
| 49 | } |
| 50 | |
| 51 | auto angle = [&face_normals](size_t f1, size_t f2) { |
| 52 | const Eigen::Ref<Vector3F>& n1 = face_normals.row(f1).transpose(); |
| 53 | const Eigen::Ref<Vector3F>& n2 = face_normals.row(f2).transpose(); |
| 54 | return atan2(n1.cross(n2).norm(), n1.dot(n2)); |
| 55 | }; |
| 56 | |
| 57 | VectorF& dihedral_angles = m_values; |
| 58 | dihedral_angles = VectorF::Zero(num_faces * vertex_per_face); |
| 59 | |
| 60 | auto compute_kernel = [&angle, &dihedral_angles, vertex_per_face]( |
| 61 | EdgeMap::const_range_type& r) { |
| 62 | for (const auto& entry : r) { |
| 63 | const auto& key = entry.first; |
| 64 | const auto& val = entry.second; |
| 65 | const size_t num_adj_faces = val.size(); |
| 66 | assert(num_adj_faces > 0); |
| 67 | Float normal_angle = 0.0; |
| 68 | if (num_adj_faces == 1) { |
| 69 | normal_angle = 0.0; |
| 70 | } else if (num_adj_faces == 2) { |
| 71 | normal_angle = angle(val[0][0], val[1][0]); |
| 72 | } else { |
| 73 | for (size_t i=0; i<num_adj_faces; i++) { |
| 74 | for (size_t j=i+1; j<num_adj_faces; j++) { |
| 75 | normal_angle = std::max(normal_angle, |
| 76 | angle(val[i][0], val[j][0])); |
| 77 | } |
nothing calls this directly
no test coverage detected