opengl - 행렬 역 정확도




math matrix (2)

몇 가지 실험을 거친 후 (매트릭스가 아닌 변환에 대해 말하면) 매트릭스의 대각선 (예 : 스케일링 계수)이 반전 값의 주요 원인임을 알 수 있습니다.

따라서 곱 p= m[0] · m[5] · m[10] · m[15] (모두! = 0 인 경우)를 결정자와 비교합니다. 이들이 0.1 < p/det < 10 과 비슷하다면 역행렬에서 어떻게 든 "신뢰할 수"있습니다. 그렇지 않으면 렌더링 전략을 변경하는 데 도움이되는 숫자 문제가 있습니다.

나는 약 5,000,000 x 1,000,000 단위의 큰 세계를 가지고 있습니다. 카메라는 물체 가까이에 있거나 세계를 볼 수있을만큼 먼 거리에있을 수 있습니다.
투영을 해제하여 월드 좌표에서 마우스 위치를 얻습니다 (Z는 깊이 버퍼에서 나옴). 문제는 행렬의 역수가 포함된다는 것입니다. 크고 작은 숫자 (예 : 원점에서 변환 및 더 많은 세계를보기 위해 크기 조정)를 동시에 사용하는 경우 계산이 불안정 해집니다.

역행렬 의 정확성을 보려고 노력하고 있습니다. 변환 행렬의 특성으로 인해 0이 아닌 것이 이상적입니다. 작은 값을 'det'한다는 것은 그 자체로는 아무것도 의미하지 않으며 매트릭스의 작은 값 때문일 수 있습니다. 그러나 숫자가 잘못되었다는 신호일 수도 있습니다.

또한 각 변환을 뒤집고 곱하여 역수를 계산할 수 있다는 것도 알고 있습니다. 더 정확한 정보를 제공합니까?

행렬이 퇴화되고 있는지, 숫자 문제가 있는지 어떻게 알 수 있습니까?


우선 4x4 동종 변환 행렬 이해를 참조하십시오.

  1. 누적 행렬의 정확도 향상 (정규화)

    변환 매트릭스의 변성 하려면 하나의 축을 기본으로 선택하십시오. 나는 일반적으로 내 앱에서 일반적으로 시야 또는 방향이므로 Z 를 선택했습니다. 그런 다음 교차 곱을 활용하여 나머지 축을 재 계산 / 정규화합니다 (이것은 서로 직각이어야하며 스케일을 사용하지 않으면 단위 크기). 이것은 직교 / 직교 정규 행렬에 대해서만 수행 할 수 있으므로 비대칭 또는 투영이 없습니다 ...

    모든 작업 후에 각 매트릭스에서 수행 된 작업 카운터를 만든 후 일부 임계 값을 초과 한 경우이를 정규화하고 카운터를 재설정하면이 작업을 수행 할 필요가 없습니다.

    이러한 행렬의 변성 감지 하려면 두 축 사이의 내적 (0 또는 그 근처)에 따라 직교성을 테스트 할 수 있습니다. 직교 정규 행렬의 경우 축 방향 벡터의 단위 크기에 대해서도 테스트 할 수 있습니다 ...

    다음은 C ++ 에서 변환 행렬 정규화가 ( 직교 정규 행렬의 경우) 어떻게 보이는지입니다.

    double reper::rep[16]; // this is my transform matrix stored as member in `reper` class
    //---------------------------------------------------------------------------
    void reper::orto(int test) // test is for overiding operation counter
            {
            double   x[3],y[3],z[3]; // space for axis direction vectors
            if ((cnt>=_reper_max_cnt)||(test)) // if operations count reached or overide
                    {
                    axisx_get(x);      // obtain axis direction vectors from matrix
                    axisy_get(y);
                    axisz_get(z);
                    vector_one(z,z);   // Z = Z / |z|
                    vector_mul(x,y,z); // X = Y x Z  ... perpendicular to y,z
                    vector_one(x,x);   // X = X / |X|
                    vector_mul(y,z,x); // Y = Z x X  ... perpendicular to z,x
                    vector_one(y,y);   // Y = Y / |Y|
                    axisx_set(x);      // copy new axis vectors into matrix
                    axisy_set(y);
                    axisz_set(z);
                    cnt=0;             // reset operation counter
                    }
            }
    
    //---------------------------------------------------------------------------
    void reper::axisx_get(double *p)
            {
            p[0]=rep[0];
            p[1]=rep[1];
            p[2]=rep[2];
            }
    //---------------------------------------------------------------------------
    void reper::axisx_set(double *p)
            {
            rep[0]=p[0];
            rep[1]=p[1];
            rep[2]=p[2];
            cnt=_reper_max_cnt; // pend normalize in next operation that needs it
            }
    //---------------------------------------------------------------------------
    void reper::axisy_get(double *p)
            {
            p[0]=rep[4];
            p[1]=rep[5];
            p[2]=rep[6];
            }
    //---------------------------------------------------------------------------
    void reper::axisy_set(double *p)
            {
            rep[4]=p[0];
            rep[5]=p[1];
            rep[6]=p[2];
            cnt=_reper_max_cnt; // pend normalize in next operation that needs it
            }
    //---------------------------------------------------------------------------
    void reper::axisz_get(double *p)
            {
            p[0]=rep[ 8];
            p[1]=rep[ 9];
            p[2]=rep[10];
            }
    //---------------------------------------------------------------------------
    void reper::axisz_set(double *p)
            {
            rep[ 8]=p[0];
            rep[ 9]=p[1];
            rep[10]=p[2];
            cnt=_reper_max_cnt; // pend normalize in next operation that needs it
            }
    //---------------------------------------------------------------------------

    벡터 연산은 다음과 같습니다.

    void  vector_one(double *c,double *a)
            {
            double l=divide(1.0,sqrt((a[0]*a[0])+(a[1]*a[1])+(a[2]*a[2])));
            c[0]=a[0]*l;
            c[1]=a[1]*l;
            c[2]=a[2]*l;
            }
    void  vector_mul(double *c,double *a,double *b)
            {
            double   q[3];
            q[0]=(a[1]*b[2])-(a[2]*b[1]);
            q[1]=(a[2]*b[0])-(a[0]*b[2]);
            q[2]=(a[0]*b[1])-(a[1]*b[0]);
            for(int i=0;i<3;i++) c[i]=q[i];
            }
  2. 비 누적 행렬의 정확도 향상

    유일한 선택은 행렬의 최소 double 정확도를 사용하는 것입니다. 가장 안전한 방법은 적어도 double 데이터 유형 (내 reper 클래스와 같은)을 기반으로 GLM 또는 자신의 행렬 수학을 사용하는 것입니다.

    저렴한 대안은 다음과 같은 배정도 기능을 사용하는 것입니다.

    glTranslated
    glRotated
    glScaled
    ...

    OpenGL 구현이 float 자를 수 있으므로 도움이되지만 안전하지 않은 경우도 있습니다. 또한 64 비트 HW 보간 기가 없으므로 파이프 라인 단계 사이의 모든 반복 결과가 float 으로 잘립니다.

    때로는 상대 참조 프레임이 도움이됩니다 (따라서 유사한 크기 값으로 작업을 유지).

    • 광선 및 타원체 교차 정확도 개선

    또한 자체 매트릭스 수학 함수를 사용하는 경우 연산 순서도 고려해야하므로 항상 가능한 한 적은 정확도를 잃게됩니다.

  3. 의사 역행렬

    경우에 따라 직교 회전 행렬의 조옮김도 역수 라는 사실을 이용할 수 있기 때문에 결정 요인이나 호너 방식 또는 가우스 제거 방법으로 역행렬의 계산을 피할 수 있습니다. 방법은 다음과 같습니다.

    void  matrix_inv(GLfloat *a,GLfloat *b) // a[16] = Inverse(b[16])
            {
            GLfloat x,y,z;
            // transpose of rotation matrix
            a[ 0]=b[ 0];
            a[ 5]=b[ 5];
            a[10]=b[10];
            x=b[1]; a[1]=b[4]; a[4]=x;
            x=b[2]; a[2]=b[8]; a[8]=x;
            x=b[6]; a[6]=b[9]; a[9]=x;
            // copy projection part
            a[ 3]=b[ 3];
            a[ 7]=b[ 7];
            a[11]=b[11];
            a[15]=b[15];
            // convert origin: new_pos = - new_rotation_matrix * old_pos
            x=(a[ 0]*b[12])+(a[ 4]*b[13])+(a[ 8]*b[14]);
            y=(a[ 1]*b[12])+(a[ 5]*b[13])+(a[ 9]*b[14]);
            z=(a[ 2]*b[12])+(a[ 6]*b[13])+(a[10]*b[14]);
            a[12]=-x;
            a[13]=-y;
            a[14]=-z;
            }

    따라서 행렬의 회전 부분이 A*inverse(A)=unit_matrix 투영은 그대로 유지되고 원점 위치가 다시 계산되므로 A*inverse(A)=unit_matrix 이 함수는 작성되어 내부 호출로 사용할 수 있습니다.

    GLfloat a[16]={values,...}
    matrix_inv(a,a);

    유효한 결과로 이어집니다. Inverse를 계산하는이 방법은 훨씬 적은 작업 (재귀 또는 축소가 없음 )을 보류하므로 더 빠르고 수치 적으로 안전합니다. 대략적으로 이것은 직교 균질 4x4 행렬에서만 작동합니다 !!! *

  4. 잘못된 역의 검출

    따라서 행렬 A 와 역 B 가 있다면

    A*B = C = ~unit_matrix
    

    두 행렬을 곱하고 단위 행렬을 확인하십시오 ...

    • C 대각선이 아닌 모든 요소의 절대 합은 0.0 가까워 야합니다.
    • C 모든 대각선 요소는 +1.0 가까워 야합니다




matrix