Calculating CSF volume for native space ventricle ROIs

This is an oddly specific, possibly vague, post that’s really just intended for my personal future reference. It’s a follow-up to this post about calculating gray matter volume from CAT12-produced ROIs and to this post about manually creating and editing ROIs with ITK-SNAP.

My goal here is to calculate total cerebrospinal fluid (CSF) volume for the left and right lateral ventricles in each subject’s native space.

ROI Creation

So, similar to my gray matter volume ROIs, I’ve saved the ventricle ROIs of interest based on the native space Neuromorphometrics atlas outputted by CAT12’s expert mode.

## ------------------------------------------------------------------------------------------------------ ##
## 		                    SAVE NATIVE SPACE ROIs FROM ATLASES				 							 
## ------------------------------------------------------------------------------------------------------ ##

# Use fslmaths to save the specific native space ROIs we're potentially interested in
if [ $save_native_ROIs == 1 ]
then  
   for SUB in "${subjects[@]}"; do
   # Set atlases 
    subjCSFpath=$CSF_folder/${SUB}/
		
   # Lateral ventricle = 25, 26; 
   fslmaths $subjCSFpath/neuromorphometrics_T1_${SUB}.nii.gz -thr 24.99999999 -uthr 25.00000001 -bin 
       $subjCSFpath/neuromorph_LLatVent_${SUB}.nii.gz
   fslmaths $subjCSFpath/neuromorphometrics_T1_${SUB}.nii.gz -thr 25.99999999 -uthr 26.00000001 -bin   
       $subjCSFpath/neuromorph_RLatVent_${SUB}.nii.gz
		
   # Note done with saving native space ROIs 
   echo 'Done saving native space ROIs for subj:' ${SUB}
   echo "Total run time so far: $((SECONDS/60)) minutes and $((SECONDS%60)) seconds";
		
done 
fi 

Next, as described here, I decided that I didn’t like the ventricle coverage of these automatically-piped-out ROIs, particularly for my older adults. So, I decided to hand correct them in ITK-SNAP.

Calculating CSF Volume in the ROIs

Now, I’m just going to go through several iterations of calculating CSF volume in different versions of my ventricle ROIs, to ensure that I know precisely what’s going on — similar to what I did in my gray matter volume calculation post, except abbreviated, because I “know” what I’m doing now. We’ll work with my subject #1001, who is a ‘healthy’ young adult.

1) CAT12’s automated calculation of left lateral ventricular volume

First, let’s start off with CAT12’s calculation based on the Neuromorphometrics atlas labels, warped to the subject’s native space. CAT12 says 6.5076 mL for left lateral ventricle. Note that for ventricles, you need to look at the the CSF spreadsheet (and not, e.g., the gray matter spreadsheet). The CSF spreadsheet is named ROI_catROI_neuromorphometrics_Vcsf.csv by default.

2) Kathleen’s calculation of CSF volume using the “whole” ROI mask

See the end of this post for more complete code to pull CSF volumes. For now, we’ll just focus on what changes between each case.

First, I’ll note that the volume per voxel will be the same in each case. The CAT12-produced CSF segments (from which you create ventricle ROIs and pull ventricular volume) all have the same voxel sizes. This is ~0.80 mm x 0.80 mm x 0.80 mm = 0.512 mm³.

Next, I’ll note that in each case I used the capital -V flag to pull the number of non-zero voxels in the ROI mask. This number will change between each example, but the code for pulling it will be the same.

The primary flag that changes between my different examples is lowercase -m (mean of all voxels in the ROI mask) versus -M (mean of all non-zero voxels in the ROI mask).

So, first I calculated ventricular volume using the “whole” ROI mask outputted by CAT12, to match CAT12’s automated calculation. That is, I used the -m flag to calculate mean intensity based on all voxels in the mask.

meanIntens=$(fslstats $p3_image -k $ROI -m);

This gave me 6.5076 mL of CSF in the left lateral ventricle, which matches the CAT12 output!

Number of voxels in neuromorph_LLatVent in 1001 is: 15781 voxels
Mean intensity for neuromorph_LLatVent in 1001 is: 0.805408

CSF for neuromorph_LLatVent in 1001 is: 6507.593547 mm3
Or: 6.50759354700000000000 mL

And a histogram like this:

You can see that there’s 1,385 voxels with intensity = 0 contributing to this calculation here, so these bring down the mean intensity a bit. But they still come out to produce the same result as the CAT12 spreadsheet because we’re also counting the voxels with intensity = 0 in the -V total voxel count — so, while we’re driving down mean intensity a little, we’re also including more total voxels in the multiplication to get total left lateral ventricle volume in mL.

3) Kathleen’s calculation of CSF volume using the “thresholded” ROI mask

This version goes through the calculation in a way that makes more logical sense to me, though the end result (volume in mL) is identical to version #2 above.

Here, we first threshold the ROI based on the native space CSF map (i.e., the p3*.nii output from CAT12). That is, we make a new ventricle mask that excludes any voxel with intensity = 0. This isn’t going to visually look too much different from the original mask. That’s because (different from the gray matter ROIs), the initial ventricle ROIs only have a relatively small number of voxels with intensity = 0 (i.e., ~1385/15781=~8.78%). Nonetheless, this cleans up the ventricle mask a little (because intensity=0 is definitely not CSF). And, it (imo) makes the volume calculation make more sense.

There’s many ways to do this step, but here’s one way that uses SPM in MATLAB.

function thresh_ROI(subjID, ROI)

% Set up SPM 
addpath /blue/rachaelseidler/share/FromExternal/Research_Projects_UF/spm12

for i=1:length(subjID)
 %% -------------------------------------------------------------------------------------------------------- %%
 %% 		                               GENERAL SET-UP											
 %% -------------------------------------------------------------------------------------------------------- %%
 % Note step
     disp(sprintf('Starting ROI thresholding for: %s',sprintf(subjID{i})))

 % Set up paths
     subjID{i}
     T1_folder = ['/blue/rachaelseidler/share/FromExternal/Research_Projects_UF/GABA_Aging_KH/Preproc_T1_ROIs/',subjID{i}]; 


 %% -------------------------------------------------------------------------------------------------------- %%
 %% 		                               THRESHOLD ROI												
 %% -------------------------------------------------------------------------------------------------------- %%
 % Select which files to input into imcalc 
   ROI
   file1 = spm_select('ExtFPList', fullfile(T1_folder),'^p3T1_.*\.nii$',1)
   file2 = spm_select('ExtFPList', fullfile(T1_folder),ROI)

 % Run imacalc to threshold ROI by retaining ROI mask only over voxels with intensity > 0
   clear matlabbatch
   matlabbatch{1}.spm.util.imcalc.input = [cellstr(file1); cellstr(file2)];
   matlabbatch{1}.spm.util.imcalc.output = [ROI{1},'_thresh0_',subjID{i}];
   matlabbatch{1}.spm.util.imcalc.outdir = cellstr(T1_folder);
   matlabbatch{1}.spm.util.imcalc.expression = 'i2.*(i1>0)';
   matlabbatch{1}.spm.util.imcalc.var = struct('name', {}, 'value', {});
   matlabbatch{1}.spm.util.imcalc.options.dmtx = 0;
   matlabbatch{1}.spm.util.imcalc.options.mask = 0;
   matlabbatch{1}.spm.util.imcalc.options.interp = 1;
   matlabbatch{1}.spm.util.imcalc.options.dtype = 4;
   spm_jobman('run',matlabbatch);
   disp(sprintf('Done ROI thresholding for: %s',sprintf(subjID{i})))
	
end
end

This now gives me a thresholded ventricle mask with all voxels where intensity=0 have been removed.

Now I can run similar code to above, except here I’ll use -M instead of -m to calculate mean intensity.

meanIntens=$(fslstats $p3_image -k $ROI -M);

This again gives us 6.5076 mL of CSF in the left lateral ventricle, which again matches the CAT12 spreadsheet output (hooray!).

Number of voxels in neuromorph_LLatVent_thresh0 in 1001 is: 14458 voxels
Mean intensity for neuromorph_LLatVent_thresh0 in 1001 is: 0.879108

CSF for neuromorph_LLatVent_thresh0 in 1001 is: 6507.593453 mm3
Or: 6.50759345300000000000 mL

And, it gives us a slightly-nicer-looking histogram, without the 1,385-voxel blip at intensity = 0-0.1.

Do note that the mean intensity is different here compared to the above case. You likely will always be concerned with mL rather than mean intensity for ventricular volume, but it is worth noting that this would be the “correct” approach compared to the above one that includes intensity=0 voxels if you wanted to actually use the mean intensity for anything.

Also, note that there are now fewer voxels. In fact, there are (15,781 total voxels) – (1385 voxels with intensity < 0.1) + (62 voxels that had 0 < intensity < 0.1) = 14,458 voxels in the new thresholded mask.

To clarify: I say above that there are ~1,385 voxels with intensity=0 in the mask (see the first histogram above) because that is the number that fell into the 0-0.1 bin. The “correct” way to say this is actually that there are 1,385 voxels with intensity < 0.1. If we threshold based on i1>0, like in the above imcalc code, we’ll still have some voxels with 0 < intensity < 0.1; in this case, we have 62 of such voxels. If we were instead to threshold at something higher, like i1>0.099, then we’d get a histogram without anything in the 0-0.1 bin. However, this introduces some calculation error compared to the CAT12 calculation, which uses a 0 threshold, not a 0.1 threshold. [Anyhow, I just reran like this to convince myself I wasn’t still including 62 intensity=0 voxels in my calculation. Do this for an understanding check, but not for a final calculation.]

TL;DR regarding CAT12 spreadsheet values

So, I’ve just shown that I understand where the CAT12 spreadsheet values are coming from. Cool. The TL;DR of this, however, is that you can go ahead and just use the CAT12 spreadsheet values without paying much attention to what’s happening under the hood, as long as you trust the ROI masks. This is where you may want to consider options for visually checking the ROI masks, and potentially correcting them if they didn’t come out well — particularly for special populations like aging brains.

As I described in more detail here, I wanted to hand correct the ventricle masks. I’ve done this and now want to compare the same two options from above on my hand-corrected masks. Primarily, I’m looking to see how the hand-corrected values compare to the uncorrected ones. And, whether I can replicate the same two methods as above on my hand-corrected masks.

3) Kathleen’s calculation of CSF volume using the “whole” hand-corrected ROI mask

First, I used the whole mask (no thresholding applied first) and the -m method of calculating mean intensity.

meanIntens=$(fslstats $p3_image -k $ROI -m);

This gave me 7.46 mL of left lateral ventricular volume.

Number of voxels in neuromorph_LLatVent_ITKcorr in 1001 is: 16485 voxels
Mean intensity for neuromorph_LLatVent_ITKcorr in 1001 is: 0.883512

CSF for neuromorph_LLatVent_ITKcorr in 1001 is: 7457.124003 mm3
Or: 7.45712400300000000000 mL

And a histogram like this, which indicates 356 voxels with 0 < intensity < 0.1.

In general, the increase in volume with hand correction makes sense (i.e., 6.51 mL to 7.46 mL). My hand-corrected masks are a bit bigger than the uncorrected ones — because the original mask was not grabbing the entire ventricle, as shown below. The red mask is the original CAT12-outputted mask based on inverse warping the Neuromorphometrics atlas, and the green is my hand-corrected mask. (Although I didn’t test this numerically… but probably should, for fun!) this underestimation effect seemed to be even worse for older adults (in this post and the image below, we’re working with a ‘healthy’ young adult).

4) Kathleen’s calculation of CSF volume using the “thresholded” hand-corrected ROI mask

Next, I’ll redo my calculations with thresholding of i1>0 first, just as I did above in #3. Then, when I use the -M flag to calculate mean intensity, I should end up with a volume that is identical to the #3 result.

meanIntens=$(fslstats $p3_image -k $ROI -m);

This is the case! I again get 7.46 mL for the (hand-corrected) left lateral ventricle.

Number of voxels in neuromorph_LLatVent_ITKcorr_thresh0 in 1001 is: 16173 voxels
Mean intensity for neuromorph_LLatVent_ITKcorr_thresh0 in 1001 is: 0.900556
CSF for neuromorph_LLatVent_ITKcorr_thresh0 in 1001 is: 7457.122400 mm3
Or: 7.45712240000000000000 mL

And now I also get a nicer-looking histogram, with a blip of 44 voxels (instead of 356) for the 0-0.1 bin.

TL;DR regarding calculating ventricular volume using hand-corrected masks

I’ll be moving forward with the volumes calculated using the last approach described here (approach #4). That is, I’ve hand corrected my ventricle masks from CAT12 (because I think the automated ones under-estimate the ventricle) >> then thresholded the masks to remove voxels with intensity=0 >> then calculated # of voxels and mean intensity based on this thresholded mask >> and used these values to calculate ventricular volume.

More complete code for pulling ROI volumes

Just refer to the FSL documentation, and the descriptions above, to understand which flags you need.

## ------------------------------------------------------------------------------------------------------ ##
## 				CALC ROI SIZE, MEAN INTENSITY & FW WITHIN MASK									
## ------------------------------------------------------------------------------------------------------ ##
# First, calc size of ROI. In MNI space, this was the same for everyone. Now, in native space, it's not. 
#   Note the capital V. This will output two numbers: 
#     The first is the number of voxels greater than 0. 
#     The second is their volume in cubic millimiters. 
#     The two numbers should be the same if the ROI is 1x1x1mm3.
# Then, use FSL to calc mean intensity of voxels within the ROI mask. 
#   Multiply: (mean intensity) x (# of voxels) x (vol per voxel). 
#     Vol per voxel = 0.8x0.8x0.8 for the CAT12 CSF images
#        (fslhd: 0.80 in each dimension; SPM display: 0.8x0.8x0.8) 
#	 So vol/voxel = 0.512 mm3 
#	 So, for native space, will need to multiply (mean intensity) x (# of voxels) x (vol per voxel)
# Note: will return values in mm3. Divide by 1000 to get values in mL. 
#	E.g., 9011.7 mm3 = 9.01 mL. 
## ------------------------------------------------------------------------------------------------------ ##

if [ $pull_ROI_vols == 1 ] 
then 
   # Delete outfile if it already exists (so that echo >> won't just append to end of an old file) 
    if [ ! -f $outfile ]; then echo "No outfile found; proceeding."; echo ""; 
    else rm $outfile; echo "Removed old outfile; starting fresh"; echo ""; fi

   # Pull CSF vol within ROI mask 
    for SUB in "${subjects[@]}"; do
	# Set ROI & volume per voxel 
	subj_ROI_folder=/"$CSF_folder"/${SUB}
	ROI=$subj_ROI_folder/$ROIname'_'${SUB}.nii*
	p3_image=$subj_ROI_folder/p3T1_${SUB}.nii
	volPerVoxel=0.512
		
	# Calc size of ROI 
	numVoxels=$(fslstats $ROI -V | awk '{print $1}')
	echo "Number of voxels in $ROIname in ${SUB} is:" $numVoxels "voxels"; 
		
	# Calc mean intensity within mask 
	meanIntens=$(fslstats $p3_image -k $ROI -M); 
	echo "Mean intensity for $ROIname in ${SUB} is:" $meanIntens;		
			
	# Multiply mean intensity by number of voxels 
	CSF=$(echo $meanIntens*$numVoxels*$volPerVoxel | bc)
	CSF_mL=$(echo $CSF/1000 | bc -l)
		
	echo "CSF for $ROIname in ${SUB} is:" $CSF  "mm3";
	echo "Or:" $CSF_mL "mL"; echo ""; 

    # Save to output file  
    # 1 mm3 = 0.001 mL; so e.g., get 4,376 mm3 = 4.376 mL 
    echo -e "${SUB}\t$meanIntens\t$numVoxels\t$CSF\t$CSF_mL" >> "$outfile"
			
done
fi 

Update: -m vs. -M if you’re using thresholded images

I realized post hoc that it doesn’t actually matter whether you use -m (all voxels) or -M (non-zero voxels only) if you’ve taken the liberty of thresholding to remove all of the zero voxels (i.e., i1>0). Both will return an identical mean (duh). The case in which it does matter to use -m instead of -M is when you haven’t thresholded to remove 0 voxels and thus you want to mean over all voxels (and you’re including voxels with intensity = 0 in the total voxel count) — as in examples #1 and #3.