Optotrak Motion Tracking System

From REALab Wiki
Jump to navigation Jump to search

Our NDI Certus Optotrak motion tracking system uses active markers to transmit data to a ceiling mounted 3 sensor camera system. It is capable of very fast motion data collection, with 1000Hz + sampling rate and .01mm resolution. (Note though that the strober box tends to heat up with higher sampling rates, so 200Hz tends to be practical and sufficient for reach trajectory studies).

Connecting the Optotrak

To turn on the optotrak, flick the switch on the cameras (located between the two plugs), and press the on button at the front of the optotrak box. (Make sure the lens caps are not covering the cameras!)

The connections are as follows:

Markers (make sure to plug them into slots 1/2, 3/4, in the right order) -> strober box -> optotrak control box (make sure to plug into the correct port, e.g, 1, 2, or 3) -> power + USB to computer.

If the system is working correctly, you should see a green light on the strober box, and a green light above the connected port on the control box.

Using NDI Software

You can use the NDI software NDI First Principles to check that the markers are working and to get a sense of the trackable volume of the cameras.

To ensure everything is connected properly: utilities > query system You should see "camera c3 - 3 sensors", if you have 3 pairs of sensors attached to the marker hub and everything is properly connected.


To get feel for range of trackable space (green light = marker seen, vs. red light = marker obstructed):

new experiment : next > next > ports (front of box) + how many markers 1 + 6 markers

cal frequency = 200hz for rapid reaching power frequency- leave alone

finish

If the system is working correctly, you should see a green light on the strober box, and a green light above the connected port on the control box.

Now you can move the markers around to test where your capture volume is - green light indicates a marker is visible, red light means it is obstructed.

Optotrack Data Collection (Matlab Code)

There are a number of matlab scripts that work with the optotrak to read in and save motion tracking data. These are located on the real room 1 PC in the Documents\MATLAB\Optotrak Matlab code directory, and have been written by Craig Chapman. These scripts also draw on additional matlab files in the \optomfiles subdirectory, which were acquired by Craig but written by someone else (so he likely can't answer all questions related to those scripts). I've summed up some of the more important scripts here with links to the scripts. The optomfiles can also be downloaded here: optomfiles.zip


Summary of selected matlab scripts:

runshapempt.m This is the main script for one of Craig's optotrack experiments. It calls on other optotrack scripts to do many things, such as establishing initializing the optotrack, setting up global and local coordinate frames, reading in marker values, and ultimately these marker values are saved in output files. The output files are then read into analysis scripts (see the next section) for statistical analysis of reaches and graphs.


Setup (lines 17-26):

    %SETUP OPTOTRAK
    %initialize optotrak
    optotrak_init

    %define coordinate frame for this experiment
    transform_info = optotrak_transform_realTime_dynamic;
    T_glob_loc = transform_info.T_glob_loc;

    %setup collection parameters for this experiment, which activates the markers
    setup_info = setup_optotrak

optotrak_init.m This script was written by Craig to set up the optitrack and is called at the beginning of the session for each subject. This script creates a local coordinate frame (dynamic and static), by looking at the frontmost marker and asking for 3 points: 0,0, then some positive x coordinate (establishes the x-axis) and then some positive y coordinate (establishes the y-axis).

T_glob_loc: This line translates the optotrak coordinates from a global optotrak coordinate frame to a local table coordinate frame.

setup optotrak.m This sets up the optotrak with values for number of markers used in the experiment, sampling frequency of the markers, etc. You can edit this script if you need to change these values.


Get start position (line 50 - 88):

    %get info about the start position - participants will have to be within
    %5cm of this point to start the first trial
    reply='N';
    while ~strcmpi(reply,'Y')

        uiwait(msgbox('Press OK to collect initial start pos data','Collect?','modal'));

        %initialize and read one frame from all markers, but call marker 1 as start
        %pos info
        puFrameNumber  = 0;
        puElements=0;
        puFlags=0;
        pDataDest=single(zeros(1,setup_info.numMarkers*3)); %nMarkers x 3 (x,y,z for each)
        returnFlags=0;

        [returnFlags, puFrameNumber,puElements,puFlags,pDataDest] = DataGetNext3D(puFrameNumber,puElements,puFlags,pDataDest);

        [new_startPos]=transform4(T_glob_loc, pDataDest(1:3))

        Question = {'Start pos data:'...
            ''...
            'Marker 1: ' int2str(new_startPos)...
            ''...
            'Accept?'...
            ''};
        ButtonName = questdlg(Question, 'Start Pos Check', 'Yes', 'No', 'Quit', 'Yes');

        switch ButtonName
            case 'Yes'
                reply = 'Y';
                startPos = new_startPos;
            case 'No'
                reply = 'N';
            case 'Quit'
                %quit opto too?
                %need to pass quit flag so this actually works?
                clear all; close all; Screen('CloseAll'); 
                disp('*** Exiting Program ***');
                return
        end
    end
    matData.startPos{1} = startPos;

datagetnext3d.m This is one of the optomfiles scripts. It gets the next available position of marker and stores it in a variable (by convention Craig calls this variable "pdtadest" - Note that the optotrack works in single precision variable types, whereas matlab likes to work in double variable types - it is important to initialize "pdtadest" as a single variable in your matlab script.) You can play around with the position of your markers and check that the values of pdtadest make sense - Note that the optotrak by convention calls any missing value a very high negative number- you can change these high negative values to "NaN" in matlab to make this clearer.

OptotrakActivateMarkers(); (line 377): this activates the markers. Make sure to deactivate the markers at the end of your session (by calling the OptotrakDeActivateMarkers() function at the end of your script, e.g., as in line 596).


Check that subject's finger is in the start position (lines 379-399):

        %wait for IR1 to be within startPos range
        inRange = 20; %say 20mm in 3D distance is within range
        inRangeFlag = 0;
        while ~inRangeFlag

            %initialize and read one frame from all markers, but call marker 1 as start
            %pos info
            puFrameNumber  = 0;
            puElements=0;
            puFlags=0;
            pDataDest=single(zeros(1,data_struct.setup_info.numMarkers*3)); %nMarkers x 3 (x,y,z for each)
            returnFlags=0;

            [returnFlags, puFrameNumber,puElements,puFlags,pDataDest] = DataGetNext3D(puFrameNumber,puElements,puFlags,pDataDest);
            curLocation = pDataDest(1:3);
            curLocation=transform4(data_struct.transform_info.T_glob_loc, curLocation);
            curRange = sqrt(sum((curLocation-matData.startPos{1}).^2));

            if curRange < inRange
                inRangeFlag = 1;
            end
        end

These lines of code get the location of the marker, changes the value to double from single, tranforms it to the right coordinate frame, then uses the resulting value to check if the marker is in the starting position.

pDataDest(1:3) denotes the xyz coordinates of only one marker (4:6 would be the next marker?)


Start recording and get onset time (lines 428-443):

        %start recording optotrak
        %initialize the file for spooling
        DataBufferInitializeFile(0,trial_fileName);

        %start collecting data
        DataBufferStart();

        screen('Flip',mainWin,0,1);
        %play a beep
        snd('Play',beep);

        targetOnTime = getsecs;
        
        %wait for movement
        [tooEarlyFlag reachStartedFlag reachOnsetPos] = detectVelOnset_opto_realTime(data_struct.setup_info,data_struct.transform_info.T_glob_loc,rxnTimeLimit,matData.startPos{1});
        reachOnset = getsecs;

The DataBuffer lines start buffering the data via the optotrak box (necessary as we're likely sampling at a very high rate)

detectVelOnset_opto_realTime.m This function determines when the subject starts their reach by storing 4 samples at a time, and calculating whether a velocity threshold is reached over these 4 samples.

reachOnset = getsecs; This is the time of reach onset. You can use this value to determine whether the subject started their reach too early (false start) or too late.


End of trial (lines 582-596):

        %spool data from buffer to file
        puRealtimeData=0;
        puSpoolComplete=0;
        puSpoolStatus=0;
        pulFramesBuffered=0;
        while (puSpoolComplete == 0)
            [puRealtimeData,puSpoolComplete,puSpoolStatus,pulFramesBuffered]=DataBufferWriteData(puRealtimeData,puSpoolComplete,puSpoolStatus,pulFramesBuffered);
            WaitSecs(.1);
            if puSpoolStatus ~= 0
                disp('Spooling Error');
                break;
            end % if
        end

        OptotrakDeActivateMarkers();

        close all;
        %can maybe get error here
        [reach_info normalizedReach] = analyzeReach(trial_fileName, data_struct.transform_info.T_glob_loc, matData.startPos, data_struct.setup_info.frameRate);

This part spools data from buffer to file, makes sure that spooling completes, then deactivate the markers.

analyzeReach.m- This function performs some pre-analysis tasks, such as filtering and filling in missing data.


At the end of each trial Craig saves the data structure in a format that will facilitate easy use of his analysis scripts.


SUMMARY

beginning of experiment:

  • initialize
  • setup frame
  • setup optotrak values

trial loop:

  • activate and start buffering data
  • check marker is at home
  • calculate reach onset time (can determine whether it's a false start or late start)
  • spool data from buffer to file
  • deactivate

example code

Here is one of Rob Whitwell's experiments which uses the above code:

OTCollectBuffer.m



Optotrak Data Analysis (Matlab Code)

You can download the matlab scripts here:

Analyzereachmatlabscripts.zip

If you followed Craig's convention for saving data, then each subject has their own folder with a .mat output file (which includes parameters, trial conditions), and also an opto folder with raw reach data (one file per trial). This raw data is read into the analysis scripts and the results of the analysis scripts get placed back into the .mat file (data_struct.fdaMat).

analyzeReach.m is the main analysis script, it loads in the data and calls other functions. The script is well documented, but I made additional notes (maybe extraneous) that might also be helpful. A detailed guide to functional data analysis of trajectory data is included in the zip file above, but here is a rough overview of the analysis process.


Load reach (lines 5 - 28):

%LOAD REACH
rawData = [];
[fid, message] = fopen(deblank(trial_fileName),'r+','l'); % little endian byte ordering
% Read the header portion of the file
[filetype,items,subitems,numframes,frequency,user_comments,sys_comments,descrip_file,cutoff,coll_time,coll_date,frame_start,extended_header,character_subitems,integer_subitems,double_subitems,item_size,padding] = read_nd_c_file_header(fid);
% Read the data portion of the file
rawData = read_nd_c_file_data(fid, items, subitems, numframes);
fclose(fid);

%TRANSFORM AND RESHAPE DATA
numMarkers = size(rawData,1);
%define value for looking for missing data - see API for details
max_neg = -3e28;

%turn missing values into nans
rawData(rawData<max_neg)=nan;

%convert the raw data (3d matrix with markers by dimensions of space x
%time) to a format suitable for transformation, then transform
%returns cell array where each cell represents the transformed x,y,z
%position data (columns) by time (number of samples) for each marker
for i = 1:numMarkers;
    data{i} = transform4(T_glob_loc, squeeze(rawData(i,:,:))' ); 
end

fopen() This function reads in a byte at a time, read in header, then read in data function

rawData matrix created, with dimensions 6x3x400. 400=data points(samples), 3=xyz values, 6=markers.

Note that the optotrak codes missing values values as -3e28, so it's helpful to first convert these to NaN.

transform4 This function transforms the coordinate frame (using the subject .mat data file which has necessary info for this transform).

So this first part makes each marker a cell, (6 cells in total if you have 6 markers), each with 3 xyz coordinates x 400 samples (200 Hz over 2 seconds in Craig's example).


Interpolation (lines 30 - 38):

%FILL MISSING DATA
toFill = 1;
fillWith = 2;
toKeep = [3 4 5 6];
bounds = {};

%fill the data
[filledData, fillInfo] = fillMissing_inpaint(data,toFill,fillWith,toKeep,startPos,bounds)
reach_info.fill_info = fillInfo;

The raw data has to be edited so that one index finger marker (e.g., marker 2) can fill in gaps from other index finger marker (e.g., marker 1). Next missing samples are interpolated. Bounds are set on this interpolation because interpolation does not work very well at the apex (inflexion) or ends of the reach.

fillMissing_inpaint This function does interpolation while preserving curvature and velocity information (techniques used here are borrowed from image processing techniques).


Filtering (lines 64-68):

%FILTER
toFilter = [1];
cutoff = 10;
reach_info.posLowpassCutoff = cutoff;
filteredData = filterData(filledData,toFilter,cutoff,frameRate);

We need to filter position data, because when look at its derivative (velocity), noisy position segments get magnified.

cutoff=10 10hz is the outer maximum that a person can move, thus this is the low pass filter value.

frameRate= 200, since we sampled at 200Hz (note: don't use "framerate" as a variable name, since it's a PTB function)


Translate and rotate (lines 77-97):

screenIR1 = filledData{2}(1,:);
screenIR2 = filledData{3}(1,:);
screenIR3 = filledData{4}(1,:);
screenIR4 = filledData{5}(1,:);
screenx = [screenIR1(1) screenIR2(1) screenIR2(1) screenIR1(1)]';
screeny = [screenIR4(2) screenIR4(2) screenIR4(2) screenIR4(2)]';
screenz = [screenIR2(3) screenIR2(3) screenIR3(3) screenIR3(3)]';

reach_info.screenPos{1} = [screenIR1;screenIR2;screenIR3;screenIR4;];

%TRANSLATE AND ROTATE - insert later once setup is known
toTranslate = [1];
%try using bottom points
if ~any(isnan([screenIR3 screenIR4]))
    orthoPoints = [screenIR3;screenIR4];
else %use the top points
    orthoPoints = [screenIR2;screenIR1];
end
originIR = 1;
[transRotData] = translateRotate(filteredData, toTranslate, orthoPoints, originIR);

These lines make everything start from the exact same origin point, and then adjusts the rotation of trajectory so that it's perpendicular to the screen (thus we figure out the screen coordinates for this).


Reach onset and offset (lines 101-106):

%DEFINE REACH ONSET AND OFFSET
toDefine = [1];
velLowpassCutoff = 8;
reach_info.velLowpassCutoff = velLowpassCutoff;
[definedReach,defined_reach_info,tangVel] = defineReach_maxY(filteredData,toDefine,velLowpassCutoff,frameRate);
reach_info.defined_reach_info = defined_reach_info;

Next, we extract only the outward movement (i.e., figure out where reach out towards the screen stops).

onset is determined as: for a consecutive number of frames hand had to be going at certain velocity and with certain acceleration

stop is determined as: maximum reach distance is the first cutoff, then go backwards from there to find the point where they went really slow (this info can be found in Craig's papers).

definedReach now just has one marker, all the redundant markers are not processed past this point.



Plotting (lines 108 - 127) - these lines plot the data for an individual trial.



Normalizing (lines 129 - 135):

%NORMALIZE YOUR REACH
toNormalize = [1];
normalizeFrames = 200;
reach_info.normalizeFrames = normalizeFrames;
normalizeType = 1;
reach_info.normalizeType = normalizeType;
normalizedReach = normalizeFDA(definedReach,toNormalize,normalizeFrames,normalizeType,frameRate);

At this point we are done for any single trial, but when you try to average trials, some reaches take up many samples from the optotrak, others take much fewer samples. Thus we have to normalize (equate) different reaches to the same number of samples before we can average the reaches together.

You can normalize in space, time, or xyz values - this script by default normalizes in space, i.e., such that the start and end point in space remain the same (but you can set the normalizeType value to correspond to normalizing in space, time, or xyz). Craig uses functional data analysis to do this normalization: first we fit a mathematical function to the data, then we can alter function to any number of points we need. Next we overfit (e.g., every 4mm) for one direction in space.

This function uses the xyz values, as well as velocity of xyz for the single marker in the definedReach variable.

Once you run the normalizeFDA function, you will have a normalized reach trajectory for every trial, so that trials can then be averaged together.

See the guide included in Analyzereachmatlabscripts.zip for more details on all this, including how to continue with the averaging.

NOTE: You should exercise caution when using these functions. Try to look at individual's reaches and make sure there's not systematic noisy data for lots of individual trials (at least for a couple of subjects).