Нам якимось чином потрібно реалізувати алогоритм зору для нашої камери.

У нас є камера, яка знаходиться у воксельній сітці і спрямована в якусь сторону, як нам зрозуміти і головне відобразити на екрані, що бачить камера? Перше шо приходить на думку - це від камери пустити промінь і шоб цей промінь міг перетнутись із якимось об'єктом і повідомити нас про це, щоб ми знали де саме відбулась колізія.

На щастя, такий алгоритм вже існує і він називається DDA(Digital Differential Analyzer) - це алгоритм трасування променя у дискретній сітці.
Лишилось тільки розібратись як він працює та імплементувати))

Як працює сам алгоритм (розглянемо 2d версію)

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

  • Кожний воксель має розмір 1 x 1, тобто займає всю клітинку.
  • Для усього проєкту буде використана правостороння система координат де Z+ йде на глядача.
    І хоча ці уточнення не виглядають важливими, ми далі побачимо, що саме такі налаштування допоможуть нам максимально спростити алгоритм.

ПОЧНІМО!

У променя є початкова і кінцева координата.

 vec2 rayStartPoint = vec2(5.6, 2.6);
 vec2 rayEndPoint = vec2(7.3, 1.3);

Задача: через які вокселі проходить промінь?

Очевидним є рішення використати цикл while (isVoxelInCell): розбити промінь на маленькі шматочки, наприклад крок 0,01 і постійно перевіряти чи не вперлись ми у воксель. Але через таку логіку буде дуже багато непотрібних перевірок, а цей алгоритм трасування променя ми будемо викликати для кожного пікселя на екрані і за його допомогою вираховувати нормалі, туман, прозорість, освітлення, тіні, колізію. чанки, LOD і т.д. Тобто він має бути досить швидким. Адже потрібно максимально зменшити кількість кроків.

Алгоритм немає очей, він не може подивитись на сітку і сказать, через які саме вокселі він проходить. Перше, що нам потрібно - це дізнатись у яку сторону дивиться промінь відносно свого початку.

Як людина, яка колись була пов'язана з армією, перше, що потрібно зробити по прибутті на нове місце - озирнутись на місцевості. Ми можемо розбити простір навколо початкової позиції променя на 4 сектори. Тобто секторів, у яких може знаходитись промінь всього 4.

360 / 4 = 90 - і це досить не погано! Ми із повної невизначенності можемо скоротити наш район пошуку, який буде, одним із 4 секторів

Маленька зноска - камера не буде передавати нам кінцеву позицію променя. вона буде передавати нормалізований вектор(вектор, який має довжину 1, його використовують для зазначення напрямку), довжину променя і початок променя. І сам нормалізований вектор буде мати вигляд типу (0.6 -0.8). це означає шо відносно його початку, яким є початкова точка (5.6, 2.6) він буде зміщений на 0.6 вправо і на -0.8 вниз.

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

Нехай початок A=(x1,y1)A = (x_1, y_1), кінець B=(x2,y2)B = (x_2, y_2).

1. Проекція вектору:

v⃗=B−A=(x2−x1,  y2−y1)\vec{v} = B - A = (x_2 - x_1,\; y_2 - y_1)

2. Довжина:

∥v⃗∥=(x2−x1)2+(y2−y1)2\lVert \vec{v} \rVert = \sqrt{(x_2 - x_1)^2 + (y_2 - y_1)^2}

3. Нормалізований вектор:

u^=v⃗∥v⃗∥=(x2−x1∥v⃗∥,  y2−y1∥v⃗∥)\hat{u} = \frac{\vec{v}}{\lVert \vec{v} \rVert} = \left( \frac{x_2 - x_1}{\lVert \vec{v} \rVert},\; \frac{y_2 - y_1}{\lVert \vec{v} \rVert} \right)

Приклад: A=(5.6,2.6)A = (5.6, 2.6), B=(7.3,1.3)B = (7.3, 1.3)

v⃗=(7.3−5.6,  1.3−2.6)=(1.7,  −1.3)∥v⃗∥=1.72+(−1.3)2=2.89+1.69=2.14u^=(1.72.14,  −1.32.14)=(0.795,  −0.607)\begin{aligned} \vec{v} &= (7.3 - 5.6,\; 1.3 - 2.6) = (1.7,\; -1.3) \\ \lVert \vec{v} \rVert &= \sqrt{1.7^2 + (-1.3)^2} = \sqrt{2.89 + 1.69} = 2.14 \\ \hat{u} &= \left( \tfrac{1.7}{2.14},\; \tfrac{-1.3}{2.14} \right) = (0.795,\; -0.607) \end{aligned} 
vec2 rayDirecion = vec2(0.795, -0.607); // відносно rayStartPoint

Це нам дає, шо вектор правіше на 0,795 і нижче на 0,607.

Це зрозуміло для нас, а як це порахувать у коді?

int getDirection(float rayCoord) { // компонента напрямку променя
    if (rayCoord > 0.0) return 1;  // праворуч / угору /
    if (rayCoord < 0.0) return -1; // ліворуч / вниз /
    return 0;                      // паралельно вісі
}

Ну шож ми звузили коло пошуку до 90 градусів)) Але як тепер точно дізнатися куди саме йде промінь?

На малюнку видно із якими двома сторонами вокселя контактує наш сектор, це права и нижня сторона (пам'ятаємо шо промінь обмежений сектором у 90 градусів) це ми вже вичислили.

Тепер ми можемо дізнатися, через яку сторону пройде наш промінь. Як? Якийх шлях коротший із двох, через ту сторону решітки і пройде промінь першим. На малюнку ми візуально бачимо, що координата х коротша ніж у і це означає, шо промінь проходить десь, через праву сторону клітинки. Але знову ж таки, це бачимо ми, компьютер не вміє бачити, йому потрібно якось це вирахувати. Для того шоб компютеру дізнатись, яка з координат менша, а яка більша, нам потрібно для початку порахувати довжини наших відрізків. Тут є одна хитрість, яку формулу використать і звучить вона так:

Потрібно брати координати тієї сторони, через яку може пройти промінь. Тобто у данному випадку для х ми використовуємо координату правої сторони бо промінь проходить крізь неї, а для y - нижньої.

const A = (5.6, 2.6); // x, y
const cellStart = (5, 2);
const cellLength = 1;

const AB = (cellStart.x + cellLength) - A.x;
AB = (5 + 1) - 5.6;
AB = 0.4;

const AC = A.y - cellStart.y;
AC = 2.6 - 2;
AC = 0.6;

0.4 < 0.6;
AB < AC; //

Давайте трохи підсумуємо:

  1. якщо х правіше точки старту то використовуємо формулу (cellStart.x + cellLength) - A.x;
  2. якщо х лівіше точки старту то використовуємо формулу A.x - cellStart.x;
  3. якщо y вище точки старту то використовуємо формулу (cellStart.y + cellLength) - A.y;
  4. якщо y нижче точки старту то використовуємо формулу A.у - cellStart.y;

Я думаю ви побачили шо формули лише дві - "рахуємо від початку клітинки або від кінця."

і ось сам код:

float getFirstBounds(int direction, float pointInCellCoord, float cellCoord) {
    if (direction > 0) return (cellCoord + 1.0) - pointInCellCoord;
    if (direction < 0) return pointInCellCoord - cellCoord;
    return 0.0;
}

Тепер, коли ми знаємо, через яку сторону клітинки пройде промінь, нам потрібно вирахувати час, за який промінь пройде від початку до вертикальної або горизонтальної сторони клітинки. Але тут ми зробимо відступ, щоб усе пояснити.

Що насправді таке нормалізований вектор?

Нормалізований вектор ще називається одиничним(тому що, його довжина дорівнює 1); Це означає, що його координати вказують лише на напрямок руху. До прикладу: вектор AB (0.894, 0.447); - це означає, що за одиницю дистанції його координата зміститься на 0.894 по х і на 0.447 по y. Маючи дві точки ми можемо побудувати пряму. Одна точка - це старт (0, 0), а інша - координати вектора (0.894б 0,447)

насправді старт це не завжди 0, 0 друга координата локальна по відношенню до старту. Тобто яким би не був старт ми завжди на 0.894 правіше і на 0.447 вище, а старт може бути любий, і в принципі його можна вважати за 0.

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

Для нормалізованого вектора довжина дорівнює 1. Отже, промінь за одиницю часу проходить рівно одну одиницю відстані.

Ця швидкість не змінюється ні за величиною, ні за напрямком: промінь не розганяється, не гальмує і не повертає. Такий рух у фізиці називають рівномірним прямолінійним.

Рівномірний прямолінійний рух — це рух уздовж прямої, під час якого тіло за будь-які рівні проміжки часу проходить однакові відстані. Швидкість такого руху стала.

Тепер, у нас є формула рівномірного лінійного руху, яка виглядає наступним чином: s⃗=v⃗⋅t\vec{s} = \vec{v} \cdot t. де:

s⃗\vec{s} - це вектор переміщення (який йде від початкової точки до кінцевої. його довжину рахують як «кінець мінус початок»:)
v⃗\vec{v} - це вектор швидкості, тобто наш нормалізований вектор.
tt - це час.

Давайте зводить усе до купи.

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

Переміщення відбувається за якийсь час. Знаючи конкретний час ми можемо вирахувати де саме буде знаходитись вектор по якійсь вісі! Наша формула s⃗=v⃗⋅t\vec{s} = \vec{v} \cdot t. може дати нам відповідь на це, але вона має і вектори і прості скалярні величини, а ми не можемо разом працювать змішуючи їх. Нам потрібно викинути вектори із розрахунку. Як? давайте рахувать для окремої вісі! Кожен вектор несе у собі дві компоненти (х, y); і шоб його розкласти ми заберемо одну із них і підставимо у формулу.

Почнемо з х Як сказано у формулі що s⃗\vec{s} рахується як кінець мінус початок і о чудо ми це вже порахували) Можливо, ви памятаєте як ми шукали, через яку сторону клітинки буде проходити промінь х там був 0.4, а у, 0.6 і промінь проходить через праву сторону клітинки, так ось це і є наша довжина переміщення променя по х. А як ми з'ясували раніше наш вектор напрямку(нормалізований) - це і є швидкість. щоб підставити у формулу на потрібно звідти витягнуть х компоненту. Але! координата вектору напрямку може мати відємний знак, якщо х знаходиться лівіше точки старту, а так як ми шукаємо tt, а tx=∣sx∣∣vx∣t_x = \frac{|s_x|}{|v_x|}. і не можу бути відємним, то якшо у одному із операндів буде -, то і у результаті буде -, нам потрібно х компоненту вектора напрямку взяти по модулю, щоб позбутися знака і тоді в результаті маємо

const tx = AB / abs(vectorDirection.x); // abs - це функція взяття модуля від числа. 
tx = 0.4 / 0.795;
tx = 0.503;

Тепер рахуємо те саме для y:

const ty = AC / abs(vectorDirection.y); // abs - це функція взяття модуля від числа. 
ty = 0.6 / 0.607; // було від'ємне
ty = 0.988;

Тут потрібно розуміть, що час, який ми знайшли по у це час, за який промінь дійшов би до горизонтальної сторони клітинки, а нас цікавить довжина відрізку від початку А.у і до ?.

Насправді, порахувать це досить просто час за який промінь пройде по х до вертикальної сторони клітинки такий самий як і для у, щоб промінь дійшов до вертикальної сторони клітинки (зауважте до вертикальної не горизонтальної), тому що промінь йде по різним вісям із різною швидкістю але за однаковий час. і користуючись тієюж формулою знайдем довжину цього відрізка

const A? = vectorDirection.y * tx; // тут без модуля бо потрібна знати куди зміститься y - вверх чи вниз. 
A? = -0.607 * 0.503;
A? = -0.305; // це зміщення по у, а не фінальна координата, щоб її знайти потрібно до початкової A.y додати зміщення.

І самий неочікуваний твіст - нам не потрібно знати довжину відрізка по у)) Ми будем далі порівнювати час. Це я написав для допитливих 😊

І код самої програми: 

float getTMax(int direction, float axisBound, float rayDirection) {
    if (direction == 0) return RAY_NO_HIT_T;
    return axisBound / abs(rayDirection);
}

Прикольно, так? стільки пояснень і так мало коду. От так воно буває коли маленька функція має під собою таке велике математичне підгрунтя.

Йдемо далі!

Нарешті ми вийшли за межі першої клітинки! Як нам рахувать це діло далі? Перша клітинка була така муторна тому, що промінь починався "десь" у її середині, а усі наступні він буде проходити від початку до кінця. Якщо подивитись на картинку то буде видно шо довжина променя при перетині різна, але, але, але! ми будемо використовувати лише одну формулу для розрахунку!

Щоб цей трюк спрацював, підійдемо до задачі з хитрістю: будемо рахувати не вхід і вихід, а лише вихід із клітинки.

Якщо глянути візуально, відстані між виходами різні. І це правда: промінь проводить у різних клітинках різний час.

Але якщо, промінь розкласти на х і у ...

Розглянемо окремо лише перетини вертикальних ліній. Вони стоять рівно через одну клітинку, тож між двома сусідніми вертикальними лініями промінь зміщується по x рівно на 1. З горизонтальними лініями так само: між двома сусідніми промінь зміщується по y рівно на 1.

Рух по x не залежить від руху по y. Швидкість стала, відстань між вертикальними лініями завжди 1, отже і час між ними завжди однаковий:

tΔx=1∣vx∣,tΔy=1∣vy∣t_{\Delta x} = \frac{1}{|v_x|}, \qquad t_{\Delta y} = \frac{1}{|v_y|}

Тобто ми маємо дві незалежні послідовності: вертикальні перетини через кожні tΔxt_{\Delta x} і горизонтальні через кожні tΔyt_{\Delta y}. Виходи з клітинок — це обидві послідовності, складені разом за часом. Тому на кожному кроці ми порівнюємо, яка лінія буде наступною: вертикальна чи горизонтальна.

У цьому алгоритмі ми оперуємо не координатами точки перетину, а часом, за який промінь дійде до горизонтальної лінії у , або до вертикальної х, іншими словами нас ціквить не "де перетин", а у яку клітинку промінь піде далі. і як було з'ясовано раніше: той часовий відрізок, який коротший - через ту сторону промінь пройде першим. 

const tDeltaX = 1 / abs(vectorDirection.x);
const tDeltaX = 1 / 0.795;
const tDeltaX = 1.257;

const tDeltaY = 1 / abs(vectorDirection.x);
const tDeltaY = 1 / 0.607;
const tDeltaY = 1.647;

Хочеш зрозуміти промінь? Стань промінем.

Все про шо ми говорили протягом цієї статі це були підготовчі роботи, ми готували усі необхідні данні для того, щоб запустити алгоритм. Давайте трохи порахуємо наш алгоритм вручну:

Стартова позиція:

const xt = 0;
const xy = 0;

Рахуємо час, для променя, і дивимось, яку лінію він перетне першу.

const xt = 0.503;
const xy = 0.988;
xt < xy -> наступна клітинка праворуч.

Тепер, коли ми знаємо, у яку клітинку піде промінь - додаємо до часу, який був менший на попередньому кроці, дельту(тобто час, за який промінь пройде до наступної клітинки) і порівнюємо час.

const xt += tDeltaХ;
const xt = 0.503 + 1.257;
const xt = 1.76;
const yt = 0.988; // оскільки y був більший у попередньому порівнянні до нього ми дельту не додаємо.
yt < xt -> наступна клітинка знизу.

Тут додаємо до уt:

const yt += tDeltaY;
const yt = 0.988 + 1.647;
const yt = 2.635;
const xt = 1.76; // оскільки х був більший у попередньому порівнянні до нього ми дельту не додаємо.
xt < yt -> наступна клітинка праворуч.

Тут додаємо до xt:

const xt += tDeltaX;
const xt = 1.76 + 1.257;
const xt = 3.017;
const yt = 1.76; // оскільки y був більший у попередньому порівнянні до нього ми дельту не додаємо.
xy < xt -> наступна клітинка знизу.

СТОП! далі наша лінія не буде сходитись із тим шо я намалював. Знаєте чому? Не тому що алгоритм не правильний, просто у мене руки із дупи ростуть і я не правильно намалював нахил лінії відносно координат, які ми рахуємо)) 🫠

Сподіваюсь, шо ідея зрозуміла.

Майже на фініші!

По суті у нас лишилось декілька питаннь:

  1. Коли промінь зрозуміє, де його кінець.
  2. Як ми знаємо чи є воксель у клітинці, яку ми пересікаємо.
  3. І що робити, коли промінь іде вздовж вісі.

Коли промінь зрозуміє, де його кінець?
Тут усе більш менш просто. Коли ми знаходимось у циклі ми заводимо ше одну змінну t і внеї додємо менший час на кожному кроці і порівнюємо із бажаним значенням. Якшо час більший за бажаний - промінь скінчився, виходимо із циклу.

 	float t = 0.0;
	float rayLength = 10.0; // максимальна довжина променю десять одиниць довжини
    for (int i = 0; i < 256; i++) {  // 256 граничне число ітерацій, щоб не було замкненог циклу.
        if (xTMax < yTMax) {
            t = xTMax; // ось тут додаєм  наш час до змінної
            xTMax += xDeltaT;
        } else {
            t = yTMax; // або тут додаєм  наш час до змінної
            yTMax += yDeltaT;
        } 

        if (t > rayLength) break; // промінь скінчився - виходимо.
    }

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

Як ми знаємо чи є воксель у клітинці, яку ми пересікаємо? У цьому самому циклі на початку ми викликаємо функцію, яка в свою чергу проходить по масиву із вокселями і перевіряє чи по координатам є воксель, якшо так - то виходимо із циклу.

float getCellStartPoint(float coord) {
    return floor(coord); // округлюємо координату клітинки вниз до цілого значення.
}

int xDirection = getDirection(rayDirection.x); // getDirection ми вже розглядали раніше
int yDirection = getDirection(rayDirection.y);

ivec2 currentCell = ivec2(getCellStartPoint(rayStartPoint.x), getCellStartPoint(rayStartPoint.y));

    for (int i = 0; i < 256; i++) {  // 256 граничне число ітерацій, щоб не було замкненог циклу.
        if (isVoxelSolid(currentCell)) { 
            hitCell = currentCell; // цю змінну ми повертаємо із функції. через out.
            return true;
        }
		if (xTMax < yTMax) {
            currentCell.x += xDirection;
        } else {
            currentCell.y += yDirection;
        } 
    }

Сама функція isVoxelSolid(currentCell) у мене поки не реалізована, але як тільки буде, я напишу про неї статтю))

І що робити, коли промінь йде вздовж вісі. Або на 0 ділити не можна!

Можливо ви помітили у коді деяких функцій 

int getDirection(float rayCoord) { // компонента напрямку променя
    if (rayCoord > 0.0) return 1;  // праворуч / угору /
    if (rayCoord < 0.0) return -1; // ліворуч / вниз /
    return 0;                      // паралельно вісі <- ОЦЮ СТРОКУ
}

Шо це за return 0 - це особливий випадок коли промінь йде паралельно вісі у цьому випадку х = 0. А оскільки на 0 ділити неможна ми заміняємо 0 на ду-у-уже велике число 

const float RAY_NO_HIT_T = 1e30; // 1 * 10 у 30 степені = нонільйон, або 1 із 30 нулями.

І тоді x у данному прикладі, завжди програє порівняння і ми завжди вибираємо y, і рухаємось завжди донизу.

P.S.

Ну щож, наче розібрали усе, що хотіли. Коли я копирсався у принципах роботи цього алгоритму, я навіть трохи приревнував геній людини, яка його створила до власної тупості)) Дуже витончено і майстерно вони виконали цю задачу, сподіваюсь, що і я колись буду взмозі написати щось настільки ж просте і геніальне одночасно!

Увесь код алгоритму

// DDA Algorithm

// Система координат: X праворуч+, Y вгору+, Z на себе+.
// Клітинка c займає відрізок [c, c+1) по будь-якій вісі,
// її початок — кут із найменшими координатами (ліво-низ-зад).
const float RAY_NO_HIT_T = 1e30;

float getCellStartPoint(float coord) {
    return floor(coord); // округлюємо вниз, щоб дізнатися початок клітинки
}

int getDirection(float rayCoord) { // компонента напрямку променя
    if (rayCoord > 0.0) return 1;  // праворуч / вгору / на себе
    if (rayCoord < 0.0) return -1; // ліворуч / вниз / від себе
    return 0;                      // паралельно вісі
}

float getFirstBounds(int direction, float pointInCellCoord, float cellCoord) {
    if (direction > 0) return (cellCoord + 1.0) - pointInCellCoord;
    if (direction < 0) return pointInCellCoord - cellCoord;
    return 0.0;
}

// t до першої межі по цій осі
float getTMax(int direction, float axisBound, float rayDirection) {
    if (direction == 0) return RAY_NO_HIT_T;
    return axisBound / abs(rayDirection);
}

// скільки t потрібно, щоб пройти цілу клітинку по цій вісі
float getDeltaT(float rayDirection, int direction) {
    if (direction == 0) return RAY_NO_HIT_T;
    return float(direction) / rayDirection;
}

// повертає true, якщо влучили; hitCell — у яку клітинку влучили
bool traceRay(vec3 rayStartPoint, float rayLength, vec3 rayDirection, out ivec3 hitCell) {
    int xDirection = getDirection(rayDirection.x);
    int yDirection = getDirection(rayDirection.y);
    int zDirection = getDirection(rayDirection.z);

    float xBound = getFirstBounds(xDirection, rayStartPoint.x, getCellStartPoint(rayStartPoint.x));
    float yBound = getFirstBounds(yDirection, rayStartPoint.y, getCellStartPoint(rayStartPoint.y));
    float zBound = getFirstBounds(zDirection, rayStartPoint.z, getCellStartPoint(rayStartPoint.z));

    float xTMax = getTMax(xDirection, xBound, rayDirection.x);
    float yTMax = getTMax(yDirection, yBound, rayDirection.y);
    float zTMax = getTMax(zDirection, zBound, rayDirection.z);

    float xDeltaT = getDeltaT(rayDirection.x, xDirection);
    float yDeltaT = getDeltaT(rayDirection.y, yDirection);
    float zDeltaT = getDeltaT(rayDirection.z, zDirection);

    float t = 0.0;
    // клітинка, де знаходиться старт
    ivec3 currentCell = ivec3(getCellStartPoint(rayStartPoint.x), getCellStartPoint(rayStartPoint.y), getCellStartPoint(rayStartPoint.z));

    for (int i = 0; i < 256; i++) {  // 256 — просто число, щоб не було нескінченного циклу
        if (isVoxelSolid(currentCell)) {
            hitCell = currentCell;
            return true;
        }

        // крокуємо по вісі, де межа найближча
        if (xTMax < yTMax && xTMax < zTMax) {
            t = xTMax;
            xTMax += xDeltaT;
            currentCell.x += xDirection;
        } else if (yTMax < zTMax) {
            t = yTMax;
            yTMax += yDeltaT;
            currentCell.y += yDirection;
        } else {
            t = zTMax;
            zTMax += zDeltaT;
            currentCell.z += zDirection;
        }

        if (t > rayLength) break; // промінь закінчився
    }

    hitCell = currentCell;
    return false;
}

/*
    Дві умови для правильної роботи: rayDirection має бути нормалізований,
    а isVoxelSolid має повертати false для клітинок за межами сітки.
*/