63. fejezet - Oriented bounding box számolás


Ez most egy elég rövid cucc lesz, mert utálom a statisztikát, meg amúgyis kettes voltam belöle (amikor épp nem buktam meg). A feladat az, hogy egy modellre kiszámoljunk egy OBB-t, de nem feltétlenül kell, hogy minimális legyen. Minden ezzel foglalkozó paper ilyen mágikus statisztikai szavakat használ, mint kovariancia mátrix, ami valami olyasmi lehet mint a szórás, csak magasabb dimenzióban (vagy valami totál más, és büntiböl kirugják a honlapomat a szerverröl). Update: kirúgták :(

Hát akkó' variancia

Mindenféle magyarázat nélkül az említett mátrixhoz kell a ponthalmaz közepe (jobb helyeken átlag). Ezt még majdnem ki is tudjuk számolni:

CODE
for( quint32 i = 0; i < desc.NumVertices(0); ++i ) { mean._x += vdata[i].x; mean._y += vdata[i].y; mean._z += vdata[i].z; } mean /= (float)desc.NumVertices(0);

A kovarianca mátrix az ettöl a ponttól vett eltérés önmagával vett tenzorszorzata, pontosabban ezeknek az átlaga. Azaz valóban egy magasabb dimenziós szóráshoz hasonló. A tenzorok a mátrixok általánosításai, a tenzorszorzás pedig olyan müvelet ami pl. két vektorból egy mátrixot csinál (n és m-edrendü tenzorból n + m-edrendüt).

CODE
for( quint32 i = 0; i < desc.NumVertices(0); ++i ) { tmp._x = vdata[i].x - mean._x; tmp._y = vdata[i].y - mean._y; tmp._z = vdata[i].z - mean._z; covariance._m11 += tmp._x * tmp._x; covariance._m21 += tmp._y * tmp._x; covariance._m31 += tmp._z * tmp._x; covariance._m12 += tmp._x * tmp._y; covariance._m22 += tmp._y * tmp._y; covariance._m32 += tmp._z * tmp._y; covariance._m13 += tmp._x * tmp._z; covariance._m23 += tmp._y * tmp._z; covariance._m33 += tmp._z * tmp._z; } covariance *= 1.0f / (float)desc.NumVertices(0);

Valami okos ember kitalálta, hogy ennek a mátrixnak a sajátvektorai megadják az OBB koordinátatengelyeit. Azt persze már nem írta le, hogy ezt hogyan kell kiszámolni, úgyhogy most fellapozzuk a lineáris algebra tankönyvet.

Egy A mátrix karakterisztikus polinomjának nevezzük az alábbi polinomot: det(A - λI) Nyilván 3x3-as mátrixra ez egy harmadfokú polinom. Ennek a polinomnak a megoldásai a mátrix sajátértékei. Hogyan is kell kiszámolni egy 3x3-as mátrix determinánsát?

CODE
| a - λ b c | | d e - λ f | = cg(e - λ) + fh(a - λ) + bd(c - λ) - cdh - bfg - (a - λ)(e - λ)(c - λ) | g h i - λ |

Rögtön mindenki át is rendezi fejben, és természetesen mindenki álmából felébresztve keni-vágja a harmadfokú egyenlet megoldóképletét is. Akik aludtak volna általános suliban, azoknak itt a kód:

CODE
int solveCubic(float out[3], float coeff[4]) { const float& a = coeff[0]; const float& b = coeff[1]; const float& c = coeff[2]; const float& d = coeff[3]; float f = (3 * c / a - (b * b) / (a * a)) / 3; float g = (2 * (b * b * b) / (a * a * a) - 9 * b * c / (a * a) + 27 * d / a) / 27; float h = (g * g) / 4 + (f * f * f) / 27; if( f == 0 && g == 0 && h == 0 ) { out[0] = out[1] = out[2] = -fCbrt(d / a); return 3; } else if( h <= 0 ) { float i = powf((g * g) / 4 - h, 0.5f); float j = fCbrt(i); float k = acosf(-(g / (2 * i))); float m = cosf(k / 3); float n = 1.732050f * sinf(k / 3); float p = -(b / (3 * a)); out[0] = 2 * j * m + p; out[1] = -j * (m + n) + p; out[2] = -j * (m - n) + p; return 3; } else if( h > 0 ) { float r = -(g / 2) + powf(h, 0.5f); float s = fCbrt(r); float t = -(g / 2) - powf(h, 0.5); float u = fCbrt(t); float p = -(b / (3 * a)); out[0] = (s + u) + p; return 1; } return 0; }

Köbgyököt már csak tudtok írni (javaslom a powf használatát némi meggondolással). Most egy érdekes dolog jön, a sajátvektorok. Azok olyan állatok, amikre Av = λv ahol λ nyilván egy sajátérték. Szemfüles matematikusaok észrevették, hogy ezekböl végtelen sok van, hiszen ha v egy sajátvektor, akkor αv is az. Nyilván nekünk ennek az ezek által generált altérnek (sajátaltér) egy bázisa kellene.

Gondoljuk meg, mi is a helyzet: det(A - λI) = 0 minden λ sajátértékre, azaz ennek a lambdás mátrixnak van legalább kettö lineárisan összefüggö sora. Vadul kezdjünk el Gauss-eliminálni (a lambdát most nem írom ki):

CODE
a b c d e f g h i a b c 0 e - db/a f - dc/a 0 h - gb/a i - gc/a a b c 0 e - db/a f - dc/a 0 0 0

A fentiekböl nyilvánvaló, hogy az utolsó sor csupa nulla. Ha esetleg a mátrix nem ilyen lenne, akkor a sorokat tetszölegesen lehet cserélgetni és ilyen lesz. Mivel az utolsó sorban minden nulla, ez a paraméter tetszöleges, legyen mondjuk K. Ekkor a (lambdához tartozó) sajátvektorok az alábbi alakúak:

CODE
[ (-c - b * (dc/a - f) / (e - db/a)) / a ] K * [ (dc/a - f) / (e - db/a) ] [ 1 ]

Höfö ellenörizni, mert én tuti nem fogom. Az implementációban egy rakat if-el találtam ki, hogy hol nulla a mátrix.

K-t akár el is lehet felejteni, mind a három sajátértékre kiszámolva ezt a vektort (és mondjuk normalizálva öket) megvan a sajátaltér bázisa. Innen már elég könnyü a dolog, ugyanis:

CODE
qVector3 end[3]; end[0] = mean + vectors[0]; end[1] = mean + vectors[1]; end[2] = mean + vectors[2]; qVector3 tmin(float_max, float_max, float_max); qVector3 tmax(-float_max, -float_max, -float_max); float t; for( quint32 i = 0; i < desc.NumVertices(0); ++i ) { tmp._x = vdata[i].x; tmp._y = vdata[i].y; tmp._z = vdata[i].z; pointDistanceToLine(t, tmp, mean, end[0]); ASGN_IF(tmin._x, >, t); ASGN_IF(tmax._x, <, t); pointDistanceToLine(t, tmp, mean, end[1]); ASGN_IF(tmin._y, >, t); ASGN_IF(tmax._y, <, t); pointDistanceToLine(t, tmp, mean, end[2]); ASGN_IF(tmin._z, >, t); ASGN_IF(tmax._z, <, t); } bb.Vertices[0] = mean + tmin._x * vectors[0] + tmin._y * vectors[1] + tmin._z * vectors[2]; bb.Vertices[1] = mean + tmin._x * vectors[0] + tmin._y * vectors[1] + tmax._z * vectors[2]; bb.Vertices[2] = mean + tmin._x * vectors[0] + tmax._y * vectors[1] + tmin._z * vectors[2]; bb.Vertices[3] = mean + tmin._x * vectors[0] + tmax._y * vectors[1] + tmax._z * vectors[2]; bb.Vertices[4] = mean + tmax._x * vectors[0] + tmin._y * vectors[1] + tmin._z * vectors[2]; bb.Vertices[5] = mean + tmax._x * vectors[0] + tmin._y * vectors[1] + tmax._z * vectors[2]; bb.Vertices[6] = mean + tmax._x * vectors[0] + tmax._y * vectors[1] + tmin._z * vectors[2]; bb.Vertices[7] = mean + tmax._x * vectors[0] + tmax._y * vectors[1] + tmax._z * vectors[2];

Mi is történik itt? Vegyünk egy vertexet és vetitsük rá a tengelyekre (nyilván külön mindegyikre) és nézzük meg, hogy így milyen messze van a középponttól. A legnagyobb ilyen érték lesz az OBB erre a tengelyre vonatkozó mérete. Segítségképpen:

CODE
#define ASGN_IF(a, op, b) if((a) op (b)) (a) = (b); // legközelebbi pont az egyenesen == a + out * (b - a) // visszatérö érték: pont-egyenes távolság float pointDistanceToLine(float& out, const qVector3& p, const qVector3& a, const qVector3& b); { qVector3 bma(b - a); qVector3 pma(p - a); float x = qVector3::Dot(bma, pma); float y = qVector3::Dot(bma, bma); if( y == 0 ) out = 0; out = x / y; pma = a + out * bma; return qVector3::Distance(p, pma); }

Ha valakinek kérdése van, az így járt.


Summarum

OBB-t számoltam és kódot nem adok mert lusta voltam tutorialt írni.


Höfö:

  • Implementáld!

back to homepage

Valid HTML 4.01 Transitional Valid CSS!