====== Transforming data across scans/sessions ======

mrLoadRet has been written to make it relatively easy to combine data sets that have been acquired with different different slice prescriptions (and/or different voxel sizes). Sometimes this is done within the same scanning session but more frequently it is because data from the same subject were acquired in separate scanning sessions.

The critical thing is that there is an alignment transform associated with every scan/run. A map that results from analyzing the data inherits the transform. Likewise there is an alignment transform associated with every ROI. These transforms can be extracted from the data structures and used to bring data and ROIs into register. The alignment transforms for each scan/run are generated by [[mrtools:mralign|mrAlign]] and are explained in detail in [[mrtools:coordinatetransforms|Coordinate Transforms]]. The alignment transforms for the ROIs are inherited from the base anatomy that is selected when the ROI is drawn/created.

mrLoadRet enables you to combine data that were acquired (typically from the same subject) on different days. Start with one session and adding/importing to it the data from another session. Use "Edit - Group - New Group" to create a new group to hold the imported data. Then use either "File - Import - Group" or "Edit - Scan - Add Scan" to copy the scans from the other session into this session.

The code below also exemplifies the conventions that we have adopted for writing mrLoadRet scripts and functions (e.g., the use of viewGet to dereference the data structures).

====== Saving and loading an analysis performed elsewhere in mrTools ======
If you performed an analysis using other tools and have an output matrix for the analysis with coordinates X*Y*Z in which each entry correspond to a value returned by the analysis, you have created something that is equivalent to what is called an 'Overlay' in mrLoadRet. The following lines show a way to save this 'Overlay' as a mrTools 'Analysis.' This will allow you to load and show the results of your analysis in mrLoadRet.

<code matlab>
% open a view 
v = newView;

% Load your X*Y*Z analysys
map = load('my-analysis.mat');

% choose on which group and scan the analysis was performed on (the analysis will be saved to that same group and scan): 
scanNum = 1; groupNum = 1;

% save the analysis file
mrDispOverlay(map,scanNum,groupNum,[],'overlayName=my-analysis','saveName=my-analysis-file-name');
</code>

Now you can start mrLoadRet and you should find a file called 'my-analysis-file-name.mat' in the Folder corresponding to group 1. Using the mrLoadRet 'File' menue you can load the analysis in the current view and explor the results. -Franco


===== Transforming the time series =====

This first example is excerpted (with only slight modification) from averageTSeries.m, analysis code for averaging time series data across runs/scans.

<code matlab>
% Create an "offscreen" view data structure (i.e., with no window, no GUI)
v = newView;

% Specify group and scans
groupName = 'Raw';
groupNum = viewGet(v,'groupNum',groupName);
scanNum1 = 1;
scanNum2 = 2;

% Also need to know the number of temporal frames in the time series. Might 
% want to check that both scans have the same number of frames but I didn't 
% do that here.
nFrames = viewGet(v,'nFrames',scanNum1);

% Extract transform from one scan to the other
scan2scan = viewGet(v,'scan2scan',scanNum1,groupNum,scanNum2,groupNum);

% swapXY is needed here, becuase of the way that warpAffine3 works.
swapXY = [0 1 0 0;1 0 0 0;0 0 1 0; 0 0 0 1];

% M is the resulting transform
M = swapXY * scan2scan * swapXY;
    
% Warp the frames
tseries = loadTSeries(v,scanNum2);
tseriesWarped = zeros(size(tseries));
for frame = 1:nFrames
  mrWaitBar(frame/nFrames,waitHandle);
  tseriesWarped(:,:,:,frame) = warpAffine3(tseries(:,:,:,frame),M,NaN,0,interpMethod,scanDims);
end  
</code>

===== Transforming a map/overlay =====

Data analysis code in mrLoadRet typically takes one or more time series as the input and either produces a time series output or it produces one or more "maps" as the output. For example, averageTSeries produces a time series of functional MRI volumes as the average of two or more time series. Likewise, the motion compensation analysis produces a time series output for each of sevearal time series inputs. corAnal, on the other hand, used for analyzing retinotopy data takes one or more time series as inputs and produces several maps (co, amp, ph) as the output. Each map is the same size as one fMRI volume (a single temporal "frame"). Such maps are often called "statistical parameter maps" in the neuroimaging field. We also call them "overlays" because they are typically rendered in pseudo-color (thresholded in some way) on top of a base anatomy.

This example corresponds to a case in which retinotopy data were acquired on one day and an event-related experiment was run on a different day. We want to analyze the event-related data only for certain voxels that correspond to particular retinotopic locations so we want to bring the two data sets into register. See above for how to import the retinotopy data into the event-related session. % In this code example, there are two groups, one for the concatenated event-relate data and the other for the retinotopy data. The radial component of the retinotopy was acquired in several runs of expanding and contracting rings, and averaged across runs. It is these average time series that were imported with the event-related data. Then corAnal was used to compute the co, amp, and ph maps. The ph map is the temporal phase of the fMRI response time series at each voxel, which corresponds to the eccentricity in the visual field that evoked the largest response in that voxel. Likewise, the angular component ph corresponds to the angular component of the visual field that evoked the largest response in each voxel. 

The end result of this code example is to computed two maps called "X" and "Y". Each of these is the size of one volume. They contain, for each voxel, the corresponding x and y locations in the visual field (in units of degrees of visual angle), in register with the event-related data. With this information in hand, you could find all voxels that correspond to a particular range of visual field locations (e.g., within 2 deg of the horizontal meridian, between 5-7 deg eccentricity to the right of fixation) and extract the event-related data for only those voxels. Or you could use [[mrtools:functionreference#mrdispoverlay|mrDispOverlay]] to add the X and Y maps along with whatever maps you compute from the event-related data so that you can readily visualize the relationship between the event-related results and visual field locations.

<code matlab>
% Create an "offscreen" view data structure (i.e., with no window, no GUI)
v = newView;

% Specify groups and scans
groupName = 'Concatenated';
groupNum = viewGet(v,'groupNum',groupName);
groupNameRet = 'Retinotopy';
groupNumRet = viewGet(v,'groupNum',groupNameRet);
angularScanNum = 1;
radialScanNum = 2;

% Get scanDims (volume size) for event-related scan
v = viewSet(v,'curGroup',groupNum);
scanNum = 1;
scanDims = viewGet(v,'scandims',scanNum);

% Get xforms from retintopy data to event-related scan
xformRadial2Scan = viewGet(v,'scan2scan',radialScanNum,groupNumRet,scanNum,groupNum);
xformAngular2Scan = viewGet(v,'scan2scan',angularScanNum,groupNumRet,scanNum,groupNum);

% Load corAnal and transform retinotopy co and ph maps
v = viewSet(v,'curGroup',groupNumRet);
scanDimsRet = viewGet(v,'scandims');
interpMethod = 'nearest';
swapXY = [0 1 0 0;1 0 0 0;0 0 1 0; 0 0 0 1];
Mradial = inv(swapXY * xformRadial2Scan * swapXY);
Mangular = inv(swapXY * xformAngular2Scan * swapXY);
v = loadAnalysis(v,'corAnal/corAnal');
co = viewGet(v,'co');
amp = viewGet(v,'amp');
ph = viewGet(v,'ph');
coRadial = warpAffine3(co.data{radialScanNum},Mradial,NaN,0,interpMethod,scanDims);
ampRadial = warpAffine3(amp.data{radialScanNum},Mradial,NaN,0,interpMethod,scanDims);
phRadial = warpAffine3(ph.data{radialScanNum},Mradial,NaN,0,interpMethod,scanDims);
coAngular = warpAffine3(co.data{angularScanNum},Mangular,NaN,0,interpMethod,scanDims);
phAngular = warpAffine3(ph.data{angularScanNum},Mangular,NaN,0,interpMethod,scanDims);

% Wrap phase of radial component by pi/2 so that 0 corresponds to the fovea.
% This is particular to how the retinotopy data were acquired in this example.
% You may need to do something differently depending on where the stimulus was
% at the beginning of each run.
phRadial = phRadial + pi/2;
id = phRadial > 2*pi;
phRadial(id) = phRadial(id) - 2*pi;

% Flip phase of angular component so that right hor meridian = 0 and upper
% vert meridian = pi/2. This is particular to how the retinotopy data were 
% acquired in this example. You may need to do something differently depending 
% on where the stimulus was at the beginning of each run, and which direction 
% it rotated.
phAngular = 2*pi - phAngular;

% Convert retinotopy measurements from temporal phase to deg visual angle
innerRet = 0.4;   % inner-most radius of expanding rings in deg of vis angle
outerRet = 11;    % outer-most radius of expanding rings
eccentricity = phRadial/(2*pi) * (outerRet-innerRet) + innerRet;
X = eccentricity .* cos(phAngular);
Y = eccentricity .* sin(phAngular);
</code>
===== Finding corresponding voxels for an ROI analysis =====

Normally, you can get the coordinates of an ROI that exists in your ROI directory using loadROITSeries. But there may be situations (like for a classification experiment) in which you want to make sure that you get exactly the same voxels from each scan taken on different sessions. First, a reminder on how you **normally** get roi coordinates:

  v = newView
  scanNum = 1;
  groupName = 'Raw';
  roiName = 'myroi';
  myroi = loadROITSeries(v,roiName,scanNum,groupName);

Then the field:

  myroi.scanCoords  

will contain the coordinates of the ROI for the scan/group that you selected. You will also have a field:

  myroi.tSeries

which will have loaded the time series for each voxel. If you just wanted to get the scan coordinates, but not  the time series, you could do:

  myroi = loadROITSeries(v,roiName,scanNum,groupName,'loadType=none');

Now, the way this all works is that the roi has its own coordinates (roi.coords), its own transformation matrix (xform) and voxel sizes (voxelSize) -- these are inherited from whatever base anatomy you originally drew the ROI on. To get coordinates for the scan, it transforms them (see: [mrtools:coordinatetransforms|coordinate transforms] for details). It does this in a way that attempts to keep the volume of the ROI the same despite possible differences in the voxel size of the ROI and then scan (it uses the function xformROIcoords).

Now because of the way the transformation works, if you look at the scanCoords for your ROI from one session, you may have a slightly different set of coordinates from a scan taken on a different day. This is normal because the head will have been in a different position and the transformations are taking care of finding the same part of cortex. The two scanCoords may not even have the same number of voxels. 

But, there are times when you might want to know for scan day 1, what are the voxels from scan day 2 which exactly **match**. That is, you want to know where exactly is ROI voxel #15 from day 1 on day 2. This is sometimes useful for a classification protocol in which you have calculated voxel weights in the classifier for each voxel and want to use those weights for analyzing the second days data. In that case, you can call loadROITSeries as follows:

  myroi = loadROITSeries(v,roiName,scanNum,groupNum,'matchScanNum=1','matchGroupNum=2','loadType=none');

This will load the roi with its scan coordinates for scanNum/groupNum but make sure the number and order of voxels match that of scan 4, group 2.