#pragma once
/**
 * @file class_a_surfacing.h
 * @brief CGM Class-A 级曲面工具 — G3 过渡、曲率匹配、高光线、反射线、等照度、诊断、保形修改
 *
 * 对标 Dassault CGM 的 Class-A 曲面功能，提供汽车/航空工业级的曲面质量控制：
 * - G3 连续性过渡（位置+切平面+曲率+曲率变化率连续）
 * - 曲率匹配优化（最小二乘调整控制点使曲率场连续）
 * - 高光线/反射线/等照度分析（工业标准曲面光顺检测）
 * - 综合曲面诊断报告（多项指标评分）
 * - 保形修改（保持特征边界的曲面变形）
 *
 * @ingroup curves
 */

#include "vde/curves/nurbs_curve.h"
#include "vde/curves/nurbs_surface.h"
#include "vde/core/point.h"
#include <vector>
#include <array>
#include <tuple>
#include <utility>
#include <string>
#include <functional>

namespace vde::curves {

using core::Point3D;
using core::Vector3D;

// ═══════════════════════════════════════════════════════════
//  G3 Blend — G3 连续性过渡曲面
// ═══════════════════════════════════════════════════════════

/// 边参数：定义过渡曲面的公共边方向与范围
struct EdgeParams {
    int edge_a = 0;          ///< surf_a 的边界索引 (0=umin, 1=umax, 2=vmin, 3=vmax)
    int edge_b = 0;          ///< surf_b 的边界索引
    double blend_width = 0.1; ///< 过渡带宽度（参数域比例）
    bool align_orientation = true; ///< 是否自动对齐曲面朝向
};

/// G3 连续性过渡结果
struct G3BlendResult {
    NurbsSurface blend_surface;           ///< 过渡曲面
    double max_position_error = 0.0;      ///< 最大位置误差
    double max_tangent_error = 0.0;       ///< 最大切平面角度误差（弧度）
    double max_curvature_error = 0.0;     ///< 最大曲率偏差
    double max_curvature_derivative_error = 0.0; ///< 最大曲率变化率偏差
    bool g0_ok = false;
    bool g1_ok = false;
    bool g2_ok = false;
    bool g3_ok = false;
};

/**
 * @brief G3 连续性过渡曲面
 *
 * 在两面公共边界处构造过渡曲面，使过渡面与两边曲面均达到 G3 连续。
 * 算法：
 *  1. 检测公共边界，提取边界曲线与跨边界导数（1阶、2阶、3阶）
 *  2. 构造 4 排过渡控制点，分别对应位置、切平面、曲率、曲率变化率约束
 *  3. 求解最小二乘系统保证 G3 连续
 *
 * @param surf_a  曲面 A
 * @param surf_b  曲面 B
 * @param params  边参数（边界索引、过渡宽度等）
 * @return G3 过渡结果（含过渡曲面与连续性验证数据）
 *
 * @note 对标 CGM G3 Filling 功能
 * @ingroup curves
 */
[[nodiscard]] G3BlendResult g3_blend(
    const NurbsSurface& surf_a,
    const NurbsSurface& surf_b,
    const EdgeParams& params);

// ═══════════════════════════════════════════════════════════
//  Curvature Matching — 曲率匹配优化
// ═══════════════════════════════════════════════════════════

/// 曲率匹配选项
struct CurvatureMatchOptions {
    double tolerance = 1e-8;       ///< 收敛容差
    int max_iterations = 50;       ///< 最大迭代次数
    int samples_per_edge = 20;     ///< 边界采样点数
    double regularization = 0.01;  ///< Tikhonov 正则化系数
    bool preserve_boundary = true; ///< 是否保持非匹配边界不变
};

/// 曲率匹配结果
struct CurvatureMatchResult {
    NurbsSurface matched_surface_b;   ///< 匹配后的曲面 B
    double initial_rms_curvature_diff = 0.0;  ///< 初始曲率 RMS 差
    double final_rms_curvature_diff = 0.0;    ///< 最终曲率 RMS 差
    int iterations = 0;                       ///< 实际迭代次数
    bool converged = false;                   ///< 是否收敛
    std::vector<Point3D> displacement;        ///< 控制点位移向量
};

/**
 * @brief 曲率匹配优化
 *
 * 调整 surf_b 的控制点使两曲面沿公共边界的曲率场匹配，达到 G2+ 连续。
 * 使用 Levenberg-Marquardt 非线性最小二乘优化，约束边界控制点位移最小。
 *
 * @param surf_a    基准曲面 A（不变）
 * @param surf_b    待匹配曲面 B（将被修改）
 * @param options   优化选项
 * @return 匹配结果（含优化后的曲面 B）
 *
 * @note 对标 CGM Match Surface 功能
 * @ingroup curves
 */
[[nodiscard]] CurvatureMatchResult curvature_matching(
    const NurbsSurface& surf_a,
    const NurbsSurface& surf_b,
    const CurvatureMatchOptions& options = {});

// ═══════════════════════════════════════════════════════════
//  Highlight Lines — 高光线分析
// ═══════════════════════════════════════════════════════════

/// 高光线参数
struct HighlightParams {
    std::vector<Vector3D> light_directions; ///< 多个光源方向（均自动归一化）
    int res_u = 50;   ///< u 方向分辨率
    int res_v = 50;   ///< v 方向分辨率
    int num_bands = 8; ///< 光带数量
    double band_width = 0.05; ///< 带宽（参数域比例）
};

/// 高光线分析结果
struct HighlightResult {
    int res_u, res_v;                                ///< 网格分辨率
    /// 高光线掩码网格 [dir_idx][u][v] ∈ {0,1}
    std::vector<std::vector<std::vector<int>>> band_mask;
    /// 各方向高光线连续性评分 [0,1]
    std::vector<double> continuity_scores;
    /// 综合评分
    double overall_score = 0.0;
    /// 不连续点列表 (u, v, dir_idx)
    std::vector<std::tuple<double, double, int>> discontinuities;
};

/**
 * @brief 高光线分析
 *
 * 模拟平行光源在曲面上的反射高光线，检测光线方向的几何不连续。
 * 高光线 ≤ C1 的断裂 = 曲面在此处仅有 C1 连续。
 * 高光线 C2 断裂 = 曲率不连续。
 *
 * @param surface   NURBS 曲面
 * @param light_dirs 光源方向列表
 * @param res_u     u 方向分辨率
 * @param res_v     v 方向分辨率
 * @param num_bands 光带数量
 * @return 高光线分析结果
 *
 * @note 对标 CGM Highlight Lines Analysis
 * @ingroup curves
 */
[[nodiscard]] HighlightResult highlight_lines(
    const NurbsSurface& surface,
    const std::vector<Vector3D>& light_dirs,
    int res_u, int res_v, int num_bands = 8);

// ═══════════════════════════════════════════════════════════
//  Reflection Lines — 反射线分析
// ═══════════════════════════════════════════════════════════

/// 反射线分析结果
struct ReflectionResult {
    int res_u, res_v;                              ///< 网格分辨率
    /// 反射线方向场 grid[u][v] = 反射光线方向
    std::vector<std::vector<Vector3D>> reflection_field;
    /// 反射线扭曲度 grid[u][v] ∈ [0,∞)
    std::vector<std::vector<double>> distortion;
    double max_distortion = 0.0;
    double mean_distortion = 0.0;
    bool is_fair = false; ///< 是否判定为光顺
};

/**
 * @brief 反射线分析
 *
 * 计算虚拟环境（点光源或平行光）在曲面上的镜面反射线模式。
 * 通过计算反射方向场的扭曲度评估曲面光顺性。
 *
 * 算法：
 *  law_of_reflection: r = 2(n·l)n - l （入射光方向 l，法线 n）
 *  扭曲度 = reflection_field 的方向导数幅值
 *
 * @param surface  NURBS 曲面
 * @param eye      观察点位置
 * @param light    光源位置（若为平行光则传入方向向量远端）
 * @param res_u    u 方向分辨率
 * @param res_v    v 方向分辨率
 * @return 反射线分析结果
 *
 * @note 对标 CGM Reflection Lines / Isophotes Analysis
 * @ingroup curves
 */
[[nodiscard]] ReflectionResult reflection_lines(
    const NurbsSurface& surface,
    const Point3D& eye,
    const Point3D& light,
    int res_u, int res_v);

// ═══════════════════════════════════════════════════════════
//  Iso-Photes — 等照度分析
// ═══════════════════════════════════════════════════════════

/// 等照度线分析结果
struct IsoPhoteResult {
    int res_u, res_v;                              ///< 网格分辨率
    /// 照度值网格 ∈ [-1, 1]（光源方向·法线的点积）
    std::vector<std::vector<double>> illumination;
    /// 等照度线密度分布（直方图）
    std::vector<double> iso_histogram;
    /// 等照度线连续性 — 梯度突变点
    std::vector<std::pair<double,double>> gradient_discontinuities;
    double grad_max = 0.0;
    double grad_mean = 0.0;
};

/**
 * @brief 等照度分析
 *
 * 计算光源照射下曲面上的等照度线（恒定 cosθ 线）。
 * 等照度线在 G1 连续处平滑，在仅有 C0 处出现转折。
 * 等照度线的质量直接反映曲面光顺度：
 *  - 间距均匀 → 曲率均匀
 *  - 无抖动 → 高阶连续性好
 *
 * @param surface  NURBS 曲面
 * @param res      分辨率（res × res 采样点）
 * @return 等照度分析结果
 *
 * @note 对标 CGM Iso-Photes Analysis
 * @ingroup curves
 */
[[nodiscard]] IsoPhoteResult iso_photes(
    const NurbsSurface& surface,
    int res);

// ═══════════════════════════════════════════════════════════
//  Surface Diagnosis — 综合曲面诊断报告
// ═══════════════════════════════════════════════════════════

/// 诊断等级
enum class DiagnosisGrade {
    A_CLASS = 0,  ///< A 级曲面 — 通过所有检测
    B_CLASS = 1,  ///< B 级曲面 — 可接受，有轻微缺陷
    C_CLASS = 2,  ///< C 级曲面 — 需要修复
    FAIL = 3      ///< 不合格
};

/// 子项诊断结果
struct DiagnosisItem {
    std::string name;        ///< 诊断项名称
    double value = 0.0;      ///< 测量值
    double threshold = 0.0;  ///< 阈值
    bool pass = false;       ///< 是否通过
    std::string description; ///< 描述
};

/// 综合曲面诊断报告
struct SurfaceDiagnosisReport {
    DiagnosisGrade grade = DiagnosisGrade::FAIL;

    /// 高斯曲率诊断
    DiagnosisItem gaussian_range;       ///< 高斯曲率范围
    DiagnosisItem gaussian_continuity;  ///< 高斯曲率连续性

    /// 平均曲率诊断
    DiagnosisItem mean_range;           ///< 平均曲率范围

    /// 高光线诊断
    DiagnosisItem highlight_score;      ///< 高光线评分

    /// 反射线诊断
    DiagnosisItem reflection_distortion; ///< 反射扭曲度

    /// 等照度诊断
    DiagnosisItem isophote_gradient;    ///< 等照度梯度

    /// 法线连续性
    DiagnosisItem normal_jump;          ///< 最大法线跳跃（角度°）

    /// 主曲率方向场
    DiagnosisItem principal_direction_flow; ///< 主方向场平滑度

    /// 可展性
    DiagnosisItem developability;       ///< 可展面检测

    double overall_score = 0.0;         ///< 综合评分 [0, 100]
    std::string recommendation;         ///< 修复建议
};

/**
 * @brief 综合曲面诊断报告
 *
 * 对单张 NURBS 曲面进行多项工业级质量检测，生成诊断报告：
 *  1. 高斯/平均曲率范围及连续性
 *  2. 高光线分析（4方向）
 *  3. 反射线扭曲度
 *  4. 等照度梯度
 *  5. 法线连续性（内部采样）
 *  6. 主曲率方向场流畅度
 *  7. 可展面检测
 *
 * @param surface  NURBS 曲面
 * @return 综合诊断报告
 *
 * @note 对标 CGM Surface Diagnosis / CATIA Generative Shape Design
 * @ingroup curves
 */
[[nodiscard]] SurfaceDiagnosisReport surface_diagnosis(
    const NurbsSurface& surface);

// ═══════════════════════════════════════════════════════════
//  Shape Modification — 保形修改
// ═══════════════════════════════════════════════════════════

/// 约束点
struct ShapeConstraint {
    Point3D target_point;     ///< 目标位置
    double u = 0.0;           ///< 曲面参数 u
    double v = 0.0;           ///< 曲面参数 v
    double weight = 1.0;      ///< 约束权重
    bool fix_position = true; ///< true=位置约束, false=仅法向约束
    Vector3D target_normal;   ///< 目标法向（仅 fix_position=false 时使用）
};

/// 保形修改结果
struct ShapeModificationResult {
    NurbsSurface modified_surface;    ///< 修改后的曲面
    std::vector<Point3D> control_displacements; ///< 控制点位移
    double max_displacement = 0.0;   ///< 最大控制点位移
    double energy_preserved_ratio = 0.0; ///< 应变能保持率 (0~1)
    double boundary_deviation = 0.0; ///< 边界偏差
    int iterations = 0;              ///< 迭代次数
};

/**
 * @brief 保形修改
 *
 * 在约束条件下修改曲面控制点，同时最小化曲面变形能量，
 * 保持曲面的整体形状特征（边界、特征线等）。
 *
 * 算法：基于薄板样条能量最小化的约束优化
 *  E = Σ (||S(u_i,v_i) - target_i||² · w_i) + λ · bending_energy(S)
 *
 * @param surface     NURBS 曲面
 * @param constraints 约束点列表
 * @return 修改结果
 *
 * @note 对标 CGM Shape Modification / Control Point Morphing
 * @ingroup curves
 */
[[nodiscard]] ShapeModificationResult shape_modification(
    const NurbsSurface& surface,
    const std::vector<ShapeConstraint>& constraints);

// ═══════════════════════════════════════════════════════════
//  G3 Blend with Constraints — 带约束的 G3 过渡
// ═══════════════════════════════════════════════════════════

/// G0/G1/G2/G3 连续性约束边参数
struct G3ConstraintEdges {
    EdgeParams g0_edge;  ///< G0 位置约束边
    EdgeParams g1_edge;  ///< G1 切平面约束边
    EdgeParams g2_edge;  ///< G2 曲率约束边
    EdgeParams g3_edge;  ///< G3 曲率变化率约束边
    bool enforce_g0 = true;
    bool enforce_g1 = true;
    bool enforce_g2 = true;
    bool enforce_g3 = true;
};

/**
 * @brief 带约束的 G3 过渡曲面
 *
 * 在两面之间构造过渡曲面，并精确满足指定的 G0/G1/G2/G3 约束。
 * 与 g3_blend 相比，此函数允许用户显式指定每条边的连续性等级。
 *
 * 算法：
 *  1. 根据各约束边采样控制点
 *  2. 构造最小二乘系统满足 G0→G3 逐级约束
 *  3. 使用 Lagrange 乘子法确保约束精确满足
 *
 * @param surf_a  曲面 A
 * @param surf_b  曲面 B
 * @param edges   显式 G0/G1/G2/G3 约束边参数
 * @return G3 过渡结果
 *
 * @ingroup curves
 */
[[nodiscard]] G3BlendResult g3_blend_with_constraints(
    const NurbsSurface& surf_a,
    const NurbsSurface& surf_b,
    const G3ConstraintEdges& edges);

// ═══════════════════════════════════════════════════════════
//  Surface Energy Minimization — 薄板能量最小化
// ═══════════════════════════════════════════════════════════

/// 能量最小化选项
struct EnergyMinimizationOptions {
    double bending_weight = 1.0;     ///< 弯曲能量权重
    double membrane_weight = 0.1;    ///< 薄膜能量权重
    double constraint_weight = 10.0;  ///< 约束满足权重
    int max_iterations = 100;         ///< 最大迭代次数
    double tolerance = 1e-6;          ///< 收敛容差
    double regularization = 1e-4;     ///< Tikhonov 正则化
    bool preserve_boundary = true;    ///< 是否保持边界不变
};

/// 能量最小化结果
struct EnergyMinimizationResult {
    NurbsSurface optimized_surface;          ///< 优化后的曲面
    double initial_bending_energy = 0.0;     ///< 初始弯曲能量
    double final_bending_energy = 0.0;       ///< 最终弯曲能量
    double initial_membrane_energy = 0.0;    ///< 初始薄膜能量
    double final_membrane_energy = 0.0;      ///< 最终薄膜能量
    int iterations = 0;                      ///< 实际迭代次数
    bool converged = false;                  ///< 是否收敛
    std::vector<Point3D> control_displacements; ///< 控制点位移
    double max_displacement = 0.0;           ///< 最大控制点位移
};

/**
 * @brief 薄板能量最小化曲面优化
 *
 * 最小化薄板弯曲能量（二阶导数平方积分）和薄膜能量（一阶导数平方积分），
 * 同时满足约束条件，实现曲面的全局光顺。
 *
 * 能量泛函：
 *  E = w_bend * ∫(κ₁² + κ₂²)dA + w_membrane * ∫(‖∂S/∂u‖² + ‖∂S/∂v‖²)dA
 *    + w_constraint * Σ‖S(u_i, v_i) - target_i‖²
 *
 * 算法：基于控制点的 Gauss-Newton 非线性最小二乘优化。
 *
 * @param surface     输入 NURBS 曲面
 * @param constraints 约束点列表（可为空，仅做光顺）
 * @param options     优化选项
 * @return 能量最小化结果
 *
 * @note 对标 ICEM Surf / Alias AutoStudio 的 Surface Fairing
 * @ingroup curves
 */
[[nodiscard]] EnergyMinimizationResult surface_energy_minimization(
    const NurbsSurface& surface,
    const std::vector<ShapeConstraint>& constraints,
    const EnergyMinimizationOptions& options = {});

// ═══════════════════════════════════════════════════════════
//  Curvature Continuity Optimization — 曲率连续性迭代精化
// ═══════════════════════════════════════════════════════════

/// 曲率连续性优化选项
struct CurvatureContinuityOptions {
    double position_tolerance = 1e-7;    ///< 位置连续性容差
    double tangent_tolerance = 1e-5;     ///< 切平面角度容差 (rad)
    double curvature_tolerance = 1e-4;   ///< 曲率偏差容差
    int max_iterations = 50;             ///< 最大迭代次数
    int boundary_samples = 30;           ///< 边界采样点数
    double step_size = 0.05;             ///< 步长因子
    bool optimize_both_surfaces = false; ///< 是否同时优化两面
};

/// 曲率连续性优化结果
struct CurvatureContinuityResult {
    NurbsSurface optimized_a;                ///< 优化后的曲面 A（若 optimize_both_surfaces=true）
    NurbsSurface optimized_b;                ///< 优化后的曲面 B
    double initial_g0_error = 0.0;           ///< 初始 G0 误差
    double final_g0_error = 0.0;             ///< 最终 G0 误差
    double initial_g1_error = 0.0;           ///< 初始 G1 误差 (rad)
    double final_g1_error = 0.0;             ///< 最终 G1 误差 (rad)
    double initial_g2_error = 0.0;           ///< 初始 G2 误差
    double final_g2_error = 0.0;             ///< 最终 G2 误差
    int iterations = 0;                      ///< 实际迭代次数
    bool converged = false;                  ///< 是否收敛
};

/**
 * @brief 曲率连续性迭代精化优化
 *
 * 在两曲面公共边界处迭代调整控制点，逐步提升 G0→G1→G2 连续性。
 * 每次迭代沿公共边界采样，计算当前连续性误差，梯度下降调整影响区域控制点。
 *
 * 算法：
 *  1. 检测公共边界
 *  2. 沿边界采样，计算各点 G0/G1/G2 误差
 *  3. 梯度下降调整边界面附近控制点
 *  4. 重新评估，调整步长（adaptive step size）
 *  5. 收敛检测
 *
 * @param surf_a  曲面 A
 * @param surf_b  曲面 B
 * @param options 优化选项
 * @return 优化结果（含优化后的曲面和收敛历史）
 *
 * @note 对标 CATIA GSD Join / ICEM Surf Match
 * @ingroup curves
 */
[[nodiscard]] CurvatureContinuityResult curvature_continuity_optimization(
    const NurbsSurface& surf_a,
    const NurbsSurface& surf_b,
    const CurvatureContinuityOptions& options = {});

// ═══════════════════════════════════════════════════════════
//  Reflection Line Discontinuity — 反射线不连续检测+修复
// ═══════════════════════════════════════════════════════════

/// 反射线不连续检测参数
struct ReflectionDiscontinuityParams {
    Vector3D light_direction;          ///< 光源方向（归一化）
    int res_u = 40;                    ///< u 方向分辨率
    int res_v = 40;                    ///< v 方向分辨率
    double discontinuity_threshold = 0.05; ///< 不连续判定阈值（方向变化率）
    int max_fix_iterations = 20;       ///< 最大修复迭代次数
    double fix_strength = 0.1;         ///< 修复强度
    bool auto_fix = true;              ///< 是否自动修复
};

/// 反射线不连续检测结果
struct ReflectionDiscontinuityResult {
    int res_u = 0, res_v = 0;          ///< 网格分辨率
    /// 反射场 grid
    std::vector<std::vector<Vector3D>> reflection_field;
    /// 不连续强度网格 ∈ [0, ∞)
    std::vector<std::vector<double>> discontinuity_map;
    /// 不连续点列表 (u, v, intensity)
    std::vector<std::tuple<double, double, double>> discontinuities;
    int total_discontinuities = 0;     ///< 不连续点总数
    double max_discontinuity = 0.0;    ///< 最大不连续强度
    double mean_discontinuity = 0.0;   ///< 平均不连续强度
    bool is_smooth = false;            ///< 是否判定为光顺
    NurbsSurface fixed_surface;         ///< 修复后的曲面（若 auto_fix=true）
    bool fixed = false;                ///< 是否已修复
    double fix_residual = 0.0;         ///< 修复残差
};

/**
 * @brief 反射线不连续检测与修复
 *
 * 计算曲面在给定光源方向下的反射线场，检测反射方向的局部突变点，
 * 并可选择性地自动修复（通过局部控制点微调消除不连续）。
 *
 * 算法：
 *  - 检测：计算反射方向场的梯度幅值，超过阈值的标记为不连续
 *  - 修复：对不连续区域的控制点施加局部高斯加权偏移，最小化反射梯度
 *
 * @param surface   NURBS 曲面
 * @param params    检测参数（光源方向、分辨率、阈值、修复选项）
 * @return 反射线不连续检测与修复结果
 *
 * @note 对标 CGM Reflection Line Analysis + Repair
 * @ingroup curves
 */
[[nodiscard]] ReflectionDiscontinuityResult reflection_line_discontinuity(
    const NurbsSurface& surface,
    const ReflectionDiscontinuityParams& params);

}  // namespace vde::curves