dwipreproc with half sphere sampling

[neurolabusc]

I tried dwipreproc with historical data where only half-sphere sampling was used and there was no reverse phase encoding. I was surprised that the script selected eddy instead of the older eddy_correct, as this is discouraged https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/eddy I think it would be great if the dwipreproc script included logic that if the user selects -rpe_none it either generates a warning or uses eddy if sampling is on the half sphere. Off the top of my head, I think you could detect a half sphere by computing the mean of bvec - if this is zero it is a whole sphere, if not normalize the length and compute find the bvec with the largest cross product. Alternatively, the documentation could be enhanced to describe warn users to avoid this situation.

[Lestropie]

It's discouraged to use half-sphere sampling at the acquisition stage, but personally I wouldn't interpret from that documentation that eddy_correct is recommended over eddy in that instance. Indeed I suspect it wouldn't be the case, and I doubt FMRIB would recommend as such; this is particularly true for higher b-values, where eddy_correct can cause more problems than it solves (since it doesn't predict the signal for a particular diffusion direction before performing the registration, so the image contrast between target and source can vary markedly).

It might be possible to automate the inclusion of --slm=linear in the case of a half-sphere acquisition as suggested in that linked page; but that can already be provided by the user, by utilising the -eddy_options command-line option in dwipreproc.

[neurolabusc]

OK, automatically detecting this or changing the documentation to suggest -eddy_options " --slm=linear " sound like reasonable solutions.

[jdtournier]

I just contacted Jesper Andersson about this to get his opinion, here's what he has to say:

Yes, it sounds like Chris is a little overcautious. It is a good idea to sample on the whole sphere, and we encourage new protocols like that. That does not mean that eddy doesn’t work when sampled on the half sphere. But there are some gotchas to be aware of.

Let’s say we have sampled on the half sphere, and let’s say that it is along the z-axis that we only have positive values for the diffusion gradients. Let us further say that the z-component causes a shear in the z-PE-plane that is proportional to that z-component. Eddy doesn’t attempt to explicitly warp to an undistorted space, because it has no knowledge of that space. Instead it warps all image towards a common space that can be thought of as an “average” space of the input images. So in our hypothetical case all images will be warped towards a point that has a z-PE-shear corresponding to the average z-component of all the scans.

So, we are now in a common space for all images (i.e. the EC-induced inter-volume variability is gone), but that space is not exactly the true space. That is the reason I recommend --slm=linear. What that does is to fit a straight line (with z-diffusion component on the x-axis) to all the estimated z-PE-shear values, and then extrapolates those to the point where the z-diffusion component is zero. And these are the values that eddy then uses. As far as we have been able to test that does a very good job also on data sampled on the half sphere.

So I think we're OK, but we could consider explicitly mentioning this in the docs - simple enough fix...

[neurolabusc]

@Lestropie and @jdtournier - thanks for the rapid and comprehensive answer. I agree the issue can be closed once noted in the documents. Certainly anyone should be encouraged to use full sphere-sampling and reverse phase-encoding. For anyone working with archival datasets, here is my Matlab script to detect half-sphere sampling:

function isHalfSphere = HalfSphere(bvec_nam)
%return true if DWI vectors sample half sphere
% bvec_nam : filename of b-vector file
%Example
% HalfSphere('dki.bvec')
if ~exist('bvec_nam','var') %file not specified
   fprintf('Select b-vec file\n');
   [A,Apth] = uigetfile({'*.bvec'},'Select b-vec file');
   if isnumeric(A), return; end; %nothing selected
   bvec_nam = strcat(Apth,char(A));
end;
isHalfSphere = false;
bvecs = load(bvec_nam); 
bvecs = unitLengthSub(bvecs)';
bvecs(~any(~isnan(bvecs')),:) = [];
if isempty(bvecs), return; end;
mn = unitLengthSub(mean(bvecs)')';
if isnan(mn), return; end;
mn = repmat(mn,size(bvecs,1),1);
minV = min(dot(bvecs',mn'));
thetaInDegrees = acosd(minV);
if thetaInDegrees < 110
    isHalfSphere = true;
    fprintf('Sampling appears to be half sphere (%g-degrees): %s\n', thetaInDegrees, bvec_nam);
end;
%end HalfSphere()

function Xu = unitLengthSub(X) 
%Use pythagorean theorem to set all vectors to have length 1
%does not require statistical toolbox 
if size(X,1) ~= 3, error('expected each column to have three rows (X,Y,Z coordinates)'); end;
Dx =  sqrt(sum((X.^2)));
Dx = [Dx; Dx; Dx];
Xu = X./Dx;
%end unitLengthSub()

[Lestropie]

I was thinking about adding a check within dwipreproc that would just look at the sum of unit vectors within each shell, and issue an appropriate warning if that sum is not sufficiently close to zero; that way the heuristic would be sensitive to half-sphere schemes where the half-sphere is oriented in any direction. However it's unclear from Jesper's description whether or not --slm=linear makes an assumption about the half-sphere being oriented along the z-axis, or whether it will help for any arbitrarily-oriented half-sphere scheme.

[jdtournier]

I was thinking about adding a check within dwipreproc that would just look at the sum of unit vectors within each shell, and issue an appropriate warning if that sum is not sufficiently close to zero; that way the heuristic would be sensitive to half-sphere schemes where the half-sphere is oriented in any direction.

I'd favour that approach if we were to go down that route. We could also extend dirstat to support reading the gradient table from images directly (I've been meaning to do that for a while...), and issue a warning about half-sphere sampling in that case also. One slight concern is about how we settle on what threshold constitutes a half-sphere scheme - should it also report on badly distributed schemes that are not strictly half-sphere, but highly asymmetric nonetheless? And if so, how bad should it be before we issue the warning...?

I've had a quick go at looking into this, running dirgen for 6 to 330 directions, then looking at the norm of the vector mean in Matlab. The straight dirgen output is not in any way optimised for spherical sampling (so it's essentially random). The half-sphere case (flipping all z<0 vectors) shows a very consistent value of 0.5. Running dirflip on the output of dirgen helps to bring that vector mean closer to zero (as it's supposed to...) - although this was only a brief run of 10k samples to keep the runtime manageable (dirflip is essentially brute-force search - might be a touch better with a longer run).

As to the --slm-linear option, from the eddy wiki page:

--slm

"Second level model" that specifies the mathematical form for how the diffusion gradients cause eddy currents. For high quality data with 60 directions, or more, sampled on the whole sphere we have not found any advantage of performing second level modelling. Hence our recommendation for such data is to use none, and that is also the default.

If the data has quite few directions and/or is has not been sampled on the whole sphere it can be advantageous to specify --slm=linear.

I'm guessing this means: assume eddy currents are linear function of gradients (?) - in which case it ought to be independent of how the directions are sampled. But what I don't get is if (as stated) there is no disadvantage to running with --slm=linear, why not run with that by default...?

[Lestropie]

Looks like that approach won't give as clean a separation as I'd expected. Maybe instead, something like: Look for a direction that, when inverted antipodally, the nearest direction (not taking into account antipodal symmetry) is more than twice the nearest direction without the antipodal flip?

But what I don't get is if (as stated) there is no disadvantage to running with --slm=linear, why not run with that by default...?

Suspect it's similar to the post-eddy alignment of shells (PEAS): There's some level of uncertainty in the estimation, so if that uncertainty exceeds the magnitude of the actual error, it's recommended to skip the estimation altogether. I tend to just trust that Jesper has a better understanding / more exhaustive testing of this stuff than I do.

[jdtournier]

Looks like that approach won't give as clean a separation as I'd expected.

?!? Looks pretty clean to me... Anything close to 0.5 is clearly half-spherical. The grey area is the 0.1 to 0.2 range - it's clearly not half-spherical, but not optimised for full spherical coverage either. This reflects the reality here...

If you were expecting the properly full-spherical coverage to be much closer to zero, then yes, it's not clean - but that's the reality with spherical sampling: uniformity is only ever approximate...

Maybe instead, something like: Look for a direction that, when inverted antipodally, the nearest direction (not taking into account antipodal symmetry) is more than twice the nearest direction without the antipodal flip?

Not sure I get this... :confused: Is this to check for half-spherical sampling specifically? If so, I think testing for norm of mean vector > 0.4 is more than sufficient...

[neurolabusc]

Here s an example of half-sphere sampling data - this may be a good example, as it is a popular DKI sequence. This particular sequence uses a twice-refocused spin-echo sequence that reduces distortion at the expense of reduced SNR (due to longer TE). These schemes were popular prior to Topup/Eddy (which allows us to have the SNR of a short TE and good undistortion).

dti.bval.txt dti.bvec.txt

[jdtournier]

Yep, that comes out bang on 0.5 for each shell:

>> norm (mean (dti_bvec(:,dti_bval==1000)'))

ans =

    0.5020

>> norm (mean (dti_bvec(:,dti_bval==2000)'))

ans =

    0.5013

Not sure whether that helps us figure out what to do though. I guess a good strategy might be to allow dirstat to take an image as input and give it further options to output specific stats so it can be used in scripts more easily. Then make dwipreproc use this check whether the input scheme is OK, and if no output a warning with a hopefully helpful hint to run with -eddy="--slm=linear " or something...?

[Lestropie]

Anything close to 0.5 is clearly half-spherical.

Well, not quite: There's a peak above 0.4 at ~ 20 directions. But then again, maybe that warning should be thrown for such a scheme? Depends on whether the warning should be issued only for something that's most definitely been designed as a hemispherical scheme vs. something that just kinda came out that way. Not that either is good from an acquisition perspective.

Not sure I get this...

Was probably not the best explanation of all time. But I have somewhere in memory a scheme that was almost hemispherical, but not quite; probably came from polyhedral tessellation, in that instead of the "border" between directions and no-directions being a great circle, it was kind of a jagged line.

Basic idea was: Imagine you have a hemispherical scheme where all directions have a non-negative Z. If you take the direction [0,0,1], the nearest direction to that is going to be only some small angle away. Whereas if you flip that to get [0,0,-1], and look for the nearest direction without anti-podal symmetry, the nearest direction will be ~ 90 degrees away.

I guess a good strategy might be to allow dirstat to take an image as input and give it further options to output specific stats so it can be used in scripts more easily.

A more modular approach rather than having something specifically in dwipreproc would be preferable. Though having dir* commands accept images might in itself be too high-level; really all we need is to be able to export the gradient directions of a specific shell to a file as [az,el] pairs.

New command dirconvert (to go from spherical to cartesian & vice-versa, and extract directions of a specific shell without the need to pipe all of the image data through dwiextract)?

[jdtournier]

Anything close to 0.5 is clearly half-spherical.

Well, not quite: There's a peak above 0.4 at ~ 20 directions. But then again, maybe that warning should be thrown for such a scheme?

Yep, it should... The blue line (dirgen only) is for essentially randomly oriented directions (i.e. no preference between +x & -x) - not surprising some of them come out very suboptimal, especially for low numbers of directions. What you should look at is the yellow line (dirgen + dirflip), which uses the exact same dirgen output and finds a combination of inversions for near-optimal full-spherical coverage. You can see that there is a solution with the exact same directions that also provides good full-sphere coverage (less than 0.1 for the peak you're referring to).

Depends on whether the warning should be issued only for something that's most definitely been designed as a hemispherical scheme vs. something that just kinda came out that way. Not that either is good from an acquisition perspective.

Yes, I don't think there's a lot of value in just detecting schemes that are strictly half-spherical. If they're biased towards one side, it's still not good, and we should warn about that.

A more modular approach rather than having something specifically in dwipreproc would be preferable. Though having dir* commands accept images might in itself be too high-level; really all we need is to be able to export the gradient directions of a specific shell to a file as [az,el] pairs.

We don't need that - dirstat will also accept Cartesian files, and full 4-column grad files (in which case you get a per-shell breakdown). All we need is to provide options to give specific metrics of interest, on the shell(s) of interest. I've already started on this on the dirstat_enhancements branch. But I do think it's worth adding support to load the directions direct from the images - it's a minor change, and will make it much easier to use as a quick check for users when they inspect their data. And it's already done anyway in that branch... 😉

[jdtournier]

@neurolabusc, @Lestropie: does #1157 address all the issues raised here? If so, let's close.

[neurolabusc]

Perfect. I am happy to see this issue closed. This solution helps others, and specifically helps me a lot (as I used to use eddy_correct for half sphere). Thanks!