diff --git a/aopy/postproc/base.py b/aopy/postproc/base.py index 78700604..3c04a420 100644 --- a/aopy/postproc/base.py +++ b/aopy/postproc/base.py @@ -586,3 +586,4 @@ def get_source_files(preproc_dir, subject, te_id, date): exp_data, exp_metadata = load_preproc_exp_data(preproc_dir, subject, te_id, date) return exp_metadata['source_files'], exp_metadata['source_dir'] + diff --git a/aopy/postproc/centerout.py b/aopy/postproc/centerout.py new file mode 100644 index 00000000..2f14a003 --- /dev/null +++ b/aopy/postproc/centerout.py @@ -0,0 +1,237 @@ +import numpy as np +from scipy.spatial.transform import Rotation as R + +rotations = dict( + yzx = np.array( # names come from rows (optitrack), but screen coords come from columns: + [[0, 1, 0, 0], # x goes into second column (y-coordinate, coming out of screen) + [0, 0, 1, 0], # y goes into third column (z-coordinate, up) + [1, 0, 0, 0], # z goes into first column (x-coordinate, right) + [0, 0, 0, 1]] + ), + zyx = np.array( + [[0, 0, 1, 0], + [0, 1, 0, 0], + [1, 0, 0, 0], + [0, 0, 0, 1]] + ), + xzy = np.array( + [[1, 0, 0, 0], + [0, 0, 1, 0], + [0, 1, 0, 0], + [0, 0, 0, 1]] + ), + xyz = np.identity(4), +) + +exp_rotations = dict( + none = np.identity(4), + about_x_90 = np.array( + [[1, 0, 0, 0], + [0, 0, 1, 0], + [0, 1, 0, 0], + [0, 0, 0, 1]] + ), + about_x_minus_90 = np.array( + [[1, 0, 0, 0], + [0, 0, -1, 0], + [0, 1, 0, 0], + [0, 0, 0, 1]] + ), + oop_xy_45 = np.array( + [[ 0.707, 0.5 , 0.5 , 0.], + [ 0. , 0.707, -0.707, 0.], + [-0.707, 0.5 , 0.5 , 0.], + [ 0., 0. , 0., 1.]] + ), + oop_xy_minus_45 = np.array( + [[ 0.707, 0.5 , -0.5 , 0.], + [ 0. , 0.707, 0.707, 0.], + [ 0.707, -0.5 , 0.5 , 0.], + [ 0., 0. , 0., 1.]] + ), + oop_xy_20 = np.array( + [[ 0.94 , 0.117, 0.321, 0.], + [-0. , 0.94 , -0.342, 0.], + [-0.342, 0.321, 0.883, 0.], + [ 0., 0. , 0., 1.]] + ), + oop_xy_minus_20 = np.array( + [[ 0.94 , 0.117, -0.321, 0.], + [-0. , 0.94 , 0.342, 0.], + [ 0.342, -0.321, 0.883, 0.], + [ 0., 0. , 0., 1.]] + )) + +def _get_mapping(exp_metadata): + ''' + Returns a mapping A that transforms centered hand coordinates to cursor coordinates through c=A*h + + Hand coordinates ordered [hx hy hz] where hx: forward/backward, hy: up/down, hz: right/left + Cursor coordinates ordered [cx cz cy] where cx: right/left, cz: in/out of screen, cy: up/down (same order as exp_data['task']['cursor']) + + Hand coordinates are centered about the hand space origin (i.e. offset already applied) + ''' + + offset = exp_metadata['offset'] + offset_arr = np.array( + [[1, 0, 0, 0], + [0, 1, 0, 0], + [0, 0, 1, 0], + [offset[0], offset[1], offset[2], 1]] + ) + scale = exp_metadata['scale'] + scale_arr = np.array( + [[scale, 0, 0, 0], + [0, scale, 0, 0], + [0, 0, scale, 0], + [0, 0, 0, 1]] + ) + rotation = exp_metadata['rotation'] # optitrack (x: forward/backward, y: up/down, z: right/left) --> screen space (x: right/left, y: forward/backward, z: up/down) + exp_rotation = exp_metadata['exp_rotation'] # out of plane perturbations applied in screen space + + if 'incremental_rotation' not in [feature.decode("utf-8") for feature in exp_metadata['features']]: + perturbation = exp_metadata['pertubation_rotation'] # in plane perturbations applied in screen space, about bmi3d y-axis (in/out of screen) + perturbation_rotation = R.from_euler('y', perturbation, degrees=True).as_matrix() + mapping = np.linalg.multi_dot((scale_arr[:3,:3], rotations[rotation][:3,:3], exp_rotations[exp_rotation][:3,:3], perturbation_rotation)).T # mapping from *centered* hand coords --> cursor coords + else: + mapping = [] + start = exp_metadata['init_rotation_y'] + stop = exp_metadata['final_rotation_y'] + step = exp_metadata['delta_rotation_y'] + for perturbation in np.arange(start, stop+step, step): + perturbation_rotation = R.from_euler('y', perturbation, degrees=True).as_matrix() + mapping.append(np.linalg.multi_dot((scale_arr[:3,:3], rotations[rotation][:3,:3], exp_rotations[exp_rotation][:3,:3], perturbation_rotation)).T) # mapping from *centered* hand coords --> cursor coords + + return mapping + + +def _transform_coords(hand_data, exp_metadata): + ''' + Transforms hand data into mapping used in center out experiment + Args: + hand_traj (3D numpy array): 3D array of hand trajectory data (n_timepoints x 3) : hand data when output from get_kinematics is in BMI3D coordinates. + exp_metadata: exp_metadata from load_preproc_exp_data(preproc_dir, subject, te_id, date) + + Returns: + transformed_hand_traj (3D numpy array) + + ''' + offset = exp_metadata['offset'] + offset_arr = np.array( + [[1, 0, 0, 0], + [0, 1, 0, 0], + [0, 0, 1, 0], + [offset[0], offset[1], offset[2], 1]] + ) + scale = exp_metadata['scale'] + scale_arr = np.array( + [[scale, 0, 0, 0], + [0, scale, 0, 0], + [0, 0, scale, 0], + [0, 0, 0, 1]] + ) + rotation = exp_metadata['rotation'] # optitrack -> screen space +# print(rotation) + exp_rotation = exp_metadata['exp_rotation'] # out of plane perturbations applied in screen space - x - right, y - towards monkey, z - up + perturbation = exp_metadata['pertubation_rotation'] # in plane perturbations applied in screen space about y axis so rotation is inplane + old = np.concatenate((np.reshape(hand_data, -1), [1])) +# print( coords.shape, old.shape, offset_arr.shape, scale_arr.shape, rotations[rotation].shape, exp_rotations[exp_rotation].shape ) + new = np.linalg.multi_dot((old, offset_arr, scale_arr, rotations[rotation], exp_rotations[exp_rotation])) + pertubation_rot = R.from_euler('y', perturbation, degrees=True) + hand_transformed = np.matmul(pertubation_rot.as_matrix(), new[0:3]) +# print(coords, new_coords) + return hand_transformed + +def get_taskspace_and_nullspace(mapping, task='2DCenterOut'): + ''' + Decomposes hand movement into components of movement in the plane of the screen and out of the plane of the screen + Args: + mapping (3 x 3 numpy array): Mapping matrix used in experiment. See :func:`get_mapping` for more details + + Returns: + t_a (3 x 3 numpy array): Matrix that projects hand movement into the task potent cursor space + n_a (3 x 3 numpy array): Matrix that projects hand movement into the null space + ''' + if task == '2DCenterOut': + a = mapping[[0,2], :] # rows are cursor and columns are hand. taking only hy and hz components of optitrack that correspond to up/down and right/left + elif task == '1DTracking': + a = mapping[2,:] + + a_plus = a.T @ np.linalg.pinv(a @ a.T) + t_a = a_plus @ a + n_a = np.eye(3) - t_a + + return t_a, n_a + +def decompose_hand_movements(hand_data, mapping, task='2DCenterOut'): + ''' + Decomposes hand movement into components of movement in the plane of the screen and out of the plane of the screen + Args: + hand_data (2D numpy array): 3D array of hand trajectory per trial (n_timepoints x 3) in optitrack space. X - forward/back, Y - up/down, Z - right/left + mapping (3 x 3 numpy array): Mapping matrix used in experiment. See :func:`_get_mapping` for more details . Rows must be cursor and columns must be hand. + + Returns: + + ''' + t_a, n_a = get_taskspace_and_nullspace(mapping, task) + hand_data = np.array(hand_data) + + hT = np.dot(t_a, hand_data.T) # 3 x n_timepoints + hN = np.dot(n_a, hand_data.T) + + hT_norm = np.mean(np.linalg.norm(hT, axis=0)) # gives the magnitude of movement in task space + hN_norm = np.mean(np.linalg.norm(hN, axis=0)) # gives the magnitude of movement in null space + + return hT_norm, hN_norm + + +def transform_optitrack2hand_coordinates(o_coords): + """ + Transforms coordinates from the Optitrack coordinates (O) to the intuitive hand coordinates for plotting (H). + + Parameters: + o_coords (numpy array): Optitrack data in Optitrack coordinates [Ox, Oy, Oz] as a numpy array (n_timepoints x 3). + + Returns: + numpy array: The transformed coordinates [Hx, Hy, Hz]. + + .. image:: _images/hand_data_coordinates.png + + """ + # Transformation matrix + T = np.array([ + [0, 0, 1], + [0, 1, 0], + [1, 0, 0] + ]) + + # Perform the matrix multiplication + h_coords = o_coords.dot(T) + + return h_coords + + +def transform_bmi3dscreen2cursor_coordinates(b_coords): + """ + Transforms coordinates from the BMI3d screen coordinates (B) to the intuitive cursor coordinates for plotting (H). + Note: Get kinematics functions output hand kinematics in optitrack coordinates & cursor kinematics in bmi3d coordinates + Parameters: + b_coords (numpy array): cursor data in BMI3d screen coordinates [Bx, By, Bz] as a numpy array (n_timepoints x 3). + + Returns: + numpy array: The transformed coordinates [Cx, Cy, Cz]. + + .. image:: _images/cursor_data_coordinates.png + + """ + # Transformation matrix + T = np.array([ + [1, 0, 0], + [0, 0, 1], + [0, 1, 0] + ]) + + # Perform the matrix multiplication + c_coords = b_coords.dot(T) + + return c_coords \ No newline at end of file diff --git a/docs/source/_images/MC_coord_defintion.png b/docs/source/_images/MC_coord_defintion.png new file mode 100644 index 00000000..aaf6343d Binary files /dev/null and b/docs/source/_images/MC_coord_defintion.png differ diff --git a/docs/source/_images/cursor_data_coordinates.png b/docs/source/_images/cursor_data_coordinates.png new file mode 100644 index 00000000..44b9eff8 Binary files /dev/null and b/docs/source/_images/cursor_data_coordinates.png differ diff --git a/docs/source/_images/hand_data_coordinates.png b/docs/source/_images/hand_data_coordinates.png new file mode 100644 index 00000000..464eb572 Binary files /dev/null and b/docs/source/_images/hand_data_coordinates.png differ diff --git a/tests/test_postproc.py b/tests/test_postproc.py index 3537c459..1f669efc 100644 --- a/tests/test_postproc.py +++ b/tests/test_postproc.py @@ -1,4 +1,5 @@ from aopy.postproc import * +from aopy.postproc.centerout import * import aopy import numpy as np import warnings @@ -248,5 +249,28 @@ def test_get_calibrated_eye_data(self): calibrated = get_calibrated_eye_data(eye_data, coefficients) np.testing.assert_array_equal(eye_data, calibrated) +class TestCenterout(unittest.TestCase): + + def test_transform_optitrack2hand_coordinates(self): + o_coords = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) + expected_h_coords = np.array([[3, 2, 1], [6, 5, 4], [9, 8, 7]]) + np.testing.assert_array_equal(transform_optitrack2hand_coordinates(o_coords), expected_h_coords) + + def test_transform_optitrack2hand_coordinates_single_point(self): + o_coords = np.array([[0, 1, 2]]) + expected_h_coords = np.array([[2, 1, 0]]) + np.testing.assert_array_equal(transform_optitrack2hand_coordinates(o_coords), expected_h_coords) + + def test_transform_bmi3dscreen2cursor_coordinates(self): + b_coords = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) + expected_c_coords = np.array([[1, 3, 2], [4, 6, 5], [7, 9, 8]]) + np.testing.assert_array_equal(transform_bmi3dscreen2cursor_coordinates(b_coords), expected_c_coords) + + def test_transform_bmi3dscreen2cursor_coordinates_single_point(self): + b_coords = np.array([[1, 0, 2]]) + expected_c_coords = np.array([[1, 2, 0]]) + np.testing.assert_array_equal(transform_bmi3dscreen2cursor_coordinates(b_coords), expected_c_coords) + + if __name__ == "__main__": unittest.main()