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検査の回数)はほとんど伸びません。セル数はレベルごとに4倍ずつ増えますが、殻をなすセルは2倍ずつしか増えないからです。

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

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

2Dで線分と軸並行ボックスなら候補軸は三つです。ボックスの $x$ 軸、ボックスの $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

一軸だけを使うと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になります。頂点を下側の端点でだけ数えるよう強制する半開区間の規則です。方向を三つ飛ばすコストも消えます。

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

表から読み取ることは二つです。第一に、軸並行の退化は判定方式そのものを変えて消すほうが、レイをもっと飛ばすより安く確実です。第二に、STLが閉じていなければ偶奇系はすべて崩れます。穴一つが内部全体を流体にしてしまいます。巻き数(winding number)や符号付き距離はこの場合でも答えを出しますが、定数コストがはるかに大きくなります。実務では偶奇で進みつつ、STLを先に閉じるほうを選びます。

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

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

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

$\mathbf{x}_f$ は流体ボクセルの中心、$\mathbf{e}_i$ は格子方向ベクトル、$\Delta x$ はボクセル寸法です。$q_i$ が方向ごとに異なるという点が肝心です。

壁の角度を0°から45°付近まで回すと、八方向の $q$ の棒が互いにずれ始めます。offsetを押していくと棒は連続的に滑るのに、ボクセルの等級ラベルは階段のように跳びます。オフセットが大きくなるとオレンジ色のBボクセルが灰色のGへ変わるところも観察ポイントです。

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

$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倍になるかも確認します。ここでずれるならソルバーを回す理由はありません。格子がすでに別の形状を抱えているからです。