c - 查找封闭的不规则三角3D曲面与笛卡尔矩形3D网格的交点

Question

我正在网上寻找一种有效的方法，可以交叉笛卡尔矩形三维网格和一个封闭的不规则三维表面，这是三角化。
该曲面表示为一组顶点和一组面。笛卡尔矩形网格存储为：

x_0, x_1, ..., x_(ni-1)
y_0, y_1, ..., y_(nj-1)
z_0, z_1, ..., z_(nk-1)

在下图中，显示了笛卡尔网格的单个单元格。此外，两个三角形的表面是示意图所示。此交叉点用红色虚线表示，实心红色圆圈表示与此特定单元格的交叉点。我的目标是找到曲面与单元边缘的交点，这些交点可以是非平面的。
我将在Matlab、C或C++中实现。

Answer 1

假设我们有一个规则的轴对齐的矩形网格，每个网格单元与单位立方体匹配（因此网格点（i，j，k）是在（i，j，k），i，j，k整数），我建议尝试一个2D三角形光栅化的3D变体。
其基本思想是绘制三角形周长，然后每个三角形与每个平面之间的交点垂直于一个轴，并在整数坐标下相交该轴线。
最后在网格单元面上显示线段，三角形通过网格单元的任何位置。每个网格单元的面上的线形成一个闭合的平面多边形。（但是，您需要连接线段并自行确定多边形的方向。）
为了只找出三角形通过的网格单元，可以使用一种简化的方法，以及一个位图（每个网格单元一位）。这个例子本质上只是三角光栅化的三维版本。
关键的观察是，如果你有一行（X0，Y0，Z0）-（X1，Y1，Z1），你可以在X轴上使用整数坐标Xi将其分割成段。
Ti＝（Xi-x0）/（x1-x0）
yi=（1-ti）Y0+ti Y1
zi=（1-ti）Z0+ti Z1
当然，沿着其他轴也是一样。
你需要做三遍，每个轴一个。如果对顶点进行排序以使坐标沿该轴不减，即P0＜P1＜P2，则一个端点处于整数坐标相交的线P0P2，而另一端点在小坐标处相交线P0P1，而在大坐标处相交P1P2线。
这些端点之间的交线垂直于一个轴，但它仍然需要被分割成不沿其他两个维度交叉整数坐标的段。幸运的是，这很简单；只需沿着这两个维度维护tj和tk，就像上面的ti一样，然后转到具有较小t值的下一个整数坐标；从0开始，到1结束。
原始三角形的边也需要绘制，只需沿所有三维分割。同样，这是直截了当的，通过保持每个轴的T，并沿着具有最小值的轴步进。下面是C99中最复杂的例子代码。
有相当多的实施细节需要考虑。
因为每个单元格与另一个单元格共享每个面，每个边与三个其他边共享，所以让我们为每个单元格（i，j，k）定义以下属性，其中i，j，k是标识该单元格的整数：
X面：x= i的细胞面，垂直于x轴
Y面：Y轴上的细胞面，垂直于Y轴
Z面：Z轴上的细胞面，垂直于Z轴
X边：从（i，j，k）到（i+1，j，k）的边
Y边：从（i，j，k）到（i，j+1，k）的边
Z边：从（i，j，k）到（i，j，k+1）的边
cell（i，j，k）的另外三个面是
（i+1，j，k）处的X面
（i，j+1，k）处的Y面
（i，j，k+1）处的Z面
类似地，每个边都是其他三个单元格的边。单元格（i，j，k）的X边也是网格单元格（i，j+1，k），（i，j，k+1）和（i，j+1，k+1）的边。单元格（i，j，k）的Y边也是网格单元格（i+1，j，k），（i，j，k+1）和（i+1，j，k+1）的边。单元格（i，j，k）的Z边也是网格单元格（i+1，j，k），（i，j+1，k）和（i+1，j+1，k）的边。
这是一张可能有用的图片。
（忽略它是左手的事实；我只是认为这样做更容易标记。）
这意味着，如果在特定网格单元面上有线段，则线段将在共享该面的两个网格单元之间共享。类似地，如果线段端点位于网格单元边上，则有四个不同的网格单元面，另一个线段端点可能位于该面上。
为了澄清这一点，下面的示例代码不仅打印坐标，还打印网格单元和线段端点所在的面/边/顶点。

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

typedef struct {
    double x;
    double y;
    double z;
} vector;

typedef struct {
    long   x;
    long   y;
    long   z;
} gridpos;

typedef enum {
    INSIDE = 0,    /* Point is inside the grid cell */
    X_FACE = 1,    /* Point is at integer X coordinate (on the YZ face) */
    Y_FACE = 2,    /* Point is at integer Y coordinate (on the XZ face) */
    Z_EDGE = 3,    /* Point is at integet X and Y coordinates (on the Z edge) */
    Z_FACE = 4,    /* Point is at integer Z coordinate (on the XY face) */
    Y_EDGE = 5,    /* Point is at integer X and Z coordinates (on the Y edge) */
    X_EDGE = 6,    /* Point is at integer Y and Z coordinates (on the X edge) */
    VERTEX = 7,    /* Point is at integer coordinates (at the grid point) */
} cellpos;

static inline cellpos cellpos_of(const vector v)
{
    return (v.x == floor(v.x))
         + (v.y == floor(v.y)) * 2
         + (v.z == floor(v.z)) * 4;
}

static const char *const face_name[8] = {
    "inside",
    "x-face",
    "y-face",
    "z-edge",
    "z-face",
    "y-edge",
    "x-edge",
    "vertex",
};

static int line_segments(const vector p0, const vector p1,
                         int (*segment)(void *custom,
                                        const gridpos src_cell, const cellpos src_face, const vector src_vec,
                                        const gridpos dst_cell, const cellpos dst_face, const vector dst_vec),
                         void *const custom)
{
    const vector  range = { p1.x - p0.x, p1.y - p0.y, p1.z - p0.z };
    const gridpos step = { (range.x < 0.0) ? -1L : (range.x > 0.0) ? +1L : 0L,
                           (range.y < 0.0) ? -1L : (range.y > 0.0) ? +1L : 0L,
                           (range.z < 0.0) ? -1L : (range.z > 0.0) ? +1L : 0L };
    const gridpos end = { floor(p1.x), floor(p1.y), floor(p1.z) };
    gridpos       prev_cell, curr_cell = { floor(p0.x), floor(p0.y), floor(p0.z) };
    vector        prev_vec,  curr_vec = p0;
    vector        curr_at = { 0.0, 0.0, 0.0 };
    vector        next_at = { (range.x != 0.0 && curr_cell.x != end.x) ? ((double)(curr_cell.x + step.x) - p0.x) / range.x : 1.0,
                              (range.y != 0.0 && curr_cell.y != end.y) ? ((double)(curr_cell.y + step.y) - p0.y) / range.y : 1.0,
                              (range.z != 0.0 && curr_cell.z != end.z) ? ((double)(curr_cell.z + step.z) - p0.z) / range.z : 1.0};
    cellpos       prev_face, curr_face;
    double        at;
    int           retval;

    curr_face = cellpos_of(p0);

    while (curr_at.x < 1.0 || curr_at.y < 1.0 || curr_at.z < 1.0) {
        prev_cell = curr_cell;
        prev_face = curr_face;
        prev_vec  = curr_vec;

        if (next_at.x < 1.0 && next_at.x <= next_at.y && next_at.x <= next_at.z) {
            /* YZ plane */
            at = next_at.x;
            curr_vec.x = round( (1.0 - at) * p0.x + at * p1.x );
            curr_vec.y =        (1.0 - at) * p0.y + at * p1.y;
            curr_vec.z =        (1.0 - at) * p0.z + at * p1.z;
        } else
        if (next_at.y < 1.0 && next_at.y < next_at.x && next_at.y <= next_at.z) {
            /* XZ plane */
            at = next_at.y;
            curr_vec.x =        (1.0 - at) * p0.x + at * p1.x;
            curr_vec.y = round( (1.0 - at) * p0.y + at * p1.y );
            curr_vec.z =        (1.0 - at) * p0.z + at * p1.z;
        } else
        if (next_at.z < 1.0 && next_at.z < next_at.x && next_at.z < next_at.y) {
            /* XY plane */
            at = next_at.z;
            curr_vec.x =        (1.0 - at) * p0.x + at * p1.x;
            curr_vec.y =        (1.0 - at) * p0.y + at * p1.y;
            curr_vec.z = round( (1.0 - at) * p0.z + at * p1.z );
        } else {
            at = 1.0;
            curr_vec = p1;
        }

        curr_face = cellpos_of(curr_vec);

        curr_cell.x = floor(curr_vec.x);
        curr_cell.y = floor(curr_vec.y);
        curr_cell.z = floor(curr_vec.z);

        retval = segment(custom,
                         prev_cell, prev_face, prev_vec,
                         curr_cell, curr_face, curr_vec);
        if (retval)
            return retval;

        if (at < 1.0) {
            curr_at = next_at;
            if (at >= next_at.x) {
                /* recalc next_at.x */
                if (curr_cell.x != end.x) {
                    next_at.x = ((double)(curr_cell.x + step.x) - p0.x) / range.x;
                    if (next_at.x > 1.0)
                        next_at.x = 1.0;
                } else
                    next_at.x = 1.0;
            }
            if (at >= next_at.y) {
                /* reclac next_at.y */
                if (curr_cell.y != end.y) {
                    next_at.y = ((double)(curr_cell.y + step.y) - p0.y) / range.y;
                    if (next_at.y > 1.0)
                        next_at.y = 1.0;
                } else
                    next_at.y = 1.0;
            }
            if (at >= next_at.z) {
                /* recalc next_at.z */
                if (curr_cell.z != end.z) {
                    next_at.z = ((double)(curr_cell.z + step.z) - p0.z) / range.z;
                    if (next_at.z > 1.0)
                        next_at.z = 1.0;
                } else
                    next_at.z = 1.0;
            }
        } else {
            curr_at.x = curr_at.y = curr_at.z = 1.0;
            next_at.x = next_at.y = next_at.z = 1.0;
        }
    }

    return 0;
}

int print_segment(void *outstream,
                  const gridpos src_cell, const cellpos src_face, const vector src_vec,
                  const gridpos dst_cell, const cellpos dst_face, const vector dst_vec)
{
    FILE *const out = outstream ? outstream : stdout;

    fprintf(out, "%.6f %.6f %.6f   %.6f %.6f %.6f   %s %ld %ld %ld   %s %ld %ld %ld\n",
                 src_vec.x, src_vec.y, src_vec.z,
                 dst_vec.x, dst_vec.y, dst_vec.z,
                 face_name[src_face], src_cell.x, src_cell.y, src_cell.z,
                 face_name[dst_face], dst_cell.x, dst_cell.y, dst_cell.z);

    fflush(out);
    return 0;
}

static int parse_vector(const char *s, vector *const v)
{
    double x, y, z;
    char   c;

    if (!s)
        return EINVAL;

    if (sscanf(s, " %lf %*[,/:;] %lf %*[,/:;] %lf %c", &x, &y, &z, &c) == 3) {
        if (v) {
            v->x = x;
            v->y = y;
            v->z = z;
        }
        return 0;
    }

    return ENOENT;
}


int main(int argc, char *argv[])
{
    vector start, end;

    if (argc != 3 || !strcmp(argv[1], "-h") || !strcmp(argv[1], "--help")) {
        fprintf(stderr, "\n");
        fprintf(stderr, "Usage: %s [ -h | --help ]\n", argv[0]);
        fprintf(stderr, "       %s x0:y0:z0 x1:y1:z1\n", argv[0]);
        fprintf(stderr, "\n");
        return EXIT_FAILURE;
    }

    if (parse_vector(argv[1], &start)) {
        fprintf(stderr, "%s: Invalid start point.\n", argv[1]);
        return EXIT_FAILURE;
    }
    if (parse_vector(argv[2], &end)) {
        fprintf(stderr, "%s: Invalid end point.\n", argv[2]);
        return EXIT_FAILURE;
    }

    if (line_segments(start, end, print_segment, stdout))
        return EXIT_FAILURE;

    return EXIT_SUCCESS;
}

程序采用两个命令行参数，即要分段的线的三维端点。如果将上面的代码编译为 example，则运行

./example 0.5/0.25/3.50 3.5/4.0/0.50

输出

0.500000 0.250000 3.500000   1.000000 0.875000 3.000000   inside 0 0 3   x-face 1 0 3
1.000000 0.875000 3.000000   1.100000 1.000000 2.900000   x-face 1 0 3   y-face 1 1 2
1.100000 1.000000 2.900000   1.900000 2.000000 2.100000   y-face 1 1 2   y-face 1 2 2
1.900000 2.000000 2.100000   2.000000 2.125000 2.000000   y-face 1 2 2   y-edge 2 2 2
2.000000 2.125000 2.000000   2.700000 3.000000 1.300000   y-edge 2 2 2   y-face 2 3 1
2.700000 3.000000 1.300000   3.000000 3.375000 1.000000   y-face 2 3 1   y-edge 3 3 1
3.000000 3.375000 1.000000   3.500000 4.000000 0.500000   y-edge 3 3 1   y-face 3 4 0

这表明线（0.5，0.25，3.50）-（3.5，4.0，0.50）被分成七段；这条特殊的线正好穿过七个网格单元。
对于光栅化情况（当您只对曲面三角形通过的网格单元感兴趣时），不需要存储线段点，只需计算所有线段点。当点位于顶点或网格单元内时，标记与该网格单元相对应的位。当一个点位于面上时，为共享该面的两个网格单元设置位。当线段端点位于边上时，请为共享该边的四个网格单元设置位。
问题？