Перспективная проекция в 3D: как реализовать на C++ и Python

[f/aspect, 0, 0, 0]
[0, f, 0, 0]
[0, 0, (far+near)/(near-far), 2*far*near/(near-far)]
[0, 0, -1, 0]

#include <cmath>
#include <iostream>

// Структура для 4×4 матрицы
struct Mat4 {
float m[4][4];

Mat4() {
// Инициализация единичной матрицы
for (int i = 0; i < 4; i++) {
for (int j = 0; j < 4; j++) {
m[i][j] = (i == j) ? 1.0f : 0.0f;
}
}
}

void print() {
for (int i = 0; i < 4; i++) {
for (int j = 0; j < 4; j++) {
std::cout << m[i][j] << " ";
}
std::cout << std::endl;
}
}
};

Mat4 createPerspectiveMatrix(float fovY, float aspectRatio, float nearPlane, float farPlane) {
Mat4 result;

// Проверка валидности входных параметров
if (nearPlane <= 0.0f || farPlane <= 0.0f || nearPlane >= farPlane || 
fovY <= 0.0f || aspectRatio <= 0.0f) {
throw std::invalid_argument("Invalid parameters for perspective projection");
}

// Преобразуем FOV из градусов в радианы
float fovRadians = fovY * 3.14159265f / 180.0f;

// Вычисляем множитель масштаба для Y
float f = 1.0f / std::tan(fovRadians / 2.0f);

// Заполняем матрицу
result.m[0][0] = f / aspectRatio; // Масштаб по X
result.m[1][1] = f; // Масштаб по Y

// Z-координата указывает на глубину в пространстве отсечения
result.m[2][2] = (farPlane + nearPlane) / (nearPlane – farPlane);
result.m[2][3] = -1.0f; // Для выполнения деления перспективы

// W-компонент для правильного вычисления глубины
result.m[3][2] = (2.0f * nearPlane * farPlane) / (nearPlane – farPlane);
result.m[3][3] = 0.0f; // Обнуляем, т.к. матрица не должна быть единичной в этой позиции

return result;
}

int main() {
try {
// Создаем матрицу перспективной проекции с FOV=60°, соотношением сторон 16:9
// и плоскостями отсечения 0.1 и 1000.0
Mat4 perspMatrix = createPerspectiveMatrix(60.0f, 16.0f/9.0f, 0.1f, 1000.0f);

std::cout << "Perspective Projection Matrix:" << std::endl;
perspMatrix.print();

// Пример применения матрицы к 3D точке
float point3D[4] = {10.0f, 5.0f, -20.0f, 1.0f}; // (x, y, z, w)
float projectedPoint[4] = {0.0f, 0.0f, 0.0f, 0.0f};

// Умножаем матрицу на точку
for (int i = 0; i < 4; i++) {
for (int j = 0; j < 4; j++) {
projectedPoint[i] += perspMatrix.m[i][j] * point3D[j];
}
}

std::cout << "\nOriginal Point: (" << point3D[0] << ", " << point3D[1] 
<< ", " << point3D[2] << ", " << point3D[3] << ")" << std::endl;

std::cout << "Projected Point before division: (" << projectedPoint[0] << ", " 
<< projectedPoint[1] << ", " << projectedPoint[2] << ", " 
<< projectedPoint[3] << ")" << std::endl;

// Выполняем деление перспективы
if (projectedPoint[3] != 0.0f) {
float w_inv = 1.0f / projectedPoint[3];
projectedPoint[0] *= w_inv;
projectedPoint[1] *= w_inv;
projectedPoint[2] *= w_inv;
projectedPoint[3] = 1.0f;
}

std::cout << "Final NDC Point: (" << projectedPoint[0] << ", " 
<< projectedPoint[1] << ", " << projectedPoint[2] << ")" << std::endl;

} catch (const std::exception& e) {
std::cerr << "Error: " << e.what() << std::endl;
return 1;
}

return 0;
}

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def create_perspective_matrix(fov_y, aspect_ratio, near_plane, far_plane):
"""
Создает матрицу перспективной проекции 4x4

Параметры:
fov_y (float): Угол обзора в вертикальной плоскости (в градусах)
aspect_ratio (float): Соотношение сторон (ширина/высота)
near_plane (float): Расстояние до ближней плоскости отсечения
far_plane (float): Расстояние до дальней плоскости отсечения

Возвращает:
numpy.ndarray: Матрица перспективной проекции 4x4
"""
# Проверка входных параметров
if near_plane <= 0 or far_plane <= 0 or near_plane >= far_plane or fov_y <= 0 or aspect_ratio <= 0:
raise ValueError("Invalid parameters for perspective projection")

# Преобразуем FOV из градусов в радианы
fov_radians = np.radians(fov_y)

# Вычисляем множитель масштаба для Y
f = 1.0 / np.tan(fov_radians / 2.0)

# Создаем матрицу
perspective = np.zeros((4, 4), dtype=np.float32)
perspective[0, 0] = f / aspect_ratio
perspective[1, 1] = f
perspective[2, 2] = (far_plane + near_plane) / (near_plane – far_plane)
perspective[2, 3] = -1.0
perspective[3, 2] = (2.0 * near_plane * far_plane) / (near_plane – far_plane)

return perspective

def apply_perspective_projection(points, perspective_matrix):
"""
Применяет перспективную проекцию к массиву точек

Параметры:
points (numpy.ndarray): Массив 3D точек формы (N, 3) или (N, 4)
perspective_matrix (numpy.ndarray): Матрица перспективной проекции 4x4

Возвращает:
numpy.ndarray: Проецированные точки в NDC форме (N, 3)
"""
# Убедимся, что у нас есть 4D точки (добавим w=1 если необходимо)
if points.shape[1] == 3:
homogeneous_points = np.hstack([points, np.ones((points.shape[0], 1), dtype=points.dtype)])
else:
homogeneous_points = points.copy()

# Применяем матрицу проекции
projected_points = np.dot(homogeneous_points, perspective_matrix.T)

# Деление перспективы
# Избегаем деления на ноль
mask = projected_points[:, 3] != 0
projected_points[mask, :3] /= projected_points[mask, 3:4]

# Возвращаем только x, y, z координаты
return projected_points[:, :3]

def create_cube():
"""
Создает куб с центром в начале координат и длиной ребра 2

Возвращает:
tuple: (вершины, ребра) где вершины – это массив 3D точек,
а ребра – список пар индексов вершин, образующих ребра куба
"""
# 8 вершин куба
vertices = np.array([
[-1, -1, -1], # 0
[ 1, -1, -1], # 1
[ 1, 1, -1], # 2
[-1, 1, -1], # 3
[-1, -1, 1], # 4
[ 1, -1, 1], # 5
[ 1, 1, 1], # 6
[-1, 1, 1] # 7
], dtype=np.float32)

# 12 ребер куба (пары индексов вершин)
edges = [
(0, 1), (1, 2), (2, 3), (3, 0), # Нижняя грань
(4, 5), (5, 6), (6, 7), (7, 4), # Верхняя грань
(0, 4), (1, 5), (2, 6), (3, 7) # Соединения граней
]

return vertices, edges

def visualize_projection(original_points, projected_points, edges):
"""
Визуализирует исходные 3D точки и их проекцию

Параметры:
original_points (numpy.ndarray): Исходные 3D точки
projected_points (numpy.ndarray): Проецированные точки в NDC
edges (list): Список пар индексов точек, образующих ребра
"""
fig = plt.figure(figsize=(15, 7))

# 3D визуализация оригинальных точек
ax1 = fig.add_subplot(121, projection='3d')
ax1.set_title('Original 3D Object')
ax1.set_xlabel('X')
ax1.set_ylabel('Y')
ax1.set_zlabel('Z')

for start, end in edges:
ax1.plot([original_points[start, 0], original_points[end, 0]],
[original_points[start, 1], original_points[end, 1]],
[original_points[start, 2], original_points[end, 2]], 'r-')

# Установим равные масштабы по осям
max_range = np.max(original_points.max(axis=0) – original_points.min(axis=0))
mid_x = (original_points[:, 0].max() + original_points[:, 0].min()) / 2
mid_y = (original_points[:, 1].max() + original_points[:, 1].min()) / 2
mid_z = (original_points[:, 2].max() + original_points[:, 2].min()) / 2
ax1.set_xlim(mid_x – max_range/2, mid_x + max_range/2)
ax1.set_ylim(mid_y – max_range/2, mid_y + max_range/2)
ax1.set_zlim(mid_z – max_range/2, mid_z + max_range/2)

# 2D визуализация проекции
ax2 = fig.add_subplot(122)
ax2.set_title('Projected Points (Perspective Projection)')
ax2.set_xlabel('X')
ax2.set_ylabel('Y')
ax2.set_aspect('equal') # Одинаковый масштаб по осям

for start, end in edges:
ax2.plot([projected_points[start, 0], projected_points[end, 0]],
[projected_points[start, 1], projected_points[end, 1]], 'b-')

# Ограничим отображение NDC пространства от -1 до 1
ax2.set_xlim(-1.5, 1.5)
ax2.set_ylim(-1.5, 1.5)
ax2.grid(True)

plt.tight_layout()
plt.show()

def main():
# Создаем куб
vertices, edges = create_cube()

# Сдвигаем куб на 5 единиц по оси Z для лучшей визуализации
vertices[:, 2] -= 5.0

# Создаем матрицу перспективной проекции
# FOV = 60 градусов, соотношение сторон = 1.0, ближняя плоскость = 1, дальняя = 100
perspective = create_perspective_matrix(60.0, 1.0, 1.0, 100.0)

# Применяем перспективную проекцию
projected_vertices = apply_perspective_projection(vertices, perspective)

# Выводим информацию
print("Original vertices:")
print(vertices)
print("\nPerspective matrix:")
print(perspective)
print("\nProjected vertices (after perspective division):")
print(projected_vertices)

# Визуализируем результаты
visualize_projection(vertices, projected_vertices, edges)

if __name__ == "__main__":
main()

// Оптимизированная версия с использованием SIMD (SSE)
#include <immintrin.h> // Для SSE/AVX инструкций

// Предположим, что наши данные выровнены по 16-байтной границе
// для эффективного использования SSE
struct alignas(16) Vec4 {
float x, y, z, w;
};

// Оптимизированное умножение вектора на матрицу с использованием SSE
Vec4 transformPointSSE(const Mat4& m, const Vec4& p) {
Vec4 result;

// Загружаем строки матрицы в SSE-регистры
__m128 row0 = _mm_load_ps(&m.m[0][0]);
__m128 row1 = _mm_load_ps(&m.m[1][0]);
__m128 row2 = _mm_load_ps(&m.m[2][0]);
__m128 row3 = _mm_load_ps(&m.m[3][0]);

// Загружаем точку в SSE-регистр
__m128 point = _mm_set_ps(p.w, p.z, p.y, p.x);

// Вычисляем скалярные произведения
__m128 res0 = _mm_dp_ps(row0, point, 0xF1);
__m128 res1 = _mm_dp_ps(row1, point, 0xF2);
__m128 res2 = _mm_dp_ps(row2, point, 0xF4);
__m128 res3 = _mm_dp_ps(row3, point, 0xF8);

// Объединяем результаты
__m128 tmp0 = _mm_or_ps(res0, res1);
__m128 tmp1 = _mm_or_ps(res2, res3);
__m128 resultVec = _mm_or_ps(tmp0, tmp1);

// Сохраняем результат
_mm_store_ps(&result.x, resultVec);

return result;
}

// Эффективная реализация для проецирования массива точек
void projectPointsBatch(const Mat4& projMatrix, 
Vec4* points, 
Vec4* results, 
size_t numPoints) {

// Оптимизация с использованием OpenMP для многопоточности
#pragma omp parallel for
for (size_t i = 0; i < numPoints; ++i) {
// Применяем матрицу к точке
results[i] = transformPointSSE(projMatrix, points[i]);

// Деление перспективы
if (results[i].w != 0.0f) {
float invW = 1.0f / results[i].w;
results[i].x *= invW;
results[i].y *= invW;
results[i].z *= invW;
results[i].w = 1.0f;
}
}
}

import numpy as np
from numba import njit, prange, float32

# Декоратор @njit компилирует функцию в машинный код
@njit(fastmath=True)
def apply_perspective_fast(points, matrix):
"""
Быстрое применение перспективной проекции с JIT-компиляцией через Numba
"""
result = np.empty((points.shape[0], 4), dtype=np.float32)

for i in range(points.shape[0]):
# Матричное умножение для каждой точки
for j in range(4):
val = 0.0
for k in range(4):
val += points[i, k] * matrix[k, j]
result[i, j] = val

# Деление перспективы
if result[i, 3] != 0:
w_inv = 1.0 / result[i, 3]
result[i, 0] *= w_inv
result[i, 1] *= w_inv
result[i, 2] *= w_inv
result[i, 3] = 1.0

return result[:, :3] # Возвращаем только x, y, z

# Параллельная версия с использованием prange
@njit(fastmath=True, parallel=True)
def apply_perspective_parallel(points, matrix):
"""
Параллельное применение перспективной проекции
"""
result = np.empty((points.shape[0], 4), dtype=np.float32)

# prange позволяет Numba распараллелить цикл
for i in prange(points.shape[0]):
# Матричное умножение для каждой точки
for j in range(4):
val = 0.0
for k in range(4):
val += points[i, k] * matrix[k, j]
result[i, j] = val

# Деление перспективы
if result[i, 3] != 0:
w_inv = 1.0 / result[i, 3]
result[i, 0] *= w_inv
result[i, 1] *= w_inv
result[i, 2] *= w_inv
result[i, 3] = 1.0

return result[:, :3] # Возвращаем только x, y, z

// Настройка перспективной камеры в Unity
Camera mainCamera = Camera.main;
mainCamera.fieldOfView = 60.0f; // FOV в градусах
mainCamera.aspect = 16.0f / 9.0f; // Соотношение сторон
mainCamera.nearClipPlane = 0.1f; // Ближняя плоскость
mainCamera.farClipPlane = 1000.0f; // Дальняя плоскость

import OpenGL.GL as gl
import OpenGL.GLU as glu
import numpy as np
import pygame

def setup_perspective(width, height, fov=45.0, near=0.1, far=100.0):
# Настройка области просмотра
gl.glViewport(0, 0, width, height)

# Переключение в режим проекции
gl.glMatrixMode(gl.GL_PROJECTION)
gl.glLoadIdentity()

# Создание перспективной проекции
aspect = width / height
glu.gluPerspective(fov, aspect, near, far)

# Переключение обратно в режим моделирования
gl.glMatrixMode(gl.GL_MODELVIEW)
gl.glLoadIdentity()

# Настройка позиции камеры
glu.gluLookAt(
5.0, 5.0, 5.0, # Позиция камеры
0.0, 0.0, 0.0, # Точка, на которую смотрит камера
0.0, 1.0, 0.0 # Направление "вверх"
)

import plotly.graph_objects as go
import numpy as np

# Создание данных для 3D-визуализации (например, функция z = sin(x) * cos(y))
x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
x_grid, y_grid = np.meshgrid(x, y)
z_grid = np.sin(x_grid) * np.cos(y_grid)

# Создание 3D-поверхности с перспективной проекцией
fig = go.Figure(data=[go.Surface(z=z_grid, x=x_grid, y=y_grid)])
fig.update_layout(
scene=dict(
xaxis=dict(range=[-5, 5]),
yaxis=dict(range=[-5, 5]),
zaxis=dict(range=[-1, 1]),
aspectratio=dict(x=1, y=1, z=0.5),
camera=dict(
eye=dict(x=1.5, y=1.5, z=1.5), # Положение камеры
up=dict(x=0, y=0, z=1) # Вектор "вверх"
)
)
)
fig.show()

import cv2
import numpy as np

# Калибровка камеры для получения внутренней матрицы
def calibrate_camera(calibration_images_path):
# Координаты точек калибровочной доски в 3D-пространстве
objp = np.zeros((6*9, 3), np.float32)
objp[:, :2] = np.mgrid[0:9, 0:6].T.reshape(-1, 2)

# Массивы для хранения точек объекта и изображения
objpoints = [] # 3D точки в реальном мире
imgpoints = [] # 2D точки в плоскости изображения

# Загружаем калибровочные изображения и находим углы шахматной доски
import glob
images = glob.glob(f'{calibration_images_path}/*.jpg')

for fname in images:
img = cv2.imread(fname)
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

# Находим углы шахматной доски
ret, corners = cv2.findChessboardCorners(gray, (9, 6), None)

if ret:
objpoints.append(objp)
imgpoints.append(corners)

# Калибруем камеру
ret, camera_matrix, dist_coeffs, rvecs, tvecs = cv2.calibrateCamera(
objpoints, imgpoints, gray.shape[::-1], None, None
)

return camera_matrix, dist_coeffs

# Использование матрицы камеры для проецирования 3D-точек на изображение
def project_points(points_3d, camera_matrix, rvec, tvec, dist_coeffs=None):
if dist_coeffs is None:
dist_coeffs = np.zeros(5)

# Проецируем 3D точки на 2D плоскость изображения
points_2d, _ = cv2.projectPoints(points_3d, rvec, tvec, camera_matrix, dist_coeffs)

return points_2d.reshape(-1, 2)

// Пример создания стереопроекционных матриц для VR (упрощенно)
void createStereoProjectionMatrices(float fov, float aspect, float nearPlane, float farPlane,
float ipd, Mat4& leftEyeProjection, Mat4& rightEyeProjection) {
// IPD – межзрачковое расстояние
float halfIPD = ipd / 2.0f;

// Создаем базовую матрицу перспективной проекции
Mat4 baseProjection = createPerspectiveMatrix(fov, aspect, nearPlane, farPlane);

// Смещаем проекцию для левого глаза
leftEyeProjection = baseProjection;
leftEyeProjection.m[0][2] = halfIPD / nearPlane; // Смещение по X

// Смещаем проекцию для правого глаза
rightEyeProjection = baseProjection;
rightEyeProjection.m[0][2] = -halfIPD / nearPlane; // Смещение в противоположную сторону
}