Meshless Deformations Based on Shape Matching

[raigan2]

Does anyone have any experience with implementing this paper? http://graphics.ethz.ch/~brunoh/downloa ... _SIG05.pdf

I've written two seperate solvers (in 2D), one based on Muller's source code ( http://graphics.ethz.ch/~mattmuel/demos ... Source.zip) which uses jacobi rotations to find the eigen-stuff, and another method which finds the eigen-stuff directly (http://www.dgp.toronto.edu/~hertzman/co ... eig2x2.cpp).

My problems are with degenerate situations: when the points are colinear, the matrix Apq looses a rank. The direct method fails in this case. After some adjustments, Muller's solver can handle this case -- the rotation matrix R ends up having some scale, so I'm "fixing" this by renormalizing. How is everyone else dealing with this?

My other problem with rank-deficient Apq is that I'm drawing the bodies by extracting a rigid transform (using center-of-mass and the rotation matrix R calculated by the shape matching solver). When Apq has linearly dependent rows, R is no longer the correct body orientation -- it aligns with the axis along which the particles are colinear. This isn't correct: imagine a shape with localspace particle positions (1,1),(0,0),(-1,-1) -- in this case, the particles are colinear along an axis which is diagonal relative to the localspace axes.. however, R aligns the localspace axes with the colinear axis.

I'm also having issues when the shape "flips" inside-out (this phenomena occurs with mass-spring systems, when one particle moves to the wrong side of the shape, reversing the clockness of the shape). Muller corrects this by negating the second column of Apq, however this isn't a correct solution in all cases; when it works, it works well (the shape is translated to prevent the vert from crossing through), however when it doesn't work, things tend to blow up (the shape is rotated by 180deg or more). My biggest problem with this is that I have no idea what to search for.. "flipping"? Inversion and inside-out also suck as search terms


I'm hoping someone has some suggestions or experience they could share; it's possible that I've missed something obvious.

[dog]

The problem with the paper is that the methods it proposes for finding the square root of a matrix are hard to make numerically stable.

Use the Denman-Beavers square root operation instead (see the link above). This will only need 8 iterations max to converge with single precision floats.

[raigan2]

Thanks!

Did you implement any sort of "inside-out prevention"? This isn't mentioned in the paper, however it's implemented in their source, and is definitely an issue that arises in practice with low particle counts.

[dog]

We had no need of inside out checks once we implemented a decent square root.

[raigan2]

I don't think it's related to the stability of the square root, I get the same results using two different solvers.

My (very weak) understanding is that if the determinant of Apq is negative, it indicates that the signed volume of the current points has a sign opposite to the signed volume of the original points, indicating that the shape has inverted itself/flipped inside-out.

Muller fixes this by negating one of the columns of Apq, however that's not a correct solution in general. Sadly I can't quite figure out what the correct solution would be..

This may only be common in 2D.

[dog]

Actually it IS related to the solver for the square root. Specifically you can get axis inversions with both the methods you used. We know because we tried them first.

[raigan2]

I see -- I'll definitely try the DB method then; I didn't realize the square root could cause the axes to flip!

[Antonio Martini]

at the following link there is a demo which implements a few deformable algorithms. In my opinion the shape matching approach is the one that behaves the worst:

http://graphics.ethz.ch/~brunoh/defcolstudio.html

some available choices for deformable bodies are:

1. FEM with implicit integration

2. Position Based Dynamics with added volume conservation constraints for tetrahedra.

3. Shape Matching

as usual what method is the best depends on what our objectives are.

cheers,
Antonio

[raigan2]

Luckily I'm not really using it for deformable stuff -- I'm abusing it to simulate rigid bodies in 2D -- "shape matching" is a type of constraint inside my Position Based Dynamics-based simulator (the "alpha" param is 1 and the "beta" is 0).

My view is that using shape-matching in this way is the explicit representation of a constraint which is often implicitly modeled by simply connecting every particle with every other particle using a distance constraint (i.e it lets each particle propagate forces to every other particle). Shape matching scales much better, and converges in a single solver iteration. Which makes it very easy to model a destructible world

Plus, unlike particle+stick networks, it should prevent against things turning inside-out.. as soon as I implement that sqrt algorithm!

[raigan2]

I've finally got around to trying the Denman-Beavers matrix square root as suggested, and I'm afraid I have to report that it is having absolutely no effect on preventing the shape from flipping "inside-out".

Using 20 or even 40 iterations doesn't make a difference; the square root is properly calculated, the problem is elsewhere. My only guess is that we're losing some info in the process since square-roots annihilate sign data.


Similar to my other two methods (polar decomposition and direct/closed-form eigendecomposition), you can prevent flipping in select cases by negating one of the columns of Apq (i.e negate the second column if particles have high varience in x, etc. as I mentioned earlier, and as seen in Muller's source), however this doesn't work in general -- you get the exact same problems as my other methods, where inside-outness will be prevented by rotating the shape 180deg, instead of by a small translation of the particles.

In fact, the DB algorithm seems to perform worse than polar decomposition, as it fails for cases where the original particle positions are colinear (some matrices become singular), while polar decomposition handles these cases.


This is quite frustrating! I'm surprised that no one else has encountered this problem. I can post code and/or demos, but I don't think that my implementation is at fault -- I've arrived at identical behaviour using three completely different implementations!

To me this suggests that it's a problem with the original algorithm, which sucks because it's an incredibly useful tool aside from this small (but critical) problem.

[dog]

The original algorithm works fine for me with Denman Beavers (well, the spirit of it does anyway) - and this is in a very aggressive simulation with extreme squashing. I suspect you may have bugs elsewhere.

I know how frustrating this sort of thing can be, but stick at it. Are you maintaining your original (un-squashed) shape in a consistent way? You can use this to avoid flips. Are you sure you are not getting degenerate calculations somewhere? Just one bad float value can cause havoc.

Hmmm, I'm not being very helpful really; but I think this is one where you'll just have to instrument your code to catch when things go bad, and look at what is causing the problem. At the very least this should let you see how you would want to change the original algorithm as you have interpreted it to avoid the error you are getting.

[raigan2]

Do you have any suggestions on how the original shape data could be used to avoid flips? I'm simply storing them as a list of 2D positions, in a coordinate system formed by the center-of-mass and the worldspace basis vectors (1,0),(0,1).


I'm skeptical that the problem is bug-related, as things works perfectly in all cases, _except_ that all solvers fail to prevent flips. Long chains of bodies, etc. are totally stable and well-behaved.

My test case for flipping is extremely simple: I'm using a shape that is a very wide slightly asymmetrical "V", such as (50,50),(70,60),(150,50). (this is in 2D)

If the middle point crosses the line formed by the other two points, the shape flips around. I can _detect_ that this has occurred by looking at the sign of the determinant of Apq, however I haven't found a way to actually fix the problem in general. Negating a column of Apq only works for certain cases.

Also, about Denman-Beavers: are you using regular matrix inversion at each iteration? Or some sort of pseudo-inverse? If the original shape is colinear (i.e (50,50),(70,50),(150,50)) the iterates become non-invertible/singular/0-determinant for me. Also, I found that even 12 iterations was far from sufficient. Here's my code:

Code: Select all

	var numits = 20;
	var Y = {a:AtA.a,b:AtA.b,c:AtA.c,d:AtA.d};
	var Z = {a:1,b:0,c:0,d:1};
	for(var n = 0; n < numits; n++)
	{
		var invY = Mat22_Inv(Y);
		var invZ = Mat22_Inv(Z);
		
		Mat22_Add_InPlace(invZ,Y);//Y = (Y + Z^-1)
		Mat22_Scale_InPlace(0.5,Y);//Y = (Y + Z^-1)/2
		
		Mat22_Add_InPlace(invY,Z);//Z = (Z + invY)
		Mat22_Scale_InPlace(0.5,Z);//Z = (Z + invY)/2
	}

It definitely produces an accurate square-root, however it doesn't fix the flipping, nor does it handle the colinear cases. In short, it's no better than the polar decomposition I was previously using. Am I missing something? My implementation may be slightly "naive"

My implementation of the algorithm is completely straightforward:
-calculate center-of-mass and Apq, given current and original particle positions
-calculate AtA = Apq^T*Apq
-calculate S = (sqrt(AtA))^-1
-R = Apq*S
-original positions q are transformed to produce target/rigid positions R*q; current positions are moved to target/rigid positions.

The calculation of S is really the only step that seems like it could be problematic. Are you explicitly calculating AtA? Reading through some docs on SVD, it seems that this step can produce some numerical problems.. but I've yet to figure out exactly how to avoid it.

Finally, since my Apq matrix is 2x2, I feel like there should be some direct solution/formula for the SVD. However, so far several searches have turned up only questionable (and ancient) usenet posts on the subject.

Thanks for your help!
raigan

[dog]

Of the top of my head you should do:

1. Integrate particles
2. Find centre of mass of particles
3. Calculate the initial transformation for the new points:
Initialize matrix to all zeroes
Loop over every particle:
Let L = local particle position relative to new centre of mass
Let Q = local particle position for the original shape
Add L * Q.x to the x-axis of your matrix
Add L * Q.y to the y-axis of your matrix etc.
Then scale the final matrix by 1.0/particle count to normalize it
4. Extract the rotation part (and this is from memory, so it might be iffy):
Let our linear transform from above be M
Let M2 = Transpose(M) * M
Let S = MatrixSquareRoot(M2)
Therefore our rotation matrix R = M * MatrixAffineInverse(S)
I then make sure that R is orthonormal (maybe you miss this step?)
5. Calculate the linear transformation and try to do volume preservation:
Transform = Transform of Original Shape * M
Calculate the determinate of this matrix.
If the determinant is sufficiently greater of less than 1 (have a margin)
then scale it according to the paper (some power thing I can't remember)
6. Then compose your final matrix (Transform from 5 scaled by beta + R scaled by 1 - beta).

etc.

My square root is something like:

Given input matrix m
Let m1 = m
Let m2 = identity
Let HalfMat = matrix with 0.5 in the diagonals
for (numIterations) do
Let invm1 = MatrixAffineInverse(m1)
Let invm2 = MatrixAffineInverse(m2)
m1 = (m1 + invm2) * HalfMat
m2 = (m2 + invm1) * HalfMat

return m1

and that's it.

Hope this helps.

(And thanks to my colleague who actually got all this working nicely)

[raigan2]

Thanks -- I'm using the exact same steps (except that I skip (5), as my beta is always 0).

I'm not checking that R is orthonormal, but if I understand correctly this wouldn't help the flipping-inside-out anyway, since in both cases R _is_ orthonormal -- in the flipped case, one of the axes has been reverse/negated.

Have you tried modeling a nearly-flat simplex with your simulator? This is the sort of shape that is very easy to turn inside-out: low particle count, and nearly degenerate/rank-deficient.

Is it possible that your simulator has the same problems as mine, but you never encounter them in practice?

Unless you take explicit steps to prevent flipping, I don't see how you can avoid it, since the problem (as far as I understand it) with flipped cases is that the best-fit transformation _is_ to flip the shape inside-out. Or rather, that you're fitting the flipped configuration, so you're screwed from the start.

Manipulating Apq to make sure the determinant is positive seems to be a solution to this problem, however it's not clear to me how you'd go about altering it to solve the general case. When the flip happens along one of the axes, you can negate one of the basis vectors described by Apq, however this isn't a general solution..

Thanks again for your help and encouragement -- it's unfortunate that for my application I can't rule-out these problematic configurations.

[dog]

It sounds like you are not doing volume preservation. This may be why I can't reproduce your problem.