четверг, 15 апреля 2010 г.

Haskell и параллелизм: Триангуляция сферы параллельным рекуррентным разбиением тетраэдра с применением плоского затенения

Надеюсь, вас не испугал заголовок статьи. Под катом приводится описание одного из методов распараллеливания вычислений в Haskell на примере изображения затененной сферы, имеются скриншоты результатов работы программы (трафика мало).

Введение

Сначала нужно решить задачу построения сферы. Для этого обычно применяют аппроксимацию.

Чтобы аппроксимировать сферическую поверхность, можно воспользоваться методом рекурсивного разбиения, а в качестве правильного многогранника выберем тетраэдр, 4 грани которого являются равносторонними треугольниками. Можно было бы использовать любой другой правильный многоугольник; например, икосаэдр отлично подходит для аппроксимации сферы.


Вершины тетраэдра

Для определения координат вершин тетраэдра можно воспользоваться данными выкладками. Предположим, что после соответствующих вычислений мы имеем следующие координаты четырех вершин тетраэдра:
(0.0, 0.0, 1.0),
(0.0, 0.942809, -0.333333),
(-0.816497, -0.471405, -0.333333),
(0.816497, -0.471405, -0.333333)


Каркас приложения

Для визуализации используем OpenGL; Haskell имеет к ней соответствующую привязку. Теперь начнем создавать каркас приложения на Haskell с ее помощью:

main = do
(progname, _) <- getArgsAndInitialize
initialDisplayMode $= [DoubleBuffered, WithDepthBuffer]
initialWindowSize $= (Size 640 640)
createWindow "Sphere approximation"
myInit
displayCallback $= display
reshapeCallback $= Just reshape
mainLoop

myInit = do
ortho (-2.0) 2.0 (-2.0) 2.0 (-10.0) 10.0
depthFunc $= Just Less

display = do
clear [ColorBuffer, DepthBuffer]
-- put your drawing code here
swapBuffers

reshape s@(Size w h) = do
viewport $= (Position 0 0, s)
postRedisplay NothingNothing


В этом наброске в функции main производится инициализация приложения, установка глобальных переменных состояния, и создается окно. В myInit устанавливается матрица ортогональности. В initialDisplayMode мы помещаем настройки, в данном случае мы указали, что хотим использовать двойную буферизацию, а также включить буфер глубины (Z-buffer); если этого не сделать, то едва ли изображение будет выглядеть так, как мы этого ожидаем. Пока ничего отображаться в окне не должно. В функции display перед началом рисования нужно очистить буферы цвета и глубины с помощью вызова clear, а после этого, поскольку мы использовали двойную буферизацию, переключить буфер с помощью swapBuffers.


Отображение тетраэдра

Теперь, после того, как у нас имеется каркас, мы попробуем что-нибудь изобразить. Изобразим тетраэдр. Для этого добавим к коду следующие строки:

type Gf3 = (GLfloat, GLfloat, GLfloat)

figCoords :: [Gf3]
figCoords = [ (0.0, 0.0, 1.0),
(0.0, 0.942809, -0.333333),
(-0.816497, -0.471405, -0.333333),
(0.816497, -0.471405, -0.333333)
]

tetrahedron coords = do
triangle (coords !! 0) (coords !! 1) (coords !! 2)
triangle (coords !! 3) (coords !! 2) (coords !! 1)
triangle (coords !! 0) (coords !! 3) (coords !! 1)
triangle (coords !! 0) (coords !! 2) (coords !! 3)

triangle (x1, y1, z1) (x2, y2, z2) (x3, y3, z3) = do
renderPrimitive Triangles $ do
vertex $ Vertex3 x1 y1 z1
vertex $ Vertex3 x2 y2 z2
vertex $ Vertex3 x3 y3 z3


а функцию display приведем к следующему виду:

display = do
clear [ColorBuffer, DepthBuffer]
tetrahedron figCoords
swapBuffers
Поскольку никакого затенения сейчас не используется, и все изображается с применением белой заливки, то вы должны увидеть белый равносторонний треугольник. Эти три вершины формируют основание тетраэдра:



Разбиение тетраэдра

Имеющийся объект — тетраэдр — не очень-то пока похож на сферу. Поэтому попробуем сделать его чуть более похожим, произведя разбивку его граней.
Идея состоит в делении каждой грани тетраэдра на 4 равных треугольника:

/\ /\
/ \ / \
/ \ /----\
/ \ / \ / \
/________\ /___\/___\

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

divideTriangle :: Gf3 -> Gf3 -> Gf3 -> Int -> IO ()
divideTriangle p1 p2 p3 0 = triangle p1 p2 p3
divideTriangle a@(x1, y1, z1) b@(x2, y2, z2) c@(x3, y3, z3) ndivisions = do
divideTriangle a v1 v2 (ndivisions - 1)
divideTriangle c v2 v3 (ndivisions - 1)
divideTriangle b v3 v1 (ndivisions - 1)
divideTriangle v1 v3 v2 (ndivisions - 1)
where v1 = normalizeCoords ((x1 + x2)/2, (y1 + y2)/2, (z1 + z2)/2)
v2 = normalizeCoords ((x1 + x3)/2, (y1 + y3)/2, (z1 + z3)/2)
v3 = normalizeCoords ((x2 + x3)/2, (y2 + y3)/2, (z2 + z3)/2)

normalizeCoords :: Gf3 -> Gf3
normalizeCoords c@(x, y, z)
| len > 0 = (x/len, y/len, z/len)
| otherwise = c
where len = sqrt (x*x + y*y + z*z)


а также модифицируем функции tetrahedron, display и main для учета количества рекуррентных разбиений следующим образом:

tetrahedron coords ndivisions = do
divideTriangle (coords !! 0) (coords !! 1) (coords !! 2) ndivisions
divideTriangle (coords !! 3) (coords !! 2) (coords !! 1) ndivisions
divideTriangle (coords !! 0) (coords !! 3) (coords !! 1) ndivisions
divideTriangle (coords !! 0) (coords !! 2) (coords !! 3) ndivisions

display nDivisions = do
clear [ColorBuffer, DepthBuffer]
tetrahedron figCoords nDivisions
swapBuffers

main = do
(progname, _) <- getArgsAndInitialize
initialDisplayMode $= [DoubleBuffered, WithDepthBuffer]
initialWindowSize $= (Size 640 640)
createWindow "Sphere approximation"
myInit
let nDivisions = 3
displayCallback $= (display nDivisions)
reshapeCallback $= Just reshape
mainLoop

После компиляции и запуска такой программы изображается что-то наподобие этого:

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


Затенение

Существует множество способов затенения. Одним из самых первых и распространенных был метод плоского затенения, где освещенность производится для полигона (треугольника) целиком. Этот метод хорошо подходит для затенения каких-либо объектов с прямыми гранями, но не совсем хорошо — для неровных поверхностей, поскольку при таком затенении на неровном объекте отчетливо становятся видны грани (Mach bands, статья "Как компьютер рисует изображения"). Тем не менее, нам не важна в данном случае фотореалистичность. Поэтому применим плоское затенение.

Для определения степени освещенности точки — в нашем случае освещается треугольник — нужно определить косинус угла между нормалью, идущей из этой точки, и вектором идущим из этой точки к точечному источнику света (статья в журнале Xakep "Освещение в играх")

Чтобы найти нормаль, можно воспользоваться векторным произведением:

crossProd :: (Num a) => (a, a, a) -> (a, a, a) -> (a, a, a) -> (a, a, a)
crossProd (x1, y1, z1) (x2, y2, z2) (x3, y3, z3) = (x, y, z)
where x = (y2 - y1) * (z3 - z1) - (z2 - z1) * (y3 - y1)
y = (z2 - z1) * (x3 - x1) - (x2 - x1) * (z3 - z1)
z = (x2 - x1) * (y3 - y1) - (y2 - y1) * (x3 - x1)

Нахождение косинуса угла между полученной нормалью и вектором к источнику света может замедлить вычисления, поэтому можно воспользоваться свойством скалярного произведения векторов:

dotProd :: (Num a) => (a, a, a) -> (a, a, a) -> a
dotProd (x1, y1, z1) (x2, y2, z2) = x1*x2 + y1*y2 + z1*z2

Функцию triangle изменим следующим образом:

triangle :: Gf3 -> Gf3 -> Gf3 -> IO ()
triangle a@(x1, y1, z1) b@(x2, y2, z2) c@(x3, y3, z3) = do
-- Detrmine a triangle normal
-- (it's just a direction, not a vector)
-- Normalized normals are happy normals ;)
let normal = normalizeCoords $ crossProd a b c
let lightPos = (1.0, 1.0, 1.0::GLfloat)
let intensity = 1.0 * dotProd normal lightPos

renderPrimitive Triangles $ do
color $ Color3 intensity intensity intensity
vertex $ Vertex3 x1 y1 z1
vertex $ Vertex3 x2 y2 z2
vertex $ Vertex3 x3 y3 z3

Здесь координаты источника света (1.0, 1.0, 1.0), коэффициент отражающей способности установлен в 1.0, интенсивность отражения изменяется для всех цветовых каналов одинаково. Теперь, если вы запустите программу, то можно будет увидеть следующее изображение (48 треугольников):

По-моему, оно походит чем-то на подсвеченную сферу. А если увеличить число разбиений до 8-1, то получится очень реалистичная сфера (49 152 треугольника):

Полный код приложения на данном этапе.

Параллелизм

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

Для Haskell имеются разнообразные библиотеки, позволяющие сделать приложение многопоточным. Мы воспользуемся вызовом forkIO, и синхронизирующей переменной MVar, которые достаточно давно доступны в Glasgow Haskell Compiler. Импортируем их, добавив следующие строки в начало исходного кода:

import Control.Concurrent
import Control.Concurrent.MVar
import Control.Monad

К текущему моменту в нашей программе можно распараллелить вызовы функции divideTriangle, осуществляемые в функции tetrahedron:

tetrahedron :: [Gf3] -> Int -> IO ()
tetrahedron coords ndivisions = do
child0 <- thread (coords !! 0) (coords !! 1) (coords !! 2) ndivisions
child1 <- thread (coords !! 3) (coords !! 2) (coords !! 1) ndivisions
child2 <- thread (coords !! 0) (coords !! 3) (coords !! 1) ndivisions
child3 <- thread (coords !! 0) (coords !! 2) (coords !! 3) ndivisions
ret <- mapM takeMVar [child0, child1, child2, child3]
return ()

thread :: Gf3 -> Gf3 -> Gf3 -> Int -> IO (MVar (IO ()))
thread p1 p2 p3 ndivisions = do
mvar <- newEmptyMVar
forkIO $ do
let r = divideTriangle p1 p2 p3 ndivisions
putMVar mvar r
return mvar

Здесь создаются 4 нити, работающие параллельно. В функции thread перед созданием нити инициализируется синхронизирующая переменная с помощью вызова newEmptyMVar. Потом создается сама нить, которая осуществляет вызов divideTriangle и помещает в MVar результат, который был возвращен функцией divideTriangle.

Само по себе содержимое MVar нас не интересуит, оно лишь служит сигналом нити, выполняющей вызов функции takeMVar, к тому, что можно продолжать дальнейшие вычисления, т. е. функция takeMVar не вернет управление, пока в MVar не будет помещено какое-либо значение. Таким образом, в списке ret будут находиться значения, которые нити поместили с помощью вызова putMVar, и программа не продолжит работу, пока все нити не поместят эти значения. Вызов return () позволяет вернуть значение, которое соответствует типу функции tetrahedron :: [Gf3] -> Int -> IO (), т. е. пустой кортеж.

Но если вы скомпилируете и запустите эту программу, то прежней сферы вы не увидите! Предлагаю вам подумать, почему так происходит (вспомните свойства языка Haskell). Полный код программы на данном этапе.

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

Заключение

В статье была рассмотрена техника распараллеливания вычислений на языке Haskell, продемонстрировано использование привязки к OpenGL. Как говорит эта статья, Haskell производит паралельные вычисления очень быстро.

Надеюсь, вам понравилось.

Литература

Эйнджел, Эдвард. Интерактивная компьютерная графика. Вводный курс на базе OpenGL, 2 изд.: Пер. с англ. — М.: Издательский дом "Вильямс", 2001. — 592 с.: ил. — Парал. тит. англ. ISBN 5-8459-0209-6 (рус.)

P. S. Возможно, кто-то посоветует способ сокрытия под кат текста так, чтобы это действительно работало в т. ч. и на Russian Lambda Planet.

2 комментария:

  1. Для сокрытия статьи под кат использую тег &lt!--more--&gt. Хотя в Russian Lambda Planet RSS отрезается даже раньше, чем тег.

    ОтветитьУдалить