Dual Contouring and Quadratic Error Function

Question

I've implemented marching cubes, dual marching cubes and adaptive marching cubes in C#, only to find out that I need dual contouring for my purposes. I've read all works about dual contouring and I get all but the core of the dual contouring itself: minimizing the quadratic error function (QEF).

Right now, I'm calculating internal voxel's vertices position simply by finding the mean between all edgePoints sharing that single vertex (3 to 6 edges) and it works well, but it obviously doesn't create the internal vertices in the right places.

Here is the piece of code I'm trying to create. Any help would be very 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

The QEF is rather difficult to understand. Hopefully I can help. The dual contouring method calculates the 'Hermite' data at each crossing point, or in other words, at each point created on the edge of the voxel the normal of the surface is known. With a point and a normal one can calculate the equation of a plane.

The QEF is the sum of the squares of the distances from the voxel's internal point to each of the planes associated with the voxel. Below is some pseudo code for calculating the 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 then to choose a point inside the voxel that minimizes the QEF. The literature suggests using the Grahm-Schmidt process to locate the optimal point but this can be complex and also can result in points that lie outside of the voxel.

Another option (hack-ish) is to create a grid of points inside the voxel and calculate the QEF for each one and pick the one with the lowest, the finer the grid the closer to the optimal point you'll come but the longer the calculations.

Answer 2

In my current implementation of dual contouring im using a very simple way to solve the QEF. Since the QEF in essence is a least squares approximation, I've found the simplest way to calculate the QEF is by calculating the pseudoinverse. 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[] { e.X, e.Y, e.Z }).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 });
    }

The function's inputs are the intersection points and normals, and the meanPoint is the average of the given intersectoin points.

Summarizing the math: this function calculates the point that lies on the intersection of all planes defined by the intersection points and normals. Since this doesn't have an exact solution, the least squares approximation is calculated, which finds the point which is 'least wrong'. Additionally, the intersection points are 'moved' so that the mean point becomes the origin. This ensures that when there are multiple solutions to the QEF, the solution closest to the mean point is chosen.