First commit

CartToSpherical.m

@@ -0,0 +1,8 @@
+% CartToSpherical.m
+%
+% Converts cartesian coordinates to spherical coordinates.
+
+function [horRads, verRads] = CartToSpherical(vec)
+    horRads = atan2(vec(3), vec(1));
+    verRads = acos(vec(2));
+end

EquirectToSpherical.m

@@ -0,0 +1,8 @@
+% EquirectToSpherical.m
+%
+% This function converts equirectangular coordinates to spherical coordinates.
+
+function [horRads, verRads] = EquirectToSpherical(xEq, yEq, widthPx, heightPx)
+    horRads = (xEq *2 * pi) / widthPx;
+    verRads = (yEq * pi) / heightPx;
+end

GetCartVectors.m

@@ -0,0 +1,56 @@
+% GetCartVectors.m
+%
+% This function returns an nx3 arrays of cartesian vectors for the eye direction
+% within the FOV of the headset, the eye direction in the world coordinates and
+% the head direction in the world.
+%
+% input:
+%   data        - ARFF data
+%   metadata    - metadata of ARFF data
+%   attributes  - attributes of ARFF data
+%
+% output:
+%   eyeFovVec   - nx3 (x,y,z) array corresponding to FOV vector of gaze centered 
+%                 around the middle of the video corresponding to (1,0,0) vector
+%   eyeHeadVec  - nx3 (x,y,z) array corresponding to the direction of the eye in 
+%                 the world. That is the head+eye position
+%   headVec     - nx3 (x,y,z) array corresponding to the direction of the head in
+%                 the world
+
+function [eyeFovVec, eyeHeadVec, headVec] = GetCartVectors(data, metadata, attributes)
+    c_xName = 'x';
+    c_yName = 'y';
+    c_xHeadName = 'x_head';
+    c_yHeadName = 'y_head';
+    c_angleHeadName = 'angle_deg_head';
+
+    xInd = GetAttPositionArff(attributes, c_xName);
+    yInd = GetAttPositionArff(attributes, c_yName);
+    xHeadInd = GetAttPositionArff(attributes, c_xHeadName);
+    yHeadInd = GetAttPositionArff(attributes, c_yHeadName);
+    angleHeadInd = GetAttPositionArff(attributes, c_angleHeadName);
+
+    eyeFovVec = zeros(size(data,1),3);
+    eyeHeadVec = zeros(size(data,1),3);
+    headVec = zeros(size(data,1),3);
+
+    for ind=1:size(data,1)
+        [horHeadRads, verheadRads] = EquirectToSpherical(data(ind, xHeadInd), data(ind, yHeadInd), metadata.width_px, metadata.height_px);
+        curHeadVec = SphericalToCart(horHeadRads, verheadRads);
+        headVec(ind,:) = curHeadVec;
+
+        angleHeadRads = data(ind, angleHeadInd) * pi / 180;
+        videoVec = [-1; 0; 0]; % middle of the video is considered its center
+
+        rot = HeadToVideoRot(curHeadVec, angleHeadRads, videoVec);
+        rot = rot';
+
+        [horGazeRads, verGazeRads] = EquirectToSpherical(data(ind, xInd), data(ind, yInd), metadata.width_px, metadata.height_px);
+        gazeVec = SphericalToCart(horGazeRads, verGazeRads);
+        eyeHeadVec(ind,:) = gazeVec;
+
+        gazeWithinVec = RotatePoint(rot, gazeVec);
+
+        eyeFovVec(ind,:) = gazeWithinVec(:);
+    end
+end

GetDispersion.m

@@ -0,0 +1,43 @@
+% GetDispersion.m
+%
+% This function returns the minimum angle (half of the total) of a cone
+% starting from the beginning of the coordinate system, which is the place in
+% the middle of the 360 degree sphere, with the cone vector denoting the
+% direction of the middle of the cone that contains all the data points.
+%
+% input:
+%   coneVec     - vector denoting the direction of the cone
+%   vecList     - nx3 list of vectors
+%
+% output:
+%   dispersion  - cone angle (degrees) that contains all vectors
+%   index       - index to vector list that gave the maximum radius
+
+function [dispersion, index] = GetDispersion(coneVec, vecList)
+    assert(size(vecList,2) == 3, 'Vector list not in correct format');
+    if (size(coneVec,1) > size(coneVec,2))
+        coneVec = coneVec';
+    end
+    assert(size(coneVec,1) == 1 && size(coneVec,2) == 3, 'Provided vector has wrong dimensions');
+    c_maxDelta = 0.001;
+
+    dispersion = 0;
+    for ind=1:size(vecList,1)
+        dotProd = sum(coneVec .* vecList(ind,:));
+
+		assert(dotProd < 1+c_maxDelta, 'Provided vectors are not normalized');
+        % account for double roundings
+        if (dotProd > 1)
+            dotProd = 1;
+        end
+
+        rads = acos(dotProd);
+        if (rads >= dispersion)
+            dispersion = rads;
+            index = ind;
+        end
+    end
+
+    dispersion = dispersion * 180 / pi;
+end
+

GetMaxDispersion.m

@@ -0,0 +1,31 @@
+% GetMaxDispersion.m
+%
+% This function gets 2 vector lists as input and returns the maximum dispersion
+% for an input vector in comparison with the second vector list.
+% 
+% input:
+%   vecList1    - list of vectors to iterate through
+%   vecList2    - (optional) list of vectors to use for dispersion
+%
+% output:
+%   maxDisp     - maximum dispersion
+%   index1      - pointer to the element in the vecList1 with maxDisp
+%   index2      - pointer to the element in the vecList2 with maxDisp
+
+function [maxDisp, index1, index2] = GetMaxDispersion(vecList1, vecList2)
+    if (nargin < 2)
+        vecList2 = vecList1;
+    end
+    maxDisp = 0;
+    index = -1;
+
+    for ind=1:size(vecList1)
+        [dispersion, tmpInd] = GetDispersion(vecList1(ind,:), vecList2);
+        if (dispersion > maxDisp)
+            maxDisp = dispersion;
+            index1 = ind;
+            index2 = tmpInd;
+        end
+    end
+
+end

GetSpeed.m

@@ -0,0 +1,38 @@
+% GetSpeed.m
+%
+% This function gets a list of vectors together with their timestamps and
+% returns a speed vector in degrees per second.
+%
+% input:
+%   vecList - nx3 (x,y,z) list of vectors. Taken from GetCartVectors.m
+%   time    - timestamp of vectors in us (microseconds)
+%   step    - (optional, default=1) distance of samples for speed calculation
+%
+% output:
+%   speed   - speed in deg/sec
+
+function speed = GetSpeed(vecList, time, step)
+    if (nargin < 3)
+        step = 1;
+    end
+    c_maxDelta = 0.001;
+    c_usToSec = 1000000;
+    
+    assert(size(time,1) == size(vecList,1), 'Provided number of vectors and timestamps should be the same')
+    speed = zeros(size(vecList,1),1);
+    for ind=step+1:size(vecList,1)
+        dotProd = sum(vecList(ind,:) .* vecList(ind-step,:));
+
+        assert(dotProd < 1+c_maxDelta, 'Provided vectors are not normalized');
+        % account for double roundings
+        if (dotProd > 1)
+            dotProd = 1;
+        end
+
+        rads = acos(dotProd);
+        degs = rads * 180 / pi;
+        elapsedTime = (time(ind) - time(ind-step)) / c_usToSec;
+        speed(ind - round(step/2)) = degs/elapsedTime;
+        assert(speed(ind) >= 0 && ~isinf(speed(ind)), ['Provided timestamps are not monotonic ' num2str(ind)]);
+    end
+end

HeadToVideoRot.m

@@ -0,0 +1,10 @@
+% HeadToVideoRot.m
+%
+% Calculates the rotation of the head vector to the provided video vector.
+
+function [rot] = HeadToVideoRot(headVec, headAngleRads, videoVec)
+    headToRef = YZXrotation(headVec, -headAngleRads);
+    videoToRef = YZXrotation(videoVec, 0);
+
+    rot = videoToRef*(headToRef');
+end

OnePointDriftCorrection.m

@@ -0,0 +1,188 @@
+% OnePointDriftCorrection.m
+%
+% This function applies drift correction for 360 degree videos. For correcting
+% the drift it has to know when the point was displayed and at which video
+% position (equirectangular) was displayed. It then calculates the mean head and
+% gaze directions during this time. It basically follows the inverse of the drift
+% correction during recording and applies the new displacements.
+%
+% NOTE: Be careful to provide a time interval in which gaze is not noisy
+%
+% input:
+%   arffFile        - path to ARFF file
+%   targetStartTime - target diplay starting time in us
+%   targetEndTime   - target diplay starting time in us
+%   targetPos       - [x, y] equirectangular position
+%   saveFile        - file to save data
+
+function OnePointDriftCorrection(arffFile, targetStartTime, targetEndTime, targetPos, saveFile)
+    c_timeName = 'time';
+	c_xName = 'x';
+    c_yName = 'y';
+    c_confName = 'confidence';
+    c_xHeadName = 'x_head';
+    c_yHeadName = 'y_head';
+    c_angleHeadName = 'angle_deg_head';
+    c_confThreshold = 0.8;
+
+    [data, metadata, attributes, relation, comments] = LoadArff(arffFile);
+
+    timeInd = GetAttPositionArff(attributes, c_timeName);
+    xInd = GetAttPositionArff(attributes, c_xName);
+    yInd = GetAttPositionArff(attributes, c_yName);
+    confInd = GetAttPositionArff(attributes, c_confName);
+    xHeadInd = GetAttPositionArff(attributes, c_xHeadName);
+    yHeadInd = GetAttPositionArff(attributes, c_yHeadName);
+
+    % convert target position to cartesian
+    [horTargetRads, verTargetRads] = EquirectToSpherical(targetPos(1), targetPos(2), metadata.width_px, metadata.height_px);
+    targetVec = SphericalToCart(horTargetRads, verTargetRads);
+
+    % find start/end ime position
+    indStart = find(data(:,timeInd) > targetStartTime);
+    assert(~isempty(indStart), 'Could not find provided start time');
+    indStart = indStart(1);
+
+    indEnd = find(data(:,timeInd) > targetEndTime);
+    assert(~isempty(indEnd), 'Could not find provided end time');
+    indEnd = indEnd(1);
+
+    % find mean head and gaze vector
+    gazeVecMean = zeros(3,1);
+    for ind=indStart:indEnd
+        if (data(ind, confInd) < c_confThreshold)
+            continue;
+        end
+
+        % convert gaze to reference vector
+        [horGazeRads, verGazeRads] = EquirectToSpherical(data(ind, xInd), data(ind, yInd), metadata.width_px, metadata.height_px);
+        gazeVec = SphericalToCart(horGazeRads, verGazeRads);
+        gazeVecMean = gazeVecMean + gazeVec;
+    end
+    gazeVecMean = gazeVecMean / norm(gazeVecMean);
+
+    [horRads, verRads] = CartToSpherical(gazeVecMean);
+    [xMean, yMean] = SphericalToEquirect(horRads, verRads, metadata.width_px, metadata.height_px);
+
+    dispX = targetPos(1) - xMean;
+    dispY = targetPos(2) - yMean;
+
+    dispXInd = GetMetaExtraPosArff(metadata, 'displacement_x');
+    newDispX = str2num(metadata.extra{dispXInd,2}) + dispX;
+    metadata.extra{dispXInd,2} = num2str(newDispX, '%.2f');
+    dispYInd = GetMetaExtraPosArff(metadata, 'displacement_y');
+    newDispY = str2num(metadata.extra{dispYInd,2}) + dispY;
+    metadata.extra{dispYInd,2} = num2str(newDispY, '%.2f');
+
+    [data(:,xInd), data(:,yInd)] = UndoLimits(data(:,xInd), data(:,yInd), data(:,confInd));
+    [data(:,xHeadInd), data(:,yHeadInd)] = UndoLimits(data(:,xHeadInd), data(:,yHeadInd), data(:,confInd));
+    for ind=1:size(data,1)
+        % Check if the limits were wrongly undone
+        angle = GetAngle(data(ind,xInd), data(ind,xHeadInd));
+
+        if (angle > 360/3)
+            [data(ind,xInd), data(ind,yInd)] = BringWithinLimits(data(ind,xInd), data(ind,yInd));
+        end
+
+        % Check if difference comes from errors during recordings
+        angle = GetAngle(data(ind,xInd), data(ind,xHeadInd));
+        if (angle > 360/3 && angle < 2*360/3) 
+            data(ind,xInd) = data(ind,xInd) - metadata.width_px/2;
+            [data(ind,xInd), data(ind,yInd)] = BringWithinLimits(data(ind,xInd), data(ind,yInd));
+        end
+
+        data(ind,xInd) = data(ind,xInd) + dispX;
+        data(ind,yInd) = data(ind,yInd) + dispY;
+        [data(ind,xInd), data(ind,yInd)] = BringWithinLimits(data(ind,xInd), data(ind,yInd));
+
+        data(ind,xHeadInd) = data(ind,xHeadInd) + dispX;
+        data(ind,yHeadInd) = data(ind,yHeadInd) + dispY;
+        [data(ind,xHeadInd), data(ind,yHeadInd)] = BringWithinLimits(data(ind,xHeadInd), data(ind,yHeadInd));
+        if(data(ind, confInd) < 0.3)
+            data(ind,xInd) = 0;
+            data(ind,yInd) = 0;
+        end
+
+    end
+
+    SaveArff(saveFile, data, metadata, attributes, relation, comments);
+
+    % Returns angle between two x coordinates
+    function angle = GetAngle(x1, x2)
+        horRads1 = EquirectToSpherical(x1, 0, metadata.width_px, 1);
+        horRads2 = EquirectToSpherical(x2, 0, metadata.width_px, 1);
+
+        angle = abs(horRads1 - horRads2) * 180 / pi;
+    end
+
+	function [x, y] = BringWithinLimits(x, y)
+        if (y < 0)
+            y = -y;
+            x = x + metadata.width_px / 2;
+        end
+
+        if (y > metadata.height_px)
+            y = 2 * metadata.height_px - y - 1;
+            x = x + metadata.width_px / 2;
+        end
+
+        if (x < 0)
+            x = ceil(abs(x)/metadata.width_px)*metadata.width_px + x - 1;
+            %x = metadata.width_px + x - 1;
+        elseif (x > metadata.width_px)
+            x = x - floor(x/metadata.width_px)*metadata.width_px;
+            %x = x - metadata.width_px;
+        end
+    end
+
+    % x,y are list vectors
+    function [xConv, yConv] = UndoLimits(x, y, confidence)
+        xConv = zeros(size(x));
+        yConv = zeros(size(y));
+        countHor = 0;
+        xConv(1) = x(1);
+        yConv(1) = y(1);
+        % remove horizontal transitions
+        for i=2:size(x,1)
+            totConf = confidence(i) + confidence(i-1);
+            if (totConf > 1.7)
+                if (x(i) - x(i-1) > metadata.width_px - 100) 
+                    countHor = countHor - 1;
+                elseif (x(i) - x(i-1) < -(metadata.width_px - 100))
+                    countHor = countHor + 1;
+                end
+            end
+
+            if (countHor > 0)
+                xConv(i) = x(i) + metadata.width_px;
+            elseif (countHor < 0)
+                xConv(i) = x(i) - metadata.width_px;
+            else
+                xConv(i) = x(i);
+            end
+            yConv(i) = y(i);
+        end
+
+        % remove vertical transitions
+        countVert = 0;
+        c_perc = 0.1;
+        for i=2:size(xConv,1)
+            diff = abs(x(i) - x(i-1));
+            yCloseToLim = y(i) > metadata.height_px - 100 | y(i) < 100;
+            
+            if (diff > (1-c_perc)*metadata.width_px/2 && diff < (1+c_perc)*metadata.width_px/2 && yCloseToLim)
+                countVert = countVert + 1;
+            end
+
+            if (mod(countVert,2) == 1)
+                xConv(i) = xConv(i) + metadata.width_px/2;
+                if (yConv(i) >= metadata.height_px/2)
+                    yConv(i) = 2*metadata.height_px  - yConv(i) - 1;
+                else
+                    yConv(i) = -yConv(i);
+                end
+            end
+
+        end
+    end
+end

Perpendicular.m

@@ -0,0 +1,24 @@
+% Perpendicular.m
+%
+% This function returns a vector perpendicular to the input vector.
+
+function [perpVec] = Perpendicular(vec)
+    range = 0.5;
+    vec = vec./norm(vec);
+
+    if (vec(1) == 0)
+        perpVec = [1 0 0]';
+        return;
+    elseif (vec(2) == 0)
+        perpVec = [0 1 0]';
+        return;
+    elseif (vec(3) == 0)
+        perpVec = [0 0 1]';
+        return;
+    end
+
+    perpVec = zeros(3,1);
+    perpVec(1) = rand - 0.5;
+    perpVec(2) = rand - 0.5;
+    perpVec(3) = -(vec(1)*perpVec(1) + vec(2)*perpVec(2))/ vec(3);
+end

Project3dVectors.m

@@ -0,0 +1,41 @@
+% Project3dVectors.m
+%
+% This function takes as input a list of Cartesian vectors and moves them to the
+% place of the equirectangular video with the least distortion. This is
+% achieved by finding the rotation of the mean vector to the center of the
+% video. Then it transforms all the vectors based on this rotation (i.e around
+% the center of the equirectangular video). In this way most traditional
+% algorithms can work directly on the transformed data.
+%
+% Note: in order for this method to work properly the vertical dispersion of the
+% provided vectors has to be relatively small (< 45 degrees) in order to be in
+% the almost linear part of the equirectangular video.
+%
+% input:
+%   vectorList  - list of vectors in spherical coordinates
+%   metadata    - metadata of the ARFF file
+%
+% output:
+%   coords      - 2D x,y coordinates of the projected vectors
+
+function [coords] = Project3dVectors(vecList, metadata)
+    meanVec = mean(vecList,1)';
+
+    videoMidVec = [-1; 0; 0];
+
+    rot = HeadToVideoRot(meanVec, 0, videoMidVec);
+    rot = rot'; % get rotation to middle of video
+
+    coords = zeros(size(vecList,1), 2);
+    maxVerRads = 0;
+    minVerRads = pi;
+    for ind=1:size(vecList,1)
+        curVec = vecList(ind,:);
+        rotVec = RotatePoint(rot, curVec);
+
+        [horRads, verRads] = CartToSpherical(rotVec);
+        maxVerRads = max(maxVerRads, verRads);
+        minVerRads = min(minVerRads, verRads);
+        [coords(ind,1), coords(ind,2)] = SphericalToEquirect(horRads, verRads, metadata.width_px, metadata.height_px);
+    end
+end

Project3dVectorsFov.m

@@ -0,0 +1,58 @@
+% Project3dVectorsFov.m
+%
+% This function projects the 3D eye FOV vectors into the FOV and stores a new
+% file where the x, y, and tilt of the are accounted for.
+%
+% input:
+%   arffFile    - file to process
+%   outputFile  - file to store converted data
+
+function Project3dVectorsFov(arffFile, outputFile)
+    [data, metadata, attributes, relation, comments] = LoadArff(arffFile);
+
+    % Proccess the FOV vectors
+    [eyeFovVec] = GetCartVectors(data, metadata, attributes); % rotation at (-1,0,0) point
+
+    widthFov = str2num(GetMetaExtraValueArff(metadata, 'fov_width_deg'));
+    widthFovRads = widthFov * pi / 180;
+    heightFov = str2num(GetMetaExtraValueArff(metadata, 'fov_height_deg'));
+    heightFovRads = heightFov * pi / 180;
+
+    widthFovPx = str2num(GetMetaExtraValueArff(metadata, 'fov_width_px'));
+    heightFovPx = str2num(GetMetaExtraValueArff(metadata, 'fov_height_px'));
+
+    xFov = zeros(size(data,1),1);
+    yFov = zeros(size(data,1),1);
+    for i=1:size(data,1)
+        [horRads, verRads] = CartToSpherical(eyeFovVec(i,:));
+        if (horRads < 0)
+            horRads = 2*pi + horRads;
+        end
+
+        horRads = horRads - pi; % reference vetor at (-1,0,0)
+        verRads = verRads - pi/2;
+
+        xFov(i) = widthFovPx * (horRads + widthFovRads / 2) / widthFovRads;
+        yFov(i) = heightFovPx * (verRads + heightFovRads / 2) / heightFovRads;
+    end
+    
+   
+    metaOut = metadata;
+    metaOut.width_px = 0;
+    metaOut.height_px = 0;
+
+    attOut = {};
+    dataOut = zeros(size(data,1),0);
+    timeInd = GetAttPositionArff(attributes, 'time');
+    [dataOut, attOut] = AddAttArff(dataOut, attOut, data(:,timeInd), attributes{timeInd,1}, attributes{timeInd,2});
+    [dataOut, attOut] = AddAttArff(dataOut, attOut, xFov, 'x', 'Numeric');
+    [dataOut, attOut] = AddAttArff(dataOut, attOut, yFov, 'y', 'Numeric');
+    confInd = GetAttPositionArff(attributes, 'confidence');
+    [dataOut, attOut] = AddAttArff(dataOut, attOut, data(:,confInd), attributes{confInd,1}, attributes{confInd,2});
+
+    relOut = 'gaze_fov';
+    commentsOut = comments;
+
+    SaveArff(outputFile, dataOut, metaOut, attOut, relOut, commentsOut);
+
+end

README.md

@@ -0,0 +1,7 @@
+
+# Matlab 360-degree Utilities
+
+This repo contains various Matlab utilities for 360-degree eye tracking data processing.
+In order for these utilities to work clone first the [matlab\_utils](https://github.com/IoAg/matlab_utils) repository which contains utiloities for processing ARFF formatted data.
+
+All files are licensed under GPLv3 unless stated otherwise.

Rodrigues.m

@@ -0,0 +1,16 @@
+% Rodrigues.m
+%
+% Implements Rodrigues' formula for creating a rotation matrix
+
+function [mat] = Rodrigues(vec, angle)
+    mat = zeros(3,3);
+	mat(1,1) = (1-cos(angle))*vec(1)*vec(1) + cos(angle);
+    mat(1,2) = (1-cos(angle))*vec(1)*vec(2) - sin(angle)*vec(3);
+    mat(1,3) = (1-cos(angle))*vec(1)*vec(3) + sin(angle)*vec(2);
+    mat(2,1) = (1-cos(angle))*vec(2)*vec(1) + sin(angle)*vec(3);
+    mat(2,2) = (1-cos(angle))*vec(2)*vec(2) + cos(angle);
+    mat(2,3) = (1-cos(angle))*vec(2)*vec(3) - sin(angle)*vec(1);
+    mat(3,1) = (1-cos(angle))*vec(3)*vec(1) - sin(angle)*vec(2);
+    mat(3,2) = (1-cos(angle))*vec(3)*vec(2) + sin(angle)*vec(1);
+    mat(3,3) = (1-cos(angle))*vec(3)*vec(3) + cos(angle);
+end

RotatePoint.m

@@ -0,0 +1,13 @@
+% RotatePoint.m
+%
+% Rotates a point by multiplying from the left with the provided rotation matrix.
+
+function rotVec = RotatePoint(rot, vec)
+    if (size(vec,2) > size(vec,1))
+        vec = vec';
+    end
+    assert(size(vec,1) == 3 && size(vec,2) == 1, 'Vector should have exactly 3 elements');
+    assert(size(rot,1) == 3 && size(rot,2) == 3, 'Rotation matrix has incorrect dimensions');
+
+    rotVec = rot * vec;
+end

RotationVecVec.m

@@ -0,0 +1,23 @@
+% RotationVecVec.m
+%
+% Calculates the rotation matrix between 2 vectors based on the Rodrigues' formula.
+
+function [rot] = RotationVecVec(vec1, vec2)
+    assert(size(vec1,1) > size(vec1,2), 'Not using column vectors');
+    assert(size(vec1,2) == 1, 'Not using column vectors');
+    assert(size(vec1,1) == size(vec2,1), 'Vector dimensions mismatch');
+    assert(size(vec1,2) == size(vec2,2), 'Vector dimensions mismatch');
+
+    vec1 = vec1 / norm(vec1);
+    vec2 = vec2 / norm(vec2);
+
+    if (sum(abs(vec1 + vec2) < 0.000001) == 3) % opposite vectors
+        midVec = Perpendicular(vec1);
+    else
+        midVec = vec1 + vec2;
+    end
+
+    midVec = midVec/norm(midVec);
+
+    rot = Rodrigues(midVec, pi);
+end

SphericalToCart.m

@@ -0,0 +1,10 @@
+% SphericalTocart.m
+%
+% This function converts spherical coordinates to cartesian coordinates.
+
+function [vec] = SphericalToCart(horRads, verRads)
+    vec = zeros(3,1);
+    vec(1) = sin(verRads) * cos(horRads);
+    vec(2) = cos(verRads);
+    vec(3) = sin(verRads) * sin(horRads);
+end

SphericalToEquirect.m

@@ -0,0 +1,27 @@
+% SphericalToEquirect.m
+%
+% This function converts spherical coordinates to equirectangular coordinates
+% (pixels).
+
+function [xEq, yEq] = SphericalToEquirect(horRads, verRads, widthPx, heightPx)
+    xEq = (horRads / (2 * pi)) * widthPx;
+    yEq = (verRads / pi) * heightPx;
+
+    if (yEq < 0)
+        yEq = -yEq;
+        xEq = xEq + widthPx/2;
+    end
+
+    if (yEq >= heightPx)
+        yEq = 2 * heightPx - yEq -1;
+        xEq = xEq + widthPx/2;
+    end
+
+    if (xEq < 0)
+        xEq = widthPx + xEq - 1;
+    end
+
+    if (xEq > widthPx)
+        xEq = xEq - widthPx;
+    end
+end

TrackingQuality.m

@@ -0,0 +1,65 @@
+% TrackingQuality.m
+%
+% This function gets as input the position of a target and a time interval
+% during which the target was fixated. It then returns the divergence of gaze
+% from the target in degrees.
+%
+% NOTE: Be careful to provide a time interval in which gaze is not noisy
+%
+% input:
+%   arffFile        - path to ARFF file
+%   targetStartTime - target diplay starting time in us
+%   targetEndTime   - target diplay starting time in us
+%   targetPos       - [x, y] equirectangular position
+%
+% output:
+%   divergence      - divergence from target in degrees
+
+function [divergence] = TrackingQuality(arffFile, targetStartTime, targetEndTime, targetPos)
+    if (targetStartTime < 0 || targetEndTime < 0)
+        divergence = -1;
+        return;
+    end
+    c_timeName = 'time';
+    c_xName = 'x';
+    c_yName = 'y';
+    c_confName = 'confidence';
+    c_confThreshold = 0.8;
+
+    [data, metadata, attributes, relation, comments] = LoadArff(arffFile);
+
+    timeInd = GetAttPositionArff(attributes, c_timeName);
+    xInd = GetAttPositionArff(attributes, c_xName);
+    yInd = GetAttPositionArff(attributes, c_yName);
+    confInd = GetAttPositionArff(attributes, c_confName);
+
+    % convert target position to cartesian
+    [horTargetRads, verTargetRads] = EquirectToSpherical(targetPos(1), targetPos(2), metadata.width_px, metadata.height_px);
+    targetVec = SphericalToCart(horTargetRads, verTargetRads);
+
+    % find start/end ime position
+    indStart = find(data(:,timeInd) > targetStartTime);
+    assert(~isempty(indStart), 'Could not find provided start time');
+    indStart = indStart(1);
+
+    indEnd = find(data(:,timeInd) > targetEndTime);
+    assert(~isempty(indEnd), 'Could not find provided end time');
+    indEnd = indEnd(1);
+
+    % find mean gaze vector
+    gazeVecMean = zeros(3,1);
+    for ind=indStart:indEnd
+        if (data(ind, confInd) < c_confThreshold)
+            continue;
+        end
+        % convert gaze to reference vector
+        [horGazeRads, verGazeRads] = EquirectToSpherical(data(ind, xInd), data(ind, yInd), metadata.width_px, metadata.height_px);
+        gazeVec = SphericalToCart(horGazeRads, verGazeRads);
+        gazeVecMean = gazeVecMean + gazeVec;
+    end
+    
+    vecNum = indEnd - indStart + 1;
+    gazeVecMean = gazeVecMean / norm(gazeVecMean);
+
+    divergence = GetDispersion(targetVec, gazeVecMean');
+end

YZXrotation.m

@@ -0,0 +1,22 @@
+% YZXrotation.m
+%
+% This function calculates the rotation matrix to the coordinate system of
+% reference.
+
+function rot = YZXrotation(vec, tiltRads)
+    theta = asin(vec(2));
+    psi = 0;
+    if (abs(theta) < pi/2 - 0.01)
+        psi = atan2(vec(3), vec(1));
+
+    rot = zeros(3);
+    rot(1,1) = cos(theta)*cos(psi);
+    rot(1,2) = -sin(theta);
+    rot(1,3) = cos(theta)*sin(psi);
+    rot(2,1) = cos(tiltRads)*sin(theta)*cos(psi) + sin(tiltRads)*sin(psi);
+    rot(2,2) = cos(tiltRads)*cos(theta);
+    rot(2,3) = cos(tiltRads)*sin(theta)*sin(psi) - sin(tiltRads)*cos(psi);
+    rot(3,1) = sin(tiltRads)*sin(theta)*cos(psi) - cos(tiltRads)*sin(psi);
+    rot(3,2) = sin(tiltRads)*cos(theta);
+    rot(3,3) = sin(tiltRads)*sin(theta)*sin(psi) + cos(tiltRads)*cos(psi);
+end