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.
|
