Построение на GPU в масштабе реального времени адаптивной модели рельефа Земли на основе эллипсоида вращения

Построение на GPU в масштабе реального времени адаптивной модели рельефа Земли на основе эллипсоида вращения

П.Ю. Тимохин, М.В. Михайлюк, А.В. Мальцев

Аннотация—В статье рассматривается задача моделирования в масштабе реального времени рельефа Земли на основе детализированных глобальных карт высот, заданных относительно эллипсоида вращения WGS-84. Предлагается технология адаптивной тесселяции треугольных патчей на графическом процессоре (GPU), которая обеспечивает построение в реальном времени сложных полигональных моделей рельефа Земли. Технология включает в себя этап разбиения видимого эллипсоида на крупные патчи (низкого уровня детализации) и этап их тесселяции до треугольников модели рельефа, выполняемой параллельно и независимо друг от друга на ядрах GPU. В работе предложены новые, распределенные методы и алгоритмы отбора патчей, необходимых для визуализации текущего кадра, повышения их уровня детализации в соответствии с картой высот и разрешением экрана и преобразования в полигоны рельефа. Новизна работы состоит в том, что тесселяция треугольных патчей эллипсоида выполняется на основе оригинальной схемы, которая позволяет эффективно локализовывать участки рельефа, требующие повышенной полигональной детализации. Это существенно уменьшает время построения модели рельефа и повышает качество создаваемых изображений.

Разработанные технология, методы и алгоритмы были реализованы в программных модулях и апробированы в системе визуализации виртуальной поверхности Земли. Моделирование и визуализация рельефа Земли выполнялись с детализированной текстурой подстилающей земной поверхности и расчетом освещенности с учетом атмосферы. Апробация проводилась в составе космической виртуальной сцены, включающей подробную полигональную модель Международной космической станции (МКС). Полученные результаты подтвердили адекватность предложенного решения поставленной задаче и его применимость для построения космических видеотренажерных комплексов, систем виртуального окружения, виртуальных лабораторий и др.

Ключевые слова—визуализация, Земля, рельеф, GPU, виртуальная модель, эллипсоид вращения, реальное время, тесселяция, распределенные вычисления.

Статья получена 20 сентября 2019.

Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 17-07-00243.

П.Ю. Тимохин, ФГУ ФНЦ НИИСИ РАН, старший научный сотрудник (e-mail: tpyu@yandex.ru).

М.В. Михайлюк, доктор физико-математических наук, профессор, ФГУ ФНЦ НИИСИ РАН, главный научный сотрудник (e-mail: mix@niisi.ras.ru).

А.В. Мальцев, кандидат физико-математических наук, ФГУ ФНЦ НИИСИ РАН, ведущий научный сотрудник (e-mail: avmaltcev@mail.ru).

I. Введение

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

Одной из важных задач космических видеотренажеров является моделирование и визуализация рельефа земной поверхности с учетом ее несферичности. Особенно это востребовано в тренировках по проведению визуально-инструментальных наблюдений Земли из космоса [4], которые включают в себя использование телекамер и фотоаппаратов с многократным увеличением [5].

Эффективным путем является моделирование рельефа Земли на основе детализированных глобальных сеток (карт) высот, заданных относительно эллипсоида вращения в международном стандарте WGS-84 [6]. Можно выделить два главных подхода к решению этой задачи.

Первый подход основан на испускании лучей из положения наблюдателя через пикселы экрана до пересечения с моделью земной поверхности [7, 8]. Преимуществом данного подхода является высокое качество синтезируемых изображений земной поверхности. Недостаток заключается в существенном падении скорости синтеза изображений с ростом числа обрабатываемых пикселов (просчитываемых лучей), что препятствует выполнению визуализации в масштабе реального времени на экранах со сверхвысоким разрешением. Недавнее появление аппаратной поддержки трассировки лучей у графических карт NVidia (линейка GeForce RTX) призвано улучшить данную ситуацию.

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

процессора (GPU). Однако, при этом качество визуализации модели рельефа ограничивается числом треугольников, которое видеокарта может визуализировать в масштабе реального времени (оно существенно меньше числа, необходимого для полного покрытия детализированной карты высот). Эффективным путем обхода этого ограничения является видозависимая адаптивная триангуляция карты высот, при которой построение треугольников выполняется только для видимых участков карты высот, а участки рельефа, на которых перепады высот не заметны на экране, представляются меньшим числом треугольников, чем остальные [9-11]. При реализации данного подхода на детализированных картах высот полигональная сложность (число треугольников) получаемой модели рельефа может существенно изменяться от кадра к кадру. В этой связи возникает задача разработки эффективных методов и алгоритмов построения в масштабе реального времени сложной адаптивной триангулированной модели земного рельефа, основанных на распределенных вычислениях и распараллеливании графических расчетов на современных многоядерных GPU.

В исследовании [12] построение полигональной модели рельефа выполняется на основе сферы, а карта высот записывается в видеопамять в виде концентрических текстур с убывающей детализацией. В работе [13] предлагается метод построения на GPU адаптивной полигональной модели рельефа на основе четырех типов предрассчитанных геометрических шаблонов. Авторы исследования [14] распараллелили построение треугольников адаптивной модели рельефа с помощью геометрического шейдера. Несмотря на перенос ряда графических вычислений на GPU, в приведенных работах ключевые расчеты, связанные с управлением тесселяцией (разбиением на треугольники) реализуются последовательно, что накладывает серьезные ограничения на размеры карт высот, которые могут быть обработаны в масштабе реального времени.

Добавление в архитектуру GPU программируемой тесселяции [15] открыло возможность управлять разбиением на треугольники параллельно на ядрах GPU (с помощью шейдеров управления и вычисления тесселяции). В работе [16] предложен метод адаптивной тесселяции на GPU рельефа с помощью параметрических четырехугольников (патчей-квадов) одинакового размера. В работе [17] для построения модели рельефа используется древовидная структура из патчей-квадов различных уровней детализации. Недостатком таких решений является образование разрывов геометрии (cracks) в местах стыковки смежных треугольников типа “сторона-вершина”. Для скрытия разрывов в [16] строятся дополнительные треугольники-заслонки, а в работе [18] предлагается использовать динамические сшивающие стрипы (ленты-треугольников). Другим недостатком патчей-квадов является склонность к образованию избыточной триангуляции в модели рельефа (если в крупный патч-квад попадает маленький участок карты высот, требующий высокой полиго-

нальной детализации, то весь патч-квад должен быть тесселирован с такой детализацией).

В данной работе предлагаются новые методы и алгоритмы адаптивной тесселяции на GPU рельефа Земли на основе эллипсоида вращения, которые обеспечивают построение сложных полигональных моделей рельефа Земли в масштабе реального времени. В отличие от моделирования рельефа на основе сферы, восстановление высот рельефа Земли относительно эллипсоида вращения является более точным, однако требует дополнительного объема вычислений, связанных с несовпадением нормали к поверхности эллипсоида и радиус-вектора. В данной работе эти вычисления выполняются параллельно и независимо друг от друга на ядрах GPU с помощью технологии шейдерной тесселяции. Новизна предлагаемого решения состоит в том, что тесселяция эллипосоида по карте высот впервые выполняется с помощью оригинальной схемы разбиения патчей-треугольников [19], которая в отличие от патчей-квадов позволяет эффективно локализовывать участки рельефа с повышенной детализацией, избегая принудительной тесселяции смежных областей. Еще одним преимуществом предложенной схемы является безразрывная стыковка патчей-треугольников, которая выполняется

автоматически во ходе тесселяции при переходе от одного уровня детализации к другому. В разделе II приводится математическая модель восстановления координат точек рельефа Земли по карте высот, заданной относительно эллипсоида вращения. В разделе III предлагается метод и алгоритм построения полигональной модели эллипсоида, используемой в качестве основы для тееселяции. В разделе IV описываются разработанные методы и алгоритмы адаптивной тесселяции полигональной модели рельефа Земли. Предлагаемые методы и алгоритмы реализованы на языках С++ и GLSL (шейдерная модель четвертого поколения) с использованием графической библиотеки OpenGL. В разделе V приводятся полученные результаты.

II. Моделирование рельефа Земли на основе карты ВЫСОТ

Рассмотрим задачу восстановления координат точек рельефа Земли по карте высот, заданной относительно эллипсоида вращения WGS-84 [6].

Определим пространственную геоцентрическую систему OrXr Yr Zr координат, начало которой

совпадает с центром эллипсоида, ось OrZr - направлена на северный полюс Земли, ось OrXr - в точку пересечения начального нулевого меридиана (в системе WGS-84 это Опорный меридиан [6]) с экватором, а ось OrYr дополняет систему до правой (рис. 1). В этой системе будут определяться координаты точек рельефа Земли.

Рассмотрим произвольную точку Ррел рельефа земной поверхности. Проведем через эту точку нормаль n к

поверхности эллипсоида (см. рис. 1). Обозначим через Рэ точку пересечения нормали n с поверхностью эллипсоида, через Q - точку пересечения меридиана, проходящего через Рэ, с экватором, а через h - длину отрезка РэРрел (высота точки Ррел над поверхностью эллипсоида).

n

Рис. 1. Высота h точки рельефа Земли относительно эллипсоида вращения

Значение высоты h записано в карте высот рельефа в

географической системе координат (/ ,ф ), в

которой /ор £ [-п/2,п/2] - угол между плоскостью экватора и нормалью n (географическая широта), а <р £[-п,П) - угол между осью 0ГХГ и вектором

OrQ (географическая долгота). Угол /геогр

положительный, если точка Рэ расположена в северном полушарии, и отрицательный, если Рэ - в южном полушарии. Угол фгеогр положителен, если Рэ находится в восточном полушарии, и отрицателен, если Рэ - в западном полушарии (см. рис. 1).

Карта высот представляет собой двухмерную текстуру, в которой тексел с индексами {0,0) хранит

значение высоты рельефа для /геогр =~П2, Ргеогр.=~л

(см. рис. 2). Для удобства работы с картой высот введем дополнительную систему координат (/’геогр,ф'геогр), в

которой /' £ [0, п] - угол между осью -0ГZr и

нормалью n, а ф’ор £ [0, 2п) - угол между осью -0ГХГ и вектором OrQ. Угол /[тгр откладывается против часовой стрелки, если смотреть из конца вектора n*OrZr, а угол р' - если смотреть из конца оси

0ГZГ на плоскость экватора. Система (/'геогр,ф'геогр) связана с системой /угеогр,фгеогр) следующим образом:

/' = / +П2, ф = ф +я.

т геогр. т геогр. I ’ т геогр. т геогр. v /

Согласно [6] координаты (хэ, уэ, гэ) рассматриваемой

точки Рэ в системе 0ГХГУГZr можно вычислить как

(х N cos /геогр. Cos Фгеогр.

Уэ = N cos /геогр. sin Фгеогр.

К х у чN {1 - e" ) sin Фгеогр. у

где N = aф - e2 sin2 /геогр , e = ^ 1 -b2/а2 - первый

эксцентриситет эллипсоида, a и b - большая и малая полуоси эллипсоида вращения.

t

У < //еогр. ^ П k 2

фгеогр.

-П я П

{0,0) П s

2

Рис. 2. Системы (/ееогр ,Фееогр) и (s, t) координат карты высот Земли

Запишем координаты точки Рэ через углы /

Фгееггр. , выразив из (1) углы /геогр. , Фгеогр.

(2):

^N ^Ф/'геегр. -П^Фф'геегр. - П

N COS

N

Р =

(х>л

Уэ =

К Zэ j

геогр. и

и подставив их в

- N sin/'геогр. cos ф'геогр^ -N sin /'геогр. sin ф'еогр.

-N(1 - e2)sin ф'

где N = aV1 - e" sin2 {/'геогр -П2) =

= а/a - e2 c°S2 /'геогр. .

Согласно [20] нормаль n к поверхности эллипсоида в точке Рэ будет совпадать с градиентом функции

2 2 2

F {х, y, z) = х +2У + ^2-1 (из уравнения эллипсоида) в

этой точке. Исходя из этого, запишем координаты (пх, n , nz) нормали n в системе 0ГХГYrZГ:

( dF (хэ , Уэ , zэ )> ( 2хэ ^ ( \

дх а2 1

dF { хэ , Уэ , zэ ) 2 Уэ 2 1

дУ а2 а а

dF { хэ , Уэ , zэ ) 2 Z а2

dz у К b2 у К b j

Рэ.

Обозначим через ne нормализованный вектор нормали n. Тогда координаты рассматриваемой точки Ррел рельефа Земли в системе 0rXrYrZr можно найти, сместив точку Рэ вдоль вектора ne на высоту h :

Ppe,= Рэ + H .

Построение полигональной модели рельефа Земли будет выполняться в два этапа. На первом этапе (на стадии подготовки данных для визуализации) построим полигональную модель эллипсоида (ПМЭ), состоящую из крупных полигонов. На втором этапе (на стадии визуализации) на основе полигонов ПМЭ, необходимых для визуализации текущего кадра, и карты высот выполним построение полигональной модели рельефа Земли. Рассмотрим эти этапы.

III. Метод построения ПМЭ

Пусть имеется карта высот рельефа размера 2т+1 х 2т текселов, где т > 8 . Чтобы построить ПМЭ, выполним следующие шаги:

а) Выберем на эллипсоиде (равномерно по углам

Угеогр. и Ргеогр.) Ппар= 2^ параллелей и ^ер. = Т- — 1

меридианов (единица вычитается во избежание повторного учета меридиана ргеогр = —п).

б) Зададим координаты (хэ, уэ, z3) позиций и

текстурные координаты

или (1, t) в зависимости от того, какому треугольнику

она будет принадлежать (см. следующий шаг). Поэтому для каждого такого пересечения зададим пару вершин с одинаковыми позициями, но разными текстурными координатами. Таким же образом зададим пересечения меридианов с Северным и Южным полюсами. Ввиду этого общее число вершин в ПМЭ будет равно n = n х (n +1) .

верш. пар. х мер. '

Пронумеруем полученные вершины как показано на рисунке 3.

\ пар. / V мер. /

n

верш.

— 1

Рис. 3. Нумерация вершин ПМЭ

в) Соединим полученные вершины в треугольники, как показано на рисунке 4. Из рисунка видно, что на Северном и Южном полюсах между соседними меридианами треугольники являются одиночными, а не парными, как в остальных случаях. Объединим каждый одиночный треугольник Северного полюса с противоположным ему одиночным треугольником Южного полюса в пару. Тогда число всех пар треугольников в ПМЭ будет равно

Пп. = ПЖрХ (nnaP.— 2).

Рис. 4. Полигональная модель эллипсоида

Пронумеруем полученные пары треугольников, как показано на рисунке 5, и зададим индексы вершин каждого треугольника. Число nUHd всех индексов будет равно 6nn .

Для выполнения описанных шагов создадим одномерные массивы Мпоз позиций (хэ, уэ, гэ) и Мтекс. текстурных координат (s, t) длиной neeрш , а также одномерный массив Мшд индексов вершин треугольников длиной nUHd. .

.

Пп. — 1

мер.

— 1

Рис. 5. Нумерация пар треугольников ПМЭ

Построение полигональной модели эллипсоида реализует Алгоритм 1.

1. Запишем в массивы Мтекс и Мпоз текстурные координаты и позиции (в Or Xr Yr Z Г) вершин ПМЭ: Вычислим шаги As меридиана и At параллели по координатным осям s и t (см. рис. 2):

As = V^ер. , At = 1 nnaP.— 1.

Цикл по i от 0 до nnae — 1 Цикл по j от 0 до nMep

Вычислим номер v (i, j)-ой вершины (см. рис. 3) и ее текстурные координаты (s, t):

V = i * (Пмер. + 1) + j , S = J'AS , t = iAt .

Если i = 0, то:

Mиим = (0,0,-b), Mme«Jv] = (s,0). Перейдем к след. j .

В противном случае, если i = ппар -1, то:

Mиим = (0,0, b), M_.[ V] = (s,1).

Перейдем к след. j .

В противном случае, если j = пмер , то:

Вычислим номер V (/',0)-й вершины:

V'= * * (Пмер+ 1).

Mn03[v] = Mn03 [v'], MmeKC[v] = (1,M_[v'].t). Перейдем к след. i .

В противном случае:

M [v] = (s, t), / = nt, d = 2ns.

Вычислим N(/'геогр) по (4).

Вычислим позицию вершины эллипсоида

Рэ (/'геогP.MеогP., N(/геогр )) согласно (3).

Mp0s [V] = Рэ .

Конец цикла.

Конец цикла.

2. Запишем в массив MUHd индексы треугольников ПМЭ:

Запишем п’ = п +1, п’ = п -1, п” = п - 2.

мер. мер. ’ пар. пар. ’ пар. пар.

Цикл по i от 0 до ппар. - 3 Цикл по j от 0 до пмер. -1 Вычислим номер t пары треугольников, соответствующей (i, j) -ой вершине: t = т + j .

Запишем i' = i +1, j' = j +1.

Если i = 0, то:

// для треугольника Южного полюса

M„Hd.[6t ] = j , MUHd[61 +1] = п'жр+ j ',

Mlmd[6t + 2] = <ер.+ j .

// для треугольника Северного полюса

MuA6t + 3] = п"парп'мер. + j ,

M„Hd.[6t + 4] = п':арп'мер + j ,

M„H,[6t + 5] = п'паРп'меР+ j .

Перейдем к след. j .

В противном случае:

MUHd.[6t] = ™'меР. + j ,

M^ t + 1] = j ' ,

M„Hd.[6t + 2] = Гп'мер+ j,

M„Hd.[6t + 3] = iWMep + j .

Mlmd[6t + 4] = Mua[6t + 2].

Mlmd[6t + 5] = Mua[6t +1].

Конец цикла.

Конец цикла.

3. Запишем массивы Mnm, MmeKC и MUHd в видеопамять с помощью технологии вершинных буферов.

Алгоритм 1. Построение ПМЭ.

В результате работы данного алгоритма получим в видеопамяти полигональную модель эллипсоида, вершины которой выровнены по узловым точкам карты высот. Далее будем выполнять построение полигональной модели рельефа (ПМР) Земли на основе полигонов ПМЭ, необходимых для визуализации текущего кадра, и карты высот. Рассмотрим это подробнее.

IV. Метод построения ПМР Земли

Чтобы построить полигональную модель рельефа Земли, в каждом кадре визуализации треугольники ПМЭ обрабатываются на запрограммированном

определенным образом (см. ниже) графическом конвейере GPU.

В результате одного прохода через конвейер:

- каждый треугольник ПМЭ, который не нужен для визуализации текущего кадра, будет удален;

- каждый оставшийся треугольник ПМЭ, имеющий размеры, достаточно маленькие для качественной аппроксимации рельефа, будет преобразован в участок ПМР;

- каждый непреобразованный треугольник ПМЭ будет разбит на треугольники ПМЭ меньшего размера (от 4 до 6 треугольников, см. раздел 3.2), как показано на рисунке 6, которые будут отправлены на новый проход через графический конвейер.

Такие проходы выполняются в цикле до тех пор, пока не будут получены все участки ПМР, необходимые для визуализации текущего кадра.

Для реализации описанной схемы построения ПМР предлагается использовать программируемую тесселя-цию [21]. В рамках этой технологии каждый треугольник ПМЭ, поступающий на графический конвейер, будет автоматически снабжен дополнительными параметрами тесселяции [21, 22], управляя которыми, мы будем либо удалять треугольник ПМЭ, либо разбивать его на треугольники по аппаратно «зашитому» алгоритму, либо оставлять без разбиения. Такой параметризованный треугольник ПМЭ назовем патчем-треугольииком.

Рассмотрим некоторый i -ый проход графического конвейера. В рамках этого прохода обработка каждого патча-треугольника будет производиться в отдельном потоке GPU с помощью разработанного комплекса шейдерных программ: шейдера управления тесселяцией, шейдера вычисления тесселяции и геометрического шейдера (см. рис. 7). Рассмотрим их более подробно.

A. Управление тесселяцией Данный шейдер работает одновременно на всех ядрах графического процессора и на каждом ядре обрабатывает свой патч-треугольник. Рассмотрим некоторый патч-треугольник  ABC, поступивший на вход шейдера управления тесселяцией (см. рис. 8). Координаты вершин  ABC заданы в объектной системе координат (Object Coordinate System, OCS [23]) полигональной модели рельефа (совпадает с системой ОГ ХГ Yr Z Г). Введем следующие обозначения:

- nA, nB , nC - нормали к поверхности эллипсоида в вершинах  ABC;

- R - участок поверхности рельефа Земли, ограниченный плоскостями {nA,nB}, {nB,nC}, {nC,nA};

-  A’B’C’ - треугольник, вершины которого лежат в угловых точках участка R ;

- hA,, hB,, hC, - высоты вершин  A'B'C' над поверхностью эллипсоида;

- d+ - наибольшее из расстояний от плоскости ABC до точек участка R рельефа (в направлении от эллипсоида); d_ - аналогичное расстояние в обратном направлении;

-  A+ B+ C+ и  A-B'_C'_ - треугольники, вершины которых лежат на прямых AA , BB , CC на расстояниях d+ и d_ от плоскости ABC'.

Рис. 7. Схема построения участка ПМР

Обработка  ABC шейдером управления тесселяцией включает в себя следующие шаги:

1. Вычисление координат вершин пятигранника ABCA+ B+ C+ (см. рис. 8), ограничивающего участок R рельефа.

2. Проверка нахождения пятигранника ABCA+ B+ C+ в поле зрения (в пирамиде видимости) виртуальной камеры.

3. Проверка видимости пятигранника ABCA+ B+ C+ по

горизонту эллипсоида.

4. Расчет параметров тесселяции  ABC.

С шагами 1 и 2 можно ознакомиться в работе [19]. Рассмотрим более подробно шаги 3 и 4.

A+

Рис. 8. Ограничивающий пятигранник участка R рельефа

Проверка видимости по горизонту. Обозначим через Ре - позицию камеры в системе OCS в рассматриваемой сцене на некотором заданном расстоянии от центра ОГ эллипсоида, а через q - вектор PeOr. Проведем из Ре произвольную касательную к

эллипсоиду и обозначим через а - угол между q и этой касательной. Будем называть а углом горизонта эллипсоида. Рассмотрим предельные случаи: а) угол ае горизонта на экваторе и б) угол ар горизонта на полюсе

эллипсоида (см. рис. 9). Из рисунка видно, что

ар<а<ае, где ар,ае е (0,л/2].

Рис. 9. Углы горизонта на экваторе и полюсе эллипсоида

Введем флаг bvjs видимости пятигранника

ABCA+ B+ C+ в пределах наименьшего угла ар горизонта (bvis = 1 означает, что пятигранник виден). Тогда, вычис-

лить значение флага bvis можно, проверив, что  ABC не лежит полностью на обратной стороне эллипсоида, а  A'+B'+C'+ - полностью в конусе с осью q симметрии и

углом 2ар раствора. Это реализует Алгоритм 2.

1. Инициализируем флаг bvis = 0.

2. Проверим, что  ABC не лежит полностью на обратной стороне эллипсоида:

Запишем r = APe, r2 = BPe, r3 = CPe.

Проверим знаки cos(rt, nA), cos(r2, nB) и cos(r3, nc): Если (r3, nA) > 0 или (r2, nB) > 0, или (r3, nC) > 0,

то h„ = 1.

3. Если bvis = 0, то проверим, что  A+ B+ C+ не лежит полностью в конусе 2аp :

Запишем p3 = PeA+ , Р2 = PeB+ , Рз = PeC+ .

Введем углы Д = (p, q), в2 = (Р2, q), А = (Рз, q),

где Рх, в2 , в 6 (0,0.

Проверим знаки cos Д , cos в2 и cos Д3:

Запишем g1 = (p, q), g2 = (Р2, q), g3 = (Рз, q).

Если gj < 0 или g2 < 0 , или g3 < 0 , то bvis = 1.

Если bvis = 0, то проверим отношения cos2 apjcos2 Д , cos2 ар! cos2 Д2 и cos2 apj cos2 Д3 :

Запишем общий сомножитель этих отношений:

d = (q, q) - b2.

Если d(p2, p2)/gf > 1 или d(Р2, Р2)/g2 > 1, или

d(Рз , Рз Vg32 > 1 , то bvs = 1.

Алгоритм 2. Проверка видимости по горизонту.

Отметим, что п. 3 Алгоритма 2 также обеспечивает корректный учет сложного неявного случая видимого рельефа, когда основание возвышенности рельефа скрыто за горизонтом, а вершины находятся выше линии горизонта.

Полученное в результате выполнения Алгоритма 2 значение флага bvis мы будем использовать при расчете значений параметров тесселяции  ABC на следующем шаге.

Расчет параметров тесселяции. Согласно технологии шейдерной тесселяции тесселируемый треугольник (в нашем случае это  ABC) считается равносторонним, а его разбиение на треугольники выполняется на основе разбиения его сторон на отрезки и построения концентрических равносторонних треугольников [22]. Управление этими геометрическими операциями осуществляется с помощью параметров тес-селяции - трех внешних и одного внутреннего уровней тесселяции.

Обозначим через loul A первый внешний уровень тесселяции  ABC - число равных отрезков, на которые будет разбита сторона BC , а через lout B и lout C - второй и третий внешние уровни тесселяции (они определены аналогично уровню 1ШЛ). Также обозначим через ln -

внутренний уровень тесселяции  ABC - число равных шагов на каждой стороне  ABC , через которые по перпендикулярам (см. рис. 10) будут построены вершины концентрических треугольников. На рисунке 10 изображены два концентрических треугольника для l n = 4 (один из них вырожден в точку).

Расчет уровней lutA,...,l n тесселяции выполняется, если bvjs = 1 (т.е. пятигранник ABCA+ B+ C+ виден). В противном случае значения всех уровней тесселяции обнуляются, что удалит  ABC из графического конвейера.

B

Рис. 10. Принцип построения концентрических треугольников

Чтобы вычислить уровни lutA,...,l n , используется

система координат стандартного (нормализованного) объема видимости (Normalized Device Coordinate System, NDCS [23]). Рассмотрим ребро A-A+ пятигранника

ABCA+ B+ C'+. Обозначим через PA, и PA, координаты его вершин в системе OCS, а через P* и Р*А, -координаты в системе NDCS. Построим проекцию r вектора PA, PA. на плоскость XNDCSYNDCS (см. рис. 11).

Рис. 11. Нормализованный объем видимости

Обозначим через qA, A, квадрат длины r, а через qB B и qcc, - квадраты аналогичных длин для ребер B—B+ и

C—C+. Введем некоторое пороговое значение qthres для qA, A, , qB в и qc, С, , при котором замена пятигранника ABCA+ B+ C+ на  A'B'C' будет не заметна на экране для наблюдателя:

qhres = (Прх max(2/ Wscr ,2I hscr ))2 ,

где wscr и hscr - ширина и высота экрана, в пикселах, а npx - пороговое число пикселов на экране (1-2 пиксела). Вычислим значения уровней loutA,...,lin тесселяции на основе сравнения qA,A, , qB, B, , qc, c, с qthres, используя разработанный Алгоритм 3.

1. Передадим в шейдер управления тесселяцией порог qthres и произведение MPMV матрицы перспективной проекции на модельно-видовую матрицу [23].

2. Вычислим значение qA,A, :

Запишем координаты PA, :

PA_ = MPMVPA

PA- = P s/p* w, pA JpA w, pA JpA w ,1).

Запишем PA аналогично PA

r = (PA+ ,x - PA x, PA y— PA y). qy a+ = lrl =(rr).

3. Вычислим qB, B, и qc, c, аналогично qA, A, .

4. Запишем флаги bA, A,, bB

bc-c+

bA- A = (4a_ A — q,hres ) , bB- B+ = (Чв'_ B+ — q,hres ) . bc-c+ = (qc-c+ — qthres) .

5. Запишем значения уровней loutA.

...,l

lout, A 2 bB

lout,B 2 bA- Abc-c+

lout ,C = 2 — bA- AbB- B+ , lin = 4 — 3bA- AbB- B+ bc-c+ .

6. Запишем флаг btsss тесселяции  ABC означает, что тесселяции  ABC не будет):

btess = 0 , если lout,A = lout,B = lout,C = 1 .

Алгоритм 3. Расчет параметров тесселяции.

( btess = 0

В результате выполнения Алгоритма 3 мы получим значения внешних уровней тесселяции ( lout, A , lout,B , loutC) равные 1 или 2, а внутреннего уровня - 1 или 4.

Отметим, что введенная в пп. 4 и 5 алгоритма замена условных переходов на вычисляемые выражения позволяет избавиться от лишних ветвлений, которые снижают эффективность выполнения алгоритмов на GPU.

Также в дальнейшем (при вычислении тесселяции  ABC) нам понадобится информация о том, какая сторона у  ABC является наибольшей в системе OCS. В общем случае ею может быть любая сторона  ABC, т.к. начало обхода его вершин не фиксировано. Чтобы задать наибольшую сторону, введем флаги bAB , bBC и bAc сторон  ABC, где только у наибольшей стороны

флаг будет равен 1, а остальные флаги равны 0. Расчет флагов bAB, bBC, bAc выполняется в системе OCS следующим образом:

1. Запишем квадраты qAB , qBC и qAC длин сторон

AB,BC и AC.

2. Запишем флаги bl, b2, b3:

b1 = (BaE > qAC ) , b2 = (BaE > qBC ),

b3 = (BbC > qAC ) .

3. Запишем флаги bAB , bBC, bAc :

bAB = b1b2 , bBC = b3 (1 — bAB ) .

byc = (1 — bAB )(1 — bBC).

После выполнения шейдера управления тесселяцией  ABC вместе с рассчитанными значениями lutA,...,lin

передается в генератор примитивов. Это

фиксированный этап графического конвейера с аппаратно «зашитым» алгоритмом. В соответствии с нашими значениями loulA,...,lin генератор примитивов

либо оставит  ABC без изменения (если все уровни тесселяции равны 1), либо разобьёт его на треугольники. Во втором случае внутри  ABC будет создан концентрический треугольник с вершиной в центре, а пространство внутри и снаружи этого треугольника будет заполнено более мелкими треугольниками (их число зависиг от значений lou,*, lou,,в, lou,,c ). при этом на сторонах  ABC могут быть добавлены

дополнительные вершины. В примере на рисунке 12 добавлены две вершины M и K, однако в данной работе рассматривается также возможная вершина J на стороне AB .

B

Рис. 12. Разбиение  ABC на треугольники в барицентрической системе координат при

lout,A = lout,в = 2 , lout,C = 1

Вершины сгенерированных треугольников заданы барицентрическими координатами в точечном базисе A , B , C (эти координаты приведены в работе [19]). Вместе с барицентрическими координатами сгенерированные вершины передаются на следующий этап графического конвейера - шейдер вычисления тесселяции.

B. Вычисление тесселяции

На данном этапе графического конвейера координаты вершин сгенерированных треугольников ПМЭ преобразуются из барицентрической системы координат в систему OCS, как показано на рис. 6. Для этого мы растянем  DEF в барицентрической системе координат, так, чтобы он совпал с  ABC, вершины G, H и I совместим с вершинами J, K и M, а вершину O сместим на середину наибольшей стороны  ABC , заданной одним из флагов bAB , bBC, bAC (см. пример на рисунке 13).

B, E

Рис. 13. Пример коррекции барицентрических координат вершин

Описанную процедуру реализует разработанный шейдер вычисления тесселяции. Данный шейдер выполняется на всех ядрах GPU и на каждом ядре обрабатывает свою вершину (одну из вершин A,...,O). Обозначим через V одну из таких вершин, поступившую на вход шейдера вычисления тесселяции. Заменим ее барицентрические координаты на координаты одной из вершин A , B , C , J , K , M и, используя их, вычислим для вершины V текстурные координаты TV, позицию PV (в системе OCS), а также идентификатор IV (см. ниже). Это реализует разработанный Алгоритм 4.

1. Передадим в шейдер вычисления тесселяции следующие параметры:

а) барицентрические координаты Pbar вершины V (из генератора примитивов);

б) барицентрические координаты Pbar,A,..., PbarO вершин A,...,O (см. работу [19]);

в) текстурные координаты TA , TB , TC и позиции PA , PB , PC (в системе OCS) вершин  ABC;

г) флаги bAB , bBC , bAC (из шейдера управления тесселяцией).

2. Запишем флаги bA...,bO вершин A,...,O (флаг bA = 1, если |p,ar -Pbar,A| <s, где s - машинная погрешность представления действительных чисел).

3. Запишем взаимоисключающие флаги b'A, b'B, b'C, bJ, b'K , b’M (только один из них будет равен 1, остальные -0):

bA = bA V bD , bB = bB V bE , bC = bC V bF ,

b'j = bJ V bG V(bO Л bAB ) , bK = bK V bH V(bO Л bBC ) , bM = bM V bI V(bO Л bAC ) .

4. Запишем новые барицентрические координаты Pbar вершины V :

р' = V b P

1 bar uU1 bar,U

U e{A, B,C, J, K, M}

5. Вычислим текстурные координаты TV и позицию PV вершины V :

T = P' T + P' T + P' T

2V 1 bar,.В A T 1 bar,.y1 B T 1 bar,.В C

Если b'A = 1 или bB = 1, или b'C = 1, то:

Pv = b’APA + b'BPB + b'cPc.

В противном случае:

Запишем координаты (у'геогр,ф'геогр) вершины V (см. раздел 1):

w' = nTV,, o' = 2nTV .

т геогр. V,t > тгеогр. V,^

Вычислим N(w'!e<J!p) согласно (4).

Вычислим PV

согласно (3).

6. Запишем целочисленный идентификатор IV вершины V (в 0-й бит запишем значение флага b'A, в 1-ый бит -флага bB и т.д.).

Алгоритм 4. Вершинная коррекция.

Результатом применения данного алгоритма ко всем исходным вершинам будут являться вершины треугольников ПМЭ в системе OCS полигональной модели рельефа. На рисунке 14 показан результат применения алгоритма к вершинам из рис. 12 (при условии, что AB - наибольшая сторона  ABC). Из рисунка 14 видно, что некоторые треугольники ( AGD ,

GOD,  EKH и др.) будут вырожденными (мы их удалим на следующем шаге, см. ниже), а вертикальный

ABG будет обеспечивать безразрывную стыковку

полученной триангуляции с соседним, более крупным треугольником ПМЭ. Как можно заметить, в общем случае в результате тесселяции  ABC мы получим до 6 невырожденных треугольников: 4 наклонных

треугольника и от 0 до 2 вертикальных треугольников. Полученные треугольники ПМЭ передаются на следующий этап графического конвейера -геометрический шейдер.

C. Построение выходной геометрии Геометрический шейдер является завершающим этапом рассматриваемого i -го прохода графического конвейера. Данный шейдер выполняется на всех ядрах GPU и на каждом ядре обрабатывает свой треугольник ПМЭ. На данном этапе из каждого невырожденного треугольника ПМЭ будет создаваться либо треугольник ПМР, либо треугольник новой ПМЭ и накапливать такие треугольники в видеопамяти. Для этого используется вершинный буфер BPMR треугольников ПМР, а также буферы BPASS0 и BPASSl треугольников

новой ПМЭ (для четного и нечетного проходов графического конвейера). Все эти вершинные буферы создаются в видеопамяти один раз на стадии загрузки системы визуализации.

Рассмотрим процесс обработки геометрическим шейдером треугольника ПМЭ на четном i -ом проходе графического конвейера (для нечетного прохода он будет аналогичен). Обозначим такой треугольник через

V0VV2. Из предыдущих этапов графического конвейера следует, что  V0V1V2 будет либо исходным  ABC, либо одним из треугольников, образованных в результате тесселяции  ABC . В первом случае мы будем создавать в видеопамяти из V0V1V2 треугольник ПМР - таковым будет  A'B'C' (см. рис. 8), а во втором случае - сам

V0VV2. Это реализует Алгоритм 5.

1. Передадим в геометрический шейдер следующие параметры:

а) флаг btess тесселяции  ABC (из шейдера управления тесселяцией);

б) текстурные координаты TV0, TV1, TV2, позиции PV0, PV1, PV2 и идентификаторы IV0, IVt, IV2 вершин ТО (из шейдера вычисления тесселяции).

2. Если b,ess = 0, то:

Вычислим текстурные координаты TA' , TB^ TC и позиции PA,, PB,, Pc вершин  ABC':

TA' — TV0 , TB' — TV1 , TC — TV2 .

Запишем высоты вершин  ABC :

hA = FMAP (TA ) , К = FMAP (TB' ) ,

hC = Fmap (Tc ) , где Fmap - функция выборки из карты высот.

Вычислим нормали nV0, nVt, nV2 к

поверхности эллипсоида в вершинах V0V1V2 по формуле (5).

Вычислим позиции PA, (PV0,nV0,hA,),

Pb' (PV1A1A') и Pc (Py2Л2,he) по формуле (6).

Сгенерируем  ABC .

Добавим текстурные координаты и позиции

вершин  ABC' в буфер BPMR.

В противном случае:

Запишем флаг bDEG вырожденности V,J^V2:

Если IV0 — IV1 или IV0 — IV2 , или IV1 — IV2 , то

bDEG — 1 (V0V1V2 вырожден), в противном случае bDEG — 0 .

Если bDEG — 0, то создадим треугольник новой

ПМЭ:

Сгенерируем  V0VV2.

Добавим текстурные координаты и позиции вершин WV, в буфер BPASS 0.

Алгоритм 5. Создание выходного треугольника i-го прохода.

Отметим, что в буфер BPASS 0 запись координат вершин ведется с начала буфера, а в буфере BPMR перед выполнением i -го прохода мы смещаем начало записи согласно количеству треугольников ПМР, накопленных за все предыдущие проходы.

После выполнения этапа геометрического шейдера к имеющимся в буфере BPMR данным от предыдущих проходов будут добавлены данные треугольников ПМР i -го прохода, а в буфер BPASS 0 будут записаны только данные выходных треугольников ПМЭ i -го прохода. Если число таких треугольников не равно 0, то мы используем буфер BPASS,0 в качестве новой ПМЭ для

(i +1) -го прохода, а для накопления выходных треугольников используем буферы BPASS t и BPMR и т.д.

Когда число выходных треугольников ПМЭ станет равным 0, цикл проходов через графический конвейер будет завершен, а в буфере BPMR будет сформирована результирующая полигональная модель рельефа Земли.

V. Результаты

На основе предложенных методов и алгоритмов был реализован программный комплекс построения и визуализации трехмерной полигональной модели рельефа Земли. Моделирование рельефа осуществлялось на основе карты высот размером 16Kx8K, тесселяция выполнялась с использованием патчей-треугольников 6 уровней детализации.

На рисунке 15 показан пример кадра полученной визуализации гор Мьянмы с высоты 400 км с четырехкратным увеличением. Визуализация построенной модели рельефа выполнялась с наложением детализированной текстуры земной поверхности и расчетом освещенности, учитывающим физические свойства атмосферы (рассеяние, поглощение). Предложенная схема тесселяции рельефа позволила значительно сократить число визуализируемых треугольников (в среднем с 3,2 млн. до 174 тысяч) без видимых потерь качества синтезируемых изображений. Тесселяция рельефа Земли выполнялась с помощью видеокарты GeForce GTX 1080 Ti на экране с разрешением Full HD, среднее время построения полигональной модели рельефа составило менее 1 мс.

VI. Заключение

В работе рассмотрена задача построения на GPU в масштабе реального времени адаптивной полигональной модели рельефа Земли на основе детализированных глобальных карт высот, заданных в стандарте WGS-84 относительно эллипсоида вращения. Предложенное решение включает новые, распределенные методы и алгоритмы построения на GPU полигональной модели рельефа с помощью технологии шейдерной тесселяции. Распределенное построение полигональной модели рельефа выполняется на основе разбиения эллипсоида вращения на патчи-треугольники, выровненные по узлам карты высот. Обработка патчей-треугольников выполняется полностью на GPU, параллельно и независимо друг от друга, с помощью разработанных шейдера управления тесселяцией, шейдера вычисления тесселяцией и геометрического шейдера. Реализованная в шейдерах оригинальная схема тесселяции эллипсоида позволила эффективно снизить число визуализируемых треугольников и сделать менее подверженными алья-сингу синтезируемые изображения Земли. Созданные распределенные методы и алгоритмы были реализованы

в программном комплексе построения и визуализации адаптивной трехмерной полигональной модели земного рельефа. Разработанный комплекс был апробирован на детализированной карте высот Земли в составе сложной космической виртуальной сцены (с подробной моделью МКС, около 2 млн. полигонов). Полученные результаты подтвердили адекватность разработанных методов и алгоритмов поставленной задаче и их применимость для космических видеотренажеров, систем виртуального окружения, образовательных приложений и др.

Библиография

[1] В.А. Кузнецов, Ю.Г. Руссу, В.П. Куприяновский, “Об использовании виртуальной и дополненной реальности,” International Journal of Open Information Technologies, т. 7, № 4, стр. 75-84, 2019.

[2] М.А. Бондаренко, В.А. Сухомлин, “Анализ алгоритмов совмещения видеоинформации в авиационных системах,” International Journal of Open Information Technologies, т. 4, № 10, стр. 76-81, 2016.

[3] М.В. Михайлюк, М.А. Торгашев, “Система визуализации “GLView” для имитационно-тренажерных комплексов подготовки космонавтов,” Пилотируемые полеты в космос, № 4(9), стр. 60-72, 2013.

[4] В.В. Коваленок, А.С. Иванченков, С.В. Авакян, “Результативность визуально-инструментальных наблюдений в долговременных пилотируемых полетах,” Пилотируемые полеты в космос, № 4(21), стр. 103-117, 2016.

[5] П.Ю. Тимохин, М.В. Михайлюк, “Метод сжатия разрядности карт высот на основе критерия визуальной значимости,” Труды НИИСИРАН, т. 7, № 1, стр. 30-35, 2017.

[6] Л.М. Бугаевский, Математическая картография. M.: «Златоуст», 1998, 400 с.

[7] P. Cozzi, K. Ring, 3D Engine Design for Virtual Globes, Boca Raton: CRC Press, 2011.

[8] P. Cozzi, D. Bagnell, “A WebGL Globe Rendering Pipeline,” GPU Pro 4. Advanced Rendering Techniques, CRC Press, 2013, pp. 39-48.

[9] T. Ulrich, “Rendering massive terrains using chunked level of detail control,” SIGGRAPH Course Notes, vol. 3, no. 5, 2002.

[10] P. Cignoni, F. Ganovelli et al, “Planet-sized batched dynamic adaptive meshes (P-BDAM),” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 2003, pp. 147-154.

[11] L.M. Hwa, M.A. Duchaineau, K.I. Joy, “Real-time optimal adaptation for planetary geometry and texture: 4-8 tile hierarchies,” IEEE Transactions on Visualization and Computer Graphics, vol. 11, no. 4, pp. 355-368, 2005.

[12] M. Clasen, H.-C. Hege, “Terrain rendering using spherical clipmaps,”

Рис. 15. Пример кадра визуализации построенной на GPU модели рельефа Земли

Eurographics/IEEE-VGTCSymposium on Visualization, 2006.

[13] Y. Livny, Z. Kogan, J. El-Sana, “Seamless patches for GPU-based terrain rendering,” The Visual Computer, vol. 25, pp. 197-208, 2008.

[14] O. Ripolles, F. Ramos et al, “Real-time tessellation of terrain on graphics hardware,” Computers & Geosciences, vol. 41, pp. 147-155, 2012.

[15] I. Cantlay, “DirectX 11 terrain tessellation,” NVIDIA WhitePaper, 2011.

[16] E. Yusov, M. Shevtsov, “High-performance terrain rendering using hardware tessellation,” Journal of WSCG, no. 19(3), pp. 85-92, 2011.

[17] B. Mistal, “GPU Terrain Subdivision and Tessellation,” GPU Pro 4. Advanced Rendering Techniques, CRC Press, pp. 3-20, 2013.

[18] L. Zhang, J. She et al, “A Multilevel Terrain Rendering Method Based on Dynamic Stitching Strips,” ISPRS International Journal of Geo-Information, vol. 8, no. 6: 255, 2019.

[19] M.V. Mikhaylyuk, P.Y. Timokhin, A.V. Maltsev, “A Method of Earth Terrain Tessellation on the GPU for Space Simulators,” Programming and Computer Software, vol. 43, no. 4, pp. 243-249, 2017.

[20] М.Я. Выгодский, Справочник по высшей математике, M.: ACT: Астрель, 2006, p. 991.

[21] M. Bailey, S. Cunningham, Graphics Shaders: Theory and Practice. Second Edition, CRC Press, 2011, 518 p.

[22] M. Segal, K. Akeley, “The OpenGL Graphics System: A Specification. Version 4.6, Core Profile,” The Khronos Group Inc., 2006-2018. URL: https://www.khronos.org/registry/OpenGL/specs/gl/glspec46.core.pdf (review date: 19.09.2019).

[23] R.J. Rost, B. Licea-Kane, OpenGL Shading Language (3rd Edition). Boston, USA: Addison Wesley Professional, 2009, 792 p.

Real-time construction on the GPU of adaptive terrestrial relief model based on the spheroid

P.Yu. Timokhin, M.V. Mikhaylyuk, A.V. Maltsev

Abstract—The paper considers a task of real-time modelling of Earth relief based on detailed global height maps specified relative to the rotational ellipsoid (spheroid) WGS-84. The technology of adaptive tessellation of triangular patches on the GPU is proposed, which provides real-time construction of complex polygonal models of the Earth's relief. The technology includes the stage of dividing the visible ellipsoid into coarse patches (low level of detail) and the stage of their tessellation into triangles of the relief model, executed in parallel and independently on the GPU cores. The paper proposes new, distributed methods and algorithms to extract the patches needed for visualization of the current frame, to increase their level of detail in accordance with the height map and screen resolution, and to transform to relief polygons. The novelty of the work is that the tessellation of triangular patches of an ellipsoid is based on the original scheme, which allows relief areas that require high level detail to be effectively localized. This significantly reduces the time of relief model construction and improves the quality of the created images.

The developed technology, methods and algorithms were implemented in program modules and tested in the visualization system of the Earth’s virtual surface. Modelling and visualization of the Earth's relief was carried out with a detailed texture of the underlying Earth's surface and the calculation of illumination taking into account the atmosphere. The approbation was carried out in a virtual space scene comprising a detailed polygonal model of the International Space Station (ISS). The obtained results confirmed the adequacy of the proposed solution to the task and its applicability for building of space video simulators, virtual environment systems, virtual laboratories, etc.

Keywords—visualization, Earth, relief, GPU, virtual model, rotational ellipsoid, real-time, tessellation, distributed computing.

References

[1] V.A. Kuznetsov, J.G. Russu, V.P. Kupriyanovsky, “On the use of virtual and augmented reality,” International Journal of Open Information Technologies, vol. 7, no. 4, pp. 75-84, 2019.

[2] M.A. Bondarenko, V.A. Sukhomlin, “Analyze of video information alignment algorithms in aviation systems,” International Journal of Open Information Technologies, vol. 4, no. 10, pp. 76-81, 2016.

[3] M.V. Mikhaylyuk, M.A. Torgashev, ““GLView” - visualization system for simulation and training complexes for cosmonaut preparing,” Manned Spaceflight, no. 4 (9), pp. 60-72, 2013.

[4] V.V. Kovalenok, A.S. Ivanchenkov, S.V. Avakyan, “Effectiveness of Visual-Instrumental Observations in Long-Term Manned Space Flights,” Manned Spaceflight, no. 4 (21), pp. 103-117, 2016.

[5] P.Y. Timokhin, M.V. Mikhaylyuk, “The method of height map bit depth compression based on visual significance criterion,” Proceedings of SRISA RAS, vol. 7, no. 1, pp. 30-35, 2017.

[6] L.M. Bugaevskiy, Mathematical cartography. M.:“Zlatoust”, 1998, 400 p.

[7] P. Cozzi, K. Ring, 3D Engine Design for Virtual Globes, Boca Raton: CRC Press, 2011.

[8] P. Cozzi, D. Bagnell, “A WebGL Globe Rendering Pipeline,” GPU Pro 4. Advanced Rendering Techniques, CRC Press, 2013, pp. 39-48.

[9] T. Ulrich, “Rendering massive terrains using chunked level of detail control,” SIGGRAPH Course Notes, vol. 3, no. 5, 2002.

[10] P. Cignoni, F. Ganovelli et al, “Planet-sized batched dynamic adaptive meshes (P-BDAM),” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 2003, pp. 147-154.

[11] L.M. Hwa, M.A. Duchaineau, K.I. Joy, “Real-time optimal adaptation for planetary geometry and texture: 4-8 tile hierarchies,” IEEE Transactions on Visualization and Computer Graphics, vol. 11, no. 4, pp. 355-368, 2005.

[12] M. Clasen, H.-C. Hege, “Terrain rendering using spherical clipmaps,” Eurographics/IEEE-VGTCSymposium on Visualization, 2006.

[13] Y. Livny, Z. Kogan, J. El-Sana, “Seamless patches for GPU-based terrain rendering,” The Visual Computer, vol. 25, pp. 197-208, 2008.

[14] O. Ripolles, F. Ramos et al, “Real-time tessellation of terrain on graphics hardware,” Computers & Geosciences, vol. 41, pp. 147-155, 2012.

[15] I. Cantlay, “DirectX 11 terrain tessellation,” NVIDIA WhitePaper, 2011.

[16] E. Yusov, M. Shevtsov, “High-performance terrain rendering using hardware tessellation,” Journal of WSCG, no. 19(3), pp. 85-92, 2011.

[17] B. Mistal, “GPU Terrain Subdivision and Tessellation,” GPU Pro 4. Advanced Rendering Techniques, CRC Press, pp. 3-20, 2013.

[18] L. Zhang, J. She et al, “A Multilevel Terrain Rendering Method Based on Dynamic Stitching Strips,” ISPRS International Journal of GeoInformation, vol. 8, no. 6: 255, 2019.

[19] M.V. Mikhaylyuk, P.Y. Timokhin, A.V. Maltsev, “A Method of Earth Terrain Tessellation on the GPU for Space Simulators,” Programming and Computer Software, vol. 43, no. 4, pp. 243-249, 2017.

[20] M.Y. Vygodskiy, Handbook of higher mathematics, M.: ACT: Astrel’, 2006, 991 p.

[21] M. Bailey, S. Cunningham, Graphics Shaders: Theory and Practice. Second Edition, CRC Press, 2011, 518 p.

[22] M. Segal, K. Akeley, “The OpenGL Graphics System: A Specification. Version 4.6, Core Profile,” The Khronos Group Inc., 2006-2018. URL: https://www.khronos.org/registry/OpenGL/specs/gl/glspec46.core.pdf (review date: 19.09.2019).

[23] R.J. Rost, B. Licea-Kane, OpenGL Shading Language (3rdEdition). Boston, USA: Addison Wesley Professional, 2009, 792 p.