円柱と光線の交差判定
レイトレーシングにおける円柱との交差判定を考える。
ウェブや本などでいくつかの記述例があるが、それは大抵、複雑であるか、軸がY軸もしくはZ軸に固定であるかのどちらかだろう。
以下ではすっきりと記述でき、if文だらけにならず、なおかつ軸を固定しない手法をあげる。
本手法は、判定時に光線を変形する。これにより非常にきれいに記述することが出来る。
変形は、ワールド座標から、円柱の軸をZ軸に一致させるよう回転し、X軸・Y軸を１／半径に、Z軸を１／円柱の軸の長さにスケーリングする。
回転マトリクスはワールド空間で変形後にX・Y・Z軸になる３つの軸を求めれば記述できる。
これは以下のように記述できる。
//円柱クラス
class cylinder{
public:
cylinder(const vector3& p0, const vector3& p1, float radius);
public:
bool test(test_info * info,const ray & r,float dist)const;//交差判定メンバ
private:
vector3 p0_;//端点
vector3 p1_;//端点
float radius_;//半径
matrix3 m_;//変形用マトリクス
};
//最小の情報を持つ軸を返す
static inline int get_axis(const vector3& v){
int plane;
float tmp[3];
for(int i = 0;i<3;i++){
tmp[i] = abs(v[i]);
}
plane=0;
if(tmp[plane]>tmp[1])plane=1;
if(tmp[plane]>tmp[2])plane=2;
return plane;
}
//円柱クラスのコンストラクタ
cylinder::cylinder(const vector3& p0, const vector3& p1, float radius):p0_(p0),p1_(p1),radius_(radius){
vector3 E(p1-p0);
//スケーリング係数
float iX,iY,iZ;
iZ = float(1)/E.length();
iX = iY = float(1)/radius;
//ワールド座標での軸
vector3 tx,ty,tz;
tz = E*iZ; //Z軸になる
tx = vector3(0,0,0);tx[get_axis(E)] = 1;//軸の最小・・・出来るだけ直交した軸からX軸・Y軸を作る。
tx = normalize(tx-dot(tz,tx)*tz);//X軸になる
ty = cross(tz,tx); //Y軸になる
//回転マトリクス
matrix3 R = matrix3(
tx[0],tx[1],tx[2],
ty[0],ty[1],ty[2],
tz[0],tz[1],tz[2]
);
//スケーリングマトリクス
matrix3 S = matrix3(
iX, 0, 0,
0,iY, 0,
0, 0,iZ
);
m_ = S*R;
}
以上で前準備が出来たので、交差判定ルーチンが書ける。
ようは、マトリクスを光線に適用してやれば、軸がZ軸に０〜１にあり、半径が１の円柱との交差判定に帰結できる。
以下交差判定ルーチン
//0からmaxまでにxが入ってるかどうか
static inline bool is_0X(float x,float max){
return (0<x&&x<max);
}
//２次元の内積
static inline float dot2(const vector3& a,const vector3& b){
return a[0]*b[0]+a[1]*b[1];
}
//交差判定ルーチン
bool cylinder::test(test_info * info,const ray & r,float dist)const{
static float EPSILON = 1e-10;//お好みで
vector3 xO = m_*(r.origin()-p0_); //変形した光線の原点
vector3 xD = m_*(r.direction()); //変形した光線の方向・・・大きさごと変形
//ここで(-1,-1,0)-(1,1,1)とのAABB交差判定ルーチンをいれてもいいかも、実質的にOBB交差判定ルーチンになる。
//２次方程式の解・・・(-b+-sqrt(b^2-a*c))/a
float D2 = dot2(xD,xD);
float O2 = dot2(xO,xO);
float OD = dot2(xO,xD);
float T = OD*OD-D2*(O2-1);
if(T<EPSILON)return false; //正面からの場合もfalseで返す。キャップが必要ならそれを考慮して書き換える。
float SQ = sqrt(T);
float iD2 = float(1)/D2;
float t; //距離
t=(-OD-SQ)*iD2;
if(is_0X(t,dist)&&is_0X(xO[2] + t*xD[2],1)){
//表面なら
//infoに対し、必要な情報（距離・位置・法線など）を記述
return true;
}else{
//キャップ（天面・底面）が必要ならここに記述
//裏面かどうか・・・裏面が必要なければ即座にfalseをかえす
t = (-OD+SQ)*iD2;
if(is_0X(t,dist)&&is_0X(xO[2] + t*xD[2],1)){
//裏面なら
//infoに対し、必要な情報（距離・位置・法線など）を記述
return true;
}
}
return false;
}また法線は以下のように計算できる。
vector3 P = r.origin()+t*r.direction();//位置
vector3 V = P-p0_;
vector3 D = p1_-p0_;
vector3 W = (dot(V,D)/dot(D,D))*D;
vector3 N = normalize(V-W); //法線ちなみに、もし上記の方法が高速に動作するならば"is_0X()"が高速に動作してるのかもしれない。0以下かどうかはその符号を調べることになるからだ。とはいってもそれはコンパイラの最適化しだい。間違っても0との比較を、『floatをintとしてとらえ最上位のビットを調べる』といったことをしてはいけない。間違いなく遅くなるから・・・(Jim Blinn先生に対する皮肉である。x87浮動小数点レジスタ->汎用レジスタの動作が必要になり、遅くなる。)
