diff --git a/ICAFIX/ReApplyFixMultiRunPipeline.sh b/ICAFIX/ReApplyFixMultiRunPipeline.sh index e9e4277..887fb9d 100755 --- a/ICAFIX/ReApplyFixMultiRunPipeline.sh +++ b/ICAFIX/ReApplyFixMultiRunPipeline.sh @@ -564,11 +564,11 @@ main() if [ -e ${ConcatNameNoExt}_hp${hp}_wf.txt ] ; then ndhpvol="`cat ${ConcatNameNoExt}_hp${hp}_wf.txt | awk '{print $1}'`" ndhpcifti="`cat ${ConcatNameNoExt}_hp${hp}_wf.txt | awk '{print $2}'`" - ndcifti="`cat ${ConcatNameNoExt}_hp${hp}_wf.txt | awk '{print $3}'`" + ndconcatvol="`cat ${ConcatNameNoExt}_hp${hp}_wf.txt | awk '{print $3}'`" elif [ -z "$p_WF" ] ; then ndhpvol="`echo $p_WF | cut -d ',' -f1`" ndhpcifti="`echo $p_WF | cut -d ',' -f2`" - ndcifti="`echo $p_WF | cut -d ',' -f3`" + ndconcatvol="`echo $p_WF | cut -d ',' -f3`" else log_Err_Abort "ERROR: cannot find Ndist used in MR-FIX. Please use option --wf to set Ndist" fi @@ -715,7 +715,7 @@ main() fi # ${hp} needs to be passed in as a string, to handle the hp=pd* case - local matlab_cmd="${ML_PATHS} functionhighpassandvariancenormalize(${tr}, '${hp}', '${fmri}', '${Caret7_Command}', '${RegString}');" + local matlab_cmd="${ML_PATHS} functionhighpassandvariancenormalize(${tr}, '${hp}', '${fmri}', '${Caret7_Command}','${RegString}', '${ndhpvol}', '${ndhpcifti}', '${ndconcatvol}');" log_Msg "Run interpreted MATLAB/Octave (${interpreter[@]}) with command..." log_Msg "${matlab_cmd}" diff --git a/ICAFIX/ReApplyFixMultiRunPipeline_RIKEN_merge.sh b/ICAFIX/ReApplyFixMultiRunPipeline_RIKEN_merge.sh index 22e1181..c3639fb 100755 --- a/ICAFIX/ReApplyFixMultiRunPipeline_RIKEN_merge.sh +++ b/ICAFIX/ReApplyFixMultiRunPipeline_RIKEN_merge.sh @@ -564,11 +564,11 @@ main() if [ -e ${ConcatNameNoExt}_hp${hp}_wf.txt ] ; then ndhpvol="`cat ${ConcatNameNoExt}_hp${hp}_wf.txt | awk '{print $1}'`" ndhpcifti="`cat ${ConcatNameNoExt}_hp${hp}_wf.txt | awk '{print $2}'`" - ndcifti="`cat ${ConcatNameNoExt}_hp${hp}_wf.txt | awk '{print $3}'`" + ndconcatvol="`cat ${ConcatNameNoExt}_hp${hp}_wf.txt | awk '{print $3}'`" elif [ -z "$p_WF" ] ; then ndhpvol="`echo $p_WF | cut -d ',' -f1`" ndhpcifti="`echo $p_WF | cut -d ',' -f2`" - ndcifti="`echo $p_WF | cut -d ',' -f3`" + ndconcatvol="`echo $p_WF | cut -d ',' -f3`" else log_Err_Abort "ERROR: cannot find Ndist used in MR-FIX. Please use option --wf to set Ndist" fi @@ -715,8 +715,8 @@ fi fi # ${hp} needs to be passed in as a string, to handle the hp=pd* case - local matlab_cmd="${ML_PATHS} functionhighpassandvariancenormalize(${tr}, '${hp}', '${fmri}', '${Caret7_Command}', '${RegString}');" - + local matlab_cmd="${ML_PATHS} functionhighpassandvariancenormalize(${tr}, '${hp}', '${fmri}', '${Caret7_Command}','${RegString}', '${ndhpvol}', '${ndhpcifti}', '${ndconcatvol}');" + log_Msg "Run interpreted MATLAB/Octave (${interpreter[@]}) with command..." log_Msg "${matlab_cmd}" @@ -1088,22 +1088,6 @@ fi Start=`echo "${Start} + ${NumTPS}" | bc -l` done - ## --------------------------------------------------------------------------- - ## Remove all the large time series files in ${ConcatFolder} - ## --------------------------------------------------------------------------- - - ## Deleting these files would save a lot of space. - ## But downstream scripts (e.g., RestingStateStats) assume they exist, and - ## if deleted they would therefore need to be re-created "on the fly" later - - # cd ${ConcatFolder} - # log_Msg "Removing large (concatenated) time series files from ${ConcatFolder}" - # $FSLDIR/bin/imrm ${concatfmri} - # $FSLDIR/bin/imrm ${concatfmri}_hp${hp} - # $FSLDIR/bin/imrm ${concatfmri}_hp${hp}_clean - # /bin/rm -f ${concatfmri}_Atlas${RegString}.dtseries.nii - # /bin/rm -f ${concatfmri}_Atlas${RegString}_hp${hp}.dtseries.nii - # /bin/rm -f ${concatfmri}_Atlas${RegString}_hp${hp}_clean.dtseries.nii cd ${DIR} log_Msg "Completed!" @@ -1113,14 +1097,6 @@ fi # "Global" processing - everything above here should be in a function # ------------------------------------------------------------------------------ -set -e # If any command exits with non-zero value, this script exits - -# Establish defaults -G_DEFAULT_REG_NAME="NONE" -G_DEFAULT_LOW_RES_MESH=32 -G_DEFAULT_MATLAB_RUN_MODE=1 # Use interpreted MATLAB -G_DEFAULT_MOTION_REGRESSION="FALSE" - # Set global variables g_script_name=$(basename "${0}") @@ -1159,9 +1135,6 @@ log_Check_Env_Var FSL_FIXDIR # Show tool versions show_tool_versions -# Establish default MATLAB run mode -G_DEFAULT_MATLAB_RUN_MODE=1 # Use interpreted MATLAB - # Establish default low res mesh for NHP if [ "$SPECIES" = "Macaque" ] ; then G_DEFAULT_LOW_RES_MESH=10 diff --git a/ICAFIX/ReApplyFixPipeline_RIKEN_merge.sh b/ICAFIX/ReApplyFixPipeline_RIKEN_merge.sh new file mode 100755 index 0000000..05264ee --- /dev/null +++ b/ICAFIX/ReApplyFixPipeline_RIKEN_merge.sh @@ -0,0 +1,750 @@ +#!/bin/bash +# +# # ReApplyFixPipeline.sh +# +# ## Copyright Notice +# +# Copyright (C) 2015-2017 The Human Connectome Project/Connectome Coordination Facility +# +# * Washington University in St. Louis +# * University of Minnesota +# * Oxford University +# +# ## Author(s) +# +# * Matthew F. Glasser, Department of Anatomy and Neurobiology, Washington University in St. Louis +# * Timothy B. Brown, Neuroinformatics Research Group, Washington University in St. Louis +# +# ## Product +# +# [Human Connectome Project][HCP] (HCP) Pipelines +# +# ## License +# +# See the [LICENSE](https://github.com/Washington-Univesity/Pipelines/blob/master/LICENSE.md) file +# +# +# [HCP]: http://www.humanconnectome.org +# + +# ------------------------------------------------------------------------------ +# Show usage information for this script +# ------------------------------------------------------------------------------ + +show_usage() +{ + cat < = user supplied value + + Note: The PARAMETERS can be specified positionally (i.e. without using the --param=value + form) by simply specifying all values on the command line in the order they are + listed below. + + e.g. ${g_script_name} /path/to/study/folder 100307 rfMRI_REST1_LR 2000 ... + + [--help] : show usage information and exit + --path= OR --study-folder= + --subject= + --fmri-name= (Do not include path, nifti extension, or the 'hp' string). + --high-pass= + [--reg-name= defaults to ${G_DEFAULT_REG_NAME}. (Use NONE for MSMSulc registration) + [--low-res-mesh=] defaults to ${G_DEFAULT_LOW_RES_MESH} + [--matlab-run-mode={0, 1, 2}] defaults to ${G_DEFAULT_MATLAB_RUN_MODE} + 0 = Use compiled MATLAB + 1 = Use interpreted MATLAB + 2 = Use interpreted Octave + [--motion-regression={TRUE, FALSE}] defaults to ${G_DEFAULT_MOTION_REGRESSION} + [--delete-intermediates={TRUE, FALSE}] defaults to ${G_DEFAULT_DELETE_INTERMEDIATES} + If TRUE, deletes the high-pass filtered files. + +EOF +} + +# Establish defaults +G_DEFAULT_REG_NAME="NONE" +if [[ $SPECIES =~ "Macaque ]] ; then + G_DEFAULT_LOW_RES_MESH=10 +elif [[ $SPECIES =~ "Marmoset" ]] ; then + G_DEFAULT_LOW_RES_MESH=4 +else + G_DEFAULT_LOW_RES_MESH=32 +fi +G_DEFAULT_MATLAB_RUN_MODE=1 # Use interpreted MATLAB +G_DEFAULT_MOTION_REGRESSION="FALSE" +G_DEFAULT_DELETE_INTERMEDIATES="FALSE" + +# ------------------------------------------------------------------------------ +# Get the command line options for this script. +# ------------------------------------------------------------------------------ + +get_options() +{ + local arguments=("$@") + + # initialize global output variables + unset p_StudyFolder # ${1} + unset p_Subject # ${2} + unset p_fMRIName # ${3} + unset p_HighPass # ${4} + unset p_RegName # ${5} + unset p_LowResMesh # ${6} + unset p_MatlabRunMode # ${7} + unset p_MotionRegression # ${8} + unset p_DeleteIntermediates # ${9} + + # set default values + p_RegName=${G_DEFAULT_REG_NAME} + p_LowResMesh=${G_DEFAULT_LOW_RES_MESH} + p_MatlabRunMode=${G_DEFAULT_MATLAB_RUN_MODE} + p_MotionRegression=${G_DEFAULT_MOTION_REGRESSION} + p_DeleteIntermediates=${G_DEFAULT_DELETE_INTERMEDIATES} + + # parse arguments + local num_args=${#arguments[@]} + local argument + local index=0 + + while [ "${index}" -lt "${num_args}" ]; do + argument=${arguments[index]} + + case ${argument} in + --help) + show_usage + exit 0 + ;; + --path=*) + p_StudyFolder=${argument#*=} + index=$(( index + 1 )) + ;; + --study-folder=*) + p_StudyFolder=${argument#*=} + index=$(( index + 1 )) + ;; + --subject=*) + p_Subject=${argument#*=} + index=$(( index + 1 )) + ;; + --fmri-name=*) + p_fMRIName=${argument#*=} + index=$(( index + 1 )) + ;; + --high-pass=*) + p_HighPass=${argument#*=} + index=$(( index + 1 )) + ;; + --reg-name=*) + p_RegName=${argument#*=} + index=$(( index + 1 )) + ;; + --low-res-mesh=*) + p_LowResMesh=${argument#*=} + index=$(( index + 1 )) + ;; + --matlab-run-mode=*) + p_MatlabRunMode=${argument#*=} + index=$(( index + 1 )) + ;; + --motion-regression=*) + p_MotionRegression=${argument#*=} + index=$(( index + 1 )) + ;; + --delete-intermediates=*) + p_DeleteIntermediates=${argument#*=} + index=$(( index + 1 )) + ;; + *) + show_usage + log_Err_Abort "unrecognized option: ${argument}" + ;; + esac + done + + local error_count=0 + + # check required parameters + if [ -z "${p_StudyFolder}" ]; then + log_Err "Study Folder (--path= or --study-folder=) required" + error_count=$(( error_count + 1 )) + else + log_Msg "Study Folder: ${p_StudyFolder}" + fi + + if [ -z "${p_Subject}" ]; then + log_Err "Subject ID (--subject=) required" + error_count=$(( error_count + 1 )) + else + log_Msg "Subject ID: ${p_Subject}" + fi + + if [ -z "${p_fMRIName}" ]; then + log_Err "fMRI Name (--fmri-name=) required" + error_count=$(( error_count + 1 )) + else + log_Msg "fMRI Name: ${p_fMRIName}" + fi + + if [ -z "${p_HighPass}" ]; then + log_Err "High Pass (--high-pass=) required" + error_count=$(( error_count + 1 )) + else + log_Msg "High Pass: ${p_HighPass}" + fi + + if [ -z "${p_RegName}" ]; then + log_Err "Reg Name (--reg-name=) required" + error_count=$(( error_count + 1 )) + else + log_Msg "Reg Name: ${p_RegName}" + fi + + if [ -z "${p_LowResMesh}" ]; then + log_Err "Low Res Mesh (--low-res-mesh=) required" + error_count=$(( error_count + 1 )) + else + log_Msg "Low Res Mesh: ${p_LowResMesh}" + fi + + if [ -z "${p_MatlabRunMode}" ]; then + log_Err "MATLAB run mode value (--matlab-run-mode=) required" + error_count=$(( error_count + 1 )) + else + case ${p_MatlabRunMode} in + 0) + log_Msg "MATLAB Run Mode: ${p_MatlabRunMode} - Use compiled MATLAB" + ;; + 1) + log_Msg "MATLAB Run Mode: ${p_MatlabRunMode} - Use interpreted MATLAB" + ;; + 2) + log_Msg "MATLAB Run Mode: ${p_MatlabRunMode} - Use interpreted Octave" + ;; + *) + log_Err "MATLAB Run Mode value must be 0, 1, or 2" + error_count=$(( error_count + 1 )) + ;; + esac + fi + + if [ -z "${p_MotionRegression}" ]; then + log_Err "motion regression setting (--motion-regression=) required" + error_count=$(( error_count + 1 )) + else + log_Msg "Motion Regression: ${p_MotionRegression}" + fi + + if [ -z "${p_DeleteIntermediates}" ]; then + log_Err "delete intermediates setting (--delete-intermediates=) required" + error_count=$(( error_count + 1 )) + else + log_Msg "Delete Intermediates: ${p_DeleteIntermediates}" + fi + + if [ ${error_count} -gt 0 ]; then + log_Err_Abort "For usage information, use --help" + fi +} + +# ------------------------------------------------------------------------------ +# Show Tool Versions +# ------------------------------------------------------------------------------ + +show_tool_versions() +{ + # Show HCP pipelines version + log_Msg "Showing HCP Pipelines version" + cat ${HCPPIPEDIR}/version.txt + + # Show wb_command version + log_Msg "Showing Connectome Workbench (wb_command) version" + ${CARET7DIR}/wb_command -version + + # Show FSL version + log_Msg "Showing FSL version" + fsl_version_get fsl_ver + log_Msg "FSL version: ${fsl_ver}" + + # Show specific FIX version, if available + if [ -f ${FSL_FIXDIR}/fixversion ]; then + fixversion=$(cat ${FSL_FIXDIR}/fixversion ) + log_Msg "FIX version: $fixversion" + fi + +} + +# ------------------------------------------------------------------------------ +# Check for whether or not we have hand reclassification files +# ------------------------------------------------------------------------------ + +have_hand_reclassification() +{ + local StudyFolder="${1}" + local Subject="${2}" + local fMRIName="${3}" + local HighPass="${4}" + + if (( HighPass >= 0 )); then + [ -e "${StudyFolder}/${Subject}/MNINonLinear/Results/${fMRIName}/${fMRIName}_hp${HighPass}.ica/HandNoise.txt" ] + else + [ -e "${StudyFolder}/${Subject}/MNINonLinear/Results/${fMRIName}/${fMRIName}.ica/HandNoise.txt" ] + fi +} + +function interpret_as_bool() +{ + case $(echo "$1" | tr '[:upper:]' '[:lower:]') in + (true | yes | 1) + echo 1 + ;; + (false | no | none | 0) + echo 0 + ;; + (*) + log_Err_Abort "error: '$1' is not valid for this argument, please use TRUE or FALSE" + ;; + esac +} + +# ------------------------------------------------------------------------------ +# Main processing of script. +# ------------------------------------------------------------------------------ + +main() +{ + log_Msg "Starting main functionality" + + # Retrieve positional parameters + local StudyFolder="${1}" + local Subject="${2}" + local fMRIName="${3}" + local HighPass="${4}" + + local RegName + if [ -z "${5}" ]; then + RegName=${G_DEFAULT_REG_NAME} + else + RegName="${5}" + fi + + local LowResMesh + if [ -z "${6}" ]; then + LowResMesh=${G_DEFAULT_LOW_RES_MESH} + else + LowResMesh="${6}" + fi + + local MatlabRunMode + if [ -z "${7}" ]; then + MatlabRunMode=${G_DEFAULT_MATLAB_RUN_MODE} + else + MatlabRunMode="${7}" + fi + + local MotionRegression + if [ -z "${8}" ]; then + MotionRegression=$(interpret_as_bool "${G_DEFAULT_MOTION_REGRESSION}") + else + MotionRegression=$(interpret_as_bool "${8}") + fi + + local DeleteIntermediates + if [ -z "${9}" ]; then + DeleteIntermediates=$(interpret_as_bool "${G_DEFAULT_DELETE_INTERMEDIATES}") + else + DeleteIntermediates=$(interpret_as_bool "${9}") + fi + + # Log values retrieved from positional parameters + log_Msg "StudyFolder: ${StudyFolder}" + log_Msg "Subject: ${Subject}" + log_Msg "fMRIName: ${fMRIName}" + log_Msg "HighPass: ${HighPass}" + log_Msg "RegName: ${RegName}" + log_Msg "LowResMesh: ${LowResMesh}" + log_Msg "MatlabRunMode: ${MatlabRunMode}" + log_Msg "MotionRegression: ${MotionRegression}" + log_Msg "DeleteIntermediates: ${DeleteIntermediates}" + + # Naming Conventions and other variables + local Caret7_Command="${CARET7DIR}/wb_command" + log_Msg "Caret7_Command: ${Caret7_Command}" + + local RegString + if [ ${RegName} != "NONE" ] ; then + RegString="_${RegName}" + else + RegString="" + fi + + if [ ! -z ${LowResMesh} ] && [ ${LowResMesh} != ${G_DEFAULT_LOW_RES_MESH} ]; then + RegString="${RegString}.${LowResMesh}k" + fi + + log_Msg "RegString: ${RegString}" + + # For INTERPRETED MODES, make sure that matlab/octave has access to the functions it needs. + # Since we are NOT using the ${FSL_FIXDIR}/call_matlab.sh script to invoke matlab (unlike 'hcp_fix') + # we need to explicitly add ${FSL_FIXDIR} (all the fix-related functions) + # and ${FSL_MATLAB_PATH} (e.g., read_avw.m, save_avw.m) to the matlab path. + # Several additional necessary environment variables (e.g., ${FSL_FIX_CIFTIRW} and ${FSL_FIX_WBC}) + # are set in ${FSL_FIXDIR}/settings.sh, which is sourced below for interpreted modes. + # Note that the ciftiopen.m, ciftisave.m functions are *appended* to the path through the ${FSL_FIX_CIFTIRW} + # environment variable within fix_3_clean.m itself. + # Note that we are NOT adding '${HCPPIPEDIR}/global/matlab' to the path (which also contains the + # CIFTI I/O functions as well), so the versions in ${FSL_FIX_CIFTIRW} will be the ones that get used + # (assuming that the user hasn't further customized their matlab path, e.g., via a startup.m addition). + export FSL_MATLAB_PATH="${FSLDIR}/etc/matlab" + local ML_PATHS="addpath('${FSL_FIXDIR}'); addpath('${FSL_MATLAB_PATH}');" + + # Some defaults + local aggressive=0 + local hp=${HighPass} + local DoVol=0 + local fixlist=".fix" + + if (( hp >= 0 )); then + hpStr="_hp${hp}" + else + hpStr="" + fi + + # fMRIName is expected to NOT include path info, or a nifti extension; make sure that is indeed the case + # (although if someone includes the hp string as part fMRIName itself, the script will still break) + fMRIName=$(basename $($FSLDIR/bin/remove_ext $fMRIName)) + + # If we have a hand classification and no regname, reapply fix to the volume as well + if have_hand_reclassification ${StudyFolder} ${Subject} ${fMRIName} ${hp} + then + fixlist="HandNoise.txt" + #TSC: if regname (which applies to the surface) isn't NONE, assume the hand classification was previously already applied to the volume data + if [[ "${RegName}" == "NONE" ]] + then + DoVol=1 + fi + fi + # WARNING: fix_3_clean doesn't actually do anything different based on the value of DoVol (its 5th argument). + # Rather, if a 5th argument is present, fix_3_clean does NOT apply cleanup to the volume, *regardless* of whether + # that 5th argument is 0 or 1 (or even a non-sensical string such as 'foo'). + # It is for that reason that the code below needs to use separate calls to fix_3_clean, with and without DoVol + # as an argument, rather than simply passing in the value of DoVol as set within this script. + # Not sure if/when this non-intuitive behavior of fix_3_clean will change, but this is accurate as of fix1.067 + # UPDATE (11/8/2019): As of FIX 1.06.12, fix_3_clean interprets its 5th argument ("DoVol") in the usual boolean + # manner. However, since we already had a work-around to this problem, we will leave the code unchanged so that + # we don't need to add a FIX version dependency to the script. + + log_Msg "Use fixlist=$fixlist" + + DIR=$(pwd) + cd ${StudyFolder}/${Subject}/MNINonLinear/Results/${fMRIName} + + # Note: fix_3_clean does NOT filter the volume (NIFTI) data -- it assumes + # that any desired filtering has already been done outside of fix. + # So here, we need to symlink to the hp-filtered volume data. + # HOWEVER, if missing, only need to generate the hp-filtered volume data if DoVol=1. + # Otherwise (if DoVol=0), the only role of filtered_func_data in fix_3_clean is to determine the TR. + # In that case, we will just symlink the *non-filtered* data to filtered_func_data + # (as a hack for fix_3_clean to determine the TR without a time-consuming filtering step + # on the volume) + + useNonFilteredAsFilteredFunc=0 + if (( $($FSLDIR/bin/imtest "${fMRIName}${hpStr}") )); then + log_Warn "Using existing $($FSLDIR/bin/imglob -extension ${fMRIName}${hpStr}) (not re-filtering)" + else # hp filtered volume file doesn't exist + if (( DoVol )); then # need to actually refilter the volume data + if (( hp > 0 )); then + tr=`$FSLDIR/bin/fslval ${fMRIName} pixdim4` + log_Msg "tr: ${tr}" + log_Msg "processing FMRI file ${fMRIName} with highpass ${hp}" + hptr=$(echo "scale = 10; $hp / (2 * $tr)" | bc -l) + + # Starting with FSL 5.0.7, 'fslmaths -bptf' no longer includes the temporal mean in its output. + # A work-around to this, which works with both the pre- and post-5.0.7 behavior is to compute + # the temporal mean, remove it, run -bptf, and then add the mean back in. + ${FSLDIR}/bin/fslmaths ${fMRIName} -Tmean ${fMRIName}${hpStr} + highpass_cmd="${FSLDIR}/bin/fslmaths ${fMRIName} -sub ${fMRIName}${hpStr} -bptf ${hptr} -1 -add ${fMRIName}${hpStr} ${fMRIName}${hpStr}" + log_Msg "highpass_cmd: ${highpass_cmd}" + ${highpass_cmd} + elif (( hp == 0 )); then + # Nothing in script currently detrends the volume if hp=0 is requested (which is the intended meaning of hp=0) + log_Err_Abort "hp = ${hp} not currently supported" + fi + else + useNonFilteredAsFilteredFunc=1 + fi + fi + + if [ ! -e ${fMRIName}${hpStr}.ica ]; then + log_Err_Abort "${fMRIName}${hpStr}.ica is expected to already exist, but does not" + fi + + cd ${fMRIName}${hpStr}.ica + + # Create symlink for filtered_func_data (per comments above) + if (( useNonFilteredAsFilteredFunc )); then + $FSLDIR/bin/imln ../${fMRIName} filtered_func_data + else + $FSLDIR/bin/imln ../${fMRIName}${hpStr} filtered_func_data + fi + + # However, hp-filtering of the *CIFTI* (dtseries) occurs within fix_3_clean. + # So here, we just create a symlink with the file name expected by + # fix_3_clean ("Atlas.dtseries.nii") to the non-filtered data. + if [ -f ../${fMRIName}_Atlas${RegString}.dtseries.nii ] ; then + log_Msg "FOUND FILE: ../${fMRIName}_Atlas${RegString}.dtseries.nii" + log_Msg "Performing imln" + + /bin/rm -f Atlas.dtseries.nii + $FSLDIR/bin/imln ../${fMRIName}_Atlas${RegString}.dtseries.nii Atlas.dtseries.nii + + log_Msg "START: Showing linked files" + ls -l ../${fMRIName}_Atlas${RegString}.dtseries.nii + ls -l Atlas.dtseries.nii + log_Msg "END: Showing linked files" + else + log_Warn "FILE NOT FOUND: ../${fMRIName}_Atlas${RegString}.dtseries.nii" + fi + + # Get Movement_Regressors.txt into the format expected by functionmotionconfounds.m + mkdir -p mc + if [ -f ../Movement_Regressors.txt ] ; then + log_Msg "Creating mc/prefiltered_func_data_mcf.par file" + cat ../Movement_Regressors.txt | awk '{ print $4 " " $5 " " $6 " " $1 " " $2 " " $3}' > mc/prefiltered_func_data_mcf.par + else + log_Err_Abort "Movement_Regressors.txt not retrieved properly." + fi + + ## --------------------------------------------------------------------------- + ## Run fix_3_clean + ## --------------------------------------------------------------------------- + + # MPH: We need to invoke fix_3_clean directly, rather than through 'fix -a ' because + # the latter does not provide access to the "DoVol" option within the former. + # (Also, 'fix -a' is hard-coded to use '.fix' as the list of noise components, although that + # could be worked around). + + export FSL_FIX_WBC="${Caret7_Command}" + # WARNING: fix_3_clean uses the environment variable FSL_FIX_WBC, but most previous + # versions of FSL_FIXDIR/settings.sh (v1.067 and earlier) have a hard-coded value for + # FSL_FIX_WBC, and don't check whether it is already defined in the environment. + # Thus, when settings.sh file gets sourced, there is a possibility that the version of + # wb_command is no longer the same as that specified by ${Caret7_Command}. So, after + # sourcing settings.sh below, we explicitly set FSL_FIX_WBC back to value of ${Caret7_Command}. + # (This may only be relevant for interpreted matlab/octave modes). + + log_Msg "Running fix_3_clean" + + case ${MatlabRunMode} in + + # See important WARNING above regarding why ${DoVol} is NOT included as an argument when DoVol=1 !! + + 0) + # Use Compiled MATLAB + + local matlab_exe="${FSL_FIXDIR}/compiled/$(uname -s)/$(uname -m)/run_fix_3_clean.sh" + + # Do NOT enclose string variables inside an additional single quote because all + # variables are already passed into the compiled binary as strings + local matlab_function_arguments=("${fixlist}" "${aggressive}" "${MotionRegression}" "${hp}") + if (( ! DoVol )); then + matlab_function_arguments+=("${DoVol}") + fi + + # fix_3_clean is part of the FIX distribution. + # If ${FSL_FIX_MCR} is already defined in the environment, use that for the MCR location. + # If not, the appropriate MCR version for use with fix_3_clean should be set in $FSL_FIXDIR/settings.sh. + if [ -z "${FSL_FIX_MCR}" ]; then + debug_disable_trap + source ${FSL_FIXDIR}/settings.sh + debug_enable_trap + export FSL_FIX_WBC="${Caret7_Command}" + # If FSL_FIX_MCR is still not defined after sourcing settings.sh, we have a problem + if [ -z "${FSL_FIX_MCR}" ]; then + log_Err_Abort "To use MATLAB run mode: ${MatlabRunMode}, the FSL_FIX_MCR environment variable must be set" + fi + fi + log_Msg "FSL_FIX_MCR: ${FSL_FIX_MCR}" + + local matlab_cmd=("${matlab_exe}" "${FSL_FIX_MCR}" "${matlab_function_arguments[@]}") + + # redirect tokens must be parsed by bash before doing variable expansion, and thus can't be inside a variable + # MPH: Going to let Compiled MATLAB use the existing stdout and stderr, rather than creating a separate log file + #local matlab_logfile=".reapplyfix.${fMRIName}${RegString}.fix_3_clean.matlab.log" + #log_Msg "Run MATLAB command: ${matlab_cmd[*]} >> ${matlab_logfile} 2>&1" + #"${matlab_cmd[@]}" >> "${matlab_logfile}" 2>&1 + log_Msg "Run compiled MATLAB: ${matlab_cmd[*]}" + "${matlab_cmd[@]}" + ;; + + 1 | 2) + # Use interpreted MATLAB or Octave + if [[ ${MatlabRunMode} == "1" ]]; then + local interpreter=(matlab -nojvm -nodisplay -nosplash) + else + local interpreter=(octave-cli -q --no-window-system) + fi + + if (( DoVol )); then + local matlab_cmd="${ML_PATHS} fix_3_clean('${fixlist}',${aggressive},${MotionRegression},${hp});" + else + local matlab_cmd="${ML_PATHS} fix_3_clean('${fixlist}',${aggressive},${MotionRegression},${hp},${DoVol});" + fi + + log_Msg "Run interpreted MATLAB/Octave (${interpreter[@]}) with command..." + log_Msg "${matlab_cmd}" + + # Use bash redirection ("here-string") to pass multiple commands into matlab + # (Necessary to protect the semicolons that separate matlab commands, which would otherwise + # get interpreted as separating different bash shell commands) + (debug_disable_trap; source "${FSL_FIXDIR}/settings.sh"; debug_enable_trap; export FSL_FIX_WBC="${Caret7_Command}"; "${interpreter[@]}" <<<"${matlab_cmd}") + ;; + + *) + # Unsupported MATLAB run mode + log_Err_Abort "Unsupported MATLAB run mode value: ${MatlabRunMode}" + ;; + esac + + log_Msg "Done running fix_3_clean" + + ## --------------------------------------------------------------------------- + ## Rename some files (relative to the default names coded in fix_3_clean) + ## --------------------------------------------------------------------------- + + # Remove any existing old versions of the cleaned data (normally they should be overwritten + # in the renaming that follows, but this ensures that any old versions don't linger) + cd ${StudyFolder}/${Subject}/MNINonLinear/Results/${fMRIName} + local fmri=${fMRIName} + if (( hp >= 0 )); then + local fmrihp=${fmri}_hp${hp} + else + local fmrihp=${fmri} + fi + + /bin/rm -f ${fmri}_Atlas${RegString}${hpStr}_clean.dtseries.nii + /bin/rm -f ${fmri}_Atlas${RegString}${hpStr}_clean_vn.dscalar.nii + + if (( DoVol )); then + $FSLDIR/bin/imrm ${fmrihp}_clean + $FSLDIR/bin/imrm ${fmrihp}_clean_vn + fi + + # Rename some of the outputs from fix_3_clean. + # Note that the variance normalization ("_vn") outputs require use of fix1.067 or later + # So check whether those files exist before moving/renaming them + if [ -f ${fmrihp}.ica/Atlas_clean.dtseries.nii ]; then + /bin/mv ${fmrihp}.ica/Atlas_clean.dtseries.nii ${fmri}_Atlas${RegString}${hpStr}_clean.dtseries.nii + else + log_Err_Abort "Something went wrong; ${fmrihp}.ica/Atlas_clean.dtseries.nii wasn't created" + fi + if [ -f ${fmrihp}.ica/Atlas_clean_vn.dscalar.nii ]; then + /bin/mv ${fmrihp}.ica/Atlas_clean_vn.dscalar.nii ${fmri}_Atlas${RegString}${hpStr}_clean_vn.dscalar.nii + fi + + if (( DoVol )); then + $FSLDIR/bin/immv ${fmrihp}.ica/filtered_func_data_clean ${fmrihp}_clean + if [ "$?" -ne "0" ]; then + log_Err_Abort "Something went wrong; ${fmrihp}.ica/filtered_func_data_clean wasn't created" + fi + if [ `$FSLDIR/bin/imtest ${fmrihp}.ica/filtered_func_data_clean_vn` = 1 ]; then + $FSLDIR/bin/immv ${fmrihp}.ica/filtered_func_data_clean_vn ${fmrihp}_clean_vn + fi + fi + log_Msg "Done renaming files" + + # Remove the 'fake-NIFTI' file created in fix_3_clean for high-pass filtering of the CIFTI (if it exists) + $FSLDIR/bin/imrm ${fmrihp}.ica/Atlas + + # Always delete things with too-generic names + $FSLDIR/bin/imrm ${fmrihp}.ica/filtered_func_data + rm -f ${fmrihp}.ica/Atlas.dtseries.nii + + # Optional deletion of highpass intermediates + if [ "${DeleteIntermediates}" == "1" ] ; then + if (( hp > 0 )); then # fix_3_clean only writes out the hp-filtered time series if hp > 0 + $FSLDIR/bin/imrm ${fmri}_hp${hp} # Explicitly use _hp${hp} here (rather than $hpStr as a safeguard against accidental deletion of the non-hp-filtered timeseries) + rm -f ${fmrihp}.ica/Atlas_hp_preclean.dtseries.nii + fi + else + #even if we don't delete it, don't leave this file with a hard to interpret name + if (( hp > 0 )); then + # 'OR' mv command with "true" to avoid returning an error code if file doesn't exist for some reason + mv -f ${fmrihp}.ica/Atlas_hp_preclean.dtseries.nii ${fmri}_Atlas_hp${hp}.dtseries.nii || true + fi + fi + + cd ${DIR} # Return to directory where script was launched + + log_Msg "Completed!" +} + +# ------------------------------------------------------------------------------ +# "Global" processing - everything above here should be in a function +# ------------------------------------------------------------------------------ + +# Set global variables +g_script_name=$(basename "${0}") + +# Allow script to return a Usage statement, before any other output +if [ "$#" = "0" ]; then + show_usage + exit 1 +fi + +# Verify that HCPPIPEDIR environment variable is set +if [ -z "${HCPPIPEDIR}" ]; then + echo "${g_script_name}: ABORTING: HCPPIPEDIR environment variable must be set" + exit 1 +fi + +# Load function libraries +source "${HCPPIPEDIR}/global/scripts/debug.shlib" "$@" # Debugging functions; also sources log.shlib +source ${HCPPIPEDIR}/global/scripts/opts.shlib # Command line option functions +source "${HCPPIPEDIR}/global/scripts/fsl_version.shlib" # Functions for getting FSL version + +opts_ShowVersionIfRequested $@ + +if opts_CheckForHelpRequest $@; then + show_usage + exit 0 +fi + +${HCPPIPEDIR}/show_version + +# Verify required environment variables are set and log value +log_Check_Env_Var HCPPIPEDIR +log_Check_Env_Var CARET7DIR +log_Check_Env_Var FSLDIR +log_Check_Env_Var FSL_FIXDIR + +# Show tool versions +show_tool_versions + +# Determine whether named or positional parameters are used and invoke the 'main' function +if [[ ${1} == --* ]]; then + # Named parameters (e.g. --parameter-name=parameter-value) are used + log_Msg "Using named parameters" + + # Get command line options + get_options "$@" + + # Invoke main functionality + # ${1} ${2} ${3} ${4} ${5} ${6} ${7} ${8} ${9} + main "${p_StudyFolder}" "${p_Subject}" "${p_fMRIName}" "${p_HighPass}" "${p_RegName}" "${p_LowResMesh}" "${p_MatlabRunMode}" "${p_MotionRegression}" "${p_DeleteIntermediates}" + +else + # Positional parameters are used + log_Msg "Using positional parameters" + main "$@" + +fi diff --git a/ICAFIX/functionhighpassandvariancenormalize.m b/ICAFIX/functionhighpassandvariancenormalize.m index 3d0b733..82d52e5 100644 --- a/ICAFIX/functionhighpassandvariancenormalize.m +++ b/ICAFIX/functionhighpassandvariancenormalize.m @@ -1,94 +1,338 @@ -function [ output_args ] = functionhighpassandvariancenormalize(TR,hp,fmri,WBC) +function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) +% +% FUNCTIONHIGHPASSANDVARIANCENORMALIZE(TR,HP,FMRI,WBC,[REGSTRING]) +% +% This function: +% (1) Detrends / high-pass filters motion confounds, NIFTI (volume) and CIFTI files +% (2) Calls icaDim.m to estimate the dimensionalities of structured and unstructured noise, +% and compute a spatial variance normalization map. +% It is written to support 'hcp_fix_multi_run', and is not intended as a general (stand-alone) +% function supporting detrending and variance normalization. +% +% TR: repetition time between frames, in sec +% HP: determines what high-pass filtering to apply to the motion confounds and data +% HP<0: No highpass applied (demeaning only). No 'hp' string will be added as part of the file names. +% HP>0: The full-width (2*sigma), in sec, to apply using 'fslmaths -bptf'. +% HP='pd#': If HP is a string, in which the first two characters are 'pd', followed by an integer value, +% then polynomial detrending is applied, with the order specified by the integer value +% embedded at the end of the string. +% The output files will include the string '_hppd#' +% HP=0: Gets interpreted as a linear (1st order) detrend +% This is consistent with the convention in FIX of using HP=0 to specify a linear detrend. +% The output file names will include the string '_hp0'. +% FMRI: The base string of the fmri time series (no extensions) +% WBC: wb_command (full path) +% varargin: +% [REGSTRING]: Additional registration-related string to add to output file names. OPTIONAL. +% [ndhpvol ndhpcifti ndvol]: number of distributions highpassed volume, highpassed cifti and concatenated volume. OPTIONAL. + +% Note: HP='pd0' would be interpreted as a true 0th order detrend, which is +% the same as demeaning. Mathematically, this is the same as the HP<0 condition, +% but the output files will be named differently (i.e., include '_hppd0'), and additional +% output files will be written relative to the HP<0 condition. + +% Authors: M. Glasser and M. Harms + +CIFTIMatlabReaderWriter=getenv('FSL_FIX_CIFTIRW'); +if (~isdeployed) + % addpath does not work and should not be used in compiled MATLAB + fprintf('Adding %s to MATLAB path\n', CIFTIMatlabReaderWriter); + addpath(CIFTIMatlabReaderWriter); +end + +%% Defaults +dovol = 1; +regstring = ''; +pdflag = false; % Polynomial detrend +pdstring = 'pd'; % Expected string at start of HP variable to request a "polynomial detrend" +hpstring = ''; + +%% Parse varargin +if length(varargin) == 1 && ~isempty(varargin{1}) + dovol = 0; %regname is only used for a new surface registration, so will never require redoing volume + regstring = varargin{1}; %this has the underscore on the front already + ndhpvol=2 + ndhpcifti=3 + ndvol=3 + if ~ischar(regstring) + error('%s: REGSTRING should be a string', mfilename); + end +elseif length(varargin) >= 3 + if ~isempty(varagin{4}) + dovol = 0 + regstring = varargin{4} + if ~ischar(regstring) + error('%s: REGSTRING should be a string', mfilename); + end + end + ndhpvol=varargin{1} + ndhpcifti=varargin{2} + ndvol=varargin{3} + if ~isint(ndhpvol) + error('%s: NDHPVOL should be an integer', mfilename); + end + if ~isint(ndhpcifti) + error('%s: NDHPCIFTI should be an integer', mfilename); + end + if ~isint(ndvol) + error('%s: NDVOL should be an integer', mfilename); + end +end + %UNTITLED4 Summary of this function goes here % Detailed explanation goes here +% Default value: ndhpvol=2;ndhpcifti=3; ndvol=3; - Takuya Hayashi, Matt Glasser +fprintf('hp=%d, ndist=%d,%d,%d',hp,ndhpvol,ndhpcifti,ndvol) + + +% Check whether polynomial detrend is being requested for the high-pass filtering. +if ischar(hp) + hp = lower(hp); + if strncmp(hp,pdstring,numel(pdstring)) + pdflag = true; + hp = str2double(hp(numel(pdstring)+1:end)); % hp is now a numeric representing the order of the polynomial detrend + elseif ~isempty(str2double(hp)) % Allow for hp to be provided as a string that contains purely numeric elements + hp = str2double(hp); + else error('%s: Invalid specification for the high-pass filter', mfilename); + end +end -cts=single(read_avw([fmri '.nii.gz'])); -ctsX=size(cts,1); ctsY=size(cts,2); ctsZ=size(cts,3); ctsT=size(cts,4); -cts=reshape(cts,ctsX*ctsY*ctsZ,ctsT); +%% Allow for compiled matlab (hp as a string is already handled above) +if (isdeployed) + TR = str2double(TR); +end + +%% Argument checking +if hp==0 && ~pdflag + warning('%s: hp=0 will be interpreted as polynomial detrending of order 1', mfilename); + pdflag=true; + hp=1; + % But in this case, preserve use of "_hp0" in the filename to indicate a linear detrend (for backwards compatibility in file naming) + hpstring = '_hp0'; +end +if pdflag + if (~isscalar(hp) || hp < 0 || hp ~= fix(hp)) + error('%s: Invalid specification for the order of the polynomial detrending', mfilename); + end +end + +%% Finish setting up hpstring to use in file names +if ~pdflag + pdstring = ''; +end +if isempty(hpstring) && hp>=0 + hpstring = ['_hp' pdstring num2str(hp)]; +end +%% Load the motion confounds, and the CIFTI (if hp>=0) (don't need either if hp<0) if hp>=0 - confounds=load([fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf.par']); + confounds=load([fmri hpstring '.ica/mc/prefiltered_func_data_mcf.par']); confounds=confounds(:,1:6); - confounds=functionnormalise([confounds [zeros(1,size(confounds,2)); confounds(2:end,:)-confounds(1:end-1,:)] ]); - confounds=functionnormalise([confounds confounds.*confounds]); + %% normalise function is in HCPPIPEDIR/global/matlab/normalise.m + confounds=normalise([confounds [zeros(1,size(confounds,2)); confounds(2:end,:)-confounds(1:end-1,:)] ]); + confounds=normalise([confounds confounds.*confounds]); - BO=ciftiopen([fmri '_Atlas.dtseries.nii'],WBC); -else - cts=demean(cts')'; + BO=ciftiopen([fmri '_Atlas' regstring '.dtseries.nii'],WBC); end -if hp==0 - save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]); - confounds=detrend(confounds); - save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf_hp'],'f',[1 1 1 TR]); - - cts=detrend(cts')'; - save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_hp' num2str(hp) '.nii.gz'],'f',[1 1 1 TR]); - call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_hp' num2str(hp) '.nii.gz']); +%% Apply hp filtering of the motion confounds, volume (if requested), and CIFTI +if pdflag % polynomial detrend case + if dovol > 0 + % Save and filter confounds, as a NIFTI, as expected by FIX + save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]); + confounds=detrendpoly(confounds,hp); + save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf_hp'],'f',[1 1 1 TR]); + + % Load volume time series and reduce to just the non-zero voxels (for memory efficiency) + % Note: Use 'range' to identify non-zero voxels (which is very memory efficient) + % rather than 'std' (which requires additional memory equal to the size of the input) + ctsfull=single(read_avw([fmri '.nii.gz'])); + ctsX=size(ctsfull,1); ctsY=size(ctsfull,2); ctsZ=size(ctsfull,3); ctsT=size(ctsfull,4); + ctsfull=reshape(ctsfull,ctsX*ctsY*ctsZ,ctsT); + ctsmask=range(ctsfull, 2) > 0; + fprintf('Non-empty voxels: %d (= %.2f%% of %d)\n', sum(ctsmask), 100*sum(ctsmask)/size(ctsfull,1), size(ctsfull,1)); + cts=ctsfull(ctsmask,:); + clear ctsfull; + + % Polynomial detrend + cts=detrendpoly(cts',hp)'; + + % Write out result, restoring to original size + ctsfull=zeros(ctsX*ctsY*ctsZ,ctsT, 'single'); + ctsfull(ctsmask,:)=cts; + save_avw(reshape(ctsfull,ctsX,ctsY,ctsZ,ctsT),[fmri hpstring '.nii.gz'],'f',[1 1 1 TR]); + clear ctsfull; + % Use -d flag in fslcpgeom (even if not technically necessary) to reduce possibility of mistakes when editing script + call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri hpstring '.nii.gz -d']); + end - BO.cdata=detrend(BO.cdata')'; - ciftisave(BO,[fmri '_Atlas_hp' num2str(hp) '.dtseries.nii'],WBC); -end -if hp>0 - save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]); - call_fsl(sprintf(['fslmaths ' fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf -bptf %f -1 ' fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf_hp'],0.5*hp/TR)); - - save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_hp' num2str(hp) '.nii.gz'],'f',[1 1 1 TR]); - call_fsl(['fslmaths ' fmri '_hp' num2str(hp) '.nii.gz -bptf ' num2str(0.5*hp/TR) ' -1 ' fmri '_hp' num2str(hp) '.nii.gz']); - cts=single(read_avw([fmri '_hp' num2str(hp) '.nii.gz'])); - cts=reshape(cts,ctsX*ctsY*ctsZ,ctsT); - call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_hp' num2str(hp) '.nii.gz']); + BO.cdata=detrendpoly(BO.cdata',hp)'; + ciftisave(BO,[fmri '_Atlas' regstring hpstring '.dtseries.nii'],WBC); + +elseif hp>0 % "fslmaths -bptf" based filtering + if dovol > 0 + % Save and filter confounds, as a NIFTI, as expected by FIX + save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]); + call_fsl(sprintf(['fslmaths ' fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf -bptf %f -1 ' fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf_hp'],0.5*hp/TR)); + + % bptf filtering; no masking here, so this is probably memory inefficient + call_fsl(['fslmaths ' fmri '.nii.gz -bptf ' num2str(0.5*hp/TR) ' -1 ' fmri hpstring '.nii.gz']); + call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri hpstring '.nii.gz -d']); + + % Load in result and reduce to just the non-zero voxels (for memory efficiency going forward) + ctsfull=single(read_avw([fmri hpstring '.nii.gz'])); + ctsX=size(ctsfull,1); ctsY=size(ctsfull,2); ctsZ=size(ctsfull,3); ctsT=size(ctsfull,4); + ctsfull=reshape(ctsfull,ctsX*ctsY*ctsZ,ctsT); + ctsmask=range(ctsfull, 2) > 0; + fprintf('Non-empty voxels: %d (= %.2f%% of %d)\n', sum(ctsmask), 100*sum(ctsmask)/size(ctsfull,1), size(ctsfull,1)); + cts=ctsfull(ctsmask,:); + clear ctsfull; + end BOdimX=size(BO.cdata,1); BOdimZnew=ceil(BOdimX/100); BOdimT=size(BO.cdata,2); save_avw(reshape([BO.cdata ; zeros(100*BOdimZnew-BOdimX,BOdimT)],10,10,BOdimZnew,BOdimT),'Atlas','f',[1 1 1 TR]); call_fsl(sprintf('fslmaths Atlas -bptf %f -1 Atlas',0.5*hp/TR)); - grot=reshape(single(read_avw('Atlas')),100*BOdimZnew,BOdimT); BO.cdata=grot(1:BOdimX,:); clear grot; - call_fsl('rm Atlas.nii.gz'); - ciftisave(BO,[fmri '_Atlas_hp' num2str(hp) '.dtseries.nii'],WBC); + grot=reshape(single(read_avw('Atlas')),100*BOdimZnew,BOdimT); BO.cdata=grot(1:BOdimX,:); clear grot; + ciftisave(BO,[fmri '_Atlas' regstring hpstring '.dtseries.nii'],WBC); + delete('Atlas.nii.gz'); + +elseif hp<0 % If no hp filtering, still need to at least demean the volumetric time series + if dovol > 0 + + % Load volume time series and reduce to just the non-zero voxels (for memory efficiency) + ctsfull=single(read_avw([fmri '.nii.gz'])); + ctsX=size(ctsfull,1); ctsY=size(ctsfull,2); ctsZ=size(ctsfull,3); ctsT=size(ctsfull,4); + ctsfull=reshape(ctsfull,ctsX*ctsY*ctsZ,ctsT); + ctsmask=range(ctsfull, 2) > 0; + fprintf('Non-empty voxels: %d (= %.2f%% of %d)\n', sum(ctsmask), 100*sum(ctsmask)/size(ctsfull,1), size(ctsfull,1)); + cts=ctsfull(ctsmask,:); + clear ctsfull; + + cts=demean(cts')'; + end + end -%Compute VN + +% NOTE: To interpret the logic of the following "Compute, Save, Apply" sections need to know that this function expects that: +% (1) individual run time series will be passed in with hp requested +% (2) concatenated (multi-run) time series will be passed in with hp<0 +%% [This is a bit of hack: would be cleaner if we disentangled the operation of this function from the specific needs of 'hcp_fix_multi_run']. + +%% Compute VN map if hp>=0 - Outcts=icaDim(cts,0,1,2,2); %Volume fits two distributions to deal with processing interpolation and multiband - OutBO=icaDim(BO.cdata,0,1,2,3); %CIFTI fits three because of the effects of volume to CIFTI mapping and regularization + if dovol > 0 + Outcts=icaDim(cts,0,1,-1,ndhpvol); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, Volume fits two distributions to deal with MNI transform + end + + OutBO=icaDim(BO.cdata,0,1,-1,ndhpcifti); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, CIFTI fits three distributions to deal with volume to CIFTI mapping else - Outcts=icaDim(cts,0,1,2,3); %Volume fits three distributions to deal with processing interpolation and multiband and concatination effects + if dovol > 0 + Outcts=icaDim(cts,0,1,-1,ndvol); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, Volume fits two distributions to deal with MNI transform + end end -%Save VN +%% Save VN map +% (Only need these for the individual runs, which, per above, are expected to be passed in with hp requested). if hp>=0 - save_avw(reshape(Outcts.noise_unst_std,ctsX,ctsY,ctsZ,1),[fmri '_vn.nii.gz'],'f',[1 1 1 1]); - call_fsl(['fslcpgeom ' fmri '_mean.nii.gz ' fmri '_vn.nii.gz -d']); -end + if dovol > 0 + fname=[fmri hpstring '_vn.nii.gz']; + vnfull=zeros(ctsX*ctsY*ctsZ,1, 'single'); + vnfull(ctsmask)=Outcts.noise_unst_std; + + save_avw(reshape(vnfull,ctsX,ctsY,ctsZ,1),fname,'f',[1 1 1 1]); + clear vnfull; + call_fsl(['fslcpgeom ' fmri '_mean.nii.gz ' fname ' -d']); + end -if hp>=0 VN=BO; VN.cdata=OutBO.noise_unst_std; - ciftisavereset(VN,[fmri '_Atlas_vn.dscalar.nii'],WBC); + fname=[fmri '_Atlas' regstring hpstring '_vn.dscalar.nii']; + disp(['saving ' fname]); + ciftisavereset(VN,fname,WBC); end -%Apply VN and Save HP_VN TCS -cts=cts./repmat(Outcts.noise_unst_std,1,ctsT); -if hp>=0 - save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_hp' num2str(hp) '_vn.nii.gz'],'f',[1 1 1 1]); - call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_hp' num2str(hp) '_vn.nii.gz -d']); -else - save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_vnf.nii.gz'],'f',[1 1 1 1]); - call_fsl(['fslmaths ' fmri '.nii.gz -sub ' fmri '.nii.gz -add ' fmri '_vnf.nii.gz ' fmri '_vnf.nii.gz']); - %call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_vnf.nii.gz -d']); +%% Apply VN and Save HP_VN TCS +% Here, NIFTI VN'ed TCS gets saved regardless of whether hp>=0 (i.e., regardless of whether individual or concatenated run) +% But CIFTI version of TCS only saved if hp>0 +if dovol > 0 + cts=cts./repmat(Outcts.noise_unst_std,1,ctsT); + % Use '_vnts' (volume normalized time series) as the suffix for the volumetric VN'ed TCS + fname=[fmri hpstring '_vnts.nii.gz']; + ctsfull=zeros(ctsX*ctsY*ctsZ,ctsT, 'single'); + ctsfull(ctsmask,:)=cts; + save_avw(reshape(ctsfull,ctsX,ctsY,ctsZ,ctsT),fname,'f',[1 1 1 1]); + clear ctsfull; + % N.B. Version of 'fslcpgeom' in FSL 6.0.0 requires a patch because it doesn't copy both the qform and sform faithfully + call_fsl(['fslcpgeom ' fmri '.nii.gz ' fname ' -d']); end - +% For CIFTI, we can use the extension to distinguish between VN maps (.dscalar) and VN'ed time series (.dtseries) if hp>=0 BO.cdata=BO.cdata./repmat(OutBO.noise_unst_std,1,size(BO.cdata,2)); - ciftisave(BO,[fmri '_Atlas_hp' num2str(hp) '_vn.dtseries.nii'],WBC); + ciftisave(BO,[fmri '_Atlas' regstring hpstring '_vn.dtseries.nii'],WBC); end %Echo Dims +%TSC: add the regstring to the output filename to avoid overwriting if hp>=0 - dlmwrite([fmri '_dims.txt'],[Outcts.calcDim OutBO.calcDim],'\t'); + if dovol > 0 + dlmwrite([fmri regstring '_dims.txt'],[Outcts.calcDim OutBO.calcDim],'\t'); + else + dlmwrite([fmri regstring '_dims.txt'],[OutBO.calcDim],'\t'); + end else + %TSC: this mode never gets called with a regstring dlmwrite([fmri '_dims.txt'],[Outcts.calcDim],'\t'); end +dlmwrite([fmri '_wf.txt'],[ndhpvol ndhpcifti ndvol],'\t'); end + +%% ---------------------------------------------- + +%% Polynomial detrending function +function Y = detrendpoly(X,p); + +% X: Input data (column major order) +% p: Order of polynomial to remove +% Y: Detrended output + +% Need to define a function to accomplish this, because MATLAB's own DETREND +% is only capable of removing a *linear* trend (i.e., "p=1" only) + + % Check data, must be in column order + [m, n] = size(X); + if (m == 1) + X = X'; + r=n; + else + r=m; + end + + if (~isscalar(p) || p < 0 || p ~= fix(p)) + error('order of polynomial (p) must be a non-negative integer'); + end + + % 5/1/2019 -- Construct the "Vandermonde matrix" (V) scaled to a maximum of 1, for better numerical properties. + % Note that Octave's DETREND function supports arbitrary polynomial orders, but computes V by taking powers + % of [1:r] (rather than [1:r]/r), which is not numerically robust as p increases. + + V = ([1 : r]'/r * ones (1, p + 1)) .^ (ones (r, 1) * [0 : p]); % "Vandermonde" design matrix + + % Cast design matrix to 'single' if the input is also 'single' (which CIFTI will be) + if strcmp(class(X),'single') + V = single(V); + end + + % Use mldivide ('\') as the linear solver, as it has the nice property of generating a warning + % if the solution is rank deficient (in MATLAB at least; Octave doesn't appear to generate a similar warning). + % [In contrast, PINV does NOT generate a warning if the singular values are less than its internal tolerance]. + % Note that even with the scaling of the Vandermonde matrix to a maximum of 1, rank deficiency starts + % becoming a problem at p=6 for data of class 'single' and 1200 time points. + % Rather than explicitly restricting the allowed order here, we'll code a restriction into the calling scripts. + + Y = X - V * (V \ X); % Remove polynomial fit + +end + diff --git a/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m b/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m new file mode 100755 index 0000000..3c0f78b --- /dev/null +++ b/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m @@ -0,0 +1,369 @@ +function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) +% +% FUNCTIONHIGHPASSANDVARIANCENORMALIZE(TR,HP,FMRI,WBC,[REGSTRING],[ndhpvol ndhpcifti ndconcatvol]) +% +% This function: +% (1) Detrends / high-pass filters motion confounds, NIFTI (volume) and CIFTI files +% (2) Calls icaDim.m to estimate the dimensionalities of structured and unstructured noise, +% and compute a spatial variance normalization map. +% It is written to support 'hcp_fix_multi_run', and is not intended as a general (stand-alone) +% function supporting detrending and variance normalization. +% +% TR: repetition time between frames, in sec +% HP: determines what high-pass filtering to apply to the motion confounds and data +% HP<0: No highpass applied (demeaning only). No 'hp' string will be added as part of the file names. +% HP>0: The full-width (2*sigma), in sec, to apply using 'fslmaths -bptf'. +% HP='pd#': If HP is a string, in which the first two characters are 'pd', followed by an integer value, +% then polynomial detrending is applied, with the order specified by the integer value +% embedded at the end of the string. +% The output files will include the string '_hppd#' +% HP=0: Gets interpreted as a linear (1st order) detrend +% This is consistent with the convention in FIX of using HP=0 to specify a linear detrend. +% The output file names will include the string '_hp0'. +% FMRI: The base string of the fmri time series (no extensions) +% WBC: wb_command (full path) +% varargin: +% [REGSTRING]: Additional registration-related string to add to output file names. OPTIONAL. +% [ndhpvol ndhpcifti ndconcatvol]: number of distributions highpassed volume, highpassed cifti and concatenated volume. OPTIONAL. +% [docifti]: Calculate cifti. If you use this option, the [REGSTRING] and [ndhpvol ndhpcifti ndconcatvol] options must be specified. +% If "" is specified for the [REGSTRING] and [ndhpvol ndhpcifti ndconcatvol] options, the default values are used. + +% Note: HP='pd0' would be interpreted as a true 0th order detrend, which is +% the same as demeaning. Mathematically, this is the same as the HP<0 condition, +% but the output files will be named differently (i.e., include '_hppd0'), and additional +% output files will be written relative to the HP<0 condition. + +% Authors: M. Glasser and M. Harms + +CIFTIMatlabReaderWriter=getenv('FSL_FIX_CIFTIRW'); +if (~isdeployed) + % addpath does not work and should not be used in compiled MATLAB + fprintf('Adding %s to MATLAB path\n', CIFTIMatlabReaderWriter); + addpath(CIFTIMatlabReaderWriter); +end + +%% Defaults +dovol = 1; +docifti = 1; +regstring = ''; +pdflag = false; % Polynomial detrend +pdstring = 'pd'; % Expected string at start of HP variable to request a "polynomial detrend" +hpstring = ''; + +%% Parse varargin +if length(varargin) == 1 && ~isempty(varargin{1}) + dovol = 0; %regname is only used for a new surface registration, so will never require redoing volume + regstring = varargin{1}; %this has the underscore on the front already + ndhpvol=2; + ndhpcifti=3; + ndconcatvol=3; + if ~ischar(regstring) + error('%s: REGSTRING should be a string', mfilename); + end +elseif length(varargin) >= 3 + if ~isempty(varargin{4}) + if ~isempty(varargin{1}) + dovol = 0; + end + regstring = varargin{1}; + if ~ischar(regstring) + error('%s: REGSTRING should be a string', mfilename); + end + ndhpvol=varargin{2}; + ndhpcifti=varargin{3}; + ndconcatvol=varargin{4}; + else + ndhpvol=varargin{1}; + ndhpcifti=varargin{2}; + ndconcatvol=varargin{3}; + if isempty(ndvol) + ndvol=2; + end + if isempty(ndvol) + ndhpcifti=3; + end + if isempty(ndconcatvol) + ndvol=3; + end + if ~isempty(varargin{5}) + docifti=varargin{5}; + end + if ~isint(ndhpvol) + error('%s: NDHPVOL should be an integer', mfilename); + end + if ~isint(ndhpcifti) + error('%s: NDHPCIFTI should be an integer', mfilename); + end + if ~isint(ndconcatvol) + error('%s: NDCONCATVOL should be an integer', mfilename); + end + if ~isint(docifti) + error('%s: DOCIFTI should be an integer', mfilename); + end +end + +%UNTITLED4 Summary of this function goes here +% Detailed explanation goes here +% Default value: ndhpvol=2;ndhpcifti=3; ndconcatvol=3; - Takuya Hayashi, Matt Glasser +fprintf('hp=%d, ndist=%d,%d,%d',hp,ndhpvol,ndhpcifti,ndconcatvol) + + +% Check whether polynomial detrend is being requested for the high-pass filtering. +if ischar(hp) + hp = lower(hp); + if strncmp(hp,pdstring,numel(pdstring)) + pdflag = true; + hp = str2double(hp(numel(pdstring)+1:end)); % hp is now a numeric representing the order of the polynomial detrend + elseif ~isempty(str2double(hp)) % Allow for hp to be provided as a string that contains purely numeric elements + hp = str2double(hp); + else error('%s: Invalid specification for the high-pass filter', mfilename); + end +end + +%% Allow for compiled matlab (hp as a string is already handled above) +if (isdeployed) + TR = str2double(TR); +end + +%% Argument checking +if hp==0 && ~pdflag + warning('%s: hp=0 will be interpreted as polynomial detrending of order 1', mfilename); + pdflag=true; + hp=1; + % But in this case, preserve use of "_hp0" in the filename to indicate a linear detrend (for backwards compatibility in file naming) + hpstring = '_hp0'; +end +if pdflag + if (~isscalar(hp) || hp < 0 || hp ~= fix(hp)) + error('%s: Invalid specification for the order of the polynomial detrending', mfilename); + end +end + +%% Finish setting up hpstring to use in file names +if ~pdflag + pdstring = ''; +end +if isempty(hpstring) && hp>=0 + hpstring = ['_hp' pdstring num2str(hp)]; +end + +%% Load the motion confounds, and the CIFTI (if hp>=0) (don't need either if hp<0) +if hp>=0 + confounds=load([fmri hpstring '.ica/mc/prefiltered_func_data_mcf.par']); + confounds=confounds(:,1:6); + %% normalise function is in HCPPIPEDIR/global/matlab/normalise.m + confounds=normalise([confounds [zeros(1,size(confounds,2)); confounds(2:end,:)-confounds(1:end-1,:)] ]); + confounds=normalise([confounds confounds.*confounds]); + if docifti > 0 + BO=ciftiopen([fmri '_Atlas' regstring '.dtseries.nii'],WBC); + end +end + +%% Apply hp filtering of the motion confounds, volume (if requested), and CIFTI +if pdflag % polynomial detrend case + if dovol > 0 + % Save and filter confounds, as a NIFTI, as expected by FIX + save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]); + confounds=detrendpoly(confounds,hp); + save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf_hp'],'f',[1 1 1 TR]); + + % Load volume time series and reduce to just the non-zero voxels (for memory efficiency) + % Note: Use 'range' to identify non-zero voxels (which is very memory efficient) + % rather than 'std' (which requires additional memory equal to the size of the input) + ctsfull=single(read_avw([fmri '.nii.gz'])); + ctsX=size(ctsfull,1); ctsY=size(ctsfull,2); ctsZ=size(ctsfull,3); ctsT=size(ctsfull,4); + ctsfull=reshape(ctsfull,ctsX*ctsY*ctsZ,ctsT); + ctsmask=range(ctsfull, 2) > 0; + fprintf('Non-empty voxels: %d (= %.2f%% of %d)\n', sum(ctsmask), 100*sum(ctsmask)/size(ctsfull,1), size(ctsfull,1)); + cts=ctsfull(ctsmask,:); + clear ctsfull; + + % Polynomial detrend + cts=detrendpoly(cts',hp)'; + + % Write out result, restoring to original size + ctsfull=zeros(ctsX*ctsY*ctsZ,ctsT, 'single'); + ctsfull(ctsmask,:)=cts; + save_avw(reshape(ctsfull,ctsX,ctsY,ctsZ,ctsT),[fmri hpstring '.nii.gz'],'f',[1 1 1 TR]); + clear ctsfull; + % Use -d flag in fslcpgeom (even if not technically necessary) to reduce possibility of mistakes when editing script + call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri hpstring '.nii.gz -d']); + end + if docifti > 0 + BO.cdata=detrendpoly(BO.cdata',hp)'; + ciftisave(BO,[fmri '_Atlas' regstring hpstring '.dtseries.nii'],WBC); + end + +elseif hp>0 % "fslmaths -bptf" based filtering + if dovol > 0 + % Save and filter confounds, as a NIFTI, as expected by FIX + save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]); + call_fsl(sprintf(['fslmaths ' fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf -bptf %f -1 ' fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf_hp'],0.5*hp/TR)); + + % bptf filtering; no masking here, so this is probably memory inefficient + call_fsl(['fslmaths ' fmri '.nii.gz -bptf ' num2str(0.5*hp/TR) ' -1 ' fmri hpstring '.nii.gz']); + call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri hpstring '.nii.gz -d']); + + % Load in result and reduce to just the non-zero voxels (for memory efficiency going forward) + ctsfull=single(read_avw([fmri hpstring '.nii.gz'])); + ctsX=size(ctsfull,1); ctsY=size(ctsfull,2); ctsZ=size(ctsfull,3); ctsT=size(ctsfull,4); + ctsfull=reshape(ctsfull,ctsX*ctsY*ctsZ,ctsT); + ctsmask=range(ctsfull, 2) > 0; + fprintf('Non-empty voxels: %d (= %.2f%% of %d)\n', sum(ctsmask), 100*sum(ctsmask)/size(ctsfull,1), size(ctsfull,1)); + cts=ctsfull(ctsmask,:); + clear ctsfull; + end + + if docifti > 0 + BOdimX=size(BO.cdata,1); BOdimZnew=ceil(BOdimX/100); BOdimT=size(BO.cdata,2); + save_avw(reshape([BO.cdata ; zeros(100*BOdimZnew-BOdimX,BOdimT)],10,10,BOdimZnew,BOdimT),'Atlas','f',[1 1 1 TR]); + call_fsl(sprintf('fslmaths Atlas -bptf %f -1 Atlas',0.5*hp/TR)); + grot=reshape(single(read_avw('Atlas')),100*BOdimZnew,BOdimT); BO.cdata=grot(1:BOdimX,:); clear grot; + ciftisave(BO,[fmri '_Atlas' regstring hpstring '.dtseries.nii'],WBC); + end + + delete('Atlas.nii.gz'); +elseif hp<0 % If no hp filtering, still need to at least demean the volumetric time series + if dovol > 0 + + % Load volume time series and reduce to just the non-zero voxels (for memory efficiency) + ctsfull=single(read_avw([fmri '.nii.gz'])); + ctsX=size(ctsfull,1); ctsY=size(ctsfull,2); ctsZ=size(ctsfull,3); ctsT=size(ctsfull,4); + ctsfull=reshape(ctsfull,ctsX*ctsY*ctsZ,ctsT); + ctsmask=range(ctsfull, 2) > 0; + fprintf('Non-empty voxels: %d (= %.2f%% of %d)\n', sum(ctsmask), 100*sum(ctsmask)/size(ctsfull,1), size(ctsfull,1)); + cts=ctsfull(ctsmask,:); + clear ctsfull; + + cts=demean(cts')'; + end + +end + + +% NOTE: To interpret the logic of the following "Compute, Save, Apply" sections need to know that this function expects that: +% (1) individual run time series will be passed in with hp requested +% (2) concatenated (multi-run) time series will be passed in with hp<0 +%% [This is a bit of hack: would be cleaner if we disentangled the operation of this function from the specific needs of 'hcp_fix_multi_run']. + +%% Compute VN map +if hp>=0 + if dovol > 0 + Outcts=icaDim(cts,0,1,-1,ndhpvol); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, Volume fits two distributions to deal with MNI transform + end + if docifti > 0 + OutBO=icaDim(BO.cdata,0,1,-1,ndhpcifti); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, CIFTI fits three distributions to deal with volume to CIFTI mapping + end +else + if dovol > 0 + Outcts=icaDim(cts,0,1,-1,ndconcatvol); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, Volume fits two distributions to deal with MNI transform + end +end + +%% Save VN map +% (Only need these for the individual runs, which, per above, are expected to be passed in with hp requested). +if hp>=0 + if dovol > 0 + fname=[fmri hpstring '_vn.nii.gz']; + vnfull=zeros(ctsX*ctsY*ctsZ,1, 'single'); + vnfull(ctsmask)=Outcts.noise_unst_std; + + save_avw(reshape(vnfull,ctsX,ctsY,ctsZ,1),fname,'f',[1 1 1 1]); + clear vnfull; + call_fsl(['fslcpgeom ' fmri '_mean.nii.gz ' fname ' -d']); + end + if docifti > 0 + VN=BO; + VN.cdata=OutBO.noise_unst_std; + end + fname=[fmri '_Atlas' regstring hpstring '_vn.dscalar.nii']; + disp(['saving ' fname]); + if docifti > 0 + ciftisavereset(VN,fname,WBC); + end +end + +%% Apply VN and Save HP_VN TCS +% Here, NIFTI VN'ed TCS gets saved regardless of whether hp>=0 (i.e., regardless of whether individual or concatenated run) +% But CIFTI version of TCS only saved if hp>0 +if dovol > 0 + cts=cts./repmat(Outcts.noise_unst_std,1,ctsT); + % Use '_vnts' (volume normalized time series) as the suffix for the volumetric VN'ed TCS + fname=[fmri hpstring '_vnts.nii.gz']; + ctsfull=zeros(ctsX*ctsY*ctsZ,ctsT, 'single'); + ctsfull(ctsmask,:)=cts; + save_avw(reshape(ctsfull,ctsX,ctsY,ctsZ,ctsT),fname,'f',[1 1 1 1]); + clear ctsfull; + % N.B. Version of 'fslcpgeom' in FSL 6.0.0 requires a patch because it doesn't copy both the qform and sform faithfully + call_fsl(['fslcpgeom ' fmri '.nii.gz ' fname ' -d']); +end +% For CIFTI, we can use the extension to distinguish between VN maps (.dscalar) and VN'ed time series (.dtseries) +if hp>=0 && docifti + BO.cdata=BO.cdata./repmat(OutBO.noise_unst_std,1,size(BO.cdata,2)); + ciftisave(BO,[fmri '_Atlas' regstring hpstring '_vn.dtseries.nii'],WBC); +end + +%Echo Dims +%TSC: add the regstring to the output filename to avoid overwriting +if hp>=0 && docifti > 0 + if dovol > 0 + dlmwrite([fmri regstring '_dims.txt'],[Outcts.calcDim OutBO.calcDim],'\t'); + else + dlmwrite([fmri regstring '_dims.txt'],[OutBO.calcDim],'\t'); + end +else + %TSC: this mode never gets called with a regstring + dlmwrite([fmri '_dims.txt'],[Outcts.calcDim],'\t'); +end + +dlmwrite([fmri '_wf.txt'],[ndhpvol ndhpcifti ndconcatvol],'\t'); +end + + +%% ---------------------------------------------- + +%% Polynomial detrending function +function Y = detrendpoly(X,p); + +% X: Input data (column major order) +% p: Order of polynomial to remove +% Y: Detrended output + +% Need to define a function to accomplish this, because MATLAB's own DETREND +% is only capable of removing a *linear* trend (i.e., "p=1" only) + + % Check data, must be in column order + [m, n] = size(X); + if (m == 1) + X = X'; + r=n; + else + r=m; + end + + if (~isscalar(p) || p < 0 || p ~= fix(p)) + error('order of polynomial (p) must be a non-negative integer'); + end + + % 5/1/2019 -- Construct the "Vandermonde matrix" (V) scaled to a maximum of 1, for better numerical properties. + % Note that Octave's DETREND function supports arbitrary polynomial orders, but computes V by taking powers + % of [1:r] (rather than [1:r]/r), which is not numerically robust as p increases. + + V = ([1 : r]'/r * ones (1, p + 1)) .^ (ones (r, 1) * [0 : p]); % "Vandermonde" design matrix + + % Cast design matrix to 'single' if the input is also 'single' (which CIFTI will be) + if strcmp(class(X),'single') + V = single(V); + end + + % Use mldivide ('\') as the linear solver, as it has the nice property of generating a warning + % if the solution is rank deficient (in MATLAB at least; Octave doesn't appear to generate a similar warning). + % [In contrast, PINV does NOT generate a warning if the singular values are less than its internal tolerance]. + % Note that even with the scaling of the Vandermonde matrix to a maximum of 1, rank deficiency starts + % becoming a problem at p=6 for data of class 'single' and 1200 time points. + % Rather than explicitly restricting the allowed order here, we'll code a restriction into the calling scripts. + + Y = X - V * (V \ X); % Remove polynomial fit + +end + diff --git a/ICAFIX/hcp_fix_multi_run_RIKEN_merge b/ICAFIX/hcp_fix_multi_run_RIKEN_merge index 15666b9..afdae20 100755 --- a/ICAFIX/hcp_fix_multi_run_RIKEN_merge +++ b/ICAFIX/hcp_fix_multi_run_RIKEN_merge @@ -79,17 +79,10 @@ e.g. hcp_fix_multi_run rfMRI_REST1_RL/rfMRI_REST1_RL.nii.gz@rfMRI_REST1_RL/rfMR (if launching the script from the '\${StudyFolder}/\${Subject}/MNINonLinear/Results' directory) Optional single character arguments: - -t : Training file (default is HCP_Style_Single_Multirun_Dedrift.RData for Human, NHPHCP_Macaque.RData for Macaque, - NHPHCP_Marmoset.RData for Marmoset) -d : Dimension for group-ICA (default is estimated by icaDIM) - -s : Species - 0: Human (default), 1: Macaque, 2: Marmoset - -r : Run mode - 0: Run all (PreICA+FIX, ICA+FIX, PostICA+FIX, default), 1: Run ICA,FIX and PostICA+FIX, - 2: Run FIX and PostICA+FIX, 3: Run PostICA+FIX) - -g : use surface geometry FIX -w ,, : Wishart filter Ndist (default is 2,3,3) - -T : FIX threshold 0 to 100 (%) (default is 10 in HUman, 40 in Macaque and Marmoset) -h : help message - -n : Do not overwrite hp if already calculated + -b : Brain mask EOF } @@ -277,34 +270,6 @@ function interpret_as_bool() ############################################################# - -# RData -if [ $Species = 0 ] ; then - if [ "$TrainingData" = "" ] ; then - TrainingData=${FSL_FIXDIR}/training_files/HCP_Style_Single_Multirun_Dedrift.RData; - fi - export SPECIES=Human - if [ "$FixThresh" = "" ] ; then - FixThresh=10 - fi -elif [ $Species = 1 ] ; then - export SPECIES=Macaque - if [ "$TrainingData" = "" ] ; then - TrainingData=${FSL_FIXDIR}/training_files/NHPHCP_Macaque.RData - fi - if [ "$FixThresh" = "" ] ; then - FixThresh=40 - fi -elif [ $Species = 2 ] ; then - export SPECIES=Marmoset - if [ "$TrainingData" = "" ] ; then - TrainingData=${FSL_FIXDIR}/training_files/NHPHCP_Marmoset.RData - fi - if [ "$FixThresh" = "" ] ; then - FixThresh=40 - fi -fi - ## --------------------------------------------------------------------------- ## Check whether FSL version is recent enough ## --------------------------------------------------------------------------- @@ -363,23 +328,17 @@ fi # Inserted by Takuya Hayashi, July 2018 shift;shift;shift;shift Species=0 -RunMode=0 WF="2,3,3" -OW="1" FixThresh="" SURFGEOM=0 +TemplateMask="" UserDim=NONE; while getopts t:d:s:r:w:nT:hg OPT do case "$OPT" in - "t" ) TrainingData="$OPTARG";; "d" ) UserDim="$OPTARG";; - "s" ) Species="$OPTARG";; - "r" ) RunMode="$OPTARG";; "w" ) WF="$OPTARG";; - "n" ) OW=0;; - "T" ) FixThresh="$OPTARG";; - "g" ) SURFGEOM=1;; + "b" ) TemplateMask="$OPTARG";; "h" ) usage_exit;; * ) usage_exit;; esac @@ -403,7 +362,6 @@ log_Msg "doMotionRegression: ${doMotionRegression}" log_Msg "TrainingData: ${TrainingData}" log_Msg "FixThresh: ${FixThresh}" log_Msg "DeleteIntermediates: ${DeleteIntermediates}" -log_Msg "SPECIES: ${SPECIES}" log_Msg "WF Ndist: ${WF}" log_Msg "UserDim: ${UserDim}" log_Msg "CARET7DIR: ${CARET7DIR}" @@ -432,179 +390,169 @@ VNCIFTISTRING="" MovementNIFTIMergeSTRING="" MovementNIFTIhpMergeSTRING="" MovementTXTMergeSTRING="" -FileListSTRING="" # Added by TH May 2019 for fmri in $fmris ; do log_Msg "Top of loop through fmris: fmri: ${fmri}" - NIFTIvolMergeSTRING=`echo "${NIFTIvolMergeSTRING}$($FSLDIR/bin/remove_ext $fmri)_demean "` - NIFTIvolhpVNMergeSTRING=`echo "${NIFTIvolhpVNMergeSTRING}$($FSLDIR/bin/remove_ext $fmri)_hp${hp}_vnts "` - SBRefVolSTRING=`echo "${SBRefVolSTRING}$($FSLDIR/bin/remove_ext $fmri)_SBRef "` - MeanVolSTRING=`echo "${MeanVolSTRING}$($FSLDIR/bin/remove_ext $fmri)_mean "` - VNVolSTRING=`echo "${VNVolSTRING}$($FSLDIR/bin/remove_ext $fmri)_hp${hp}_vn "` - CIFTIMergeSTRING=`echo "${CIFTIMergeSTRING} -cifti $($FSLDIR/bin/remove_ext $fmri)_Atlas_demean.dtseries.nii"` - CIFTIhpVNMergeSTRING=`echo "${CIFTIhpVNMergeSTRING} -cifti $($FSLDIR/bin/remove_ext $fmri)_Atlas_hp${hp}_vn.dtseries.nii"` - MeanCIFTISTRING=`echo "${MeanCIFTISTRING} -cifti $($FSLDIR/bin/remove_ext $fmri)_Atlas_mean.dscalar.nii "` - VNCIFTISTRING=`echo "${VNCIFTISTRING} -cifti $($FSLDIR/bin/remove_ext $fmri)_Atlas_hp${hp}_vn.dscalar.nii "` - MovementNIFTIMergeSTRING=`echo "${MovementNIFTIMergeSTRING}$($FSLDIR/bin/remove_ext $fmri)_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf.nii.gz "` - MovementNIFTIhpMergeSTRING=`echo "${MovementNIFTIhpMergeSTRING}$($FSLDIR/bin/remove_ext $fmri)_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf_hp.nii.gz "` + fmriNoExt=$($FSLDIR/bin/remove_ext $fmri) # $fmriNoExt still includes leading directory components + + # Create necessary strings for merging across runs + # N.B. Some of these files don't exist yet, and are about to get created + NIFTIvolMergeSTRING+="${fmriNoExt}_demean " + NIFTIvolhpVNMergeSTRING+="${fmriNoExt}_hp${hp}_vnts " #These are the individual run, VN'ed *time series* + SBRefVolSTRING+="${fmriNoExt}_SBRef " + MeanVolSTRING+="${fmriNoExt}_mean " + VNVolSTRING+="${fmriNoExt}_hp${hp}_vn " #These are the individual run, VN'ed NIFTI *maps* (created by functionhighpassandvariancenormalize) + CIFTIMergeSTRING+="-cifti ${fmriNoExt}_Atlas_demean.dtseries.nii " + CIFTIhpVNMergeSTRING+="-cifti ${fmriNoExt}_Atlas_hp${hp}_vn.dtseries.nii " + MeanCIFTISTRING+="-cifti ${fmriNoExt}_Atlas_mean.dscalar.nii " + VNCIFTISTRING+="-cifti ${fmriNoExt}_Atlas_hp${hp}_vn.dscalar.nii " #These are the individual run, VN'ed CIFTI *maps* (created by functionhighpassandvariancenormalize) + MovementNIFTIMergeSTRING+="${fmriNoExt}_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf " + MovementNIFTIhpMergeSTRING+="${fmriNoExt}_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf_hp " cd `dirname $fmri` - log_Msg "pwd: "$(pwd) - fmri=`basename $fmri` - fmri=`$FSLDIR/bin/imglob $fmri` - [ `imtest $fmri` != 1 ] && echo No valid 4D_FMRI input file specified && exit 1 - fmri_orig=$fmri - FileListSTRING="${FileListSTRING}${fmri}@" # Added by TH May 2019 - - mkdir -p ${fmri}_hp$hp.ica/mc + log_Debug_Msg "pwd: "$(pwd) + fmri=`basename $fmri` # After this, $fmri no longer includes the leading directory components + fmri=`$FSLDIR/bin/imglob $fmri` # After this, $fmri will no longer have an extension (if there was one initially) + if [ `$FSLDIR/bin/imtest $fmri` != 1 ]; then + log_Err_Abort "Invalid 4D_FMRI input file specified: ${fmri}" + fi + + # Get Movement_Regressors.txt into the format expected by functionhighpassandvariancenormalize.m + mkdir -p ${fmri}_hp${hp}.ica/mc if [ -f Movement_Regressors.txt ] ; then - log_Msg "About to create ${fmri}_hp$hp.ica/mc/prefiltered_func_data_mcf.par file" - cat Movement_Regressors.txt | awk '{ print $4 " " $5 " " $6 " " $1 " " $2 " " $3}' > ${fmri}_hp$hp.ica/mc/prefiltered_func_data_mcf.par + log_Debug_Msg "About to create ${fmri}_hp${hp}.ica/mc/prefiltered_func_data_mcf.par file" + cat Movement_Regressors.txt | awk '{ print $4 " " $5 " " $6 " " $1 " " $2 " " $3}' > ${fmri}_hp${hp}.ica/mc/prefiltered_func_data_mcf.par else - log_Error "Movement_Regressors.txt not retrieved properly." - exit -1 + log_Err_Abort "Movement_Regressors.txt not retrieved properly." fi - tr=`$FSLDIR/bin/fslval $fmri pixdim4` - log_Msg "tr: $tr" - log_Msg "processing FMRI file $fmri with highpass $hp" + #Demean movement regressors + demeanMovementRegressors Movement_Regressors.txt Movement_Regressors_demean.txt + MovementTXTMergeSTRING+="$(pwd)/Movement_Regressors_demean.txt " - if [ $OW = 1 ] ; then imrm ${fmri}_demean;fi - if [ ! -e ${fmri}_demean.nii.gz ] ; then - log_Msg "processing FMRI mean and demean" - #Demean volume, movement regressors, and data + #Demean volumes ${FSLDIR}/bin/fslmaths $fmri -Tmean ${fmri}_mean ${FSLDIR}/bin/fslmaths $fmri -sub ${fmri}_mean ${fmri}_demean - fi - - demeanMovementRegressors Movement_Regressors.txt Movement_Regressors_demean.txt - MovementTXTMergeSTRING=`echo "${MovementTXTMergeSTRING}$(pwd)/Movement_Regressors_demean.txt "` - - if [ $OW = 1 ] ; then rm -rf $($FSLDIR/bin/remove_ext $fmri)_Atlas_demean.dtseries.nii ; fi - if [ ! -e $($FSLDIR/bin/remove_ext $fmri)_Atlas_demean.dtseries.nii ] ; then - log_Msg "Calculating mean and demeaned fmri Atlas" - ${FSL_FIX_WBC} -cifti-reduce $($FSLDIR/bin/remove_ext $fmri)_Atlas.dtseries.nii MEAN $($FSLDIR/bin/remove_ext $fmri)_Atlas_mean.dscalar.nii - ${FSL_FIX_WBC} -cifti-math "TCS - MEAN" $($FSLDIR/bin/remove_ext $fmri)_Atlas_demean.dtseries.nii -var TCS $($FSLDIR/bin/remove_ext $fmri)_Atlas.dtseries.nii -var MEAN $($FSLDIR/bin/remove_ext $fmri)_Atlas_mean.dscalar.nii -select 1 1 -repeat - fi - #Highpass Filter and Variance Normalize Data - #if [ $hp -gt 0 ] ; then - # log_Msg "running highpass" - # hptr=`echo "10 k $hp 2 / $tr / p" | dc -` - # log_Msg "hptr $hptr" - # ${FSLDIR}/bin/fslmaths $fmri -sub ${fmri}_mean -bptf $hptr -1 ${fmri}_hp$hp - #fi - - #log_Msg "functionmotionconfounds log file is to be named: .fix.functionmotionconfounds.log instead of .fix.log" - #cd ${fmri}_hp$hp.ica - #${FSL_FIXDIR}/call_matlab.sh -l .fix.functionmotionconfounds.log -f functionmotionconfounds $tr $hp - #cd .. - - #${FSL_FIX_WBC} -cifti-convert -to-nifti ${fmri}_Atlas.dtseries.nii ${fmri}_Atlas_FAKENIFTI.nii.gz - #${FSLDIR}/bin/fslmaths ${fmri}_Atlas_FAKENIFTI.nii.gz -bptf $hptr -1 ${fmri}_Atlas_hp$hp_FAKENIFTI.nii.gz - #${FSL_FIX_WBC} -cifti-convert -from-nifti ${fmri}_Atlas_hp$hp_FAKENIFTI.nii.gz ${fmri}_Atlas.dtseries.nii ${fmri}_Atlas_hp$hp.dtseries.nii - #$FSLDIR/bin/imrm ${fmri}_Atlas_FAKENIFTI ${fmri}_Atlas_hp$hp_FAKENIFTI + #Demean CIFTI + ${Caret7_Command} -cifti-reduce ${fmri}_Atlas.dtseries.nii MEAN ${fmri}_Atlas_mean.dscalar.nii + ${Caret7_Command} -cifti-math "TCS - MEAN" ${fmri}_Atlas_demean.dtseries.nii -var TCS ${fmri}_Atlas.dtseries.nii -var MEAN ${fmri}_Atlas_mean.dscalar.nii -select 1 1 -repeat + + ## "1st pass" VN on the individual runs; high-pass gets done here as well + tr=`$FSLDIR/bin/fslval $fmri pixdim4` #No checking currently that TR is same across runs + log_Msg "tr: ${tr}" + log_Msg "processing FMRI file ${fmri} with highpass ${hp}" - #cd ${fmri}_hp$hp.ica + ## MPH: Log output to existing stdout/stderr, rather than a separate .fix.functionhighpassandvariancenormalize.log file #if [ -e .fix.functionhighpassandvariancenormalize.log ] ; then - # rm .fix.functionhighpassandvariancenormalize.log + # /bin/rm .fix.functionhighpassandvariancenormalize.log #fi - if [ -e .fix.functionhighpassandvariancenormalize_riken.log ] ; then - rm .fix.functionhighpassandvariancenormalize_riken.log - fi - if [ $OW = 1 ] ; then - if [ -e ${fmri}_hp${hp}.nii.gz ] ; then - log_Msg "Found and removing ${fmri}_hp${hp}" - imrm ${fmri}_hp${hp} - fi - if [ -e ${fmri}_Atlas_hp${hp}.dtseries.nii ] ; then - log_Msg "Found and removing ${fmri}_Atlas_hp${hp}.dtseries.nii" - rm -rf ${fmri}_Atlas_hp${hp}.dtseries.nii - fi - fi - if [ -e ${fmri}_hp${hp}.nii.gz ] ; then - log_Msg "Found ${fmri}_hp${hp} and confirming it is already filtered" - v=`fslstats ${fmri}_hp${hp} -M` - if [ `echo "$v > 0.000001" | bc` -eq 1 ] || [ `echo "$v < -0.000001" | bc` -eq 1 ] ; then - log_Msg "Found ${fmri}_hp${hp} is not yet filtered" - imrm ${fmri}_hp${hp} + ## We will use the call_matlab.sh script available in ${FSL_FIXDIR}. + ## But the arguments need to be customized depending on the matlab mode, because: + ## (1) the -r argument does not support the use of compiled matlab, but is needed for interpreted matlab/octave + ## to allow addition of 'addpath' as a command, so that functionhighpassandvariancenormalize.m can be found within matlab/octave. + ## (2) for compiled matlab, we need to use the -c and -b flags to specify the location of the MCR and the compiled binary. + ## (We cannot use the defaults in ${FSL_FIXDIR}/settings.sh here, since functionhighpassandvariancenormalize.m is an HCP + ## specific script that is not part of the FIX distribution). + + # For INTERPRETED MODES, make sure that matlab/octave has access to the functions it needs. + # normalise.m (needed by functionhighpassandvariancenormalize.m) is in '${HCPPIPEDIR}/global/matlab' + # Note that the call_matlab.sh script itself adds ${FSL_FIXDIR} (all the fix-related functions) + # and ${FSL_MATLAB_PATH} (e.g., read_avw.m, save_avw.m) to the matlab path. + # Note that fix_3_clean.m *appends* ${FSL_FIX_CIFTIRW} (defined in ${FSL_FIXDIR}/settings.sh) + # to the matlab path (i.e., the ciftiopen.m, ciftisave.m functions). + # BUT the CIFTI I/O functions are also now included within '${HCPPIPEDIR}/global/matlab' as well. + # AND since 'addpath' *prepends* to the matlab path, the versions in '${HCPPIPEDIR}/global/matlab' will therefore take precedence! + ML_PATHS="addpath('${HCPPIPEDIR}/global/matlab'); addpath('${this_script_dir}/scripts');" + + case ${FSL_FIX_MATLAB_MODE} in + 0) + # Use compiled Matlab + if [ -z "${MATLAB_COMPILER_RUNTIME}" ]; then + log_Err_Abort "To use MATLAB run mode: ${FSL_FIX_MATLAB_MODE}, the MATLAB_COMPILER_RUNTIME environment variable must be set" else - log_Msg "Found ${fmri}_hp${hp} is already filtered" + log_Msg "MATLAB_COMPILER_RUNTIME: ${MATLAB_COMPILER_RUNTIME}" fi - fi - if [ -e ${fmri}_Atlas_hp${hp}.dtseries.nii ] ; then - log_Msg "Found and using ${fmri}_Atlas_hp${hp}.dtseries.nii" - fi - - log_Msg "${FSL_FIXDIR}/call_matlab.sh -l .fix.functionhighpassandvariancenormalize_riken.log -f functionhighpassandvariancenormalize_riken $tr $hp $fmri ${FSL_FIX_WBC} $ndhpvol $ndhpcifti $ndcifti" - #${FSL_FIXDIR}/call_matlab.sh -l .fix.functionhighpassandvariancenormalize_riken.log -f functionhighpassandvariancenormalize_riken $tr $hp $fmri ${FSL_FIX_WBC} - ${FSL_FIXDIR}/call_matlab.sh -l .fix.functionhighpassandvariancenormalize_riken.log -f functionhighpassandvariancenormalize_riken $tr $hp $fmri ${FSL_FIX_WBC} $ndhpvol $ndhpcifti $ndcifti + CSUBDIR=${HCPPIPEDIR}/ICAFIX/scripts/Compiled_functionhighpassandvariancenormalize + matlab_cmd="${FSL_FIXDIR}/call_matlab.sh -c ${MATLAB_COMPILER_RUNTIME} -b ${CSUBDIR} -f functionhighpassandvariancenormalize ${tr} ${hp} ${fmri} ${Caret7_Command}" + log_Msg "Run compiled MATLAB with command..." + log_Msg "${matlab_cmd}" + ${matlab_cmd} + ;; + + 1 | 2) + # Use interpreted MATLAB or Octave + # ${hp} needs to be passed in as a string, to handle the hp=pd* case + matlab_cmd="${ML_PATHS} functionhighpassandvariancenormalize(${tr}, '${hp}', '${fmri}', '${Caret7_Command}', '${ndhpvol}, '${ndhpcifti}', '${ndcifti});" + log_Msg "Run interpreted MATLAB/Octave with command..." + log_Msg "${FSL_FIXDIR}/call_matlab.sh -r ${matlab_cmd}" + ${FSL_FIXDIR}/call_matlab.sh -r "${matlab_cmd}" # Note: the semicolon within ${matlab_cmd} needs to be protected with quotation marks + ;; + *) + # Unsupported MATLAB run mode + log_Err_Abort "Unsupported MATLAB run mode value (FSL_FIX_MATLAB_MODE) in $FSL_FIXDIR/settings.sh: ${FSL_FIX_MATLAB_MODE}" + ;; + esac + + log_Msg "Dims: $(cat ${fmri}_dims.txt)" + + # Demean the movement regressors (in the 'fake-NIFTI' format returned by functionhighpassandvariancenormalize) + fslmaths ${fmri}_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf -Tmean ${fmri}_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf_mean + fslmaths ${fmri}_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf -sub ${fmri}_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf_mean ${fmri}_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf + $FSLDIR/bin/imrm ${fmri}_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf_mean - #log_Msg "matlab -nodesktop -nosplash -r \"addpath('$FSL_FIXDIR'); addpath('${FSLDIR}/etc/matlab'); functionhighpassandvariancenormalize_riken($tr,$hp,'$fmri','$FSL_FIX_WBC',$ndhpvol,$ndhpcifti,$ndcifti);\"" - #matlab -nodesktop -nosplash -r "addpath('$FSL_FIXDIR'); addpath('${FSLDIR}/etc/matlab'); functionhighpassandvariancenormalize_riken($tr,$hp,'$fmri','$FSL_FIX_WBC',$ndhpvol,$ndhpcifti,$ndcifti);" >> .fix.functionhighpassandvariancenormalize_riken.log 2>&1 + cd ${DIR} # Return to directory where script was launched - #cd .. - if [ "$(cat ${fmri}_dims.txt)" != "" ] ; then - log_Msg "Dims: $(cat ${fmri}_dims.txt)" - else - log_Msg "Dims was not properly estimated by IcaDim" - fi - - fslmaths $(pwd)/${fmri}_hp$hp.ica/mc/prefiltered_func_data_mcf_conf.nii.gz -Tmean $(pwd)/${fmri}_hp$hp.ica/mc/prefiltered_func_data_mcf_conf_mean.nii.gz - fslmaths $(pwd)/${fmri}_hp$hp.ica/mc/prefiltered_func_data_mcf_conf.nii.gz -sub $(pwd)/${fmri}_hp$hp.ica/mc/prefiltered_func_data_mcf_conf_mean.nii.gz $(pwd)/${fmri}_hp$hp.ica/mc/prefiltered_func_data_mcf_conf.nii.gz - $FSLDIR/bin/imrm $(pwd)/${fmri}_hp$hp.ica/mc/prefiltered_func_data_mcf_conf_mean.nii.gz - - fmri=${fmri}_hp$hp - #cd ${fmri}.ica - # Per https://github.com/Washington-University/Pipelines/issues/60, the following line doesn't appear to be necessary - #$FSLDIR/bin/imln ../$fmri filtered_func_data - #cd .. log_Msg "Bottom of loop through fmris: fmri: ${fmri}" -done -#AlreadyHP="-1" #Don't run highpass on concatenated data +done ### END LOOP (for fmri in $fmris; do) + +## --------------------------------------------------------------------------- +## Concatenate the individual runs and create necessary files +## --------------------------------------------------------------------------- #Make Concatenated Folder ConcatFolder=`dirname ${ConcatName}` log_Msg "Making concatenated folder: ${ConcatFolder}" if [ ! -e ${ConcatFolder} ] ; then - mkdir -p ${ConcatFolder} + mkdir ${ConcatFolder} else - if [ $OW = 1 ] ; then - rm -r ${ConcatFolder} - fi - mkdir -p ${ConcatFolder} + /bin/rm -r ${ConcatFolder} + log_Warn "Previous ${ConcatFolder} directory removed" + mkdir ${ConcatFolder} fi ConcatNameNoExt=$($FSLDIR/bin/remove_ext $ConcatName) # No extension, but still includes the directory path -if [ -e `remove_ext ${ConcatName}`_fmris.txt ] ; then rm `remove_ext ${ConcatName}`_fmris.txt;fi # Add by TH May 2019 -echo ${FileListSTRING} >> `remove_ext ${ConcatName}`_fmris.txt # Add by TH May 2019 - -fslmerge -tr `remove_ext ${ConcatName}`_demean ${NIFTIvolMergeSTRING} $tr -fslmerge -tr `remove_ext ${ConcatName}`_hp${hp}_vnts ${NIFTIvolhpVNMergeSTRING} $tr -fslmerge -t `remove_ext ${ConcatName}`_SBRef ${SBRefVolSTRING} -fslmerge -t `remove_ext ${ConcatName}`_mean ${MeanVolSTRING} -fslmerge -t `remove_ext ${ConcatName}`_hp${hp}_vn ${VNVolSTRING} -fslmaths `remove_ext ${ConcatName}`_SBRef -Tmean `remove_ext ${ConcatName}`_SBRef -fslmaths `remove_ext ${ConcatName}`_mean -Tmean `remove_ext ${ConcatName}`_mean -fslmaths `remove_ext ${ConcatName}`_hp${hp}_vn -Tmean `remove_ext ${ConcatName}`_hp${hp}_vn -fslmaths `remove_ext ${ConcatName}`_hp${hp}_vnts -mul `remove_ext ${ConcatName}`_hp${hp}_vn `remove_ext ${ConcatName}`_hp${hp} -fslmaths `remove_ext ${ConcatName}`_demean -add `remove_ext ${ConcatName}`_mean `remove_ext ${ConcatName}` -fslmaths `remove_ext ${ConcatName}`_SBRef -bin `remove_ext ${ConcatName}`_brain_mask # Inserted to create mask to be used in melodic for suppressing memory error - Takuya Hayashi - +# Merge volumes from the individual runs +fslmerge -tr ${ConcatNameNoExt}_demean ${NIFTIvolMergeSTRING} $tr +fslmerge -tr ${ConcatNameNoExt}_hp${hp}_vnts ${NIFTIvolhpVNMergeSTRING} $tr +fslmerge -t ${ConcatNameNoExt}_SBRef ${SBRefVolSTRING} +fslmerge -t ${ConcatNameNoExt}_mean ${MeanVolSTRING} +fslmerge -t ${ConcatNameNoExt}_hp${hp}_vn ${VNVolSTRING} +# Average across runs +fslmaths ${ConcatNameNoExt}_SBRef -Tmean ${ConcatNameNoExt}_SBRef +fslmaths ${ConcatNameNoExt}_mean -Tmean ${ConcatNameNoExt}_mean # "Grand" mean across runs +fslmaths ${ConcatNameNoExt}_demean -add ${ConcatNameNoExt}_mean ${ConcatNameNoExt} + # Preceding line adds back in the "grand" mean + # Resulting file not used below, but want this concatenated version (without HP or VN) to exist +fslmaths ${ConcatNameNoExt}_hp${hp}_vn -Tmean ${ConcatNameNoExt}_hp${hp}_vn # Mean VN map across the individual runs +fslmaths ${ConcatNameNoExt}_hp${hp}_vnts -mul ${ConcatNameNoExt}_hp${hp}_vn ${ConcatNameNoExt}_hp${hp} + # Preceding line restores the mean VN map + # The resulting TCS becomes the input to a 2nd pass of VN, which then becomes the input to melodic +fslmaths ${ConcatNameNoExt}_SBRef -bin ${ConcatNameNoExt}_brain_mask + # Preceding line creates mask to be used in melodic for suppressing memory error - Takuya Hayashi # Same thing for the CIFTI -${FSL_FIX_WBC} -cifti-merge `remove_ext ${ConcatName}`_Atlas_demean.dtseries.nii ${CIFTIMergeSTRING} -${FSL_FIX_WBC} -cifti-average `remove_ext ${ConcatName}`_Atlas_mean.dscalar.nii ${MeanCIFTISTRING} -${FSL_FIX_WBC} -cifti-average `remove_ext ${ConcatName}`_Atlas_hp${hp}_vn.dscalar.nii ${VNCIFTISTRING} -${FSL_FIX_WBC} -cifti-math "TCS + MEAN" `remove_ext ${ConcatName}`_Atlas.dtseries.nii -var TCS `remove_ext ${ConcatName}`_Atlas_demean.dtseries.nii -var MEAN `remove_ext ${ConcatName}`_Atlas_mean.dscalar.nii -select 1 1 -repeat -${FSL_FIX_WBC} -cifti-merge `remove_ext ${ConcatName}`_Atlas_hp${hp}_vn.dtseries.nii ${CIFTIhpVNMergeSTRING} -${FSL_FIX_WBC} -cifti-math "TCS * VN" `remove_ext ${ConcatName}`_Atlas_hp${hp}.dtseries.nii -var TCS `remove_ext ${ConcatName}`_Atlas_hp${hp}_vn.dtseries.nii -var VN `remove_ext ${ConcatName}`_Atlas_hp${hp}_vn.dscalar.nii -select 1 1 -repeat - +${Caret7_Command} -cifti-merge ${ConcatNameNoExt}_Atlas_demean.dtseries.nii ${CIFTIMergeSTRING} +${Caret7_Command} -cifti-average ${ConcatNameNoExt}_Atlas_mean.dscalar.nii ${MeanCIFTISTRING} +${Caret7_Command} -cifti-math "TCS + MEAN" ${ConcatNameNoExt}_Atlas.dtseries.nii -var TCS ${ConcatNameNoExt}_Atlas_demean.dtseries.nii -var MEAN ${ConcatNameNoExt}_Atlas_mean.dscalar.nii -select 1 1 -repeat +${Caret7_Command} -cifti-merge ${ConcatNameNoExt}_Atlas_hp${hp}_vn.dtseries.nii ${CIFTIhpVNMergeSTRING} +${Caret7_Command} -cifti-average ${ConcatNameNoExt}_Atlas_hp${hp}_vn.dscalar.nii ${VNCIFTISTRING} +${Caret7_Command} -cifti-math "TCS * VN" ${ConcatNameNoExt}_Atlas_hp${hp}.dtseries.nii -var TCS ${ConcatNameNoExt}_Atlas_hp${hp}_vn.dtseries.nii -var VN ${ConcatNameNoExt}_Atlas_hp${hp}_vn.dscalar.nii -select 1 1 -repeat # At this point the concatenated VN'ed time series (both volume and CIFTI, following the "1st pass" VN) can be deleted log_Msg "Removing the concatenated VN'ed time series" @@ -638,102 +586,120 @@ done ## Prepare for melodic on concatenated file ## --------------------------------------------------------------------------- -#ConcatFolder=`dirname ${ConcatName}` -cd ${ConcatFolder} - -#Concat fMRI is intensity normalized (not variance normalized) and has no overall mean but is detrended and within run demeaned and variance normalized for both volume and CIFTI -#concatfmri=`basename $(remove_ext ${ConcatName})`_hp${hp} -#concat_fmri_orig=`basename $(remove_ext ${ConcatName})` - -mkdir -p ${concatfmri}.ica -log_Msg "About to put contents of ${MovementTXTMergeSTRING} in Movement_Regressors_demean.txt file" -cat ${MovementTXTMergeSTRING} > Movement_Regressors_demean.txt -mkdir -p ${concatfmri}.ica/mc -fslmerge -tr ${concatfmri}.ica/mc/prefiltered_func_data_mcf_conf_hp ${MovementNIFTIhpMergeSTRING} $tr -fslmerge -tr ${concatfmri}.ica/mc/prefiltered_func_data_mcf_conf ${MovementNIFTIMergeSTRING} $tr -# Added the following line - needed for FIX feature extraction - Takuya Hayashi July 2018 -cat Movement_Regressors_demean.txt | awk '{ print $4 " " $5 " " $6 " " $1 " " $2 " " $3}' > ${concatfmri}.ica/mc/prefiltered_func_data_mcf.par - -cd `dirname $fmri` -} - -RunICA () { -AlreadyHP="-1" #Don't run highpass on concatenated data +# concatfmri is intensity normalized and has no overall mean. +# It is also not variance normalized (the mean VN map *across runs* was "restored" above). +# But is detrended and within run demeaned and variance normalized for both volume and CIFTI +concatfmri=`basename ${ConcatNameNoExt}` # Directory path is now removed +concat_fmri_orig=`basename $(remove_ext ${ConcatName})` +concatfmrihp=${concatfmri}_hp${hp} + +#MPH: Need to be in original directory here, for the ${MovementNIFTI*STRING} variables to work +#if user did not use absolute paths for input arguments +mkdir -p ${ConcatFolder}/${concatfmrihp}.ica +log_Debug_Msg "About to put contents of ${MovementTXTMergeSTRING} in Movement_Regressors_demean.txt file" +cat ${MovementTXTMergeSTRING} > ${ConcatFolder}/Movement_Regressors_demean.txt +mkdir ${ConcatFolder}/${concatfmrihp}.ica/mc +fslmerge -tr ${ConcatFolder}/${concatfmrihp}.ica/mc/prefiltered_func_data_mcf_conf_hp ${MovementNIFTIhpMergeSTRING} $tr +fslmerge -tr ${ConcatFolder}/${concatfmrihp}.ica/mc/prefiltered_func_data_mcf_conf ${MovementNIFTIMergeSTRING} $tr + +#MPH: Now we can switch to ${ConcatFolder} cd ${ConcatFolder} -if [ "$tr" = "" ] ; then - fmri=`echo ${fmris} | tr ' ' '\n' | tail -1` - tr=`$FSLDIR/bin/fslval $fmri pixdim4` - log_Msg "tr: $tr" -fi - -log_Msg "running MELODIC" -log_Msg "About to run melodic: Contents of ${concatfmri}.ica follow" -if [ "${DEBUG}" = "TRUE" ] ; then - ls -lRa `remove_ext ${concatfmri}`.ica -fi - -log_Msg "Running melodic located at: ${FSLDIR}/bin/melodic" -log_Msg "Beginning of melodic version log, help, and checksum" -if [ "${DEBUG}" = "TRUE" ] ; then - log_Msg "${FSLDIR}/bin/melodic --version" - ${FSLDIR}/bin/melodic --version - log_Msg "${FSLDIR}/bin/melodic --help" - ${FSLDIR}/bin/melodic --help - log_Msg "md5sum ${FSLDIR}/bin/melodic" - md5sum ${FSLDIR}/bin/melodic -fi -log_Msg "End of melodic version log, help, and checksum" -###Run on concatenated file -- +## "2nd pass" VN +## This is solely so that we can apply our own variance normalization as part of the spatial ICA, +## rather than using the VN algorithm internal to melodic +## (The resulting ${concatfmrihp}_vnts time series is NOT the one that gets cleaned by FIX). +## See comments above about regarding the use of call_matlab here. + +AlreadyHP="-1" #Don't run highpass on concatenated data (since hp was done on the individual runs) + +case ${FSL_FIX_MATLAB_MODE} in + 0) + # Use compiled Matlab + CSUBDIR=${HCPPIPEDIR}/ICAFIX/scripts/Compiled_functionhighpassandvariancenormalize + matlab_cmd="${FSL_FIXDIR}/call_matlab.sh -c ${MATLAB_COMPILER_RUNTIME} -b ${CSUBDIR} -f functionhighpassandvariancenormalize ${tr} ${AlreadyHP} ${concatfmrihp} ${Caret7_Command}" + log_Msg "Run compiled MATLAB with command..." + log_Msg "${matlab_cmd}" + ${matlab_cmd} + ;; + + 1 | 2) + # Use interpreted MATLAB or Octave + # ${hp} needs to be passed in as a string, to handle the hp=pd* case + matlab_cmd="${ML_PATHS} functionhighpassandvariancenormalize(${tr}, '${hp}', '${fmri}', '${Caret7_Command}', '${ndhpvol}, '${ndhpcifti}', '${ndcifti});" + log_Msg "Run interpreted MATLAB/Octave with command..." + log_Msg "${FSL_FIXDIR}/call_matlab.sh -r ${matlab_cmd}" + ${FSL_FIXDIR}/call_matlab.sh -r "${matlab_cmd}" # Note: the semicolon within ${matlab_cmd} needs to be protected with quotation marks + ;; + *) + # Unsupported MATLAB run mode + log_Err_Abort "Unsupported MATLAB run mode value (FSL_FIX_MATLAB_MODE) in $FSL_FIXDIR/settings.sh: ${FSL_FIX_MATLAB_MODE}" + ;; +esac # Takuya Hayashi inserted, July 2018 - #${FSL_FIXDIR}/call_matlab.sh -l .fix.functionhighpassandvariancenormalize.log -f functionhighpassandvariancenormalize $tr $AlreadyHP $concatfmri ${FSL_FIX_WBC} - if [ -e ${concatfmri}_WF.txt ] ; then rm ${concatfmri}_WF.txt;fi # Add by TH Aug 2019 echo $ndhpvol $ndhpcifti $ndcifti > ${concatfmri}_WF.txt -log_Msg "${FSL_FIXDIR}/call_matlab.sh -l .fix.functionhighpassandvariancenormalize_riken.log -f functionhighpassandvariancenormalize_riken ${tr} ${AlreadyHP} ${concatfmri} ${FSL_FIX_WBC} ${ndhpvol} ${ndhpcifti} ${ndcifti}" - - ${FSL_FIXDIR}/call_matlab.sh -l .fix.functionhighpassandvariancenormalize_riken.log -f functionhighpassandvariancenormalize_riken ${tr} ${AlreadyHP} ${concatfmri} ${FSL_FIX_WBC} ${ndhpvol} ${ndhpcifti} ${ndcifti} - if [ "$UserDim" != "NONE" ] ; then Dim=$UserDim log_Msg "Melodic use Dim=${Dim} user specified" else - Dim=`cat ${concatfmri}_dims.txt` - log_Msg "Melodic use Dim=${Dim} estimated by IcaDim" -fi + Dim=`cat ${concatfmrihp}_dims.txt` -#Need to decide if using Wishart-based dim estimate or melodic-based estimate. Also need to decide if turning off melodic VN and feeding in ours --vn --dim=${Dim} + log_Msg "Dim: ${Dim}" +fi -# Original melodic command without logging of command -# ${FSLDIR}/bin/melodic -i $concatfmri -o ${concatfmri}.ica/filtered_func_data.ica -d -250 --nobet --report --Oall --tr=$tr -# -# melodic command with verbose output and debugging output, used when trying to determine why melodic was crashing -# melodic_cmd="${FSLDIR}/bin/melodic -i $concatfmri -o ${concatfmri}.ica/filtered_func_data.ica --nobet --report --Oall --tr=$tr --verbose --debug" +## --------------------------------------------------------------------------- +## Run melodic on concatenated file +## --------------------------------------------------------------------------- -#melodic_cmd="${FSLDIR}/bin/melodic -i $concatfmri -o ${concatfmri}.ica/filtered_func_data.ica --nobet --report --Oall --tr=$tr" +log_Debug_Msg "About to run melodic: Contents of ${concatfmrihp}.ica follow" +if [ ! -z "${log_debugOn}" ] ; then + ls -lRa ${concatfmrihp}.ica +fi -MELODIC="${FSLDIR}/bin/melodic" -MELODIC=/media/myelin/brainmappers/Connectome_Project/MultiRunICA+FIXTest/melodic -MELODIC=/mnt/pub/devel/FSL-devel/MELODIC/melodic_20170712 +#grab melodic from $PATH by default, don't hardcode it with respect to $FSLDIR +#we need to do "if which ..." because the script currently uses ERR trap +if which melodic &> /dev/null +then + MELODIC=$(which melodic 2> /dev/null) +else + #if it isn't even in $PATH, fall back on FSLDIR + MELODIC="${FSLDIR}/bin/melodic" +fi -# Added "-m `remove_ext ${ConcatName}`_brain_mask" to avoid memory issue - Takuya Hayahsi July 2018 -melodic_cmd="${MELODIC} -i ${concatfmri}_vnts -o ${concatfmri}.ica/filtered_func_data.ica --nobet --report --Oall --tr=$tr --vn --dim=${Dim} --verbose --debug -m `remove_ext ${ConcatName}`_brain_mask" +log_Msg "Running MELODIC located at: $MELODIC" +log_Debug_Msg "Beginning of melodic version log, help, and checksum" +if [ ! -z "${log_debugOn}" ] ; then + log_Debug_Msg "$MELODIC --version" + $MELODIC --version + log_Debug_Msg "$MELODIC --help" + $MELODIC --help + log_Debug_Msg "md5sum $MELODIC" + md5sum $MELODIC +fi +log_Debug_Msg "End of melodic version log, help, and checksum" + +# Use the Wishart-based dim estimate instead of the melodic-based estimate (--dim=${Dim}) +# Also turn off melodic VN and feed in the input that has already been VN'ed instead (--vn flag, which turns OFF melodic's vn) +# Added "-m ${concatfmri}_brain_mask" to avoid memory issue - Takuya Hayahsi July 2018 +melodic_cmd="${MELODIC} -i ${concatfmrihp}_vnts -o ${concatfmrihp}.ica/filtered_func_data.ica --nobet --report --Oall --tr=${tr} --vn --dim=${Dim} -m ${concatfmri}_brain_mask" +if [ ! -z "${log_debugOn}" ] ; then + melodic_cmd="${melodic_cmd} --verbose --debug" +fi log_Msg "melodic_cmd: ${melodic_cmd}" ${melodic_cmd} return_code=$? log_Msg "melodic has been run: return_code = ${return_code}" -log_Msg "melodic has been run: Contents of ${concatfmri}.ica follow" -if [ "${DEBUG}" = "TRUE" ] ; then - ls -lRa `remove_ext ${concatfmri}`.ica +log_Debug_Msg "melodic has been run: Contents of ${concatfmrihp}.ica follow" +if [ ! -z "${log_debugOn}" ] ; then + ls -lRa ${concatfmrihp}.ica fi if [ "${return_code}" -ne "0" ] ; then - log_Error "melodic has returned a non-zero code" - log_Error "Exiting this script with -1 return value." - exit -1 + log_Err_Abort "melodic has returned a non-zero code" fi # At this point (following melodic), ${concatfmrihp}_vnts is no longer necessary @@ -748,13 +714,24 @@ $FSLDIR/bin/imrm ${concatfmrihp}.ica/filtered_func_data.ica/concat_data ## Housekeeping related to files expected for FIX ## --------------------------------------------------------------------------- -cd `remove_ext ${concatfmri}`.ica +cd ${concatfmrihp}.ica -$FSLDIR/bin/imln ../${concatfmri} filtered_func_data -$FSLDIR/bin/imln filtered_func_data.ica/mask mask +# This is the concated volume time series from the 1st pass VN, with requested +# hp-filtering applied and with the mean VN map multiplied back in +$FSLDIR/bin/imln ../${concatfmrihp} filtered_func_data + +# This is the concated CIFTI time series from the 1st pass VN, with requested +# hp-filtering applied and with the mean VN map multiplied back in +# Unlike 'hcp_fix' (single-run), here we symlink to the hp-filtered CIFTI and +# use "AlreadyHP=-1" to skip any additional filtering in fix_3_clean. +if [ -f ../${concatfmri}_Atlas_hp${hp}.dtseries.nii ] ; then + $FSLDIR/bin/imln ../${concatfmri}_Atlas_hp${hp}.dtseries.nii Atlas.dtseries.nii +fi -if [ `$FSLDIR/bin/imtest ../${concat_fmri_orig}_SBRef` = 1 ] ; then - $FSLDIR/bin/imln ../${concat_fmri_orig}_SBRef mean_func +# Other necessary files (N.B. Movement regressors already prep'ed above) +$FSLDIR/bin/imln filtered_func_data.ica/mask mask +if [ `$FSLDIR/bin/imtest ../${concatfmri}_SBRef` = 1 ] ; then + $FSLDIR/bin/imln ../${concatfmri}_SBRef mean_func else $FSLDIR/bin/imln filtered_func_data.ica/mean mean_func fi @@ -767,7 +744,7 @@ mkdir -p reg cd reg i_am_at=`pwd` -log_Msg "current folder ${i_am_at}" +log_Debug_Msg "current folder ${i_am_at}" $FSLDIR/bin/imln ../../../../T1w_restore_brain highres $FSLDIR/bin/imln ../../../../wmparc wmparc @@ -776,24 +753,26 @@ $FSLDIR/bin/makerot --theta=0 > highres2example_func.mat if [ `$FSLDIR/bin/imtest ../../../../T2w` = 1 ] ; then $FSLDIR/bin/fslmaths ../../../../T1w -div ../../../../T2w veins -odt float # Added for NHP - if [ $SPECIES = Human ] ; then - $FSLDIR/bin/flirt -in ${FSL_FIXDIR}/mask_files/hcp_0.7mm_brain_mask -ref veins -out veinbrainmask -applyxfm + if [[ $TemplateMask != "" ]] ; then + $FSLDIR/bin/fslmaths $TemplateMask -bin veinbrainmask else - $FSLDIR/bin/fslmaths $TemplateMask -bin veinbrainmask + $FSLDIR/bin/flirt -in ../../../../brainmask_fs -ref veins -out veinbrainmask -applyxfm fi - #$FSLDIR/bin/flirt -in ${FSL_FIXDIR}/mask_files/hcp_0.7mm_brain_mask -ref veins -out veinbrainmask -applyxfm $FSLDIR/bin/fslmaths veinbrainmask -bin veinbrainmask $FSLDIR/bin/fslmaths veins -div `$FSLDIR/bin/fslstats veins -k veinbrainmask -P 50` -mul 2.18 -thr 10 -min 50 -div 50 veins $FSLDIR/bin/flirt -in veins -ref example_func -applyxfm -init highres2example_func.mat -out veins_exf $FSLDIR/bin/fslmaths veins_exf -mas example_func veins_exf fi + +# Return to ${ConcatFolder} +# Do not use 'cd ${ConcatFolder}', because ${ConcatFolder} may not be an absolute path cd ../.. -} -RunFix () { -AlreadyHP="-1" #Don't run highpass on concatenated data -cd ${ConcatFolder} -log_Msg "running FIX" +## --------------------------------------------------------------------------- +## Actually run FIX +## --------------------------------------------------------------------------- + +log_Msg "Running FIX" # Changes to handle user specified training data file if [ "X${TrainingData}" != X ]; then @@ -809,66 +788,109 @@ if [ "X${TrainingData}" != X ]; then if [ ! -f "${TrainingData}" ]; then TrainingData=${FSL_FIXDIR}/training_files/${TrainingData} fi - log_Msg "User has specified a training data file: ${TrainingData}" # finally, if the TrainingData file is not found, report an error and get out of here if [ ! -f "${TrainingData}" ]; then - log_Error "FIX training data not found: ${TrainingData}" - exit -1 + log_Err_Abort "FIX training data not found: ${TrainingData}" fi # added for fix_0c_reg_masks - Takuya Hayashi cp ${FSLDIR}/etc/flirtsch/ident.mat ${concatfmri}.ica/reg/standard2example_func.mat else - log_Msg "User has NOT specified a training data file" - log_Msg "using training data ${FSL_FIXDIR}/training_files/HCP_hp2000.RData" - training_hp="2000" + # User has not specified a training data file + #TSC: so, let's look for it and give a useful error, rather than pretending they didn't want to do what they said + #TSC: if you WANT to use, say, hp150 with hp2000 training, it should be explicitly requested by the user + #TSC: exception: we will recommend polynomial detrending for performance, and hp 2000 is already trained, and largely indistinguishable + log_Msg "training data file not specified" + training_hp="${hp}" + if [[ ! -f "${FSL_FIXDIR}/training_files/HCP_hp${hp}.RData" ]]; then + if [[ "${pdFlag}" == "TRUE" && "${hp}" != "pd0" ]]; then # "hp=pd0" would be a demeaning only, so doesn't qualify here + #hack: hp 2000 is close enough to a polynomial detrend, but the latter is far faster + #we know the hp 2000 training data exists, so use it if we don't see an explicit hppd# set + training_hp=2000 + else + log_Err_Abort "no standard training data found for specified high pass (${hp}), please specify training data manually or use a standard high pass setting" + fi + fi TrainingData=${FSL_FIXDIR}/training_files/HCP_hp${training_hp}.RData fi +log_Msg "using training data file: ${TrainingData}" - +# set up fix command +# use array for whitespace safety, even if the rest of the script isn't if [[ ${doMotionRegression} == "TRUE" ]]; then - if [ $SURFGEOM = 0 ] ; then - fix_cmd=("${FSL_FIXDIR}/fix" "${concatfmri}.ica" "${TrainingData}" "${FixThresh}" -m -h "${AlreadyHP}") - else - fix_cmd=("${FSL_FIXDIR}/fix_hcp" "${concatfmri}.ica" "${TrainingData}" "${FixThresh}" -m -h "${AlreadyHP}") - fi - + fix_cmd=("${FSL_FIXDIR}/fix" "${concatfmrihp}.ica" "${TrainingData}" "${FixThresh}" -m -h "${AlreadyHP}") else #-h is actually a subargument to -m, and will cause problems if specified without (or even not directly following) -m # Update (11/8/2019): No longer true as of FIX 1.06.12, but since filtering of the CIFTI in MR-FIX occurs in # functionhighpassandvariancenormalize it doesn't matter, so we leave the code unchanged so that we don't # need to add a FIX version dependency to the script - fix_cmd=("${FSL_FIXDIR}/fix" "${concatfmri}.ica" "${TrainingData}" "${FixThresh}") + fix_cmd=("${FSL_FIXDIR}/fix" "${concatfmrihp}.ica" "${TrainingData}" "${FixThresh}") fi - log_Msg "fix_cmd: ${fix_cmd[*]}" +## MPH: The 'fix' script itself will continue to log to its own custom files +## Alert user to where those are +log_Msg "Check ${concatfmrihp}.ica/fix/logMatlab.txt for log output from feature extraction" +log_Msg "Check ${concatfmrihp}.ica/.fix_2b_predict.log for log output from component classification" +log_Msg "Check ${concatfmrihp}.ica/.fix.log for log output from cleanup stage" "${fix_cmd[@]}" - return_code=$? -if [ "${return_code}" -ne "0" ] ; then +if [ "${return_code}" -ne "0" ]; then log_Err_Abort "return_code from fix_cmd: ${return_code}" fi log_Msg "Done running FIX" -} +## --------------------------------------------------------------------------- +## Rename some files (relative to the default names coded in fix_3_clean) +## --------------------------------------------------------------------------- -RunPostFix () { -cd ${ConcatFolder} +$FSLDIR/bin/immv ${concatfmrihp}.ica/filtered_func_data_clean ${concatfmrihp}_clean +if [ "$?" -ne "0" ]; then + log_Err_Abort "Something went wrong; ${concatfmrihp}.ica/filtered_func_data_clean wasn't created" +fi +if [ -f ${concatfmrihp}.ica/Atlas_clean.dtseries.nii ]; then + /bin/mv ${concatfmrihp}.ica/Atlas_clean.dtseries.nii ${concatfmri}_Atlas_hp${hp}_clean.dtseries.nii +else + log_Err_Abort "Something went wrong; ${concatfmrihp}.ica/Atlas_clean.dtseries.nii wasn't created" +fi -log_Msg "Moving ${concatfmri}.ica/filtered_func_data_clean to ${concatfmri}_clean" -$FSLDIR/bin/immv ${concatfmri}.ica/filtered_func_data_clean ${concatfmri}_clean -$FSLDIR/bin/immv ${concatfmri}.ica/filtered_func_data_clean_vn ${concatfmri}_clean_vn +# The variance normalization ("_vn") outputs of fix (fix_3_clean) require use of fix1.067 or later +# So check whether those files exist before moving/renaming them +if [ `$FSLDIR/bin/imtest ${concatfmrihp}.ica/filtered_func_data_clean_vn` = 1 ]; then + $FSLDIR/bin/immv ${concatfmrihp}.ica/filtered_func_data_clean_vn ${concatfmrihp}_clean_vn +fi +if [ -f ${concatfmrihp}.ica/Atlas_clean_vn.dscalar.nii ]; then + /bin/mv ${concatfmrihp}.ica/Atlas_clean_vn.dscalar.nii ${concatfmri}_Atlas_hp${hp}_clean_vn.dscalar.nii + if [ -f ${concatfmri}_Atlas_hp${hp}_clean.dtseries.nii ] ; then + ${FSL_FIX_WBC} -cifti-math "TCS/VN" ${concat_fmri_orig}_Atlas_hp${hp}_clean_vn.dtseries.nii -var TCS ${concat_fmri_orig}_Atlas_hp${hp}_clean.dtseries.nii -var VN ${concat_fmri_orig}_Atlas_hp${hp}_clean_vn.dscalar.nii -select 1 1 -repeat + fi +fi +log_Msg "Done renaming files" + +# Remove the 'fake-NIFTI' file created in fix_3_clean for high-pass filtering of the CIFTI (if it exists) +$FSLDIR/bin/imrm ${concatfmrihp}.ica/Atlas + +# Always delete things with too-generic names +$FSLDIR/bin/imrm ${concatfmrihp}.ica/filtered_func_data +rm -f ${concatfmrihp}.ica/Atlas.dtseries.nii + +# Optional deletion of highpass intermediates +# (hp<0 not supported in this script currently, so no need to condition on value of hp) +if [ "${DeleteIntermediates}" == "TRUE" ] +then + $FSLDIR/bin/imrm ${concatfmri} ${concatfmrihp} + rm -f ${concatfmri}_Atlas.dtseries.nii ${concatfmri}_Atlas_hp${hp}.dtseries.nii +fi -log_Msg "Checking for existence of ${concatfmri}.ica/Atlas_clean.dtseries.nii" -if [ -f ${concatfmri}.ica/Atlas_clean.dtseries.nii ] ; then - log_Msg "Moving ${concatfmri}.ica/Atlas_clean.dtseries.nii to ${concat_fmri_orig}_Atlas_hp${hp}_clean.dtseries.nii" - /bin/mv ${concatfmri}.ica/Atlas_clean.dtseries.nii ${concat_fmri_orig}_Atlas_hp${hp}_clean.dtseries.nii - /bin/mv ${concatfmri}.ica/Atlas_clean_vn.dscalar.nii ${concat_fmri_orig}_Atlas_hp${hp}_clean_vn.dscalar.nii - /bin/mv ${concatfmri}.ica/Atlas_unst.dtseries.nii ${concat_fmri_orig}_Atlas_hp${hp}_unst.dtseries.nii - ${FSL_FIX_WBC} -cifti-math "TCS/VN" ${concat_fmri_orig}_Atlas_hp${hp}_clean_vn.dtseries.nii -var TCS ${concat_fmri_orig}_Atlas_hp${hp}_clean.dtseries.nii -var VN ${concat_fmri_orig}_Atlas_hp${hp}_clean_vn.dscalar.nii -select 1 1 -repeat +# Convert sICA+FIX cleaned movement regressors to text +if [ -f ${concatfmrihp}.ica/mc/prefiltered_func_data_mcf_conf_hp_clean.nii.gz ] ; then + fslmeants -i ${concatfmrihp}.ica/mc/prefiltered_func_data_mcf_conf_hp_clean.nii.gz -o Movement_Regressors_hp${hp}_clean.txt --showall + # Strip header lines included as part of '--showall' flag + nVols=$($FSLDIR/bin/fslnvols ${concatfmrihp}.ica/mc/prefiltered_func_data_mcf_conf_hp_clean.nii.gz) + # Execute tail in a subshell, so we can successfully overwrite file with same name + echo "$(tail -n ${nVols} Movement_Regressors_hp${hp}_clean.txt)" > Movement_Regressors_hp${hp}_clean.txt fi #Cleaned volume and CIFTI have no mean and are intensity normalized by multiplying in the overall variance normalization @@ -914,18 +936,6 @@ for fmri in $fmris ; do Start=`echo "${Start} + ${NumTPS}" | bc -l` done -log_Msg "Finished" - -} - -if [ "$RunMode" = 0 ] ; then - RunPreFix;RunICA;RunFix;RunPostFix -elif [ "$RunMode" = 1 ] ; then - RunICA; RunFix;RunPostFix -elif [ "$RunMode" = 2 ] ; then - RunFix;RunPostFix -elif [ "$RunMode" = 3 ] ; then - RunPostFix -fi cd ${DIR} +log_Msg "Completed!"