level set method 水平集方法基本问题
2014-03-13 13:55
232 查看
版里搜了下,呵呵,发现level set method应该是从图像识别领域借过来的。基本思想: 令(n+1)维函数 t = phi(x(t)) 等于0,即水平集。
基本问题:
1 初始化
定义 level set 函数 phi(x) 为符号距离函数,界面线上距离为0。一般为隐函数。
1.1 符号判定, 界面闭合线内/外
1.1.1 2D留数算法: 将边界C离散成n段,c_1, c_2, .. c_n 为结点, oc_i 表示某点o 指向c_i的向量
sum ( atan(( oc_i x oc_(i+1)) / (oc_i . oc_(i+1))), i = 1: n)
if sum == 0 ---> 点o 在边界C外
if sum == 2 * pi ---> 点o在边界C内
1.1.2 交点判定算法: 给定已知符号的点P0, 待定符号的点Px 与 P0 之间连线穿过边界C,
如果穿过次数是偶数表示Px与P0符号相同;
如果是奇数,则符号相反。
实现:
将边界C离散成n段,初始定义所有网格点(2D面)符号为正(C外),当离散的边界线段,与某条垂向网格线相交,则将该网格线上,在边界线段上侧的点全部变负号。
1.2 距离计算
1.2.1 边界C可由参数方程表示,先搜索距离边界线最近的网格点,计算距离;
1.2.2 边界C分段连续,直接计算各个端点到网格点的距离
2 演化方程
2.1 边界线随时间的变化,采用PDE描述 phi(x(t), t) = 0
d phi(x)/ dt + grad_phi (a. n + b.t) = 0
其中n是法向单位向量,t是切向单位向量。a是法向速度,b是切向速度。n影响边界线拓扑。
<- - n = grad_ phi / norm( grad_phi)
--> d phi(x)/ dt = a norm(grad_phi) = 0
2.2 演化方程的数值解
演化方程具有 Hamilton-Jacobi形式, d phi /dt + H(grad_phi) = 0 。 这类方程不存在经典解,只有弱解。解C0连续但不保证C1连续,即使初值C1连续。守恒型双曲PDE也有类似的解性质。比如超声速压缩流动,给定连续的初状态值,仍然可能演化出不连续的间断激波。所以演化方程的高阶精度数值解可借鉴守恒型双曲方程的ENO, WENO等解方案。LSM演化方程实际上是就是一个包含源项的守恒型双曲方程。
常用一阶精度的演化方程常采用迎风格式,Godunov格式, Lax_Friedrich格式等。
2.2.1 2D Godunov格式
norm(grad_phi) ~= ( max(-D_{+x},D_{-x}, 0) ^2 + max(-D_{+y}, D_{-y}, 0)^2 ) ^.5;
D_{+x} forward差分, D_{-x} backward差分。同理 D_{+y}, D_{-y}。 物理意义上要求保证信息必须从已知传向未知,所以对于每个当前要求的phi, 必须选择传播正方向上的差分格式。
一个简单的判断:当地(local) phi_x (i,j) 大于还是小于0. 如果 大于0,表示phi(i,j)沿x轴正方向传播,显然先要经过phi(i-1, j),故应该选用 D_{-x}差分。而且此时还满足
abs( D_{-x}) > abs(D_{+x}), 所以才有统一的上式样子。同理,如果local phi_x(i,j) 小于0, 则表示 phi(i,j )沿x轴负方向传播,即先经过 phi(i+1, j), 故采用 D_{+x}。同理y方向。
MATALB CODE here
2.2.2 ENO格式
参见 力学所李新亮老师CFD讲义, 激波平滑。
2.2 法向速度延扩
从边界线到全水平集面,任意垂直边界线的直线上速度恒定 --》 所以水平集面上任意一点的速度,可由两条垂直边界线的直线l1, l2唯一确定
2.3 边界设定
布置ghost point, 速度边界Ok, 位移边界可能出问题 -- > 重初始化(重新选择一个level set funtion)
假定周期边界 --》 保证单元边界符号连续
3 窄带level set method
只跟新距离边界线一定范围的网格点,远场点保持不变。节省计算/存储。
预告:下篇介绍 reinitialization in level set method
基本问题:
1 初始化
定义 level set 函数 phi(x) 为符号距离函数,界面线上距离为0。一般为隐函数。
1.1 符号判定, 界面闭合线内/外
1.1.1 2D留数算法: 将边界C离散成n段,c_1, c_2, .. c_n 为结点, oc_i 表示某点o 指向c_i的向量
sum ( atan(( oc_i x oc_(i+1)) / (oc_i . oc_(i+1))), i = 1: n)
if sum == 0 ---> 点o 在边界C外
if sum == 2 * pi ---> 点o在边界C内
1.1.2 交点判定算法: 给定已知符号的点P0, 待定符号的点Px 与 P0 之间连线穿过边界C,
如果穿过次数是偶数表示Px与P0符号相同;
如果是奇数,则符号相反。
实现:
将边界C离散成n段,初始定义所有网格点(2D面)符号为正(C外),当离散的边界线段,与某条垂向网格线相交,则将该网格线上,在边界线段上侧的点全部变负号。
1.2 距离计算
1.2.1 边界C可由参数方程表示,先搜索距离边界线最近的网格点,计算距离;
1.2.2 边界C分段连续,直接计算各个端点到网格点的距离
2 演化方程
2.1 边界线随时间的变化,采用PDE描述 phi(x(t), t) = 0
d phi(x)/ dt + grad_phi (a. n + b.t) = 0
其中n是法向单位向量,t是切向单位向量。a是法向速度,b是切向速度。n影响边界线拓扑。
<- - n = grad_ phi / norm( grad_phi)
--> d phi(x)/ dt = a norm(grad_phi) = 0
2.2 演化方程的数值解
演化方程具有 Hamilton-Jacobi形式, d phi /dt + H(grad_phi) = 0 。 这类方程不存在经典解,只有弱解。解C0连续但不保证C1连续,即使初值C1连续。守恒型双曲PDE也有类似的解性质。比如超声速压缩流动,给定连续的初状态值,仍然可能演化出不连续的间断激波。所以演化方程的高阶精度数值解可借鉴守恒型双曲方程的ENO, WENO等解方案。LSM演化方程实际上是就是一个包含源项的守恒型双曲方程。
常用一阶精度的演化方程常采用迎风格式,Godunov格式, Lax_Friedrich格式等。
2.2.1 2D Godunov格式
norm(grad_phi) ~= ( max(-D_{+x},D_{-x}, 0) ^2 + max(-D_{+y}, D_{-y}, 0)^2 ) ^.5;
D_{+x} forward差分, D_{-x} backward差分。同理 D_{+y}, D_{-y}。 物理意义上要求保证信息必须从已知传向未知,所以对于每个当前要求的phi, 必须选择传播正方向上的差分格式。
一个简单的判断:当地(local) phi_x (i,j) 大于还是小于0. 如果 大于0,表示phi(i,j)沿x轴正方向传播,显然先要经过phi(i-1, j),故应该选用 D_{-x}差分。而且此时还满足
abs( D_{-x}) > abs(D_{+x}), 所以才有统一的上式样子。同理,如果local phi_x(i,j) 小于0, 则表示 phi(i,j )沿x轴负方向传播,即先经过 phi(i+1, j), 故采用 D_{+x}。同理y方向。
MATALB CODE here
2.2.2 ENO格式
参见 力学所李新亮老师CFD讲义, 激波平滑。
2.2 法向速度延扩
从边界线到全水平集面,任意垂直边界线的直线上速度恒定 --》 所以水平集面上任意一点的速度,可由两条垂直边界线的直线l1, l2唯一确定
2.3 边界设定
布置ghost point, 速度边界Ok, 位移边界可能出问题 -- > 重初始化(重新选择一个level set funtion)
假定周期边界 --》 保证单元边界符号连续
3 窄带level set method
只跟新距离边界线一定范围的网格点,远场点保持不变。节省计算/存储。
预告:下篇介绍 reinitialization in level set method
相关文章推荐
- 水平集方法 level set method
- level set method reinitialization (水平集方法的重初始化)
- level-set method vs. immersed boundary method (水平集方法与沉浸元方法比较)
- ldd3 misc-progs/setlevel.c编译问题的解决
- 解决textview.setMovementMethod后Listview的OnItemClickListener点击失效问题
- method setOnKeyListener的问题
- 解决问题: .Net异常 Method not found: Void System.Web.UI.WebControls.BaseDataList.set_Caption(System.String).
- Particle Level Set + Lattice Boltzmann Method
- The method setClass(Context, Class<?>) in the type Intent is not applicable for the问题
- implementation of the particle levelset method
- 论文阅读:Touching Tooth Segmentation from CT Image Sequences Using Coupled Level Set Method
- C/C++基础:set的基本操作总结及注意细节问题
- 解决listview中的textview设置了setMovementMethod导致onItemclick无效的问题
- 自定义颜色的setTextColor问题
- 如何回答关于List和Set的问题
- 分组背包问题Matlab实现——之基本背包问题
- hadoop+eclipse 调试程序出现的问题warn: no job jar file set 以及 点击run on hadoop 无反应
- response.setHeader(...)下载中文文件名乱码问题(转于九万里——城市夜空‘s Blog)
- AdaultBird--我的黑马程序员之路!Chapter3---多线程同步,同步唤醒和等待,守护线程,String类,基本数据类型包装类,集合框架中的List和Set