Add assessment of spatial coverage within both fMRIVolume and fMRISurface #138
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -23,8 +23,8 @@ Usage() { | |
| echo " --motionmatdir=<input motion correcton matrix directory>" | ||
| echo " --motionmatprefix=<input motion correcton matrix filename prefix>" | ||
| echo " --ofmri=<input fMRI 4D image>" | ||
| echo " --freesurferbrainmask=<input FreeSurfer brain mask, nifti format in atlas (MNI152) space>" | ||
|
There was a problem hiding this comment. I don't know if this is correct. There was a problem hiding this comment. Here are the inputs to --freesurferbrainmask, --biasfield (and --t1) are all from There was a problem hiding this comment. I was not claiming it was wrong, just that I didn't know, and highlighting it so someone who did know would double-check. |
||
| echo " --biasfield=<input biasfield image, in atlas (MNI152) space>" | ||
| echo " --gdfield=<input warpfield for gradient non-linearity correction>" | ||
| echo " --scoutin=<input scout image (EPI pre-sat, before gradient non-linearity distortion correction)>" | ||
| echo " --scoutgdcin=<input scout gradient nonlinearity distortion corrected image (EPI pre-sat)>" | ||
|
|
@@ -139,7 +139,7 @@ ${FSLDIR}/bin/fslmaths ${WD}/${BiasFieldFile}.${FinalfMRIResolution} -thr 0.1 ${ | |
| # NB: warpfield resolution is 10mm, so 1mm to fMRIres downsample loses no precision | ||
| ${FSLDIR}/bin/convertwarp --relout --rel --warp1=${fMRIToStructuralInput} --warp2=${StructuralToStandard} --ref=${WD}/${T1wImageFile}.${FinalfMRIResolution} --out=${OutputTransform} | ||
|
|
||
| # Add stuff for estimating RMS motion | ||
| invwarp -w ${OutputTransform} -o ${OutputInvTransform} -r ${ScoutInputgdc} | ||
| applywarp --rel --interp=nn -i ${FreeSurferBrainMask}.nii.gz -r ${ScoutInputgdc} -w ${OutputInvTransform} -o ${ScoutInputgdc}_mask.nii.gz | ||
| if [ -e ${fMRIFolder}/Movement_RelativeRMS.txt ] ; then | ||
|
|
@@ -154,46 +154,58 @@ fi | |
| if [ -e ${fMRIFolder}/Movement_AbsoluteRMS_mean.txt ] ; then | ||
| /bin/rm ${fMRIFolder}/Movement_AbsoluteRMS_mean.txt | ||
| fi | ||
|
|
||
| # Copy some files from ${WD} (--workingdir arg) to ${fMRIFolder} (--fmrifolder arg) | ||
| ${FSLDIR}/bin/imcp ${WD}/${T1wImageFile}.${FinalfMRIResolution} ${fMRIFolder}/${T1wImageFile}.${FinalfMRIResolution} | ||
| ${FSLDIR}/bin/imcp ${WD}/${FreeSurferBrainMaskFile}.${FinalfMRIResolution} ${fMRIFolder}/${FreeSurferBrainMaskFile}.${FinalfMRIResolution} | ||
| ${FSLDIR}/bin/imcp ${WD}/${BiasFieldFile}.${FinalfMRIResolution} ${fMRIFolder}/${BiasFieldFile}.${FinalfMRIResolution} | ||
|
|
||
| mkdir -p ${WD}/prevols | ||
| mkdir -p ${WD}/postvols | ||
|
|
||
| # Apply combined transformations to fMRI in a one-step resampling | ||
| # (combines gradient non-linearity distortion, motion correction, and registration to atlas (MNI152) space, but keeping fMRI resolution) | ||
| ${FSLDIR}/bin/fslsplit ${InputfMRI} ${WD}/prevols/vol -t | ||
| FrameMergeSTRING="" | ||
| FrameMergeSTRINGII="" | ||
| for ((k=0; k < $NumFrames; k++)); do | ||
| vnum=`${FSLDIR}/bin/zeropad $k 4` | ||
|
|
||
| # Add stuff for estimating RMS motion | ||
| rmsdiff ${MotionMatrixFolder}/${MotionMatrixPrefix}${vnum} ${MotionMatrixFolder}/${MotionMatrixPrefix}0000 ${ScoutInputgdc} ${ScoutInputgdc}_mask.nii.gz | tail -n 1 >> ${fMRIFolder}/Movement_AbsoluteRMS.txt | ||
| if [ $k -eq 0 ] ; then | ||
| echo "0" >> ${fMRIFolder}/Movement_RelativeRMS.txt | ||
| else | ||
| rmsdiff ${MotionMatrixFolder}/${MotionMatrixPrefix}${vnum} $prevmatrix ${ScoutInputgdc} ${ScoutInputgdc}_mask.nii.gz | tail -n 1 >> ${fMRIFolder}/Movement_RelativeRMS.txt | ||
| fi | ||
| prevmatrix="${MotionMatrixFolder}/${MotionMatrixPrefix}${vnum}" | ||
|
|
||
| # Combine GCD with motion correction | ||
| ${FSLDIR}/bin/convertwarp --relout --rel --ref=${WD}/prevols/vol${vnum}.nii.gz --warp1=${GradientDistortionField} --postmat=${MotionMatrixFolder}/${MotionMatrixPrefix}${vnum} --out=${MotionMatrixFolder}/${MotionMatrixPrefix}${vnum}_gdc_warp.nii.gz | ||
| # Add in the warp to MNI152 | ||
| ${FSLDIR}/bin/convertwarp --relout --rel --ref=${WD}/${T1wImageFile}.${FinalfMRIResolution} --warp1=${MotionMatrixFolder}/${MotionMatrixPrefix}${vnum}_gdc_warp.nii.gz --warp2=${OutputTransform} --out=${MotionMatrixFolder}/${MotionMatrixPrefix}${vnum}_all_warp.nii.gz | ||
| # Apply one-step warp, using spline interpolation | ||
| ${FSLDIR}/bin/applywarp --rel --interp=spline --in=${WD}/prevols/vol${vnum}.nii.gz --warp=${MotionMatrixFolder}/${MotionMatrixPrefix}${vnum}_all_warp.nii.gz --ref=${WD}/${T1wImageFile}.${FinalfMRIResolution} --out=${WD}/postvols/vol${vnum}.nii.gz | ||
|
|
||
| # Generate a mask for keeping track of spatial coverage (use nearest neighbor interpolation here) | ||
| ${FSLDIR}/bin/fslmaths ${WD}/prevols/vol${vnum}.nii.gz -mul 0 -add 1 ${WD}/prevols/vol${vnum}_mask.nii.gz | ||
| ${FSLDIR}/bin/applywarp --rel --interp=nn --in=${WD}/prevols/vol${vnum}_mask.nii.gz --warp=${MotionMatrixFolder}/${MotionMatrixPrefix}${vnum}_all_warp.nii.gz --ref=${WD}/${T1wImageFile}.${FinalfMRIResolution} --out=${WD}/postvols/vol${vnum}_mask.nii.gz | ||
|
|
||
| # Create strings for merging | ||
| FrameMergeSTRING+="${WD}/postvols/vol${vnum}.nii.gz " | ||
| FrameMergeSTRINGII+="${WD}/postvols/vol${vnum}_mask.nii.gz " | ||
|
|
||
| #Do Basic Cleanup | ||
| rm ${MotionMatrixFolder}/${MotionMatrixPrefix}${vnum}_gdc_warp.nii.gz | ||
| rm ${MotionMatrixFolder}/${MotionMatrixPrefix}${vnum}_all_warp.nii.gz | ||
| done | ||
|
|
||
| # Merge together results and restore the TR (saved beforehand) | ||
| ${FSLDIR}/bin/fslmerge -tr ${OutputfMRI} $FrameMergeSTRING $TR_vol | ||
| ${FSLDIR}/bin/fslmerge -tr ${OutputfMRI}_mask $FrameMergeSTRINGII $TR_vol | ||
|
|
||
| # Generate a spatial coverage mask that captures the voxels that have data available at *ALL* time points | ||
| # (gets applied in IntensityNormalization.sh; so don't change name here without changing it in that script as well). | ||
| fslmaths ${OutputfMRI}_mask -Tmin ${OutputfMRI}_mask | ||
|
|
||
| # Combine transformations: gradient non-linearity distortion + fMRI_dc to standard | ||
|
|
@@ -214,12 +226,11 @@ ${FSLDIR}/bin/fslmaths ${WD}/gdc_dc_jacobian -Tmean ${WD}/gdc_dc_jacobian | |
| #and resample it to MNI space | ||
| ${FSLDIR}/bin/applywarp --rel --interp=spline -i ${WD}/gdc_dc_jacobian -r ${WD}/${T1wImageFile}.${FinalfMRIResolution} -w ${StructuralToStandard} -o ${JacobianOut} | ||
|
|
||
| # Compute average motion across frames | ||
| cat ${fMRIFolder}/Movement_RelativeRMS.txt | awk '{ sum += $1} END { print sum / NR }' >> ${fMRIFolder}/Movement_RelativeRMS_mean.txt | ||
| cat ${fMRIFolder}/Movement_AbsoluteRMS.txt | awk '{ sum += $1} END { print sum / NR }' >> ${fMRIFolder}/Movement_AbsoluteRMS_mean.txt | ||
|
|
||
| # Do Basic Cleanup | ||
| rm -r ${WD}/postvols | ||
| rm -r ${WD}/prevols | ||
|
|
||
|
|
||
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm not particularly fond of having an internal aspect of the data determine whether an output file gets generated. If we added some automated fix for the zeros (more dilation), having the mask of zeros exist regardless of whether it is empty might make the code simpler.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I did that here because due to the dilation, zeros in the CIFTI should be rare and the *_nonzero.stats.txt is written out regardless, so that provides the "evidence" that the check was performed on a particular set of processed data. That said, if you and @glasserm prefer that we write out the *_nonzero.dscalar.nii file regardless, I can do that.
I don't see how adding more dilation impacts whether or not we should write out the map of nonzero grayordinates. At some point, we can't just keep dilating, and frankly it seems like we already have a lot of dilation built in -- if I'm reading
RibbonVolumeToSurfaceMappingcorrectly, we already dilate the "native" mesh data by 10 mm, and then dilate by 30 mm after we resample to the 32k_fsLR mesh.Ideally, it would be great to have a dscalar that captures the grayordinates that contain "dilated" values themselves, but I don't know how easy that would be to implement...
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I prefer Mike's solution