实现斐波那契球面与Bowyer-Watson算法时遭遇段错误求助

问题背景 #

调试生成斐波那契球面并通过Bowyer-Watson算法做三角剖分的代码时，持续出现段错误，在线调试工具也无法解决，以下是原代码的问题分析及修复方案。

原代码中的关键错误 #

• 数组越界：add_vertex里的edge edges[3]固定大小为3，但处理多个符合条件的三角形时，添加的边数会远超3，直接触发段错误，需改为动态分配。
• 循环变量重定义：add_vertex中判断边唯一性的嵌套循环里，内外层都用了i作为循环变量，导致变量覆盖，逻辑混乱，需修改内层变量名。
• 代码位置错误：add_vertex中创建新三角形、重分配数组的代码被放在边遍历循环内，导致每处理一条边就执行一次数组操作并返回，逻辑完全错误，需移到循环外。
• 顶点赋值覆盖：create_fibonacci_sphere中给point_list[n]连续两次赋值，第二次覆盖了第一次的y坐标，导致数据丢失，需修正赋值逻辑。
• 语法结构错误：triangulate中移除超级三角形的for循环缺少闭合}，create_fibonacci_sphere最后一行printf在函数外，均会导致编译或运行错误。
• 浮点数比较不当：vertex_equals用==直接比较浮点数，因精度问题会导致顶点匹配错误，需改用差值极小值判断。
• 参数类型不匹配：triangulate调用时传入size_t类型的arr_size，但函数参数是int，需显式转换。

修复后的完整代码 #

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

// 顶点结构体
typedef struct Vertex
{
    double x;
    double y;
} vertex;

// 比较两个顶点（考虑浮点数精度）
int vertex_equals(vertex v1, vertex v2)
{
    const double eps = 1e-9;
    return fabs(v1.x - v2.x) < eps && fabs(v1.y - v2.y) < eps;
}

// 边结构体
typedef struct Edge
{
    vertex v0;
    vertex v1;
} edge;

// 比较两条边（无向边）
int edge_equals(edge e1, edge e2)
{
    return (vertex_equals(e1.v0, e2.v0) && vertex_equals(e1.v1, e2.v1)) ||
           (vertex_equals(e1.v0, e2.v1) && vertex_equals(e1.v1, e2.v0));
}

// 外接圆结构体
typedef struct Circumcircle
{
    vertex c;
    double r;
} circum_circ;

// 三角形结构体
typedef struct Triangle
{
    vertex v0;
    vertex v1;
    vertex v2;
    circum_circ circumcirc;
} triangle;

// 计算三角形的外接圆
circum_circ calc_circum_circ(vertex v0, vertex v1, vertex v2)
{
    double A = v1.x - v0.x;
    double B = v1.y - v0.y;
    double C = v2.x - v0.x;
    double D = v2.y - v0.y;

    double E = A * (v0.x + v1.x) + B * (v0.y + v1.y);
    double F = C * (v0.x + v2.x) + D * (v0.y + v2.y);

    double G = 2.0 * (A * (v2.y - v1.y) - B * (v2.x - v1.x));

    double dx, dy;
    vertex center;
    double radius;

    // 共线点处理
    if (round(fabs(G)) == 0)
    {
        double minx = fmin(v0.x, fmin(v1.x, v2.x));
        double miny = fmin(v0.y, fmin(v1.y, v2.y));
        double maxx = fmax(v0.x, fmax(v1.x, v2.x));
        double maxy = fmax(v0.y, fmax(v1.y, v2.y));

        center.x = (minx + maxx) / 2;
        center.y = (miny + maxy) / 2;

        dx = center.x - minx;
        dy = center.y - miny;
    }
    else
    {
        double cx = (D * E - B * F) / G;
        double cy = (A * F - C * E) / G;

        center.x = cx;
        center.y = cy;

        dx = center.x - v0.x;
        dy = center.y - v0.y;
    }
    radius = sqrt(dx * dx + dy * dy);

    circum_circ circle = {center, radius};
    return circle;
}

// 创建三角形
triangle make_triangle(vertex v0, vertex v1, vertex v2)
{
    triangle t;
    t.v0 = v0;
    t.v1 = v1;
    t.v2 = v2;
    t.circumcirc = calc_circum_circ(v0, v1, v2);
    return t;
}

// 判断顶点是否在三角形的外接圆内
int in_circumcircle(triangle t, vertex v)
{
    double dx = t.circumcirc.c.x - v.x;
    double dy = t.circumcirc.c.y - v.y;
    // 用平方比较避免开根号，提升性能同时避免精度问题
    return (dx * dx + dy * dy) <= (t.circumcirc.r * t.circumcirc.r + 1e-9);
}

// 创建包含所有点的超级三角形
triangle make_super_triangle(vertex* vertices, int num_vertices)
{
    double minx = INFINITY, miny = INFINITY, maxx = -INFINITY, maxy = -INFINITY;
    for(int i=0; i<num_vertices; i++)
    {
        if(vertices[i].x < minx) minx = vertices[i].x;
        if(vertices[i].y < miny) miny = vertices[i].y;
        if(vertices[i].x > maxx) maxx = vertices[i].x;
        if(vertices[i].y > maxy) maxy = vertices[i].y;
    }

    double dx = (maxx - minx) * 10;
    double dy = (maxy - miny) * 10;

    vertex v0 = {minx - dx, miny - dy * 3};
    vertex v1 = {minx - dx, maxy + dy};
    vertex v2 = {maxx + dx * 3, maxy + dy};

    return make_triangle(v0, v1, v2);
}

// 添加顶点并更新三角剖分结果
triangle* add_vertex(vertex v, triangle* triangles, int* num_triangles)
{
    // 动态分配边数组，避免越界
    edge* edges = malloc(sizeof(edge) * (*num_triangles) * 3);
    if (!edges) return NULL;
    int num_edges = 0;

    // 移除所有外接圆包含当前顶点的三角形，并收集它们的边
    for(int i = 0; i < *num_triangles; i++)
    {
        if(in_circumcircle(triangles[i], v))
        {
            edges[num_edges++] = (edge){triangles[i].v0, triangles[i].v1};
            edges[num_edges++] = (edge){triangles[i].v1, triangles[i].v2};
            edges[num_edges++] = (edge){triangles[i].v2, triangles[i].v0};
            // 用最后一个三角形覆盖当前位置，减少数组移动
            triangles[i] = triangles[--(*num_triangles)];
            i--;
        }
    }

    // 筛选出唯一的边
    edge* unique_edges = malloc(sizeof(edge) * num_edges);
    if (!unique_edges) { free(edges); return NULL; }
    int num_unique_edges = 0;

    for(int i = 0; i < num_edges; i++)
    {
        int is_unique = 1;
        // 内层循环变量改用k，避免和外层i冲突
        for(int k = 0; k < num_edges; k++)
        {
            if(i != k && edge_equals(edges[i], edges[k]))
            {
                is_unique = 0;
                break;
            }
        }
        if(is_unique)
        {
            unique_edges[num_unique_edges++] = edges[i];
        }
    }

    // 创建新的三角形
    triangle* new_triangles = malloc(sizeof(triangle) * num_unique_edges);
    if (!new_triangles) { free(edges); free(unique_edges); return NULL; }

    for(int i = 0; i < num_unique_edges; i++)
    {
        new_triangles[i] = make_triangle(unique_edges[i].v0, unique_edges[i].v1, v);
    }

    // 扩展三角形数组
    triangle* result = realloc(triangles, sizeof(triangle) * (*num_triangles + num_unique_edges));
    if(result == NULL)
    {
        free(new_triangles);
        free(edges);
        free(unique_edges);
        return NULL;
    }

    triangles = result;
    for(int i=0; i < num_unique_edges; i++)
    {
        triangles[(*num_triangles)++] = new_triangles[i];
    }

    // 释放临时分配的内存
    free(new_triangles);
    free(edges);
    free(unique_edges);
    return triangles;
}

// 执行Delaunay三角剖分
triangle* triangulate(vertex* vertices, int num_vertices, int* num_triangles)
{
    triangle st = make_super_triangle(vertices, num_vertices);

    triangle* triangles = malloc(sizeof(triangle));
    if(triangles == NULL)
    {
        return NULL;
    }

    triangles[0] = st;
    *num_triangles = 1;

    // 逐个添加顶点
    for(int i=0; i < num_vertices; i++)
    {
        triangles = add_vertex(vertices[i], triangles, num_triangles);
        if(triangles == NULL)
        {
            free(triangles);
            return NULL;
        }
    }

    // 移除包含超级三角形顶点的三角形
    for(int i = 0; i < *num_triangles; i++)
    {
        if (vertex_equals(triangles[i].v0, st.v0) || vertex_equals(triangles[i].v0, st.v1) || vertex_equals(triangles[i].v0, st.v2) ||
            vertex_equals(triangles[i].v1, st.v0) || vertex_equals(triangles[i].v1, st.v1) || vertex_equals(triangles[i].v1, st.v2) ||
            vertex_equals(triangles[i].v2, st.v0) || vertex_equals(triangles[i].v2, st.v1) || vertex_equals(triangles[i].v2, st.v2)) {
            triangles[i] = triangles[--(*num_triangles)];
            i--;
        }
    }

    // 调整数组大小以匹配实际三角形数量
    triangle* result = realloc(triangles, sizeof(triangle) * (*num_triangles));
    if(result == NULL)
    {
        free(triangles);
        return NULL;
    }
    triangles = result;

    return triangles;
}

void create_fibonacci_sphere(int radius, int num_of_points, float points[][3], size_t arr_size)
{
    double golden_ratio = (1 + sqrt(5)) / 2;
    double angle_increment = 2 * M_PI * golden_ratio;

    // 生成斐波那契球面点
    for (int i = 0; i < num_of_points; ++i)
    {
        double t = (double)i / num_of_points;
        double angle1 = acos(1 - 2 * t);
        double angle2 = angle_increment * i;

        points[i][0] = radius * sin(angle1) * cos(angle2);
        points[i][1] = radius * sin(angle1) * sin(angle2);
        points[i][2] = radius * cos(angle1);
    }

    // 转换为2D顶点（这里选择x-z平面，若需要x-y平面可修改第二个参数）
    vertex point_list[num_of_points];
    for (int n = 0; n < num_of_points; n++)
    {
        point_list[n] = (vertex) {points[n][0], points[n][2]};
        // 若需要x-y平面，改为：point_list[n] = (vertex) {points[n][0], points[n][1]};
    }

    int num_triangles;
    // 显式转换size_t为int
    triangle* triangles = triangulate(point_list, (int)arr_size, &num_triangles);
    if (triangles)
    {
        // 打印三角剖分结果
        for (int i = 0; i < num_triangles; ++i)
        {
            printf("%f %f %f %f %f %f\n", triangles[i].v0.x, triangles[i].v0.y,
                   triangles[i].v1.x, triangles[i].v1.y,
                   triangles[i].v2.x, triangles[i].v2.y);
        }
        free(triangles);
    }

    printf("%d points generated on Fibonacci sphere with radius %d.\n", num_of_points, radius);
}

int main(void)
{
    int radius = 50;
    int num_of_points = 1000;
    float vertices[num_of_points][3];
    create_fibonacci_sphere(radius, num_of_points, vertices, sizeof(vertices)/sizeof(vertices[0]));

    printf("this will print if there are no errors.\n");
    return 0;
}

额外说明 #

• 修复了所有导致段错误的核心问题，包括数组越界、逻辑错误、语法错误等。
• 优化了浮点数比较逻辑，避免精度问题引发的判断错误。
• 增加了内存分配失败的检查，提升代码稳定性。
• 注释了顶点平面选择的位置，可根据需求切换x-y或x-z平面。

内容的提问来源于stack exchange，提问作者Chillzy