Dual Contouring - Recherche du point caractéristique, les normales désactivées

Je suis ce tutoriel pour implémenter Dual Contouring http://www.sandboxie.com/misc/isosurf/isosurfaces.html

Ma source de données est une grille 16x16x16; Je traverse cette grille de bas en haut, de gauche à droite, de près à loin.

Pour chaque index de ma grille, je crée une structure cubique:

public Cube(int x, int y, int z, Func<int, int, int, IsoData> d, float isoLevel) {
            this.pos = new Vector3(x,y,z);
            //only create vertices need for edges
            Vector3[] v = new Vector3[4];
            v[0] = new Vector3 (x + 1, y + 1, z);
            v[1] = new Vector3 (x + 1, y, z + 1);
            v[2] = new Vector3 (x + 1, y + 1, z + 1);
            v[3] = new Vector3 (x, y + 1, z + 1);
            //create edges from vertices
            this.edges = new Edge[3];
            edges[0] = new Edge (v[1], v[2], d, isoLevel);
            edges[1] = new Edge (v[2], v[3], d, isoLevel);
            edges[2] = new Edge (v[0], v[2], d, isoLevel);
        }

En raison de la façon dont je traverse la grille, je n'ai qu'à regarder 4 sommets et 3 arêtes. Dans cette image, les sommets 2, 5, 6, 7 correspondent à mes sommets 0, 1, 2, 3 et les bords 5, 6, 10 correspondent à mes bords 0, 1, 2.

Un bord ressemble à ceci:

    public Edge(Vector3 p0, Vector3 p1, Func<int, int, int, IsoData> d, float isoLevel) {
        //get density values for edge vertices, save in vector , d = density function, data.z = isolevel 
        this.data = new Vector3(d ((int)p0.x, (int)p0.y, (int)p0.z).Value, d ((int)p1.x, (int)p1.y, (int)p1.z).Value, isoLevel);
        //get intersection point
        this.mid = LerpByDensity(p0,p1,data);
        //calculate normals by gradient of surface
        Vector3 n0 = new Vector3(d((int)(p0.x+1),   (int)p0.y,      (int)p0.z       ).Value - data.x,
                                 d((int)p0.x,       (int)(p0.y+1),  (int)p0.z       ).Value - data.x,
                                 d((int)p0.x,       (int)p0.y,      (int)(p0.z+1)   ).Value - data.x);

        Vector3 n1 = new Vector3(d((int)(p1.x+1),   (int)p1.y,      (int)p1.z       ).Value - data.y,
                                 d((int)p1.x,       (int)(p1.y+1),  (int)p1.z       ).Value - data.y,
                                 d((int)p1.x,       (int)p1.y,      (int)(p1.z+1)   ).Value - data.y);
        //calculate normal by averaging normal of edge vertices
        this.normal = LerpByDensity(n0,n1,data);
    }

Je vérifie ensuite tous les bords pour un changement de signe, s'il y en a un, je trouve les cubes environnants et j'obtiens le point caractéristique de ces cubes.

Maintenant, cela fonctionne si je définis le point de fonctionnalité au centre du cube, j'obtiens alors le look minecraft en bloc. Mais ce n'est pas ce que je veux.

Pour trouver le point de fonctionnalité, je voulais le faire comme dans cet article: https://gamedev.stackexchange.com/a/83757/49583

Fondamentalement, vous démarrez le sommet au centre de la cellule. Ensuite, vous faites la moyenne de tous les vecteurs pris du sommet à chaque plan et déplacez le sommet le long de cette résultante, et répétez cette étape un nombre fixe de fois. J'ai trouvé que le déplacer ~ 70% le long de la résultante se stabiliserait dans le moins d'itérations.

J'ai donc eu une classe d'avion:

private class Plane {

        public Vector3 normal;
        public float distance;

        public Plane(Vector3 point, Vector3 normal) {
            this.normal = Vector3.Normalize(normal);
            this.distance = -Vector3.Dot(normal,point);
        }

        public float Distance(Vector3 point) {
            return Vector3.Dot(this.normal, point) + this.distance;
        }

        public Vector3 ShortestDistanceVector(Vector3 point) {
            return this.normal * Distance(point);
        }
 }

et une fonction pour obtenir le point caractéristique, où je crée 3 plans, un pour chaque bord et moyenne la distance au centre:

 public Vector3 FeaturePoint {
            get {
                Vector3 c = Center;
 //                 return c; //minecraft style

                Plane p0 = new Plane(edges[0].mid,edges[0].normal);
                Plane p1 = new Plane(edges[1].mid,edges[1].normal);
                Plane p2 = new Plane(edges[2].mid,edges[2].normal);

                int iterations = 5;
                for(int i = 0; i < iterations; i++) {
                    Vector3 v0 = p0.ShortestDistanceVector(c);
                    Vector3 v1 = p1.ShortestDistanceVector(c);
                    Vector3 v2 = p2.ShortestDistanceVector(c);
                    Vector3 avg = (v0+v1+v2)/3;
                    c += avg * 0.7f;
                }

                return c;
            }
        }

Mais ça ne marche pas, les sommets sont partout. Où est l'erreur? Puis-je réellement calculer la normale du bord en faisant la moyenne de la normale des sommets du bord? Je ne peux pas obtenir la densité au milieu du bord, car je n'ai qu'une grille entière comme source de données ...

Edit: j'ai également trouvé ici http://www.mathsisfun.com/algebra/systems-linear-equations-matrices.html que je peux utiliser des matrices pour calculer l'intersection des 3 plans, du moins c'est comme ça que je l'ai compris, donc J'ai créé cette méthode

 public static Vector3 GetIntersection(Plane p0, Plane p1, Plane p2) {              
            Vector3 b = new Vector3(-p0.distance, -p1.distance, -p2.distance);

            Matrix4x4 A = new Matrix4x4 ();
            A.SetRow (0, new Vector4 (p0.normal.x, p0.normal.y, p0.normal.z, 0));
            A.SetRow (1, new Vector4 (p1.normal.x, p1.normal.y, p1.normal.z, 0));
            A.SetRow (2, new Vector4 (p2.normal.x, p2.normal.y, p2.normal.z, 0));
            A.SetRow (3, new Vector4 (0, 0, 0, 1));

            Matrix4x4 Ainv = Matrix4x4.Inverse(A);

            Vector3 result = Ainv * b;
            return result;
        }

qui avec ces données

        Plane p0 = new Plane (new Vector3 (2, 0, 0), new Vector3 (1, 0, 0));
        Plane p1 = new Plane (new Vector3 (0, 2, 0), new Vector3 (0, 1, 0));
        Plane p2 = new Plane (new Vector3 (0, 0, 2), new Vector3 (0, 0, 1));

        Vector3 cq = Plane.GetIntersection (p0, p1, p2);

calcule une intersection à (2.0, 2.0, 2.0), donc je suppose que cela fonctionne correctement. Pourtant, pas les bons sommets. Je pense vraiment que ce sont mes normales.

Tout d'abord, vos normales devraient être parfaitement correctes si elles sont calculées par des différences avant / arrière / centrales. Le problème est que vous avez déplacé votre point central vers la mauvaise direction dans votre fonction FeaturePoint, ce qui vous éloigne davantage du minimum.

Vector3 c = Center;
Plane p0 = new Plane(edges[0].mid,edges[0].normal);
Plane p1 = new Plane(edges[1].mid,edges[1].normal);
Plane p2 = new Plane(edges[2].mid,edges[2].normal);

int iterations = 5;
for(int i = 0; i < iterations; i++) {
    Vector3 v0 = p0.GetDistanceToPoint(c) * edges[0].normal;
    Vector3 v1 = p1.GetDistanceToPoint(c) * edges[1].normal;
    Vector3 v2 = p2.GetDistanceToPoint(c) * edges[2].normal;
    Vector3 avg = (v0+v1+v2)/3;
    c -= avg * 0.7f; // Error was here!
}
return c;

Cela s'est produit parce que votre code ne converge pas contre un point et saute donc hors de votre boîte de voxel. Je ne sais pas si le code de Est-ce que quelqu'un peut expliquer le double contour? était destiné à utiliser une approche de projection où le point est projeté sur l'avion à travers:

distance = Vector3.Dot(point - origin, normal);
projectedPoint = point - distance * normal;

mais c'est la même méthode. Si vous réécrivez la projection dans votre code d'origine, cela se traduit par:

    Vector3 v0 = c - p0.GetDistanceToPoint(c) * edges[0].normal;
    Vector3 v1 = c - p1.GetDistanceToPoint(c) * edges[1].normal;
    Vector3 v2 = c - p2.GetDistanceToPoint(c) * edges[2].normal;
    c = (v0+v1+v2)/3;

qui peut être réécrit en:

    Vector3 v0 = p0.GetDistanceToPoint(c) * edges[0].normal;
    Vector3 v1 = p1.GetDistanceToPoint(c) * edges[1].normal;
    Vector3 v2 = p2.GetDistanceToPoint(c) * edges[2].normal;
    c = c - (v0+v1+v2)/3;

et donc aboutir au premier code. En projetant le point sur trois plans non plans, il converge lentement vers un minimum car vous minimisez la distance de chaque plan à votre point comme indiqué dans l'image.

Les points rouges indiquent le point caractéristique, les lignes bleues les normales et le point violet le point projeté sur le plan. Vous n'avez pas non plus besoin d'utiliser le facteur 0,7 car il devrait converger plus rapidement sans lui. Si vous utilisez cette méthode, méfiez-vous, que l'algorithme peut ne pas fonctionner si vous avez des plans sans intersection.