Algorithm inorder to Convert a Triangular Facet to Voxelised Model

Question

I have the following vertices: (x1,y1,z1),(x2,y2,z2) and (x3,y3,z3). I have followed the below approach:

1. Calculating the bounding box of the facet.

   The length of the bounding box is Maximum(xmax-xmin, ymax-ymin,zmax-zmin) where xmax,xmin,....,zmin can be calculated by looping through all the vertices of the triangle.

2. Now,considering a particular resolution size(h), I divide my bounding box into grids.

   For example, if my bounding box is of length 10, No. of grids will be 1000(10*10*10).

3. I now consider vertex1,vertex2 and vertex2,vertex3 and vertex1,vertex3 and apply Bresenham algorithm.

   Now,I use the following Bresenham algorithm from this link.

   function visitAll(gx0, gy0, gz0, gx1, gy1, gz1, visitor) {
   
       var gx0idx = Math.floor(gx0);
       var gy0idx = Math.floor(gy0);
       var gz0idx = Math.floor(gz0);
   
       var gx1idx = Math.floor(gx1);
       var gy1idx = Math.floor(gy1);
       var gz1idx = Math.floor(gz1);
   
       var sx = gx1idx > gx0idx ? 1 : gx1idx < gx0idx ? -1 : 0;
       var sy = gy1idx > gy0idx ? 1 : gy1idx < gy0idx ? -1 : 0;
       var sz = gz1idx > gz0idx ? 1 : gz1idx < gz0idx ? -1 : 0;
   
       var gx = gx0idx;
       var gy = gy0idx;
       var gz = gz0idx;
   
       //Planes for each axis that we will next cross
       var gxp = gx0idx + (gx1idx > gx0idx ? 1 : 0);
       var gyp = gy0idx + (gy1idx > gy0idx ? 1 : 0);
       var gzp = gz0idx + (gz1idx > gz0idx ? 1 : 0);
   
       //Only used for multiplying up the error margins
       var vx = gx1 === gx0 ? 1 : gx1 - gx0;
       var vy = gy1 === gy0 ? 1 : gy1 - gy0;
       var vz = gz1 === gz0 ? 1 : gz1 - gz0;
   
       //Error is normalized to vx * vy * vz so we only have to multiply up
       var vxvy = vx * vy;
       var vxvz = vx * vz;
       var vyvz = vy * vz;
   
       //Error from the next plane accumulators, scaled up by vx*vy*vz
       // gx0 + vx * rx === gxp
       // vx * rx === gxp - gx0
       // rx === (gxp - gx0) / vx
       var errx = (gxp - gx0) * vyvz;
       var erry = (gyp - gy0) * vxvz;
       var errz = (gzp - gz0) * vxvy;
   
       var derrx = sx * vyvz;
       var derry = sy * vxvz;
       var derrz = sz * vxvy;
   
       do {
           visitor(gx, gy, gz);
   
           if (gx === gx1idx && gy === gy1idx && gz === gz1idx) break;
   
           //Which plane do we cross first?
           var xr = Math.abs(errx);
           var yr = Math.abs(erry);
           var zr = Math.abs(errz);
   
           if (sx !== 0 && (sy === 0 || xr < yr) && (sz === 0 || xr < zr)) {
               gx += sx;
               errx += derrx;
           }
           else if (sy !== 0 && (sz === 0 || yr < zr)) {
               gy += sy;
               erry += derry;
           }
           else if (sz !== 0) {
               gz += sz;
               errz += derrz;
           }
   
       } while (true);
   }
   

4. From the above step, I obtain some sample points.I loop through each sample and obtain the grid,in which this sampling point lies.

So How do I fill my inside triangular region with grids(or voxels) ?

I would be really glad,if someone can correct my approach and also provide an efficient algorithm for filling voxels inside my triangular facet.