Updated cortical thickness and curvature regression

PostFreeSurfer/PostFreeSurferPipeline.sh

@@ -70,6 +70,7 @@ opts_AddOptional '--inflatescale' 'InflateExtraScale' 'number' "surface inflatio
 opts_AddOptional '--processing-mode' 'ProcessingMode' 'HCPStyleData|LegacyStyleData' "disable some HCP preprocessing requirements to allow processing of data that doesn't meet HCP acquisition guidelines - don't use this if you don't need to" 'HCPStyleData'
 opts_AddOptional '--structural-qc' 'QCMode' 'yes|no|only' "whether to run structural QC, default 'yes'" 'yes'
 opts_AddOptional '--use-ind-mean' 'UseIndMean' 'YES or NO' "whether to use the mean of the subject's myelin map as reference map's myelin map mean, defaults to 'YES'" 'YES'
+opts_AddOptional '--metric-regression' 'MetricReg' 'OLD, NEW or BOTH' "whether to use the updated curvature-thickness regression, defaults to 'BOTH'" 'BOTH'
 
 opts_ParseArguments "$@"
 
@@ -281,6 +282,7 @@ if ((doProcessing)); then
     argList+=("$CorrectionSigma")
     argList+=("$RegName")                                  # ${39}
     argList+=("$UseIndMean")
+    argList+=("$MetricReg")                                # ${41}
     "$PipelineScripts"/CreateMyelinMaps.sh "${argList[@]}"
 fi
 

PostFreeSurfer/scripts/CreateMyelinMaps.sh

@@ -93,6 +93,7 @@ ReferenceMyelinMaps="${37}"
 CorrectionSigma="${38}"
 RegName="${39}"
 UseIndMean="${40}"
+MetricReg="${41}"
 
 log_Msg "RegName: ${RegName}"
 
@@ -171,7 +172,6 @@ if [ "${T2wPresent}" = "YES" ] ; then
   ${CARET7DIR}/wb_command -add-to-spec-file "$T1wFolder"/"$NativeFolder"/"$Subject".native.wb.spec INVALID "$T1wFolder"/T1wDividedByT2w_ribbon.nii.gz
 fi
 
-
 MapListFunc="corrThickness@shape"
 if [ "${T2wPresent}" = "YES" ] ; then
   MapListFunc+=" MyelinMap@func SmoothedMyelinMap@func"
@@ -191,7 +191,9 @@ for Hemisphere in L R ; do
     RegSphere="${AtlasSpaceFolder}/${NativeFolder}/${Subject}.${Hemisphere}.sphere.reg.reg_LR.native.surf.gii"
   fi
 
-  ${CARET7DIR}/wb_command -metric-regression "$AtlasSpaceFolder"/"$NativeFolder"/"$Subject"."$Hemisphere".thickness.native.shape.gii "$AtlasSpaceFolder"/"$NativeFolder"/"$Subject"."$Hemisphere".corrThickness.native.shape.gii -roi "$AtlasSpaceFolder"/"$NativeFolder"/"$Subject"."$Hemisphere".roi.native.shape.gii -remove "$AtlasSpaceFolder"/"$NativeFolder"/"$Subject"."$Hemisphere".curvature.native.shape.gii
+  if [ "${MetricReg}" = "OLD" ] || [ "${MetricReg}" = "BOTH" ]; then
+     ${CARET7DIR}/wb_command -metric-regression "$AtlasSpaceFolder"/"$NativeFolder"/"$Subject"."$Hemisphere".thickness.native.shape.gii "$AtlasSpaceFolder"/"$NativeFolder"/"$Subject"."$Hemisphere".corrThickness.native.shape.gii -roi "$AtlasSpaceFolder"/"$NativeFolder"/"$Subject"."$Hemisphere".roi.native.shape.gii -remove "$AtlasSpaceFolder"/"$NativeFolder"/"$Subject"."$Hemisphere".curvature.native.shape.gii
+  fi
 
   #Reduce memory usage by smoothing on downsampled mesh
   LowResMesh="${LowResMeshesArray[0]}"
@@ -217,6 +219,10 @@ for Hemisphere in L R ; do
   done
 done
 
+if [ "${MetricReg}" = "NEW" ] || [ "${MetricReg}" = "BOTH" ]; then
+  $HCPPIPEDIR/global/scripts/CorrThick.sh --subject-dir="$StudyFolder" --subject="$Subject" --regnames="$RegName"
+fi
+
 LowResMeshList=""
 for LowResMesh in "${LowResMeshesArray[@]}" ; do
   LowResMeshList+="${AtlasSpaceFolder}/fsaverage_LR${LowResMesh}k@${LowResMesh}k_fs_LR@atlasroi "

global/scripts/CorrThick.sh

@@ -0,0 +1,139 @@
+#!/bin/bash
+set -eu
+
+pipedirguessed=0
+if [[ ""$HCPPIPEDIR:-"" == "" ]]
+then
+	pipedirguessed=1
+	#fix this if the script is more than one level below HCPPIPEDIR
+	export HCPPIPEDIR="$(dirname -- "$0")/.."
+fi
+
+source "$HCPPIPEDIR/global/scripts/newopts.shlib" "$@"
+source "$HCPPIPEDIR/global/scripts/debug.shlib" "$@"
+
+#description of script/command
+opts_SetScriptDescription "Run curvature-corrected cortical thickness python script and save curvatures, regression coefficients, and curvature-corrected thickness"
+
+opts_AddMandatory '--subject-dir' 'SubjectDir' 'path' "folder containing all subjects"
+opts_AddMandatory '--subject' 'Subject' 'subject ID' "subject-id"
+opts_AddOptional '--regnames' 'RegNamesStr' 'my reg str' "set the desired registration name(s) separated by @, 'string' 'RegName@RegName@[email protected].' default MSMSulc" "MSMSulc"
+opts_AddOptional '--hemi' 'Hemi' 'hemisphere' "provide hemisphere for regression calculation, L=Left, R=Right, or B=Both, default 'B'" "B"
+opts_AddOptional '--surf' 'Surface' 'surface' "provide surface for regression calculation, white or midthickness, default midthickness" "midthickness"
+opts_AddOptional '--patch-size' 'PatchSize' 'distance' "provide patch kernel size in millimeters FWHM for regression, default 6" "6"
+opts_AddOptional '--surf-smooth' 'SurfSmooth' 'distance' "provide surface smoothing in millimeters FWHM, default 2.14" "2.14"
+opts_AddOptional '--metric-smooth' 'MetricSmooth' 'distance' "provide metric smoothing in millimeters FWHM, default 2.52" "2.52"
+opts_AddOptional '--skip-computation' 'SkipCompute' 'YES or NO' "whether or not to compute the curvature-corrected cortical thickness, if it is already available, but just to resample it to 164k and 32k, defaults to 'NO'" "NO"
+
+opts_ParseArguments "$@"
+
+if ((pipedirguessed))
+then
+	log_Err_Abort "HCPPIPEDIR is not set, you must first source your edited copy of Examples/Scripts/SetUpHCPPipeline.sh"
+fi
+
+#display the parsed/default values
+opts_ShowValues
+
+#sanity check boolean strings and convert to 1 and 0
+SkC=$(opts_StringToBool "$SkipCompute")
+
+#set paths
+NonlinearFolder="$SubjectDir"/"$Subject"/MNINonLinear
+NativeFolder="$NonlinearFolder"/Native
+T1wFolder="$SubjectDir"/"$Subject"/T1w/Native
+
+#Make intermediate directory to save intermediate files
+mkdir -p "$NativeFolder"/CorrThick
+
+#Loop through left and right hemispheres
+if [[ "$Hemi" == *B* ]]
+then
+	Hemi="L R"
+fi
+
+RegNames=`echo "$RegNamesStr" | sed s/"@"/" "/g`
+
+LowResMesh="32"
+HighResMesh="164"
+MapListFunc="MRcorrThickness MRcorrThickness_intercept MRcorrThickness_normcoeffs MRcorrThickness_curvs MRcorrThickness_coeffs"
+
+#Generate MRcorrThickness in Native Space
+if ((! SkC)); then
+	for Hemisphere in $Hemi ; do
+		if [[ "$Hemisphere" == "L" ]] ; then
+			Structure="CORTEX_LEFT"
+		elif [[ "$Hemisphere" == "R" ]] ; then
+			Structure="CORTEX_RIGHT"
+		fi
+		(
+			cd "$HCPPIPEDIR"/global/scripts/CorrThick
+			python CorrThick.py "$SubjectDir" "$Subject" "$Structure" "$Hemisphere" "$Surface" "$PatchSize" "$SurfSmooth" "$MetricSmooth"
+		)
+	done	
+fi
+
+#Set the Color Palette(s) and Resample to HighResMesh and LowResMesh
+for Hemisphere in $Hemi ; do
+
+	for Map in $MapListFunc ; do
+		if [[ "$Map" == MRcorrThickness || "$Map" == MRcorrThickness_intercept ]] ; then
+			wb_command -metric-palette "$NativeFolder"/"$Subject"."$Hemisphere"."$Map".native.shape.gii MODE_AUTO_SCALE_PERCENTAGE -pos-percent 4 96 -interpolate true -palette-name videen_style -disp-pos true -disp-neg false -disp-zero false -normalization NORMALIZATION_SELECTED_MAP_DATA
+		elif [[ "$Map" == MRcorrThickness_normcoeffs ]] ; then
+			wb_command -metric-palette "$NativeFolder"/"$Subject"."$Hemisphere"."$Map".native.shape.gii MODE_AUTO_SCALE_ABSOLUTE_PERCENTAGE -interpolate true -palette-name ROY-BIG-BL -disp-pos true -disp-neg true -disp-zero false -normalization NORMALIZATION_ALL_MAP_DATA
+		else
+			wb_command -metric-palette "$NativeFolder"/"$Subject"."$Hemisphere"."$Map".native.shape.gii MODE_AUTO_SCALE_ABSOLUTE_PERCENTAGE -interpolate true -palette-name ROY-BIG-BL -disp-pos true -disp-neg true -disp-zero false -normalization NORMALIZATION_SELECTED_MAP_DATA
+		fi
+	done 
+
+	for RegName in $RegNames ; do
+		RegSphere=""$NativeFolder"/"$Subject"."$Hemisphere".sphere."$RegName".native.surf.gii"
+
+		if [[ "$RegName" == "MSMSulc" ]] ; then
+			RegString=""
+		else
+			RegString="_$RegName"
+		fi
+		for Map in $MapListFunc ; do
+			wb_command -metric-resample "$NativeFolder"/"$Subject"."$Hemisphere"."$Map".native.shape.gii "$RegSphere" "$NonlinearFolder"/"$Subject"."$Hemisphere".sphere."$HighResMesh"k_fs_LR.surf.gii ADAP_BARY_AREA "$NonlinearFolder"/"$Subject"."$Hemisphere"."$Map""$RegString"."$HighResMesh"k_fs_LR.shape.gii -area-surfs "$T1wFolder"/"$Subject"."$Hemisphere".midthickness.native.surf.gii "$NonlinearFolder"/"$Subject"."$Hemisphere".midthickness."$HighResMesh"k_fs_LR.surf.gii -current-roi "$NativeFolder"/"$Subject"."$Hemisphere".roi.native.shape.gii
+			wb_command -metric-mask "$NonlinearFolder"/"$Subject"."$Hemisphere"."$Map""$RegString"."$HighResMesh"k_fs_LR.shape.gii "$NonlinearFolder"/"$Subject"."$Hemisphere".atlasroi."$HighResMesh"k_fs_LR.shape.gii "$NonlinearFolder"/"$Subject"."$Hemisphere"."$Map""$RegString"."$HighResMesh"k_fs_LR.shape.gii
+			wb_command -metric-resample "$NativeFolder"/"$Subject"."$Hemisphere"."$Map".native.shape.gii "$RegSphere" "$NonlinearFolder"/fsaverage_LR"$LowResMesh"k/"$Subject"."$Hemisphere".sphere."$LowResMesh"k_fs_LR.surf.gii ADAP_BARY_AREA "$NonlinearFolder"/fsaverage_LR"$LowResMesh"k/"$Subject"."$Hemisphere"."$Map""$RegString"."$LowResMesh"k_fs_LR.shape.gii -area-surfs "$T1wFolder"/"$Subject"."$Hemisphere".midthickness.native.surf.gii "$NonlinearFolder"/fsaverage_LR"$LowResMesh"k/"$Subject"."$Hemisphere".midthickness."$LowResMesh"k_fs_LR.surf.gii -current-roi "$NativeFolder"/"$Subject"."$Hemisphere".roi.native.shape.gii
+			wb_command -metric-mask "$NonlinearFolder"/fsaverage_LR"$LowResMesh"k/"$Subject"."$Hemisphere"."$Map""$RegString"."$LowResMesh"k_fs_LR.shape.gii "$NonlinearFolder"/fsaverage_LR"$LowResMesh"k/"$Subject"."$Hemisphere".atlasroi."$LowResMesh"k_fs_LR.shape.gii "$NonlinearFolder"/fsaverage_LR"$LowResMesh"k/"$Subject"."$Hemisphere"."$Map""$RegString"."$LowResMesh"k_fs_LR.shape.gii
+		done
+	done
+done
+
+#Create CIFTI Files in Native, HighResMesh, and LowResMesh
+if [[ "$Hemi" == *L* && "$Hemi" == *R* ]] ; then
+	for STRING in "$NativeFolder"@native@roi "$NonlinearFolder"@"$HighResMesh"k_fs_LR@atlasroi "$NonlinearFolder/fsaverage_LR"$LowResMesh"k@"$LowResMesh"k_fs_LR@atlasroi" ; do
+		Folder=`echo $STRING | cut -d "@" -f 1`
+		Mesh=`echo $STRING | cut -d "@" -f 2`
+		ROI=`echo $STRING | cut -d "@" -f 3`
+
+		for RegName in $RegNames ; do
+	  		if [[ "$RegName" == "MSMSulc" || "$Folder" == "$NativeFolder" ]] ; then
+				RegString=""
+			else
+				RegString="_$RegName"
+			fi
+
+			for Map in $MapListFunc ; do
+				wb_command -cifti-create-dense-scalar "$Folder"/"$Subject"."$Map""$RegString"."$Mesh".dscalar.nii -left-metric "$Folder"/"$Subject".L."$Map""$RegString"."$Mesh".shape.gii -roi-left "$Folder"/"$Subject".L."$ROI"."$Mesh".shape.gii -right-metric "$Folder"/"$Subject".R."$Map""$RegString"."$Mesh".shape.gii -roi-right "$Folder"/"$Subject".R."$ROI"."$Mesh".shape.gii
+				if [[ "$Map" == MRcorrThickness || "$Map" == MRcorrThickness_intercept ]] ; then
+					wb_command -cifti-palette "$Folder"/"$Subject"."$Map""$RegString"."$Mesh".dscalar.nii MODE_AUTO_SCALE_PERCENTAGE "$Folder"/"$Subject"."$Map""$RegString"."$Mesh".dscalar.nii -pos-percent 4 96 -interpolate true -palette-name videen_style -disp-pos true -disp-neg false -disp-zero false -normalization NORMALIZATION_SELECTED_MAP_DATA
+				elif [[ "$Map" == MRcorrThickness_normcoeffs ]] ; then
+					wb_command -cifti-palette "$Folder"/"$Subject"."$Map""$RegString"."$Mesh".dscalar.nii MODE_AUTO_SCALE_ABSOLUTE_PERCENTAGE "$Folder"/"$Subject"."$Map""$RegString"."$Mesh".dscalar.nii -interpolate true -palette-name ROY-BIG-BL -disp-pos true -disp-neg true -disp-zero false -normalization NORMALIZATION_ALL_MAP_DATA
+				else
+					wb_command -cifti-palette "$Folder"/"$Subject"."$Map""$RegString"."$Mesh".dscalar.nii MODE_AUTO_SCALE_ABSOLUTE_PERCENTAGE "$Folder"/"$Subject"."$Map""$RegString"."$Mesh".dscalar.nii -interpolate true -palette-name ROY-BIG-BL -disp-pos true -disp-neg true -disp-zero false -normalization NORMALIZATION_SELECTED_MAP_DATA
+				fi
+			done
+			if [[ "$Folder" == "$NativeFolder" ]] ; then
+				break
+			fi
+		done
+	done
+fi
+
+#Remove preliminary directory and all its contents
+rm -rf "$NativeFolder"/CorrThick
+

global/scripts/CorrThick/CorrThick.py

@@ -0,0 +1,42 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Aug 20 14:44:13 2024
+
+@author: brainmappers
+"""
+import sys
+
+subjects_dir=str(sys.argv[1])
+subject=str(sys.argv[2])
+structure=str(sys.argv[3])
+hemi=str(sys.argv[4])
+surface=str(sys.argv[5])
+number=float(sys.argv[6])
+iteration=str(sys.argv[7])
+smooth=str(sys.argv[8])
+
+import neighbor_info
+import curvature
+import metric_regression
+import roi
+import wb
+
+################################################################################
+mesh='164k'
+wb.wb_metric_resample_to_164k(subjects_dir,subject,hemi,surface,mesh) 
+wb.wb_surf_resample_to_164k(subjects_dir,subject,hemi,surface,mesh) 
+x,y,z,a,b,c = wb.wb_taubin(subjects_dir,subject,hemi,surface,mesh,iteration)
+neighbor_info.neighbor_info(a,b,c,x,y,z,subjects_dir,subject,hemi,surface,mesh)
+curvature.Gaussian_curvature(x,y,z,subjects_dir,subject,hemi,surface)
+curvature.mean_curvature(x,y,z,subjects_dir,subject,hemi,surface)
+curvature.k1_k2_SI(x,y,z,subjects_dir,subject,hemi,surface)
+wb.wb_smooth(subjects_dir,subject,hemi,surface,mesh,smooth)
+wb.wb_rois(subjects_dir,subject,hemi,surface,mesh,number)
+rois,weights=roi.roi(subjects_dir,subject,hemi,surface,mesh,number)
+metric_regression.metric_regression(subjects_dir,subject,hemi,surface,mesh,rois,weights) 
+wb.wb_metric_resample_to_native(subjects_dir,subject,hemi,surface) 
+wb.wb_structure(subjects_dir,subject,hemi,surface,structure)
+wb.wb_set_map_names(subjects_dir,subject,hemi)
+
+