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 ECAT format images
10

11
An ECAT format image consists of:
12

13
* a *main header*;
14
* at least one *matrix list* (mlist);
15

16
ECAT thinks of memory locations in terms of *blocks*.  One block is 512
17
bytes.  Thus block 1 starts at 0 bytes, block 2 at 512 bytes, and so on.
18

19
The matrix list is an array with one row per frame in the data.
20

21
Columns in the matrix list are:
22

23
* 0: Matrix identifier (frame number)
24
* 1: matrix data start block number (subheader followed by image data)
25
* 2: Last block number of matrix (image) data
26
* 3: Matrix status
27

28
    * 1: hxists - rw
29
    * 2: exists - ro
30
    * 3: matrix deleted
31

32
There is one sub-header for each image frame (or matrix in the terminology
33
above).  A sub-header can also be called an *image header*.  The sub-header is
34
one block (512 bytes), and the frame (image) data follows.
35

36
There is very little documentation of the ECAT format, and many of the comments
37
in this code come from a combination of trial and error and wild speculation.
38

39
XMedcon can read and write ECAT 6 format, and read ECAT 7 format: see
40
http://xmedcon.sourceforge.net and the ECAT files in the source of XMedCon,
41
currently ``libs/tpc/*ecat*`` and ``source/m-ecat*``.  Unfortunately XMedCon is
42
GPL and some of the header files are adapted from CTI files (called CTI code
43
below).  It's not clear what the licenses are for these files.
44
"""
45

46 33
import warnings
47 33
from numbers import Integral
48

49 33
import numpy as np
50

51 33
from .volumeutils import (native_code, swapped_code, make_dt_codes,
52
                          array_from_file)
53 33
from .spatialimages import SpatialImage
54 33
from .arraywriters import make_array_writer
55 33
from .wrapstruct import WrapStruct
56 33
from .fileslice import canonical_slicers, predict_shape, slice2outax
57 33
from .deprecated import deprecate_with_version
58

59 33
BLOCK_SIZE = 512
60

61 33
main_header_dtd = [
62
    ('magic_number', '14S'),
63
    ('original_filename', '32S'),
64
    ('sw_version', np.uint16),
65
    ('system_type', np.uint16),
66
    ('file_type', np.uint16),
67
    ('serial_number', '10S'),
68
    ('scan_start_time', np.uint32),
69
    ('isotope_name', '8S'),
70
    ('isotope_halflife', np.float32),
71
    ('radiopharmaceutical', '32S'),
72
    ('gantry_tilt', np.float32),
73
    ('gantry_rotation', np.float32),
74
    ('bed_elevation', np.float32),
75
    ('intrinsic_tilt', np.float32),
76
    ('wobble_speed', np.uint16),
77
    ('transm_source_type', np.uint16),
78
    ('distance_scanned', np.float32),
79
    ('transaxial_fov', np.float32),
80
    ('angular_compression', np.uint16),
81
    ('coin_samp_mode', np.uint16),
82
    ('axial_samp_mode', np.uint16),
83
    ('ecat_calibration_factor', np.float32),
84
    ('calibration_unitS', np.uint16),
85
    ('calibration_units_type', np.uint16),
86
    ('compression_code', np.uint16),
87
    ('study_type', '12S'),
88
    ('patient_id', '16S'),
89
    ('patient_name', '32S'),
90
    ('patient_sex', '1S'),
91
    ('patient_dexterity', '1S'),
92
    ('patient_age', np.float32),
93
    ('patient_height', np.float32),
94
    ('patient_weight', np.float32),
95
    ('patient_birth_date', np.uint32),
96
    ('physician_name', '32S'),
97
    ('operator_name', '32S'),
98
    ('study_description', '32S'),
99
    ('acquisition_type', np.uint16),
100
    ('patient_orientation', np.uint16),
101
    ('facility_name', '20S'),
102
    ('num_planes', np.uint16),
103
    ('num_frames', np.uint16),
104
    ('num_gates', np.uint16),
105
    ('num_bed_pos', np.uint16),
106
    ('init_bed_position', np.float32),
107
    ('bed_position', '15f'),
108
    ('plane_separation', np.float32),
109
    ('lwr_sctr_thres', np.uint16),
110
    ('lwr_true_thres', np.uint16),
111
    ('upr_true_thres', np.uint16),
112
    ('user_process_code', '10S'),
113
    ('acquisition_mode', np.uint16),
114
    ('bin_size', np.float32),
115
    ('branching_fraction', np.float32),
116
    ('dose_start_time', np.uint32),
117
    ('dosage', np.float32),
118
    ('well_counter_corr_factor', np.float32),
119
    ('data_units', '32S'),
120
    ('septa_state', np.uint16),
121
    ('fill', '12S')
122
]
123 33
hdr_dtype = np.dtype(main_header_dtd)
124

125

126 33
subheader_dtd = [
127
    ('data_type', np.uint16),
128
    ('num_dimensions', np.uint16),
129
    ('x_dimension', np.uint16),
130
    ('y_dimension', np.uint16),
131
    ('z_dimension', np.uint16),
132
    ('x_offset', np.float32),
133
    ('y_offset', np.float32),
134
    ('z_offset', np.float32),
135
    ('recon_zoom', np.float32),
136
    ('scale_factor', np.float32),
137
    ('image_min', np.int16),
138
    ('image_max', np.int16),
139
    ('x_pixel_size', np.float32),
140
    ('y_pixel_size', np.float32),
141
    ('z_pixel_size', np.float32),
142
    ('frame_duration', np.uint32),
143
    ('frame_start_time', np.uint32),
144
    ('filter_code', np.uint16),
145
    ('x_resolution', np.float32),
146
    ('y_resolution', np.float32),
147
    ('z_resolution', np.float32),
148
    ('num_r_elements', np.float32),
149
    ('num_angles', np.float32),
150
    ('z_rotation_angle', np.float32),
151
    ('decay_corr_fctr', np.float32),
152
    ('corrections_applied', np.uint32),
153
    ('gate_duration', np.uint32),
154
    ('r_wave_offset', np.uint32),
155
    ('num_accepted_beats', np.uint32),
156
    ('filter_cutoff_frequency', np.float32),
157
    ('filter_resolution', np.float32),
158
    ('filter_ramp_slope', np.float32),
159
    ('filter_order', np.uint16),
160
    ('filter_scatter_fraction', np.float32),
161
    ('filter_scatter_slope', np.float32),
162
    ('annotation', '40S'),
163
    ('mt_1_1', np.float32),
164
    ('mt_1_2', np.float32),
165
    ('mt_1_3', np.float32),
166
    ('mt_2_1', np.float32),
167
    ('mt_2_2', np.float32),
168
    ('mt_2_3', np.float32),
169
    ('mt_3_1', np.float32),
170
    ('mt_3_2', np.float32),
171
    ('mt_3_3', np.float32),
172
    ('rfilter_cutoff', np.float32),
173
    ('rfilter_resolution', np.float32),
174
    ('rfilter_code', np.uint16),
175
    ('rfilter_order', np.uint16),
176
    ('zfilter_cutoff', np.float32),
177
    ('zfilter_resolution', np.float32),
178
    ('zfilter_code', np.uint16),
179
    ('zfilter_order', np.uint16),
180
    ('mt_4_1', np.float32),
181
    ('mt_4_2', np.float32),
182
    ('mt_4_3', np.float32),
183
    ('scatter_type', np.uint16),
184
    ('recon_type', np.uint16),
185
    ('recon_views', np.uint16),
186
    ('fill', '174S'),
187
    ('fill2', '96S')]
188 33
subhdr_dtype = np.dtype(subheader_dtd)
189

190
# Ecat Data Types
191
# See:
192
# http://www.turkupetcentre.net/software/libdoc/libtpcimgio/ecat7_8h_source.html#l00060
193
# and:
194
# http://www.turkupetcentre.net/software/libdoc/libtpcimgio/ecat7r_8c_source.html#l00717
195 33
_dtdefs = (  # code, name, equivalent dtype
196
    (1, 'ECAT7_BYTE', np.uint8),
197
    # Byte signed? https://github.com/nipy/nibabel/pull/302/files#r28275780
198
    (2, 'ECAT7_VAXI2', np.int16),
199
    (3, 'ECAT7_VAXI4', np.int32),
200
    (4, 'ECAT7_VAXR4', np.float32),
201
    (5, 'ECAT7_IEEER4', np.float32),
202
    (6, 'ECAT7_SUNI2', np.int16),
203
    (7, 'ECAT7_SUNI4', np.int32))
204 33
data_type_codes = make_dt_codes(_dtdefs)
205

206

207
# Matrix File Types
208 33
ft_defs = (  # code, name
209
    (0, 'ECAT7_UNKNOWN'),
210
    (1, 'ECAT7_2DSCAN'),
211
    (2, 'ECAT7_IMAGE16'),
212
    (3, 'ECAT7_ATTEN'),
213
    (4, 'ECAT7_2DNORM'),
214
    (5, 'ECAT7_POLARMAP'),
215
    (6, 'ECAT7_VOLUME8'),
216
    (7, 'ECAT7_VOLUME16'),
217
    (8, 'ECAT7_PROJ'),
218
    (9, 'ECAT7_PROJ16'),
219
    (10, 'ECAT7_IMAGE8'),
220
    (11, 'ECAT7_3DSCAN'),
221
    (12, 'ECAT7_3DSCAN8'),
222
    (13, 'ECAT7_3DNORM'),
223
    (14, 'ECAT7_3DSCANFIT'))
224 33
file_type_codes = dict(ft_defs)
225

226 33
patient_orient_defs = (  # code, description
227
    (0, 'ECAT7_Feet_First_Prone'),
228
    (1, 'ECAT7_Head_First_Prone'),
229
    (2, 'ECAT7_Feet_First_Supine'),
230
    (3, 'ECAT7_Head_First_Supine'),
231
    (4, 'ECAT7_Feet_First_Decubitus_Right'),
232
    (5, 'ECAT7_Head_First_Decubitus_Right'),
233
    (6, 'ECAT7_Feet_First_Decubitus_Left'),
234
    (7, 'ECAT7_Head_First_Decubitus_Left'),
235
    (8, 'ECAT7_Unknown_Orientation'))
236 33
patient_orient_codes = dict(patient_orient_defs)
237

238
# Indexes from the patient_orient_defs structure defined above for the
239
# neurological and radiological viewing conventions
240 33
patient_orient_radiological = [0, 2, 4, 6]
241 33
patient_orient_neurological = [1, 3, 5, 7]
242

243

244 33
class EcatHeader(WrapStruct):
245
    """Class for basic Ecat PET header
246

247
    Sub-parts of standard Ecat File
248

249
    * main header
250
    * matrix list
251
      which lists the information for each frame collected (can have 1 to many
252
      frames)
253
    * subheaders specific to each frame with possibly-variable sized data
254
      blocks
255

256
    This just reads the main Ecat Header, it does not load the data or read the
257
    mlist or any sub headers
258
    """
259 33
    template_dtype = hdr_dtype
260 33
    _ft_codes = file_type_codes
261 33
    _patient_orient_codes = patient_orient_codes
262

263 33
    def __init__(self,
264
                 binaryblock=None,
265
                 endianness=None,
266
                 check=True):
267
        """Initialize Ecat header from bytes object
268

269
        Parameters
270
        ----------
271
        binaryblock : {None, bytes} optional
272
            binary block to set into header, By default, None in which case we
273
            insert default empty header block
274
        endianness : {None, '<', '>', other endian code}, optional
275
            endian code of binary block, If None, guess endianness
276
            from the data
277
        check : {True, False}, optional
278
            Whether to check and fix header for errors.  No checks currently
279
            implemented, so value has no effect.
280
        """
281 33
        super(EcatHeader, self).__init__(binaryblock, endianness, check)
282

283 33
    @classmethod
284 5
    def guessed_endian(klass, hdr):
285
        """Guess endian from MAGIC NUMBER value of header data
286
        """
287 33
        if not hdr['sw_version'] == 74:
288 33
            return swapped_code
289
        else:
290 33
            return native_code
291

292 33
    @classmethod
293 33
    def default_structarr(klass, endianness=None):
294
        """ Return header data for empty header with given endianness
295
        """
296 33
        hdr_data = super(EcatHeader, klass).default_structarr(endianness)
297 33
        hdr_data['magic_number'] = 'MATRIX72'
298 33
        hdr_data['sw_version'] = 74
299 33
        hdr_data['num_frames'] = 0
300 33
        hdr_data['file_type'] = 0  # Unknown
301 33
        hdr_data['ecat_calibration_factor'] = 1.0  # scale factor
302 33
        return hdr_data
303

304 33
    def get_data_dtype(self):
305
        """ Get numpy dtype for data from header"""
306 33
        raise NotImplementedError("dtype is only valid from subheaders")
307

308 33
    def get_patient_orient(self):
309
        """ gets orientation of patient based on code stored
310
        in header, not always reliable
311
        """
312 33
        code = self._structarr['patient_orientation'].item()
313 33
        if code not in self._patient_orient_codes:
314 0
            raise KeyError('Ecat Orientation CODE %d not recognized' % code)
315 33
        return self._patient_orient_codes[code]
316

317 33
    def get_filetype(self):
318
        """ Type of ECAT Matrix File from code stored in header"""
319 33
        code = self._structarr['file_type'].item()
320 33
        if code not in self._ft_codes:
321 0
            raise KeyError('Ecat Filetype CODE %d not recognized' % code)
322 33
        return self._ft_codes[code]
323

324 33
    @classmethod
325 5
    def _get_checks(klass):
326
        """ Return sequence of check functions for this class """
327 33
        return ()
328

329

330 33
def read_mlist(fileobj, endianness):
331
    """ read (nframes, 4) matrix list array from `fileobj`
332

333
    Parameters
334
    ----------
335
    fileobj : file-like
336
        an open file-like object implementing ``seek`` and ``read``
337

338
    Returns
339
    -------
340
    mlist : (nframes, 4) ndarray
341
        matrix list is an array with ``nframes`` rows and columns:
342

343
        * 0: Matrix identifier (frame number)
344
        * 1: matrix data start block number (subheader followed by image data)
345
        * 2: Last block number of matrix (image) data
346
        * 3: Matrix status
347

348
            * 1: hxists - rw
349
            * 2: exists - ro
350
            * 3: matrix deleted
351

352
    Notes
353
    -----
354
    A block is 512 bytes.
355

356
    ``block_no`` in the code below is 1-based.  block 1 is the main header,
357
    and the mlist blocks start at block number 2.
358

359
    The 512 bytes in an mlist block contain 32 rows of the int32 (nframes,
360
    4) mlist matrix.
361

362
    The first row of these 32 looks like a special row.  The 4 values appear
363
    to be (respectively):
364

365
    * not sure - maybe negative number of mlist rows (out of 31) that are
366
      blank and not used in this block.  Called `nfree` but unused in CTI
367
      code;
368
    * block_no - of next set of mlist entries or 2 if no more entries. We also
369
      allow 1 or 0 to signal no more entries;
370
    * <no idea>.  Called `prvblk` in CTI code, so maybe previous block no;
371
    * n_rows - number of mlist rows in this block (between ?0 and 31) (called
372
      `nused` in CTI code).
373
    """
374 33
    dt = np.dtype(np.int32)
375 33
    if endianness is not native_code:
376 33
        dt = dt.newbyteorder(endianness)
377 33
    mlists = []
378 33
    mlist_index = 0
379 33
    mlist_block_no = 2  # 1-based indexing, block with first mlist
380 28
    while True:
381
        # Read block containing mlist entries
382 33
        fileobj.seek((mlist_block_no - 1) * BLOCK_SIZE)  # fix 1-based indexing
383 33
        dat = fileobj.read(BLOCK_SIZE)
384 33
        rows = np.ndarray(shape=(32, 4), dtype=dt, buffer=dat)
385
        # First row special, points to next mlist entries if present
386 33
        n_unused, mlist_block_no, _, n_rows = rows[0]
387 33
        if not (n_unused + n_rows) == 31:  # Some error condition here?
388 0
            mlist = []
389 0
            return mlist
390
        # Use all but first housekeeping row
391 33
        mlists.append(rows[1:n_rows + 1])
392 33
        mlist_index += n_rows
393 33
        if mlist_block_no <= 2:  # should block_no in (1, 2) be an error?
394 33
            break
395 33
    return np.row_stack(mlists)
396

397

398 33
def get_frame_order(mlist):
399
    """Returns the order of the frames stored in the file
400
    Sometimes Frames are not stored in the file in
401
    chronological order, this can be used to extract frames
402
    in correct order
403

404
    Returns
405
    -------
406
    id_dict: dict mapping frame number -> [mlist_row, mlist_id]
407

408
    (where mlist id is value in the first column of the mlist matrix )
409

410
    Examples
411
    --------
412
    >>> import os
413
    >>> import nibabel as nib
414
    >>> nibabel_dir = os.path.dirname(nib.__file__)
415
    >>> from nibabel import ecat
416
    >>> ecat_file = os.path.join(nibabel_dir,'tests','data','tinypet.v')
417
    >>> img = ecat.load(ecat_file)
418
    >>> mlist = img.get_mlist()
419
    >>> get_frame_order(mlist)
420
    {0: [0, 16842758]}
421
    """
422 33
    ids = mlist[:, 0].copy()
423 33
    n_valid = np.sum(ids > 0)
424 33
    ids[ids <= 0] = ids.max() + 1  # put invalid frames at end after sort
425 33
    valid_order = np.argsort(ids)
426 33
    if not all(valid_order == sorted(valid_order)):
427
        # raise UserWarning if Frames stored out of order
428 33
        warnings.warn_explicit(f'Frames stored out of order; true order = {valid_order}\n'
429
                               'frames will be accessed in order STORED, NOT true order',
430
                               UserWarning, 'ecat', 0)
431 33
    id_dict = {}
432 33
    for i in range(n_valid):
433 33
        id_dict[i] = [valid_order[i], ids[valid_order[i]]]
434 33
    return id_dict
435

436

437 33
def get_series_framenumbers(mlist):
438
    """ Returns framenumber of data as it was collected,
439
    as part of a series; not just the order of how it was
440
    stored in this or across other files
441

442
    For example, if the data is split between multiple files
443
    this should give you the true location of this frame as
444
    collected in the series
445
    (Frames are numbered starting at ONE (1) not Zero)
446

447
    Returns
448
    -------
449
    frame_dict: dict mapping order_stored -> frame in series
450
            where frame in series counts from 1; [1,2,3,4...]
451

452
    Examples
453
    --------
454
    >>> import os
455
    >>> import nibabel as nib
456
    >>> nibabel_dir = os.path.dirname(nib.__file__)
457
    >>> from nibabel import ecat
458
    >>> ecat_file = os.path.join(nibabel_dir,'tests','data','tinypet.v')
459
    >>> img = ecat.load(ecat_file)
460
    >>> mlist = img.get_mlist()
461
    >>> get_series_framenumbers(mlist)
462
    {0: 1}
463
    """
464 33
    nframes = len(mlist)
465 33
    frames_order = get_frame_order(mlist)
466 33
    mlist_nframes = len(frames_order)
467 33
    trueframenumbers = np.arange(nframes - mlist_nframes, nframes)
468 33
    frame_dict = {}
469 33
    for frame_stored, (true_order, _) in frames_order.items():
470
        # frame as stored in file -> true number in series
471 33
        try:
472 33
            frame_dict[frame_stored] = trueframenumbers[true_order] + 1
473 33
        except IndexError:
474 33
            raise IOError('Error in header or mlist order unknown')
475 33
    return frame_dict
476

477

478 33
def read_subheaders(fileobj, mlist, endianness):
479
    """ Retrieve all subheaders and return list of subheader recarrays
480

481
    Parameters
482
    ----------
483
    fileobj : file-like
484
        implementing ``read`` and ``seek``
485
    mlist : (nframes, 4) ndarray
486
        Columns are:
487
        * 0 - Matrix identifier.
488
        * 1 - subheader block number
489
        * 2 - Last block number of matrix data block.
490
        * 3 - Matrix status
491
    endianness : {'<', '>'}
492
        little / big endian code
493

494
    Returns
495
    -------
496
    subheaders : list
497
        List of subheader structured arrays
498
    """
499 33
    subheaders = []
500 33
    dt = subhdr_dtype
501 33
    if endianness is not native_code:
502 33
        dt = dt.newbyteorder(endianness)
503 33
    for mat_id, sh_blkno, sh_last_blkno, mat_stat in mlist:
504 33
        if sh_blkno == 0:
505 0
            break
506 33
        offset = (sh_blkno - 1) * BLOCK_SIZE
507 33
        fileobj.seek(offset)
508 33
        tmpdat = fileobj.read(BLOCK_SIZE)
509 33
        sh = np.ndarray(shape=(), dtype=dt, buffer=tmpdat)
510 33
        subheaders.append(sh)
511 33
    return subheaders
512

513

514 33
class EcatSubHeader(object):
515

516 33
    _subhdrdtype = subhdr_dtype
517 33
    _data_type_codes = data_type_codes
518

519 33
    def __init__(self, hdr, mlist, fileobj):
520
        """parses the subheaders in the ecat (.v) file
521
        there is one subheader for each frame in the ecat file
522

523
        Parameters
524
        -----------
525
        hdr : EcatHeader
526
            ECAT main header
527
        mlist : array shape (N, 4)
528
            Matrix list
529
        fileobj : ECAT file <filename>.v  fileholder or file object
530
                  with read, seek methods
531
        """
532 33
        self._header = hdr
533 33
        self.endianness = hdr.endianness
534 33
        self._mlist = mlist
535 33
        self.fileobj = fileobj
536 33
        self.subheaders = read_subheaders(fileobj, mlist, hdr.endianness)
537

538 33
    def get_shape(self, frame=0):
539
        """ returns shape of given frame"""
540 33
        subhdr = self.subheaders[frame]
541 33
        x = subhdr['x_dimension'].item()
542 33
        y = subhdr['y_dimension'].item()
543 33
        z = subhdr['z_dimension'].item()
544 33
        return x, y, z
545

546 33
    def get_nframes(self):
547
        """returns number of frames"""
548 33
        framed = get_frame_order(self._mlist)
549 33
        return len(framed)
550

551 33
    def _check_affines(self):
552
        """checks if all affines are equal across frames"""
553 33
        nframes = self.get_nframes()
554 33
        if nframes == 1:
555 33
            return True
556 33
        affs = [self.get_frame_affine(i) for i in range(nframes)]
557 33
        if affs:
558 27
            i = iter(affs)
559 27
            first = next(i)
560 33
            for item in i:
561 33
                if not np.allclose(first, item):
562 0
                    return False
563 27
        return True
564

565 33
    def get_frame_affine(self, frame=0):
566
        """returns best affine for given frame of data"""
567 33
        subhdr = self.subheaders[frame]
568 33
        x_off = subhdr['x_offset']
569 33
        y_off = subhdr['y_offset']
570 33
        z_off = subhdr['z_offset']
571

572 33
        zooms = self.get_zooms(frame=frame)
573

574 33
        dims = self.get_shape(frame)
575
        # get translations from center of image
576 33
        origin_offset = (np.array(dims) - 1) / 2.0
577 33
        aff = np.diag(zooms)
578 33
        aff[:3, -1] = -origin_offset * zooms[:-1] + np.array([x_off, y_off,
579
                                                              z_off])
580 33
        return aff
581

582 33
    def get_zooms(self, frame=0):
583
        """returns zooms  ...pixdims"""
584 33
        subhdr = self.subheaders[frame]
585 33
        x_zoom = subhdr['x_pixel_size'] * 10
586 33
        y_zoom = subhdr['y_pixel_size'] * 10
587 33
        z_zoom = subhdr['z_pixel_size'] * 10
588 33
        return (x_zoom, y_zoom, z_zoom, 1)
589

590 33
    def _get_data_dtype(self, frame):
591 33
        dtcode = self.subheaders[frame]['data_type'].item()
592 33
        return self._data_type_codes.dtype[dtcode]
593

594 33
    def _get_frame_offset(self, frame=0):
595 33
        return int(self._mlist[frame][1] * BLOCK_SIZE)
596

597 33
    def _get_oriented_data(self, raw_data, orientation=None):
598
        """
599
        Get data oriented following ``patient_orientation`` header field. If
600
        the ``orientation`` parameter is given, return data according to this
601
        orientation.
602

603
        :param raw_data: Numpy array containing the raw data
604
        :param orientation: None (default), 'neurological' or 'radiological'
605
        :rtype: Numpy array containing the oriented data
606
        """
607 33
        if orientation is None:
608 33
            orientation = self._header['patient_orientation']
609 33
        elif orientation == 'neurological':
610 0
            orientation = patient_orient_neurological[0]
611 33
        elif orientation == 'radiological':
612 0
            orientation = patient_orient_radiological[0]
613
        else:
614 0
            raise ValueError('orientation should be None,\
615
                neurological or radiological')
616

617 33
        if orientation in patient_orient_neurological:
618 0
            raw_data = raw_data[::-1, ::-1, ::-1]
619 33
        elif orientation in patient_orient_radiological:
620 0
            raw_data = raw_data[::, ::-1, ::-1]
621

622 33
        return raw_data
623

624 33
    def raw_data_from_fileobj(self, frame=0, orientation=None):
625
        """
626
        Get raw data from file object.
627

628
        :param frame: Time frame index from where to fetch data
629
        :param orientation: None (default), 'neurological' or 'radiological'
630
        :rtype: Numpy array containing (possibly oriented) raw data
631

632
        .. seealso:: data_from_fileobj
633
        """
634 33
        dtype = self._get_data_dtype(frame)
635 33
        if self._header.endianness is not native_code:
636 33
            dtype = dtype.newbyteorder(self._header.endianness)
637 33
        shape = self.get_shape(frame)
638 33
        offset = self._get_frame_offset(frame)
639 33
        fid_obj = self.fileobj
640 33
        raw_data = array_from_file(shape, dtype, fid_obj, offset=offset)
641 33
        raw_data = self._get_oriented_data(raw_data, orientation)
642 33
        return raw_data
643

644 33
    def data_from_fileobj(self, frame=0, orientation=None):
645
        """
646
        Read scaled data from file for a given frame
647

648
        :param frame: Time frame index from where to fetch data
649
        :param orientation: None (default), 'neurological' or 'radiological'
650
        :rtype: Numpy array containing (possibly oriented) raw data
651

652
        .. seealso:: raw_data_from_fileobj
653
        """
654 33
        header = self._header
655 33
        subhdr = self.subheaders[frame]
656 33
        raw_data = self.raw_data_from_fileobj(frame, orientation)
657
        # Scale factors have to be set to scalars to force scalar upcasting
658 33
        data = raw_data * header['ecat_calibration_factor'].item()
659 33
        data = data * subhdr['scale_factor'].item()
660 33
        return data
661

662

663 33
class EcatImageArrayProxy(object):
664
    """ Ecat implemention of array proxy protocol
665

666
    The array proxy allows us to freeze the passed fileobj and
667
    header such that it returns the expected data array.
668
    """
669

670 33
    def __init__(self, subheader):
671 33
        self._subheader = subheader
672 33
        self._data = None
673 33
        x, y, z = subheader.get_shape()
674 33
        nframes = subheader.get_nframes()
675 33
        self._shape = (x, y, z, nframes)
676

677 33
    @property
678 5
    def shape(self):
679 33
        return self._shape
680

681 33
    @property
682 5
    def ndim(self):
683 33
        return len(self.shape)
684

685 33
    @property
686 5
    def is_proxy(self):
687 33
        return True
688

689 33
    def __array__(self, dtype=None):
690
        """ Read of data from file
691

692
        This reads ALL FRAMES into one array, can be memory expensive.
693

694
        If you want to read only some slices, use the slicing syntax
695
        (``__getitem__``) below, or ``subheader.data_from_fileobj(frame)``
696

697
        Parameters
698
        ----------
699
        dtype : numpy dtype specifier, optional
700
            A numpy dtype specifier specifying the type of the returned array.
701

702
        Returns
703
        -------
704
        array
705
            Scaled image data with type `dtype`.
706
        """
707
        # dtype=None is interpreted as float64
708 33
        data = np.empty(self.shape)
709 33
        frame_mapping = get_frame_order(self._subheader._mlist)
710 33
        for i in sorted(frame_mapping):
711 33
            data[:, :, :, i] = self._subheader.data_from_fileobj(
712
                frame_mapping[i][0])
713 33
        if dtype is not None:
714 33
            data = data.astype(dtype, copy=False)
715 33
        return data
716

717 33
    def __getitem__(self, sliceobj):
718
        """ Return slice `sliceobj` from ECAT data, optimizing if possible
719
        """
720 33
        sliceobj = canonical_slicers(sliceobj, self.shape)
721
        # Indices into sliceobj referring to image axes
722 33
        ax_inds = [i for i, obj in enumerate(sliceobj) if obj is not None]
723 33
        assert len(ax_inds) == len(self.shape)
724 33
        frame_mapping = get_frame_order(self._subheader._mlist)
725
        # Analyze index for 4th axis
726 33
        slice3 = sliceobj[ax_inds[3]]
727
        # We will load volume by volume.  Make slicer into volume by dropping
728
        # index over the volume axis
729 33
        in_slicer = sliceobj[:ax_inds[3]] + sliceobj[ax_inds[3] + 1:]
730
        # int index for 4th axis, load one slice
731 33
        if isinstance(slice3, Integral):
732 33
            data = self._subheader.data_from_fileobj(frame_mapping[slice3][0])
733 33
            return data[in_slicer]
734
        # slice axis for 4th axis, we will iterate over slices
735 33
        out_shape = predict_shape(sliceobj, self.shape)
736 33
        out_data = np.empty(out_shape)
737
        # Slice into output data with out_slicer
738 33
        out_slicer = [slice(None)] * len(out_shape)
739
        # Work out axis corresponding to volume in output
740 33
        in2out_ind = slice2outax(len(self.shape), sliceobj)[3]
741
        # Iterate over specified 4th axis indices
742 33
        for i in list(range(self.shape[3]))[slice3]:
743 33
            data = self._subheader.data_from_fileobj(
744
                frame_mapping[i][0])
745 33
            out_slicer[in2out_ind] = i
746 33
            out_data[tuple(out_slicer)] = data[in_slicer]
747 33
        return out_data
748

749

750 33
class EcatImage(SpatialImage):
751
    """ Class returns a list of Ecat images, with one image(hdr/data) per frame
752
    """
753 33
    _header = EcatHeader
754 33
    header_class = _header
755 33
    valid_exts = ('.v',)
756 33
    _subheader = EcatSubHeader
757 33
    files_types = (('image', '.v'), ('header', '.v'))
758

759 33
    ImageArrayProxy = EcatImageArrayProxy
760

761 33
    def __init__(self, dataobj, affine, header,
762
                 subheader, mlist,
763
                 extra=None, file_map=None):
764
        """ Initialize Image
765

766
        The image is a combination of
767
        (array, affine matrix, header, subheader, mlist)
768
        with optional meta data in `extra`, and filename / file-like objects
769
        contained in the `file_map`.
770

771
        Parameters
772
        ----------
773
        dataobj : array-like
774
            image data
775
        affine : None or (4,4) array-like
776
            homogeneous affine giving relationship between voxel coords and
777
            world coords.
778
        header : None or header instance
779
            meta data for this image format
780
        subheader : None or subheader instance
781
            meta data for each sub-image for frame in the image
782
        mlist : None or array
783
            Matrix list array giving offset and order of data in file
784
        extra : None or mapping, optional
785
            metadata associated with this image that cannot be
786
            stored in header or subheader
787
        file_map : mapping, optional
788
            mapping giving file information for this image format
789

790
        Examples
791
        --------
792
        >>> import os
793
        >>> import nibabel as nib
794
        >>> nibabel_dir = os.path.dirname(nib.__file__)
795
        >>> from nibabel import ecat
796
        >>> ecat_file = os.path.join(nibabel_dir,'tests','data','tinypet.v')
797
        >>> img = ecat.load(ecat_file)
798
        >>> frame0 = img.get_frame(0)
799
        >>> frame0.shape == (10, 10, 3)
800
        True
801
        >>> data4d = img.get_fdata()
802
        >>> data4d.shape == (10, 10, 3, 1)
803
        True
804
        """
805 33
        self._subheader = subheader
806 33
        self._mlist = mlist
807 33
        self._dataobj = dataobj
808 33
        if affine is not None:
809
            # Check that affine is array-like 4,4.  Maybe this is too strict at
810
            # this abstract level, but so far I think all image formats we know
811
            # do need 4,4.
812 33
            affine = np.array(affine, dtype=np.float64, copy=True)
813 33
            if not affine.shape == (4, 4):
814 0
                raise ValueError('Affine should be shape 4,4')
815 33
        self._affine = affine
816 33
        if extra is None:
817 33
            extra = {}
818 33
        self.extra = extra
819 33
        self._header = header
820 33
        if file_map is None:
821 33
            file_map = self.__class__.make_file_map()
822 33
        self.file_map = file_map
823 33
        self._data_cache = None
824 33
        self._fdata_cache = None
825

826 33
    @property
827 5
    def affine(self):
828 33
        if not self._subheader._check_affines():
829 0
            warnings.warn('Affines different across frames, loading affine '
830
                          'from FIRST frame', UserWarning)
831 33
        return self._affine
832

833 33
    def get_frame_affine(self, frame):
834
        """returns 4X4 affine"""
835 0
        return self._subheader.get_frame_affine(frame=frame)
836

837 33
    def get_frame(self, frame, orientation=None):
838
        """
839
        Get full volume for a time frame
840

841
        :param frame: Time frame index from where to fetch data
842
        :param orientation: None (default), 'neurological' or 'radiological'
843
        :rtype: Numpy array containing (possibly oriented) raw data
844
        """
845 33
        return self._subheader.data_from_fileobj(frame, orientation)
846

847 33
    def get_data_dtype(self, frame):
848 27
        subhdr = self._subheader
849 27
        dt = subhdr._get_data_dtype(frame)
850 27
        return dt
851

852 33
    @property
853 5
    def shape(self):
854 33
        x, y, z = self._subheader.get_shape()
855 33
        nframes = self._subheader.get_nframes()
856 33
        return(x, y, z, nframes)
857

858 33
    def get_mlist(self):
859
        """ get access to the mlist
860
        """
861 33
        return self._mlist
862

863 33
    def get_subheaders(self):
864
        """get access to subheaders"""
865 33
        return self._subheader
866

867 33
    @classmethod
868 33
    @deprecate_with_version('from_filespec class method is deprecated.\n'
869
                            'Please use the ``from_file_map`` class method '
870
                            'instead.',
871
                            '2.1', '4.0')
872 5
    def from_filespec(klass, filespec):
873 33
        return klass.from_filename(filespec)
874

875 33
    @staticmethod
876 5
    def _get_fileholders(file_map):
877
        """ returns files specific to header and image of the image
878
        for ecat .v this is the same image file
879

880
        Returns
881
        -------
882
        header : file holding header data
883
        image : file holding image data
884
        """
885 33
        return file_map['header'], file_map['image']
886

887 33
    @classmethod
888 33
    def from_file_map(klass, file_map, *, mmap=True, keep_file_open=None):
889
        """class method to create image from mapping
890
        specified in file_map
891
        """
892 33
        hdr_file, img_file = klass._get_fileholders(file_map)
893
        # note header and image are in same file
894 33
        hdr_fid = hdr_file.get_prepare_fileobj(mode='rb')
895 33
        header = klass._header.from_fileobj(hdr_fid)
896 33
        hdr_copy = header.copy()
897
        # LOAD MLIST
898 33
        mlist = np.zeros((header['num_frames'], 4), dtype=np.int32)
899 33
        mlist_data = read_mlist(hdr_fid, hdr_copy.endianness)
900 33
        mlist[:len(mlist_data)] = mlist_data
901
        # LOAD SUBHEADERS
902 33
        subheaders = klass._subheader(hdr_copy, mlist, hdr_fid)
903
        # LOAD DATA
904
        # Class level ImageArrayProxy
905 33
        data = klass.ImageArrayProxy(subheaders)
906
        # Get affine
907 33
        if not subheaders._check_affines():
908 0
            warnings.warn('Affines different across frames, loading affine '
909
                          'from FIRST frame', UserWarning)
910 33
        aff = subheaders.get_frame_affine()
911 33
        img = klass(data, aff, header, subheaders, mlist,
912
                    extra=None, file_map=file_map)
913 33
        return img
914

915 33
    def _get_empty_dir(self):
916
        """
917
        Get empty directory entry of the form
918
        [numAvail, nextDir, previousDir, numUsed]
919
        """
920 33
        return np.array([31, 2, 0, 0], dtype=np.int32)
921

922 33
    def _write_data(self, data, stream, pos, dtype=None, endianness=None):
923
        """
924
        Write data to ``stream`` using an array_writer
925

926
        :param data: Numpy array containing the dat
927
        :param stream: The file-like object to write the data to
928
        :param pos: The position in the stream to write the data to
929
        :param endianness: Endianness code of the data to write
930
        """
931 33
        if dtype is None:
932 33
            dtype = data.dtype
933

934 33
        if endianness is None:
935 33
            endianness = native_code
936

937 33
        stream.seek(pos)
938 33
        make_array_writer(data.newbyteorder(endianness),
939
                          dtype).to_fileobj(stream)
940

941 33
    def to_file_map(self, file_map=None):
942
        """ Write ECAT7 image to `file_map` or contained ``self.file_map``
943

944
        The format consist of:
945

946
        - A main header (512L) with dictionary entries in the form
947
            [numAvail, nextDir, previousDir, numUsed]
948
        - For every frame (3D volume in 4D data)
949
          - A subheader (size = frame_offset)
950
          - Frame data (3D volume)
951
        """
952 33
        if file_map is None:
953 33
            file_map = self.file_map
954

955
        # It appears to be necessary to load the data before saving even if the
956
        # data itself is not used.
957 33
        self.get_fdata()
958 33
        hdr = self.header
959 33
        mlist = self._mlist
960 33
        subheaders = self.get_subheaders()
961 33
        dir_pos = 512
962 33
        entry_pos = dir_pos + 16  # 528
963 33
        current_dir = self._get_empty_dir()
964

965 33
        hdr_fh, img_fh = self._get_fileholders(file_map)
966 33
        hdrf = hdr_fh.get_prepare_fileobj(mode='wb')
967 33
        imgf = hdrf
968

969
        # Write main header
970 33
        hdr.write_to(hdrf)
971

972
        # Write every frames
973 33
        for index in range(0, self.header['num_frames']):
974
            # Move to subheader offset
975 33
            frame_offset = subheaders._get_frame_offset(index) - 512
976 33
            imgf.seek(frame_offset)
977

978
            # Write subheader
979 33
            subhdr = subheaders.subheaders[index]
980 33
            imgf.write(subhdr.tobytes())
981

982
            # Seek to the next image block
983 33
            pos = imgf.tell()
984 33
            imgf.seek(pos + 2)
985

986
            # Get frame
987 33
            image = self._subheader.raw_data_from_fileobj(index)
988

989
            # Write frame images
990 33
            self._write_data(image, imgf, pos + 2, endianness='>')
991

992
            # Move to dictionnary offset and write dictionnary entry
993 33
            self._write_data(mlist[index], imgf, entry_pos, endianness='>')
994

995 33
            entry_pos = entry_pos + 16
996

997 33
            current_dir[0] = current_dir[0] - 1
998 33
            current_dir[3] = current_dir[3] + 1
999

1000
            # Create a new directory is previous one is full
1001 33
            if current_dir[0] == 0:
1002
                # self._write_dir(current_dir, imgf, dir_pos)
1003 0
                self._write_data(current_dir, imgf, dir_pos)
1004 0
                current_dir = self._get_empty_dir()
1005 0
                current_dir[3] = dir_pos / 512
1006 0
                dir_pos = mlist[index][2] + 1
1007 0
                entry_pos = dir_pos + 16
1008

1009 33
        tmp_avail = current_dir[0]
1010 33
        tmp_used = current_dir[3]
1011

1012
        # Fill directory with empty data until directory is full
1013 33
        while current_dir[0] > 0:
1014 33
            entry_pos = dir_pos + 16 + (16 * current_dir[3])
1015 33
            self._write_data(np.zeros(4, dtype=np.int32), imgf, entry_pos)
1016 33
            current_dir[0] = current_dir[0] - 1
1017 33
            current_dir[3] = current_dir[3] + 1
1018

1019 33
        current_dir[0] = tmp_avail
1020 33
        current_dir[3] = tmp_used
1021

1022
        # Write directory index
1023 33
        self._write_data(current_dir, imgf, dir_pos, endianness='>')
1024

1025 33
    @classmethod
1026 5
    def from_image(klass, img):
1027 0
        raise NotImplementedError("Ecat images can only be generated "
1028
                                  "from file objects")
1029

1030 33
    @classmethod
1031 5
    def load(klass, filespec):
1032 33
        return klass.from_filename(filespec)
1033

1034

1035 33
load = EcatImage.load

Read our documentation on viewing source code .

Loading