'벡터'에 해당되는 글 2건

  1. 2009.02.03 벡터는 벡터일 뿐
  2. 2008.06.06 SSE의 오묘함 2
dev.log2009. 2. 3. 22:19
벡터는 벡터일 뿐, 3차원 벡터는 좌표 3개가 아니다.
3D게임 프로그래밍을 시작하는, 특히 2D게임을 만들다가 3D로 전향하려는 사람들에게서 잘 드러나는 것 같은데, 벡터를 벡터로 보지 않고 계산해야 할 좌표 3개로 보는 경향이 있는 것 같다. 사실, 벡터를 조작하여 얻고자 하는 무언가는 대부분의 경우 주어진 벡터로부터 특정한 방향을 가진 벡터, 혹은 주어진 벡터의 길이와 관계된 것일 경우가 많다. 이런 것을 구하기 위해서는 벡터를 좌표 3개로 분할할 필요가 없는데도, 습관처럼 좌표를 쪼개서 계산한다면, 벡터 연산에 익숙하지 않기 때문일 것 같다.

3차원 공간상의 벡터에 정의된 연산은 대략 반전, 덧셈, 뺄셈, 실수 곱셈, 실수 나눗셈, 그리고 내적과 외적이 거의 전부다. 이를 적절히 이용하면 간략한 표현으로 원하는 결과를 얻을 수 있는데, 좌표로 쪼개서 계산하면 고려해야할 예외적 상황이 증가하여 오히려 코드가 복잡해지는 결과가 초래되곤 한다.

오늘 본 문제는 '어떤 점에서 선분에 내린 수선의 발이 그 선분을 분할하는 비율'이다. 조금 더 상세히 정의하자면 선분의 시작점 s, 끝점 t가 주어졌을때, 임의의 점 p에서 내린 수선의 발 f에 대해서 |sf|/|st|의 값을 구하라는 것이다. 이런 문제는, 고등학교 수학시간에도 나오는 문제이지만, 이를 좌표로 쪼개서 풀려고 다음과 같이 하는 것을 보았다. (정확히 기억하고 있진 않을 수도 있지만)
  1. 각 벡터 p, s, t를 한 평면, 이를테면 xy평면에 투영하여 p', s', t'으로 한다.
  2. xy평면에서 선분s't'의 방향은 ( tx-sx, ty-sy )이다.
  3. xy평면에서 선분s't'과 직교하는 선분의 방향 d는 ( -ty+sy, tx-sx )이다.
  4. p'의 s'에 대한 변위 o를 구한다. 이때 o( p'x-s'x, p'y-s'y )이다.
  5. 점 o를 지나고 d와 평행한 직선과 선분 s't'의 교점 q를 연립방정식으로 구하면 그 점이 수선의 발.
  6. qx/(tx-sx)가 구하고자 하는 값이다.

대충 돌아가기는 하겠지만, 몇가지 예외사항이 존재한다. 우선 과정1에서 선분st가 xy평면에 수직일 경우 다른 평면을 골라 투영해야 한다. 또한 6단계에서 선분s't'이 y축에 평행하면 0/0의 꼴이 되어 부정이 되므로 역시 다른 좌표축을 골라야 한다. 여러가지 복잡한 예외를 고려해야 하기 때문에 결과적으로는 코드 길이가 원래보다 2배 정도 길어지게 마련이다. 하지만, 이는 내적의 성질을 이용해서 다음과 같은 간략한 방법으로 구할 수 있다.

이렇게 하면 선분의 방향이 어떠하든 일반적인 해를 구할 수 있다. 내적은 두 벡터의 길이의 곱에 사잇각의 코사인을 곱한 것과 같다는 성질을 이용한 것인데, 벡터의 길이에 코사인을 곱하면 다른 쪽 벡터에 투영된 길이가 나온다는 점을 생각하면 되겠다.

물론, 이런 방법으로는 길이의 비만 구할 수 있기 때문에 수선의 발 자체가 필요한 경우에는 따로 계산을 해야 하겠지만, 그때에도 위 결과로부터 직접 구할 수 있다. 선분 ts와 선분 fs의 길이의 비를 이미 구했으므로, 선분 ts를 이 비율대로 곱하기만 하면 수선의 발이 직접 구해진다. 따라서 수선의 발 f는 다음과 같이 구할 수 있다.

벡터에 정의된 연산이 몇개가 안되므로, 성급하게 벡터를 좌표로 쪼개려고 하는 것보다는, 벡터 연산으로 요리조리 조작해 보는 편이 더 일반적으로 유효한 결과를 가져온다. 3D프로그래밍의 시작은 벡터를 좌표의 집합이 아닌, 벡터 그 자체로 보는 것이라고나.


Posted by uhm
dev.log2008. 6. 6. 09:11
요즘은 SSE 수학 클래스 만들기에 열을 올리고 있다. 아직 행렬 클래스는 (노느라) 바뻐서 못만들고는 있지만 벡터 클래스는 요모조모 손봐가면서 테스트중. 전에 쓴 대로, SSE는 확실히 성능향상에 도움은 된다. 그런데 그게 좀 오묘하다. 명령의 순서가 희한하게 걸리면 오히려 속도가 느려진다.

다음은 주어진 연산순서를 각각 1000만번씩 반복했을 때 걸린 시간을 10번씩 측정한 결과이다. (CPU - Core2Duo E6600, Memory - DDR2 PC2-5300 6GB, 32bit application on 64bit Windows)

연산순서 기본구현
SSE구현
cross->length->dot 0.307732
0.307663
0.307266
0.307756
0.308705
0.307553
0.307968
0.307418
0.307018
0.307521
0.173833
0.173594
0.17392
0.17435
0.173717
0.17394
0.173853
0.173623
0.173899
0.173708
cross->dot->length 0.308549
0.307164
0.307153
0.307491
0.307836
0.307346
0.307413
0.308006
0.307355
0.307362
0.366581
0.366742
0.366583
0.367308
0.366557
0.366436
0.36658
0.366613
0.36664
0.366343

첫번째 결과는 당연히 SSE가 빠르겠거니.. 하는 추정에 부합하는데, 두번째 결과를 보면 연산 순서만 바뀌었을 뿐인데 SSE구현쪽의 결과가 훨씬 더 느리다. 반면, 기본 구현은 연산 순서랑 관계없이 거의 일정한 속도를 보장해 준다. 내 생각으로는 SSE명령을 수행 후에 메모리에 연산 결과를 쓰면서 레지스터중 일부를 초기화하고 메모리에 억세스하는 비용이 크기 때문일거 같은데, 어셈코드를 봐도 사실 정확한 건 잘 모르겠다. (공부하자)

들쭉날쭉한 속도를 개선하기 위해 요모조모로 테스트해 보다가 내적의 문제점을 깨달았다. 내적은, SSE로 구현할 때 애로사항이 많다. SSE의 기본은 (사실은 SSE가 아니라 그 근간인 SIMD라고 해야 맞지만) 부동소수 4개를 각기 따로따로 별도의 데이터 패쓰를 거쳐서 처리한다는, 일종의 data parallelism 인데, 내적은 따로따로 처리되어야 하는 데이터 3개, 즉 xyz좌표를 필수적으로 한군데에서 처리하도록 해야 한다. SSE의 설계 이념과 매우 동떨어진 연산이 될 수 밖에 없다. 물론 SSE3에서는 패킹된 데이터4개를 하나로 더해주는 명령이 추가되었지만, 아직 SSE3를 지원하지 않는 CPU가 많이 사용되고 있으므로 적용하기는 약간 무리다. 따라서 지금으로썬 내적의 SSE 구현이 병목이므로 내적을 SSE를 쓰지 않은 일반 구현으로 바꾸는 것이 최선일듯.

SSE로 구현한 벡터클래스의 내적을 일반구현으로 바꿔서 테스트해봤다.다음이 그 결과다.
연산순서 기본구현 SSE구현 
cross->length->dot 0.31281
0.328362
0.326691
0.312755
0.312775
0.325598
0.318856
0.31096
0.325269
0.308923
0.163027
0.1661
0.163232
0.163836
0.163102
0.162626
0.163001
0.169913
0.163181
0.162604
cross->dot->length 0.307485
0.307058
0.306469
0.307107
0.308057
0.306862
0.308106
0.307052
0.307446
0.306855
0.161829
0.161945
0.162188
0.16885
0.163519
0.161771
0.161899
0.162052
0.16204
0.162045
결과가 아주 이쁘장하다. 흡족하다. 저정도면 어느 상황에서도 비교적 빠른 속도를 보장해주는 벡터클래스가 되지 않을까.

그래서 결과적으로 완성된 코드는 다음과 같다.
#if !defined( __VECTOR3_H__ )
#define __VECTOR3_H__

#include <cmath>


#pragma intrinsic( sqrt )
// Internal Combustion Engine
namespace ice
{
template <class scalar_t, bool use_implicit = true > struct vector3_t
{
    typedef scalar_t scalar_type;
    scalar_type x, y, z;

public:
    vector3_t() {}
    vector3_t( scalar_type a, scalar_type b, scalar_type c ) : x(a), y(b), z(c) {}
    vector3_t( const scalar_type* xyz ) : x(xyz[0]), y(xyz[1]), z(xyz[2]) {}

    scalar_type operator[] ( int i ) const { return *(&x + i); }
    scalar_type& operator[] ( int i ) { return *(&x +i); }

    vector3_t& operator+= ( const vector3_t& v )
    {
        x += v.x; y += v.y; z += v.z;
        return *this;
    }
    vector3_t operator+ ( const vector3_t& v ) const
    {
        return vector3( x+v.x, y+v.y, z+v.z );
    }
    vector3_t& operator-= ( const vector3_t& v )
    {
        x -= v.x; y -= v.y; z -= v.z;
        return *this;
    }
    vector3_t operator- ( const vector3_t& v ) const
    {
        return vector3( x-v.x, y-v.y, z-v.z );
    }
    vector3_t& operator*= ( scalar_type s )
    {
        x *= s; y *= s; z *= s;
        return *this;
    }
    vector3_t operator* ( scalar_type s ) const
    {
        return vector3( x*s, y*s, z*s );
    }
    vector3_t& operator/= ( scalar_type s )
    {
        x /= s; y /= s; z /= s;
        return *this;
    }
    vector3_t operator/ ( scalar_type s ) const
    {
        return vector3( x/s, y/s, z/s );
    }

    static scalar_type dot( const vector3_t& v1, const vector3_t& v2 )
    {
        return v1.x*v2.x + v1.y*v2.y + v1.z*v2.z;
    }
    static vector3_t cross( const vector3_t& v1, const vector3_t& v2 )
    {
        return vector3_t(  v1.y * v2.z - v1.z * v2.y,
                                v1.z * v2.x - v1.x * v2.z,
                                v1.x * v2.y - v1.y * v2.x );
    }

    scalar_type squared() const { return x*x + y*y + z*z; }
    scalar_type length() const { return std::sqrt(squared()); }
    void normalize()
    {
        *this /= length();
    }
    vector3_t unit() const
    {
        return *this/length();
    }
};
}

#if defined( _USE_SSE_ )
#pragma pack( push, 16 )

#include <intrin.h>
#define ALIGN_MMX __declspec(align(16))

namespace ice
{

template<> struct ALIGN_MMX vector3_t<float, false>
{
    typedef float scalar_type;
    scalar_type x, y, z;

    vector3_t() {}
    vector3_t( scalar_type _a, scalar_type _b, scalar_type _c )
    {
        _mm_store_ps( &x, _mm_set_ps( 0, _c, _b, _a ) );
    }
    vector3_t( const scalar_type* xyz )
    {
        _mm_store_ps( &x, _mm_load_ps( xyz ) );
    }
    vector3_t( const vector3_t& v )
    {
        _mm_store_ps( &x, _mm_set_ps( 0, v.z, v.y, v.x ) );
    }

    scalar_type operator[] ( int i ) const { return *(&x + i); }
    scalar_type& operator[] ( int i ) { return *(&x +i); }

    vector3_t& operator+= ( const vector3_t& v )
    {
        _mm_store_ps( &x, _mm_add_ps( _mm_load_ps(&x), _mm_load_ps(&v.x) ) );
        return *this;
    }
    vector3_t operator+ ( const vector3_t& v ) const
    {
        return vector3_t( _mm_add_ps( _mm_load_ps(&x), _mm_load_ps(&v.x) ) );
    }
    vector3_t& operator-= ( const vector3_t& v )
    {
        _mm_store_ps( &x, _mm_sub_ps( _mm_load_ps(&x), _mm_load_ps(&v.x) ) );
        return *this;
    }
    vector3_t operator- ( const vector3_t& v ) const
    {
        return vector3_t( _mm_sub_ps( _mm_load_ps(&x), _mm_load_ps(&v.x) ) );
    }
    vector3_t& operator*= ( scalar_type s )
    {
        _mm_store_ps( &x, _mm_mul_ps( _mm_load_ps(&x), _mm_set1_ps( s ) ) );
        return *this;
    }
    vector3_t operator* ( scalar_type s ) const
    {
        return vector3_t( _mm_mul_ps( _mm_load_ps(&x), _mm_set1_ps( s ) ) );
    }
    vector3_t& operator/= ( scalar_type s )
    {
        _mm_store_ps( &x, _mm_div_ps( _mm_load_ps(&x), _mm_set1_ps( s ) ) );
        return *this;
    }
    vector3_t operator/ ( scalar_type s ) const
    {
        return vector3_t( _mm_div_ps( _mm_load_ps(&x), _mm_set1_ps( s ) ) );
    }

    static scalar_type dot( const vector3_t& v1, const vector3_t& v2 )
    {
        return v1.x*v2.x + v1.y*v2.y + v1.z*v2.z;
    }
    static vector3_t cross( const vector3_t& v1, const vector3_t& v2 )
    {
        //v1.y * v2.z - v1.z * v2.y
        //v1.z * v2.x - v1.x * v2.z
        //v1.x * v2.y - v1.y * v2.x
        __m128 a = _mm_shuffle_ps( *(__m128*)&v1, *(__m128*)&v1,
_MM_SHUFFLE( 3, 1, 2, 0 ) );
        __m128 c = _mm_shuffle_ps( *(__m128*)&v2, *(__m128*)&v2,
_MM_SHUFFLE( 3, 2, 0, 1 ) );
        __m128 b = _mm_shuffle_ps( *(__m128*)&v1, *(__m128*)&v1,
_MM_SHUFFLE( 3, 2, 0, 1 ) );
        __m128 d = _mm_shuffle_ps( *(__m128*)&v2, *(__m128*)&v2,
_MM_SHUFFLE( 3, 1, 2, 0 ) );
        return vector3_t( _mm_sub_ps( _mm_mul_ps( a, c ), _mm_mul_ps(b,d) ) );
    }

    scalar_type squared() const { return dot(*this,*this); }
    scalar_type length()  const
    {
        return _mm_sqrt_ss( _mm_set_ss( squared() ) ).m128_f32[0];
    }
    void normalize()
    {
        _mm_store_ps( &x, _mm_div_ps( _mm_load_ps(&x),
_mm_sqrt_ps( _mm_set1_ps( squared() ) ) ) );
    }
    vector3_t unit() const
    {
        return vector3_t( _mm_div_ps( _mm_load_ps(&x),
_mm_sqrt_ps( _mm_set1_ps( squared() ) ) ) );
    }

private:
    vector3_t( __m128 v )
    {
        _mm_store_ps( &x, v );
    }
};

typedef vector3_t<float, false>  vector3;
}
#pragma pack ( pop )
#else

namespace ice
{
typedef vector3_t<float>  vector3;
}

#endif // #if defined( _USE_SSE_ )

#endif // #if defined( _VECTOR3_H_ )

Posted by uhm