Shapematching algorithm

Hi Sweenie,

I'm not going to write any pseudo-code but I'll try to explain what they are doing at each step in that section.

There are two sets of points: one set describing the orginal object (x^0) and one set describing the deformed object (x). Ignoring w_i, eq. (5) describes the squared distance between the points of the original object and the deformed object. I guess the authors figure this is a good way to measure how close/similar the two shapes are. Note that working with the square of a distance is often gives a qualitatively similar result to using the actual distance but it means that you don't have to call nasty square root functions. This could also be interpreted as using "least squares minimization". Now t^0 is the center of mass (CoM) of the original shape and t is the CoM of the deformed shape (eq. (6) gives those equation). Subtracting t^0 from x^0 moves the object so that it's CoM is at the origin. Subtracting t from x does the same. So the CoM's of the two objects are both at the origin and all that remains is to "rotate" the original object so that it matches the deformed object.

The points in these two sets are probably not going to match up perfectly. So the authors try to find the transformation that minimizes the difference between the shapes. So they do a co-ordinate transform to get p_i and q_i (tidies up the notation a little) and get sum_i [ m_i (A q_i - p_i)^2 ]. This is just eq. (5) reworked a little. The problem now is to find the matrix A that minimizes that formula. They do this using calculus.

Unfortunately, things do get complicated at this point. They treat each element of the matrix A as a variable. All other parts of the problem are constant. They differentiate the equation given above with respect to every element in A. This results in the equation

0 = sum_i [ m_i (A q_i - p_i) transpose(q_i) ].

Hopefully, you can see how that arises. Once you have this equation you need to solve for A. I won't go through the steps involve but you will get the equation

A sum_i [ m_i q_i transpose(q_i) ] = sum_i [ m_i p_i transpose(q_i) ].

Each of the sumations represents a matrix. So its convenient to define

A_qq = sum_i [ m_i q_i transpose(q_i) ]
A_pq = sum_i [ m_i p_i transpose(q_i) ]

They do not seem to define A_qq and A_pq explicitly in their article and they should. So our equation becomes

A A_qq = A_pq.

Assuming A_qq has an inverse, the solution is

A = A_pq inverse(A_qq).

This is almost what they have in their article but they've made a type and left of the inverse symbol (actually, it is hard to know if they have since they didn't define A_qq or A_pq).

Since A_qq is symmetric, it is not a rotation and the rotation of the shape must be in A_pq. But there could also be some scaling in A_pq. So the last part of the derviation is where they extract only the part of A_pq that rotates the object and that is their solution to the problem.

I hope that helped, let me know if there is anything you would like me to explain further.


-Josh

Hi Josh.

Well, it brought me a little bit closer to understanding, but I'm afraid I'm still lost.

Does _i define each point or each component of a 3D-vector(x, y & z)?
I suppose it doesn't define a 3dvector since I only thought one could transpose a matrix.

Even though the solution is just infront of me I simply don't know how to build the rotation matrix.
Just to let you know what level my equationskill is at I will show you the type of equations I can solve...
x + 4 = 10
x = 10 - 4
(please don't laugh) :)

So I guess I should stay in the shallow part of the pool.
Anyway, thank you very much for trying to explain.
I'm sure someone else will find it useful though.

ps.

Since I din't want to give up on this that easy, I tried a method that fits my level of math.
I took the angles between the points and CoM of the original shape and the angles between the points and CoM of the deformed shape.
When I added all angle difference for each point and used the averaged angle difference to rotate the original shape.

It worked somewhat well, a maximum of 3-4 degrees from the optimal rotation.
To verify my estimated rotation I did a full 360 turn of the original shape(in 1 degree steps) and logged the sum of all distances between each rotated original point and deformed point to see visually what the optimal rotation should be.

But the body is jittering too much and is very unstable so I guess this amateur method din't work. :)

Hi Sweenie,

sorry about the notation, I didn't want the explanation to become to turgid [smile]

"sum_i" means that this is a summation over the index "i". To be more carefully there really ought to be an upper limit on the index but we can ignore that i think. So "sum_i x_i" means sum up the x's whatever they might be. In the article you're reading the x's are vectors and x_i represents the ith vector in the list of vectors. The way that the authors have formulated their problem, the vectors could have any dimension and it will still work.

You can transpose both a matrix and a vector. Usually, the column vector is the default vector. If you transpose it you get a row vectors and vice versa. Now if you have a pair of vectors, p and q, and suppose that p' is the transpose of p, there are two similar products that we can form

p'q : the inner product (a scalar)
p q' : the outer product (a matrix)

So the product of the two vectors p and q in your article produce matrices because they are outer products.

Hopefully that sheds a little more light on the issue [smile]

Good luck!


-Josh

So, if I want the outer product of p q' would it be like this?

[p.x*q.x p.x*q.y p.x*q.z]
[p.y*q.x p.y*q.y p.y*q.z]
[p.z*q.x p.z*q.y p.z*q.z]

If m_i was a scalar defining the squared distance between each point on the original and deformed shape, would the matrix look like this?

[m*p.x*q.x m*p.x*q.y m*p.x*q.z]
[m*p.y*q.x m*p.y*q.y m*p.y*q.z]
[m*p.z*q.x m*p.z*q.y m*p.z*q.z]

When you put the sum_i before the left bracked, does that mean that i'm supposed to add the matrices together?

I thought that I was supposed to sum all the components of p and q together, but that produced a zero vector for p if the points of the original shape defined a perfect square.

hi Sweenie,

You have the outer product correct. Just to be very careful about the scalar multiplication

Quote:

[m*p.x*q.x m*p.x*q.y m*p.x*q.z] [m*p.y*q.x m*p.y*q.y m*p.y*q.z] [m*p.z*q.x m*p.z*q.y m*p.z*q.z]

This is a scalar, m, times the matrix you obtained from the outer product of p and q.

OK, so the summation notation can be a little confusing because people are sometimes sloppy with how they use it. One thing to remember about summation is that it is just adding things up. It's that simple. The summation symbol is short hand so that we don't have to write out a big list of numbers to add up,

sum(i=1:n) a(i) = a(1) + a(2) + ... + a(n)

here the summation goes from 1 to n in the index i. (I'm changing my notation of the summation so that I hope it'll be clearer).

Now, suppose you find something like

S = sum(i=1:n) p(i) * q(i)

this is equal to

S = p(1)*q(1) + p(2)*q(2) + ... + p(n)*q(n).

Now, I hope, the order of the opertations becomes clearer; You first multiply all those p's and q's together and then add up the result. It doesn't matter if the p's and q's are scalars or matrices or functions, the order of the operations will be like this. If you are confused, write a small part of the series out like I did above and it usually becomes clear.

Good luck!


-Josh

hmm, ok.

But if I write it out as a serie I'll get this...

A_pq = m(1)*p(1)*q(1)' + m(2)*p(2)*q(2)' ...

I read that as...
multiply the scalar m with the vector p then multiply the result with the transposed vector q which produces a 3x3 matrix.
Then do the same for the next set/index and add the resulting matrix to the previous one until all points are processed.
I thought that was what I was doing in the second matrix in the post above.

One thing though, I've learned from some math-sites that when multiplying two 3d-vectors, you can use the DotProduct(innerproduct) which returns a scalar or the CrossProduct(outerproduct) which returns a 3d-vector.
But what is the difference between the outerproduct which returns a vector and the outerproduct which returns a 3x3 matrix?

I'm sure you feel like beating me senseless with a mathsbook right now, but please wait until I got this thing working. ;)

Quote:

A_pq = m(1)*p(1)*q(1)' + m(2)*p(2)*q(2)' ...

OK, the point I was making about scalar multiplication by m is that there is a distinction between

m p_1 q_1 + m p_2 q_2 + ... + m p_n q_n

and

m_1 p_1 q_1 + m_2 p_2 q_2 + ... + m_n p_n q_n

The value of each m_i may be different in the second series. In the first series, all of the terms are multiplied by the same scalar so it doesn't have a subscript because it does change. When you said that you multiplied a matrix by m_i and I only see an m in the matrix, I was concerned that you might have the wrong idea. Now, if your sum is written as

sum(i=1:n) m_i p_i q_i

then the series you expanded and gave in your last post is 100% correct.

Quote:

One thing though, I've learned from some math-sites that when multiplying two 3d-vectors, you can use the DotProduct(innerproduct) which returns a scalar or the CrossProduct(outerproduct) which returns a 3d-vector. But what is the difference between the outerproduct which returns a vector and the outerproduct which returns a 3x3 matrix?

Yep, this is a common source of confusion. The inner product is also called the scalar product and the dot product. It probably has other names as well. However, the outer product is not the same as the cross product. The outer product is also called the dyadic product or the tensor product... annoying isn't it [wink] If you find a website that tells you that the outer product is the same as a cross product, find a different site. I suggest looking here for reliable definitions.

Quote:

I'm sure you feel like beating me senseless with a mathsbook right now, but please wait until I got this thing working. ;)

hehe, no worries Sweenie [smile]

aah! I'm an idiot.
I've just realized that m_i isn't the squared distance between the points.
It's simply the mass of each point which I could define as 1.0 for all points if I want to remove it from the equation.

I misread the line in your first post...

Quote:

Ignoring w_i, eq. (5) describes the squared distance between the points of the original object and the deformed object

Since they said in the paper that m_i = w_i I somehow got the idea that m was defining the squared distance. :/

Now, since I thought that solving A_pq was enough, I used that matrix to transform a vector defined by Orig.Shape CoM and Orig.Shape Point 0.
Using the untransformed vector and the transformed vector I built a rotation matrix based on these two vectors.
Since my function to build this matrix normalized the vectors I thought that the scaling in A_pq wasn't going to interfere.
But it seems it does, because I get very weird results.

Well, I guess I have to do the polar decomposition after all.
Sounds easy enough. *shivers* ;)
[Edit]
In the paper they say they Diagonalize the matrix A'pq A_pq using a bunch of Jacobi rotations... ???
Now, I've found some sourcecode that does a single value decomposition of a transformationsmatrix, it seems to me that this extracts the translation, rotation and scaling into separate matrices.
Would that give me the same result?
Ok, A_pq is a 3x3 matrix, but couldn't I just "upgrade" it to a 4x4 matrix?



Thanks for you help Josh. :)

[Edited by - Sweenie on August 31, 2005 3:01:53 AM]

Quote:

Original post by Sweenie Well, I guess I have to do the polar decomposition after all. Sounds easy enough. *shivers* ;) [Edit] In the paper they say they Diagonalize the matrix A'pq A_pq using a bunch of Jacobi rotations... ??? Now, I've found some sourcecode that does a single value decomposition of a transformationsmatrix, it seems to me that this extracts the translation, rotation and scaling into separate matrices. Would that give me the same result?

Hmm... I didn't think you could use a singular value decomposition in that way but after a little digging around on google it seems you can!

However, I think you'll find that Jacobi rotations will prove simpler. Also search for "Givens rotations". The idea is simply to apply a rotation to the matrix that sets specific elements or sets of elements to zero. So you just go through the matrix setting all the off-diagonal elements to zero, thus diagonlizing it.

Ok, A_pq is a 3x3 matrix, but couldn't I just "upgrade" it to a 4x4 matrix?
quote]

I'm not sure what you mean here. IIRC, the analysis they used made no assumptions about the dimension of the space they were using, so I guess you could apply their result to a 4x4 matrix.

details of this thread aside: it seems like quite an awesome technique for games, i hope to see it in a game/demo soon.

Sweenie: learning more about linear algebra can be a great help, here and with many other problems. i dont know what would be a good place for this on the net, maybe someone else has a pointer?

Quote:

Sweenie: learning more about linear algebra can be a great help, here and with many other problems.

Indeed it would.
I suppose this kind of math is what people learn at college, which I never been to.
I've made several attempts to buy a mathsbook about linear algebra suitable for beginners, but it's seems to exist thousands of different books out there and since most of them are quite expensive I didn't want to buy one that turns out to be too easy or too hard.
I've borrowed a book from a friend once, but it required some pre-skills and thus they lost me at page 3(not the copyright notice).

So what I know today is basically what I managed to read on the net and by browsing through sourcecodes.
For example I know what a transformationmatrix does and how to use it, but I don't know why it works.

It's pretty frustrating to read all these exciting and awesome papers and when you reach the "goodies"(algorithms) all you see is letters and symbols nicely wrapped inside a bunch of brackets and parentheses.

Anyway, with some effort and alot of help from the "symbologist" jjd(Josh), I've almost managed to get the Shapematching algo working.

I'm having a litte problem with the Givens rotations.

Following a tutorial I read I managed to zero one of the matrix elements, but I'm failing with the others.
Let me explain:

Let's say I define my "undecomposed??" matrix like this.

A =
[a b c]
[d e f]
[g h i]

Now I want to zero the g element.
First I compose a new Given rotationmatrix which I will call G1
scalar ca = atan(g/d) ;cosangle
[ 1 0 0 ]
[ 0 cos(ca) sin(ca)]
[ 0 -sin(ca) cos(ca)]

Now multiplying G1 with A (G1*A) I get the A-matrix with element g zeroed.
I've verified the matrix-result with the one in the tutorial and it came out right.

Now I thought that to zero element d I just let ca be equal to atan(d/a).
But that didn't work.

What could I be doing wrong?
[Edit]
Wait! I think I found a pattern...
if a,d & g is a vector x,y,z
Then g/d = z/y and the rotationmatrix is a rotation around the x-axis.
Now if d/a = y/x then the rotationmatrix should be around the z-axis.
[ cos(ca) sin(ca) 0 ]
[ -sin(ca) cos(ca) 0 ]
[ 0 0 1 ]

Now I actually get both g & d zeroed.

Can anyone confirm that I'm going in the right direction?



[Edited by - Sweenie on September 1, 2005 3:30:54 AM]

Quote:

Original post by Sweenie I'm having a litte problem with the Givens rotations. Following a tutorial I read I managed to zero one of the matrix elements, but I'm failing with the others. Let me explain: Let's say I define my "undecomposed??" matrix like this. A = [a b c] [d e f] [g h i] Now I want to zero the g element. First I compose a new Given rotationmatrix which I will call G1 scalar ca = atan(g/d) ;cosangle [ 1 0 0 ] [ 0 cos(ca) sin(ca)] [ 0 -sin(ca) cos(ca)] Now multiplying G1 with A (G1*A) I get the A-matrix with element g zeroed. I've verified the matrix-result with the one in the tutorial and it came out right. Now I thought that to zero element d I just let ca be equal to atan(d/a). But that didn't work. What could I be doing wrong? [Edit] Wait! I think I found a pattern... if a,d & g is a vector x,y,z Then g/d = z/y and the rotationmatrix is a rotation around the x-axis. Now if d/a = y/x then the rotationmatrix should be around the z-axis. [ cos(ca) sin(ca) 0 ] [ -sin(ca) cos(ca) 0 ] [ 0 0 1 ] Now I actually get both g & d zeroed. Can anyone confirm that I'm going in the right direction?

It looks like you've got it! :)

Just a couple of points I'll make. The first is just something you might not have picked up in your readings, but it is pretty important in linear algebra. The order in which matrices are multiplied affects the results, i.e. AB != BA (in general). So when you want to multiply a maxtrix A by a matrix B you could mean either AB (post-multiplication by B) or BA (pre-multiplication by B). In technical-speak, matrices do not "commute" under multiplication.

The other point is that you can exchange the 3 trig functions in your calculation with 1 square root. Your rotation matrix is something like

Q =
[1 0 0]
[0 c s]
[0 -s c]

where c = cos(theta) and s = sin(theta). But you only need the values c and s, not necessarily theta. If I (pre) multiply your original A matrix with Q I get

A(3,1) = g * c - d * s

where I'm using fortran indexing (i.e. the first index is 1, not 0). We want A(3,1) to equal zero, so

g * c = d * s,

or

t * c = s

where t = g / d (assuming d != 0). What else do we know about the problem? Well, our transformation, although it may not seem obvious, are dealing with right-angled triangles, so Pythagoras theorem applies. Using that,

cos(theta) = d / sqrt(d*d + g*g)
sin(theta) = g / sqrt(d*d + g*g)

where d is the "length" of the adjacent side of the triangle and g is the "length" of the opposite side. So, we don't need to explicitly determine theta. Here's some matlab code that performs a Givens rotation on a vector

[a,b] --> [*,0]




-Josh

Ah, great! :)
Speed improvements are always welcome.

Just a couple of minutes ago, I struck gold!
I found an italian site with LOADS of math functions nicely explained and presented in Ms Excel addons.
This makes everything so much easier since I can try out my matrices in Excel and verify that my values are correct.
I just downloaded and Excel Addon that contains functions for Matrices, including Jacobi rotations amongst others, and it's totally free. :)

Here is the link... (And it's all in english, except the introtext on the mainpage)
Foxes XNumbers 4.4 & Matrix 1.9

How silly of me! I should have pointed you in the direction of Octave. It is a free matlab clone. It is pretty good for the odd matrix calculation and comes with numerous functions that you'll find useful.


-Josh

Yes!

I finally made it! :)

See for yourself...
---> Demo
Run the exe and click anywhere on screen to attract the body towards the mousepointer.
Press o, m or d to toggle rendering of the original, matched and deformed shape.

@jjd.
Thanks again for all your help.

@Anonymous Poster
Thanks for sharing your sourcecode.
It really helped me puzzle together all the pieces.

Quote:

If you get it working, can you tell me what I'm doing wrong?

To tell you the truth, I have no idea.
In fact, I actually used your sqrt routine except for the jacobian rotation part which is based on the code from the MATRIX lib mentiond a few posts above.
Maybe it's the last step where you use the beta-value and the A_qq matrix.
I never do that part, since I only wanted to simulate rigid bodies so far.

If you download the ziparchive above, you'll find the sourcecode for the demo.
It's written in Blitzmax(BASIC), but it shouldn't be any problem for a C++ coder.(I guess BASIC is like pseudocode for c++ people ;) )