查找封闭的不规则三角化 3D 曲面与笛卡尔矩形 3D 网格的交点
Finding the intersection of a closed irregular triangulated 3D surface with a Cartesian rectangular 3D grid
我正在在线搜索一种有效的方法，该方法可以将笛卡尔矩形 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++ 中实现。
将问题分解成更小的步骤。
Finding the intersection of a line segment with a triangle is easy.
一旦你实现了这个，只需执行一个嵌套循环来检查网格中的线与表面三角形的每个组合之间的交集。
假设我们有一个规则的轴对齐矩形网格，每个网格单元匹配单位立方体（因此网格点（i，j,k) 在 (i,j,k), i,j,k 整数), 我建议尝试2D 三角形光栅化的 3D 变体。
基本思路是画出三角形的周长，然后画出三角形与垂直于轴并与该轴相交于整数坐标的每个平面的每个交点。
无论三角形穿过网格单元的什么地方，您最终都会在网格单元面上得到线段。每个网格单元面上的线形成一个封闭的平面多边形。 （但是，您需要自己连接线段并确定多边形的方向。）
为了仅找出三角形穿过的网格单元，可以使用简化的方法和位图（每个网格单元一位）。这种情况本质上只是三角形光栅化的 3D 版本。
关键的观察是，如果你有一行 (X0,Y 0,Z0)-(X1,Y1,Z1 ), 你可以沿着 xxi 将它分割成整数坐标段]轴使用
ti = (xi - X0) / (X1 - X0)
yi = (1 - ti) Y0 + ti Y1
zi = (1 - ti) Z0 + ti Z1
当然，其他轴也类似。
您需要进行三遍，每个轴各一遍。如果对顶点进行排序，使坐标沿该轴不递减，即 p0 ≤ p 1≤p2,一端点在整数坐标交线p0p2，另一个端点与线p相交0p1 在小坐标，线 p 1p2 大坐标.
这些端点之间的交线垂直于一个轴，但仍需要将其拆分为不沿其他两个维度与整数坐标相交的线段。幸运的是，这很简单；沿着这两个保持 tj 和 tk尺寸就像上面的 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 edge: 来自(i,j,[=60=的边]k) 到 (i+1,j,k)
- Y edge: 来自(i,j,[=60=的边]k) 到 (i,j+1,k)
- Z edge: 来自(i,j,[=60=的边]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).
这是一张可能有帮助的图片。
（忽略它是左撇子的事实；我只是觉得这样标记会更容易。）
这意味着如果您在特定网格单元面上有一条线段，则该线段在共享该面[的两个网格单元之间共享=218=]。类似地，如果线段端点在网格单元 edge 上，则有四个不同的网格单元 faces 另一个线段端点可能在。
为了澄清这一点，我下面的示例代码不仅打印坐标，还打印网格单元和 face/edge/vertex 线段端点打开。
#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) 被分成七段；这条特定的线恰好穿过七个网格单元。
对于光栅化的情况——当你只对表面三角形穿过哪些网格单元感兴趣时——你不需要存储线段点，只需要计算它们。当点位于顶点或网格单元内时，标记对应于该网格单元的位。当一个点在一个面上时，为共享该面的两个网格单元设置位。当线段端点位于一条边上时，为共享该边的四个网格单元设置位。
有问题吗？
按照 here 在每个曲面三角形边上的说明使用线-平面相交。您将使用六个平面，每个网格单元面一个。
我正在在线搜索一种有效的方法，该方法可以将笛卡尔矩形 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++ 中实现。
将问题分解成更小的步骤。
Finding the intersection of a line segment with a triangle is easy.
一旦你实现了这个，只需执行一个嵌套循环来检查网格中的线与表面三角形的每个组合之间的交集。
假设我们有一个规则的轴对齐矩形网格，每个网格单元匹配单位立方体（因此网格点（i，j,k) 在 (i,j,k), i,j,k 整数), 我建议尝试2D 三角形光栅化的 3D 变体。
基本思路是画出三角形的周长，然后画出三角形与垂直于轴并与该轴相交于整数坐标的每个平面的每个交点。
无论三角形穿过网格单元的什么地方，您最终都会在网格单元面上得到线段。每个网格单元面上的线形成一个封闭的平面多边形。 （但是，您需要自己连接线段并确定多边形的方向。）
为了仅找出三角形穿过的网格单元，可以使用简化的方法和位图（每个网格单元一位）。这种情况本质上只是三角形光栅化的 3D 版本。
关键的观察是，如果你有一行 (X0,Y 0,Z0)-(X1,Y1,Z1 ), 你可以沿着 xxi 将它分割成整数坐标段]轴使用
ti = (xi - X0) / (X1 - X0)
yi = (1 - ti) Y0 + ti Y1
zi = (1 - ti) Z0 + ti Z1
当然，其他轴也类似。
您需要进行三遍，每个轴各一遍。如果对顶点进行排序，使坐标沿该轴不递减，即 p0 ≤ p 1≤p2,一端点在整数坐标交线p0p2，另一个端点与线p相交0p1 在小坐标，线 p 1p2 大坐标.
这些端点之间的交线垂直于一个轴，但仍需要将其拆分为不沿其他两个维度与整数坐标相交的线段。幸运的是，这很简单；沿着这两个保持 tj 和 tk尺寸就像上面的 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 edge: 来自(i,j,[=60=的边]k) 到 (i+1,j,k)
- Y edge: 来自(i,j,[=60=的边]k) 到 (i,j+1,k)
- Z edge: 来自(i,j,[=60=的边]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).
这是一张可能有帮助的图片。
（忽略它是左撇子的事实；我只是觉得这样标记会更容易。）
这意味着如果您在特定网格单元面上有一条线段，则该线段在共享该面[的两个网格单元之间共享=218=]。类似地，如果线段端点在网格单元 edge 上，则有四个不同的网格单元 faces 另一个线段端点可能在。
为了澄清这一点，我下面的示例代码不仅打印坐标，还打印网格单元和 face/edge/vertex 线段端点打开。
#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) 被分成七段；这条特定的线恰好穿过七个网格单元。
对于光栅化的情况——当你只对表面三角形穿过哪些网格单元感兴趣时——你不需要存储线段点，只需要计算它们。当点位于顶点或网格单元内时，标记对应于该网格单元的位。当一个点在一个面上时，为共享该面的两个网格单元设置位。当线段端点位于一条边上时，为共享该边的四个网格单元设置位。
有问题吗？
按照 here 在每个曲面三角形边上的说明使用线-平面相交。您将使用六个平面，每个网格单元面一个。