21#include <lagrange/Mesh.h>
22#include <lagrange/common.h>
29template <
typename MeshType>
32 typename MeshType::Index seed_facet_id,
33 const Eigen::Matrix<typename MeshType::Scalar, 3, 1>& bc,
34 typename MeshType::Scalar radius = 0.0)
36 using Scalar =
typename MeshType::Scalar;
37 if (mesh.get_dim() != 3) {
38 throw std::runtime_error(
"Input mesh must be 3D mesh.");
40 if (mesh.get_vertex_per_facet() != 3) {
41 throw std::runtime_error(
"Input mesh is not triangle mesh");
43 if (!mesh.is_connectivity_initialized()) {
44 mesh.initialize_connectivity();
47 using Index =
typename MeshType::Index;
48 using FacetType = Eigen::Matrix<Index, 3, 1>;
49 using VertexType =
typename MeshType::VertexType;
50 using AttributeArray =
typename MeshType::AttributeArray;
53 radius = std::numeric_limits<Scalar>::max();
55 const Index num_facets = mesh.get_num_facets();
56 const Index num_vertices = mesh.get_num_vertices();
57 const auto& vertices = mesh.get_vertices();
58 const auto& facets = mesh.get_facets();
60 const FacetType seed_facet = facets.row(seed_facet_id);
61 const VertexType seed_point = vertices.row(seed_facet[0]) * bc[0] +
62 vertices.row(seed_facet[1]) * bc[1] +
63 vertices.row(seed_facet[2]) * bc[2];
65 using Entry = std::pair<Index, Scalar>;
66 auto comp = [](
const Entry& e1,
const Entry& e2) {
return e1.second > e2.second; };
67 std::priority_queue<Entry, std::vector<Entry>,
decltype(comp)> Q(comp);
68 AttributeArray dist(num_vertices, 1);
70 Q.push({seed_facet[0], (vertices.row(seed_facet[0]) - seed_point).norm()});
71 Q.push({seed_facet[1], (vertices.row(seed_facet[1]) - seed_point).norm()});
72 Q.push({seed_facet[2], (vertices.row(seed_facet[2]) - seed_point).norm()});
74 std::list<Entry> involved_vts;
76 Entry entry = Q.top();
77 involved_vts.push_back(entry);
79 Scalar curr_dist = dist(entry.first, 0);
81 if (curr_dist < 0 || entry.second < curr_dist) {
82 dist(entry.first, 0) = entry.second;
83 const auto& adj_vertices = mesh.get_vertices_adjacent_to_vertex(entry.first);
84 for (
const auto& vj : adj_vertices) {
85 Scalar d = entry.second + (vertices.row(vj) - vertices.row(entry.first)).norm();
87 Scalar next_dist = dist(vj, 0);
88 if (next_dist < 0 || d < next_dist) {
96 mesh.add_vertex_attribute(
"dijkstra_distance");
97 mesh.import_vertex_attribute(
"dijkstra_distance", dist);
std::optional< std::vector< Index > > compute_dijkstra_distance(SurfaceMesh< Scalar, Index > &mesh, const DijkstraDistanceOptions< Scalar, Index > &options={})
Computes dijkstra distance from a seed facet.
Definition: compute_dijkstra_distance.cpp:24
#define la_runtime_assert(...)
Runtime assertion check.
Definition: assert.h:169
Main namespace for Lagrange.
Definition: AABBIGL.h:30