Double contour and quadratic error function

Question

I implemented marching cubes, double marching cubes and adaptive marching cubes in C #, only to find out that I need double circuits for my purposes. I read all the works on double contouring, and I get everything except the core of the double contour: minimizing the quadratic error function (QEF).

Right now, I am calculating the internal position of the voxel versels by simply finding the average between all the boundary points that use this single vertex (3 to 6 edges) and it works well, but it obviously does not create internal vertices on the right places.

Here is the code snippet I'm trying to create. Any help would be much appreciated

/// <summary> /// ORIGINAL WORK: Dual Contouring of Hermite Data by Tao Ju (remember me of a MechCommander 2 character) /// 2.3 Representing and minimizing QEFs /// The function E[x] can be expressed as the inner /// product (Ax-b)T (Ax-b) where A is a matrix whose rows are the /// normals ni and b is a vector whose entries are ni*pi. <------------ (dot product?)> /// Typically, the quadratic function E[x] is expanded into the form /// E[x] = xT AT Ax - 2xT AT b + bT b (2) /// where the matrix AT A is a symmetric 3x3 matrix, AT b is a column /// vector of length three and bT b is a scalar. The advantage of this expansion /// is that only the matrices AT A, AT b and bT b need be stored /// (10 floats), as opposed to storing the matrices A and b. Furthermore, /// a minimizing value ˆ x for E[x] can be computed by solving /// the normal equations AT Aˆ x = AT b. /// </summary> public Vector3 GetMinimumError(Vector3 p0, Vector3 p1, Vector3 p2, Vector3 n0, Vector3 n1, Vector3 n2) { //so, here we are. I'm creating a vector to store the final value. Vector3 position = Vector3.Zero; //Values of b are supposed to b (:P) three floats. The only way i know to find a float value //by multiplying 2 vectors is to use dot product. Vector3 b = new Vector3( Vector3.Dot(p0, n0), Vector3.Dot(p1, n1), Vector3.Dot(p2, n2)); //What the transpose of a vector is supposed to be? //I don't know, but i think should be the vector itself :) float bTb = Vector3.Dot(b, b); //i create a square matrix 3x3, so i can use c# matrix transformation libraries. //i know i will probably have to build bigger matrix later on, but it should fit for now Matrix A = new Matrix( n0.X, n0.Y, n0.Z, 0, n1.X, n1.Y, n1.Z, 0, n2.X, n2.Y, n2.Z, 0, 0, 0, 0, 0); //easy Matrix AT = Matrix.Transpose(A); //EASY Matrix ATA = Matrix.Multiply(AT, A); //Another intuition. Hope makes sense... Vector3 ATb = Vector3.Transform(b, AT); //... // some cool stuff about solving // the normal equations AT Aˆ x = AT b //... return position; //profit! } 

Answer 1

QEF is pretty hard to understand. Hope I can help. The double-contour method calculates Hermite data at each intersection, or, in other words, at each point created on the edge of the voxel, a normal surface is known. With a point and a normal, you can calculate the plane equation.

QEF is the sum of the squared distances from the inner voxel point to each of the planes associated with the voxel. Below is some pseudo code for calculating QEF.

 double get_QEF(Point3d point, Voxel3d voxel) { double QEF = 0.0; foreach(plane in voxel.planes) { double dist_to_plane = plane.distance(point); QEF += dist_to_plane*dist_to_plane; } return(QEF); } 

The goal is to select a point within the voxel that minimizes QEF. It is suggested in the literature to use the Graham-Schmidt process to determine the optimal point, but this can be complex and can also lead to points that lie outside the voxel.

Another option (hack-ish) is to create a grid of points inside the voxel and calculate the QEF for each of them and choose the one with the lowest, the more precisely the grid is closer to the optimal point that you will arrive, but the more calculations.

Answer 2

In my current dual-circuit im implementation, using a very simple way to solve QEF. Since QEF is essentially a least squares approximation, I found the easiest way to calculate QEF by calculating a pseudo-inverse. This pseudoinverse can be calculated using any algebraic library in your language.

This is the code I'm using:

  public static Vector<float> CalculateCubeQEF(Vector3[] normals, Vector3[] positions, Vector3 meanPoint) { var A = DenseMatrix.OfRowArrays(normals.Select(e => new[] { eX, eY, eZ }).ToArray()); var b = DenseVector.OfArray(normals.Zip(positions.Select(p => p - meanPoint), Vector3.Dot).ToArray()); var pseudo = PseudoInverse(A); var leastsquares = pseudo.Multiply(b); return leastsquares + DenseVector.OfArray(new[] { meanPoint.X, meanPoint.Y, meanPoint.Z }); } 

Function inputs are the intersection and normal points, and the midpoint is the average of these intersection points.

To summarize the math: this function calculates a point lying at the intersection of all planes defined by intersection points and normals. Since this does not have an exact solution, the least squares approximation is calculated, which finds the “least erroneous” point. In addition, the intersection points “move”, so that the midpoint becomes the beginning. This ensures that if there are multiple solutions for QEF, the solution closest to the midpoint is selected.