Funzione di doppio errore di contornitura e quadratica

Question

Ho implementato cubetti in marcia, cubi a doppia marcia e cubetti di marchingegno adattivi in ​​C#, solo per scoprire che per i miei scopi ho bisogno di un doppio contorno. Ho letto tutti i lavori sul dual contouring e ottengo tutto tranne il nucleo del dual contouring stesso: minimizzando la funzione di errore quadratico (QEF).

In questo momento, sto calcolando la posizione dei vertici interni del voxel semplicemente trovando la media tra tutti i edgePoints che condividono quel singolo vertice (da 3 a 6 spigoli) e funziona bene, ma ovviamente non crea i vertici interni nel posti giusti

Ecco il codice che sto cercando di creare. Qualsiasi aiuto sarebbe molto apprezzato

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

Il QEF è piuttosto difficile da capire. Spero di poterti aiutare. Il duplice metodo di contouring calcola i dati di 'Hermite' ad ogni punto di incrocio, o in altre parole, in ogni punto creato sul bordo del voxel è nota la normale della superficie. Con un punto e uno normale è possibile calcolare l'equazione di un piano.

Il QEF è la somma dei quadrati delle distanze dal punto interno del voxel a ciascuno dei piani associati al voxel. Di seguito è riportato un pseudo codice per il calcolo del 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); 
} 

L'obiettivo è quindi di scegliere un punto all'interno del voxel che riduce al minimo il QEF. La letteratura suggerisce di utilizzare il processo di Grahm-Schmidt per individuare il punto ottimale, ma questo può essere complesso e può anche portare a punti che si trovano al di fuori del voxel.

Un'altra opzione (hack-ish) è creare una griglia di punti all'interno del voxel e calcolare il QEF per ognuno e scegliere quello con il minimo, più la griglia è fine il più vicino al punto ottimale che verrai ma più lunghi sono i calcoli.

Answer 2

Nella mia attuale implementazione della doppia contornatura im utilizzando un modo molto semplice per risolvere il QEF. Essendo il QEF in sostanza un'approssimazione dei minimi quadrati, ho trovato che il modo più semplice per calcolare il QEF è il calcolo dello pseudoinverso. Questo pseudoinverso può essere calcolato utilizzando qualsiasi libreria algebrica nella tua lingua.

Questo è il codice che sto utilizzando:

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

ingressi della funzione sono i punti di intersezione e normali, e la meanPoint è la media della proposta intersectoin punti.

Riassunto dei calcoli matematici: questa funzione calcola il punto che si trova sull'intersezione di tutti i piani definiti dai punti di intersezione e dalle normali. Poiché questo non ha una soluzione esatta, viene calcolata l'approssimazione dei minimi quadrati, che trova il punto che è 'meno sbagliato'. Inoltre, i punti di intersezione vengono "spostati" in modo che il punto medio diventi l'origine. Questo assicura che quando ci sono più soluzioni al QEF, viene scelta la soluzione più vicina al punto medio.