이 글에서, 나는 널리 사용되는 공식의 유도 뒤에 있는 사고 프로세스를 설명하려고 한다: 그것은 두 개의 3D 벡터 사이의 회전을 나타내는 쿼터니언을 찾는 것이다. 정말 새로운 것은 없지만, 희망스럽게도, 몇 가지 아이디어들이 다른 때에 재사용될 수 있다.
Note : 여기에서 보여지는 과정은 목적에서 불완전하다. production code에서 사용될 수 있는 버전에 대해, 대신에 다음의 자료를 보아라, Quaternion from two vectors: the final version.
Naive method
회전은 회전축과 한 각을 사용하여 가장 잘 시각화된다. degenerate cases를 제외하고서, 그 회전축은 두 개의 원래 벡터들의 외적을 연산하여 얻어질 수 있다:
그러고나서, 그 각도는 외적 그리고/또는 내적의 성질들을 사용하여 얻어질 수 있다:
외적과 내적공식.. 외적은 cosTheta, 내적은 sinTheta
θ는 항상 0과 π사이이기 떄문에(저 기호가 뭔지모르겠다, 여튼 각도 제한이 있다), 우리는 오직 내적만을 신경쓴다. 이것은 우리에게 쿼터니언을 만드는 어떤 명백한 코드를 준다 (명확성을 위해 θ = 0인 corner cases를 생략하고):
quat quat::fromtwovectors(vec3 u, vec3 v) { float cos_theta = dot(normalize(u), normalize(v)); float angle = acos(cos_theta); vec3 w = normalize(cross(u,v)); return quat::fromaxisangle(angle, w); }
이것은 단순하지만, 작동한다. "quaternion from two vectors"에 대해 구글링 하는 것은 이 포럼 포스트와 이 질문을 보여주는데, 여기에서 변형 이론이 제안된다.
Looking under the hood
axis와 angle로부터 쿼터니언을 구성할 때 무슨 일이 발생하는지를 봐보자:
이것은 quat::fromaxisangle에 대한 코드가 어느정도 이렇게 보인다는 것을 의미한다:
quat quat::fromaxisangle(float angle, vec3 axis) { float half_sin = sin(0.5f * angle); float half_cos = cos(0.5f * angle); return quat(half_cos, half_sin * axis.x, half_sin * axis.y, half_sin * axis.z); }
Avoiding trigonometry
만약 너가 Inigo Quilze의 최근 자료 avoiding trigonometry를 읽었다면, 너는 아마도 우리가 cos(θ)로부터 θ를 연산하고, 그러고나서 sin(θ/2)와 cos(θ/2)를 연산했다는 사실에 눈살을 찌뿌렸을 것이다.
사실, 그것을 하는 좀 더 간단한 방법이 있다; precalculus의 반각공식이 우리에게 다음을 말해준다:
이것은 우리가 우리의 쿼터니언을 만드는 코드를 간단하게 만든다:
quat quat::fromtwovectors(vec3 u, vec3 v) { float cos_theta = dot(normalize(u), normalize(v)); float half_cos = sqrt(0.5f * (1.f + cos_theta)); float half_sin = sqrt(0.5f * (1.f - cos_theta)); vec3 w = normalize(cross(u,v)); return quat(half_cos, half_sin * w.x, half_sin * w.y, half_sin * w.z); }
이것은 꽤 좋다. 잘 알려진 삼각법 공식을 사용하여, 우리는모든 삼각기하 함수 호출을 제거했다!
Avoiding square roots
우리가 조금 더 잘할 수 있다는 것이 발생한다. 우리가 세 개의 벡터들을 표준화시킨다는 것에 주목해라: u, v 그리고 cross(u,v). 그것은 세 개의 square roots이다. 요지는, 우리가 이미 이 공식을 통해서 w의 norm을 안다는 것이다:
글고 우리는 또 다시 precalculus로부터 sin(θ)를 안다:
또한, sqrt(a) * sqrt(b) = sqrt(ab)라는 사실을 이용하여, 우리는 더 적은 square root를 수행한다.
그러므로 우리는 다음의 성능향상 개선을 만들어낼 수 있다:
glm::quat fromtwovectors(const glm::vec3 &u, const glm::vec3& v) { float norm_u_norm_v = glm::sqrt(glm::dot(u,u) * glm::dot(v, v)); float cos_theta = glm::dot(u, v) / norm_u_norm_v; float half_cos = sqrt(0.5f * (1.f + cos_theta)); float half_sin = sqrt(0.5f * (1.f - cos_theta)); glm::vec3 w = cross(u, v) * (1.f / (norm_u_norm_v * 2.f * half_sin * half_cos); return quat(half_cos, half_sin * w.x, half_sin * w.y, half_sin * w.z); }
오 잠깐! 우리는 w를 연산하기 위해 sin(θ/2)로 나누고{vec3 w를 계산할 때 }, 그러고나서 또 다시 sin(θ/2)를 곱한다{return quat할 때}. 이것은 우리가 그 변수가 필요없다는 것을 의미하고, 우리는 더 간단하게 할 수 있다:
glm::quat fromtwovectors(const glm::vec3 &u, const glm::vec3& v) { float norm_u_norm_v = glm::sqrt(glm::dot(u,u) * glm::dot(v, v)); float cos_theta = glm::dot(u, v) / norm_u_norm_v; float half_cos = sqrt(0.5f * (1.f + cos_theta)); glm::vec3 w = cross(u, v) * (1.f / (norm_u_norm_v * 2.f * half_cos); return quat(half_cos, w.x, w.y, w.z); }
이것은 OgreVector3.h에서 Ogre3D engine에 의해서 사용되는 코드이다. 그것들이 최종 결과에 대해 부가적인 표준화 단계를 수행하는 것을 제외하고. 이것은 수학적으로 쓸모없지만, numerical 안정성 문제 때문에, 그럼에도 불구하고 그렇게 하는 것은 아마도 안전하다.
이 마지막 표준화 단계는 실제로 그 코드를 간결하게 하는 한 기회이다.
Improving on Ogre3D
우리는 두 개의 square roots와 4개의 나눗셈으로 내려왔고, 더해서 꽤 몇 개의 곱셈과 더하기가 있다. 우리가 작동시키는 플랫폼에 따라서, 더 간단하게 하고 성능을 개선시키는 것이 가능하다. 예를들어, 많은 SIMD 아키텍쳐에서, 쿼터니언을 표준화 시키는 것은 매우 빠를 수 있다.
이것은 만약 우리가 그 쿼터니언에 2.f * half_cos로 모든 컴포넌트에 곱한다면 얻는 코드이고, normalize()가 그 나머지 일을 하게 한다:
glm::quat fromtwovectors(const glm::vec3 &u, const glm::vec3& v) { float norm_u_norm_v = glm::sqrt(glm::dot(u,u) * glm::dot(v, v)); float cos_theta = glm::dot(u, v) / norm_u_norm_v; float half_cos = sqrt(0.5f * (1.f + cos_theta)); glm::vec3 w = cross(u, v) * (1.f / norm_u_norm_v); return normalize(quat(2 * half_cos * half_cos, w.x,
w.y, w.z)); }
이제 half_cos는 오직그것의 제곱 형태로만 나타나고, 그것은 square root로부터 나오기 때문에, 우리는 간단히 그 square root를 제거할 수 있다:
glm::quat fromtwovectors(const glm::vec3 &u, const glm::vec3& v) { float norm_u_norm_v = glm::sqrt(glm::dot(u,u) * glm::dot(v, v)); float cos_theta = glm::dot(u, v) / norm_u_norm_v; glm::vec3 w = cross(u, v) * (1.f / (norm_u_norm_v); return normalize(quat(1.f + cos_theta, w.x, w.y, w.z)); }
1.f + cos_theta는 아까 반각공식에서 제곱하고 * 2하면 저렇게 나와서 코드가 저렇게 된다.
}
그리고 같은 이유로, 우리는 모든 쿼터니언 컴포넌트를 norm_u_norm_v로 곱할 수 있다:
glm::quat fromtwovectors(const glm::vec3 &u, const glm::vec3& v) { float norm_u_norm_v = glm::sqrt(glm::dot(u,u) * glm::dot(v, v)); glm::vec3 w = cross(u, v); return normalize(quat(norm_u_norm_v + dot(u, v), w.x, w.y, w.z)); }
우리는 여전히 두 개의 square roots와 4개의 나누기 연산을 하는데, 거기에서 어떤 것들은 normalize()에 숨겨져 있지만, 그 코드는 이제 상당히 더 짧다.
Final form
만약 u와 v가 단위 벡터라고 강요될 수 있다면, norm_u_norm_v는 생략되어질 수 있고, 간단히 1.0f로 대체될 수 있다:
glm::quat fromtwovectors(const glm::vec3 &u, const glm::vec3& v) { glm::vec3 w = cross(u, v); return normalize(quat(1.f + dot(u, v), w.x, w.y, w.z)); }
아름답지 않냐? 우리가 시작한 sin(), cos() 그리고 acos() 어지러운 것들을 고려하면?
이 알고리즘은 인터넷 전반에서 발견될 수 있지만, 나는 누가 처음에 이것을 만들어냈는지를 모른다. 또한, 많은 3D 엔진들 (publicly하게 이용가능하고 다소 private하기도 한) 은 그것으로부터 이득을 얻는다.
Update(06/01/2014)
아래의 코멘트에서, Michael Norel은 non-unit 버전에 대해 다음의 개선을 제안한다.
d = dot(u, v)와 w = cross(u, v)의 값들이 어쨋든 계산되기 때문에, sqlength(u) * sqlength(v)의 값은 다른 방식으로 연산될 수 있다. 즉, d * d + sqlength(w). 다음의 코드는 적어도 세 번의 곱을 덜 한다:
glm::quat fromtwovectors(const glm::vec3 &u, const glm::vec3& v) { glm::vec3 w = cross(u, v); glm::quat q = glm::quat(glm::dot(u,v), w.x, w.y, w.z); q.w += length(q); return normalize(q); }
또한, Marc B. Reynolds는 unit version에서, 최종 표준화 factor는 sqrt((1 + dot(u, v))^2 + sqlength(cross(u,v)))이고, 이것은 sqrt(2 +2dot(u,v))로 줄어든다, sin^2+ cos^2 항등식 떄문에. 그것은 다음의 가급적 개선된 버전으로 이끈다:
glm::quat fromtwovectors(const glm::vec3 &u, const glm::vec3& v) { float m = sqrt(2.f + 2.f * dot(u,v)); glm::vec3 w = (1.f / m) * glm::cross(u,v); return glm::quat(0.5f * m, w.x, w.y, w.z); }
Final Version
나는 방향들 중의 하나에서 다른 것으로 transform시킬 두 개의 임의의 방향 벡터에서 쿼터니언을 구성하는 첫 번째 글을 보였다. 그 자료는 신중히 두 벡터가 서로를 마주보는 특별 케이스를 빠뜨렸다. 그래서 나는 어떻게 orthogonal vector를 빠르게 고르는지를 설명하는 두 번째 article을 썼었다.
두 글로부터의 결과를 합치는 것은 간단할 것이다. 여기에 그 결과인 최종 함수가 있다:
/* Build a unit quaternion representing the rotation * from u to v. The input vectors need not be normalised. */ quat quat::fromtwovectors(vec3 u, vec3 v) { float norm_u_norm_v = sqrt(dot(u, u) * dot(v, v)); float real_part = norm_u_norm_v + dot(u, v); vec3 w; if (real_part < 1.e-6f * norm_u_norm_v) { /* If u and v are exactly opposite, rotate 180 degrees * around an arbitrary orthogonal axis. Axis normalisation * can happen later, when we normalise the quaternion. */ real_part = 0.0f; w = abs(u.x) > abs(u.z) ? vec3(-u.y, u.x, 0.f) : vec3(0.f, -u.z, u.y); } else { /* Otherwise, build quaternion the standard way. */ w = cross(u, v); } return normalize(quat(real_part, w.x, w.y, w.z)); }
----------------------------------------------------------------------------------------------
두 번째 Article 바로 번역
On picking an orthogonal vector (and combing coconuts)
다음의 문제를 고려해라: 3D 공간에서 non-zero vector v가 주어지고, v에 orthogonal한 벡터 w를 찾아라 (즉, dot(v, w) = 0인).
처음 보기에, 이것은 쉬운 것처럼 보인다. 모든 가능한 v에 대해, 한 평면에 놓여있는 무한한 orthogonal vectors들이 있다:
우리는 하나를 고를 필요가 있다. 그러나, 그것은 보기와는 다르게 사소하지 않다.
First attempts
댓글 없음:
댓글 쓰기