gazesim/code/sim.py

610 lines
24 KiB
Python

from __future__ import division
'''
This is the core part of the simulation framework.
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.
(notice that you need to disable plotting code so you can visualize and manipulate
the visualization otherwise after closing the plots the app stops)
'''
try:
import visual as vs
from visual import vector as vv
except ImportError:
from vector import Vector as vv
# from vector import Vector as vector
pass
from vector import Vector as v
import numpy as np
import cv2, cv
from geom import *
## 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
def log(self, msg):
if self._log: print msg
################################################################################
## Eye Specific Variables and Methods
################################################################################
# difference between center of cornea sphera and the plane containing limbus (based on anatomical data)
cornea_limbus_offset = 5.25
# difference between centers of spheres representing sclera and cornea (based on anatomical data)
sclera_cornea_center_offset = 4.7
# (based on anatomical data)
sclera_radius = 11.5
# (based on anatomical data)
cornea_radius = 7.8
# *unit for measurement is millimeter
# This direction should always be towards the target 3D point
# this is just the default value, adjust_eye method updates this vector
pupil_direction = v(1, 0, 0)
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
cornea_pos = sclera_pos + pupil_direction * sclera_cornea_center_offset
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
scene_eye_distance_factor = 1.0
def moveEyeToSceneCamera(self):
'''
distance_factor is between 0 and 1, controls where on the
difference vector between scene camera and eyeball to place
the eyeball (1 exactly on the target position, 0 does not change anything)
'''
if self.place_eyeball_on_scene_camera:
_diff = v(self.scene_camera.t) - self.sclera_pos
_diff = _diff * self.scene_eye_distance_factor
self.sclera_pos = self.sclera_pos + _diff
eye_camera_pos = v(self.eye_camera.t) + _diff
x, y, z = eye_camera_pos.x, eye_camera_pos.y, eye_camera_pos.z
self.eye_camera.setT((x, y, z))
self.recomputeEyeInner()
def setEyeRelativeToSceneCamera(self, xyz):
'''
places the eye in xyz relative to scene camera as the origin
(moves eye camera accordingly)
xyz: Vector object
'''
dest = v(self.scene_camera.t) + xyz # new dest originating from scene camera
_diff = dest - self.sclera_pos # difference from current pos to dest
self.sclera_pos = self.sclera_pos + _diff
eye_camera_pos = v(self.eye_camera.t) + _diff
x, y, z = eye_camera_pos.x, eye_camera_pos.y, eye_camera_pos.z
self.eye_camera.setT((x, y, z))
self.recomputeEyeInner()
def get_target_dir(self):
'''
Returns the direction of target from the eye
'''
return (self.target_position - self.sclera_pos).norm()
def adjust_eye(self):
'''
Adjusts eye in a way that it gazes at the target
'''
self.pupil_direction = self.get_target_dir()
cornea_pos = self.sclera_pos + self.pupil_direction * self.sclera_cornea_center_offset
# self.eye_inner.pos = cornea_pos
## Old approach which considers pupil center as the point positioned on the surface of
## the cornea sphere in direction of the gaze
# self.pupil_position = cornea_pos + self.pupil_direction * self.cornea_radius
## Now we consider pupil center on the same ray but closer to the cornea sphere's center
## their displacement is defined by average anatomical data and pupil is consider as the
## center of the limbus circle
self.pupil_position = cornea_pos + self.pupil_direction * self.cornea_limbus_offset
################################################################################
################################################################################
## Calibration and Test Data
num_calibration = 25 # calibration points
num_test = 16
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
calibration_grid = True
calibration_random_depth = False
test_grid = False
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
'''
self.calibration_points = generatePoints(25, self.min_xyz, self.max_xyz,
grid=self.calibration_grid,
randomZ=self.calibration_random_depth,
depth=self.calibration_depth, offset=ofr,
xoffset=self.scene_camera.t[0],
yoffset=self.scene_camera.t[1],
zoffset=self.scene_camera.t[2])
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
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
self.test_points = generatePoints(16, min_xyz, max_xyz,
grid=self.test_grid,
randomZ=self.test_random_depth,
randFixedZ=self.test_random_fixed_depth,
depth=self.test_depth, offset=ofr,
xoffset=self.scene_camera.t[0],
yoffset=self.scene_camera.t[1],
zoffset=self.scene_camera.t[2])
def setCalibrationDepth(self, depth):
'''
Sets calibration depth wrt to scene camera
'''
self.calibration_random_depth = False
if isinstance(depth, list):
# adding scene camera depth to each calibration depth
self.calibration_depth = map(lambda d: d+self.scene_camera.t[2], depth)
else:
self.calibration_depth = self.scene_camera.t[2] + depth
self.reset()
def setTestDepth(self, depth):
'''
Sets test depth wrt to scene camera
'''
self.test_random_depth = False
self.test_random_fixed_depth = False
self.test_depth = self.scene_camera.t[2] + depth
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))
self.eye_camera.setT((self.sclera_pos + v(0, 0,
self.sclera_cornea_center_offset + \
self.cornea_radius + \
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))
# Scene Camera
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
################################################################################
################################################################################
## Simulations Steps
def init(self):
self.setupCameras()
self.reset()
# self.moveEyeToSceneCamera() # Uncomment to place eyeball center on scene camera
def reset(self, ofr=0.5):
self.pupil_positions = []
self.c = 0
self.calibration_points = []
self.test_points = []
self.generateCalibrationPts(ofr = ofr)
self.generateTestPts(ofr = ofr)
# these might change if points are generated on a grid
self.num_calibration = len(self.calibration_points)
self.num_test = len(self.test_points)
if self.place_eyeball_on_scene_camera:
self.moveEyeToSceneCamera()
def gazeAtNextPoint(self, calibration = True):
'''
Gazes at the next Calibration/Test target point and records both pupil and
target positions.
'''
if calibration:
if self.c >= self.num_calibration: # finished processing calibration points
self.log('Processed all calibration points...')
return False
else:
self.target_position = v(self.calibration_points[self.c])
self.active_gaze_point = self.target_position
else:
if self.c >= self.num_test: # finished processing calibration points
self.log('Processed all test points...')
return False
else:
self.target_position = v(self.test_points[self.c])
self.active_gaze_point = self.target_position
self.c+=1
self.adjust_eye()
self.pupil_positions.append(self.pupil_position)
return True
def processData(self, calibration = True):
proceed = self.gazeAtNextPoint(calibration)
while proceed:
proceed = self.gazeAtNextPoint(calibration)
tr_target_projections_3d = None
tr_pupil_projections_3d = None
def computeGazeMap(self):
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 = []
hw = self.scene_camera.image_width/2.
for i in xrange(len(targets)):
p = targets[i]
if abs(p[0])>hw or abs(p[1])>hw:
invalid.append(i)
for i in invalid[::-1]:
targets = targets[:i] + targets[i+1:]
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)))
# 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)
# NOTE: calibration points are in world CCS
self.tr_3d_pupil_locations = self.pupil_positions[:]
self.tr_pupil_locations = pupil_locations[:] # for future use
self.tr_target_projections = targets[:]
# 3D position of target point projections
base = v(self.scene_camera.t) + v(self.scene_camera.direction) * self.scene_camera.f
pts = map(lambda p: np.array((p[0], p[1], 0)), targets) # adding extra z = 0
self.tr_target_projections_3d = map(lambda p:p+np.array([base.x, base.y, base.z]), pts)
# 3D position of pupil point projections
base = v(self.eye_camera.t) + v(self.eye_camera.direction) * self.eye_camera.f
pts = map(lambda p: np.array((p[0], p[1], 0)), pupil_locations) # adding extra z = 0
# TODO: apply back transformation to world CS
# temporary solution with our setting: x <> -x
pts = map(lambda p: np.array((-p[0], p[1], 0)), pupil_locations)
self.tr_pupil_projections_3d = map(lambda p:p+np.array([base.x, base.y, base.z]), pts)
# Normalizing points
self.log('Normalizing points...')
targets = self.scene_camera.getNormalizedPts(targets)
pupil_locations = self.eye_camera.getNormalizedPts(pupil_locations)
# Computing polynomial map
self.log('Computing polynomial map...')
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))]))
self.log('Successfully computed map...')
# Converting these to numpy arrays to support nice operands
targets = np.array(targets)
pupil_locations = np.array(pupil_locations)
self.log('Scene Camera Image Width: %s' % self.scene_camera.image_width)
self.c=0
self.pupil_positions = []
te_target_projections_3d = None
te_pupil_projections_3d = None
def projectAndMap(self):
targets_3d = map(lambda vec: v(vec), self.test_points)
# Processing Projections
targets = map(lambda p: v(p), self.test_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 = []
hw = self.scene_camera.image_width/2.
for i in xrange(len(targets)):
p = targets[i]
if abs(p[0])>hw or abs(p[1])>hw:
invalid.append(i)
for i in invalid[::-1]:
targets = targets[:i] + targets[i+1:]
pupil_locations = pupil_locations[:i] + pupil_locations[i+1:]
self.test_points = self.test_points[:i] + self.test_points[i+1:]
self.log('Removed %s invalid test points...' % (old-len(targets)))
self.te_3d_pupil_locations = self.pupil_positions[:]
self.te_pupil_locations = pupil_locations[:] # for future use
self.te_target_projections = targets[:]
# 3D position of target point projections
base = v(self.scene_camera.t) + v(self.scene_camera.direction) * self.scene_camera.f
pts = map(lambda p: np.array((p[0], p[1], 0)), targets) # adding extra z = 0
self.te_target_projections_3d = map(lambda p:p+np.array([base.x, base.y, base.z]), pts)
# 3D position of pupil point projections
base = v(self.eye_camera.t) + v(self.eye_camera.direction) * self.eye_camera.f
pts = map(lambda p: np.array((p[0], p[1], 0)), pupil_locations) # adding extra z = 0
# TODO: apply back transformation to world CS
# temporary solution: x <> -x
pts = map(lambda p: np.array((-p[0], p[1], 0)), pupil_locations)
self.te_pupil_projections_3d = map(lambda p:p+np.array([base.x, base.y, base.z]), pts)
if not len(targets): # all points are invalid, terminate simulation!
return -1
# Normalizing points
self.log('Normalizing points...')
targets = np.array(self.scene_camera.getNormalizedPts(targets))
pupil_locations = np.array(self.eye_camera.getNormalizedPts(pupil_locations))
# Compute mapped points
self.log('Computing %s mapped points...' % len(pupil_locations))
mapped_targets = np.array(map(self.mp, np.array(pupil_locations)))
self.log('Computed %s mapped points...' % len(mapped_targets))
C = v(self.scene_camera.t)
# Denormalize points and add camera plane's z to each
base = C + v(self.scene_camera.direction) * self.scene_camera.f
targets = self.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 = self.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))
N = len(targets_3d)
AE = [getAngularDiff(v(targets[i]), v(mapped_targets[i]), C) for i in xrange(N)]
AAE = np.mean(AE)
AE_STD = np.std(AE)
self.log('AE = %s, AAE = %s' % (AE, AAE))
targets_3d = map(lambda vec: vec-v(self.scene_camera.t), targets_3d) # originating from scene camera
mapped_targets = map(lambda vec: vec-v(self.scene_camera.t), mapped_targets)
mapped_targets = map(lambda pt: pt/pt[2], mapped_targets) # making z = 1
mapped_3d = map(lambda mt: v(v(mt[0])*mt[1].z), zip(mapped_targets, targets_3d))
# Computing physical distance error (in meters)
PHE = list((u-v).mag/1000. for u,v in zip(targets_3d, mapped_3d))
N = len(targets_3d)
APHE = np.mean(PHE)
PHE_STD = np.std(PHE)
PHE_m, PHE_M = min(PHE), max(PHE)
return AAE, PHE, AE_STD, PHE_STD
################################################################################
################################################################################
def __init__(self, log=True):
self._log=log
self.init()
def runCalibration(self):
self.reset()
self.processData(calibration = True)
self.computeGazeMap()
def runTest(self):
self.processData(calibration = False)
return self.projectAndMap()
def runSimulation(self):
'''
Runs a single simulation with default parameters and logging enabled
'''
self.runCalibration()
return self.runTest()
def computeGazeMapFromRealData(self, pupil_locations, targets, scene_camera_dimensions):
'''
pupil_locations are gaze positions obtained from pupil tracker. they are assumed to be normalized
targets are 2D projections of markers on scene camera. they are not assumed to be normalized.
(would be normalized here)
scene_camera_dimensions should be width and height of scene_camera images
'''
# Normalizing points
self.log('Normalizing target points...')
px, py = np.array(scene_camera_dimensions)/2.
w, h = np.array(scene_camera_dimensions)
targets = map(lambda p:np.array([(p[0] - px + w/2)/w,
(p[1] - py + h/2)/h]), targets)
# Computing polynomial map
self.log('Computing polynomial map...')
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))]))
self.log('Successfully computed map...')
def projectAndMapRealData(self, pupil_locations, targets, scene_camera_dimensions, targets_3d, camera_matrix, dist_coeffs):
ret = []
# Normalizing points
# self.log('Normalizing points...')
# px, py = np.array(scene_camera_dimensions)/2.
# w, h = np.array(scene_camera_dimensions)
# targets = map(lambda p:np.array([(p[0] - px + w/2)/w,
# (p[1] - py + h/2)/h]), targets)
# Compute mapped points
self.log('Computing %s mapped points...' % len(pupil_locations))
mapped_targets = np.array(map(self.mp, np.array(pupil_locations)))
self.log('Computed %s mapped points...' % len(mapped_targets))
self.log('Denormalizing estimated points...')
px, py = np.array(scene_camera_dimensions)/2.
w, h = np.array(scene_camera_dimensions)
offset = np.array([w/2-px, h/2-py])
mapped_targets = map(lambda p:(p*np.array([w, h]))-offset, mapped_targets)
N = len(targets)
PE = list(np.linalg.norm(b-a) for a,b in zip(targets, mapped_targets))
APE = sum(PE)/N
VAR = sum((pe - APE)**2 for pe in PE)/N
STD = np.sqrt(VAR)
m, M = min(PE), max(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=dist_coeffs)
mapped_3d = map(lambda p: [p[0], p[1], 1], mapped_3d[0])
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)
AAE = sum(AE)/N
VAR = sum((ae - AAE)**2 for ae in AE)/N
STD = np.sqrt(VAR)
m, M = min(AE), max(AE)
ret.extend([AAE, VAR, STD, m, M])
# multiplying corresponding ray with depth of ground truth target to get intersection
# with its depth plane
targets_3d = map(lambda vec: v(vec), targets_3d)
mapped_3d = map(lambda mt: v(v(mt[0])*mt[1].z), zip(mapped_3d, targets_3d))
# Computing physical distance error (in meters)
PHE = list((u-v).mag for u,v in zip(targets_3d, mapped_3d))
N = len(targets_3d)
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)
ret.extend([APHE, PHE_VAR, PHE_STD, PHE_m, PHE_M])
return ret
def perform2D2DCalibrationOnReadData(self, cp, ct, scene_camera_dimensions):
'''
cp, ct are calibration pupil points (normalized) and target points respectively
'''
self.computeGazeMapFromRealData(cp, ct, scene_camera_dimensions)
def run2D2DTestOnRealData(self, tp, tt, scene_camera_dimensions, targets_3d, camera_matrix, dist_coeffs):
'''
tp, tt are test pupil points (normalized) and target points respectively
'''
return self.projectAndMapRealData(tp, tt, scene_camera_dimensions, targets_3d, camera_matrix, dist_coeffs)
def runCustomSimulation(self, args):
self.calibration_grid, self.calibration_random_depth, self.test_grid, self.test_random_depth, self.test_random_fixed_depth = args
return self.runSimulation()
display_pupil_position = True
display_calibration_points = True
display_test_points = False
display_eye_camera = True
display_scene_camera = True
calibration = False #
display_test_point_rays = False #
display_active_gaze_point = False
display_calibration_point_rays = True
display_test_point_rays = True
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
# converts a vector object into a visual vector object
v_to_vv = lambda vec: vv(vec.x, vec.y, vec.z)
## World Axes
# drawAxes(None, vs.color.white, 13, (0, 0, 0))
## Eyeball component (outer sphere)
eye_outer = vs.sphere(pos=v_to_vv(self.sclera_pos),
radius=self.sclera_radius,
color=vs.color.red)
## Eyeball component (inner sphere)
eye_inner = vs.sphere(pos=v_to_vv(self.cornea_pos),
radius = self.cornea_radius,
color=vs.color.white)
self.recomputeEyeInner()
eye_inner.pos = v_to_vv(self.cornea_pos)
## Pupil position
if self.display_pupil_position:
vs.points(pos=[v_to_vv(self.pupil_position)], size=5, color=vs.color.black)
## calibration points
if self.display_calibration_points:
vs.points(pos=self.calibration_points, size=10, color=vs.color.yellow)
## Test points
if self.display_test_points:
vs.points(pos=self.test_points, size=10, color=vs.color.blue)
## Eye camera
if self.display_eye_camera:
drawCameraFrame(self.eye_camera)
## Scene camera
if self.display_scene_camera:
drawCameraFrame(self.scene_camera)
# Display rays from scene camera towards calibration points
if self.display_calibration_point_rays:
# Cast rays from scene camera towards calibration points
for point in self.calibration_points:
diff = vv(point) - vv(self.scene_camera.t)
drawLine(None, vv(self.scene_camera.t), diff.mag, diff.norm())
# Display rays from scene camera towards test points
if self.display_test_point_rays:
# Cast rays from scene camera towards calibration points
for point in self.test_points:
diff = vv(point) - vv(self.scene_camera.t)
drawLine(None, vv(self.scene_camera.t), diff.mag, diff.norm())
# active gaze point
if self.display_active_gaze_point and self.active_gaze_point:
vs.points(pos=[v_to_vv(self.active_gaze_point)], size=10, color=vs.color.red)
if self.tr_target_projections_3d:
vs.points(pos=self.tr_target_projections_3d, size=7, color=vs.color.red)
if self.tr_pupil_projections_3d:
vs.points(pos=self.tr_pupil_projections_3d, size=7, color=vs.color.red)
if self.te_target_projections_3d:
vs.points(pos=self.te_target_projections_3d, size=7, color=vs.color.blue)
if self.rays:
for ray in self.rays:
drawLine(None, tuple(ray.position), ray.length, tuple(ray.direction), ray.color)