ImplicitSurfaceNetwork/surface_integral/src/SurfaceIntegrator.cpp

#include "SurfaceIntegrator.hpp"
#include "quadrature.hpp"
#include <Eigen/Geometry> // For vector and cross product operations
#include <cmath>          // For math functions like sqrt
#include <set>

namespace internal
{

// Constructor 1: Initialize only with a reference to the surface
SurfaceAreaCalculator::SurfaceAreaCalculator(const subface& surface) : m_surface(surface) {}

// Constructor 2: Initialize with surface and u-breaks (e.g., trimming curves)
SurfaceAreaCalculator::SurfaceAreaCalculator(const subface&               surface,
                                             const stl_vector_mp<double>& u_breaks,
                                             double                       umin,
                                             double                       umax,
                                             double                       vmin,
                                             double                       vmax)
    : m_surface(surface), m_u_breaks(u_breaks), Umin(umin), Umax(umax), Vmin(vmin), Vmax(vmax)
{
}

SurfaceAreaCalculator::SurfaceAreaCalculator(const subface& surface, const parametric_plane& uv_plane)
    : m_surface(surface), m_uv_plane(uv_plane), Umin(0.0), Umax(0.0), Vmin(0.0), Vmax(0.0)
{
    if (!uv_plane.chain_vertices.empty()) {
        // 初始化为第一个点的坐标
        double min_u = uv_plane.chain_vertices[0].x();
        double max_u = uv_plane.chain_vertices[0].x();
        double min_v = uv_plane.chain_vertices[0].y();
        double max_v = uv_plane.chain_vertices[0].y();

        // 遍历所有链顶点
        for (const auto& pt : uv_plane.chain_vertices) {
            double u = pt.x();
            double v = pt.y();

            if (u < min_u) min_u = u;
            if (u > max_u) max_u = u;
            if (v < min_v) min_v = v;
            if (v > max_v) max_v = v;
        }

        Umin = min_u;
        Umax = max_u;
        Vmin = min_v;
        Vmax = max_v;
    } else {
        // 没有顶点时使用默认范围 [0, 1] × [0, 1]
        Umin = 0.0;
        Umax = 1.0;
        Vmin = 0.0;
        Vmax = 1.0;
    }

    std::set<uint32_t> unique_vertex_indices;

    // 插入所有类型的顶点索引
    unique_vertex_indices.insert(uv_plane.singularity_vertices.begin(), uv_plane.singularity_vertices.end());
    unique_vertex_indices.insert(uv_plane.polar_vertices.begin(), uv_plane.polar_vertices.end());
    unique_vertex_indices.insert(uv_plane.parallel_start_vertices.begin(), uv_plane.parallel_start_vertices.end());

    std::set<double> unique_u_values;

    for (uint32_t idx : unique_vertex_indices) {
        if (idx < uv_plane.chain_vertices.size()) {
            double u = uv_plane.chain_vertices[idx].x();
            unique_u_values.insert(u);
        }
    }

    // 转换为 vector 并排序（set 默认是有序的）
    m_u_breaks = stl_vector_mp<double>(unique_u_values.begin(), unique_u_values.end());
}

// Set u-breaks (optional trimming or partitioning lines)
void SurfaceAreaCalculator::set_ubreaks(const stl_vector_mp<double>& u_breaks) { m_u_breaks = u_breaks; }

// Main entry point to compute surface area
template <typename Func>
double SurfaceAreaCalculator::calculate(Func&& func, int gauss_order) const
{
    // 在u方向进行高斯积分
    auto u_integrand = [&](double u) {
        // 对每个u，找到v方向的精确交点
        std::vector<double> v_breaks = find_vertical_intersections(u);

        // 在v方向进行高斯积分
        auto v_integrand = [&](double v) {
            // 判断点是否在有效域内
            if (IsPointInsideself(u, v, self.outerEdges, self.innerEdges, self.Umin, self.Umax, self.Vmin, self.Vmax)) {
                try {
                    gp_Pnt p;
                    gp_Vec dU, dV;
                    surface->D1(u, v, p, dU, dV);

                    const double jacobian = dU.Crossed(dV).Magnitude();
                    return func(u, v, p, dU, dV) * jacobian;
                } catch (...) {
                    return 0.0; // 跳过奇异点
                }
            }
            return 0.0;         // 点不在有效域内
        };

        double v_integral = 0.0;
        for (size_t i = 0; i < v_breaks.size() - 1; ++i) {
            double a = v_breaks[i];
            double b = v_breaks[i + 1];

            // 检查区间中点是否有效
            double mid_v = (a + b) / 2.0;
            if (IsPointInsideself(u, mid_v, self.outerEdges, self.innerEdges, self.Umin, self.Umax, self.Vmin, self.Vmax)) {
                v_integral += integrate_1D(a, b, v_integrand, gauss_order);
            } else {
                std::cout << "uv out of domain: (" << u << "," << mid_v << ")" << std::endl;
            }
        }

        return v_integral;
    };

    // 在u方向积分
    double integral = 0.0;
    for (size_t i = 0; i < u_breaks.size() - 1; ++i) {
        double a = u_breaks[i];
        double b = u_breaks[i + 1];

        // 检查区间中点是否有效
        double mid_u           = (a + b) / 2.0;
        auto   v_intersections = find_vertical_intersections(mid_u, self.outerEdges, self.Vmin, self.Vmax);

        if (!v_intersections.empty()) { // 确保该u区间有有效区域

            integral += integrate_1D(a, b, u_integrand, gauss_order, is_u_near_singularity(mid_u));
        }
    }

    return integral;
}

// 直线u=u_val与边界的交点
std::vector<double> SurfaceAreaCalculator::find_vertical_intersections(double u_val)
{
    std::vector<double>  intersections;
    std::vector<int32_t> uPositionFlags;
    uPositionFlags.reserve(m_uv_plane.chain_vertices.size());

    std::transform(m_uv_plane.chain_vertices.begin(),
                   m_uv_plane.chain_vertices.end(),
                   std::back_inserter(uPositionFlags),
                   [&](const auto& currentVertex) {
                       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_val);
                    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());
                }
            }
        }
    }

    // 去重排序
    sort_and_unique_with_tol(intersections);
    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 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++; }
                }
            }
            /*
            case v1.x() == v2.x() == u, do not count as intersection.but will cout in next iteration
            */
        }
    }
    return intersection_count % 2 == 1; // in domain
}

bool is_u_near_singularity(double u, double tol = 1e-6)
{
    for (auto idx : m_uv_plane.singularity_vertices) {
        double singular_u = m_uv_plane.chain_vertices[idx].x();
        if (std::abs(u - singular_u) < tol) { return true; }
    }
    // 可扩展：判断是否靠近极点、极性顶点等
    for (auto idx : m_uv_plane.polar_vertices) {
        double polar_u = m_uv_plane.chain_vertices[idx].x();
        if (std::abs(u - polar_u) < tol) { return true; }
    }
    return false;
}

void SurfaceAreaCalculator::sort_and_unique_with_tol(std::vector<double>& vec, double epsilon)
{
    if (vec.empty()) return;

    std::sort(vec.begin(), vec.end());

    size_t write_index = 0;
    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];
        }
    }

    vec.resize(write_index + 1);
}

// 牛顿法求解器
double newton_method(const std::function<equation_intermediate_t(double)>& F,
                     double                                                v_initial,
                     double                                                tolerance,
                     int                                                   max_iterations)
{
    double v = v_initial;

    for (int i = 0; i < max_iterations; ++i) {
        equation_intermediate_t res = F(v);
        double                  f   = res.f;
        double                  df  = res.df;

        if (std::abs(f) < tolerance) {
            std::cout << "✅ Converged in " << i + 1 << " iterations. v = " << v << std::endl;
            return v;
        }

        if (std::abs(df) < 1e-10) {
            std::cerr << "⚠️ Derivative near zero. No convergence." << std::endl;
            return v;
        }

        v -= f / df;

        std::cout << "Iteration " << i + 1 << ": v = " << v << ", f = " << f << std::endl;
    }

    std::cerr << "❌ Did not converge within " << max_iterations << " iterations." << std::endl;
    return v;
}


} // namespace internal