extensions for 2D convex shape

[MichaelRiss]

Hello Sébastien,

I tried to play around with 2D convex shapes and noticed that they are not completely supported yet. Then I tried to patch in what's missing and to my surprise succeeded in getting a basic 2D convex demo up and running. Maybe you want to take a look at it, comment and - if everything is ok - merge?
This pull request contains the parts of nphysics.

Note: I also created "demo2Dconvex" branches in my github repositories. They have the dependencies properly set so that they compile and can be tested.

Michi

[sebcrozet]

I'd prefer not to allow non_snake_case as this is non-idiomatic. Instead the "convention" I tend to use is to append the dimension number without the D. For example: triangle_unit_inertia_tensor_wrt_point2.

[MichaelRiss]

Is it possible that crates.io needs a while to update to the newest version of ncollide?
Travis-CI seems to pull in the old version of ncollide without triangle_center.

[MichaelRiss]

And yes, I will squash the commits in the end. ;)
Just tell me when you are ok with the current state.

To speed things up I squashed my master branch and added another branch fullHistory with all the individual commits.

[sebcrozet]

Last thing before a merge: I don't think your computations for the 2D triangle moment of inertia is correct. I think it should be bh³/36 + d^2 where b is the triangle basis length, h its height, and d the distance from the triangle centroid and the point relative to which we are computing the inertia momentum.

[MichaelRiss]

Indeed, I also think my calculation is wrong. I'm still trying to wrap my head around the math.

[MichaelRiss]

This is how far I came: After some digging I found a planar moment of inertia cheat sheet:
http://www.bau.uni-siegen.de/subdomains/bauinformatik/lehre/tm2/arbeitsblaetter/arbeitsblatt_08_flaechentraegheitsmomente_bsp.pdf (it's in German but the formulas should be readable).
I found your suggested formula there (without the term for the parallel axis theorem), but it's for I_y - the moment of inertia around the y-axis (careful, the axis orientation is a bit weird). I think we rather need I_p - the polar moment of inertia for rotations around the S point.

I then started to develop an algorithm but soon ended up with a lot of trigonometric function calls - and that's not the type of operations you would typically see in such a code. Lots of refinement and optimization would be necessary to get this approach into shape.
Then I grew impatient as this is something that for sure was solved before. Looking around how other people tackled the problem, I came across the Box2D project and its solution is very elegant and efficient (in Box2D/Collision/Shapes/b2PolygonShape.cpp). Therefore, I simply adapted their approach until we find something more elegant or efficient (which may be hard 😄).

What do you think about the general approach?
The demo looks ok with the new way to calculate the inertia, but in the long run I would like to add some tests that verify the calculation with well known shapes.

[MichaelRiss]

I'm a bit stuck right now with the tests. Using the test crate requires rust nightly due to the experimental feature status and ncollide does not compile with rust nightly (it seems it has problems with inferring types). Looking at the errors it's not clear to me whether changes in ncollide are needed or whether this is a compiler regression.

Could be this one:
rust-lang/rust#26279

[MichaelRiss]

It seems like the test crate is not needed at the moment. So I commented it out to get the tests to work with stable rust. I can now test square, rectangle and triangle inertia calculation - and they all pass.

Note: I had to break your design by introducing an additional function convex_hull_unit_angular_inertia2. The problem was that your function returns a matrix with the same dimensions as the input coordinate vectors. This is correct in 3D but not in 2D, there the result is a single value (Mat1). Please comment on that, is there something that I missed about your design?

Now that the tests pass and the demo works, I'm happy with my changes and done for the moment.

Please take a look at it, comment and merge if you like it.

[sebcrozet]

Note: I had to break your design by introducing an additional function convex_hull_unit_angular_inertia2. The problem was that your function returns a matrix with the same dimensions as the input coordinate vectors. This is correct in 3D but not in 2D, there the result is a single value (Mat1). Please comment on that, is there something that I missed about your design?

The dimensions of the input point P and of the inertia tensor matrix I are actually not forced to be equal. One method you could use to merge the functions ending with 2 with the other ones is to switch on the point dimension and return the appropriate result for the selected one. And if you need to create a matrix by hand, you can do it by indexing without knowing the actual type of I:

let mut i: I = na::zero();
i[(0, 0)] = na::cast(42.0); // the 2D inertia momentum. 

In your case, the following should work:

pub fn convex_hull_unit_angular_inertia<P, I>(dim: usize, points: &[P]) -> I
    where P: Point,
          I: Zero +
             Add<I, Output = I> +
             Mul<<P::Vect as Vect>::Scalar, Output = I> +
             IndexMut<(usize, usize), Output = <P::Vect as Vect>::Scalar>,
          P::Vect: Outer,
          <P::Vect as Outer>::OuterProductType: EigenQR<<P::Vect as Vect>::Scalar, P::Vect> +
                                                Mul<P, Output = P> +
                                                Add<<P::Vect as Outer>::OuterProductType, Output = <P::Vect as Outer>::OuterProductType> +
                                                Zero + Copy {
    assert!(dim == 2 || dim == 3);

    match dim {
        2 => {
            let convex_mesh = transformation::convex_hull2(points);

            unsafe {
                let (area, _, i): (_, _, <P::Vect as Vect>::Scalar) = convex_mesh_mass_properties2(&convex_mesh, na::one());

                let mut tensor: I = na::zero();
                tensor[(0, 0)] = i * (na::one::<<P::Vect as Vect>::Scalar>() / area);

                tensor
            }
        }
        3 => {
            let convex_mesh = transformation::convex_hull3(points);
            unsafe {
                let (vol, _, i): (_, _, I) = convex_mesh_mass_properties(&convex_mesh, na::one());

                i * (na::one::<<P::Vect as Vect>::Scalar>() / vol)
            }
        }
        _ => {
            unimplemented!()
        }
    }
}

And you cloud do the same to get rid of convex_mesh_mass_properties2 as well.

[MichaelRiss]

Wow! I still got a lot to learn about rust. Thanks for the hint!

Now I could realize your first suggestion but with the second convex_mesh_mass_properties I ran into a problem with the input parameters. The 3D version uses TriMesh and the 2D version Polyline.
Any idea how to untangle this?

[sebcrozet]

You're right, I forgot about the TriMesh/Polyline. I don't think there is a way to untangle this without introducing some weird traits.

Okay, so as soon as you remove the commented code (that you left to show me what could go wrong), and squash the commits, I will merge this!

[MichaelRiss]

Done.

[sebcrozet]

Thanks for all this work and your tenacity!

[MichaelRiss]

Darn! I just realized that the picking in the demo is broken. I had by mistake (C)SFML 2.3 installed instead of 2.1. The demos compiled with 2.3 but picking didn't work at all, so I didn't noticed the regression with convex shapes. I'm looking at the code now looking for a fix, but it looks deeply mathematical and complex. I hope I get somewhere.
But it was by no chance intentional to deliver broken code - I'm deeply sorry!

[sebcrozet]

If you are talking about a panic like "Attempted to translate the identity matrix" whenever you try to grab a convex object, don't worry about that. It's a bug independent from your work (originating from this line)!

I did not notice this either because I was running SFML 2.3 as well, but now that you spotted it, I should be able to fix this.

[sebcrozet]

This might be a bit tricky to fix actually. I need to modify a couple of algorithms to make it work without adding an overhead wrt. the current (bugged anyway) implementation.

I opened the issue #35 for this.

[MichaelRiss]

Yes, I'm referring to the 'Attempted to translate the identity matrix.' panic. Still no solution from my side so far. :-/

[sebcrozet]

I created the issue #35 on the wrong repository. It is now there: dimforge/ncollide#87, with a comment about how I am going to fix it.