Level Set 方法

Level Set方法是一种处理几何体形状变化的计算机图形学算法。

与直观的基于“点、线、面”的处理几何元素不同，Level Set方法基于隐式方法表达和处理几何图形，将$\mathbb{R}$维空间内的一个几何图形，表示为$\mathbb{R}+1$维空间中的一个等值线/等值面。

这里所说的”几何图形“，特指一般的、不规则的、难以使用参数方程来描述的几何图形。对于可以用解析方法简单描述的图形，直接用解析方法来处理即可。然而，工程中接触的图形，多难以使用解析方法描述，或解析方法过于复杂而失去了通用性。

In [1]:
import matplotlib.pyplot as pyplot
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pylab as pylab
%matplotlib notebook
In [2]:
import numpy as np
x_length=1
y_length=1

x_dim=100
y_dim=100

circle_x = 0.4
circle_y = 0.6
circle_r = 0.3
In [3]:
"""
一个圆形，以及其在三维空间中的等值面
圆形是一个二维图形
其在Level Set方法中，表现为三维空间内某个曲面与 z=0 截面的交线

Level Set方法译为水平集方法
"""

theta = np.linspace(0,2*np.pi,360)
x=circle_r * np.cos(theta) + circle_x
y=circle_r * np.sin(theta) + circle_y


gx, gy = np.meshgrid(np.linspace(0,x_length,x_dim),np.linspace(0,y_length,y_dim))

# signed Distance Field (SDF) 有向距离场
sdf_circle = np.sqrt((gx-circle_x)**2 + (gy-circle_y)**2) - circle_r



fig=pyplot.figure()
ax = fig.add_subplot(111, projection='3d')

ax.plot(x,y,zs=0,zdir='z')

ax.plot_surface(X=gx,Y=gy,Z=sdf_circle,rstride=5, cstride=5)
Out[3]:
<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7f36f9f66860>

用等值面来表示几何图形，有以下优点：

  • 通用性好：无论何种图形，映射到高维空间内都是一组仅包含标量的点，故对于各种图形，我们仅需要开发一套算法
  • 适应力强：图形的拓扑结构变化（分离/合并/特征的消失与出现），反映在高维空间内，仅仅是不同的点具有了不同的数值，不需要特别处理

缺点只有一个：

  • 计算量大：由于需要处理的数据高了一个维度，数据量大增。同时，需要迭代求解一组偏微分方程。

水平集方法

水平集方法是一种隐式描述移动界面的欧拉方法， 最早由 Osher 和 Sethian 在1988年引入。

水平集方程的推导

在区域$\Omega$内，考虑随时间$t$演化的界面$\Gamma=\Gamma(x,t)$，其中$x$表示$\Gamma$上的一点。移动界面可以用一个函数$\Phi(x,t)$的零水平集表示, 即$\Gamma=\{x|\Phi(x,t)=0\}$。给定$\Omega$上的速度向量场$\vec{V}(x,t)$,则方程$(1)$描述了$\Phi$的水平集以速度$\vec{V}(x,t)$在区域$\Omega$中传播的过程，该方程就称为水平集方程。

$${\Phi}_t + V \cdot \nabla\Phi = 0 \tag{1}$$

二维情况

假设我们关注$\Phi(x,y,t)=0$所表示的移动界面。这意味着界面上的点在移动过程中，始终符合$\Phi(x,y,t)=0$.如此则有：

$$\frac{d}{dt}{\Phi\left(x(t),y(t),t\right)}=0$$

应用链式求导法则：

$${\Phi}_t+{{\Phi}_x}{x_t} + {{\Phi}_y}{y_t} = 0$$$${\Phi}_t+ (x_t,y_t)\cdot ({\Phi}_x, {\Phi}_y) = 0$$$${\Phi}_t+ \vec{V}\cdot{\nabla{\Phi}}=0 \tag{2}$$

其中$\vec{V}=(x_y,y_t)$.这就是我们需要求解的水平及方程。方程$(2)$描述了水平集在外部速度场中的移动，建立了隐式表达与显式表达之间的桥梁。

进一步看，水平集曲面的法向量可以如下表示：

$$\vec{N} = \frac{\nabla\Phi(x,y,t)}{\left|\nabla\Phi(x,y,t)\right|}$$

观察方程$(2)$，$\vec{V}\cdot{\nabla{\Phi}}$这一项可以理解为：这些点在水平集曲面法向上的分量才会影响曲面的形状。切向上面的速度只是在曲面内部平移。

法向上的分量可以表示为：

$${V_n} = \vec{V}\cdot{\vec{N}} = \vec{V}{\cdot} \frac{\nabla\Phi(x,y,t)}{\left|\nabla\Phi(x,y,t)\right|}$$

继续变换：

$${V_n} {\left|\nabla\Phi(x,y,t)\right|} = \vec{V}{\cdot} \nabla\Phi(x,y,t)$$

将这一项代入方程$(2)$，得到用于计算的方程：

$${\Phi}_t + {V_n} {\left|\nabla\Phi(x,y,t)\right|}=0$$

这是一个标量方程，很容易用数学方法求解。

${\left|\nabla\Phi(x,y,t)\right|}$一项，在标准的有向距离场中应该恒等于1，但非均匀速度场下，每一步之后该特征就会被破坏，需要进行重新初始化操作。如果速度场式均匀的，则方程变为：

$${\Phi}_t + {V_n}=0$$$$\Phi = C + t{V_n}$$

有解析解。实际上Level Set方法退化为最小距离场方法。

Level Set 方法原始文献：

Osher, Stanley, and James A Sethian. “Fronts Propagating with Curvature-Dependent Speed: Algorithms Based on Hamilton-Jacobi Formulations.” Journal of Computational Physics 79, no. 1: 12–49.

本文参考了：

http://weihuayi.github.io/numerical/Levelset-method.html

Level Set Methom (LSM) 的两个变种：最小距离场法 Minimum Distance Method (MDM) 与快速推进法 Fast Marching Method (FMM).

  • MDM：假设速度场$\vec{V}$是均匀的，即在水平集曲面在所有点处的法向速度相同。MDM是在此前提下对LSM的简化。
  • FMM：速度场$\vec{V}$恒为负（几何体在扩张）前提下对LSM的优化。
In [4]:
# 下面，对一个圆形和一个多变形的扩张，使用LSM进行仿真

# 工具函数们
import sdf

# 一个多边形
polygon_array = np.asarray([[0.1, 0.1],
                            [0.2, 0.4],
                            [0.5, 0.8],
                            [0.5, 0.5],
                            [0.6, 0.1]])

# 计算SDF
gx, gy = np.meshgrid(np.linspace(0,x_length,x_dim),np.linspace(0,y_length,y_dim))

# 计算一条边的SDF
num_edges = polygon_array.shape[0]
df = sdf.generate_df_for_single_edge(
    gx,gy, # 网格点坐标
    polygon_array[num_edges - 1, 0], polygon_array[num_edges - 1, 1], # 边（线段）的起点坐标
    polygon_array[0, 0], polygon_array[0, 1] # 边（线段）的终点坐标
)
# 循环计算其他边
for i in range(num_edges - 1):
    df_temp = sdf.generate_df_for_single_edge(gx, gy,
                                              polygon_array[i, 0], polygon_array[i, 1],
                                              polygon_array[i + 1, 0], polygon_array[i + 1, 1]
                                             )
    # 对两个（有向）距离场取最小值，实质上是几何上的的布尔并(OR)运算
    df = np.minimum(df, df_temp)

# 计算符号场
sf = sdf.generate_sf_for_polygon(gx, gy, polygon_array)

# 有向距离场
sdf_final = sf*df

# 与圆做布尔并
circle_r = 0.1 # 圆调整一下，否则扩张过程看起来不明显
circle_x = 0.8
circle_y = 0.8
sdf_circle = np.sqrt((gx-circle_x)**2 + (gy-circle_y)**2) - circle_r
sdf_final = np.minimum(sdf_circle, sdf_final)

# 查看SDF

fig=pyplot.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X=gx,Y=gy,Z=sdf_final,rstride=5, cstride=5)
Out[4]:
<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7f36f9e32ef0>
In [5]:
# 以及当前的等值线

fig=pyplot.figure()
ax = fig.add_subplot(111)
ax.contour(gx, gy, sdf_final, [0], colors='k')
Out[5]:
<matplotlib.contour.QuadContourSet at 0x7f36f9dc1240>
In [6]:
# 计算数据

steps = 20

Vn = -1 # 负的速度是扩张

CFL = 0.8
dx = np.sqrt((x_length/(x_dim-1))**2 + (y_length/(y_dim-1))**2)

# 速度场
v = np.zeros(gx.shape) + dx * Vn * CFL

# 分配一个保存数据用的list
data = [sdf_final]

# 计算
for s in range(steps):
    # 假装这里进行了重新初始化（本例实际上不需要）
    
    # 假装计算 $\nabla\Phi$（本例中恒等于1）（参考资料：WENO5，Gudunov格式）
    nabla = np.ones(gx.shape)
    
    # 时间推进（最好用TVD-RK，本例只用一阶）
    sdf_final = sdf_final + v * nabla
    
    # 保存中间数据
    data.append(sdf_final)
In [7]:
from matplotlib import animation
# 显示数据
fig = pyplot.figure()
ax = pyplot.axes()
curve = ax.contour(gx, gy, data[0], [0], colors='k')


def animate(d):
    ax.clear()
    curve = ax.contour(gx, gy, d, [0], colors='k')
    return curve,

anim = animation.FuncAnimation(fig, animate,frames = data, interval=200, blit=True)

anim.save('result.mp4', fps=20, extra_args=['-vcodec', 'libx264'])