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

Автор: | 05.02.2018

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

Введение

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

изображения запишутся в виде:

Алгоритм и тестовые примеры

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

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

Задача 3D реконструкция на основе уравнений (1) была протестирована на следующем примере (рис.5).

Были подготовлены изображения разноцветного 3D-треугольника с помощью специально разработанного для этой цели WinAPI OpenGL приложения (Программа photo3). Изображения создаются под заданными углами к вспомогательной СК относительно главной СК и сохраняются через клавишу Prtsc в буферной памяти, а затем обрабатываются в графическом редакторе PAINT и сохраняются как BMP файл (24-розрядный рисунок).

Для тестирования задачи была разработана программа (Программа angle3), которая загружает файлы 2-х изображений и по ним реконструирует 3D объект в виде облака точек.

                    

 

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

  • перевести координаты точки из пространства фото (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).

В случае, когда фотографии получены с обычного фотоаппарата, взаимное положение камер может быть определено на глаз. Ниже приводится результат (рис. 10 -12) использования программы 3D реконструкции детской игрушки по 2 фотографиям.  Фотографии были получены под разными углами поворота фотоаппарата.  Углы в этом случае определялись с погрешностью в пределах 5 градусов. Полученный результат свидетельствует о том, что даже при такой невысокой точности определения положения камеры был получен относительно неплохой результат.

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

Выводы

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


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

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

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

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

Программа photo3

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 gm = 30.0, bt = 0.0, al = 45.0;
 gm = 3.14*gm / 180;
 bt = 3.14*bt / 180;
 al = 3.14*al / 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

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 ** Engine::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();
 void SetTranslationMatrix_5();
void SetTranslationMatrix_6();
 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(); // hi=- 45 Z
 CurrentMatrix.MultiplyMatrices(Tr); 
 Tr.SetTranslationMatrix_5(); // fi=-30 X
 CurrentMatrix.MultiplyMatrices(Tr); 
 Tr.SetTranslationMatrix_6(); // psi=15 y
 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 gm = 30.0, bt = 25.0, al = 0.0; //вариант проходящий!!! только через направление взгляда
//double gm = 30.0, bt = 0.0, al = -45.0; //вариант проходящий, но dmin>10
 gm = 3.14*gm / 180;
 bt = 3.14*bt / 180;
 al = 3.14*al / 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(){ // hi=45 Z
 SetUnit();
 double hi = 15.0;
 hi = 3.14*hi / 180;
data[0][0] = cos(hi); // a // a c p 0 //
 data[0][1] = sin(hi); // 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(hi); // b 
 data[1][1] = cos(hi); // d
 //data[1][2] = -1.0; // q
 //data[1][3] = 0.1; // 0
//data[2][0] = 0.707; // h угол 45
 //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_5(){ // fi=30 X
 SetUnit();
 double fi = 25.0;
 fi = 3.14*fi / 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(fi); // d 
 data[1][2] = sin(fi); // q 
 //data[1][3] = 0.1; // 0 угол 30
//data[2][0] = 0.5; // h
 data[2][1] = -sin(fi); // f 
 data[2][2] = cos(fi); // 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(){ // psi=-15 y
 SetUnit();
 double psi = -30.0;
 psi = 3.14*psi / 180;
 data[0][0] = cos(psi); // a // a c p 0 //
 //data[0][1] = -0.707; // c // b d q 0 //
 data[0][2] = -sin(psi); // 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(psi); // h угол -15
 //data[2][1] = 1.0; // f
 data[2][2] = cos(psi); // 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;
}

Тестовые рисунки размещаются в папке Data.

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

3D-реконструкция по изображениям — Станислав Протасов

3D-cканирование с помощью фотокамеры — Виталий Окулов

123D Catch получение 3D модели из фотографий

Обучение по использованию 3д-сканера

Говорящие руки. Технология 3D-сканирования

Топ-10 приложений для 3D-принтеров

Faceworx + MeshMixer = 3d модель по фото

 

 

Авторы: Николай Свирневский и Антон Дядюк

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

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