Алгоритм 3D-реконструкции под Android device (Algorithm of 3D reconstruction for Android device)

Автор: | 05.02.2018

Введение
Формализация задачи
Определение соответствия между координатами точек изображений
Алгоритм 3D-реконструкции
Тестирование алгоритма
Выводы
Программные реализации
Полезные ссылки

Введение

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

Трехмерная реконструкция (англ. 3D reconstruction) – это процесс получения 3D объектов на основе изображений.  На вход алгоритма обработки подается набор из нескольких изображений (два или более), результатом работы которого является 3D фотография (рис.1).

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

 Формализация задачи

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

Оси системы координат (СК) устройства (рис. 3) исходят из центра устройства (камеры). Вращение определяется углами Эйлера (alpha  beta  gamma), которые представляют собой разницу в градусах между координатами устройства и земными координатами.

Угол alpha=0°, когда ось Y устройства направлена на север. Углы beta=0° и gamma=0°, когда плоскость устройства параллельна земной поверхности. При повороте устройства вокруг соответствующей оси (см. рис.3) значение углов alpha и gamma  увеличивается в направлении против часовой стрелки,  bета– по часовой стрелке.

На рисунке 4 представлена схема для определения математической зависимости между координатами точек 2-х изображений. Точки описаны в плоскости XY (xy) главного и вспомогательного изображений, ось Z (z) направлена перпендикулярно к плоскости камеры и проходит через ее центр. Главным считается то изображение, фото которого получено первым. Значение углов определяется разницей (дельта) углов вспомогательного положения камеры относительно главного:

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

Определение соответствия между координатами точек изображений

Положение вспомогательной СК xyz относительно главной XYZ задается композицией 3-х матричных преобразований – вращения относительно оси Z, вращения относительно оси X и вращения относительно оси Y:

Матрицы записаны так, что положительное значение углов alpha и  gamma  — в направлении против часовой стрелки, bета – по часовой стрелке.

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

Результат перемножения 3-х обратных матриц и матрицы ортогонального преобразования дает следующий результат:

Преобразование координат точек из главной СК  во вспомогательную СК   в матричном представлении имеет вид:

В итоге, уравнения преобразования точек пространственного объекта в  точки текущего изображения запишутся в виде:

Алгоритм 3D-реконструкции

Система уравнений (1) позволяет найти координату Z каждой точки (X,Y) на главном изображении, если известны углы и координаты соответствующей точки (x,y) на вспомогательном изображении. Соответствие между точками может устанавливаться по их цвету.

Краткое описание алгоритма:

  • Для каждой точки главного изображения (X,Y) определяем соответствующую точку на вспомогательном изображении (x,y), осуществляя итерацию координаты Z
  • Соответствующая точка найдена, если совпадают 3 цвета точек (RGB) на обоих изображениях.
  • Вместе с соответствующей точкой определяется и координата Z.

Детальное описание алгоритма:

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

  • перевести координаты точки из пространства фото (I,J) в пространство логических координат (XY)
  • задать координату Z (Z=Z+dZ)
  • выполнить пересчет координат точки из главной СК (X,Y,Z) во вспомогательную СК (x,y)
  • перевести координаты точки из пространства логических координат (x,y) в пространства фото (I1,J1)
  • проверить на соответствие по цвету точки главного изображения (I,J) и вспомогательного (I1,J1)

Цифровое изображение хранится в виде прямоугольной матрицы, элементы  которой несут информацию о цвете элементарных участков изображения (пикселов), а номера   строки I и столбца J элемента определяют его положение в матрице.

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

Левая СК принята при записи изображений в файл во всех форматах и используется в большинстве программ по обработке изображений. В фотограмметрии традиционно применяется правая СК для анализа снимка, и многие современные цифровые фотограмметрические системы используют именно эту систему координат. Для перехода из одной СК в другую достаточно выполнить преобразование J = H — J.

В логической системе координат (рис.8) вершины нормализуются, т.е. удовлетворяют следующим условиям:

  • центр координат перенесен в центр экрана;
  • направление оси Y вверх;
  • координаты (от -1 до +1) соотнесены с шириной (Width) и высотой окна (Height) .

Ниже приведены соотношения для перехода из СК цифрового изображения в нормализованную СК и обратно.

X= 2 * I / Width  - 1
Y= 1 -  2 * J / Height         для левой СК       Y= 2 * J / Height -1
0 < Z < 1
I1 = Width * (x + 1) / 2
J1 = - Height * (y - 1) / 2    для левой СК       J1= Height * (y - 1) / 2

Определяются разницы сигналов между предполагаемыми соответствующими точками:

dR = abs(rgb1[I1][ J1].rgbRed - rgb[I][ J].rgbRed);
dG = abs(rgb1[I1][ J1].rgbGreen - rgb[I][ J].rgbGreen);
dB = abs(rgb1[I1][ J1].rgbBlue - rgb[I][ J].rgbBlue);

Соответствие между точками установлено, если выполняется комбинированное условие:

(dR < dmin и dG < dmin и dB < dmin)

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

Рассмотренный выше алгоритм 3D реконструкции по изображениям применим для случая, когда известно взаимное положение камер в системе 3D видеонаблюдения (рис.9).

Тестирование алгоритма

Тест №1

Алгоритм решения задачи 3D реконструкции  тестируется на следующем примере (рис.5).

Тестирование проводится в 2 этапа:

  1. Для тестирования генерируются два изображения цветных треугольников при помощи WinAPI OpenGL приложения photo3. Изображения создаются под заданными в программе углами alpha, bета и  gamma вспомогательной СК относительно главной СК. Для первого фото углы alpha=0, bета=0 и  gamma=0  Каждое изображение сохраняется через клавишу  prtsc в буферной памяти. Затем обрабатывается в графическом редакторе Paint — полностью вырезается квадратная часть окна, остаются только белый фон и треугольник. Сохраняются оба изображения как 24-разрядный BMP-файл размером 500*500 пикселов.
  2.  WinAPI приложение  angle3 загружает подготовленные на 1-м этапе изображения и  по ним реконструирует 3D объект в виде облака точек.

                    

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

Тест №2

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

Ниже приводится результат (рис. 10 -12) использования программы 3D реконструкции детской игрушки по 2 фотографиям.  Фотографии были получены под углами поворота камеры: alpha = 0.0, beta = 30.0 и  gamma = 20.0.  Углы определялись с погрешностью в пределах 5 градусов.

Тестовые фото игрушки «Незнайка» можно скачать здесь Data.

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

Тест №3

Для дальнейшего тестирования задачи было разработано приложение под мобильное устройство с операционной системой  Android, которая загружает фотографии и по ним реконструирует 3D объект в виде облака точек. Результаты приведены (рис. 13 — 14).

Программная реализация Android приложения на этом сайте не приводится. Есть только лишь заготовка для такого приложения (см. здесь). Можете попробовать самостоятельно реализовать такое приложение.

Выводы

  1. В работе описан и программно реализован (на С++ и Java) под операционную систему Android алгоритм 3D-реконструкции по 2-м изображениям. Он построен на основе решения одной из множества задач аффинных преобразований в однородных координатах. Алгоритм описан только для простейшего случая ортогонального проецирования, когда взаимное положение камер определяется только 3-мя угловыми параметрами. Это обусловлено тем, что пока в мобильных устройствах достаточно точно определяется только лишь ориентация устройства.
  2. Алгоритм может быть обобщен для случая, когда взаимное положение камер определяется большим количеством параметров, включая аппарат центрального проецирования и учитывая не только вращение, но и перемещение устройства в пространстве. Реализация такого алгоритма будет возможна, когда мобильные устройства смогут с достаточной точностью определять параметры перемещения и расстояния в пространстве. Недостаток параметров также может быть компенсирован за счет совершенствования алгоритма, например – путем нахождения на обоих фото особых точек и установления соответствия между ними на изображениях.
  3. Даже аппарат центрального проецирования есть всего лишь идеализированной моделью процесса формирования изображения в цифровой камере. Для того, чтобы максимально приблизить модель к реальному процессу, необходимо решить ряд задач предварительной обработки изображений, отделения фона, учесть потери качества, связанные с сжатием изображений, т.е., – решить комплекс задач компьютерного зрения.
  4. Программа, в которой реализован алгоритм, восстанавливает по фотографиям лишь облако точек пространственного объекта. Однако, она может быть легко модифицирована под поверхностную модель, например, методом триангуляции Делоне. Триангуляция – это разбиение геометрического объекта на треугольники, которые строятся по ближним точкам. При этом  цвет треугольника интерполируется по цвету этих точек (такие возможности обеспечивает графическая библиотека opengl, встроенная  в android).

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

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

Для тестирования алгоритма восстановления  3D-объекта по изображениям были разработаны 2 программы, которые обеспечивают:

  • генерацию рисунков 3D объектов (программа photo3)
  • загрузку  файлов изображений  и обработка точек изображений для определения координаты Z (программа angle3).

Для запуска приложения используйте Visual Studio не старше 17 версии. В 19 версии появляются ошибки, которые не нашел способ исправить.

Если возникают ошибки из-за преобразования символов, то в окне свойств проекта установите «Использовать набор символов Юникода»:

Для английской версии VS:

При возникновении ошибки типа error C4996: ‘fopen’: This function or variable may be unsafe   в свойствах проекта отключите предупреждения, генерируемые _CRT_SECURE_NO_DEPRECATE. Как это сделать, смотри здесь.

Программа photo3 для генерации изображений треугольников

Создается проект  WinAPI  приложения (File>New>Project>Visual C++>Other>Empty Project>OK). В проект добавляются (Add>New Item…) файлы с кодом, который приводится ниже.  В директорию, где расположены исходные файлы, забрасываются  файлы библиотеки  OpenGL.

Файлы библиотеки  OpenGL можно получить из скаченного архива glutdlls37beta.

Header Files

api.h

#ifndef __API__
#define __API__
#pragma comment( lib, "glut32.lib" ) // Подключается библиотека glut32.lib
#include <windows.h> // Заголовочный файл для использования функций Windows
#include "glut.h" // Заголовочный файл библиотеки glut32
#endif
// Файлы glut.h, glut32.lib размещаются в директории проекта, 
// где и файлы-исходники (.cpp) и заголовочные файлы (.h)

engine.h

#include "api.h"
#ifndef __ENGINE
#define __ENGINE
class Engine {
 GLsizei Height, Width;
 GLvoid DrawModel(GLvoid); // Отрисовка модели
public:
 GLvoid Resize(GLsizei width, GLsizei height); // Функция, вызываемая при изменении размеров окна
 GLvoid Init(GLvoid); // Функция, для задания начальных параметров
 GLvoid Draw(GLvoid); // Отрисовка (render) сцены 
 int SetWindowPixelFormat(HDC hDC); //функция, которая устанавливает параметры контекста воспроизведения OpenGL 
};
#endif

Source Files

main.cpp

#include "engine.h"
HWND hWnd;
HGLRC hGLRC;
HDC hDC;
Engine *engine;
LRESULT CALLBACK WindowFunc(HWND hWnd,UINT msg,WPARAM wParam,LPARAM lParam)
{
float pos[4] = {3,3,3,1};
float dir[3] = {-1,-1,-1};
PAINTSTRUCT ps;
switch(msg)
{
case WM_CREATE:
engine = new Engine();
hDC = GetDC(hWnd); // получаем контекст устройства нашего окна
engine->SetWindowPixelFormat(hDC); // устанавливаем параметры контекста воспроизведения OpenGL
hGLRC = wglCreateContext(hDC); // создаем контекст воспроизведения OpenGL
wglMakeCurrent(hDC, hGLRC); // делаем его текущим
engine->Init();
break;
case WM_DESTROY:
if (hGLRC) 
{ // удаляем созданный выше контекст воспроизведения OpenGL
wglMakeCurrent(NULL, NULL);
wglDeleteContext(hGLRC);
} // освобождаем контекст устройства нашего окна
ReleaseDC(hWnd, hDC);
PostQuitMessage(0);
break;
case WM_PAINT:
BeginPaint(hWnd, &ps);
engine->Draw();
EndPaint(hWnd, &ps);
break;
case WM_SIZE:
engine->Resize( LOWORD(lParam), HIWORD(lParam) );
break; 
case WM_ERASEBKGND:
 return 1;
break;
default:
return DefWindowProc(hWnd,msg,wParam,lParam);
}
return 0;
}
int WINAPI WinMain(HINSTANCE hThisInst,HINSTANCE hPrevInst,LPSTR str,int nWinMode)
{
MSG msg;
WNDCLASS wcl;
wcl.hInstance=hThisInst;
wcl.lpszClassName = (LPCWSTR) "OpenGLWinClass";
wcl.lpfnWndProc = WindowFunc;
wcl.style = CS_OWNDC | CS_HREDRAW | CS_VREDRAW;
wcl.hIcon = NULL;
wcl.hCursor = LoadCursor(NULL,IDC_ARROW);
wcl.lpszMenuName = NULL;
wcl.cbClsExtra = 0;
wcl.cbWndExtra = 0;
wcl.hbrBackground = (HBRUSH)GetStockObject(WHITE_BRUSH);
RegisterClass(&wcl);
hWnd = CreateWindow((LPCWSTR) "OpenGLWinClass",L"Win32 OpenGL Template", 
WS_OVERLAPPEDWINDOW | WS_CLIPCHILDREN | WS_CLIPSIBLINGS,
100,
100,
516,
540,
HWND_DESKTOP, NULL,
hThisInst, NULL);
ShowWindow(hWnd,nWinMode);
UpdateWindow(hWnd);
while(1)
{
while( PeekMessage(&msg,NULL,0,0,PM_NOREMOVE) ) 
if(GetMessage(&msg,NULL,0,0))
{ 
TranslateMessage(&msg);
DispatchMessage(&msg);
}
else
return 0;
//display(); // функция display вызывается в бесконечном цикле, причем, когда в очереди сообщений окна нет ничего.
}
return 0;
}

engine.cpp

#include "engine.h"
#include <math.h>
GLvoid Engine::Resize(GLsizei width, GLsizei height){
 if (height == 0) 
 {
 height = 1; 
 }
 glViewport(0, 0, width, height); // Устанавливается область просмотра
 Height = height;
 Width = width; 
}
GLvoid Engine::Init(GLvoid){
 glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
 glEnable(GL_DEPTH_TEST); // Включается тест глубины
}
GLvoid Engine::Draw(GLvoid) 
{
 //double al = 0.0, bt = 0.0, gm = 0.0;
 double al = 0.0, bt = 30.0, gm = 20.0;
 
 al = 3.14*al / 180;
 bt = 3.14*bt / 180;
 gm = 3.14*gm / 180;
 glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); // Очищается буфер кадра и буфер глубины
 // Текущая матрица сбрасывается на единичную
 glMatrixMode(GL_PROJECTION); // Действия будут производиться с матрицей проекции
 glLoadIdentity();
 glOrtho(-1, 1, -1, 1, 3, 10); // Устанавливается ортогональная проекция 
 glMatrixMode(GL_MODELVIEW); // Действия будут производиться с матрицей модели
 glLoadIdentity(); 
 gluLookAt(5 * cos(bt)*sin(gm), 5 * sin(bt), 5 * cos(bt)*cos(gm), 0, 0, 0, -sin(al), cos(al), 0); //Определяется вид ( Mview )
 //gluLookAt (0, 0, 5, 0, 0, 0, 0, 1, 0); //Определяется вид ( Mview )
//рисуется модель 
 DrawModel(); 
SwapBuffers(wglGetCurrentDC());
}
GLvoid Engine::DrawModel(GLvoid) // Отрисовка модели
{
 glBegin(GL_TRIANGLES);
 glColor3f(1, 0, 0); // set vertex color to red
 glVertex3f(1, 0, 0);
 glColor3f(0, 1, 0); // set vertex color to green
 glVertex3f(0, 1, 0);
 glColor3f(0, 0, 1); // set vertex color to blue
 glVertex3f(0, 0, 1);
 glEnd();
 }
// код функции, которая устанавливает параметры контекста воспроизведения OpenGL. 
int Engine::SetWindowPixelFormat(HDC hDC)
{
int m_GLPixelIndex;
PIXELFORMATDESCRIPTOR pfd;
pfd.nSize = sizeof(PIXELFORMATDESCRIPTOR);
pfd.nVersion = 1;
pfd.dwFlags = PFD_DRAW_TO_WINDOW | 
PFD_SUPPORT_OPENGL | 
PFD_DOUBLEBUFFER;
pfd.iPixelType = PFD_TYPE_RGBA;
pfd.cColorBits = 32;
pfd.cRedBits = 8;
pfd.cRedShift = 16;
pfd.cGreenBits = 8;
pfd.cGreenShift = 8;
pfd.cBlueBits = 8;
pfd.cBlueShift = 0;
pfd.cAlphaBits = 0;
pfd.cAlphaShift = 0;
pfd.cAccumBits = 64; 
pfd.cAccumRedBits = 16;
pfd.cAccumGreenBits = 16;
pfd.cAccumBlueBits = 16;
pfd.cAccumAlphaBits = 0;
pfd.cDepthBits = 32;
pfd.cStencilBits = 8;
pfd.cAuxBuffers = 0;
pfd.iLayerType = PFD_MAIN_PLANE;
pfd.bReserved = 0;
pfd.dwLayerMask = 0;
pfd.dwVisibleMask = 0;
pfd.dwDamageMask = 0;
// Передается на рассмотрение системе, выбранный нами формат пикселей. Функция просматривает в контексте устройства - hdc наиболее подходящий формат пикселей и выбирает его
m_GLPixelIndex = ChoosePixelFormat( hDC, &pfd);
if(m_GLPixelIndex==0) // Let's choose a default index.
{
m_GLPixelIndex = 1; 
if(DescribePixelFormat(hDC,m_GLPixelIndex,sizeof(PIXELFORMATDESCRIPTOR),&pfd)==0)
return 0;
}
// установить формат пикселей в контексте устройства
if (SetPixelFormat( hDC, m_GLPixelIndex, &pfd)==FALSE)
return 0;
return 1;
}

Программа angle3

Создается проект  WinAPI  приложения (File>New>Project>Visual C++>Other>Empty Project>OK). В проект добавляются (Add>New Item…) файлы с кодом, который приводится ниже.

Тестовые рисунки размещаются в папке Data проекта приложения — там же, где и исходные файлы проекта. Они загружаются из функции WndProc в файле main.cpp.

Начальное положение ЛСК xyz относительно ГСК  XYZ  определяется через матрицы поворота относительно осей ГСК —  gamma (Y),  beta (X) и alpha (Z) устанавливаются в функциях класса Matrix:

  • SetTranslationMatrix_4();     // gamma (Y)
  • SetTranslationMatrix_5();     // beta (X)
  • SetTranslationMatrix_6();     // alpha (Z)

В функции Engine::GetZ (в файле engine.cpp) непосредственно реализуется алгоритм восстановления координаты Z. Здесь задаются углы alpha (Z),  beta (X) и gamma (Y), которые задействованы в системе уравнений (1). На точность 3D-реконструкции влияют 2 параметра:

  • dz = 0.01;  —  шаг итерации по координате Z;
  • dmin  —  точность совпадения по цветовым параметрам. Для теста c использованием  фото  треугольников значение этого параметра от 1 до 5.  Для теста c использованием   фото игрушки «Незнайка» значение этого параметра значительно больше (dmin = 15), поскольку углы были определены приближенно.

Значения углов должны соответствовать тем углам, под которым было получено 2-е изображение в программе photo3. Только в этом случае представляется возможность восстановить 3D модель треугольника по 2-м изображениям.

Тестовое фото triangle.bmp — при alpha=0, beta=30, gamma=20:

Тестовое фото triangle1.bmp — при alpha=0, beta=30, gamma=20:

Тестовые фото игрушки «Незнайка» можно скачать здесь Data .  Второе фото получено под углами gamma=30, beta=30, alpha=0.

После запуска программы последовательно нажимаются клавиши 1,2,3 и 4. В результате появляются

  1. Изображение 1-й картинки
  2. Изображение 2-й картинки
  3. Точки объекта при Z=0
  4. Точки объекта  с восстановленной координатой Z.

 

Header Files

action.h

#include "matrix.h"
#include "vec.h"
class Action {
public:
 vec old_mouse;
 Matrix CurrentMatrix;
 void InitAction(double x, double y);
 void Rotate(double x, double y);
 void Translate(double x, double y);
 void Transform_0();
 void Transform_3();
 void Transform_4();
 void Transform_5();
};

engine.h

#include <windows.h>
#include "action.h"
#include "viewport.h"
#include <stdio.h>
class Engine{
 //Matrix current_rot;
 Action *action;
 unsigned short read_u16(FILE *fp);
 unsigned int read_u32(FILE *fp);
 int read_s32(FILE *fp);
public:
 Viewport viewport;
 void Draw(HDC hdc);
 void SetAction(Action *_action);
 int LoadBMP(HWND hWnd, HBITMAP hBmp, HDC hCmpDC, BITMAP bm, LPCWSTR fileName);
 RGBQUAD ** MatrixBMP(const char *fname, int &mx, int &my);
 void Filter(RGBQUAD **rgb, RGBQUAD **rgb1, int mx, int my);
 void GetZ(RGBQUAD **rgb, RGBQUAD **rgb1, int mx, int my);
 void ImagePixel(HDC hdc, RGBQUAD **rgb,int mx, int my);
};

geometry.h

#ifndef _POINT 
 struct _Point
 {
 double x, y, z;
 };
 #define _POINT
#endif

matrix.h

#include "geometry.h"
#ifndef _MATRIX
class Matrix{
 double data[4][4];
public:
 Matrix();
 void SetUnit();
 void SetRotationMatrix(double alpha);
 void SetRotationMatrixbySinCos(double sinalpha, double cosalpha);
 void SetTranslationMatrix(double tx, double ty);
 void SetTranslationMatrix_0();
 void SetTranslationMatrix_3();
void SetTranslationMatrix_4();     // gamma (Y)
void SetTranslationMatrix_5();     // beta (X)
void SetTranslationMatrix_6();     // alpha (Z)
 void SetTranslationMatrix_7();
 void MultiplyMatrices(Matrix &right);
 void ApplyMatrixtoPoint(_Point &point);
};
#define _MATRIX
#endif
void SetW(bool _w);

vec.h

#include "math.h"
typedef double vec_float;
class vec
{
public:
 vec_float x, y;
 vec() {}
 vec(vec_float xx, vec_float yy)
 {
 x = xx;
 y = yy;
 }
 vec(const vec& vector){
 x = vector.x;
 y = vector.y;
 }
 inline void set(vec_float xx, vec_float yy)
 {
 x = xx;
 y = yy;
 }
 inline vec operator + (vec t) // сложение
 {
 return vec(x + t.x, y + t.y);
 }
 inline vec operator - (vec t) // вычитание
 {
 return vec(x - t.x, y - t.y);
 }
 inline vec operator * (vec_float t) // произведение на число
 {
 return vec(x * t, y * t);
 }
 inline vec_float operator * (vec t) // скалярное произведение
 {
 return x * t.x + y * t.y;
 }
 inline vec_float operator ^ (vec t) // длина результата векторного произведения с учетом направления
 {
 return x * t.y - y * t.x;
 }
 inline vec_float length() // длина вектора
 {
 return sqrt(x * x + y * y);
 //return x * x + y * y;
 }
 inline vec unit() // нормализация вектора
 {
 vec_float l = length();
 if (l == 0.0f) return vec(0.0f, 0.0f);
 return vec(x / l, y / l);
 }
 inline bool zero() // определяет нулевой ли вектор
 {
 return (x == 0.0f) && (y == 0.0f);
 }
 inline bool equals(vec t) // проверяет вектора на точное совпадение
 {
 return (x == t.x) && (y == t.y);
 }
};

viewport.h

#include "geometry.h"
#ifndef _VIEWPORT
class Viewport{
 int Margin;
 int Height, Width;
public:
 Viewport();
 void SetWindowSize(int _Width, int _Height);
 _Point T(_Point point);
 _Point T_inv(_Point point);
 void SetMargin(int _Margin = 0);
};
#define _VIEWPORT
#endif

Source Files

main.cpp

#include "engine.h"
const double PI = 3.141592653;
LRESULT CALLBACK WndProc(HWND, UINT, WPARAM, LPARAM);
char szClassName[] = "CG6";
char szWindowCaption[] = "CG #6 Mouse Tracking and Rotation";
////////////////////////////////////////////////////////////////////////////////////////////////////
int WINAPI WinMain(HINSTANCE hInstance, HINSTANCE hPrevInstance, LPSTR lpCmdLine, int nCmdShow)
{
 HWND hWnd; 
 MSG lpMsg;
 WNDCLASS wc;
 // Заполняем структуру класса окна
 wc.style = CS_HREDRAW | CS_VREDRAW;
 wc.lpfnWndProc = WndProc;
 wc.cbClsExtra = 0;
 wc.cbWndExtra = 0;
 wc.hInstance = hInstance;
 wc.hIcon = NULL;
 wc.hCursor = LoadCursor(NULL, IDC_ARROW);
 wc.hbrBackground = (HBRUSH)GetStockObject(WHITE_BRUSH);
 wc.lpszMenuName = NULL;
 wc.lpszClassName = (LPCWSTR)szClassName;
// Регистрируем класс окна
 if (!RegisterClass(&wc))
 {
 MessageBox(NULL, L"Cannot register class", L"Error", MB_OK);
 return 0;
 }
 // Создаем основное окно приложения
 hWnd = CreateWindow( 
 (LPCWSTR)szClassName, // Имя класса 
 L"Шаблон WinAPI приложения", // Текст заголовка 
 WS_OVERLAPPEDWINDOW, // Стиль окна 
 100, 100, // Позиция левого верхнего угла 
 516, 540, // Ширина и высота окна 
 (HWND) NULL, // Указатель на родительское окно NULL 
 (HMENU) NULL, // Используется меню класса окна 
 (HINSTANCE)hInstance, // Указатель на текущее приложение
 NULL ); // Передается в качестве lParam в событие WM_CREATE
 if (!hWnd) 
 {
 MessageBox(NULL, L"Cannot create main window", L"Error", MB_OK);
 return 0;
 }
// Показываем наше окно
 ShowWindow(hWnd, nCmdShow); 
 UpdateWindow(hWnd);
// Выполняем цикл обработки сообщений до закрытия приложения
 while (GetMessage(&lpMsg, NULL, 0, 0)) {
 TranslateMessage(&lpMsg);
 DispatchMessage(&lpMsg);
 }
return (lpMsg.wParam);
}

////////////////////////////////////////////////////////////////////////////////////////////////////
LRESULT CALLBACK WndProc(HWND hWnd, UINT messg, WPARAM wParam, LPARAM lParam)
{
 PAINTSTRUCT ps;
 static RECT Rect;
 static HDC hdc, hCmpDC;
 static HBITMAP hBmp;
 static Action *action;
 static Engine *engine;
 static BITMAP bm;
static bool r1 /*= false*/;
static bool r2 /*= false*/;
static bool r3 /*= false*/;
static bool r4 /*= false*/;
static bool w /*= false*/;
static int mx;
static int my;
static int mx1;
static int my1;
static int im;
static RGBQUAD **v1;
static RGBQUAD **v2;
switch (messg)
 {
 case WM_CREATE:
 engine = new Engine(); 
 action = new Action(); //вызыв. констр. Matrix и иниц. матрица (станов. единичной)
 engine->SetAction(action); // иниц.закрытой перем engine->action = action;
 r1 = true;
 r2 = true;
 r3 = true;
 r4 = true;
 w = true;
 SetW(w);
 im = 0;
break;
 case WM_PAINT:
 GetClientRect(hWnd, &Rect);
 hdc = BeginPaint(hWnd, &ps);
 SetBkColor(hdc, 0xEECCCC);
 // Создание нового контекста для двойной буфферизации
 hCmpDC = CreateCompatibleDC(hdc);
 hBmp = CreateCompatibleBitmap(hdc, Rect.right - Rect.left, Rect.bottom - Rect.top);
 SelectObject(hCmpDC, hBmp);
 // Закраска фоновым цветом
 LOGBRUSH br;
 br.lbStyle = BS_SOLID;
 //br.lbColor = 0xEECCCC;
 br.lbColor = 0xFFFFFF;
 HBRUSH brush;
 brush = CreateBrushIndirect(&br);
 FillRect(hCmpDC, &Rect, brush);
 DeleteObject(brush);
// Загрузка картинок
 if (im == 1){ engine->LoadBMP(hWnd, hBmp, hCmpDC, bm, TEXT("Data/triangle.bmp"));}
 if (im == 2){engine->LoadBMP(hWnd, hBmp, hCmpDC, bm, TEXT("Data/triangle1.bmp"));} 
 // Рисование
 if (im >= 3)
 {
 engine->Draw(hCmpDC);
 engine->ImagePixel(hCmpDC, v1, mx, my); 
 }
 // Вывод на экран
 SetStretchBltMode(hdc, COLORONCOLOR);
 BitBlt(hdc, 0, 0, Rect.right - Rect.left, Rect.bottom - Rect.top, hCmpDC, 0, 0, SRCCOPY);
DeleteDC(hCmpDC);
 DeleteObject(hBmp);
 hCmpDC = NULL;
 EndPaint(hWnd, &ps);
 break;
case WM_SIZE: // какое раньше событие WM_PAINT или WM_SIZE
 GetClientRect(hWnd, &Rect);
 engine->viewport.SetWindowSize(Rect.right - Rect.left, Rect.bottom - Rect.top);
 // x = false;
 break; //определение закрытых перем. viewport.Height, viewport.Width;
case WM_ERASEBKGND:
 return 1;
 break;
case WM_LBUTTONDOWN:
 _Point mouse_point;
 mouse_point.x = LOWORD(lParam); //Сохраним координаты курсора мыши в системе окна
 mouse_point.y = HIWORD(lParam);
 //преобразуем координаты курсора в логические
 mouse_point = engine->viewport.T_inv(mouse_point);
 // запоминаем координаты курсора в объекте old_mouse
 action->InitAction(mouse_point.x, mouse_point.y);
 //InvalidateRect(hWnd, NULL, FALSE);
 break;
case WM_MOUSEMOVE:
 if (UINT(wParam) & MK_LBUTTON) {
 _Point mouse_point;
 mouse_point.x = LOWORD(lParam); // Координаты в системе окна
 mouse_point.y = HIWORD(lParam);
 mouse_point = engine->viewport.T_inv(mouse_point);// Логические координаты 
 if (UINT(wParam) & MK_CONTROL) {
//Перемещает систему координат, опираясь на текущее и предыдущее положения мыши
 action->Translate(mouse_point.x, mouse_point.y);
 } else {
 action->Rotate(mouse_point.x, mouse_point.y);
 }
 InvalidateRect(hWnd, NULL, FALSE);
 }
 break;
 case WM_KEYDOWN:
 int KeyPressed;
 KeyPressed = int(wParam);
 if (KeyPressed == int('W')) {w = (w == true) ? false : true; SetW(w); }
 if (KeyPressed == int('0')) {action->Transform_0();}
 if (KeyPressed == int('1'))
 {
 // Считывание и формирование массива 1-го изображения (только 1 раз)
 if (r1 == true) {
 //int mx, my; // для размеров изображения
 v1 = engine->MatrixBMP("Data/triangle.bmp", mx, my); // выделяем память и читаем туда файл
 r1 = false;
 }
 // Загрузка 1-й картинки
 im = 1;
 }
 if (KeyPressed == int('2'))
 {
 // Считывание и формирование массива 2-го изображения (только 1 раз)
 if (r2 == true) {
 v2 = engine->MatrixBMP("Data/triangle1.bmp", mx1, my1); // выделяем память и читаем туда файл
 r2 = false;
 }
 // Загрузка 2-й картинки
 im = 2;
 }
 if (KeyPressed == int('3'))
 {
 if (r3 == true) {
 engine->Filter(v1, v2, mx, my); // фильтрация точек изображений и присвоение точкам 1-го z=0;
 action->Transform_4(); // Подготовка соответствующей матрицы преобразований
 r3 = false; 
 } 
 im = 3;
 }
 if (KeyPressed == int('4'))
 { 
 // определение соответствующих точек 2-х изображений и нахождение для каждой точки 3-й координаты 
 if (r4 == true) {
 engine->GetZ(v1, v2, mx, my); //Только один раз определяем координату Z для точек главного изображения
 delete v2; // массив точек 2-го изображения можно удалить
 r4 = false;
 }
 im = 4;
 }
 InvalidateRect(hWnd, NULL, FALSE);
 break;
case WM_DESTROY:
 PostQuitMessage(0);
 break;
default:
 return (DefWindowProc(hWnd, messg, wParam, lParam));
 }
return (0);
}

action.cpp

#include "action.h"
/*
 Сигнализирует о начале некоторого действия
*/
void Action::InitAction(double x, double y){
 old_mouse.set(x, y);
}
/*
 Поворачивает систему координат, опираясь на текущее
 и предыдущее положения мыши
*/
void Action::Rotate(double x, double y){
 vec new_mouse(x, y);
 vec_float sina, cosa;
 sina = old_mouse.unit() ^ new_mouse.unit();
 cosa = old_mouse.unit() * new_mouse.unit(); 
 Matrix Rot;
 Rot.SetRotationMatrixbySinCos(sina, cosa);
 CurrentMatrix.MultiplyMatrices(Rot);
 old_mouse = new_mouse;
}
/*
 Перемещает систему координат, опираясь на текущее
 и предыдущее положения мыши
*/
void Action::Translate(double x, double y){
 vec new_mouse(x, y);
 vec delta = new_mouse - old_mouse;
 Matrix Tr;
 Tr.SetTranslationMatrix(delta.x, delta.y);
 CurrentMatrix.MultiplyMatrices(Tr);
 old_mouse = new_mouse;
}
void Action::Transform_0(){
 Matrix Tr;
 Tr.SetTranslationMatrix_0();
 CurrentMatrix.MultiplyMatrices(Tr);
}
void Action::Transform_3(){
 Matrix Tr;
 //Tr.SetTranslationMatrix_3();
 //CurrentMatrix.MultiplyMatrices(Tr);
 Tr.SetUnit();
}
void Action::Transform_4(){
 Matrix Tr;
 Tr.SetTranslationMatrix_4();       //  gamma (Y)
 CurrentMatrix.MultiplyMatrices(Tr); 
 Tr.SetTranslationMatrix_5();       //  beta (X)
 CurrentMatrix.MultiplyMatrices(Tr); 
 Tr.SetTranslationMatrix_6();       //  alpha (Z)
 CurrentMatrix.MultiplyMatrices(Tr); 
}
void Action::Transform_5(){
 Matrix Tr;
 Tr.SetTranslationMatrix_7();
 CurrentMatrix.MultiplyMatrices(Tr);
 Tr.SetTranslationMatrix_6();
 CurrentMatrix.MultiplyMatrices(Tr);
}

engine.cpp

#include <windows.h>
#include "geometry.h"
#include "matrix.h"
#include "engine.h"
/*
 Привязываем к движку объект, который отвечает за действия пользователя
*/
void Engine::SetAction(Action *_action)
{
 action = _action;
}
/*
 Выводит графику на контекст hdc
*/
void Engine::Draw(HDC hdc){
 //_Point p[8];
 _Point p[12];
 //p[1].x = 0.5;
 //p[1].y = 0.0;
 // p[1].z = 0.0;
 // p[2].x = 0.0;
 // p[2].y = 0.5;
 // p[2].z = 0.0; // 0.5 = 1
 //p[3].x = 0.0;
 //p[3].y = 0.0;
 // p[3].z = 0.5;
 p[1].x = 1.0;
 p[1].y = 0.0;
 p[1].z = 0.0;
 p[2].x = 0.0;
 p[2].y = 1.0;
 p[2].z = 0.0; // 0.5 = 1
 p[3].x = 0.0;
 p[3].y = 0.0;
 p[3].z = 1.0; 
 // локальная система координат
p[4].x = 0.0;
 p[4].y = 0.0;
 p[4].z = 0.0;
p[5].x = 0.5;
 p[5].y = 0.0;
 p[5].z = 0.0;
p[6].x = 0.0;
 p[6].y = 0.5;
 p[6].z = 0.0;
p[7].x = 0.0;
 p[7].y = 0.0;
 p[7].z = 0.5;
for (int i = 1; i <8; i++){
 //Вращение и перемещение осуществляется в логических координатах
 action->CurrentMatrix.ApplyMatrixtoPoint(p[i]);
 // Переход из логических в оконные координаты.
 p[i] = viewport.T(p[i]);
 }
// треугольник
MoveToEx(hdc, p[1].x, p[1].y, NULL);
 LineTo(hdc, p[2].x, p[2].y);
 LineTo(hdc, p[3].x, p[3].y);
 LineTo(hdc, p[1].x, p[1].y);
// оси координат
int i = 4;
for (int j = 1; j <= 3; j++){
MoveToEx(hdc, p[i].x, p[i].y, NULL);
LineTo(hdc, p[i+j].x, p[i+j].y);
}
// глобальная система координат
p[8].x = 0.0;
p[8].y = 0.0;
p[8].z = 0.0;
p[9].x = 1.0;
p[9].y = 0.0;
p[9].z = 0.0;
p[10].x = 0.0;
p[10].y = 1.0;
p[10].z = 0.0;
p[11].x = 0.0;
p[11].y = 0.0;
p[11].z = 1.0;
for (int i = 1; i <12; i++){
 // Переход из логических в оконные координаты.
 p[i] = viewport.T(p[i]);
}
// оси координат
i = 8;
for (int j = 1; j <= 3; j++){
 MoveToEx(hdc, p[i].x, p[i].y, NULL);
 LineTo(hdc, p[i+j].x, p[i+j].y);
}
}
int Engine::LoadBMP(HWND hWnd, HBITMAP hBmp, HDC hCmpDC, BITMAP bm, LPCWSTR fileName){
 // Загрузка картинки
hBmp = (HBITMAP)LoadImage(NULL, fileName, IMAGE_BITMAP, 0, 0, LR_LOADFROMFILE | LR_CREATEDIBSECTION);
if (hBmp == NULL)
 {
 MessageBoxW(hWnd, TEXT("Файл не найден"), TEXT("Загрузка изображения"), MB_OK | MB_ICONHAND);
 DestroyWindow(hWnd);
 return 1;
 }
GetObject(hBmp, sizeof(bm), &bm);
 SelectObject(hCmpDC, hBmp);
 ReleaseDC(hWnd, hCmpDC);
 return 1;
}
void Engine::Filter(RGBQUAD **rgb, RGBQUAD **rgb1, int mx, int my){
 int i, j;
 double g[3][3]; // box filter
double div = 1.0; // коэффициент нормирования
 g[0][0] = 0.0; g[1][0] = 0.0; g[2][0] = 0.0;
 g[0][1] = 0.0; g[1][1] = 1.0; g[2][1] = 0.0;
 g[0][2] = 0.0; g[1][2] = 0.0; g[2][2] = 0.0;
 //double div = 9.0; // коэффициент нормирования
 //g[0][0] = 1.0; g[1][0] = 1.0; g[2][0] = 1.0;
 //g[0][1] = 1.0; g[1][1] = 1.0; g[2][1] = 1.0;
 //g[0][2] = 1.0; g[1][2] = 1.0; g[2][2] = 1.0;
//double div = 18.0; // коэффициент нормирования
//g[0][0] = 1.0; g[1][0] = 2.0; g[2][0] = 2.0;
 //g[0][1] = 2.0; g[1][1] = 4.0; g[2][1] = 2.0;
 //g[0][2] = 2.0; g[1][2] = 2.0; g[2][2] = 1.0;
//double div = 1.0; // коэффициент нормирования
 //g[0][0] = -1.0; g[1][0] = -1.0; g[2][0] = -1.0;
 //g[0][1] = -1.0; g[1][1] = 9.0; g[2][1] = -1.0;
 //g[0][2] = -1.0; g[1][2] = -1.0; g[2][2] = -1.0;
 RGBQUAD **rgbtmp = new RGBQUAD*[mx];
// Усреднение (сглаживание) цветов пикселов в соответствии с ближайшими пикселами 1-го рисунка
for (i = 0; i < mx; i++) {
 rgbtmp[i] = new RGBQUAD[my];
 }
for (j = 1; j < my - 1; j++) {
 for (i = 1; i < mx - 1; i++) {
 rgbtmp[i][j].rgbRed = (rgb[i - 1][j + 1].rgbRed*g[0][0] + rgb[i][j + 1].rgbRed*g[1][0] + rgb[i + 1][j + 1].rgbRed*g[2][0] +
 rgb[i - 1][j].rgbRed*g[0][1] + rgb[i][j].rgbRed*g[1][1] + rgb[i + 1][j].rgbRed*g[2][1] +
 rgb[i - 1][j - 1].rgbRed*g[0][2] + rgb[i][j - 1].rgbRed*g[1][2] + rgb[i + 1][j - 1].rgbRed*g[2][2]) / div;
 rgbtmp[i][j].rgbGreen = (rgb[i - 1][j + 1].rgbGreen*g[0][0] + rgb[i][j + 1].rgbGreen*g[1][0] + rgb[i + 1][j + 1].rgbGreen*g[2][0] +
 rgb[i - 1][j].rgbGreen*g[0][1] + rgb[i][j].rgbGreen*g[1][1] + rgb[i + 1][j].rgbGreen*g[2][1] +
 rgb[i - 1][j - 1].rgbGreen*g[0][2] + rgb[i][j - 1].rgbGreen*g[1][2] + rgb[i + 1][j - 1].rgbGreen*g[2][2]) / div;
 rgbtmp[i][j].rgbBlue = (rgb[i - 1][j + 1].rgbBlue*g[0][0] + rgb[i][j + 1].rgbBlue*g[1][0] + rgb[i + 1][j + 1].rgbBlue*g[2][0] +
 rgb[i - 1][j].rgbBlue*g[0][1] + rgb[i][j].rgbBlue*g[1][1] + rgb[i + 1][j].rgbBlue*g[2][1] +
 rgb[i - 1][j - 1].rgbBlue*g[0][2] + rgb[i][j - 1].rgbBlue*g[1][2] + rgb[i + 1][j - 1].rgbBlue*g[2][2]) / div;
 }
}
for (i = 1; i < mx - 1; i++) {
 rgb[i] = rgbtmp[i]; 
 }
// Усреднение (сглаживание) цветов пикселов в соответствии с ближайшими пикселами 2-го рисунка
 for (i = 0; i < mx; i++) {
 rgbtmp[i] = new RGBQUAD[my];
 }
for (j = 1; j < my; j++) {
 for (i = 1; i < mx; i++) {
 rgb[i][j].rgbReserved = 128; // z=0 rgb[i][j].rgbReserved используется для хранения координаты Z (0-255) 
 }
 }
for (j = 1; j < my - 1; j++) {
 for (i = 1; i < mx - 1; i++) {
 rgbtmp[i][j].rgbRed = (rgb1[i - 1][j + 1].rgbRed*g[0][0] + rgb1[i][j + 1].rgbRed*g[1][0] + rgb1[i + 1][j + 1].rgbRed*g[2][0] +
 rgb1[i - 1][j].rgbRed*g[0][1] + rgb1[i][j].rgbRed*g[1][1] + rgb1[i + 1][j].rgbRed*g[2][1] +
 rgb1[i - 1][j - 1].rgbRed*g[0][2] + rgb1[i][j - 1].rgbRed*g[1][2] + rgb1[i + 1][j - 1].rgbRed*g[2][2]) / div;
 rgbtmp[i][j].rgbGreen = (rgb1[i - 1][j + 1].rgbGreen*g[0][0] + rgb1[i][j + 1].rgbGreen*g[1][0] + rgb1[i + 1][j + 1].rgbGreen*g[2][0] +
 rgb1[i - 1][j].rgbGreen*g[0][1] + rgb1[i][j].rgbGreen*g[1][1] + rgb1[i + 1][j].rgbGreen*g[2][1] +
 rgb1[i - 1][j - 1].rgbGreen*g[0][2] + rgb1[i][j - 1].rgbGreen*g[1][2] + rgb1[i + 1][j - 1].rgbGreen*g[2][2]) / div;
 rgbtmp[i][j].rgbBlue = (rgb1[i - 1][j + 1].rgbBlue*g[0][0] + rgb1[i][j + 1].rgbBlue*g[1][0] + rgb1[i + 1][j + 1].rgbBlue*g[2][0] +
 rgb1[i - 1][j].rgbBlue*g[0][1] + rgb1[i][j].rgbBlue*g[1][1] + rgb1[i + 1][j].rgbBlue*g[2][1] +
 rgb1[i - 1][j - 1].rgbBlue*g[0][2] + rgb1[i][j - 1].rgbBlue*g[1][2] + rgb1[i + 1][j - 1].rgbBlue*g[2][2]) / div;
 }
 }
for (i = 1; i < mx - 1; i++) {
 rgb1[i] = rgbtmp[i];
 }
}
void Engine::GetZ(RGBQUAD **rgb, RGBQUAD **rgb1, int mx, int my){
 _Point p;
 _Point p1;
 double al = 0.0, bt = 30.0, gm = 20.0; 
 al = 3.14*al / 180;
 bt = 3.14*bt / 180;
 gm = 3.14*gm / 180;
double dz=0.01, pxy;
 int i, j, i1, j1, SUM, dmin;
dmin = 1.0;
 int dR, dG, dB;
for (j = 1; j < my; j++) {
 for (i = 1; i < mx; i++) {
 SUM = rgb[i][j].rgbRed + rgb[i][j].rgbGreen + rgb[i][j].rgbBlue;
 if (SUM < 765) { // Отфильтровываем светлые точки (255,255,255)
 // переходим от координат изображения (снизу вверх) к оконным координатам (сверху вниз)
 p.x = i;
 p.y = my - j;
rgb[i][j].rgbReserved = 127; // z=0 rgb[i][j].rgbReserved используется для хранения координаты Z (0-255) 
 p.z = mx *rgb[i][j].rgbReserved / 255;
 //Выполняется преобразование оконных координат в логические координаты x[-1, 1] y[-1, 1] z[-1, 1]
 p = viewport.T_inv(p);
 if (p.x*p.x + p.y*p.y + p.z*p.z < 1.0) // Отфильтровываем точки за пределами сферы единичного радиуса
 do {
p.z = p.z + dz;
 if (p.x*p.x + p.y*p.y + p.z*p.z >= 1.0){ p.z = 1.1; break; }; // Отфильтровываем точки за пределами сферы единичного радиуса
 p1.x = p.x*(cos(gm)*cos(al) - sin(gm)*sin(bt)*sin(al)) + p.y*cos(bt)*sin(al) - p.z*(sin(gm)*cos(al) + cos(gm)*sin(bt)*sin(al));
 p1.y = -p.x*(cos(gm)*sin(al) + sin(gm)*sin(bt)*cos(al)) + p.y*cos(bt)*cos(al) + p.z*(sin(gm)*sin(al) - cos(gm)*sin(bt)*cos(al));//p1.z = 0.0; в данном случае координата z не имеет значения
 p1 = viewport.T(p1); // Переход из логических в оконные координаты.
 // Переход от оконных координат к пикселам изображения (снизу вверх)
 i1 = p1.x;
 j1 = my - p1.y;
 dR = abs(rgb1[i1][j1].rgbRed - rgb[i][j].rgbRed);
 dG = abs(rgb1[i1][j1].rgbGreen - rgb[i][j].rgbGreen);
 dB = abs(rgb1[i1][j1].rgbBlue - rgb[i][j].rgbBlue);
} while ((dR > dmin || dG > dmin || dB > dmin) && p.z <= 1.0);
 if (p.z <= 1) {
 p = viewport.T(p); // Переход из логических в оконные координаты (главное изображение)
 rgb[i][j].rgbReserved = p.z * 255 / mx; // запоминаем координату z в байтах (0-255)
 }
 else
 {
 rgb[i][j].rgbRed = 255; rgb[i][j].rgbGreen = 255; rgb[i][j].rgbBlue = 255;
 };
 }
}
}
}void Engine::ImagePixel(HDC hdc, RGBQUAD **rgb, int mx, int my) {
COLORREF clr;
 _Point p;
 int i, j, SUM;
 //int dR, dG, dB;
 //выводим результат
 for (j = 0; j < my; j++) {
 for (i = 0; i < mx; i++) {
 SUM = rgb[i][j].rgbRed + rgb[i][j].rgbGreen + rgb[i][j].rgbBlue;
 if (SUM < 755) { // Отфильтровываем светлые точки (255,255,255)
 p.x = i;
 p.y = my - j;
//rgb[i][j].rgbReserved = 127; // z=0
 p.z = mx *rgb[i][j].rgbReserved / 255;
 //if (p.x == mx / 2 && p.y == mx / 2) { p.z = mx; }; // рисуется точка
 //Выполняется преобразование из координат экрана в логические координаты x[-1, 1] y[-1, 1] z[-1, 1]
 p = viewport.T_inv(p);
// 
 //p.z = 0.5;
 // Блок выполняется при каждой перерисовке
 //Вращение и перемещение осуществляется в логических координатах
 action->CurrentMatrix.ApplyMatrixtoPoint(p);
 // Переход из логических в оконные координаты.
 p = viewport.T(p);
 // Рисуем точку
 clr = RGB(rgb[i][j].rgbRed, rgb[i][j].rgbGreen, rgb[i][j].rgbBlue);
 SetPixel(hdc, p.x, p.y, clr);
 }
 }
 } 
}
RGBQUAD ** Engine::MatrixBMP(const char *fname, int &mx, int &my)
{
 FILE * pFile = fopen(fname, "rb");
// считываем заголовок файла
 BITMAPFILEHEADER header;
header.bfType = read_u16(pFile);
 header.bfSize = read_u32(pFile);
 header.bfReserved1 = read_u16(pFile);
 header.bfReserved2 = read_u16(pFile);
 header.bfOffBits = read_u32(pFile);
// считываем заголовок изображения
 BITMAPINFOHEADER bmiHeader;
bmiHeader.biSize = read_u32(pFile);
 bmiHeader.biWidth = read_s32(pFile);
 bmiHeader.biHeight = read_s32(pFile);
 bmiHeader.biPlanes = read_u16(pFile);
 bmiHeader.biBitCount = read_u16(pFile);
 bmiHeader.biCompression = read_u32(pFile);
 bmiHeader.biSizeImage = read_u32(pFile);
 bmiHeader.biXPelsPerMeter = read_s32(pFile);
 bmiHeader.biYPelsPerMeter = read_s32(pFile);
 bmiHeader.biClrUsed = read_u32(pFile);
 bmiHeader.biClrImportant = read_u32(pFile);
 RGBQUAD **rgb = new RGBQUAD*[bmiHeader.biWidth];
 for (int i = 0; i < bmiHeader.biWidth; i++) {
 rgb[i] = new RGBQUAD[bmiHeader.biHeight];
 }
int kr = (int)bmiHeader.biWidth * 3 % 4;
 if (kr != 0) {kr = 4 - kr; 
}
for (int j = 0; j < bmiHeader.biHeight; j++) {
 for (int i = 0; i < bmiHeader.biWidth; i++) {
 rgb[i][j].rgbBlue = getc(pFile);
 rgb[i][j].rgbGreen = getc(pFile);
 rgb[i][j].rgbRed = getc(pFile);
 }
// пропускаем последние 1-3 байта в конце строки для обеспечения кратности 4
 for (int i = 0; i < kr; i++) {
 getc(pFile);
 }
}
// Передаем значения ширины и высоты глобальным переменным
 mx = bmiHeader.biWidth;
 my = bmiHeader.biHeight;
fclose(pFile);
 return rgb;
}
unsigned short Engine::read_u16(FILE *fp)
{
 unsigned char b0, b1;
b0 = getc(fp);
 b1 = getc(fp);
return ((b1 << 8) | b0);
}
unsigned int Engine::read_u32(FILE *fp)
{
 unsigned char b0, b1, b2, b3;
b0 = getc(fp);
 b1 = getc(fp);
 b2 = getc(fp);
 b3 = getc(fp);
return ((((((b3 << 8) | b2) << 8) | b1) << 8) | b0);
}
int Engine::read_s32(FILE *fp)
{
 unsigned char b0, b1, b2, b3;
b0 = getc(fp);
 b1 = getc(fp);
 b2 = getc(fp);
 b3 = getc(fp);
return ((int)(((((b3 << 8) | b2) << 8) | b1) << 8) | b0);
}

matrix.cpp

#include <math.h>
#include <memory.h>
#include "matrix.h"
#include "geometry.h"
bool w;
void SetW(bool _w) { w = _w; }
/*
 Конструктор, инициализирует матрицу единичной (матрица тождественного преобразования)
*/
Matrix::Matrix() {
 SetUnit();
}
/*
 Текущая матрица становится единичной
*/
void Matrix::SetUnit() {
 //memset(data, sizeof(data), 0);
 for (int i = 0; i < 4; i++) {
 for (int j = 0; j < 4; j++) {
 data[i][j] = 0.0;
 }
 }
 for (int i = 0; i < 4; i++) {
 data[i][i] = 1.0;
 }
}
/*
 Устанавливает текущую матрицу в качестве матрицы вращения на угол alpha,
 заданный косинусом и синусом
*/
void Matrix::SetRotationMatrixbySinCos(double sinalpha, double cosalpha) {
 SetUnit();
 /*data[0][0] = cosalpha;
 data[0][2] = sinalpha;
 data[2][0] = -sinalpha;
 data[2][2] = cosalpha;*/
 data[0][0] = cosalpha; // a // a c p 0 //
 //data[0][1] = 1.0; // c // b d q 0 //
 data[0][2] = sinalpha; // p // h f r 0//
 //data[0][3] = 1.0; // 0 // m n l 1 //
 //data[1][0] = 0.1; // b 
 //data[1][1] = 0.0; // d
 //data[1][2] = 0.1; // q
 //data[1][3] = 0.1; // 0
 data[2][0] = -sinalpha; // h
 //data[2][1] = 0.25; // f
 data[2][2] = cosalpha; // r
 //data[2][3] = 1.0; // 0 //w= data[2][3] * _z + 1 data[2][3] = 1/F F =1
 //data[3][0] = 0.25; // m
 //data[3][1] = 0.25; // n
 //data[3][2] = 1.0; // l
 //data[3][3] = 1.0; // 1
}
/*
Устанавливает текущую матрицу в качестве матрицы вращения на угол alpha
*/
void Matrix::SetRotationMatrix(double alpha) {
 SetRotationMatrixbySinCos(sin(alpha), cos(alpha));
}
/*
 Устанавливает текущую матрицу в качестве матрицы параллельного переноса на вектор (tx, ty)
*/
void Matrix::SetTranslationMatrix(double tx, double ty) {
 SetUnit();
 double txx = tx * tx;
 double txy = tx * ty;
 double tyy = ty * ty;
 double tsx = sqrt(1 - txx);
 double tsy = sqrt(1 - tyy);
 data[0][0] = tsx; // a // a c p 0 //
 data[0][1] = -txy; // c // b d q 0 //
 data[0][2] = tx * tsy; // p // h f r 0 //
 //data[0][3] = 1.0; // 0 // m n l 1 //
 //data[1][0] = 0.0; // b 
 data[1][1] = tsy; // d 
 data[1][2] = ty; // q 
 //data[1][3] = 0.0; // 0
 data[2][0] = -tx; // h
 data[2][1] = -tsx * ty; // f 
 data[2][2] = tsx * tsy; // r 
 //data[2][3] = 0.0; // 0
 //data[3][0] = 0.25; // m
 //data[3][1] = 0.25; // n
 //data[3][2] = 1.0; // l
 //data[3][3] = 1.0; // 1
}
void Matrix::SetTranslationMatrix_0() {
 SetUnit();
 //data[0][0] = 0.0; // a // a c p 0 //
 //data[0][1] = 1.0; // c // b d q 0 //
 //data[0][2] = 0.0; // p // h f r 0//
 //data[0][3] = 1.0; // 0 // m n l 1 //
 //data[1][0] = 0.1; // b 
 //data[1][1] = 0.0; // d
 //data[1][2] = 0.1; // q
 //data[1][3] = 0.1; // 0
 //data[2][0] = 0.5; // h
 //data[2][1] = 0.25; // f
 data[2][2] = 0.0; // r
 //data[2][3] = -1.0; // 0 //w= data[2][3] * _z + 1 data[2][3] = 1/F F =1
 //data[3][0] = 0.25; // m
 //data[3][1] = 0.25; // n
 //data[3][2] = 1.0; // l
 //data[3][3] = 1.0; // 1
}
void Matrix::SetTranslationMatrix_3() {
 SetUnit(); // единичная матрица
}

void Matrix::SetTranslationMatrix_4() { //gm (Y)
 SetUnit();
 double gm = 20.0;
 gm = 3.14*gm / 180;
 data[0][0] = cos(gm); // a     // a c p 0 //
 //data[0][1] = -0.707; // c    // b d q 0 //
 data[0][2] = sin(gm); // p     // h f r 0 //
 //data[0][3] = 1.0; // 0       // m n l 1 //
 //data[1][0] = 0.707; // b 
 //data[1][1] = 0.707; // d
 //data[1][2] = -1.0; // q
 //data[1][3] = 0.1; // 0
 data[2][0] = -sin(gm); // h      угол 20 (Y)
 //data[2][1] = 1.0; // f
 data[2][2] = cos(gm); // r
 //data[2][3] = 0.0; // 0
 //data[3][0] = 0.25; // m
 //data[3][1] = 0.25; // n
 //data[3][2] = 1.0; // l
 //data[3][3] = 1.0; // 1
}

void Matrix::SetTranslationMatrix_5() { // bt (X) 
 SetUnit();
 double bt = 30.0;
 bt = 3.14*bt / 180;
 //data[0][0] = 0.0; // a     // a c p 0 //
 //data[0][1] = 1.0; // c     // b d q 0 //
 //data[0][2] = 0.0; // p     // h f r 0 //
 //data[0][3] = 1.0; // 0     // m n l 1 //
 //data[1][0] = 0.1; // b 
 data[1][1] = cos(bt); // d 
 data[1][2] = sin(bt); // q 
 //data[1][3] = 0.1; // 0      угол 30 (X)
 //data[2][0] = 0.5; // h
 data[2][1] = -sin(bt); // f 
 data[2][2] = cos(bt); // r 
 //data[2][3] = 0.0; // 0
 //data[3][0] = 0.25; // m
 //data[3][1] = 0.25; // n
 //data[3][2] = 1.0; // l
 //data[3][3] = 1.0; // 1
}

void Matrix::SetTranslationMatrix_6() { // al (Z)
 SetUnit();
 double al = 0.0;
 al = 3.14*al / 180;
 data[0][0] = cos(al); // a    // a c p 0 //
 data[0][1] = -sin(al); // c   // b d q 0 //
 //data[0][2] = -0.707; // p   // h f r 0 //
 //data[0][3] = 1.0; // 0      // m n l 1 //
 data[1][0] = sin(al); // b 
 data[1][1] = cos(al); // d
 //data[1][2] = -1.0; // q
 //data[1][3] = 0.1; // 0 угол 30 (Z)
 //data[2][0] = 0.707; // h 
 //data[2][1] = 1.0; // f
 //data[2][2] = 0.707; // r
 //data[2][3] = 0.0; // 0
 //data[3][0] = 0.25; // m
 //data[3][1] = 0.25; // n
 //data[3][2] = 1.0; // l
 //data[3][3] = 1.0; // 1
}

void Matrix::SetTranslationMatrix_7() {
 SetUnit();
 //data[0][0] = 0.0; // a // a c p 0 //
 //data[0][1] = 1.0; // c // b d q 0 //
 //data[0][2] = 0.0; // p // h f r 0//
 //data[0][3] = 1.0; // 0 // m n l 1 //
 //data[1][0] = 0.1; // b 
 data[1][1] = 0.866; // d 
 data[1][2] = -0.5; // q 
 //data[1][3] = 0.1; // 0 угол 30
 //data[2][0] = 0.5; // h
 data[2][1] = 0.5; // f 
 data[2][2] = 0.866; // r 
 //data[2][3] = 0.0; // 0
 //data[3][0] = 0.25; // m
 //data[3][1] = 0.25; // n
 //data[3][2] = 1.0; // l
 //data[3][3] = 1.0; // 1
}
/*
 Умножает текущую матрицу на матрицу, переданную в качестве параметра
*/
void Matrix::MultiplyMatrices(Matrix &right) {
 double temp[4][4];
 double val;
 memcpy(temp, data, sizeof(data));
 for (int i = 0; i < 4; i++) {
 for (int j = 0; j < 4; j++) {
 val = 0;
 for (int v = 0; v < 4; v++) {
 /*val += temp[i][v] * right.data[v][j];*/
 if (w == true) val += temp[i][v] * right.data[v][j]; // глобальная система координат
 else val += right.data[i][v] * temp[v][j]; // локальная система координат
 }
 data[i][j] = val;
 }
 }
}
/*
 Перемножает точку, переданную в качестве параметра на текущую матрицу.
 При этом последний столбец матрицы не учитывается
*/
void Matrix::ApplyMatrixtoPoint(_Point &point) {
 double _x, _y, _z, w;
 _x = point.x;
 _y = point.y;
 _z = point.z;
 point.x = _x * data[0][0] + _y * data[1][0] + _z * data[2][0] + data[3][0];
 point.y = _x * data[0][1] + _y * data[1][1] + _z * data[2][1] + data[3][1];
 point.z = _x * data[0][2] + _y * data[1][2] + _z * data[2][2] + data[3][2];
 w = data[2][3] * _z + 1;
 point.x = point.x / w;
 point.y = point.y / w;
}

viewport.cpp

#include "viewport.h"
/*
 Задает размер окна
*/
void Viewport::SetWindowSize(int _Width, int _Height){
 Width = _Width;
 Height = _Height;
}
/*
 Выполняет преобразование из координат [-1, 1] x [-1, 1]
 в координаты экрана, с учетом отступов
*/
_Point Viewport::T(_Point point){
 _Point TPoint;
 TPoint.x = Margin + (1.0 / 2) * (point.x + 1) * (Width - 2 * Margin);
 TPoint.y = Margin + (-1.0 / 2) * (point.y - 1) * (Height - 2 * Margin);
 TPoint.z = Margin + (1.0 / 2) * (point.z + 1) * (Width - 2 * Margin);
 return TPoint;
}
/*
 Выполняет преобразование из координат экрана, в координаты
 [-1, 1] x [-1, 1] с учетом отступов
*/
_Point Viewport::T_inv(_Point point){
 _Point TPoint;
 TPoint.x = double(point.x - Margin) / (1.0 / 2) / (Width - 2 * Margin) - 1;
 TPoint.y = double(point.y - Margin) / (-1.0 / 2) / (Height - 2 * Margin) + 1;
 TPoint.z = double(point.z - Margin) / (1.0 / 2) / (Width - 2 * Margin) - 1;
 return TPoint;
}
/*
 Конструктор, выставляет отступы по умолчанию
*/
Viewport::Viewport(){
 SetMargin();
}
/*
 Устанавливает отступ по краям экрана
*/
void Viewport::SetMargin(int _Margin){
 Margin = _Margin;
}

 

Полезные ссылки:

 

Автор: Николай Свирневский

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *