I am trying to calculate normalized cross cross-correlations (NormCrossCor) the way the Reproducible Brain Charts (RBC) did, using nilearn.
For instance, the RBC set a criterion for NormCrossCor at >=.8 when they analyzed the dCCNP data (see Developmental Chinese Color Nest Project - RBC ).
As shown below, I used the np.corrcoef() function to calculate NormCorr between
a) the masked subject’s T1 MRI in MNI space and
b) the masked MNI152NLin2009cAsym template.
However, when I tried this approach, T1 MRI from most dCCNP participants did not meet the NormCrossCor >= 0.8 criteria. When the RBC people did, only one participant failed based on either NormCrossCor >= 0.8 or median FD <=.2.
What am I doing wrong?
What I did:
I have downloaded the Chinese Colornest Project (CCNP) bids data from the RepoBrainChart github ( GitHub - ReproBrainChart/CCNP_BIDS: RBC Repository of CCNP BIDS data · GitHub ).
I have processed the bids data using the fMRIprep 25.2.5 container and used I am using nilearn 0.13.1 and nibabel 5.4.0:
apptainer run --cleanenv ${bindpath} \
${image_path} \
${bids_path} \
${out_path} \
participant \
--participant_label ${participant_id} \
--omp-nthreads 4 \
--mem_mb 24000 \
--skip_bids_validation \
--level full \
--fs-license-file ${license_path} \
-w ${work_dir_path} \
--clean-workdir \
--cifti-output \
--notrack \
--fs-subjects-dir ${fs_dir_path}
All of the participants completed successfully.
Here is my snippet of Python code to calculate NCC..
# Get the exact T1w template fMRIPrep uses
# 'desc-preproc' or no description is standard for the T1 image
t1_path = tf.get('MNI152NLin2009cAsym', desc=None, resolution=1, suffix='T1w')
# Get the corresponding brain mask
mask_path = tf.get('MNI152NLin2009cAsym', resolution=1, desc='brain', suffix='mask')
# Load them into nilearn
mni_template_img = image.load_img(t1_path)
mni_brain_mask_img = image.load_img(mask_path)
def calculate_ncc(sub_mni_path, sub_mni_mask_path, mni_template_img, mni_brain_mask_img):
"""
Calculates the Normalized Cross-Correlation (NCC) between a subject
scan and a MNI152NLin2009cAsym template use by fMRIprep.
Standardized by using Z-score normalized strictly within the brain mask
to prevent background voxels from inflating the correlation score. As the
images are defaced we have combined the MNI152NLin2009cAsymtemplate mask
and the subject mask, to create a combined mask.
Pearson correlation mathematically equals Z-scored NCC, ensuring alignment
is evaluated purely on tissue contrast patterns regardless of global scanner
brightness.
"""
# Validation: Ensure the subject file exists before processing
if not sub_mni_path or not os.path.exists(sub_mni_path):
return 0.0
# Load the subject image (standard NIfTI format)
sub_img = image.load_img(sub_mni_path)
# Load the subject mask
sub_mask_img = image.load_img(sub_mni_mask_path)
# Handle Spatial Resampling Safely (Checks shape and spatial grid affine)
if (sub_img.shape != template_img.shape) or not np.allclose(sub_img.affine, template_img.affine):
sub_img = image.resample_to_img(sub_img, template_img, interpolation='linear')
# Always force the subject mask onto the template mask's grid to ensure perfect array alignment
if (sub_mask_img.shape != mask_img.shape) or not np.allclose(sub_mask_img.affine, mask_img.affine):
sub_mask_img = image.resample_to_img(sub_mask_img, mask_img, interpolation='nearest')
# Extract Data Arrays
sub_data = sub_img.get_fdata()
# Get template
temp_data = template_img.get_fdata()
# Generate Boolean Masks (Thresholding handles interpolation rounding)
mask_data = mask_img.get_fdata() > 0.5
sub_mask_data = sub_mask_img.get_fdata() > 0.5
# Mask Intersection
combined_mask = mask_data & sub_mask_data
# Extract 1D vectors strictly within the brain tissue
sub_vector = sub_data[combined_mask]
temp_vector = temp_data[combined_mask]
# Check for dead signals / uniform flat backgrounds
if len(sub_vector) == 0 or np.std(sub_vector) == 0 or np.std(temp_vector) == 0:
return 0.0
# Calculate Pearson Correlation (Equivalent to Z-scored NCC)
correlation = np.corrcoef(sub_vector, temp_vector)[0, 1]
# Handle NaNs just in case
return float(correlation) if not np.isnan(correlation) else 0.0