from __future__ import division import os, sys # import seaborn from pylab import rcParams import cv2 import numpy as np from numpy import linalg as LA from time import time from itertools import combinations import matplotlib.pyplot as plt from matplotlib.pyplot import legend import matplotlib.patches as mpatches from minimize import findInitialW, _q, g, minimizeEnergy, g3D3D, findW3D3D from minimize import g as gaze_ray from geom import getSphericalCoords, getAngularDiff from vector import Vector as v from sim import GazeSimulation from recording.util.tools import is_outlier from recording.tracker import readCameraParams from geom import getRotationMatrix import pdb PARTICIPANTS = ['p10', 'p16', 'p13', 'p24', 'p5', 'p14', 'p26', 'p12', 'p20', 'p7', 'p15', 'p11', 'p21', 'p25'] class Experiment: def performExperiment(self): print 'Running Experiment...' start = time() self.__run__() self.runningTime = time() - start print 'Running Time:', self.runningTime def __run__(self): pass class Parallax2Dto2DMapping(Experiment): def __run__(self): sim = GazeSimulation(log = False) sim.place_eyeball_on_scene_camera = False sim.setEyeRelativeToSceneCamera(v(-65, -33, -73)) sim.setCalibrationDepth(1 * 1000) # mm, wrt scene camera sim.setTestDepth(1.5 * 1000) sim.calibration_grid = True sim.calibration_random_depth = False sim.test_grid = True sim.test_random_depth = False sim.test_random_fixed_depth = False sim.reset() print 'scene_camera', sim.scene_camera.t print 'calibration' print len(sim.calibration_points) print min(np.array(sim.calibration_points)[:,0]) - sim.scene_camera.t[0], max(np.array(sim.calibration_points)[:,0]) - sim.scene_camera.t[0] print min(np.array(sim.calibration_points)[:,1]) - sim.scene_camera.t[1], max(np.array(sim.calibration_points)[:,1]) - sim.scene_camera.t[1] print 'depths', set(np.array(sim.calibration_points)[:,2]) print len(sim.test_points) print min(np.array(sim.test_points)[:,0]) - sim.scene_camera.t[0], max(np.array(sim.test_points)[:,0]) - sim.scene_camera.t[0] print min(np.array(sim.test_points)[:,1]) - sim.scene_camera.t[1], max(np.array(sim.test_points)[:,1]) - sim.scene_camera.t[1] print 'depths', set(np.array(sim.test_points)[:,2]) plt.ylabel('Y (mm)') plt.xlabel('X (mm)') plt.plot(np.array(sim.calibration_points)[:,0], np.array(sim.calibration_points)[:,1], 'bo') plt.plot(np.array(sim.test_points)[:,0], np.array(sim.test_points)[:,1], 'ro') plt.show() # self.sim = sim # visualize the setting in windows class Parallax2Dto3DMapping(Experiment): ''' IMPORTANT! In all experiments, scene camera's rvec = (0, 0, 0) i.e. the corresponding rotation matrix is the identity matrix therefore I have not included the dot production with this rotation matrix to convert points in world coordinates into scene camera coordinates. however, one should know that if the scene camera is rotated differently this transformation is essential. ''' def __run__(self): sim = GazeSimulation(log = False) sim.place_eyeball_on_scene_camera = False sim.setEyeRelativeToSceneCamera(v(-65, -33, -73)) # based on actual measurements in our study sim.setCalibrationDepth(1 * 1000) # mm, wrt scene camera sim.setTestDepth(1.5 * 1000) sim.calibration_grid = True sim.calibration_random_depth = False sim.test_grid = True sim.test_random_depth = False sim.test_random_fixed_depth = False depths = map(lambda d:d*1000, [1, 1.25, 1.5, 1.75, 2.0]) print '> Computing results for multiple calibration depths...' results, results_std = [], [] for num_of_calibration_depths in xrange(1, 6): # from 1 calibration depths to 5 print '> Considering only %s calibration depth(s)...' %num_of_calibration_depths sim.reset() aae_2ds_aae = [] aae_2ds_phe = [] aae_3ds_aae = [] aae_3ds_phe = [] aae_3D3Ds = [] # angular error for calibs in combinations(depths, num_of_calibration_depths): # Now calibs is a set of depths from each of which we need calibration data print 'Current calibration depths', calibs calibs = list(calibs) cp, ct = [], [] sim.reset() sim.setCalibrationDepth(calibs) # Perform calibration sim.runCalibration() cp, ct, p3d = sim.tr_pupil_locations, sim.calibration_points, sim.tr_3d_pupil_locations # target positions are computed relative to the scene CCS ti = map(lambda target: v(target) - v(sim.scene_camera.t), ct) # Computing pupil pose for each gaze ni = map(lambda p: (v(p)-v(sim.sclera_pos)).norm(), p3d) # ground truth gaze vectors w, e, w0 = minimizeEnergy(cp, ti) e = v(e) # transforming pupil pose to eye camera CS eyeR = np.array(sim.eye_camera.R[:3]) ni = map(lambda pose: eyeR.dot(np.array(pose)), ni) R, e3d3d = minimizeEnergy(ni, ti, pose_given=True) # R = LA.inv(R) e3d3d = v(e3d3d) # Now we have calibration data from multiple depths, we can test on all depths for test_depth in depths: sim.setTestDepth(test_depth) aae_2d_aae, aae_2d_phe, aae_2d_std, _ = sim.runTest() # last one is PHE std aae_2ds_aae.append((aae_2d_aae, aae_2d_std)) aae_2ds_phe.append(aae_2d_phe) # Fetching test points t, p, p3d = sim.test_points, sim.te_pupil_locations, sim.te_3d_pupil_locations t = map(lambda target: v(target) - v(sim.scene_camera.t), t) # target coords in scene CCS # 3D3D t_3d3d = t[:] ni = map(lambda p: v(v(p)-v(sim.sclera_pos)).norm(), p3d) # ground truth gaze vectors # transforming pupil pose to eye camera CS ni = map(lambda r: v(eyeR.dot(np.array(r))), ni) # applying estimated rotation to pose vector in eye camera coordinates (Rn) # R is estimated rotation between scene camera and eye coordinate system (not eye camera!) # in other words, R is the rotation part of e Rni = map(lambda n: v(R.dot(np.array(n))), ni) # now ready to compare Rn with t-e # Intersecting gaze rays originating from the eye with the planes defined by each # target. then we can simply compute angular error between each intersection and # the corresponding 3D target gis = map(lambda vec: v(vec), Rni) # gaze rays originating from eyeball # we multiply g such that it hits t's z-plane i.e. multiply all coordinates by factor (t.z-e.z)/g.z # then we add e to the final g so that it originates from scene camera. now both g and t are in the # same coordinate system and originate from the same point, so we can compare them gprimes = map(lambda tg: v(((tg[0].z - e3d3d.z)/tg[1].z)*tg[1] + e3d3d), zip(t_3d3d, gis)) AE = list(np.degrees(np.arctan((v(p[0]).cross(p[1])/(v(p[0]).dot(p[1]))).mag)) for p in zip(gprimes, t_3d3d)) N = len(t) AAE = np.mean(AE) STD = np.std(AE) m, M = min(AE), max(AE) aae_3D3Ds.append((AAE, STD)) qi = map(_q, p) # computing feature vectors from raw pupil coordinates in 2D # computing unit gaze vectors corresponding to pupil positions # here we use the computed mapping matrix w gis = map(lambda q: g(q, w), qi) # Intersecting gaze rays originating from the eye with the planes defined by each # target. then we can simply compute angular error between each intersection and # the corresponding 3D target t = map(lambda vec: v(vec), t) gis = map(lambda vec: v(vec), gis) gprimes = map(lambda tg: v(((tg[0].z - e.z)/tg[1].z)*tg[1] + e), zip(t, gis)) AE = list(np.degrees(np.arctan((v(p[0]).cross(p[1])/(v(p[0]).dot(p[1]))).mag)) for p in zip(gprimes, t)) N = len(t) AAE = np.mean(AE) STD = np.std(AE) m, M = min(AE), max(AE) # Computing physical distance error (in meters) PHE = list((u-v).mag/1000 for u,v in zip(t, gprimes)) N = len(t) APHE = np.mean(PHE) PHE_STD = np.std(PHE) PHE_m, PHE_M = min(PHE), max(PHE) aae_3ds_aae.append((AAE, STD)) aae_3ds_phe.append((PHE, PHE_STD)) # results only contains AAE results.append([np.mean(np.array(aae_2ds_aae)[:,0]), np.mean(np.array(aae_3ds_aae)[:,0]), np.mean(np.array(aae_3D3Ds)[:,0])]) results_std.append([np.std(np.array(aae_2ds_aae)[:,0]), np.std(np.array(aae_3ds_aae)[:,0]), np.std(np.array(aae_3D3Ds)[:,0])]) ###################################################################################################### # plot code based on EffectNumberofClusters.py mean2D2D = [res[0] for res in results] mean2D3D = [res[1] for res in results] mean3D3D = [res[2] for res in results] std2D2D = [res[0] for res in results_std] std2D3D = [res[1] for res in results_std] std3D3D = [res[2] for res in results_std] N = 5 ind = np.asarray([0.25,1.25,2.25,3.25,4.25]) width = 0.5 # the width of the bars x2 = [0.45,1.45,2.45,3.45,4.45] x4 = [0.55,1.55,2.55,3.55,4.55] x6 = [0.50,1.50,2.50,3.50,4.50] fig = plt.figure(figsize=(14.0, 10.0)) ax = fig.add_subplot(111) rects1 = ax.errorbar(x2, mean2D2D,yerr=[std2D2D,std2D2D],fmt='o',color='red',ecolor='red',lw=3, capsize=5, capthick=2) plt.plot(x2, mean2D2D, marker="o", linestyle='-',lw=3,color='red',label = r'2D-to-2D') rects2 =ax.errorbar(x4, mean2D3D,yerr=[std2D3D,std2D3D],fmt='o',color='blue',ecolor='blue',lw=3, capsize=5, capthick=2) plt.plot(x4, mean2D3D, marker="o", linestyle='-',lw=3,color='blue', label = r'2D-to-3D') rects3 =ax.errorbar(x6, mean3D3D,yerr=[std3D3D,std3D3D],fmt='o',color='blue',ecolor='blue',lw=3, capsize=5, capthick=2) plt.plot(x6, mean3D3D, marker="o", linestyle='-',lw=3,color='blue', label = r'3D-to-3D') legend(fontsize=20,loc='upper right') ax.set_ylabel(r'Angular Error',fontsize=22) ax.set_xlabel(r'Number of Calibration Depths',fontsize=22) ax.set_xticks(ind+0.25) ax.set_xticklabels( ('D1', 'D2', 'D3','D4', 'D5') ,fontsize=18) TOPICs = [0.0,0.5,1.5,2.5,3.5,4.5,5.0] print TOPICs LABELs = ["",r'1',r'2', r'3', r'4', r'5', ""] plt.xticks(TOPICs, LABELs,fontsize=18) TOPICS = [0.5,1,1.5,2,2.5,3,3.5,4,4.5,5] print TOPICS LABELS = [r'0.5', r'1',r'1.5', r'2',r'2.5', r'3',r'3.5', r'4',r'4.5',r'5'] plt.yticks(TOPICS, LABELS,fontsize=18) def autolabel(rects): # attach some text labels for rect in rects: height = rect.get_height() ax.text(0.26+rect.get_x()+rect.get_width()/2., height +0.35, "%.2f"%float(height), ha='center', va='bottom',fontweight='bold',fontsize=13.5) left = 0.1 # the left side of the subplots of the figure right = 0.975 # the right side of the subplots of the figure bottom = 0.075 # the bottom of the subplots of the figure top = 0.925 # the top of the subplots of the figure wspace = 0.2 # the amount of width reserved for blank space between subplots hspace = 0.4 # the amount of height reserved for white space between subplots plt.subplots_adjust(left=left, bottom=bottom, right=right, top=top, wspace=wspace, hspace=hspace) plt.show() ###################################################################################################### self.sim = sim class Parallax3Dto3DMapping(Experiment): # GT pupil pose instead of estimating the pose def __run__(self): sim = GazeSimulation(log = False) sim.place_eyeball_on_scene_camera = False sim.setEyeRelativeToSceneCamera(v(-65, -33, -73)) sim.setCalibrationDepth(1 * 1000) # mm, wrt scene camera sim.setTestDepth(1.5 * 1000) sim.calibration_grid = True sim.calibration_random_depth = False sim.test_grid = True sim.test_random_depth = False sim.test_random_fixed_depth = False cdepths = [1, 1.5, 2.0] tdepths = np.linspace(1, 2, 5) results = [] for i, cdepth in enumerate(cdepths): sim.setCalibrationDepth(cdepth * 1000) # Performing calibration sim.runCalibration() aae_2Ds = [] aae_3Ds = [] aae_3D3Ds = [] std_2Ds = [] std_3Ds = [] std_3D3Ds = [] # Fetching calibration points t, p3d, p2d = sim.calibration_points, sim.tr_3d_pupil_locations, sim.tr_pupil_locations # target positions are computed relative to the scene CCS ti = map(lambda target: v(target) - v(sim.scene_camera.t), t) # Computing pupil pose for each gaze ni = map(lambda p: (v(p)-v(sim.sclera_pos)).norm(), p3d) # ground truth gaze vectors eye_scene_diff = v(sim.sclera_pos) - v(sim.scene_camera.t) e = np.array(eye_scene_diff) # eyeball coords in scene CCS e_org = v(e[:]) # transforming pose vector to eye camera coordinates eyeR = np.array(sim.eye_camera.R[:3]) # result is identical to cv2.Rodrigues(np.array([0., np.pi, 0.]))[0] ni = map(lambda pose: eyeR.dot(np.array(pose)), ni) R, e3d3d = minimizeEnergy(ni, ti, pose_given=True) e3d3d = v(e3d3d) # R = LA.inv(R) w2d3d, e2d3d, w0 = minimizeEnergy(p2d, ti) e2d3d = v(e2d3d) print print R print 'e3d3d', e3d3d, 'distance =', (e3d3d-e_org).mag print 'e2d3d', e2d3d, 'distance =', (e2d3d-e_org).mag print 'real e', e_org print for j, tdepth in enumerate(tdepths): sim.setTestDepth(tdepth * 1000) aae_2d, _, std_2d, _ = sim.runTest() aae_2Ds.append(aae_2d) std_2Ds.append(std_2d) # Fetching test points t, p3d, p2d = sim.test_points, sim.te_3d_pupil_locations, sim.te_pupil_locations ti = map(lambda target: v(target) - v(sim.scene_camera.t), t) # target coords in scene CCS ni = map(lambda p: v(v(p)-v(sim.sclera_pos)).norm(), p3d) # ground truth gaze vectors # transforming pose vector to eye camera coordinates ni = map(lambda r: v(eyeR.dot(np.array(r))), ni) # applying estimated rotation to pose vector in eye camera coordinates (Rn) ni = map(lambda n: v(R.dot(np.array(n))), ni) # now ready to compare Rn with t-e # 3D to 3D # Intersecting gaze rays originating from the eye with the planes defined by each # target. then we can simply compute angular error between each intersection and # the corresponding 3D target gis = map(lambda vec: v(vec), ni) gprimes = map(lambda tg: v(((tg[0].z - e3d3d.z)/tg[1].z)*tg[1] + e3d3d), zip(ti, gis)) AE = list(np.degrees(np.arctan((v(p[0]).cross(p[1])/(v(p[0]).dot(p[1]))).mag)) for p in zip(gprimes, ti)) # AE = list(np.degrees(np.arccos(v(p[0]).dot(p[1])/v(p[0]).mag/v(p[1]).mag)) for p in zip(gprimes, ti)) N = len(ti) AAE = np.mean(AE) STD_3D3D = np.std(AE) m, M = min(AE), max(AE) aae_3D3Ds.append(AAE) std_3D3Ds.append(STD_3D3D) # 2D to 3D qi = map(_q, p2d) # computing feature vectors from raw pupil coordinates in 2D # computing unit gaze vectors corresponding to pupil positions # here we use the computed mapping matrix w gis = map(lambda q: gaze_ray(q, w2d3d), qi) # Intersecting gaze rays originating from the eye with the planes defined by each # target. then we can simply compute angular error between each intersection and # the corresponding 3D target gis = map(lambda vec: v(vec), gis) gprimes = map(lambda tg: v(((tg[0].z - e2d3d.z)/tg[1].z)*tg[1] + e2d3d), zip(ti, gis)) AE = list(np.degrees(np.arctan((v(p[0]).cross(p[1])/(v(p[0]).dot(p[1]))).mag)) for p in zip(gprimes, ti)) N = len(t) AAE = np.mean(AE) STD_2D3D = np.std(AE) m, M = min(AE), max(AE) # Computing physical distance error (in meters) PHE = list((u-v).mag for u,v in zip(ti, gprimes)) N = len(ti) APHE = np.mean(PHE) PHE_STD = np.std(PHE) PHE_m, PHE_M = min(PHE), max(PHE) aae_3Ds.append(AAE) std_3Ds.append(STD_2D3D) print 'depth', cdepth, 'finished.' results.append([aae_2Ds, aae_3Ds, aae_3D3Ds, std_2Ds, std_3Ds, std_3D3Ds]) ###################################################################################################### ## Plotting part clrs = ['r', 'g', 'b', 'k', 'o'] colors = ['blue', 'orange', 'red', 'black', 'orange'] patches = [] fig = plt.figure(figsize=(14.0, 10.0)) ax = fig.add_subplot(111) dmap = {0: '1', 1:'3', 2:'5'} x1 = [0.375,1.375,2.375,3.375,4.375] x2 = [0.425,1.425,2.425,3.425,4.425] x3 = [0.475,1.475,2.475,3.475,4.475] x4 = [0.525,1.525,2.525,3.525,4.525] x5 = [0.575,1.575,2.575,3.575,4.575] x6 = [0.625,1.625,2.625,3.625,4.625] xrange_2d2d = [x1, x3, x5] xrange_2d3d = [x2, x4, x6] for i in [0, 1, 2]: cdepth_results = results[i] _xrange = xrange_2d2d[i] aae_2d2d = cdepth_results[0] std_2d2d = cdepth_results[3] rects1 = ax.errorbar(_xrange, aae_2d2d,yerr=[std_2d2d, std_2d2d],fmt='o',color=colors[i],ecolor=colors[i],lw=3, capsize=5, capthick=2) plt.plot(_xrange, aae_2d2d, marker="o", linestyle='-',lw=3,color=colors[i],label = '2D-to-2D Calibration Depth ' + dmap[i]) for i in [0, 1, 2]: cdepth_results = results[i] _xrange = xrange_2d3d[i] aae_2d3d = cdepth_results[1] std_2d3d = cdepth_results[4] rects2 = ax.errorbar(_xrange, aae_2d3d,yerr=[std_2d3d, std_2d3d],fmt='o',color=colors[i],ecolor=colors[i],lw=3, capsize=5, capthick=2) plt.plot(_xrange, aae_2d3d, marker="o", linestyle='--',lw=3,color=colors[i],label = '2D-to-3D Calibration Depth ' + dmap[i]) for i in [0, 1, 2]: cdepth_results = results[i] _xrange = xrange_2d2d[i] aae_3d3d = cdepth_results[2] std_3d3d = cdepth_results[5] rects3 = ax.errorbar(_xrange, aae_3d3d,yerr=[std_3d3d, std_3d3d],fmt='o',color=colors[i],ecolor=colors[i],lw=3, capsize=5, capthick=2) plt.plot(_xrange, aae_3d3d, marker="o", linestyle='-.',lw=3,color=colors[i],label = '3D-to-3D Calibration Depth ' + dmap[i]) ax.set_ylabel(r'\textbf{Angular Error}',fontsize=22) ax.set_xlabel(r'\textbf{Depth}',fontsize=22) TOPICS = [-0.2, 0, 0.2, 0.4,0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4]#,110]#,120] LABELS = [r'', r'0', r'0.2',r'0.4',r'0.6', r'0.8', r'1.0', r'1.2', r'1.4', r'1.6', r'1.8', r'2.0', r'2.2', r'2.4']#, ""]#, ""] plt.yticks(TOPICS, LABELS,fontsize=18) plt.xlabel('Depth') plt.ylabel('Angular Error') plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, ncol=3, mode="expand", borderaxespad=0., fontsize=18) TOPICs = [0.0,0.5,1.5,2.5,3.5,4.5,5.0] LABELs = ['', 'D1 - 1m', 'D2 - 1.25m', 'D3 - 1.5m', 'D4 - 1.75m', 'D5 - 2.0m', ''] plt.xticks(TOPICs, LABELs,fontsize=18) left = 0.1 # the left side of the subplots of the figure right = 0.975 # the right side of the subplots of the figure bottom = 0.075 # the bottom of the subplots of the figure top = 0.925 # the top of the subplots of the figure wspace = 0.2 # the amount of width reserved for blank space between subplots hspace = 0.4 # the amount of height reserved for white space between subplots plt.subplots_adjust(left=left, bottom=bottom, right=right, top=top, wspace=wspace, hspace=hspace) plt.show() ###################################################################################################### # self.sim = sim # This should be root path to the participants folder (including subfolders for each participant) ROOT_DATA_DIR = '/home/mmbrian/HiWi/etra2016_mohsen/code/recording/data/participants' class Parallax2Dto2DRealData(Experiment): ''' First it runs single calibration depth vs test depth over all combinations of depths per participant (25 cases) Then we perform estimation with data from multiple calibration depths for each test depth ''' def __run__(self): sim = GazeSimulation(log = False) camera_matrix, dist_coeffs = readCameraParams() root_result_path = 'results/2d2d/' if not os.path.exists(root_result_path): os.makedirs(root_result_path) errors = [] for d1 in os.listdir(ROOT_DATA_DIR): if d1.startswith('p'): # every participant if not d1 in PARTICIPANTS: continue participant_label = d1 participant_results = [] d2 = os.path.join(ROOT_DATA_DIR, d1) # .../pi/ d2 = os.path.join(d2, os.listdir(d2)[0]) # .../pi/../ print '> Processing participant', d1 participant_experiment = {} for d3 in os.listdir(d2): # every recording d4 = os.path.join(d2, d3) # .../pi/../00X/ intervals_dir = os.path.join(d4, 'gaze_intervals.npy') if not os.path.isfile(intervals_dir): # Participant not completely processed print '> One or more recordings need processing...' break else: p = np.load(os.path.join(d4, 'p.npy')) break_participant = False for point in p: if np.isnan(point[0]): # For this gaze position, not input with nonzero confidence exist break_participant = True break if break_participant: print '> One or more recordings miss gaze coords...' break # print p t2d = np.load(os.path.join(d4, 't2d.npy')) t3d = np.load(os.path.join(d4, 't3d.npy')) participant_experiment[d3] = [p, t2d, t3d] if len(participant_experiment) == 10: print '> All recordings processed' keys = sorted(participant_experiment.keys()) depths = zip(keys[::2], keys[1::2]) for calib_depth in depths: cp, ct = participant_experiment[calib_depth[0]][0], participant_experiment[calib_depth[0]][1] cdepth_value = getDepth(calib_depth) # Perform calibration sim.perform2D2DCalibrationOnReadData(cp, ct, (1280, 720)) for test_depth in depths: tdepth_value = getDepth(test_depth) tp, tt, tt3d = participant_experiment[test_depth[1]][0], participant_experiment[test_depth[1]][1], participant_experiment[test_depth[1]][2] error = sim.run2D2DTestOnRealData(tp, tt, (1280, 720), tt3d, camera_matrix, dist_coeffs) participant_results.append([cdepth_value, tdepth_value] + error) print len(participant_results), 'combinations processed...' np.save('results/2d2d/%s_2d2d_all.npy' % participant_label, np.array(participant_results)) np.savetxt('results/2d2d/%s_2d2d_all.csv' % participant_label, np.array(participant_results), delimiter=",") print '> Computing results for multiple calibration depths...' for num_of_calibration_depths in xrange(2, 6): # from 2 calibration depths to 5 participant_results = [] print '> Computing results for combining %s calibration depths...' % num_of_calibration_depths for calibs in combinations(depths, num_of_calibration_depths): # Now calibs is a set of depths, from each of those we need calibration data cp, ct = [], [] calib_depths_label = [] for calib in calibs: if len(cp): cp = np.concatenate((cp, participant_experiment[calib[0]][0]), axis=0) ct = np.concatenate((ct, participant_experiment[calib[0]][1]), axis=0) else: cp = participant_experiment[calib[0]][0] ct = participant_experiment[calib[0]][1] calib_depths_label.append(getDepth(calib)) # Perform calibration sim.perform2D2DCalibrationOnReadData(cp, ct, (1280, 720)) # Now we have calibration data from multiple depths, we can test on all depths for test_depth in depths: tdepth_value = getDepth(test_depth) tp, tt, tt3d = participant_experiment[test_depth[1]][0], participant_experiment[test_depth[1]][1], participant_experiment[test_depth[1]][2] error = sim.run2D2DTestOnRealData(tp, tt, (1280, 720), tt3d, camera_matrix, dist_coeffs) participant_results.append(calib_depths_label + [tdepth_value] + error) print len(participant_results), 'combinations processed...' result_path = 'results/2d2d/%s_calibration_depths/' % num_of_calibration_depths if not os.path.exists(result_path): os.makedirs(result_path) np.save(result_path + '%s.npy' % participant_label, np.array(participant_results)) np.savetxt(result_path + '%s.csv' % participant_label, np.array(participant_results), delimiter=",") print 'done.' def getDepth(depth_experiments): _map = {'000': 1, '002': 1.25, '004': 1.5, '006': 1.75, '008': 2.0} return _map[depth_experiments[0]] class Parallax2Dto3DRealData(Experiment): def __run__(self): sim = GazeSimulation(log = False) aae_3ds = [] # Path to whatever directory you want the results to be stored in root_result_path = '/home/mmbrian/3D_Gaze_Tracking/work/results/2D3D/' if not os.path.exists(root_result_path): os.makedirs(root_result_path) for d1 in os.listdir(ROOT_DATA_DIR): if d1.startswith('p'): # every participant if not d1 in PARTICIPANTS: continue participant_label = d1 participant_results = [] d2 = os.path.join(ROOT_DATA_DIR, d1) # .../pi/ d2 = os.path.join(d2, os.listdir(d2)[0]) # .../pi/../ print '> Processing participant', d1 participant_experiment = {} for d3 in os.listdir(d2): # every recording d4 = os.path.join(d2, d3) # .../pi/../00X/ intervals_dir = os.path.join(d4, 'gaze_intervals.npy') if not os.path.isfile(intervals_dir): # Participant not completely processed print '> One or more recordings need processing...' break else: p = np.load(os.path.join(d4, 'p.npy')) # p = np.load(os.path.join(d4, 'p_mean.npy')) break_participant = False for point in p: if np.isnan(point[0]): # For this gaze position, not input with nonzero confidence exist break_participant = True break if break_participant: print '> One or more recordings miss gaze coords...' break # print p t2d = np.load(os.path.join(d4, 't2d.npy')) t3d = np.load(os.path.join(d4, 't3d.npy')) # t2d = np.load(os.path.join(d4, 't2d_mean.npy')) # t3d = np.load(os.path.join(d4, 't3d_mean.npy')) participant_experiment[d3] = [p, t2d, t3d] if len(participant_experiment) == 10: print '> All recordings processed' keys = sorted(participant_experiment.keys()) depths = zip(keys[::2], keys[1::2]) for calib_depth in depths: cp, ct3d = participant_experiment[calib_depth[0]][0], participant_experiment[calib_depth[0]][2] cdepth_value = getDepth(calib_depth) # Performing calibration w, e, w0 = minimizeEnergy(cp, ct3d) e = v(e) for test_depth in depths: tdepth_value = getDepth(test_depth) tp, tt3d = participant_experiment[test_depth[1]][0], participant_experiment[test_depth[1]][2] t, p = tt3d, tp qi = map(_q, p) # computing feature vectors from raw pupil coordinates in 2D # computing unit gaze vectors corresponding to pupil positions # here we use the computed mapping matrix w gis = map(lambda q: g(q, w), qi) # Intersecting gaze rays originating from the eye with the planes defined by each # target. then we can simply compute angular error between each intersection and # the corresponding 3D target t = map(lambda vec: v(vec), t) gis = map(lambda vec: v(vec), gis) gprimes = map(lambda tg: v(((tg[0].z - e.z)/tg[1].z)*tg[1] + e), zip(t, gis)) AE = list(np.degrees(np.arctan((v(p[0]).cross(p[1])/(v(p[0]).dot(p[1]))).mag)) for p in zip(gprimes, t)) N = len(t) AAE = sum(AE)/N VAR = sum((ae - AAE)**2 for ae in AE)/N STD = np.sqrt(VAR) m, M = min(AE), max(AE) # Computing physical distance error (in meters) PHE = list((u-v).mag for u,v in zip(t, gprimes)) N = len(t) APHE = sum(PHE)/N PHE_VAR = sum((phe - APHE)**2 for phe in PHE)/N PHE_STD = np.sqrt(VAR) PHE_m, PHE_M = min(PHE), max(PHE) participant_results.append([cdepth_value, tdepth_value] + [AAE, VAR, STD, m, M, APHE, PHE_VAR, PHE_STD, PHE_m, PHE_M]) print len(participant_results), 'combinations processed...' np.save(os.path.join(root_result_path, '%s_2d3d_all.npy' % participant_label), np.array(participant_results)) np.savetxt(os.path.join(root_result_path, '%s_2d3d_all.csv' % participant_label), np.array(participant_results), delimiter=",") print '> Computing results for multiple calibration depths...' for num_of_calibration_depths in xrange(2, 6): # from 2 calibration depths to 5 participant_results = [] print '> Computing results for combining %s calibration depths...' % num_of_calibration_depths for calibs in combinations(depths, num_of_calibration_depths): # Now calibs is a set of depths, from each of those we need calibration data cp, ct3d = [], [] calib_depths_label = [] for calib in calibs: if len(cp): cp = np.concatenate((cp, participant_experiment[calib[0]][0]), axis=0) ct3d = np.concatenate((ct3d, participant_experiment[calib[0]][2]), axis=0) else: cp = participant_experiment[calib[0]][0] ct3d = participant_experiment[calib[0]][2] calib_depths_label.append(getDepth(calib)) # Performing calibration w, e, w0 = minimizeEnergy(cp, ct3d) e = v(e) # Now we have calibration data from multiple depths, we can test on all depths for test_depth in depths: tdepth_value = getDepth(test_depth) tp, tt3d = participant_experiment[test_depth[1]][0], participant_experiment[test_depth[1]][2] t, p = tt3d, tp qi = map(_q, p) gis = map(lambda q: g(q, w), qi) t = map(lambda vec: v(vec), t) gis = map(lambda vec: v(vec), gis) gprimes = map(lambda tg: v(((tg[0].z - e.z)/tg[1].z)*tg[1] + e), zip(t, gis)) AE = list(np.degrees(np.arctan((v(p[0]).cross(p[1])/(v(p[0]).dot(p[1]))).mag)) for p in zip(gprimes, t)) N = len(t) AAE = sum(AE)/N VAR = sum((ae - AAE)**2 for ae in AE)/N STD = np.sqrt(VAR) m, M = min(AE), max(AE) # Computing physical distance error (in meters) PHE = list((u-v).mag for u,v in zip(t, gprimes)) N = len(t) APHE = sum(PHE)/N PHE_VAR = sum((phe - APHE)**2 for phe in PHE)/N PHE_STD = np.sqrt(VAR) PHE_m, PHE_M = min(PHE), max(PHE) participant_results.append(calib_depths_label + [tdepth_value] + [AAE, VAR, STD, m, M, APHE, PHE_VAR, PHE_STD, PHE_m, PHE_M]) print len(participant_results), 'combinations processed...' result_path = os.path.join(root_result_path, '%s_calibration_depths/' % num_of_calibration_depths) if not os.path.exists(result_path): os.makedirs(result_path) np.save(result_path + '%s.npy' % participant_label, np.array(participant_results)) np.savetxt(result_path + '%s.csv' % participant_label, np.array(participant_results), delimiter=",") print 'done.' # class Parallax3Dto3DRealData(Experiment): # def __run__(self): # sim = GazeSimulation(log = False) # aae_3ds = [] # root_result_path = '/home/mmbrian/3D_Gaze_Tracking/work/results/3D3D/' # root_pose_path = '/home/mmbrian/3D_Gaze_Tracking/work/Marker_Eye_Images/ImagesUndist/' # root_data_path = '/home/mmbrian/HiWi/etra2016_mohsen/code/recording/data/participants/' # take_only_nearest_neighbor_for_calibration = True # participants = ['p14'] # if not os.path.exists(root_result_path): # os.makedirs(root_result_path) # for d1 in os.listdir(root_data_path): # if d1.startswith('p'): # every participant # if not d1 in participants: # # if not d1 in PARTICIPANTS: # continue # participant_label = d1 # participant_results = [] # d2 = os.path.join(root_data_path, d1) # .../pi/ # d2 = os.path.join(d2, os.listdir(d2)[0]) # .../pi/../ # print '> Processing participant', d1 # participant_experiment = {} # for d3 in os.listdir(d2): # every recording # d4 = os.path.join(d2, d3) # .../pi/../00X/ # # pose_info = np.loadtxt(open(os.path.join(root_pose_path+d1+'/'+d3, "null_pupils.csv"),"rb"),delimiter=";") # pose_info = np.loadtxt(open(os.path.join(root_pose_path+d1+'/'+d3, "simple_pupils.csv"),"rb"),delimiter=";") # frames_numbers = pose_info[:, 0] # pose_estimates = pose_info[:,4:7] # pose_info = dict(zip(frames_numbers, pose_estimates)) # p_frames = np.load(os.path.join(d4, 'p_frames.npy')) # # print d4 # # Fetching pose information for every target # poses = [] # for target in p_frames: # pose = [] # for fn in target: # all frames corresponding to this pupil # # first fn corresponds to the nearest neighbor # # for test use all correspondents from these 3 or 2 estimates # # i.e. each pose-marker creates a correspondence so 3*16=48 correspondents for test # # for calibration compare two cases, one similar to above take all the 75 correspondents # # and the other taking only the pose corresponding to nearest neighbor which results in # # the same number of correspondents as target markers # try: # pose.append(pose_info[fn]) # except KeyError, err: # print err # poses.append(pose) # t2d = np.load(os.path.join(d4, 't2d.npy')) # t3d = np.load(os.path.join(d4, 't3d.npy')) # participant_experiment[d3] = [poses, t2d, t3d] # keys = sorted(participant_experiment.keys()) # depths = zip(keys[::2], keys[1::2]) # for calib_depth in depths: # pose_data, ct3d = participant_experiment[calib_depth[0]][0], participant_experiment[calib_depth[0]][2] # cdepth_value = getDepth(calib_depth) # if take_only_nearest_neighbor_for_calibration: # pose = np.array(list(p[0] for p in pose_data)) # calib_3ds = ct3d[:] # else: # calib_3ds = [] # pose = [] # for i, p3d in enumerate(ct3d): # for p in pose_data[i]: # pose.append(p) # calib_3ds.append(p3d) # # Performing calibration # # First we convert gaze rays to actual pupil pose in our right hand coordinate system # # _pose = [(np.arctan(g.x/g.z), np.arctan(g.y/g.z)) for g in map(v, pose)] # _pose = map(v, pose) # print '> Running tests for calibration depth', cdepth_value # if any(g.z == 0 for g in _pose): # print 'Calibration is flawed' # # print pose # else: # print 'Calibration data is okay' # # print [g.mag for g in map(v, pose)] # # w, e, w0 = minimizeEnergy(_pose, calib_3ds, pose_given=True) # R, e = minimizeEnergy(pose, calib_3ds, pose_given=True) # # R = LA.inv(R) # print 'R', R # print 'e', e # e = v(e) # for test_depth in depths: # tdepth_value = getDepth(test_depth) # tpose_data, tt3d = participant_experiment[test_depth[1]][0], participant_experiment[test_depth[1]][2] # test_3ds = [] # tpose = [] # for i, p3d in enumerate(tt3d): # for p in tpose_data[i]: # tpose.append(p) # test_3ds.append(p3d) # # applying estimated rotation to bring pose vectors to scene camera coordinates # tpose = map(lambda p: v(R.dot(np.array(p))), tpose) # if any(g.z == 0 for g in map(v, tpose)): # print 'Test depth', tdepth_value, 'is flawed' # gis = map(lambda vec: v(vec), tpose) # t = map(lambda vec: v(vec), test_3ds) # gprimes = map(lambda tg: v(((tg[0].z - e.z)/tg[1].z)*tg[1] + e), zip(t, gis)) # # AE = list(np.degrees(np.arctan((v(p[0]).cross(p[1])/(v(p[0]).dot(p[1]))).mag)) for p in zip(gprimes, t)) # AE = list(np.degrees(np.arccos(v(p[0]).dot(p[1])/v(p[0]).mag/v(p[1]).mag)) for p in zip(gprimes, t)) # AAE = np.mean(AE) # STD = np.std(AE) # m, M = min(AE), max(AE) # # Computing physical distance error (in meters) # PHE = list((u-v).mag for u,v in zip(t, gprimes)) # APHE = np.mean(PHE) # PHE_STD = np.std(PHE) # PHE_m, PHE_M = min(PHE), max(PHE) # print 'Calibration', cdepth_value, 'Test', tdepth_value, AAE, 'degrees', APHE, 'meters' # participant_results.append([cdepth_value, tdepth_value] + [AAE, STD, m, M, APHE, PHE_STD, PHE_m, PHE_M]) # print len(participant_results), 'combinations processed...' # np.save(os.path.join(root_result_path, '%s_3d3d_all.npy' % participant_label), np.array(participant_results)) # np.savetxt(os.path.join(root_result_path, '%s_3d3d_all.csv' % participant_label), np.array(participant_results), delimiter=",") def main(): if len(sys.argv) <= 1: print 'Please select a mode.' return mode = sys.argv[1] if mode == 'pts': # Performs an experiment by fixing calibration depth and testing for different test depths # to investigate parallax error in 2D to 2D mapping ex = Parallax2Dto2DMapping() ex.performExperiment() if mode == '2d3d': # This also does 3D gaze estimation and plots estimation results for both 3D and 2D estimation ex = Parallax2Dto3DMapping() ex.performExperiment() if mode == '2d2d_2d3d': ex = Parallax3Dto3DMapping() ex.performExperiment() # ex = Parallax2Dto2DRealData() # ex.performExperiment() # ex = Parallax2Dto3DRealData() # ex.performExperiment() if mode == '3d3d_real': ex = Parallax3Dto3DRealData() ex.performExperiment() if __name__ == '__main__': main()