1
# emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*-
2
# vi: set ft=python sts=4 ts=4 sw=4 et:
3
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
4
#
5
#   See COPYING file distributed along with the NiBabel package for the
6
#   copyright and license terms.
7
#
8
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
9 33
""" Read / write access to SPM2 version of analyze image format """
10 33
import numpy as np
11

12 33
from . import spm99analyze as spm99  # module import
13

14 33
image_dimension_dtd = spm99.image_dimension_dtd[:]
15 33
image_dimension_dtd[
16
    image_dimension_dtd.index(('funused2', 'f4'))
17
] = ('scl_inter', 'f4')
18

19
# Full header numpy dtype combined across sub-fields
20 33
header_dtype = np.dtype(spm99.header_key_dtd +
21
                        image_dimension_dtd +
22
                        spm99.data_history_dtd)
23

24

25 33
class Spm2AnalyzeHeader(spm99.Spm99AnalyzeHeader):
26
    """ Class for SPM2 variant of basic Analyze header
27

28
    SPM2 variant adds the following to basic Analyze format:
29

30
    * voxel origin;
31
    * slope scaling of data;
32
    * reading - but not writing - intercept of data.
33
    """
34

35
    # Copies of module level definitions
36 33
    template_dtype = header_dtype
37

38 33
    def get_slope_inter(self):
39
        """ Get data scaling (slope) and intercept from header data
40

41
        Uses the algorithm from SPM2 spm_vol_ana.m by John Ashburner
42

43
        Parameters
44
        ----------
45
        self : header
46
           Mapping with fields:
47
           * scl_slope - slope
48
           * scl_inter - possible intercept (SPM2 use - shared by nifti)
49
           * glmax - the (recorded) maximum value in the data (unscaled)
50
           * glmin - recorded minimum unscaled value
51
           * cal_max - the calibrated (scaled) maximum value in the dataset
52
           * cal_min - ditto minimum value
53

54
        Returns
55
        -------
56
        scl_slope : None or float
57
            slope.  None if there is no valid scaling from these fields
58
        scl_inter : None or float
59
            intercept.  Also None if there is no valid slope, intercept
60

61
        Examples
62
        --------
63
        >>> fields = {'scl_slope': 1, 'scl_inter': 0, 'glmax': 0, 'glmin': 0,
64
        ...           'cal_max': 0, 'cal_min': 0}
65
        >>> hdr = Spm2AnalyzeHeader()
66
        >>> for key, value in fields.items():
67
        ...     hdr[key] = value
68
        >>> hdr.get_slope_inter()
69
        (1.0, 0.0)
70
        >>> hdr['scl_inter'] = 0.5
71
        >>> hdr.get_slope_inter()
72
        (1.0, 0.5)
73
        >>> hdr['scl_inter'] = np.nan
74
        >>> hdr.get_slope_inter()
75
        (1.0, 0.0)
76

77
        If 'scl_slope' is 0, nan or inf, cannot use 'scl_slope'.
78
        Without valid information in the gl / cal fields, we cannot get
79
        scaling, and return None
80

81
        >>> hdr['scl_slope'] = 0
82
        >>> hdr.get_slope_inter()
83
        (None, None)
84
        >>> hdr['scl_slope'] = np.nan
85
        >>> hdr.get_slope_inter()
86
        (None, None)
87

88
        Valid information in the gl AND cal fields are needed
89

90
        >>> hdr['cal_max'] = 0.8
91
        >>> hdr['cal_min'] = 0.2
92
        >>> hdr.get_slope_inter()
93
        (None, None)
94
        >>> hdr['glmax'] = 110
95
        >>> hdr['glmin'] = 10
96
        >>> np.allclose(hdr.get_slope_inter(), [0.6/100, 0.2-0.6/100*10])
97
        True
98
        """
99
        # get scaling factor from 'scl_slope' (funused1)
100 33
        slope = float(self['scl_slope'])
101 33
        if np.isfinite(slope) and slope:
102
            # try to get offset from scl_inter
103 33
            inter = float(self['scl_inter'])
104 33
            if not np.isfinite(inter):
105 33
                inter = 0.0
106 33
            return slope, inter
107
        # no non-zero and finite scaling, try gl/cal fields
108 33
        unscaled_range = self['glmax'] - self['glmin']
109 33
        scaled_range = self['cal_max'] - self['cal_min']
110 33
        if unscaled_range and scaled_range:
111 33
            slope = float(scaled_range) / unscaled_range
112 33
            inter = self['cal_min'] - slope * self['glmin']
113 33
            return slope, inter
114 33
        return None, None
115

116 33
    @classmethod
117 5
    def may_contain_header(klass, binaryblock):
118 33
        if len(binaryblock) < klass.sizeof_hdr:
119 0
            return False
120

121 33
        hdr_struct = np.ndarray(shape=(), dtype=header_dtype,
122
                                buffer=binaryblock[:klass.sizeof_hdr])
123 33
        bs_hdr_struct = hdr_struct.byteswap()
124 33
        return (binaryblock[344:348] not in (b'ni1\x00', b'n+1\x00') and
125
                348 in (hdr_struct['sizeof_hdr'], bs_hdr_struct['sizeof_hdr']))
126

127

128 33
class Spm2AnalyzeImage(spm99.Spm99AnalyzeImage):
129
    """ Class for SPM2 variant of basic Analyze image
130
    """
131 33
    header_class = Spm2AnalyzeHeader
132

133

134 33
load = Spm2AnalyzeImage.load
135 33
save = Spm2AnalyzeImage.instance_to_filename

Read our documentation on viewing source code .

Loading