Mesh Voxelizer in C#

Author: Adok

A mesh voxelizer is an algorithm that converts arbitrary three-dimensional objects, usually specified as polygons, into voxels (vector objects), that is a collection of small rectangular objects. Such an algorithm is sometimes needed for processing and visualizing 3D objects. Some theoretical background on how a mesh voxelizer works is provided by the paper "An Accurate Method for Voxelizing Polygon Meshes". In this article I would like to present you my solution in C#.

My solution uses a simple data structure called Vert that can be defined as:

    public class Vert
    {
        public double x, y, z;
    }

Moreover, it uses Verts, which is a list of elements of the Vert data type:

using Verts = System.Collections.Generic.List<Vert>;

The core of the mesh voxelizer is the following method:

        public bool IsVoxelOnFace(double[] vc, double L, double Lx, double Ly, double Lz)
        {
            if (this.Verts.Count < 3) return false;

            // 1. representation of vertices
            double rc = Rc(L, 26);

            foreach (Vert vx in this.Verts)
            {
                if ((vc[0] - vx.x) * (vc[0] - vx.x) +
                    (vc[1] - vx.y) * (vc[1] - vx.y) +
                    (vc[2] - vx.z) * (vc[2] - vx.z) <= rc * rc)
                {
                    return true;
                }
            }

            // 2. representation of edges
            int stop_count = this.Verts.Count;
            if (!this.IsSorted) this.Sort();
            for (int i = 0; i < stop_count; i++)
            {
                Vert vx1 = this.Verts.Item[i];
                Vert vx2 = this.Verts.Item[(i + 1) % this.Verts.Count];

                double xa = vx1.x;
                double ya = vx1.y;
                double za = vx1.z;
                double xb = vx2.x;
                double yb = vx2.y;
                double zb = vx2.z;
                double xp = vc[0];
                double yp = vc[1];
                double zp = vc[2];

                if (Math.Min(xa, xb) <= xp + rc &&
                    Math.Min(ya, yb) <= yp + rc &&
                    Math.Min(za, zb) <= zp + rc &&
                    Math.Max(xa, xb) >= xp - rc &&
                    Math.Max(ya, yb) >= yp - rc &&
                    Math.Max(za, zb) >= zp - rc)
                {
                    double[] w = new double[3];
                    double[] v1 = new double[3] { xp - xa, yp - ya, zp - za };
                    double[] v2 = new double[3] { xb - xa, yb - ya, zb - za };
                    w[0] = v1[1] * v2[2] - v1[2] * v2[1];
                    w[1] = v1[2] * v2[0] - v1[0] * v2[2];
                    w[2] = v1[0] * v2[1] - v1[1] * v2[0];

                    double div = v2[0] * v2[0] + v2[1] * v2[1] + v2[2] * v2[2];
                    double d = (w[0] * w[0] + w[1] * w[1] + w[2] * w[2]) / (div != 0 ? div : 1);

                    if (d <= rc * rc) return true;
                }
            }

            // 3. representation of body
            double[] vc_proj = Apply_Rot_Matrix(ToOrigin(vc));
            return bboxtest(vc_proj, Math.Abs(t(Lz, 26)), Lx / sqrt2, Ly / sqrt2) &&
                (wn_PnPoly(vc_proj, Lx, Ly));
        }

It checks whether a given voxel (the coordinates of the midpoint are provided in the array vc) hits a polygon. The parameter L is the approximate size of the voxel, assuming it is quadratic. To support non-quadratic voxels, also the parameters Lx, Ly and Lz can be specified.

Now this method relies on several other methods which I will present you here now:

        public double t26(double L)
        {
            return Rc26(L); 
        }

        public double Rc26(double L)
        {
            return L / 2 * 1.732;
        }

        public double t6(double L)
        {
            return Rc6(L);

        }

        public double Rc6(double L)
        {
            return L / 2;
        }

        public double t(double L, int which = 6)
        {
            switch (which)
            {
                case 6: return t6(L);
                case 26: return t26(L);
                default: return 0;
            }
        }

        public double Rc(double L, int which = 6)
        {
            switch (which)
            {
                case 6: return Rc6(L);
                case 26: return Rc26(L);
                default: return 0;
            }
        }

As you see these methods are parameterized so several variants of it can be used. Basically they are to determine the distance from the midpoint to the corners of a voxel.

        public void Sort()
        {
            bool[] plane = this.Plane;

            Verts copy = this.Verts.Clone();
            if (plane[2] == true)
            {
                copy.Item.Sort(
                    delegate(Vert v1, Vert v2)
                    {
                        return v1.x.CompareTo(v2.x);
                    }
                    );
            }
            else if (plane[1] == true)
            {
                copy.Item.Sort(
                    delegate(Vert v1, Vert v2)
                    {
                        return v1.z.CompareTo(v2.z);
                    }
                    );
            }
            else
            {
                copy.Item.Sort(
                    delegate(Vert v1, Vert v2)
                    {
                        return v1.y.CompareTo(v2.y);
                    }
                    );
            }

            Verts half1 = copy.Clone();
            Verts half2 = new Verts();

            if (plane[2] == true)
            {
                for (int i = 0; i < half1.Item.Count - 2; )
                {
                    Vert a = half1.Item[i];
                    Vert b = half1.Item[i + 1];
                    Vert c = half1.Item[i + 2];

                    if ((b.y - a.y) * (c.x - b.x) > (c.y - b.y) * (b.x - a.x))
                    {
                        half2.Add(b);
                        half1.Remove(i + 1);
                    }
                    else
                    {
                        i++;
                    }
                }
            }
            else if (plane[1] == true)
            {
                for (int i = 0; i < half1.Item.Count - 2; )
                {
                    Vert a = half1.Item[i];
                    Vert b = half1.Item[i + 1];
                    Vert c = half1.Item[i + 2];

                    if ((b.x - a.x) * (c.z - b.z) > (c.x - b.x) * (b.z - a.z))
                    {
                        half2.Add(b);
                        half1.Remove(i + 1);
                    }
                    else
                    {
                        i++;
                    }
                }
            }
            else if (plane[0] == true)
            {
                for (int i = 0; i < half1.Item.Count - 2; )
                {
                    Vert a = half1.Item[i];
                    Vert b = half1.Item[i + 1];
                    Vert c = half1.Item[i + 2];

                    if ((b.z - a.z) * (c.y - b.y) > (c.z - b.z) * (b.y - a.y))
                    {
                        half2.Add(b);
                        half1.Remove(i + 1);
                    }
                    else
                    {
                        i++;
                    }
                }
            }
            else
            {
                double[,] rot_matrix = Rot_Matrix;

                Verts verts_test = new Verts();
                for (int i = 0; i < half1.Item.Count; i++)
                {
                    Vert a = ToOrigin(half1.Item[i]);
                    Vert new_vert = new Vert();
                    new_vert.x = Apply_Rot_Matrix_X(a);
                    new_vert.y = Apply_Rot_Matrix_Y(a);
                    verts_test.Add(new_vert);
                }

                for (int i = 0; i < verts_test.Item.Count - 2; )
                {
                    Vert a = verts_test.Item[i];
                    Vert b = verts_test.Item[i + 1];
                    Vert c = verts_test.Item[i + 2];

                    if ((b.y - a.y) * (c.x - b.x) > (c.y - b.y) * (b.x - a.x))
                    {
                        half2.Add(half1.Item[i + 1]);
                        half1.Remove(i + 1);
                        verts_test.Remove(i + 1);
                    }
                    else
                    {
                        i++;
                    }
                }
            }

            for (int i = half2.Count; i > 0; i--)
            {
                half1.Add(half2.Item[i - 1]);
            }

            this.Verts = half1;

            IsSorted = true;
        }

This method is for sorting the vertices. It itself depends on some auxiliary methods to apply the rotation matrix:

        public double Apply_Rot_Matrix(double[] v, int t)
        {
            double[,] rot_matrix = Rot_Matrix;
            return rot_matrix[t, 0] * v[0] + rot_matrix[t, 1] * v[1] + rot_matrix[t, 2] * v[2];
        }

        public double Apply_Rot_Matrix_X(double[] v)
        {
            return Apply_Rot_Matrix(v, 0);
        }

        public double Apply_Rot_Matrix_Y(double[] v)
        {
            return Apply_Rot_Matrix(v, 1);
        }

        public double Apply_Rot_Matrix_Z(double[] v)
        {
            return Apply_Rot_Matrix(v, 2);
        }

        public double Apply_Rot_Matrix(Vert v, int t)
        {
            double[,] rot_matrix = Rot_Matrix;
            return rot_matrix[t, 0] * v.x + rot_matrix[t, 1] * v.y + rot_matrix[t, 2] * v.z;
        }

        public double Apply_Rot_Matrix_X(Vert v)
        {
            return Apply_Rot_Matrix(v, 0);
        }

        public double Apply_Rot_Matrix_Y(Vert v)
        {
            return Apply_Rot_Matrix(v, 1);
        }

        public double Apply_Rot_Matrix_Z(Vert v)
        {
            return Apply_Rot_Matrix(v, 2);
        }

        public Vert Apply_Rot_Matrix(Vert v)
        {
            Vert v_new = new Vert();

            double[,] rot_matrix = Rot_Matrix;
            v_new.x = rot_matrix[0, 0] * v.x + rot_matrix[0, 1] * v.y + rot_matrix[0, 2] * v.z;
            v_new.y = rot_matrix[1, 0] * v.x + rot_matrix[1, 1] * v.y + rot_matrix[1, 2] * v.z;
            v_new.z = rot_matrix[2, 0] * v.x + rot_matrix[2, 1] * v.y + rot_matrix[2, 2] * v.z;

            return v_new;
        }

        public double[] Apply_Rot_Matrix(double[] v)
        {
            double[] v_new = new double[3];

            double[,] rot_matrix = Rot_Matrix;
            v_new[0] = rot_matrix[0, 0] * v[0] + rot_matrix[0, 1] * v[1] + rot_matrix[0, 2] * v[2];
            v_new[1] = rot_matrix[1, 0] * v[0] + rot_matrix[1, 1] * v[1] + rot_matrix[1, 2] * v[2];
            v_new[2] = rot_matrix[2, 0] * v[0] + rot_matrix[2, 1] * v[1] + rot_matrix[2, 2] * v[2];

            return v_new;
        }

        private double[,] _rot_matrix;
        private bool _rot_matrix_defined = false;

        public double[,] Rot_Matrix
        {
            get
            {
                if (_rot_matrix_defined) return _rot_matrix;

                double[] v1 = NormalVector_normalized;
                double[] v2 = new double[3] { 0, 0, 1 };
                double[] v = new double[3];

                v1 = NormalVector_normalized;

                v[0] = v1[1] * v2[2] - v1[2] * v2[1];
                v[1] = v1[2] * v2[0] - v1[0] * v2[2];
                v[2] = v1[0] * v2[1] - v1[1] * v2[0];

                double v_abs = Math.Sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
                if (v_abs != 0)
                {
                    v[0] /= v_abs;
                    v[1] /= v_abs;
                    v[2] /= v_abs;
                }

                double cos_angle = v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
                double angle = Math.Acos(cos_angle);
                double sin_angle = Math.Sin(angle);

                double[,] rot_matrix = new double[3, 3];

                double cos_angle_inv = 1 - cos_angle;
                double a0 = v[0] * cos_angle_inv;
                double a1 = v[1] * cos_angle_inv;
                double b0 = v[0] * sin_angle;
                double b1 = v[1] * sin_angle;
                double b2 = v[2] * sin_angle;
                double c10 = v[1] * a0;
                double c20 = v[2] * a0;
                double c21 = v[2] * a1;

                rot_matrix[0, 0] = v[0] * a0 + cos_angle;
                rot_matrix[0, 1] = c10 - b2;
                rot_matrix[0, 2] = c20 + b1;
                rot_matrix[1, 0] = c10 + b2;
                rot_matrix[1, 1] = v[1] * a1 + cos_angle;
                rot_matrix[1, 2] = c21 - b0;
                rot_matrix[2, 0] = c20 - b1;
                rot_matrix[2, 1] = c21 + b0;
                rot_matrix[2, 2] = v[2] * v[2] * cos_angle_inv + cos_angle;

                _rot_matrix = rot_matrix;
                _rot_matrix_defined = true;

                return rot_matrix;
            }
        }

These auxiliary methods in return make use of a data structure that stores the normal vector of a face:

        private double[] _normalVector;
        private bool _normalVector_defined = false;

        public double[] NormalVector
        {
            get
            {
                if (_normalVector_defined) return _normalVector;

                double[] v1 = new double[3];
                double[] v2 = new double[3];
                double[] v = new double[3];

                v1[0] = this.Verts.Item[1].x - this.Verts.Item[0].x;
                v1[1] = this.Verts.Item[1].y - this.Verts.Item[0].y;
                v1[2] = this.Verts.Item[1].z - this.Verts.Item[0].z;

                // to avoid parallel lines
                for (int i = 2; i < this.Verts.Count; i++)
                {
                    v2[0] = this.Verts.Item[i].x - this.Verts.Item[i - 1].x;
                    v2[1] = this.Verts.Item[i].y - this.Verts.Item[i - 1].y;
                    v2[2] = this.Verts.Item[i].z - this.Verts.Item[i - 1].z;

                    if ((Math.Abs(v1[0] / v2[0]) != Math.Abs(v1[1] / v2[1])) || (Math.Abs(v1[0] / v2[0]) != Math.Abs(v1[2] / v2[2]))) break;
                }

                v[0] = v1[1] * v2[2] - v1[2] * v2[1];
                v[1] = v1[2] * v2[0] - v1[0] * v2[2];
                v[2] = v1[0] * v2[1] - v1[1] * v2[0];

                _normalVector = v;
                _normalVector_defined = true;

                return v;
            }
        }

        private double[] _normalVector_normalized;
        private bool _normalVector_normalized_defined = false;

        public double[] NormalVector_normalized
        {
            get
            {
                if (_normalVector_normalized_defined) return _normalVector_normalized;

                double[] v = NormalVector;

                double a = Math.Sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);

                if (a == 0) return v;

                v[0] /= a;
                v[1] /= a;
                v[2] /= a;

                _normalVector_normalized = v;
                _normalVector_normalized_defined = true;

                return v;
            }
        }

Moreover, this method is used:

        public Vert ToOrigin(Vert v)
        {
            double[] n = NormalVector_normalized_times_lambda;
            Vert v_new = new Vert();

            v_new.x = v.x - n[0];
            v_new.y = v.y - n[1];
            v_new.z = v.z - n[2];

            return v_new;
        }

        public double[] ToOrigin(double[] v)
        {
            double[] n = NormalVector_normalized_times_lambda;
            double[] v_new = new double[3];

            v_new[0] = v[0] - n[0];
            v_new[1] = v[1] - n[1];
            v_new[2] = v[2] - n[2];

            return v_new;
        }

Next, bboxtest and Verts_proj, for testing whether a vertex is within a given bounding box:

        public bool bboxtest(double[] vc_proj, double zdist, double Lx_divided_sqrt2, double Ly_divided_sqrt2)
        {
            if (this.Verts_proj == null) return false;

            Dbl_3D_BBox bbox = this.Verts_proj.BoundingBox();

            if (bbox == null) return false;

            if (vc_proj[2] < bbox.zmin - zdist) return false;
            if (vc_proj[2] > bbox.zmax + zdist) return false;
            if (vc_proj[0] < bbox.xmin - Lx_divided_sqrt2) return false;
            if (vc_proj[0] > bbox.xmax + Lx_divided_sqrt2) return false;
            if (vc_proj[1] < bbox.ymin - Ly_divided_sqrt2) return false;
            if (vc_proj[1] > bbox.ymax + Ly_divided_sqrt2) return false;

            return true;
        }

        public Verts _verts_proj;
        public bool _verts_proj_defined = false;

        public Verts Verts_proj
        {
            get
            {
                if (this.Verts.Count < 3) return null;

                if (_verts_proj_defined) return _verts_proj;

                if (!this.IsSorted) this.Sort();

                Verts verts_proj = new Verts();

                for (int i = 0; i < this.Verts.Count; i++)
                {
                    Vert v = this.Verts.Item[i];
                    verts_proj.Add(Apply_Rot_Matrix(ToOrigin(v)));
                }

                _verts_proj_defined = true;
                _verts_proj = verts_proj;

                return verts_proj;
            }
        } 

Finally, we need to check whether a vertex is inside or outside of a polygon, which we do using these methods:

        public double isLeft(Vert P0, Vert P1, double x, double y)
        {
            return ( (P1.x - P0.x) * (y - P0.y) - (x - P0.x) * (P1.y - P0.y) );
        }

        public bool wn_PnPoly(double[] vc, double Lx, double Ly)
        {
            int wn = 0;
            int n = this.Verts_proj.Count;
            List<Vert> V = this.Verts_proj.Item;

            double x_min = vc[0] - Lx / sqrt2;
            double y_min = vc[1] - Ly / sqrt2;
            double x_max = vc[0] + Lx / sqrt2;
            double y_max = vc[1] + Ly / sqrt2;

            if (Double.IsNaN(x_min) || Double.IsNaN(x_max) || Double.IsNaN(y_min) || Double.IsNaN(y_max)) return false;

            for (double x = x_min; x <= x_max; x = x_min == x_max ? x_max + 1 : x + (x_max - x_min) / 2)
            {
                for (double y = y_min; y <= y_max; y = y_min == y_max ? y_max + 1 : y + (y_max - y_min) / 2)
                {
                    wn = 0;
                    for (int i = 0; i < n; i++)
                    {
                        int j = (i + 1) % n;
                        if (V[i].y <= y)
                        {
                            if (V[j].y > y)
                                if (isLeft(V[i], V[j], x, y) > 0)
                                    ++wn;
                        }
                        else
                        {
                            if (V[j].y <= y)
                                if (isLeft(V[i], V[j], x, y) < 0)
                                    --wn;
                        }
                    }
                    if (wn != 0) return true;
                }
            }
            return wn != 0;
        }

Oh, and sqrt2 is defined as:

public static double sqrt2 = 1.4142;

So to use these methods in a mesh voxelizer, you would run a for loop that iterates through all the voxels and then calls IsVoxelOnFace to check whether the voxel hits the polygon. Do that for all polygons you have and you will get a nice voxel representation of your mesh.

Note that this code is for the CPU and could be further optimized by rewriting it for GPU.

I hope this helps you write your own mesh voxelizer.