A couple questions on 3d math

Sage
Posts: 1,199
Joined: 2004.10
Post: #1
I've got a couple questions on 3d math, specifically about inverting 4x4 matrices, and creating planes from triangles. In both cases, I have code that works, I'm really just asking because I want to understand -- better -- what's going on.

First, regarding matrices. I've been using a nice set of 3d math classes/methods which I adapted from demos from humus.ca. His mat4 class has an "inverse" method, which returns the inverse of a mat4. I've used his inverse code to inverse the modelview for certain kinds of particle systems as well as for inversing the transform on a model for the purposes of stencil shadowing. Unfortunately, while implementing my own version of gluUnProject I saw that MESA uses a different approach to inverting a 4x4 matrix ( they use the approach described in the matrix faq which I remember using in high school in my calculus classes ). I want to know why there are these two approaches ( I'll post the two versions below ) and which is "correct" and why.

Code:
// this one's from MESA's gluUnProject and works for gluUnproject
bool mat4::inverse( mat4 &inv ) const
{
    int i, j, k, swap;
    double t;
    float temp[4][4];
    
    for (i=0; i<4; i++)
    {
        for (j=0; j<4; j++)
        {
            temp[i][j] = this->mat[i*4+j];
        }
    }

    inv.identity();
    
    for (i = 0; i < 4; i++)
    {
        /*
         ** Look for largest element in column
         */
        swap = i;
        for (j = i + 1; j < 4; j++)
        {
            if (fabs(temp[j][i]) > fabs(temp[i][i]))
            {
                swap = j;
            }
        }
        
        if (swap != i) {
            /*
             ** Swap rows.
             */
            for (k = 0; k < 4; k++)
            {
                t = temp[i][k];
                temp[i][k] = temp[swap][k];
                temp[swap][k] = t;
                
                t = inv[i*4+k];
                inv[i*4+k] = inv[swap*4+k];
                inv[swap*4+k] = t;
            }
        }
        
        if ( fabsf( temp[i][i] ) < EPSILON )
        {
            /*
             ** No non-zero pivot.  The matrix is singular, which shouldn't
             ** happen.  This means the user gave us a bad matrix.
             */
            return false;
        }
        
        t = temp[i][i];
        for (k = 0; k < 4; k++)
        {
            temp[i][k] /= t;
            inv[i*4+k] /= t;
        }
        for (j = 0; j < 4; j++)
        {
            if (j != i)
            {
                t = temp[j][i];
                for (k = 0; k < 4; k++)
                {
                    temp[j][k] -= temp[i][k]*t;
                    inv[j*4+k] -= inv[i*4+k]*t;
                }
            }
        }
    }
    
    return true;
}

// this one's from Humus.ca and works for everything else
mat4 mat4::inverse() const {
    mat4 ret;
    float det;
    det = mat[0] * mat[5] * mat[10];
    det += mat[4] * mat[9] * mat[2];
    det += mat[8] * mat[1] * mat[6];
    det -= mat[8] * mat[5] * mat[2];
    det -= mat[4] * mat[1] * mat[10];
    det -= mat[0] * mat[9] * mat[6];
    det = 1.0f / det;
    ret[0] =  (mat[5] * mat[10] - mat[9] * mat[6]) * det;
    ret[1] = -(mat[1] * mat[10] - mat[9] * mat[2]) * det;
    ret[2] =  (mat[1] * mat[6] - mat[5] * mat[2]) * det;
    ret[3] = 0.0f;
    ret[4] = -(mat[4] * mat[10] - mat[8] * mat[6]) * det;
    ret[5] =  (mat[0] * mat[10] - mat[8] * mat[2]) * det;
    ret[6] = -(mat[0] * mat[6] - mat[4] * mat[2]) * det;
    ret[7] = 0.0f;
    ret[8] =  (mat[4] * mat[9] - mat[8] * mat[5]) * det;
    ret[9] = -(mat[0] * mat[9] - mat[8] * mat[1]) * det;
    ret[10] =  (mat[0] * mat[5] - mat[4] * mat[1]) * det;
    ret[11] = 0.0f;
    ret[12] = -(mat[12] * ret[0] + mat[13] * ret[4] + mat[14] * ret[8]);
    ret[13] = -(mat[12] * ret[1] + mat[13] * ret[5] + mat[14] * ret[9]);
    ret[14] = -(mat[12] * ret[2] + mat[13] * ret[6] + mat[14] * ret[10]);
    ret[15] = 1.0f;
    return ret;
}

OK, second question. I have code for generating a plane equation from a given CCW-wound triangle. I got it from nVIDIA's robust stencil shadows demo, IIRC. It's worked for me thus far, but over the last few days I've been implementing a render-to-texture variant of my user interface code for GUIs which are rendered on surfaces in the game world. I'm using gluUnProject to intersect a ray following the mouse cursor with the plane of the user interface to allow mouse clicks in the "world" to be forwarded to the in-game user interface.

Now, the way I'm doing it is to create a plane for the surface from a CCW triangle in world space ( part of the quad the UI is rendered on ). I do ray intersections with this plane ( and two orthogonal planes ) to get mouse coordinates on the surface of the GUI. The oddity is that I have to invert the fourth component of the plane ( where the plane is a*x + b*y + c*z = d ) to get the intersection. So, here's the code I've used all this time, and I'm hoping that perhaps it's wrong in the first place, and that the d component shouldn't be negated:

Code:
void plane::set( const vec3 &A, const vec3 &B, const vec3 &C )
{
    float p0, p1, p2, p3;

    p0 = A.y*(B.z-C.z) + B.y*(C.z-A.z) + C.y*(A.z-B.z);
    p1 = A.z*(B.x-C.x) + B.z*(C.x-A.x) + C.z*(A.x-B.x);
    p2 = A.x*(B.y-C.y) + B.x*(C.y-A.y) + C.x*(A.y-B.y);
    
    // the NVIDIA code says this should be negative, as is shown here
    p3 = -( A.x*(B.y*C.z - C.y*B.z) +
            B.x*(C.y*A.z - A.y*C.z) +
            C.x*(A.y*B.z - B.y*A.z) );    

    float rmag = 1.0f / sqrtf( p0*p0 +  p1*p1 + p2*p2 );
    
    a = p0 * rmag;
    b = p1 * rmag;
    c = p2 * rmag;
    d = p3 * rmag;
}

When I use the code above, it doesn't work. But when I negate 'd', it works correctly. So is the code above the right way, or the wrong way?

Thanks,
Quote this message in a reply
Member
Posts: 509
Joined: 2002.05
Post: #2
TomorrowPlusX Wrote:Now, the way I'm doing it is to create a plane for the surface from a CCW triangle in world space ( part of the quad the UI is rendered on ). I do ray intersections with this plane ( and two orthogonal planes ) to get mouse coordinates on the surface of the GUI. The oddity is that I have to invert the fourth component of the plane ( where the plane is a*x + b*y + c*z = d ) to get the intersection. So, here's the code I've used all this time, and I'm hoping that perhaps it's wrong in the first place, and that the d component shouldn't be negated:

Code:
void plane::set( const vec3 &A, const vec3 &B, const vec3 &C )
{
    float p0, p1, p2, p3;

    p0 = A.y*(B.z-C.z) + B.y*(C.z-A.z) + C.y*(A.z-B.z);
    p1 = A.z*(B.x-C.x) + B.z*(C.x-A.x) + C.z*(A.x-B.x);
    p2 = A.x*(B.y-C.y) + B.x*(C.y-A.y) + C.x*(A.y-B.y);
    
    // the NVIDIA code says this should be negative, as is shown here
    p3 = -( A.x*(B.y*C.z - C.y*B.z) +
            B.x*(C.y*A.z - A.y*C.z) +
            C.x*(A.y*B.z - B.y*A.z) );    

    float rmag = 1.0f / sqrtf( p0*p0 +  p1*p1 + p2*p2 );
    
    a = p0 * rmag;
    b = p1 * rmag;
    c = p2 * rmag;
    d = p3 * rmag;
}

When I use the code above, it doesn't work. But when I negate 'd', it works correctly. So is the code above the right way, or the wrong way?

Thanks,

I think its because you are using a*x + b*y + c*z = d as the equation of a plane, where it can also be written as a*x + b*y + c*z - d = 0. I learned it the second way in Calculus 3 today (isn't that odd that I read this thread two hours after that class?)
Quote this message in a reply
Sage
Posts: 1,403
Joined: 2005.07
Post: #3
well those two equations are equivalent, I think the error comes from one part of the algorithm being based on maths using ax+by+cz+d=0 and other parts using ax+by+cz=d.

as for the matrices, there many ways to code somthing like this (although I dont follow what is happening in the first one) as long as both find the inverse, they are both correct. I would definitly use somthing more like the second one, although you may want to put div by zero checks in.

Sir, e^iπ + 1 = 0, hence God exists; reply!
Quote this message in a reply
Member
Posts: 161
Joined: 2005.07
Post: #4
An explanation on how to invert matrices in OpenGL can be found here:
http://gpwiki.org/index.php/MathGem:Fast..._Inversion
Quote this message in a reply
Sage
Posts: 1,199
Joined: 2004.10
Post: #5
imikedaman Wrote:An explanation on how to invert matrices in OpenGL can be found here:
http://gpwiki.org/index.php/MathGem:Fast..._Inversion

That's cleared it up for me, actually. Thanks! What that tells me is that the simpler matrix inversion is good for matrices which only represent translation & rotation, but not scaling. And that corresponds to my usage. The more complex one works for all matrices, and that explains why it's used in gluUnProject, since that method needs to be able to invert the modelview * projection, which will include scaling.

And regarding the plane math, it certainly makes sense. I'm still a little fuzzy as to why some code wants one and some code the other, but I accept it.

Thanks, folks!
Quote this message in a reply
DoG
Moderator
Posts: 869
Joined: 2003.01
Post: #6
The MESA matrix inversion uses Gaussian elemination with pivoting, the humus version uses adjoint/determinant to determine the inverse. Both work on all invertible matrices, but they do react differently when the matrix becomes non-invertible, either due to numerical error or it being a truly non-invertible matrix.

With Gaussian elimination, it can be determined when things go tits up in most cases, while using the latter method all you get is NaN all over the place.
Quote this message in a reply
Post Reply