How to handle division by zero (or epsilon)

[jcohenadad]

hey team! what should be the strategy for handling division by zero (or epsilon) when computing MTR (this applies to other methods too, e.g. MTsat)?

• detect small divisor, mask out those elements during division and replace those elements in the output array by X
• divide everything and handle the output (look for aberrant values and replace by X)
  • one problem with that is that when saving niftis with a very large range (i.e. including e.g. -/+ 10^15), the precision gets hurt in an uncontrolled manner (requires to take care of the output dtype, pushing physical space (e.g. with float64) unnecessarily).
  • another problem (minor), is that division by zero might be handled differently across platforms, potentially yielding differences in precision (see my comment below) and/or warning messages.

What should be X: zeros? nan? inf?

• i see two problems with replacing by zero:
  • users might blindly use the output image to aggregate metrics, and include zeros in their calculations.
  • a mask could be provided alongside the output image, but some uninformed users will not consider the mask and go on with their processing. I've seen that happening.
  • some users will want to resample their output image to another space, yielding critical bias when interpolating zeroed-values (at least, nan/inf would output an error, depending on the software).

I've been trying to handle this at both levels (during the division process, and at the output level), to provide maximum flexibility. See an example here (feedback welcome :-)).

[mathieuboudreau]

Slightly related, but during my PhD I always just handled the outputs (second option) for all methods where voxels were independently processed (except for the case of the divisor being equal exactly to 0, obviously). This let me handle functionality “spurious values” for all my methods the same way i.e. I could reuse masking functions (filter out values beyond a range), and then save masks of those voxels. I would set the values to 0 (as I would with any other mask), but provide the accompanying "outlier" mask & values map for further quality assurance if desired.

That was my logic behind it, anyways.

[mathieuboudreau]

Doing it after also lets you handle your outlier mask with more representative values (i.e. NaN for true divisions by 0 vs the true values of small divisors.), things like that. But this was just my own preference, there might be a better workflow that someone else is using.

[mathieuboudreau]

@spinicist, care to chime in on what you do in QUIT?

[spinicist]

Sadly, I currently do very little in this situation.

The models which require non-linear fitting have the ability to detect bad input data or fitting failures and inform the user (a message is printed for each voxel that fails). Most of the time the problem is a bad mask that includes background.

I think the only straight division in QUIT is in the MTSat program. This does absolutely nothing about bad values. Because MTSat is a simple calculation and not a non-linear fit, it uses the ITK ParallelizeImageRegion function and so it would be quite tedious to add per-voxel error reporting (I'd need mutexes or similar to avoid race conditions in the print statements). I suppose a simple solution could be a flag that I set on encountering bad values and then I print a single message at the end (I quite like this idea). It would also be fairly straightforward to make it detect NaN or Inf and clamp those a particular value, but currently I don't do that.

I'm not sure there is a particularly good way to solve this beyond reporting suspicious values during processing. In my experience a division by zero (or small value) is either a bug or the result of a bad mask that included background regions - hence if the program reports a bad value people should check their masks and reprocess?