Search code examples
cmatlabgeometryintersectiontriangulation

Finding the intersection of a closed irregular triangulated 3D surface with a Cartesian rectangular 3D grid

I am searching online for an efficient method that can intersect a Cartesian rectangular 3D grid with a close irregular 3D surface which is triangulated.

This surface is represented is a set of vertices, V, and a set of faces, F. The Cartesian rectangular grid is stored as:

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

In the figure below, a single cell of the Cartesian grid is shown. In addition, two triangles of the surface is schematically shown. This intersection is shown by a dotted red lines with the solid red circles the intersection points with this particular cell. My goal is to find the points of intersection of the surface with the edges of the cells, which can be non-planar.

I will implement in either MATLAB, C, or C++.

enter image description here

Solution

  • Assuming we have a regular axis-aligned rectangular grid, with each grid cell matching the unit cube (and thus grid point (i,j,k) is at (i,j,k), with i,j,k integers), I would suggest trying a 3D variant of 2D triangle rasterization.

    The basic idea is to draw the triangle perimeter, then every intersection between the triangle and each plane perpendicular to an axis and intersecting that axis at integer coordinates.

    You end up with line segments on grid cell faces, wherever the triangle passes through a grid cell. The lines on the faces of each grid cell form a closed planar polygon. (However, you'll need to connect the line segments and orient the polygon yourself.)

    For finding out only the grid cells the triangle passes through, a simplified approach can be used, and a bitmap (one bit per grid cell). This case is essentially just a 3D version of triangle rasterization.

    The key observation is that if you have a line (X0,Y0,Z0)-(X1,Y1,Z1), you can split it into segments at integer coordinates xi along the x axis using

    ti = (xi - X0) / (X1 - X0)

    yi = (1 - ti) Y0 + ti Y1

    zi = (1 - ti) Z0 + ti Z1

    Similarly along the other axes, of course.

    You'll need to do three passes, one along each axis. If you sort the vertices so that the coordinates are nondecreasing along that axis, i.e. p0 ≤ p1 ≤ p2, one endpoint is at integer coordinates intersecting line p0p2, and the other endpoint intersects line p0p1 at small coordinates, and line p1p2 at large coordinates.

    The intersection line between those endpoints is perpendicular to one axis, but it still needs to be split into segments that do not cross integer coordinates along the other two dimensions. This is fortunately simple; just maintain tj and tk along those two dimensions just like ti above, and step to the next integer coordinate that has the smaller t value; start at 0, and end at 1.

    The edges of the original triangle also need to be drawn, just split along all three dimensions. Again, this is straightforward, by maintaining the t for each axis, and stepping along the axis with the smallest value. I have example code in C99 for this, the most complicated case, below.

    There are quite a few implementation details to consider.

    Because each cell shares each face with another cell, and each edge with three other edges, let's define the following properties for each cell (i,j,k), where i,j,k are integers identifying the cell:

    • X face: The cell face at x=i, perpendicular to the x axis
    • Y face: The cell face at y=j, perpendicular to the y axis
    • Z face: The cell face at z=k, perpendicular to the z axis
    • X edge: The edge from (i,j,k) to (i+1,j,k)
    • Y edge: The edge from (i,j,k) to (i,j+1,k)
    • Z edge: The edge from (i,j,k) to (i,j,k+1)

    The other three faces for cell (i,j,k) are

    • the X face at (i+1,j,k)
    • the Y face at (i,j+1,k)
    • the Z face at (i,j,k+1)

    Similarly, each edge is an edge for three other cells. The X edge of cell (i,j,k) is also an edge for grid cells (i,j+1,k), (i,j,k+1), and (i,j+1,k+1). The Y edge of cell (i,j,k) is also an edge for grid cells (i+1,j,k), (i,j,k+1), and (i+1,j,k+1). The Z edge of cell (i,j,k) is also an edge for grid cells (i+1,j,k), (i,j+1,k), and (i+1,j+1,k).

    Here is an image that might help. Illustration of a regular rectangular axis-aligned unit grid cell

    (Ignore the fact that it's left-handed; I just thought it'd be easier to label this way.)

    This means that if you have a line segment on a specific grid cell face, the line segment is shared between the two grid cells sharing that face. Similarly, if a line segment endpoint is on a grid cell edge, there are four different grid cell faces the other line segment endpoint could be on.

    To clarify this, my example code below prints not only the coordinates, but the grid cell and the face/edge/vertex the line segment endpoint is on.

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

    The program takes two command-line parameters, the 3D endpoints for the line to be segmented. If you compile the above to say example, then running

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

    outputs

    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

    which shows that line (0.5, 0.25, 3.50) - (3.5, 4.0, 0.50) gets split into seven segments; this particular line passing through exactly seven grid cells.

    For the rasterization case -- when you are only interested in which grid cells the surface triangles pass through --, you do not need to store the line segment points, only compute them all. When a point is at a vertex or inside a grid cell, mark the bit corresponding to that grid cell. When a point is at a face, set the bit for the two grid cells that share that face. When a line segment endpoint is at an edge, set the bit for the four grid cells that share that edge.

    Questions?