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 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.

Thursday, November 9, 2017

not all MB4 fMRI scans are free of gradients and banding

I was encouraged by the images in the previous post: it looked like the gradient and banding ("washing out") artifacts were not visible in the MB4 people. I need to revise that a bit: I do have bands and gradients in at least some MB4 people, though to a much lesser degree than in our MB8 datasets. How much this affects our task analyses is still to be determined.

On the left is the same MB4 person and run from the previous post (27 September 2017), voxelwise standard deviation over time of a task fMRI run, after going through the HCP minimal preprocessing pipelines. I was glad to see that the vessels were brightest (as they should be), though concerned about the frontal dropout. The person on the right is another person scanned at MB4 (doing the same task; same encoding, scanner, etc.); same preprocessing and color scaling for both.

The vessels clearly don't stand out as much in the person on the right. It's hard to tell in the above image, but there's a gradient in the standard deviation and tSNR images, with better signal on the edges than in the center of the brain. Below is another run from the person on the right, tSNR calculated on the images ready to go into afni for GLM calculation (so through the HCP minimal preprocessing pipelines, plus smoothing and voxelwise normalizing). This tSNR image is shown at six different color scalings; it's easier to see interactively (download the NIfTI here), but hopefully it's clear that the darker colors (lower tSNR) spreads from the center of the brain to the edges, rather than uniformly.

Here is the same person and run again, with the standard deviation (first two) and tSNR (right pane) of the raw (no preprocessing, first and lower) and minimally preprocessed (right two) images. I marked a banding distortion with green lines, as well as the frontal dropout. The banding is perfectly horizontal in the raw image, at 1/4 of the way up the image (see larger image), which makes sense, since this is an MB4 acquisition. I included the larger view of this raw image since all three banding artifacts are somewhat visible; in our dataset the inferior band is generally the most prominent.

The banding and gradient artifacts are certainly less visually prominent in our MB4 than our MB8 images, but they are present in some MB4 people. I haven't systematically evaluated (and probably won't be able to soon) all of our participants, so don't have a sense of how often this occurs, or how much it impacts detection of task BOLD (which is of course the key question).

Below the jump: movement regressors for the two runs in the top image; the person with the banding and gradients had very low motion; less than the person from the September post.