Extend the checking of the subcortical segmentation

[mharms]

With these changes, it will now report separately for each structure, and include additional info regarding overlap with the Atlas ROIs.

Example output file (from the same subject that prompted us to add this report in the first place):

Structure,nMissingBrainMaskFromAtlas,nMissingROIFromAtlas,nROIOverlapAtlas,nROIOutsideAtlas
ALL,409,,,
CEREBELLUM_LEFT,155,1304,7405,957
THALAMUS_LEFT,0,217,1071,203
CAUDATE_LEFT,0,303,425,96
PUTAMEN_LEFT,0,606,454,120
PALLIDUM_LEFT,0,129,168,110
BRAIN_STEM,19,468,3004,183
HIPPOCAMPUS_LEFT,28,317,447,139
AMYGDALA_LEFT,44,163,152,87
ACCUMBENS_LEFT,0,89,46,32
DIENCEPHALON_VENTRAL_LEFT,0,150,556,98
CEREBELLUM_RIGHT,81,1279,7865,914
THALAMUS_RIGHT,0,268,980,149
CAUDATE_RIGHT,0,332,423,114
PUTAMEN_RIGHT,0,359,651,102
PALLIDUM_RIGHT,0,93,167,95
HIPPOCAMPUS_RIGHT,19,390,405,92
AMYGDALA_RIGHT,61,175,157,34
ACCUMBENS_RIGHT,0,107,33,40
DIENCEPHALON_VENTRAL_RIGHT,2,210,502,66

Note that the sum of nROIOverlapAtlas and nROIOutsideAtlas is the number of voxels in the subject-specific ROI. And the sum of nMissingROIFromAtlas and nROIOverlapAtlas is the number of voxels in the Atlas ROIs.

Having this info readily available will be valuable I think. For example, in the example above, we can quickly discern that a number of structures had rather poor overlap with the Atlas ROIs (e.g., nMissingROIFromAtlas > nROIOverlapAtlas), such as PUTAMEN_LEFT, both AMYGDALA, and both ACCUMBENS.

[glasserm]

This may be a bit of information overload. The original purpose of the feature was to assess overlap between the brain mask and the atlas ROIs, not the individual structure ROIs and the atlas ROIs, which are expected to have numerous mismatches in all subjects. It is very surprising that the first column is non-zero for more than the brain stem and possibly the cerebellum if this is the dilation subject. How is this occurring?

[mharms]

So, this quantifies the mismatch the individual structure ROIs and the Atlas ROIs, which IMO is the more relevant info anyway, since it provides info on the extent to which the CIFTI voxels are actually informed by "real" underlying data, rather than simply filled in via the dilation process during CIFTI creation.

The first column is non-zero for more than just the brain stem and cerebellum because other structures (e.g., hipp, amyg) also form the border with the brain mask.

[glasserm]

I guess I don't see how this additional info would be useful for evaluating data usability. Do you envision another use for it?

[mharms]

Sure. One could set a threshold of the minimum amount of overlap with the Atlas that one is willing to accept to include that structure from that subject as part of an analysis. Having this available gives us a mechanism to get a handle on how pervasive of an issue this is.

[glasserm]

I am not following. As I have said, it is EXPECTED that there will not be perfect overlap between the atlas and individual ROIs. This is because although nonlinear volume registration does a decent job with subcortical structures, it does not do a perfect job and so there will be some differences in overlap. The purpose of the CIFTI approach to subcortical structures is to use the data from the individual's ROI to populate the data from the atlas ROI, rather than tolerating residual partial voluming of structures after a volume-based alignment. Assuming proper FOV coverage of the structures (and no artifacts like that brain stem case), the subject is usable even if there are many non-zero values in your table. It is not clear to me how the changes have enabled us to better discriminate usable cases and those with FOV issues or major artifacts. The fact that some atlas subcortical structures end up near the edge of the brain mask without major artifacts or FOV issues actually means that the extant table (and feature) is less useful than I thought it would be. Small differences between the brain mask and atlas ROIs are also expected if the atlas ROIs border the brain mask edge in some places. I guess what is really needed for determining data usability is something that detects large artifacts or FOV mismatches that occur over multiple contiguous slices. Thus, there needs to be some mm threshold built into these measures to ensure that small, expected variances are not being counted. Perhaps if you dilated the brain mask 2x the voxel size and measured atlas ROI voxels outside of that, that would be a useful measure.

[coalsont]

I would consider reordering the columns, having the number that is expected to be big when things are good buried in the middle is confusing. Maybe overlap, brainmask, missing from atlas, discarded from individual?

[mharms]

Can we please focus on the accuracy of the code? The changes are simply quantifying the degree of overlap, so that we have a precise definition of it. I never said that the overlap was expected to be perfect, but I think we can agree that it also shouldn't be "awful". Currently, we have no way to quantify this. My additions provide it.

[coalsont]

@glasserm I think part of what Mike wants is to quantify how much of the data in the atlas cifti space is just dilation. A few voxels isn't a problem, yes, but until now we have basically been hiding this fact from the user entirely.

[glasserm]

It isn't clear why that needs to be a part of the pipelines versus some script Mike runs on some data to test things (I write such scripts all the time and only rarely do they become pipeline additions). I did think the original purpose of the code was sensible--to identify unusable subjects with major issues, but I don't think the implementation currently will do a good job at that. I suggested above an approach to improve upon this.

[coalsont]

We haven't yet decided what the threshold for unusable for this issue is, even in our own data. Dilating would only hide some of the non-overlap, making the measure less accurate and sensitive. Different users may want a different threshold for usability, or may even want to try it as a confound regressor.

"It's not zero" isn't the only way to handle a measurement.

[glasserm]

The code above will not help in making a threshold determination either. I don't think the overall set of ideas is mature enough to be included in the HCP Pipelines at this time, as it does not serve the originally intended purpose (flag bad subjects) and the purpose it appears to be attempting to serve is really better done with a separate experimental script. Alternatively, a 2 voxels away threshold is reasonable and I think would be useful for the brain mask, if there was a desire to implement something now. I am less convinced of the utility of the individual versus atlas ROI analysis, but I suppose a 2 voxels away threshold could be used there.

[mharms]

I believe that the current file provides all the information necessary to evaluate the overlap thoroughly. It also appropriately exposes the overlap to users that may care about this issue.

[glasserm]

It doesn't say how far the voxels are away, and that is CRITICAL to interpreting the information.

[mharms]

I can add code that additionally adds nMissing2xDilROIFromAtlas -- the number of voxels missing from the Atlas ROIs after a 2x dilation of the subject specific ROIs. Is that a reasonable compromise?

I don't see value in a 2x dilation on the brain mask. In fact, once I got into code, and thought about it more, I don't see much value in the nMissingBrainMaskFromAtlas variable at all. Given the parcel-constrained nature by which the CIFTI is created, it is irrelevant if the Atlas ROI voxels have spatial coverage by subject-specific voxels labelled with anything other than their specific associated Atlas structure. So, I'd actually be in favor of removing that variable -- I left it in just to maintain the calculation that you had previously added.

[glasserm]

If we have 2x dilation for both the ROIs and the brain mask, I guess the code is okay for inclusion, as it at least could be used for the original purpose.

[mharms]

What is the purpose of comparing to the brain mask at all? It serves no purpose.

[mharms]

I'm about ready to test a revision. Do I have your blessing to drop the brainmask comparisons? If not, please provide a rationale as to what purpose they serve.

[glasserm]

The brainmask comparisons are the only numbers I would look at from your code and would address whether a subject has issues or not regarding usability. I don't intend to accept this pull request without them. I have explained their purpose repeatedly above.

[mharms]

No you haven't. Not in the face of my point that "Given the parcel-constrained nature by which the CIFTI is created, it is irrelevant if the Atlas ROI voxels have spatial coverage by subject-specific voxels labelled with anything other than their specific associated Atlas structure". If half of the Atlas ROI for hippocampus is "covered" by voxels labelled as amygdala in the subject, those are dilation voxels for the purpose of setting the grayordinate values.

[glasserm]

The brainmask flags overall artifact or FOV issues that would warrant closer inspection. Segmentation or registration variances are a different class of issues that will likely be much more common and much less concerning.

[mharms]

It baffles me how 50% of a structure being based on dilated values is "less concerning", but in the interest of making everyone satisfied, I will dilate the brainmask and include your additional requested variable. We can each then focus on what we think is most relevant.

[mharms]

Updated output. First column of numbers is the one that you requested.

Structure,nMissing2xDilBrainMaskFromAtlas,nMissingBrainMaskFromAtlas,nMissing2xDilROIFromAtlas,nMissingROIFromAtlas,nROIOverlapAtlas,nROIOutsideAtlas
ALL,8,409,,,,
CEREBELLUM_LEFT,1,155,1304,8707,7405,957
THALAMUS_LEFT,0,0,217,1288,1071,203
CAUDATE_LEFT,0,0,303,724,425,96
PUTAMEN_LEFT,0,0,606,775,454,120
PALLIDUM_LEFT,0,0,129,281,168,110
BRAIN_STEM,0,19,468,3471,3004,183
HIPPOCAMPUS_LEFT,0,28,317,749,447,139
AMYGDALA_LEFT,1,44,163,284,152,87
ACCUMBENS_LEFT,0,0,89,118,46,32
DIENCEPHALON_VENTRAL_LEFT,0,0,150,706,556,98
CEREBELLUM_RIGHT,1,81,1279,9141,7865,914
THALAMUS_RIGHT,0,0,268,1241,980,149
CAUDATE_RIGHT,0,0,332,744,423,114
PUTAMEN_RIGHT,0,0,359,969,651,102
PALLIDUM_RIGHT,0,0,93,256,167,95
HIPPOCAMPUS_RIGHT,0,19,390,777,405,92
AMYGDALA_RIGHT,5,61,175,295,157,34
ACCUMBENS_RIGHT,0,0,107,124,33,40
DIENCEPHALON_VENTRAL_RIGHT,0,2,210,710,502,66

[coalsont]

Can "nROIOverlapAtlas" be changed to be last, as it looks like the only number where larger is better?

I find it odd that most of the structures have similar nMissingROIFromAtlas to their overlap (always larger in this output), is this a particularly bad (or small) subject, were the atlas ROIs made using union rather than >50% of subjects, or something else (bug in the code)?

[mharms]

I did it in that order because it seemed more logical to cover all the variables related to missing or overlap first. nROIOutsideAtlas is a bit of a different beast, thence it is at the end . Also, the sum of nMissingROIFromAtlas,nROIOverlapAtlas is equal to the number of voxels in the Atlas ROIs, so it makes sense to me to keep those variables consecutive. i.e., I focused more on the logical structure of the variable ordering. I'm not sure if either organization scheme is a priori better.

The structurals in this test subject weren't so great.

[glasserm]

Do these data exist somewhere I could look at them?

[glasserm]

It does seem in this subject that a lot of the anterior temporal lobes was missed by the FreeSurfer segmentation, though apparently in this run the brainstem was less of an issue. The dilated brain mask might be useful to keep around for any follow up checking of things, rather than having to infer what it would have looked like from the non-dilated mask.

[coalsont]

Thinking about this, the expected nonoverlap doesn't really seem dependent on the voxel size, it depends more on how well the structurals are registered and segmented. If someone were to make a 4mm grayordinates space (because they hate being able to see any detail, presumably), this would dilate 8mm, while if resolution improves to 1.2mm, this would only dilate by 2.4mm.

Of course, in nonhuman subjects, the dilation distance would need to be scaled by brain size, and using the voxel resolution sort of does that.

Another option would be to actually measure the distance of all the non-overlap voxels from the ROI, and average them, but this would need a new wb_command implemented.

[glasserm]

I think it is okay because we are talking about specific standard grayordinates spaces.

[mharms]

@coalsont Are you still reviewing?

[coalsont]

I don't have any more comments on the current PR code.

[mharms]

So, is this ready to merge?

[mharms]

@coalsont @glasserm Anything else? Can I merge this?

[glasserm]

OK