nipy/nibabel — pull request #222 by matthew-brett (2014-01-29 23:33:08+00:00)
Title: BF: remove intercept only iuint-iuint scaling

Short version: remove conditional branch trying to preserve precision on
writing arrays by keeping working type as an integer type.

The advantages to removal are:

* fixing a bug for explicit scaling in `array_to_file`
* making the code simpler to think about

The disadvantages are:

* minor loss in precision for int64 / int64 output types.
* float output type rather than integer output type when reading images
where this scaling trick could be used.

Long version:

``int`` means signed integer type
``uint`` means unsigned integer type
``iuint`` refers to set of all signed and unsigned integer types

When writing data to file, in `volumeutils.array_to_file`.

If:

* input type is iuint
* output type is iuint
* slope is 1
* integer is not 0 and is an integer

Then:

* Working type calculated by `volumeutils._inter_type`
* intercept cast to working type

Here is the relevant code from 1.3.0::

# (u)int to (u)int with inter alone - select precision
if (slope == 1 and inter != 0 and
in_dtype.kind in 'iu' and out_dtype.kind in 'iu' and
inter == np.round(inter)): # (u)int to (u)int offset only scaling
# Does range of in type minus inter fit in out type? If so, use that as
# working type. Otherwise use biggest float for max integer precision
inter = inter.astype(_inter_type(in_dtype, -inter, out_dtype))

The problem arises from the sign flip in the line above. Let's call the
internal intercept value in `_inter_type` : `minter` (minus inter in this
case). Define:

* `in_min = np.iinfo(in_dtype).min`
* `in_max = np.iinfo(in_dtype).max`

Looking at `_inter_type`, we will get a returned uint type if all of `in_min +
minter`, `in_max + minter`, `in_min`, `in_max`, `minter` are > 0. If
`_inter_type` returns uint as working type and `inter` is negative, then the
last line above will cast a negative number into a positive one and cause a
mess.

The problem then arises only when:

* input type is uint (in_min = 0)
* `inter` is negative
* output_type is larger than input type (`in_max + minter` > `in_max`)

See appended script for demonstration.

This is clearly a problem for `array_to_file` - but is it a problem for the
image reading / writing prior to the recent pull request allowing arbitrary
scaling?

Then the question is: is there any calculated scaling that results in:

* uint input type
* negative intercept
* larger iuint output type?

But - if the input type is uint and smaller then the output type, then it can
and will cast to the output type, so calculating scaling between these types
will never generate a negative intercept. See the assert in appended script
for demonstration.

The problem doesn't arise for ``volumeutils.apply_read_scaling``, even though
we are using `_inter_type`. This is because we pass `inter` to `_inter_type`
rather than `-inter`. Because we are checking `minter` in `_inter_type`, and
this will be negative, we can't get a uint return type for a negative
intercept.