3D triangle rasterization into voxel grid

Prologue:

This is Q&A to cover 3D triangle rasterization into voxel grid I was asked for in order to solve a different problem related to material erosion/removal during manufactoring process simulation. The main Idea behind this question is how to port scanline based 2D triangle rasterization like this into 3D voxels.

So the question is How to efficiently rasterize 3D triangle into uniform voxel grid represented as 3D array with Voxels represented as floating point value (density or color if you like).

Using C++ without any additional libs for the rasterization itself.

Solution

So first updated algorithm:

definitions

Lets the triangle be given by its 3 int 3D vertexes and one float color:
```
void triangle(int x0,int y0,int z0,int x1,int y1,int z1,int x2,int y2,int z2,float c);
```
The target voxel grid is 3D array like this:
```
float map[vxs][vys][vzs];
```
where vxs,vys,vzs is the resolution of grid. Scanline algorithm requires left and right buffers however in 3D we can not use single axis for it so instead I chose the major axis (axis which has biggest side of AABB BBOX of input triangle) so buffer size vsz should be max of the 3 resolutions. As its just 1D array its not that memory hungry so I choose this instead for convenience:
```
const int vsz=vxs|vys|vzs;
int lx[vsz],ly[vsz],lz[vsz]; // left buffers for filling
int rx[vsz],ry[vsz],rz[vsz]; // right buffers for filling
```

compute AABB and chose major axis

 int X0,Y0,Z0,X1,Y1,Z1,dx,dy,dz;
 // BBOX
            X0=x0;            X1=x0;
 if (X0>x1) X0=x1; if (X1<x1) X1=x1;
 if (X0>x2) X0=x2; if (X1<x2) X1=x2;
            Y0=y0;            Y1=y0;
 if (Y0>y1) Y0=y1; if (Y1<y1) Y1=y1;
 if (Y0>y2) Y0=y2; if (Y1<y2) Y1=y2;
            Z0=z0;            Z1=z0;
 if (Z0>z1) Z0=z1; if (Z1<z1) Z1=z1;
 if (Z0>z2) Z0=z2; if (Z1<z2) Z1=z2;
 dx=X1-X0;
 dy=Y1-Y0;
 dz=Z1-Z0;

now select major axis from the value of dx,dy,dz (biggest one). So for example if dx is biggest then from now on the left right buffers index will be x coordinate.

rasterize triangle edges into left/right buffers

now to chose if the edge should go to the left or right buffer we can use the sign of the major coordiante difference. So chose which one it is first and then order the edge vertexes by major coordinate (to avoid holes on joined triangles).

For the rasterization of edge (line) simply use integer version of DDA algortihm which is easily portable to higher dimensions and its also faster than Bresenham on modern CPUs.

Now instead of rendering voxels (x,y,z) of the line into map we simply render it into left/right buffers. So if the major axis is x and target buffer is left it would be:
```
ly[x]=y;
lz[x]=z;
```
As there are 3 possible major axises we need 3 versions of this function... Also there are 3 cases of the target buffer (left, right, both) but for that you might use pointers so non need to have each code twice...

Also in case the edges might be out of voxel grid range clipping is needed to add. I used simple if condition inside DDA iteration however to improve speed its better to clip the line before that...
fill in the triangle using lines

simply go through the left/right buffers and join the vertexes at the same index they represent with line. The line itself might be simplified as we knwo the major axis will not change in it however for simplicity I used normal 3D DDA line (the same as #3 just directly renders into map instead of left/right buffers).

Putting all together:

// voxel grid resolution
const int vxs=100;
const int vys=100;
const int vzs=100;
// helper buffers size should be max(vxs,vys,vzs)
const int vsz=vxs|vys|vzs;
float map[vxs][vys][vzs];       // voxel space map
int lx[vsz],ly[vsz],lz[vsz];    // left buffers for filling
int rx[vsz],ry[vsz],rz[vsz];    // right buffers for filling

void line(int x0,int y0,int z0,int x1,int y1,int z1,float c) // rencer line
    {
    int i,n,cx,cy,cz,sx,sy,sz;
    // line DDA parameters
    x1-=x0; sx=0; if (x1>0) sx=+1; if (x1<0) { sx=-1; x1=-x1; } if (x1) x1++;           n=x1;
    y1-=y0; sy=0; if (y1>0) sy=+1; if (y1<0) { sy=-1; y1=-y1; } if (y1) y1++; if (n<y1) n=y1;
    z1-=z0; sz=0; if (z1>0) sz=+1; if (z1<0) { sz=-1; z1=-z1; } if (z1) z1++; if (n<z1) n=z1;
    // single pixel (not a line)
    if (!n)
        {
        if ((x0>=0)&&(x0<vxs)&&(y0>=0)&&(y0<vys)&&(z0>=0)&&(z0<vzs)) map[x0][y0][z0]=c;
        return;
        }
    // ND DDA algo i is parameter
    for (cx=cy=cz=n,i=0;i<n;i++)
        {
        if ((x0>=0)&&(x0<vxs)&&(y0>=0)&&(y0<vys)&&(z0>=0)&&(z0<vzs)) map[x0][y0][z0]=c;
        cx-=x1; if (cx<=0){ cx+=n; x0+=sx; }
        cy-=y1; if (cy<=0){ cy+=n; y0+=sy; }
        cz-=z1; if (cz<=0){ cz+=n; z0+=sz; }
        }
    }

void _plinex(int x0,int y0,int z0,int x1,int y1,int z1) // helper edge render
    {
    int i,n,cx,cy,cz,sx,sy,sz,*by,*bz;
    // target buffer & order points by x
    if (x1>=x0)
        {
        i=x0; x0=x1; x0=i;
        i=y0; y0=y1; y0=i;
        i=z0; z0=z1; z0=i; by=ly; bz=lz;
        } else {           by=ry; bz=rz; }
    // line DDA parameters
    x1-=x0; sx=0; if (x1>0) sx=+1; if (x1<0) { sx=-1; x1=-x1; } if (x1) x1++;           n=x1;
    y1-=y0; sy=0; if (y1>0) sy=+1; if (y1<0) { sy=-1; y1=-y1; } if (y1) y1++; if (n<y1) n=y1;
    z1-=z0; sz=0; if (z1>0) sz=+1; if (z1<0) { sz=-1; z1=-z1; } if (z1) z1++; if (n<z1) n=z1;
    // single pixel (not a line)
    if (!n)
        {
        if ((x0>=0)&&(x0<vxs))
            {
            ly[x0]=y0; lz[x0]=z0;
            ry[x0]=y0; rz[x0]=z0;
            }
        return;
        }
    // ND DDA algo i is parameter
    for (cx=cy=cz=n,i=0;i<n;i++)
        {
        if ((x0>=0)&&(x0<vxs)){ by[x0]=y0; bz[x0]=z0; }
        cx-=x1; if (cx<=0){ cx+=n; x0+=sx; }
        cy-=y1; if (cy<=0){ cy+=n; y0+=sy; }
        cz-=z1; if (cz<=0){ cz+=n; z0+=sz; }
        }
    }

void _pliney(int x0,int y0,int z0,int x1,int y1,int z1) // helper edge render
    {
    int i,n,cx,cy,cz,sx,sy,sz,*bx,*bz;
    // target buffer & order points by y
    if (y1>=y0)
        {
        i=x0; x0=x1; x0=i;
        i=y0; y0=y1; y0=i;
        i=z0; z0=z1; z0=i; bx=lx; bz=lz;
        } else {           bx=rx; bz=rz; }
    // line DDA parameters
    x1-=x0; sx=0; if (x1>0) sx=+1; if (x1<0) { sx=-1; x1=-x1; } if (x1) x1++;           n=x1;
    y1-=y0; sy=0; if (y1>0) sy=+1; if (y1<0) { sy=-1; y1=-y1; } if (y1) y1++; if (n<y1) n=y1;
    z1-=z0; sz=0; if (z1>0) sz=+1; if (z1<0) { sz=-1; z1=-z1; } if (z1) z1++; if (n<z1) n=z1;
    // single pixel (not a line)
    if (!n)
        {
        if ((y0>=0)&&(y0<vys))
            {
            lx[y0]=x0; lz[y0]=z0;
            rx[y0]=x0; rz[y0]=z0;
            }
        return;
        }
    // ND DDA algo i is parameter
    for (cx=cy=cz=n,i=0;i<n;i++)
        {
        if ((y0>=0)&&(y0<vys)){ bx[y0]=x0; bz[y0]=z0; }
        cx-=x1; if (cx<=0){ cx+=n; x0+=sx; }
        cy-=y1; if (cy<=0){ cy+=n; y0+=sy; }
        cz-=z1; if (cz<=0){ cz+=n; z0+=sz; }
        }
    }

void _plinez(int x0,int y0,int z0,int x1,int y1,int z1) // helper edge render
    {
    int i,n,cx,cy,cz,sx,sy,sz,*bx,*by;
    // target buffer & order points by z
    if (z1>=z0)
        {
        i=x0; x0=x1; x0=i;
        i=y0; y0=y1; y0=i;
        i=z0; z0=z1; z0=i; bx=lx; by=ly;
        } else {           bx=rx; by=ry; }
    // line DDA parameters
    x1-=x0; sx=0; if (x1>0) sx=+1; if (x1<0) { sx=-1; x1=-x1; } if (x1) x1++;           n=x1;
    y1-=y0; sy=0; if (y1>0) sy=+1; if (y1<0) { sy=-1; y1=-y1; } if (y1) y1++; if (n<y1) n=y1;
    z1-=z0; sz=0; if (z1>0) sz=+1; if (z1<0) { sz=-1; z1=-z1; } if (z1) z1++; if (n<z1) n=z1;
    // single pixel (not a line)
    if (!n)
        {
        if ((z0>=0)&&(z0<vzs))
            {
            lx[z0]=x0; ly[z0]=y0;
            rx[z0]=x0; ry[z0]=y0;
            }
        return;
        }
    // ND DDA algo i is parameter
    for (cx=cy=cz=n,i=0;i<n;i++)
        {
        if ((z0>=0)&&(z0<vzs)){ bx[z0]=x0; by[z0]=y0; }
        cx-=x1; if (cx<=0){ cx+=n; x0+=sx; }
        cy-=y1; if (cy<=0){ cy+=n; y0+=sy; }
        cz-=z1; if (cz<=0){ cz+=n; z0+=sz; }
        }
    }

void triangle(int x0,int y0,int z0,int x1,int y1,int z1,int x2,int y2,int z2,float c) // render triangle
    {
    int X0,Y0,Z0,X1,Y1,Z1,dx,dy,dz,x,y,z;
    // BBOX
               X0=x0;            X1=x0;
    if (X0>x1) X0=x1; if (X1<x1) X1=x1;
    if (X0>x2) X0=x2; if (X1<x2) X1=x2;
               Y0=y0;            Y1=y0;
    if (Y0>y1) Y0=y1; if (Y1<y1) Y1=y1;
    if (Y0>y2) Y0=y2; if (Y1<y2) Y1=y2;
               Z0=z0;            Z1=z0;
    if (Z0>z1) Z0=z1; if (Z1<z1) Z1=z1;
    if (Z0>z2) Z0=z2; if (Z1<z2) Z1=z2;
    dx=X1-X0;
    dy=Y1-Y0;
    dz=Z1-Z0;
         if ((dx>=dy)&&(dx>=dz))    // x is major axis
        {
        // render circumference into left/right buffers
        _plinex(x0,y0,z0,x1,y1,z1);
        _plinex(x1,y1,z1,x2,y2,z2);
        _plinex(x2,y2,z2,x0,y0,z0);
        // fill the triangle
        if (X0<0) X0=0; if (X1>=vxs) X1=vxs-1;
        for (x=X0;x<=X1;x++)
            {
            y0=ly[x]; z0=lz[x];
            y1=ry[x]; z1=rz[x];
            line(x,y0,z0,x,y1,z1,c);
            }
        }
    else if ((dy>=dx)&&(dy>=dz))    // y is major axis
        {
        // render circumference into left/right buffers
        _pliney(x0,y0,z0,x1,y1,z1);
        _pliney(x1,y1,z1,x2,y2,z2);
        _pliney(x2,y2,z2,x0,y0,z0);
        // fill the triangle
        if (Y0<0) Y0=0; if (Y1>=vys) Y1=vys-1;
        for (y=Y0;y<=Y1;y++)
            {
            x0=lx[y]; z0=lz[y];
            x1=rx[y]; z1=rz[y];
            line(x0,y,z0,x1,y,z1,c);
            }
        }
    else if ((dz>=dx)&&(dz>=dy))    // z is major axis
        {
        // render circumference into left/right buffers
        _plinez(x0,y0,z0,x1,y1,z1);
        _plinez(x1,y1,z1,x2,y2,z2);
        _plinez(x2,y2,z2,x0,y0,z0);
        // fill the triangle
        if (Z0<0) Z0=0; if (Z1>=vzs) Z1=vzs-1;
        for (z=Z0;z<=Z1;z++)
            {
            x0=lx[z]; y0=ly[z];
            x1=rx[z]; y1=ry[z];
            line(x0,y0,z,x1,y1,z,c);
            }
        }
    }

The edge functions _pline can be speeded up, if you realize the major axis is the same as iteration axis i you can have one computation less. Also the 3 incarnations can be merged to single one (just by reordering the x,y,z order of called parameters and passing the left/right buffers pointers. And the moving the clipping outside DDA iteration hould boost performance quite a lot too.

Here preview of triangle intersecting cylinder:

preview1

And here preview of material removal using this triangle rasterization:

preview2

The orange triangle is not voxel one its just overlay, the voxel one is not rendered directly but substracted from the voxel grid instead to remove material during animation ...