Funkcja podwójnego konturowania i kwadratowego błędu

Question

Wprowadziłem kostki maszerujące, podwójne kostki maszerujące i kostki marszu adaptacyjnego w języku C#, tylko po to, aby przekonać się, że do moich celów potrzebuję podwójnego konturowania. Przeczytałem wszystkie prace dotyczące podwójnego konturowania i otrzymuję wszystko poza rdzeniem samego podwójnego konturowania: minimalizowanie kwadratowej funkcji błędu (QEF).

W tej chwili obliczam pozycję wierzchołków wewnętrznych wokseli, po prostu znajdując średnią między wszystkimi krawędziami dzielącymi ten pojedynczy wierzchołek (od 3 do 6 krawędzi) i działa dobrze, ale oczywiście nie tworzy wewnętrznych wierzchołków w odpowiednie miejsca.

Oto fragment kodu, który próbuję utworzyć. Każda pomoc byłaby bardzo cenna.

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

Trudno jest zrozumieć QEF. Mam nadzieję, że mogę pomóc. Metoda podwójnego konturowania oblicza dane "Hermite" w każdym punkcie przecięcia, lub innymi słowy, w każdym punkcie utworzonym na krawędzi woksela, gdzie znana jest normalna powierzchnia. Z punktem i normalnym można obliczyć równanie płaszczyzny.

QEF to suma kwadratów odległości od punktu wewnętrznego woksela do każdej z płaszczyzn związanych z wokselem. Poniżej znajduje się pseudo kod do obliczania 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); 
} 

Celem jest wybranie punktu wewnątrz woksela, który minimalizuje QEF. Literatura sugeruje użycie procesu Grahm-Schmidta do zlokalizowania optymalnego punktu, ale może to być złożone, a także może prowadzić do punktów leżących poza wokselem.

Inną opcją (hack-ish) jest utworzenie siatki punktów wewnątrz woksela i obliczyć QEF dla każdego z nich i wybrać ten z najniższą, im mniejszą siatkę, tym bliżej punktu optymalnego, który przyjdziesz ale im dłuższe są obliczenia.

Answer 2

W mojej obecnej implementacji podwójnego konturowania im za pomocą bardzo prostego sposobu rozwiązania QEF. Ponieważ w skrócie QEF jest przybliżeniem najmniejszych kwadratów, znalazłem najprostszy sposób obliczenia QEF przez obliczenie pseudo-odwrotności. Ta pseudo-odwrotność może być obliczona za pomocą dowolnej biblioteki algebraicznej w twoim języku.

Jest to kod używam:

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

wejścia funkcji są punkty przecięcia i normalne, a meanPoint jest średnią dany intersectoin punktów.

Podsumowanie matematyki: ta funkcja oblicza punkt leżący na przecięciu wszystkich płaszczyzn zdefiniowanych przez punkty przecięcia i linie normalne. Ponieważ nie ma to dokładnego rozwiązania, obliczane jest przybliżenie najmniejszych kwadratów, które określa punkt, który "najmniej się myli". Dodatkowo punkty przecięcia są "przesuwane" tak, że punkt środkowy staje się punktem początkowym. Zapewnia to, że gdy istnieje wiele rozwiązań do QEF, wybierane jest rozwiązanie najbliższe średniej.