ZCWang
/
ImplicitSurfaceNetwork
1
0
Fork 0
Code Issues 1 Pull Requests Projects Releases Wiki Activity
Browse Source

feat: enhance integrator initialization and improve v intersection finding logic 

- 1.  apply bbox of chain to speed up intersection finding
- 2.  Replace ray-based inside-test with connectivity-driven approach using
precomputed chain group map. Two chains are considered co-located in the
domain if they are mutually connected (bidirectional membership in
chain_group_map).
feat-integrator
mckay 2 weeks ago
parent
190dfc417f
commit
d3367a08b4
2 changed files with 102 additions and 68 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 20
      surface_integral/interface/subface_integrator.hpp
  2. 150
      surface_integral/src/subface_integrator.cpp

20
surface_integral/interface/subface_integrator.hpp
View File

@ -26,12 +26,7 @@ public:
* throughout the lifetime of this integrator.
*/
integrator_t(const stl_vector_mp<object_with_index_mapping<subface>>& surfaces,
const flat_hash_map_mp<uint32_t, parametric_plane_t>& uv_planes)
: m_subfaces(surfaces), m_uv_planes(uv_planes)
{
std::cout << "Integrator initialized with " << surfaces.size() << " subfaces." << std::endl;
std::cout << "m_uv_planes size: " << m_uv_planes.size() << std::endl;
}
const flat_hash_map_mp<uint32_t, parametric_plane_t>& uv_planes);
/// Default destructor
~integrator_t() = default;
@ -61,14 +56,19 @@ private:
/// Non-owning reference to the map of parametric planes (ID -> parametric_plane_t)
const flat_hash_map_mp<uint32_t, parametric_plane_t>& m_uv_planes;
// Precomputed bounding boxes for each chain in each parametric plane
flat_hash_map_mp<uint32_t, stl_vector_mp<Eigen::AlignedBox2d>> m_chain_bboxes_hash{};
stl_vector_mp<double> compute_u_breaks(const parametric_plane_t& param_plane, double u_min, double u_max) const;
stl_vector_mp<double> find_v_intersections_at_u(const subface& subface,
void find_v_intersections_at_u(const subface& subface,
const parametric_plane_t& param_plane,
double u_val,
bool if_add_boundary=true,
double v_min=0.0,
double v_max=0.0) const;
double v_min,
double v_max,
stl_vector_mp<double>& intersections,
stl_vector_mp<uint16_t>& intersected_chains
) const;
bool is_point_inside_domain(const subface& subface, const parametric_plane_t& param_plane, double u, double v) const;

150
surface_integral/src/subface_integrator.cpp
View File

@ -8,6 +8,37 @@
namespace internal
{
integrator_t(const stl_vector_mp<object_with_index_mapping<subface>>& surfaces,
const flat_hash_map_mp<uint32_t, parametric_plane_t>& uv_planes)
: m_subfaces(surfaces), m_uv_planes(uv_planes)
{
std::cout << "Integrator initialized with " << surfaces.size() << " subfaces." << std::endl;
std::cout << "m_uv_planes size: " << m_uv_planes.size() << std::endl;
// calculate chain bbox
for (int subface_index = 0; subface_index < m_subfaces.size(); ++subface_index) {
if (m_uv_planes.find(subface_index) != m_uv_planes.end()) {
const auto& param_plane = m_uv_planes.at(subface_index);
auto& chain_bboxs = m_chain_bboxes_hash[subface_index];
chain_bboxes.clear();
chain_bboxs.reverse(chain_vertices.size());
for (const auto& chain: param_plane.chain_vertices){
Eigen::AlignedBox2d bbox;
bbox.setEmpty();
if (!chain.empty()) {
for (const auto& pt : chain) {
bbox.extend(pt);
}
}
chain_bboxes.push_back(bbox);
}
}
}
m_chain_bboxes_hash
}
double integrator_t::calculate(
int gauss_order,
double (*func)(double u, double v, const Eigen::Vector3d& p, const Eigen::Vector3d& dU, const Eigen::Vector3d& dV)) const
@ -40,33 +71,32 @@ double integrator_t::calculate_one_subface(
// Gaussian integration in u direction
auto u_integrand = [&](double u) {
// Find exact v intersections for each u
stl_vector_mp<double> v_breaks = find_v_intersections_at_u(subface, param_plane, u, true,v_min, v_max);
stl_vector_mp<double> v_intersections;
stl_vector_mp<uint16_t> intersected_chains;
find_v_intersections_at_u(subface, param_plane, u, v_min, v_max, v_intersections, intersected_chains);
// Gaussian integration in v direction
auto v_integrand = [&](double v) {
// Check if point is inside domain
if (is_point_inside_domain(subface, param_plane, u, v)) {
try {
// Get evaluators for both directions
auto eval_du = subface.fetch_curve_constraint_evaluator(parameter_v_t{}, v); // Fix v, vary u → ∂/∂u
auto eval_dv = subface.fetch_curve_constraint_evaluator(parameter_u_t{}, u); // Fix u, vary v → ∂/∂v
auto res_u = eval_du(u); // f(u,v), grad_f = ∂r/∂u
auto res_v = eval_dv(v); // f(u,v), grad_f = ∂r/∂v
Eigen::Vector3d p = res_u.f.template head<3>(); // Point (x,y,z)
Eigen::Vector3d dU = res_u.grad_f.template head<3>(); // ∂r/∂u
Eigen::Vector3d dV = res_v.grad_f.template head<3>(); // ∂r/∂v
// Area element: ||dU × dV||
Eigen::Vector3d cross = dU.cross(dV);
double jacobian = cross.norm(); // Jacobian (area scaling factor)
return func(u, v, p, dU, dV) * jacobian;
} catch (...) {
return 0.0; // Skip singular points
}
try {
// Get evaluators for both directions
auto eval_du = subface.fetch_curve_constraint_evaluator(parameter_v_t{}, v); // Fix v, vary u → ∂/∂u
auto eval_dv = subface.fetch_curve_constraint_evaluator(parameter_u_t{}, u); // Fix u, vary v → ∂/∂v
auto res_u = eval_du(u); // f(u,v), grad_f = ∂r/∂u
auto res_v = eval_dv(v); // f(u,v), grad_f = ∂r/∂v
Eigen::Vector3d p = res_u.f.template head<3>(); // Point (x,y,z)
Eigen::Vector3d dU = res_u.grad_f.template head<3>(); // ∂r/∂u
Eigen::Vector3d dV = res_v.grad_f.template head<3>(); // ∂r/∂v
// Area element: ||dU × dV||
Eigen::Vector3d cross = dU.cross(dV);
double jacobian = cross.norm(); // Jacobian (area scaling factor)
return func(u, v, p, dU, dV) * jacobian;
} catch (...) {
return 0.0; // Skip singular points
}
return 0.0; // Point not in domain
};
double v_integral = 0.0;
@ -76,7 +106,7 @@ double integrator_t::calculate_one_subface(
// Check midpoint validity
double mid_v = (a + b) / 2.0;
if (is_point_inside_domain(subface, param_plane, u, mid_v)) {
if (is_edge_inside_domain(subface, param_plane, u, mid_v)) {
v_integral += integrate_1D(a, b, v_integrand, gauss_order);
} else {
std::cout << "uv out of domain: (" << u << "," << mid_v << ")" << std::endl;
@ -94,7 +124,9 @@ double integrator_t::calculate_one_subface(
double b = u_breaks[i + 1];
double mid_u = (a + b) / 2.0;
auto v_intersections = find_v_intersections_at_u(subface, param_plane, mid_u,true, v_min, v_max);
stl_vector_mp<double> v_intersections;
stl_vector_mp<uint16_t> intersected_chains;
find_v_intersections_at_u(subface, param_plane, mid_u, v_min, v_max, v_intersections, intersected_chains);
if (!v_intersections.empty()) {
integral += integrate_1D(a, b, u_integrand, gauss_order, is_u_near_singularity(mid_u));
@ -133,28 +165,37 @@ stl_vector_mp<double> integrator_t::compute_u_breaks(const parametric_plane_t& p
return stl_vector_mp<double>(break_set.begin(), break_set.end());
}
stl_vector_mp<double> integrator_t::find_v_intersections_at_u(const subface& subface,
void integrator_t::find_v_intersections_at_u(const subface& subface,
const parametric_plane_t& param_plane,
const stl_vector_mp<Eigen::AlignedBox2d>& chain_bboxs,
double u_val,
bool if_add_boundary,
double v_min,
double v_max
double v_max,
stl_vector_mp<double>& intersections,
stl_vector_mp<uint16_t>& intersected_chains
) const
{
stl_vector_mp<double> intersections;
if (if_add_boundary) {
intersections.push_back(v_min);
intersections.push_back(v_max);
}
intersections.clear();
intersected_chains.clear();
uint16_t chain_idx =0;
// Iterate over each boundary chain
for (const auto& chain : param_plane.chain_vertices) {
const size_t n_vertices = chain.size()-1;
if (n_vertices < 2){
std::cout << "Warning: Degenerate chain with less than 2 vertices." << std::endl;
chain_idx++;
continue; // Skip degenerate chains
}
// filter bbox
if (chain_bbox[chain_idx].isEmpty()
|| (chain_bbox[chain_idx].x()<u_val-1e-8 || chain_bbox[chain_idx].x()>u_val+1e-8)) {
chain_idx++;
continue;
}
// Iterate over each edge in the chain (including closing edge if desired)
for (size_t i = 0; i < n_vertices; ++i) {
const Eigen::Vector2d& v1 = chain[i]; // Current vertex: (u1, v1)
@ -182,6 +223,8 @@ stl_vector_mp<double> integrator_t::find_v_intersections_at_u(const subface&
if (pos1 == 0 && pos2 == 0) {
intersections.push_back(v1_val);
intersections.push_back(v2_val);
intersected_chains.push_back(chain_idx);
intersected_chains.push_back(chain_idx);
}
// Case 3: One endpoint on u = u_val or segment crosses u = u_val
else if (std::abs(u1 - u2) < eps) {
@ -189,6 +232,8 @@ stl_vector_mp<double> integrator_t::find_v_intersections_at_u(const subface&
if (pos1 == 0 || pos2 == 0) {
intersections.push_back(v1_val);
intersections.push_back(v2_val);
intersected_chains.push_back(chain_idx);
intersected_chains.push_back(chain_idx);
}
} else {
// General case: non-vertical segment crossing u = u_val
@ -211,48 +256,37 @@ stl_vector_mp<double> integrator_t::find_v_intersections_at_u(const subface&
try {
double v_solution = newton_method(target_function, v_initial);
intersections.push_back(v_solution);
intersected_chains.push_back(chain_idx);
} catch (...) {
// If Newton's method fails (e.g., divergence), fall back to linear approximation
intersections.push_back(v_initial);
intersected_chains.push_back(chain_idx);
}
}
}
chain_idx ++;
}
// Final step: sort and remove duplicates within tolerance
if (!intersections.empty()) { sort_and_unique_with_tol(intersections, 1e-8); }
//if (!intersections.empty()) { sort_and_unique_with_tol(intersections, 1e-8); }
return intersections;
}
/*
point (u, v) is inside the domain by ray-casting algorithm
To determine whether a point (u, v) is inside or outside a domain by counting the intersections of a vertical ray starting
from the point and extending upwards,
NOTE: when v_intersect - v < threshold, further checks are required to accurately determine if the intersection point lies
precisely on the boundary segment.
*/
bool integrator_t::is_point_inside_domain(const subface& subface,
const parametric_plane_t& param_plane,
double u,
double v) const
bool is_edge_inside_domain(const stl_vector_mp<stl_vector_mp<uint16_t>> & chain_group_indices,
uint16_t chain_idx1,
uint16_t chain_idx2) const
{
auto intersections = find_v_intersections_at_u(subface, param_plane, u, false);
const double tol_near = 1e-8;
const double tol_above = 1e-12;
uint32_t count = 0;
for (double v_int : intersections) {
double diff = v_int - v;
if (diff > tol_above) {
count++;
} else if (std::abs(diff) < tol_near) {
return true; // on boundary → inside
}
const auto& group1 = chain_group_indices[chain_idx1];
const auto& group2 = chain_group_indices[chain_idx2];
if (group1.find(chain_idx2) != group1.end() && group2.find(chain_idx1) != group2.end()) {
return true;
}
return (count % 2) == 1;
return false;
}
bool integrator_t::is_u_near_singularity(double u, double tol) const { return false; }
void integrator_t::sort_and_unique_with_tol(stl_vector_mp<double>& vec, double epsilon) const

Write Preview
Loading…
Cancel
Save