Skip to content

Create a nifti of within subject change (2 timepoints, output from dbm)

lizhep edited this page Aug 5, 2024 · 9 revisions

Goal: Create a nifti that represents the difference between two scans (usually different timepoints). One reason you may want to do this if you want to use a multivariate technique such as non-negative matrix factorization to examine variation and don't have a longitudinal approach, but want to capture change between two time points. To create a single image per subject representing this change, the following script warps two images together and creates a single output per subject.

For this, I used files on Niagara. Create a script in the directory where you ran the 2level, for example named within_subject_difference.sh. The files used here are outputs from Gabe's two-level dbm.

Currently, the parts of the script below only relevant for absolute jacobians are commented out.

Lines to edit to make the script run for you:

  • path to the 2level run on niagara
  • the time point files

(note: here the script is using the first two time points, however you can use any two timepoints if you have more than 2)

  • time point 1 = ./output/$subject/subject*lsq601Warp.nii.gz
  • time point 2 = ./output/$subject/subject*lsq611Warp.nii.gz
  • time point 3 = ./output/$subject/subject*lsq621Warp.nii.gz
  • time point n = ./output/$subject/subject*lsq6n1Warp.nii.gz

this is useful is you didn't explicitly hard code time point in your mri file name - however this workaround will not work for the section of the code that uses the absolute jacobians - you'll likely need to rename the mri files.

  • if you do change the time points that you will use, you should change everything that has "1_to_2" and "2_to_1" with the appropriate numbers
module load cobralab

for file in /home/m/mchakrav/lanicupo/scratch/dbm_20200807/output/subject*; do subject=$(basename $file);

if [ -s ./output/$subject/subject*lsq611Warp.nii.gz ]; then
#First create jacobian between the two timepoints
        #Merge two warps to create one between the timepoints
	antsApplyTransforms --verbose -d 3 -r ./output/$subject/*_template0.nii.gz \
	-t ./output/$subject/subject*lsq601InverseWarp.nii.gz \
	-t ./output/$subject/subject*lsq611Warp.nii.gz \
	-o [ warp_2_to_1.nii.gz,1 ] --verbose
	mv warp_2_to_1.nii.gz ./output/$subject/
	echo "Moving warp_2_to_1.nii.gz to subject folder"

#	#Create a jacobian from the merged warps
	echo $subject
	CreateJacobianDeterminantImage 3 ./output/$subject/warp_2_to_1.nii.gz \
	./output/$subject/jacobian_1_to_2.nii.gz 1 1

#Second create jacobian of either nlin+affine (absolutes) or just nlin (relative)
        #FOR ABSOLUTES
        #subtract the affines of time one and two
        #echo ./output/jacobians/groupwise/$subject*_mri_2_lsq6_affine.nii.gz
        #ImageMath 3 ./output/$subject/affine_jacobian_1_to_2.nii.gz - \
        #./output/jacobians/groupwise/$subject*_mri_2_lsq6_affine.nii.gz \
        #./output/jacobians/groupwise/$subject*_mri_1_lsq6_affine.nii.gz


        #FOR RELATIVES
        #Get ants transform from deformation field
        ANTSUseDeformationFieldToGetAffineTransform ./output/$subject/warp_2_to_1.nii.gz 0.25 \
                affine ./output/$subject/delin_subject.mat \
                ./output/$subject/*_otsumask.nii.gz

        #Create composite warp
        #output from step above is transform
        antsApplyTransforms -d 3 --verbose \
                -t [./output/$subject/delin_subject.mat,1] \
                -r ./output/$subject/${subject}_template0.nii.gz \
                -o [./output/$subject/delin_subject.nii.gz,1]

        #Create a jacobian determinant of composite warp (output from step above)
        CreateJacobianDeterminantImage 3 ./output/$subject/delin_subject.nii.gz \
                ./output/$subject/delin_jacobian.nii.gz 1 1

#Third combine the two jacobians created above
        #For the ABSOLUTES, combine the jacobians and the affines
      #combine the affines between time 1 and 2 and the jacobians between 1 and 2
        #ImageMath 3 ./output/$subject/absolute_jacobian_1_to_2.nii.gz + ./output/$subject/jacobian_1_to_2.nii.gz ./output/$subject/affine_jacobian_1_to_2.nii.gz

        #for the RELATIVES, combine the jacobians and the nlin
        ImageMath 3 ./output/$subject/relative_jacobian_1_to_2.nii.gz + ./output/$subject/jacobian_1_to_2.nii.gz ./output/$subject/delin_jacobian.nii.gz

#Fourth put them in common space
        antsApplyTransforms -d 3 -i ./output/$subject/relative_jacobian_1_to_2.nii.gz \
        -t ./output/secondlevel/secondlevel_*${subject}*_template*1Warp.nii.gz \
        -t ./output/secondlevel/secondlevel_*${subject}*_template*GenericAffine.mat \
        -o ./output/$subject/relative_1_to_2_jacobian_in_common_space.nii.gz -r ./output/secondlevel/secondlevel_template0.nii.gz

#Fifth blur them
        SmoothImage 3 ./output/$subject/relative_1_to_2_jacobian_in_common_space.nii.gz 0.08492569002 ./output/$subject/relative_common_smooth.nii.gz 1 0
        else
	echo "$subject has only one scan"

        fi
done

Note: This is currently only available for relative jacobians. Create a script for example named within_subject_difference.sh and copy the following:

#!/usr/bin/env bash
set -euo pipefail


## The script should be run as within-subchange.sh [path to dbm output][path to dbm input txt] [path to output for within subject change niftis]

# Name of output directory from DBM 
output_dbm=$1
# Name of output directory for nifti with within subject change  
output_within_sub=$3

# Create directory if it does not exist
if [[ -d ${output_within_sub} ]]
then
    echo "${output_within_sub} exists on your filesystem."
else 
    echo "${output_within_sub} does not exists on your filesystem. Creating it."
    mkdir ${output_within_sub}
fi

# For smoothing calculations 
calc () { awk "BEGIN{ print $* }" ;}
fwhm=4vox #Use 4vox as in DBM pipeline
sigma_num=$(calc "$(echo ${fwhm} | grep -o -E '^[0-9]+')/(2*sqrt(2*log(2)))")


# Loop through subjects in input txt 
i=1
while read -r subject_scans; do
    IFS=',' read -r -a scans <<<${subject_scans}
    if [[ ${#scans[@]} -gt 1 ]] && [[ ${#scans[@]} -lt 3 ]]; then
        echo 'Subject has two scans, continue'
        # Verify if directory exists:
        if [[ -d ${output_within_sub}/subject_${i} ]]; then
                echo "${output_within_sub}/subject_${i} exists on your filesystem."
                if [ "$(ls -A ${output_within_sub}/subject_${i})" ]; then
                        echo "${output_within_sub}/subject_${i} is not empty, deleting its contents and creating new one."
                        rm -r ${output_within_sub}/subject_${i}
                        mkdir ${output_within_sub}/subject_${i}
                fi
        else 
                echo "${output_within_sub}/subject_${i} does not exists on your filesystem. Creating it."
                mkdir ${output_within_sub}/subject_${i}
        fi
        scan_name=$(basename -s .nii ${scans[0]})
#First create jacobian between the two timepoints
        antsApplyTransforms --verbose -d 3 -r ${output_dbm}/firstlevel/subject_${i}/final/average/template_sharpen_shapeupdate.nii.gz \
        -t ${output_dbm}/firstlevel/subject_${i}/final/transforms/$(basename -s .nii ${scans[0]})_1InverseWarp.nii.gz \
        -t ${output_dbm}/firstlevel/subject_${i}/final/transforms/$(basename -s .nii ${scans[1]})_1Warp.nii.gz \
        -o [ ${output_within_sub}/subject_${i}/warp_2_to_1.nii.gz,1 ] --verbose
        echo "warp_2_to_1.nii.gz in subject folder"
        CreateJacobianDeterminantImage 3 ${output_within_sub}/subject_${i}/warp_2_to_1.nii.gz \
        ${output_within_sub}/subject_${i}/jacobian_1_to_2.nii.gz 1 1
#Second create jacobian of either nlin+affine (absolutes) or just nlin (relative)
#FOR RELATIVES
        #Get ants transform from deformation field
        ANTSUseDeformationFieldToGetAffineTransform ${output_within_sub}/subject_${i}/warp_2_to_1.nii.gz 0.25 \
                affine ${output_within_sub}/subject_${i}/delin_subject.mat \
                ${output_dbm}/firstlevel/subject_${i}/final/average/mask_shapeupdate.nii.gz

        #Create composite warp
        #output from step above is transform
        antsApplyTransforms -d 3 --verbose \
                -t [${output_within_sub}/subject_${i}/delin_subject.mat,1] \
                -r ${output_dbm}/firstlevel/subject_${i}/final/average/template_sharpen_shapeupdate.nii.gz \
                -o [${output_within_sub}/subject_${i}/delin_subject.nii.gz,1]

        #Create a jacobian determinant of composite warp (output from step above)
        CreateJacobianDeterminantImage 3 ${output_within_sub}/subject_${i}/delin_subject.nii.gz \
                ${output_within_sub}/subject_${i}/delin_jacobian.nii.gz 1 1

#Third combine the two jacobians created above
        #For the RELATIVES, combine the jacobians and the nlin
        ImageMath 3 ${output_within_sub}/subject_${i}/relative_jacobian_1_to_2.nii.gz + ${output_within_sub}/subject_${i}/jacobian_1_to_2.nii.gz ${output_within_sub}/subject_${i}/delin_jacobian.nii.gz

#Fourth put them in common space
        antsApplyTransforms -d 3 -i ${output_within_sub}/subject_${i}/relative_jacobian_1_to_2.nii.gz \
        -t ${output_dbm}/secondlevel/final/transforms/subject_${i}_*1Warp.nii.gz \
        -t ${output_dbm}/secondlevel/final/transforms/subject_${i}_*GenericAffine.mat \
        -o ${output_within_sub}/subject_${i}/relative_1_to_2_jacobian_in_common_space.nii.gz \
        -r ${output_dbm}/secondlevel/final/average/template_sharpen_shapeupdate.nii.gz

# Fifth put them in target space
        if [[ -d ${output_dbm}/secondlevel/final-target ]] && [[ -z "$(ls ${output_dbm}/secondlevel/final-target)" ]]; then
                echo 'No final-target transforms, not registering output to final-target'
        else
                antsApplyTransforms -d 3 --verbose -n Linear \
                -i ${output_within_sub}/subject_${i}/relative_1_to_2_jacobian_in_common_space.nii.gz \
                -o ${output_within_sub}/subject_${i}/relative_1_to_2_jacobian_in_final_target.nii.gz \
                -t ${output_dbm}/secondlevel/final-target/to_target_1Warp.nii.gz \
                -t ${output_dbm}/secondlevel/final-target/to_target_0GenericAffine.mat \
                -r ${output_dbm}/secondlevel/final-target/final_target.nii.gz    
        fi
#Fifth blur them
        SmoothImage 3 ${output_within_sub}/subject_${i}/relative_1_to_2_jacobian_in_common_space.nii.gz ${sigma_num} ${output_within_sub}/subject_${i}/${scan_name%%_*}_relative_common_smooth.nii.gz 0 0
# If fifth section was uncommented, then uncomment this too:
        SmoothImage 3 ${output_within_sub}/subject_${i}/relative_1_to_2_jacobian_in_final_target.nii.gz ${sigma_num} ${output_within_sub}/subject_${i}/${scan_name%%_*}_relative_final-target_smooth.nii.gz 0 0
    else
        echo 'Subject has one scan, do not include'
    fi
    ((++i))
done < $2

Then you can run the script like this:

./within-subchange.sh [path to dbm output] [path to dbm input txt] [path to output for within subject change niftis]
Clone this wiki locally