Sign in to follow this  
adamv215

Help with Gram-Schmidt Orthonormalization code

Recommended Posts

adamv215    100
Hi everyone. I took this Orthonomalization code from the JSci library and adapted it, but I think there's wonky going on. I pass in [0, 2, 0], [0, 0, 2] and a random third vector. As I unerstand G-S orthonormalization, it should create a set of vectors with the third vector being perpendicular to the first two. So the result should be either [1, 0, 0] or [-1, 0, 0]. But instead I get [0, 0, 0]. I wonder if anybody could catch any mistakes I was making in my code:
public Vector3d[] orthonormalize(Vector3d[] vecs) {
    for(int i = 0; i < vecs.length; i++) {
        for(int j = 0; j < i; j++) {
            Vector3d temp = new Vector3d(vecs[j]);  //Sanity check to make sure nothing is being overwritten
            double dot = temp.dot(vecs[i]);
            Vector3d temp2 = new Vector3d(vecs[j]); //Sanity check to make sure nothing is being overwritten
            temp2.scale(dot);
            vecs[i].sub(temp2);
            vecs[i].normalize();
        }
    }
Also, this is my first post, so I apologize if it is in the wrong forum. Thanks!! p.s. I realize I can just use the cross-product to find a normal but will eventually need to run this code in more than three dimensions.

Share this post


Link to post
Share on other sites
adamv215    100
So far as I can tell. Scale just multiplies each element by a double. Sub is just vector subtraction and Normalization is just Normalization too. Hmmm...

Share this post


Link to post
Share on other sites
Halifax2    295
Are you sure that's correct? I thought that you were supposed to subtract the sum of all projections to get the perpendicular, not subtract the projection, normalize, subtract the projection from the modified vector, normalize, ad infinitum. For example, the algorithm is (in pseudo-code):

out[0] = normalize(orig[0])
for( i = 1; i < orig.length; ++ i )
{
Vector3D accum = Vector3D(0,0,0)

// Calculate the sum of all projections
for( j = 0; j < i - 1; ++ j )
{
accum += project(vecs[j], vecs[i])
// e.g. project vecs[i] onto vecs[j]
}

// Subtract the projections from the original vector and normalize
out[i] = normalize(orig[i] - accum);
}




Your code appears to be doing the following:

for( i = 0; i < orig.length; ++ i )
{
for( j = 0; j < i; ++ j )
{
vecs[i] -= normalize(project(vecs[j], vecs[i]));
}
}




Essentially the error appears to be in the fact that you aren't projecting onto the original vector. You project onto the modified vector.

I may be incorrect in some of my assumptions, so please disagree and explain.

Share this post


Link to post
Share on other sites
Sneftel    1788
Your pseudocode of adamv's method isn't quite correct, since it shows normalization before subtraction rather than after subtraction, but your description of adamv's method is correct. It doesn't really matter which you do. The important invariant is that after step (i,j), i is perpendicular to all vectors 0..j. This is maintained in either case.

Share this post


Link to post
Share on other sites
Halifax2    295
Your right, it was incorrect; I wasn't sure if that invariant was maintained or not, I couldn't think it through. Given that it is maintained though and the code above is correct, then the only other guess I have is that the third vector is either linearly dependent to the first or second vector.

Share this post


Link to post
Share on other sites
adamv215    100
Sneftel's first comment may have been closer to the problem. Probably because I'm assigning random doubles to the third vector when I do the first step of the normal calculation (x*x+y*y+z*z) the java debugger reports a value of "infinity." While I'm not sure the Gram-Schmidt code works, I'll take this particular issue over to the general programming forum.

Share this post


Link to post
Share on other sites
Halifax2    295
Yeah, Sneftel is probably the closet. Just to test the algorithm I wrote up a small sample program for myself (I was bored anyways); you can notice in the demo vectors included the numerical instability of the gram-schmidt orthonormalization e.g. the vectors aren't fully orthogonal, the angle between them evaluates to something like 1.1E-16.

If you want to take away the numerical instabilities, then I would recommend (once you get it working) you read the wikipedia page, specifically "Numerical Stability."

At any rate, I hope you get it working!

EDIT:

I like this Pennsylvania representation. :)

Share this post


Link to post
Share on other sites
grhodes_at_work    1385
Keep in mind that this forum is focused on math and physics for game development. Now, of course Gram-Schmidt orthonormalization absolutely has its place in game development. It is also a common academic topic. If you are using this for game development, it would be appropriate for you to say what you are using it for. If you are not using it for game development, it is still appropriate for you to tell us what you are working on, and also to review the Forum FAQ, which outlines the forum's policies for academic discussions of math and physics.

Share this post


Link to post
Share on other sites
adamv215    100
So noted. Thank you for taking the time to inform me of this, I will read the FAQ and abide by it in the future. My explanation follows:

I have been interested in mesh generation in a long time, especially for a combat space sim. I have thought that this was the logical next step from games like Elite, and I found a neat little paper on using Vonoroi Diagrams for mesh generation.

From what I have gathered, Vonoroi Diagrams are created by projecting the points into a Dim+1 space and calculating a Convex Hull. But qhull is so huge and a lot of what was there was unnecessary to me (and I prefer Java) that I found myself struggling to use the library. I then found a pretty well written Quickhull implementation in Java in 3D could be expanded into higher dimensions pretty easily. But I found myself struggling to compute the Normal Vector of the hyperplanes I was generating. I tried Gaussian Elimination, but sometimes this required division by zero, so I moved to Gram-Schmidt instead.

Thanks,
adamv215

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this