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### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 9import warnings 10 11import numpy as np 12 13from .volumeutils import (native_code, swapped_code, make_dt_codes, 14 array_from_file) 15from .spatialimages import SpatialImage, ImageDataError 16from .arraywriters import make_array_writer 17 18 19MAINHDRSZ = 502 20main_header_dtd = [ 21 ('magic_number', '14S'), 22 ('original_filename', '32S'), 23 ('sw_version', np.uint16), 24 ('system_type', np.uint16), 25 ('file_type', np.uint16), 26 ('serial_number', '10S'), 27 ('scan_start_time',np.uint32), 28 ('isotope_name', '8S'), 29 ('isotope_halflife', np.float32), 30 ('radiopharmaceutical','32S'), 31 ('gantry_tilt', np.float32), 32 ('gantry_rotation',np.float32), 33 ('bed_elevation',np.float32), 34 ('intrinsic_tilt', np.float32), 35 ('wobble_speed',np.uint16), 36 ('transm_source_type',np.uint16), 37 ('distance_scanned',np.float32), 38 ('transaxial_fov',np.float32), 39 ('angular_compression', np.uint16), 40 ('coin_samp_mode',np.uint16), 41 ('axial_samp_mode',np.uint16), 42 ('ecat_calibration_factor',np.float32), 43 ('calibration_unitS', np.uint16), 44 ('calibration_units_type',np.uint16), 45 ('compression_code',np.uint16), 46 ('study_type','12S'), 47 ('patient_id','16S'), 48 ('patient_name','32S'), 49 ('patient_sex','1S'), 50 ('patient_dexterity','1S'), 51 ('patient_age',np.float32), 52 ('patient_height',np.float32), 53 ('patient_weight',np.float32), 54 ('patient_birth_date',np.uint32), 55 ('physician_name','32S'), 56 ('operator_name','32S'), 57 ('study_description','32S'), 58 ('acquisition_type',np.uint16), 59 ('patient_orientation',np.uint16), 60 ('facility_name', '20S'), 61 ('num_planes',np.uint16), 62 ('num_frames',np.uint16), 63 ('num_gates',np.uint16), 64 ('num_bed_pos',np.uint16), 65 ('init_bed_position',np.float32), 66 ('bed_position','15f'), 67 ('plane_separation',np.float32), 68 ('lwr_sctr_thres',np.uint16), 69 ('lwr_true_thres',np.uint16), 70 ('upr_true_thres',np.uint16), 71 ('user_process_code','10S'), 72 ('acquisition_mode',np.uint16), 73 ('bin_size',np.float32), 74 ('branching_fraction',np.float32), 75 ('dose_start_time',np.uint32), 76 ('dosage',np.float32), 77 ('well_counter_corr_factor', np.float32), 78 ('data_units', '32S'), 79 ('septa_state',np.uint16), 80 ('fill', '12S') 81 ] 82hdr_dtype = np.dtype(main_header_dtd) 83 84 85subheader_dtd = [ 86 ('data_type', np.uint16), 87 ('num_dimensions', np.uint16), 88 ('x_dimension', np.uint16), 89 ('y_dimension', np.uint16), 90 ('z_dimension', np.uint16), 91 ('x_offset', np.float32), 92 ('y_offset', np.float32), 93 ('z_offset', np.float32), 94 ('recon_zoom', np.float32), 95 ('scale_factor', np.float32), 96 ('image_min', np.int16), 97 ('image_max', np.int16), 98 ('x_pixel_size', np.float32), 99 ('y_pixel_size', np.float32), 100 ('z_pixel_size', np.float32), 101 ('frame_duration', np.uint32), 102 ('frame_start_time', np.uint32), 103 ('filter_code', np.uint16), 104 ('x_resolution', np.float32), 105 ('y_resolution', np.float32), 106 ('z_resolution', np.float32), 107 ('num_r_elements', np.float32), 108 ('num_angles', np.float32), 109 ('z_rotation_angle', np.float32), 110 ('decay_corr_fctr', np.float32), 111 ('corrections_applied', np.uint32), 112 ('gate_duration', np.uint32), 113 ('r_wave_offset', np.uint32), 114 ('num_accepted_beats', np.uint32), 115 ('filter_cutoff_frequency', np.float32), 116 ('filter_resolution', np.float32), 117 ('filter_ramp_slope', np.float32), 118 ('filter_order', np.uint16), 119 ('filter_scatter_fraction', np.float32), 120 ('filter_scatter_slope', np.float32), 121 ('annotation', '40S'), 122 ('mt_1_1', np.float32), 123 ('mt_1_2', np.float32), 124 ('mt_1_3', np.float32), 125 ('mt_2_1', np.float32), 126 ('mt_2_2', np.float32), 127 ('mt_2_3', np.float32), 128 ('mt_3_1', np.float32), 129 ('mt_3_2', np.float32), 130 ('mt_3_3', np.float32), 131 ('rfilter_cutoff', np.float32), 132 ('rfilter_resolution', np.float32), 133 ('rfilter_code', np.uint16), 134 ('rfilter_order', np.uint16), 135 ('zfilter_cutoff', np.float32), 136 ('zfilter_resolution',np.float32), 137 ('zfilter_code', np.uint16), 138 ('zfilter_order', np.uint16), 139 ('mt_4_1', np.float32), 140 ('mt_4_2', np.float32), 141 ('mt_4_3', np.float32), 142 ('scatter_type', np.uint16), 143 ('recon_type', np.uint16), 144 ('recon_views', np.uint16), 145 ('fill', '174S'), 146 ('fill2', '96S')] 147subhdr_dtype = np.dtype(subheader_dtd) 148 149# Ecat Data Types 150_dtdefs = ( # code, name, equivalent dtype 151 (1, 'ECAT7_BYTE', np.uint8), 152 (2, 'ECAT7_VAXI2', np.int16), 153 (3, 'ECAT7_VAXI4', np.float32), 154 (4, 'ECAT7_VAXR4', np.float32), 155 (5, 'ECAT7_IEEER4', np.float32), 156 (6, 'ECAT7_SUNI2', np.uint16), 157 (7, 'ECAT7_SUNI4', np.int32)) 158data_type_codes = make_dt_codes(_dtdefs) 159 160 161# Matrix File Types 162ft_defs = ( # code, name 163 (0, 'ECAT7_UNKNOWN'), 164 (1, 'ECAT7_2DSCAN'), 165 (2, 'ECAT7_IMAGE16'), 166 (3, 'ECAT7_ATTEN'), 167 (4, 'ECAT7_2DNORM'), 168 (5, 'ECAT7_POLARMAP'), 169 (6, 'ECAT7_VOLUME8'), 170 (7, 'ECAT7_VOLUME16'), 171 (8, 'ECAT7_PROJ'), 172 (9, 'ECAT7_PROJ16'), 173 (10, 'ECAT7_IMAGE8'), 174 (11, 'ECAT7_3DSCAN'), 175 (12, 'ECAT7_3DSCAN8'), 176 (13, 'ECAT7_3DNORM'), 177 (14, 'ECAT7_3DSCANFIT')) 178 179patient_orient_defs = ( #code, description 180 (0, 'ECAT7_Feet_First_Prone'), 181 (1, 'ECAT7_Head_First_Prone'), 182 (2, 'ECAT7_Feet_First_Supine'), 183 (3, 'ECAT7_Head_First_Supine'), 184 (4, 'ECAT7_Feet_First_Decubitus_Right'), 185 (5, 'ECAT7_Head_First_Decubitus_Right'), 186 (6, 'ECAT7_Feet_First_Decubitus_Left'), 187 (7, 'ECAT7_Head_First_Decubitus_Left'), 188 (8, 'ECAT7_Unknown_Orientation')) 189 190#Indexes from the patient_orient_defs structure defined above for the 191#neurological and radiological viewing conventions 192patient_orient_radiological = [0, 2, 4, 6] 193patient_orient_neurological = [1, 3, 5, 7] 194 195class EcatHeader(object): 196 """Class for basic Ecat PET header 197 Sub-parts of standard Ecat File 198 main header 199 matrix list 200 which lists the information for each 201 frame collected (can have 1 to many frames) 202 subheaders specific to each frame 203 with possibly-variable sized data blocks 204 205 This just reads the main Ecat Header, 206 it does not load the data 207 or read the mlist or any sub headers 208 209 """ 210 211 _dtype = hdr_dtype 212 _ft_defs = ft_defs 213 _patient_orient_defs = patient_orient_defs 214 215 def __init__(self, 216 fileobj=None, 217 endianness=None): 218 """Initialize Ecat header from file object 219 220 Parameters 221 ---------- 222 fileobj : {None, string} optional 223 binary block to set into header, By default, None 224 in which case we insert default empty header block 225 endianness : {None, '<', '>', other endian code}, optional 226 endian code of binary block, If None, guess endianness 227 from the data 228 """ 229 if fileobj is None: 230 self._header_data = self._empty_headerdata(endianness) 231 return 232 233 hdr = np.ndarray(shape=(), 234 dtype=self._dtype, 235 buffer=fileobj) 236 if endianness is None: 237 endianness = self._guess_endian(hdr) 238 239 if endianness != native_code: 240 dt = self._dtype.newbyteorder(endianness) 241 hdr = np.ndarray(shape=(), 242 dtype=dt, 243 buffer=fileobj) 244 self._header_data = hdr.copy() 245 246 return 247 248 def get_header(self): 249 """returns header """ 250 return self 251 252 @property 253 def binaryblock(self): 254 return self._header_data.tostring() 255 256 @property 257 def endianness(self): 258 if self._header_data.dtype.isnative: 259 return native_code 260 return swapped_code 261 262 263 def _guess_endian(self, hdr): 264 """Guess endian from MAGIC NUMBER value of header data 265 """ 266 if not hdr['sw_version'] == 74: 267 return swapped_code 268 else: 269 return native_code 270 271 @classmethod 272 def from_fileobj(klass, fileobj, endianness=None): 273 """Return /read header with given or guessed endian code 274 275 Parameters 276 ---------- 277 fileobj : file-like object 278 Needs to implement ``read`` method 279 endianness : None or endian code, optional 280 Code specifying endianness of data to be read 281 282 Returns 283 ------- 284 hdr : EcatHeader object 285 EcatHeader object initialized from data in file object 286 287 Examples 288 -------- 289 290 291 """ 292 raw_str = fileobj.read(klass._dtype.itemsize) 293 return klass(raw_str, endianness) 294 295 def write_to(self, fileobj): 296 fileobj.write(self.binaryblock) 297 298 def _empty_headerdata(self,endianness=None): 299 """Return header data for empty header with given endianness""" 300 #hdr_data = super(EcatHeader, self)._empty_headerdata(endianness) 301 dt = self._dtype 302 if not endianness is None: 303 dt = dt.newbyteorder(endianness) 304 hdr_data = np.zeros((), dtype=dt) 305 hdr_data['magic_number'] = 'MATRIX72' 306 hdr_data['sw_version'] = 74 307 hdr_data['num_frames']= 0 308 hdr_data['file_type'] = 0 # Unknown 309 hdr_data['ecat_calibration_factor'] = 1.0 # scale factor 310 return hdr_data 311 312 313 def get_data_dtype(self): 314 """ Get numpy dtype for data from header""" 315 raise NotImplementedError("dtype is only valid from subheaders") 316 317 318 def copy(self): 319 return self.__class__( 320 self.binaryblock, 321 self.endianness) 322 323 324 def __eq__(self, other): 325 """ checks for equality between two headers""" 326 self_end = self.endianness 327 self_bb = self.binaryblock 328 if self_end == other.endianness: 329 return self_bb == other.binaryblock 330 other_bb = other._header_data.byteswap().tostring() 331 return self_bb == other_bb 332 333 def __ne__(self, other): 334 ''' equality between two headers defined by ``header_data`` 335 336 For examples, see ``__eq__`` method docstring 337 ''' 338 return not self == other 339 340 def __getitem__(self, item): 341 ''' Return values from header data 342 343 Examples 344 -------- 345 >>> hdr = EcatHeader() 346 >>> hdr['magic_number'] #23dt next : bytes 347 'MATRIX72' 348 ''' 349 return self._header_data[item].item() 350 351 def __setitem__(self, item, value): 352 ''' Set values in header data 353 354 Examples 355 -------- 356 >>> hdr = EcatHeader() 357 >>> hdr['num_frames'] = 2 358 >>> hdr['num_frames'] 359 2 360 ''' 361 self._header_data[item] = value 362 363 def get_patient_orient(self): 364 """ gets orientation of patient based on code stored 365 in header, not always reliable""" 366 orient_code = dict(self._patient_orient_defs) 367 code = self._header_data['patient_orientation'].item() 368 if not orient_code.has_key(code): 369 raise KeyError('Ecat Orientation CODE %d not recognized'%code) 370 return orient_code[code] 371 372 def get_filetype(self): 373 """ gets type of ECAT Matrix File from 374 code stored in header""" 375 ft_codes = dict(self._ft_defs) 376 code = self._header_data['file_type'].item() 377 if not ft_codes.has_key(code): 378 raise KeyError('Ecat Filetype CODE %d not recognized'%code) 379 return ft_codes[code] 380 381 def __iter__(self): 382 return iter(self.keys()) 383 384 def keys(self): 385 ''' Return keys from header data''' 386 return list(self._dtype.names) 387 388 def values(self): 389 ''' Return values from header data''' 390 data = self._header_data 391 return [data[key] for key in self._dtype.names] 392 393 def items(self): 394 ''' Return items from header data''' 395 return zip(self.keys(), self.values()) 396 397class EcatMlist(object): 398 399 def __init__(self,fileobj, hdr): 400 """ gets list of frames and subheaders in pet file 401 402 Parameters 403 ----------- 404 fileobj : ECAT file <filename>.v fileholder or file object 405 with read, seek methods 406 407 Returns 408 ------- 409 mlist : numpy recarray nframes X 4 columns 410 1 - Matrix identifier. 411 2 - subheader record number 412 3 - Last record number of matrix data block. 413 4 - Matrix status: 414 1 - exists - rw 415 2 - exists - ro 416 3 - matrix deleted 417 """ 418 self.hdr = hdr 419 self._mlist = self.get_mlist(fileobj) 420 421 def get_mlist(self, fileobj): 422 fileobj.seek(512) 423 dat=fileobj.read(128*32) 424 425 dt = np.dtype([('matlist',np.int32)]) 426 if not self.hdr.endianness is native_code: 427 dt = dt.newbyteorder(self.hdr.endianness) 428 nframes = self.hdr['num_frames'] 429 mlist = np.zeros((nframes,4), dtype='uint32') 430 record_count = 0 431 done = False 432 433 while not done: #mats['matlist'][0,1] == 2: 434 435 mats = np.recarray(shape=(32,4), dtype=dt, buf=dat) 436 if not (mats['matlist'][0,0] + mats['matlist'][0,3]) == 31: 437 mlist = [] 438 return mlist 439 440 nrecords = mats['matlist'][0,3] 441 mlist[record_count:nrecords+record_count,:] = mats['matlist'][1:nrecords+1,:] 442 record_count+= nrecords 443 if mats['matlist'][0,1] == 2 or mats['matlist'][0,1] == 0: 444 done = True 445 else: 446 # Find next subheader 447 tmp = int(mats['matlist'][0,1]-1)#cast to int 448 fileobj.seek(0) 449 fileobj.seek(tmp*512) 450 dat = fileobj.read(128*32) 451 452 return mlist 453 454 def get_frame_order(self): 455 """Returns the order of the frames stored in the file 456 Sometimes Frames are not stored in the file in 457 chronological order, this can be used to extract frames 458 in correct order 459 460 Returns 461 ------- 462 id_dict: dict mapping frame number -> [mlist_row, mlist_id] 463 464 (where mlist id is value in the first column of the mlist matrix ) 465 466 Examples 467 -------- 468 >>> import os 469 >>> import nibabel as nib 470 >>> nibabel_dir = os.path.dirname(nib.__file__) 471 >>> from nibabel import ecat 472 >>> ecat_file = os.path.join(nibabel_dir,'tests','data','tinypet.v') 473 >>> img = ecat.load(ecat_file) 474 >>> mlist = img.get_mlist() 475 >>> mlist.get_frame_order() 476 {0: [0, 16842758]} 477 """ 478 mlist = self._mlist 479 ids = mlist[:, 0].copy() 480 n_valid = np.sum(ids > 0) 481 ids[ids <=0] = ids.max() + 1 # put invalid frames at end after sort 482 valid_order = np.argsort(ids) 483 if not all(valid_order == sorted(valid_order)): 484 #raise UserWarning if Frames stored out of order 485 warnings.warn_explicit('Frames stored out of order;'\ 486 'true order = %s\n'\ 487 'frames will be accessed in order '\ 488 'STORED, NOT true order'%(valid_order), 489 UserWarning,'ecat', 0) 490 id_dict = {} 491 for i in range(n_valid): 492 id_dict[i] = [valid_order[i], ids[valid_order[i]]] 493 494 return id_dict 495 496 def get_series_framenumbers(self): 497 """ Returns framenumber of data as it was collected, 498 as part of a series; not just the order of how it was 499 stored in this or across other files 500 501 For example, if the data is split between multiple files 502 this should give you the true location of this frame as 503 collected in the series 504 (Frames are numbered starting at ONE (1) not Zero) 505 506 Returns 507 ------- 508 frame_dict: dict mapping order_stored -> frame in series 509 where frame in series counts from 1; [1,2,3,4...] 510 511 Examples 512 -------- 513 >>> import os 514 >>> import nibabel as nib 515 >>> nibabel_dir = os.path.dirname(nib.__file__) 516 >>> from nibabel import ecat 517 >>> ecat_file = os.path.join(nibabel_dir,'tests','data','tinypet.v') 518 >>> img = ecat.load(ecat_file) 519 >>> mlist = img.get_mlist() 520 >>> mlist.get_series_framenumbers() 521 {0: 1} 522 523 524 525 """ 526 frames_order = self.get_frame_order() 527 nframes = self.hdr['num_frames'] 528 mlist_nframes = len(frames_order) 529 trueframenumbers = np.arange(nframes - mlist_nframes, nframes) 530 frame_dict = {} 531 try: 532 for frame_stored, (true_order, _) in frames_order.items(): 533 #frame as stored in file -> true number in series 534 frame_dict[frame_stored] = trueframenumbers[true_order]+1 535 return frame_dict 536 except: 537 raise IOError('Error in header or mlist order unknown') 538 539class EcatSubHeader(object): 540 541 _subhdrdtype = subhdr_dtype 542 _data_type_codes = data_type_codes 543 544 def __init__(self, hdr, mlist, fileobj): 545 """parses the subheaders in the ecat (.v) file 546 there is one subheader for each frame in the ecat file 547 548 Parameters 549 ----------- 550 hdr : EcatHeader 551 552 mlist : EcatMlist 553 554 fileobj : ECAT file <filename>.v fileholder or file object 555 with read, seek methods 556 557 558 """ 559 self._header = hdr 560 self.endianness = hdr.endianness 561 self._mlist = mlist 562 self.fileobj = fileobj 563 self.subheaders = self._get_subheaders() 564 565 def _get_subheaders(self): 566 """retreive all subheaders and return list of subheader recarrays 567 """ 568 subheaders = [] 569 header = self._header 570 endianness = self.endianness 571 dt = self._subhdrdtype 572 if not self.endianness is native_code: 573 dt = self._subhdrdtype.newbyteorder(self.endianness) 574 if self._header['num_frames'] > 1: 575 for item in self._mlist._mlist: 576 if item[1] == 0: 577 break 578 self.fileobj.seek(0) 579 offset = (int(item[1])-1)*512 580 self.fileobj.seek(offset) 581 tmpdat = self.fileobj.read(512) 582 sh = (np.recarray(shape=(), dtype=dt, 583 buf=tmpdat)) 584 subheaders.append(sh.copy()) 585 else: 586 self.fileobj.seek(0) 587 offset = (int(self._mlist._mlist[0][1])-1)*512 588 self.fileobj.seek(offset) 589 tmpdat = self.fileobj.read(512) 590 sh = (np.recarray(shape=(), dtype=dt, 591 buf=tmpdat)) 592 subheaders.append(sh) 593 return subheaders 594 595 def get_shape(self, frame=0): 596 """ returns shape of given frame""" 597 subhdr = self.subheaders[frame] 598 x = subhdr['x_dimension'].item() 599 y = subhdr['y_dimension'].item() 600 z = subhdr['z_dimension'].item() 601 return (x,y,z) 602 603 def get_nframes(self): 604 """returns number of frames""" 605 mlist = self._mlist 606 framed = mlist.get_frame_order() 607 return len(framed) 608 609 610 def _check_affines(self): 611 """checks if all affines are equal across frames""" 612 nframes = self.get_nframes() 613 if nframes == 1: 614 return True 615 affs = [self.get_frame_affine(i) for i in range(nframes)] 616 if affs: 617 i = iter(affs) 618 first = i.next() 619 for item in i: 620 if not np.all(first == item): 621 return False 622 return True 623 624 def get_frame_affine(self,frame=0): 625 """returns best affine for given frame of data""" 626 subhdr = self.subheaders[frame] 627 x_off = subhdr['x_offset'] 628 y_off = subhdr['y_offset'] 629 z_off = subhdr['z_offset'] 630 631 zooms = self.get_zooms(frame=frame) 632 633 dims = self.get_shape(frame) 634 # get translations from center of image 635 origin_offset = (np.array(dims)-1) / 2.0 636 aff = np.diag(zooms) 637 aff[:3,-1] = -origin_offset * zooms[:-1] + np.array([x_off,y_off,z_off]) 638 return aff 639 640 def get_zooms(self,frame=0): 641 """returns zooms ...pixdims""" 642 subhdr = self.subheaders[frame] 643 x_zoom = subhdr['x_pixel_size'] * 10 644 y_zoom = subhdr['y_pixel_size'] * 10 645 z_zoom = subhdr['z_pixel_size'] * 10 646 return (x_zoom, y_zoom, z_zoom, 1) 647 648 649 def _get_data_dtype(self, frame): 650 dtcode = self.subheaders[frame]['data_type'].item() 651 return self._data_type_codes.dtype[dtcode] 652 653 def _get_frame_offset(self, frame=0): 654 mlist = self._mlist._mlist 655 offset = (mlist[frame][1]) * 512 656 return int(offset) 657 658 def _get_oriented_data(self, raw_data, orientation=None): 659 ''' 660 Get data oriented following ``patient_orientation`` header field. If the 661 ``orientation`` parameter is given, return data according to this 662 orientation. 663 664 :param raw_data: Numpy array containing the raw data 665 :param orientation: None (default), 'neurological' or 'radiological' 666 :rtype: Numpy array containing the oriented data 667 ''' 668 if orientation is None: 669 orientation = self._header['patient_orientation'] 670 elif orientation == 'neurological': 671 orientation = patient_orient_neurological[0] 672 elif orientation == 'radiological': 673 orientation = patient_orient_radiological[0] 674 else: 675 raise ValueError('orientation should be None,\ 676 neurological or radiological') 677 678 if orientation in patient_orient_neurological: 679 raw_data = raw_data[::-1, ::-1, ::-1] 680 elif orientation in patient_orient_radiological: 681 raw_data = raw_data[::, ::-1, ::-1] 682 683 return raw_data 684 685 def raw_data_from_fileobj(self, frame=0, orientation=None): 686 ''' 687 Get raw data from file object. 688 689 :param frame: Time frame index from where to fetch data 690 :param orientation: None (default), 'neurological' or 'radiological' 691 :rtype: Numpy array containing (possibly oriented) raw data 692 693 .. seealso:: data_from_fileobj 694 ''' 695 dtype = self._get_data_dtype(frame) 696 if not self._header.endianness is native_code: 697 dtype=dtype.newbyteorder(self._header.endianness) 698 shape = self.get_shape(frame) 699 offset = self._get_frame_offset(frame) 700 fid_obj = self.fileobj 701 raw_data = array_from_file(shape, dtype, fid_obj, offset=offset) 702 raw_data = self._get_oriented_data(raw_data, orientation) 703 return raw_data 704 705 def data_from_fileobj(self, frame=0, orientation=None): 706 ''' 707 Read scaled data from file for a given frame 708 709 :param frame: Time frame index from where to fetch data 710 :param orientation: None (default), 'neurological' or 'radiological' 711 :rtype: Numpy array containing (possibly oriented) raw data 712 713 .. seealso:: raw_data_from_fileobj 714 ''' 715 header = self._header 716 subhdr = self.subheaders[frame] 717 raw_data = self.raw_data_from_fileobj(frame, orientation) 718 data = raw_data * header['ecat_calibration_factor'] 719 data = data * subhdr['scale_factor'] 720 return data 721 722 723 724 725class EcatImage(SpatialImage): 726 """This class returns a list of Ecat images, 727 with one image(hdr/data) per frame 728 """ 729 _header = EcatHeader 730 header_class = _header 731 _subheader = EcatSubHeader 732 _mlist = EcatMlist 733 files_types = (('image', '.v'), ('header', '.v')) 734 735 736 class ImageArrayProxy(object): 737 ''' Ecat implemention of array proxy protocol 738 739 The array proxy allows us to freeze the passed fileobj and 740 header such that it returns the expected data array. 741 ''' 742 def __init__(self, subheader): 743 self._subheader = subheader 744 self._data = None 745 x, y, z = subheader.get_shape() 746 nframes = subheader.get_nframes() 747 self.shape = (x, y, z, nframes) 748 749 def __array__(self): 750 ''' Cached read of data from file 751 This reads ALL FRAMES into one array, can be memory expensive 752 use subheader.data_from_fileobj(frame) for less memory intensive 753 reads 754 ''' 755 if self._data is None: 756 self._data = np.empty(self.shape) 757 frame_mapping = self._subheader._mlist.get_frame_order() 758 for i in sorted(frame_mapping): 759 self._data[:,:,:,i] = self._subheader.data_from_fileobj(frame_mapping[i][0]) 760 return self._data 761 762 def __init__(self, data, affine, header, 763 subheader, mlist , 764 extra = None, file_map = None): 765 """ Initialize Image 766 767 The image is a combination of 768 (array, affine matrix, header, subheader, mlist) 769 with optional meta data in `extra`, and filename / file-like objects 770 contained in the `file_map`. 771 772 Parameters 773 ---------- 774 data : None or array-like 775 image data 776 affine : None or (4,4) array-like 777 homogeneous affine giving relationship between voxel coords and 778 world coords. 779 header : None or header instance 780 meta data for this image format 781 subheader : None or subheader instance 782 meta data for each sub-image for frame in the image 783 mlist : None or mlist instance 784 meta data with array giving offset and order of data in file 785 extra : None or mapping, optional 786 metadata associated with this image that cannot be 787 stored in header or subheader 788 file_map : mapping, optional 789 mapping giving file information for this image format 790 791 Examples 792 -------- 793 >>> import os 794 >>> import nibabel as nib 795 >>> nibabel_dir = os.path.dirname(nib.__file__) 796 >>> from nibabel import ecat 797 >>> ecat_file = os.path.join(nibabel_dir,'tests','data','tinypet.v') 798 >>> img = ecat.load(ecat_file) 799 >>> frame0 = img.get_frame(0) 800 >>> frame0.shape == (10, 10, 3) 801 True 802 >>> data4d = img.get_data() 803 >>> data4d.shape == (10, 10, 3, 1) 804 True 805 """ 806 self._subheader = subheader 807 self._mlist = mlist 808 self._data = data 809 if not affine is None: 810 # Check that affine is array-like 4,4. Maybe this is too strict at 811 # this abstract level, but so far I think all image formats we know 812 # do need 4,4. 813 affine = np.asarray(affine) 814 if not affine.shape == (4,4): 815 raise ValueError('Affine should be shape 4,4') 816 self._affine = affine 817 if extra is None: 818 extra = {} 819 self.extra = extra 820 self._header = header 821 if file_map is None: 822 file_map = self.__class__.make_file_map() 823 self.file_map = file_map 824 825 def _set_header(self, header): 826 self._header = header 827 828 def get_data(self): 829 """returns scaled data for all frames in a numpy array 830 returns as a 4D array """ 831 if self._data is None: 832 raise ImageDataError('No data in this image') 833 return np.asanyarray(self._data) 834 835 def get_affine(self): 836 if not self._subheader._check_affines(): 837 warnings.warn('Affines different across frames, loading affine from FIRST frame', 838 UserWarning ) 839 return self._affine 840 841 def get_frame_affine(self, frame): 842 """returns 4X4 affine""" 843 return self._subheader.get_frame_affine(frame=frame) 844 845 def get_frame(self,frame, orientation=None): 846 ''' 847 Get full volume for a time frame 848 849 :param frame: Time frame index from where to fetch data 850 :param orientation: None (default), 'neurological' or 'radiological' 851 :rtype: Numpy array containing (possibly oriented) raw data 852 ''' 853 return self._subheader.data_from_fileobj(frame, orientation) 854 855 def get_data_dtype(self,frame): 856 subhdr = self._subheader 857 dt = subhdr._get_data_dtype(frame) 858 return dt 859 860 @property 861 def shape(self): 862 x,y,z = self._subheader.get_shape() 863 nframes = self._subheader.get_nframes() 864 return(x, y, z, nframes) 865 866 def get_mlist(self): 867 """ get access to the mlist """ 868 return self._mlist 869 870 def get_subheaders(self): 871 """get access to subheaders""" 872 return self._subheader 873 874 @classmethod 875 def from_filespec(klass, filespec): 876 return klass.from_filename(filespec) 877 878 879 @staticmethod 880 def _get_fileholders(file_map): 881 """ returns files specific to header and image of the image 882 for ecat .v this is the same image file 883 884 Returns 885 ------- 886 header : file holding header data 887 image : file holding image data 888 """ 889 return file_map['header'], file_map['image'] 890 891 @classmethod 892 def from_file_map(klass, file_map): 893 """class method to create image from mapping 894 specified in file_map""" 895 hdr_file, img_file = klass._get_fileholders(file_map) 896 #note header and image are in same file 897 hdr_fid = hdr_file.get_prepare_fileobj(mode = 'rb') 898 header = klass._header.from_fileobj(hdr_fid) 899 hdr_copy = header.copy() 900 ### LOAD MLIST 901 mlist = klass._mlist(hdr_fid, hdr_copy) 902 ### LOAD SUBHEADERS 903 subheaders = klass._subheader(hdr_copy, 904 mlist, 905 hdr_fid) 906 ### LOAD DATA 907 ## Class level ImageArrayProxy 908 data = klass.ImageArrayProxy(subheaders) 909 910 ## Get affine 911 if not subheaders._check_affines(): 912 warnings.warn('Affines different across frames, loading affine from FIRST frame', 913 UserWarning ) 914 aff = subheaders.get_frame_affine() 915 img = klass(data, aff, header, subheaders, mlist, extra=None, file_map = file_map) 916 return img 917 918 def _get_empty_dir(self): 919 ''' 920 Get empty directory entry of the form 921 [numAvail, nextDir, previousDir, numUsed] 922 ''' 923 return np.array([31, 2, 0, 0], dtype=np.uint32) 924 925 def _write_data(self, data, stream, pos, dtype=None, endianness=None): 926 ''' 927 Write data to ``stream`` using an array_writer 928 929 :param data: Numpy array containing the dat 930 :param stream: The file-like object to write the data to 931 :param pos: The position in the stream to write the data to 932 :param endianness: Endianness code of the data to write 933 ''' 934 if dtype is None: 935 dtype = data.dtype 936 937 if endianness is None: 938 endianness = native_code 939 940 stream.seek(pos) 941 writer = make_array_writer( 942 data.newbyteorder(endianness), 943 dtype).to_fileobj(stream) 944 945 def to_file_map(self, file_map=None): 946 ''' Write ECAT7 image to `file_map` or contained ``self.file_map`` 947 948 The format consist of: 949 950 - A main header (512L) with dictionary entries in the form 951 [numAvail, nextDir, previousDir, numUsed] 952 - For every frame (3D volume in 4D data) 953 - A subheader (size = frame_offset) 954 - Frame data (3D volume) 955 ''' 956 if file_map is None: 957 file_map = self.file_map 958 959 data = self.get_data() 960 hdr = self.get_header() 961 mlist = self.get_mlist()._mlist 962 subheaders = self.get_subheaders() 963 dir_pos = 512L 964 entry_pos = dir_pos + 16L #528L 965 current_dir = self._get_empty_dir() 966 967 hdr_fh, img_fh = self._get_fileholders(file_map) 968 hdrf = hdr_fh.get_prepare_fileobj(mode='wb') 969 imgf = hdrf 970 971 #Write main header 972 hdr.write_to(hdrf) 973 974 #Write every frames 975 for index in xrange(0, self.get_header()['num_frames']): 976 #Move to subheader offset 977 frame_offset = subheaders._get_frame_offset(index) - 512 978 imgf.seek(frame_offset) 979 980 #Write subheader 981 subhdr = subheaders.subheaders[index] 982 imgf.write(subhdr.tostring()) 983 984 #Seek to the next image block 985 pos = imgf.tell() 986 imgf.seek(pos + 2) 987 988 #Get frame and its data type 989 image = self._subheader.raw_data_from_fileobj(index) 990 dtype = image.dtype 991 992 #Write frame images 993 self._write_data(image, imgf, pos+2, endianness='>') 994 995 #Move to dictionnary offset and write dictionnary entry 996 self._write_data(mlist[index], imgf, entry_pos, 997 np.uint32, endianness='>') 998 999 entry_pos = entry_pos + 16L 1000 1001 current_dir[0] = current_dir[0] - 1 1002 current_dir[3] = current_dir[3] + 1 1003 1004 #Create a new directory is previous one is full 1005 if current_dir[0] == 0: 1006 #self._write_dir(current_dir, imgf, dir_pos) 1007 self._write_data(current_dir, imgf, dir_pos) 1008 current_dir = self._get_empty_dir() 1009 current_dir[3] = dir_pos / 512L 1010 dir_pos = mlist[index][2] + 1 1011 entry_pos = dir_pos + 16L 1012 1013 tmp_avail = current_dir[0] 1014 tmp_used = current_dir[3] 1015 1016 #Fill directory with empty data until directory is full 1017 while current_dir[0] > 0: 1018 entry_pos = dir_pos + 16L + (16L * current_dir[3]) 1019 self._write_data(np.array([0,0,0,0]), imgf, entry_pos, np.uint32) 1020 current_dir[0] = current_dir[0] - 1 1021 current_dir[3] = current_dir[3] + 1 1022 1023 current_dir[0] = tmp_avail 1024 current_dir[3] = tmp_used 1025 1026 #Write directory index 1027 self._write_data(current_dir, imgf, dir_pos, endianness='>') 1028 1029 1030 @classmethod 1031 def from_image(klass, img): 1032 raise NotImplementedError("Ecat images can only be generated "\ 1033 "from file objects") 1034 1035 @classmethod 1036 def load(klass, filespec): 1037 return klass.from_filename(filespec) 1038 1039 1040load = EcatImage.load 1041