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

i searching online efficient method can intersect cartesian rectangular 3d grid close irregular 3d surface triangulated.

this surface represented set of vertices, v, , set of faces, f. cartesian rectangular grid stored as:

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

in figure below, single cell of cartesian grid shown. in addition, 2 triangles of surface schematically shown. intersection shown dotted red lines solid red circles intersection points particular cell. goal find points of intersection of surface edges of cells, can non-planar.

i implement in either matlab, c, or c++.

assuming have regular axis-aligned rectangular grid, each grid cell matching unit cube (and grid point (i,j,k) @ (i,j,k), i,j,k integers), suggest trying 3d variant of 2d triangle rasterization.

the basic idea draw triangle perimeter, every intersection between triangle , each plane perpendicular axis , intersecting axis @ integer coordinates.

you end line segments on grid cell faces, wherever triangle passes through grid cell. lines on faces of each grid cell form closed planar polygon. (however, you'll need connect line segments , orient polygon yourself.)

for finding out grid cells triangle passes through, simplified approach can used, , bitmap (one bit per grid cell). case 3d version of triangle rasterization.

the key observation if have line (x₀,y₀,z₀)-(x₁,y₁,z₁), can split segments @ integer coordinates x_i along x axis using

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

y_i = (1 - t_i) y₀ + t_i y₁

z_i = (1 - t_i) z₀ + t_i z₁

similarly along other axes, of course.

you'll need 3 passes, 1 along each axis. if sort vertices coordinates nondecreasing along axis, i.e. p₀ ≤ p₁ ≤ p₂, 1 endpoint @ integer coordinates intersecting line p₀p₂, , other endpoint intersects line p₀p₁ @ small coordinates, , line p₁p₂ @ large coordinates.

the intersection line between endpoints perpendicular 1 axis, still needs split segments not cross integer coordinates along other 2 dimensions. fortunately simple; maintain t_j , t_k along 2 dimensions t_i above, , step next integer coordinate has smaller t value; start @ 0, , end @ 1.

the edges of original triangle need drawn, split along 3 dimensions. again, straightforward, maintaining t each axis, , stepping along axis smallest value. have example code in c99 this, complicated case, below.

there quite few implementation details consider.

because each cell shares each face cell, , each edge 3 other edges, let's define following properties each cell (i,j,k), i,j,k integers identifying cell:

• x face: cell face @ x=i, perpendicular x axis
• y face: cell face @ y=j, perpendicular y axis
• z face: cell face @ z=k, perpendicular z axis
• x edge: edge (i,j,k) (i+1,j,k)
• y edge: edge (i,j,k) (i,j+1,k)
• z edge: edge (i,j,k) (i,j,k+1)

the other 3 faces cell (i,j,k) are

• the x face @ (i+1,j,k)
• the y face @ (i,j+1,k)
• the z face @ (i,j,k+1)

similarly, each edge edge 3 other cells. x edge of cell (i,j,k) edge grid cells (i,j+1,k), (i,j,k+1), , (i,j+1,k+1). y edge of cell (i,j,k) edge grid cells (i+1,j,k), (i,j,k+1), , (i+1,j,k+1). z edge of cell (i,j,k) edge grid cells (i+1,j,k), (i,j+1,k), , (i+1,j+1,k).

here image might help.

(ignore fact it's left-handed; thought it'd easier label way.)

this means if have line segment on specific grid cell face, line segment shared between 2 grid cells sharing face. similarly, if line segment endpoint on grid cell edge, there 4 different grid cell faces other line segment endpoint on.

to clarify this, example code below prints not coordinates, grid cell , face/edge/vertex line segment endpoint 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 inside grid cell */     x_face = 1,    /* point @ integer x coordinate (on yz face) */     y_face = 2,    /* point @ integer y coordinate (on xz face) */     z_edge = 3,    /* point @ integet x , y coordinates (on z edge) */     z_face = 4,    /* point @ integer z coordinate (on xy face) */     y_edge = 5,    /* point @ integer x , z coordinates (on y edge) */     x_edge = 6,    /* point @ integer y , z coordinates (on x edge) */     vertex = 7,    /* point @ integer coordinates (at 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 */             @ = next_at.x;             curr_vec.x = round( (1.0 - at) * p0.x + @ * p1.x );             curr_vec.y =        (1.0 - at) * p0.y + @ * p1.y;             curr_vec.z =        (1.0 - at) * p0.z + @ * p1.z;         } else         if (next_at.y < 1.0 && next_at.y < next_at.x && next_at.y <= next_at.z) {             /* xz plane */             @ = next_at.y;             curr_vec.x =        (1.0 - at) * p0.x + @ * p1.x;             curr_vec.y = round( (1.0 - at) * p0.y + @ * p1.y );             curr_vec.z =        (1.0 - at) * p0.z + @ * p1.z;         } else         if (next_at.z < 1.0 && next_at.z < next_at.x && next_at.z < next_at.y) {             /* xy plane */             @ = next_at.z;             curr_vec.x =        (1.0 - at) * p0.x + @ * p1.x;             curr_vec.y =        (1.0 - at) * p0.y + @ * p1.y;             curr_vec.z = round( (1.0 - at) * p0.z + @ * p1.z );         } else {             @ = 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 2 command-line parameters, 3d endpoints line segmented. if compile above example, 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 line (0.5, 0.25, 3.50) - (3.5, 4.0, 0.50) gets split 7 segments; particular line passing through 7 grid cells.

for rasterization case -- when interested in grid cells surface triangles pass through --, not need store line segment points, compute them all. when point @ vertex or inside grid cell, mark bit corresponding grid cell. when point @ face, set bit 2 grid cells share face. when line segment endpoint @ edge, set bit 4 grid cells share edge.

questions?