Функция двойного контура и квадратичной ошибки

Question

Я реализовал марширующие кубы, двойные марширующие кубы и адаптивные марширующие кубы в C# только для того, чтобы узнать, что мне нужно двойное контурирование для моих целей. Я прочитал все работы о двойном контуре, и я получил все, кроме сути самого двойного контура: минимизация функции квадратичной ошибки (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!
    }

Answer 1

QEF довольно сложно понять. Надеюсь, я смогу помочь. Метод двойного контура вычисляет данные "Эрмита" в каждой точке пересечения или, другими словами, в каждой точке, созданной на краю вокселя, известна нормаль поверхности. С точкой и нормалью можно вычислить уравнение плоскости.

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. В литературе предлагается использовать процесс Грамма-Шмидта, чтобы найти оптимальную точку, но это может быть сложным, а также может привести к точкам, которые лежат за пределами вокселя.

Другой вариант (hack-ish) состоит в том, чтобы создать сетку точек внутри вокселя и рассчитать QEF для каждого и выбрать ту, которая имеет наименьшее значение. Чем тоньше сетка, тем ближе к оптимальной точке вы придете, но чем дольше расчеты.

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 выбирается решение, наиболее близкое к средней точке.