Enable Creation of Cleaned Movement Regressors

[glasserm]

Regress sICA noise components out of movement regressors to enable their use on sICA+FIX cleaned fMRI data. The files are present in text and .nii.gz formats and for both the concatinated and single runs

[glasserm]

fix_3_clean.m code:

%%%% compute movement regressors with noise components removed
if domot == 0
  confounds = functionmotionconfounds(TR,hp);
  if aggressive == 1
    % aggressively regress out noise ICA components from movement regressors
    betaconfounds = pinv(ICA(:,DDremove),1e-6) * confounds;                              % beta for confounds (bad only)
    confounds = confounds - (ICA(:,DDremove) * betaconfounds);    % cleanup
  else
    % non-aggressively regress out noise ICA components from movement regressors
    betaconfounds = pinv(ICA,1e-6) * confounds;                              % beta for confounds (good *and* bad)
    confounds = confounds - (ICA(:,DDremove) * betaconfounds(DDremove,:));    % cleanup
  end
  save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),'mc/prefiltered_func_data_mcf_conf_hp_clean','f',[1 1 1 TR]);
end

[coalsont]

FYI, tail can take a filename as argument, and can be specified as "skip n lines" using +.

[glasserm]

I am comfortable with the existing usage.

[mharms]

Given that ${ConcatFolder}/Movement_Regressors_hp${hp}_clean.txt has already been created, a seemingly much simpler way to accomplish this would be:

outFile=$(dirname $fmri)/Movement_Regressors_hp${hp}_clean.txt
sed -n "${Start},${Stop}p" ${ConcatFolder}/Movement_Regressors_hp${hp}_clean.txt > ${outFile}

[glasserm]

I wanted the intermediate NIFTI file in place with the current implementation.

[mharms]

Why do we want a NIFTI version of the movement regressors lying around? If it serves some purpose, can you add a comment explaining what its value is?

[mharms]

1. Should we also be deriving a _dt version of the cleaned motion regressors?
2. Once we have hcp_fix_multi_run ready, we need to transfer the relevant code to hcp_fix as well (for the single-run condition).
3. There is an oddity in that the "cleaned" version of the Movement_Regressors will be nTP x 24, whereas the original is nTP x 12. Not sure if we want to do anything about that, but raising the awareness of it...

[glasserm]

Overall many of the above changes fall in the category of probably accomplishing the same thing with different code, which I am not really sure is necessary (as you may see, I generally don't make these kind of comments on pull requests as it seems like it just generates extra work for the person implementing the request requiring more code changes and testing for the same result).

1. These data are already detrended as stated in the name.
2. hcp_fix does not really need this code does it? This was to allow combined processing of resting state and task data with MR+FIX but differential treatment of movement regressors. This situation would not arise for SR+FIX. Honestly, at this point SR+FIX is kind of deprecated in favor of MR+FIX, and I'm not sure it is worth spending a lot of time on it if we don't plan to recommend people use it.
3. That is true, but I expect we need to do it this way as it is probably not identical to compute derivatives after regression instead of before.

[coalsont]

Confusing code is bad code in that it makes it more challenging to work on the code later, especially if it is someone else that has to work on it. The suggested changes should be fairly easy to make and test. Part of the point of code review is to find simpler and less confusing code to do the same thing.

[mharms]

Code readability and easy interpretability should be valued as well.

1. I don't see _dt in the file name.
2. Adding the analog to hcp_fix would allow one to do the "protected" motion regression within the context of SR-FIX. It should be trivial to add, so why not make that option available (e.g., for people processing long single-run data).
3. Computing the backward differences might be fine, but the squared terms would indeed not be the same.

[glasserm]

Well I responded to each comment above. I don't agree with any of them in this case. Both of you at times don't see that the costs of changing and testing stuff just for the sake of it are real. I likely also don't always see the costs of code that is quickly written being confusing to others. We do need to balance these things, but in this case the decisions I made above were well considered.

[glasserm]

1. _hp${Highpass}_clean = detrend
2. Because everything takes time and this isn't a priority for me.
3. The user would be able to compute derivatives after the fact or use the ones with regressed out timeseries, so this should work for both users.

[mharms]

@glasserm In this case, both @coalsont and myself are indicating that we don't have the right balance. The usage statement for fslmeants clearly indicates that the mask is an optional argument. So that particular code is completely unnecessary (and outright to confusing to anyone that might be trying to understand why it is there).

[glasserm]

Here is the usage of fslmeants:

fslmeants -i filtered_func_data -o meants.txt -m my_mask
fslmeants -i filtered_func_data -m my_mask
fslmeants -i filtered_func_data -c 24 19 10

As you can see,
fslmeants -i filtered_func_data --showall
is not one of the given choices. I interpreted "optional" in this case to mean that either -m or -c was required to be specified. I could be wrong, but elected to implement what I knew would work, rather than what might not work.

[coalsont]

I tested --showall without -m for you, since you were unwilling. Unless the output format is different, it seems to have worked. None of the help examples use --showall, but you were willing to try that option and see what happened...

I don't think changing the cat ... | tail is important, but the other lines definitely are more confusing than what they are doing, and I think they either deserve edits to make them simpler, or edits to add comments explaining why the code isn't simpler.

[mharms]

Unrelated, but Sri just coincidentally brought this to my attention, as he was attempting to ascertain some timing info from the log files.

Can we move this bit of code in hcp_fix_multi_run

#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

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"

down to just after here:

## ---------------------------------------------------------------------------
## Run melodic on concatenated file
## ---------------------------------------------------------------------------

so that the message about "Running MELODIC" (and version info if log_debugOn exists) actually occurs right before we run melodic.

thx

[glasserm]

I tested it without the mask and it worked, so I removed those lines. As for this other change, feel free to make it yourself and it does not require a pull request. I am going to merge this unless there is any last comment this evening.