Hi. This might be a stupid request, but I've seen many helpful people on this board so I thought I give it a try anyway. In the paper "Meshless Deformations Based On Shapematching" (see link below) in chapter 3.3 they describe a method to obtain the optimal rotation for a shape. I fully understand the concept but I'm totally lost when it comes to the equation. If someone have the time to "translate" that part into some kind of c++ pseudocode I would be very grateful. Here is the paper. http://graphics.ethz.ch/~mattmuel/publications/MeshlessDeformations_SIG05.pdf Thanks!

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.