As-Rigid-As Possible Surface Modeling

son·2023년 1월 16일
0

ARAP은 메쉬 변형 알고리즘이다. CGAL에도 구현되어있고, python으로도 Open3D 통해서 사용할 수 있다. Vertex 수천개 정도 스케일에서는 꽤 빠르게 돌아가는 편이다.

ARAP으로 변형한 Armadillo. 오른쪽 그림에서 오른쪽 다리의 빨간 점들이 fixed point이고 왼손의 노란 점이 handle이다.

배경

두 surface 사이의 차이는 shell energy라는 것을 통해 측정할 수 있는데

Es(S,S)=ΩksIIF2+kbIIIIF2E_s(\mathcal{S},\mathcal{S}')=\int_{\Omega}k_s\|\mathbf{I}'-\mathbf{I}\|^2_F+k_b\|\mathbf{II}'-\mathbf{II}\|^2_F

I\mathbf{I}II\mathbf{II}는 각각 first, second fundamental form이다. I\mathbf{I}은 surface의 strechting과 관련 있고, II\mathbf{II}는 bending과 관련 있다. 따라서 만약 surface S\mathcal{S}S\mathcal{S}'사이의 deform이 rigid하면 shell energy는 0이다. Non-rigid 한 deform을 찾을 때에도, shell energy가 0은 아니더라도 최소화 되도록 찾으면 (as-rigid-as possible 하게), local한 shape이 유지될 수 있다.

Cell Ci\mathcal{C}_i (여기선 vertex ii의 one-ring neighborhood로 정의)와 deform된 버전 Ci\mathcal{C}_i'에 대해 deformation이 rigid하다면 rotation matrix Ri\mathbf{R}_i

pipj=Ri(pipj),jN(i)\mathbf{p}'_i-\mathbf{p}'_j=\mathbf{R}_i(\mathbf{p}_i-\mathbf{p}_j), \forall j\in \mathcal{N}(i)

로 나타낼 수 있다. pi\mathbf{p}_i는 vertex 좌표. N(i)\mathcal{N}(i)는 vertex ii의 neighborhood (참고로,
jN(i)\forall j\in \mathcal{N}(i)에 대해 iN(j)i\in\mathcal{N}(j)이므로 cell들은 overlapping되어 있다.). Deformation이 non-rigid면 deform을 제일 잘 approximation하는 Ri\mathbf{R}_i는 아래의 weighted square를 minimize하는 것이다.

E(Ci,Ci)=jN(i)wij(pipj)Ri(pipj)2E(\mathcal C_i,\mathcal C'_i) = \sum_{j\in\mathcal N(i)}w_{ij}\|( \mathbf{p}'_i-\mathbf{p}'_j)-\mathbf{R}_i(\mathbf{p}_i-\mathbf{p}_j)\|^2

pipj\mathbf{p}_i-\mathbf{p}_jeij\mathbf{e}_{ij}로 적자. 그럼 위 식은 다시 쓰면

jwij(eijRieij)T(eijRieij)=jwij(eijTeij2eijTRieij+eijTeij)\sum_j w_{ij}(\mathbf{e}'_{ij}-\mathbf{R}_i\mathbf{e}_{ij})^T(\mathbf{e}'_{ij}-\mathbf{R}_i\mathbf{e}_{ij})=\sum_j w_{ij}(\mathbf{e}'^T_{ij}\mathbf{e}'_{ij}-2\mathbf{e}'^T_{ij}\mathbf{R}_i\mathbf{e}_{ij}+\mathbf{e}^T_{ij}\mathbf{e}_{ij})

나머지 항은 상수니까 Ri\mathbf{R}_i을 포함하는 term만 minimize할 수 있는데,

arg minRij2wijeijTRieij=arg maxRijwijeijTRieij\argmin_{\mathbf{R}_{i}} \sum_j -2w_{ij}\mathbf{e}'^T_{ij}\mathbf{R}_i\mathbf{e}_{ij} = \argmax_{\mathbf{R}_{i}} \sum_j w_{ij}\mathbf{e}'^T_{ij}\mathbf{R}_i\mathbf{e}_{ij}

Weighted inner product를 trace를 사용한 형태로 변형시키면 (Tr(yxT)=xTyTr(\mathbf{y}\mathbf{x}^T)=\mathbf{x}^T\mathbf{y}),

=arg maxRiTr(jwijRieijeijT)=arg maxRiTr(RijwijeijeijT)=\argmax_{\mathbf{R}_{i}} Tr\left(\sum_j w_{ij}\mathbf{R}_i\mathbf{e}_{ij}\mathbf{e}'^T_{ij}\right)=\argmax_{\mathbf{R}_{i}} Tr\left(\mathbf{R}_i\sum_j w_{ij}\mathbf{e}_{ij}\mathbf{e}'^T_{ij}\right)

이다.

Si=jN(i)wijeijeijT=PiDiPiT\mathbf{S}_i=\sum_{j\in\mathcal{N}(i)} w_{ij}\mathbf{e}_{ij}\mathbf{e}'^T_{ij}=\mathbf{P}_i\mathbf{D}_i\mathbf{P}'^T_i

Si\mathbf{S}_ieij\mathbf{e}_{ij}로 구성된 3×N(i)3\times|\mathcal{N}(i)| 행렬 Pi\mathbf{P}_iwijw_{ij}로 구성된 N(i)×N(i)|\mathcal{N}(i)|\times|\mathcal{N}(i)| 대각행렬 Di\mathbf{D}_i를 이용해서 쓸 수 있다. 만약 행렬 M\mathbf{M}이 PSD이면 어느 orthogonal 행렬 R\mathbf{R}에 대해서도 Tr(M)Tr(RM)Tr(\mathbf{M})\geq Tr(\mathbf{RM})이므로, RiSi\mathbf{R}_i \mathbf{S}_i를 symmetric PSD로 만드는 행렬 Ri\mathbf{R}_iTr(RiSi)Tr(\mathbf{R}_i\mathbf{S}_i)를 maximize한다. 그러한 Ri\mathbf{R}_i 는 SVD해서 (Si=UiΣiViT\mathbf{S}_i = \mathbf{U}_i \mathbf{\Sigma}_i\mathbf{V}^T_i)

Ri=ViUiT\mathbf{R}_i=\mathbf{V}_i \mathbf{U}^T_i

를 통해서 구할 수 있다.

구체적인 메소드

앞에서 정의했던 cell에 대한 에너지를 전체 surface에 대해 합하면

E(S)=i=1nwiE(Ci,Ci)=iwijN(i)wij(pipj)Ri(pipj)2E(\mathcal{S}')=\sum^n_{i=1}w_i E(\mathcal{C}_i, \mathcal{C}'_i)=\sum_i w_i\sum_{j\in\mathcal N(i)}w_{ij}\|( \mathbf{p}'_i-\mathbf{p}'_j)-\mathbf{R}_i(\mathbf{p}_i-\mathbf{p}_j)\|^2

이 때, weight wiw_iwijw_{ij}를 어떻게 골라야 할까? Deformation energy가 최대한 mesh-independent 하기 위해서 (특히 이 경우 non-uniform한 cell의 영향을 줄이기 위해서) cotangent weight을 사용한다.

wij=(cotαij+cotβij)/2w_{ij}=(\cot \alpha_{ij}+\cot \beta_{ij})/2

α\alphaβ\betae\mathbf{e}에 마주보고 있는 angle이다. Computing Discrete Minimal Surfaces and Their Conjugates을 참조하면 Dirichlet energy란 함수가 얼마만큼 변화하는가에 대한 measure이다(ED(f)=12Ωf2E_{D}(f)=\frac{1}{2}\int_{\Omega}|\nabla f|^2). 이런 surface에서 angle-preserving 한 mapping을 구하려면, 얼마만큼 Cauchy-Riemann equation을 만족하느냐를 나타내는 conformal energy라는 것을 사용하면 되는데, 이 conformal energy를 minimize하는 것이 Dirichlet energy 빼기 mapping area를 minimize하는 것과 동치라고 한다(https://math.stackexchange.com/questions/3213598/what-are-the-use-cases-of-the-dirichlet-energy-in-computer-vision 참조). 아무튼 계산해보면 triangle간의 linear map에 대한 Dirichlet energy가

ED(f)=14i=1edges(cotαi+cotβi)ei2E_D(f)=\frac{1}{4}\sum_{i=1}^{edges}(\cot\alpha_i + \cot\beta_i)|\mathbf{e}_i|^2

가 된다 한다. 이 내용을 ARAP에서 이용한 것이다. Cell energy가 cell area에 비례하는데 Voronoi area를 이용해서 normalize 해줄 수 있다.

우리가 ARAP을 통해 하고 싶은 것은 sparse한 user input을 constraint으로 두고, 나머지에 해당하는 vertex의 위치를 적절히 결정하고 싶은 것이다. 따라서 우리는 문제를 p\mathbf{p}’에 대해 풀어야한다. 또 R\mathbf{R} 역시 찾아야한다. 그래서 두 step으로 나눠서 한번은 R\mathbf{R}을 고정하고 p\mathbf p’를 찾고, 한번은 p\mathbf p’를 고정하고 R\mathbf R을 찾는다.

먼저 각 i에 대해 Ri\mathbf{R}_i을 찾는 것은 서로 독립적이므로 앞서 언급한 Si\mathbf{S}_i의 SVD를 이용한 방식으로 찾을 수 있다. (참고로, ARAP 외에 다른 deform 방식 중에는 이웃한 Ri\mathbf{R}_i 이 서로 비슷한 값을 갖도록 유도하는 것들도 있다.)

pi\mathbf{p}'_i는 partial derivate를 이용해서 찾을 수 있다.

E(S)pi=pi(jN(i)wij(pipj)Ri(pipj)2+jN(i)wji(pjpi)Rj(pjpi)2)=jN(i)2wij((pipj)Ri(pipj))jN(i)2wji((pjpi)Rj(pjpi))\frac{\partial E(\mathcal{S}')}{\partial\mathbf{p}'_i}=\frac{\partial }{\partial\mathbf{p}'_i}\left( \sum_{j\in\mathcal{N}(i)} w_{ij}\|( \mathbf{p}'_i-\mathbf{p}'_j)-\mathbf{R}_i(\mathbf{p}_i-\mathbf{p}_j)\|^2 + \sum_{j\in\mathcal{N}(i)} w_{ji}\|( \mathbf{p}'_j-\mathbf{p}'_i)-\mathbf{R}_j(\mathbf{p}_j-\mathbf{p}_i)\|^2\right) \\ = \sum_{j\in\mathcal{N}(i)} 2w_{ij}\left(( \mathbf{p}'_i-\mathbf{p}'_j)-\mathbf{R}_i(\mathbf{p}_i-\mathbf{p}_j)\right) -\sum_{j\in\mathcal{N}(i)} 2w_{ji}\left(( \mathbf{p}'_j-\mathbf{p}'_i)-\mathbf{R}_j(\mathbf{p}_j-\mathbf{p}_i)\right)

이 때, wij=wjiw_{ij}=w_{ji} 이므로

E(S)pi=jN(i)4wij((pipj)12(Ri+Rj)(pipj))\frac{\partial E(\mathcal{S}')}{\partial\mathbf{p}'_i} = \sum_{j\in\mathcal{N}(i)} 4w_{ij}\left(( \mathbf{p}'_i-\mathbf{p}'_j)-\frac{1}{2}(\mathbf{R}_i+\mathbf{R}_j)(\mathbf{p}_i-\mathbf{p}_j)\right)

derivative가 0이 되려면

jN(i)wij(pipj)=jN(i)wij2(Ri+Rj)(pipj)\sum_{j\in\mathcal{N}(i)} w_{ij}( \mathbf{p}'_i-\mathbf{p}'_j)=\sum_{j\in\mathcal{N}(i)} \frac{w_{ij}}{2}(\mathbf{R}_i+\mathbf{R}_j)(\mathbf{p}_i-\mathbf{p}_j)

왼쪽 term은 discrete Laplace-Beltrami operator가 p\mathbf p'에 적용된 형태임을 볼 수 있다. 따라서 Lp=b\mathbf{Lp}'=\mathbf{b} 꼴이고 Cholesky decomposition을 이용해서 빠르게 풀 수 있다. (https://en.wikipedia.org/wiki/Cholesky_decomposition 의 application 항목 첫번째 문단 참조)

앞서 이웃한 Ri\mathbf{R}_i 이 서로 비슷한 값을 갖도록 유도하는 것들도 있다고 했는데, 이 term을 ARAP에 더한 것이 Smoothed Rotation Enhanced As-Rigid-As Possible (SR_ARAP) Deformation 이고 cell energy term에 RiRjF2\|\mathbf{R}_i-\mathbf{R}j\|^2_{F}을 더한다.

0개의 댓글