복셀 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배로 가는지도 확인한다. 여기서 어긋나면 솔버를 돌릴 이유가 없다. 격자가 이미 다른 형상을 담고 있기 때문이다.