Monday, May 14, 2018

mean pattern subtraction confusion

A grad student brought several papers warning against subtracting the mean pattern (and other types of normalization) before correlational analyses (including, but not only, RSA) to my attention. I was puzzled: Pearson correlation ignores additive and  multiplicative transformations, so how could it be affected by subtracting the mean? Reading closely, my confusion was from what exactly was meant by "subtracting the mean pattern"; it is still the case that not all forms of "normalization" before correlational MVPA are problematic.

The key is where and how the mean-subtraction and/or normalization is done. Using the row- and column-scaling terminology from Pereira, Mitchell, and Botvinick (2009) (datasets arranged with voxels in columns, examples (trials, frames, whatever) in rows): the mean pattern subtraction warned against in Walther et al. (2016) and Garrido et al. (2013) is a particular form of column-scaling, not row-scaling. (A different presentation of some of these same ideas is in Hebart and Baker (2018), Figure 5.)

Here's my illustration of two different types of subtracting the mean; code to generate these figures is below the jump. The original pane shows the "activity patterns" for each of five trials in a four-voxel ROI. The appearance of these five trials after each has been individual normalized (row-scaled) is in the "trial normalized" pane, while the appearance after the mean pattern was subtracted is at right (c.f. Figure 1D in Walther et al. (2016) and Figure 1 in Garrido et al. (2013)).

Note that Trial2 is Trial1+15; Trial3 is Trial1*2; Trial4 is the mirror image of Trial1; Trial5 has a shape similar to Trial4 but positioned near Trial1. (Example adapted from Figure 8.1 of Cluster analysis for researchers (1984) by H. Charles Romesburg.)

And here are matrix views of the Pearson correlation and Euclidean distance between each pair of trials in each version of the dataset:

The mean pattern subtraction did not change the patterns' appearance as much as the trial normalization did: the patterns are centered on 0, but the extents are the same (e.g., voxel2 ranges from 0 to 80 in the original dataset, -40 to 40 after mean pattern subtraction). The row (trial) normalized image has a much different appearance: centered on zero, but the patterns for trials 1, 2, and 3 are now identical, and the range is much less (e.g., voxel2 is now a bit less than -1 to a bit more than 1). Accordingly, the trial-normalized similarity matrix is identical to the original when calculated with correlation, but different when calculated with Euclidean distance; the mean pattern subtracted similarity matrix is identical to the original when calculated with Euclidean distance, but different when calculated with Pearson correlation.

This is another case where lack of clear terminology can lead to confusion; "mean pattern subtraction" and "subtracting the mean from each pattern" are quite different procedures and have quite different effects, depending on your distance metric (correlation or Euclidean).

Sunday, April 15, 2018

more crescents: with sequence variations

Back in January I posted about the "crescent" artifacts that show up in the PA runs of some people. (In a few people reversed crescents appear in the AP runs as well.) The working hypothesis is still that these are N/2 ghosts and perhaps related to insufficient fat suppression. I am still extremely interested in any advice people have for avoiding these, especially as I now have evidence that task signal is noticeably reduced in the "crescent" areas (details soon, hopefully).

We have begun a set of pilot scans to see if some parameter combinations produce better subcortical and frontal signal in a reward task. So far all scans have been on a Siemens Prisma, 64 32 channel head coil, CMRR MB4 sequences; we'll be doing the tests on a Siemens Vida as well in a few weeks. So far, scans are acquired with 2.4 or 3.0 mm isotropic voxels, either "flat" (AC-PC aligned as usual) or "tilted" (30 degrees off AC-PC); more acquisition details below the jump.

Of the two pilot people (so far), one has the crescent artifact and the other does not, and the appearance of the crescents in the different acquisitions is interesting. All of these images are voxelwise standard deviation, calculated over the entire run (no censoring, but extremely low motion), and on raw images (preprocessing is in progress).

First, here are sagittal views of a flat (left, scan 15) and tilted (right, scan 37) 2.4 mm isotropic run. The crescent artifact is visible in both; I marked the approximate ends with green arrows (click to enlarge). The multiband slice boundary is visible in both a fourth of the way up the image (red arrows).

Here are axial slices of the set of runs we have so far for this person. All are with the same color scaling (0 to 200); brighter is higher standard deviation. These are raw images, so the slice appearance varies quite a bit between the "flat" and "tilted" runs.
The "crescents" are visible in all PA runs, though perhaps easiest to spot with the 2p4 (2.4 mm isotropic) voxels. The slices in which the crescents appear varies between the tilted and flat acquisitions (e.g., k=31-46 for run 15_2p4flat_PA; 22-31 for run 37_2p4tilt_PA). It will be easier to compare the crescent locations after preprocessing.

The 3p0 (3.0 mm isotropic voxels) images are generally more uniform and dark than the 2p4 runs, likely reflecting improved signal-to-noise, particularly in the middle of the brain. While the large vessels are brightest in all runs (as they should be), the runs with 2.4 mm voxels (2p4) generally have a "starburst" type effect (brighter in the center, darker towards the edges), which is worrying, particularly since we want good signal in reward areas.

I will share other observations on this blog as the piloting and analyses progress. Please contact me if you'd like to run your own analyses; we'd be happy to share and are very interested in others' thoughts.

UPDATE 18 April 2018: I've wondered before if head size was a factor in which people have the crescent artifact, using the total intracranial volume measurement produced by freesurfer as a proxy. I don't have those measurements yet, but they kindly allowed me to measure their heads as if fitting them for hats, and they were nearly identical: about 58 cm for the person without the crescent artifact, and about 57 for the person shown in this person (with the artifact). Both people have a normal healthy body size; the person without the artifact was a bit shorter (around 5'2") than the other (around 5'7"). So, at least for these two pilots, external head size doesn't seem to matter for the artifact.

more acquisition parameters below the jump

Thursday, February 8, 2018

Connectome Workbench: making a surface version of a volumeric image

This is a tutorial for making a surface version of a volumetric (NIfTI) image (e.g., a ROI) for visualizing with the Connectome Workbench. This post replaces a tutorial I wrote back in 2013, and assumes a bit of familiarity with Workbench. If you're using Workbench for the first time I suggest that you complete the tutorial in this post before trying this one.

First, a warning: as advertised, the steps in this post will make a surface version of a volumetric image. However, the surface version will be an approximation, and likely only suitable for visualization purposes (e.g., making an illustration of a set of ROIs for a talk). If you have an entire dataset that you want to prepare for surface analysis (e.g., running GLMs), you need different procedures (e.g., SUMA, FreeSurfer). Again, I suggest the directions in this post (wb_command -volume-to-surface-mapping) be used cautiously, for quick visualizations, and accompanied by careful confirmation that the mapping produced a reasonable result.

needed files

Before we can make a surface version of a volumetric image, we need to know what it's aligned to, so that we can pick the proper corresponding surface template. Recall that gifti surface files (*.surf.gii) are sort of the underlay anatomy for surface images (e.g., in my Getting Started post on Workbench we load surf.gii files to get a blank brain), so we'll need gifti surface files that will work with our volumetric image (and to serve as the underlay when we're ready to plot the converted volumetric ROI).

For this demo, we'll use fakeBrain.nii.gz (it should let you download without signing in; this is the same NIfTI as shown in other posts), which is aligned to MNI. One MNI dataset with the necessary surface files is the HCP 1200 Subjects Group Average Data; this post describes the files and gives the direct ConnectomeDB download link.

The HCP 1200 Subject Group Average download contains multiple *.surf.gii for each hemisphere, including midthickness, pial, and very_inflated. We can use any of these for visualization in Workbench, but which we pick for the volume to surface conversion does make a difference in what the resulting surface image will look like. It seems best to start with the midthickness surface for the conversion, then try others if the projection seems off.

using wb_command

The wb_command -volume-to-surface-mapping function does the conversion. wb_command.exe (on my Windows machine; file extension may vary) should be in the same directory as wb_view.exe, which you use to start the Workbench GUI. Don't double-click wb_command.exe - it's a command line program. Instead, open up a command prompt and navigate to the directory containing wb_command.exe (on my machine, /bin_windows64/). If you  type wb_command at the prompt it should print out some version and help information; if not, check if you're in the correct directory, and try ./wb_command if you're on a linux-type system.

Now we're ready: we give the function our input NIfTI (fakeBrain.nii.gz), our surface gifti (, the output file we want it to make (demoL.shape.gii), and the options for it to use (-trilinear). Since surface gifti files are just for one hemisphere, we have to do the command twice, once for each. (I included the full path to each file below; update for your machine.)

 wb_command -volume-to-surface-mapping d:/temp/fakeBrain.nii.gz d:/Workbench/HCP_S1200_GroupAvg_v1/ d:/temp/demoL.shape.gii -trilinear  
 wb_command -volume-to-surface-mapping d:/temp/fakeBrain.nii.gz d:/Workbench/HCP_S1200_GroupAvg_v1/ d:/temp/demoR.shape.gii -trilinear  

We now have demoL.shape.gii and demoR.shape.gii, surface versions of fakeBrain.nii.gz, which can be viewed in Workbench (or other gifti-aware programs). Check the surface projection carefully: does the ROI align properly with the anatomy? If not, try a different .surf.gii or fitting option (e.g., -enclosing) in the wb_command call; these can make a big difference.

Below the jump are the surface images from the above commands, plotted on the S1200 midthickness in Workbench.

Wednesday, January 17, 2018

Holy crescents, Batman!

Quite a few of the posts over the last year or so have arisen from things that catch my eye as I review the SMS/MB4 images we're collecting in our ongoing project, and this is another. For quick comparison, I make (with knitr; we may give mriqc a try) files showing slices from mean, standard deviation, and tSNR images for participants, runs, and sessions.

Some participants have obvious bright crescent-shaped artifacts in their standard deviation images (the examples above are from two people; both calculated from non-censored frames, after completing the HCP Minimal Preprocessing pipeline). Looking over people and runs (some participants have completed 6 imaging sessions, over months), people have the crescents or not - their presence doesn't vary much with session (scanning day), task, or movement level (apparent or real).

They do, however, vary with encoding direction: appearing in PA phase encoding runs only. Plus, they seem to vary with subject head size, more likely in small-headed people (large-headed people seem more likely to have "ripples", but that's an artifact for another day).

All that (and thanks to conversations with practiCal fMRI and @DataLoreNeuro) gave a hint: these crescents appear to be N/2 ghost artifacts.

Playing with the contrast and looking outside the brain has convinced me that the crescents do align with the edges of ghost artifacts, which I tried to show above. These are from a raw image (the HCP Minimal Preprocessing pipelines mask the brain), so it's hard to see; I can share example NIfTIs if anyone is interested.

So, why do we have the bright ghosts, what should we do about it, and what does that mean for analysis of images we've already collected? Suggestions are welcome! For analysis of existing images, I suspect that these will hurt our signal quality a little: we want the task runs to be comparable, but they're not in people with the crescent: voxels within the crescent areas have quite different tSNR in the PA and AP runs.

Holy crescents, Batman! (We've been watching the 1966 Batman TV series.)

Wednesday, January 10, 2018

afni: 3dTstat with not-censored timepoints only

In a few recent posts I've shown images of the mean and standard deviation (calculated across time for each voxel), for QC tests. These are easy to calculate in afni (example here), but the 3dTstat command I used includes all timepoints (TRs), unless you specify otherwise. As described previously, we've been using a threshold of FD > 0.9 for censoring high-motion frames before doing GLMs. Thus, I wanted to calculate the mean and standard deviation images only including frames that were not marked for censoring (i.e., restrict the frames used by 3dTstat). This was a bit of a headache to code up, so R and afni code are after the jump, in the hopes it will be useful for others.

Wednesday, November 29, 2017

assigning arbitrary values to surface parcels: method 2

The previous post describes a method for assigning arbitrary values to surface MMP parcels via GIfTI files. Tim Coalson kindly pointed me to another method, which I'll demonstrate here. Both methods work, but one or the other might be easier in particular situations.

In the previous post I used matlab and ran wb_command from the command prompt (on Windows, by opening a command window in the directory with wb_command.exe, then using full paths to the input and output files). Here, I use R, and call wb_command from within R using its system() function. You may need to adjust system for other operating systems, or simply replace it with print and copy-paste the full line to the command prompt.

 rm(list=ls());  # clear R's memory  
 setwd("d:/Workbench/workbench-v1.2.3/bin_windows64/"); # location of wb_command.exe, to call wb_command functions via system() within R  
 local.path <- "d:/temp/";  # local directory for reading and writing files  
 s1200.path <- "d:/Workbench/HCP_S1200_GroupAvg_v1/"; # HCP S1200 Group Average Data Release  
 template.fname <- "MMPtemplate.pscalar.nii";  # filename for the template we'll make  
 # make a template pscalar.nii from the MMP atlas  
 system(paste0("wb_command -cifti-parcellate ", s1200.path, "S1200.thickness_MSMAll.32k_fs_LR.dscalar.nii ",   
        s1200.path, "Q1-Q6_RelatedValidation210.CorticalAreas_dil_Final_Final_Areas_Group_Colors.32k_fs_LR.dlabel.nii COLUMN ",  
        local.path, template.fname));  
 if (!file.exists(paste0(local.path, template.fname))) { stop(paste0("missing: ", local.path, template.fname)); }  
 # you can make a text version of the template, which has 360 rows, but the values in each row aren't the parcel numbers.  
 # system(paste0("wb_command -cifti-convert -to-text ", local.path, template.fname, " ", local.path, "temp.txt"));  
 # the text file needs to be arranged with the 180 right hemisphere parcels in the first 180 rows (in order),   
 # then the 180 parcels for the left hemisphere.  
 # build a text file with the values to plot.   
 out.vec <- rep(0,360);  
 out.vec[c(1,10,15)] <- 1;  # right hemisphere parcels given value 1  
 out.vec[c(180+1,180+10,180+15)] <- 2;  # corresponding left hemisphere parcels given value 2  
 write.table(out.vec, paste0(local.path, "plotDemo.txt"), col.names=FALSE, row.names=FALSE);  
 # create a CIFTI from the text file for viewing in Workbench  
 system(paste0("wb_command -cifti-convert -from-text ", local.path, "plotDemo.txt ", local.path, template.fname, " ", local.path, "plotDemo.pscalar.nii"));  
 if (!file.exists(paste0(local.path, "plotDemo.pscalar.nii"))) { stop(paste("missing:", local.path, "plotDemo.pscalar.nii")); }  

And here is the resulting image, with parcels 1, 10, and 15 plotted in pink on the right hemisphere (1), and yellow on the left (2). Instructions for generating this view are below the jump.

Monday, November 27, 2017

assigning arbitrary values to surface parcels: method 1

This tutorial describes a method for plotting arbitrary per-parcel values (such as from an analysis) on the surface. For example, let's say I want to display MMP parcels 1, 10, and 15 (only) in red, or (more usefully) to assign continuous numbers to each parcel, and then display parcels with larger numbers in hotter colors.

This post describes a method using matlab and creating GIfTI files; see the next post for a method using R, wb_command functions, and creating a CIFTI file. Both methods work, but one or the other might be more convenient in particular situations.

I assume that you have a surface version of the parcellation to use as a template. For example, the MMP parcellation was released in CIFTI format as part of the HCP 1200 release, and the Gordon (2014) parcellation can be downloaded here.

I'll be using the MMP in this example; if you want to follow along, download a copy of the S1200 Group Average Data Release; I put mine at d:/Workbench/HCP_S1200_GroupAvg_v1/. The MMP template is named Q1-Q6_RelatedValidation210.CorticalAreas_dil_Final_Final_Areas_Group_Colors.32k_fs_LR.dlabel.nii. (If you're not sure how to use the files in the S1200 release, try this tutorial to get started.)

I have the values to assign to the parcels in a text file with 180 lines (one line for each MMP parcel). For this tutorial, let's do the simple example of assigning a color to parcels 1, 10, and 15 only. An easy way to do this is to make a text file with 1s on these rows and 0s on all the other rows. I prefer R, but since the GIfTI library is in matlab, here's matlab code for making the little text file:
 out = repelem(0, 180);  
 out(1) = 1;  
 out(10) = 1;  
 out(15) = 1;  
 csvwrite('D:\temp\parcelNums.csv', out')  

The MMP template is in CIFTI format, but we can extract GIfTI files for each hemisphere using wb_command cifti-separate:
 wb_command -cifti-separate D:\Workbench\HCP_S1200_GroupAvg_v1\Q1-Q6_RelatedValidation210.CorticalAreas_dil_Final_Final_Areas_Group_Colors.32k_fs_LR.dlabel.nii COLUMN -metric CORTEX_LEFT D:\temp\mmpL.func.gii  
 wb_command -cifti-separate D:\Workbench\HCP_S1200_GroupAvg_v1\Q1-Q6_RelatedValidation210.CorticalAreas_dil_Final_Final_Areas_Group_Colors.32k_fs_LR.dlabel.nii COLUMN -metric CORTEX_RIGHT D:\temp\mmpR.func.gii  

This matlab code reads in the text file with the new parcel values and the GIfTI templates, then writes GIfTI files with the new parcel values substituted for the parcel label numbers:
 addpath 'C:/Program Files/MATLAB/gifti-1.6';  % so matlab can find the library  
 mmpL = gifti('D:\temp\mmpL.func.gii');  % load the left side gifti MMP atlas  
 mmpR = gifti('D:\temp\mmpR.func.gii');  % and the right side MMP atlas  
 newvals = csvread('D:\temp\parcelNums.csv');  % 180 integers; new value for each parcel  
 Lout = mmpL;  % output gifti  
 Lout.cdata(:,1) = repelem(0, size(mmpL.cdata,1));  % replace the values with zeros  
 Rout = mmpR;  
 Rout.cdata(:,1) = repelem(0, size(mmpR.cdata,1));   
 for i=1:180  % i = 1;  
   inds = find(mmpR.cdata == i);  % find vertices for parcel i  
   Rout.cdata(inds,1) = newvals(i);  % put the new value in for this parcel's vertices  
   inds = find(mmpL.cdata == (i+180)); % in MMP, left hemisphere vertices are 181:360  
   Lout.cdata(inds,1) = newvals(i);    
 save(Lout,'D:\temp\plotL.func.gii','Base64Binary');  % save the gifti  

You can now to plot these GIfTIs in Workbench (see this post if you're not sure how); I plotted them on the S1200 Group Average (MNI) anatomy:
I clicked to put a marker in the left visual parcel. The value at this vertex is 1, as assigned (green lines). I loaded in the MMP atlas as well (blue lines), so it tells me (correctly!) that the marker is in the L_V1_ROI.