1001 lines
30 KiB
Tcsh
Executable file
1001 lines
30 KiB
Tcsh
Executable file
#!/bin/tcsh -f
|
|
# fs-synthmorph-reg - sources
|
|
if(-e $FREESURFER_HOME/sources.csh) then
|
|
source $FREESURFER_HOME/sources.csh
|
|
endif
|
|
|
|
if($?FSMNI152DIR == 0) setenv FSMNI152DIR $FREESURFER/average/mni_icbm152_nlin_asym_09c
|
|
unsetenv FS_LOCAL_PYTHONPATH
|
|
|
|
set VERSION = '$Id$';
|
|
set scriptname = `basename $0`
|
|
|
|
set outdir = ();
|
|
set subject = ();
|
|
set m3z = ();
|
|
set m3zinv = ()
|
|
set invol = ()
|
|
set targvol = ()
|
|
set lta2 = ()
|
|
set vgthresh = 1e-5
|
|
set StripInput = 0
|
|
set StripTarget = 0
|
|
set ReRun = 0
|
|
|
|
set threads = 1
|
|
set ForceUpdate = 0
|
|
set antsreg = antsRegistrationSyNQuick.sh
|
|
set UseQuick = 0
|
|
set CropInput = 1
|
|
set CropTarget = 1
|
|
set MNITarget = 1
|
|
# Synthmorph uses 1mm internally, so useing target res=1mm does not
|
|
# slow things down or take more mem and it is probably gives a little
|
|
# bit better registration. ANTs speed and performance will be
|
|
# affected.
|
|
set MNITargetRes = 1.0mm
|
|
set MNIOutputRes = 1.0mm
|
|
set DoCBIG = 0
|
|
set DoTest = 0
|
|
set dim = 3
|
|
set UseAnts = 0 # 0 means use synthmorph
|
|
set ComputeInverse = 1
|
|
set PitStr = ""
|
|
set AffineOnly = 0
|
|
|
|
set tmpdir = ();
|
|
set cleanup = 1;
|
|
set LF = ();
|
|
|
|
set inputargs = ($argv);
|
|
set PrintHelp = 0;
|
|
if($#argv == 0) goto usage_exit;
|
|
set n = `echo $argv | grep -e -help | wc -l`
|
|
if($n != 0) then
|
|
set PrintHelp = 1;
|
|
goto usage_exit;
|
|
endif
|
|
set n = `echo $argv | grep -e -version | wc -l`
|
|
if($n != 0) then
|
|
echo $VERSION
|
|
exit 0;
|
|
endif
|
|
goto parse_args;
|
|
parse_args_return:
|
|
goto check_params;
|
|
check_params_return:
|
|
|
|
set StartTime = `date`;
|
|
set tSecStart = `date '+%s'`;
|
|
set year = `date +%Y`
|
|
set month = `date +%m`
|
|
set day = `date +%d`
|
|
set hour = `date +%H`
|
|
set min = `date +%M`
|
|
|
|
if($#outdir) then
|
|
mkdir -p $outdir/log
|
|
else
|
|
set outdir = `dirname $m3z`
|
|
mkdir -p $outdir
|
|
endif
|
|
pushd $outdir > /dev/null
|
|
set outdir = `pwd`;
|
|
popd > /dev/null
|
|
|
|
if($#tmpdir == 0) set tmpdir = $outdir/tmp
|
|
mkdir -p $tmpdir
|
|
|
|
# Set up log file
|
|
if($#LF == 0) set LF = $outdir/log/fs-synthmorph-reg.Y$year.M$month.D$day.H$hour.M$min.log
|
|
if($LF != /dev/null) rm -f $LF
|
|
echo "Log file for fs-synthmorph-reg" >> $LF
|
|
date | tee -a $LF
|
|
echo "" | tee -a $LF
|
|
echo "setenv SUBJECTS_DIR $SUBJECTS_DIR" | tee -a $LF
|
|
echo "cd `pwd`" | tee -a $LF
|
|
echo $0 $inputargs | tee -a $LF
|
|
ls -l $0 | tee -a $LF
|
|
echo "" | tee -a $LF
|
|
cat $FREESURFER_HOME/build-stamp.txt | tee -a $LF
|
|
echo $VERSION | tee -a $LF
|
|
uname -a | tee -a $LF
|
|
echo "pid $$" | tee -a $LF
|
|
if($?PBS_JOBID) then
|
|
echo "pbsjob $PBS_JOBID" >> $LF
|
|
endif
|
|
if($?SLURM_JOB_ID) then
|
|
echo SLURM_JOB_ID $SLURM_JOB_ID >> $LF
|
|
endif
|
|
|
|
#========================================================
|
|
# Note: Ants might not be deterministc with threads
|
|
setenv OMP_NUM_THREADS $threads
|
|
setenv ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS $threads
|
|
|
|
if($StripInput) then
|
|
#set involstripped = $tmpdir/invol.stripped.nii.gz # breaks header
|
|
set involstripped = $tmpdir/invol.stripped.mgz
|
|
set ud = `UpdateNeeded $involstripped $invol`
|
|
if($ud || $ForceUpdate) then
|
|
set cmd = (mri_synthstrip -i $invol -o $involstripped -t $threads)
|
|
echo $cmd | tee -a $LF
|
|
$cmd | tee -a $LF
|
|
if($status) exit 1
|
|
endif
|
|
set invol = $involstripped
|
|
endif
|
|
if($StripTarget) then
|
|
#set targvolstripped = $tmpdir/targ.stripped.nii.gz # breaks header
|
|
set targvolstripped = $tmpdir/targ.stripped.mgz
|
|
set ud = `UpdateNeeded $targvolstripped $targvol`
|
|
if($ud || $ForceUpdate) then
|
|
set cmd = (mri_synthstrip -i $targvol -o $targvolstripped -t $threads)
|
|
echo $cmd | tee -a $LF
|
|
$cmd | tee -a $LF
|
|
if($status) exit 1
|
|
endif
|
|
set targvol = $targvolstripped
|
|
endif
|
|
|
|
if($CropInput) then # Crop for speed
|
|
# Input volume
|
|
set involcrop = $tmpdir/invol.crop.nii.gz
|
|
set ud = `UpdateNeeded $involcrop $invol`
|
|
if($ud || $ForceUpdate) then
|
|
set cmd = (mri_mask -bb 3 $invol $invol $involcrop)
|
|
echo $cmd | tee -a $LF
|
|
$cmd | tee -a $LF
|
|
if($status) exit 1
|
|
endif
|
|
# To get back to the full FoV
|
|
set regcroptoinvol = $tmpdir/reg.crop-to-invol.lta
|
|
set ud = `UpdateNeeded $regcroptoinvol $involcrop $invol`
|
|
if($ud || $ForceUpdate) then
|
|
set cmd = (lta_convert --inlta identity.nofile --src $involcrop \
|
|
--trg $invol --outlta $regcroptoinvol)
|
|
echo $cmd | tee -a $LF
|
|
$cmd |& tee -a $LF
|
|
if($status) exit 1
|
|
endif
|
|
else
|
|
set involcrop = $invol
|
|
set regcroptoinvol = ()
|
|
endif
|
|
|
|
if($CropTarget) then
|
|
if(! $MNITarget) then
|
|
# Target vol
|
|
set targvolcrop = $tmpdir/targvol.crop.nii.gz
|
|
set ud = `UpdateNeeded $targvolcrop $targvol`
|
|
if($ud || $ForceUpdate) then
|
|
set cmd = (mri_mask -crop 3 $targvol $targvol $targvolcrop)
|
|
echo $cmd | tee -a $LF
|
|
$cmd | tee -a $LF
|
|
if($status) exit 1
|
|
endif
|
|
# To get back to the full FoV
|
|
set regcroptotarg = $tmpdir/reg.crop-to-targ.lta
|
|
set ud = `UpdateNeeded $regcroptotarg $targvolcrop $targvol`
|
|
if($ud || $ForceUpdate) then
|
|
set cmd = (lta_convert --inlta identity.nofile --src $targvolcrop \
|
|
--trg $targvol --outlta $regcroptotarg)
|
|
echo $cmd | tee -a $LF
|
|
$cmd |& tee -a $LF
|
|
if($status) exit 1
|
|
endif
|
|
else
|
|
set targvolcrop = $FSMNI152DIR/reg-targets/mni152.${MNITargetRes}$PitStr.cropped.nii.gz
|
|
set regcroptotarg = $FSMNI152DIR/reg-targets/reg.${MNITargetRes}.cropped.to.${MNIOutputRes}.lta
|
|
# Note that lta does not have PitStr
|
|
endif
|
|
else
|
|
set targvolcrop = $targvol
|
|
set regcroptotarg = ()
|
|
endif
|
|
|
|
if($UseAnts) then
|
|
# Compute both the affine and the warp with ANTs. The Affine goes from
|
|
# the input croped space to the target cropped space. The warp goes
|
|
# from target cropped to target cropped.
|
|
set warp = $tmpdir/reg.1Warp.nii.gz
|
|
set affinemat = $tmpdir/reg.0GenericAffine.mat
|
|
set ud = `UpdateNeeded $warp $involcrop $targvol`
|
|
if($ud || $ForceUpdate) then
|
|
date | tee -a $LF
|
|
pushd $tmpdir # this program can crash if stuff in the current dir
|
|
if(! $DoCBIG) then
|
|
set cmd = ($antsreg -d $dim -m $involcrop -f $targvolcrop -o reg. -n $threads)
|
|
else
|
|
set cmd = (antsRegistration --dimensionality $dim \
|
|
--output [reg.,reg.Warped.nii.gz,reg.InverseWarped.nii.gz] \
|
|
--initial-moving-transform [ $involcrop, $targvolcrop, 1 ] \
|
|
--collapse-output-transforms 1 \
|
|
--use-histogram-matching 1 \
|
|
--use-estimate-learning-rate-once 1 \
|
|
--metric mattes[ $fixed, $moving, 1, 32, regular, 0.3] \
|
|
--transform affine[ 0.1 ] \
|
|
--convergence [ 100x100x200, 1.e-8, 20 ] \
|
|
--smoothing-sigmas 4x2x1vox \
|
|
--shrink-factors 3x2x1 -l 1 \
|
|
-metric cc[ $fixed, $moving, 1, 4] \
|
|
--transform SyN[ .20, 3, 0] \
|
|
--convergence [ 100x100x50, 0, 5 ] \
|
|
--smoothing-sigmas 1x0.5x0vox \
|
|
--shrink-factors 4x2x1)
|
|
endif
|
|
echo "\n\n"| tee -a $LF
|
|
echo $cmd | tee -a $LF
|
|
fs_time $cmd |& tee -a $LF
|
|
if($status) goto error_exit
|
|
echo "\n\n"| tee -a $LF
|
|
if($CropInput || $CropTarget) then
|
|
# Somehow ANTs will change the geom slightly, so map it back to the
|
|
# the target space. This is a bit of a hack that seems to work, but
|
|
# it is not principled (ie, it was derived by trial and error, so
|
|
# there is a danger that it will not work all the time).
|
|
set cmd = (mri_vol2vol --regheader --targ $targvolcrop --mov $warp --o $warp)
|
|
echo "\n\n"| tee -a $LF
|
|
echo $cmd | tee -a $LF
|
|
fs_time $cmd |& tee -a $LF
|
|
if($status) goto error_exit
|
|
echo "\n\n"| tee -a $LF
|
|
endif
|
|
popd
|
|
endif
|
|
|
|
# Convert the ANTs affine to txt file (maps involcrop to the input of the targvolcrop space)
|
|
set affinetxt = $tmpdir/reg.0GenericAffine.txt
|
|
set ud = `UpdateNeeded $affinetxt $affinemat`
|
|
if($ud || $ForceUpdate) then
|
|
date | tee -a $LF
|
|
set cmd = (ConvertTransformFile $dim $affinemat $affinetxt --hm --ras)
|
|
echo "\n\n"| tee -a $LF
|
|
echo $cmd | tee -a $LF
|
|
fs_time $cmd |& tee -a $LF
|
|
if($status) goto error_exit
|
|
echo "\n\n"| tee -a $LF
|
|
endif
|
|
|
|
# Convert the ANTs affine to an LTA
|
|
set affinelta = $tmpdir/reg.0GenericAffine.lta
|
|
set ud = `UpdateNeeded $affinelta $affinetxt`
|
|
if($ud || $ForceUpdate) then
|
|
date | tee -a $LF
|
|
set cmd = (lta_convert --src $involcrop --trg $targvolcrop --outlta $affinelta)
|
|
if($dim == 2) set cmd = ($cmd --inniftyreg2d $affinetxt);
|
|
if($dim == 3) set cmd = ($cmd --inniftyreg $affinetxt);
|
|
echo "\n\n"| tee -a $LF
|
|
echo $cmd | tee -a $LF
|
|
fs_time $cmd |& tee -a $LF
|
|
if($status) goto error_exit
|
|
echo "\n\n"| tee -a $LF
|
|
endif
|
|
|
|
# Create one lta that maps from the uncropped input vol space to the taret space
|
|
if($CropInput || $CropTarget) then
|
|
set regcroptargtocropinvol = $tmpdir/reg.croptarg_to_cropinvol.lta
|
|
set ud = `UpdateNeeded $regcroptargtocropinvol $affinelta $regcroptoinvol`
|
|
if($ud || $ForceUpdate) then
|
|
date | tee -a $LF
|
|
set cmd = (mri_concatenate_lta -invert1 -invertout $affinelta $regcroptoinvol $regcroptargtocropinvol)
|
|
echo "\n\n"| tee -a $LF
|
|
echo $cmd | tee -a $LF
|
|
fs_time $cmd |& tee -a $LF
|
|
if($status) goto error_exit
|
|
echo "\n\n"| tee -a $LF
|
|
endif
|
|
else
|
|
set regcroptargtocropinvol = $affinelta
|
|
endif
|
|
|
|
# Concatenate the affines and the warp (ANTS)
|
|
set ud = `UpdateNeeded $m3z $warp $targvolcrop $regcroptargtocropinvol $regcroptotarg `
|
|
if($ud || $ForceUpdate) then
|
|
date | tee -a $LF
|
|
set cmd = (mri_warp_convert --initk $warp --insrcgeom $targvolcrop \
|
|
--outm3z $m3z --vg-thresh $vgthresh)
|
|
if($CropInput) set cmd = ($cmd --lta1-inv $regcroptargtocropinvol )
|
|
if($#regcroptotarg) set cmd = ($cmd --lta2 $regcroptotarg)
|
|
echo "\n\n"| tee -a $LF
|
|
echo $cmd | tee -a $LF
|
|
fs_time $cmd |& tee -a $LF
|
|
if($status) goto error_exit
|
|
echo "\n\n"| tee -a $LF
|
|
endif
|
|
|
|
else # End ANTS ================= else Use synthmorph
|
|
|
|
set gpuopt = ""
|
|
# Compute both the affine, maps from involcrop to targvolcrop
|
|
set affinelta = $outdir/aff.lta
|
|
set ud = `UpdateNeeded $affinelta $involcrop $targvolcrop`
|
|
if($ud || $ForceUpdate) then
|
|
date | tee -a $LF
|
|
set cmd = (mri_synthmorph -m affine -t $affinelta $involcrop $targvolcrop -j $threads $gpuopt)
|
|
echo "\n\n"| tee -a $LF
|
|
echo $cmd | tee -a $LF
|
|
fs_time $cmd |& tee -a $LF
|
|
if($status) goto error_exit
|
|
endif
|
|
|
|
# Create one lta that maps from the uncropped input vol space to the (cropped) target space
|
|
set regtargtoinvol = $outdir/reg.targ_to_invol.lta
|
|
set reginvoltotarg = $outdir/reg.invol_to_targ.lta
|
|
if($CropInput || $CropTarget) then
|
|
set reginvoltocroptarg = $tmpdir/reg.invol_to_croptarg.lta
|
|
set ud = `UpdateNeeded $reginvoltocroptarg $affinelta $regcroptoinvol`
|
|
if($ud || $ForceUpdate) then
|
|
date | tee -a $LF
|
|
set cmd = (mri_concatenate_lta -invert1 -invertout $affinelta $regcroptoinvol $reginvoltocroptarg)
|
|
echo "\n\n"| tee -a $LF
|
|
echo $cmd | tee -a $LF
|
|
fs_time $cmd |& tee -a $LF
|
|
if($status) goto error_exit
|
|
echo "\n\n"| tee -a $LF
|
|
endif
|
|
if($MNITarget) then
|
|
# Create one lta that maps from the uncropped input vol space to the uncropped target space
|
|
# This produces a linear tranform that can be used to map from target to input volume just
|
|
# like the nonlinear transform
|
|
set ud = `UpdateNeeded $regtargtoinvol $reginvoltocroptarg`
|
|
if($ud || $ForceUpdate) then
|
|
set regtargtocropped = $FSMNI152DIR/reg-targets/reg.$MNIOutputRes.to.$MNITargetRes.cropped.lta
|
|
set cmd = (mri_concatenate_lta -invert2 $regtargtocropped $reginvoltocroptarg $regtargtoinvol)
|
|
echo $cmd | tee -a $LF
|
|
fs_time $cmd |& tee -a $LF
|
|
if($status) goto error_exit
|
|
echo "\n\n"| tee -a $LF
|
|
endif
|
|
endif
|
|
# Compute the inverse: invol-to-targ
|
|
set ud = `UpdateNeeded $reginvoltotarg $regtargtoinvol`
|
|
if($ud || $ForceUpdate) then
|
|
set cmd = (lta_convert --invert --inlta $regtargtoinvol --outlta $reginvoltotarg)
|
|
echo $cmd | tee -a $LF
|
|
fs_time $cmd |& tee -a $LF
|
|
if($status) goto error_exit
|
|
echo "\n\n"| tee -a $LF
|
|
endif
|
|
else
|
|
pushd $tmpdir
|
|
ln -fs aff.lta reg.invol_to_targ.lta
|
|
popd
|
|
set reginvoltocroptarg = $affinelta
|
|
set reginvoltotarg = $affinelta
|
|
# Compute target to invol
|
|
set ud = `UpdateNeeded $regtargtoinvol $reginvoltotarg`
|
|
if($ud || $ForceUpdate) then
|
|
set cmd = (lta_convert --invert --inlta $reginvoltotarg --outlta $regtargtoinvol)
|
|
echo $cmd | tee -a $LF
|
|
fs_time $cmd |& tee -a $LF
|
|
if($status) goto error_exit
|
|
echo "\n\n"| tee -a $LF
|
|
endif
|
|
|
|
endif
|
|
|
|
if($MNITarget) then
|
|
# Create an xfm file which can be used as the taliarch.xfm
|
|
set xfm = $outdir/reg.invol_to_targ.xfm
|
|
set ud = `UpdateNeeded $xfm $reginvoltotarg`
|
|
if($ud) then
|
|
set cmd = (lta_convert --inlta $reginvoltotarg --outmni $xfm)
|
|
echo $cmd | tee -a $LF
|
|
$cmd |& tee -a $LF
|
|
if($status) goto error_exit
|
|
echo "\n\n"| tee -a $LF
|
|
endif
|
|
endif
|
|
|
|
echo "" | tee -a $LF
|
|
echo "To check affine registration" | tee -a $LF
|
|
echo "tkregisterfv --mov $invol --targ $targvol --reg $reginvoltotarg" | tee -a $LF
|
|
echo "" | tee -a $LF
|
|
|
|
if($AffineOnly) then
|
|
echo "AffineOnly specified, so exiting now" | tee -a $LF
|
|
goto done
|
|
endif
|
|
|
|
# Now compute the nonlinear part
|
|
set deform = $tmpdir/deform.mgz
|
|
set smvol = $tmpdir/synthmorph.out.mgz
|
|
set ud = `UpdateNeeded $deform $affinelta $involcrop $targvolcrop`
|
|
if($ud) then
|
|
set cmd = (mri_synthmorph -m deform -t $deform -i $affinelta $involcrop $targvolcrop -j $threads $gpuopt)
|
|
if($DoTest) set cmd = ($cmd -o $smvol); # can be a pain when you want to rerun to create test output
|
|
echo "\n\n" | tee -a $LF
|
|
echo $cmd | tee -a $LF
|
|
fs_time $cmd |& tee -a $LF
|
|
if($status) goto error_exit
|
|
endif
|
|
|
|
# Compute warp to the full FoV of the input (synthmorph)
|
|
set ud = `UpdateNeeded $m3z $deform $targvolcrop $reginvoltocroptarg $regcroptotarg`
|
|
if($ud || $ForceUpdate) then
|
|
# This call is not intuitive to me. It appears that the insrcgeom
|
|
# will set the vox2ras for the image-side of the warp. The warp
|
|
# maps from target voxel to cropped input RAS. The insrcgeom
|
|
# indicates how to go from that RAS to a CRS in the cropped input,
|
|
# then the lta1 indicates how to go from the cropped CRS to the
|
|
# uncropped CRS.
|
|
set cmd = (mri_warp_convert --inras $deform --insrcgeom $involcrop\
|
|
--outm3z $m3z --vg-thresh $vgthresh)
|
|
if($CropInput) set cmd = ($cmd --lta1-inv $regcroptoinvol)
|
|
if($#regcroptotarg) set cmd = ($cmd --lta2 $regcroptotarg)
|
|
echo "\n\n$cmd" | tee -a $LF
|
|
fs_time $cmd |& tee -a $LF
|
|
if($status) goto error_exit
|
|
endif
|
|
|
|
endif # Synthmorph
|
|
|
|
if($ComputeInverse) then
|
|
set ud = `UpdateNeeded $m3zinv $m3z`
|
|
if($ud || $ForceUpdate) then
|
|
set cmd = (mri_ca_register -invert-and-save $m3z $m3zinv)
|
|
echo "\n\n$cmd" | tee -a $LF
|
|
fs_time $cmd |& tee -a $LF
|
|
if($status) goto error_exit
|
|
endif
|
|
endif
|
|
|
|
# Test against output generated by the registration program
|
|
if($DoTest) then
|
|
if($UseAnts) then
|
|
set mvol = $tmpdir/reg.Warped.nii.gz
|
|
else
|
|
set mvol = $smvol
|
|
endif
|
|
if($CropInput || $CropTarget) then
|
|
# Ants will output to the cropped target space, need to resample
|
|
# to the full uncropped space.
|
|
set mvol2 = $tmpdir/morph.out.nii.gz
|
|
set ud = `UpdateNeeded $mvol2 $mvol $targvol`
|
|
if($ud || $ForceUpdate) then
|
|
set cmd = (mri_vol2vol --regheader --mov $mvol --targ $targvol --o $mvol2)
|
|
echo $cmd | tee -a $LF
|
|
fs_time $cmd |& tee -a $LF
|
|
if($status) exit 1
|
|
endif
|
|
set mvol = $mvol2
|
|
endif
|
|
|
|
set ud = `UpdateNeeded $testvol $invol $m3z $mvol`
|
|
if($ud || $ForceUpdate) then
|
|
set cmd = (mri_convert -rt nearest $invol -at $m3z $testvol)
|
|
echo $cmd | tee -a $LF
|
|
fs_time $cmd |& tee -a $LF
|
|
if($status) exit 1
|
|
set cmd = (mri_diff --po $testvol $mvol)
|
|
echo $cmd | tee -a $LF
|
|
$cmd | tee $tmpdir/test.diff.dat |& tee -a $LF
|
|
endif
|
|
echo "tkmeditfv -f $targvol -aux $testvol" |& tee -a $LF
|
|
endif
|
|
|
|
#========================================================
|
|
|
|
done:
|
|
|
|
# Cleanup
|
|
if($cleanup) rm -rf $tmpdir
|
|
|
|
# Done
|
|
echo " " |& tee -a $LF
|
|
set tSecEnd = `date '+%s'`;
|
|
@ tSecRun = $tSecEnd - $tSecStart;
|
|
set tRunMin = `echo $tSecRun/60|bc -l`
|
|
set tRunMin = `printf %5.2f $tRunMin`
|
|
set tRunHours = `echo $tSecRun/3600|bc -l`
|
|
set tRunHours = `printf %5.2f $tRunHours`
|
|
echo "Started at $StartTime " |& tee -a $LF
|
|
echo "Ended at `date`" |& tee -a $LF
|
|
echo "Fs-Synthmorph-Reg-Run-Time-Sec $tSecRun" |& tee -a $LF
|
|
echo "Fs-Synthmorph-Reg-Run-Time-Min $tRunMin" |& tee -a $LF
|
|
echo "Fs-Synthmorph-Reg-Run-Time-Hours $tRunHours" |& tee -a $LF
|
|
echo " " |& tee -a $LF
|
|
echo "fs-synthmorph-reg Done" |& tee -a $LF
|
|
exit 0
|
|
|
|
###############################################
|
|
|
|
############--------------##################
|
|
error_exit:
|
|
echo "ERROR:"
|
|
|
|
exit 1;
|
|
###############################################
|
|
|
|
############--------------##################
|
|
parse_args:
|
|
set cmdline = ($argv);
|
|
while( $#argv != 0 )
|
|
|
|
set flag = $argv[1]; shift;
|
|
|
|
switch($flag)
|
|
|
|
case "--i":
|
|
if($#argv < 1) goto arg1err;
|
|
set invol = $argv[1]; shift;
|
|
if(! -e $invol) then
|
|
echo "ERROR: cannot find $invol"
|
|
exit 1
|
|
endif
|
|
set invol = `getfullpath $invol`
|
|
breaksw
|
|
|
|
case "--t":
|
|
if($#argv < 1) goto arg1err;
|
|
set targvol = $argv[1]; shift;
|
|
if(! -e $targvol) then
|
|
echo "ERROR: cannot find $targvol"
|
|
exit 1
|
|
endif
|
|
set targvol = `getfullpath $targvol`
|
|
set MNITarget = 0
|
|
set Crop = 1
|
|
breaksw
|
|
|
|
case "--warp":
|
|
case "--m3z":
|
|
if($#argv < 1) goto arg1err;
|
|
set m3z = $argv[1]; shift;
|
|
breaksw
|
|
|
|
case "--affine-only":
|
|
set AffineOnly = 1
|
|
breaksw
|
|
|
|
case "--o":
|
|
if($#argv < 1) goto arg1err;
|
|
set outdir = $argv[1]; shift;
|
|
breaksw
|
|
|
|
case "--s":
|
|
if($#argv < 1) goto arg1err;
|
|
set subject = $argv[1]; shift;
|
|
breaksw
|
|
|
|
case "--sd":
|
|
if($#argv < 1) goto arg1err;
|
|
setenv SUBJECTS_DIR $argv[1]; shift;
|
|
breaksw
|
|
|
|
case "--threads":
|
|
if($#argv < 1) goto arg1err;
|
|
set threads = $argv[1]; shift;
|
|
breaksw
|
|
|
|
case "--vg-thresh":
|
|
if($#argv < 1) goto arg1err;
|
|
set vgthresh = $argv[1]; shift;
|
|
breaksw
|
|
|
|
case "--synthmorph":
|
|
set UseAnts = 0
|
|
breaksw
|
|
case "--ants":
|
|
set UseAnts = 1
|
|
breaksw
|
|
|
|
case "--test":
|
|
set DoTest = 1
|
|
breaksw
|
|
case "--no-test":
|
|
set DoTest = 0
|
|
breaksw
|
|
|
|
case "--2d":
|
|
set dim = 2
|
|
set DoTest = 0
|
|
breaksw
|
|
|
|
case "--quick":
|
|
set antsreg = antsRegistrationSyNQuick.sh
|
|
set UseQuick = 1
|
|
breaksw
|
|
case "--no-quick":
|
|
case "--syn":
|
|
set antsreg = antsRegistrationSyN.sh
|
|
set UseQuick = 0
|
|
breaksw
|
|
|
|
case "--crop":
|
|
set CropInput = 1
|
|
set CropTarget = 1
|
|
breaksw
|
|
case "--no-crop":
|
|
set CropInput = 0
|
|
set CropTarget = 0
|
|
breaksw
|
|
|
|
case "--strip":
|
|
set StripInput = 1
|
|
set StripTarget = 1
|
|
breaksw
|
|
case "--no-strip":
|
|
set StripInput = 0
|
|
set StripTarget = 0
|
|
breaksw
|
|
case "--strip-input":
|
|
set StripInput = 1
|
|
breaksw
|
|
case "--no-strip-input":
|
|
set StripInput = 0
|
|
breaksw
|
|
case "--strip-target":
|
|
set StripTarget = 1
|
|
breaksw
|
|
case "--no-strip-target":
|
|
set StripTarget = 0
|
|
breaksw
|
|
|
|
case "--inv":
|
|
set ComputeInverse = 1
|
|
breaksw
|
|
case "--no-inv":
|
|
set ComputeInverse = 0
|
|
breaksw
|
|
|
|
case "--mni":
|
|
set MNITarget = 1
|
|
breaksw
|
|
case "--no-nmni":
|
|
set MNITarget = 0
|
|
breaksw
|
|
|
|
case "--mni-res":
|
|
case "--mni-int-res":
|
|
case "--mni-targ-res":
|
|
if($#argv < 1) goto arg1err;
|
|
set MNITargetRes = $argv[1]; shift;
|
|
if($MNITargetRes != 1.0mm && $MNITargetRes != 1.5mm && $MNITargetRes != 2.0mm) then
|
|
echo "ERROR: --mni-res must be 1.0mm, 1.5mm, or 2.0mm"
|
|
exit 1
|
|
endif
|
|
set MNITarget = 1
|
|
breaksw
|
|
|
|
case "--mni-output-res":
|
|
case "--mni-out-res":
|
|
if($#argv < 1) goto arg1err;
|
|
set MNIOutputRes = $argv[1]; shift;
|
|
if($MNIOutputRes != 1.0mm && $MNIOutputRes != 1.5mm && $MNIOutputRes != 2.0mm && $MNIOutputRes != conformed) then
|
|
echo "ERROR: --mni-out-res must be 1.0mm, 1.5mm, 2.0mm, or conformed"
|
|
exit 1
|
|
endif
|
|
set MNITarget = 1
|
|
breaksw
|
|
|
|
case "--mni-1":
|
|
set MNIOutputRes = 1.0mm
|
|
set MNITarget = 1
|
|
breaksw
|
|
|
|
case "--mni-1.5":
|
|
set MNIOutputRes = 1.5mm
|
|
set MNITarget = 1
|
|
breaksw
|
|
|
|
case "--mni-2":
|
|
set MNIOutputRes = 2.0mm
|
|
set MNITarget = 1
|
|
breaksw
|
|
|
|
case "--mni-conformed":
|
|
set MNIOutputRes = conformed
|
|
set MNITarget = 1
|
|
breaksw
|
|
|
|
case "--pituitary":
|
|
set PitStr = ".pit"
|
|
breaksw
|
|
|
|
case "--cbig":
|
|
set DoCBig = 1
|
|
breaksw
|
|
case "--no-cbig":
|
|
set DoCBig = 0
|
|
breaksw
|
|
|
|
case "--force":
|
|
set ForceUpdate = 1
|
|
breaksw
|
|
case "--rerun":
|
|
# Prevents it from bailing out on the first check. This can be
|
|
# handy when debugging and you want it to get into the interior
|
|
# of the code.
|
|
set ReRun = 1
|
|
breaksw
|
|
|
|
case "--log":
|
|
if($#argv < 1) goto arg1err;
|
|
set LF = $argv[1]; shift;
|
|
breaksw
|
|
|
|
case "--nolog":
|
|
case "--no-log":
|
|
set LF = /dev/null
|
|
breaksw
|
|
|
|
case "--tmp":
|
|
case "--tmpdir":
|
|
if($#argv < 1) goto arg1err;
|
|
set tmpdir = $argv[1]; shift;
|
|
set cleanup = 0;
|
|
breaksw
|
|
|
|
case "--nocleanup":
|
|
set cleanup = 0;
|
|
breaksw
|
|
|
|
case "--cleanup":
|
|
set cleanup = 1;
|
|
breaksw
|
|
|
|
case "--debug":
|
|
set verbose = 1;
|
|
set echo = 1;
|
|
breaksw
|
|
|
|
default:
|
|
echo ERROR: Flag $flag unrecognized.
|
|
echo $cmdline
|
|
exit 1
|
|
breaksw
|
|
endsw
|
|
|
|
end
|
|
|
|
goto parse_args_return;
|
|
############--------------##################
|
|
|
|
############--------------##################
|
|
check_params:
|
|
|
|
if($UseAnts == -1) then
|
|
echo "ERROR: must spec either --ants or --synthmorph"
|
|
exit 1
|
|
endif
|
|
|
|
# This will keep the temp stuff from being deleted
|
|
if($#outdir && $#tmpdir == 0) then
|
|
set tmpdir = $outdir
|
|
set cleanup = 0
|
|
endif
|
|
|
|
if($#subject) then
|
|
set sd = $SUBJECTS_DIR/$subject
|
|
if(! -e $sd) then
|
|
echo "ERROR: cannot find $subject"
|
|
exit 1;
|
|
endif
|
|
if($#invol == 0) set invol = $sd/mri/norm.mgz
|
|
if($#outdir == 0 && $MNITarget) then
|
|
if($UseAnts) then
|
|
set outdir = $sd/mri/transforms/ants.warp.$MNITargetRes.$MNIOutputRes$PitStr
|
|
else
|
|
set outdir = $sd/mri/transforms/synthmorph.$MNITargetRes.$MNIOutputRes$PitStr
|
|
endif
|
|
#don't create particular files in the transform dir
|
|
#set m3z = $outdir.nii.gz
|
|
#set m3zinv = $outdir.inv.nii.gz # Just in case the invert is requested
|
|
endif
|
|
endif
|
|
|
|
# This is too complicated now as it is not just setting the target volume to
|
|
# the cropped; have to set the reg-crop-to-target too
|
|
#if($MNITarget) set CropTarget = 0;
|
|
|
|
if(! $MNITarget) then
|
|
# Could get this to work, but not really needed
|
|
set CropTarget = 0;
|
|
set CropInput = 0;
|
|
endif
|
|
|
|
if($#outdir == 0) then
|
|
echo "ERROR: must spec output dir"
|
|
exit 1
|
|
endif
|
|
if($#invol == 0) then
|
|
echo "ERROR: must spec input volume"
|
|
exit 1
|
|
endif
|
|
foreach f ($invol $targvol)
|
|
if(! -e $f) then
|
|
echo "ERROR: cannot find $f"
|
|
exit 1
|
|
endif
|
|
end
|
|
|
|
mkdir -p $outdir
|
|
set outdir = `getfullpath $outdir`
|
|
if($#m3z == 0) then
|
|
if($MNITarget) then
|
|
# Lose ants vs synthseg distinction in the output
|
|
set m3z = $outdir/warp.to.mni152.${MNITargetRes}.${MNIOutputRes}.nii.gz
|
|
set m3zinv = $outdir/warp.to.mni152.${MNITargetRes}.${MNIOutputRes}.inv.nii.gz
|
|
else
|
|
set m3z = $outdir/warp.nii.gz
|
|
set m3zinv = $outdir/warp.inv.nii.gz
|
|
endif
|
|
endif
|
|
|
|
if($ComputeInverse && $#m3zinv == 0) then
|
|
set stem = `fname2stem $m3z`
|
|
set ext = `fname2ext $m3z`
|
|
set m3zinv = $stem.inv.$ext
|
|
endif
|
|
|
|
if($#targvol && $MNITarget) then
|
|
echo "ERROR: cannot specify both target volume and mni target"
|
|
exit 1
|
|
endif
|
|
|
|
# MNI is already skull stripped
|
|
if($MNITarget) set targvol = $FSMNI152DIR/reg-targets/mni152.${MNIOutputRes}${PitStr}.nii.gz
|
|
|
|
if($UseQuick) then
|
|
set antsreg = antsRegistrationSyNQuick.sh
|
|
else
|
|
set antsreg = antsRegistrationSyN.sh
|
|
endif
|
|
|
|
set testvol = ()
|
|
if($DoTest) set testvol = $outdir/test.nii.gz
|
|
|
|
set ud1 = `UpdateNeeded $m3z $targvol $invol`
|
|
set ud2 = 0
|
|
if($ComputeInverse) set ud2 = `UpdateNeeded $m3zinv $targvol $invol`
|
|
if($ReRun == 0 && $ud1 == 0 && $ud2 == 0 && $ForceUpdate == 0) then
|
|
echo "fs-synthmorph-reg: update not needed"
|
|
exit 0
|
|
endif
|
|
|
|
goto check_params_return;
|
|
############--------------##################
|
|
|
|
############--------------##################
|
|
arg1err:
|
|
echo "ERROR: flag $flag requires one argument"
|
|
exit 1
|
|
############--------------##################
|
|
arg2err:
|
|
echo "ERROR: flag $flag requires two arguments"
|
|
exit 1
|
|
############--------------##################
|
|
|
|
############--------------##################
|
|
usage_exit:
|
|
echo ""
|
|
echo "fs-synthmorph-reg -- frontend for running mri_synthmorph (24GB)"
|
|
echo " --i invol"
|
|
echo " --t targetvol"
|
|
echo " --warp warp.{m3z,mgz,nii.gz} : output warp file (or save in outdir)"
|
|
echo " --o outdir"
|
|
echo " --s subject (invol=norm.mgz, warp=mri/transforms/synthmorph.1.0mm.1.0mm.nii.gz targ=MNI152.1mm)"
|
|
echo " --threads threads"
|
|
echo " --mni-targ-res 1.0mm, 1.5mm, 2.0mm (default is 1.0mm)"
|
|
echo " --mni-out-res 1.0mm, 1.5mm, 2.0mm, conformed (default is 1.0mm)"
|
|
echo " --no-inv : do not compute warp inverse (computed by default)"
|
|
echo " --affine-only : only perform the affine registration part"
|
|
echo " --strip : skull strip input and target"
|
|
echo " --strip-input : skull strip input"
|
|
echo " --strip-target : skull strip target"
|
|
echo " --crop/--no-crop (default is $CropInput, but no-crop if not registering to mni)"
|
|
echo " --pituitary : select the MNI target volume that does not mask out the pituitary"
|
|
echo " --test : resample the input volume to the output space and store as test.nii.gz"
|
|
echo " --vg-thresh vgthresh : threshold for testing diffs in volume geom"
|
|
echo " --ants : use ANTs registration intstead of synthmorph"
|
|
echo " --force : regenerate all output "
|
|
echo " --rerun : for debugging only, when you want it to get into the inteiror of the code"
|
|
echo " --tmp tmpdir"
|
|
echo " "
|
|
|
|
|
|
if(! $PrintHelp) exit 1;
|
|
echo $VERSION
|
|
cat $0 | awk 'BEGIN{prt=0}{if(prt) print $0; if($1 == "BEGINHELP") prt = 1 }'
|
|
exit 1;
|
|
|
|
#---- Everything below here is printed out as part of help -----#
|
|
BEGINHELP
|
|
|
|
This is a frontend for running mri_synthmorph, especially for
|
|
registering to mni152 space. For typical usage, it will consume less
|
|
than 15GB of memory. For affine only it requires less than 7GB.
|
|
|
|
Simple usage:
|
|
|
|
fs-synthmorph-reg --i inputvol.nii.gz --o synthmorphdir
|
|
|
|
This will register to the mni152 saving output files in synthmorphdir. The nonlinear
|
|
warp will be outdir/warp.to.mni152.1.0mm.1.0mm.nii.gz which can be applied like:
|
|
|
|
mri_convert inputvol.nii.gz -at warp.to.mni152.1.0mm.1.0mm.nii.gz inputvol.mni152.nii.gz
|
|
|
|
which can be checked with
|
|
|
|
freeview inputvol.mni152.nii.gz $FREESURFER/average/mni_icbm152_nlin_asym_09c/reg-targets/mni152.1.0mm.nii.gz
|
|
|
|
If you have created the warp from a FreeSurfer subject space (either with --s or passing a
|
|
conformed volume), then you can apply the warp to the surfaces, eg,
|
|
|
|
mris_apply_reg --warp lh.white warp.to.mni152.1.0mm.1.0mm.inv.nii.gz lh.warp.mni152
|
|
|
|
which you can check with
|
|
|
|
freeview inputvol.mni152.nii.gz $FREESURFER/average/mni_icbm152_nlin_asym_09c/reg-targets/mni152.1.0mm.nii.gz -f lh.warp.mni152
|
|
|
|
Skull stripping: mri_synthmorph is quite robust to the presence or
|
|
absence of skull stripping, but there will be some differences, so, to
|
|
be on the safe side, we recommend that images be skull stripped prior
|
|
to registration. When registering to the mni152, it will automatically
|
|
use a skull stripped volume. If you want fs-synthmorph-reg to skull
|
|
strip for you, then add --strip (uses mri_synthstrip). Note that if
|
|
you pass a FS subject (--s), it will use the norm.mgz volume, which is
|
|
already skull stripped. Note that ANTs performance and speed will be
|
|
very sensitive to stripping.
|
|
|
|
Cropping: internally, fs-synthmorph-reg will crop the input to a
|
|
minimal window around non-zero voxels in the image. In the end, this
|
|
probably does not have much of an effect for synthmorph as it will
|
|
sample to 256^3 internally. For ANTs, this can speed up the program
|
|
dramatically and reduce its memory footprint. fs-synthmorph-reg will
|
|
manage all of the transforms so that they all reference the uncropped
|
|
space, ie, the user is completely insulated from the cropping process
|
|
and its implications. To turn off cropping, use --no-crop. If you are
|
|
not stripping, it probably does not make sense to crop.
|
|
|
|
Target and Output MNI Resolution: the first 1.0mm in the warp file
|
|
name means that the registration is done to a version of the mni152
|
|
with an isotropic resolution of 1.0mm (this will be how finely the
|
|
warp field is sampled). The second 1.0mm means that the warp file will
|
|
actually warp to the the mni152 with an isotropic resolution of 1.0mm.
|
|
You can change the resolution of the warp field with the
|
|
--mni-targ-res flag (options are 1.0mm, 1.5mm, 2.0mm), but,
|
|
internally, synthmorph uses 1mm, so this makes sense. You can change
|
|
the resolution of the output in several ways. From this program, you
|
|
can specify the output resolution with --mni-out-res (options are
|
|
1.0mm, 1.5mm, 2.0mm). You can also create an LTA file to the desired
|
|
output resolution and then use mri_warp_convert to change the warp.
|
|
Or you can create an LTA file to the desired resolution, and then use
|
|
mri_vol2vol --gcam to apply the warp and the LTA (instead of
|
|
mri_convert, which will only the the warp).
|
|
|
|
Affine and affine-only: by default, this script will perform a
|
|
non-linear registration to the target space. The first step of this
|
|
will be to perform an affine registration. Sometimes, this is useful
|
|
too. This affine registration is stored in reg.targ_to_invol.lta. If
|
|
you ONLY want the affine registration, then you can specify
|
|
--affine-only.
|
|
|
|
Arbitrary Targets: you can use fs-synthmorph-reg to register to target
|
|
volumes other than the mni152. To do this, simply specify --t
|
|
targetvol. You can have fs-synthmorph-reg skull strip the target by
|
|
adding --strip-target. By default, the target will be cropped unless
|
|
--no-crop (see Cropping above).
|
|
|
|
Please cite SynthMorph:
|
|
Anatomy-specific acquisition-agnostic affine registration learned from
|
|
fictitious images
|
|
Hoffmann M, Hoopes A, Fischl B*, Dalca AV* (*equal contribution)
|
|
SPIE Medical Imaging: Image Processing, 12464, p 1246402, 2023
|
|
https://doi.org/10.1117/12.2653251
|
|
https://synthmorph.io/#papers (PDF)
|
|
|
|
SynthMorph: learning contrast-invariant registration without acquired images
|
|
Hoffmann M, Billot B, Greve DN, Iglesias JE, Fischl B, Dalca AV
|
|
IEEE Transactions on Medical Imaging, 41 (3), 543-558, 2022
|
|
https://doi.org/10.1109/TMI.2021.3116879
|
|
|
|
If you use SynthStrip in your analysis, please cite:
|
|
SynthStrip: Skull-Stripping for Any Brain Image
|
|
A Hoopes, JS Mora, AV Dalca, B Fischl, M Hoffmann
|
|
NeuroImage 206 (2022), 119474
|
|
https://doi.org/10.1016/j.neuroimage.2022.119474
|
|
Website: https://synthstrip.io
|
|
|