CGAL生成融合双椭球SurfaceMesh

 生成曲面如下图所示

#include <iostream>
#include <vector>
#include <cmath>
#include<fstream>

#include<glm/glm.hpp>

#include <vtkSmartPointer.h>
#include <vtkPoints.h>
#include <vtkPolyData.h>
#include <vtkDelaunay3D.h>

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/Point_set_3/IO.h>
#include <CGAL/remove_outliers.h>
#include <CGAL/grid_simplify_point_set.h>
#include <CGAL/jet_smooth_point_set.h>
#include <CGAL/jet_estimate_normals.h>
#include <CGAL/mst_orient_normals.h>
#include <CGAL/poisson_surface_reconstruction.h>
#include <CGAL/Advancing_front_surface_reconstruction.h>
#include <CGAL/Scale_space_surface_reconstruction_3.h>
#include <CGAL/Scale_space_reconstruction_3/Jet_smoother.h>
#include <CGAL/Scale_space_reconstruction_3/Advancing_front_mesher.h>
#include <CGAL/Polygon_mesh_processing/manifoldness.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <CGAL/Polygon_mesh_processing/smooth_shape.h>
#include <cstdlib>
#include <CGAL/IO/OBJ.h>
#include <CGAL/IO/OFF.h>

#include <CGAL/Alpha_shape_3.h>
#include <CGAL/Alpha_shape_vertex_base_3.h>
#include <CGAL/Alpha_shape_cell_base_3.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Polygon_mesh_processing/connected_components.h>
#include <CGAL/Polygon_mesh_processing/repair.h>

#include <iostream>
#include <Eigen/Dense>
#include <cmath>

typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point_3;
typedef Kernel::FT FT;
typedef CGAL::Point_set_3<Point_3> Point_set;
typedef CGAL::Surface_mesh<Point_3> Surface_mesh;

// x表示沿长轴正向的距离
glm::dvec3 getEllipsoidPoint(double a, double b, double c, double x)
{
    // phi = 90 方位角是90°
    //double x = a * sin(theta);
    double theta = asin(x / a);
    double y = 0;
    double z = c * cos(theta);
    return glm::dvec3(x, y, z);
}

double calcYinXoY(double x, double a, double b)
{
    double phi = acos(x / a);
    return (b * sin(phi));
}

// 已知YOZ平面内，Y方向偏移，求高度Z值
double calcZinYoZ(double y, double b, double c)
{
    double theta = asin(y / b);
    return (c * cos(theta));
}

// 已知XoY平面内，Y方向偏移，求X值
double calcXinXoY(double y, double a, double b)
{
    double phi = asin(y / b);
    return (a * cos(phi));
}

// 不能全部的点z坐标都乘，越接近中线ypos，乘的数要越接近1
void scaleZ(std::vector<glm::dvec3>& data, double ypos, double scale, bool yneg = true)
{
    for (auto& d : data)
    {
        if (yneg && d.y < 0)
        {

        d.z *= (-d.y / ypos * scale + (ypos + d.y) / ypos);
        }
    }
}

// 不能全部的点z坐标都乘，越接近中线ypos，乘的数要越接近1
void scaleX(std::vector<glm::dvec3>& data, double xpos, double scale, double a, double b, bool yneg = true)
{
    for (auto& d : data)
    {
        
        if (yneg && d.y < 0)
        {
            auto yscale = scale * fabs(d[1]) / b + (b - fabs(d[1])) / b;
            if (fabs(d.x) < xpos)
                d.x *= yscale;
            else 
               d.x *= (a - fabs(d.x)) / (a-xpos) * yscale + (fabs(d.x) - xpos)/ (a-xpos);
        }
    }
}

void scaleY(std::vector<glm::dvec3>& data, double scale, bool neg = true)
{
    for (auto& d : data)
    {
        if (neg && d.y < 0)
        {
            d.y *= scale;
        }
        else if (!neg && d.y > 0)
        {
            d.y *= scale;
        }
    }
}

void clipY(std::vector<glm::dvec3>& data, double thresh_y, bool neg = true)
{
    auto it = std::remove_if(data.begin(), data.end(), [&](const glm::dvec3& point) {
        if (neg && point.y < -thresh_y)
            return true;
        if (!neg && point.y > thresh_y)
            return true;
        return false;
        });

    // 删除移动到末尾的元素
    data.erase(it, data.end());
}

void mirrorY(std::vector<glm::dvec3>& data, double y_pos, bool neg = true)
{
    std::vector<glm::dvec3> tmp;
    for (auto d : data)
    {
        if(neg)
            tmp.push_back(glm::dvec3(d.x, d.y + 2.0 * (-y_pos - d.y), d.z));
        else
            tmp.push_back(glm::dvec3(d.x, d.y + 2.0 * (y_pos - d.y), d.z));
    }
    data.insert(data.end(), tmp.begin(), tmp.end());
}

// 生成椭球表面点云数据
std::vector<glm::dvec3> generateEllipsoidPoints(double a, double b, double c, int thetaSteps = 100, int phiSteps = 150) {
    std::vector<glm::dvec3> points;

    // 越远离XOY平面，采样点数可以越少

    // 角度步长
    double thetaStep = M_PI / thetaSteps;
    double phiStep = 2 * M_PI / phiSteps;

    // 生成点云
    for (int i = 0; i <= thetaSteps; ++i) {
        double theta = i * thetaStep; // 极角

        int n_phi = round(9.0 * phiSteps / 10.0 * (double)(thetaSteps / 2.0 - fabs(thetaSteps / 2.0 - i)) / thetaSteps * 2.0 + phiSteps / 10.0);
        phiStep = 2 * M_PI / n_phi;
        for (int j = 0; j < n_phi; ++j) {
            double phi = j * phiStep; // 方位角

            // 计算点的坐标
            double x = a * sin(theta) * cos(phi);
            double y = b * sin(theta) * sin(phi);
            double z = c * cos(theta);

            // 添加到点云
            points.push_back({ x, y, z });
        }
    }

    return points;
}

// 保存点云数据到 txt 文件
void savePointCloudToTxt(const std::vector<glm::dvec3>& points, const std::string& filename) {
    // 打开文件
    std::ofstream outFile(filename);
    if (!outFile.is_open()) {
        std::cerr << "无法打开文件: " << filename << std::endl;
        return;
    }

    // 写入点云数据
    for (const auto& point : points) {
        outFile << point.x << " " << point.y << " " << point.z << "\n";
    }

    // 关闭文件
    outFile.close();
    std::cout << "点云数据已保存到: " << filename << std::endl;
}

typedef std::array<std::size_t, 3> Facet;
struct Construct {
    Surface_mesh& mesh;

    template < typename PointIterator>
    Construct(Surface_mesh& mesh, PointIterator b, PointIterator e)
        : mesh(mesh)
    {
        for (; b != e; ++b) {
            boost::graph_traits<Surface_mesh>::vertex_descriptor v;
            v = add_vertex(mesh);
            mesh.point(v) = *b;
        }
    }

    Construct& operator=(const Facet f)
    {
        typedef boost::graph_traits<Surface_mesh>::vertex_descriptor vertex_descriptor;
        typedef boost::graph_traits<Surface_mesh>::vertices_size_type size_type;
        mesh.add_face(vertex_descriptor(static_cast<size_type>(f[0])),
            vertex_descriptor(static_cast<size_type>(f[1])),
            vertex_descriptor(static_cast<size_type>(f[2])));
        return *this;
    }

    Construct&
        operator*() { return *this; }

    Construct&
        operator++() { return *this; }

    Construct
        operator++(int) { return *this; }

};

// 假设输入是 std::vector<glm::dvec3> points
void reconstructSurface(std::vector<glm::dvec3>& points) {
    // 读取点云数据
    std::vector<Point_3> point_set;
    namespace PMP = CGAL::Polygon_mesh_processing;
    for (int i = 0; i < points.size(); i++)
    {
        point_set.push_back(Point_3(points[i].x, points[i].y, points[i].z));
    }

    Surface_mesh m;
    Construct construct(m, point_set.begin(), point_set.end());

    CGAL::advancing_front_surface_reconstruction(point_set.begin(),
        point_set.end(),
        construct);

    int nb_iterations = 15; // 平滑迭代次数
    CGAL::Polygon_mesh_processing::smooth_shape(m,
        nb_iterations,
        CGAL::Polygon_mesh_processing::parameters::vertex_point_map(m.points())
    );

    printf("reconstruct done!\n");

    // 保存网格到OFF文件
    std::ofstream off_file("output.off");
    if (off_file)
    {
        if (!CGAL::IO::write_OFF(off_file, m))
        {
            std::cerr << "Error: could not write to output.off" << std::endl;
            return ;
        }
        std::cout << "Mesh saved to output.off" << std::endl;
    }
    else
    {
        std::cerr << "Error: could not open output.off for writing" << std::endl;
        return ;
    }



}

struct GeoParams
{
    double d = 0.0;    // 针距
    double time = 0.0;
    double power = 0.0;

    double fbMid = 0.0;
    double fbMax = 0.0;
    double width = 0.0;
    double udMid = 0.0;
    double udMax = 0.0;
    double forthMid = 0.0;
    double forthMax = 0.0;
};

int main() {
    GeoParams gps;

    // 双针
    gps.d = 25;
    gps.fbMid = 23;
    gps.fbMax = 36;
    gps.width = 50;
    gps.udMid = 16;
    gps.udMax = 23;

    double a = gps.fbMax / 2.0;               // X 轴方向 长轴方向
    double b = gps.width / 2.0 - gps.d / 2.0; // Y 轴方向 俯视图 短轴方向
    double c = gps.udMax / 2.0;               // Z 轴方向 高度方向 厚度方向
    
    // 生成椭球表面点云数据
    std::vector<glm::dvec3> ellipsoidPoints = generateEllipsoidPoints(a, b, c);
    auto y0 = calcYinXoY(gps.fbMid / 2.0, a, b);
    auto y1 = gps.d / 2.0;
    auto scale = y1 / y0;
    scaleY(ellipsoidPoints, scale);

    auto z0 = calcZinYoZ(gps.d / 2.0 / scale, b, c);
    auto z1 = gps.udMid / 2.0;
    auto scale_z = z1 / z0;
    // 不能全部的点z坐标都乘，越接近中线，乘的数要越接近1
    scaleZ(ellipsoidPoints, gps.d / 2.0, scale_z);

    auto x0 = calcXinXoY(gps.d / 2.0, a, b);
    auto x1 = a - gps.forthMax + gps.forthMid;
    auto scale_x = x1 / x0;
    // 不能全部的点x坐标都乘，越接近长轴顶端，乘的数要越接近1
    scaleX(ellipsoidPoints, x0, scale_x, a);

    ellipsoidPoints.push_back(glm::dvec3(0, 0, c));
    ellipsoidPoints.push_back(glm::dvec3(0, 0, -c));
    
    
    clipY(ellipsoidPoints, gps.d / 2.0);
    mirrorY(ellipsoidPoints, gps.d / 2.0);
    reconstructSurface(ellipsoidPoints);
    return 0;
}