CG技術の実装と数理： 実装部分の紹介
研究集会CG技術の実装と数理で実装したDecoupling Noises and Features via Weighted l1-analysis Compressed Sensingの技術を簡単にまとめます．
この論文では，ノイズ付のメッシュに対して，なるべくシャープな特徴を残した状態で平滑化することを試みています．特徴とノイズを分けるのは簡単ではありませんが，L1解析を使うことによりシャープな特徴を解析できるというのが基本アイデアになります． 発表資料，概要は，CG技術の実装と数理概要をご参照ください．
提案手法では，2ステップの平滑化を提案していて， - Phase I: 異なるノイズ設定に対して，最適な平滑化メッシュを計算 * 離散Laplacian平滑化モデル * GCVによる最適な平滑化 - Phase II: シャープさを考慮しつつ平滑化 * L1解析によるシャープな特徴の抽出 * シャープさを考慮しつつ平滑化 の2ステップを繰り返しながら，段々とシャープな特徴を増やしていき平滑化を進めていくという流れになります．
Phase I: 異なるノイズ設定に対して，最適な平滑化メッシュを計算
離散Laplacian平滑化モデル
ノイズ付の頂点列$P$を頂点制約に，平滑化後のLaplacian項 $L S$が0になるようにして($L$はLaplacian行列)，平滑化後の頂点$S$を求めます．
$\mathrm{arg min}_S ||P - S||^2 + \lambda || L S ||^2 $
$\lambda$の値を変えて平滑化の度合いを調整できますが，入力ノイズの設定が変わってしまった場合，新しい$\lambda$に調整してやる必要があります．
以下，平滑化部分のPythonコードを抜粋した物です．
def laplacianSmoothing(P, L, lmbd):
M = L.T.dot(L)
A = scipy.sparse.csc_matrix(np.identity(M.shape[0])) + lmbd * M
lin_solve = lin_solve_op(A, P)
P_new = lin_solve(A, P)
return P_new
lin_solve_op(A, P)は，入力行列の種類，サイズに応じて適切なsolverを返す関数です．内部実装では，scipy.sparse.linalg.spsolveをPの各列に適用しているだけです．
GCVによる最適な平滑化
この論文では，$\lambda$の調整処理を自動化するため，Generalized Cross Validation (GCV)という手法を使って最適な$\lambda$を計算しています．参考までですが，この手法はスプライン曲線の平滑化パラメータを自動計算する際によく使われているようです．
詳しい手法は論文を参照していただくとよいと思いますが，GCVでは，ノイズ付頂点列$P$，ラプラシアン行列$L$，$\lambda$を入力とする評価関数$GCV$を定義し，評価関数が最小=最適な平滑化としています．$P$，$L$は固定されているため，最適な$\lambda$を求めてあげれば良いのですが，関数に$I + \lambda L^t L$の逆行列を含んでおり，最小値を式から求めるのは少々困難です．そこで，この論文では，直線探索を使って評価関数$GCV$を最小化する$\lambda$を求めています．
以下，GCVのPythonコードを抜粋した物です．
def gcvSolve(P, L, lmbd_min = 0.05, lmbd_max = 1000.0):
eigenValues = computeEigenValues(L.T.dot(L))
lmbd_opt = scipy.optimize.fminbound( gcvFunc(eigenValues, P, M), lmbd_min, lmbd_max )
return lmbd_opt
今回の実装では，scipy.optimize.fminboundを使って評価関数を最小にする$\lambda$を計算しています．この論文では，高速化のため，$M = L^t L$の固有値を利用して評価関数を近似しているので，gcvFuncにeigenValues(固有値)が含まれる形になっています．
実際にやってみると，下図のような形で最適な平滑化パラメータを求めることができます．
Phase II: シャープさを考慮しつつ平滑化
L1解析によるシャープな特徴の抽出
この論文では，ノイズ付メッシュからシャープな特徴を抽出するために重み付のL1解析を使っています．各頂点の重みは，平滑化メッシュ$S$のLaplacian強度$| L_i S |$から計算され，平滑化メッシュにおいて曲率が高い部分を優先するように調整しています．今回の実装では，Alternating Direction Method of Multipliers (ADMM)の手法を参考にし，L1解析部分を実装しました．
def solve_admm(A, b, W, lmbd = 1.0, rho = 1.0, maxiter = 200, shrinkage, error_check):
mat_prod = mat_prod_op(A)
dot_prod = dot_prod_op(A, b)
m,n = A.shape
l = b.shape[1]
x = numpy.random.normal(size = (n, l)) * 0.001
z = numpy.random.normal(size = (m, l)) * 0.001
a = numpy.random.normal(size = (m, l)) * 0.001
AtA = mat_prod(A.T, A)
Atb = dot_prod(A.T, b)
WtW = mat_prod(W.T, W)
soft_thre = lmbd / rho
A_sub = AtA + rho * WtW
lusolve = lu_solve_op(A_sub, Atb)
for i in range(maxiter):
b_sub = Atb + rho * dot_prod(W.T, z - a)
x = lusolve(b_sub)
if isNormalize:
x = normalizeVectors(x)
z = shrinkage((dot_prod(W, x) + a ), soft_thre)
a = a + dot_prod(W, x) - z
if errorCheck(x-z):
break
return x, z
ただし，実験したところでは，メッシュ解像度の影響を受けてしまい，論文ほどの精度は得られませんでした．
以下は，比較的上手くいった高解像度メッシュ(頂点数9000)の結果ですが，ノイズが強くなってくると，影響を受けて周りの頂点を誤検出します．
シャープさを考慮しつつ平滑化
シャープさを考慮しつつ平滑化するステップでは，解析したシャープな特徴を基にLaplacian行列を修正することでシャープな特徴が保たれるようにします． - コーナー部分: $L(s_i) = 0$ - クリース部分: $L(s_i) = s_j + s_k - 2 s_i$ ($s_j$，$s_k$は$s_i$につながるクリースの端点)
以下，Laplacian行列の修正部分のPythonコードになります．
def updateLaplacianMatrix(L, creaseIDs, cornerIDs, VCList):
Ldense = L.todense()
w_crease = 10.0
sharpIDs = set(creaseIDs) | set(cornerIDs)
for p in creaseIDs:
Ldense[p,:] = 0.0
w_sum = 0.0
for q in VCList[p]:
if q in sharpIDs:
Ldense[p, q] = -w_crease
w_sum += w_crease
Ldense[p, p] = w_sum
for p in cornerIDs:
Ldense[p,:] = 0.0
L_new = scipy.sparse.csc_matrix(Ldense)
return L_new
論文中では，クリースについても同じ重みを使っていましたが，あまり良い結果が得られなかったので，クリース部分のつながりを優先するようにw_creaseの重みを調整しています．
今回の実装では，L1解析部分は上手くいきませんでしたが，シャープな特徴が事前に分かっている場合には，良好な結果が得られました．
参考文献
[1] Decoupling Noises and Features via Weighted l1-analysis Compressed Sensing: http://staff.ustc.edu.cn/~lgliu/Projects/2014_DecouplingNoise/default.htm [2] MATLAB scripts for alternating direction method of multipliers: https://web.stanford.edu/~boyd/papers/admm/