diff --git a/unityvr/analysis/posAnalysis.py b/unityvr/analysis/posAnalysis.py index 9d25848..5b292a0 100644 --- a/unityvr/analysis/posAnalysis.py +++ b/unityvr/analysis/posAnalysis.py @@ -14,7 +14,7 @@ ##functions to process posDf dataframe #obtain the position dataframe with derived quantities -def position(uvrDat, derive = True, rotate_by = None, filter_date = '2021-09-08', plot = False, plotsave=False, saveDir=None): +def position(uvrDat, derive = True, rotate_by = None, filter_date = '2021-09-08', plot = False, plotsave=False, saveDir=None, computeVel=False, **computeVelocitiesKwargs): ## input arguments # set derive = True if you want to compute derived quantities (ds, s, dTh (change in angle), radangle (angle in radians(-pi,pi))) # rotate_by: angle (degrees) by which to rotate the trajectory to ensure the bright part of the panorama is at 180 degree heading. @@ -22,7 +22,7 @@ def position(uvrDat, derive = True, rotate_by = None, filter_date = '2021-09-08' # converting from Unity's left handed angle convention to right handed convention was implemented after a certain # date in the preproc.py file - posDf = uvrDat.posDf + posDf = uvrDat.posDf.copy() #NOT INPLACE #angle correction #this is required only for data that was preprocessed before the filter_date @@ -33,9 +33,10 @@ def position(uvrDat, derive = True, rotate_by = None, filter_date = '2021-09-08' #rotate if rotate_by is not None: + + #rotate the trajectory posDf['x'], posDf['y'] = rotation_deg(posDf['x'],posDf['y'],rotate_by) - posDf['dx'], posDf['dy'] = rotation_deg(posDf['dx'],posDf['dy'],rotate_by) - posDf['dxattempt'], posDf['dyattempt'] = rotation_deg(posDf['dxattempt'],posDf['dyattempt'],rotate_by) + posDf['angle'] = (posDf['angle']+rotate_by)%360 uvrDat.metadata['rotated_by'] = (uvrDat.metadata['rotated_by']+rotate_by)%360 if ('rotated_by' in uvrDat.metadata) else (rotate_by%360) @@ -43,10 +44,9 @@ def position(uvrDat, derive = True, rotate_by = None, filter_date = '2021-09-08' posDf.dc2cm = 10 if derive: - posDf['ds'] = np.sqrt(posDf['dx']**2+posDf['dy']**2) - posDf['s'] = np.cumsum(posDf['ds']) - posDf['dTh'] = (np.diff(posDf['angle'],prepend=posDf['angle'].iloc[0]) + 180)%360 - 180 - posDf['radangle'] = ((posDf['angle']+180)%360-180)*np.pi/180 + posDf = posDerive(posDf) + if computeVel: + posDf = computeVelocities(posDf,**computeVelocitiesKwargs) if plot: fig, ax = viz.plotTrajwithParameterandCondition(posDf, figsize=(10,5), parameter='angle') @@ -119,17 +119,16 @@ def rotation_deg(x,y,theta): # Add derrived quantities: Velocities -def computeVelocities(posDf, convf = 10, window=7, order=3): - #convf: conversion from dc (to cm if 10) +def computeVelocities(inDf, window=7, order=3): + posDf = inDf.copy() #window and order for filter - from scipy.signal import savgol_filter # add derrived parameter to positional dataframe - posDf['vT'] = np.hypot(np.gradient(posDf.x.values*convf), np.gradient(posDf.y.values*convf))*(1/posDf.dt) - posDf['vR'] = np.gradient(np.unwrap(posDf.angle.values)) + posDf['vT'] = np.hypot(np.gradient(posDf.x.values), np.gradient(posDf.y.values))*(1/posDf.dt) #translational velocity in decimeter/second + posDf['vR'] = np.gradient(np.unwrap(posDf.angle.values)) #rotational velocity in degrees/second - posDf['vT_filt'] = savgol_filter(posDf.vT, window, order) - posDf['vR_filt'] = savgol_filter(posDf.vR, window, order) + posDf['vT_filt'] = savgolFilterInterpolated(posDf, 'vT', window, order) #savgol filtered velocities + posDf['vR_filt'] = savgolFilterInterpolated(posDf, 'vR', window, order) return posDf def getTimeDf(uvrDat, trialDir, posDf = None, imaging = False, rate = 9.5509): @@ -159,4 +158,38 @@ def getTimeDf(uvrDat, trialDir, posDf = None, imaging = False, rate = 9.5509): timeDf['count'] = np.arange(1,len(timeDf['time'])+1,1) timeDf = carryAttrs(timeDf,posDf) - return timeDf \ No newline at end of file + return timeDf + +def posDerive(inDf): + posDf = inDf.copy() #NOT INPLACE + posDf['dx'] = np.diff(posDf['x'], prepend=0) #allocentric translation vector x component + posDf['dy'] = np.diff(posDf['y'], prepend=0) #allocentric translation vector y component + posDf['ds'] = np.sqrt(posDf['dx']**2+posDf['dy']**2) #magnitude of the translation vector + posDf['s'] = np.cumsum(posDf['ds']) #integrated pathlength from start + posDf['dTh'] = (np.diff(posDf['angle'],prepend=posDf['angle'].iloc[0]) + 180)%360 - 180 + posDf['radangle'] = ((posDf['angle']+180)%360-180)*np.pi/180 + + #derive forward and side velocities: does not depend on rotation of trajectory and angle + if 'dx_ft' not in posDf: + xy = np.diff(posDf[['x', 'y']], axis=0) + rotation_mats = np.array([np.array([[np.cos(theta), -np.sin(theta)], + [np.sin(theta), np.cos(theta)]]).T for theta in np.deg2rad(posDf['angle'])[:-1]]) + posDf['dx_ft'], posDf['dy_ft'] = np.vstack([[0,0],np.einsum('ijk,ik->ij', rotation_mats, xy)]).T #in egocentric frame, units: decimeters, #dx_ft is forward motion (forward is positive), dy_ft is side motion (rightward is positive) + return posDf + +def savgolFilterInterpolated(posDf, value, window, order, timeStr = 'time', kind='linear'): + from scipy.signal import savgol_filter + + df = posDf.copy() + + #filtering with a savgol filter assumes that the data is evenly spaced + #interpolating the data to be evenly spaced + interpolated = sp.interpolate.interp1d(df[timeStr], df[value], kind=kind, bounds_error=False, fill_value='extrapolate') + + equalTime = np.linspace(df[timeStr].min(),df[timeStr].max(),2*len(df[timeStr])) + + filt = savgol_filter(interpolated(equalTime), window, order) + + filtInterpolated = sp.interpolate.interp1d(equalTime, filt, kind='nearest', bounds_error=False, fill_value='extrapolate') + + return filtInterpolated(df.time) \ No newline at end of file diff --git a/unityvr/preproc/logproc.py b/unityvr/preproc/logproc.py index e9cf95f..2baa686 100644 --- a/unityvr/preproc/logproc.py +++ b/unityvr/preproc/logproc.py @@ -11,8 +11,9 @@ objDfCols = ['name','collider','px','py','pz','rx','ry','rz','sx','sy','sz'] posDfCols = ['frame','time','x','y','angle'] -ftDfCols = ['frame','ficTracTReadMs','ficTracTWriteMs','dx','dy','dz'] +ftDfCols = ['frame','ficTracTReadMs','ficTracTWriteMs','wx_ft','wy_ft','wz_ft'] dtDfCols = ['frame','time','dt'] +tempDfCols = ['frame','time','temperature'] nidDfCols = ['frame','time','dt','pdsig','imgfsig'] texDfCols = ['frame','time','xtex','ytex'] vidDfCols = ['frame','time','img','duration'] @@ -266,10 +267,10 @@ def posDfFromLog(dat): 'x': match['worldPosition']['x'], 'y': match['worldPosition']['z'], #axes are named differently in Unity 'angle': (-match['worldRotationDegs']['y'])%360, #flip due to left handed convention in Unity - 'dx':match['actualTranslation']['x'], - 'dy':match['actualTranslation']['z'], - 'dxattempt': match['attemptedTranslation']['x'], - 'dyattempt': match['attemptedTranslation']['z'] + 'dx_ft':match['actualTranslation']['x'], #forward motion + 'dy_ft':match['actualTranslation']['z'], #right sideslip + 'dxattempt_ft': match['attemptedTranslation']['x'], + 'dyattempt_ft': match['attemptedTranslation']['z'] } entries[entry] = pd.Series(framedat).to_frame().T print('correcting for Unity angle convention.') @@ -288,9 +289,9 @@ def ftDfFromLog(dat): framedat = {'frame': match['frame'], 'ficTracTReadMs': match['ficTracTimestampReadMs'], 'ficTracTWriteMs': match['ficTracTimestampWriteMs'], - 'dx': match['ficTracDeltaRotationVectorLab']['x'], - 'dy': match['ficTracDeltaRotationVectorLab']['y'], - 'dz': match['ficTracDeltaRotationVectorLab']['z']} + 'wx_ft': match['ficTracDeltaRotationVectorLab']['x'], + 'wy_ft': match['ficTracDeltaRotationVectorLab']['y'], + 'wz_ft': match['ficTracDeltaRotationVectorLab']['z']} entries[entry] = pd.Series(framedat).to_frame().T if len(entries) > 0: