# Dual Contouring - Finding the feature point, normals off

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.

• Unity already has a Plane structure defined (see here), which has the methods you gave already defined (except the shortest vector method, which you can add to the Plane structure using C# extension methods). You can use the GetDistanceToPoint method instead of your Distance method. – EvilTak Nov 16 '15 at 10:16
• Thanks for your comment, I replaced my implementation with the Unity implementation and using this function private Vector3 shortestDistanceVector(Plane p, Vector3 point) { return p.GetDistanceToPoint(point) * p.normal; } I also get only random vertices. I suspect my normals are totally off. I also added an edit, where I tried a second method, maybe you can have a look at it and tell me what I did wrong there. – ElDuderino Nov 16 '15 at 10:49
• Can I actually calculate the edge normal by averaging the normal of the edge vertices? - I may be mistaken, but I think I've seen advice elsewhere saying never to interpolate in order to get normals - they just don't interpolate well. Calculate per face, it's safer. Really, you should first be constructing a minimum test case to ensure your normals calculation is correct. Then move on with this. – Engineer Nov 16 '15 at 10:54
• But I only get the faces after I have the normals, I need the normals to create the planes and get the vertices for the faces from them. And as said with my current structure I can index into my data only at the edge vertices. Or what faces are your talking about? – ElDuderino Nov 16 '15 at 11:14
• @ElDuderino Faces like faces (or triangles) of a mesh, but I don't know how you can get that from your data. If you can generate triangles instead of edges, then normal calculation becomes really easy. – EvilTak Nov 16 '15 at 11:29

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.

• Hey, awesome to get an answer after 2 years :) I never found a solution, so I stopped this project, but I'll revisit it with this knowledge and let you know how it went. Have a +1 til then. – ElDuderino May 30 '17 at 19:16
• Awesome! I'm glad I could help you. Let me know if it works for you. – Tim Rolff May 30 '17 at 20:17