minimal surface - 演算法筆記

polyline smoothing🚧

polyline smoothing

「平滑化」。當折線太亂，我們可以移動頂點、截彎取直。

我使用兩種函式庫做特徵分解。左邊jsfeat，右邊math.js。

演算法（midpoint smoothing）

每個點移至兩個鄰點的中點。

x̕ᵢ = ½ (xᵢ₋₁ + xᵢ₊₁)

X座標、Y座標、Z座標，互不干涉，分頭計算。

兩端點保持不動，避免折線縮緊。

實施越多次，折線越平滑。

演算法（Laplacian smoothing）

梯散：一個點到鄰點的中點的差距。梯散的本質：凹陷程度。

∆xᵢ = ½ [(xᵢ₋₁ - xᵢ) + (xᵢ₊₁ - xᵢ)] = ½ (xᵢ₋₁ + xᵢ₊₁) - xᵢ

平滑：一個點移至鄰點的中點。平滑的本質：填坑。

x̕ᵢ = xᵢ + ∆xᵢ

midpoint smoothing與Laplacian smoothing等價。

演算法（diffusion smoothing）

朝著中點移動，自訂比例0≤α≤1。0不動、1抵達中點。

x̕ᵢ = xᵢ + α ∆xᵢ   (0 ≤ α ≤ 1)

避免折線前後震盪、來回抖動。

引入頻譜（spectral analysis）

化作矩陣。折線平滑K次，即是矩陣的K次方。

特徵分解。折線平滑K次，即是特徵值的K次方。

化作矩陣O(N²)，特徵分解O(N³T)，折線平滑K次O(N² + NlogK)。適用於K很大的情況。

N = 5
[ x̕₀ ]   [ 1 0 0 0 0 ] [ x₀ ]
[ x̕₁ ]   [ ½ 0 ½ 0 0 ] [ x₁ ]
[ x̕₂ ] = [ 0 ½ 0 ½ 0 ] [ x₂ ]
[ x̕₃ ]   [ 0 0 ½ 0 ½ ] [ x₃ ]
[ x̕₄ ]   [ 0 0 0 0 1 ] [ x₄ ]

未釘住端點、釘住端點、封閉折線，矩陣稍微不同。

開放折線        開放折線        封閉折線
未釘住端點      釘住端點        左右循環
[ 0 1 0 0 0 ]   [ 1 0 0 0 0 ]   [ 0 ½ 0 0 ½ ]
[ ½ 0 ½ 0 0 ]   [ ½ 0 ½ 0 0 ]   [ ½ 0 ½ 0 0 ]
[ 0 ½ 0 ½ 0 ]   [ 0 ½ 0 ½ 0 ]   [ 0 ½ 0 ½ 0 ]
[ 0 0 ½ 0 ½ ]   [ 0 0 ½ 0 ½ ]   [ 0 0 ½ 0 ½ ]
[ 0 0 0 1 0 ]   [ 0 0 0 0 1 ]   [ ½ 0 0 ½ 0 ]

順帶一提，封閉折線的情況下，恰好是對稱矩陣。對稱矩陣擁有正規正交特徵基底，容易求得特徵分解。

順帶一提，這些特殊矩陣，數學家已經推導出特徵向量暨特徵值的公式，不必實施特徵分解演算法。詳情請見後面小節的超連結。

演算法（spectral smoothing）

特徵向量恰是cosine波，特徵值大小恰好符合頻率高低。

低頻是主體，高頻是細節。保留低頻的特徵向量的特徵值，其餘特徵值歸零。保留主體、消除細節。

數學式子

每個點移至兩個鄰點的中點。

x̕ᵢ = ½ (xᵢ₋₁ + xᵢ₊₁)

數學式子改寫成Laplace operator。

x̕ᵢ = xᵢ + ∆xᵢ

數學式子改寫成Laplacian matrix。

x̕ = D⁻¹Ax = (I - D⁻¹L)x = x - D⁻¹Lx

    [ 0 1 0 0 0 ]       [ 1 0 0 0 0 ]       [  1 -1  0  0  0 ]
    [ 1 0 1 0 0 ]       [ 0 2 0 0 0 ]       [ -1  1 -1  0  0 ]
A = [ 0 1 0 1 0 ]   D = [ 0 0 2 0 0 ]   L = [  0 -1  1 -1  0 ]
    [ 0 0 1 0 1 ]       [ 0 0 0 2 0 ]       [  0  0 -1  1 -1 ]
    [ 0 0 0 1 0 ]       [ 0 0 0 0 1 ]       [  0  0  0  1 -1 ]
 adjacency matrix       degree matrix         Laplacian matrix

數學式子改寫成normalized Laplacian matrix。

x̕ = Âx = (I - L̂)x = x - L̂x    where Â = D⁻¹A , L̂ = D⁻¹L

    [ 0 1 0 0 0 ]       [  1 -1  0  0  0 ]
    [ ½ 0 ½ 0 0 ]       [ -½  ½ -½  0  0 ]
Â = [ 0 ½ 0 ½ 0 ]   L̂ = [  0 -½  ½ -½  0 ]
    [ 0 0 ½ 0 ½ ]       [  0  0 -½  ½ -½ ]
    [ 0 0 0 1 0 ]       [  0  0  0  1 -1 ]
       normalized               normalized
 adjacency matrix         Laplacian matrix

未釘住端點、釘住端點、封閉折線，矩陣稍微不同。

開放折線（未釘住端點）
    [ 0 1 0 0 0 ]       [ 1 0 0 0 0 ]       [  1 -1  0  0  0 ]
    [ 1 0 1 0 0 ]       [ 0 2 0 0 0 ]       [ -1  1 -1  0  0 ]
A = [ 0 1 0 1 0 ]   D = [ 0 0 2 0 0 ]   L = [  0 -1  1 -1  0 ]
    [ 0 0 1 0 1 ]       [ 0 0 0 2 0 ]       [  0  0 -1  1 -1 ]
    [ 0 0 0 1 0 ]       [ 0 0 0 0 1 ]       [  0  0  0  1 -1 ]
 adjacency matrix       degree matrix         Laplacian matrix

開放折線（釘住端點）
    [ 1 0 0 0 0 ]       [ 2 0 0 0 0 ]       [  0  0  0  0  0 ]
    [ 1 0 1 0 0 ]       [ 0 2 0 0 0 ]       [ -1  1 -1  0  0 ]
A = [ 0 1 0 1 0 ]   D = [ 0 0 2 0 0 ]   L = [  0 -1  1 -1  0 ]
    [ 0 0 1 0 1 ]       [ 0 0 0 2 0 ]       [  0  0 -1  1 -1 ]
    [ 0 0 0 0 1 ]       [ 0 0 0 0 2 ]       [  0  0  0  0  0 ]
 adjacency matrix       degree matrix         Laplacian matrix

封閉折線
    [ 0 1 0 0 1 ]       [ 2 0 0 0 0 ]       [  1 -1  0  0  1 ]
    [ 1 0 1 0 0 ]       [ 0 2 0 0 0 ]       [ -1  1 -1  0  0 ]
A = [ 0 1 0 1 0 ]   D = [ 0 0 2 0 0 ]   L = [  0 -1  1 -1  0 ]
    [ 0 0 1 0 1 ]       [ 0 0 0 2 0 ]       [  0  0 -1  1 -1 ]
    [ 1 0 0 1 0 ]       [ 0 0 0 0 2 ]       [  1  0  0 -1  1 ]
 adjacency matrix       degree matrix         Laplacian matrix

數學式子

化作矩陣。折線平滑K次，即是矩陣的K次方。

x̕ = Âᵏx = (I - L̂)ᵏx

特徵分解。折線平滑K次，即是特徵值的K次方。

x̕ = EΛᵏE⁻¹x = E(I - Γ)ᵏE⁻¹x     where Â = EΛE⁻¹ , L̂ = EΓE⁻¹

自訂比例0≤α≤1。避免折線前後震盪、來回抖動。

x̕ = (I - αL̂)ᵏx = E(I - αΓ)ᵏE⁻¹x     where L̂ = EΓE⁻¹

L̂的特徵向量恰是cosine波，特徵值大小恰好符合頻率高低。

特徵向量暨特徵值公式（矩陣並未正規化）
https://math.stackexchange.com/questions/3354537/
https://www.cs.yale.edu/homes/spielman/561/2009/lect02-09.pdf

開放折線（未釘住端點）：path graph
開放折線（釘住端點）：查無資料
封閉折線：ring graph

低頻是主體，高頻是細節。保留L̃較大的特徵值暨特徵向量，其餘特徵值歸零。導致(I - αL̂)僅減去細節、未減去主體。

梯度、梯散

章節標題本來是梯散。為了方便講解，我利用了平滑化來引入梯散。

梯度、梯散是微分幾何的基本運算。梯散擁有許多意義，諸如離散曲率。梯散也有許多功能，例如平滑化、測地距離。

mesh smoothing (vertex)🚧

vertex smoothing

當網格太亂，我們可以移動頂點、截彎取直。

基本原理仍是梯散，方法如同先前章節。但是本章節打算深入介紹更多知識。

two perspective:  midpoint / laplacian
two domain:       spatial analysis / spectral analysis
two formation:    quadratic form optimization / linear equation
two method:       iteration / resolution
two objective:    minimize edge length / minimize area
two area formula: circumcenter / incenter
two location:     midpoint / perpendicular foot

兩種觀點：走往中點、加上梯散

x̕ᵢ = sum xⱼ / sum j   midpoint
x̕ᵢ = xᵢ + ∆xᵢ         laplacian

兩種空間：空域、頻域

計算梯散，可以在空域或頻域。實務上大家使用空域。

curve:    t->x(t) t->y(t) t->z(t) t=0⋯1     3 function
polyline: i->x[i] i->y[i] i->z[i] i=0⋯n-1   3 array
surface:  uv->x(u,v) uv->y(u,v) uv->z(u,v)  three 2d function
mesh:     uv merge to i

eigenfunction of laplace operator of curve:    t->f(i)   infinite functions
                                                         (plot as nodal sets)
eigenvector   of laplace operator of polyline: i->f(i)   n vectors
                                                         (sin/cos wave)

eigenfunction/vector plot on curve/polyline as various colors. [-1,+1]->[0,255]
nodal set of eigenfunction -----> meaningless?
nodal set of Fiedler vector ----> look likes rings

第二小特徵值的特徵向量（Fiedler vector），如同套環。

兩種形式：二次最佳化、一次方程組

兩種算法：疊代、求解

兩種元件：邊長、面積

laplacian solver

1. fixed point iterator: explicit / implicit
2. sparse matrix solver

smoothing by minimize x*Lx  (with boundary condition: two endpoint)
smoothing by solve Lx = 0
smoothing by iterate x = x + Lx (fixed point iteration of linear function)

discrete continuous
x*Lx     integral || grad x(u,v) ||^2 dudv
Lx       laplace x

polyline smoothing geometry meaning
                combinatorial laplacian           cotangent laplacian
quadratic form  sum of pairwise squared distance  sum of triangle area
linear form     pit                               proj foot of minimal surface

兩種面積公式：外心、內心

兩種方向：中點、垂足

三維的情況下，L改成L_cot，中點變垂足。

沿著垂直方向移動，更能保留折線原始造型。

x' = x - L̂_cotx;

填坑。將凹陷程度添加回去。

坑底：點
地面：「鄰點的平均數」或者「鄰點的最小曲面的垂足（cotangent Laplacian）」
凹陷程度：坑底到地面的差距

等價於Laplace equation求解，不動點遞推法Richardson iteration。逐步更新點座標，同步地、漸漸地平滑。

外心面積公式

《A Laplacian for Nonmanifold Triangle Meshes》

cotangent laplacian: obuse angle
《Discrete Differential-Geometry Operators for Triangulated 2-Manifolds》
AMixed = 0
For each triangle T from the 1-ring neighborhood of x
If T is non-obtuse,
    // Voronoi safe
    // Add Voronoi formula (see Section 3.3)
    AMixed+ = Voronoi region of x in T
Else
    // Voronoi inappropriate
    // Add either area(T)/4 or area(T)/2
    If the angle of T at x is obtuse
        AMixed+ = area(T)/2
    Else
        AMixed+ = area(T)/4

mesh smoothing (face)🚧

face normal smoothing

面法向量平滑化，再以面法向量推導點座標。

一、法向量平滑化：diffusion filter。

n̕ᵢ = nᵢ + α ∆nᵢ

二、頂點位置最佳化：「邊向量」與「鄰面的法向量」盡量垂直。

min   sum       sum    ‖normal(f) ∙ (pᵢ - pⱼ)‖²
 p  j∈adj(i) f∈face(ij)

二甲、頂點位置公式解（最小平方法）。

《Linear Anisotropic Mesh Filters》

二乙、頂點更新位置（梯度下降法）：「鄰邊向量」投影到「鄰面法向量」，加總之後，作為凹陷程度。

《Fuzzy Vector Median-Based Surface Smoothing》

p̕ᵢ = pᵢ + α  sum       sum    [normal(f) ∙ (pᵢ - pⱼ)] normal(f)
           j∈adj(i) f∈face(ij)

p̕ᵢ = pᵢ + α  sum       sum    proj_normal(f)(pᵢ - pⱼ)
           j∈adj(i) f∈face(ij)

其他領域的平滑化演算法

折線平滑化（座標值XYZ）、訊號平滑化（強度值），演算法互通，可以互相使用。請見本站文件「signal filtering」。

網格平滑化（座標值XYZ）、圖片平滑化（像素值RGB），演算法互通，可以互相使用。請見本站文件「image filtering」。