From 67123bb970fb49442a63354f6c8925e0ec4f34df Mon Sep 17 00:00:00 2001
From: mohsen-mansouryar <mohsen@mpi-inf.mpg.de>
Date: Thu, 28 Apr 2016 18:14:38 +0200
Subject: [PATCH] updated code with more comments

---
 code/experiment.py             | 441 ---------------------------------
 code/geom.py                   |  61 +++--
 code/minimize.py               |   5 +
 code/parallax_2D3D_3Cdepths.py |   4 +-
 code/parallax_analysis.py      | 353 ++++++++++----------------
 code/recording_experiment.py   |   7 +-
 code/short_sim.py              | 358 --------------------------
 code/sim.py                    | 106 ++++----
 code/svis.py                   |   9 +-
 code/util.py                   |  28 +--
 code/vector.py                 |   6 +
 11 files changed, 262 insertions(+), 1116 deletions(-)
 delete mode 100644 code/experiment.py
 delete mode 100644 code/short_sim.py

diff --git a/code/experiment.py b/code/experiment.py
deleted file mode 100644
index a314a03..0000000
--- a/code/experiment.py
+++ /dev/null
@@ -1,441 +0,0 @@
-from __future__ import division
-import numpy as np
-from numpy import linalg as LA
-
-from time import time
-
-from sklearn import linear_model
-
-from sim import GazeSimulation
-from svis import Ray
-from minimize import findInitialW, _q, g, minimizeEnergy, minimizeEnergyGivenE, minimizeUsingScaleVector
-
-from geom import getSphericalCoords, getAngularDiff
-
-import visual as vs
-# from visual import vector as v # for vector operations
-from vector import Vector as v # for vector operations
-
-class Experiment:
-
-	def performExperiment(self):
-		print 'Running Experiment...'
-		start = time()
-		self.__run__()
-		self.runningTime = time() - start
-		print 'Running Time:', self.runningTime
-		if self.minimizationTime:
-			print 'Minimization Time:', self.minimizationTime
-
-	def __run__(self):
-		pass
-
-
-class Experiment1(Experiment):
-
-	def __run__(self):
-		'''
-		Simply evaluates 2D gaze estimation (using pupil calibration method)
-		'''
-		sim = GazeSimulation(log = False)
-		## Uncomment for shifting eyeball to scene camera coordinates
-		sim.scene_eye_distance_factor = 0.6 # controls how much to move eyeball towards scene camera [0-1]
-		sim.place_eyeball_on_scene_camera = False
-
-		# calibration_grid, calibration_random_depth,	test_grid, test_random_depth, test_random_fixed_depth
-		experiments = [
-			[True, 				False,				False,  		True,				False],
-			[True, 				False,				False,  		False,				False],
-			[True, 				False,				False,  		False,				True],
-			[False, 			True,				False,  		True,				False],
-			[False, 			True,				False,  		False,				False],
-			[False, 			True,				False,  		False,				True],
-			[False, 			False,				False,  		True,				False],
-			[False, 			False,				False,  		False,				False],
-			[False, 			False,				False,  		False,				True],
-			# [True, 				False,				True,  			False,				False],
-		]
-
-		nsim = 100
-		for experiment in experiments:
-			res = [sim.runCustomSimulation(experiment) for n in xrange(nsim)]
-			res = filter(lambda aae: aae>=0, res)
-			M = max(res)
-			m = min(res)
-			aae = sum(res) / len(res)
-			print 'AAE: %2.5f - Min AAE: %2.5f | Max AAE: %2.5f' % (aae, m, M)
-
-		eye_scene_diff = v(sim.sclera_pos) - v(sim.scene_camera.t)
-		print "Eyeball-SceneCamera distance was", eye_scene_diff.mag
-
-		self.sim = sim
-		self.minimizationTime = -1
-
-
-class Experiment1p5(Experiment):
-
-	def __run__(self):
-		'''
-		Read 3D gaze estimation experiment
-		'''
-		sim = GazeSimulation(log = False)
-		## Uncomment for shifting eyeball to scene camera coordinates
-		sim.place_eyeball_on_scene_camera = True
-		sim.scene_eye_distance_factor = 0.6
-		# sim.scene_eye_distance_factor = 1.0
-
-		# 	calibration_grid, calibration_random_depth,	test_grid, test_random_depth, test_random_fixed_depth
-		# -> grid setting
-		# con = [True, 				False,				True,  		False,				True]	
-		# -> general setting with grid learning
-		# con = [True, 				False,				False,  		True,				False]
-		con = [True, 				False,				False,  		False,				False]
-		# -> general setting with random learning
-		# con = [False, 				True,				False,  		True,				False]	
-		sim.num_calibration = 36
-		sim.num_test = 25
-		aae_2d = sim.runCustomSimulation(con)
-
-		## Fetching the calibration points
-		t, p = sim.calibration_points, sim.tr_pupil_locations
-		# target positions are computed relative to the scene CCS
-		t = map(lambda target: v(target) - v(sim.scene_camera.t), t) # verified
-		# p = sim.eye_camera.getNormalizedPts(p)
-
-		start = time()
-
-		eye_scene_diff = v(sim.sclera_pos) - v(sim.scene_camera.t)
-		e = np.array(eye_scene_diff) # eyeball coords in scene CCS
-		print 'Real e', e
-		e_org = e[:]
-		print 'Eye-Scene distance', eye_scene_diff.mag
-		# w, w0 = minimizeEnergyGivenE(p, t, e)
-
-		## Finding the optimal transformation matrix by minimizing the nonlinear energy
-		# w0 is the initial w by solving the leastsq with e=(0,0,0)
-		# w is by solving the leastsq again optimizing for both e and w
-		w, e, w0 = minimizeEnergy(p, t)
-		
-		# here were are initializing minimization with an e close to the ground
-		# truth e by adding random perturbations to it
-		# w, e, w0 = minimizeEnergy(p, t, e + np.random.rand(1, 3)[0])
-
-		print 'Estimated e', e, '( Distance:', LA.norm(e - e_org) , ')'
-		self.minimizationTime = time() - start
-		print
-
-		sim.w = w
-
-		## Fetching test points
-		t, p = sim.test_points, sim.te_pupil_locations
-		t = map(lambda target: v(target) - v(sim.scene_camera.t), t) # target coords in scene CCS
-
-		# closest point distance to scene camera
-		cDist = min(v(pt).mag for pt in t)
-		# farthest point distance to scene camera
-		fDist = max(v(pt).mag for pt in t)
-		# average point distance to scene camera
-		avgDist = sum(v(pt).mag for pt in t)/len(t) 
-
-
-		# p = sim.eye_camera.getNormalizedPts(p)
-		# print
-		# print p
-		# print
-
-
-		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) 
-		gis0 = map(lambda q: g(q, w0), qi) 
-
-		# now we can compare unit gaze vectors with their corresponding gaze rays t
-		# normalizing gaze rays first
-		t = map(lambda vec: v(vec).norm(), t)
-		# TODO: compare spherical coordinates instead
-
-
-		AE = list(np.degrees(np.arctan((v(p[0]).cross(p[1])/(v(p[0]).dot(p[1]))).mag)) for p in zip(gis, t))
-		N = len(t)
-		AAE = sum(AE)/N
-		VAR = sum((ae - AAE)**2 for ae in AE)/N
-		print 'AAE:', AAE, '\nVariance:', VAR, 'STD:', np.sqrt(VAR), '\nMin:', min(AE), 'Max:', max(AE), '(N=' + str(N) + ')'
-		print 'Target Distances: m=%s M=%s Avg=%s' % (cDist, fDist, avgDist)
-
-		AE0 = list(np.degrees(np.arctan((v(p[0]).cross(p[1])/(v(p[0]).dot(p[1]))).mag)) for p in zip(gis0, t))
-		AAE0 = sum(AE0)/N
-		print 'AAE (only optimizing W for e=(0,0,0)):', AAE0
-
-		print 'AAE (2D-to-2D mapping):', aae_2d
-		print ('Improvement (2D-2D vs 2D-3D): %s' % round((aae_2d - AAE)*100/aae_2d, 2)) + '%'
-		print ('Improvement (2D-2D vs 2D-2D): %s' % round((aae_2d - AAE0)*100/aae_2d, 2)) + '%'
-
-
-		sim.display_test_points = True
-		sim.display_calibration_points = False
-
-		## Visualizing gaze vector
-		ray_scale_factor = 120
-		ground_truth_gaze_rays = []
-		for pt in t:
-			ground_truth_gaze_rays.append(Ray(sim.scene_camera.t, ray_scale_factor, v(pt), vs.color.green))
-		sim.rays = ground_truth_gaze_rays
-
-		estimated_gaze_rays = []
-		for pt in gis:
-			estimated_gaze_rays.append(Ray(sim.scene_camera.t, ray_scale_factor, v(pt), vs.color.red))
-		sim.rays.extend(estimated_gaze_rays)
-
-		AE = []
-		for pair in zip(t, gis):
-			gt_r, e_r = pair
-			base = v(sim.scene_camera.t)
-			gt_ray_endpoint = base + (v(gt_r) * ray_scale_factor)
-			e_ray_endpoint = base + (v(e_r) * ray_scale_factor)
-			AE.append(getAngularDiff(gt_ray_endpoint, e_ray_endpoint, base))
-
-			diff = gt_ray_endpoint - e_ray_endpoint
-			sim.rays.append(Ray(e_ray_endpoint, diff.mag, diff.norm(), vs.color.orange))
-		
-		# Recomputing AAE (using law of cosines) sum(AE)/N
-		assert(round(sum(AE)/N - AAE) == 0)
-
-		self.sim = sim
-		print
-
-def p3dEx():
-	ex = Experiment1p5()
-	ex.performExperiment()
-	return ex
-
-def experiment2(s = 1):
-	sim = GazeSimulation(log = False)
-	sim.place_eyeball_on_scene_camera = True
-
-	# 2/3 for calibration, 1/3 for test
-	# sim.num_calibration = 66
-	# sim.num_test = 33
-	# sim.reset()
-
-	# 	calibration_grid, calibration_random_depth,	test_grid, test_random_depth, test_random_fixed_depth
-	# con = [True, 				False,				False,  		False,				True]
-	# con = [False, 				True,				False,  		True,				False]
-	con = [True, 				False,				True,  		False,				True]	
-	# sim.runSimulation()
-	# sim.num_calibration = 36
-	sim.num_calibration = 25
-	sim.num_test = 25
-	sim.runCustomSimulation(con)
-
-	# retrieving calibration points (t for 3D target positions and p for 2D pupil locations)
-	t, p = sim.calibration_points, sim.tr_pupil_locations
-	# t = sim.tr_target_projections # projections of target points
-	# target positions are computed relative to the scene CCS
-	t = map(lambda target: v(target) - v(sim.scene_camera.t), t) # verified
-	
-	# p = sim.eye_camera.getNormalizedPts(p)
-	
-	# computing the initial mapping matrix w by solving the energy minimization 
-	# for e=(0, 0, 0)
-	# w = findInitialW(p, t, True)
-	# print w
-	# w, e = minimizeEnergy(p, t)
-	# print w
-	w = np.array([[s, 0],
-				  [0, s]])
-	# w = minimizeUsingScaleVector(p, t)
-
-	# e = np.array(v(sim.sclera_pos) - v(sim.scene_camera.t)) # eyeball coords in scene CCS
-	# w = minimizeEnergyGivenE(p, t, e)
-
-	t, p = sim.test_points, sim.te_pupil_locations
-	t = map(lambda target: v(target) - v(sim.scene_camera.t), t) # target coords in scene CCS
-	
-	# TODO, is this necessary? (due to eye camera rotation)
-	p = map(lambda pt: np.array([-pt[0], pt[1]]), p)
-
-	# p = sim.eye_camera.getNormalizedPts(p)
-
-	# this w can be used for 2D-2D mapping (apparently!) testing it
-	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) 
-
-	# now we can compare unit gaze vectors with their corresponding gaze rays t
-	# normalizing gaze rays
-	t = map(lambda vec: v(vec).norm(), t)
-	SAE = sum(np.degrees(np.arctan((v(p[0]).cross(p[1])/(v(p[0]).dot(p[1]))).mag)) for p in zip(gis, t))
-	# return SAE / len(t)
-
-	sim.display_test_points = True
-	sim.display_calibration_points = False
-	## Visualizing gaze vector
-	ray_scale_factor = 120
-	ground_truth_gaze_rays = []
-	for pt in t:
-		ground_truth_gaze_rays.append(Ray(sim.scene_camera.t, ray_scale_factor, v(pt), vs.color.green))
-	sim.rays = ground_truth_gaze_rays
-
-	estimated_gaze_rays = []
-	for pt in gis:
-		estimated_gaze_rays.append(Ray(sim.scene_camera.t, ray_scale_factor, v(pt), vs.color.red))
-	sim.rays.extend(estimated_gaze_rays)
-
-	AE = []
-	for pair in zip(t, gis):
-		gt_r, e_r = pair
-		base = v(sim.scene_camera.t)
-		gt_ray_endpoint = base + (v(gt_r) * ray_scale_factor)
-		e_ray_endpoint = base + (v(e_r) * ray_scale_factor)
-		AE.append(getAngularDiff(gt_ray_endpoint, e_ray_endpoint, base))
-
-		diff = gt_ray_endpoint - e_ray_endpoint
-		sim.rays.append(Ray(gt_ray_endpoint, diff.mag, -diff.norm(), vs.color.orange))
-
-	return sim
-	# print e
-
-def experiment3():
-	sim = GazeSimulation(log = False)
-	# 2/3 for calibration, 1/3 for test
-	# sim.num_calibration = 66
-	# sim.num_test = 33
-	# sim.reset()
-
-	# 	calibration_grid, calibration_random_depth,	test_grid, test_random_depth, test_random_fixed_depth
-	con = [True, 				False,				False,  		True,				False]
-	# con = [False, 				True,				False,  		True,				False]
-	# sim.runSimulation()
-	sim.runCustomSimulation(con)
-	
-	# retrieving calibration points (t for 3D target positions and p for 2D pupil locations)
-	t, p = sim.calibration_points, sim.tr_pupil_locations
-	# t = sim.tr_target_projections # projections of target points
-	# target positions are computed relative to the scene CCS
-	t = map(lambda target: v(target) - v(sim.scene_camera.t), t)
-	
-	# p = sim.eye_camera.getNormalizedPts(p)
-	
-	# Applying regression using LARS (Least Angle Regression)
-	qi = map(_q, p)
-	clf = linear_model.Lars(n_nonzero_coefs=np.inf) # n_nonzero_coefs=1
-	t = map(lambda vec: v(vec).norm(), t)
-	clf.fit(qi, t)
-
-	t, p = sim.test_points, sim.te_pupil_locations
-	t = map(lambda target: v(target) - v(sim.scene_camera.t), t)
-	# p = sim.eye_camera.getNormalizedPts(p)
-
-	# this w can be used for 2D-2D mapping (apparently!) testing it
-	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: clf.predict(q)[0], qi) 
-
-	# now we can compare unit gaze vectors with their corresponding gaze rays t
-	# normalizing gaze rays
-	t = map(lambda vec: v(vec).norm(), t)
-	SAE = sum(np.degrees(np.arctan((v(p[0]).cross(p[1])/(v(p[0]).dot(p[1]))).mag)) for p in zip(gis, t))
-	print SAE / len(t)
-	# print e
-
-def experimentGT():
-	## Here we are using the ground truth gaze vectors (from the eye) and we compare
-	## them with ground truth vectors initiating from the scene camera
-	## In case eyeball and scene camera are located at the same point, the AAE should
-	## be reported as zero, this is only to make sure we're doing the right thing
-
-	sim = GazeSimulation(log = False)
-	# Uncomment to shift eyeball to scene camera
-	# sim.place_eyeball_on_scene_camera = True
-
-	# 2/3 for calibration, 1/3 for test
-	# sim.num_calibration = 66
-	# sim.num_test = 33
-	# sim.reset()
-
-	# 	calibration_grid, calibration_random_depth,	test_grid, test_random_depth, test_random_fixed_depth
-	con = [True, 				False,				False,  		True,				False]
-	# con = [False, 				True,				False,  		True,				False]
-	# sim.runSimulation()
-	
-	sim.runCustomSimulation(con)
-	
-	# retrieving calibration points (t for 3D target positions and p for 2D pupil locations)
-	t, p, tp = sim.calibration_points, sim.tr_pupil_locations, sim.tr_target_projections
-	# t = sim.tr_target_projections # projections of target points
-	# target positions are computed relative to the scene CCS
-	
-	t_org = t[:]
-	t = map(lambda target: v(target) - v(sim.scene_camera.t), t)
-	
-	rays_from_scene = []
-	for pt in t:
-		rays_from_scene.append(Ray(sim.scene_camera.t, v(pt).mag, v(pt).norm(), vs.color.red))
-	sim.rays = rays_from_scene
-
-	
-	rays_from_scene_to_target_projections = []
-	base = v(sim.scene_camera.t) + v(sim.scene_camera.direction) * sim.scene_camera.f
-	for pp in tp:
-		# print pp
-		# since we wanna draw rays starting from scene camera,
-		# ray is relative to scene CCS
-		ray = base + v(pp[0], pp[1], 0) - v(sim.scene_camera.t) 
-		rays_from_scene_to_target_projections.append(Ray(v(sim.scene_camera.t), ray.mag * 4, ray.norm(), vs.color.orange))
-	sim.rays.extend(rays_from_scene_to_target_projections)
-
-	t = map(lambda vec: v(vec).norm(), t)
-	
-	alpha_GT = []
-	_deg = lambda triple: (np.degrees(triple[0]), np.degrees(triple[1]), triple[2])
-	print 
-	for pt in t:
-		alpha_GT.append(_deg(getSphericalCoords(pt[0], pt[1], pt[2])))
-		# print _deg(getSphericalCoords(pt[0], pt[1], pt[2])), pt[2]
-	print 
-
-	# starting from eyeball
-	base = v(sim.sclera_pos)
-	# direction vector from eyeball towards eye camera
-	d = -v(sim.eye_camera.direction).norm()
-	# calculating distance between eyeball and image plane of the eye camera
-	dist = (v(sim.eye_camera.t) - v(sim.sclera_pos)).mag - sim.eye_camera.f
-
-	p = map(lambda p: np.array((-p[0], p[1])), p) # -x due to orientation
-	sf = 4
-	p = map(lambda p: p.dot(np.array([[sf, 0], 
-									  [0 , sf]])), p)
-	gis = map(lambda p: (d*dist + v(p[0], p[1], 0)).norm(), p)
-
-	alpha_test = []
-
-	rays_from_eye = []
-	for pt in gis:
-		alpha_test.append(_deg(getSphericalCoords(pt[0], pt[1], pt[2])))
-		rays_from_eye.append(Ray(sim.sclera_pos, 200, pt, vs.color.yellow))
-	sim.rays.extend(rays_from_eye)
-
-	# gis = map(lambda t: (v(t) - v(sim.sclera_pos)).norm(), t_org)
-	
-	# now we can compare unit gaze vectors with their corresponding gaze rays t
-	# normalizing gaze rays
-	
-	# for _t, _gis in zip(t, gis): print _t, _gis
-	SAE = sum(np.degrees(np.arctan((v(p[0]).cross(p[1])/(v(p[0]).dot(p[1]))).mag)) for p in zip(gis, t))
-	print SAE / len(t)
-	return sim
-
-
-
-# experiment1()
-
-# for s in np.arange(1, 20, 0.25):
-# 	print experiment2(s), s
-
-# experiment3()
-
-# experimentGT()
-# experiment1p5()
\ No newline at end of file
diff --git a/code/geom.py b/code/geom.py
index debff31..ea47e04 100644
--- a/code/geom.py
+++ b/code/geom.py
@@ -2,9 +2,13 @@ from __future__ import division
 import numpy as np
 from random import sample
 
-def_fov = np.pi * 2./3
-
 def generatePoints(n, min_xyz, max_xyz, grid=False, randomZ=True, randFixedZ=False, depth=None, offset=0.5, xoffset=0, yoffset=0, zoffset=0):
+  '''
+  Generates 3D points either on a grid or randomly. the depth of each point can also be selected
+  at random, given fixed depth for all points, or a random fixed depth for the points.
+
+  UPDATE: the order of the points in the end is also randomly shuffled.
+  '''
   if randFixedZ: # means all points have the same z but z is chosen at random between max and min z
     z = min_xyz[2] + np.random.random() * (max_xyz[2] - min_xyz[2])
   else: # same depth
@@ -65,7 +69,7 @@ def getSphericalCoords(x, y, z):
       /
      Z   
 
-    with a CounterClockwise rotation of the axis vectors
+    with a CCW rotation of the axis vectors
     '''
     r = np.sqrt(x*x + y*y + z*z)
     teta = np.arctan(x/z)
@@ -96,6 +100,10 @@ def getRotationMatrixFromAngles(r):
     [R_11 R_12 R_13]
     [R_11 R_12 R_13]
     [0    0    0   ] 
+
+    This is convenient since we can insert a translation column
+    and get a full transformation matrix for vectors in homogeneous
+    coordinates
     '''
     cos = map(np.cos, r)
     sin = map(np.sin, r)
@@ -155,20 +163,22 @@ def getRotationMatrix(a, b):
         return np.concatenate((np.eye(3) + vx, [[0, 0, 0]]))
     return np.concatenate((np.eye(3) + vx + vx.dot(vx)*((1-c)/s/s), [[0, 0, 0]]))
 
-
-
+# Default FOV angle for a camera
+def_fov = np.pi * 2./3
 
 class PinholeCamera:
     '''
-    Models a basic Pinhole Camera with 9 degrees of freedom
+    Models a Pinhole Camera with 9 degrees of freedom
     '''
-    # Intrinsic parameters
+    ################################################################################################
+    ## Intrinsic parameters
     f = 1 # focal length
     p = (0, 0) # position of principal point in the image plane
+    
     # Extrinsic parameters
     # this rotation corresponds to a camera setting pointing towards (0, 0, -1)
     r = (0, 0, 0) # rotations in x, y', and z'' planes respectively 
-    t = (0, 0, 0) # camera center translation w.r.t world coordinate system (with no rotation)
+    t = (0, 0, 0) # camera center translation w.r.t world coordinate system
     #
     # Using the above parameters we can construct the camera matrix
     #
@@ -178,15 +188,15 @@ class PinholeCamera:
     #
     # and thus we have x = PX for every point X in the word coordinate system 
     # and its corresponding projection in the camera image plane x
-    # NOTE: points are assumed to be represented by homogeneous vectors
-    #
+    # NOTE: points are assumed to be represented in homogeneous coordinates
+
     # Other parameters
     label = ''
-
-    direction = (0, 0, 1) # camera direction
+    direction = (0, 0, 1) # camera direction vector (viewpoint direction)
     fov = def_fov # field of view (both horizontal and vertical)
     image_width = 2*f*np.tan(fov/2.)
-    ################################################
+    ################################################################################################
+    
     def __init__(self, label, f = 1, r = (0, 0, 0), t = (0, 0, 0), direction = (0, 0, 1), fov=def_fov):
         self.label = label
         self.f = f
@@ -197,8 +207,11 @@ class PinholeCamera:
         self.recomputeCameraMatrix(True, True)
 
     def recomputeCameraMatrix(self, changedRotation, changedIntrinsic):
+        '''
+        Updates camera matrix
+        '''
         if changedRotation:
-            # # Computing rotation matrix using elemental rotations
+            # Computing rotation matrix using elemental rotations
             self.R = getRotationMatrixFromAngles(self.r)
             # by default if rotation is 0 then camera optical axis points to positive Z
             self.direction = np.array([0, 0, 1, 0]).dot(self.R) 
@@ -208,8 +221,8 @@ class PinholeCamera:
         self.Rt = np.concatenate((self.R, np.array([[_t[0], _t[1], _t[2], 1]]).T), axis=1)
         # instead of the above, we could also represent translation matrix as [I|-C] and
         #                              [R -RC]
-        # then compute Rt as R[I|-C] = [0   1]                   [R]      [-RC]
-        # but we're basically do the same thing by concatenating [0] with [  1]
+        # then compute Rt as R[I|-C] = [0   1]                [R]      [-RC]
+        # but we basically do the same thing by concatenating [0] with [  1]
 
         if changedIntrinsic:
             # Computing intrinsic matrix
@@ -218,10 +231,11 @@ class PinholeCamera:
                 [f, 0, px, 0],
                 [0, f, py, 0],
                 [0, 0, 1,  0]])
+
         # Full Camera Projection Matrix
         self.P = self.K.dot(self.Rt)
         
-    ################################################
+    ################################################################################################
     ## Intrinsic parameter setters
     def setF(self, f, auto_adjust = True):
         self.f = f
@@ -235,7 +249,8 @@ class PinholeCamera:
     def setFOV(self, fov):
         self.fov = fov
         self.image_width = 2*self.f*np.tan(fov/2.)
-    ################################################
+    
+    ################################################################################################
     ## Extrinsic parameter setters
     def setT(self, t, auto_adjust = True):
         self.t = t
@@ -245,11 +260,13 @@ class PinholeCamera:
       self.r = r
       if auto_adjust:
           self.recomputeCameraMatrix(True, False)
-    ################################################
+    
+    ################################################################################################
+    ## Main methods
     def project(self, p):
         '''
         Computes projection of a point p in world coordinate system
-        using its homogeneous vector coordinates [x, y, z, 1]
+        using its homogeneous coordinates [x, y, z, 1]
         '''
         if len(p) < 4: p = (p[0], p[1], p[2], 1)
         projection = self.P.dot(np.array(p))
@@ -272,6 +289,10 @@ class PinholeCamera:
         return map(lambda p:(p*self.image_width)-offset, pts)
 
 class Camera(PinholeCamera):
+    '''
+    A wrapper around abstract Pinhole camera object with further
+    methods for visualization.
+    '''
     default_radius = 2.0
 
     def updateShape(self):
diff --git a/code/minimize.py b/code/minimize.py
index a9d1608..625d175 100644
--- a/code/minimize.py
+++ b/code/minimize.py
@@ -6,6 +6,11 @@ from sklearn.preprocessing import PolynomialFeatures as P
 from vector import Vector as v
 import cv2
 
+'''
+This module contains all the methods used in solving the minimization problems
+in our paper.
+'''
+
 M = 7 # this is the size of the feature vector
 _q = lambda p: np.array([p[0], p[1], p[0]*p[0], p[1]*p[1], p[0]*p[1], p[0]*p[1]*p[0]*p[1], 1])
 
diff --git a/code/parallax_2D3D_3Cdepths.py b/code/parallax_2D3D_3Cdepths.py
index 47fe325..e376617 100644
--- a/code/parallax_2D3D_3Cdepths.py
+++ b/code/parallax_2D3D_3Cdepths.py
@@ -164,6 +164,8 @@ class Parallax2Dto3DMapping(Experiment):
 
 				results.append([calibs, aae_2ds_aae, aae_3ds_aae, aae_3D3Ds])
 
+		######################################################################################################
+		## Plotting part
 		plt.ylabel('Angular Error')
 		plt.xlabel('Depth')
 
@@ -223,7 +225,6 @@ class Parallax2Dto3DMapping(Experiment):
 
 		plt.yticks(TOPICS, LABELS,fontsize=18)
 
-		# plt.legend(fontsize=20)
 		plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3,
            ncol=3, mode="expand", borderaxespad=0., fontsize=15)
 		plt.yticks(fontsize='18')
@@ -242,6 +243,7 @@ class Parallax2Dto3DMapping(Experiment):
 
 		plt.subplots_adjust(left=left, bottom=bottom, right=right, top=top, wspace=wspace, hspace=hspace)
 		plt.show()
+		######################################################################################################
 
 if __name__ == '__main__':
 	# This also does 3D gaze estimation and plots estimation results for both 3D and 2D estimation
diff --git a/code/parallax_analysis.py b/code/parallax_analysis.py
index eb956cb..2812aa3 100644
--- a/code/parallax_analysis.py
+++ b/code/parallax_analysis.py
@@ -60,17 +60,13 @@ class Parallax2Dto2DMapping(Experiment):
 		print 'scene_camera', sim.scene_camera.t
 		print 'calibration'
 		print len(sim.calibration_points)
-		# print 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 'test'
 		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])
-
-		
+		print 'depths', set(np.array(sim.test_points)[:,2])	
 		
 		plt.ylabel('Y (mm)')
 		plt.xlabel('X (mm)')
@@ -84,18 +80,17 @@ class Parallax2Dto2DMapping(Experiment):
 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 differentl7y
-	this transformation is essential. I would add the corresponding computations later on.
+	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))
-		# sim.setEyeRelativeToSceneCamera(v(-65, -33, 0)) # assuming eyeball and scene camera are coplanar i.e. e = (e.x, e.y, 0)
+		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)
@@ -213,49 +208,8 @@ class Parallax2Dto3DMapping(Experiment):
 			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])])
 		
-		# Old plot code
 		######################################################################################################
-		# plt.ylabel('Angular Error')
-		# plt.xlabel('Depth')
-
-		# fig = plt.figure(figsize=(14.0, 10.0))
-		# ax = fig.add_subplot(111)
-
-		# clrs = ['b', 'r', 'orange']
-
-
-		# _xrange = [0.5,1.5,2.5,3.5,4.5]
-		# ax.plot(_xrange, [res[0] for res in results], 'r', label='2D-to-2D', marker="o", linestyle='-',lw=3)
-		# ax.plot(_xrange, [res[1] for res in results], 'b', label='2D-to-3D', marker="o", linestyle='-',lw=3)
-		# ax.plot(_xrange, [res[2] for res in results], 'g', label='3D-to-3D', marker="o", linestyle='-',lw=3)
-
-		# ax.set_ylabel(r'Angular Error',fontsize=22, fontweight='bold')
-		# ax.set_xlabel(r'Number of Calibration Depths',fontsize=22, fontweight='bold')
-		
-		# plt.legend(fontsize=20)
-		# # plt.legend(loc="upper left", ncol=3, title=r"$d_c$")
-		# # plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3,
-		# # 				ncol=5, mode="expand", borderaxespad=0., fontsize=20)
-		# # plt.xticks(fontsize='18')
-		# plt.yticks(fontsize='18')
-
-
-		# TOPICs = [0.0,0.5,1.5,2.5,3.5,4.5,5.0]
-		# LABELs = ['', '1', '2', '3', '4', '5', '']
-		# # ax.set_xticklabels(LABELs,fontsize=18)
-		# 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()
-		######################################################################################################
-		# New plot code based on EffectNumberofClusters.py
+		# 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]
@@ -263,36 +217,17 @@ class Parallax2Dto3DMapping(Experiment):
 		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
-
-		#    x1 = [0.4,1.4,2.4,3.4,4.4]
+		width = 0.5 # the width of the bars
 		x2 = [0.45,1.45,2.45,3.45,4.45]
-		#    x3 = [0.5,1.5,2.5,3.5,4.5]
 		x4 = [0.55,1.55,2.55,3.55,4.55]
-		#    x5 = [0.6,1.6,2.6,3.6,4.6]
 		x6 = [0.50,1.50,2.50,3.50,4.50]
 
 		fig = plt.figure(figsize=(14.0, 10.0))
-
 		ax = fig.add_subplot(111)
 
-		#    print mean2D2D
-		#    print mean2D3D
-
-		#    ax.axhline(linewidth=2, y = np.mean(mean2D2D),color='r')
-		#    ax.axhline(linewidth=2, y = np.mean(mean2D3D),color='blue')
-
-		#    ax.axhline(linewidth=2, y = minvaluevalue,color='black')
-
-		#    ax.text(0.98, Participantmeanvalue+0.5, "Mean %.2f" % Participantmeanvalue,fontsize=12, fontweight='bold',color='r')
-		#    ax.text(0.98, maxvaluevalue+0.5, "Maximum %.2f" % maxvaluevalue,fontsize=12, fontweight='bold',color='black')
-		#    ax.text(0.98, minvaluevalue+0.5, "Minimum %.2f" % minvaluevalue,fontsize=12, fontweight='bold', color='black')
-
-		#    rects1 = ax.bar(ind, Participantmean,width, color='r',edgecolor='black',)#, hatch='//')
 		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')
 
@@ -304,35 +239,20 @@ class Parallax2Dto3DMapping(Experiment):
 
 		legend(fontsize=20,loc='upper right')
 
-		#    rects3 = ax.errorbar(x3, meanC3,yerr=[stdC3,stdC3],fmt='o',color='black',ecolor='black',lw=3, capsize=5, capthick=2)
-		#    plt.plot(x3, meanC3, marker="o", linestyle='-',lw=3,color='black')
-		#    
-		#    rects4 =ax.errorbar(x4, meanC4,yerr=[stdC4,stdC4],fmt='o',color='green',ecolor='green',lw=3, capsize=5, capthick=2)
-		#    plt.plot(x4, meanC4, marker="o", linestyle='-',lw=3,color='green')
-		#    
-		#    rects5 =ax.errorbar(x5, meanC5,yerr=[stdC5,stdC5],fmt='o',color='orange',ecolor='orange',lw=3, capsize=5, capthick=2)
-		#    plt.plot(x5, meanC5, marker="o", linestyle='-',lw=3,color='orange')
-
-
 		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]#,110]#,120]
+		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', ""]#, ""]#, ""]
-
-		#    fig.canvas.set_window_title('Distance Error Correlation')
+		LABELs = ["",r'1',r'2', r'3', r'4', r'5', ""]
 		plt.xticks(TOPICs, LABELs,fontsize=18)
 
-		#    legend([rects1,rects2], [r'\LARGE\textbf{2D2D}', r'\LARGE\textbf{2D3D}'], loc='lower right') 
-
-		TOPICS = [0.5,1,1.5,2,2.5,3,3.5,4,4.5,5]#,110]#,120]
+		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']#, ""]#, ""]
+		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']
 
-		#    fig.canvas.set_window_title('Accuracy - Activity Statistics')
 		plt.yticks(TOPICS, LABELS,fontsize=18)
 
 		def autolabel(rects):
@@ -342,8 +262,6 @@ class Parallax2Dto3DMapping(Experiment):
 				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)
 
-			#    autolabel(rects1)
-
 
 		left  = 0.1  # the left side of the subplots of the figure
 		right = 0.975    # the right side of the subplots of the figure
@@ -351,19 +269,15 @@ class Parallax2Dto3DMapping(Experiment):
 		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()
 		######################################################################################################
 
 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.setEyeRelativeToSceneCamera(v(-65, -33, 0)) # assuming eyeball and scene camera are coplanar i.e. e = (e.x, e.y, 0)
 		sim.setCalibrationDepth(1 * 1000) # mm, wrt scene camera
 		sim.setTestDepth(1.5 * 1000)
 
@@ -465,8 +379,7 @@ class Parallax3Dto3DMapping(Experiment): # GT pupil pose instead of estimating t
 				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))
-				# 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(t)
 				AAE = np.mean(AE)
 				STD_2D3D = np.std(AE)
@@ -478,15 +391,14 @@ class Parallax3Dto3DMapping(Experiment): # GT pupil pose instead of estimating t
 				PHE_STD = np.std(PHE)
 				PHE_m, PHE_M = min(PHE), max(PHE)
 
-				# aae_3Ds.append((AAE, STD, PHE, PHE_STD))
 				aae_3Ds.append(AAE)
 				std_3Ds.append(STD_2D3D)
-				# break
 
 			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 = []
@@ -532,7 +444,6 @@ class Parallax3Dto3DMapping(Experiment): # GT pupil pose instead of estimating t
 
 		ax.set_ylabel(r'\textbf{Angular Error}',fontsize=22)
 		ax.set_xlabel(r'\textbf{Depth}',fontsize=22)
-		# ax.set_ylim((0, 2.4))
 
 		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']#, ""]#, ""]
@@ -558,9 +469,10 @@ class Parallax3Dto3DMapping(Experiment): # GT pupil pose instead of estimating t
 
 		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):
 	'''
@@ -662,8 +574,6 @@ class Parallax2Dto2DRealData(Experiment):
 
 
 		print 'done.'
-		# plt.plot(tdrange, aaes)
-		# plt.show()
 
 def getDepth(depth_experiments):
 	_map = {'000': 1, '002': 1.25, '004': 1.5, '006': 1.75, '008': 2.0}
@@ -673,6 +583,7 @@ 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)
@@ -823,136 +734,136 @@ class Parallax2Dto3DRealData(Experiment):
 
 
 
-class Parallax3Dto3DRealData(Experiment):
-	def __run__(self):
-		sim = GazeSimulation(log = False)
-		aae_3ds = []
+# 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
+# 		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']
+# 		participants = ['p14']
 
-		if not os.path.exists(root_result_path):
-			os.makedirs(root_result_path)
+# 		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/
+# 		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)
+# 					# 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]
+# 					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])
+# 				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
+# 				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]
+# 					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)
+# 						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)
+# 						# 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'
+# 						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))
+# 						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)
+# 						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)
+# 						# 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 '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=",")
+# 				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():
@@ -971,8 +882,6 @@ def main():
 		ex = Parallax2Dto3DMapping()
 		ex.performExperiment()
 
-	# ex = Parallax2Dto3DMappingEye()
-	# ex.performExperiment()i
 	if mode == '2d2d_2d3d':
 		ex = Parallax3Dto3DMapping()
 		ex.performExperiment()
diff --git a/code/recording_experiment.py b/code/recording_experiment.py
index 31de8fc..2807e10 100644
--- a/code/recording_experiment.py
+++ b/code/recording_experiment.py
@@ -8,12 +8,15 @@ from time import time
 from geom import getSphericalCoords, getAngularDiff
 from recording.tracker import Marker
 
-# from visual import vector as v # for vector operations
-from vector import Vector as v
+try:
+	from visual import vector as v
+except ImportError:
+	from vector import Vector as v
 
 DATA_DIR = './recording/data/'
 # DATA_DIR = '.\\recording\\data\\'
 
+# Update these accordingly for testing out your own data
 EYE_CAMERA_IMAGE_WIDTH = 640
 EYE_CAMERA_IMAGE_HEIGHT = 360
 
diff --git a/code/short_sim.py b/code/short_sim.py
deleted file mode 100644
index bd934da..0000000
--- a/code/short_sim.py
+++ /dev/null
@@ -1,358 +0,0 @@
-from __future__ import division
-
-import visual as vs
-import wx 
-
-import numpy as np
-import pylab as plt
-
-import math 
-from random import random as rnd
-
-from util import *
-from geom import *
-
-from pupil.calibrate import get_map_from_cloud
-mp = None # mapping obtained from training set
-
-display_width = 800
-display_height = 600
-
-panel_width = 224
-win = vs.window(width=display_width + panel_width, height=display_height, title='Gaze Simulation')
-scene = vs.display(window=win, width=display_width, height=display_height)
-
-################################################################################
-## World Axes
-drawAxes(None, vs.color.white, 13, (0, 0, 0))
-################################################################################
-################################################################################
-## Eye Representation
-eye_camera_pos = (-15, 20, 0)
-scene_camera_pos = (0, 20, 0)
-
-# This is gonna be the offset distance between the two spheres which are modelling an eye
-eye_ball_offset = 4
-
-eye_outer_pos = vs.vector((-15, 20, 15))
-
-##_diff = vs.vector(scene_camera_pos) - eye_outer_pos
-##eye_outer_pos = eye_outer_pos + _diff
-##eye_camera_pos = vs.vector(eye_camera_pos) + _diff
-##eye_camera_pos = (eye_camera_pos.x, eye_camera_pos.y, eye_camera_pos.z)
-
-
-eye_outer_radius = 5.5
-eye_outer = vs.sphere(pos=eye_outer_pos,
-                      radius=eye_outer_radius, 
-                      color=vs.color.red)
-
-# This direction should always be towards the target 3D point
-# this is just the default value
-pupil_direction = vs.vector(1, 0, 0)
-pupil_position = None # computed whenever target is changed
-pupil_positions = [] # 3D positions of the pupil to be projected onto the eye camera
-
-eye_inner_radius = 3
-eye_inner_pos = eye_outer_pos + pupil_direction * eye_ball_offset
-eye_inner = vs.sphere(pos=eye_inner_pos,
-                      radius = eye_inner_radius,
-                      color=vs.color.white)
-
-
-target_position = vs.vector(0, 0, 0)
-def get_target_dir():
-  '''
-  Returns the direction of target from the eye
-  '''
-  return (target_position - eye_outer_pos).norm()
-
-def adjust_eye():
-  '''
-  Adjusts eye in a way that it gazes at the target
-  '''
-  global pupil_direction, pupil_position
-  pupil_direction = get_target_dir()
-  eye_inner.pos = eye_outer_pos + pupil_direction * eye_ball_offset
-  pupil_position = eye_inner.pos + pupil_direction * eye_inner_radius
-
-adjust_eye()
-
-def drawPupilPoint():
-  vs.points(pos=[pupil_position], size=10, color=vs.color.blue)
-################################################################################
-################################################################################
-## 3D Points for Training and Test
-num_training = 25
-num_test = 25
-
-min_xyz = -25, -5, -110
-max_xyz = 45, 65, -30
-
-training_points = generatePoints(num_training, min_xyz, max_xyz, True, False)
-test_points = generatePoints(num_test, min_xyz, max_xyz, False, True)
-
-
-def drawTrainingPoints():
-  vs.points(pos=training_points, size=10, color=vs.color.yellow)
-def drawTestPoints():
-  vs.points(pos=test_points, size=10, color=vs.color.blue)
-
-drawTrainingPoints()
-# drawTestPoints()
-################################################################################
-
-################################################################################
-## Camera Object
-class Camera(PinholeCamera):
-  default_radius = 2.0
-
-  def drawCameraFrame(self):  # create frame and draw its contents
-    c = vs.vector(self.t)
-    self.center = vs.sphere(pos=c,
-                            radius=Camera.default_radius,
-                            color=vs.color.green)
-    self.dir = vs.arrow(pos=c,
-                        axis=vs.vector(self.direction) * self.f,
-                        shaftwidth=1.0)
-    self.img_plane = vs.box(pos=c + self.dir.axis,
-                            length=self.image_width, 
-                            width=0.5,
-                            height=self.image_width,
-                            color=vs.color.white,
-                            opacity=0.5)
-
-  def updateShape(self):
-    c = vs.vector(self.t)
-    # Update camera center
-    self.center.pos = c
-    # Update the arrow
-    self.dir.pos = c
-    self.dir.axis = vs.vector(self.direction) * self.f
-    # Update image plane
-    self.img_plane.pos = c + self.dir.axis
-    self.img_plane.length = self.image_width
-    self.img_plane.height = self.image_width
-
-# Eye Camera
-eye_camera = Camera(label='Eye Camera', f = 5)
-eye_camera.setR((0, 0, 0))
-eye_camera.setT(v(eye_camera_pos))# + v(30, 0, 0))
-# eye_camera.setDir((0, 0, 1))
-eye_camera.drawCameraFrame()
-# Scene Camera
-scene_camera = Camera(label='Scene Camera', f = 25, fov=math.pi/3)
-scene_camera.setT(v(scene_camera_pos) + v(10, 5, -5))
-scene_camera.setR((0, np.pi, 0))
-# scene_camera.setDir((0, 0, -1))
-scene_camera.drawCameraFrame()
-
-# Cast rays from scene camera towards training points
-for point in training_points:
-  diff = vs.vector(point) - vs.vector(scene_camera.t)
-  drawLine(None, vs.vector(scene_camera.t), diff.mag, diff.norm())
-
-################################################################################
-## Widgets
-pan = win.panel   # addr of wx window object 
-pan.SetSize((display_width + panel_width, display_height))
-
-pan_minx = display_width + 7
-wx.StaticText(pan, pos=(pan_minx, 7),
-    label = "Select an item..." )
-
-def hReset(evt):
-  pass
-
-c = 0
-def gazeAtNextPoint(evt):
-  global c, target_position, pupil_positions
-  drawTrainingPoints()
-  if c >= num_training:
-    print 'Processed all training points...'
-    return
-
-  print 'Shifted gaze to point %d = %s...' % (c, training_points[c])
-
-  vs.points(pos=[training_points[c]], size=10, color=vs.color.red)
-
-  target_position = vs.vector(training_points[c])
-  c+=1
-  adjust_eye()
-  # drawPupilPoint()
-  pupil_positions.append(pupil_position)
-
-def computeMap(evt):
-  global mp, c, pupil_positions
-  # Processing Projections
-  targets = map(lambda p:vs.vector(p), training_points)
-  targets = map(lambda p: scene_camera.project((p.x, p.y, p.z, 1)), targets)
-  targets = map(lambda p: (-p[0], p[1]), targets) # x values are in a different CS, hence there are inverted
-  # TODO: instead of inverting which only holds for this specific setting, rotate using scene_camera.R
-  
-  pupil_locations = map(lambda p: eye_camera.project((p.x, p.y, p.z, 1)), pupil_positions)
-  # eye_camera is using the same world CS, so no reversing is required
-
-  # filtering out points outside of the image plane
-  old = len(targets)
-  invalid = []
-  for i in xrange(len(targets)):
-    p = targets[i]
-    if abs(p[0])>scene_camera.image_width/2. or abs(p[1])>scene_camera.image_width/2.:
-      invalid.append(i)
-  for i in invalid[::-1]:
-    targets = targets[:i] + targets[i+1:]
-    pupil_locations = pupil_locations[:i] + pupil_locations[i+1:]
-  print 'Removed %s invalid target points...' % (old-len(targets))
-  # Displaying target point projections
-  base = scene_camera.center.pos + scene_camera.dir.axis
-  pts = map(lambda p: np.array((p[0], p[1], 0)), targets)
-  vs.points(pos=map(lambda p:p+np.array([base.x, base.y, base.z]), pts), size=7, color=vs.color.red)
-  print 'Eye Camera Image Width:', eye_camera.image_width
-  # Displaying pupil point projections
-  base = eye_camera.center.pos + eye_camera.dir.axis
-  # drawLine(None, eye_camera.t, eye_camera.f * 3, eye_camera.direction, vs.color.yellow)
-  pts = map(lambda p: np.array((p[0], p[1], 0)), pupil_locations)
-  vs.points(pos=map(lambda p:p+np.array([base.x, base.y, base.z]), pts), size=3, color=vs.color.red)
-  # normalizing points
-  print 'Normalizing points...'
-  targets = scene_camera.getNormalizedPts(targets)
-  pupil_locations = eye_camera.getNormalizedPts(pupil_locations)
-  # computing polynomial map
-  print 'Computing polynomial map...'
-  mp = get_map_from_cloud(np.array([(pupil_locations[i][0], pupil_locations[i][1],
-                            targets[i][0], targets[i][1]) for i in xrange(len(targets))]))
-  # from sklearn.neighbors import KNeighborsRegressor as knn
-  # ngh = knn(n_neighbors = 1)
-  # ngh.fit(pupil_locations, targets)
-  # mp = ngh.predict
-  print 'Successfully computed map...'
-
-  # Converting these to numpy arrays to support nice operands
-  targets = np.array(targets)
-  pupil_locations = np.array(pupil_locations)
-
-  print 'Scene Camera Image Width:', scene_camera.image_width
-  plt.subplot(211)
-  plt.title('Target Points')
-  plt.plot(targets[:, 0], targets[:, 1], 'ro')
-  plt.axis([0, 1, 0, 1])
-  plt.xlabel('X')
-  plt.ylabel('Y')
-
-  labels = ['p{0}'.format(i) for i in xrange(len(targets))]
-  for label, x, y in zip(labels, targets[:, 0], targets[:, 1]):
-    plt.annotate(
-        label, 
-        xy = (x, y), xytext = (-13, 13),
-        textcoords = 'offset points', ha = 'right', va = 'top',
-        bbox = dict(boxstyle = 'round,pad=0.25', fc = 'yellow', alpha = 0.25),
-        arrowprops = dict(arrowstyle = '->', connectionstyle = 'arc3,rad=0'))
-
-  plt.subplot(212)
-  plt.title('Pupil Points')
-  plt.plot(pupil_locations[:, 0], pupil_locations[:, 1], 'bo')
-  plt.axis([0, 1, 0, 1])
-  plt.xlabel('X')
-  plt.ylabel('Y')
-
-  labels = ['p{0}'.format(i) for i in xrange(len(pupil_locations))]
-  for label, x, y in zip(labels, pupil_locations[:, 0], pupil_locations[:, 1]):
-    plt.annotate(
-        label, 
-        xy = (x, y), xytext = (-13, 13),
-        textcoords = 'offset points', ha = 'right', va = 'top',
-        bbox = dict(boxstyle = 'round,pad=0.25', fc = 'yellow', alpha = 0.25),
-        arrowprops = dict(arrowstyle = '->', connectionstyle = 'arc3,rad=0'))
-
-  c=0
-  pupil_positions = []
-  plt.show()
-
-def showTestPts(evt):
-  global c, target_position, pupil_positions
-  drawTestPoints()
-  if c >= num_test:
-    print 'Processed all test points...'
-    return
-
-  vs.points(pos=[test_points[c]], size=10, color=vs.color.red)  
-
-  print 'Shifted gaze to point %d = %s...' % (c, test_points[c])
-  target_position = vs.vector(test_points[c])
-  adjust_eye()
-  pupil_positions.append(pupil_position)
-  c+=1
-
-def projectAndMap(evt):
-  # Processing Projections
-  targets = map(lambda p:vs.vector(p), test_points)
-  targets = map(lambda p: scene_camera.project((p.x, p.y, p.z, 1)), targets)
-  targets = map(lambda p: (-p[0], p[1]), targets) # due to orientation of SC
-  pupil_locations = map(lambda p: eye_camera.project((p.x, p.y, p.z, 1)), pupil_positions)
-  # filtering out points outside of the image plane
-  old = len(targets)
-  invalid = []
-  for i in xrange(len(targets)):
-    p = targets[i]
-    if abs(p[0])>scene_camera.image_width/2. or abs(p[1])>scene_camera.image_width/2.:
-      invalid.append(i)
-  for i in invalid[::-1]:
-    targets = targets[:i] + targets[i+1:]
-    pupil_locations = pupil_locations[:i] + pupil_locations[i+1:]
-  print 'Removed %s invalid target points...' % (old-len(targets))
-  # normalizing points
-  print 'Normalizing points...'
-  targets = np.array(scene_camera.getNormalizedPts(targets))
-  pupil_locations = np.array(eye_camera.getNormalizedPts(pupil_locations))
-  # Compute mapped points
-  print pupil_locations
-  print 'Computing %s mapped points...' % len(pupil_locations)
-  mapped_targets = np.array(map(mp, np.array(pupil_locations)))
-  print 'Computed %s mapped points...' % len(mapped_targets)
-  print mapped_targets
-  
-  C = scene_camera.center.pos
-  # ssd = sum((vs.vector(mapped_targets[i]) - vs.vector(targets[i])).mag**2 for i in xrange(len(mapped_targets)))
-  # Computing Angular Error (Sum of Angular Differences)
-
-  # denormalize points and add camera plane's z to each 
-  base = scene_camera.center.pos + scene_camera.dir.axis
-
-  targets = scene_camera.getDenormalizedPts(targets)
-  targets = np.array(map(lambda p: np.concatenate((p + np.array([base.x, base.y]), np.array([base.z])), axis=0), targets))
-
-  mapped_targets = scene_camera.getDenormalizedPts(mapped_targets)
-  mapped_targets = np.array(map(lambda p: np.concatenate((p + np.array([base.x, base.y]), np.array([base.z])), axis=0), mapped_targets))
-
-  AE = list(getAngularDiff(vs.vector(targets[i]), vs.vector(mapped_targets[i]), C) for i in xrange(len(mapped_targets)))
-  N = len(mapped_targets)
-  AAE = sum(AE)/N
-  VAR = sum((ae - AAE)**2 for ae in AE)/N
-  print 'AAE = %s, Variance = %s, STD = %s' % (AAE, VAR, np.sqrt(VAR))
-
-
-  vs.points(pos=mapped_targets, size=2, color=vs.color.yellow)
-
-  plt.subplot(111)
-  plt.plot(targets[:, 0], targets[:, 1], 'bo', label='Projections')
-  plt.plot(mapped_targets[:, 0], mapped_targets[:, 1], 'ro', label='Mapped Estimations')
-
-  for _mp, tp in zip(mapped_targets, targets):
-    coords = zip(*(_mp, tp))
-    plt.plot(coords[0], coords[1], 'k-')
-
-  plt.legend()
-  plt.show()
-
-y = 25
-Button('Reset', pan_minx, y, hReset, (panel_width - 28, 40), pan)
-y += 40 + 7
-Button('Gaze at Next Point', pan_minx, y, gazeAtNextPoint, (panel_width - 28, 40), pan)
-y += 40 + 7
-Button('Compute Map', pan_minx, y, computeMap, (panel_width - 28, 40), pan)
-y += 40 + 7
-Button('Generate Test Points', pan_minx, y, showTestPts, (panel_width - 28, 40), pan)
-y += 40 + 7
-Button('Project, Map, Plot!', pan_minx, y, projectAndMap, (panel_width - 28, 40), pan)
-################################################################################
diff --git a/code/sim.py b/code/sim.py
index 92adfcc..9ac7661 100644
--- a/code/sim.py
+++ b/code/sim.py
@@ -1,29 +1,36 @@
 from __future__ import division
+
 '''
 This is the core part of the simulation framework.
-I am basically decoupling logic from visualization, this module covers the logic.
+This module mainly covers the logic meaning that you can create 
+a GazeSimulation object, specify the environment and and then 
+simply run the experiment with no visualization to get the data.
+
+However, you can as well visualize the environment after the simulation
+environment is set up and you ran the experiment. in order to do so, 
+after you have a GazeSimulation object, simply call visualize() method on it.
 '''
 
-# import visual as vs
-# from visual import vector as v # for vector operations
-from vector import Vector as v # for vector operations
+try:
+	import visual as vs
+	from visual import vector as v
+except ImportError:
+	from vector import Vector as vector
 import numpy as np
 import cv2, cv
 
-# from sklearn.neighbors import KNeighborsRegressor as knn
-# ## (u, v, u^2, v^2, uv, u^2*v^2, 1)
-# _q = lambda p: (p[0], p[1], p[0]*p[0], p[1]*p[1], p[0]*p[1], p[0]*p[1]*p[0]*p[1], 1)
-# _qi = lambda q: (q[0], q[1])
-
-# from util import *
 from geom import *
-# from svis import Visualizable # Uncomment if you want visualization
+## Uncomment the following import if you want visualization
+# from svis import Visualizable 
 
 # 2D to 2D calibration method from the pupil project
 from pupil.calibrate import get_map_from_cloud
 
+# a factor to convert from centimeter to millimeter later when
+# generating data (unit for our measurements is millimeter)
 point_scale_factor = 10
 
+## Uncomment the following if you want visualization (and comment the line after)
 # class GazeSimulation(Visualizable):
 class GazeSimulation:
 	_log = True
@@ -43,9 +50,9 @@ class GazeSimulation:
 	# *unit for measurement is millimeter
 
 	# This direction should always be towards the target 3D point
-	# this is just the default value
+	# this is just the default value, adjust_eye method updates this vector
 	pupil_direction = v(1, 0, 0)
-	pupil_position = None # computed whenever target is changed
+	pupil_position = None # 3D point which is computed whenever target is changed
 	pupil_positions = [] # 3D positions of the pupil to be projected onto the eye camera's image plane
 
 	sclera_pos = v(-15, 20, -25) # wrt world coordinate system
@@ -53,6 +60,9 @@ class GazeSimulation:
 
 	target_position = v(0, 0, 0)
 	def recomputeEyeInner(self):
+		'''
+		This follows our eye model to compute the position of cornea at every time.
+		'''
 		self.cornea_pos = self.sclera_pos + self.pupil_direction * self.sclera_cornea_center_offset
 
 	place_eyeball_on_scene_camera = False
@@ -109,18 +119,15 @@ class GazeSimulation:
 		self.pupil_position = cornea_pos + self.pupil_direction * self.cornea_limbus_offset
 	################################################################################  
 	################################################################################  
-	## calibration and Test Data
+	## Calibration and Test Data
 	num_calibration = 25 # calibration points
 	num_test = 16
-	# TODO: consider the least number of calibration points required, also investigate if accuracy is
-	# significantly improved with more points (my guess is that it will not)
 
-	# TODO: considering measurement unit is mm, adjust these parameters accordingly to simulate more
-	# realistic settings
-	# TODO: also consider the range of eye movement (50x50 degrees) to compute realistic boundaries
-	# for position of gaze targets
 	min_xyz = -55 * point_scale_factor, -45 * point_scale_factor, 1000
 	max_xyz = 55 * point_scale_factor, 45 * point_scale_factor, 2000
+
+	# These are considered if we have 25 calibration points and 16 test points,
+	# Changing those parameters require to recompute these values
 	grid_x_width = 22 * point_scale_factor
 	grid_y_width = 18 * point_scale_factor
 	calibration_depth, test_depth = None, None
@@ -131,9 +138,11 @@ class GazeSimulation:
 	test_random_depth = False
 	test_random_fixed_depth = False
 	def generateCalibrationPts(self, ofr = 0):
-		# For now we're considering 25 calibration points (5x5 grid) and 16 test points (inner 4x4 grid)
-		# as per the real world experiment. these parameters could be given though, they're hard coded for
-		# simplicity
+		'''
+		For now we're considering 25 calibration points (5x5 grid) and 16 test points (inner 4x4 grid)
+		as per the real world experiment. these parameters could be given though, they're hard coded for
+		simplicity
+		'''
 		self.calibration_points = generatePoints(25, self.min_xyz, self.max_xyz, 
 												grid=self.calibration_grid, 
 												randomZ=self.calibration_random_depth,
@@ -141,11 +150,14 @@ class GazeSimulation:
 												xoffset=self.scene_camera.t[0],
 												yoffset=self.scene_camera.t[1],
 												zoffset=self.scene_camera.t[2])
-		# Test point for checking spherical coordinates
-		# self.calibration_points.append(np.array(v(self.scene_camera.t)+120*v(self.scene_camera.direction) + v(-10, 0, 0)))
-	def generateTestPts(self, ofr = 0): # TODO
+		
+	def generateTestPts(self, ofr = 0): 
+		'''
+		Same as above, the values used here only apply to our setup, simply replace them with your
+		own values or introduce new variables to compute them at runtime. (due to our tight time schedule
+		this part of the code uses hardcoded material)
+		'''
 		min_xyz, max_xyz = self.min_xyz, self.max_xyz
-		# if self.num_calibration == 25 and self.num_test == 16 and self.test_grid and self.calibration_grid:
 		min_xyz = -55 * point_scale_factor+self.grid_x_width/2., -45 * point_scale_factor+self.grid_y_width/2., 1000
 		max_xyz = 55 * point_scale_factor-self.grid_x_width/2., 45 * point_scale_factor-self.grid_y_width/2., 2000
 		# The above two lines are also currently hard coded only to match the inner grid with our setting
@@ -158,10 +170,10 @@ class GazeSimulation:
 											xoffset=self.scene_camera.t[0],
 											yoffset=self.scene_camera.t[1],
 											zoffset=self.scene_camera.t[2])
-		# self.test_points = map(lambda p: v(p)-v(self.scene_camera.t), self.test_points)
+
 	def setCalibrationDepth(self, depth):
 		'''
-		sets calibration depth wrt to scene camera
+		Sets calibration depth wrt to scene camera
 		'''
 		self.calibration_random_depth = False
 		if isinstance(depth, list): 
@@ -173,7 +185,7 @@ class GazeSimulation:
 
 	def setTestDepth(self, depth):
 		'''
-		sets test depth wrt to scene camera
+		Sets test depth wrt to scene camera
 		'''
 		self.test_random_depth = False
 		self.test_random_fixed_depth = False
@@ -181,29 +193,26 @@ class GazeSimulation:
 		self.reset()
 		
 
+	# temporary variable to store active gaze point during the experiment
 	active_gaze_point = None
 	################################################################################  
 	################################################################################  
 	## Camera Setup
 	def setupCameras(self):
 		# Eye Camera
-		self.eye_camera = Camera(label='Eye Camera', f = 13, fov=np.radians(50)) # TODO: f and fov must be realistic
-		# self.eye_camera.setT((-15, 20, 0))
+		self.eye_camera = Camera(label='Eye Camera', f = 13, fov=np.radians(50)) 
 		self.eye_camera.setT((self.sclera_pos + v(0, 0, 
 												self.sclera_cornea_center_offset + \
 												self.cornea_radius + \
-												45)).tuple()) # 35mm in front of the surface of the eye
+												45)).tuple()) # 45mm in front of the surface of the eye
 		# rotating around y axis, camera now pointing towards negative Z
 		# with this rotation, also x becomes inverted
 		self.eye_camera.setR((0, np.pi, 0)) 
-		# self.eye_camera.setDir((0, 0, -1)) # deprecated
 
 		# Scene Camera
-		self.scene_camera = Camera(label='Scene Camera', f = 50, fov=np.radians(100)) # TODO: f and fov must be realistic
-		# self.scene_camera.setT((0, 20, 0)) 
-		self.scene_camera.setT((self.sclera_pos + v(-25, 25, 25)).tuple()) # sclera_pos is v(-15, 20, -25)
+		self.scene_camera = Camera(label='Scene Camera', f = 50, fov=np.radians(100))
+		self.scene_camera.setT((self.sclera_pos + v(-25, 25, 25)).tuple()) # sclera_pos is by default v(-15, 20, -25)
 		self.scene_camera.setR((0, 0, 0)) # camera pointing towards positive Z
-		# self.scene_camera.setDir((0, 0, 1)) # deprecated
 	################################################################################  
 	################################################################################  
 	## Simulations Steps
@@ -229,7 +238,7 @@ class GazeSimulation:
 
 	def gazeAtNextPoint(self, calibration = True):
 		'''
-		Gazes at the next calibration/Test target point and records both pupil and 
+		Gazes at the next Calibration/Test target point and records both pupil and 
 		target positions.
 		'''
 		if calibration:
@@ -263,6 +272,7 @@ class GazeSimulation:
 		targets = map(lambda p:v(p), self.calibration_points)
 		targets = map(lambda p: self.scene_camera.project((p.x, p.y, p.z, 1)), targets)
 		pupil_locations = map(lambda p: self.eye_camera.project((p.x, p.y, p.z, 1)), self.pupil_positions)
+		
 		# Filtering out points outside of the image plane
 		old = len(targets)
 		invalid = []
@@ -276,8 +286,7 @@ class GazeSimulation:
 			pupil_locations = pupil_locations[:i] + pupil_locations[i+1:]
 			self.calibration_points = self.calibration_points[:i] + self.calibration_points[i+1:]
 		self.log('Removed %s invalid calibration points...' % (old-len(targets)))
-		# print 'Removed %s invalid target points...' % (old-len(targets))
-
+		
 		# to get calibration points p and t, refer to self.calibration_points and 
 		# self.pupil_locations (also target_projections has the corresponding projection
 		# points for the targets)
@@ -309,10 +318,7 @@ class GazeSimulation:
 		
 		self.mp = get_map_from_cloud(np.array([(pupil_locations[i][0], pupil_locations[i][1],
 		                targets[i][0], targets[i][1]) for i in xrange(len(targets))]))
-		# ngh = knn(n_neighbors = 1)
-		# ngh.fit(map(_q, pupil_locations), map(_q, targets))
-		# self.mp = lambda p: _qi(ngh.predict(_q(p))[0])
-
+		
 		self.log('Successfully computed map...')
 		# Converting these to numpy arrays to support nice operands
 		targets = np.array(targets)
@@ -401,7 +407,6 @@ class GazeSimulation:
 		APHE = np.mean(PHE)
 		PHE_STD = np.std(PHE)
 		PHE_m, PHE_M = min(PHE), max(PHE)
-		# ret.extend([APHE, PHE_VAR, PHE_STD, PHE_m, PHE_M])
 		return AAE, PHE, AE_STD, PHE_STD
 
 	################################################################################  
@@ -417,7 +422,6 @@ class GazeSimulation:
 
 	def runTest(self):
 		self.processData(calibration = False)
-		# AAE = self.projectAndMap()
 		return self.projectAndMap()
 
 	def runSimulation(self):
@@ -473,15 +477,11 @@ class GazeSimulation:
 		VAR = sum((pe - APE)**2 for pe in PE)/N
 		STD = np.sqrt(VAR)
 		m, M = min(PE), max(PE)
-
-		# ret.append(['pixel_error: APE, VAR, STD, min, MAX, data', APE, VAR, STD, m, M, PE])
 		ret.extend([APE, VAR, STD, m, M])
 
 		self.log('Computing 3D back projection of estimated points...')
-		# mapped_3d = cv2.undistortPoints(np.array([mapped_targets]), cameraMatrix=camera_matrix, distCoeffs=np.array([0, 0, 0, 0, 0]))
 		mapped_3d = cv2.undistortPoints(np.array([mapped_targets]), cameraMatrix=camera_matrix, distCoeffs=dist_coeffs)
 		mapped_3d = map(lambda p: [p[0], p[1], 1], mapped_3d[0])
-		# print v(mapped_3d[0]).mag, np.sqrt(mapped_3d[0][0]**2 + mapped_3d[0][1]**2)
 
 		AE = list(np.degrees(np.arctan((v(p[0]).cross(v(p[1]))/(v(p[0]).dot(v(p[1])))).mag)) for p in zip(mapped_3d, targets_3d))
 		N = len(targets_3d)
@@ -534,6 +534,10 @@ class GazeSimulation:
 	rays = []
 
 	def draw(self):
+		'''
+		This visualizes the simulation environment. be careful about the lines you have to
+		uncomment at the top of this module so that this part works.
+		'''
 		from util import drawAxes, drawLine
 		from svis import drawCameraFrame
 		## World Axes
diff --git a/code/svis.py b/code/svis.py
index c86d34b..0eeb854 100644
--- a/code/svis.py
+++ b/code/svis.py
@@ -1,9 +1,12 @@
 '''
 This is the visualization backend for bare_sim module
 '''
-import visual as vs
-# from visual import vector as v
-from vector import Vector as v
+try:
+	import visual as vs
+	from visual import vector as v
+except ImportError:
+	from vector import Vector as v
+
 from geom import Camera
 
 from util import drawAxes
diff --git a/code/util.py b/code/util.py
index e70d43e..848d09d 100644
--- a/code/util.py
+++ b/code/util.py
@@ -1,12 +1,12 @@
-import visual as vs
-# from visual import vector as v
-from vector import Vector as v
+try:
+  import visual as vs
+  from visual import vector as v
+except ImportError:
+  from vector import Vector as vector
+
 import wx
 
 def drawAxes(frame, color, size, position):
-    # directions = (vs.vector(size, 0, 0),
-    #               vs.vector(0, size, 0),
-    #               vs.vector(0, 0, size))
     directions = (vs.vector(size, 0, 0),
                   vs.vector(0, size, 0),
                   vs.vector(0, 0, size))
@@ -21,24 +21,16 @@ def drawAxes(frame, color, size, position):
                  text=labels[i],
                  pos=origin+directions[i])
 
-# def wCoords(x, y, z):
-#   return wCoords(v(x, y, z))
-# def wCoords(p):
-#   return v(-p.x, p.y, -p.z)  
-
 def drawLine(frame, position, length, direction, color = vs.color.white):
 	'''
 	Draws a straight line from given position, with given length and direction
 	(all are relative to given frame)
 	'''
-	return vs.curve(frame=frame, pos=[v(position), v(position)+v(direction).norm()*length], color = color)
+	return vs.curve(
+            frame=frame, 
+            pos=[v(position), v(position)+v(direction).norm()*length], 
+            color = color)
 	
-# def reDrawLine(line, position, length, direction):
-# 	line.pos = [position, position+direction.norm()*length]
-
-def Button(label, x, y, func, size, pan):
-   btn = wx.Button(pan, label=label, pos=(x, y), size=size)
-   btn.Bind(wx.EVT_BUTTON, func)	
 
 
 
diff --git a/code/vector.py b/code/vector.py
index a7cd8e3..293330e 100644
--- a/code/vector.py
+++ b/code/vector.py
@@ -1,6 +1,12 @@
 from math import sqrt
 from numpy import array
 
+'''
+Vector object has the same functionality as VPython's vector. 
+this module is used when only vector processing is required and
+VPython is not available (e.g. non-Windows platforms)
+'''
+
 class Vector:
 	x, y, z = (0.0, 0.0, 0.0)