# 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`

.