Курсовая работа: Реализация алгоритма обратной трассировки лучей для моделей с большим числом полигонов

Факультет информатики и систем управления

Кафедра "Программное обеспечение ЭВМ и информационные технологии"

Курсовой проект

по машинной графике

Расчетно-пояснительная записка

Тема:

"Реализация алгоритма обратной трассировки лучей для моделей с большим числом полигонов"

Москва 2004

Оглавление

1. Введение

2. Конструкторская часть

2.1 Обоснование использованных алгоритмов

2.2 Структура данных

2.2.1 Источники света

2.2.2 Объекты для визуализации

2.2.3 Текстуры

2.3 Алгоритм обратной трассировки лучей

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

2.3.2 Математическая основа обратной трассировки лучей

2.3.3 Составление матрицы

2.3.4 Программная реализация

2.3.5 Определение пересечения луча с треугольником

2.3.4 Формирование отраженного луча

2.3.5 Формирование преломленного луча

2.4 Оболочки

2.4.1 Алгоритм построения иерархических оболочек

2.4.2 Алгоритм обхода оболочек в трассировке лучей

2.5 Текстурирование

2.5.1. Процедуры для работы с текстурами

2.5.2. Собственно текстурирование

2.6 Закраска Фонга

2.7 Освещение

2.7.1. Модель освещения Уиттеда

2.7.2 Диффузное отражение

2.7.3 Зеркальное отражение

2.7.4 Фоновая освещенность

2.7.5 Прозрачность

2.7.6 Процедуры расчета освещенности

3. Технологическая часть

3.1 Выбор языка программирования и обоснование выбора

3.2 Модульная структура программы

3.3 Интерфейс программы

4. Экспериментально-исследовательская часть

Тест № 1

Тест № 2

Тест № 3

Заключение

Список литературы

1. Введение

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

В моей курсовой работе для рендеринга применяется алгоритм обратной трассировки. На сегодняшний день он считается одним из лучших для формирования реалистических изображений. Его используют большинство трехмерных графических редакторов. Применяется так же алгоритм z-буфера. Но эффекты отражения и преломления, там достигаются искусственными методами и считаются не совсем точно. Подобным методом пользуются в программах, где крайне важна скорость. Примером таких программ служат компьютерные игры. Почти все они построены на использовании алгоритма z-буфера.

Целью моей курсовой было так же сделать алгоритм обратной трассировки как можно более быстрым. Для этого я применил метод иерархических оболочек. Его применение позволяет сделать время рендеринга, пропорциональным логарифму от числа объектов, а не числу объектов. Добиться с помощью этого реального времени, конечно, не удастся, но делает время ожидания приемлемым, равным порядка 5-30 секунд для 30000 треугольников на сцене.

Модуль Engine программы, может быть использован отдельно, в других программах Delphi. С помощью всего нескольких функций пользователь сможет задать сцену любой сложности и произвести рендеринг сцены. Модуль содержит функции для:

управления камерой

управления источниками света

задания объектов на сцене.

поворота объектов

рендеринга сцены

вывода изображения в задаваемое окно

По использованию модуль Engine очень похож на модуль OpenGL.

2. Конструкторская часть

2.1 Обоснование использованных алгоритмов

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

Достоинством алгоритма является то, что он не требователен к памяти, в отличие от алгоритма z-буфера. А недостатком является то, что работает он сравнительно долго и не позволяет строить изображения в реальном времени.

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

Для сглаживания изображения применен алгоритм закраски Фонга. Он является самым затратным по времени. Метод Гуро, например, быстрее Фонга примерно в 5 раз. Но время его выполнения от общего времени рендеринга не превышает 3 процентов. Зато он дает великолепные результаты. В частности, блики выглядят куда реалистичнее, чем если использовать метод Гуро.

2.2 Структура данных

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

2.2.1 Источники света

Источники света не имеют никаких геометрических размеров, они являются точечными и не рисуются при рендеринге. Информация об источниках света хранится в массиве Svet. В i-ом элементе массива хранится информация об i-ом источнике света. Элемент массива представляет собой запись:

TLight=record

tip: integer;

lim: real;

Center: TPoint;

R,G,B: real;

DirX,DirY,DirZ: real;

end;

Поле tip содержит информацию о типе источника. Если оно равно 1, то источник светит во все стороны. Если оно равно 2, то источник светит внутри конуса, направляющая которого DirX, DirY, DirZ, а угол при вершине равен 2*Lim. Угол измеряется в радианах. Если тип источника - 3, то источник также светит в конусе, но по мере отклонения от образующей его интенсивность уменьшается и на угле Lim равна нулю.

Поле Center содержит координаты источника в глобальной системе координат.

Поля R,G,B содержат интенсивность источника по красной, зеленой и синей компоненте. Они могут принимать значения от 0 до 1.

Если источник первого типа, то нет необходимости вводить поля DirX, DirY, DirZ и Lim, так как они не требуются для расчета интенсивности.

2.2.2 Объекты для визуализации

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

Objects (массив объектов), Vse (массив треугольников), Toch (массив точек).

Массив Objects

Элемент массива представляет собой запись:

TObj=record

StartT,EndT: integer;

StartG,EndG: integer;

XC,YC,ZC,R: real;

nnn,NPr: real;

end;

StartT, EndT соответствуют индексам в массиве точек. Они указывают, что точки с номером, большим или равным StartT и меньшим или равным EndT, принадлежат данному объекту.

StartG, EndG соответствуют индексам в массиве треугольников. Они указывают, что треугольники с номером, большим или равным StartG и меньшим или равным EndG, принадлежат данному объекту.

В NPr содержится показатель преломления данного объекта.

В nnn содержится коэффициент затухания света в данном объекте.

Массив Toch

Элемент массива представляет собой запись:

TApex=record

X,Y,Z: real;

nx,ny,nz: real;

end;

Поля X,Y,Z содержат координаты точки.

Поля nx, ny, nz содержат значение нормали в данной точке. Эти поля используются при закраске по методу Фонга.

Массив Vse

Массив содержит полную информацию обо всех треугольниках сцены.

Элемент массива представляет собой запись:

TGran=record

Nom: array [1. .3] of integer;

ColorR,ColorG,ColorB: Byte;

KOt,KPr,KRas,KDif,KBlik: real;

Tek: array [1. .3] of array [1. .2] of integer;

TNom: integer;

PaintType: boolean;

XC,YC,ZC,R: real;

O: integer;

p: real;

end;

Массив Nom содержит номера точек, которые являются вершинами треугольника.

ColorR, ColorG, ColorB содержат цвет треугольника.

Поля KOt, KPr, KRas, KDif, KBlik, содержат оптические коэффициенты поверхности треугольника.

O - номер объекта, которому принадлежит данный треугольник.

XC, YC, ZC, R - координаты центра и радиус сферической оболочки треугольника.

PaintType - способ закраски треугольника.

TNom - номер текстуры, которая наложена на треугольник.

Массив Tek содержит текстурные координаты, каждой вершины треугольника.

Запись треугольника не содержит координат вершин, она содержит ссылки на вершины. Таким образом, сразу несколько треугольников, могут ссылаться на одну и ту же вершину.

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

Это упрощает изменение координат вершины. Если бы координаты вершин хранились непосредственно в записи треугольника, то пришлось бы изменить все записи треугольников, включающие данную вершину. В данном случае надо изменить координаты только одной точки.

Очень удобно поворачивать, объекты. Если необходимо повернуть заданный объект, достаточно пробежать по всем его точкам и изменить их координаты.

Настройки цвета, коэффициентов и сглаживания у каждого треугольника свои, а не одинаковые у всех треугольников объекта. Это дает большую свободу в формировании сцены.

Коэффициенты затухания и преломления задаются в записи объекта, так как они характеризуют весь объект целиком.

Данная организация не обладает избыточностью, все треугольники и вершины описываются ровно один раз.

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

Иерархические оболочки

Для хранения иерархических оболочек используется массив Shapes. Он состоит из записей:

TShape=record

tip: integer;

S: integer;

G: TSpisok;

S1: integer;

G1: TSpisok;

Low: array [1. .8] of integer;

NLow: integer;

XC,YC,ZC,R: real;

end;

Первый элемент в массиве Shapes соответствует оболочке, включающей все треугольники сцены.

Поле tip принимает два значения: 1 и 2. Если у оболочки нет подоболочек, то tip равен 2, в противном случае равен 1.

G - это указатель на список треугольников, принадлежащих данной оболочке, S - их число.

G1 - это указатель на список треугольников, которые принадлежат оболочке и очень большие, S1 - их число.

Low - массив содержащий номера подоболочек в массиве Shapes, Nlow - число этих подоболочек.

XC, YC, ZC - координаты центра этой оболочки.

R - радиус оболочки.

Таким образом, в памяти разворачивается дерево. Из корня (т.е. из первой оболочки) его легко обойти. Проще всего это можно сделать, пользуясь рекурсивным алгоритмом.

2.2.3 Текстуры

Информация о текстурах хранится в массиве Tex. Для каждой текстуры хранятся ее размеры (lx, ly) и указатель на область памяти, куда загружена текстура (PT).

TTex=record

lx,ly: integer;

PT: PRGBI;

end

2.3 Алгоритм обратной трассировки лучей

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

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

Рассмотрим, как формируется изображение. Изображение получается из-за того, что свет попадает в камеру. Выпустим из источников света множество лучей. Назовем их первичными лучами. Часть этих лучей улетит в свободное пространство, а часть попадет на объекты. На них лучи могут преломиться, отразится. При этом часть энергии луча поглотится. Преломленные и отраженные лучи образуют множество вторичных лучей. Далее эти лучи опять же преломятся и отразятся и образуют новое поколение лучей. В конечном итоге часть лучей попадет в камеру и сформирует изображение.

Существуют алгоритмы, работающие по такому алгоритму. Но они крайне неэффективны, так как большинство лучей, исходящих из источника, не попадают в камеру. А приемлемая картинка получается, если трассировать большое число лучей, что займет очень много времени. Данный алгоритм называется прямой трассировкой лучей.

Метод обратной трассировки лучей позволяет значительно сократить перебор световых лучей. Этот метод разработали в 80-х годах Уиттед и Кэй. В этом методе отслеживаются лучи не от источников, а из камеры. Таким образом, трассируется определенное число лучей, равное разрешению картинки.

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

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

Если объект прозрачный, то необходимо построить вторичный луч такой, чтобы при преломлении он давал исходный луч. Некоторые тела могут, обладать свойством диффузного преломления. При этом образуется не один, а множество преломленных лучей. Как и в случае отражения, я этим пренебрегаю.

Таким образом, первичный луч, найдя пересечение с объектом, делится в общем случае на два луча (отраженный и преломленный). Далее эти два луча делятся еще на два и так далее.

Рис 1.

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

Если поколение луча равно максимальной глубине рекурсии, то возвращаем среднюю яркость по всем составляющим. Если нет, то идем дальше

Определяем ближайший треугольник, с которым пересекается луч.

Если такого треугольника нет, возвращаем цвет фона, если есть, идем дальше.

Если поверхность, с которой было найдено пересечение, отражает, то формируем отраженный луч и вызываем рекурсивно процедуру Ray с поколением луча, увеличенным на 1.

Если поверхность, с которой было найдено пересечение, преломляет, то формируем преломленный луч и вызываем рекурсивно процедуру Ray с поколением луча, увеличенным на 1.

Определяем итоговую освещенность пиксела, учитывая расположение источников, свойства материала, а так же интенсивности отраженного и преломленного луча.

Я уже рассмотрели ряд ограничений метода трассировки, когда говорили о диффузном преломлении и о неровном зеркале. Рассмотрим и некоторые другие.

Освещать сцену могут только специальные объекты - источники света. Они точечные и не могут поглощать, преломлять и отражать свет.

Свойства отражающей поверхности состоят из двух компонент - диффузной и зеркальной.

При диффузном отражении учитываются только лучи от источников света. Если источник освещает точку, через зеркало (зайчиком), то считается, что точка не освещена.

Зеркальность тоже делится на две составляющие.

reflection - учитывает отражение от других объектов (не источников света)

specular - учитывает блики от источников света

В трассировке не учитываются зависимости от длины волны света:

коэффициента преломления

коэффициента поглощения

коэффициента отражения

Так как я не моделирую диффузное отражение и преломление, то не смогу получить фоновую подсветку. Поэтому вводим минимальную фоновую освещенность. Часто она позволяет просто значительно улучшить качество изображения.

Алгоритм трассировки позволяет рисовать очень качественные тени. Это не потребует большой переделки алгоритма. В него придется кое-что добавить. При расчете освещенности точки необходимо пустить в каждый из источников света "Теневой фронт". "Теневой фронт" - это луч, с помощью которого проверяется, лежит ли что-нибудь между точкой и источником. Если между ними лежит непрозрачный объект, то точка находится в тени. Это значит, что данный источник, не делает свой вклад в итоговую освещенность точки. Если лежит прозрачный объект, то интенсивность источника уменьшается. Прорисовка теней является очень затратной по времени. Так что, в некоторых ситуациях их отключают.

В моей программе есть возможность включить сглаживание изображения. Сглаживание заключается в том, что для определения цвета пиксела. пускается не один луч, а четыре и определяется среднее значение цвета у этих лучей. Если необходимо найти цвет пиксела (i,j), то пускаются 4 луча в точки экранной плоскости с координатами (i-0.25,j-0.25), (i-0.25,j+0.25), (i+0.25,j-0.25), (i+0.25,j+0.25).

2.3.2 Математическая основа обратной трассировки лучей

Координаты всех объектов сцены определены в некой глобальной системе координат (в том числе и камеры). После формирования первичного луча создадим подсистему, у которой центр совпадает с точкой выхода луча и ось OZ направлена по лучу. Вычислим матрицу перехода из первой системы координат во вторую. Это позволит просто искать пересечения с треугольником, со сферой, векторы преломления и отражения. При необходимости переводим координаты нужных объектов в новую систему координат и работаем уже в ней. Если необходимо построить вторичный луч, создаем еще одну систему координат, связанную с вторичным лучом, и считаем матрицу для перехода из 2 системы в 3. Чтобы получить матрицу перехода из 1 в 3 необходимо умножить матрицу перехода из 2 в 3 умножить на матрицу перехода из 1 в 2. Таким образом, мы работаем все время в какой-то подсистеме. Нам не надо переводить никакие координаты обратно в глобальную систему. Поэтому не надо и составлять обратную матрицу.

2.3.3 Составление матрицы

Составление матрицы преобразования из текущей системы координат в систему координат, центр которой находится в точке (x, y, z) и ось OZ которой направлена по (dx, dy, dz). Для такого преобразования необходимо:

совершить сдвиг в точку (x, y, z)

совершить поворот вокруг OZ

совершить поворот вокруг OX

1. Матрица сдвига: .

2. Необходимо совершить поворот вокруг оси OZ по часовой стрелке на угол α.

.

Матрица поворота, таким образом, будет равна:

3. Необходимо совершить поворот вокруг оси OX по часовой стрелке на угол β.

.

Матрица поворота, таким образом, будет равна:

Умножив M₃ на M₂ , а результат на M₁ получим искомую матрицу перехода:

2.3.4 Программная реализация

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

TMatrix=Array [0. .11] of real

Это массив из 12 элементов типа real. Он представляет собой последовательно записанные три верхние строчки матрицы. Я не включил последнюю строчку, так как она одинаковая у всех матриц преобразования и равна (0, 0, 0,1).

Для операций над матрицами в программе предусмотрены следующие процедуры:

1. Procedure MatrixAB (var Res: TMatrix; const M1,M2: TMatrix)

Процедура умножает матрицу M1 на матрицу M2 и помещает результат в Res.

2. Procedure ShiftMatrix (var M: TMatrix; Z: real)

Создает матрицу перехода из текущей системы координат в систему координат, сдвинутую по оси OZ на z.

Процедура перемещает систему координат, задаваемую матрицей M, по оси OZ на z.

3. Procedure SetMatrix (var M: TMatrix; dx,dy,dz,x,y,z: real) overload

Создает матрицу перехода из текущей системы координат в систему координат, находящуюся в точке (x,y,z), ось OZ которой направлена по вектору (dx,dy,dz).

4. Procedure SetMatrix (var M: TMatrix; dx,dy,dz: real) overload

Создает матрицу перехода из текущей системы координат в систему координат, ось OZ которой направлена по вектору (dx,dy,dz).

Преобразование координат

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

Для преобразования точки из одной системы координат в другую в программе существует процедура Trans (const M: TMatrix; var F: TPoint; const V: TPoint).

В V содержатся координаты точки, координаты которой надо преобразовать.

В F содержатся результат.

M - матрица преобразования.

В процедуру Ray передается только матрица перехода из глобальной системы координат в систему, связанную с лучом (M_i ).

Процедура находит координаты вторичного луча в новой системе координат.

Составляет матрицу перехода из текущей системы в систему, связанную с лучом (L_{i +1} ).

Умножает матрицы M_{i +1} =L_{i +1} M_i

Вызывает рекурсивно Ray с параметром M_{i +1}

2.3. 5 Определение пересечения луча с треугольником

Преобразуем все вершины треугольника в локальную систему координат, связанную с лучом. Луч в этой системе координат имеет координаты (0, 0,1). После этого задача сводится плоской. Необходимо определить, лежит ли точка (0, 0) внутри треугольника в проекции на плоскость OXY.

Преобразуем координату x вершин треугольника в локальную систему координат.

Проверяем, лежат ли точки треугольника все справа от нуля или все слева. Если да, то пересечения с треугольником нет. Если нет, то:

Преобразуем координату у вершин треугольника в локальную систему координат.

Проверяем, лежат ли точки треугольника все сверху от нуля или все снизу. Если да, то пересечения с треугольником нет. Если нет, то:

Необходимо, чтобы при обходе треугольника по часовой стрелке точка (0,0) лежала справа от каждой стороны (либо наоборот). Это можно установить, проверив одного ли знака три векторных произведения:

, , .

Если не одного знака, то пересечения нет, Если одного, то:

Преобразуем координату z вершин треугольника в локальную систему координат.

Определяем нормаль к треугольнику, для этого умножим векторно два вектора

и

В случае если, NZ равно нулю, луч параллелен оси OZ и соответственно лучу. Значит, пересечения нет. В другом случае пересечение есть.

Находим координаты пересечения.

Составим уравнение плоскости, в которой находится треугольник:

,

подставим x=0 и y=0,

2.3.4 Формирование отраженного луча

Обозначим отраженный луч через R, вектор, направленный против падающего луча - S, вектор нормали - N. Рассмотрим единичные векторы этих векторов R₁ , S₁ , N₁ . Так как все три вектора находятся в одной плоскости, то можно записать R₁ +S₁ =N^’ . Длина вектора N^’ равна 2cosθ. Так как векторы N^’ и N₁ сонаправленные, то можно записать:

N’=N₁ 2cosθ.

Таким образом

.

Поставим условие, что падающий и отраженный лучи имеют одинаковую длину.

Так как падающий луч в локальной системе координат имеет координаты (0, 0,1). То вектор S будет иметь координаты (0, 0, - 1). Подставим его координаты в выражение для отраженного луча. Получим

2.3.5 Формирование преломленного луча

Обозначим преломленный луч, имеющий единичную длину, через T₁ . Единичный вектор нормали - через N₁ . А вектор, направленный противоположно падающему лучу - через S₁ . Разложим вектор S₁ на A и N_s , а вектор T₁ на B и N_T . Вектор равен

.

Найдем вектор N_T . Этот вектор противоположен по направлению вектору нормали, а длина его равна

.

Таким образом

.

Для того, что бы определить cosα₂ , запишем закон преломления

.

Воспользуемся тождеством

Получим

Значение для cosα₁ можно выразить через скалярное произведение единичных векторов S₁ и N₁ . С учетом этого можно записать выражение для вектора

N_T :

Найдем вектор B. Он располагается на одной прямой с вектором A, причем

.

Найдем отношение длин векторов A и B:

. Отсюда

.

Если подкоренное выражение отрицательно, то это соответствует полному внутреннему отражению. В этом случае преломленный луч создаваться не будет. Учтем то, что вектор S равен (0, 0,1).

2.4 Оболочки

В алгоритме трассировки лучей от 70 до 90 процентов временных затрат занимает процедура определения пересечений луча с объектами сцены. Если перебирать все объекты сцены, то время будет пропорционально Cn, где С - количество пикселей, а n - число объектов на сцене. Улучшить алгоритм можно, если каким-нибудь образом попробовать сократить число перебираемых объектов. Очень простой, но эффективной является идея иерархических оболочек. Предположим, что нам надо изобразить содержимое полки с посудой. Проведем мысленно большую сферу вокруг полки, так чтобы она включала и полку, и посуду на ней. Затем сферу вокруг каждого предмета посуды. Теперь представим себе процесс рендеринга. Проверяем, пересекается ли луч с самой большой сферой. Если нет, то это значит, что луч не пересекает внутренние оболочки, переходим к другому лучу. Если пересекает, то смотрим пересечение луча с подсферами данной сферы. Как видно из такого примера, многие лучи могут совершить всего одну проверку и отсеяться.

При построении оболочек необходимо, чтобы главная оболочка целиком включала все ее подоболочки, иначе не будет работать правило: "Если луч не пересекает главную оболочку, то он не пересекает и все ее подоболочки". Так же при построении оболочек желательно, чтобы оболочки, имеющие один порядок вложенности, пересекались по как можно меньшему объему. Это улучшит эффективность алгоритма. Метод оболочек помогает сделать время рендеринга пропорциональным Clog (n).

2.4.1 Алгоритм построения иерархических оболочек

Для начала строятся сферические оболочки вокруг всех треугольников сцены. Этим занимается процедура Obol1. Сферы выбираются таким образом, чтобы все точки треугольника лежали внутри сферы. Для построения оболочки вокруг треугольника процедура:

Находит минимальный параллелепипед, в котором находится весь треугольник.

Определяет середину параллелепипеда.

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

Далее вызывается процедура Obol2. Она и занимается собственно построением иерархии. Процедура создает одну оболочку, в которую записывает номера всех треугольников, сцены и вызывает рекурсивную процедуру step с параметром 1. Это значит, что step обработает первую оболочку.

Процедура step:

Определяет центр и радиус сферы, которая включает в себя все треугольники данной оболочки.

Перемещаем мысленно систему координат в центр этой сферы. Создаем 8 оболочек, каждая оболочка отвечает за свой октант в перенесенной системе координат.

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

Если какая-то оболочка оказалась пуста (т.е. если в нее не попал не один треугольник), то она уничтожается.

Устанавливается связь: в запись оболочки добавляются номера ее подоболочек.

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

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

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

2.4.2 Алгоритм обхода оболочек в трассировке лучей

Устанавливаем текущей первую оболочку.

Для трассируемого луча проверяется, пересекается ли он с текущей оболочкой. Если нет, то выводим цвет фона. Если да, то идем на п.3.

Пересекаем луч с большими треугольниками данной оболочки.

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

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

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

2.5 Текстурирование

2.5.1. Процедуры для работы с текстурами

Для работы с текстурой предусмотрены 3 функции:

1. AddTeksture

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

имела размеры, равные степени двойки

была в формате BMP

имела глубину цвета равную 24.

Функция возвратит true при успешной загрузке текстуры и false при неудаче.

2. CloseTex

Удаляет все ранее загруженные текстуры из памяти.

3. GetTexPoint

Эта функция позволяет получить 3 компоненты цвета точки текстуры. Входными параметрами для нее являются:

X и Y - координаты точки текстуры

Nom - номер текстуры

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

Третья функция используется непосредственно на этапе текстурирования.

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

2.5.2. Собственно текстурирование

Поставим в соответствие каждому треугольнику формулы преобразования координат точек треугольника в текстурные координаты:

Определим коэффициенты a,b,c,d,e,f.

Поставим в соответствие каждой вершине треугольника нужную текстурную координату.

Мы получили две системы линейных уравнений. Системы будут иметь единственное решение в том случае, если определитель системы не равен нулю. Определим его значение.

Это условие соответствует тому, что три точки треугольника лежат на одной прямой. Если это так, то треугольник не текстурируем.

В противном случае определяем коэффициенты. Точка пересечения треугольника и луча имеет во вспомогательной системе координат нулевые координаты X и Y.

Поэтому X_T =C и Y_T =F. Имеет смысл искать только коэффициенты С и F.

Далее вызывается, описанная выше процедура GetTexPoint, с текстурными координатами round (C) и round (F). Получаем цвет нужной нам точки треугольника.

2.6 Закраска Фонга

В программе есть возможность сгладить выборочно нужные треугольники. Для этого в атрибутах треугольника есть флаг PaintType. Если он равен True, совершается сглаживание Фонга. Если равен False, то треугольник не сглаживается. Для удобства в программе введены две константы Const_Paint_Fong и Const_Paint_Flat, равные соответственно True и False. Наличие такого флажка, делает возможным строить практически любые по форме тела.

Закраска Фонга заключается в следующем:

Определяются нормали в вершинах грани

Определяются внешние нормали у всех граней, содержащих данную вершину

Нормаль в вершине равна среднему значению нормалей прилежащих граней

Билинейной интерполяцией вычисляется нормаль в каждом пикселе.

В данной программе первый шаг алгоритма не осуществляется. Необходимо еще при моделировании сцены определить значение нормалей в вершинах. Второй шаг алгоритма использует метод использованный в текстурировании. Мы ставим в соответствие каждому треугольнику формулы преобразования координат точек треугольника в x,y,z компоненты нормали:

Опуская вычисления, определим с, f, e.

Вычисленные значения являются значениями нормали в точке пересечения луча и треугольника.

Закраска Фонга является очень трудоемкой, она намного медленнее закраски Гуро, так как приходится интерполировать три величины и приходится каждый раз совершать пересчет освещенности. Но она является наиболее качественной, так как дает реалистичные блики и качественные плавные изменения яркости на гладких телах. К недостаткам метода можно отнести то, что при создании последовательных кадров закраска может заметно меняться от кадра к кадру.

2.7 Освещение

2.7.1. Модель освещения Уиттеда

В данной программе я применил модель освещения Уиттеда. Она задается следующей формулой:

K_a - коэффициент рассеянного отражения.

K_d - коэффициент диффузного отражения.

K_s - коэффициент зеркальности.

K_r - коэффициент отражения.

K_t - коэффициент преломления.

I_a - интенсивность фонового освещения.

I_d - интенсивность, учитываемая для диффузного рассеивания.

I_s - интенсивность, учитываемая для зеркальности.

I_r - интенсивность излучения, приходящего по отраженному лучу.

I_t - интенсивность излучения, приходящего по преломленному лучу.

С - цвет поверхности.

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

При расчете освещенности, необходимо посчитать не одну формулу, а три (по каждой компоненте света: красной, синей и зеленой). При этом вместо С используется нужная компонента цвета поверхности, вместо I_r и I_t нужные компоненты излучения, приходящего по отраженному и преломленному лучу. Эти значения будут возвращены процедурами Ray, вызванными для преломленного и отраженного лучей.

2.7.2 Диффузное отражение

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

, где θ - угол между нормалью и направлением на источник.

Косинус угла между двумя векторами можно посчитать по формуле:

.

Соответственно cos (θ) можно посчитать по формуле

.

Так как нормаль - это единичный вектор, то

Поэтому, для расчета диффузной интенсивности в формуле Уиттеда необходимо просуммировать все Icosθ для каждого, видимого из данной точки, источника.

2.7.3 Зеркальное отражение

Падающий луч, попадая на слегка шероховатую поверхность реального зеркала порождает, не один отраженный луч, а несколько лучей, рассеиваемых по различным направлениям. Зона рассеивания зависит от качества полировки и может быть описана некоторым законом распределения. Как правило, форма зоны рассеивания симметрична относительно линии идеального зеркально отраженного луча. К числу простейших, но достаточно часто используемых, относится эмпирическая модель распределения Фонга, согласно которой интенсивность зеркально отраженного излучения пропорциональна cos^p α. α - угол отклонения от линии идеально отраженного луча. Показатель p находится в диапазоне от 1 до 200 и зависит от качества полировки.

.

α - угол между вектором наблюдения и вектором отражения луча из данного источника (R и S).

Но можно поступить наоборот. Можно найти угол между отражением S и L, он будет равен первому. Но в данной процедуре выгодно искать угол между S и L. Это дает экономию в расчетах. Так как во вспомогательной системе координат падающий луч (S) и ось OZ совпадают. Поэтому отраженный луч (W) ищется очень просто:

W_x =-N_x *N_y

W_y =-N_y *N_z

W_z =0,5-N_z ²

Поэтому

Для расчета зеркальности в формуле Уиттеда необходимо просуммировать все Icos^p α для каждого, видимого из данной точки, источника.

2.7.4 Фоновая освещенность

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

Полностью формула Уиттеда будет выглядеть:

Интенсивности излучений, приходящих по отраженному и преломленному лучу умножаются на , где R - расстояние, пройденное лучом в среде, а δ - коэффициент затухания света в среде. Таким образом, я предусматриваю затухание света в среде.

2.7.5 Прозрачность

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

Простая модель.

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

Модель преломления.

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

1-K_d -K_s -K_a (Part).

Если коэффициент преломления объекта, в который входит луч больше 10000, то формируется только отраженный луч, поскольку преломленный луч в такой оптически плотной среде будет иметь очень малую интенсивность. Коэффициент отражения в этом случае равен Part, вся энергия пойдет на отражение.

Если коэффициент преломления объекта меньше 10000, то определяем, соответствует ли угол падения луча углу полного внутреннего отражения.

Если это угол полного внутреннего отражения, то преломленный луч не формируется, формируется только отраженный луч, с коэффициентом отражения, равным Part.

Если это не угол полного внутреннего отражения, то формируются и преломленный, и отраженный луч. Зависимость между энергией преломленного и отраженного луча выражает формула Френеля:

Она определяет долю отраженной энергии. Доля пропущенной энергии равна 1-F. Соответственно, коэффициент отражения устанавливается равным Part*F, а преломления Part* (1-F).

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

2.7.6 Процедуры расчета освещенности

Расчет диффузной и зеркальной составляющих осуществляет процедура Elluminaty.

Elluminaty (const M: TMatrix; MaxN: integer; NX, NY, NZ, RX, RY, RZ: real; BodyR, BodyG, BodyB: byte; var DifR, DifG, DifB, BlikR, BlikG, BlikB: real)

Параметры функции:

BodyR, BodyG, BodyB - компоненты цвета поверхности

DifR, DifG, DifB - возвращаемые функцией компоненты диффузного света

BlikR, BlikG, BlikB - возвращаемые функцией компоненты зеркального света

NX, NY, NZ - нормаль к поверхности, посчитанная в системе координат, где Z направлен по ходу луча, падающего на поверхность.

RX, RY, RZ - отражение падающего на поверхность луча (вектора наблюдателя).

MaxN - номер треугольника, на котором находится рассматриваемый пиксел.

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

Процедура Elluminaty перебирает все источники света и определяет вклад каждого источника в итоговую освещенность точки.

Для каждого источника она определяет его координаты во вспомогательной системе координат (T. x, T. y, T. z). Cоставляет матрицу преобразования к системе координат, в которой ось OZ совпадает с направлением на источник.

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

Функция Per2 возвращает коэффициент, равный доле света, который дошел от источника до точки. Она перебирает треугольники, лежащие между источником и рассматриваемой точкой. Если встречается непрозрачный треугольник, то функция завершается со значением 0. При встрече с прозрачным объектом результирующий коэффициент умножается на прозрачность треугольника.

3. Технологическая часть

3.1 Выбор языка программирования и обоснование выбора

Для разработки программы использовался язык Object Pascal. В качестве среды разработки использовалась среда Borland Delphi 6. Выбор именно Delphi 6 был обусловлен тем, что данная среда обладает богатыми возможностями для отладки программ, а так же тем, что она обладает большим собранием визуальных компонент, позволяющим сделать качественный интерфейс.

3.2 Модульная структура программы

Программа написана с помощью структурного подхода.

В ней содержится 6 модулей.4 модуля относятся к модулям форм. Это модули

Main, Model, Options, Figure. Главной формой в программе является Main. В поле этой формы выводится изображение, так же с помощью нее осуществляется вызов, других форм. Форма Model служит для преобразований над сценой, камерой, источниками света. С помощью формы Options осуществляется настройка опций рендеринга. Форма Figure, служит для параметрического задания первой модели. Так как главная форма управляет всеми остальными формами, то модули этих форм должны быть подключены к Main. Но помимо этого все подформы должны иметь доступ к главной форме. Поэтому я установил двустороннюю связь между Main и модулями подформ.

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

Модуль DataControl подключается к тем модулям, в которых осуществляется ввод данных с клавиатуры. Модуль содержит функции для проверки правильности ввода данных. Это функции IsReal и IsInt. Они проверяют, является ли переданная им строка числом с плавающей точкой или целым. Поскольку ввод данных осуществляется во всех формах кроме Main, то модуль подключен ко всем ним.

Полная структура связей выглядит так:

3.3 Интерфейс программы

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

Меню главной формы состоит из пунктов: "Преобразования", "Фигура", "Опции", "Рендеринг".

При выборе пункта "Рендеринг", произойдет просчет изображения с текущими параметрами источников света, камеры и вывод его на экран.

При выборе пункта "Фигура", можно выбрать из развернувшегося списка сцену, которую необходимо изобразить.

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

При выборе пункта "Опции" на экране, так же появится форма, в которой можно менять параметры рендеринга.

Форма делится на три части:

верхняя часть отвечает за преобразования над сценой.

средняя отвечает за преобразования над камерой.

нижняя позволяет управлять источниками освещения.

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

Сцена.

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

Повороты возможны вокруг осей OX, OY, OZ. Для каждой оси выделена своя строчка с полем для ввода угла и двумя кнопками. Для поворота необходимо ввести в поле, соответствующее нужной оси, угол поворота и нажать на одну из двух кнопок. При нажатии на "+" поворот осуществится на положительный угол, при нажатии на "-" - на отрицательный угол.

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

Камера.

Для камеры тоже существует переключатель между режимами сдвига и поворота. Повороту соответствует поднятый флажок.

Для поворота камеры вверх или вниз необходимо ввести в первое поле угол поворота и нажать "+" для поворота вверх или "-" для поворота вниз.

Для поворота камеры влево или вправо необходимо ввести угол поворота во второе поле и нажать "+" для поворота вправо или "-" для поворота влево.

Для поворота камеры по часовой стрелке или против часовой действуем аналогично, для поворота по часовой надо нажать "+", против часовой "-".

Перемещению камеры соответствует опущенный флажок.

Все аналогично операциям со сценой. Первая строчка позволяет двигаться вверх-вниз, вторая - вправо-влево, а третья - взад-вперед.

Источники света.

Управлять источниками света позволяет третья часть формы. На форме есть две кнопки со стрелками. Эти кнопки предназначены для того, чтобы листать источники света. По мере просмотра источников можно видеть все их характеристики (на форме расположены поля, в которых высвечиваются координаты, тип, интенсивность по трем составляющим и другие характеристики источника). Все параметры источника можно менять. Чтобы это сделать необходимо:

Найти с помощью стрелок нужный источник.

В поля характеристик источника ввести нужные значения.

Нажать кнопку изменить.

Обязательными параметрами для всех источников являются координаты и интенсивность по 3 составляющим. Если источник является обычным (т.е. светящим во все стороны), то этих параметров достаточно. Если же источник светит в конусе, то необходимо ввести направление и угол, в каком светит источник. Если источник второго типа, то можно поставить галочку в пункте затухание. В этом случае интенсивность источника будет максимальна на образующей конуса и постепенно снижается ближе к краям. На краях интенсивность равна нулю.

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

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

Включено, сглаживание или нет.

Включены тени или нет.

Включено ли качественное моделирование эффекта преломления

Глубина трассировки

Разрешение получаемого изображения

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

Все эти настройки можно менять. Для этого нужно ввести в поля новые значения и нажать "OK".

4. Экспериментально-исследовательская часть

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

Тест № 1

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

Сразу видно, что при пустой сцене рендеринг занимает 1.43 секунды. Это время, затраченное на перебор всех пикселов и вызов функций обработки пикселов. Оно не включает время на поиск пересечения луча с треугольником. График же похож на график логарифма. Вид графика подтверждается теорией. Алгоритм не перебирает все треугольники, а делает всего несколько сравнений со сферическими оболочками, число которых равно примерно логарифму по основанию 8 от числа треугольников.

Тест № 2

Рассмотрим теперь сцену, в которой полученное изображение занимает не весь экран. Для этого я буду удалять камеру от сцены. При очень сильном удалении сцена будет занимать всего 1 пиксел. Сцена будет использоваться одна и та же и будет состоять из 40000 треугольников. Проведем зависимость времени рендеринга от площади, занимаемой изображением на экране. Картинка будет иметь размеры 800х600.

График представляет собой линейную зависимость. Если луч не пересекает сцену, то после сравнения с самой большой оболочкой для многих лучей сразу же будет установлено, что они ничего не пересекают. Поэтому время рендеринга будет определяться временем на обработку лучей, пересекающих сцену. А поскольку время (t) обработки каждого из лучей в среднем одинаковое. То время будет носить зависимость n*t, где n - число пикселей. А это линейная функция от n.

Тест № 3

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

Видно, что при хаотичном разбросе треугольников, алгоритм работает куда медленнее. Улучшение есть, но оно очень мало. При хаотичном разбросе количество оболочек очень большое и поиск по ним становится, сравним по скорости с перебором объектов. По сравнению с первым тестом получилось замедление примерно в 100 раз. Алгоритм хорошо справляется с реальными объектами, так как они чаще всего хорошо сгруппированы.

Заключение

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

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

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

Помимо этого в программе:

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

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

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

Существует множество вариантов модернизации программы.

Модернизация алгоритма трассировки:

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

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

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

редактировать геометрию объектов на сцене

редактировать свойства материалов

добавлять и удалять объекты сцены.

управлять источниками света.

загружать и сохранять сцены. Для этого понадобится введение специального формата.

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

Список литературы

1. Роджерс Д. Алгоритмические основы машинной графики: пер. с англ. - М.: Мир, 1989.

2. Порев В.Н. Компьютерная графика 2004.

3. Тихомиров Ю. Программирование трехмерной графики. - СПб.: ИРМ - Санкт-Петербург, 1998.

4. Гантмахер Ф.Р. Теория матриц. - М.: Наука, 1967.