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""" Read / write access to SPM2 version of analyze image format """
10import numpy as np
11
12from . import spm99analyze as spm99  # module import
13
14image_dimension_dtd = spm99.image_dimension_dtd[:]
15image_dimension_dtd[
16    image_dimension_dtd.index(('funused2', 'f4'))
17] = ('scl_inter', 'f4')
18
19# Full header numpy dtype combined across sub-fields
20header_dtype = np.dtype(spm99.header_key_dtd +
21                        image_dimension_dtd +
22                        spm99.data_history_dtd)
23
24
25class 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    template_dtype = header_dtype
37
38    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        slope = float(self['scl_slope'])
101        if np.isfinite(slope) and slope:
102            # try to get offset from scl_inter
103            inter = float(self['scl_inter'])
104            if not np.isfinite(inter):
105                inter = 0.0
106            return slope, inter
107        # no non-zero and finite scaling, try gl/cal fields
108        unscaled_range = self['glmax'] - self['glmin']
109        scaled_range = self['cal_max'] - self['cal_min']
110        if unscaled_range and scaled_range:
111            slope = float(scaled_range) / unscaled_range
112            inter = self['cal_min'] - slope * self['glmin']
113            return slope, inter
114        return None, None
115
116    @classmethod
117    def may_contain_header(klass, binaryblock):
118        if len(binaryblock) < klass.sizeof_hdr:
119            return False
120
121        hdr_struct = np.ndarray(shape=(), dtype=header_dtype,
122                                buffer=binaryblock[:klass.sizeof_hdr])
123        bs_hdr_struct = hdr_struct.byteswap()
124        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
128class Spm2AnalyzeImage(spm99.Spm99AnalyzeImage):
129    """ Class for SPM2 variant of basic Analyze image
130    """
131    header_class = Spm2AnalyzeHeader
132
133
134load = Spm2AnalyzeImage.load
135save = Spm2AnalyzeImage.instance_to_filename
136