GPU上でのvoxel構築手法(2)

以前の記事に引き続きvoxelの構築手法です。6-separatingの場合の具体的な実装を考えてみたいと思います。

以前の記事で説明した三角形との交差判定について

三角形とvoxelとの交差判定を説明した部分が間違ってました。参照した論文を読んだ際に、私が間違えて欠落させた条件があります。

  • 26-separatingの場合は先の条件に加えてAABBでの交差テストが必要です。
  • 6-separationgの場合は1Dと2Dのbounding boxのテストの際に平面の交差テストが必要です。

と記述されていました。

6-separatingの場合に追加する必要のある交差判定計算

参照した論文では省略されていて具体的には言及されていませんが、追加する必要のあると思われる処理を考えてみました。エッジフラグの判定だけで不十分と思われるケースを補完するには、各投影面でvoxelの中点を結んだ「ひし形」と、45度回転した軸でAABBの交差テストを行う必要があると考えられます。

voxelのサイズを1とすると、交差判定のひし形は45度回転すると0.5*sqrt(2.0)のAABBになります。そこでvoxelの中心から45度回転させて判定します。voxelの中心からの相対座標を計算すれば、三角形の頂点の回転は下記の計算で行えます。加えてAABBの判定は、回転した3頂点のX,Yそれぞれの最大最小値を算出して、1Dの重複テストをかければわかります。

const float2 v1, v2, v3;
float rf = 1.0 / sqrt(2.0);
float2 rv1, rv2, rv3;

rv1.x = v1.x * rf - v1.y * rf;
rv1.y = v1.x * rf + v1.y * rf;
.
.
.

float minX =  min(min(rv1.x, rv2.x), rv3.x);
float maxX  = max(max(rv1.x, rv2.x), rv3.x);
if ((maxX + 0.25*sqrt(2.0)) * (minX - 0.25*sqrt(2.0)) < 0.0)
  discard;
.
.
.

この交差判定は大小関係の判定のみで、実際の座標値は必要ないため以下の処理で問題ないと思われます。

rv1.x = v1.x - v1.y;
rv1.y = v1.x + v1.y;

float minX =  min(min(rv1.x, rv2.x), rv3.x);
float maxX  = max(max(rv1.x, rv2.x), rv3.x);
.
.
if ((maxX + 0.5) * (minX - 0.5) >= 0.0)
  discard;

6-separatingの場合のエッジフラグの判定


voxelの中心点を原点として考えます。この場合voxelの中心点におけるエッジフラグは、エッジ上の頂点(v)と法線(n)から、

dot(-v[0], n[0])

が正か負かで判定できます。
次にvoxelの中心からのオフセットの方法を考えます。たとえばX方向に+0.5オフセットした場所でのエッジフラグの判定は、

dot((0.5,0.0)-v[0], n[0])

なので、

dot((0.5, 0.0), n[0]) - dot(v[0], n[0])

と計算できます。

判定対象の4点のうち、エッジフラグの値が最大値となる点のオフセット方向は、エッジの法線方向に存在する点です。
6-separatingの場合はX,Yいずれかのオフセットのみなので、

max(abs(n[0].x), abs(n[0].y)) * 0.5 - dot(v[0], n[0]) 

で計算できると思います。

geometry shaderとpixel shaderにおける計算の仕分け

以上の処理をgeometry shaderとpixel shaderで分担して行います。

geometry shaderの処理
エッジフラグのテストの際に使用するオフセット値は、以下のようにXYZ各軸へ投影した際に使えるように値を計算しておきます。
各投影面におけるエッジの法線ベクトルは2D投影したエッジのベクトルを90度回転したものです。
三角形の面法線の向きをチェックしてどちらに回転させるかを決めます。エッジの法線ベクトルはgeometry shaderから渡すには数が多すぎです。オフセットの値のみ渡すことにします。

const float3 planeNormal;
const float3 e[3];

float3 eOfs[3];

// XY plane projection
float2 eNrm[3];
float isFront = sign(planeNormal.z);
eNrm[0] = float2(e[0].y, -e[0].x) * isFront;
eNrm[1] = float2(e[1].y, -e[1].x) * isFront;
eNrm[2] = float2(e[2].y, -e[2].x) * isFront;

const float2 an[3];
an[0] = abs(eNrm[0]);
an[1] = abs(eNrm[1]);
an[2] = abs(eNrm[2]);

eOfs[0].x = max(an[0].x, an[0].y)) * 0.5;
eOfs[0].y = max(an[1].x, an[1].y)) * 0.5;
eOfs[0].z = max(an[2].x, an[2].y)) * 0.5;
// YZ plane projection
// ZX plane projection

回転とAABBの判定部分もgeometry shaderで行える部分があります。
先ほどはvoxelの中心からのオフセット値で計算しましたが、voxelの中心値を後ほど回転して判定すればよいので、下記のように投影して回転した最大最小値を計算しておきます。

float4  prjMinmax[3];

// XY plane projection
float3 rpX = float3(vPos[0].x - vPos[0].y, vPos[1].x - vPos[1].y, vPos[2].x - vPos[2].y);
float3 rpY = float3(vPos[0].x + vPos[0].y, vPos[1].x + vPos[1].y, vPos[2].x + vPos[2].y);
prjMinmax[0].x =  min(rpX);
prjMinmax[0].y  = max(rpX);
prjMinmax[0].z =  min(rpY);
prjMinmax[0].w  = max(rpY);

// YZ plane projection
// ZX plane projection

ここまでの計算を行うとgeometry shaderから pixel shaderに渡すパラメーターは

const float3 vPos[3];
const flaot3 eOfs[3];
const float4 prjMinMax[3];

と、なかなかの量になります。これだけ長い計算を、並列性が比較的低いgeomerty shaderで行うのが得策かどうかわかりません。
pixel shader側に引き渡すパラメータも多いです。実際には、これに加えてvoxelに書き込む情報も必要になります。

pixel shader側の処理
pixel shaderではまず回転とAABBのテストを行います。

const float3 voxelCenter;

// XY plane projection
flaot2 rVoxelCenter = float2(voxelCenter.x - voxelCenter.y, coxelCenter.x + voxelCenter.y);
if ((prjMinMax[0].y -(voxelCenter.x-0.5)) * (prjMinMax[0].x - (voxelCenter.x+0.5)) >= 0.0 ||
    (prjMinMax[0].w -(voxelCenter.y-0.5)) * (prjMinMax[0].z - (voxelCenter.y+0.5)) >= 0.0)
discard;

// ZY plane projection
// XZ plane projection

続いてエッジフラグのチェックを行います。

const float planeNormal;
const float3 e[3];

float3 oPos[0] = vPos[0] - voxelCenter;
float3 oPos[1] = vPos[1] - voxelCenter;
float3 oPos[2] = vPos[2] - voxelCenter;

// XY plane projection
float2 eNrm[3];
float isFront = sign(planeNormal.z);
eNrm[0] = float2(e[0].y, -e[0].x) * isFront;
eNrm[1] = float2(e[1].y, -e[1].x) * isFront;
eNrm[2] = float2(e[2].y, -e[2].x) * isFront;

if (eOfs[0].x - dot(oPos[0].xy, eNrm[0].xy) < 0 ||
    eOfs[0].y - dot(oPos[1].xy, eNrm[1].xy) < 0 ||
    eOfs[0].z - dot(oPos[2].xy, eNrm[2].xy) < 0)
discard;

// ZY plane projection
// XZ plane projection

6-separatingの交差判定計算は長い

実際に書いてみると、とても長いシェーダーになりそうです。下図は6-separatingの判定領域になります。(注:6-separatingの判定領域の3次元表現は出来ないようです。参考までに図は残しておきます。)

なぜ6-separatingなのか

6-separatingはvoxelにおける最小の空間分割表現であると思われます。つまりポリゴンをvoxelizeしたときに、充填されるvoxel数が最小になるはずです。voxelの数が後々の処理で重荷になる可能性があることを考えると、voxelizeされるvoxel数は必要かつ最小であることが望ましいと思われます。その一方で、26-separatingに比べると処理は複雑になります。voxelizeしたvoxelの使用用途によっては正確な6-separatingである必要はなく「6-separating以上」であれば十分なケースも考えられます。このように条件を緩和すればこのあたりの処理はもっと早くすることが可能だと思います。次は違うvoxelizeの方法を検討してみたいと思います。と同時に、そろそろ実装して数字を取ってみる必要がある気がしてきました。

コメントを残す

以下に詳細を記入するか、アイコンをクリックしてログインしてください。

WordPress.com ロゴ

WordPress.com アカウントを使ってコメントしています。 ログアウト / 変更 )

Twitter 画像

Twitter アカウントを使ってコメントしています。 ログアウト / 変更 )

Facebook の写真

Facebook アカウントを使ってコメントしています。 ログアウト / 変更 )

Google+ フォト

Google+ アカウントを使ってコメントしています。 ログアウト / 変更 )

%s と連携中