1. 程式人生 > >《實時碰撞檢測演算法技術》讀書筆記(六):最近點計算(下)

《實時碰撞檢測演算法技術》讀書筆記(六):最近點計算(下)

點至3D矩形的最近點

實際上等同於計算OBB上的最近點,其中3D矩形可看做是z向為0OBB

struct Rect {
    Point c;
    Vector u[2];
    float e[2];
}

z軸為0並重寫函式ClosestPtPointOBB()

//Given point p, return point q on (or in) Rect r, closest to p
void ClosestPtPointRect(Point p, Rect r, Point &q)
{
    Vector d = p - r.c;
    //Start result at center of rect; make steps from there
    q = r.c;
    //For each rect axis...
    for(int i = 0; i < 2; i++) {
         //...project d onto that axis to get the distance
        //along the axis of d from the rect center
        float dist = Dot(d, r.u[i]);
        //if distance farther than the rect extents, clamp to the box
        if(dist > r.e[i]) dist = r.e[i];
        if(disr < -r.e[i]) dist = -r.e[i];
        //Step that distance along the axis to get world coordinate
        q += dist * r.u[i];
    }
}

點至三角形的最近點

給定三角形ABC以及一點P,且令ABC上距P的最近點為Q。一種解決方案是:如果P正交投影后位於三角形內部,則投影點即為最近點Q,如果P的正交投影位於三角形ABC的外部,則最近點位於三角形的某一邊上。此時則遍歷每一條邊(ABBCCA),計算並返回P的最近點。但明顯此種方法並不高效。一種較好的演算法是:考查P位於三角形的哪一個Voronoi特徵域中。一旦確定,只須將P正交投影至該特徵域中即可獲取最近點Q

P位於包含頂點AVoronoi域由穿越A的兩個平面的負半空間的交集構成,其法線分別為B-AC-A


確定P位於某一邊的Voronoi域,可以採用:計算PABC

上正交投影的質心座標。

若點P不屬於任一頂點或邊的Voronoi域,則Q必位於三角形ABC的內部,實際上為正交投影R


P點在ABC上正交投影的質心座標為三角形RABRBCRCA的有符號面積與三角形ABC有符號面積之比。令n為三角形ABC的法線,且對於某t值,有R = P - tn,則R = uA + vB + wC的質心座標(u,v,w)為:u = rbc/abc, v = rca/abc以及w = rab/abc。通過某些向量計算,可以達到一定程度上的簡化目的。如,表示式rab可以按如下方式進行簡化:

Vector n = Cross(b - a, c - a);

float rab = Dot(n, Cross(a - r, b - r));    //proportional to signed area of RAB

float rbc = Dot(n, Cross(b - r, c - r));    //proportional to signed area of RBC

float rca = Dot(n, Cross(c - r, a - r));    //proportional to signed area of RCA

float abc = rab + rbc + rca;           //proportional to signed of ABC

(上面程式碼中,n應該為垂直於三角形的向量(法線)。Cross代表叉乘,設α為ab夾角,則Cross(a,b)的數學意義是absinα*n。而三角形中,若已知兩邊及它們的夾角α,其面積則是0.5*absinα。即三角形ABC的面積=1/2*abs(AB×AC)


R的質心座標可以直接通過P獲得,而無須計算R

於是,首先判斷P點是否在三角形3個頂點域內,是則取該頂點為最近點,否則繼續判斷是否在邊域內,若在邊域內,取投影在邊上的點,否則位於三角形內部,取正交投影座標R

對下面部分程式碼的解釋:

點投影在直線上:P = A + s*AB, s = snom/(snom + sdenom)(snom + sdenom)AB長度,snom指投影到AB上的長度APsnom負表示PAAB角度大於90度)。

面與面之間夾角:float vc = Dot(n, Cross(a - p, b - p))vc負表示面之間夾角大於90度。

Point ClosestPtPointTriangle(Point p, Point a, Point b, Point c)
{
    Vector ab = b - a;
    Vector ac = c - a;
    Vector bc = c - b;
 
    //Compute parametric position s for projection P’ for P on AB.
    //P’ = A + s*AB, s = snom/(snom + sdenom)
    float snom = Dot(p - a, ab), sdenom = Dot(p - b, a - b);
 
    //Compute parametric position t for projection P’ of P on AC,
    //P’ = A + t*AC, t = tnom/(tnom + tdenom)
    float tnom = Dot(p - a, ac), tdenom = Dot(p - c, a - c);
 
    if(snom <= 0.0f && tnom <= 0.0f) return a;    //Vertex region early out
 
    //Compute parametric position u for projection P’ of P on BC,
    //P’ = B + u*BC, u = unom/(unom + udenom)
    float unom = Dot(p - b, bc), udenom = Dot(p - c, b - c);
 
    if(sdenom <= 0.0f && unom <= 0.0f) return b;      //Vertex region early out
    if(tdenom <= 0.0f && udenom <= 0.0f) return c;    //Vertex region early out
 
    //P is outside (or on) AB if the triple scalar product [N PA PB] <= 0
    Vector n = Cross(b - a, c - a);
    float vc = Dot(n, Cross(a - p, b - p));
    //If P outside AB and within feature region of AB,
    //return projection of P onto AB
    if(vc <= 0.0f && snom >= 0.0f && sdenom >= 0.0f)
        return a + snom / (snom + sdenom) * ab;
 
    //P is outside (or on) BC if the triple scalar product [N PB PC] <= 0
    float va = Dot(n, Cross(b - p, c - p));
    //If P outside BC and within feature region of BC,
    //return projection of P onto BC
    if(va <= 0.0f && unom >= 0.0f && udenom >= 0.0f)
        return b + unom / (unom + udenom) * bc;
 
    //P is outside (or on) CA if the triple scalar product [N PC PA] <= 0
    float vb = Dot(n, Cross(c - p, a - p));
    //If P outside CA and within feature region of CA,
    //return projection of P onto CA
    if(vb <= 0.0f && tnom >= 0.0f && tdenom >= 0.0f)
        return a + tnom / (tnom + tdenom) * ac;
 
    //P must project inside face region. Compute Q using barycentric coordinates
    float u = va / (va + vb + vc);
    float v = vb / (va + vb + vc);
    float w = 1.0f - u - v;    // = vc / (va + vb + vc)
    return u * a + v * b + w * c;
}

以上程式碼呼叫了4次叉積運算。與點積運算相比,叉積運算通常成本較高。因此,採用某種更加經濟的表示式是值得考慮的。根據拉格朗日恆等式:


可表示3個標量三重積:

Vector n = Cross(b - a, c - a);

float va = Dot(n, Cross(b - p, c - p));

float vb = Dot(n, Cross(c - p, a - p));

float vc = Dot(n, Cross(a - p, b - p));

並利用6個點積計算對以下內容進行表達:

float d1 = Dot(b - a, p - a);

float d2 = Dot(c - a, p - a);

float d3 = Dot(b - a, p - b);

float d4 = Dot(c - a, p - b);

float d5 = Dot(b - a, p - c);

float d6 = Dot(c - a, p - c);

最終有:

float va = d3 * d6 - d5 * d4;

float vb = d5 * d2 - d1 * d6;

float vc = d1 * d4 - d3 * d2;

實際上,d1~d6也可以用於計算snomsdenomtnomtdenomunom以及udenom

float snom = d1;

float sdenom = -d3;

float tnom = d2;

float tdenom = -d6;

float unom = d4 - d3;

float udenom = d5 - d6;

其中向量n不再需要。

點到四面體的最近點

給定點P,當前問題為:計算四面體ABCD上(或之內)距離點P的最近點Q。一種較為直接的演算法是:針對四面體中的每一個面,呼叫函式ClosestPointTriangle()並計算最近點。在多個計算結果中,返回距點P最近的點Q。獨立於距離測試,還需要考查點P是否位於四面體中的面內。若P位於四面體中某一面內,則其即為最近點。


//Test if point and d lie on opposite side of plane through abc
int PointOutsideOfPlane(Point p, Point a, Point b, Point c, Point d)
{
    float signp = Dot(p - a, Cross(b - a, c - a));    //[AP AB AC]
    float signd = Dot(d - a, Cross(b - a, c - a));    //[AD AB AC]
    //Points on opposite sides if expression signs are opposite
    return signp * signd < 0.0f;
}

上訴演算法執行良好且易於實現。然而該演算法仍有改進餘地。首先,假定頂點PVoronoi特徵域存在,則QP在該正交域上的正交投影。針對四面體,需要考查總計14個特徵域:4個頂點域、6個邊域以及4個面域。若P不屬於任何一個特徵域,則必定位於四面體內部。

與前述測試方法相比,該演算法似乎無明顯的效能改觀。然而,Voronoi域間卻可以共享多數計算過程而無須重複計算。

點到凸多面體的最近點

一種易於實現的、複雜度為O(n)的演算法是:計算多面體上每個面中距P的最近點,並返回一個最近結果。同時,還要兼顧測試點P是否位於H各面的背面,即P位於H內部這一情況。為了加速測試過程,面距離測試計算規定:點P位於各面的正面。

對於較大的多面體,一種較為快速的預計算方法為:以層次結構的方式構建多面體的各部分。利用這種預構造的層級方案將使最近點計算的時間複雜度達到log級別。如Dobkin-Kirkpatrick層次結構演算法。

此外還有一些高效計算多面體最近點的其他一些演算法如GJK演算法等。

兩直線間的最近點

對於2D空間,兩直線不平行則相交,相交時最近點為交點,平行時的最近點為線上任意點。而在3D空間中則未必。故在3D空間中一個較好方案是:若兩直線之間彼此並未達到足夠接近,則設為不相交(未必平行)。基於此,相交測試轉換為兩直線間之間的距離,並檢測該距離是否小於給定的閥值。

給定L1由頂點Q1P1連成,L2由頂點Q2P2連成。


針對st的幾組計算結果,L1(s)lL2(t)分別對應直線上的最近點,且v(s, t) = L1(s) - L2(t)。當v具有最小長度值時,兩頂點則為最近點。其實即為計算st滿足以下垂直約束條件:


在此省略推導給出結果:

設a = d1·d1, b = d1·d2, c = d1·r, e = d2·d2, f = d2·r, d = ae - b^2

s = (bf - ce) / d     t = (af - bc) / d

注意:當d=0時兩直線平行且需要另行處理。


兩線段上的最近點


不難看出,在某些情況下,若直線間某一最近點位於相關線段的外部延長線上,該點是可以擷取至相應線段的最近端點處。

若直線間的最近點皆位於各自線段的外部延長線上,則上述擷取操作需要重複計算兩次。

給定直線上L2一點S2(t) = P2 + td2,L1上最近點L1(s)為:

s = (S2(t) - P1)·d1/d1·d1 = (P2 + td2 - P1)·d1/d1·d1

類似地,給定直線S1上一點S1(s) = P1 + sd1,直線L2上最近點L2(t)為:

t = (S1(s) - P2)·d2/d2·d2 = (P1 + sd1 - P2)·d2/d2·d2

可以使用高斯消元法求解線性方程組。

s = (bt - c)/a

 t = (bs + f)/e

其中,a = d1·d1, b = d1·r, e = d2·d2, f = d2·r

2D線段相交計算

測試兩線段ABCD是否相交,可以首先計算二者(延長線)的交點,且檢視該交點是否位於各自線段的包圍盒中。需注意的是,需要單獨處理兩線段平行這一特例。

計算兩線段(延長線)交點時,顯式定義直線L1的方程L1(t) = A + t(B - A),並隱式定義另一直線方程n·(X - C) = 0,其中垂直於CD。利用A+tB - A)帶入方程n·(X - C) = 0並求解t


則交點P = L1(t)可將t值代入上述顯式方程得到。另外,可通過檢測0<=t<=1從而確定P是否位於AB包圍盒內。

另一種2D線段相交測試方案為:兩線段相交僅當兩端點彼此位於線段兩側。可以通過判斷三角形ABDABC的環繞方向加以解決,而三角形的有符號面積可以確定其環繞方向(|a × b| = |a||b|sino,設點a(x1,y1), b(x2, y2)夾角為o|a||b|sino = |x1y2 - x2y1|,可得三角形有符號面積為0.5*x1y2 - x2y1)。 兩個非零向量ab平行,當且僅當|a × b|0)。對於某些應用,可以在執行線段相交測試前線處理線段的AABB盒相交測試。對於3D空間內成本高昂的線段相交測試,這一點尤為重要。


//Return 2 times the signed triangle area. The result is positive if
//abc is ccw, negative if abc is cw, zero if abc is degenerate.
float Signed2DTriArea(Point a, Point b, Point c) {
    return (a.x - c.x) * (b.y - c.y) - (a.y - c.y) * (b.x - c.x);
}
 

線段和三角形最近點


在此不討論平行情況。

最近點計算涉及下列物件:

線段PQ與三角形邊AB

線段PQ與三角形邊BC

線段PQ與三角形邊CA

線段端點P與三角形平面(P的投影位於三角形ABC內)

線段端點Q與三角形平面(Q的投影位於三角形ABC內)

上述計算返回具有唯一最小距離的最近點作為結果。

在某些情況下,可以適當調整測試的數量以減少計算量。如,若線段兩端點投影均位於三角形內部,則無需再執行線段-邊的最近點測試,因為最近點必然產生於線段的某一投影點中;若只是線段的一個端點投影位於三角形內部,則只須執行一次相應的線段-邊測試;當線段的投影均位於三角形外部,也只須執行一次相應的線段-邊測試。對於後兩種情況,可以通過計算線段端點所處的Voronoi域確定相關的線段-邊測試。

兩個三角形之間的最近點計算

兩個三角形T1T2之間的最近點計算一般可認為:相應最近點位於某一個三角形的邊上。因此,其結果轉化成可對三角形中6條邊之間的可能組合,執行線段-三角形的最近點計算。則包含最小(平方)距離的一組點即為最近點集合。

但線段-三角形之間距離測試成本較高。一種較好T1T2間最近點集測試方案為:一組最近點產生於兩個三角形的相應邊上,或產生於某一頂點與另一個三角形的內部點之間。問題轉化成:計算兩三角形間一對邊的最近距離,以及頂點與三角形之間的最近點(相應頂點的投影位於三角形內部)。綜上,需要執行6次頂點-三角形測試以及9次邊-邊測試。最終從多組最近點中,選擇具有最小距離的一組頂點作為最終三角形間的最近點。


(三角形平行以及重合的情況需要剔除)