Merge multi-shell CSD into core

[jdtournier]

The actual multi-shell dwi2fod version now works, as far as I can tell, and is probably ready to be integrated into the master branch. What is still needed are the tools supporting it, in particular the estimation of the various responses. In the first instance, it would be good to have at least something that can work on any given dataset. What would be needed to get all this ready for release?

[bjeurissen]

Estimation of the various responses is indeed something that needs to be fixed.

Currently my pipeline is as follows:

1) The act_anat_prepare_fsl script is used to identify voxels belonging to the different tissue classes (based on anatomical scan).

2) The dwi2tensor and tensor2metric commands are used to obtain FA estimates to exclude CSF and GM voxels that are too anisotropic and to exclude WM voxels that are too isotropic. This results in a response voxel mask for each tissue type.

3) The old estimate_response command from the 0.2 branch is used to estimate the responses per tissue and per shell and with the right lmax. Each shell is selected using the dwiextract command with the -shell option. One response text file is saved per tissue type, each containing one line per shell.

4) The new multishell command takes a list of lmax'es (one per tissue type), and response text file/tissue odf output pairs per tissue type.

1) and 2) are probably fine, especially for now.

4) is in decent shape as well (it just lacks support for the -shell option, e.g. to exclude certain shells).

3) is currently the messy part because:

• it relies on the 0.2 branch
• it reorients the DWI signals per shell instead of for all the shells at the same time

The easiest fix (for the time being) would be for dwi2response to support:

• input of the final response mask (no iterations)
• isotropic (lmax=0) reponses
• responses per shell (meaning also accepting multiple lmaxes)

[jdtournier]

OK, good to know. For 3), you can probably use a combination of amp2sh and sh2response - the latter needs the SH fit to the DW signal, the mask within which to estimate the response, and the directions of the fibres within each voxel. I'd coded this up to provide a bit more flexibility for something I'm doing here, but with a view to provide a more flexible approach to the whole response estimation issue. Basically, using this allows you to use a Rician-bias corrected SH fit if you can get it (I'd modified amp2sh to do that a while back, but can't find the code any more...), and you can chose to use a tensor fit, a FOD peak fit, or the FOD mean lobe direction as your input direction, etc.

But this doesn't solve the issue of getting the response over all shells directly. I guess you can script that as-is for now, and we'll see whether there's a simple way of supporting multi-shell explicitly later.

In the meantime, I'll go modify the lmax detection code to allow lmax=0...

[bjeurissen]

sh2response is a great solution! Didn't know about that command...

Will prepare a script that handles the different shells.

[dchristiaens]

I had extended the 0.2 estimate_response to multi-shell data long time ago, in which a single tensor was fitted to the full voxel data, yielding a single rotation for all shells. I remember it was fairly straightforward and I don't mind sharing the code, only I don't have the time to port this to dwi2response currently.

The main difference is in the way the SF directions are estimated. dwi2response takes the iterative approach from Tax et al. (2014). Would you want to extend this to MSMT-CSD at each iteration, or rather SSST-CSD of the highest b-value shell only?

[jdtournier]

@dchristiaens - I think I'd advise against investing too much time into this at this stage. Using the sh2response approach, we can feed in the directions computed using dwi2tensor | tensor2metric -vec -mod none as-is, which would accomplish exactly what you have...

[dchristiaens]

Just came to the same conclusion checking the sh2response docs. Works for me if Ben prepares a script.

[jdtournier]

OK, managed to pull the Rician bias corrected version of amp2sh out of git - even though I'd deleted the branch. You gotta love git - tip of the day: git fsck --lost-found to find your long-lost discarded branches... :)

[bjeurissen]

Okay, the current pipeline functions as in the example below:

• dwi.mif contains 4 shells, mask*.mif contains the response voxel masks for each tissue type
• estimate the multi-shell responses for tissues 1 to 3 (supply lmax for each shell/tissue combination): msdwi2response dwi.mif mask1.mif 0,0,0,0 response1.txt msdwi2response dwi.mif mask2.mif 0,0,0,0 response2.txt msdwi2response dwi.mif mask3.mif 0,4,6,8 response3.txt
• get the ODF for each tissue (supply lmax for each tissue ODF): msdwi2fod -mask mask.mif -lmax 0,0,8 dwi.mif response1.txt odf1.mif response2.txt odf2.mif response3.txt odf3.mif

The mask*.mif files can be obtained by your method of choice. The main thing I am not happy about is the name of the command msdwi2fod. msdwi might suggest that it only works for multi-shell, while it works just as well for single shell. fod, suggests that it provides fiber orientation distribution functions, but in fact it provides just general odfs of whatever you put in the responses.

[jdtournier]

The amount of time we spend naming stuff...

Mostly agree with you on this one. Ideally, this should be included in dwi2fod... Not too worried about the fod part, it is really designed to get the FOD out, even if it is now capable of getting other tissue types out at the same time.

So maybe we should just do the right thing from the outset and merge it into dwi2fod directly. The algorithms are completely different between the regular CSD and the multi-tissue version, so we should just be able to lump them together without conflict. Looking at the msdwi2fod syntax, I think it already is entirely compatible with the dwi2fod syntax... We only need to check how many responses have been supplied.

That said, we should also provide a means to override that, for those rare cases where we'll want to use the same algorithm to compare between single-tissue and multi-tissue on an equal footing.

How does that sound? Probably a bit more work upfront, but ultimately the right place, I feel... Comments...?

[Lestropie]

Personally I'd like to add a number of things to dwi2response, which will likely ultimately result in the removal of sh2response. It should also not be too difficult to provide multi-shell response estimation in a script using Ben's existing approach for now, until an alternative mechanism (i.e. not dependent on the ACT image) is found.

So what needs to be added to dwi2response is:

• No iterative optimisation; receive a mask of single-fibre voxels from the user.
• When using the above, have the option of either doing an iterative response function calculation based on a heuristic initialisation, so that CSD can be run to get the fibre orientations, or the user provides a peak direction image (like is currently done in sh2response).
• Rician bias correction (not necessarily required in this context, but an advantage of having everything done together).

If this can be done, the multi-shell dwi2response script would proceed as follows:

1. Select non-partial-volumed voxels for each tissue type using the ACT image.
2. Feed the WM mask from step 1 to dwi2response using the highest b-value shell - use the iterative optimisation.
3. Use the resulting WM SF mask from dwi2response, re-running dwi2response once for each shell to get the WM response function for each shell.
4. For GM and CSF, remove voxels from relevant mask where FA is too high.
5. For GM and CSF, simply calculate the mean DWI intensity for each shell, and run mrstats using the masks from step 4 to get the mean intensity in each shell for each of GM and CSF. That gives you the lmax=0 'response function' for each tissue / shell combination.
6. Concatenate the three results to produce the multi-shell response function text file.

[jdtournier]

If the estimation of the response function(s) was a problem with an obvious solution, I'd tend to agree that a single dwi2response executable was the right thing to do. But if anything, what the past few months have pretty clearly demonstrated is that there is no single best way of handling this. Sure, we can keep adding checks and balances every time someone produces data that causes issues, but then the onus is on us to fix it - there is no scope for individuals not steeped in the MRtrix way to debug and/or amend the procedure. And it also causes issues when trying to do something that is clearly related, but not what it was intended for (the multi shell responses, for instance). It also means that the estimation might not so easily benefit from changes elsewhere that might help get better results - for example, it would be trivial to take advantage of the new Rician bias correction in amp2sh if that was all scripted . Basically anyone with some scripting knowledge could add in support for that feature - when it's an executable, one of us will need to do a lot more work to add in the feature and test it...

Given that the current procedure is a linear sequence of stages, I don't see there being any advantage to having one executable perform all of these otherwise unrelated steps. There is no reason that I can see why we couldn't break up the various steps into dedicated executables, and string them up in a script. That would provide a means for users to actually get their hands dirty and try to fix it without a massive learning curve...

For the current version of dwi2response, I think the only step that isn't currently already available as part of some other app is the FMLS segmentation to get the directions. Personally I think sh2peaks would be perfectly adequate, but if you feel very strongly about it, there's no reason not to code up a similar shsegment app (or whatever we might want to call it) that could be substituted in its place.

Generally speaking, the issues we've been having recently have really made me think that there is value in the old tenet of Unix: each app should do one job, and do it well...

[Lestropie]

I don't see that the -rician option should be mutually exclusive between amp2sh and dwi2response. If someone's intending to use RCSD, it would make a lot more sense to add a -rician option to dwi2response than to have a script call dwi2response / amp2sh / dwi2tensor / tensor2metric / sh2response.

And I definitely don't agree with removing dwi2response and replacing it with a script. Its presence doesn't prevent people from experimenting with other solutions, and most users just want the damn response function. It's much easier to just pull out the relevant information from the FOD segmentation within an executable, than to try to output all of that data to images then try to string it all together in a script (e.g. there's a heuristic in there for rejecting FOD lobes based on the relative size of the negative lobes; that'd be annoying to try to do in a script).

MSMT CSD is tougher and definitely needs a script (unless an automated solution is found later where you simply feed in the DWI and an optional mask and it spits out response functions based on data-driven tissue segmentation and/or more complex heuristics). But I'm not convinced that sh2response is necessary for this script (requiring explicit conversion from amplitudes to SH and calculation of directions beforehand) if you could instead handle that within dwi2response by giving it a command-line option instructing it to just use the provided mask and not do any SF voxel selection optimisation. It's still producing a response function based on DWI data, it's only the voxel selection that is different.

If you still want to go that way, getting directions for sh2response without using the tensor model would be best handled using fod2fixel to produce a fixel image, then fixel2voxel to extract the XYZ direction of each fixel (not yet implemented but not hard), then mrconvert -coord 3 0:2. But you'd still have had to derive a response function in order to do CSD in order to get those directions.

[jdtournier]

I don't see that the -rician option should be mutually exclusive between amp2sh and dwi2response.

It's not. The point is that there's no need to have the option if it can be provided by some other command. If we had a script, it would be trivial for any of us or our users to modify any one stage of the pipeline. With a big executable, it's down to the developers familiar with the app (i.e. you) to implement even relatively trivial changes. So in this case, if dwi2response was already a script, you wouldn't need to implement the -rician option, it would basically already be there.

And I definitely don't agree with removing dwi2response and replacing it with a script. Its presence doesn't prevent people from experimenting with other solutions, and most users just want the damn response function.

Well, before I added sh2response, there was no alternative within MRtrix, so no, people couldn't experiment very easily. And if they did want to experiment, they'd have to know that these other commands could be used for that, and they'd essentially need to re-implement the whole pipeline from scratch - they wouldn't have an existing script to start from, which a pretty big barrier to entry. And while I appreciate most users won't care much how their response function is obtained, they might want to have a go a debugging it if it happens not to work on their data, which is a lot easier when it's part of a script that they can trivially modify.

It's much easier to just pull out the relevant information from the FOD segmentation within an executable, than to try to output all of that data to images then try to string it all together in a script

I'm not sure why you'd say that. I had a quick look through the code for dwi2response, and it seems relatively straightforward to break up into separate parts. The main thing would be to provide a replacement for sh2peaks, which would provide the mean directions for the two largest FOD lobes within each voxel, sorted and weighted by the FOD integral (or whatever metric you deem most appropriate). And there'd be nothing stopping you from adding in all the heuristics you think are necessary to make sure the output is sensible, I don't see that it would be in any way dependent on any other stages in the pipeline (?).

For reference, this is more or less what I used for the number of directions for HARDI paper:

# start with sharpest possible response function (i.e. SH decomposition of a flat disc).
$ echo "1,-1,1,-1,1" > init_response.txt

# initial estimate of FOD, within brain mask:
$ dwi2fod dwi.mif -mask mask.mif init_response.txt fod.mif

# extract 2 largest peaks into a 4D image split over volumes:
$ sh2peaks fod.mif -mask mask.mif -num 2 peak-[].mif

# compute ratio of peak squared amplitudes (admittedly not pretty, but clear enough):
# (note this computes the ratio peak1/(peak2+1) to prevent divide-by-zero if peak2 is zero)
$ mrcalc peak-0.mif 2 -pow peak-1.mif 2 -pow peak-2.mif 2 -pow -add -add peak-3.mif 2 -pow peak-4.mif 2 -pow peak-5.mif 2 -pow -add -add 1 -add -div peak_ratio.mif

# get mask of 300 voxels with largest peak1/peak2 ratio:
$ mrthreshold peak_ratio.mif -top 300 single-fibre.mif

# update estimate of response function within that mask:
$ amp2sh dwi.mif -mask single-fibre.mif - | sh2response - single-fibre.mif peak-[0:2].mif response.txt

With the data I had at the time, this converged immediately - running this through a second iteration produced the same mask again. The big difference between this and what is currently in dwi2response, aside from using the FOD segmentation rather than the peak amplitude, is that the initial response function is the sharpest possible. This might seem like the wrong thing to do, but if you think about it, what this effectively does is low-pass filter the FOD - deconvolving a high-frequency kernel from a given input will attenuate the high-frequencies in the output compared to using to a low-frequency kernel.

In any event, I don't think there's any need to do anything about dwi2response at this stage - you've clearly invested a lot of effort into it, and you don't want that going to waste. However, moving forward with the multi-shell stuff, I would rather use a heavily scripted approach tying together small fast executables, which would give us and our users lots of flexibility to try out different things. Obviously, if some method comes up in the future that can only be done as a single integrated approach, we'd need to re-assess - but then it would probably need to be written from scratch anyway, irrespective of whether the current implementation is scripted or not.

[Lestropie]

It's not. The point is that there's no need to have the option if it can be provided by some other command.

I'd equally argue that that sort of thing should be a modular functor, as it will likely be used in other areas. And doing that correction with 'some other command' rather than within dwi2response would mean writing three images to file (shell-of-interest DWIs in SH, voxel mask, directions from FODs) that were already there in dwi2response in RAM ready to go. But six of one, half a dozen of the other.

Well, before I added sh2response, there was no alternative within MRtrix, so no, people couldn't experiment very easily.

That's only because when we were filling feature holes in MRtrix3 that were present in 0.2 we omitted estimate_response. Giving dwi2response an option for the user to specify the SF mask, and not do any subsequent SF voxel selection optimisation, would cover off that functionality, with less required pre-processing steps than sh2response, albeit missing out on the added flexibility (assuming you want the flexibility to e.g. define your fibre directions based on something other than the FODs). Again, six of one...

I had a quick look through the code for dwi2response, and it seems relatively straightforward to break up into separate parts.

There are definitely separate components, and I think Chantal's own implementation was a script; but I could see from the outset that there would be considerable disk overhead and scripting awkwardness in doing so: there were enough steps that could conceivably be done using existing commands but this particular use case was so specific that handling it all in RAM seemed sensible, and was then in no way restrictive in terms of heuristic design as I developed it. I also didn't know from Chantal's paper that it was going to cause so much grief.

Trouble with interpreting the results of FOD segmentation is the amount of higher-order information, which is a lot easier to get at from the actual segmented FOD class rather than having to export it to an image. The main one is the utilisation of negative lobe peak amplitudes and integrals for positive lobe rejection (i.e. a positive lobe has to be noticeably larger than the negative ringing in order to be retained); though in RF estimation this isn't as crucial. There's also a heuristic in there looking at fibre dispersion (ratio of integral to peak amplitude), which I guess you could get at by writing to a fixel image and storing the peak amplitude as the 'value' parameter.

I think writing FOD segmentation results to a fixel image using fod2fixel in this context is more making use of existing modular tools, as opposed to creating a dedicated command for running FOD segmentation and getting directions / amplitudes of just the two largest fixels purely for the sake of RF estimation; that'd be a very specific command prerequisite for a supposedly flexible script. If fixel2voxel is given additional -split options to extract individual directions / sizes / values, that should be able to serve the intended purpose.

For reference, this is more or less what I used for the number of directions for HARDI paper:

I suspect Siemens would be interested in that; they moaned about the execution time of dwi2response just like everyone else has.

However, moving forward with the multi-shell stuff, I would rather use a heavily scripted approach tying together small fast executables, which would give us and our users lots of flexibility to try out different things.

I've never been against using a script for multi-shell RF estimation at this point in time. I'm just not convinced that such a script has prerequisite commands that either don't exist yet or that need to reside separately to existing commands. Assuming there is no Rician bias correction yet in MSMT-CSD, the existing dwi2response could be run directly on the WM voxels identified by the 5TT image re-gridding to provide the WM SF voxels (then feed each shell back in independently using that SF mask to get the per-shell response functions), and I don't see any need to run sh2response / dwi2response for getting the GM / CSF responses (it's lmax=0 so just select your exemplar voxels using your heuristic of choice then take the mean value across DWI volumes & voxels).

For RCSD, I think adding a -rician option to dwi2response would make a lot more sense for most users than having to switch between a binary command and a (likely very short) script depending on which algorithm they will be using. Multi-shell responses, fair enough, it's fundamentally different and we don't have a black-box solution so a script makes sense.

[jdtournier]

@> I'd equally argue that that sort of thing should be a modular functor, as it will likely be used in other areas.

Fully agree.

doing that correction with 'some other command' rather than within dwi2response would mean writing three images to file

It sounds like you're worried about performance here. MRtrix is very heavily optimised for this kind of thing - this is why we use memory-mapping everywhere we can - if used properly, it reduces the overhead of loading the data to essentially nothing (although that's no longer strictly true given the delayed write-back we added a while back to avoid issues running on distributed filesystems, but it's still a remarkably quick operation). And with modern OS's with decent memory caching, you won't incur any disk IO related penalty (although the discussion around issue #199 is clearly relevant here). Finally, the major bottleneck in all of this is the CSD itself, everything else is trivial in comparison.

And I'd argue having these files on disk is a good thing: it allows us to inspect the output and debug - and most important allows users to debug their own issues.

Giving dwi2response an option for the user to specify the SF mask, and not do any subsequent SF voxel selection optimisation, would cover off that functionality, with less required pre-processing steps than sh2response, albeit missing out on the added flexibility (assuming you want the flexibility to e.g. define your fibre directions based on something other than the FODs).

I think one of the issues with the multi-shell stuff is we do need the flexibility to supply both the mask and the directions. I'm not against adding that into dwi2reponse if you feel like it, I just don't think it's necessary given that alternatives already exist.

I think writing FOD segmentation results to a fixel image using fod2fixel in this context is more making use of existing modular tools, as opposed to creating a dedicated command for running FOD segmentation

Absolutely, you don't want to make your life any harder than it needs to be...

the existing dwi2response could be run directly on the WM voxels identified by the 5TT image re-gridding to provide the WM SF voxels

As I mentioned above, unless we can provide both the mask and the directions, we can't really run dwi2response directly. I mean, we could, but there's no guarantee that the directions assumed for the fibre axis will match across shells, which would not be great.

For RCSD, I think adding a -rician option to dwi2response would make a lot more sense for most users than having to switch between a binary command and a (likely very short) script depending on which algorithm they will be using.

Sure, there's no problem with doing that, but it might not be that simple. The way I do Rician bias correction at the moment relies on having a noise map available (since a per-voxel noise estimate can be a little unstable, especially if included in the ML estimation). I've yet to commit the command to derive one (it's sitting on some other branch)... We could simply rely on users providing a precomputed map, or compute it within dwi2response if needed. Either way, it's starting to make dwi2response sound like a bit of a Swiss army knife...

Anyway, clearly we have different opinions on the matter, and that's fine. We just need to come up with a way forward... I'm personally strongly in favour of scripting the whole lot as much as possible, because this allows a much wider group of people to contribute to improving whatever strategy we devise - and makes for much easier debugging, crucially not necessarily requiring our input since others may also get familiar enough with the procedure to assist. I don't think it'll impact performance to any great extent. Maybe the script will not be very pretty, but it'll be a lot more intelligible to most users than the original C++...

For now, let's fix up whatever issues remain for dwi2response, and we'll start a fresh script for the multi-shell stuff (whether it uses dwi2response or not). We'll worry about what happens down the track later...

[Lestropie]

I think one of the issues with the multi-shell stuff is we do need the flexibility to supply both the mask and the directions. ... there's no guarantee that the directions assumed for the fibre axis will match across shells, which would not be great.

Yes, since you want to enforce the same fibre direction across shells, it's fairly equally clunky whether that feeds to dwi2response or sh2response. Unless you fed the whole volume to dwi2response and assumed that SF voxels and directions are calculated using the largest b-value shell with subsequent RF calculation for all shells; OK for a script, but a pretty heavy assumption for a command.

As I mentioned above, unless we can provide both the mask and the directions, we can't really run dwi2response directly.

What I was meaning there is that just because you have multi-shell, doesn't mean you have to revert to a different heuristic for selecting single-fibre voxels. Once you have a WM mask, you can still use dwi2response to determine a SF voxel mask. How to then get the response functions from that mask is in some respects a separate step.

The way I do Rician bias correction at the moment relies on having a noise map available (since a per-voxel noise estimate can be a little unstable, especially if included in the ML estimation).

I think I'm in preference of the pre-calculation method. But if the user needs a noise map for response function estimation for RCSD, presumably they will be using the same map for FOD estimation. So you wouldn't want to wrap that noise map estimation step in a script where intermediate images are discarded; and if you have to provide such an image as input to an RF estimation script, it's no more complex than providing that image to dwi2response using the -rician option.

So... here's a separate idea. Suspect you won't like it because it goes against the Unix tenet, but here goes anyway. Give up on the idea of having a single command to go straight from raw DWI to a response function in all cases; it's clearly not adequate as a black box even before MSMT. Remove the response function output from dwi2response and make the SF mask the default output; rename to dwi2sf or dwirecursivesf (so it's clear that it's specifically the Tax et al. method). Modify sh2response to operate on DWIs instead (so include shell selection, and duplicate the -rician capability of amp2sh) and call that dwi2response. It's combining what you're currently doing with two steps (amp2sh | sh2response) into one, which admittedly isn't quite as flexible, but you're unlikely to be estimating response functions from SH data that didn't come from DWI. Plus it may well enable extra cleverness, e.g. if you're doing Rician bias correction, you could also weight the response function fit by the inverse of the noise variance in each SF voxel.

[jdtournier]

OK, I think we're coming to a consensus...

Yes, essentially the main thing is to divorce the single-fibre mask estimation from the response function estimation itself. But not just the mask, also the directions within them, so we can feed that information to the response estimation stage and guarantee the same directions are used across shells.

Your suggestion of removing the response estimation proper from dwi2response is one idea, but it seems weird not to have as the output the very thing that is optimised as part of that process. I think it would make more sense for it to provide as an optional output the final single-fibre mask and directions, so that the information can be used in a separate pipeline, or inspected for debugging. Regardless, the point is you could then use dwi2response to produce the mask and directions and feed that information through to the multi-shell response estimation (or a different response function estimator for single shell if you felt like it), and that does indeed makes good sense.

As to merging amp2sh | sh2response into a single command, I don't see the rationale. Why would you do invest the time and effort into this when you have the functionality already? This is all going to end up in a script anyway, it's not going to be any messier for the end user. And you can still add in your suggestion of estimating a weighted response function using the noise estimate, you'd just need to pass that information to sh2response, which would be a trivial change...

[Lestropie]

dwi2response can already optionally output the SF mask, and could trivially be made to output directions for those voxels as well (note though that it will only have directions for voxels in the SF mask). It might just seem unusual to a standard user to have dwi2response (or whatever it's called) doing both a SF mask optimisation and RF estimation, and another command sh2response that just does RF estimation, but with more functionality. Why doesn't dwi2response also have the extra functionality? Why doesn't dwi2response just feed data to sh2response so that the extra functionality is accessible? Oh wait, I can do it if I trigger these optional outputs, and generate two response functions instead of one along the way... Removing the RF output from dwi2response (or making it a command-line option, with the SF mask the default argument output) may actually mean that the process makes more sense for the common user. Also, if dwi2response stays as it is, people will keep complaining about it not giving good RFs, whereas if they're forced to do two explicit steps, they're more likely to appreciate that finding good SF voxels isn't necessarily easy, and/or do manual SF mask editing or try to come up with a better estimation mechanism themselves.

RE amp2sh | sh2response, another way of putting it: If neither dwi2response nor sh2response existed, and you needed to write a command from scratch to do RF estimation, would you have it receive as input a DWI, or an SH image (which almost nobody deals with unless it's an FOD image)? The former makes a lot more sense to me, and I reckon anyone digging around in there would see an explicit SH conversion and turn their head sideways.

[jdtournier]

Yes, both good points. That said, until there is a universally accepted solution to this problem, I don't really think we should be thinking too hard about the commands themselves (in terms of their discoverability to users). If the recommended procedure is to use a script anyway, this is not an issue I'd be too worried about at this stage... We need a way for others to experiment with their response function estimation, and it might be messy for them if they want to delve into the script, but it's not something I'd be too worried about as long as each step is sensible once you see it in context.

And even if a single solution comes up that we all agree on, I'd still be of the opinion that it should be a script unless it is clear that the interdependencies between the various steps make that impractical and/or grossly inefficient. I really think there is value is providing users with a simple means of understanding the inner workings of the commands they're using, without having to look at the actual C++ code. I really think you'll also come to appreciate that when someone else can debug their own issues without requiring you to spend your time inspecting their data and fixing up your commands. Especially when the data they have is a little bit exotic and the procedure requires a tweak specifically for these data, which you don't want to include the main command since it'll only ever be used for these data... With a script, there is more scope for this to handled by the users without your involvement.

[jdtournier]

msdwi2fod now merged to updated_syntax - closing.