Confusion on rotation of LCS in 3D

[doxley]

I have spent much more time than I'd care to admit trying to get my head around rotation in 3D with scripting. I've searched and read everything I could find including the very helpful description by @DeepSOIC (https://forum.freecadweb.org/viewtopic.php?t=16110). My basic problem seems to revolve about confusion of operations in the LCS vs thinking in the global system.

I've been working on this somehwat contrived example just to try to work it out and I'd appreciate some guidance. I assume an arbitrary line in 3D with an LCS attached to its endpoint (LastParameter). I'd like to rotate the LCS such that the z-axis is along the extension of the line. Secondarily, I'd then like to rotate the x-axis to lie in a plane parallel to the global xy-plane (not considered in example).

The angles below correspond to the projection of the line onto one of the xz, yz, or a plane along line perpendicular to the xy plane. Ideally, these are extrinsic Euler angles, but apparently not. Alternatively, the problem is that I'm rotating in the local coordinate system and don't know how to do the rotations relative to the global coordinate system. I'm also not clear on what order I would need (if there is one for these angles).

In any case, I'd appreciate any suggestions that can let me get back on track. Also, does FreeCAD have a more straightforward way to achieve what I'm trying to do or is it necessary to do it "by hand" as my code is attempting.

Thanks,
--Don

Code: Select all

def rotate_lcs_to_line(line):
    """Align LCS such that z-axis is along line, positive from start to end."""
    line_start = line.Shape.Edges[0].Vertexes[0]
    line_end =line.Shape.Edges[0].Vertexes[1]
    lcs = FreeCAD.ActiveDocument.addObject('PartDesign::CoordinateSystem', "pos2_lcs")
    lcs.Placement = line.Shape.Placement
    lcs.Support = [(line, 'Vertex2')]
    lcs.MapMode = 'Translate'
    del_x = line_end.X - line_start.X
    del_y = line_end.Y - line_start.Y
    del_z = line_end.Z - line_start.Z
    x_angle = np.arccos(del_z / np.sqrt(del_y ** 2 + del_z ** 2)) 
    y_angle = np.arccos(del_z / np.sqrt(del_x ** 2 + del_z ** 2))
    z_angle = np.pi /2 - np.arcsin(del_z / np.sqrt(del_x **2 + del_y **2 + del_z **2))
    # Rotations are not independent, thus any combo fails
    lcs.Placement.Matrix.rotateY(y_angle)
    lcs.Placement.Matrix.rotateX(x_angle)
    lcs.Placement.Matrix.rotateZ(z_angle)
    return
        
start = FreeCAD.Vector(15, 40, 50)
end = FreeCAD.Vector(35, 60, 80)
ls = Part.makeLine(start, end)
s = FreeCAD.ActiveDocument.addObject("Part::Feature", "Line")
s.Shape = ls
rotate_lcs_to_line(s)
    

Code: Select all

OS: Ubuntu 22.04.1 LTS (ubuntu:GNOME/ubuntu)
Word size of FreeCAD: 64-bit
Version: 0.21.
Build type: Release
Branch: unknown
Hash: 1faf86c3da12c1712e1d5ec015721d9aa02f8672
Python 3.10.4, Qt 5.15.3, Coin 4.0.0, Vtk 9.1.0, OCC 7.5.1
Locale: English/United States (en_US)
Installed mods: 
  * Curves 0.5.2
  * Assembly4 0.12.3
  * Manipulator 1.4.9
  * fasteners 0.3.50
  

[onekk]

Hello.

Try to read this.


https://wiki.freecadweb.org/Sandbox:Edw ... #Rotations


Although link could be correct, it is quite old, and probably form 2016 something is changed in FC methods.

A side note, something in your formatting is wrong, so your post is difficult to read.

@edwilliams16 do you have some hints?

I have banged my head, but I know that a simple solution will exist, but my math skills are very limited.

This answer could be very useful for many other people, and if added to your sandbox, will result in a big improvement.

Regards

Carlo D.

[doxley]

Thx, @onekk.

I'm trying to work my way through the sandbox reference - looks very good. I've also edited the post to hopefully fix the formatting problem.
Thx,
--Don

[edwilliams16]

Just placing it directly seems simpler than figuring out the right attachment offset.

Code: Select all

doc = App.ActiveDocument
obj = doc.getObject("BaseFeature")
shp = obj.Shape
n1 = shp.Curve.Direction
zhat = App.Vector(0,0,1)
a =shp.Vertexes[1].Point
rot = App.Rotation(n1.cross(zhat), zhat, n1, 'ZXY')
pl = App.Placement(a, rot)
lcs = doc.getObject('Body').newObject('PartDesign::CoordinateSystem','Local_CS')
lcs.Placement = pl
doc.recompute()

[doxley]

Apologies for slow reply, I've been out of touch for a while.

Thanks @edwilliams16. This is exactly what I was looking for. I think part of my failing was a failure to understand the meaning/use of the cross product and thus my failure to exploit it.

I agree with @onekk that adding something like this explicitly to your sandbox would be valuable.

Thanks,
--Don

[onekk]

Just placing it directly seems simpler than figuring out the right attachment offset.
...

Hello, as usual I'm bothering you.

It will be a good thing that if possible you will make some "graphical explanation" of the geometrical meaning of cross (and dot) product, for us "math impaired" people.

I have not found a "decent" introductory explanation, and from some other posts, I'm not the only one.

I've done some progress in "filling gaps" of my math background, but the "whole" vector math is very difficult to "contestualize" in FreeCAD, as examples are done using some different "terminology".

Kind Regards and sorry for bothering.

Carlo D.

[edwilliams16]

@onekk It's very hard to give intuitive explanations of vector operations that are suitable for all audiences. In my opinion, 3blue1brown (Grant Sanderson) is the master of mathematical exposition.
His take on dot products is
https://www.youtube.com/watch?v=LyGKycYT2v0

and on cross products.
https://www.youtube.com/watch?v=eu6i7WJeinw
https://www.youtube.com/watch?v=BaM7OCEm3G0

Depending on your mathematical background, you may need to refer back to earlier videos of his.

[onekk]

@onekk It's very hard to give intuitive explanations of vector operations that are suitable for all audiences...

Thanks I will try to follow the explanation.

My math background is not very high, mostly statistics and some economic, but not far beyound derivative, series and limits, trigonometry and vector math was only "cited" in my math textbooks, and were not very used in exams tests...

Regards.

Carlo D.

[doxley]

Thanks @edwilliams16. The references by 3blue1brown (Grant Sanderson) were REALLY GOOD. I think I would be well advised to work my way through his full course!

One thought. I've been playing with your suggested solution trying to get a better picture in my head of the Rotation [rot = App.Rotation(n1.cross(zhat), zhat, n1, 'ZXY') ]. While the end result is clear, the movement of the axes takes some thinking. I was wondering if a "movie" in which two objects in arbitrary locations/orientations where rotated as is done with my line and LCS would help visualization. The notion would be to create some intermediate steps rotating the Z axis into place followed by some additional steps rotating the X axis into place. While there is, of course, no intermediate steps in the actual Rotation, I think it might make it easier for the user to grasp what is really happening. Also, rather than using zhat for the Y vector, make it something silly (e.g., (1000, -35, 2.6)..) to emphasize that Y is ignored in this sequence.

Just a thought...
--Don

[edwilliams16]

Thanks @edwilliams16. The references by 3blue1brown (Grant Sanderson) were REALLY GOOD. I think I would be well advised to work my way through his full course!

One thought. I've been playing with your suggested solution trying to get a better picture in my head of the Rotation [rot = App.Rotation(n1.cross(zhat), zhat, n1, 'ZXY') ]. While the end result is clear, the movement of the axes takes some thinking. I was wondering if a "movie" in which two objects in arbitrary locations/orientations where rotated as is done with my line and LCS would help visualization. The notion would be to create some intermediate steps rotating the Z axis into place followed by some additional steps rotating the X axis into place. While there is, of course, no intermediate steps in the actual Rotation, I think it might make it easier for the user to grasp what is really happening. Also, rather than using zhat for the Y vector, make it something silly (e.g., (1000, -35, 2.6)..) to emphasize that Y is ignored in this sequence.

Just a thought...
--Don

You could make a movie using the slerp() method https://wiki.freecadweb.org/Sandbox:Edw ... lation..29