Create coordinate system with vector components

Jimmy He
Jimmy He Member, Employee Posts: 24
10 Comments 5 Likes First Answer First Anniversary
✭✭✭✭

If I know the vector components of the principal coordinate axes, how should I go about creating a coordinate system using Mechanical scripting?

Comments

  • Jimmy He
    Jimmy He Member, Employee Posts: 24
    10 Comments 5 Likes First Answer First Anniversary
    ✭✭✭✭

    Below is a code snippet to achieve this:
    `

    import math
    
    ##########################################################
    # Created by P. Thieffry
    # Last update: 2021/04/07
    def dotProduct(v1,v2):
        return v2[0]*v1[0]+v2[1]*v1[1]+v2[2]*v1[2]
    
    def TransformationMatrix(origCS,targCS):
        # compute transformation matrix to align origCS to targCS
        # All coordinates are expressed in the global CS
        #
        tMat=[[dotProduct(origCS[0],targCS[0]),dotProduct(origCS[0],targCS[1]),dotProduct(origCS[0],targCS[2])],
        [dotProduct(origCS[1],targCS[0]),dotProduct(origCS[1],targCS[1]),dotProduct(origCS[1],targCS[2])],
        [dotProduct(origCS[2],targCS[0]),dotProduct(origCS[2],targCS[1]),dotProduct(origCS[2],targCS[2])]]
        return tMat
    ##########################################################
    
    ##########################################################
    # Modified from: https://learnopencv.com/rotation-matrix-to-euler-angles/
    def rotationMatrixToEulerAngles_ZYX( R ) :
        sy = math.sqrt(R[0][0] * R[0][0] +  R[1][0] * R[1][0])
        singular = sy < 1e-6
        factor = 180./math.pi
        if  not singular :
            x = math.atan2(R[2][1] , R[2][2])*factor
            y = math.atan2(-R[2][0], sy)*factor
            z = math.atan2(R[1][0], R[0][0])*factor
        else :
            x = math.atan2(-R[1][2], R[1][1])*factor
            y = math.atan2(-R[2][0], sy)*factor
            z = 0
        return [z, y, x]
    ##########################################################
    
    def norm( v ):
        return ( v[0]**2 + v[1]**2 + v[2]**2 )**0.5
    
    def unit_vec( v ):
        n = norm( v )
        return [ _/n for _ in v ]
    
    def CreateCS( origin , xyz , name , unit='mm' ):
        # Create coordinate
        testCS = ExtAPI.DataModel.Project.Model.CoordinateSystems.AddCoordinateSystem()
        testCS.OriginDefineBy= CoordinateSystemAlignmentType.Fixed
        testCS.OriginX = Quantity(origin[0],unit)
        testCS.OriginY = Quantity(origin[1],unit)
        testCS.OriginZ = Quantity(origin[2],unit)
    
        # Get rotations
        origCS = [[1,0,0],[0,1,0],[0,0,1]]
        tMat = TransformationMatrix(origCS,xyz)
        angs = rotationMatrixToEulerAngles_ZYX(tMat)
    
        ##########################################################
        # Created by M.H. Pernelle
        # Last update: 2021/06
        # Perform transformation
        testCS.PrimaryAxisDefineBy = CoordinateSystemAlignmentType.GlobalX
        # Rotation Z
        testCS.RotateZ()
        testCS.SetTransformationValue(1,angs[0])
        # Rotation Y
        testCS.RotateY()
        testCS.SetTransformationValue(2,angs[1])
        # Rotation X
        testCS.RotateX()
        testCS.SetTransformationValue(3,angs[2])
        ##########################################################
    
        # Set name
        testCS.Name = name
    
    # Example
    xyz = [[0.5,0.5,-0.70711],[-0.14645,0.85355,0.5],[0.85355,-0.14645,0.5]]
    CreateCS( [1,2,3] , xyz , 'my_cs' , unit='m' )
    

    `

  • Mike.Thompson
    Mike.Thompson Member, Employee Posts: 382
    25 Answers 100 Comments Second Anniversary 25 Likes
    ✭✭✭✭

    @Jimmy He ,
    If you have any time/desire I am hoping to expand this solution to have a function that creates a CS with inputs of XYZ origin and a Z vector definition. X and Y can be arbitrary, but the Z direction should be from a unit vector.

    Thanks for any help with this.

  • Jimmy He
    Jimmy He Member, Employee Posts: 24
    10 Comments 5 Likes First Answer First Anniversary
    ✭✭✭✭
    edited January 2024

    Hi @Mike.Thompson

    Please see the updated solution here. The user can specify the unit of the origin XYZ. Angle unit is detected in the code and appropriate conversion is performed.

    import math
    
    ##########################################################
    # Created by P. Thieffry
    # Last update: 2021/04/07
    def dotProduct(v1,v2):
        return v2[0]*v1[0]+v2[1]*v1[1]+v2[2]*v1[2]
    
    def TransformationMatrix(origCS,targCS):
        # compute transformation matrix to align origCS to targCS
        # All coordinates are expressed in the global CS
        #
        tMat=[[dotProduct(origCS[0],targCS[0]),dotProduct(origCS[0],targCS[1]),dotProduct(origCS[0],targCS[2])],
        [dotProduct(origCS[1],targCS[0]),dotProduct(origCS[1],targCS[1]),dotProduct(origCS[1],targCS[2])],
        [dotProduct(origCS[2],targCS[0]),dotProduct(origCS[2],targCS[1]),dotProduct(origCS[2],targCS[2])]]
        return tMat
    ##########################################################
    
    ##########################################################
    # Modified from: https://learnopencv.com/rotation-matrix-to-euler-angles/
    def rotationMatrixToEulerAngles_ZYX( R ) :
        sy = math.sqrt(R[0][0] * R[0][0] +  R[1][0] * R[1][0])
        singular = sy < 1e-6
        if  not singular :
            x = math.atan2(R[2][1] , R[2][2])
            y = math.atan2(-R[2][0], sy)
            z = math.atan2(R[1][0], R[0][0])
        else :
            x = math.atan2(-R[1][2], R[1][1])
            y = math.atan2(-R[2][0], sy)
            z = 0
        return [z, y, x]
    ##########################################################
    
    def norm( v ):
        return ( v[0]**2 + v[1]**2 + v[2]**2 )**0.5
    
    def unit_vec( v ):
        n = norm( v )
        return [ _/n for _ in v ]
    
    def cross(a, b):
        result = [a[1]*b[2] - a[2]*b[1],
                a[2]*b[0] - a[0]*b[2],
                a[0]*b[1] - a[1]*b[0]]
        return result
    
    def CreateCS( origin , z , name , unit='mm' ):
        angle_unit = DataModel.CurrentUnitFromQuantityName("Angle")
        convert_factor = Quantity(1,'rad').ConvertUnit( angle_unit ).Value
    
        # Create coordinate
        testCS = ExtAPI.DataModel.Project.Model.CoordinateSystems.AddCoordinateSystem()
        testCS.OriginDefineBy= CoordinateSystemAlignmentType.Fixed
        testCS.OriginX = Quantity(origin[0],unit)
        testCS.OriginY = Quantity(origin[1],unit)
        testCS.OriginZ = Quantity(origin[2],unit)
    
        z_unit = unit_vec(z)
        x_unit = unit_vec([ z_unit[1] , -z_unit[0] , 0 ])
        y_unit = unit_vec(cross(z_unit, x_unit))
        xyz = [ x_unit , y_unit , z_unit ]
    
        # Get rotations
        origCS = [[1,0,0],[0,1,0],[0,0,1]]
        tMat = TransformationMatrix(origCS,xyz)
        angs = rotationMatrixToEulerAngles_ZYX(tMat)
    
        ##########################################################
        # Created by M.H. Pernelle
        # Last update: 2021/06
        # Perform transformation
        testCS.PrimaryAxisDefineBy = CoordinateSystemAlignmentType.GlobalX
        # Rotation Z
        testCS.RotateZ()
        testCS.SetTransformationValue(1,angs[0]*convert_factor)
        # Rotation Y
        testCS.RotateY()
        testCS.SetTransformationValue(2,angs[1]*convert_factor)
        # Rotation X
        testCS.RotateX()
        testCS.SetTransformationValue(3,angs[2]*convert_factor)
        ##########################################################
    
        # Set name
        testCS.Name = name
    
    # Example
    z = [1.,1.,0.]
    CreateCS( [0,0,0] , z , 'my_cs' , unit='m' )