Cortical thickness values skew to 0

[Jiajia-yj]

@ntustison I used the antspynet.cortical_thickness(...) and the dkt = antspynet.desikan_killiany_tourville_labeling(...) to extract the ROI values from T1 weighted MRI, However, the median of the ROI values is lower than the normal range. I'm wondering if I need to set some threshold for the voxel values when get the ROI values?
The code I use is just the example code:
'''
import ants
import antspynet
import antstorch

t1 = ants.image_read(antspynet.get_antsxnet_data('mprage_hippmapp3r'))
kk = antspynet.cortical_thickness(t1, verbose=True)
ants.plot(t1, overlay=kk['thickness_image'], overlay_alpha=0.75)
'''
However, the cortical thickness value's Peak frequency is around 0.2-0.4 mm, which is very small and skewed to 0. This also happens when I preprocess the raw image from ADNI dataset. The code I used to preprocess the ADNI images is:
'''
# 2. Load MNI152 template and brain mask
# (Consider using res-01 if your raw data is 1mm isotropic)
# ----------------------------------------------------------------
template_path = "tpl-MNI152NLin2009cAsym_res-02_T1w.nii.gz"
mask_path = "tpl-MNI152NLin2009cAsym_res-02_desc-brain_mask.nii.gz"

template  = ants.image_read(template_path)
mni_mask  = ants.image_read(mask_path)
mni_brain = template * mni_mask          

# ----------------------------------------------------------------
# 3. Load raw T1 image and reorient to LPI
# ----------------------------------------------------------------
raw     = ants.image_read(input_path)
img_lpi = ants.reorient_image2(raw, orientation="LPI")

# ----------------------------------------------------------------
# 4. Cortical thickness estimation in NATIVE space  (ANTsPyNet)
# ----------------------------------------------------------------
print("  [Stage 4] Cortical thickness estimation (native space)...")
kk = antspynet.cortical_thickness(img_lpi, verbose=False)

thickness_native = kk["thickness_image"]
seg_native       = kk["segmentation_image"]   

# Binary brain mask derived from segmentation
brain_mask_native = ants.threshold_image(seg_native, low_thresh=1,
                                         high_thresh=None,
                                         inval=1, outval=0)

# Skull-stripped native-space T1
t1_brain_native = img_lpi * brain_mask_native

'''
Also, by using the DKT function, the ROI values are around 0.8, while the values performed by the platform (like ADNI) have a median around 2-3 mm, the code is listed below:
'''
# ----------------------------------------------------------------
# 5. DKT parcellation in NATIVE space (ANTsPyNet)
# ----------------------------------------------------------------
print(" [Stage 5] DKT parcellation (native space)...")
dkt = antspynet.desikan_killiany_tourville_labeling(
img_lpi,
do_lobar_parcellation=True,
version=1,
verbose=False)
dkt_labels = dkt["segmentation_image"]

# ----------------------------------------------------------------
# 6. Extract per-ROI cortical thickness statistics (FIXED ZERO-SKEW)
# ----------------------------------------------------------------
print("  [Stage 6] Computing cortical thickness per DKT ROI...")

# FIX: Intersect DKT labels with an explicit mask of where thickness > 0.
# This keeps non-cortical edge voxels from diluting your Mean/Sigma metrics.
thickness_positive_mask = ants.threshold_image(thickness_native, low_thresh=0.001,
                                               high_thresh=None, inval=1, outval=0)
valid_dkt_labels = dkt_labels * thickness_positive_mask

stats = ants.label_stats(thickness_native, valid_dkt_labels)
df    = pd.DataFrame(stats)
df    = df[df["LabelValue"] > 0].reset_index(drop=True)

csv_path = out("dkt_cortical_thickness.csv")
df.to_csv(csv_path, index=False)
print(f"  Saved → {csv_path}")

'''

[ntustison]

Can you post the original image and the cortical thickness image?

[Jiajia-yj]

Cortical thickness image: 002_S_0295.nii.gz
The raw image is too big, I cannot upload to here.
While I can share the running results for the example code you posted:
'''
import ants
import antspynet
import antstorch

t1 = ants.image_read(antspynet.get_antsxnet_data('mprage_hippmapp3r'))
kk = antspynet.cortical_thickness(t1, verbose=True)
ants.plot(t1, overlay=kk['thickness_image'], overlay_alpha=0.75)
'''

We can see the cortical thickness values peaks around 0.2-0.4

[ntustison]

That's not a good image to use. What are the results for the MNI template?

[Jiajia-yj]

The results for the MNI template are below, I also attached the raw image, and the code I use.

process_with_dkt.py

7505_2020-12-14_I1387299.nii.gz
7505_2020-12-14_I1387299_roi_stats.csv

raw_7505_2020-12-14_I1387299.nii.gz

[Jiajia-yj]

For the MNI image, I registered the CT to the MNI 2mm template, and the cortical thickness image overlaid with the template is :

From the plot, we can see that the cortical thickness map is very sparse and that some regions are missing cortical thickness values. Could this indicate a problem with my pipeline?

[ntustison]

The raw image is not good in that there's no gray matter/white matter contrast. Try the attached image. I just ran it with both kk variants and compared the results with the same processing I did a year ago and they're the the same.

defaced_mprage.nii.gz

[Jiajia-yj]

Thanks @ntustison, I used defaced_mprage.nii.gz you provided. Since you point out the raw image quality is important, do you have any suggestion or criterions (any info or image description?) that can check the images in batch? I will need to preprocess thousands of images, it will be impossible for me to visual check all images.

Back to the image defaced_mprage.nii.gz process:
The distribution of non-zero cortical thickness voxel values looks correct (what's your opinion?), but the DKT ROI still have some vary low ROIs which is concerned. Should I add a threshold in DKT ROI statistics step?

Meanwhile, I also attached the CT registed to MNI 2mm template image below: # CT within the naive space # CT registered to MNI space

The code I use to process the MRI is below:
'''
from templateflow import api as tflow
import ants

#1. Get the T1w Template (193, 229, 193)
template_path = tflow.get('MNI152NLin2009cAsym',
resolution=2,
desc=None,
suffix='T1w',
extension='.nii.gz')

#2. Get the Mask (Matching the above)
mask_path = tflow.get('MNI152NLin2009cAsym',
resolution=2,
desc='brain',
suffix='mask',
extension='.nii.gz')

template = ants.image_read(str(template_path))
mask = ants.image_read(str(mask_path))

img = ants.image_read("defaced_mprage.nii.gz")

#3. get the crotical thickness map

kk = antspynet.cortical_thickness(img, verbose=True)
cortical_thickness = kk['thickness_image']
ants.plot(img, overlay=kk['thickness_image'], overlay_alpha=0.75)

#4. dkt parcellation
dkt = antspynet.desikan_killiany_tourville_labeling(img, do_lobar_parcellation=True, version=1, verbose=True)
dkt_labels = dkt['segmentation_image']
ants.plot(img, overlay=dkt_labels, overlay_alpha=0.75)

#label stats
stats = ants.label_stats(cortical_thickness, dkt_labels)
df = pd.DataFrame(stats)
df = df[df['LabelValue'] > 0]
print(df.head())
print(df.shape)

#5. register to MNI space
mni_t1 = template
mni_mask = mask > 0
mni_brain = mni_t1 * mni_mask

moving_brain = kk["segmentation_image"]
moving_mask = moving_brain > 0
reg = ants.registration(
fixed=mni_brain,
moving=moving_brain,
type_of_transform="SyNRA",
fixed_mask=mni_mask,
moving_mask=moving_mask,
verbose=True,
)
ct_mni = ants.apply_transforms(
fixed=mni_brain,
moving=cortical_thickness,
transformlist=reg['fwdtransforms'],
interpolation='nearestNeighbor'
)
ct_mni = ct_mni * mni_mask
ct_mni_smoothed = ants.smooth_image(ct_mni, sigma=1.0)
ants.plot(ct_mni_smoothed, title="Cortical Thickness in MNI Space ( Smoothed)")

'''

The results images are attached here:
1_cortical_thickness.nii.gz
2_dkt_labels.nii.gz
3_final_ct_mni_2mm.nii.gz
4_final_ct_mni_2mm_smoothed.nii.gz

[ntustison]

Okay, thanks for sending. The cortical thickness and dkt images look correct. Also, you don't need to send any histograms as they don't really help diagnose any potential issues.

Since you point out the raw image quality is important, do you have any suggestion or criterions (any info or image description?) that can check the images in batch? I will need to preprocess thousands of images, it will be impossible for me to visual check all images.

Not really. We haven't had any need with the thousands of images that we've processed. The raw image that you sent was particularly poor. Basically, if you're unable to see the gray/white matter contrast, the pipeline isn't going to find it either.

but the DKT ROI still have some vary low ROIs which is concerned. Should I add a threshold in DKT ROI statistics step?

I don't think that's going to help. If you have an alternative parcellation scheme, you're certainly welcome to use that. However, be aware that this is the approach that we've used for thousands of images and a significant number of papers after working on this pipeline for a long time (prior to 2014). And it's probably a little premature at this point to ask about changes. So my advice would be to process a significant number of subjects and actually looking at the overall trend before trying to change the pipeline.

I will mention, though, that there's a modified thickness pipeline that we developed more recently (along with version 1 of the DKT labeling) that some people prefer which allows you to perform the segmentation step using the following combinations :

• t1 only
• t1, t2
• t1, fa
• t1, t2, fa

Instead of antspynet.cortical_thickness(ț1, ...), one specifies

• t1 only -- antspynet.cortical_thickness([t1, None, None], ...)
• t1, t2 -- antspynet.cortical_thickness([t1, t2, None], ...)
• t1, fa -- antspynet.cortical_thickness([t1, None, fa], ...)
• t1, t2, fa -- antspynet.cortical_thickness([t1, t2, fa], ...)

[Jiajia-yj]

Hi @ntustison, the code is running great and giving me great data. Now I met another prolem that I have a batch of raw T1W images needed to be preprocessed and extract the cortical thickness value for the downstream voxel-wise analysis. So I need to register all images to MNI space. My question is, do you recommend that I register the raw images to MNI templates first and then use antspynet.cortical_thickness() function to extract the CT or extract the CT first then register the CT to MNI space?
Thanks for you help.

[ntustison]

"extract the CT first then register the CT to MNI space". <----This

Registration first could effect the most accurate thickness measurement.

[Jiajia-yj]

Thanks @ntustison, just some more details I might need to ask. Below is my code about register the CT to MNI space. What I did is to wrap the CT map to MNI brain mask, and applied it to CT and then smooth it.
Is this step correct? Or should I warp the brain from the raw image to the MNI brain mask instead? Also does the "SyNRA" transform and "nearestNeighbor" interpolation recommend? And is the smoothness step necessary?

banner("STAGE 4", "Registration & Saving Transforms")

mni_t1 = template
mni_mask = mask > 0
mni_brain = mni_t1 * mni_mask

moving_brain = kk["segmentation_image"]
moving_mask = moving_brain > 0

# Run the registration
reg = ants.registration(
    fixed=mni_brain,
    moving=moving_brain,
    type_of_transform="SyNRA",
    fixed_mask=mni_mask,
    moving_mask=moving_mask,
    verbose=True,
)


# Apply to Cortical Thickness
ct_mni = ants.apply_transforms(
    fixed=mni_brain,
    moving=cortical_thickness,
    transformlist=reg['fwdtransforms'],
    interpolation='nearestNeighbor'
)

# Apply MNI mask and Smooth
ct_mni = ct_mni * mni_mask
ct_mni_smoothed = ants.smooth_image(ct_mni, sigma=1.0)

Thank you so much for your help.

[ntustison]

Is this step correct? Or should I warp the brain from the raw image to the MNI brain mask instead?

It looks correct but it's usually pretty easy to verify by looking at the images to see if the results are expected.

Also does the "SyNRA" transform and "nearestNeighbor" interpolation recommend?

I prefer the family of antsRegistrationSyN* set of transforms but that will probably work. You should look at the set of transform types for more information.

Nearest neighbor is not recommended. You should use "linear" or another interpolator that respects real values.

And is the smoothness step necessary?

Depends on what you're trying to do. If you're going to average the thickness values in a parcellated region or regions, then I would say no but if you're going to perform something like voxel-based morphometry with the set of aligned images then, yes, smoothing is probably important.