diff --git a/examples/human_pose_3d/cameras.py b/examples/human_pose_3d/cameras.py new file mode 100755 index 000000000..c80500d5e --- /dev/null +++ b/examples/human_pose_3d/cameras.py @@ -0,0 +1,210 @@ +"""Utilities to deal with the cameras of human3.6m""" + +from xml.dom import minidom +import math +import numpy as np + +CAMERA_ID_TO_NAME = { + 1: "54138969", + 2: "55011271", + 3: "58860488", + 4: "60457274", +} + + +def eulerAnglesToRotationMatrix(theta): + R_x = np.array([[1, 0, 0], + + [0, math.cos(theta[0]), -math.sin(theta[0])], + + [0, math.sin(theta[0]), math.cos(theta[0])] + + ]) + R_y = np.array([[math.cos(theta[1]), 0, math.sin(theta[1])], + + [0, 1, 0], + + [-math.sin(theta[1]), 0, math.cos(theta[1])] + + ]) + + R_z = np.array([[math.cos(theta[2]), -math.sin(theta[2]), 0], + + [math.sin(theta[2]), math.cos(theta[2]), 0], + + [0, 0, 1] + + ]) + + R = np.dot(R_x, np.dot(R_y, R_z)) + return R + + +def project_point_radial( P, R, T, f, c, k, p ): + """ + Project points from 3d to 2d using camera parameters + including radial and tangential distortion + + Args + P: Nx3 points in world coordinates + R: 3x3 Camera rotation matrix + T: 3x1 Camera translation parameters + f: (scalar) Camera focal length + c: 2x1 Camera center + k: 3x1 Camera radial distortion coefficients + p: 2x1 Camera tangential distortion coefficients + Returns + Proj: Nx2 points in pixel space + D: 1xN depth of each point in camera space + radial: 1xN radial distortion per point + tan: 1xN tangential distortion per point + r2: 1xN squared radius of the projected points before distortion + """ + + # P is a matrix of 3-dimensional points + assert len(P.shape) == 2 + assert P.shape[1] == 3 + + # print(f"\n****************In camera loop {P[:32].shape} {P[:32]}*****************\n") + + N = P.shape[0] + X = R.dot( P.T - T ) # rotate and translate + XX = X[:2,:] / X[2,:] + r2 = XX[0,:]**2 + XX[1,:]**2 + + radial = 1 + np.einsum( 'ij,ij->j', np.tile(k,(1, N)), np.array([r2, r2**2, r2**3]) ) + tan = p[0]*XX[1,:] + p[1]*XX[0,:] + + XXX = XX * np.tile(radial+tan,(2,1)) + np.outer(np.array([p[1], p[0]]).reshape(-1), r2 ) + Proj = (f * XXX) + c + Proj = Proj.T + + D = X[2,] + + return Proj, D, radial, tan, r2 + + +def world_to_camera_frame(P, R, T): + """ + Convert points from world to camera coordinates + + Args + P: Nx3 3d points in world coordinates + R: 3x3 Camera rotation matrix + T: 3x1 Camera translation parameters + Returns + X_cam: Nx3 3d points in camera coordinates + """ + + assert len(P.shape) == 2 + assert P.shape[1] == 3 + + X_cam = R.dot( P.T - T ) # rotate and translate + + + return X_cam.T + + +def camera_to_world_frame(P, R, T): + """Inverse of world_to_camera_frame + + Args + P: Nx3 points in camera coordinates + R: 3x3 Camera rotation matrix + T: 3x1 Camera translation parameters + Returns + X_cam: Nx3 points in world coordinates + """ + + assert len(P.shape) == 2 + assert P.shape[1] == 3 + # print(f"\n P shape: {P.shape}") + # print(f"\n R in camera_to_world_frame shape: {R.shape} {type(R)}") + # print(f"\n T shape: {T.shape}") + # print(f"\n P.T shape: {(P.T).shape}") + # print(f"\n R.T shape: {(R.T).shape}") + # print(f"\n R.T.dot( P.T ) shape: {(R.T.dot( P.T )).shape}") + X_cam = R.T.dot( P.T ) + T # rotate and translate + + return X_cam.T + + +def load_camera_params(w0, subject, camera): + """Load h36m camera parameters + + Args + w0: 300-long array read from XML metadata + subect: int subject id + camera: int camera id + Returns + R: 3x3 Camera rotation matrix + T: 3x1 Camera translation parameters + f: (scalar) Camera focal length + c: 2x1 Camera center + k: 3x1 Camera radial distortion coefficients + p: 2x1 Camera tangential distortion coefficients + name: String with camera id + """ + + # Get the 15 numbers for this subject and camera + w1 = np.zeros(15) + start = 6 * ((camera-1)*11 + (subject-1)) + w1[:6] = w0[start:start+6] + w1[6:] = w0[(265+(camera-1)*9 - 1): (264+camera*9)] + + def rotationMatrix(r): + R1, R2, R3 = [np.zeros((3, 3)) for _ in range(3)] + + # [1 0 0; 0 cos(obj.Params(1)) -sin(obj.Params(1)); 0 sin(obj.Params(1)) cos(obj.Params(1))] + R1[0:] = [1, 0, 0] + R1[1:] = [0, np.cos(r[0]), -np.sin(r[0])] + R1[2:] = [0, np.sin(r[0]), np.cos(r[0])] + + # [cos(obj.Params(2)) 0 sin(obj.Params(2)); 0 1 0; -sin(obj.Params(2)) 0 cos(obj.Params(2))] + R2[0:] = [ np.cos(r[1]), 0, np.sin(r[1])] + R2[1:] = [0, 1, 0] + R2[2:] = [-np.sin(r[1]), 0, np.cos(r[1])] + + # [cos(obj.Params(3)) -sin(obj.Params(3)) 0; sin(obj.Params(3)) cos(obj.Params(3)) 0; 0 0 1];% + R3[0:] = [np.cos(r[2]), -np.sin(r[2]), 0] + R3[1:] = [np.sin(r[2]), np.cos(r[2]), 0] + R3[2:] = [0, 0, 1] + + return (R1.dot(R2).dot(R3)) + + R = rotationMatrix(w1) + # print(f"\n R in load_camera_params: {R} {type(R)}") + T = w1[3:6][:, np.newaxis] + f = w1[6:8][:, np.newaxis] + c = w1[8:10][:, np.newaxis] + k = w1[10:13][:, np.newaxis] + p = w1[13:15][:, np.newaxis] + name = CAMERA_ID_TO_NAME[camera] + + return R, T, f, c, k, p, name + + +def load_cameras(bpath, subjects=[1,5,6,7,8,9,11]): + """Loads the cameras of h36m + + Args + bpath: path to xml file with h36m camera data + subjects: List of ints representing the subject IDs for which cameras are requested + Returns + rcams: dictionary of 4 tuples per subject ID containing its camera parameters for the 4 h36m cams + """ + rcams = {} + + xmldoc = minidom.parse(bpath) + string_of_numbers = xmldoc.getElementsByTagName('w0')[0].firstChild.data[1:-1] + + # Parse into floats + w0 = np.array(list(map(float, string_of_numbers.split(" ")))) + + assert len(w0) == 300 + + for s in subjects: + for c in range(4): # There are 4 cameras in human3.6m + rcams[(s, c+1)] = load_camera_params(w0, s, c+1) + + return rcams diff --git a/examples/human_pose_3d/data_utils.py b/examples/human_pose_3d/data_utils.py new file mode 100755 index 000000000..d519970a3 --- /dev/null +++ b/examples/human_pose_3d/data_utils.py @@ -0,0 +1,485 @@ +"""Utility functions for dealing with human3.6m data.""" + +import copy +import os +import glob + +import numpy as np +import cdflib + +import cameras + +# Human3.6m IDs for training and testing +TRAIN_SUBJECTS = [1, 5, 6, 7, 8] +TEST_SUBJECTS = [9, 11] + +coco_part_labels = [ + 'nose', 'eye_l', 'eye_r', 'ear_l', 'ear_r', + 'sho_l', 'sho_r', 'elb_l', 'elb_r', 'wri_l', 'wri_r', + 'hip_l', 'hip_r', 'kne_l', 'kne_r', 'ank_l', 'ank_r' +] + +# Joints in H3.6M -- data has 32 joints, but only 17 that move; these are the indices. +H36M_NAMES = [''] * 32 +H36M_NAMES[0] = 'Hip' +H36M_NAMES[1] = 'RHip' +H36M_NAMES[2] = 'RKnee' +H36M_NAMES[3] = 'RFoot' +H36M_NAMES[6] = 'LHip' +H36M_NAMES[7] = 'LKnee' +H36M_NAMES[8] = 'LFoot' +H36M_NAMES[12] = 'Spine' +H36M_NAMES[13] = 'Thorax' +H36M_NAMES[14] = 'Neck/Nose' +H36M_NAMES[15] = 'Head' +H36M_NAMES[17] = 'LShoulder' +H36M_NAMES[18] = 'LElbow' +H36M_NAMES[19] = 'LWrist' +H36M_NAMES[25] = 'RShoulder' +H36M_NAMES[26] = 'RElbow' +H36M_NAMES[27] = 'RWrist' + +# Joints in COCO, 2D poses from HigherHRNet --> data has 17 joints; these are the indices. Hip,12,14,16,11,13,15,Spine,Thorax,0,Head,5,7,9,6,8,10 +# to make compatible with Human3.6M, Nose -> Neck/Nose; Ankle -> Foot +COCO_NAMES = [''] * 17 +COCO_NAMES[0] = 'Head' # Nose renamed as head +COCO_NAMES[1] = 'Thorax' +COCO_NAMES[2] = 'Spine' +COCO_NAMES[4] = 'Hip' +COCO_NAMES[5] = 'LShoulder' +COCO_NAMES[6] = 'RShoulder' +COCO_NAMES[7] = 'LElbow' +COCO_NAMES[8] = 'RElbow' +COCO_NAMES[9] = 'LWrist' +COCO_NAMES[10] = 'RWrist' +COCO_NAMES[11] = 'LHip' +COCO_NAMES[12] = 'RHip' +COCO_NAMES[13] = 'LKnee' +COCO_NAMES[14] = 'RKnee' +COCO_NAMES[15] = 'LFoot' +COCO_NAMES[16] = 'RFoot' + + +def filter_moving_joints_3d(poses3d): + """ + Selects 16 moving joints (Neck/Nose excluded) from 32 predicted joints in 3d + + Args + poses3d: Nx96 points in camera coordinates + Returns + p3d: Nx48 points (moving joints) + """ + reshaped = np.reshape(poses3d, [poses3d.shape[0], 32, 3]) + idx = [0,1,2,3,6,7,8,12,13,15,17,18,19,25,26,27] + p3d = reshaped[:,idx,:] + p3d = p3d.reshape(poses3d.shape[0],-1) + return p3d + + +def preprocess_2d_data(poses_2d): + """Preprocesses 2d detections by creating some extra joints + and converting COCO joint order to H36M. + + Args + poses_2d: list of 2d detections obtained from HigherHRNet + Returns + poses: nx32 np array with 2d poses + """ + # Permutation that goes from COCO detections to H36M ordering. + COCO_TO_GT_PERM = np.array([COCO_NAMES.index(h) for h in H36M_NAMES if h != '' and h in COCO_NAMES]) + assert np.all(COCO_TO_GT_PERM == np.array([4, 12, 14, 16, 11, 13, 15, 2, 1, 0, 5, 7, 9, 6, 8, 10])) + + poses_2d = np.array(poses_2d) + poses_2d = poses_2d[:, :, :2] # Take x and y coord from the HigherHRNet output + poses_2d = np.reshape(poses_2d, + (poses_2d.shape[0], -1)) # reshape to make it compatible with input the model expects + + # make Thorax, mid-point of shoulders i.e. COCO_NAMES[5] & COCO_NAMES[6] + poses_2d[:, 2:4] = (poses_2d[:, 10:12] + poses_2d[:, 12:14]) / 2 + + # make Hip, mid-point of hips i.e. COCO_NAMES[11] & COCO_NAMES[12] + poses_2d[:, 8:10] = (poses_2d[:, 22:24] + poses_2d[:, 24:26]) / 2 + + # make Spine, mid-point of thorax and hip i.e. COCO_NAMES[1] & COCO_NAMES[4] + poses_2d[:, 4:6] = (poses_2d[:, 2:4] + poses_2d[:, 8:10]) / 2 + + # Reshape into (n, 17, 2) matrix + poses_2d = np.reshape(poses_2d, [poses_2d.shape[0], 17, 2]) + + # Permute the loaded data to make it compatible with H36M + poses = poses_2d[:, COCO_TO_GT_PERM, :] + + # Reshape back into n x (32*2) matrix + poses = np.reshape(poses, [poses.shape[0], -1]) + return poses + + +def load_joints_2d(joints): + # Permutation that goes from COCO detections to H36M ordering. + COCO_TO_GT_PERM = np.array([COCO_NAMES.index(h) for h in H36M_NAMES if h != '' and h in COCO_NAMES]) + assert np.all(COCO_TO_GT_PERM == np.array([4, 12, 14, 16, 11, 13, 15, 2, 1, 0, 5, 7, 9, 6, 8, 10])) + + poses_2d = np.array(joints) + poses_2d = np.reshape(poses_2d, (poses_2d.shape[0], -1)) + + # make Thorax, mid-point of shoulders + poses_2d[:, 2:4] = (poses_2d[:, 10:12] + poses_2d[:, 12:14]) / 2 + + # make Hip, mid-point of hips + poses_2d[:, 8:10] = (poses_2d[:, 22:24] + poses_2d[:, 24:26]) / 2 + + # make Spine, mid-point of thorax and hip + poses_2d[:, 4:6] = (poses_2d[:, 2:4] + poses_2d[:, 8:10]) / 2 + + # Reshape into (n, 17, 2) matrix + poses_2d = np.reshape(poses_2d, [poses_2d.shape[0], 17, 2]) + + # Permute the loaded data to make it compatible with H36M + poses = poses_2d[:, COCO_TO_GT_PERM, :] + + # Reshape into nx32 matrix + poses = np.reshape(poses, [poses.shape[0], -1]) + return poses + + +def load_params(): + """Loads normalization statistics: mean and stdev, dimensions used and ignored from npy files + + Returns + data_mean: nxd np array with the mean of the data + data_std: nxd np array with the standard deviation of the data + dim_to_use: nxd np array of dimensions used in the model + dim_to_ignore: nxd np array of dimensions not used in the model + """ + path_prefix = os.path.dirname(os.path.abspath(__file__)) + data_mean_2d = np.load(os.path.join(path_prefix, 'files/data_mean_2d.npy')) + data_std_2d = np.load(os.path.join(path_prefix, 'files/data_std_2d.npy')) + dim_to_use_2d = np.load(os.path.join(path_prefix, 'files/dim_to_use_2d.npy')) + dim_to_ignore_2d = np.load(os.path.join(path_prefix, 'files/dim_to_ignore_2d.npy')) + data_mean_3d = np.load(os.path.join(path_prefix, 'files/data_mean_3d.npy')) + data_std_3d = np.load(os.path.join(path_prefix, 'files/data_std_3d.npy')) + dim_to_use_3d = np.load(os.path.join(path_prefix, 'files/dim_to_use_3d.npy')) + dim_to_ignore_3d = np.load(os.path.join(path_prefix, 'files/dim_to_ignore_3d.npy')) + + return data_mean_2d, data_std_2d, dim_to_use_2d, dim_to_ignore_2d, data_mean_3d, data_std_3d, dim_to_use_3d, dim_to_ignore_3d + + +def normalize_data(data, data_mean, data_std, dim_to_use): + """Normalizes a dictionary of poses + + Args + data: dictionary where values are + data_mean: np vector with the mean of the data + data_std: np vector with the standard deviation of the data + dim_to_use: list of dimensions to keep in the data + Returns + data_out: dictionary with same keys as data, but values have been normalized + """ + data_out = {} + + for key in data.keys(): + data[key] = data[key][:, dim_to_use] + mu = data_mean[dim_to_use] + stddev = data_std[dim_to_use] + data_out[key] = np.divide((data[key] - mu), stddev) + + return data_out + + +def unNormalizeData(normalized_data, data_mean, data_std, dimensions_to_ignore): + """Un-normalizes a matrix whose mean has been substracted and that has been divided by + standard deviation. Some dimensions might also be missing + + Args + normalized_data: nxd matrix to unnormalize + data_mean: nxd np array with the mean of the data + data_std: nxd np array with the standard deviation of the data + dimensions_to_ignore: nxd np array of dimensions that were removed from the original data + Returns + orig_data: the input normalized_data, but unnormalized + """ + T = normalized_data.shape[0] # Batch size + D = data_mean.shape[0] # Dimensionality; 2d data: 64; 3d data: 96 (32 joints) + + orig_data = np.zeros((T, D), dtype=np.float32) + dimensions_to_use = np.array([dim for dim in range(D) + if dim not in dimensions_to_ignore]) + + orig_data[:, dimensions_to_use] = normalized_data + + stdMat = data_std.reshape((1, D)) + stdMat = np.repeat(stdMat, T, axis=0) + meanMat = data_mean.reshape((1, D)) + meanMat = np.repeat(meanMat, T, axis=0) + orig_data = np.multiply(orig_data, stdMat) + meanMat + return orig_data + + +def normalization_stats(complete_data, dim, predict_14=False): + """Computes normalization statistics: mean and stdev, dimensions used and ignored + + Args + complete_data: nxd np array with poses + dim. integer={2,3} dimensionality of the data + predict_14. boolean. Whether to use only 14 joints + Returns + data_mean: np vector with the mean of the data + data_std: np vector with the standard deviation of the data + dimensions_to_ignore: list of dimensions not used in the model + dimensions_to_use: list of dimensions used in the model + """ + if not dim in [2, 3]: + raise ValueError('dim must be 2 or 3') + + data_mean = np.mean(complete_data, axis=0) + data_std = np.std(complete_data, axis=0) + + # Encodes which 17 (or 14) 2d-3d pairs we are predicting + dimensions_to_ignore = [] + if dim == 2: + dimensions_to_use = np.where(np.array([x != '' and x != 'Neck/Nose' for x in H36M_NAMES]))[0] + dimensions_to_use = np.sort(np.hstack((dimensions_to_use * 2, dimensions_to_use * 2 + 1))) + dimensions_to_ignore = np.delete(np.arange(len(H36M_NAMES) * 2), dimensions_to_use) + else: # dim == 3 + dimensions_to_use = np.where(np.array([x != '' for x in H36M_NAMES]))[0] + dimensions_to_use = np.delete(dimensions_to_use, [0, 7, 9] if predict_14 else 0) + + dimensions_to_use = np.sort(np.hstack((dimensions_to_use * 3, + dimensions_to_use * 3 + 1, + dimensions_to_use * 3 + 2))) + dimensions_to_ignore = np.delete(np.arange(len(H36M_NAMES) * 3), dimensions_to_use) + + return data_mean, data_std, dimensions_to_ignore, dimensions_to_use + + +def define_actions(action): + """Given an action string, returns a list of corresponding actions. + + Args + action: String. either "all" or one of the h36m actions + Returns + actions: List of strings. Actions to use. + Raises + ValueError: if the action is not a valid action in Human 3.6M + """ + actions = ["Directions", "Discussion", "Eating", "Greeting", + "Phoning", "Photo", "Posing", "Purchases", + "Sitting", "SittingDown", "Smoking", "Waiting", + "WalkDog", "Walking", "WalkTogether"] + + if action == "All" or action == "all": + return actions + + if not action in actions: + raise ValueError("Unrecognized action: %s" % action) + + return [action] + + +def load_data(bpath, subjects, actions, dim=3): + """Loads 2d ground truth from disk, and puts it in an easy-to-acess dictionary + + Args + bpath: String. Path where to load the data from + subjects: List of integers. Subjects whose data will be loaded + actions: List of strings. The actions to load + dim: Integer={2,3}. Load 2 or 3-dimensional data + Returns: + data: Dictionary with keys k=(subject, action, seqname) + values v=(nx(32*2) matrix of 2d ground truth) + There will be 2 entries per subject/action if loading 3d data + There will be 8 entries per subject/action if loading 2d data + """ + + if not dim in [2, 3]: + raise ValueError('dim must be 2 or 3') + + data = {} + + for subj in subjects: + for action in actions: + dpath = os.path.join(bpath, 'S{0}'.format(subj), 'MyPoseFeatures/D{0}_Positions'.format(dim), + '{0}*.cdf'.format(action)) + # #print( dpath ) + fnames = glob.glob(dpath) + loaded_seqs = 0 + + for fname in fnames: + seqname = os.path.basename(fname) + + # This rule makes sure SittingDown is not loaded when Sitting is requested + if action == "Sitting" and seqname.startswith("SittingDown"): + continue + + # This rule makes sure that WalkDog and WalkTogeter are not loaded when + # Walking is requested. + if seqname.startswith(action): + # #print( fname ) + loaded_seqs = loaded_seqs + 1 + + cdf_file = cdflib.CDF(fname) + poses = cdf_file.varget("Pose").squeeze() + cdf_file.close() + + data[(subj, action, seqname)] = poses + + if dim == 2: + assert loaded_seqs == 8, "Expecting 8 sequences, found {0} instead".format(loaded_seqs) + else: + assert loaded_seqs == 2, "Expecting 2 sequences, found {0} instead".format(loaded_seqs) + + return data + + +def transform_world_to_camera(poses_set, cams, ncams=4): + """Project 3d poses from world coordinate to camera coordinate system + + Args + poses_set: dictionary with 3d poses + cams: dictionary with cameras + ncams: number of cameras per subject + Return: + t3d_camera: dictionary with 3d poses in camera coordinate + """ + t3d_camera = {} + for t3dk in sorted(poses_set.keys()): + + subj, action, seqname = t3dk + t3d_world = poses_set[t3dk] + + for c in range(ncams): + R, T, _, _, _, _, name = cams[(subj, c + 1)] + camera_coord = cameras.world_to_camera_frame(np.reshape(t3d_world, [-1, 3]), R, T) + camera_coord = np.reshape(camera_coord, [-1, len(H36M_NAMES) * 3]) + + sname = seqname[:-3] + name + ".h5" # e.g.: Waiting 1.58860488.h5 + t3d_camera[(subj, action, sname)] = camera_coord + + return t3d_camera + + +def read_3d_data(actions, data_dir, camera_frame, rcams, predict_14=False): + """Loads 3d poses, zero-centres and normalizes them + + Args + actions: list of strings. Actions to load + data_dir: string. Directory where the data can be loaded from + camera_frame: boolean. Whether to convert the data to camera coordinates + rcams: dictionary with camera parameters + predict_14: boolean. Whether to predict only 14 joints + Returns + train_set: dictionary with loaded 3d poses for training + test_set: dictionary with loaded 3d poses for testing + data_mean: vector with the mean of the 3d training data + data_std: vector with the standard deviation of the 3d training data + dim_to_ignore: list with the dimensions to not predict + dim_to_use: list with the dimensions to predict + train_root_positions: dictionary with the 3d positions of the root in train + test_root_positions: dictionary with the 3d positions of the root in test + """ + # Load 3d data + train_set = load_data(data_dir, TRAIN_SUBJECTS, actions, dim=3) + test_set = load_data(data_dir, TEST_SUBJECTS, actions, dim=3) + + if camera_frame: + train_set = transform_world_to_camera(train_set, rcams) + test_set = transform_world_to_camera(test_set, rcams) + + # Apply 3d post-processing (centering around root) + train_set, train_root_positions = postprocess_3d(train_set) + test_set, test_root_positions = postprocess_3d(test_set) + + # Compute normalization statistics + complete_train = copy.deepcopy(np.vstack(list(train_set.values()))) + data_mean, data_std, dim_to_ignore, dim_to_use = normalization_stats(complete_train, dim=3, predict_14=predict_14) + + # Divide every dimension independently + train_set = normalize_data(train_set, data_mean, data_std, dim_to_use) + test_set = normalize_data(test_set, data_mean, data_std, dim_to_use) + + return train_set, test_set, data_mean, data_std, dim_to_ignore, dim_to_use, train_root_positions, test_root_positions + + +def postprocess_3d(poses_set): + """Center 3d points around root + + Args + poses_set: dictionary with 3d data + Returns + poses_set: dictionary with 3d data centred around root (center hip) joint + root_positions: dictionary with the original 3d position of each pose + """ + root_positions = {} + for k in poses_set.keys(): + # Keep track of the global position + root_positions[k] = copy.deepcopy(poses_set[k][:, :3]) + + # Remove the root from the 3d position, so that other joints get equal weight and model doesnt focus just on the root + poses = poses_set[k] + poses = poses - np.tile(poses[:, :3], [1, len(H36M_NAMES)]) + poses_set[k] = poses + + return poses_set, root_positions + + +def create_2d_data(actions, data_dir, rcams): + """Creates 2d poses by projecting 3d poses with the corresponding camera + parameters. Also normalizes the 2d poses + + Args + actions: list of strings. Actions to load + data_dir: string. Directory where the data can be loaded from + rcams: dictionary with camera parameters + Returns + train_set: dictionary with projected 2d poses for training + test_set: dictionary with projected 2d poses for testing + data_mean: vector with the mean of the 2d training data + data_std: vector with the standard deviation of the 2d training data + dim_to_ignore: list with the dimensions to not predict + dim_to_use: list with the dimensions to predict + """ + + # Load 3d data + train_set = load_data(data_dir, TRAIN_SUBJECTS, actions, dim=3) + test_set = load_data(data_dir, TEST_SUBJECTS, actions, dim=3) + + # Create 2d data by projecting with camera parameters + train_set = project_to_cameras(train_set, rcams) + test_set = project_to_cameras(test_set, rcams) + + # Compute normalization statistics. + complete_train = copy.deepcopy(np.vstack(list(train_set.values()))) + data_mean, data_std, dim_to_ignore, dim_to_use = normalization_stats(complete_train, dim=2) + + # Divide every dimension independently + train_set = normalize_data(train_set, data_mean, data_std, dim_to_use) + test_set = normalize_data(test_set, data_mean, data_std, dim_to_use) + + return train_set, test_set, data_mean, data_std, dim_to_ignore, dim_to_use + + +def project_to_cameras(poses_set, cams, ncams=4): + """Project 3d poses using camera parameters + + Args + poses_set: dictionary with 3d poses + cams: dictionary with camera parameters + ncams: number of cameras per subject + Returns + t2d: dictionary with 2d poses + """ + t2d = {} + + for t3dk in sorted(poses_set.keys()): + subj, a, seqname = t3dk + t3d = poses_set[t3dk] + + for cam in range(ncams): + R, T, f, c, k, p, name = cams[(subj, cam + 1)] + pts2d, _, _, _, _ = cameras.project_point_radial(np.reshape(t3d, [-1, 3]), R, T, f, c, k, p) + + pts2d = np.reshape(pts2d, [-1, len(H36M_NAMES) * 2]) + sname = seqname[:-3] + name + ".h5" # e.g.: Waiting 1.58860488.h5 + t2d[(subj, a, sname)] = pts2d + + return t2d diff --git a/examples/human_pose_3d/global_pose.py b/examples/human_pose_3d/global_pose.py new file mode 100755 index 000000000..818ade9c8 --- /dev/null +++ b/examples/human_pose_3d/global_pose.py @@ -0,0 +1,171 @@ +"""Predicting 3d poses from 2d joints""" +import os +import pickle +import numpy as np +from scipy.optimize import * +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +import copy +os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2' +import tensorflow as tf +from linear_model import mse_loss +import data_utils +import viz +import helper_functions + +from tensorflow.keras.utils import get_file +from paz.applications import HigherHRNetHumanPose2D +from paz.backend.image import load_image, show_image +from paz.backend.camera import Camera + +def joints_2d_from_image(): + # URL = ('https://github.com/oarriaga/altamira-data/releases/download' + # '/v0.10/single_person_test_pose.png') + # filename = os.path.basename(URL) + # fullpath = get_file(filename, URL, cache_subdir='paz/tests') + path = '/home/dfki.uni-bremen.de/kshinde/Downloads/test_image.jpg' + image = load_image(path) + H, W = image.shape[:2] + detect = HigherHRNetHumanPose2D() + inferences = detect(image) + # image = inferences['image'] + # show_image(image) + return inferences['keypoints'], H, W + + +def optimize_trans(initial_root_translation, poses3d, Ki, f, img_center): + """Optimization function to minimize the distance between 2d poses and projection of 3d poses + + Args + initial_root_translation: initial guess of absolute position of root joint in camera space + poses3d: 3d poses to be optimized + Ki: 2d poses + f: focal length + img_center: principal point of the camera + Returns + person_sum: sum of L2 distances between each joint per person + """ + # add root translation to poses3d + initial_root_translation = np.reshape(initial_root_translation, (-1, 3)) + new_poses3d = poses3d + np.tile(initial_root_translation, (1, 16)) + + # Project all poses translation 3D to 2D + ppts = helper_functions.proj_3d_to_2d(new_poses3d.reshape((-1, 3)), f, img_center) + ppts = ppts.reshape((Ki.shape[0],-1,2)) + Ki = Ki.reshape((Ki.shape[0],-1,2)) + person_sum = 0 + + for i in range(Ki.shape[0]): + person_sum += np.sum(np.linalg.norm(Ki[i] - ppts[i], axis=1)) + + return person_sum + + +def predict_3d_poses(): + """Predicts 3d human pose for each person from the multi-human 2d poses obtained from HigherHRNet""" + + poses_2d, img_h, img_w = joints_2d_from_image() + + # Load 2d and 3d normalization stats for H36M dataset + data_mean_2d, data_std_2d, dim_to_use_2d, dim_to_ignore_2d, data_mean_3d, \ + data_std_3d, dim_to_use_3d, dim_to_ignore_3d = data_utils.load_params() + print("\n==> done loading normalization stats.") + + poses_2d = data_utils.load_joints_2d(poses_2d) + print(f"poses_2d : {poses_2d} {poses_2d.shape}") + + # Normalize 2d poses + mu = data_mean_2d[dim_to_use_2d] + stddev = data_std_2d[dim_to_use_2d] + enc_in = np.divide((poses_2d - mu), stddev) + + # load the model + model_path = '/home/dfki.uni-bremen.de/kshinde/Projects/models/baseline_model' + # latter part added because custom loss is defined, is a TF bug + model = tf.keras.models.load_model(model_path, custom_objects={'mse_loss': mse_loss}) + print("\n==> Model loaded!") + + # pass 2d poses and get predictions + poses3d = model.predict(enc_in) + + # denormalize + poses3d = data_utils.unNormalizeData(poses3d, data_mean_3d, data_std_3d, dim_to_ignore_3d) + poses3d_copy = poses3d.copy() + + return poses_2d, poses3d, poses3d_copy, img_h, img_w + + +def translate_root(): + """Finds the optimal translation of root joint for each person to give a good enough estimate + of the global human pose in camera coordinates""" + + poses_2d, poses3d, poses3d_copy, img_h, img_w = predict_3d_poses() + + p3d_16 = data_utils.filter_moving_joints_3d(poses3d) + Ki = poses_2d.astype(np.float32) # 2d poses + + # get human root joint in 2d + root_2d = poses_2d[:, :2] + + # FIXME: Get the intrinsics (image center and focal length) information from the image + cam = Camera() + cam.intrinsics_from_HFOV(HFOV=70, image_shape=[img_h, img_w]) + f = cam.intrinsics[0, 0] + img_center = np.array([[cam.intrinsics[0, 2], cam.intrinsics[1, 2]]]) + + s2d = helper_functions.s2d(poses_2d) + s3d = helper_functions.s3d(p3d_16) + + initial_root_translation = helper_functions.init_translation(f, root_2d, img_center, s2d, s3d) + initial_root_translation = initial_root_translation.flatten() + + root_translation = least_squares(optimize_trans, initial_root_translation, verbose=0, args=(p3d_16, Ki, f, img_center)) + + print(f"\nOPTIMIZATION result : {root_translation}\n{root_translation.x}\n{root_translation.x.shape}") + + root_translation = np.reshape(root_translation.x, (-1, 3)) + root_translation = np.tile(root_translation, (1, 32)) + + # Get global pose Pg = Pi + t* + new_ppts = np.zeros(shape=(poses_2d.shape[0], 64)) + for i in range(poses3d.shape[0]): + poses3d[i] = poses3d[i] + root_translation[i] + ppts = helper_functions.proj_3d_to_2d(poses3d[i].reshape((-1, 3)), f, img_center) + new_ppts[i] = np.reshape(ppts, [1, 64]) + print(f"\nposes3d after optimization {poses3d} {poses3d.shape}") + print(f"\nRoots after optimization {poses3d[:,:3]} {poses3d[:,:3].shape}") + + visualize(poses_2d, p3d_16, poses3d, new_ppts) + + +def visualize(poses_2d, poses3d, ps3d, ppts): + # 1080p = 1,920 x 1,080 + fig = plt.figure(figsize=(19.2, 10.8)) + gs1 = gridspec.GridSpec(1,4) + gs1.update(wspace=-0.00, hspace=0.05) # set the spacing between axes. + plt.axis('off') + + ax = plt.subplot(gs1[0]) + viz.show2Dpose(poses_2d, ax, add_labels=True) + ax.invert_yaxis() + ax.title.set_text('HRNet 2D poses') + + ax1 = plt.subplot(gs1[1], projection='3d') + ax1.view_init(-90, -90) + viz.show3Dpose(poses3d, ax1, add_labels=True) + ax1.title.set_text('Baseline prediction') + + ax2 = plt.subplot(gs1[2], projection='3d') + ax2.view_init(-90, -90) + viz.show3Dpose(ps3d, ax2, add_labels=True) + ax2.title.set_text('Optimized 3D poses') + + ax3 = plt.subplot(gs1[3]) + viz.show2Dpose(ppts, ax3, add_labels=True) + ax3.invert_yaxis() + ax3.title.set_text('2D projection of optimized poses') + plt.show() + + +if __name__ == "__main__": + translate_root() diff --git a/examples/human_pose_3d/helper_functions.py b/examples/human_pose_3d/helper_functions.py new file mode 100644 index 000000000..c7170e4f5 --- /dev/null +++ b/examples/human_pose_3d/helper_functions.py @@ -0,0 +1,114 @@ +import numpy as np + + +def proj_3d_to_2d(P, f, c): + """ + Project points in camera frame from 3d to 2d using intrinsic matrix of the camera + + Args + P: Nx3 points in camera coordinates + f: (scalar) Camera focal length + c: 2x1 image center + Returns + Nx2 points in pixel space + """ + assert len(P.shape) == 2 + assert P.shape[1] == 3 + + z = P[:, 2] + x = (f / z) * P[:, 0] + c[0, 0] + y = (f / z) * P[:, 1] + c[0, 1] + return np.column_stack((x, y)) + + +def s2d(poses2d): + """Computes sum of bone lengths in 2d + + Args + poses2d: np array of poses in 2d + + Returns + sum_bl: sum of length of all bones in the 2d skeleton + """ + assert poses2d[0].shape == (32,), "channels should have 32 entries, it has %d instead" % poses2d[0].shape + sum_bl = np.zeros(poses2d.shape[0]) # sum of bone lengths, each entry is for each person + poses2d = np.reshape(poses2d, (poses2d.shape[0], 16, -1)) + + start_joints = np.array([1, 2, 3, 1, 5, 6, 1, 8, 9, 9, 11, 12, 9, 14, 15]) - 1 + end_joints = np.array([2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]) - 1 + + for idx, person in enumerate(poses2d): + for i in np.arange(len(start_joints)): + bone_length = np.linalg.norm(person[start_joints[i]] - person[end_joints[i]]) + sum_bl[idx] += bone_length + return sum_bl + + +def s3d(poses3d): + """Computes sum of bone lengths in 3d + + Args + poses3d: np array of predicted poses in 3d + + Returns + sum_bl: sum of length of all bones in the 3d skeleton + """ + sum_bl = np.zeros(poses3d.shape[0]) # sum of bone lengths, each entry is for each person + poses3d = np.reshape(poses3d, (poses3d.shape[0], 16, -1)) + + # start_joints = np.array([1, 2, 3, 1, 7, 8, 1, 13, 14, 15, 14, 18, 19, 14, 26, 27]) - 1 + # end_joints = np.array([2, 3, 4, 7, 8, 9, 13, 14, 15, 16, 18, 19, 20, 26, 27, 28]) - 1 + #TODO: CHECK THIS PART THE INDICES!!!! + start_joints = np.array([1, 2, 3, 1, 5, 6, 1, 8, 9, 9, 11, 12, 9, 14, 15]) - 1 + end_joints = np.array([2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]) - 1 + + for idx, person in enumerate(poses3d): + for i in np.arange(len(start_joints)): + bone_length = np.linalg.norm(person[start_joints[i]] - person[end_joints[i]]) + sum_bl[idx] += bone_length + + return sum_bl + + +def init_translation(f, root_2d, img_center, s2d, s3d): + """Computes initial 3d translation of root joint + + Args + f: focal length of the camera in pixels + root_2d: 2d root joint from HigherHRNet + img_center: center of the image (or principal point) + s2d: sum of bone lengths of 2d skeleton + s3d: sum of bone lengths of 3d skeleton (or can be of it's orthographic projection) + + Returns + Array of initial estimate of the global position of the root joint in 3d + """ + ratio = s3d / s2d + ox, oy = img_center[0] + Z = f * ratio # depth coord + X = (root_2d[:, 0] - ox) * ratio # horz coord + Y = (root_2d[:, 1] - oy) * ratio # vert coord + return np.column_stack((X, Y, Z)) + + +def orthographic_proj(P): + """Computes orthographic projection of 3d pose + + Args + P: 3d pose + + Returns + Array containing the X and Y coordinate of 3d pose + """ + assert len(P.shape) == 2 + assert P.shape[1] == 3 + + x = P[:, 0] + y = P[:, 1] + return np.column_stack((x, y)) + + +def pix_to_cam(point, f, c): + point = (point - c) / f + z = 2.8 # f in mm for ZED camera, replace with focal length of the camera used + return np.insert(point, 2, z, axis=1) diff --git a/examples/human_pose_3d/linear_model.py b/examples/human_pose_3d/linear_model.py new file mode 100644 index 000000000..4213dca25 --- /dev/null +++ b/examples/human_pose_3d/linear_model.py @@ -0,0 +1,154 @@ +"""Simple model to regress 3d human poses from 2d joint locations""" + +import numpy as np +import tensorflow as tf +from tensorflow.keras.models import Model +from tensorflow.keras.layers import ( + BatchNormalization, + ReLU, + Dense, + Dropout, + Input +) + + +def two_linear(xin, linear_size, residual, dropout_keep_prob, max_norm, batch_norm, idx): + """Make a bi-linear block with optional residual connection + + Args + xin: the batch that enters the block + linear_size: integer. The size of the linear units + residual: boolean. Whether to add a residual connection + dropout_keep_prob: float [0,1]. Probability of dropping something out + max_norm: boolean. Whether to clip weights to 1-norm + batch_norm: boolean. Whether to do batch normalization + idx: integer. Number of layer (for naming/scoping) + Returns + y: the batch after it leaves the block + """ + initializer = tf.keras.initializers.HeNormal() + # Linear 1 + if max_norm: + y = Dense(linear_size, use_bias=True, kernel_initializer=initializer, bias_initializer=initializer, + kernel_constraint=tf.keras.constraints.MaxNorm(max_value=1), name="linear2_" + str(idx))(xin) + else: + y = Dense(linear_size, use_bias=True, kernel_initializer=initializer, bias_initializer=initializer, + name="linear2_" + str(idx))(xin) + + if batch_norm: + y = BatchNormalization(name="batch_normalization1" + str(idx))(y) # , training=isTraining + + y = ReLU()(y) + y = Dropout(dropout_keep_prob)(y) + + # Linear 2 + if max_norm: + y = Dense(linear_size, use_bias=True, kernel_initializer=initializer, bias_initializer=initializer, + kernel_constraint=tf.keras.constraints.MaxNorm(max_value=1), name="linear3_" + str(idx))(y) + else: + y = Dense(linear_size, use_bias=True, kernel_initializer=initializer, bias_initializer=initializer, + name="linear3_" + str(idx))(y) + + if batch_norm: + y = BatchNormalization(name="batch_normalization2" + str(idx))(y) + + y = ReLU()(y) + y = Dropout(dropout_keep_prob)(y) + + # Residual every 2 blocks + y = (xin + y) if residual else y + + return y + + +def get_all_batches(data_x, data_y, camera_frame, batch_size): + """Obtain a list of all the batches, randomly permuted + + Args + data_x: dictionary with 2d inputs + data_y: dictionary with 3d expected outputs + camera_frame: whether the 3d data is in camera coordinates + + Returns + encoder_inputs: list of 2d batches + decoder_outputs: list of 3d batches + """ + input_size = 16*2 + output_size = 16*3 + # Figure out how many frames we have + n = 0 + for key2d in data_x.keys(): + n2d, _ = data_x[key2d].shape + n = n + n2d + + encoder_inputs = np.zeros((n, input_size), dtype=float) + decoder_outputs = np.zeros((n, output_size), dtype=float) + + # Put all the data into big arrays + idx = 0 + for key2d in data_x.keys(): + (subj, b, fname) = key2d + # keys should be the same if 3d is in camera coordinates + key3d = key2d if (camera_frame) else (subj, b, '{0}.h5'.format(fname.split('.')[0])) + key3d = (subj, b, fname[:-3]) if fname.endswith('-sh') and camera_frame else key3d + + n2d, _ = data_x[key2d].shape + encoder_inputs[idx:idx + n2d, :] = data_x[key2d] + decoder_outputs[idx:idx + n2d, :] = data_y[key3d] + idx = idx + n2d + + # Make the number of examples a multiple of the batch size + n_extra = n % batch_size + if n_extra > 0: # Otherwise examples are already a multiple of batch size + encoder_inputs = encoder_inputs[:-n_extra, :] + decoder_outputs = decoder_outputs[:-n_extra, :] + + return encoder_inputs, decoder_outputs + + +def mse_loss(y_true, y_pred): + loss = tf.reduce_mean(tf.square(y_true - y_pred)) + return loss + + +def LinearModel( + linear_size, + num_layers, + residual, + max_norm, + batch_norm, + dropout_keep_prob, + predict_14, + input_shape=(32,) + ): + HUMAN_2D_SIZE = 16 * 2 + HUMAN_3D_SIZE = 14 * 3 if predict_14 else 16 * 3 + + inputs = Input(shape=input_shape) + initializer = tf.keras.initializers.HeNormal() + if max_norm: + y3 = Dense(linear_size, use_bias=True, kernel_initializer=initializer, bias_initializer=initializer, + kernel_constraint=tf.keras.constraints.MaxNorm(max_value=1), name="linear1_")(inputs) + else: + y3 = Dense(linear_size, use_bias=True, kernel_initializer=initializer, bias_initializer=initializer, + name="linear1_")(inputs) + if batch_norm: + y3 = BatchNormalization(name="batch_normalization")(y3) + y3 = ReLU()(y3) + y3 = Dropout(dropout_keep_prob)(y3) + + # === Create multiple bi-linear layers === + for idx in range(num_layers): + y3 = two_linear(y3, linear_size, residual, dropout_keep_prob, max_norm, batch_norm, idx) + + # === Last linear layer has HUMAN_3D_SIZE in output === + if max_norm: + y = Dense(HUMAN_3D_SIZE, use_bias=True, kernel_initializer=initializer, bias_initializer=initializer, + kernel_constraint=tf.keras.constraints.MaxNorm(max_value=1), name="linear4_")(y3) + else: + y = Dense(HUMAN_3D_SIZE, use_bias=True, kernel_initializer=initializer, bias_initializer=initializer, + name="linear4_")(y3) + # === End linear model === + model = Model(inputs, outputs=y, name='LinearModel') + return model + diff --git a/examples/human_pose_3d/procrustes.py b/examples/human_pose_3d/procrustes.py new file mode 100644 index 000000000..59baae78e --- /dev/null +++ b/examples/human_pose_3d/procrustes.py @@ -0,0 +1,64 @@ +import numpy as np + + +def compute_similarity_transform(X, Y, compute_optimal_scale=False): + """ + A port of MATLAB's `procrustes` function to Numpy. + Adapted from http://stackoverflow.com/a/18927641/1884420 + + Args + X: array NxM of targets, with N number of points and M point dimensionality + Y: array NxM of inputs + compute_optimal_scale: whether we compute optimal scale or force it to be 1 + + Returns: + d: squared error after transformation + Z: transformed Y + T: computed rotation + b: scaling + c: translation + """ + + muX = X.mean(0) + muY = Y.mean(0) + + X0 = X - muX + Y0 = Y - muY + + ssX = (X0**2.).sum() + ssY = (Y0**2.).sum() + + # centred Frobenius norm + normX = np.sqrt(ssX) + normY = np.sqrt(ssY) + + # scale to equal (unit) norm + X0 = X0 / normX + Y0 = Y0 / normY + + # optimum rotation matrix of Y + A = np.dot(X0.T, Y0) + U,s,Vt = np.linalg.svd(A,full_matrices=False) + V = Vt.T + T = np.dot(V, U.T) + + # Make sure we have a rotation + detT = np.linalg.det(T) + V[:,-1] *= np.sign( detT ) + s[-1] *= np.sign( detT ) + T = np.dot(V, U.T) + + traceTA = s.sum() + + if compute_optimal_scale: # Compute optimum scaling of Y. + b = traceTA * normX / normY + d = 1 - traceTA**2 + Z = normX*traceTA*np.dot(Y0, T) + muX + else: # If no scaling allowed + b = 1 + d = 1 + ssY/ssX - 2 * traceTA * normY / normX + Z = normY*np.dot(Y0, T) + muX + + c = muX - b*np.dot(muY, T) + + return d, Z, T, b, c diff --git a/examples/human_pose_3d/train.py b/examples/human_pose_3d/train.py new file mode 100644 index 000000000..efea21a66 --- /dev/null +++ b/examples/human_pose_3d/train.py @@ -0,0 +1,322 @@ +"""Predicting 3d poses from 2d joints""" +import os +import sys +import time +import numpy as np + +os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2' +import tensorflow as tf +from tensorflow.keras.optimizers import Adam +from tensorflow.keras.callbacks import Callback, TensorBoard, LearningRateScheduler +from tensorflow.python.keras.callbacks import CallbackList +from absl import app, flags, logging +from absl.flags import FLAGS +from datetime import datetime +import cameras +import data_utils +from linear_model import LinearModel +from linear_model import mse_loss, get_all_batches +import procrustes + +flags.DEFINE_float("learning_rate", 1e-3, "Learning rate") +flags.DEFINE_float("dropout", 1, "Dropout keep probability. 1 means no dropout") +flags.DEFINE_integer("batch_size", 64, "Batch size to use during training") +flags.DEFINE_integer("epochs", 200, "How many epochs we should train for") +flags.DEFINE_boolean("camera_frame", False, "Convert 3d poses to camera coordinates") +flags.DEFINE_boolean("max_norm", False, "Apply maxnorm constraint to the weights") +flags.DEFINE_boolean("batch_norm", False, "Use batch_normalization") + +# Data loading +flags.DEFINE_boolean("predict_14", False, "predict 14 joints") +flags.DEFINE_string("action", "All", "The action to train on. 'All' means all the actions") + +# Architecture +flags.DEFINE_integer("linear_size", 1024, "Size of each model layer.") +flags.DEFINE_integer("num_layers", 2, "Number of layers in the model.") +flags.DEFINE_boolean("residual", False, "Whether to add a residual connection every 2 layers") + +# Evaluation +flags.DEFINE_boolean("procrustes", False, "Apply procrustes analysis at test time") +flags.DEFINE_boolean("evaluateActionWise", False, "The dataset to use either h36m or heva") + +# Directories +flags.DEFINE_string("cameras_path", "data/h36m/metadata.xml", "File with h36m metadata, including cameras") +flags.DEFINE_string("data_dir", "SCRATCH/3d-pose-baseline/data/h36m/", "Data directory") +flags.DEFINE_string("train_dir", "SCRATCH/3d-pose-baseline/experiments", "Training directory.") + +# Train or load +flags.DEFINE_boolean("sample", False, "Set to True for sampling.") +flags.DEFINE_boolean("use_cpu", False, "Whether to use the CPU") +flags.DEFINE_integer("load", 0, "Try to load a previous checkpoint.") + +# Misc +flags.DEFINE_boolean("use_fp16", False, "Train using fp16 instead of fp32.") + + +class CustomCallbackList(CallbackList): + """This Class avoids the warning printed when callback takes more time as compared to training""" + + def _call_batch_hook(self, mode, hook, batch, logs=None): + """Helper function for all batch_{begin | end} methods.""" + if not self.callbacks: + return + hook_name = 'on_{mode}_batch_{hook}'.format(mode=mode, hook=hook) + + logs = logs or {} + for callback in self.callbacks: + batch_hook = getattr(callback, hook_name) + batch_hook(batch, logs) + + +tf.keras.callbacks.CallbackList = CustomCallbackList # tf.python.keras.callbacks.CallbackList + + +def get_train_dir(): + return os.path.join(FLAGS.train_dir, + FLAGS.action, + 'dropout_{0}'.format(FLAGS.dropout), + 'epochs_{0}'.format(FLAGS.epochs) if FLAGS.epochs > 0 else '', + 'lr_{0}'.format(FLAGS.learning_rate), + 'residual' if FLAGS.residual else 'not_residual', + 'depth_{0}'.format(FLAGS.num_layers), + 'linear_size{0}'.format(FLAGS.linear_size), + 'batch_size_{0}'.format(FLAGS.batch_size), + 'procrustes' if FLAGS.procrustes else 'no_procrustes', + 'maxnorm' if FLAGS.max_norm else 'no_maxnorm', + 'batch_normalization' if FLAGS.batch_norm else 'no_batch_normalization', + 'predict_14' if FLAGS.predict_14 else 'predict_17') + + +def denormalize(enc_in, dec_out, poses3d, data_mean_2d, data_std_2d, + dim_to_ignore_2d, data_mean_3d, data_std_3d, dim_to_ignore_3d): + """ + Function that denormalizes the inputs + + Args + data_mean_3d: the mean of the training data in 3d + data_std_3d: the standard deviation of the training data in 3d + dim_to_use_3d: out of all the 96 dimensions that represent a 3d body in h36m, compute results for this subset + dim_to_ignore_3d: complelment of the above + data_mean_2d: mean of the training data in 2d + data_std_2d: standard deviation of the training data in 2d + dim_to_use_2d: out of the 64 dimensions that represent a body in 2d in h35m, use this subset + dim_to_ignore_2d: complement of the above + encoder_inputs: input for the network + decoder_outputs: expected output for the network + + Returns + enc_in: denormalized encoder inputs + dec_out: adenormalized decoder outputs + poses3d: denormalized 3D poses + """ + # denormalize + enc_in = data_utils.unNormalizeData(enc_in, data_mean_2d, data_std_2d, dim_to_ignore_2d) + dec_out = data_utils.unNormalizeData(dec_out, data_mean_3d, data_std_3d, dim_to_ignore_3d) + poses3d = data_utils.unNormalizeData(poses3d, data_mean_3d, data_std_3d, dim_to_ignore_3d) + + return enc_in, dec_out, poses3d + + +def calculate_err(poses3d, dec_out, n_joints, all_dists): + """ + Compute Euclidean distance error per joint + + Args + poses3d: predicted 3d poses + dec_out: expected output for the network + n_joints: number of joints + all_dists: empty list + Returns + all_dists: list of L2 distance + """ + # Compute Euclidean distance error per joint + sqerr = (poses3d - dec_out) ** 2 # Squared error between prediction and expected output + dists = np.zeros((sqerr.shape[0], n_joints)) # Array with L2 error per joint in mm + dist_idx = 0 + for k in np.arange(0, n_joints * 3, 3): + # Sum across X,Y, and Z dimenstions to obtain L2 distance + dists[:, dist_idx] = np.sqrt(np.sum(sqerr[:, k:k + 3], axis=1)) + dist_idx = dist_idx + 1 + + all_dists.append(dists) + + return all_dists + + +def lr_exp_decay(epoch): + n_batch = 24371 + decay_rate = 0.96 + decay_steps = 100000 + p = ((epoch + 1) * n_batch) / decay_steps + lr = tf.multiply(FLAGS.learning_rate, tf.pow(decay_rate, p)) + tf.summary.scalar('learning_rate', data=lr, step=epoch) + return lr + + +class ValCallback(Callback): + def on_test_batch_end(self, batch, logs=None): + if (batch + 1) % 1000 == 603: + print("...Evaluating: batch {}".format(batch + 1)) + + +class EvaluateEpoch(Callback): + def __init__(self, enc_inputs_val, dec_outputs_val, model, data_mean_2d, data_std_2d, dim_to_ignore_2d, + data_mean_3d, data_std_3d, dim_to_use_3d, dim_to_ignore_3d, logdir): + super(EvaluateEpoch, self).__init__() + self.enc_in = enc_inputs_val + self.dec_out = dec_outputs_val + self.model = model + self.dm_2d = data_mean_2d + self.ds_2d = data_std_2d + self.dti_2d = dim_to_ignore_2d + self.dm_3d = data_mean_3d + self.ds_3d = data_std_3d + self.dtu_3d = dim_to_use_3d + self.dti_3d = dim_to_ignore_3d + self.logdir = logdir + + def on_epoch_end(self, epoch, logs=None): + n_joints = 17 + all_dists, start_time, loss = [], time.time(), 0. + + print("==> Working on test epoch {0}".format(epoch + 1)) + + dp = 1.0 # dropout keep probability is always 1 at test time + poses3d = self.model.predict(x=self.enc_in) + poses3d = np.asarray(poses3d) + step_loss = mse_loss(self.dec_out, poses3d) + loss += step_loss + + # denormalize + enc_in, dec_out, poses3d = denormalize(self.enc_in, self.dec_out, poses3d, self.dm_2d, self.ds_2d, + self.dti_2d, self.dm_3d, self.ds_3d, self.dti_3d) + + # Keep only the relevant dimensions + dtu3d = np.hstack((np.arange(3), self.dtu_3d)) + dec_out = dec_out[:, dtu3d] + poses3d = poses3d[:, dtu3d] + + if FLAGS.procrustes: + # Apply per-frame procrustes alignment if asked to do so + for j in range(FLAGS.batch_size): + gt = np.reshape(dec_out[j, :], [-1, 3]) + out = np.reshape(poses3d[j, :], [-1, 3]) + _, Z, T, b, c = procrustes.compute_similarity_transform(gt, out, compute_optimal_scale=True) + out = (b * out.dot(T)) + c + + poses3d[j, :] = np.reshape(out, [-1, 17 * 3]) + + all_dists = calculate_err(poses3d, dec_out, n_joints, all_dists) + step_time = (time.time() - start_time) + all_dists = np.vstack(all_dists) + + # Error per joint and total for all passed batches + joint_err = np.mean(all_dists, axis=0) + total_err = np.mean(all_dists) + + print("=============================\n" + "Step-time (s): %.4f\n" + "Val loss avg: %.4f\n" + "Val error avg (mm): %.2f\n" + "=============================" % (step_time, loss, total_err)) + + for i in range(n_joints): + # 6 spaces, right-aligned, 5 decimal places + print("Error in joint {0:02d} (mm): {1:>5.2f}".format(i + 1, joint_err[i])) + print("=============================") + + # Saving the evaluation results + filename = os.path.join(self.logdir, 'Val_log.txt') + if not os.path.isdir(os.path.dirname(filename)): + os.makedirs(os.path.dirname(filename)) + with open(filename, 'a') as eval_log_file: + eval_log_file.write('Epoch: {}, total_err: {}, joint_err: {}, loss: {}, step_time: {} s\n'. + format(str(epoch), total_err, joint_err, loss, step_time)) + + +def train(): + """Train a linear model for 3d pose estimation""" + train_dir = get_train_dir() + print(f"==> train_dir {train_dir}") + + # Logs dir for train and test runs + path_prefix = os.path.dirname(os.path.abspath(__file__)) + logdir = os.path.join(path_prefix, 'logs/scalars/', datetime.now().strftime("%Y%m%d-%H%M%S")) + save_dir = os.path.join(path_prefix, 'saved_model/') + + if not os.path.exists(save_dir): + os.system('mkdir -p {}'.format(save_dir)) + file_writer = tf.summary.create_file_writer(logdir + "/metrics") + file_writer.set_as_default() + + actions = data_utils.define_actions(FLAGS.action) + + # Load camera parameters + SUBJECT_IDS = [1, 5, 6, 7, 8, 9, 11] + this_file = os.path.dirname(os.path.realpath(__file__)) + rcams = cameras.load_cameras(os.path.join(this_file, "..", FLAGS.cameras_path), SUBJECT_IDS) + + # Load 3d data and load (or create) 2d projections + train_set_3d, test_set_3d, data_mean_3d, data_std_3d, dim_to_ignore_3d, dim_to_use_3d, train_root_positions, \ + test_root_positions = data_utils.read_3d_data(actions, FLAGS.data_dir, FLAGS.camera_frame, rcams, FLAGS.predict_14) + + # Read groundtruth 2D projections + train_set_2d, test_set_2d, data_mean_2d, data_std_2d, dim_to_ignore_2d, dim_to_use_2d = data_utils.create_2d_data( + actions, FLAGS.data_dir, rcams) + print("\n==> done reading and normalizing data.") + + # === Create the model === + print("\n==> Creating %d bi-layers of %d units." % (FLAGS.num_layers, FLAGS.linear_size)) + model = LinearModel( + FLAGS.linear_size, + FLAGS.num_layers, + FLAGS.residual, + FLAGS.max_norm, + FLAGS.batch_norm, + FLAGS.dropout, + FLAGS.predict_14, + ) + print("\n==> Model created!") + + # Define loss function (criterion) and optimizer and compile model + + lr_schedule = tf.keras.optimizers.schedules.ExponentialDecay( + FLAGS.learning_rate, + decay_steps=100000, + decay_rate=0.96, + ) + optimizer = Adam(learning_rate=FLAGS.learning_rate) # lr_schedule + criterion = mse_loss + model.compile(optimizer, criterion) # model.compile(opt, loss) is syntax + + encoder_inputs, decoder_outputs = get_all_batches(train_set_2d, train_set_3d, FLAGS.camera_frame, + FLAGS.batch_size) + + enc_inputs_val, dec_outputs_val = get_all_batches(test_set_2d, test_set_3d, FLAGS.camera_frame, + FLAGS.batch_size) + + # Callbacks + callbacks = [LearningRateScheduler(lr_exp_decay, verbose=1), + EvaluateEpoch(enc_inputs_val, dec_outputs_val, model, data_mean_2d, data_std_2d, dim_to_ignore_2d, + data_mean_3d, data_std_3d, dim_to_use_3d, dim_to_ignore_3d, logdir), + ValCallback(), + TensorBoard(log_dir=logdir)] # TqdmCallback(verbose=2) + start_time_epoch = time.time() + # == Train == + history = model.fit(x=encoder_inputs, y=decoder_outputs, batch_size=FLAGS.batch_size, epochs=FLAGS.epochs, + verbose=2, callbacks=callbacks, validation_data=(enc_inputs_val, dec_outputs_val), shuffle=True) + print("\ncompleted in {0:.2f} s".format((time.time() - start_time_epoch))) + # Save the model + print("\n==> Saving the model... ", end="") + start_time = time.time() + model.save(save_dir + 'baseline_model') + print("\ndone in {0:.2f} ms".format(1000 * (time.time() - start_time))) + + # Reset global time and loss + step_time, loss = 0, 0 + + sys.stdout.flush() + + +if __name__ == "__main__": + train() diff --git a/examples/human_pose_3d/viz.py b/examples/human_pose_3d/viz.py new file mode 100644 index 000000000..cb643f405 --- /dev/null +++ b/examples/human_pose_3d/viz.py @@ -0,0 +1,120 @@ +"""Functions to visualize human poses""" + +import numpy as np +import data_utils + + +def show3Dpose(channels, ax, lcolor="#3498db", rcolor="#e74c3c", add_labels=False): # blue, orange + """Visualize a 3d skeleton + + Args + channels: 48x1 vector. The pose to plot. + ax: matplotlib 3d axis to draw on + lcolor: color for left part of the body + rcolor: color for right part of the body + add_labels: whether to add coordinate labels + Returns + Nothing. Draws on ax. + """ + + # assert channels.shape[1] == 16 * 3, "channels should have 48 entries, it has %d instead" % channels.size + + if channels.shape[1] == 48: + vals = np.reshape(channels, (channels.shape[0], 16, -1)) + # print(f"\nvals 3d {vals} {vals.shape}", flush=True) + I = np.array([1, 2, 3, 1, 5, 6, 1, 8, 9, 9, 11, 12, 9, 14, 15]) - 1 # start points + J = np.array([2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]) - 1 # end points + + else: + # print(f"ELSE loop p3d shape is : {channels.shape}", flush=True) + vals = np.reshape(channels, (channels.shape[0], 32, -1)) + I = np.array([1, 2, 3, 1, 7, 8, 1, 13, 14, 14, 18, 19, 14, 26, 27]) - 1 # start points + J = np.array([2, 3, 4, 7, 8, 9, 13, 14, 16, 18, 19, 20, 26, 27, 28]) - 1 # end points + LR = np.array([1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1], dtype=bool) + + def plot_person(person, lcolor="#3498db", rcolor="#e74c3c"): + for i in np.arange(len(I)): + x, y, z = [np.array([person[I[i], j], person[J[i], j]]) for j in range(3)] + ax.plot(x, y, z, lw=2, c=lcolor if LR[i] else rcolor) + + i = 0 + # Now we plot all persons one by one in same plot + for person in vals: + if i == 1: + plot_person(person, lcolor="#9b59b6", rcolor="#2ecc71") + else: + plot_person(person) + i += 1 + + RADIUS = 750 # space around the subject + xroot, yroot, zroot = vals[0, 0, 0], vals[0, 0, 1], vals[0, 0, 2] + ax.set_xlim3d([-RADIUS + xroot, RADIUS + xroot]) + ax.set_zlim3d([-RADIUS + zroot, RADIUS + zroot]) + ax.set_ylim3d([-RADIUS + yroot, RADIUS + yroot]) + + if add_labels: + ax.set_xlabel("x") + ax.set_ylabel("y") + ax.set_zlabel("z") + + # Get rid of the panes (actually, make them white) + white = (1.0, 1.0, 1.0, 0.0) + ax.w_xaxis.set_pane_color(white) + ax.w_zaxis.set_pane_color(white) + # Keep y (xz) pane + + +def show2Dpose(channels, ax, lcolor="#3498db", rcolor="#e74c3c", add_labels=False): + """Visualize a 2d skeleton + + Args + channels: nx64 vector. n is num persons detected to plot. + ax: matplotlib axis to draw on + lcolor: color for left part of the body + rcolor: color for right part of the body + add_labels: whether to add coordinate labels + Returns + Nothing. Draws on ax. + """ + + if channels.shape[1] == 32: + vals = np.reshape(channels, (channels.shape[0], 16, -1)) + # print(f"\nvals 2d f loop {vals} {vals.shape}") + I = np.array([1, 2, 3, 1, 5, 6, 1, 8, 9, 9, 11, 12, 9, 14, 15]) - 1 # start points + J = np.array([2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]) - 1 # end points + + else: + # print(f"ELSE loop p2d shape is : {channels.shape}") + vals = np.reshape(channels, (channels.shape[0], 32, -1)) + # print(f"\nvals 2d else loop {vals} {vals.shape}") + I = np.array([1, 2, 3, 1, 7, 8, 1, 13, 14, 14, 18, 19, 14, 26, 27]) - 1 # start points + J = np.array([2, 3, 4, 7, 8, 9, 13, 14, 16, 18, 19, 20, 26, 27, 28]) - 1 # end points + LR = np.array([1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1], dtype=bool) + + def plot_person(person): + for i in np.arange(len(I)): + x, y = [np.array([person[I[i], j], person[J[i], j]]) for j in range(2)] + ax.plot(x, y, lw=2, c=lcolor if LR[i] else rcolor) + + # Now we plot all persons one by one in same plot + for person in vals: + plot_person(person) + + # Get rid of the ticks + # ax.set_xticks([]) + # ax.set_yticks([]) + # + # # Get rid of tick labels + # ax.get_xaxis().set_ticklabels([]) + # ax.get_yaxis().set_ticklabels([]) + + RADIUS = 350 # space around the subject + xroot, yroot = vals[0, 0, 0], vals[0, 0, 1] + ax.set_xlim([-RADIUS + xroot, RADIUS + xroot]) + ax.set_ylim([-RADIUS + yroot, RADIUS + yroot]) + + if add_labels: + ax.set_xlabel("x") + ax.set_ylabel("z") + + ax.set_aspect('equal')