I would personally encourage to adopt the lower threshold for the template fixel masks that I have suggested as a docs change in this pull request, since 0.2 has been found to be consistently inappropriate in a lot of studies now, and this impact is very important on CFE, which affects other fixels even if they are in the mask at the higher threshold.

I don't see this change in the diff...

Also: This is still being applied as a peak amplitude threshold... has an integral threshold been tried instead?

The reason that I retained mtbin before, but made it provide the warning/error, was to in the mean time alarm users that they should probably switch.

I'm still personally for stripping it of its functionality for the sake of being safe for good FBA science, ...

(& various points from @jdtournier)

Consider instead having mtbin throw an exception immediately, unless an option -override is provided at the command-line. That gives users feedback of a very high priority, will stall any automated analysis, but still allows access to the old command if users wish to test the differences themselves.

I would however fix the basis function error in mtbin. That's fundamentally a bug, which would otherwise be fixed with a software update if the higher-level issues with the command were not also addressed. Leaving mtbin in place so that its output could be compared to mtnormalise, but leaving an outright bug in place, would be disingenuous, since it would not prove that a complete overhaul of the command was necessary over and above simply fixing the basis function bug. If one wants access to the current mtbin, including the basis function bug, it can always be accessed by checking out an older code version.

To be honest, I'm not sure labelling it as 3.0_RC2 achieves this enough, but I dont' think we want to label it as 3.1_RC1 either, when we haven't even released 3.0 yet...

3.0_RC2 makes sense to me. Along with an Announcement post.

Happy for the tag to come along with a few other changes as well (like the fixelstats -mask one).

That's probably the only one that's likely to go with it. We haven't even created a development branch yet for accumulating updates, so adding anything else would involve finishing and combining multiple pull requests; and I don't think it's worth delaying this trying to bundle other things to go with it.

After a long and intense strategic discussion, he convinced me (and at that point, I agreed) of trying to get rid of mtbin and making that very clear to users ASAP, so they would just as ASAP change their practices to do FBA with the new method.

If the issues are of the severity that you're describing, I'd have preferred that an announcement be sent out before a replacement was ready. If results generated using this command had the potential to be worse than what would have been produced using what the recommended pipeline was before mtbin was made available, then our recommendation should have reverted back to that pipeline, along with informing users of the issue and that developing a replacement would be pursued with relative urgency. They could then choose whether to continue using the method with stringent manual checking, perform FBA processing using the former pipeline, or wait for the update. As things have played out, mtbin is still sitting there on master with no warnings, which strikes me as inappropriate if the problems are indeed so severe. This also would have somewhat reduced the time pressure on getting the replacement working. Then again, maybe there was an attempt to do this but there was inadequate communication response; I don't know.

Command name (comments)

mtbin has never really worked for me, even if the acronym is OK in hindsight.

I don't see that it's any different to tcksift. If a method is conceptually unique enough to warrant an acronym, it's probably going to be difficult to encapsulate the operation of the method within a command name operator postfix of 10 characters or less without using said acronym.

I'm thinking of candidates like fodbiascorrect or equivalent, ...

... in terms of the name I personally think mtnormalise is a prime candidate ...

These demonstrate the issue above. For instance, we have dwibiascorrect that does bias field correction, and dwinormalise that does intensity normalisation. This command does both; so neither "operator" is sufficiently descriptive. I appreciate that the emphasis on the bias field correction aspect of this method may be reduced if dwibiascorrect is recommended to appear in the pipeline and/or a command-line option is added to this command to disable the bias field estimation, but it's still a fundamental part of the default operation of the command.

Command name (proposed)

Having gone through the command and compared it to the mtbin code, here's my unashamed opinion. Not looking to trigger or berate anyone, just laying down my reasoning of the situation. It's a difficult circumstance to deal with no matter who says what, so let's just get my bit over and done with.

If we consider those features that are common between this command and mtbin:

• Bias field correction;
• Combined with intensity normalisation;
• Of multi-tissue decomposed data;
• Based on objective of unity sum of tissue volume fractions in each voxel;
• Using 3D polynomial basis functions;
• And a single global normalisation factor per tissue;
• Incorporating voxel mask outlier rejection based on inter-quartile range.

Those features that differ between the two commands:

• Objective function (i.e. using a log-basis);
• Arrangement of various algorithm steps, i.e. ordering, nested iteration;
• Removal of stray degree-of-freedom in global tissue normalisation factors.

So while I appreciate that you want to emphasize these differences in order to justify the importance of the changes / the magnitude of your contribution / the claim of command authorship, ultimately this command has more in common with "the MTBIN method" than it has differences. Therefore, to demand absolute separation of this command from mtbin (both the command and the ISMRM abstract) is disrespecting the underlying intellectual property that belongs to another developer, yet is completely fundamental to its operation. This is most certainly not the kind of IP relationship we want occurring between developers (and the negative perception of such was compounded by pushing the change without independent review).

OK, that's the nasty bit done. Now onto some constructive input.

I advocate the following compromise:

Rename this new command to mtbin2.

This:

• Retains the reference to the fundamental concept from which the method originates;
• Concisely describes what the command provides to the user (within the limitations of the acronym at least);
• Is logically consistent for any users already using the mtbin command;
• Clearly indicates that the internal approach to solving the problem has changed significantly enough to justify a version increment (just like tcksift / tcksift2);
• Justifies the command authorship as it stands now;
• Is compatible with retaining the old command;
• Does not conflict with any existing command (i.e. as would occur if the new command were named mtnormalise).

Finally, some code snippets I didn't change but wanted to leave a comment on, putting them here so that maybe they provide a little cooldown period from the above:

summed_log_values.emplace_back (summed_log.value());

emplace_back() is intended for constructing objects in-place within the container; it's not just a "new" push_back(). E.g. Check my usage here.

const size_t max_inner_iter = DEFAULT_INNER_MAXITER_VALUE;

Provide access to this as a command-line option?

auto summed_log = ImageType::scratch (header_3D);
vector<float> summed_log_values;

These tripped me up initially. If I'm going to sum a set of values, then take the logarithm of the result, I would intuitively call that "log_summed". "summed_log" may be more consistent with what I see people doing with file names in their pipelines, i.e. appending the most recent processing step to the name. Not sure if this is a variable personal preference thing, I just found it weird.

INFO ("iteration: " + str(iter));

I wouldn't run INFO level outputs just spitting out iteration numbers, since it's not really giving much "information" to the user. So I've dropped them to DEBUG; but I reckon you could include INFO level outputs at the end of execution providing e.g. the final multiplicative factor per tissue.

outlier_rejection (3.f);
outlier_rejection(1.5f);

Worth having as command-line options?

output_headers[j].keyval()["lognorm_scale"] = str(lognorm_scale);

Is this going to need to be regressed against? Or should it be incorporated into scaling the output images?

See below an introduction into / overview of the structure of mtnormalise, which we claim to be "pretty damn robust", which means I haven't been able to make it fail at its task over a lot of datasets (I will detail this later). For people concerned with the difference relative to the "old" mtbin, I will of course assume you're fully up to date with the exact structure of that one. So here goes my fancy sketch:

Users INput a number of tissues (e.g. directly from MT-CSD) and a brain mask.

1. "Non-physical" voxels, defined as the sum of all tissues being zero or less, are detected and excluded from the mask (that the user provided at the input )forever. We don't want to (even aim to) fit zero, because if we would hypothetically succeed, we'd be dividing by zero when we do the correction, which gets us already into one of the problems of mtbin. When I refer to the "basic mask" beyond this point, this refers to the user input mask without these non-physical voxels. There's only 1 "basic mask", it never changes beyond this initial step.
2. The summed tissue image intensities within the basic mask are subjected to a weak outlier detection. Outlier detection is always in the log-domain, where we compute the quartiles and apply generous Tukey-style fences (Q1-3xIQR and Q3+3xIQR). Whenever mtbin computed outliers (but it never did this generous one at the start, only another one later on), it didn't do this in the log-domain, resulting almost always in a lower Tukey fence that was negative, which made that lower fence useless: it didn't even reject negative/zero values, let alone any very small ones (i.e. positive but close to zero). When I refer to the "(current) working mask" beyond this point, this refers to the "basic mask" without the current outlier voxels. The current working mask changes whenever outlier detection is performed. Outlier detection is always performed in the log-domain and using all current data in the basic mask. So whenever this happens, the outlier status of any voxel in the basic mask can change: voxels that weren't outliers before can become outliers, but voxels that were outliers before can also become non-outliers ("inliers"? ).
3. We compute the balance of the tissue types. In a sort of perfect world, this wouldn't be needed, but in our world we do this to essentially account for the intensity scale (local T2) of where we once happened to sample our response functions. In particular the WM one: we typically get it in areas that happen to have longer T2 than most of the other WM (e.g. CST, splenium, ..), so out of the box the sum of tissue types image typically shows a lot of the WM areas being lower intensity relative to GM/CSF. Tissue balancing means solving for N (equal to number of tissues) weights that result in a weighted sum of tissue maps being as constant over all voxels in the current working mask as possible. The weights output by this step are guaranteed (by our implementation) to have a geometric mean of 1 (equivalent: log-domain weights have an arithmetic mean of 0). So with multiplicative effects in the back of your mind, this operation doesn't affect the global intensity level of the weighted sum of tissues as a whole. mtbin had a mixed meaning at this step, since it didn't guarantee such a geometric mean of 1 (or any other mean of 1, for that matter). It mixed the balancing with a correction for the global intensity level. It solved for N degrees of freedom in this kind of step, while we solve for N-1 degrees of freedom.
4. The current weighted sum of tissues image is subjected to outlier detection with normal Tukey fences (Q1-1.5xIQR and Q3+1.5xIQR). Of course in the log-domain, using all current data in the basic mask. This results in a new current working mask. If the ("new")current working mask is different from the previous working mask, we go back to step 3. Steps 3 and 4 iterate in a tight loop to co-optimise the tissue balance and outlier definitions. This is all about getting the data right for input into the actual (one, and only) normalisation step that this command is all about. These steps are computationally cheaper compared to the normalisation step that will follow, so it's worth spending some tight looping on. Furthermore, there's only so much that can change here between these two steps; whereas the following intensity normalisation is much more flexible (polynomial in log-domain). Especially before the first time we will be heading into the normalisation step, we want to do an effort to not give it rather unbalanced tissues and/or a bad definition of outliers (either not rejecting important outliers, or rejecting too much valuable data. Should the subsequent normalisation go wonky because we don't do these steps decently enough, then any steps thereafter may become even more ill informed. So tight loop it is. We iterate until exact convergence (which is exactly established when the current working mask is not changing any more) or until a safety cap in iterations. This has been set at 7, which is a decent balance. Even if we reach 7, the result becomes already very stable, and we will revisit this tight loop many more times to give the algo further chances to fine-tune all the little details. We do need a cap, because due to the distribution of values being discrete samples, there exist a risk that steps 3 and 4 start to oscillate between 2 or more states. The risk is minor, but I searched and searched and searched until I stumbled upon a subject that could break things without such a cap within a particular iteration of the algo, where it kept on oscillating. mtbin doesn't do a tight loop.
5. After exiting the dashed box of co-optimising steps 3 and 4, we fit the normalisation field represented by a polynomial basis to the current weighted sum of tissues in the log-domain for all voxels in the current working mask with the goal of the normalisation field correcting towards a given value (the well-known 0.282... by default, user can change this for their own unusual uses). This finishes one "iteration", which started with the steps 3-4 tight looping dashed box. We go back to this dashed box, which now applies the current normalisation field for the purpose of establishing the tissue balance (3) and outlier detection (4). The default number of "iterations" (the big loop) is set to 15, the user can change this. In practice, something around 10 is quite safe for both clinical resolutions (i.e. numbers of voxels to which everything is fit) as well as upsampled or HCP-range resolutions to ensure reasonable convergence, but since the algorithm is relatively fast anyway, we play it extra safe and have set the default to 15. There's no need to waste time on computing an overall convergence criterium, and any criterium you can think of still has weaknesses that may catch you unaware in specific settings. So just a fixed number of iterations that is on the safe side, was in our experience the best guarantee for the "ultimate robustness" we are so proud of. At the more or less similar step "5" in mtbin, it works of course not in the log-domain, but it also does a re-normalisation towards an arithmetic average of 1 bias-style of field. It's very hard to explain without sketches and 5 hands (ask @rtabbara about how we communicated this among ourselves if you're interested), but this mtbin practice essentially spreads out the normalisation across the balancing step and the polynomial fitting step. This is annoying to get a grip on what information is fitted where, as when the polynomial fit changes to accommodate a local intensity change, this influences the rest of the values too, which then has to be countered by the balancing-and-somewhat-global-intensity-level fitting step of mtbin. Parts of the normalisation are bouncing back and forth between these steps. Furthermore mtbin uses the tissue balance-and-somewhat-intensity-level weights as its convergence criterium for the whole algorithm. This is unreliable both ways! The algorithm can appear to be quite far from converging since this bouncing back and forth of information is happening without improving fits a lot or even at all. The "at all" may surprise you, but this is a further consequence of the faulty arithmetic: solving least squares in the non-log-domain (additive mindset goes here), but then dividing by this solution (multiplicative mindset goes here) is inconsistent with each other. You could do mtbin's sort of biasy-normalisation step twice in a row (second time on the corrected output of the first) without any other actions in between, and still get a spatially varying result the second time. mtnormalise's normalisation step doesn't suffer from this, as it fits with an additive mindset in the log-domain, and corrects with an additive mindset in the log-domain (or multiplicative in the non-log-domain, whichever domain suits your mind the best to visualise this mentally). So that's mtbin potentially failing to detect convergence for quite a while. The other way around is possible too, since mtbin does its own "kind-of similar (but each very different as should be clear by now) to steps 3-4-5" steps in one big loop (no tight looping any 2 steps). So even though mtbin's balancy-globally-somethings weights may at one point not vary a lot or almost at all between 2 of its bigger iterations, this doesn't guarantee that the current sort of mutual state between mtbin's outliers and biasy-normalisation field hasn't changed. It's not because it (again "sort of") global intensity level has been settled in it's own kind-of step 3, that the spatially varying aspect still isn't changing due to the outlier definition changing in turn. In fact, these are probably the 2 steps that risk influencing (as in "co-optimise") each other the most. This is why getting the outlier detection stable is quite important with respect to the normalisation field. This is one of the primary reasons for the tight loop in mtnormalise. I hope I got that all across in this bulk of text. You would've certainly gotten it if you were here with us for several hours at the white board and with the data that at times got us back to the white board...

I also hope the above makes it clear it's not easy (in practice, for me) to talk about both algorithms with respect to each other, since it's not just the order or looping structure of some boxes that is different, or even on top of that the log-domain in certain crucial definitions of those boxes, but also in a way the meaning of what our boxes are after, versus mtbin. To get a grip on making this optimisation stable and reliable for many different datasets (i.e. essentially initialisations), we consciously found our way towards this clear separation of "tissue balance" in step 3, and "normalisation" (but then all of it) in step 5. This means that "switching off the bias field correction", where you would/may be thinking of the bias field and global normalisation as something separate in your mind, becomes in our case just setting the polynomial basis for step 5 to an order of zero; i.e. just fitting the constant term. All the rest of the theory above holds as usual, the normalisation just doesn't vary spatially. In practice, you'll see and should at this point understand that such an order zero step 5 then can't influence the outlier definition and tissue balance, if these would've already converged among themselves in the tight loop. However, since this may not have happened exactly (yet probably still close to) during the first capped 7 tight loops, then still the next 7 ones will in principle just continue among themselves, not being influenced at that stage by anything the normalisation step 5 delivered. Each time step 5 is encountered in such a zero order scenario, the normalisation that it fits may change still due to the current state of the 3-4 tight loop. So co-optimising 3-4 still influences the final main 5 output, even though 5 doesn't influence 3-4 in this scenario. In this case, the algorithm thus essentially reduces to the 3-4 loop without a cap, and delivers the 5-fit that belongs with that. There is still an overall cap of 7x15 runs of 3 and 4 of course. But in practice, you'll suddenly see a 3-4 loop below 7 iterations, and then nothing changes any more and everything has exactly converged. Setting the order to anything above zero (order three is a wise default) allows spatial variance, and is virtually impossible to exactly converge, due to the continuous nature of step 3 and step 5 co-optimising until numerical precision can't tell the changes any more. But, and this is the most important but(t) end of the business here, you'll see that this algorithm and all of the above specific choices lead to virtual total convergence. I will detail later on (but I'm close to finishing for the night now) what this was all tested on, but you can easily give this a go yourself on any dataset you've got lying around, which will probably give more insight via experience than all of my words here can. What you can do is: run it on a dataset of your liking (maybe not neonatal to begin with, but it'd be cool to give it a go on that too if you're keen), and specify the -info option for some insights into the run as it happens. Be sure to also specify the -check_mask and -check_norm options, which give you the (final) "inlier" mask of the run and the actual normalisation field (your input data divided by the normalisation field is your output data). Watch an interesting dynamic unfold (add some national geographic music if desired). But then the most interesting bit maybe: give the outputted corrected tissues as the input to a new run of mtnormalise, all the above options for insight added again. And take a look at the normalisation field (within your brain mask, but even typically far beyond that, often easily the whole field of view) of that second run if that is very close to essentially constant 1, you'll know the first run had essentially converged, and an equally intensive complete second run wasn't able to do anything of significance beyond that. Also take a look at the inlier mask (i.e. the final working mask): you'll see that what is typically removed as outliers in a healthy developed human is areas affected by potential other not fully corrected artefacts (more extreme effects of ringing in some cases, hyper or hypo intense regions due to EPI distortions at extremities, ...), some non-brain areas if some bits of junk were accidentally included in the brain mask, and areas of low signal (due to iron, e.g. dentate nucleus, etc..) or high/long T2 relative to the bulk of a given tissue type's other T2 (e.g. core CST bits, maybe bits of splenium).

Finally to add some insight to how step 2 came to be added to our design at some stage (you can track this down in the commits at some point): before that one, the chicken and egg problem that kept us up all night was the order of 3 and 4. Once the thing's running smoothly, that order doesn't really matter, it's then just about getting the pair 3-4 right, and iterate optimising that pair with 5. But at the start, we face a dilemma. In a decently behaved scenario, no worries. But:

• if (some) extreme outliers in the data somewhere, this risks getting the tissue balance completely wrong if 3 goes first, with the subsequent risk of one or more tissue types almost entire "area" becoming an outlier subsequently, and then all risks going wonky. Note this is probably quite a rare scenario, but we didn't stop at dreaming those up too.
• if there is one (or more) tissue types very unbalanced at the input (due to very weird calibration of responses for a very weird scenario), this risks getting the outliers wrong if 4 goes first, resulting in the same kind of risk to wonky behaviour as in the aforementioned case, but now due to the opposite reason at the input. Again a very rare scenario (or at least very non-standard) indeed, but well, you never know.

So both scenarios are rare, and it's a bit hard to guess which one is more likely when we integrate over all use cases done by actual users from now until the end of times. So it still kept us up at night, and deciding the order of 3-4, while irrelevant for most scenarios, became a small obsession (from my part mostly, @rtabbara is probably more reasonable on this front 😉 ) to get right. Adding step 2 is the compromise between both: we do detect outliers first, but we do it much more relaxed in that step. So the tissue types have to be unbalanced up too extremely silly extents for the above horror scenario to unfold. That's all the reasoning there is to that specific one, which, as mentioned, was a final touch to the algorithm when the rest was established at that point.

Ok, the above documents mtnormalise and our journey to it better than mtbin ever was. I'm going home for the day now, but I'll list some of what we did for testing later (but it'll be a more concise overview then) and conjure up another example on an actual dataset that shows the essence of mtbin failing (i.e. the thing bugging us so much we felt we had to embark on this crazy journey). I will then also reply to some of the other things that have been commented in this pull request.

OK, I just rebased this to merge onto the dev branch. There'll be a few more issues going into it soon.

As to this branch, I guess the only thing left to discuss is the name - feels like we agree on just about everything else.

• I like @Lestropie's idea of an -override option to allow mtbin to proceed, but with a default behaviour of throwing an exception. That really should cover all bases.

• we'll need to sort out an announcement message to go with the new updated tag name - we'll get the process underway on the staff forum when things have settled down a bit.

• naming-wise, I appreciate @Lestropie's point about a lot of commands already explicitly being named as acronyms, but the distinction here is that many of these are named after algorithms that have been published under these acronyms. The tcksift / tcksift2 commands are a case in point: they implement the algorithm published under that exact acronym, it would make no sense to name them anything else. In contrast, we never proposed a BIN algorithm, and there's no mention of such in Dave's abstract. So in this case, the 'bias-field and intensity normalisation' acronym comes entirely from internal decisions, rather then referring to published literature. So I don't see that there's any reason to use that name, and I would still rather use a more descriptive term here. The only issues with the proposed mtnormalise as far as I'm concerned is that there is already an existing command with that name (although I don't expect anyone is using it), and that it breaks the logical connection that would have existed with a transition from mtbin to mtbin2 - but that really only impacts users already familiar with mtbin, and there are most likely not too numerous yet given how recently it was released. So I'd vote mtnormalise, but I'm admittedly not all that fussed about it - mtbin2 would be fine by me if that was the consensus.

We have done various testing using:

• 40 subjects' 3-tissue results, mixed healthy subjects and several variants of MS
• 97 subjects' 3-tissue results, mixed healthy subjects and cases of epilepsy, many related to DEPDC5
• several randomly selected elderly peoples' 3-tissue results, some with MCI and AD; all with WMH lesions, some with large amounts of very large lesions, some with massive calcifications of the choroid plexus, many with visible atrophy, sometimes extreme atrophy and massive ventricles
• a young child's 3-tissue result, both original and upsampled resolutions
• HCP data (multi-tissue result)
• "low quality" data (low b-value, limited number of directions) (3-tissue result)
• 2-tissue WM-CSF results of some of the above, just to see where I could push it. Works perfectly by the way, almost always fitting a normalisation field that is very similar in contrast to the one produced when informed by 3-tissue CSD. Single-tissue results (WM only) don't do well of course, we're missing too much of the (signal) picture there obviously.

For most, both mtbin as well as a few stages of development of mtnormalise have been run. We did mrstats on the results, checking for values in realistic ranges. For the first 2 groups, I directly compared the spatially variant difference between the spatial contrasts of the actual "bias" fields, by dividing one by the other. This was just looking at contrast (and variance of that contrast), as mtbin's output tries to filter out the global part from the bias field that is output via an option, whereas mtnormalise outputs the complete normalisation field, that includes this scale in it (i.e. it directly relates input to output).
Next, I'll briefly produce some output on a typical mtbin ultimate fail case (there's also less ultimate fails, "average" fails one could call them, but I won't be reproducing all of that again). The ultimate fail case is probably the best one to keep here for future reference, should anyone want to re-understand the reason of this effort.

Ultimate fail case. I picked one of the 40 subjects of the above first bullet for this. There were 3 such "ultimate fail" cases among those 40.

Input:
WM

GM

CSF

Everything outside the mask is black, inside the mask is jet colour scale. As you can see, there is a little bit of crap attached to the mask at the front. This happens every once in a while; but the bit was really small in this case.

mtbin run:

thijs@bri-pub-59:~/Documents/mtbin_fail$ mtbin fod_wm.mif binwm.mif fod_gm.mif bingm.mif fod_csf.mif bincsf.mif -mask mask.mif -info
mtbin: [WARNING] This command is deprecated and may produce inappropriate results in several cases. Please use the new mtnormalise.
mtbin: [.   ] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] opening image "fod_wm.mif"...
mtbin: [INFO] opening image "fod_wm.mif"...
mtbin: [INFO] opening image "fod_gm.mif"...
mtbin: [INFO] opening image "fod_gm.mif"...
mtbin: [INFO] opening image "fod_csf.mif"...
mtbin: [INFO] opening image "fod_csf.mif"...
mtbin: [INFO] opening image "mask.mif"...
mtbin: [  . ] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] iteration: 1
mtbin: [ .  ] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] scale factors: 1.7575459653980563 1.2945077816176735 1.6570568353748945
mtbin: [ .  ] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] iteration: 2
mtbin: [  . ] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] scale factors: 1.8697080019987342 1.2526105713984046 1.5752244097784691
mtbin: [ .  ] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] percentage change in estimated scale factors: 4.8522324151506906
mtbin: [  . ] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] iteration: 3
mtbin: [.   ] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] scale factors: 1.8879995524327986 1.2414655745914718 1.5499852244482419
mtbin: [ .  ] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] percentage change in estimated scale factors: 1.1567705867484146
mtbin: [INFO] rejecting outliers
mtbin: [  . ] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] iteration: 4
mtbin: [   .] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] scale factors: 1.8968117770775175 1.2633873872754533 1.5427502634419703
mtbin: [  . ] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] percentage change in estimated scale factors: 0.89977550588256838
mtbin: [INFO] rejecting outliers
mtbin: [.   ] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] iteration: 5
mtbin: [ .  ] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] scale factors: 1.8820840428746743 1.2666131760713555 1.5488073825435695
mtbin: [.   ] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] percentage change in estimated scale factors: 0.47479786711318139
mtbin: [INFO] rejecting outliers
mtbin: [   .] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] iteration: 6
mtbin: [ .  ] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] scale factors: 1.8758530823174706 1.2677816438422407 1.5541932766251769
mtbin: [   .] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] percentage change in estimated scale factors: 0.25702100915106807
mtbin: [INFO] rejecting outliers
mtbin: [.   ] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] iteration: 7
mtbin: [ .  ] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] scale factors: 1.8735276544870516 1.2681914502510261 1.5560461165498927
mtbin: [.   ] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] percentage change in estimated scale factors: 0.091835546910303559
mtbin: [INFO] creating image "binwm.mif"...
mtbin: [INFO] creating image "bingm.mif"...
mtbin: [INFO] creating image "bincsf.mif"...
mtbin: [done] performing intensity normalisation and bias field correction...

Output:
WM

GM

CSF

Look at the range for insight into the problem at the crosshairs.

Maybe we could iron that out, one might think? Let's try and feed these outputs into mtbin again:

thijs@bri-pub-59:~/Documents/mtbin_fail$ mtbin binwm.mif binwm2.mif bingm.mif bingm2.mif bincsf.mif bincsf2.mif -mask mask.mif -info
mtbin: [WARNING] This command is deprecated and may produce inappropriate results in several cases. Please use the new mtnormalise.
mtbin: [.   ] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] opening image "binwm.mif"...
mtbin: [INFO] opening image "binwm.mif"...
mtbin: [INFO] opening image "bingm.mif"...
mtbin: [INFO] opening image "bingm.mif"...
mtbin: [INFO] opening image "bincsf.mif"...
mtbin: [INFO] opening image "bincsf.mif"...
mtbin: [INFO] opening image "mask.mif"...
mtbin: [   .] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] iteration: 1
mtbin: [  . ] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] scale factors:    1.4802299640746948 -0.063235442264358788 -0.044279212053310965
mtbin: [ .  ] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] iteration: 2
mtbin: [ .  ] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] scale factors:   1.5043068175025223 -0.12061755527460417 -0.14954650969024424
mtbin: [.   ] performing intensity normalisation and bias field correction...... 
mtbin: [INFO] percentage change in estimated scale factors: -108.95076241046257
mtbin: [INFO] creating image "binwm2.mif"...
mtbin: [INFO] creating image "bingm2.mif"...
mtbin: [INFO] creating image "bincsf2.mif"...
mtbin: [done] performing intensity normalisation and bias field correction...

Yeah, no that is. This independently also shows mtbin can't handle outliers at the intput correctly. I could show the results, but they're not very interesting: all the outputs are only filled with NaNs.

Next up, mtnormalise on the same original inputs:

thijs@bri-pub-59:~/Documents/mtbin_fail$ mtnormalise fod_wm.mif normwm.mif fod_gm.mif normgm.mif fod_csf.mif normcsf.mif -mask mask.mif -info
mtnormalise: [.   ] loading input images... 
mtnormalise: [INFO] opening image "fod_wm.mif"...
mtnormalise: [INFO] opening image "fod_gm.mif"...
mtnormalise: [INFO] opening image "fod_csf.mif"...
mtnormalise: [INFO] opening image "mask.mif"...
mtnormalise: [done] loading input images
mtnormalise: [INFO] Iteration: 1
mtnormalise: [INFO] Balance factors (1):  1.1269888270383537 0.83458131093149701  1.0631920838588858
mtnormalise: [INFO] Balance factors (2):  1.1302228674059285 0.85288434429789217  1.0373988626651387
mtnormalise: [INFO] Balance factors (3):  1.1318211499078685 0.85719413069624539  1.0307254657088498
mtnormalise: [INFO] Balance factors (4):   1.132379985930005 0.85819950072061013  1.0290099106320219
mtnormalise: [INFO] Balance factors (5):  1.1325249810738918 0.85840185774959132  1.0286356237439809
mtnormalise: [INFO] Balance factors (6): 1.1325792406890618 0.8584543080499426 1.0285234987419798
mtnormalise: [INFO] Balance factors (7): 1.1325929315584391 0.8584708076848907   1.02849129809756
mtnormalise: [  7%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 2
mtnormalise: [INFO] Balance factors (1):   1.189519082947619 0.83580783339757081  1.0058243592202578
mtnormalise: [INFO] Balance factors (2):  1.1886476197216203 0.83561858207473494  1.0067897512421515
mtnormalise: [INFO] Balance factors (3):  1.1884634705477113 0.83556803140952518  1.0070066693396178
mtnormalise: [INFO] Balance factors (4):  1.1884252634695098 0.83554681899202565  1.0070646101019718
mtnormalise: [INFO] Balance factors (5):  1.1884175964463488 0.83554592377551717  1.0070721861256011
mtnormalise: [INFO] Balance factors (6):  1.1884149882788633 0.83554631188080863   1.007073928528333
mtnormalise: [ 13%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 3
mtnormalise: [INFO] Balance factors (1):  1.2090475033398744 0.83468189433174322 0.99091326547982295
mtnormalise: [INFO] Balance factors (2):  1.2087350200479985 0.83251386046932874 0.99375064224789567
mtnormalise: [INFO] Balance factors (3):  1.2086625182978181 0.83212525255797698  0.9942743683453853
mtnormalise: [INFO] Balance factors (4):  1.2086564552478103 0.83205824492942537 0.99435942765899143
mtnormalise: [INFO] Balance factors (5):  1.2086571504831494 0.83204635838830765 0.99437306101367684
mtnormalise: [INFO] Balance factors (6):  1.2086560776783439 0.83204611908353054 0.99437422961246202
mtnormalise: [INFO] Balance factors (7):  1.2086533781067477 0.83204400658859967 0.99437897522809326
mtnormalise: [ 20%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 4
mtnormalise: [INFO] Balance factors (1):  1.2184027209129167 0.83336970873486638 0.98485303841823424
mtnormalise: [INFO] Balance factors (2):  1.2182993136491886  0.8328864751036873 0.98550808297715209
mtnormalise: [INFO] Balance factors (3):  1.2182862309769167   0.832791435342223 0.98563113522808043
mtnormalise: [INFO] Balance factors (4):  1.2182925748411273 0.83276733872201092 0.98565452254993668
mtnormalise: [INFO] Balance factors (5):  1.2182956806284044 0.83276460755591797 0.98565524241347091
mtnormalise: [INFO] Balance factors (6):  1.2182983504896874 0.83276306039989334 0.98565491358607527
mtnormalise: [INFO] Balance factors (7):  1.2183031317952797  0.8327601763522845 0.98565445886797276
mtnormalise: [ 27%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 5
mtnormalise: [INFO] Balance factors (1):  1.2227214440811764 0.83397267268925035 0.98066494451266284
mtnormalise: [INFO] Balance factors (2):  1.2225272066619066 0.83388167171442973 0.98092779074452063
mtnormalise: [INFO] Balance factors (3):    1.22247697345839 0.83386824121334902 0.98098389810114606
mtnormalise: [INFO] Balance factors (4):  1.2224696638697587 0.83386603327920883 0.98099236125149825
mtnormalise: [INFO] Balance factors (5):  1.2224584177912003 0.83386920976458945 0.98099764898796527
mtnormalise: [INFO] Balance factors (6):  1.2224512288265836 0.83386965663452239 0.98100289229884696
mtnormalise: [INFO] Balance factors (7):  1.2224487270254882 0.83386856356171291 0.98100618591504096
mtnormalise: [ 33%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 6
mtnormalise: [INFO] Balance factors (1):  1.2245163590437684 0.83456386354203682 0.97853380262124257
mtnormalise: [INFO] Balance factors (2):  1.2244056351185353 0.83447689182653495 0.97872428715799964
mtnormalise: [INFO] Balance factors (3):   1.224394237060368  0.8344609937445937 0.97875204498792623
mtnormalise: [INFO] Balance factors (4):  1.2243927643036279 0.83445914021400724 0.97875539631832853
mtnormalise: [ 40%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 7
mtnormalise: [INFO] Balance factors (1):  1.2253990816148794 0.83480589634295954 0.97754541169176445
mtnormalise: [INFO] Balance factors (2):  1.2253170399018469 0.83482801444638066 0.97758496271064765
mtnormalise: [INFO] Balance factors (3):  1.2252950572841252  0.8348338886446145 0.97759562247191811
mtnormalise: [INFO] Balance factors (4):   1.225292342857758 0.83483296412246799 0.97759887079313634
mtnormalise: [INFO] Balance factors (5):  1.2252934583496049 0.83483184566417323 0.97759929052725891
mtnormalise: [INFO] Balance factors (6):  1.2252949255898011 0.83483158046210271    0.97759843044729
mtnormalise: [ 47%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 8
mtnormalise: [INFO] Balance factors (1):  1.2257837234164528 0.83497010609315991 0.97704647616654072
mtnormalise: [INFO] Balance factors (2):   1.225750362439149 0.83497203975047074  0.9770708054851992
mtnormalise: [INFO] Balance factors (3):  1.2257394892765092 0.83497274588621773 0.97707864646820231
mtnormalise: [INFO] Balance factors (4):   1.225736807239481 0.83497427534039859  0.9770789946597529
mtnormalise: [INFO] Balance factors (5):  1.2257301896214523 0.83497878037322271 0.97707899808332421
mtnormalise: [ 53%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 9
mtnormalise: [INFO] Balance factors (1):  1.2259900791111886 0.83504533890846699 0.97679401056865123
mtnormalise: [INFO] Balance factors (2):  1.2259909721429323 0.83501525240982388 0.97682849397034566
mtnormalise: [INFO] Balance factors (3):  1.2259899425834846 0.83501072262069875 0.97683461341999678
mtnormalise: [INFO] Balance factors (4):  1.2259880848662925 0.83500820228226003 0.97683904202271676
mtnormalise: [INFO] Balance factors (5):  1.2259891318288583 0.83500617969533952 0.97684057396618962
mtnormalise: [ 60%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 10
mtnormalise: [INFO] Balance factors (1):   1.226124239522588 0.83503260996137918 0.97670201972269999
mtnormalise: [INFO] Balance factors (2):  1.2261396166638276 0.83501532025075009 0.97670999400567193
mtnormalise: [ 67%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 11
mtnormalise: [INFO] Balance factors (1): 1.2262180084811556 0.8350208610216211 0.9766410726415935
mtnormalise: [INFO] Balance factors (2):  1.2262061250255982  0.8350165031538872 0.97665563454441184
mtnormalise: [INFO] Balance factors (3):   1.226200222704223 0.83501964353846392 0.97665666260743356
mtnormalise: [INFO] Balance factors (4):  1.2261977532838548 0.83502150904643213 0.97665644754376979
mtnormalise: [INFO] Balance factors (5):  1.2261950458787421 0.83502302710216947 0.97665682842800339
mtnormalise: [ 73%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 12
mtnormalise: [INFO] Balance factors (1):  1.2262470345934766 0.83502717022846329 0.97661057585141031
mtnormalise: [INFO] Balance factors (2): 1.2262289713515675 0.8350223142642299 0.9766306414694621
mtnormalise: [ 80%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 13
mtnormalise: [INFO] Balance factors (1):   1.226255452141173 0.83501758405893756 0.97661508358263949
mtnormalise: [INFO] Balance factors (2):  1.2262411846326819 0.83502051005848921  0.9766230244506735
mtnormalise: [INFO] Balance factors (3):  1.2262379748212491 0.83502195015901415 0.97662389655397019
mtnormalise: [ 87%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 14
mtnormalise: [INFO] Balance factors (1):  1.2262532024969077  0.8350164013229906 0.97661825855095852
mtnormalise: [INFO] Balance factors (2):  1.2262556960950668 0.83501608283547835 0.97661664508791846
mtnormalise: [ 93%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 15
mtnormalise: [INFO] Balance factors (1):  1.2262663602035184 0.83501636748421404 0.97660781911844519
mtnormalise: [INFO] Balance factors (2):  1.2262501899812361 0.83501894657178843 0.97661768092889978
mtnormalise: [INFO] Balance factors (3):  1.2262494086349061 0.83501894335555538 0.97661830697555219
mtnormalise: [100%] performing log-domain intensity normalisation
mtnormalise: [.   ] writing output images... 
mtnormalise: [INFO] creating image "normwm.mif"...
mtnormalise: [ .  ] writing output images... 
mtnormalise: [INFO] creating image "normgm.mif"...
mtnormalise: [   .] writing output images... 
mtnormalise: [INFO] creating image "normcsf.mif"...
mtnormalise: [done] writing output images

Output:
WM

GM

CSF

All nice and good. We may wonder if we really converged as good as it seems. Let's throw the outputs at mtnormalise again and see what happens:

thijs@bri-pub-59:~/Documents/mtbin_fail$ mtnormalise normwm.mif normwm2.mif normgm.mif normgm2.mif normcsf.mif normcsf2.mif -mask mask.mif -check_norm checknorm2.mif -info
mtnormalise: [.   ] loading input images... 
mtnormalise: [INFO] opening image "normwm.mif"...
mtnormalise: [INFO] opening image "normgm.mif"...
mtnormalise: [INFO] opening image "normcsf.mif"...
mtnormalise: [INFO] opening image "mask.mif"...
mtnormalise: [done] loading input images
mtnormalise: [INFO] Iteration: 1
mtnormalise: [INFO] Balance factors (1):  1.2237739278136768 0.82050066976746128  0.9959094737156714
mtnormalise: [INFO] Balance factors (2):  1.2257145712891209 0.83249291833359973 0.98000908812744281
mtnormalise: [INFO] Balance factors (3):  1.2261619963657453 0.83458941458993841  0.9771905904550684
mtnormalise: [INFO] Balance factors (4):  1.2262260132714449 0.83494538097694648 0.97672298606466479
mtnormalise: [INFO] Balance factors (5):  1.2262463093283977  0.8350063269265936 0.97663553148925997
mtnormalise: [INFO] Balance factors (6):  1.2262537229259864 0.83501552795893264 0.97661886553813115
mtnormalise: [INFO] Balance factors (7):  1.2262555482368245 0.83501615484556668 0.97661667862400914
mtnormalise: [  7%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 2
mtnormalise: [INFO] Balance factors (1):  1.2262610988437703 0.83501380510768686 0.97661500621697661
mtnormalise: [INFO] Balance factors (2): 1.2262546106836345 0.8350113610405886 0.9766230320785998
mtnormalise: [INFO] Balance factors (3):  1.2262533893195933 0.83501104839232843 0.97662437047917972
mtnormalise: [ 13%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 3
mtnormalise: [INFO] Balance factors (1):  1.2262562362164768 0.83500958317172336   0.976623816844851
mtnormalise: [INFO] Balance factors (2):  1.2262542731185193 0.83500464344311065 0.97663115784190646
mtnormalise: [ 20%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 4
mtnormalise: [INFO] Balance factors (1):  1.2262557573489861 0.83500342658315774 0.97663139900497697
mtnormalise: [INFO] Balance factors (2):  1.2262599859949519 0.83500023957075176 0.97663175875503683
mtnormalise: [ 27%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 5
mtnormalise: [INFO] Balance factors (1):  1.2262628376731874 0.83499971926498384 0.97663009615292562
mtnormalise: [ 33%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 6
mtnormalise: [INFO] Balance factors (1):  1.2262634541779471 0.83499967639536876 0.97662965529252066
mtnormalise: [ 40%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 7
mtnormalise: [INFO] Balance factors (1): 1.2262635910829818 0.8349996754276654 0.9766295473894725
mtnormalise: [ 47%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 8
mtnormalise: [INFO] Balance factors (1):   1.226263621773559 0.83499967626895777 0.97662952196267638
mtnormalise: [ 53%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 9
mtnormalise: [INFO] Balance factors (1):  1.2262636287348627  0.8349996765700366 0.97662951606635839
mtnormalise: [ 60%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 10
mtnormalise: [INFO] Balance factors (1):  1.2262636303949705 0.83499967660579266 0.97662951470238424
mtnormalise: [ 67%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 11
mtnormalise: [INFO] Balance factors (1):  1.2262636307923214 0.83499967661374208 0.97662951437662404
mtnormalise: [ 73%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 12
mtnormalise: [INFO] Balance factors (1):  1.2262636308675581 0.83499967661522545 0.97662951431496892
mtnormalise: [ 80%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 13
mtnormalise: [INFO] Balance factors (1):  1.2262636308849626 0.83499967661103891 0.97662951430600464
mtnormalise: [ 87%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 14
mtnormalise: [INFO] Balance factors (1):  1.2262636308882184 0.83499967661029362 0.97662951430428313
mtnormalise: [ 93%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 15
mtnormalise: [INFO] Balance factors (1):  1.2262636308886687 0.83499967661064367 0.97662951430351486
mtnormalise: [100%] performing log-domain intensity normalisation
mtnormalise: [INFO] creating image "checknorm2.mif"...
mtnormalise: [.   ] writing output images... 
mtnormalise: [INFO] creating image "normwm2.mif"...
mtnormalise: [ .  ] writing output images... 
mtnormalise: [INFO] creating image "normgm2.mif"...
mtnormalise: [  . ] writing output images... 
mtnormalise: [INFO] creating image "normcsf2.mif"...
mtnormalise: [done] writing output images

It definitely seems to ultimately settle, with from outer iteration 5 onwards literally no more changes in the outlier mask. Note that in between each of those lines, both the spatially variant normalisation as well as subsequently the rebalancing are given a go (see diagram in other post), and still this doesn't cause changes in the outlier mask. The normalisation field as well as the rebalancing are still co-optimising though, but the changes are in the realm of numerical insignificance. Interesting observation: why does the algorithm still seem to go so happily in those first 4 iterations? Well, simple explanation really: the initial weaker outlier rejection sets the initial rebalancing off to a slightly different start, compared to what it was "used to" and had settled on during the last iterations of the previous run. But it quickly finds its way again during iteration 1 itself, and the subsequent iterations are just for obsession's sake. To convince us of that, I've output the actual normalisation field of this run as well.

Here's the stats on that one:

thijs@bri-pub-59:~/Documents/mtbin_fail$ mrstats checknorm2.mif -mask mask.mif 
      volume       mean     median      stdev        min        max      count
       [ 0 ]   0.999989 0.999992967 6.86213e-05   0.998896    1.00039     792008

Beautiful. And just because we're overly confident, let's check that for the entire field of view, even spatially quite far beyond the actual mask:

thijs@bri-pub-59:~/Documents/mtbin_fail$ mrstats checknorm2.mif 
      volume       mean     median      stdev        min        max      count
       [ 0 ]   0.999849 0.999973893 0.000519785   0.995579    1.00179    4423680

Booyah. Note this is far from trivial, and actually quite exceptional: we all know how a polynomial behaves beyond the domain of data that it was fit to. The edge voxels, especially those in that crappy extra bit to the mask, risk having all the influence to still have this polynomial go to very low or high values by the time it reaches the edge of the field of view. But it doesn't; it's marvellously flat.

Finally, just for fun, let's try and give mtnormalise that first output from mtbin, the one that really tripped up mtbin on its second run. We're not expecting mtnormalise to iron out those extreme values again of course. On the contrary, it should avoid modelling those and tripping over them. They will remain in the result though, but we hope they don't affect the rest of the process too much. So here goes:

thijs@bri-pub-59:~/Documents/mtbin_fail$ mtnormalise binwm.mif binnormwm.mif bingm.mif binnormgm.mif bincsf.mif binnormcsf.mif -mask mask.mif -check_norm checkbinnorm.mif -info
mtnormalise: [.   ] loading input images... 
mtnormalise: [INFO] opening image "binwm.mif"...
mtnormalise: [INFO] opening image "bingm.mif"...
mtnormalise: [INFO] opening image "bincsf.mif"...
mtnormalise: [INFO] opening image "mask.mif"...
mtnormalise: [done] loading input images
mtnormalise: [INFO] Iteration: 1
mtnormalise: [INFO] Balance factors (1):  1.2070962510994618 0.80854625008047976  1.0245973578333272
mtnormalise: [INFO] Balance factors (2):  1.2118917428063674 0.81934780480316716  1.0070890619198065
mtnormalise: [INFO] Balance factors (3):  1.2128426508061805 0.82115370149778588  1.0040863986504014
mtnormalise: [INFO] Balance factors (4):  1.2130369508556027 0.82144374255020602  1.0035710945748428
mtnormalise: [INFO] Balance factors (5):   1.213071079314588 0.82147463960051914  1.0035051152301673
mtnormalise: [INFO] Balance factors (6):  1.2130776400884289 0.82148215316943263  1.0034905095360391
mtnormalise: [INFO] Balance factors (7):  1.2130771015159643 0.82148307230614359  1.0034898322771468
mtnormalise: [  7%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 2
mtnormalise: [INFO] Balance factors (1):  1.2232968427959132 0.82370735374196868 0.99241930074854945
mtnormalise: [INFO] Balance factors (2):  1.2215742577226556 0.82504625680864907  0.9922059555657825
mtnormalise: [INFO] Balance factors (3):   1.221194896453287 0.82527544104788075 0.99223855446785247
mtnormalise: [INFO] Balance factors (4):  1.2211061666752951 0.82532791606023848 0.99224756202815101
mtnormalise: [INFO] Balance factors (5):  1.2210819049399824 0.82534571450270722  0.9922458789647759
mtnormalise: [INFO] Balance factors (6):  1.2210801904210222 0.82534737969880723 0.99224527024911613
mtnormalise: [INFO] Balance factors (7):  1.2210773769433891 0.82534891818400347 0.99224570688465985
mtnormalise: [ 13%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 3
mtnormalise: [INFO] Balance factors (1):  1.2251096161156256 0.82826067533902314 0.98550313262162603
mtnormalise: [INFO] Balance factors (2):    1.22475513567079 0.82858098626971921 0.98540728264610478
mtnormalise: [INFO] Balance factors (3):   1.224671058640491 0.82864492902010389 0.98539888895466488
mtnormalise: [INFO] Balance factors (4):  1.2246585400718677 0.82865323849571437 0.98539908041683788
mtnormalise: [ 20%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 4
mtnormalise: [INFO] Balance factors (1):  1.2263728607934417 0.83067234011510338 0.98162976667701618
mtnormalise: [INFO] Balance factors (2):  1.2263980368494447 0.83076543766418487 0.98149961382917095
mtnormalise: [INFO] Balance factors (3):  1.2264022823761807 0.83077484451876227 0.98148510262892907
mtnormalise: [ 27%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 5
mtnormalise: [INFO] Balance factors (1):  1.2270299854985087 0.83218025099967663 0.97932630293678147
mtnormalise: [INFO] Balance factors (2):  1.2271915963222202 0.83218431920901448 0.97919254699397795
mtnormalise: [INFO] Balance factors (3):  1.2272260639880685 0.83218847457283751 0.97916015629455222
mtnormalise: [INFO] Balance factors (4):  1.2272297142025377  0.8321868807048467 0.97915911928341115
mtnormalise: [ 33%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 6
mtnormalise: [INFO] Balance factors (1):  1.2274110026970686 0.83318153091048341 0.97784575230206605
mtnormalise: [INFO] Balance factors (2):  1.2275214191775266 0.83317488392942185 0.97776559478498271
mtnormalise: [INFO] Balance factors (3):   1.227550932017839 0.83317641501552497 0.97774029054900302
mtnormalise: [INFO] Balance factors (4):   1.227561020616657 0.83318322481113682 0.97772426385153199
mtnormalise: [INFO] Balance factors (5):  1.2275638506982198   0.833185193647487 0.97771969938125103
mtnormalise: [ 40%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 7
mtnormalise: [INFO] Balance factors (1):  1.2275170803292466 0.83384972873070196 0.97697773037626867
mtnormalise: [INFO] Balance factors (2):  1.2276648837485291 0.83382000049678484 0.97689493622386048
mtnormalise: [INFO] Balance factors (3):   1.227695674649093 0.83381333893114717 0.97687823995407308
mtnormalise: [INFO] Balance factors (4):  1.2277047040844111 0.83380988842787163 0.97687509780375903
mtnormalise: [ 47%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 8
mtnormalise: [INFO] Balance factors (1):  1.2276177525737155 0.83425081539557078 0.97642794453073611
mtnormalise: [INFO] Balance factors (2):  1.2276604892783305  0.8342855359023299 0.97635331896761124
mtnormalise: [INFO] Balance factors (3):  1.2276693501121132 0.83429112371855296 0.97633973277758379
mtnormalise: [INFO] Balance factors (4):  1.2276740317737236 0.83429518539948977 0.97633125637212637
mtnormalise: [INFO] Balance factors (5):  1.2276733396442783 0.83429723324950655 0.97632941031625553
mtnormalise: [ 53%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 9
mtnormalise: [INFO] Balance factors (1):  1.2275855060028729 0.83454835163490593 0.97610546480547233
mtnormalise: [INFO] Balance factors (2):  1.2276384873754862 0.83449788593676344 0.97612236565372368
mtnormalise: [INFO] Balance factors (3):  1.2276414002074907  0.8344929649480769 0.97612580576325247
mtnormalise: [INFO] Balance factors (4):  1.2276424422390577 0.83449416945270127 0.97612356828697233
mtnormalise: [ 60%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 10
mtnormalise: [INFO] Balance factors (1):  1.2275618224385711 0.83464064033041929 0.97601636401460767
mtnormalise: [INFO] Balance factors (2):  1.2276127647748549 0.83462160024411325 0.97599812697049126
mtnormalise: [INFO] Balance factors (3):  1.2276148080825051 0.83462140582561239 0.97599672981780294
mtnormalise: [ 67%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 11
mtnormalise: [INFO] Balance factors (1):  1.2275569482966961 0.83471357162057691  0.9759349617344123
mtnormalise: [INFO] Balance factors (2):  1.2275901378579466 0.83470506309014203 0.97591852389394851
mtnormalise: [INFO] Balance factors (3):  1.2275987876749435  0.8347013393026137 0.97591600120604638
mtnormalise: [ 73%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 12
mtnormalise: [INFO] Balance factors (1):  1.2275627799401321  0.8347598892108804 0.97587617488018297
mtnormalise: [INFO] Balance factors (2):  1.2275719785839339 0.83475493249336641 0.97587465692702924
mtnormalise: [INFO] Balance factors (3):  1.2275735918357749 0.83475570065647731 0.97587247642959596
mtnormalise: [ 80%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 13
mtnormalise: [INFO] Balance factors (1):  1.2275500083753368 0.83480853770490371 0.97582945819610845
mtnormalise: [INFO] Balance factors (2):  1.2275739178541132 0.83479773625648446 0.97582307796063128
mtnormalise: [INFO] Balance factors (3):  1.2275762229561777 0.83479733180667826 0.97582171836659082
mtnormalise: [INFO] Balance factors (4):  1.2275819419170997 0.83479377984197334 0.97582132429138246
mtnormalise: [INFO] Balance factors (5):  1.2275762229561777 0.83479733180667826 0.97582171836659082
mtnormalise: [INFO] Balance factors (6):  1.2275819419170997 0.83479377984197334 0.97582132429138246
mtnormalise: [INFO] Balance factors (7):  1.2275762229561777 0.83479733180667826 0.97582171836659082
mtnormalise: [ 87%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 14
mtnormalise: [INFO] Balance factors (1):  1.2275649230210419 0.83483217269045518 0.97578997564851011
mtnormalise: [INFO] Balance factors (2):  1.2275706862581834  0.8348315470461154 0.97578612575809753
mtnormalise: [INFO] Balance factors (3):  1.2275772802640423 0.83482777764171179 0.97578529009789183
mtnormalise: [ 93%] performing log-domain intensity normalisation... 
mtnormalise: [INFO] Iteration: 15
mtnormalise: [INFO] Balance factors (1):  1.2275698467879521 0.83485016159715353 0.97576503604901133
mtnormalise: [INFO] Balance factors (2):   1.227578072278104 0.83485122078134344 0.97575725989861617
mtnormalise: [100%] performing log-domain intensity normalisation
mtnormalise: [INFO] creating image "checkbinnorm.mif"...
mtnormalise: [.   ] writing output images... 
mtnormalise: [INFO] creating image "binnormwm.mif"...
mtnormalise: [  . ] writing output images... 
mtnormalise: [INFO] creating image "binnormgm.mif"...
mtnormalise: [.   ] writing output images... 
mtnormalise: [INFO] creating image "binnormcsf.mif"...
mtnormalise: [done] writing output images

Excellent. It even settles on close to the same balance coefficients that it did for the input data that was not first mangled by mtbin. The stats on this normalisation field (just within the mask):

thijs@bri-pub-59:~/Documents/mtbin_fail$ mrstats checkbinnorm.mif -mask mask.mif 
      volume       mean     median      stdev        min        max      count
       [ 0 ]   0.982526 0.975953877  0.0718988    0.76743    2.59254     792008

Still a significant correction versus the original output of mtbin (stdev is probably the best insight here). This is what it looks like:

That's still quite a significant relative impact in the brain area itself as well, with relative differences between core brain regions that exceed 10% easily. One of the main things I worry about with respect to mtbin still being the current advertised step on master. That and those extreme outliers that mtbin introduces in these ultimate fail cases, as they trip up the (current, on master) population_template script; but I've described that elsewhere and in many, many emails. Using mtnormalise outputs, population_template shouldn't trip up any more (as this was mostly due to the extreme value of these outliers, not so much those extra bits to the mask); but it's still a wise idea regardless to make population_template a bit more robust to outliers. Thanks to @maxpietsch for working on that! (sorry for not commenting on those things yet, but this thing here has sucked up all my time sadly)

I will do those replies to other comments later, as I need to catch a train before they stop running today. In the mean time, just to quickly add the gist on the name thing: I'm very much for mtnormalise, and very much against mtbin2. I'll add detail later (probably tomorrow).

I would personally encourage to adopt the lower threshold for the template fixel masks that I have suggested as a docs change in this pull request, since 0.2 has been found to be consistently inappropriate in a lot of studies now, and this impact is very important on CFE, which affects other fixels even if they are in the mask at the higher threshold.

I don't see this change in the diff...

Also: This is still being applied as a peak amplitude threshold... has an integral threshold been tried instead?

Change is at: https://github.com/MRtrix3/mrtrix3/pull/1036/files#diff-eb4aa7a31236ad6b638bffa6a4485951R129
Has only been changed in the MT pipeline, because a safer threshold may still be needed for ST-CSD. For MT-CSD, my consistent experience until now has been an actual need to set a lower threshold, since otherwise essential fixels are missing from the analysis mask, leading even to loss in power in other fixels. While reproducing a certain study's result, this heavily impacted the final result of the study. It's also quite safe to lower this threshold, since MT-CSD average templates typically look amazingly clean, with almost no to literally completely no spurious peaks in the core bits of the template at all. In a recent study with a template of 40 subjects, we could even lower this threshold to 0.5 and get a pretty awesome looking (and very complete) fixel analysis mask.

Peaks vs integral: in these scenarios would practically make no difference at all (apart from maybe a few random differences at the far extents of the WM (essentially already in the GM), again due to the template being that clean. But now thinking of it, actually "peaks" is also not that bad of a choice: we're essentially after things that are quite certain to be genuine fixels. Peaks have the actual benefit of (especially in a template) being a measure that also factors in certainty, as peaks of subjects that were less certain or less well aligned across the template subjects will be lower in peak amplitude. I know it of course also factors in dispersion, but even that could in a way be seen as less "certainty of fixel". There's definitely arguments against this reasoning too; but so I feel peaks isn't really that bad per se for this particular step. Overall, I wouldn't be too fussed about it personally. Peaks does a pretty good job in dwi2response tournier as well.

Consider instead having mtbin throw an exception immediately, unless an option -override is provided at the command-line.

This has now been implemented, with a lot of warnings here and there to make any user think about it carefully.

Code snippets here and there

All good.

Added options

Outputting the weighted sum image doesn't really help much for better understanding in any scenario I could think of. I don't want to confuse users with this. But users that do want to see this (and then probably understand the algorithm very well, so they can interpret the weighted sum image in context), can simply mrcalc a weighted sum of the outputs of the command, using the weights provided by the -info output. I think that's sufficient for users who have such a keen technical interest in this. The other outputs (inlier / final processing mask and final normalisation field) as it stands are more directly useful and interpretable. The mask is an easy check to spot nothing has gone too wonky, and the normalisation field directly relates input to output in one summary image (and directly shows exactly what correction has actually been done). The lognorm_scale entry that is written, is essentially the geometric average of the final normalisation field in the final processing mask (so without outlier areas). Since we reason (working in the log-domain) that we're looking at log-normal distributions, the geometric average is also known as the/a scale parameter of the log-normal distribution. Note that mtnormalise doesn't try to "factor this out of the bias field" like mtbin used to somewhat do. The normalisation field still directly relates input to output, and users are free to compute any other metric on the normalisation field within the mask to get a sense of what kind of "global normalisation/scaling" has been detected and corrected for in their data. We do write this suggestion as a number to the header, so it's directly accessible over any bigger group of datasets: see the MT FBA pipeline instructions on what to potentially use it for (akin to what mtbin also provided in a way in another kind of parameter that was written to the header).

Allowing changing the outlier fences is also not really needed I reckon; unless a user would at some point come up with a compelling case of why they'd need to, I wouldn't confuse them with it. The numbers chosen here (1.5 and 3.0) are pretty standardised definitions for outliers and extreme outliers (used in e.g SPSS and other packages for boxplots too). They're also non-parametric, so quite safe to use. It doesn't really matter how much the incoming bias fields vary spatially, eventually the algorithm converges and the outliers are only relative to the corrected image at that stage, and end up indicating outliers that could also not be fit by the polynomial.

Finally,

Name

I should rectify this: there has been no intent of abusing the name for IP purposes. In fact, I think it would be unwise to decide command names via IP whatevers: the command names' primary purpose should be to indicate what the command does, rather than err towards an encrypted/obfuscated acronym with version numbers encoded. As @jdtournier mentioned, there wasn't even a published acronym for the method here to begin with; and if I ever write/publish something about this one, I won't be calling it anything "BIN2". If anything, this is a great opportunity to choose a more clear name. The reason why I was to insisting on not going under the mtbin label, is because the last thing I personally want to be known for, is introducing extra confusion. After all the confusion and problems with mtbin, I'm only left dreaming of a future where things are offered in a simple, clear and robust way. mtbin2 would only add confusion to a name that a lot of people already weren't grasping to begin with. Just my 2 cents there. IP-wise: I have no problems acknowledging the inspiration here; I've naturally added Dave to the command author list too. I don't really care actually who gets added to that: anyone's welcome as far as I'm concerned. As long as the thing works and is clear to users. The old mtbin was actually not even referring to the abstract. So even as to documenting what the thing actually did, I don't really think the change to the new command is a regression of anything on any front...

Only annoying thing about the name change right now, is that mtnormalise used to exist as a command, and has a similar user interface for the basic use case (i.e. not asking for any of the specialised outputs). But this is only a temporary annoyance now. Better to make the change sooner than later, and at least we'll end up with one single, simple and clear name that leads users towards this functionality for the purpose that it is intended for.

Not wanting to add a nasty bit of my own, but just being very objective here by the way: tcksift and tcksift2 are in fact the outliers here (so not really the model, I reckon). There are literally no other commands that have an acronym of a method in their name, only of data types (unless I'm overlooking some; apologies if that's the case). And definitely no other ones that use the command name for version numbering. Even more, one could reason: the "2" is used in many command names for an entirely different purpose, that is standardised across command names. So inspired by consistency, we could be discussing these two command names. But please, for heaven's and our time's sake, let's not go there and use the practical (rather than theoretical) bits of our minds here: I'm pretty sure most people already know by now what "sift" stands for, so I'm not asking here to still change these to e.g. tckfilter or tck2fodfit or something. But if we would be travelling several years back in time, I don't think that would've been such crazy proposals to be fair. Even nowadays, one could tcksift2 tracks with an ODF that would e.g. come from NODDI or something (not encouraging such practises though; although I know at least one person who has tckgened on such a NODDI ODF), so that's not per se so "SD informed" if one would really be picky.

Well, just my 2 cents. But for all intents and purposes, I think the command names should just be as clear and simple as possible. "normalise" still only has 9 characters, and "mt" is rather short too. There's a lot of much longer MRtrix command names out there.

Alright, that's it then, finally. I will merge into the tag branch, so don't panic please. 😅 If there's still discussion to be done, there's plenty of opportunities still. ☺️

Not wanting to add a nasty bit of my own, but just being very objective here by the way

But please, for heaven's and our time's sake, let's not go there

🙄

let's not go there

As in: actually considering changing those names (not bringing up the former facts). I don't think these kinds of "in principle" arguments should drive name changes. For me, a command name is part of the user interface, which should be designed just for the purpose of interfacing with the user. Communicating history can be done in manuals, appendices, ... but should not get in the way of the primary goal of the software: to be used.

OK, this branch is already merged with tag_3.0_RC2, but there is one pending issue that hasn't been settled yet, which is the name of the command. I think there's a general consensus towards mtnormalise -
it's what's been merged at the moment anyway. Is everyone is happy with this? If so, this should tie off the discussion around these changes, we can work on drafting a release message for the new tag.

If so, this should tie off the discussion around these changes, we can work on drafting a release message for the new tag.

👍