I have a set of 3D marker data attached on a body that rotates around a semi-constrained joint. I would like to find the center of rotation in three-dimensional coordinate system (x,y,z). Do you know any matlab codes for this calculation? 
Thanks a bunch!



 

"Serkan" wrote in message <j9vf76$e5v$1@newscl01ah.mathworks.com>...
> I have a set of 3D marker data attached on a body that rotates around a semi-constrained joint. I would like to find the center of rotation in three-dimensional coordinate system (x,y,z). Do you know any matlab codes for this calculation? 
> Thanks a bunch!
- - - - - - - - -
  Perhaps we can rephrase the statement of your problem.  You have a set of 3D cartesian coordinates of your "marker data", P, and at a later moment another set, Q, of coordinates of these corresponding points after a period of rotation about some line.  You would like to know what that line of rotation is.  Do I state things correctly?

  It can be shown that any displacement of a rigid 3D body can be accomplished by a rotation about some axis combined with a translation along that axis - a "screw" type displacement.  Presumably in your case there is known to be no such translational component along the axis, only a pure rotation about the axis through some angle.

  It can also be shown that the best fit between the two corresponding point sets in a least squares sense can be achieved by translating each of their centroids to the origin and then rotating one of them via a rotation matrix about the origin to best fit the other.  This latter can be accomplished using matlab's 'svd' (singular value decomposition) function utilizing the Kabsch algorithm.  See for example:

 http://en.wikipedia.org/wiki/Kabsch_algorithm

The equivalent can also be done using the Statistics Toolbox function 'procrustes' appropriately.  See its documentation.

  Your next task is to determine from this rotation matrix the corresponding axis of rotation and the angle of rotation.  For this you can use the technique shown in:

 http://en.wikipedia.org/wiki/Rotation_matrix

using matlab's 'eig' function.

  Finally, from the axis of rotation relative to the centroid of one of the sets, the rotation angle, and the translation between the two centroids, it is a problem in solid geometry to determine where the parallel axis of "pure" rotation is located that requires no translation, and which is presumably what you are seeking.  I'll let you solve that part on your own.  You can probably make good use of matlab's cross product function, 'cross', for this purpose.

Roger Stafford

Thank very much Roger for your kind response. I will take a look at the techniques and resources you referred. I hope that you would be available to answer my questions if I had further along the road. 
Kind Regards...

"Roger Stafford" wrote in message <ja149d$2k4$1@newscl01ah.mathworks.com>...
>
>   It can also be shown that the best fit between the two corresponding point sets in a least squares sense can be achieved by translating each of their centroids to the origin and then rotating one of them via a rotation matrix about the origin to best fit the other.  This latter can be accomplished using matlab's 'svd' (singular value decomposition) function utilizing the Kabsch algorithm.  See for example:
> 
>  http://en.wikipedia.org/wiki/Kabsch_algorithm
> 
> The equivalent can also be done using the Statistics Toolbox function 'procrustes' appropriately.  See its documentation.
> 
>   Your next task is to determine from this rotation matrix the corresponding axis of rotation and the angle of rotation.  For this you can use the technique shown in:
=========

@Serkan, all of the tasks that Roger has described here, including the extraction of the rotation axis, can be done virtually simultaneously using this tool from the File Exchange

http://www.mathworks.com/matlabcentral/fileexchange/26186-absolute-orientation-horns-method

To illustrate this, consider the following data, in which I rotate the columns of A about the axis that is parallel to [1 2 5] and translated by [10,20,30] to produce data set B.

>> A=rand(3,5); 
>> B=AxelRot(A,20,[1 2 5],[ 10,20,30]);
>> 
>> reg=absor(A,B);
>> axisOfrotation=reg.q(2:end)/reg.q(2)

axisOfrotation =

    1.0000
    2.0000
    5.0000

>> axisShift=pinv(eye(3)-reg.R)*reg.t

axisShift =

    3.3333
    6.6667
   -3.3333



A few notes:

1. The AxelRot function is my own utility function, but is available here

http://www.mathworks.com/matlabcentral/fileexchange/30864-3d-rotation-about-shifted-axis


2. The axisShift is not uniquely defined. Both it and [10;20;30] however will satisfy the equation

 x-reg.R*x=reg.t

The following verifies this:

>> (eye(3)-reg.R)*axisShift - reg.t

ans =

  1.0e-014 *

   -0.0888
   -0.1443
   -0.5634

>> (eye(3)-reg.R)*[10;20;30] - reg.t

ans =

  1.0e-013 *

   -0.1377
    0.3242
   -0.1349

"Roger Stafford" wrote in message <ja149d$2k4$1@newscl01ah.mathworks.com>...
>
>   It can be shown that any displacement of a rigid 3D body can be accomplished by a rotation about some axis combined with a translation along that axis - a "screw" type displacement. 
================

I don't see how that can be true, Roger.
If I start with an arbitrary roto-translation

y=R*x+t

then your comment here seems to say that there is a 2nd equivalent decomposition

y=Ra*x+ta

where Ra and ta are a rotation/translation about/along some axis. But the representation of a rototranslation is unique, so that would imply R=Ra, t=ta implying in turn that all rotations are rotations about an axis and all translations are translations along an axis. Obviously, there's something I'm missing here...

 

"Matt J" wrote in message <ja1a92$n0j$1@newscl01ah.mathworks.com>...
> "Roger Stafford" wrote in message <ja149d$2k4$1@newscl01ah.mathworks.com>...
> >   It can be shown that any displacement of a rigid 3D body can be accomplished by a rotation about some axis combined with a translation along that axis - a "screw" type displacement. 
> 
> I don't see how that can be true, Roger.
> If I start with an arbitrary roto-translation
> 
> y=R*x+t
> 
> then your comment here seems to say that there is a 2nd equivalent decomposition
> 
> y=Ra*x+ta
> 
> where Ra and ta are a rotation/translation about/along some axis. But the representation of a rototranslation is unique, so that would imply R=Ra, t=ta implying in turn that all rotations are rotations about an axis and all translations are translations along an axis. Obviously, there's something I'm missing here...
- - - - - - - - - -
  Yes, the statement I made is quite true, Matt.  It is a fundamental theorem in kinematics due to Chasles.  See the site:

 http://en.wikipedia.org/wiki/Chasles'_theorem

It is also stated in the book "Analytical Dynamics" by E. T. Whittaker where it says, "the most general displacement of a rigid body can be obtained by first translating the body, and then rotating it about a line."  Also it says, "This combination of a translation and a rotation round a line parallel to the direction of translation is called a screw."

  Remember that writing the homogeneous expression R*x with R a "rotation matrix" involves only a rotation about a line passing through the coordinate origin.  Chasles' theorem speaks of a rotation about a line that doesn't necessarily pass through the origin, together with a translation parallel to that same line.

Roger Stafford
 

"Roger Stafford" wrote in message <ja1f47$96s$1@newscl01ah.mathworks.com>...
> "Matt J" wrote in message <ja1a92$n0j$1@newscl01ah.mathworks.com>...
> > I don't see how that can be true, Roger.
>   Yes, the statement I made is quite true, Matt.  It is a fundamental theorem in 
- - - - - - - -
  Afterword:  A "screw" type transformation of the Chasles type could be written as:

 y = R*(x-a)+b

where the axis of (homogeneous) rotation matrix R is parallel to the vector b-a, and a and b are points along the (offset) axis of rotation.

Roger Stafford

Thank you guys!

"Serkan" wrote in message <ja3dld$hrd$1@newscl01ah.mathworks.com>...
> Thank you guys!

Hi Serkan,

could you find a suitable solution for finding the rotation center?

It would be kind if you could post any hints about the algo to develop

Thanks