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

Question

我正在在线搜索一种有效的方法，该方法可以将笛卡尔矩形 3D 网格与一个封闭的不规则 3D 表面相交，该表面是三角形的。

这个表面由一组顶点V和一组面F表示。笛卡尔矩形网格存储为：

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 变体。

基本思想是绘制三角形周长，然后是三角形与垂直于轴的每个平面之间的每个交点，并在整数坐标处与该轴相交。

在三角形通过网格单元的任何地方，您最终都会在网格单元面上看到线段。每个网格单元面上的线形成一个封闭的平面多边形。 （但是，您需要自己连接线段并确定多边形的方向。）

为了仅找出三角形通过的网格单元，可以使用简化的方法和位图（每个网格单元一个位）。这种情况本质上只是三角形光栅化的 3D 版本。

关键的观察是，如果你有一条线 (X₀,Y₀,Z₀)-(X₁,Y₁ ,Z₁)，您可以将其拆分为整数坐标 x_i 沿 x 轴使用

t_i = (x_i - X₀) / ( X₁ - X₀)

y_i = (1 - t_i) Y₀ + t_iY₁

z_i = (1 - t_i) Z₀ + t_iZ₁

当然，沿其他轴也是如此。

您需要进行三遍，沿每个轴进行一次。如果您对顶点进行排序，使坐标沿该轴不递减，即 p₀ ≤ p₁ ≤ p₂，一个端点在整数坐标与p₀p2，另一个端点在小坐标处与线p₀p₁相交，和线 p₁p₂ 在大坐标。

这些端点之间的相交线垂直于一个轴，但仍需要将其分割成不与其他两个维度上的整数坐标相交的线段。幸运的是，这很简单；就像 ti，并步进到具有较小 t 值的下一个整数坐标；从 0 开始，到 1 结束。

原来三角形的边也需要画出来，沿三个维度分割即可。同样，这很简单，通过保持每个轴的 t 并沿着具有最小值的轴步进。我有 C99 中的示例代码，这是最复杂的情​​况，如下所示。

有很多实现细节需要考虑。

因为每个单元格与另一个单元格共享每个面，并且每个边与其他三个边共享，所以让我们为每个单元格定义以下属性（i,j,k)，其中 i,j,k 是标识单元格的整数：

• X 面：x=i 处的单元面，垂直于 x 轴
• Y 面：y=j 处的单元面，垂直于 y 轴
• Z 面：z=k 处的单元面，垂直于 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)

单元格的其他三个面（i,j,k）是

• X 面在 (i+1,j,k)
• Y 面在 (i,j+1,k)
• Z 面在 (i,j,k+1)

同样，每条边都是其他三个单元格的边。单元格(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)。

这是一张可能有帮助的图片。

（忽略它是左撇子的事实；我只是认为这样标记会更容易。）

这意味着如果您在特定网格单元face 上有一条线段，则该线段将在共享该face 的两个网格单元之间共享。类似地，如果线段端点位于网格单元边上，则另一个线段端点可能位于四个不同的网格单元面上。

为了澄清这一点，下面的示例代码不仅打印坐标，还打印网格单元和线段端点所在的面/边/顶点。

#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;
}

程序采用两个命令行参数，即要分割的线的 3D 端点。如果你编译上面说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) 被分成七段；这条特殊的线正好穿过七个网格单元。

对于光栅化情况——当您只对曲面三角形通过哪些网格单元感兴趣时——您不需要存储线段点，只需计算它们。当一个点位于顶点或网格单元内时，标记对应于该网格单元的位。当一个点在一个面上时，为共享该面的两个网格单元设置该位。当线段端点位于边缘时，为共享该边缘的四个网格单元设置该位。

问题？

Answer 2

将问题分解成更小的步骤。

Finding the intersection of a line segment with a triangle is easy.

实现之后，只需执行一个嵌套循环，检查网格中的每个线组合与曲面三角形之间的交点。

Answer 3

按照here 的说明在每个曲面三角形边上使用线平面相交。您将使用六个平面，每个网格单元面一个。