1,681個のうち27個が内と外を取り違えた — STLボクセル化の偶奇判定とリンク切断

レイの偶奇判定は頂点ひとつで裏返ります。方向を増やすのは応急処置で、本当の解は区間を半開区間で閉じることです。

入力はSTLファイル一つとボクセル寸法一つだけです#

格子ボルツマン法(LBM)のソルバーに形状を渡すとき、手元にあるのは二つだけです。三角形の表面リストであるSTLファイル一つと、ボクセル一辺の長さ一つ。ところが出力すべきものはずっと多くなります。ボクセルごとに流体・固体・境界のどれなのか、境界なら隣接セルとの結合線のうちどれが壁に切られたのか、切られたなら壁までの距離比がいくつなのか。

この記事はその変換を四つの段階に割ります。オクトリーで候補を絞り、分離軸定理で交差を判定し、レイの偶奇で内外を分け、結合線の上に交点を打ちます。各段階で実際に壊れる場所がどこかも一緒に見ていきます。最後の段階で出てくる $q$ 値ひとつが、境界条件の精度を決めます。

三角形すべてをボクセルごとに検査すると値が積で膨らみます

もっとも素朴なやり方は、ボクセル一つごとに三角形すべてとの交差を調べることです。ボクセルが $N^3$ 個、三角形が $M$ 個なら検査回数は $N^3 M$ になります。 $N = 256$ 、 $M = 200{,}000$ なら $3.4 \times 10^{12}$ 回です。一日では終わりません。

オクトリー(octree、空間を8分割で再帰分割する木)がこの積を断ち切ります。ルートノード一つから始めて、そのノードに掛かる三角形のリストを作ります。リストが空でなければ子ノードを8個作り、各子は親のリストだけを再検査します。リストが空になればそこで止めます。

止まったノードが重要です。その中には表面が一つもないという意味なので、そのノード全体が流体か、全体が固体かのどちらかです。判定は一度で済みます。検査コストが体積ではなく表面積に付く構造へ変わります。

オクトリーのレベルスライダーを3から7へ上げると、オレンジの境界殻と緑のSAT棒が別々に動きます。

octree level5 (32²)

depth 0tree 0brute 0B 0

Raise the octree level and watch two things move in opposite directions: the orange boundary shell gets thinner in space but larger in count, while the green bar barely grows. Only the boxes the surface actually touches are ever split, so the tree pays for the surface, not for the volume.

殻は薄くなって個数は増えますが、検査回数の棒はほとんど止まります。セルはレベルごとに4倍、殻のセルは2倍です。

分離軸定理は軸を三つ見るだけです

ノードと三角形が重なるかどうかの判定には分離軸定理(Separating Axis Theorem, SAT)を使います。凸な二つの物体が交わらないなら、両者を分け隔てる軸が必ず一つ存在します。裏返して言えば、候補となる軸をすべて調べて一つも分け隔てられなければ、二つの物体は重なっています。

2Dで線分と軸並行ボックスなら候補軸は三つです。ボックスの $x$ 軸、ボックスの $y$ 軸、そして線分の法線。法線軸の判定はこう書けます。

\left| \mathbf{n} \cdot (\mathbf{p}_0 - \mathbf{c}) \right| > h_x |n_x| + h_y |n_y|

$\mathbf{n}$ は線分の法線、 $\mathbf{p}_0$ は線分の端点の一つ、 $\mathbf{c}$ はボックス中心、 $h_x, h_y$ はボックスの半辺です。不等式が成り立てばその軸が両者を分け隔てたという意味で、検査は即座に終わります。

3Dの三角形とボックスなら候補軸は13個です。ボックス面の法線が3個、三角形面の法線が1個、そして両者の辺方向を外積した軸が9個。13回の内積比較で片が付くという点が、ここでSATを使う理由です。早期終了がよく効くので平均コストは13よりずっと低くなります。

交差したノードは境界(B)ボクセルになります。このとき、そのノードに掛かる三角形のアドレスも一緒に保存しておきます。後で結合線と表面の交点を打つときに、そのリストをもう一度使うからです。

内と外を分けるのは交差回数の偶奇です

残るのは表面が触れなかったノードたちです。これらは全部が流体か全部が固体ですが、どちらなのかは形状全体を見ないと分かりません。

古典的な答えはジョルダン曲線定理です。点から任意の方向へ半直線を飛ばし、表面と交差した回数を数えます。奇数なら内部、偶数なら外部。表面が閉じてさえいれば方向はどれでも構いません。実装も短く済みます。三角形一つにつきレイ-三角形交差が一つあれば足ります。

問題はレイが三角形の辺や頂点をちょうど通るときです。その地点は二つの三角形が共有するので、交差が2回数えられることがあります。偶奇が裏返ります。しかもこの状況は例外的ではありません。ボクセル中心は格子の上に規則的に並び、CADから出たSTLの頂点も格子座標にぴったり合っている場合が多いのです。二つの規則格子が出会うと、軸並行のレイは頂点を頻繁に貫きます。

原典の文書が「 $x, y, z$ の三方向のうち一方向でも奇数なら内部」と書いているのは、この危険への備えです。一方向が失敗しても残りが拾うという理屈です。実際にどれだけ拾うのか数えてみました。

Pythonで裏返ったボクセルを数えてみました#

2Dに落として、ひし形を一つ41×41の格子へ入れました。頂点が $(\pm 1, 0)$ 、 $(0, \pm 1)$ なので、格子の中心線 $y = 0$ と $x = 0$ の上にちょうど載ります。閉区間による判定(両端点を含む)と半開区間による判定(片端だけ含む)を、同じ格子の上で比べました。

# ひし形(頂点が格子の中心線上にちょうど載るように配置)
def diamond_poly(r=1.0):
    return [(r, 0.0), (0.0, r), (-r, 0.0), (0.0, -r)]
 
def edges_of(poly):
    return [(poly[i], poly[(i + 1) % len(poly)]) for i in range(len(poly))]
 
# よく使われる「閉区間」判定 — 頂点を二度数えてしまう
def naive_crossings(px, py, poly, axis):
    n = 0
    for (x1, y1), (x2, y2) in edges_of(poly):
        if axis == 'x':
            a, b, c1, c2 = y1, y2, x1, x2
            p, q = py, px
        else:
            a, b, c1, c2 = x1, x2, y1, y2
            p, q = px, py
        if a == b:
            continue
        if min(a, b) <= p <= max(a, b):          # 両端を含む -> 頂点が重複
            t = (p - a) / (b - a)
            if c1 + t * (c2 - c1) > q:
                n += 1
    return n
 
# 半開区間による判定 — 頂点をちょうど一度だけ数える
def halfopen_crossings(px, py, poly, axis):
    n = 0
    for (x1, y1), (x2, y2) in edges_of(poly):
        if axis == 'x':
            a, b, c1, c2 = y1, y2, x1, x2
            p, q = py, px
        else:
            a, b, c1, c2 = x1, x2, y1, y2
            p, q = px, py
        if (a > p) != (b > p):                    # [a, b) の半開区間
            t = (p - a) / (b - a)
            if c1 + t * (c2 - c1) > q:
                n += 1
    return n
 
def cell_centers(n, lo=-1.5, hi=1.5):
    h = (hi - lo) / n
    return [lo + (i + 0.5) * h for i in range(n)], h
 
def truth_inside(px, py, r=1.0):
    return abs(px) + abs(py) < r                  # ひし形の解析的な判定
 
def sweep_axes(n=41):
    xs, h = cell_centers(n)
    poly = diamond_poly()
    bad = {'x-only': 0, 'y-only': 0, 'x-or-y': 0, 'half-open': 0}
    for py in xs:
        for px in xs:
            ref = truth_inside(px, py)
            ox = naive_crossings(px, py, poly, 'x') % 2 == 1
            oy = naive_crossings(px, py, poly, 'y') % 2 == 1
            hx = halfopen_crossings(px, py, poly, 'x') % 2 == 1
            bad['x-only'] += (ox != ref)
            bad['y-only'] += (oy != ref)
            bad['x-or-y'] += ((ox or oy) != ref)
            bad['half-open'] += (hx != ref)
    return bad, len(xs) ** 2, h
 
bad, total, h = sweep_axes(41)
print(f"grid 41x41, voxel size h = {h:.5f}, cells tested = {total}")
for k, v in bad.items():
    print(f"  {k:10s} misclassified {v:4d}  ({100*v/total:.2f}%)")
 
poly = diamond_poly()
for (px, py, tag) in [(0.0, 0.0, 'center'), (0.0, 0.9146, 'just above'), (-1.2, 0.0, 'outside left')]:
    cx = naive_crossings(px, py, poly, 'x')
    cy = naive_crossings(px, py, poly, 'y')
    hx = halfopen_crossings(px, py, poly, 'x')
    print(f"{tag:12s} ({px:+.4f},{py:+.4f})  closed x={cx} y={cy} | half-open x={hx} | truth={'IN' if truth_inside(px,py) else 'OUT'}")

grid 41x41, voxel size h = 0.07317, cells tested = 1681
  x-only     misclassified   27  (1.61%)
  y-only     misclassified   27  (1.61%)
  x-or-y     misclassified    1  (0.06%)
  half-open  misclassified    0  (0.00%)
center       (+0.0000,+0.0000)  closed x=2 y=2 | half-open x=1 | truth=IN
just above   (+0.0000,+0.9146)  closed x=1 y=2 | half-open x=1 | truth=IN
outside left (-1.2000,+0.0000)  closed x=4 y=0 | half-open x=2 | truth=OUT

wall angle63°offset0.40

F 0FB 0B 0G 0

Drag the offset and watch the q bars slide continuously while the class labels jump in steps. Turn the wall to a diagonal and the eight q values stop agreeing with each other — that spread is exactly what a halfway bounce-back throws away. Push the wall far enough and orange B voxels turn grey: no fluid link left, nothing to stream into.

壁角度を0°から45°へ回すと、八本の $q$ 棒が互いにずれます。offsetを押すと棒は滑り、等級ラベルは階段状に跳びます。オフセットが大きくなるとオレンジのBが灰色のGへ変わります。

一軸だけを使うと1,681個のうち27個が裏返ります。すべて $y = 0$ の行にある内部ボクセルです。レイが頂点 $(1, 0)$ を通ったせいで交差が1ではなく2と数えられ、偶奇が裏返って「外部」になりました。

二軸をORで束ねると誤分類は27個から1個へ落ちます。原典の規則は実際に働いています。残った1個は原点 $(0,0)$ です。 $x$ のレイも $y$ のレイもそれぞれ頂点を貫いたので、両方とも偶数を返しました。方向を増やす方式の限界がここで露わになります。3Dなら $z$ 軸がこの点を拾ってくれますが、三軸が同時に失敗する形状を作るのも難しくありません。

最後の行が本当の解です。区間を min <= p <= max から (a > p) != (b > p) へ変えると誤分類が0になります。頂点を下側の端点でだけ数えるよう強制する半開区間の規則です。方向を三つ飛ばすコストも消えます。

一軸の誤分類27個は、1,681セルの1.61%です。格子は41×41、区間は $[-1.5, 1.5]$ 、ボクセル寸法は $h = 0.07317$ です。裏返ったセルはすべて $y = 0$ の行の内部です。レイが頂点 $(1, 0)$ を通り、交差が1ではなく2になり、偶奇が外部へひっくり返りました。二軸のORは27を1へ落とします。残る1個、0.06%は原点 $(0,0)$ です。閉区間では $x=2$ 、 $y=2$ で両方とも偶数です。半開区間は $x=1$ になり、奇数へ戻ります。すぐ上の $(0, 0.9146)$ は閉区間の $x=1$ なので一軸はすでに正しく、 $y=2$ だけが偶数です。ORがここを拾います。左の外側 $(-1.2, 0)$ は真値が外部です。閉区間の $x=4$ 、 $y=0$ なので判定は合いますが、交差は2ではなく4に膨らんでいます。半開区間は $x=2$ に減らしても偶数のままです。偶奇だけは通った点が、交差の個数はすでに二倍です。

補間式へ $q = 0.5$ を入れると $(1-2q)$ は0です。後ろの項が消えて $2q/(1+q) = 1$ になります。half-wayが式の中から戻ってきます。1.61%を0.00%へ消す前に付けた $q$ は、壁のない27セルに壁を立てた値です。D2Q9のリンクは8本、D3Q27は26本です。一セルが裏返ると、そこからぶら下がるリンクの $q$ がまとめて別の壁を指します。消したBの中心を交点にすると $q = 1$ です。定義は $0 \le q_i < 1$ なのに、削除直後のリンクだけがその端に触れます。方向を増やして残した1セルと、区間を閉じて消した0セルでは費用が違います。前者は軸ごとに $O(M)$ を積み、後者は比較が一行です。

四つの判定を同じ表に並べます

手法	コスト	非閉曲面	軸並行の退化	副次的な利得
レイの偶奇、閉区間	$O(M)$ / 点	即座に破綻	裏返る (1.61%)	なし
レイの偶奇、半開区間	$O(M)$ / 点	即座に破綻	なし (0.00%)	なし
多軸のOR投票	$3 \times O(M)$	即座に破綻	ほぼなし (0.06%)	なし
符号付き距離 / 巻き数	$O(M)$ / 点、定数が大きい	耐える	なし	壁までの距離が付いてくる

表が示す差は費用より、壊れ方です。軸並行の退化は判定区間を変えて消すほうが、レイを増やすより安いです。STLが閉じていなければ、偶奇系は穴一つで内部全体を流体にします。巻き数(winding number)や符号付き距離はそこでも答えを出しますが、定数コストがずっと大きいです。頂点を外すレイは、連続衝突で壁を通り抜ける場合と同じ場所です。実務では偶奇で進み、STLは先に閉じます。

リンクが切れた場所に $q$ が生まれます#

ここまで来るとボクセルごとにF(流体)・B(境界)・S(固体)が付きます。ところがLBMにはもう一段階必要です。LBMはボクセル中心の値だけを使うのではなく、隣へ結合線に沿って分布関数を流し込むからです。D2Q9なら8本、D3Q27なら26本。壁が切るのはボクセルではなく、この結合線です。

境界ボクセルの結合線ごとに表面との交点 $I$ を打ちます。交点が複数あればボクセル中心にもっとも近いものを使います。その地点までの距離比が $q$ です。

q_i = \frac{\left| \mathbf{x}_I - \mathbf{x}_f \right|}{\left| \mathbf{e}_i \right| \Delta x}, \qquad 0 \le q_i < 1

$\mathbf{x}_f$ は流体ボクセルの中心、 $\mathbf{e}_i$ は格子方向ベクトル、 $\Delta x$ はボクセル寸法です。 $q_i$ は方向ごとに別の数です。出力直後の角度ダイヤルが、その八つの数を棒で描きます。

$q$ を無視してすべて $0.5$ に置くのが標準のhalf-way bounce-backです。実装はもっとも短く済みますが、壁の位置が実際の表面ではなく格子にスナップされます。曲面では階段状の誤差が残り、収束次数が2次から1次へ落ちます。 $q$ を使う補間bounce-backは、壁から返ってくる値をこう作ります。

f_{\bar{i}}(\mathbf{x}_f, t + \Delta t) = \frac{1}{1 + q} \left[ 2 q \, f_i^{\star}(\mathbf{x}_f) + (1 - 2q) \, f_i^{\star}(\mathbf{x}_f - \mathbf{e}_i \Delta t) \right]

$f^\star$ は衝突直後の分布関数、 $\bar{i}$ は $i$ の反対方向です。 $q = 0.5$ を入れると第二項が消えてhalf-wayの規則へ戻ります。つまり $q$ はhalf-wayを特殊な場合として含む一般化です。この補間を使うには、前の段階がリンクごとに $q$ を保存しておく必要があります。だからSATの段階で三角形のアドレスを捨てずに持ってきたのです。

流体の隣接がない境界ボクセルをなぜ消すのか

原典の手続きには、目に留まりにくい規則がもう一つあります。BボクセルのうちFボクセルと結ばれた線が一本もないものは、そのBを消します。消した場所はゴースト(G)ボクセルになります。

理由は計算量ではなく定義です。bounce-backは流体から来た分布関数を送り返す操作です。入ってくる分布関数がなければ、送り返すものもありません。そんなボクセルに境界条件を掛けると、初期化されたままのゴミ値が毎ステップのストリーミングに乗って流れ出します。

消す順序も重要です。Bを先に消してしまうと、そのボクセルと結ばれていた隣接セルの立場ではリンクが切れたことになります。原典の文書はこの場合、削除したBの中心点をそのまま交点 $I$ とせよと書いています。 $q = 1$ のリンクになるわけです。この処理がないと、削除の直後に未定義のリンクが残ります。

最後に、Fボクセルのうち交点 $I$ を一つでも持つものはFBへ昇格します。実際の計算で境界処理のループが回る対象がこのFB集合です。Fは純粋なストリーミング、FBはストリーミング + 補間bounce-back、Bは値の供給者、Gはそもそも外れます。四つの等級がそれぞれ別のカーネルに乗ります。等級をあらかじめ分けておく理由は、壁ノードが失う分布関数の本数を数えるときも同じです。分岐を毎ステップ判定せずリストとして持っておくのが、LBMの基本戦略です。

ボクセル一つが自分の名前を得るまで

STL一つが格子になる過程を辿り直すと、判定は四回ありました。オクトリーのノードが表面に触れるか(SATの13軸)、触れなかったノードは内か外か(レイの偶奇)、結合線が切られるか(符号の変化)、切られたボクセルが流体とつながっているか(リンクの本数)。

もっとも静かに間違うのは二番目です。一番目と三番目は間違えると絵が目立って壊れますが、偶奇が裏返ったボクセルは形状の内部に一列に埋まっているので、コンターでは見つけにくいのです。流量が数パーセントずれたまま検証を通ってしまうこともあります。

だから新しい形状を入れて最初にやることは、FとSの個数を解析的な体積と突き合わせることです。ボクセル寸法を半分に縮めて $F$ の個数が8倍になるかも確認します。ここでずれるならソルバーを回す理由はありません。格子がすでに別の形状を抱えているからです。