Implicit Functions & Signed Distance Functions

Reference: Stanley Osher and Ronald Fedkiw, Level Set Methods and Dynamic Implicit Surfaces, Springer, 2003, Ch. 1-2

레벨셋(level set) 방법을 공부하기 전에 가장 먼저 짚고 넘어가야 하는 것은 “인터페이스(interface)를 어떻게 표현할 것인가”라는 질문입니다. 물이 튀는 표면, 두 영역을 나누는 경계선, 움직이는 곡선 같은 것들을 컴퓨터 안에서 어떻게 다뤄야 할지 정하는 문제인데, 이 선택이 이후 등장하는 모든 수치해석 기법의 난이도를 결정합니다. 이 챕터는 그 표현 방법을 크게 두 가지, 즉 explicit 표현과 implicit 표현으로 나누고, 그중에서도 왜 implicit function이, 그리고 더 나아가 signed distance function이 좋은 선택인지를 설명합니다.

Introduction: 인터페이스를 어떻게 표현할 것인가

어떤 영역을 두 부분으로 나누는 경계, 즉 인터페이스를 표현하는 가장 직관적인 방법은 그 경계를 이루는 점들을 직접 나열하는 것입니다. 예를 들어 1차원에서 실수 전체를 $(-\infty,-1)$, $(-1,1)$, $(1,\infty)$ 세 구간으로 나눈다면, 경계는 $\partial\Omega={-1,1}$이라는 두 점으로 명시적으로 적을 수 있습니다. 이런 방식을 explicit interface representation이라고 부릅니다. 2차원에서는 경계가 곡선이므로 점을 나열하는 대신 매개변수 $s$에 대한 벡터 함수 $\vec{x}(s)$로 곡선을 표현해야 하고, 3차원에서는 곡면을 삼각형 메쉬 등으로 나타내야 합니다.

문제는 차원이 올라갈수록 explicit 표현이 점점 다루기 어려워진다는 데 있습니다. 특히 3차원 곡면의 경우 단순히 점들의 위치뿐 아니라 그 점들이 어떻게 서로 연결되어 있는지(connectivity)까지 알아야 하는데, 이 연결 관계를 정확히 계산하는 것 자체가 까다로운 문제입니다. 더 심각한 것은 표면이 움직이면서 서로 합쳐지거나(merge) 갈라지는(pinch) 경우입니다. 수영장에서 물이 튀는 장면을 떠올려보면, 표면의 연결 관계가 매 순간 바뀌기 때문에 explicit 표현에서는 이 연결 관계를 매번 다시 계산해줘야 합니다. 2차원에서는 그나마 매개변수 구간을 합치거나 나누는 정도로 처리할 수 있지만, 3차원에서는 이 “표면 수술(interface surgery)”이 훨씬 복잡해지고, 구멍이 뚫리는 등의 부정확한 결과가 나오기도 합니다.

이 문제를 피하기 위한 대안이 implicit interface representation입니다. 경계를 직접 나열하는 대신, 어떤 함수 $\phi(\vec{x})$를 하나 정의하고 그 함수의 등고선(isocontour) 중 하나, 보통 $\phi(\vec{x})=0$인 지점들의 집합을 경계로 삼는 방식입니다. 예를 들어 $\phi(x)=x^2-1$이라는 함수를 생각하면, $\phi(x)=0$을 만족하는 점은 정확히 $x=\pm 1$이므로 이것이 곧 인터페이스가 됩니다.

여기서 핵심은 $\phi$가 $n$차원 공간 전체에서 정의되는 반면, 인터페이스는 그보다 한 차원 낮은 $n-1$차원 집합이라는 점입니다. 얼핏 보면 낭비처럼 보입니다. 1차원의 점 두 개를 표현하기 위해 1차원 함수 전체를 들고 오고, 2차원의 곡선 하나를 표현하기 위해 2차원 함수 전체를 쓰는 셈이니까요. 하지만 저자들은 이 “한 차원 높여서 표현하기(embedding)”가 오히려 강력한 도구들을 제공한다고 설명합니다. 가장 대표적인 것이 연결 관계를 신경 쓸 필요가 없다는 점입니다. 3차원 표면을 다룰 때도 균일한 3차원 Cartesian grid 하나만 있으면 되고, 표면이 합쳐지거나 갈라지는 상황도 $\phi$ 값의 부호가 자연스럽게 바뀌는 것으로 처리되기 때문에 별도의 “수술”이 필요 없습니다. 이것이 이 챕터, 그리고 결국 레벨셋 방법 전체가 implicit 표현을 채택하는 근본적인 이유입니다.

Implicit Function으로 인터페이스 다루기

부호로 안과 밖을 구분하기

일단 $\phi(\vec{x})=0$을 인터페이스로 정하고 나면, 아주 자연스러운 성질 하나가 따라옵니다. 임의의 점 $\vec{x}_0$이 내부에 있는지 외부에 있는지를 $\phi(\vec{x}_0)$의 부호만 보고 판단할 수 있다는 것입니다. 논문에서는 $\Omega^- $(내부)를 $\phi<0$인 영역, $\Omega^+$(외부)를 $\phi>0$인 영역으로 정의합니다. explicit 표현에서는 어떤 점이 인터페이스 안쪽에 있는지를 판단하려면 그 점에서 바깥으로 광선을 쏘아 인터페이스와 몇 번 교차하는지(홀수면 내부, 짝수면 외부)를 세는 식의 번거로운 계산이 필요했는데, implicit 표현에서는 그냥 $\phi$ 값 하나를 계산하면 끝입니다. 이산화된 grid에서는 물론 $\phi$ 값이 격자점에서만 알려져 있으므로 보간(interpolation)을 거쳐야 하지만, 이 보간 오차는 인터페이스의 위치를 살짝 이동시키는 정도의 영향만 주고 격자점 수를 늘리면 줄어들기 때문에 대부분의 수치해석 방법이 감내할 수 있는 수준입니다.

Gradient와 법선 벡터

인터페이스 위의 한 점에서 표면이 어느 방향을 향하고 있는지, 즉 법선(normal) 벡터를 구하는 문제를 생각해보겠습니다. 여기서 사용하는 것이 $\phi$의 gradient입니다.

이 gradient는 $\phi$의 등고선에 항상 수직이고, $\phi$ 값이 증가하는 방향을 가리킵니다. 즉, 인터페이스 위의 점 $\vec{x}_0$에서 $\nabla\phi$를 계산하면 그 방향이 바로 바깥쪽을 향하는 법선 방향과 일치합니다. 다만 크기까지 법선 벡터가 되려면 정규화가 필요하므로, 단위 법선은 다음과 같이 정의됩니다.

이 식을 볼 때 주의할 점은, $\nabla\phi$가 인터페이스 위에서만 법선 방향을 나타낸다는 사실입니다. 그런데 저자들은 이 식을 인터페이스 위의 점뿐 아니라 도메인 전체에 대해 적용해서 $\vec{N}$이라는 함수를 정의합니다. 즉, 인터페이스에서 정의되는 개념(법선)을 한 차원 높은 도메인 전체로 “임베딩”하는 것입니다. 이는 앞서 설명한 implicit 표현의 철학과 정확히 같은 맥락으로, 인터페이스에서만 의미를 갖는 정보를 전체 grid 위에서 다룰 수 있게 만들어 계산을 단순화합니다.

실용적으로는 이 gradient를 유한차분(finite difference)으로 근사해야 합니다. 논문은 세 가지 기본 공식을 제시합니다.

forward/backward는 1차 정확도, central은 2차 정확도를 가집니다. 다만 이 세 공식은 모두 분모에 $\lvert\nabla\phi\rvert$가 등장하는 구조와 맞물려 있어서, $\nabla\phi=0$이 되는 특이점(예를 들어 $\phi(x)=x^2-1$에서 $x=0$)에서는 법선이 정의되지 않는다는 문제가 생깁니다. 이 챕터가 제시하는 처리 방식이 흥미로운데, 그냥 임의의 방향을 하나 골라서 법선으로 정해버리면 된다는 것입니다. 언뜻 무책임해 보이지만, 저자들은 이것이 실은 $\phi$ 값에 반올림 오차 수준의 아주 작은 섭동(perturbation)을 주는 것과 동등하다고 설명합니다. 수치해석 방법이 작은 섭동에 안정적이라면(대부분의 well-posed 문제가 그렇습니다), 이런 특이점에서의 임의 선택은 전체 계산에 해를 끼치지 않습니다. 반대로 분모가 0에 가까워질 때 작은 $\epsilon$을 더해서 나누는 흔한 트릭은 오히려 $\lvert\vec{N}\rvert\neq 1$인 이상한 법선을 만들어내므로 피해야 한다고 지적합니다.

Curvature

인터페이스가 얼마나 휘어 있는지를 나타내는 curvature $\kappa$는 법선의 divergence로 정의됩니다.

여기서 $\kappa>0$이면 볼록(convex), $\kappa<0$이면 오목(concave), $\kappa=0$이면 평평한 영역입니다. 법선의 각 성분을 따로 미분해서 계산할 수도 있지만, 식 (1.2)의 $\vec{N}=\nabla\phi/\lvert\nabla\phi\rvert$를 대입하면 $\phi$의 1차, 2차 미분만으로 곡률을 직접 계산하는 식을 얻을 수 있습니다.

식만 보면 복잡해 보이지만, 결국 이 식이 하는 일은 법선을 각 방향으로 미분해서 더하는 계산을 $\phi$의 2차 미분들로 직접 풀어쓴 것뿐입니다. $\phi_{xx}$ 같은 2차 미분은 $D_x^+D_x^-\phi = (\phi_{i+1}-2\phi_i+\phi_{i-1})/\Delta x^2$ 형태의 2차 정확도 중심차분으로 근사합니다.

한 가지 실용적으로 중요한 포인트는, $\phi$가 원점처럼 모든 편미분이 0이 되는 지점 근처에서는 $\kappa\to\infty$로 발산할 수 있다는 것입니다. 예를 들어 원 $\phi=x^2+y^2-1$의 경우 원점에서 $\kappa$가 정의되지 않습니다. 하지만 저자들은 이것이 “임베딩”으로 인한 인공적인 특이점일 뿐, 실제 인터페이스(원 자체)에서는 아무 문제가 없다고 지적합니다. 게다가 Cartesian grid의 최소 해상도가 $\Delta x$인 이상, 그보다 작은 반지름의 곡률을 정확히 표현하는 것 자체가 의미가 없으므로, 실제 구현에서는 $-1/\Delta x\le\kappa\le 1/\Delta x$로 값을 clamp해서 사용합니다.

CSG: Boolean 연산

implicit 표현의 또 다른 장점은 여러 도형의 합집합, 교집합, 여집합 같은 연산(constructive solid geometry, CSG)을 함수 연산만으로 처리할 수 있다는 점입니다. 두 implicit function $\phi_1,\phi_2$가 있을 때,

• 합집합(union): $\phi=\min(\phi_1,\phi_2)$
• 교집합(intersection): $\phi=\max(\phi_1,\phi_2)$
• 여집합(complement): $\phi=-\phi_1$
• 차집합($\phi_1$에서 $\phi_2$의 내부를 뺀 영역): $\phi=\max(\phi_1,-\phi_2)$

직관적으로 보면, $\phi<0$이 “내부”라는 약속을 그대로 이용하는 것입니다. 두 영역 중 하나라도 내부($\phi<0$)이면 합집합의 내부이므로 둘 중 더 작은(더 음수인) 값을 취하는 min이 union이 되고, 두 영역 모두 내부여야 교집합의 내부이므로 max가 intersection이 되는 식입니다. explicit 표현으로 이런 연산을 직접 하려면 곡선/곡면끼리의 교차를 다시 계산해야 하지만, implicit 표현에서는 grid 위의 값끼리 min/max만 취하면 되므로 CAD 같은 응용에서 특히 유용합니다.

Heaviside 함수와 Delta 함수, 그리고 적분

지금까지는 인터페이스의 안/밖을 구분하는 이야기였다면, 이제는 그 위에서 어떻게 적분을 계산할지의 문제입니다. 내부 영역 $\Omega^-$에서 함수 $f$를 적분하려면, 먼저 내부인지 아닌지를 나타내는 characteristic function이 필요합니다.

그런데 이 $\chi^-$는 다차원 변수 $\vec{x}$의 함수입니다. 이를 1차원 변수 $\phi$만의 함수로 바꿔 쓰면 다루기 편해지는데, 그래서 등장하는 것이 Heaviside 함수입니다.

$\chi^+(\vec{x})=H(\phi(\vec{x}))$, $\chi^-(\vec{x})=1-H(\phi(\vec{x}))$로 쓸 수 있으므로, 내부 영역에서의 적분은

로 표현됩니다. 여기서 핵심은 적분 범위가 여전히 도메인 전체 $\Omega$라는 점입니다. 원래는 내부 영역의 정확한 경계를 알아야 그 위에서 적분할 수 있는데, Heaviside 함수를 곱해두면 외부 영역의 기여가 자동으로 0이 되어 사라지므로, 굳이 내부/외부를 나눌 필요 없이 전체 도메인에 대해 적분하면 됩니다.

경계(인터페이스) 자체에서의 적분, 즉 표면적분은 Dirac delta 함수로 표현됩니다. $\delta(\phi)=H’(\phi)$로 정의하면, 표면적분은

가 됩니다. $\delta(\phi)$가 $\phi=0$을 제외한 모든 곳에서 0이기 때문에, 이 적분은 사실상 인터페이스 위에서만 값을 갖고 나머지 영역에서는 자동으로 걸러집니다. 즉, 내부/외부/경계를 각각 따로 다루지 않고 전체 도메인에 대한 하나의 적분식으로 표현할 수 있다는 것이 이 절 전체의 요점입니다. $\lvert\nabla\phi\rvert$ 항이 붙는 이유는, $\delta(\phi(\vec{x}))$가 $\phi$ 공간에서의 델타함수를 $\vec{x}$ 공간으로 변환할 때 생기는 야코비안 보정 때문입니다(연쇄법칙으로 $\nabla H(\phi(\vec{x}))=H’(\phi)\nabla\phi$가 되고, 법선 방향 $\vec{N}=\nabla\phi/\lvert\nabla\phi\rvert$과 내적하면 $\lvert\nabla\phi\rvert$가 남습니다).

실제 계산에서는 $\delta(\phi)$가 측도 0인 집합에서만 0이 아니므로 grid 샘플링만으로는 적분을 근사하기 어렵습니다. 그래서 폭 $\epsilon$만큼 “뭉갠(smeared-out)” 버전의 Heaviside/delta 함수를 사용합니다.

전형적으로 $\epsilon=1.5\Delta x$ 정도를 사용하며, 이는 signed distance function으로 정규화된 $\phi$ 기준으로 인터페이스 폭이 grid 셀 3개 정도가 되게 하는 값입니다. 다만 이 뭉갠 버전은 $O(\Delta x)$ 수준의 1차 정확도만 보장하므로, 더 정밀한 결과가 필요하면 marching cubes 같은 등고선 추출 알고리즘을 별도로 써야 한다고 저자들은 언급합니다.

Signed Distance Function: 더 좋은 φ를 고르는 법

여기까지는 $\phi$가 그저 $\phi=0$에서 인터페이스를 지나가는 어떤 함수이기만 하면 된다는 전제로 이야기를 진행했습니다. 하지만 $\phi$를 어떻게 고르느냐에 따라 수치적 안정성과 정확도가 크게 달라집니다. 예를 들어 $\phi(x)=x^2-1$처럼 gradient가 지나치게 가파르거나 평평해지는 함수는 보간이나 미분 근사에서 오차가 커지기 쉽습니다. 저자들이 제안하는 좋은 선택이 바로 signed distance function입니다.

거리 함수에서 signed distance function으로

먼저 (부호 없는) 거리 함수는 다음과 같이 정의됩니다.

즉 한 점 $\vec{x}$에서 인터페이스 위의 가장 가까운 점까지의 유클리드 거리입니다. 이 함수의 중요한 성질은 $\lvert\nabla d\rvert=1$이라는 것인데, 직관적으로는 “인터페이스에 두 배 가까워지면 $d$ 값도 정확히 절반이 된다”는 의미입니다. 기하학적으로 보면, $\vec{x}$에서 가장 가까운 인터페이스 점 $\vec{x}_C$까지 잇는 직선 경로가 거리 함수 $d$가 가장 가파르게 감소하는 방향(steepest descent)이 되고, 이 경로를 따라가면 이동 거리와 $d$ 값의 변화량이 정확히 1:1로 대응하기 때문입니다. (단, 인터페이스 위에서 가장 가까운 점이 유일하지 않은 등거리(equidistant) 점들에서는 이 성질이 깨집니다.)

이 거리 함수에 부호를 붙인 것이 signed distance function입니다.

일반 implicit function과 signed distance function을 가르는 결정적 조건은 $\lvert\nabla\phi(\vec{x})\rvert=1$입니다. 예를 들어 1차원에서 $\phi(x)=x^2-1$ 대신 $\phi(x)=\lvert x\rvert-1$을 쓰면 같은 인터페이스($x=\pm1$)를 나타내면서도 $\lvert\nabla\phi\rvert=1$을 만족하는 signed distance function이 됩니다. 2차원 원의 경우도 $\phi=x^2+y^2-1$ 대신 $\phi=\sqrt{x^2+y^2}-1$을 쓰면 됩니다.

거리 함수와 signed distance function의 차이는 미분 가능성에서 드러납니다. (부호 없는) 거리 함수는 인터페이스에서 $d=0$이 국소 최솟값이 되므로 그 지점에 뾰족한 꺾임(kink)이 생기고, 이 근방에서 미분을 근사하기가 까다롭습니다. 반면 signed distance function은 인터페이스를 가로지르며 단조(monotonic)롭게 변하기 때문에 인터페이스 근방에서도 훨씬 안정적으로 미분할 수 있습니다. 이것이 레벨셋 방법이 부호 없는 거리 함수 대신 signed distance function을 선호하는 이유입니다.

signed distance function은 또한 가장 가까운 인터페이스 점을 직접 계산하는 데도 쓸 수 있습니다.

$\phi(\vec{x})$가 이미 부호와 크기를 모두 담은 “인터페이스까지의 거리”이므로, 현재 위치에서 법선 방향으로 그만큼 되돌아가면 인터페이스 위의 가장 가까운 점에 도달한다는, 상당히 직관적인 식입니다. 다만 이 역시 등거리점 근처에서는 근사적으로만 성립합니다.

논문은 이 대목에서 일반적인 수치해석 원칙 하나를 짚고 넘어갑니다. “일반적으로만(generally) 참인 식이라도, 그것이 실패하는 방식이 우아하다면(graceful) 수치적 근사에 사용해도 괜찮다”는 것입니다. 즉 $\lvert\nabla\phi\rvert=1$이나 $\vec{x}_C=\vec{x}-\phi\vec{N}$ 같은 식들이 특이점에서 정확히 성립하지 않더라도, 그 실패가 전체 계산을 무너뜨리지 않는 한 실용적으로 사용할 수 있다는 관점입니다. 앞서 법선 계산에서 분모가 0일 때 임의의 방향을 선택해도 괜찮다고 한 것과 같은 맥락입니다.

Kink는 왜 실제로 문제가 되지 않는가

$\phi(x)=\lvert x\rvert-1$은 $x=0$에서 미분 불가능한 꺮임을 가집니다. 이론적으로는 문제처럼 보이지만, 저자들은 실제 Cartesian grid에서는 표본점(sample point)들만 존재하므로 이 꺾임 자체를 “해상”할 방법이 없다고 지적합니다. grid 위에서 이 꺾임은 수치적으로 자연스럽게 뭉개지고, 그 지점에서의 유한차분 값은 $[-1,1]$ 구간의 어떤 유한한 값이 될 뿐입니다. 즉 최악의 경우라 해봐야 gradient가 0이 되는 것이고, 이는 앞서 다룬 법선의 특이점 처리(임의 방향 선택)로 이미 대응 가능한 상황입니다.

더 일반적으로는, 1차원의 꺾임이 2차원에서는 꺾임의 선(line of kinks)으로, 3차원에서는 꺾임의 평면(plane of kinks)으로 확장될 수 있습니다. 예컨대 3차원에서 $\phi(\vec{x})=\lvert x\rvert-1$은 두 평면 $x=\pm1$을 나타내는데, $x=0$인 평면 전체가 꺾임입니다. 이 역시 grid 위에서는 수치적으로 뭉개져서 별도의 특수 처리 없이 넘어갈 수 있다고 설명합니다. 다만 뭉개진 꺾임 근처에서는 국소적으로 $\lvert\nabla\phi\rvert\neq 1$이 되므로, signed distance function이라는 가정에 의존하는 공식들을 쓸 때는 주의가 필요합니다.

Signed Distance Function에서 더 단순해지는 공식들

$\lvert\nabla\phi\rvert=1$이 (거의) 항상 성립한다는 사실은 앞서 다룬 법선과 곡률 공식을 정규화 상수 없이 단순화시켜줍니다.

여기서 $\Delta\phi=\phi_{xx}+\phi_{yy}+\phi_{zz}$는 $\phi$의 Laplacian입니다. 일반 implicit function에서 곡률을 구하려면 앞서 본 것처럼 1차, 2차 미분이 뒤섞인 복잡한 분수식 (1.8)을 계산해야 했는데, signed distance function을 쓰면 이것이 단순히 Laplacian 하나로 줄어듭니다. 표면적분에 등장하던 $\lvert\nabla\phi\rvert$ 항도 1이므로 생략할 수 있어서, $\int_\Omega f(\vec{x})\,\delta(\phi(\vec{x}))\,d\vec{x}$처럼 더 단순한 형태가 됩니다.

다만 저자들은 이 단순화된 공식들을 과신하지 말라고 경고합니다. 실제 수치 계산에서는 뭉개진 꺾임 등의 이유로 $\lvert\nabla\phi\rvert=1$이 정확히 성립하지 않는 경우가 많기 때문에, 식 (2.5)와 (2.6)이 일반적으로만 참이라는 사실을 항상 염두에 두어야 합니다. 오히려 어떤 상황에서는 법선이나 곡률 자체보다 $\nabla\phi$나 $\Delta\phi$를 직접 사용하는 편이 더 적절할 수 있다고 언급합니다.

Comment

이 두 챕터를 다 읽고 나서 기억해야 할 것은 결국 하나입니다. “인터페이스를 명시적으로 들고 다니지 말고, 한 차원 높은 공간에서 스칼라 함수의 등고선으로 표현하자”는 결정 하나가 연결 관계 추적, 안/밖 판정, 법선·곡률 계산, 적분까지 이후에 등장하는 거의 모든 문제를 재구성한다는 점입니다. explicit 표현이 차원마다 별도의 도구(1차원은 점, 2차원은 매개변수 곡선, 3차원은 삼각형 메쉬와 연결 관계)를 필요로 하는 반면, implicit 표현은 차원에 상관없이 “Cartesian grid 위에 스칼라 값 저장하기”라는 동일한 틀로 문제를 통일합니다. 이 통일성이야말로 저자들이 반복해서 강조하는 implicit 표현의 진짜 가치라고 생각합니다.

signed distance function 챕터에서 특히 인상 깊은 부분은 “일반적으로만 참인 식(generally true equation)”을 다루는 태도입니다. $\lvert\nabla\phi\rvert=1$이나 법선 계산 공식이 특이점에서 엄밀하게는 깨지는데도, 저자들은 이를 억지로 고치려 하지 않고 “그 실패가 우아하게(gracefully) 일어나는가”만을 기준으로 삼습니다. 이는 이론적 엄밀함보다 수치적 견고함을 우선하는 실용적 태도이며, 이후 챕터들에서 다룰 (아마도) upwind scheme이나 재초기화(reinitialization) 같은 기법들이 왜 그런 형태를 취하는지 이해하는 데도 도움이 되는 사고방식이라고 봅니다.

다만 이 챕터 자체는 순수하게 정의와 성질을 나열하는 도입부이므로, 여기서 소개된 signed distance function이 실제로 어떻게 유지·재계산되는지(즉 $\phi$가 시간에 따라 움직이면서도 계속 $\lvert\nabla\phi\rvert=1$을 근사적으로 유지하도록 만드는 절차)는 다뤄지지 않습니다. 이 부분은 책의 뒷장에서 다뤄질 내용으로 보이며, 지금 단계에서는 “왜 signed distance function이라는 특정한 $\phi$를 선호하는가”에 대한 동기만 확실히 챙기고 넘어가면 충분합니다. implicit function과 레벨셋 방법을 처음 접하는 사람이라면, 이 챕터의 내용만으로도 이후 등장하는 PDE 기반 인터페이스 진화 방법들을 이해하기 위한 언어(부호, gradient, 법선, 곡률, Heaviside/delta 함수)를 충분히 갖출 수 있습니다.