iDevGames Forums
A couple questions on 3d math - Printable Version

+- iDevGames Forums (http://www.idevgames.com/forums)
+-- Forum: Development Zone (/forum-3.html)
+--- Forum: Graphics & Audio Programming (/forum-9.html)
+--- Thread: A couple questions on 3d math (/thread-3872.html)



A couple questions on 3d math - TomorrowPlusX - Sep 18, 2006 06:53 AM

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,


A couple questions on 3d math - Jake - Sep 18, 2006 09:53 AM

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?)


A couple questions on 3d math - unknown - Sep 18, 2006 10:24 AM

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.


A couple questions on 3d math - imikedaman - Sep 18, 2006 11:21 AM

An explanation on how to invert matrices in OpenGL can be found here:
http://gpwiki.org/index.php/MathGem:Fast_Matrix_Inversion


A couple questions on 3d math - TomorrowPlusX - Sep 18, 2006 01:45 PM

imikedaman Wrote:An explanation on how to invert matrices in OpenGL can be found here:
http://gpwiki.org/index.php/MathGem:Fast_Matrix_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!


A couple questions on 3d math - DoG - Sep 19, 2006 02:38 PM

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.