물리 기반 시뮬레이션/기초

[ 물리 기반 시뮬레이션 ] 03. 위치 기반 역학 (Position Based Dynamics)

coding-l7 2023. 7. 3. 18:01

📌 위치 기반 시뮬레이션 (Position Based Simulation)

물리 기반 시뮬레이션 분야는 강체의 역학, 유체의 흐름과 같은 물리 현상의 시뮬레이션을 위한 새로운 방법을 찾는 데에 관심이 있고, 컴퓨터 그래픽스에서는 속도와 안정성이 중요하다.

하지만, 기존의 힘을 속도와 가속도를 적분하는 방식은 계산이 복잡하다.

➡️ 따라서, 위치에 직접적으로 작용하여 위치를 제어하는 위치 기반 시뮬레이션이 등장했다.

 

위치 기반 시뮬레이션의 주요 기능과 장점은 다음과 같다:

✔️ 불안정성 제거 (명시적 적분 제어)

✔️ 정점의 위치를 직접 조작 가능

✔️ 일반적인 제약 조건을 처리

✔️ 간단한 구현 (명시적 위치 기반 Solver)

📌 알고리즘 (Algorithm)

먼저, 알고리즘을 이야기하기 전에 동적 객체를 표현하기 위해 다음과 같은 조건을 고려해야 한다.

  • 동적 객체를 나타낼 때 $N$ 개의 정점과 $M$ 개의 제약 조건으로 나타낸다.
  • 각각의 정점 $i \in [1,...,N]$ 는 질량 $m_i$, 위치 $x_i$ , 속도 $v_i$ 를 가진다.
  • 각각의 제약조건 $j \in [1,...,M]$  는 다음과 같이 구성되어 있다.
    • 카디널리티 $n_j$
    • 함수 $C_j$
    • 강성 매개변수 $k_j \in [0,1]$
    • $\text{equality}$ 또는 $\text{inequality}$ 유형 중 하나

이때, 유형이 $\text{equality}$인 제약 조건 $j$ 는 $C_j = 0$일 때 만족하며, $\text{inequality}$의 경우 $C_j \geq 0$ 일 때 만족으로 간주한다.

이는 시간 단계 $\Delta t$ 에 기반하여 동적 객체를 다음과 같이 시뮬레이션한다:

🧩 의사 코드 (pseudo-code)

1) 	모든 정점 i에 대해 반복
2) 	  xi = x0i, vi = v0i, wi = 1/mi //초기화
3) 	반복 종료
4) 	시뮬레이션 동안 반복
5)		모든 정점 i에 대해 반복 
6)			vi ← vi +∆twifext(xi) //외부 힘 계산 ex) 중력
7)		dampVelocities(v1,...,vN) //감쇠
8)		모든 정점 i에 대해 반복 
9)			pi ← xi +∆tvi 
10)		모든 정점 i에 대해 반복 
11)			generateCollisionConstraints(xi → pi) //충돌 제약 조건 생성
12)		solverIterations 횟수만큼 반복
13)			projectConstraints(C1,...,CM+Mcoll ,p1,...,pN) //제약 조건 투영
14)		반복 종료
15)		모든 정점 i에 대해 반복
16)			vi ← (pi −xi)/∆t //각 입자의 속도와 위치를 예측값 p를 이용해 업데이트
17)			xi ← pi
18)		endfor
19)		velocityUpdate(v1,...,vN)
20)	반복 종료

각 줄의 설명은 다음과 같다:

  • 1-3 : 각 정점의 위치, 속도, 역질량 값을 초기화
  • 5-6 : 각 정점을 외부 힘에 대하여 속도를 계산 (여기에서는 $v_i = v_i + \Delta tg$ 로 중력을 표현)
  • 7 : 전역 감쇠를 추가 (아래의 내용 중 감쇠 참고)
  • 8-9 : 각 정점의 새 위치에 대한 예측값 $p_i$ 를 오일러 방법을 이용해 계산
  • 10-11 : 시간 단계마다 변하는 충돌 제약 조건을 생성
  • 12-14 : 위에서 계산한 예측값 $p_i$ 가 제약 조건을 만족하도록 반복적으로 조정, 고정된 제약 조건과 충돌 제약 조건을 모두 고려한다. (Gauss-Seidel 방법과 Newton-Raphson 방법) (아래의 내용 중 $\text{solver}$ 참고)
  • 15-18 : 정점의 위치를 최적화된 예측값으로 이동, 이를 바탕으로 속도를 계산 (Verlet 적분)
  • 19-20 : 충돌하는 정점의 속도를 마찰 및 반발 계수에 따라 조정

이때, $\text{iterator}$를 이용한 반복에서 반복 횟수를 증가시키면, 충돌 감지에서 $\text{solver}$로 병목 현상이 이동한다.

📌 Solver

$\text{solver}$의 입력은 $M+M_{coll}$ 개의 제약 조건과 $p_1,...,p_N$ 의 새로운 위치 예측값이다. 이러한 $\text{solver}$는 $p$ 값을 수정하여 모든 제약 조건을 만족하도록 시도하며, 이는 비선형 방정식으로 표현된다.

 

예를 들어, 간단한 거리 제약 조건의 경우 다음과 같이 나타낼 수 있다:

제약 조건을 만족시키기 위해 Gaussian-Seidel 반복법을 사용하며, 각 제약 조건을 독립적으로 한 번에 하나씩 해결한다. 그러나 GS와 달리 제약 조건을 해결하기 위해 비선형 연산이 요구된다. 또한, 단순히 제약 조건들을 순회하며 입자를 제약 조건이 유효한 위치로 투영하고, 위치에 대한 수정을 즉시 반영한다.

📌 제약 조건 투영 (Constraints Projection)

제약 조건 투영은 각 정점에 대한 새로운 위치 값 $p$ 를 변화시켜 제약 조건을 만족시키는 과정이다.

시뮬레이션 동안 반복되는 과정에서 점의 이동과 관련된 중요한 문제로는 선형 운동량과 각 운동량의 보존이 있다.

 

선형 운동량과 각 운동량은 각각 다음과 같이 표현된다:

선형 운동량 P, 각 운동량 L

 

다음과 같은 식을 만족할 경우 선형 운동량과 각 운동량이 보존된다:

선형 운동량 보존
각 운동량 보존

이때, $\Delta p_i$ 는 점 $i$ 의 이동량을 나타내며 $r_i$ 는 $p_i$ 와 임의의 공통 회전 중심 사이의 거리이다.

또한, 선형 운동량 보존의 경우 질량 중심을 보존하는 것을 의미한다.

 

운동량 보존은 주로 객체의 구조나 모양을 유지하는 데에 필요한 내부 제약 조건에서 만족되어야 하며, 투영이 위의 조건을 하나라도 만족하지 못할 경우, 외부 힘과 같은 역할을 하는 유령 힘이 발생한다.

 

내부 제약 조건의 경우, 위에서의 운동량 보존에 의해 $C$ 는 이동 및 회전과 관련이 없음을 의미한다. 즉, 점들의 회전이나 이동은 제약 함수의 값에 영향을 미치지 않는다. 따라서, 제약 조건의 최대 변화량을 나타내는 기울기 $\nabla _pC$ 는 물체의 이동 및 회전을 방향을 의미하는 강체 모드에 수직이다.

 

따라서, $\Delta p$ 를 $\nabla_pC$ 의 방향으로 선택하면 모든 질량이 동일한 경우 선형 운동량과 각운동량이 자동으로 보존된다. 그러므로, 현재 위치 $p$ 가 주어졌을 때, 제약 조건을 만족하는 $\Delta p$  즉, $C(p+\Delta p) = 0$ 을 만족하는 보정 $\Delta p$ 를 찾아야 한다. 이 식은 다음과 같이 근사할 수 있다:

3)

또한, $\Delta p$ 를 $\nabla_pC$ 의 방향으로 제한하기 위해 스칼라 $\lambda$를 선택하여 다음과 같이 나타낸다:

4)

4) 식을 위의 3) 근사식에 대입하여 스칼라 $\lambda$ 값에 대해 계산한 뒤, 4)식에 대입하면 다음과 같은 식을 얻을 수 있다:

마지막의 $\Delta p$ 와 관련된 식은 비선형 방정식의 해를 구하는 방법인 Newton-Raphson 방법의 형태를 가진다.

또한, 이를 각 정점 $p_i$ 에 대한 식으로 나타내면 다음과 같이 나타낼 수 있다:

이때, 스케일링 팩터 $s$ 는 다음과 같다:

만약, 각 정점이 개별적으로 질량을 가진다면, $w_i = 1/m_i$로 가중치를 적용하며, 식은 다음과 같다:

📐 거리 제약 조건 (Distance Constraints)

위와 같은 거리 제약 조건에 대한 예를 들 수 있다. 먼저, 각각의 점에 대한 도함수는 다음과 같다.

이때 $n$ 은 점 $p_2$ 에서 $p_1$ 방향을 가리키는 크기 1의 벡터이다. 또한 이를 통해 스케일링 팩터 $s$ 를 구할 수 있다.

최종적으로, $p_1$ 과 $p_2$ 의 변화율은 다음과 같다.

위에서는 제약의 유형($\text{equality, inequality}$) 과 강성 매개변수 $k$ 를 고려하지 않았다.

먼저, 제약의 유형이 $\text{equality}$ 인 경우에 항상 투영을 수행하며, $\text{inequality}$ 이 경우에는 제약 함수 $C$ 가 $C < 0$ 일 경우에 투영을 수행한다.

 

또한, 강성 매개변수 $k$ 를 적용하는 방법으로는 여러 가지가 있다.

먼저, 가장 간단한 방법으로는 보정값 $\Delta p$ 에 $k \in [0, 1]$ 을 곱하는 것이다. 하지만, $\text{solver}$ 의 반복 루프에서 $k$ 의 효과가 비선형적으로 나타날 수 있다. 예를 들어 단일 거리 제약 조건에서 $n_s$ 만큼의 $\text{solver}$ 를 반복한다면, 오차가 $\Delta p(1-k)^{n_s}$ 가 된다.

따라서, 선형적인 관계를 얻기 위해 보정값을 직접적으로 $k$ 가 아닌, $k' = 1-{(1-k)}^{1/n_s}$를 곱한다. 이를 이용하면 오차는 $\Delta p(1-k')^{n_s} = \Delta p(1-k)$ 가 되어, $k$ 에 선형적으로 의존하고 $n_s$ 에 독립적인 선형 관계를 갖는다.

📌 충돌 감지 및 반응 (Collision Detection & Response)

위치 기반 접근 방식의 장점 중 하나는 충돌 반응이 간단하게 구현된다는 점이다.

위의 알고리즘에서 $M_{coll}$ 개의 충돌 제약 조건이 생성되는데, 이는 각 시간 단계마다 처음부터 생성된다.

충돌 제약 조건의 개수 $M_{coll}$ 은 충돌하는 점의 수에 따라 달라지며, 연속적인 충돌 및 정적인 충돌 모두 처리할 수 있다.

연속적인 충돌 처리

연속적인 충돌 처리를 위해 각 정점 $i$ 에 대해 $x_i \rightarrow p_i$ 방향의 벡터를 테스트한다. 이 벡터가 객체 내부로 들어가면 해당 위치에서 진입점 $q_c$ 와 표면 법선 $n_c$ 을 계산한다. 또한, $\text{inequailty}$ 제약 조건 $C(p) = (p-q_c) \cdot n_c$ 와 강성 $k = 1$ 을 제약 조건 목록에 추가한다.

정적인 충돌 처리

만약, 벡터 $x_i \rightarrow p_i$ 가 객체 내부에 완전히 위치한다면, 연속적인 충돌 감지는 실패한 것이며, 이 경우 정적인 충돌 처리로 넘어간다.

정적인 충돌 처리를 위해 $p_i$ 에 가장 가까운 표면 점 $q_s$ 와 해당 위치에서의 표면 법선 $n_s$ 를 계산한다. 또한, 부등식 제약 조건 $C(p) = (p-q_s) \cdot n_s$ 와 강성 $k = 1$ 을 제약 조건 목록에 추가한다.

 

➡️ 충돌 제약 조건 생성은 $\text{solver}$ 루프 밖에서 수행되며, 이는 시뮬레이션의 속도를 빠르게 해준다. 

 

위에서 이야기한 충돌 처리는 정적 객체와의 충돌에 한정된다. 즉, 충돌 상대에게는 충격이 전달되지 않는다. 따라서, 두 개의 동적인 객체의 충돌 처리는 두 개의 객체를 동시에 시뮬레이션해야 한다.

 

예를 들어, 한 객체의 한 점 $q$ 가 다른 객체의 삼각형 $p_1, p_2, p_3$ 을 통과한다면, 다음과 같은 제약 함수로 표현되는 부등식 제약 조건을 삽입할 수 있다.

이는 점 $q$ 를 삼각형의 올바른 면(면에 대하여 원래 있던 방향)에 위치시키며, 이 제약 함수는 강체 모드와 독립적이기 때문에 선형 운동량과 각 운동량을 보존한다.

📌 감쇠 (Damping)

위의 알고리즘에서 언급된 내용으로, 감쇠는 새로운 위치를 계산하기 전에 속도를 감쇠 시킨다.

이에 대한 과정은 다음과 같은 의사 코드로 나타낼 수 있다:

(1) xcm = (∑i ximi)/(∑i mi)
(2) vcm = (∑i vimi)/(∑i mi)
(3) L = ∑i ri ×(mivi)
(4) I = ∑i r˜ir˜Ti mi
(5) ω = I^(−1)L
(6) 모든 정점 i에 대해서 반복
(7) 	∆vi = vcm + ω × ri − vi
(8) 	vi ← vi +kdamping∆vi
(9) 반복 종료

이때, $r_i = x_i - x_{cm}$ 이고, $ \tilde{r}_i$ 는 $\tilde{r}_iv = r_i \times v$ 성질을 갖는 3x3 행렬이다. 또한, $k_{damping}$ 은 감쇠 계수로 $[0, 1]$ 범위를 가진다.

  • 1-5 줄에서는 global linear 속도 $v_{cm}$과 각속도 $ω$ 를 계산한다.
  • 6-9 줄에서는 각 정점에 대해 개별 속도 $v_i$ 의 변화 $\Delta v_i$ 를 전역 운동 $v_cm + ω × r_i$ 와의 차이로 감쇠시킨다. ➡️ 따라서 $k_{damping} = 1$ 인 경우 전역 운동만 유지되며, 정점의 집합은 강체처럼 작동한다.

📌 연결 (Attachments)

위치 기반 접근 방식에서 정적 또는 동적인 객체에 정점을 연결하는 것은 매우 간단하다. 정점의 위치는 단순히 정적 대상 위치로 설정되거나 동적 객체의 위치와 일치하도록 각 시간 단계마다 업데이트된다. 이 정점을 포함하는 다른 제약 조건들이 이를 움직이지 않도록 하기 위해 역질량 $w_i$ 를 0으로, 즉, 정점의 질량을 무한대로 설정한다.

참고 자료 : 

Position Based Dynamics - Matthias Müller, Bruno Heidelberger, Marcus Hennix, John Ratcliff

선형 운동량 : 

https://ko.wikipedia.org/wiki/운동량

각 운동량 : 

https://ko.wikipedia.org/wiki/각운동량

가우스-자이델(Gauss-Seidel) 법 :

https://study2give.tistory.com/entry/가우스자이델

뉴턴-랩슨(Newton-Raphson) 법 : 

https://evening-ds.tistory.com/38