mat4 - opengl projection matrix




Précision inverse de la matrice (2)

Après quelques expériences, je constate que (en parlant de transformations, pas de matrice), la diagonale (facteurs d’échelle) de la matrice ( m , avant d’inverser) est la principale responsable de la valeur déterminante.

Je compare donc le produit p= m[0] · m[5] · m[10] · m[15] (si tous sont! = 0) avec le déterminant. Si elles sont similaires 0.1 < p/det < 10 je peux "faire confiance" en quelque sorte à la matrice inverse. Sinon, j'ai des problèmes numériques qui conseillent de changer de stratégie de rendu.

J'ai un grand monde, environ 5 000 000 x 1 000 000 d'unités. La caméra peut être proche d'un objet ou assez loin pour voir le monde entier.
J'obtiens la position de la souris en coordonnées du monde en effectuant une projection (Z provient du tampon de profondeur). Le problème est que cela implique une matrice inverse . Lorsque vous utilisez des grands et des petits nombres (par exemple, en vous éloignant de l'origine et en les redimensionnant pour voir plus de monde) simultanément, les calculs deviennent instables.

En essayant de voir la précision de cette matrice inverse, je regarde le déterminant. Idéalement, il ne sera jamais nul, en raison de la nature des matrices de transformation. Je sais qu'être «dét» une petite valeur ne signifie rien en soi, cela peut être dû à de petites valeurs dans la matrice. Mais cela peut aussi indiquer que les chiffres sont erronés.

Je sais aussi que je peux calculer l'inverse en inversant chaque transformation et en les multipliant. Fournit-il plus de précision?

Comment savoir si ma matrice est en train de dégénérer, de subir des problèmes numériques?


pour commencer, voir Comprendre les matrices de transformation homogène 4x4

  1. Amélioration de la précision des matrices cumulatives (normalisation)

    Pour éviter la dégénérescence de la matrice de transformation, sélectionnez un axe principal. En général, je choisis Z comme il est généralement question d’affichage ou de direction dans mes applications. Exploitez ensuite le produit croisé pour recalculer / normaliser le reste des axes (qui devrait être perpendiculaire l’un à l’autre et, à moins que l’échelle ne soit utilisée, indiquer également la taille de l’unité). Cela ne peut être fait que pour les matrices orthogonales / orthonormales, donc pas de biais ou de projections ...

    Vous n'avez pas besoin de faire cela après chaque opération, mais vous devez simplement créer un compteur d'opérations effectuées sur chaque matrice. Si un seuil est franchi, réglez-le et réinitialisez le compteur.

    Pour détecter la dégénérescence de telles matrices, vous pouvez tester l'orthogonalité par produit scalaire entre deux axes (il doit être égal à zéro ou très proche de celui-ci). Pour les matrices orthonormales, vous pouvez également tester la taille d'unité des vecteurs de direction d'axe ...

    Voici à quoi ressemble ma normalisation de matrice de transformation (pour les matrices orthonormales ) en 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
            }
    //---------------------------------------------------------------------------

    Les opérations vectorielles ressemblent à ceci:

    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. Amélioration de la précision pour les matrices non cumulatives

    Votre seul choix est d'utiliser au moins la double précision de vos matrices. Le plus sûr est d'utiliser GLM ou votre propre matrice mathématique basée au moins sur le type de données double (comme ma classe de reper ).

    Une alternative bon marché utilise des fonctions à double précision comme

    glTranslated
    glRotated
    glScaled
    ...

    ce qui, dans certains cas, aide mais n’est pas sûr car l’implémentation OpenGL peut le tronquer en float . De plus, il n'existe pas encore d'interpolateur matériel 64 bits, de sorte que tous les résultats itérés entre les étapes de pipeline sont tronqués en float s.

    Parfois, le cadre de référence relatif est utile (conservez donc des opérations sur des valeurs d'amplitude similaires), par exemple, voir:

    • amélioration de la précision des intersections entre rayons et ellipsoïdes

    Aussi, si vous utilisez vos propres fonctions mathématiques matricielles, vous devez également tenir compte de l'ordre des opérations afin de toujours perdre le moins de précision possible.

  3. Matrice pseudo-inverse

    Dans certains cas, vous pouvez éviter de calculer la matrice inverse à l'aide de déterminants, d'un schéma de Horner ou d'une méthode d'élimination de Gauss, car vous pouvez également exploiter le fait que Transpose de la matrice de rotation orthogonale est également son inverse . Voici comment cela se fait:

    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;
            }

    Ainsi, la partie rotationnelle de la matrice est transposée, la projection reste A*inverse(A)=unit_matrix et la position d'origine est recalculée de sorte A*inverse(A)=unit_matrix Cette fonction est écrite de sorte qu'elle peut être utilisée sur place afin d'appeler

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

    conduire à des résultats valables aussi. Cette méthode de calcul Inverse est plus rapide et plus sûre numériquement, car elle dépend beaucoup moins d’opérations (pas de récursion ni de réductions, pas de divisions ). Bien sûr, cela ne fonctionne que pour les matrices 4x4 orthogonales homogènes !!! *

  4. Détection de l'inverse erroné

    Donc si vous avez la matrice A et son inverse B alors:

    A*B = C = ~unit_matrix
    

    Donc, multipliez les deux matrices et vérifiez l'unité de matrice ...

    • abs total de tous les éléments non diagonaux de C doit être proche de 0.0
    • tous les éléments diagonaux de C doivent être proches de +1.0






matrix