Matrix4::invert breaks because of low precision in approx_eq

[jarrett]

When inverting a matrix, we check that the determinant is non-zero using approx_eq. If the determinant is approximately zero, the inversion fails.

But approx_eq only goes to a precision of 10^-5.

That's problematic, because projection matrices often have very small determinants. This is especially true when the camera is zoomed out very far. In that case, we're transforming a large volume of world space to the clip space volume. Which means that the determinant tends to be very small. (This is because the determinant is the ratio of the two volumes.)

[jarrett]

Admittedly, this can be worked around by building a separate ortho projection with very close near and far planes. This keeps the viewing volume small, resulting in a larger determinant. This works just fine for unprojecting a point in screen space. But it does mean we need two different versions of the projection matrix: one for projecting, and one for inverting. So the workaround is less than ideal.

[kvark]

One way to work around it is to only change your camera position/zoom as long as the projection matrix is invertible. Takes a bit of extra logic in your code, but hey - you've got to support near-zero determinant somehow.

[jarrett]

Well, even perfectly reasonable zoom levels yield determinants small enough to trigger approx_eq. For example, it's not unusual to have 1000 x 1000 x 100 world units on screen at once. But that's a determinant of 10^-8, which is well below the threshold for approx_eq.

Given that the precision of approx_eq is 10^-5, I think the largest possible viewing volume in world units is:

46 x 46 x 46 = 97,336

That seems to like a very small viewing volume for most 3d applications.

[kvark]

Just to clarify, the problem is only relevant with orthographic projections. Instead of expanding the projection to cover the world units, you can just scale the world down.

Another way to go would be to have a specialized invert variant for matrices that bear no rotation or perspective, thus working for the orthographic projections only. This specialized version would be highly optimized and not compute the determinant anywhere, but rather work out the equations in place.

[jarrett]

Just to clarify, the problem is only relevant with orthographic projections.

Why would that be true? The viewing volume for a perspective camera could also be large. That would result in a tiny determinant, correct?

Instead of expanding the projection to cover the world units, you can just scale the world down.

Scaling the world down is exactly what the camera's model-view-projection matrix does. It scales from a large volume, e.g. 1000 x 1000 x 100, to clip space, which is 2 x 2 x 2. And that scaling is itself the problem, because it's what makes the determinant so small.

I think the only way to solve this without workarounds is to increase the precision of the degeneracy check. Would it be a problem to do that?

[kvark]

Why would that be true? The viewing volume for a perspective camera could also be large. That would result in a tiny determinant, correct?

No, because a perspective projection fits the view not by scaling but by projecting. The determinant is going to be decent.

Scaling the world down is exactly what the camera's model-view-projection matrix does.

Scaling something is not a exactly a problem, inverting large scaling IS a problem. So, if your scale ends up in the Model part of the M-V-P, then it doesn't participate in the inverse, hence everything is good.

I think the only way to solve this without workarounds is to increase the precision of the degeneracy check. Would it be a problem to do that?

It's just a half-measure, it doesn't actually "solve" the problem.

[jarrett]

But to transform from screen space to world space, you need to multiply by inv(projection x model x view). So that big scale factor is necessarily built into the matrix that goes from screen to world.

To put it another way: if you only use the inverse projection matrix to transform from screen to world, you'll get the wrong world coorda, because you will have neglected to account for the model-view transform.

On Apr 24, 2015, at 5:07 PM, Dzmitry Malyshau notifications@github.com wrote:

Why would that be true? The viewing volume for a perspective camera could also be large. That would result in a tiny determinant, correct?

No, because a perspective projection fits the view not by scaling but by projecting. The determinant is going to be decent.

Scaling the world down is exactly what the camera's model-view-projection matrix does.

Scaling something is not a exactly a problem, inverting large scaling IS a problem. So, if your scale ends up in the Model part of the M-V-P, then it doesn't participate in the inverse, hence everything is good.

—
Reply to this email directly or view it on GitHub.

[kvark]

@jarrett You don't need inv(MVP) for the regular rendering, so I assume you are talking about re-projection from the screen space, used for mouse cursor tracing and such. This narrows down the problematic case even more. You can work around the inverse limitations for this case by just doing the inverse of the projection manually:

X = inv(P * V * M);
X = inv(V*M) * inv(P);

Basically, you can unproject a point manually (since orthographic projection is rather simple), then multiply by the computed inv(V*M) matrix that doesn't have the problem.

If you do that, you could contribute the unprojection method back into cgmath::Ortho.

[jarrett]

so I assume you are talking about re-projection from the screen space, used for mouse cursor tracing and such.

Yes, that is exactly right.

But here's the problem: One of the two inversions, inv(V*M) or inv(P), has to fail. If I scale in V * M, then inv(V * M) fails. If instead I use a large projection window in P, then inv(P) fails. In other words, there will always be at least one matrix that transforms from large world coordinates to small screen coordinates, and this matrix will necessarily have a small determinant.

Anyhow, is there some reason why the precision of the floating point comparison can't be increased? Do you know why 10^-5 was chosen as the epsilon?

[kvark]

One of the two inversions, inv(V*M) or inv(P), has to fail.

No. What I'm saying is - instead of computing inv(P) and then transforming a point, you can just manually un-project the point, involving no matrix inversion.

Let's say your orthographic projection is: w = P(q) = scale*q + offset, then you can un-project q without inversing any matrix: q = (w - offset) / scale.

Do you know why 10^-5 was chosen as the epsilon?

I bet it was an arbitrary decision. But whatever you want to change it to would be no less arbitrary, and that's why I don't like it as a solution.

[jarrett]

I bet it was an arbitrary decision. But whatever you want to change it to would be no less arbitrary, and that's why I don't like it as a solution.

That is a good point. So instead, how about this: An unsafe invert variant, maybe called unsafe_invert, that skips the determinant? Because there are situations where the programmer knows for sure the matrix is invertible, and computing the determinant is just a waste, or worse, a source of errors.

If you're OK with that idea, I'll implement it.

[kvark]

@jarrett sounds good!

[jarrett]

Here's a first attempt:

jarrett@6a17b8a

I'm interested in hearing your thoughts.

[jarrett]

I just pushed a new revision with the latest from bjs/cgmath-rs merged in:

https://github.com/jarrett/cgmath-rs

Any interest in getting these features into the crate?

[brendanzab]

I'll look at this after work - sorry for the delay. I definitely think there should be a way to sidestep the calculations, if the programmer thinks it is ok.

I'm also looking at improving approx_eq to handle small floating point values, and moving it to a separate crate. These are the articles I'm looking at:

• http://floating-point-gui.de/errors/comparison/
• http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
• https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/

[FredrikNoren]

Any updates on this?

[FredrikNoren]

I just noticed I don't need to use unsafe in my "hacked" mat4 invert (stolen from @jarrett I think):

pub fn mat4_invert(mat: Matrix4<f32>) -> Matrix4<f32> {
    let det = mat.determinant();
    let one: f32 = one();
    let inv_det = one / det;
    let t = mat.transpose();
    let cf = |i, j| {
        let mat = match i {
            0 => Matrix3::from_cols(t.y.truncate_n(j),
                                    t.z.truncate_n(j),
                                    t.w.truncate_n(j)),
            1 => Matrix3::from_cols(t.x.truncate_n(j),
                                    t.z.truncate_n(j),
                                    t.w.truncate_n(j)),
            2 => Matrix3::from_cols(t.x.truncate_n(j),
                                    t.y.truncate_n(j),
                                    t.w.truncate_n(j)),
            3 => Matrix3::from_cols(t.x.truncate_n(j),
                                    t.y.truncate_n(j),
                                    t.z.truncate_n(j)),
            _ => panic!("out of range")
        };
        let sign = if (i+j) & 1 == 1 {-one} else {one};
        mat.determinant() * sign * inv_det
    };

    Matrix4::new(cf(0, 0), cf(0, 1), cf(0, 2), cf(0, 3),
                 cf(1, 0), cf(1, 1), cf(1, 2), cf(1, 3),
                 cf(2, 0), cf(2, 1), cf(2, 2), cf(2, 3),
                 cf(3, 0), cf(3, 1), cf(3, 2), cf(3, 3))
}

Not sure what they changed but that compiles in latest rust nightly

[Andlon]

I just got bitten by this behavior. I actually suggest to not try to do an approximate comparison for the determinant, and instead simply check for exact equality with zero. Why? Simply put, in imprecise arithmetic, the determinant is a completely unreliable way to check if a matrix is invertible or not. Consider the following example.

Let I denote the 3x3 identity matrix, and let A = 1e-6 * I. Then A is easily invertible, with inverse 1e6 * I. However, it's determinant is 1e-18, which is below double-precision machine epsilon (1.11e-16), and would fail the current ULPS-based tests in cgmath.

Having matrices like this is also not at all uncommon. I just encountered this issue when dealing with inertia tensors and objects of very varying sizes.

[kvark]

@Andlon we only have so much precision. One way to improve it is to choose the scale where most numbers would be around e0. It's not reasonable to deal with 1e-6 values and expect things to work as well as around 1e0. At the end of the day, cgmath allows you to use f64.

[Andlon]

@kvark: I agree that it's not reasonable for things to work as well near 1e0 as near 1e-6, but I think it's reasonable to expect that things should work near 1e-6, and in fact, with all due respect, I consider this a bug. In terms of scale, your suggestion assumes that there is a single uniform scale where most numbers are in fact near e0, which may not be the case.

If you really want to determine if the inversion is successful, performing an LU(P) decomposition can give you a much better indication whether the matrix is non-singular to working precision or not, and a more accurate inverse to boot. However, if you want to use an explicit inverse, I suggest not using the determinant as an indicator, because it will very easily break when it should not. With the current API of cgmath, this also means that there is no way to override this behavior for the user.

I'd be happy to contribute some time to improve the situation and send a PR, if we could agree on desired behavior!

EDIT: Oh, and I ran into the issues when using Matrix3<f64>, so there's not more precision to be gained there, unfortunately!

[kvark]

@Andlon good points!
Sounds like all we need is a precise_invert() method using LU(P) decomposition?

[Andlon]

I think actually we should introduce LU decompositions into cgmath as a separate entity, because you should actually never invert a matrix. Rather, you're better off solving the system Ax = b. If you need to solve many such systems for the same A, an LU decomposition lets you accomplish this as efficiently as inverting the matrix and then repeatedly multiplying the inverse with the vector, and with much greater precision. However, in computer graphics, it has some merits to work with inverses, because most literature is presented this way, and inverses are easier to compose. Moreover, the disadvantages of explicit inverses are not so pronounced for such small dimensions, although they are still very real. In any case, it may be useful to users to have LU decomposition available.

If it turns out that inversion by LU decomposition is as fast as explicit inverses, we could just swap the current implementation. The potential gains in accuracy are truly enormous.

I'd be up for working on this, but at the moment I'm still slightly undecided between whether to use cgmath or nalgebra going further. I've used cgmath so far, and have generally been happy, but I miss the equivalent of nalgebra's UnitQuat. I'd have to make a decision here before I sink time into something like this!

[jarrett]

I'm not familiar with LU decompositions, so if you all believe it'll solve the problem, I'll defer to you.

Whatever you do, I'd just recommend documenting it in such a way that a typical programmer can understand. I think most of us know what it means to invert a matrix and when to do it. As Andlon said, many references on computer graphics teach us to use inversions. But I've never read one that mentions LU decomposition.

Here's one possible approach to documentation: Add a heading in the readme called "Inversion." Under it, explain that you recommend LU decompositions as a fast alternative and present some realistic sample code. E.g. an example of unprojecting mouse coordinates to a ray in world space.

[Andlon]

@jarrett: excellent point!

[brendanzab]

LU decomopositions would be great. I would also note that it might make sense to re-look at our current inversion code to ensure it is up to scratch. It's something I'm not too well versed in.

[brendanzab]

Maybe open an issue about UnitQuat and its rationale?

[tomaka]

For anyone having this issue I have a no-invert-fail branch in my fork: https://github.com/tomaka/cgmath-rs/tree/no-invert-fail