integration-of-complex-surfaces
/
ImplicitSurfaceNetwork
mirror of http://118.195.195.192:3000/ZCWang/ImplicitSurfaceNetwork.git
2
0
Fork 1
Code Issues Projects Releases Wiki Activity
Browse Source

Refactor SurfaceAreaCalculator: streamline constructors, improve function naming, and enhance integration methods

wch
mckay 3 months ago
parent
ada2bdd02a
commit
d1b82fdcd7
2 changed files with 147 additions and 158 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 20
      surface_integral/interface/SurfaceIntegrator.hpp
  2. 245
      surface_integral/src/SurfaceIntegrator.cpp

20
surface_integral/interface/SurfaceIntegrator.hpp
View File

@ -1,7 +1,7 @@
#pragma once #pragma once
#include <primitive_process/interface/base/subface.hpp> #include <primitive_process/interface/base/subface.hpp>
#include "surface_integral/include/quadrature.hpp"
namespace internal namespace internal
{ {
@ -19,8 +19,7 @@ class SurfaceAreaCalculator
public: public:
// 构造函数，接受对 subface 的引用 // 构造函数，接受对 subface 的引用
explicit SurfaceAreaCalculator(const subface& surface); explicit SurfaceAreaCalculator(const subface& surface);
[[deprecated("Use calculate_new() instead")]] [[deprecated("Use calculate_new() instead")]] SurfaceAreaCalculator(const subface& surface,
SurfaceAreaCalculator(const subface& surface,
const stl_vector_map<double>& u_breaks, const stl_vector_map<double>& u_breaks,
double umin, double umin,
double umax, double umax,
@ -45,17 +44,14 @@ private:
double Vmin = 0.0; double Vmin = 0.0;
double Vmax = 1.0; double Vmax = 1.0;
// 私有辅助函数（可扩展） // 私有辅助函数
template <typename Func>
double SurfaceAreaCalculator::GaussIntegrate1D(double a, double b, Func&& func, int q) const;
// 直线u=u_val与边界的交点 // 直线u=u_val与边界的交点
std::vector<double> FindVerticalIntersectionsOCCT(double u_val, std::vector<double> find_vertical_intersections(double u_val,
const std::vector<Handle(Geom2d_Curve)>& edges, const std::vector<Handle(Geom2d_Curve)>& edges,
double v_min, double v_min,
double v_max); double v_max);
bool IsPointInsideDomain(double u, bool is_point_inside_domain(double u,
double v, double v,
const std::vector<Handle(Geom2d_Curve)>& outerEdges, const std::vector<Handle(Geom2d_Curve)>& outerEdges,
const std::vector<Handle(Geom2d_Curve)>& innerEdges, const std::vector<Handle(Geom2d_Curve)>& innerEdges,
@ -63,7 +59,11 @@ private:
double u_max, double u_max,
double v_min, double v_min,
double v_max); double v_max);
double integrate_over_uv(int num_samples) const; double newton_method(const std::function<equation_intermediate_t(double)>& F,
double v_initial,
double tolerance = 1e-8,
int max_iterations = 100);
void sort_and_unique_with_tol(std::vector<double>& vec, double epsilon = 1e-8);
}; };
} // namespace internal } // namespace internal

245
surface_integral/src/SurfaceIntegrator.cpp
View File

@ -20,9 +20,8 @@ SurfaceAreaCalculator::SurfaceAreaCalculator(const subface& surfa
} }
SurfaceAreaCalculator::SurfaceAreaCalculator(const subface& surface, const parametric_plane& uv_plane) SurfaceAreaCalculator::SurfaceAreaCalculator(const subface& surface, const parametric_plane& uv_plane)
: m_surface(surface), m_uv_plane(uv_plane), : m_surface(surface), m_uv_plane(uv_plane), Umin(0.0), Umax(0.0), Vmin(0.0), Vmax(0.0)
Umin(0.0), Umax(0.0), Vmin(0.0), Vmax(0.0) { {
if (!uv_plane.chain_vertices.empty()) { if (!uv_plane.chain_vertices.empty()) {
// 初始化为第一个点的坐标 // 初始化为第一个点的坐标
double min_u = uv_plane.chain_vertices[0].x(); double min_u = uv_plane.chain_vertices[0].x();
@ -83,28 +82,12 @@ double SurfaceAreaCalculator::calculate(Func&& func, int gauss_order) const
// 在u方向进行高斯积分 // 在u方向进行高斯积分
auto u_integrand = [&](double u) { auto u_integrand = [&](double u) {
// 对每个u，找到v方向的精确交点 // 对每个u，找到v方向的精确交点
std::vector<double> v_breaks = {self.Vmin, self.Vmax}; std::vector<double> v_breaks = find_vertical_intersections(u);
auto outer_intersections = FindVerticalIntersectionsOCCT(u, self.outerEdges, self.Vmin, self.Vmax);
v_breaks.insert(v_breaks.end(), outer_intersections.begin(), outer_intersections.end());
auto inner_intersections = FindVerticalIntersectionsOCCT(u, self.innerEdges, self.Vmin, self.Vmax);
v_breaks.insert(v_breaks.end(), inner_intersections.begin(), inner_intersections.end());
// 排序并去重
sort_and_unique_with_tol(v_breaks, 1e-6);
// 在v方向进行高斯积分 // 在v方向进行高斯积分
auto v_integrand = [&](double v) { auto v_integrand = [&](double v) {
// 判断点是否在有效域内 // 判断点是否在有效域内
if (IsPointInsideself(u, if (IsPointInsideself(u, v, self.outerEdges, self.innerEdges, self.Umin, self.Umax, self.Vmin, self.Vmax)) {
v,
self.outerEdges,
self.innerEdges,
self.Umin,
self.Umax,
self.Vmin,
self.Vmax)) {
try { try {
gp_Pnt p; gp_Pnt p;
gp_Vec dU, dV; gp_Vec dU, dV;
@ -126,15 +109,8 @@ double SurfaceAreaCalculator::calculate(Func&& func, int gauss_order) const
// 检查区间中点是否有效 // 检查区间中点是否有效
double mid_v = (a + b) / 2.0; double mid_v = (a + b) / 2.0;
if (IsPointInsideself(u, if (IsPointInsideself(u, mid_v, self.outerEdges, self.innerEdges, self.Umin, self.Umax, self.Vmin, self.Vmax)) {
mid_v, v_integral += gauss_integrate_1D(a, b, v_integrand, gauss_order);
self.outerEdges,
self.innerEdges,
self.Umin,
self.Umax,
self.Vmin,
self.Vmax)) {
v_integral += GaussIntegrate1D(a, b, v_integrand, gauss_order);
} else { } else {
std::cout << "uv out of domain: (" << u << "," << mid_v << ")" << std::endl; std::cout << "uv out of domain: (" << u << "," << mid_v << ")" << std::endl;
} }
@ -151,11 +127,11 @@ double SurfaceAreaCalculator::calculate(Func&& func, int gauss_order) const
// 检查区间中点是否有效 // 检查区间中点是否有效
double mid_u = (a + b) / 2.0; double mid_u = (a + b) / 2.0;
auto v_intersections = FindVerticalIntersectionsOCCT(mid_u, self.outerEdges, self.Vmin, self.Vmax); auto v_intersections = find_vertical_intersections(mid_u, self.outerEdges, self.Vmin, self.Vmax);
if (!v_intersections.empty()) { // 确保该u区间有有效区域 if (!v_intersections.empty()) { // 确保该u区间有有效区域
double integralp = GaussIntegrate1D(a, b, u_integrand, gauss_order); double integralp = gauss_integrate_1D(a, b, u_integrand, gauss_order);
integral += integralp; integral += integralp;
std::cout << "integral " << i << ": " << integralp << std::endl; std::cout << "integral " << i << ": " << integralp << std::endl;
} }
@ -164,131 +140,144 @@ double SurfaceAreaCalculator::calculate(Func&& func, int gauss_order) const
return integral; return integral;
} }
// 在区间[a,b]上进行高斯积分采样
template<typename Func>
double SurfaceAreaCalculator::GaussIntegrate1D(double a, double b, Func&& func, int q) const {
double sum = 0.0;
for (int i = 0; i < q; ++i) {
double x=a+(b-a)*GaussQuad::x(q, i);
double w=GaussQuad::w(q, i);
sum += w * func(x);
}
return sum * (b-a);
}
// 直线u=u_val与边界的交点 // 直线u=u_val与边界的交点
std::vector<double> SurfaceAreaCalculator::FindVerticalIntersectionsOCCT( std::vector<double> SurfaceAreaCalculator::find_vertical_intersections(double u_val)
double u_val,
const std::vector<Handle(Geom2d_Curve)>& edges,
double v_min,
double v_max)
{ {
std::vector<double> intersections; std::vector<double> intersections;
std::vector<int32_t> uPositionFlags;
// 创建垂直线(u=u_val) uPositionFlags.reserve(m_uv_plane.chain_vertices.size());
gp_Lin2d vertical_line(gp_Pnt2d(u_val, v_min), gp_Dir2d(0, 1));
Handle(Geom2d_Line) hVerticalLine = new Geom2d_Line(vertical_line); std::transform(m_uv_plane.chain_vertices.begin(),
m_uv_plane.chain_vertices.end(),
for (const auto& edge : edges) { std::back_inserter(uPositionFlags),
if (edge.IsNull()) { [&](const auto& currentVertex) {
continue; double uDifference = currentVertex.x() - u_val;
if (uDifference < 0) return -1; // 在参考值左侧
if (uDifference > 0) return 1; // 在参考值右侧
return 0; // 等于参考值
});
uint32_t group_idx = 0;
for (uint32_t element_idx = 0; element_idx < m_uv_plane.chains.index_group.size() - 1; element_idx++) {
if (element_idx > m_uv_plane.chains.start_indices[group_idx + 1]) group_idx++;
if (element_idx + 1 < m_uv_plane.chains.start_indices[group_idx + 1]) {
uint32_t vertex_idx1 = m_uv_plane.chains.index_group[element_idx];
uint32_t vertex_idx2 = m_uv_plane.chains.index_group[element_idx + 1];
if (uPositionFlags[vertex_idx1] * uPositionFlags[vertex_idx2] <= 0) {
auto v1 = m_uv_plane.chain_vertices[vertex_idx1], v2 = m_uv_plane.chain_vertices[vertex_idx2];
if (v1.x() != v2.x()) {
// "The line segment is vertical (u₁ == u₂), so there is no unique v value corresponding to the given u."
double v_initial = v1.y() + (v2.y() - v1.y()) * (u_val - v1.x()) / (v2.x() - v1.x());
auto curve_evaluator = m_surface.fetch_curve_constraint_evaluator(u);
auto solver_evaluator = m_surface.fetch_solver_evaluator();
auto target_function = [&](double v) -> equation_intermediate_t {
constraint_curve_intermediate temp_res = curve_evaluator(v);
return solver_evaluator(std::move(temp_res));
};
double v_solution = newton_method(target_function, v_initial);
intersections.push_back(v_solution);
} else {
intersections.push_back(v1.y());
intersections.push_back(v2.y());
} }
// 使用Geom2dAPI_InterCurveCurve求交
Geom2dAPI_InterCurveCurve intersector(edge, hVerticalLine);
if (intersector.NbPoints() > 0) {
for (int i = 1; i <= intersector.NbPoints(); ++i) {
gp_Pnt2d pt = intersector.Point(i);
double v = pt.Y();
if (v >= v_min && v <= v_max) {
intersections.push_back(v);
} }
} }
} }
// 补充采样法确保可靠性 // 去重排序
Geom2dAdaptor_Curve adaptor(edge); sort_and_unique_with_tol(intersections);
double first = adaptor.FirstParameter(); return intersections;
double last = adaptor.LastParameter();
const int samples = 20;
gp_Pnt2d prev_p;
adaptor.D0(first, prev_p);
double prev_v = prev_p.Y();
for (int i = 1; i <= samples; ++i) {
double param = first + i * (last - first) / samples;
gp_Pnt2d curr_p;
adaptor.D0(param, curr_p);
double curr_v = curr_p.Y();
if ((prev_v - u_val) * (curr_v - u_val) <= 0) {
// 牛顿迭代法
double t = param - (last - first)/samples;
for (int iter = 0; iter < 5; ++iter) {
gp_Pnt2d p;
gp_Vec2d deriv;
adaptor.D1(t, p, deriv);
if (std::abs(p.X() - u_val) < Precision::Confusion()) {
if (p.Y() >= v_min && p.Y() <= v_max) {
intersections.push_back(p.Y());
} }
break;
/*
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 is_point_inside_domain(double u, double v)
{
uint32_t group_idx = 0, intersection_count = 0;
for (uint32_t element_idx = 0; element_idx < m_uv_plane.chains.index_group.size() - 1; element_idx++) {
if (element_idx > m_uv_plane.chains.start_indices[group_idx + 1]) group_idx++;
if (element_idx + 1 < m_uv_plane.chains.start_indices[group_idx + 1]) {
uint32_t vertex_idx1 = m_uv_plane.chains.index_group[element_idx];
uint32_t vertex_idx2 = m_uv_plane.chains.index_group[element_idx + 1];
auto v1 = m_uv_plane.chain_vertices[vertex_idx1], v2 = m_uv_plane.chain_vertices[vertex_idx2];
if ((v1.x() <= u && v2.x() > u) || (v2.x() < u && v1.x() >= u)) {
double v_intersected = v1.y() + (v2.y() - v1.y()) * (u_val - v1.x()) / (v2.x() - v1.x());
if (v_interdected - v >= 1e-6) { intersection_count++; }
else if (std::abs(v_intersected - v) < 1e-6) {
auto curve_evaluator = m_surface.fetch_curve_constraint_evaluator(u);
auto solver_evaluator = m_surface.fetch_solver_evaluator();
auto target_function = [&](double v) -> equation_intermediate_t {
constraint_curve_intermediate temp_res = curve_evaluator(v);
return solver_evaluator(std::move(temp_res));
};
double v_solution = newton_method(target_function, v_initial);
if (std::abs(v_solution - v) > 0) {
intersection_count++;
} }
t -= (p.X() - u_val) / deriv.X();
} }
} }
prev_v = curr_v; /*
case v1.x() == v2.x() == u, do not count as intersection.but will cout in next iteration
*/
} }
} }
return intersection_count % 2 == 1; // in domain
// 去重排序
std::sort(intersections.begin(), intersections.end());
intersections.erase(std::unique(intersections.begin(), intersections.end()), intersections.end());
return intersections;
} }
// 牛顿法求解器
// Private method: Numerical integration over the UV parameter domain double SurfaceAreaCalculator::newton_method(const std::function<equation_intermediate_t(double)>& F,
double SurfaceAreaCalculator::integrate_over_uv(int num_samples) const double v_initial,
double tolerance,
int max_iterations)
{ {
double area = 0.0; double v = v_initial;
double du = 1.0 / num_samples; for (int i = 0; i < max_iterations; ++i) {
double dv = 1.0 / num_samples; equation_intermediate_t res = F(v);
double f = res.f;
double df = res.df;
// Fetch evaluators for point and derivative calculation if (std::abs(f) < tolerance) {
auto point_func = m_surface.fetch_point_by_param_evaluator(); std::cout << "✅ Converged in " << i + 1 << " iterations. v = " << v << std::endl;
auto solver_func = m_surface.fetch_solver_evaluator(); return v;
}
for (int i = 0; i < num_samples; ++i) { if (std::abs(df) < 1e-10) {
double u = (i + 0.5) * du; // Midpoint sampling in u-direction std::cerr << "⚠️ Derivative near zero. No convergence." << std::endl;
return v;
}
for (int j = 0; j < num_samples; ++j) { v -= f / df;
double v = (j + 0.5) * dv; // Midpoint sampling in v-direction
// Evaluate the surface point at (u, v) std::cout << "Iteration " << i + 1 << ": v = " << v << ", f = " << f << std::endl;
Eigen::Vector4d pt = point_func(Eigen::Vector2d(u, v)); }
// Get constraint and solve intermediate result std::cerr << "❌ Did not converge within " << max_iterations << " iterations." << std::endl;
auto constraint = m_surface.fetch_curve_constraint_evaluator(internal::parameter_v_t{}, v)(u); return v;
auto result = std::get<internal::parametric_equation_intermediate>(solver_func(std::move(constraint))); }
// Extract partial derivatives ∂S/∂u and ∂S/∂v void SurfaceAreaCalculator::sort_and_unique_with_tol(std::vector<double>& vec, double epsilon)
Eigen::Vector3d dSdu = result.grad_f.col(0); {
Eigen::Vector3d dSdv = result.grad_f.col(1); if (vec.empty()) return;
// Compute the differential area element (cross product norm) std::sort(vec.begin(), vec.end());
double dA = dSdu.cross(dSdv).norm();
// Accumulate area using numerical integration size_t write_index = 0;
area += dA * du * dv; for (size_t read_index = 1; read_index < vec.size(); ++read_index) {
if (std::fabs(vec[read_index] - vec[write_index]) > epsilon) {
++write_index;
vec[write_index] = vec[read_index];
} }
} }
return area; vec.resize(write_index + 1);
} }
} // namespace internal } // namespace internal
Write Preview
Loading…
Cancel
Save