Двойной контур - Поиск характерной точки, нормали выключены

Я следую этому руководству по реализации Dual Contouring http://www.sandboxie.com/misc/isosurf/isosurfaces.html

Мой источник данных - сетка 16x16x16; Я пересекаю эту решетку снизу вверх, слева направо, почти далеко.

Для каждого индекса моей сетки я создаю структуру куба:

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

Из-за того, как я пересекаю сетку, мне нужно смотреть только на 4 вершины и 3 ребра. На этой картинке вершины 2, 5, 6, 7 соответствуют моим вершинам 0, 1, 2, 3, а ребра 5, 6, 10 соответствуют моим ребрам 0, 1, 2.

Край выглядит так:

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

Затем я проверяю все ребра на предмет смены знака, если есть такой, я нахожу окружающие кубы и получаю характерную точку этих кубов.

Теперь это работает, если я установлю характерную точку в центре куба, то получу блочный вид майнкрафта. Но это не то, что я хочу.

Чтобы найти особенность, я хотел сделать это, как в этом посте: https://gamedev.stackexchange.com/a/83757/49583

По сути, вы начинаете вершину в центре ячейки. Затем вы усредняете все векторы, взятые из вершины в каждую плоскость, перемещаете вершину вдоль полученного результата и повторяете этот шаг фиксированное число раз. Я обнаружил, что перемещение его на ~ 70% по результату стабилизируется за наименьшее количество итераций.

Итак, я получил класс самолета:

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

и функция для получения характерной точки, где я создаю 3 плоскости, по одной на каждое ребро и усредняю ​​расстояние до центра:

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

Но это не работает, вершины повсюду. Где ошибка? Могу ли я на самом деле вычислить нормаль ребра путем усреднения нормали вершин ребра? Я не могу получить плотность в средней точке края, поскольку у меня есть только целочисленная сетка в качестве источника данных ...

Редактировать: я также нашел здесь http://www.mathsisfun.com/algebra/systems-linear-equations-matrices.html, что я могу использовать матрицы для вычисления пересечения 3-х плоскостей, по крайней мере, так я понял, поэтому Я создал этот метод

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

который с этими данными

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

вычисляет пересечение в (2.0, 2.0, 2.0), поэтому я предполагаю, что оно работает правильно. Тем не менее, не верные вершины. Я действительно думаю, что это мои нормальные.

Прежде всего, ваши нормальные показатели должны быть полностью нормальными, если они рассчитываются с помощью различий вперед / назад / центральности. Проблема в том, что вы переместили свою центральную точку в неправильном направлении в вашей функции FeaturePoint, что привело к дальнейшему отклонению от минимума.

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;

Это произошло потому, что ваш код не сходится с точкой и, следовательно, выпрыгивает из окна вашего вокселя. Я не знаю, если код из Кто-то может объяснить двойной контур? предполагалось использовать проекционный подход, при котором точка проецируется на плоскость через:

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

но это тот же метод. Если вы переписываете проекцию в исходный код, это приводит к:

    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;

который можно переписать на:

    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;

и, следовательно, результат в первом коде. Проецируя точку на три неплоские плоскости, она медленно сходится к минимуму, потому что вы минимизируете расстояние от каждой плоскости до вашей точки, как показано на рисунке.

Красные точки обозначают характерную точку, синие линии - нормали, а фиолетовую - точку, спроецированную на плоскость. Вам также не нужно использовать коэффициент 0,7, потому что без него он должен сходиться быстрее. Если вы используете этот метод, будьте осторожны, что алгоритм может не работать, если у вас непересекающиеся плоскости.