Построение изолиний поверхности

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

Задача.

Поверхность – двумерный массив размером nx на ny точек, в каждой из которых содержится Z отметка высоты. Пусть, поверхность хранится в массиве вида:

std::vector<double> surf;

К элементу с номером i, j можно обратиться следующим образом: surf[j*nx + i];

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

Элементарные элементы из которых состоит поверхность – ячейки – прямоугольники, образованные четыремя соседними вершинами. Эти ячейки не подходят для трассировки, т.к. через каждый прямоугольник можно провести две изолинии, притом, несколькими способами.

Треугольник является такой фигурой, через которую можно провести только одну изолинию. Это означает что если плоскость, описываемая уравнением Z = const пересекает какую-либо сторону треугольника, она гарантировано пересечет одну из двух оставшихся сторон. Это легко доказывается.

Поэтому для трассировки изолиний необходимо разбить существующие прямоугольники на треугольники следующим образом:



Подготовка.

Трассировка изолиний выполняется с заданным шагом dz в интервале изменения высот поверхности. Далее приведено описание построения изолинии для одной величины Z.

Первым этапом является проверка всех существующих граней на пересечение с заданной величиной Z.

Каждую грань легко задать одним узлом и направлением. Всего существует три направления: горизонтальное, вертикальное, диагональное. Поэтому, для хранения отмеченных граней понадобится массив такого же размера, как и масив поверхности nx на ny. Каждая ячейка массива будет хранить unsigned char число, являющееся битовым полем. В битовом поле закодированны три направления: 001 – горизонтальное, 010 – вертикальное, 100 – диагональное. Если ячейка массива содержит 0 – значит из этого узла не выходит ни одного ребра, пересекающегося с заданной величиной Z.

std::vector<unsigned char> edges; // bit





// так можно отметить вертикальную грань


if ((z1 < curz) && (z >= curz))


edges[y * nx + x] |= 2;

Итак, в цикле перебираем все точки поверхности и заносим информацию об отмеченных ребрах с специальный массив.



Трассировка.

Для того что бы трассировать линии удобно применить такой контейнер как очередь (FILO – First In Last Out). На каждом шаге цикла из очереди будет доставаться одна грань, рассчитываться координаты точки пересечения этой грани с плоскостью Z = Z0, анализироваться треугольник, которому принадлежит грань и добавляться новые грани в конец очереди.

Элементы хранимые в очереди – грани треугольников – отрезки соединяющие две соседние точки в массиве поверхности. Грань может быть представлена такой структурой данных:

struct Edge


{


int knot1;


int knot2;


};

Каждый узел поверхности однозначно определяется целым числом, значение которого можно вычислить как j*nx+i, где i и j – координаты узла в массиве.
Так как каждая грань принадлежит двум треугольникам, удобно задавать ее координатами начала и конца. Это позволяет определить направление вектора-грани, что дает информацию о треугольнике, к которому принадлежит грань. Будем считать что все треугольники обходятся против часовой стрелки, и что бы узнать третью координату треугольника, достаточно обойти его против часовой стрелки по направлению вектора-грани.

Первый этап

Для начала надо найти хотя бы одну отмеченную грань. Последовательно перебираем массив, где хранятся метки граней и ищем значение не равное 0. Это означает, что из данный узел является началом хотя бы одной грани (горизонтальной, вертикальной, диагональной), пересекающей поверхность Z = Z0;

После того как узел найден, проверяем какая грань, выходящая из него, отмечена и добавляем эту грань в очередь. Причем добавляется прямое и инвертированное направление грани: (knot1, knot2) и (knot2, knot1). Это позволяет начать трассировку линий в двух направлениях.

Втоой этап.

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

а) Достаем грань из очереди и находим третью вершину, вместе с которой она образует треугольник. Третью вершину легко найти зная координаты вершин грани (knot1, knot2).
б) Проверяем две оставшиеся грани реугольника, если какая-либо грань пересекает плоскоссть Z = Z0 (проверяем массив с отмечеными гранями), добавляем ее в очередь.
в) Считаем точные координаты точки пересечения исходной грани треугольника с горизонтальной плоскостью. Это просто сделать, зная значения высот узлов грани Z1, Z2. и координаты точек (x1, y1) (x2, y2)
г) Добавляем полученную координату к текущей изолинии. Убираем метку этой грани из массива отмеченных граней и удяляем грань из очереди.

Повторяем цикл, пока очередь не исчерпается.

Стоит заметить, что самая первая грань в очереди добавляется два раза, в прямом и инвертированном направлении. Это необходимо для того, что бы начать построение изолинии сразу в двух направлениях. Если изолиния замкнутая, к таком способу можно не прибегать, но если изолиния разомкнутая и заканчивается на границах, требуется отследить ее в обе стороны, пока изолиния не дойдет до границы.. Для того что бы знать, в начало или в конец изолинии добавлять новую точку, можно хранить вместе с каждой отмеченной гранью переменную int direction, которая может принимать значения +1 или -1.

struct Edge


{





int knot1;






int knot2;










int direction;



};

Я использовал слово “грань”, которое в общем случае обозначает поверхность объемного многогранница для удобства, т.к. словосочетание “сторона треугольника” кажется мне не ёмким