spinalcordtoolbox icon indicating copy to clipboard operation
spinalcordtoolbox copied to clipboard

Averaging of motion parameters should use unsigned displacement

Open jcohenadad opened this issue 1 year ago • 16 comments

Feature requested by a user: https://forum.spinalcordmri.org/t/computation-of-motion-parameters-sct-fmri-moco/1200/1

In sct_fmri_moco and sct_dmri_moco, there is a feature to average motion parameters (signed translations) across time points, for each slice. However, this averaging is done on the signed translations:

https://github.com/spinalcordtoolbox/spinalcordtoolbox/blob/29476b99e05d59a2b42ed7d20d4f40a58a53c497/spinalcordtoolbox/moco.py#L425-L426

~The problem is that, if for time point i Tx=-0.5mm and for time point i+1 Tx=0.5mm, the average would be 0mm.~

Proposed solution:

Compute the root mean square displacement, or, as suggested by the user, perform the average on the absolute (ie: unsigned) displacement values.

2024-01-22 EDIT: Sorry, I misunderstood the actual issue here, and my explanation above is wrong (I struck through the wrong part). What the output moco_params.tsv file is, is an average across SLICES, for each time point (and not an average across time point, for each slice).

jcohenadad avatar Jan 17 '24 13:01 jcohenadad

Mathematically, and thinking about algorithms that might use these generated values, I think signed values for X and Y (which is the current behaviour) makes the most sense.

It sounds like you want to answer the question "for each slice, how much motion was there on average?", in which case separating X and Y is probably not useful. I would suggest adding a third column (or a new file) for the slice-wise average of the magnitude, sqrt(x*x + y*y), which is a meaningful and always positive number (assuming that X and Y are both in mm, otherwise we should convert to mm first).

mguaypaq avatar Jan 17 '24 17:01 mguaypaq

Should we get rid of the X and Y column then, and only keep 'sqrt(xx + yy)'?

jcohenadad avatar Jan 17 '24 17:01 jcohenadad

These parameters are mainly used to remove motion-related fluctuation in the signal. The intent is to identify the amount of motion in a specific volume, and thus the parameters should be in my opinion sensitive to the magnitude of the motion. The rms, absolute mean or the maximum of the slicewise motions may be better options than the current mean calculation. However, I agree that using both x and y (as suggested by @mguaypaq) make probably more sense. I would suggest to keep at least one metric for each X and Y direction (not sure yet which one) and add the magnitude ('sqrt(xx + yy)') as a third metric in the outputs.

CarolineLndl avatar Jan 17 '24 23:01 CarolineLndl

I agree with Caroline to keep 2 separate X and Y metrics, and adding a third composite one as suggested. Better to have more variables available than too few.

combesaje avatar Jan 18 '24 14:01 combesaje

@CarolineLndl @combesaje are the current mean(X) and mean(Y) metrics of any use at the moment? If so for what

@mguaypaq I'm not sure I see your point about the value of the averaged signed X and Y. I really don't see how it would be useful.

My concern in keeping them, is that they might be misused (as might have been the case already by previous researchers).

jcohenadad avatar Jan 19 '24 02:01 jcohenadad

Actually what I've used are the slicewise moco_X and moco_Y files for denoising fMRI (also potentially QC and scrubbing), not the slice-average in the .tsv file so maybe you're right that they have little use. Generally I like having access to more granular data but since those params can be calculated by the user from the moco niftis I don't feel strongly about it.

combesaje avatar Jan 19 '24 11:01 combesaje

Like @CarolineLndl, I've been using the volume-wise averages of x and y motion correction as confounds in the fMRI analysis, to account for motion detected in volumes. I now see that RMS would be a better metric for this.

Yet X and Y data might still be of use to better understand movement dynamics having taken place while scanning, like the participant shifting position / sliding, which could be seen through framewise displacement computation for instance ?

RaphaSchl avatar Jan 19 '24 18:01 RaphaSchl

Yet X and Y data might still be of use to better understand movement dynamics having taken place while scanning, like the participant shifting position / sliding, which could be seen through framewise displacement computation for instance ?

Just to clarify: we do not propose to get rid of the volume-wise and slice-wise X and Y translation parameters. Instead, we propose to get rid of the average X and Y parameters, and replace them with ('sqrt(xx + yy)') (calculated for each slice).

jcohenadad avatar Jan 19 '24 19:01 jcohenadad

Ok I have done some tested. I think mean slice-wise magnitude is the best option. As you can see after the timeserie cleaning (3rd graph) it is really similar to using x and y rms or absolute mean. However, this will reduce the model's DoF (only one regressor compared to two) and probably explain the movement better (you can see it on figure 1 and 2 black line). I my opinion we should remove the actual mean from the output and only keep the magnitude. The correlation between several timeseries (which is the main metric used) should not be modified that much using either magnitude, rms or abs mean and a little for the mean (in case of large movements). Signal amplitude should have been impacted a little more. Motion_params_21012024.pdf

CarolineLndl avatar Jan 21 '24 16:01 CarolineLndl

Great! It sounds like we're converging on a consensus for replacing the two slice-wise averages (in X and Y) with slice-wise mean magnitude (that is, sqrt(x*x+y*y)). I'll implement that.

mguaypaq avatar Jan 22 '24 21:01 mguaypaq

OK, first of all, I apologize but I misunderstood what moco_params.tsv was showing in the first place. I edited my comment https://github.com/spinalcordtoolbox/spinalcordtoolbox/issues/4344#issue-2086237735.

Now, I have a bit of trouble understanding @CarolineLndl experiment here https://github.com/spinalcordtoolbox/spinalcordtoolbox/issues/4344#issuecomment-1902691319 and its conclusion.

My questions/comments:

  • On plots 1 and 2, what does 'motion_magn' mean, and how is it computed?
  • And to be sure: is 'motion_x_mean' what is currently output by moco_params.tsv? I.e. average of signed Tx across all slices, for each time point?
  • And what about motion_x_rms: is it sqrt(Tx**Tx)?
  • Now the plot 3 about the signal time series: is this showing the time series of a single voxel at a specific location that is supposed to correlate with subject motion? (eg: at an interface between two structures)? If so, then what are rms, mean and mean_abs?
  • What we currently propose to do is sqrt(x*x+y*y), which I interpret as:

$$\sqrt {\sum_{slice=1}^{n} T_x[slice]^{2} + T_y[slice]^{2} } $$

What I find problematic with this approach, is that it scales with the number of slices, making it difficult to come up with some sort of 'standard' threshold for 'what is a lot of motion, and what is not'. Shouldn't we be normalizing it with the number of slices?

jcohenadad avatar Jan 22 '24 22:01 jcohenadad

Sorry for the lack of explanation:

  • motion_magn is the mean magnitude across slices (average of magnitude across all slices for each time point)
  • motion_x_mean is the current moco_params.tsv (average of signed Tx across all slices for each time point)
  • motion_x_rms is the rms across slices (np.sqrt(np.mean(array_X**2)) across all slices for each time point)
  • Plot 3 is a signal from a random voxel of interest (in the cord) for which I regressed out the motion parameter (either rms, mean of mean_abs). The volume wise motion parameters are one of the regressors used during the denoising step, to remove the contribution of the motion from the signal.

You can, in addition to the previous step, decide to regress out specific outlier volumes (the one with large motion) but in that case people calculate the DVARS (RMS intensity difference of volume N to volume N+1). From the dvars you can calculate your "threshold" for large motion. But this step is already implemented in fsl I'm not sure we need to add the DVARS calculation in the sct (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FSLMotionOutliers). In sum, what we need here are volume-based parameters that can be regressed out from the signal, and thus remove the contribution of the motion in the signal. I hope my explanations are more clear ...

CarolineLndl avatar Jan 22 '24 23:01 CarolineLndl

Thank you for the clarification @CarolineLndl, additional questions:

motion_magn is the mean magnitude across slices (average of magnitude across all slices for each time point)

By 'magnitude' do you mean the unsigned Tx value? So, for each time point, something like this:

$$\sum_{slice=1}^{n} \left| T_x[slice] \right| $$

In sum, what we need here are volume-based parameters that can be regressed out from the signal, and thus remove the contribution of the motion in the signal. I hope my explanations are more clear ...

yup, cristal clear! 💎

jcohenadad avatar Jan 23 '24 02:01 jcohenadad

For the magnitude, I think the previous formula was ok sqrt(xx+yy) $$\sqrt { \sum_{slice=1}^n T_x[slice]^{2} + T_y[slice]^{2} }$$

CarolineLndl avatar Jan 29 '24 20:01 CarolineLndl

Hi everyone,

I agree with Julien that the value of each volume should represent the average over the slices rather than the sum only.

I would personally go for the average of the unsigned Tx/Ty values, as in: $$\dfrac{\sum\limits_{slice=1}^{n} |Tx[slice]|}{n}$$ Indeed, I would assume that the motion in x and y may differently impact the signals.

Otherwise, if the consensus is more in favor of a single summary value, the previous formula (sqrt(x*x+y*y)) would be fine I suppose, provided that it is normalized by the number of slices.

Thanks a lot!

NawalK avatar Jan 30 '24 12:01 NawalK

A slight correction: the formula I had in mind has the square root in a different place: $$\frac{1}{n} \sum_{slice=1}^n \sqrt{T_x[slice]^2 + T_y[slice]^2}$$

In words: for each time point and each slice, we have a translation vector with X and Y coordinates. We take the Euclidean length of this translation vector to get a number (which is physically meaningul: if X and Y were in millimeters, then this length is also in millimeters). We then average this number across all slices to get a single number, which represents the average amount of motion at a particular time point, regardless of the direction of the motion. This is also correctly normalized so that we can ignore the number of slices.

I propose this because it's the natural way to generalize from this formula for one less dimension (meaning, only X, no Y coordinates): $$\frac{1}{n} \sum_{slice=1}^n |T_x[slice]|$$ In two dimensions, the absolute value becomes Euclidean distance.

mguaypaq avatar Mar 07 '24 19:03 mguaypaq

As a small heads-up: SCT v6.3 has been released, with the mean magnitude moco_params.tsv behavior now included. :relaxed:

joshuacwnewton avatar Apr 26 '24 12:04 joshuacwnewton