thickinthehead fails on large images ***plus discussion about ANTs thickness measures***

[binarybottle]

After fixing a bug in issue #149, the thickinthehead function fails when running on large image volumes. Inside this functon, Mindboggle calls the ANTs ResampleImageBySpacing command to resample segmented and labeled images to twice their resolution for extracting more precise inner and outer cortical edges. If an image starts at, for example::

dimensions: 320 x 320 x 224
voxel sizes: 0.800000, 0.800000, 0.799866

then after resampling, it becomes::

dimensions: 639 x 640 x 448
voxel sizes: 0.400000, 0.400000, 0.399933

The subsequent command is ANTs ImageMath, which Kills the process when it encounters such a large image:

ImageMath 3 inner_edge.nii.gz MD noncortex.nii.gz 1 Killed

I was able to confirm that Mindboggle's thickinthehead produces reasonable results if run on a downsampled version of this scan with maximum 256 voxels in any dimension.

I am not sure how to resolve this other than to skip the thickinthehead calculation for image volumes larger than 256 voxels in any dimension, or to downsample the image before running thickinthehead.

[satra]

does it fail because it runs out of memory? would be good to know: @ntustison any thoughts here?

[binarybottle]

I think it's a memory issue, because the process is Killed for two different images when they get to about 500MB on my laptop. Do you think the best way to resolve this issue is to warn, raise an error, downsample to a size that works (256 mm^3 with 1mm^3 voxel size), or break up the volume somehow into pieces (like explode_scalars but for volumes)?

[satra]

@binarybottle - unless you have surfaces or lines that you are mapping on to voxels, taking two volumes and simply upsampling them will not necessarily get your more precision. what's perhaps more useful is if you have partial volume info to take into account. so i think you can perhaps do all operations at native volume resolution instead of upsampling.

freesurfer has mri_segstats. perhaps ants has something similar that can take image + label (+ partial volume map) and return volume stats per label.

[ntustison]

freesurfer has mri_segstats. perhaps ants has something similar that can take image + label (+ partial volume map) and return volume stats per label.

LabelGeometryMeasures.

I'm somewhat surprised that your problem is due to a morphological operation. Can you post the data so I can run that particular command on the input data?

[binarybottle]

Thank you, @satra and @ntustison!

@satra -- Thank you for suggesting mri_segstats (LabelGeometryMeasures), which would be good to include in the pipeline for additional output. But this doesn't help with the thickinthehead processing. The reason I upsample is so that when I do morphological operations on the segmented images I can extract an outer and an inner layer of voxels for the gray matter for each region, from which I estimate an area of the region. As I mentioned at the top of issue #149, the thickinthehead measure is computed simply as follows for each labeled cortical region:

volume = (number of voxels) x (product of voxel dimensions)
area = (outer edge volume + inner edge volume) / 2
thickinthehead = volume / area

@ntustison -- I have uploaded mindboggle123 (antsCorticalThickness, recon-all, and mindboggle) output for the example RU brain image I used in the above discussion to Open Science Framework: https://osf.io/36gdy/ Please go to examples -> output -> issue149_thickinthehead_RU5739870

[ntustison]

@binarybottle Your thickinthehead measurement sounds a lot like what we put together in ImageMath in discussions with you (called LabelThickness). If that's the case, I thought we also discussed why that isn't the best approach for calculating digital surface areas. Don't you remember these plots?

[binarybottle]

Wow -- I completely forgot about our LabelThickness discussion back in 2013! I just dredged up our extensive email exchange with Brian Avants when I first wrote thickinthehead and compared its results to FreeSurfer and antsCT (9/2013), when you and Brian created LabelThickness as a version of thickinthehead in ANTs (10/2013), and your Christmas gift two months later of LabelThickness2, which I thought you said had resolved the concern related to sampling resolutions in your plot:

12/25/2013:

First, I put ATITH in antsCorticalThickness.sh. You just need to specify the cortical label image using the ‘-r’ option (‘r’ being for the phonetic beginning of Arno’s name—‘a' was already used). This is an optional argument. If the label image isn’t specified or doesn’t exist, we simply go with our global thickness prior.

Second, currently the “LabelThickness” option is the one utilized in antsCorticalThickness.sh. However, I wanted to make sure that “LabelThickness2” was working so I used the Kirby data set to calculate ATITH for all the labels using both “LabelThickness” and “LabelThickness2”. I also did this for different spacing resamplings. So for each data set, I resampled the cortical thickness DKT label image using isotropic sampling resolutions of (0.5, 0.6, 0.7, 0.8, 0.9, 1) mm. I then calculated ATITH and spit that out to a .csv file. I then plotted the results by averaging overall each subjects for each label at each sampling resolution. Both plots for “LabelThickness” and “Label Thickness2” are attached (in addition to the .csv files and R script). LabelThickness2 seems to be working and considering the behavior, we might want to switch out the option in antsCorticalThickness (albeit modified by a scalar factor).

And I forgot that all I need to do to compute LabelThickness2 is to add the -r flag when running antsCorticalThickness.sh!

Your plot of LabelThickness2 across sampling resolutions (12/25/2013): labelThickness2.pdf

Your reproducibility plots and table for antsCT, Freesurfer, and ANTs thickinthehead (ATITH) thickness values (12/26/2013): reproducibilityAntsCTBA.pdf reproducibilityAtithBA.pdf reproducibilityFreesurferBA.pdf reproducibilityTables.pdf

LabelThickness2 computes thickness as: thicknessprior = volume[label] / surface[label] * 2.0 * volumeelement; Where a label's edge and area are computed as:

      for( unsigned int ii = 0; ii < Gsz; ii++ )  {
              if( GHood.GetPixel( ii ) == 0 ) {  isedge = true; }
      }
      if( isedge ) {
              surface[(unsigned long) label] = surface[(unsigned long) label] + 1;
      }

Please excuse my ignorance of the algorithm and of C++, but does this surface computation tally the number of voxels on the pial surface of the label region, on the pial and gray-white surfaces of the label region, or on every surface of the label region including edges abutting other labels? And it computes surface area at the native resolution of the image, yes? If so and if I understand correctly, Mindboggle's thickinthehead does three things differently: it (1) doubles the resolution for a thinner surface, (2) takes the average of the inner and outer surface areas to account for sharp curvature at the gray-white boundary vs. shallow curvature at the pial boundary, and (3) does not multiply by a factor of 2. Why is volumeelement necessary if both the volume and surface are in units of voxels of the same size? Shouldn't both (or neither) be multiplied by the volume of a voxel?

[ntustison]

Yep, all of that material is still in the KellyKapowski GitHub repo.

LabelThickness2 computes thickness as: ...

That's not LabelThickness2. That code snippet is part of the function above it LabelThickness. This is LabelThickness2 which uses this digitally-based surface area estimation technique.

[binarybottle]

Thank you for clarifying! Now it's in two repos for good measure!

I ran the following and get thickness measures of around 1-2mm per region: ImageMath 3 thickness_prior.nii.gz LabelThickness2 cbic_labels.nii.gz

So the code in LabelThickness2, unlike LabelThickness, estimates surface based on the Crofton formula?: https://github.com/midas-journal/midas-journal-852/blob/master/perf3D.cxx#L221-L223 Which surface exactly does LabelThickness2 estimate -- pial, gm/wm + pial, or entire periphery including label boundaries?

[ntustison]

Calculation is generic in that it is performed simply based on the labeled objects in your image.

[binarybottle]

So if I have an image volume in which the whole brain has non-zero values and the outside of the brain has zero-value voxels, will LabelThickness2 compute thickness for a given labeled region based on the surface area of the region in contact with zero-value voxels only?

[ntustison]

That's not quite accurate. For computing surface area of a given object (with a single object being defined by a single label), both LabelThickness2 and LabelGeometryMeasures calculate surface areas by defining "inside" vs. "outside" as "label" vs. "not label" (not as non-zero vs. zero). For example, your DKT labeling protocol, will give you surface areas for each ROI.

[binarybottle]

I see. So if I understand correctly, I have two issues with computing inside vs. outside for a given cortical region:

(1) If LabelThickness2 includes surfaces abutting other cortical regions in its estimate of cortical label surface area, then it will overestimate surface area and underestimate thickness.

(2) If LabelThickness2 includes both pial/non-brain and gray/white surfaces, then it makes sense to divide the surface area by 2 (like thickinthehead does to get the average of the outer and inner boundaries as an estimate of the midline surface area). This would bring the LabelThickness2 measures to the 4mm range. Is this correct?

[ntustison]

Both statements are correct. For (1), though, I would say that your measurement already begins with incredibly simplistic assumptions (i.e., thickness = volume / surface area) and since those measurements are being applied to the cortex, you're probably already over/underestimating your measurements anyway (not to mention the averaging over the entire region). What you're losing in the precise definition of "thickness" is probably more than made up for in the accuracy of the measurement. I would also say that the utility of a measurement, such as cortical thickness, as a biomarker, is not how well it outputs numbers which conform to accepted values but rather how discriminative it is in distinguishing groups and individuals. So, for (2), sure go ahead and divide by 2 (which will admittedly help in presenting results) but that isn't going to do anything for e.g., age prediction. We talked about this previously in the context of the 2014 paper, your recent ATITH paper, and our longitudinal cortical thickness paper. The latter is posted on biorxiv but you might want to look at our updated version which is ready for resubmission.

[binarybottle]

(1) I completely agree that thickinthead makes simplistic assumptions, and I know that most of the surface is spread over the folds and not along the borders between adjacent cortical regions, but including these borders will result in a worse estimate of surface area for smaller regions with proportionately larger borders, and I don't feel comfortable leaving these borders in.

(2) I also agree that if prediction is the goal, then whatever works, works, and that multiplying by 2 won't make a difference. I do think that absolute values that reflect the range of expected anatomical measures do matter for many neuroanatomists, however, so I would prefer to take the average of the outer and inner surface areas rather than lump them together. I apologize for not remembering our discussions about this -- at least my poor memory forces me to revisit things and perhaps experience them anew...

(3) Would it be possible to create a LabelThickness3 that takes into account (1) and (2) above and would be closer to Mindboggle's thickinthehead measure?

(4) Great paper! Thank you for sharing! I know you're probably too far along, but since you've already run the ADNI data through FS and ANTs, could we run the FS and ANTs outputs through Mindboggle as well and compare thickinthehead values computed by Mindboggle?

[binarybottle]

Back to the main program for a minute...

I am considering resampling the segmented image so that it has isometric 1mm^3 voxels, which would be consistent with the original intended input for Mindboggle's thickinthehead algorithm (see https://github.com/nipy/mindboggle/issues/149). I don't trust that a surface estimate will be reliable across different voxel dimensions, but there's this hopeful article from 2015: "Comparison of several surface area estimators for 3D binary images": https://popups.uliege.be/0351-580X/index.php?id=3680&file=1&pid=3676

I tried resampling a segmented image with FreeSurfer and with ANTs and viewing in freeview, and the two results appear different!:

mri_convert -i ru_seg.nii.gz -o ru_seg_1x1x1fs.nii.gz -vs 1 1 1 -rt nearest

ResampleImageBySpacing 3 ru_seg.nii.gz ru_seg_1x1x1.nii.gz 1 1 1 0 0 1

[satra]

@binarybottle - these days people are acquiring anything from 0.5 to 1mm in vivo. wouldn't resampling to 1 reduce the clarity these scans provide? why not leave every file at native resolution for calculation?

in general, other approaches may be useful in the volumetric context as well: https://riojournal.com/article/12346/

[binarybottle]

Thank you, @satra. While the effects of erosion and dilation are benefitted by upsampling, I have removed upsampling while revising the algorithm (see issue #149). While I'm not satisfied with the revision, it is at least less susceptible to artifacts due to voxel dimensions, so I'll test it out before deciding whether to overhaul or remove it entirely. It also doesn't run into occasional memory issues, so I can close this issue.