В Лекция 4, при примерите как да се ползват новите Vector и Matrix класове, имаше слайд ("Нетривиален пример") как да изчислим ортонормиран базис около зададен вектор.

Кодът там ползва random вектор, за да създаде първия перпендикулярен на v вектор (a), но има защити, които проверяват дали fabs(dot(randomVector, v)) > 0.9. По време на лекцията обясних,
че това е нужно, защото не е добра идея да подаваме на векторното произведение почти еднакви вектори, тъй като полученият трети вектор ще страда от намалена прецизност. Този вид проблеми се обобщават с термина Numerical stability и се анализират при проектиране на алгоритми, които ползват floating point числа.

Една от причините в курса по 3D графика да ползваме Vector-и от double-и е за да не се налага да мислим особено много за този род проблеми, тъй като там навлизаме в доста дълбоки води.
Ако се цели да се направи най-бързият рейтрейсър възможен, всички матрици и вектори биха ползвали float-и, а алгоритмите ще са старателно написани по най-стабилният (числено) начин, а това не би било особено fun занимание ако целта ни не е да ставаме експерти по numerical stability.

Все пак, дори и за сметките да ползваме double, floating point проблемите си стоят, просто изскачат по-трудно. Ще дам два примера

10 ротации с 1/10 оборот

В branch "et04-fp-stability" съм качил следното парче код:

const double N = 10;
Vector a(0.9372438, -0.2811731, 0.2061936);
Vector b = a;
Matrix m = rotationAroundY(toRadians(360 / N));
printf("Iter b                        dot(a, b) a^b                       len(a^b) norm(a^b)                dot(a, nc) dot(b, nc)\n");
for (int i = 1; i <= int(N); i++) {
	b *= m; // rotate by 1/N turn around Y
	Vector c = a^b;
	Vector nc = c;
	nc.normalize();
	printf("%4d ", i);
	printVec(b);
	printf("%8.5f  ", dot(a, b));
	printVec(c);
	printf(" %.5f  ", c.length());
	printVec(nc);
	printf("%9.6f  %9.6f\n", dot(a, nc), dot(b, nc));
}
return 0;

Какво се случва там?

• a е единичен вектор, който е само леко отместен от +X;
• b в началото е копие на a;
• на всяка итерация на for-цикъла, b се завърта с 1/10 оборот около оста Y (т.е. с по 36° на всяка итерация), и се смятат:
  • скаларното произведение dot(a, b);
  • векторното произведение a^b;
  • резултатният вектор a^b се нормира (и се нарича nc);
  • изчисляват се скаларните произведения между nc и {a, b};
• докато се върти около Y, векторът b се движи по повърхността на един широк конус (с връх 0, 0, 0 и ъгъл надолу към -Y);
• на 10тата итерация, b вече е завъртян на 360° и би следвало да съвпада с a.

Чистата математика би предсказала:

• a^b да е ненулев вектор за всички итерации без 10тата;
• Резултатите да са недефинирани за 10тата итерация. Би следвало a^b да е нулев вектор, който няма как да бъде нормиран;
• dot(a, nc) и dot(b, nc) следва да са нули винаги (резултатният вектор трябва да е перпендикулярен спрямо двата входни);

Реално се случва следното:

На проблемната 10-та итерация наистина a^b изглежда да е нулев вектор, и дължината му е не по-голяма от 10^-5, обаче векторът дефакто не е нулев и може да бъде нормиран, като даже излиза в правилната полуполовина и дори е прилично перпендикулярен на a и b! Разликата до прав ъгъл е под 2°.

Примерът горе е показателен за следните неща:

• Всяка манипулация на точки в пространството (в случая ротация) леко ги деформира. Например: единичен вектор след 1000 ротации може да изисква ново нормиране;
• 10 ротации от по 1/10 би следвало да резултират в идентитет, но това не е така;
• Точки, които би трябвало да съвпадат, всъщност "почти съвпадат", и в известен смисъл това е по-опасно от пълното съвпадение. При идентични вектори, горния код би крашнал, заради опита да нормираме нулев вектор. Тук тази математически невалидна операция тихичко минава и даже тихичко дава почти вярни резултати, а това би могло да ни докара много трудни за хващане бъгове

Light visibility biasing

В кода на visible() функцията при изчисление на сенките за Lambert шейдъра приложихме една промяна, за да избегнем визуалните проблеми, които се виждаха в началото (виж видеото от лекцията в 1:48:08). Този подход се нарича най-общо се нарича biasing и решава следният проблем:

• повечето Geometry::intersect реализации не правят твърда разлика дали обектът се намира пред или зад нас, т.е. подобно на Sphere::intersect връщат една или повече пресечни точки във вид на параметър p, където пресечната точка ip се смята като ray.start + ray.dir * p. Съответно, ако има пресечни точки "наобратно" зад ray.start, те ще се появят като решения с p ≤ 0
• когато трасираме лъч от камерата към геометриите, самата камера е невидима и можем да предполагаме, че в близост до нея няма нежелани геометрии; т.е. околността на ray.start е "чиста"
• когато трасираме visibility лъчите от пресечните точки към точката на светлината, вече ray.start тръгва от някоя геометрия и околността е "нечиста"
  • т.е. е възможно да получим тези "фантомни" пресечни точки с p ≅ 0
• възможно е заради грешки от закръгленията, досущ като в предния пример, да ни се получи параметър p = много малко, но положително число. Тогава ще си мислим, че нещо засенчва пресечната точка ip.
• нещо повече, възможно е действително геометрията да се самозасенчва. Представете си, че шейдваме точка във вътрешната извивка на един тор. Лъчът в посока към светлината спокойно може да има 4 решения на параметричното уравнение ray.start + ray.dir * p:
  1. отрицателно от близката външна стена на тора, зад нас
  2. ≅ 0, за шейдваната точка
  3. положително, от отсрещната вътрешна стена
  4. още по-голямо, от отсрещната външна стена

  Нас реално ни интересува №3, но ако стратегията на Torus::intersect е да ни върне "най-малкото положително решение", то би могло да ни върне №2.

• в лекция 4 решихме този проблем, като не трасираме visibility лъчи от геометрията към лампата, а от лампата към геометрията, и освен това сравняваме разстоянията чрез епсилон сравнение. Това, на което разчитаме, без да го заявяваме явно, е че околността на точката на лампата ще е "чиста";
• Това е вярно само за Omni Light, обаче. Като реализираме Rectangular light, например, тя е видима и си е геометрия като всяка друга, т.е. трасирането на visibility rays и в двете посоки е потенциално проблематично

Истинското решение на този проблем изисква или по някакъв начин да можем да викаме трасиране на лъч с игнориране на някакъв списък от околности около пресечни точки (което понякога има смисъл, особено за триъгълни мрежи: там можем да игнорираме триъгълника, от който "тръгваме"), или да отместим началото и края на лъча, така че да са мъничко, но измеримо отместени от повърхностите на геометриите. Именно това се нарича biasing, и изглежда така:

• Искаме да проверим за visibility между intersectionInfo.ip и lightPos
• Отместваме intersectionInfo.ip по протежение на нормала на пресечната точка: tempA = intersectionInfo.ip + intersectionInfo.norm * 1e-6
• Отместваме lightPos по протежение на посоката на светлината (което ще е добре дефинирано за светлини като RectLight, както ще видите по-натам): tempB = lightPos + lightDir * 1e-6
• Трасираме лъча между tempA и tempB, и проверяваме за пресечни точки на разстояние по-малко от distance(intersectionInfo.ip, lightPos). Ако не намерим такива - значи шейдваната точка е осветена от тази светлина.