Setting joint hierarchy rotations to a new bind pose with quaterions

Question

Edit: Adding a more concrete example and some images. Normally you can take the HIERARCHY section of a BVH file, plot the joints with their positions relative their parents and display the bind pose (See first image)

With the 'flipped' skeleton, rendering the same thing results in image2. However, if I then apply the rotations from the 0th frame to the flipped version, I get image 1 again.

What I want to do is rebase all of the frames rotations in the 'flipped' animation so they correspond to the 'normal' skeleton's default joint positions and then use the 'normal' rotation instead of the conditional to render the frames.

The 'flipped' skeletons are generated by Mixamo/Rokoko, but if you import a BVH file into blender and then export it, it will be in the 'normal' skeleton format.

These two lines are a rotation from the Mixamo formatted BVH followed by the rotation generated by Blender relative the unfolded skeleton (both represented in degrees in YZX rotate order)

-0.000000 -0.152100 12.667808 
-0.2075    2.8384   -0.0478 

tl;dr: I want to do this via a script rather than using Blender.

Developing a set of opensource BVH (Biovision Heuristics) animation tools and I've run into a snag I can't seem to solve. I'd like to convert a 'flipped' skeleton to the more common format.

The BVH whitepaper is not too specific and has resulted in 2 different implementations.

I can render both but I'd really like to translate one to the other and do away with case-based rendering.

The following is a small recursive function that builds the world positions for the joints with the more-common and flipped implementations. The only difference is whether we use the cumulative rotation or its conjugate.

(Note: I for clarity I omitted a bunch of if joint.parent is None: logic where the position | rotation is defaulted to 0,0,0)

self and parent are type Joint. w_position is world position of the joint (parent relative position + world position of parent) position is the parent-relative position of the joint. joint.rotation is a parent-relative rotation rotation (the parameter) is the cumulative rotation of the joint's ancestors (world relative rotation)

    def compute_frame_wpos( self, frame_no : int, flipped=False, rotation : glm.quat = glm.quat(glm.vec3(0,0,0)) ):
        '''This function sets the world position of the joints for a given      
          frame.  Recursively calls its children, passing them the cumulative   
          rotations of their ancestors'''            
                          
        if flipped:
            self.w_position = self.parent.w_position + \                        
                              self.position * rotation 
        else:
            self.w_position = self.parent.w_position + \                            
                          self.position * glm.conjugate(rotation)   
                              
        #Recurse                                                                
        child_rotation = self.frames[frame_no].rotation * rotation              
        for child in self.children:
            child.compute_frame_wpos( frame_no, flipped, child_rotation ) 

The initial (parent relative) positions of the two skeletons are different, the flipped one rendering as folded up like a pretzel for the bind pose. What I would like to do is use the 0th frame (the resting pose) in place of the bind pose and correct all the rotations.

I then created an unfolded pose for the flipped skeleton by calling compute_frame_wpos for the 0th (resting) frame, then applying the following recursion to build new local positions from world positions.

    def _make_relative_from_world( self ):
        '''A recursive call to set the joint's parent relative positions from the world positions.'''
          
        self.position = self.w_position - self.parent.w_position

        for child in joint.children:
            child._make_relative_from_world(  )

Finally, I got stuck, trying to fix the rotations. In my calling code, I loop through the frames like this:

Skeleton is an encapsulating class containing the joint hierarchy and an alias list. skeleton contains the original skeleton with its imported parent-relative joint positions. ref_skeleton contains the new unfolded skeleton with the joint positions and world_positions described above.

for frame_no in range(0, self.skeleton.num_frames):
    self.skeleton.get_root().compute_frame_wpos( frame_no, True )
    self._reset_resting_rotation( ref_skel, frame_no )

My naive implementation employs a shorest-arc from the frame's world position for a joint to the world position for the new pose.

    def _reset_frame_rotation( self, reference_skeleton : Skeleton, frame_no : int ):
        '''Recursively change frame's rotations to be relative the reference skeleton'''

        #Get the joint matching this one in the reference_skeleton
        ref_joint = reference_skeleton.joints[ self.name ]

        frame_pos = self.w_position - self.parent.w_position

        self.frames[frame_no].rotation = shortest_arc( frame_pos, ref_joint.position )

        for child in self.children:
            child._reset_frame_rotation( reference_skeleton, frame_no ) 

Pretty sure the shortest_arc implementation is correct but I'll include it, just in case:

def shortest_arc(point_a : glm.vec3, point_b : glm.vec3) -> glm.quat:
    '''Given the vectors a and b return a quaternion representing the shortest arc between'''                     
                                                                                
    result = glm.quat(1.0,0.0,0.0,0.0)                                          
                                                                                
    ab_dot = glm.dot(point_a,point_b)  #dot product                             
    cross = glm.cross(point_a,point_b) #cross product                           
    cross_sqr = glm.dot(cross,cross)  #Squared length of the crossproduct       
                                                                                
    #Test if args have sufficient magnitude                                     
    if ab_dot*ab_dot+cross_sqr != 0.0:                                          
        #Test if the arguments are (anti)parllel                                
        if cross_sqr > 0.0:                                                     
            sqr = math.sqrt(ab_dot * ab_dot + cross_sqr) + ab_dot               
                                                                                
            #Inverted magnitude of the quaternion.                              
            mag = 1.0 / math.sqrt(cross_sqr + sqr * sqr)                        
            result = glm.quat(sqr * mag, cross[0] * mag, cross[1] * mag, cross[2] * mag) 
        elif ab_dot < 0.0: #Test if the angle is > than PI/2 (anti parallel)    
            #The arguments are anti-parallel, we have to choose an axis
            point_c = point_a - point_b
            
            #The length is projected on the XY plane                            
            mag = math.sqrt( point_c[0] * point_c[0] + point_c[1] * point_c[1] )
            
            #F32 magnitude threshhold. Probably not needed.                     
            f32_mag_thresh = 0.0000001
            if mag > f32_mag_thresh:
                result = glm.quat( 0.0, -point_c[1]/mag, point_c[0]/mag, 0.0 )  
            else:
                result = glm.quat( 0.0, 1.0, 0.0, 0.0 )                         
    return result

In writing this, I realize I probably need to get the world rotations of the hierarchy from both the reference and final skeleton and subtract their difference from my final frame rotation and will try that and update if it works.