Dual Contouring - Finding the feature point, normals off

Question

I am following this tutorial to implement Dual Contouring http://www.sandboxie.com/misc/isosurf/isosurfaces.html

My data source is a grid 16x16x16; I traverse this grid bottom to top, left to right, near to far.

For each index of my grid, I create a cube structure :

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

Due to how I traverse the grid, I only need to look at 4 vertice and 3 edges. In this picure, the vertices 2, 5, 6, 7 correspond to my vertices 0, 1, 2, 3, and the edges 5, 6, 10 correspond to my edges 0, 1, 2.

An edge looks like this :

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

I then check all the edges for a sign change, if there is one I find the surrounding cubes and get the feature point of those cubes.

Now that works if I set the feature point to the cube center, I then get the blocky minecraft look. But that's not what I want.

To find the feature point, I wanted to do it as in this post : https://gamedev.stackexchange.com/a/83757/49583

Basically, you start the vertex in the centre of the cell. Then you average all the vectors taken from the vertex to each plane and move the vertex along that resultant, and repeat this step a fixed number of times. I found moving it ~70% along the resultant would stabilize in the least amount of iterations.

So I got a Plane class :

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

and a function to get the feature point, where I create 3 planes, one for each edge and average the distance to the center :

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

But it's not working, the vertices are all over the place. Where is the error? Can I actually calculate the edge normal by averaging the normal of the edge vertices? I cannot get the density at the edge midpoint, as I only have an integer grid as datasource ...

Edit : I also found here http://www.mathsisfun.com/algebra/systems-linear-equations-matrices.html that I can use matrices to compute the intersection of the 3 planes, at least that's how I understood it, so I created this method

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

which with this data

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

calculates an intersection at (2.0, 2.0, 2.0), so I assume it works correct. Still, not the correct vertices. I really think it is my normals.

Answer 1

First of all your normals should be totaly fine if they are calculated through forward- / backward- / central- differences. The problem is, that you moved your center point torwards the wrong direction in your FeaturePoint function which results in going further away from the 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;

This happend because your code does not converge against a point and therfore jump out of your voxel box. I don't know if the code from Can someone explain dual contouring? was intended to use a projection approach where the point gets projected onto the plane through:

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

but it's the same method. If you rewrite the projection into your original code this results in:

    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;

which can be rewritten to:

    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;

and therefore result in the first code. By projecting the point on three non planar planes it slowly converges towards a minimum because you minimize the distance from each plane to your point as show in the picture.

The red dots denote the feature point, the blue lines the normals and the purple point the point projected on the plane. You also don't need to use the factor 0.7 because it should converge faster without it. If you using this method beware, that the algorithm may not be working if you have non intersecting planes.