###########################
Scalefactors and intercepts
###########################

SPM Analyze and nifti1 images have *scalefactors*.  nifti1 images also have
*intercepts*.  If ``A`` is an array in memory, and ``S`` is the array that will
be written to disk, then::

    R = (A - intercept) / scalefactor

and ``R == S`` if ``R`` is already the data dtype we need to write.

If we load the image from disk, we exactly recover ``S`` (and ``R``).  To get
something approximating ``A`` (say ``Aprime``) we apply the intercept and
scalefactor::

    Aprime = (S * scalefactor) + intercept

In a perfect world ``A`` would be exactly the same as ``Aprime``.  However
``scalefactor`` and ``intercept`` are floating point values.  With floating
point, if ``r = (a - b) / c; p = (r * c) + b`` it is not necessarily true that
``p == a``. For example:

>>> import numpy as np
>>> a = 10
>>> b = np.e
>>> c = np.pi
>>> r = (a - b) / c
>>> p = (r * c) + b
>>> p == a
False

So there will be some error in this reconstruction, even when ``R`` is the same
type as ``S``.

More common is the situation where ``R`` is a different type from ``S``.  If
``R`` is of type ``r_dtype``, ``S`` is of type ``s_dtype`` and
``cast_function(R, dtype)`` is some function that casts ``R`` to the desired
type ``dtype``, then::

    R = (A - intercept) / scalefactor
    S = cast_function(R, s_dtype)
    R_prime = cast_function(S, r_dtype)
    A_prime = (R_prime * scalefactor) + intercept

The type of ``R`` will depend on what numpy did for upcasting ``A, intercept,
scalefactor``.

In order that ``cast_function(S, r_dtype)`` can best reverse ``cast_function(R,
s_dtype)``, the second needs to know the type of ``R``, which is not stored. The
type of ``R`` depends on the types of ``A`` and of ``intercept, scalefactor``.
We don't know the type of ``A`` because it is not stored.

``R`` is likely to be a floating point type because of the application of
scalefactor and intercept. If ``(intercept, scalefactor)`` are not the identity
(0, 1), then we can ensure that ``R`` is at minimum the type of the ``intercept,
scalefactor`` by making these be at least 1D arrays, so that floating point
types will upcast in ``R = (A - intercept) / scalefactor``.

The cast of ``R`` to ``S`` and back to ``R_prime`` can lose resolution if the
types of ``R`` and ``S`` have different resolution.

Our job is to select:

* scalefactor
* intercept
* ``cast_function``

such that we minimize some measure of difference between ``A`` and
``A_prime``.