双重轮廓和二次误差函数

Question

我已经用C#实现了 marching cubes、dual marching cubes 和自适应 marching cubes，但我发现我需要 dual contouring 来达到我的目的。我已经阅读了所有关于 dual contouring 的资料，但是我还没有完全理解 dual contouring 的核心：最小化二次误差函数（QEF）。

目前，我通过找到共享单个顶点的所有边缘点（3 到 6 条边）之间的平均值来计算内部体素的顶点位置，这样做效果不错，但显然不能在正确的位置创建内部顶点。

这里是我正在尝试创建的代码片段。非常感谢任何帮助。

/// <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!
    }

- Fario Bostanzo

---

这里有一个不错的解释链接，并且还有Python实现。 - Fedor

Answer 1

对于理解QEF来说可能有些困难，但我希望能够帮助。双重轮廓法通过计算在每个交点处的'Hermite'数据（或者换句话说，在已知表面法线的每个沿体素边缘创建的点处）来进行计算。使用一个点和一个法线可以计算一个平面的方程。

QEF是将体素内部点到与体素相关联的每个平面的距离平方和。下面是一些伪代码用于计算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); 
}

目标是选择一个点在体素内，使QEF最小化。文献建议使用Grahm-Schmidt过程来定位最优点，但这可能很复杂，并且可能导致点位于体素外部。

另一种选择（hack-ish）是在体素内创建一个点网格，为每个点计算QEF，并选择具有最低QEF的点，网格越细，接近最优点的距离就越近，但计算时间会更长。

- GCB613

Answer 2

在我的双轮廓实现中，我使用了一种非常简单的方法来解决QEF问题。由于QEF本质上是最小二乘逼近，我发现计算QEF的最简单方法是通过计算伪逆来进行。这个伪逆可以使用您语言中的任何代数库来计算。

以下是我使用的代码:

    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 });
    }

该函数的输入是交点和法线，meanPoint是给定交点的平均值。

简单概括一下数学：该函数计算位于由交点和法线定义的所有平面的交点上的点。由于这没有确切的解决方案，因此计算最小二乘逼近，找到“最少错误”的点。此外，交点“移动”，使得平均点成为原点。这确保了当QEF有多个解时，选择离平均点最近的解。

- MHGameWork