Программные продукты и системы

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

Существующие алгоритмы компьютерной визуализации линий уровня, применяемые в различных пакетах прикладных программ (MatLab, MatCAD, Beta-Soft и т.д.), не могут использоваться в отечественных программных комплексах, поскольку, во-первых, алгоритмы их реализации недоступны, во-вторых, невозможно их контролировать и управлять ими, в-третьих, нельзя масштабировать область визуализации в пропорциях, адекватных реальным.

Анализ существующих алгоритмов визуализации скалярных полей в виде линий уровня показал, что наиболее пригодным для наших целей является алгоритм обработки полигональной сетки изоповерхности трехмерных скалярных полей Marching Cubes. Данный алгоритм впервые был предложен в 1987 году Вильямом Лоренсеном и Харви Клайном на конференции SIGGRAPH [1]. Он позволяет обрабатывать полигональные сетки изоповерхности трехмерного скалярного поля, называемого еще сеткой вокселей.

Аналогом данного алгоритма для двухмерного скалярного поля стал алгоритм Marching Squares. Описанию применения данного алгоритма для визуализации скалярных двухмерных полей в виде изолиний и встраиванию его в программный комплекс для моделирования температурных полей электронных модулей [2] посвящена данная статья.

Описание алгоритма Marching Squares

Введем следующие понятия и обозначения:

- линия уровня (D) – множество точек P на плоскости, в которых функция U(P) принимает одинаковые значения, U(P)=const;

-      битовая маска – двухмерная матрица Mi,j (i=1, ..., N, j=1, ..., M), в каждой ячейке принимающая значение 0 или 1;

-      матрица значений – двухмерная матрица Ii,j (i=1, ..., N, j=1, ..., M), содержащая исходную регулярную сетку;

-      Imax – максимальное значение в матрице значений, Imax=max(Ii,j);

-      Imin – минимальное значение в матрице значений, Imin=min(Ii,j).

Входными данными для алгоритма является регулярная сетка, представленная матрицей значений размерностью N´M, для которой строится битовая маска по следующему правилу:

Mi,j=

где L - уровень контура, расположенный в диапазоне от Imin до Imax.

Затем битовая маска покрывается контурной сеткой, состоящей из контурных ячеек (рис. 1). Каждая контурная ячейка представляет прямоугольник, охватывающий блок битовой сетки размером 2´2. Отметим также, что размер контурной сетки на одну ячейку меньше в каждом направлении.

При последовательном просмотре контурной сетки вершины каждой контурной ячейки, лежащие ниже линии уровня, помечаются «-», а лежащие выше – «+». Таким образом, для каждой контурной ячейки может быть составлено 16 комбинаций расположения вершин относительно линии уровня (рис. 2).

Для обеспечения быстрого нахождения нужной комбинации контурной ячейки применяется бинарная индексация вершин. Для этого каждую вершину нумеруют двоичными числами 0001 (1), 0010 (2), 0100 (4) и 1000 (8) (рис. 2). Затем, обходя все вершины по часовой стрелке с наиболее значащего бита в левом верхнем углу и применяя побитовую операцию «ИЛИ» и сдвиг влево, находят индекс искомой комбинации K по следующей формуле:

,

где V8, V4, V2, V1 – значения битовой маски при соответствующих вершинах контурной ячейки.

Для полученной комбинации контурной ячейки применяется линейная интерполяция, чтобы найти точное положение контурной линии относительно края ячейки (рис. 2).

Реализация алгоритма Marching Squares на языке Pascal

Алгоритм Marching Squares для заданного уровня контура позволяет строить линии уровня, последовательно обходя все ячейки контурной сетки. На уровне контура могут быть найдены несколько линий уровня. При последовательном обходе всех ячеек контурной сетки среди найденных точек линий уровня невозможно однозначно сопоставить точку и соответствующую ей линию уровня. Чтобы устранить этот недостаток, необходимо выполнить следующие действия:

-      найти контурную ячейку Vi,j, в которой проходит линия уровня Di, просмотрев контурную сетку;

-      определить следующую контурную ячейку, в которой линия Di имеет продолжение;

-      обойти последовательно по часовой стрелке контурные ячейки, в которых проходит линия уровня Di, пока не будет найдена контурная ячейка Vi,j.

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

Приведем описание класса, реализующего алгоритм Marching Squares:

TDouble1DArray = array of Double;

TDouble2DArray = array of array of Double;

TBit2DArray = array of array of Boolean;

TAlgrotihm = class(TObject)

private

      Data: TDouble2DArray; // массив входных данных

X_Size,Y_Size : integer; // размеры массива входных данных

IsolineArray: TIsolineArray; // список линий уровня

Mask: TBit2DArray; // битовая маска

BitField: TBit2DArray;

Section: Double;

Contours : TDouble1DArray;

Direction: TPoint;

N, W, S, E : TPoint;

Origin : TPoint

SectionIndex : integer;

Xnode,Ynode : Integer;

CountSection: Integer;

W, H: Integer;

procedure SearchIsolines;

procedure Check(x,y: integer);

function GetDirection (index:integer;var Isoline: TIsoline): TPoint;

public

procedure CalculateIsoline(const aData: TDoub­le2DArray;

aW,aW:Integer);

end;

Для хранения матрицы значений и ее размеров используются переменные Data, X_Size и Y_Size соответственно. Переменная Mask хранит битовую маску; битовое поле BitField необходимо для отметки ячеек контурной сетки, которые были просмотрены; SectionIndex, Section и CountSec­tion хранят номер уровня контура, его значение и количество уровней соответственно. Переменная Direction определяет направление окончания обхода контурной линии; переменные N, W, S, E означают направления движения при обходе контурной линии; Origin хранит координаты первой найденной контурной ячейки, в которой проходят линии уровня; XNode, YNode – координаты текущей просматриваемой контурной ячейки по осям X и Y соответственно. Переменная IsolineArray определяется типом данных TIsolineArray = array of TList, представляющим собой одномерный массив, каждый элемент которого хранит список линий уровня. Для хранения одной линии уровня введем класс TIsoline, который будет хранить список точек линии уровня, ее значение, направление градиента и обеспечивать возможность добавления новой точки в список точек линии. Направление градиента принимает истинное значение, если градиент направлен во внутреннюю область, образуемую линией уровня. Функция CalculateIsoline вычисляет точки линий уровня. На вход функции подаются массив данных, длина и высота области вывода изображения. Функции SearchIsolines, Check, GetDirection являются вспомогательными и обеспечивают нахождение линий уровня для заданного значения уровня контура. Алгоритм реализован на языке Pascal и встраивается в программный пакет [2].

Методика заливки областей

Для описания настоящей методики введем следующие обозначения. Тип данных Node описывает узел дерева со следующими переменными: Parent – ссылка на узел-предок; Children – список всех узлов-потомков; Isoline – ссылка на линию уровня, связанную с текущим узлом; Region – регион, границей которого являются линия уровня Isoline. Переменная RootList обозначает список, содержащий ссылки на корни root для всех деревьев.

Описанный ранее класс TAlgorithm хранит найденные линии уровня в массиве IsolineArray. Для каждого элемента массива IsolineArrayх[i], хранящего список линий уровня, создается дерево линий уровня, состоящее из единственного узла – корня (root), который добавляется в список всех корней деревьев (RootList). Далее для каждой линии уровня из списка линий уровня создаются узел Node и регион Region, хранящийся внутри этого узла. Созданные узлы вставляются в дерево по вхождении регионов так, что все узлы-потомки, идущие от некоторого узла-предка N вниз, имеют регионы, входящие в регион этого узла-предка N. На рисунке 3 изображено дерево, в котором регион узла N5 лежит внутри региона узла N4, при этом регионы узлов N3, N4 входят в регион узла N0.

Рассмотрим подробнее вставку узла в дерево. Сначала Parent инициализируется значением ссылки на корень дерева Root. Затем для каждого дочернего узла (Child) производятся следующие проверки. Если регион вставляемого узла Node не пересекается ни с одним из регионов дочерних узлов, значит, данный узел является дочерним для Parent. Если регион Node пересекает один из регионов Child, проверяется, является ли регион дочернего узла внутренним для региона узла Node. Если является, значит, регион вставляемого узла содержит внутри себя регион узла Child. Следовательно, узел Node должен располагаться между ним и узлом Parent, являясь новым родительским узлом для Child (рис. 4). Если же регион Node расположен внутри региона Child, значение Parent принимается равным Child, и проверки повторяются.

После упорядочения всех узлов во всех деревьях RootList вычисляется регион для корневого узла в каждом дереве. Для этого совершается обход дерева от его узлов к корню по следующему правилу: регион текущего узла объединяется с регионами дочерних узлов, если градиент линии уровня для текущего узла направлен внутрь региона, или вычитается в противном случае.

Таким образом, после обхода всех узлов корневой узел дерева содержит в себе область, в которой значения скалярной величины больше значения уровня контура (Section) для данного дерева. Чтобы получить регион со значением, равным значению Section, из региона корневого узла списка RootList для индекса I необходимо вычитать регион корня с индексом I+1.

В результате операций вычитания RootList будет содержать список, где каждый элемент, являясь корнем дерева с порядковым номером N, содержит регион областей со значением step N.

В заключение необходимо отметить, что разработанный для визуализации скалярных полей в виде линий уровня алгоритм показал свою работоспособность, он встроен в программный комплекс для моделирования температурных полей электронных модулей [2]. В отличие от существующих алгоритмов, которые используются в зарубежных программных средствах, предназначенных для проектирования технических систем, он позволяет обрабатывать массивы данных с большей скоростью, легко надстраивается, имеет динамически изменяющуюся цветовую палитру по выбору пользователя, а также легко адаптируется к любому программному комплексу.

Автор выражает признательность научному руководителю профессору А.Г. Мадере (НИИСИ РАН, г. Москва) за постановку задачи и обсуждение результатов.

Литература

1. Lorensen W.E., Cline H.E. Marching Cubes: A high resolution 3D surface construction algorithm // Computer Graphics. Vol. 21. №. 4. July 1987.

2. Мадера А.Г., Кандалов П.И. Моделирование трехмерных температурных полей в электронных модулях // Программные продукты и системы. 2010. № 2 (90). С. 29–33.