diff --git a/squat_analysis/function/handler.py b/squat_analysis/function/handler.py index 724d83b..cd8f4bd 100644 --- a/squat_analysis/function/handler.py +++ b/squat_analysis/function/handler.py @@ -70,82 +70,45 @@ def handler(event, context): lowpass_cutoff_frequency_for_coordinate_values=filter_frequency, n_repetitions=n_repetitions) squat_events = squat.get_squat_events() - - max_knee_flexion_angle_r_mean, max_knee_flexion_angle_r_std, _ = squat.compute_peak_angle('knee_angle_r') - max_knee_flexion_angle_l_mean, max_knee_flexion_angle_l_std, _ = squat.compute_peak_angle('knee_angle_l') - max_knee_flexion_angle_mean_mean = np.round(np.mean(np.array([max_knee_flexion_angle_r_mean, max_knee_flexion_angle_l_mean]))) - max_knee_flexion_angle_mean_std = np.round(np.mean(np.array([max_knee_flexion_angle_r_std, max_knee_flexion_angle_l_std]))) - - max_hip_flexion_angle_r_mean, max_hip_flexion_angle_r_std, _ = squat.compute_peak_angle('hip_flexion_r') - max_hip_flexion_angle_l_mean, max_hip_flexion_angle_l_std, _ = squat.compute_peak_angle('hip_flexion_l') - max_hip_flexion_angle_mean_mean = np.round(np.mean(np.array([max_hip_flexion_angle_r_mean, max_hip_flexion_angle_l_mean]))) - max_hip_flexion_angle_mean_std = np.round(np.mean(np.array([max_hip_flexion_angle_r_std, max_hip_flexion_angle_l_std]))) - - max_hip_adduction_angle_r_mean, max_hip_adduction_angle_r_std, _ = squat.compute_peak_angle('hip_adduction_r') - max_hip_adduction_angle_l_mean, max_hip_adduction_angle_l_std, _ = squat.compute_peak_angle('hip_adduction_l') - max_hip_adduction_angle_mean_mean = np.round(np.mean(np.array([max_hip_adduction_angle_r_mean, max_hip_adduction_angle_l_mean]))) - max_hip_adduction_angle_mean_std = np.round(np.mean(np.array([max_hip_adduction_angle_r_std, max_hip_adduction_angle_l_std]))) - rom_knee_flexion_angle_r_mean, rom_knee_flexion_angle_r_std, _ = squat.compute_range_of_motion('knee_angle_r') - rom_knee_flexion_angle_l_mean, rom_knee_flexion_angle_l_std, _ = squat.compute_range_of_motion('knee_angle_l') - rom_knee_flexion_angle_mean_mean = np.round(np.mean(np.array([rom_knee_flexion_angle_r_mean, rom_knee_flexion_angle_l_mean]))) - rom_knee_flexion_angle_mean_std = np.round(np.mean(np.array([rom_knee_flexion_angle_r_std, rom_knee_flexion_angle_l_std]))) - + # Detect squat type + eventTypes = squat.squatEvents['eventTypes'] + + # Return squat type (none detected, more than one type, or all same type). + squat_type = '' + unique_types = set(eventTypes) + if len(unique_types) < 1: + squat_type = 'No squats detected' + + elif len(unique_types) > 1: + squat_type = 'Mixed squat types detected' + + else: + if eventTypes[0] == 'double_leg': + squat_type = 'Double leg squats' + elif eventTypes[0] == 'single_leg_l': + squat_type = 'Single leg squats (left)' + elif eventTypes[0] == 'single_leg_r': + squat_type = 'Single leg squats (right)' + + # Compute metrics. + max_trunk_lean_ground_mean, max_trunk_lean_ground_std, max_trunk_lean_ground_units = squat.compute_trunk_lean_relative_to_ground() + max_trunk_flexion_mean, max_trunk_flexion_std, max_trunk_flexion_units = squat.compute_trunk_flexion_relative_to_ground() + squat_depth_mean, squat_depth_std, squat_depth_units = squat.compute_squat_depth() + + # Store metrics dictionary. squat_scalars = {} - squat_scalars['peak_knee_flexion_angle_mean'] = {'value': max_knee_flexion_angle_mean_mean} - squat_scalars['peak_knee_flexion_angle_mean']['label'] = 'Mean peak knee flexion angle (deg)' - squat_scalars['peak_knee_flexion_angle_mean']['colors'] = ["red", "yellow", "green"] - peak_knee_flexion_angle_threshold = 100 - squat_scalars['peak_knee_flexion_angle_mean']['min_limit'] = float(np.round(0.90*peak_knee_flexion_angle_threshold)) - squat_scalars['peak_knee_flexion_angle_mean']['max_limit'] = float(peak_knee_flexion_angle_threshold) - - squat_scalars['peak_knee_flexion_angle_std'] = {'value': max_knee_flexion_angle_mean_std} - squat_scalars['peak_knee_flexion_angle_std']['label'] = 'Std peak knee flexion angle (deg)' - squat_scalars['peak_knee_flexion_angle_std']['colors'] = ["green", "yellow", "red"] - std_threshold_min = 2 - std_threshold_max = 4 - squat_scalars['peak_knee_flexion_angle_std']['min_limit'] = float(std_threshold_min) - squat_scalars['peak_knee_flexion_angle_std']['max_limit'] = float(std_threshold_max) - - squat_scalars['peak_hip_flexion_angle_mean'] = {'value': max_hip_flexion_angle_mean_mean} - squat_scalars['peak_hip_flexion_angle_mean']['label'] = 'Mean peak hip flexion angle (deg)' - squat_scalars['peak_hip_flexion_angle_mean']['colors'] = ["red", "yellow", "green"] - peak_hip_flexion_angle_threshold = 100 - squat_scalars['peak_hip_flexion_angle_mean']['min_limit'] = float(np.round(0.90*peak_hip_flexion_angle_threshold)) - squat_scalars['peak_hip_flexion_angle_mean']['max_limit'] = float(peak_hip_flexion_angle_threshold) - - squat_scalars['peak_hip_flexion_angle_std'] = {'value': max_hip_flexion_angle_mean_std} - squat_scalars['peak_hip_flexion_angle_std']['label'] = 'Std peak hip flexion angle (deg)' - squat_scalars['peak_hip_flexion_angle_std']['colors'] = ["green", "yellow", "red"] - squat_scalars['peak_hip_flexion_angle_std']['min_limit'] = float(std_threshold_min) - squat_scalars['peak_hip_flexion_angle_std']['max_limit'] = float(std_threshold_max) - - squat_scalars['peak_knee_adduction_angle_mean'] = {'value': max_hip_adduction_angle_mean_mean} - squat_scalars['peak_knee_adduction_angle_mean']['label'] = 'Mean peak knee adduction angle (deg)' - squat_scalars['peak_knee_adduction_angle_mean']['colors'] = ["red", "green", "red"] - knee_adduction_angle_threshold = 5 - squat_scalars['peak_knee_adduction_angle_mean']['min_limit'] = float(-knee_adduction_angle_threshold) - squat_scalars['peak_knee_adduction_angle_mean']['max_limit'] = float(knee_adduction_angle_threshold) - - squat_scalars['peak_knee_adduction_angle_std'] = {'value': max_hip_adduction_angle_mean_std} - squat_scalars['peak_knee_adduction_angle_std']['label'] = 'Std peak knee adduction angle (deg)' - squat_scalars['peak_knee_adduction_angle_std']['colors'] = ["green", "yellow", "red"] - squat_scalars['peak_knee_adduction_angle_std']['min_limit'] = float(std_threshold_min) - squat_scalars['peak_knee_adduction_angle_std']['max_limit'] = float(std_threshold_max) - - squat_scalars['rom_knee_flexion_angle_mean'] = {'value': rom_knee_flexion_angle_mean_mean} - squat_scalars['rom_knee_flexion_angle_mean']['label'] = 'Mean range of motion knee flexion angle (deg)' - squat_scalars['rom_knee_flexion_angle_mean']['colors'] = ["red", "yellow", "green"] - rom_knee_flexion_angle_threshold_min = 85 - rom_knee_flexion_angle_threshold_max = 115 - squat_scalars['rom_knee_flexion_angle_mean']['min_limit'] = float(rom_knee_flexion_angle_threshold_min) - squat_scalars['rom_knee_flexion_angle_mean']['max_limit'] = float(rom_knee_flexion_angle_threshold_max) - - squat_scalars['rom_knee_flexion_angle_std'] = {'value': rom_knee_flexion_angle_mean_std} - squat_scalars['rom_knee_flexion_angle_std']['label'] = 'Std range of motion knee flexion angle (deg)' - squat_scalars['rom_knee_flexion_angle_std']['colors'] = ["green", "yellow", "red"] - squat_scalars['rom_knee_flexion_angle_std']['min_limit'] = float(std_threshold_min) - squat_scalars['rom_knee_flexion_angle_std']['max_limit'] = float(std_threshold_max) + squat_scalars['max_trunk_lean_ground'] = {'value': np.round(max_trunk_lean_ground_mean, 2), + 'std': np.round(max_trunk_lean_ground_std, 2), + 'label': 'Mean max trunk lean (deg)'} + + squat_scalars['max_trunk_flexion'] = {'value': np.round(max_trunk_flexion_mean, 2), + 'std': np.round(max_trunk_flexion_std, 2), + 'label': 'Mean max trunk flexion (deg)'} + + squat_scalars['squat_depth'] = {'value': np.round(squat_depth_mean, 2), + 'std': np.round(squat_depth_std, 2), + 'label': 'Mean squat depth (m)'} # %% Return indices for visualizer and line curve plot. # %% Create json for deployement. @@ -178,13 +141,14 @@ def handler(event, context): y_axes.remove('mtp_angle_r') y_axes.remove('mtp_angle_l') - # Create results dictionnary. + # Create results dictionary. results = { 'indices': times, 'metrics': squat_scalars, 'datasets': datasets, 'x_axis': 'time', - 'y_axis': y_axes} + 'y_axis': y_axes, + 'squat_type': squat_type} return { 'statusCode': 200, diff --git a/squat_analysis/function/squat_analysis.py b/squat_analysis/function/squat_analysis.py index 75d66c5..eeb4608 100644 --- a/squat_analysis/function/squat_analysis.py +++ b/squat_analysis/function/squat_analysis.py @@ -23,11 +23,11 @@ import numpy as np import pandas as pd -from scipy.signal import find_peaks +from scipy.signal import find_peaks, peak_widths from matplotlib import pyplot as plt from utilsKinematics import kinematics - +import opensim class squat_analysis(kinematics): @@ -53,28 +53,42 @@ def __init__(self, session_dir, trial_name, n_repetitions=-1, # Coordinate values. self.coordinateValues = self.get_coordinate_values() + + # Body transforms + self.bodyTransformDict = self.get_body_transform_dict() # Making sure time vectors of marker and coordinate data are the same. if not np.allclose(self.markerDict['time'], self.coordinateValues['time'], atol=.001, rtol=0): raise Exception('Time vectors of marker and coordinate data are not the same.') - # Trim marker data and coordinate values. + if not np.allclose(self.bodyTransformDict['time'], self.coordinateValues['time'], atol=.001, rtol=0): + raise Exception('Time vectors of body transofrms and coordinate data are not the same.') + + # Trim marker, body transforms, and coordinate data. if self.trimming_start > 0: self.idx_trim_start = np.where(np.round(self.markerDict['time'] - self.trimming_start,6) <= 0)[0][-1] - self.markerDict['time'] = self.markerDict['time'][self.idx_trim_start:,] + self.markerDict['time'] = self.markerDict['time'][self.idx_trim_start:] for marker in self.markerDict['markers']: self.markerDict['markers'][marker] = self.markerDict['markers'][marker][self.idx_trim_start:,:] + self.bodyTransformDict['time'] = self.bodyTransformDict['time'][self.idx_trim_start:] + for body in self.bodyTransformDict['body_transforms']: + self.bodyTransformDict['body_transforms'][body] = \ + self.bodyTransformDict['body_transforms'][body][self.idx_trim_start:] self.coordinateValues = self.coordinateValues.iloc[self.idx_trim_start:] if self.trimming_end > 0: self.idx_trim_end = np.where(np.round(self.markerDict['time'],6) <= np.round(self.markerDict['time'][-1] - self.trimming_end,6))[0][-1] + 1 - self.markerDict['time'] = self.markerDict['time'][:self.idx_trim_end,] + self.markerDict['time'] = self.markerDict['time'][:self.idx_trim_end] for marker in self.markerDict['markers']: self.markerDict['markers'][marker] = self.markerDict['markers'][marker][:self.idx_trim_end,:] + self.bodyTransformDict['time'] = self.bodyTransformDict['time'][self.idx_trim_end:] + for body in self.bodyTransformDict['body_transforms']: + self.bodyTransformDict['body_transforms'][body] = \ + self.bodyTransformDict['body_transforms'][body][self.idx_trim_end:] self.coordinateValues = self.coordinateValues.iloc[:self.idx_trim_end] # Segment squat repetitions. - self.squatEvents = self.segment_squat(n_repetitions=n_repetitions) + self.squatEvents = self.segment_squat(n_repetitions=n_repetitions, visualizeSegmentation=True) self.nRepetitions = np.shape(self.squatEvents['eventIdxs'])[0] # Initialize variables to be lazy loaded. @@ -111,7 +125,7 @@ def compute_scalars(self, scalarNames, return_all=False): nonexistant_methods = [entry for entry in scalarNames if 'compute_' + entry not in method_names] if len(nonexistant_methods) > 0: - raise Exception(str(['compute_' + a for a in nonexistant_methods]) + ' does not exist in gait_analysis class.') + raise Exception(str(['compute_' + a for a in nonexistant_methods]) + ' does not exist in squat_analysis class.') scalarDict = {} for scalarName in scalarNames: @@ -121,37 +135,37 @@ def compute_scalars(self, scalarNames, return_all=False): scalarDict[scalarName]['units']) = thisFunction(return_all=return_all) return scalarDict + - def segment_squat(self, n_repetitions=-1, height_value=0.2, visualizeSegmentation=False): + def segment_squat(self, n_repetitions=-1, peak_proportion_threshold=0.7, + peak_width_rel_height=0.95, peak_distance_sec=0.5, + toe_height_threshold=0.05, + visualizeSegmentation=False): pelvis_ty = self.coordinateValues['pelvis_ty'].to_numpy() dt = np.mean(np.diff(self.time)) # Identify minimums. pelvSignal = np.array(-pelvis_ty - np.min(-pelvis_ty)) - pelvSignalPos = np.array(pelvis_ty - np.min(pelvis_ty)) - idxMinPelvTy,_ = find_peaks(pelvSignal, distance=.7/dt, height=height_value) - - # Find the max adjacent to all of the minimums. - minIdxOld = 0 + pelvRange = np.abs(np.max(pelvis_ty) - np.min(pelvis_ty)) + peakThreshold = peak_proportion_threshold * pelvRange + idxMinPelvTy,_ = find_peaks(pelvSignal, prominence=peakThreshold, + distance=peak_distance_sec/dt) + peakWidths = peak_widths(pelvSignal, idxMinPelvTy, + rel_height=peak_width_rel_height) + + # Store start and end indices. startEndIdxs = [] - for i, minIdx in enumerate(idxMinPelvTy): - if i < len(idxMinPelvTy) - 1: - nextIdx = idxMinPelvTy[i+1] - else: - nextIdx = len(pelvSignalPos) - startIdx = np.argmax(pelvSignalPos[minIdxOld:minIdx]) + minIdxOld - endIdx = np.argmax(pelvSignalPos[minIdx:nextIdx]) + minIdx - startEndIdxs.append([startIdx,endIdx]) - minIdxOld = np.copy(minIdx) + for start_ips, end_ips in zip(peakWidths[2], peakWidths[3]): + startEndIdxs.append([start_ips, end_ips]) + startEndIdxs = np.rint(startEndIdxs).astype(int) # Limit the number of repetitions. if n_repetitions != -1: startEndIdxs = startEndIdxs[-n_repetitions:] - # Extract events: start and end of each repetition. - eventIdxs = np.array(startEndIdxs) - eventTimes = self.time[eventIdxs] + # Store start and end event times + eventTimes = self.time[startEndIdxs] if visualizeSegmentation: plt.figure() @@ -162,16 +176,35 @@ def segment_squat(self, n_repetitions=-1, height_value=0.2, visualizeSegmentatio label='Start/End rep') if c_v == 0: plt.legend() + plt.hlines(-peakWidths[1], peakWidths[2], peakWidths[3], color='k', + label='peak start/end') plt.xlabel('Frames') plt.ylabel('Position [m]') plt.title('Vertical pelvis position') plt.draw() + + # Detect squat type (double leg, single leg right, single leg left) + # Use toe markers + eventTypes = [] + markersDict = self.markerDict['markers'] + + for eventIdx in startEndIdxs: + lToeYMean = np.mean(markersDict['L_toe_study'][eventIdx[0]:eventIdx[1]+1, 1]) + rToeYMean = np.mean(markersDict['r_toe_study'][eventIdx[0]:eventIdx[1]+1, 1]) + if lToeYMean - rToeYMean > toe_height_threshold: + eventTypes.append('single_leg_r') + elif rToeYMean - lToeYMean > toe_height_threshold: + eventTypes.append('single_leg_l') + else: + eventTypes.append('double_leg') + # Output. squatEvents = { 'eventIdxs': startEndIdxs, 'eventTimes': eventTimes, - 'eventNames':['repStart','repEnd']} + 'eventNames':['repStart','repEnd'], + 'eventTypes': eventTypes} return squatEvents @@ -253,6 +286,147 @@ def compute_range_of_motion(self, coordinate, return_all=False): return ranges_of_motion, units else: return range_of_motion_mean, range_of_motion_std, units + + def compute_squat_depth(self, return_all=False): + pelvis_ty = self.coordinateValues['pelvis_ty'].to_numpy() + + squat_depths = np.zeros((self.nRepetitions)) + for i in range(self.nRepetitions): + rep_range = self.squatEvents['eventIdxs'][i] + + pelvis_ty_range = pelvis_ty[rep_range[0]:rep_range[1]+1] + squat_depths[i] = np.max(pelvis_ty_range) - np.min(pelvis_ty_range) + + # Average across all squats. + squat_depths_mean = np.mean(squat_depths) + squat_depths_std = np.std(squat_depths) + + # Define units. + units = 'm' + + if return_all: + return squat_depths, units + else: + return squat_depths_mean, squat_depths_std, units + + def compute_trunk_lean_relative_to_pelvis(self, return_all=False): + torso_transforms = self.bodyTransformDict['body_transforms']['torso'] + pelvis_transforms = self.bodyTransformDict['body_transforms']['pelvis'] + + max_trunk_leans = np.zeros((self.nRepetitions)) + for i in range(self.nRepetitions): + rep_range = self.squatEvents['eventIdxs'][i] + + torso_transforms_range = torso_transforms[rep_range[0]:rep_range[1]+1] + pelvis_transforms_range = pelvis_transforms[rep_range[0]:rep_range[1]+1] + + trunk_leans_range = np.zeros(len(torso_transforms_range)) + for j in range(len(torso_transforms_range)): + y_axis = opensim.Vec3(0, 1, 0) + torso_y_in_ground = torso_transforms_range[j].xformFrameVecToBase(y_axis).to_numpy() + + z_axis = opensim.Vec3(0, 0, 1) + pelvis_z_in_ground = pelvis_transforms_range[j].xformFrameVecToBase(z_axis).to_numpy() + + acos_deg = np.rad2deg(np.arccos(np.dot(torso_y_in_ground, pelvis_z_in_ground))) + trunk_leans_range[j] = 90 - acos_deg + + # using absolute value for now. perhaps keep track of both positive + # and negative trunk lean angles (to try to detect a bad side?) + max_trunk_leans[i] = np.max(np.abs(trunk_leans_range)) + units = 'deg' + + if return_all: + return max_trunk_leans, units + + else: + trunk_lean_mean = np.mean(max_trunk_leans) + trunk_lean_std = np.std(max_trunk_leans) + return trunk_lean_mean, trunk_lean_std, units + + def compute_trunk_lean_relative_to_ground(self, return_all=False): + torso_transforms = self.bodyTransformDict['body_transforms']['torso'] + pelvis_transforms = self.bodyTransformDict['body_transforms']['pelvis'] + + max_trunk_leans = np.zeros((self.nRepetitions)) + for i in range(self.nRepetitions): + rep_range = self.squatEvents['eventIdxs'][i] + + torso_transforms_range = torso_transforms[rep_range[0]:rep_range[1]+1] + pelvis_transforms_range = pelvis_transforms[rep_range[0]:rep_range[1]+1] + + trunk_leans_range = np.zeros(len(torso_transforms_range)) + for j in range(len(torso_transforms_range)): + y_axis = opensim.Vec3(0, 1, 0) + torso_y_in_ground = torso_transforms_range[j].xformFrameVecToBase(y_axis).to_numpy() + + # find the heading of the pelvis (in the ground x-z plane) + x_axis = opensim.Vec3(1, 0, 0) + pelvis_x_in_ground = pelvis_transforms_range[j].xformFrameVecToBase(x_axis).to_numpy() + pelvis_x_in_ground_xz_plane = pelvis_x_in_ground + pelvis_x_in_ground_xz_plane[1] = 0 + pelvis_x_in_ground_xz_plane /= np.linalg.norm(pelvis_x_in_ground_xz_plane) + + # find the z-axis for comparison to the torso (z-axis in ground + # that is rotated to the heading of the pelvis) + rotated_z_axis = np.cross(pelvis_x_in_ground_xz_plane, np.array([0, 1, 0])) + + acos_deg = np.rad2deg(np.arccos(np.dot(torso_y_in_ground, rotated_z_axis))) + trunk_leans_range[j] = 90 - acos_deg + + # using absolute value for now. perhaps keep track of both positive + # and negative trunk lean angles (to try to detect a bad side?) + max_trunk_leans[i] = np.max(np.abs(trunk_leans_range)) + + units = 'deg' + + if return_all: + return max_trunk_leans, units + + else: + trunk_lean_mean = np.mean(max_trunk_leans) + trunk_lean_std = np.std(max_trunk_leans) + return trunk_lean_mean, trunk_lean_std, units + + def compute_trunk_flexion_relative_to_ground(self, return_all=False): + torso_transforms = self.bodyTransformDict['body_transforms']['torso'] + pelvis_transforms = self.bodyTransformDict['body_transforms']['pelvis'] + + max_trunk_flexions = np.zeros((self.nRepetitions)) + for i in range(self.nRepetitions): + rep_range = self.squatEvents['eventIdxs'][i] + + torso_transforms_range = torso_transforms[rep_range[0]:rep_range[1]+1] + pelvis_transforms_range = pelvis_transforms[rep_range[0]:rep_range[1]+1] + + trunk_flexions_range = np.zeros(len(torso_transforms_range)) + for j in range(len(torso_transforms_range)): + y_axis = opensim.Vec3(0, 1, 0) + torso_y_in_ground = torso_transforms_range[j].xformFrameVecToBase(y_axis).to_numpy() + + # find the heading of the pelvis (in the ground x-z plane) + x_axis = opensim.Vec3(1, 0, 0) + pelvis_x_in_ground = pelvis_transforms_range[j].xformFrameVecToBase(x_axis).to_numpy() + pelvis_x_in_ground_xz_plane = pelvis_x_in_ground + pelvis_x_in_ground_xz_plane[1] = 0 + pelvis_x_in_ground_xz_plane /= np.linalg.norm(pelvis_x_in_ground_xz_plane) + + acos_deg = np.rad2deg(np.arccos(np.dot(torso_y_in_ground, pelvis_x_in_ground_xz_plane))) + trunk_flexions_range[j] = 90 - acos_deg + + # using absolute value for now. perhaps keep track of both positive + # and negative trunk lean angles (to try to detect a bad side?) + max_trunk_flexions[i] = np.max(trunk_flexions_range) + + units = 'deg' + + if return_all: + return max_trunk_flexions, units + + else: + max_trunk_flexions_mean = np.mean(max_trunk_flexions) + max_trunk_flexions_std = np.std(max_trunk_flexions) + return max_trunk_flexions_mean, max_trunk_flexions_std, units def get_coordinates_segmented(self): @@ -325,3 +499,4 @@ def get_center_of_mass_segmented_normalized_time(self): comValuesTimeNormalized['indiv'] = [pd.DataFrame(data=d, columns=colNames) for d in comValuesSegmentedNorm] return comValuesTimeNormalized + \ No newline at end of file