1
$\begingroup$

I'm interested in using Solve[] to find the centre point of rotation (as an $x,y,z$ coordinate point) and axis of rotation (as a vector) from a TransformationFunction obtained by using the Mathematica function FindGeometricTransform. eg:

pts1 = {{-21.365, -1.61273, 2.41973}, {-41.0366, -4.33682, 4.78811},
        {-18.1104, -20.673, 7.53}, {-19.804, 3.79904, 21.6102}};

pts2 = {{-17.9409, -3.2446, -7.46078}, {-35.9907, -7.76684, -14.7927},
        {-14.3971, -22.658, -4.21113}, {-25.7926, -1.61099, 10.8609}};

FindGeometricTransform gives the transformation:

t = FindGeometricTransform[pts2, pts1][[2]]

I posted a similar thread at Mathforum under the topic: How to calculate the 3-D centre point of rotation given the angle of rotation

and got the answer:

Then Solve will give you the fixed point, i.e. the centre of rotation:

 {x, y, z} /. Solve[t[{x, y, z}] == {x, y, z}, {x, y, z}][[1]]

But when I try and use this get the following message:

During evaluation of In[928]:= RowReduce::luc: Result for RowReduce of badly conditioned matrix {{-0.438289,-0.152175,-0.114141,1.80057},{-0.196422,-0.022041,0.0708135,0.32083},{-0.12289,0.142966,0.458519,0.443714}} may contain significant numerical errors. >>

Out[928]= {-5.22037*10^14, 1.37944*10^15, -3.42994*10^14}

the output is also wrong as the answer should be close to my guesstimate of:

approxCenterOfRotation = {0.05418732005730931`,1.3533759077820666`, -0.1590650885642857`}

However, for some reason I cannot reply to the post on Math Forum to get say that this solution didn't work.

This problem can be visualised with the following (NB: I calculated the axisOfRotation vector using a different method [i.e. not using Solve], but am still interested in a method using Solve to calculate this if possible...): 

axisOfRotation(*normalised*) = {0.347494, -0.904472, 0.247341}

viewVector = axisOfRotation*100

Show[
     Graphics3D[{Black, PointSize[0.01], Point[approxCenterOfRotation]}],
     ListPointPlot3D[{pts1, pts2, t /@ pts1},PlotStyle -> PointSize[0.01]],
     Graphics3D[{Blue,Line[{{pts1[[1]], pts1[[2]]}, {pts1[[1]], pts1[[3]]}, {pts1[[1]],
  pts1[[4]]}}]}],
     Graphics3D[{Purple,Line[{{pts2[[1]], pts2[[2]]}, {pts2[[1]], pts2[[3]]}, {pts2[[1]],pts2[[4]]}}]}],
     Graphics3D[{Brown,Line[{{(t /@ pts1)[[1]], (t /@ pts1)[[2]]}, {(t /@ pts1)[[1]], (t /@pts1)[[3]]}, {(t /@ pts1)[[1]], (t /@ pts1)[[4]]}}]}],
     Graphics3D[{Orange, Thick, Dashed,Line[{(approxCenterOfRotation +(axisOfRotation*20)),(approxCenterOfRotation + (axisOfRotation*-20))}]}],
         PlotRange -> All, BoxRatios -> Automatic, AspectRatio -> Automatic,
         Axes -> True, ImageSize -> 700, ViewPoint -> viewVector]

Any suggestions ?

Thanks for your solution Sjoerd C. de Vries, but I'm still having problems using it for a similar case: 

pts1 = {{-4.1703933347009725`, 
    1.4117161073762858`, -1.9926556291047952`}, {-22.189270524165632`,
     6.791481789144898`, -8.802808787459966`}, {-10.82311078603393`, 
-17.22865684593722`, 0.8844690943141205`}, {-9.74367022928643`, 
    6.269145198021164`, 16.59077697323315`}};

pts2 = {{4.170393334702319`, -1.4117161073741804`, 
    1.9926556291087356`}, {-10.018009366114274`, 
    2.942544547202209`, -11.413672837383066`}, {-3.8405873454718273`, 
-19.553858138726234`, 4.578554129848027`}, {-7.427598618275446`, 
    5.792664273078895`, 16.607151377933118`}};

t = FindGeometricTransform[pts2, pts1][[2]]

centreOfRotation = t[{0, 0, 0}]

rotMat = {t[[1, 1, {1, 2, 3}]], t[[1, 2, {1, 2, 3}]], 
  t[[1, 3, {1, 2, 3}]]}

{rotVec} = NullSpace[rotMat - Transpose[rotMat]];
rotVec

viewVector = rotVec*100

Show[
 Graphics3D[{Orange, PointSize[0.02], Point[centreOfRotation]}],

 Graphics3D[{Green, PointSize[0.025], Point[pts1[[4]]]}],
 Graphics3D[{Green, PointSize[0.025], Point[pts2[[4]]]}],

 ListPointPlot3D[{pts1, pts2, t /@ pts1}, 
  PlotStyle -> PointSize[0.01]], 
 Graphics3D[{Blue, 
   Line[{{pts1[[1]], pts1[[2]]}, {pts1[[1]], pts1[[3]]}, {pts1[[1]], 
      pts1[[4]]}}]}], 
 Graphics3D[{Purple, 
   Line[{{pts2[[1]], pts2[[2]]}, {pts2[[1]], pts2[[3]]}, {pts2[[1]], 
      pts2[[4]]}}]}], 
 Graphics3D[{Brown, 
   Line[{{(t /@ pts1)[[1]], (t /@ pts1)[[2]]}, {(t /@ pts1)[[1]], (t /@
          pts1)[[3]]}, {(t /@ pts1)[[1]], (t /@ pts1)[[4]]}}]}],
 Graphics3D[{Orange, Thick, Dashed, 
   Line[{(centreOfRotation + (rotVec*20)), (centreOfRotation + 
(rotVec*-20))}]}], PlotRange -> All, BoxRatios -> Automatic, 
 AspectRatio -> Automatic, Axes -> True, AxesLabel -> {"X", "Y", "Z"},
  ImageSize -> 700, ViewPoint -> viewVector]

The orange point is the point calculated by centreOfRotation = t[{0, 0, 0}], but looking at the graphic, the center point of rotation should lie closer to the two green points. Also, the axis of rotation (rotVec) looks incorrect. Any further suggestions?

Improve this question
$\endgroup$
1
  • $\begingroup$ sorry, there was an extra \ in the code... third line of code from the bottom: (Graphics3D[{Orange, Thick, Dashed,Line[{(approxCenterOfRotation +(axisOfRotation*20)),(approxCenterOfRotation + (axisOfRotation-20))}]}],) should read: ( Graphics3D[{Orange, Thick, Dashed,Line[{(approxCenterOfRotation +(axisOfRotation*20)),(approxCenterOfRotation + (axisOfRotation*-20))}]}],*) $\endgroup$
    – dime
    Jan 8, 2014 at 23:01

2 Answers 2

Reset to default
1
$\begingroup$

A rotation in 3D does not have a centre point, technically -- just an axis. The reason that Solve complains is that there should be a one-parameter family of solutions. But since the TransformationFunction is in terms of approximate reals, it is slightly off. Things that should cancel out exactly and produce a null space of dimension one do not. But they are close, which is what Mathematica is warning about in its messages. (It doesn't know whether they should or should not cancel.)

One way to get around it is to recognize that the rank of your system should be two and use only two (nontrivial) equations of the system. (Update: Some systems have a null space spanned by one of the standard basis vectors. In such cases, the corresponding variable is missing from the system. The updated code solves for variables that are actually present in the system. Also, small numbers are chopped.) 

axisRule = 
 With[{eqns = DeleteCases[Thread@Chop[newt[{x, y, z}] == 1. {x, y, z}], True]},
  First @ Solve[
    eqns[[1 ;; 2]], 
    DeleteDuplicates@Cases[eqns, _Symbol, Infinity][[1 ;; 2]]
  ]]
(*
   {x -> -0.687317 + 1.522 z, y -> 12.3478 - 4.02178 z}
*)

We can use the answer to parametrize the axis:

axis[z_] := Evaluate[{x, y, z} /. axisRule]

And plot it:

Graphics3D[{PointSize[Large],
  Darker@Green, Point[pts1],
  Red, Point[pts2],
  Thick, Blue, Line[{axis[-5], axis[12]}],
  {Thin, Black, Dashed, Line[Transpose[{pts1, pts2}]],
   PointSize[Medium], Point[Flatten[NestList[t, t@pts2, 9], 1]]}
  }]

Mathematica graphics

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thanks Michael E2. Your solution works well for one of my problems. In this problem the second set of points was created by rotating (alone) the first set of points. For the other problem the second set of points was created by rotating AND translating the first set of points. Your solution did not work for this problem. Was this because there was translation as well? I think so. However, I'm not sure how there is no center point of rotation? We can use RotationTransform[theta,w,p], p is the center of rotation. Is there a way to get from the resulting TransformationFunction back to theta,w,p? $\endgroup$
    – dime
    Jan 10, 2014 at 4:42
  • $\begingroup$ simple example: pt1={{0,0,0},{1,0,0},{0,1,0},{0,0,1}}; w={0,1,0};p={-0.25,-0.25,-0.25}; simpt=RotationTransform[0.35,w,p] pt2 = simpt[pt1] Graphics3D[{PointSize[0.02], Darker@Green, Point[pt1], Darker@Red,Point[pt2],PointSize[0.03],Black,Point[p], Dashed,Thick,Orange,Line[{p-w,p+w}], Thin,Green,Line[{{pt1[[1]],pt1[[2]]},{pt1[[1]],pt1[[3]]},{pt1[[1]],pt1[[4]]}}],Thin,Red,Line[{{pt2[[1]],pt2[[2]]},{pt2[[1]],pt2[[3]]},{pt2[[1]], pt2[[4]]}}]}] newt=FindGeometricTransform[pt1, pt2][[2]] then using axisRule: axis[1]={-0.83,1.7*10^15,1}, but how can this be plotted easily? $\endgroup$
    – dime
    Jan 10, 2014 at 5:06
  • $\begingroup$ @dime See if the update helps. Some linear algebra might help with other edge cases. $\endgroup$
    – Michael E2
    Jan 10, 2014 at 6:00
0
$\begingroup$

My feeling is that the transformation has actually a center of rotation close to infinity and that what the method returns is just affected by numerical errors.

To see this you may do the following:

rpts1 = Rationalize[pts1, 10^-20];
rpts2 = Rationalize[pts2, 10^-20];

t = FindGeometricTransform[rpts2, rpts1][[2]];

rt = Rationalize[t, 10^-20];

fixedPoint = {x, y, z} /. Solve[rt[{x, y, z}] == {x, y, z}, {x, y, z}];

fixedPoint // N

(-7.47793*10^(14) -2.10147*10^(15) -1.07537*10^(14)

which are very big numbers.

Actually, rationalizing the two points pts1 and pts2 seems to have no effect on the transformation t, since probably the algorithm is purely numeric. But rationalizing the transformation shows the big values in the fixed point of the transformation.

Improve this answer
$\endgroup$
1
  • $\begingroup$ I don't think this is correct. See my answer and see what you think. $\endgroup$
    – Michael E2
    Jan 9, 2014 at 22:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.