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 ECAT format images 10 11An ECAT format image consists of: 12 13* a *main header*; 14* at least one *matrix list* (mlist); 15 16ECAT thinks of memory locations in terms of *blocks*. One block is 512 17bytes. Thus block 1 starts at 0 bytes, block 2 at 512 bytes, and so on. 18 19The matrix list is an array with one row per frame in the data. 20 21Columns 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 32There is one sub-header for each image frame (or matrix in the terminology 33above). A sub-header can also be called an *image header*. The sub-header is 34one block (512 bytes), and the frame (image) data follows. 35 36There is very little documentation of the ECAT format, and many of the comments 37in this code come from a combination of trial and error and wild speculation. 38 39XMedcon can read and write ECAT 6 format, and read ECAT 7 format: see 40http://xmedcon.sourceforge.net and the ECAT files in the source of XMedCon, 41currently ``libs/tpc/*ecat*`` and ``source/m-ecat*``. Unfortunately XMedCon is 42GPL and some of the header files are adapted from CTI files (called CTI code 43below). It's not clear what the licenses are for these files. 44""" 45 46import warnings 47from numbers import Integral 48 49import numpy as np 50 51from .volumeutils import (native_code, swapped_code, make_dt_codes, 52 array_from_file) 53from .spatialimages import SpatialImage 54from .arraywriters import make_array_writer 55from .wrapstruct import WrapStruct 56from .fileslice import canonical_slicers, predict_shape, slice2outax 57from .deprecated import deprecate_with_version 58 59BLOCK_SIZE = 512 60 61main_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] 123hdr_dtype = np.dtype(main_header_dtd) 124 125 126subheader_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')] 188subhdr_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_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)) 204data_type_codes = make_dt_codes(_dtdefs) 205 206 207# Matrix File Types 208ft_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')) 224file_type_codes = dict(ft_defs) 225 226patient_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')) 236patient_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 240patient_orient_radiological = [0, 2, 4, 6] 241patient_orient_neurological = [1, 3, 5, 7] 242 243 244class 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 template_dtype = hdr_dtype 260 _ft_codes = file_type_codes 261 _patient_orient_codes = patient_orient_codes 262 263 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 super(EcatHeader, self).__init__(binaryblock, endianness, check) 282 283 @classmethod 284 def guessed_endian(klass, hdr): 285 """Guess endian from MAGIC NUMBER value of header data 286 """ 287 if not hdr['sw_version'] == 74: 288 return swapped_code 289 else: 290 return native_code 291 292 @classmethod 293 def default_structarr(klass, endianness=None): 294 """ Return header data for empty header with given endianness 295 """ 296 hdr_data = super(EcatHeader, klass).default_structarr(endianness) 297 hdr_data['magic_number'] = 'MATRIX72' 298 hdr_data['sw_version'] = 74 299 hdr_data['num_frames'] = 0 300 hdr_data['file_type'] = 0 # Unknown 301 hdr_data['ecat_calibration_factor'] = 1.0 # scale factor 302 return hdr_data 303 304 def get_data_dtype(self): 305 """ Get numpy dtype for data from header""" 306 raise NotImplementedError("dtype is only valid from subheaders") 307 308 def get_patient_orient(self): 309 """ gets orientation of patient based on code stored 310 in header, not always reliable 311 """ 312 code = self._structarr['patient_orientation'].item() 313 if code not in self._patient_orient_codes: 314 raise KeyError('Ecat Orientation CODE %d not recognized' % code) 315 return self._patient_orient_codes[code] 316 317 def get_filetype(self): 318 """ Type of ECAT Matrix File from code stored in header""" 319 code = self._structarr['file_type'].item() 320 if code not in self._ft_codes: 321 raise KeyError('Ecat Filetype CODE %d not recognized' % code) 322 return self._ft_codes[code] 323 324 @classmethod 325 def _get_checks(klass): 326 """ Return sequence of check functions for this class """ 327 return () 328 329 330def 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 dt = np.dtype(np.int32) 375 if endianness is not native_code: 376 dt = dt.newbyteorder(endianness) 377 mlists = [] 378 mlist_index = 0 379 mlist_block_no = 2 # 1-based indexing, block with first mlist 380 while True: 381 # Read block containing mlist entries 382 fileobj.seek((mlist_block_no - 1) * BLOCK_SIZE) # fix 1-based indexing 383 dat = fileobj.read(BLOCK_SIZE) 384 rows = np.ndarray(shape=(32, 4), dtype=dt, buffer=dat) 385 # First row special, points to next mlist entries if present 386 n_unused, mlist_block_no, _, n_rows = rows[0] 387 if not (n_unused + n_rows) == 31: # Some error condition here? 388 mlist = [] 389 return mlist 390 # Use all but first housekeeping row 391 mlists.append(rows[1:n_rows + 1]) 392 mlist_index += n_rows 393 if mlist_block_no <= 2: # should block_no in (1, 2) be an error? 394 break 395 return np.row_stack(mlists) 396 397 398def 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 ids = mlist[:, 0].copy() 423 n_valid = np.sum(ids > 0) 424 ids[ids <= 0] = ids.max() + 1 # put invalid frames at end after sort 425 valid_order = np.argsort(ids) 426 if not all(valid_order == sorted(valid_order)): 427 # raise UserWarning if Frames stored out of order 428 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 id_dict = {} 432 for i in range(n_valid): 433 id_dict[i] = [valid_order[i], ids[valid_order[i]]] 434 return id_dict 435 436 437def 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 nframes = len(mlist) 465 frames_order = get_frame_order(mlist) 466 mlist_nframes = len(frames_order) 467 trueframenumbers = np.arange(nframes - mlist_nframes, nframes) 468 frame_dict = {} 469 for frame_stored, (true_order, _) in frames_order.items(): 470 # frame as stored in file -> true number in series 471 try: 472 frame_dict[frame_stored] = trueframenumbers[true_order] + 1 473 except IndexError: 474 raise IOError('Error in header or mlist order unknown') 475 return frame_dict 476 477 478def 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 subheaders = [] 500 dt = subhdr_dtype 501 if endianness is not native_code: 502 dt = dt.newbyteorder(endianness) 503 for mat_id, sh_blkno, sh_last_blkno, mat_stat in mlist: 504 if sh_blkno == 0: 505 break 506 offset = (sh_blkno - 1) * BLOCK_SIZE 507 fileobj.seek(offset) 508 tmpdat = fileobj.read(BLOCK_SIZE) 509 sh = np.ndarray(shape=(), dtype=dt, buffer=tmpdat) 510 subheaders.append(sh) 511 return subheaders 512 513 514class EcatSubHeader(object): 515 516 _subhdrdtype = subhdr_dtype 517 _data_type_codes = data_type_codes 518 519 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 self._header = hdr 533 self.endianness = hdr.endianness 534 self._mlist = mlist 535 self.fileobj = fileobj 536 self.subheaders = read_subheaders(fileobj, mlist, hdr.endianness) 537 538 def get_shape(self, frame=0): 539 """ returns shape of given frame""" 540 subhdr = self.subheaders[frame] 541 x = subhdr['x_dimension'].item() 542 y = subhdr['y_dimension'].item() 543 z = subhdr['z_dimension'].item() 544 return x, y, z 545 546 def get_nframes(self): 547 """returns number of frames""" 548 framed = get_frame_order(self._mlist) 549 return len(framed) 550 551 def _check_affines(self): 552 """checks if all affines are equal across frames""" 553 nframes = self.get_nframes() 554 if nframes == 1: 555 return True 556 affs = [self.get_frame_affine(i) for i in range(nframes)] 557 if affs: 558 i = iter(affs) 559 first = next(i) 560 for item in i: 561 if not np.allclose(first, item): 562 return False 563 return True 564 565 def get_frame_affine(self, frame=0): 566 """returns best affine for given frame of data""" 567 subhdr = self.subheaders[frame] 568 x_off = subhdr['x_offset'] 569 y_off = subhdr['y_offset'] 570 z_off = subhdr['z_offset'] 571 572 zooms = self.get_zooms(frame=frame) 573 574 dims = self.get_shape(frame) 575 # get translations from center of image 576 origin_offset = (np.array(dims) - 1) / 2.0 577 aff = np.diag(zooms) 578 aff[:3, -1] = -origin_offset * zooms[:-1] + np.array([x_off, y_off, 579 z_off]) 580 return aff 581 582 def get_zooms(self, frame=0): 583 """returns zooms ...pixdims""" 584 subhdr = self.subheaders[frame] 585 x_zoom = subhdr['x_pixel_size'] * 10 586 y_zoom = subhdr['y_pixel_size'] * 10 587 z_zoom = subhdr['z_pixel_size'] * 10 588 return (x_zoom, y_zoom, z_zoom, 1) 589 590 def _get_data_dtype(self, frame): 591 dtcode = self.subheaders[frame]['data_type'].item() 592 return self._data_type_codes.dtype[dtcode] 593 594 def _get_frame_offset(self, frame=0): 595 return int(self._mlist[frame][1] * BLOCK_SIZE) 596 597 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 if orientation is None: 608 orientation = self._header['patient_orientation'] 609 elif orientation == 'neurological': 610 orientation = patient_orient_neurological[0] 611 elif orientation == 'radiological': 612 orientation = patient_orient_radiological[0] 613 else: 614 raise ValueError('orientation should be None,\ 615 neurological or radiological') 616 617 if orientation in patient_orient_neurological: 618 raw_data = raw_data[::-1, ::-1, ::-1] 619 elif orientation in patient_orient_radiological: 620 raw_data = raw_data[::, ::-1, ::-1] 621 622 return raw_data 623 624 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 dtype = self._get_data_dtype(frame) 635 if self._header.endianness is not native_code: 636 dtype = dtype.newbyteorder(self._header.endianness) 637 shape = self.get_shape(frame) 638 offset = self._get_frame_offset(frame) 639 fid_obj = self.fileobj 640 raw_data = array_from_file(shape, dtype, fid_obj, offset=offset) 641 raw_data = self._get_oriented_data(raw_data, orientation) 642 return raw_data 643 644 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 header = self._header 655 subhdr = self.subheaders[frame] 656 raw_data = self.raw_data_from_fileobj(frame, orientation) 657 # Scale factors have to be set to scalars to force scalar upcasting 658 data = raw_data * header['ecat_calibration_factor'].item() 659 data = data * subhdr['scale_factor'].item() 660 return data 661 662 663class 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 def __init__(self, subheader): 671 self._subheader = subheader 672 self._data = None 673 x, y, z = subheader.get_shape() 674 nframes = subheader.get_nframes() 675 self._shape = (x, y, z, nframes) 676 677 @property 678 def shape(self): 679 return self._shape 680 681 @property 682 def ndim(self): 683 return len(self.shape) 684 685 @property 686 def is_proxy(self): 687 return True 688 689 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 data = np.empty(self.shape) 709 frame_mapping = get_frame_order(self._subheader._mlist) 710 for i in sorted(frame_mapping): 711 data[:, :, :, i] = self._subheader.data_from_fileobj( 712 frame_mapping[i][0]) 713 if dtype is not None: 714 data = data.astype(dtype, copy=False) 715 return data 716 717 def __getitem__(self, sliceobj): 718 """ Return slice `sliceobj` from ECAT data, optimizing if possible 719 """ 720 sliceobj = canonical_slicers(sliceobj, self.shape) 721 # Indices into sliceobj referring to image axes 722 ax_inds = [i for i, obj in enumerate(sliceobj) if obj is not None] 723 assert len(ax_inds) == len(self.shape) 724 frame_mapping = get_frame_order(self._subheader._mlist) 725 # Analyze index for 4th axis 726 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 in_slicer = sliceobj[:ax_inds[3]] + sliceobj[ax_inds[3] + 1:] 730 # int index for 4th axis, load one slice 731 if isinstance(slice3, Integral): 732 data = self._subheader.data_from_fileobj(frame_mapping[slice3][0]) 733 return data[in_slicer] 734 # slice axis for 4th axis, we will iterate over slices 735 out_shape = predict_shape(sliceobj, self.shape) 736 out_data = np.empty(out_shape) 737 # Slice into output data with out_slicer 738 out_slicer = [slice(None)] * len(out_shape) 739 # Work out axis corresponding to volume in output 740 in2out_ind = slice2outax(len(self.shape), sliceobj)[3] 741 # Iterate over specified 4th axis indices 742 for i in list(range(self.shape[3]))[slice3]: 743 data = self._subheader.data_from_fileobj( 744 frame_mapping[i][0]) 745 out_slicer[in2out_ind] = i 746 out_data[tuple(out_slicer)] = data[in_slicer] 747 return out_data 748 749 750class EcatImage(SpatialImage): 751 """ Class returns a list of Ecat images, with one image(hdr/data) per frame 752 """ 753 _header = EcatHeader 754 header_class = _header 755 valid_exts = ('.v',) 756 _subheader = EcatSubHeader 757 files_types = (('image', '.v'), ('header', '.v')) 758 759 ImageArrayProxy = EcatImageArrayProxy 760 761 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 self._subheader = subheader 806 self._mlist = mlist 807 self._dataobj = dataobj 808 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 affine = np.array(affine, dtype=np.float64, copy=True) 813 if not affine.shape == (4, 4): 814 raise ValueError('Affine should be shape 4,4') 815 self._affine = affine 816 if extra is None: 817 extra = {} 818 self.extra = extra 819 self._header = header 820 if file_map is None: 821 file_map = self.__class__.make_file_map() 822 self.file_map = file_map 823 self._data_cache = None 824 self._fdata_cache = None 825 826 @property 827 def affine(self): 828 if not self._subheader._check_affines(): 829 warnings.warn('Affines different across frames, loading affine ' 830 'from FIRST frame', UserWarning) 831 return self._affine 832 833 def get_frame_affine(self, frame): 834 """returns 4X4 affine""" 835 return self._subheader.get_frame_affine(frame=frame) 836 837 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 return self._subheader.data_from_fileobj(frame, orientation) 846 847 def get_data_dtype(self, frame): 848 subhdr = self._subheader 849 dt = subhdr._get_data_dtype(frame) 850 return dt 851 852 @property 853 def shape(self): 854 x, y, z = self._subheader.get_shape() 855 nframes = self._subheader.get_nframes() 856 return(x, y, z, nframes) 857 858 def get_mlist(self): 859 """ get access to the mlist 860 """ 861 return self._mlist 862 863 def get_subheaders(self): 864 """get access to subheaders""" 865 return self._subheader 866 867 @classmethod 868 @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 def from_filespec(klass, filespec): 873 return klass.from_filename(filespec) 874 875 @staticmethod 876 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 return file_map['header'], file_map['image'] 886 887 @classmethod 888 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 hdr_file, img_file = klass._get_fileholders(file_map) 893 # note header and image are in same file 894 hdr_fid = hdr_file.get_prepare_fileobj(mode='rb') 895 header = klass._header.from_fileobj(hdr_fid) 896 hdr_copy = header.copy() 897 # LOAD MLIST 898 mlist = np.zeros((header['num_frames'], 4), dtype=np.int32) 899 mlist_data = read_mlist(hdr_fid, hdr_copy.endianness) 900 mlist[:len(mlist_data)] = mlist_data 901 # LOAD SUBHEADERS 902 subheaders = klass._subheader(hdr_copy, mlist, hdr_fid) 903 # LOAD DATA 904 # Class level ImageArrayProxy 905 data = klass.ImageArrayProxy(subheaders) 906 # Get affine 907 if not subheaders._check_affines(): 908 warnings.warn('Affines different across frames, loading affine ' 909 'from FIRST frame', UserWarning) 910 aff = subheaders.get_frame_affine() 911 img = klass(data, aff, header, subheaders, mlist, 912 extra=None, file_map=file_map) 913 return img 914 915 def _get_empty_dir(self): 916 """ 917 Get empty directory entry of the form 918 [numAvail, nextDir, previousDir, numUsed] 919 """ 920 return np.array([31, 2, 0, 0], dtype=np.int32) 921 922 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 if dtype is None: 932 dtype = data.dtype 933 934 if endianness is None: 935 endianness = native_code 936 937 stream.seek(pos) 938 make_array_writer(data.newbyteorder(endianness), 939 dtype).to_fileobj(stream) 940 941 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 if file_map is None: 953 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 self.get_fdata() 958 hdr = self.header 959 mlist = self._mlist 960 subheaders = self.get_subheaders() 961 dir_pos = 512 962 entry_pos = dir_pos + 16 # 528 963 current_dir = self._get_empty_dir() 964 965 hdr_fh, img_fh = self._get_fileholders(file_map) 966 hdrf = hdr_fh.get_prepare_fileobj(mode='wb') 967 imgf = hdrf 968 969 # Write main header 970 hdr.write_to(hdrf) 971 972 # Write every frames 973 for index in range(0, self.header['num_frames']): 974 # Move to subheader offset 975 frame_offset = subheaders._get_frame_offset(index) - 512 976 imgf.seek(frame_offset) 977 978 # Write subheader 979 subhdr = subheaders.subheaders[index] 980 imgf.write(subhdr.tobytes()) 981 982 # Seek to the next image block 983 pos = imgf.tell() 984 imgf.seek(pos + 2) 985 986 # Get frame 987 image = self._subheader.raw_data_from_fileobj(index) 988 989 # Write frame images 990 self._write_data(image, imgf, pos + 2, endianness='>') 991 992 # Move to dictionnary offset and write dictionnary entry 993 self._write_data(mlist[index], imgf, entry_pos, endianness='>') 994 995 entry_pos = entry_pos + 16 996 997 current_dir[0] = current_dir[0] - 1 998 current_dir[3] = current_dir[3] + 1 999 1000 # Create a new directory is previous one is full 1001 if current_dir[0] == 0: 1002 # self._write_dir(current_dir, imgf, dir_pos) 1003 self._write_data(current_dir, imgf, dir_pos) 1004 current_dir = self._get_empty_dir() 1005 current_dir[3] = dir_pos / 512 1006 dir_pos = mlist[index][2] + 1 1007 entry_pos = dir_pos + 16 1008 1009 tmp_avail = current_dir[0] 1010 tmp_used = current_dir[3] 1011 1012 # Fill directory with empty data until directory is full 1013 while current_dir[0] > 0: 1014 entry_pos = dir_pos + 16 + (16 * current_dir[3]) 1015 self._write_data(np.zeros(4, dtype=np.int32), imgf, entry_pos) 1016 current_dir[0] = current_dir[0] - 1 1017 current_dir[3] = current_dir[3] + 1 1018 1019 current_dir[0] = tmp_avail 1020 current_dir[3] = tmp_used 1021 1022 # Write directory index 1023 self._write_data(current_dir, imgf, dir_pos, endianness='>') 1024 1025 @classmethod 1026 def from_image(klass, img): 1027 raise NotImplementedError("Ecat images can only be generated " 1028 "from file objects") 1029 1030 @classmethod 1031 def load(klass, filespec): 1032 return klass.from_filename(filespec) 1033 1034 1035load = EcatImage.load 1036