1""" Classes to wrap DICOM objects and files 2 3The wrappers encapsulate the capabilities of the different DICOM 4formats. 5 6They also allow dictionary-like access to named fields. 7 8For calculated attributes, we return None where needed data is missing. 9It seemed strange to raise an error during attribute processing, other 10than an AttributeError - breaking the 'properties manifesto'. So, any 11processing that needs to raise an error, should be in a method, rather 12than in a property, or property-like thing. 13""" 14 15import operator 16import warnings 17 18import numpy as np 19 20from . import csareader as csar 21from .dwiparams import B2q, nearest_pos_semi_def, q2bg 22from ..openers import ImageOpener 23from ..onetime import auto_attr as one_time 24from ..pydicom_compat import tag_for_keyword, Sequence 25from ..deprecated import deprecate_with_version 26 27 28class WrapperError(Exception): 29 pass 30 31 32class WrapperPrecisionError(WrapperError): 33 pass 34 35 36def wrapper_from_file(file_like, *args, **kwargs): 37 r""" Create DICOM wrapper from `file_like` object 38 39 Parameters 40 ---------- 41 file_like : object 42 filename string or file-like object, pointing to a valid DICOM 43 file readable by ``pydicom`` 44 \*args : positional 45 args to ``dicom.read_file`` command. 46 \*\*kwargs : keyword 47 args to ``dicom.read_file`` command. ``force=True`` might be a 48 likely keyword argument. 49 50 Returns 51 ------- 52 dcm_w : ``dicomwrappers.Wrapper`` or subclass 53 DICOM wrapper corresponding to DICOM data type 54 """ 55 from ..pydicom_compat import read_file 56 57 with ImageOpener(file_like) as fobj: 58 dcm_data = read_file(fobj, *args, **kwargs) 59 return wrapper_from_data(dcm_data) 60 61 62def wrapper_from_data(dcm_data): 63 """ Create DICOM wrapper from DICOM data object 64 65 Parameters 66 ---------- 67 dcm_data : ``dicom.dataset.Dataset`` instance or similar 68 Object allowing attribute access, with DICOM attributes. 69 Probably a dataset as read by ``pydicom``. 70 71 Returns 72 ------- 73 dcm_w : ``dicomwrappers.Wrapper`` or subclass 74 DICOM wrapper corresponding to DICOM data type 75 """ 76 sop_class = dcm_data.get('SOPClassUID') 77 # try to detect what type of dicom object to wrap 78 if sop_class == '1.2.840.10008.5.1.4.1.1.4.1': # Enhanced MR Image Storage 79 # currently only Philips is using Enhanced Multiframe DICOM 80 return MultiframeWrapper(dcm_data) 81 # Check for Siemens DICOM format types 82 # Only Siemens will have data for the CSA header 83 try: 84 csa = csar.get_csa_header(dcm_data) 85 except csar.CSAReadError as e: 86 warnings.warn('Error while attempting to read CSA header: ' + 87 str(e.args) + 88 '\n Ignoring Siemens private (CSA) header info.') 89 csa = None 90 if csa is None: 91 return Wrapper(dcm_data) 92 if csar.is_mosaic(csa): 93 # Mosaic is a "tiled" image 94 return MosaicWrapper(dcm_data, csa) 95 # Assume data is in a single slice format per file 96 return SiemensWrapper(dcm_data, csa) 97 98 99class Wrapper(object): 100 """ Class to wrap general DICOM files 101 102 Methods: 103 104 * get_affine() (deprecated, use affine property instead) 105 * get_data() 106 * get_pixel_array() 107 * is_same_series(other) 108 * __getitem__ : return attributes from `dcm_data` 109 * get(key[, default]) - as usual given __getitem__ above 110 111 Attributes and things that look like attributes: 112 113 * affine : (4, 4) array 114 * dcm_data : object 115 * image_shape : tuple 116 * image_orient_patient : (3,2) array 117 * slice_normal : (3,) array 118 * rotation_matrix : (3,3) array 119 * voxel_sizes : tuple length 3 120 * image_position : sequence length 3 121 * slice_indicator : float 122 * series_signature : tuple 123 """ 124 is_csa = False 125 is_mosaic = False 126 is_multiframe = False 127 b_matrix = None 128 q_vector = None 129 b_value = None 130 b_vector = None 131 132 def __init__(self, dcm_data): 133 """ Initialize wrapper 134 135 Parameters 136 ---------- 137 dcm_data : object 138 object should allow 'get' and '__getitem__' access. Usually this 139 will be a ``dicom.dataset.Dataset`` object resulting from reading a 140 DICOM file, but a dictionary should also work. 141 """ 142 self.dcm_data = dcm_data 143 144 @one_time 145 def image_shape(self): 146 """ The array shape as it will be returned by ``get_data()`` 147 """ 148 shape = (self.get('Rows'), self.get('Columns')) 149 if None in shape: 150 return None 151 return shape 152 153 @one_time 154 def image_orient_patient(self): 155 """ Note that this is _not_ LR flipped """ 156 iop = self.get('ImageOrientationPatient') 157 if iop is None: 158 return None 159 # Values are python Decimals in pydicom 0.9.7 160 iop = np.array(list(map(float, iop))) 161 return np.array(iop).reshape(2, 3).T 162 163 @one_time 164 def slice_normal(self): 165 iop = self.image_orient_patient 166 if iop is None: 167 return None 168 # iop[:, 0] is column index cosine, iop[:, 1] is row index cosine 169 return np.cross(iop[:, 1], iop[:, 0]) 170 171 @one_time 172 def rotation_matrix(self): 173 """ Return rotation matrix between array indices and mm 174 175 Note that we swap the two columns of the 'ImageOrientPatient' 176 when we create the rotation matrix. This is takes into account 177 the slightly odd ij transpose construction of the DICOM 178 orientation fields - see doc/theory/dicom_orientaiton.rst. 179 """ 180 iop = self.image_orient_patient 181 s_norm = self.slice_normal 182 if iop is None or s_norm is None: 183 return None 184 R = np.eye(3) 185 # np.fliplr(iop) gives matrix F in 186 # doc/theory/dicom_orientation.rst The fliplr accounts for the 187 # fact that the first column in ``iop`` refers to changes in 188 # column index, and the second to changes in row index. 189 R[:, :2] = np.fliplr(iop) 190 R[:, 2] = s_norm 191 # check this is in fact a rotation matrix. Error comes from compromise 192 # motivated in ``doc/source/notebooks/ata_error.ipynb``, and from 193 # discussion at https://github.com/nipy/nibabel/pull/156 194 if not np.allclose(np.eye(3), np.dot(R, R.T), atol=5e-5): 195 raise WrapperPrecisionError('Rotation matrix not nearly ' 196 'orthogonal') 197 return R 198 199 @one_time 200 def voxel_sizes(self): 201 """ voxel sizes for array as returned by ``get_data()`` 202 """ 203 # pix space gives (row_spacing, column_spacing). That is, the 204 # mm you move when moving from one row to the next, and the mm 205 # you move when moving from one column to the next 206 pix_space = self.get('PixelSpacing') 207 if pix_space is None: 208 return None 209 zs = self.get('SpacingBetweenSlices') 210 if zs is None: 211 zs = self.get('SliceThickness') 212 if zs is None or zs == '': 213 zs = 1 214 # Protect from python decimals in pydicom 0.9.7 215 zs = float(zs) 216 pix_space = list(map(float, pix_space)) 217 return tuple(pix_space + [zs]) 218 219 @one_time 220 def image_position(self): 221 """ Return position of first voxel in data block 222 223 Parameters 224 ---------- 225 None 226 227 Returns 228 ------- 229 img_pos : (3,) array 230 position in mm of voxel (0,0) in image array 231 """ 232 ipp = self.get('ImagePositionPatient') 233 if ipp is None: 234 return None 235 # Values are python Decimals in pydicom 0.9.7 236 return np.array(list(map(float, ipp))) 237 238 @one_time 239 def slice_indicator(self): 240 """ A number that is higher for higher slices in Z 241 242 Comparing this number between two adjacent slices should give a 243 difference equal to the voxel size in Z. 244 245 See doc/theory/dicom_orientation for description 246 """ 247 ipp = self.image_position 248 s_norm = self.slice_normal 249 if ipp is None or s_norm is None: 250 return None 251 return np.inner(ipp, s_norm) 252 253 @one_time 254 def instance_number(self): 255 """ Just because we use this a lot for sorting """ 256 return self.get('InstanceNumber') 257 258 @one_time 259 def series_signature(self): 260 """ Signature for matching slices into series 261 262 We use `signature` in ``self.is_same_series(other)``. 263 264 Returns 265 ------- 266 signature : dict 267 with values of 2-element sequences, where first element is 268 value, and second element is function to compare this value 269 with another. This allows us to pass things like arrays, 270 that might need to be ``allclose`` instead of equal 271 """ 272 # dictionary with value, comparison func tuple 273 signature = {} 274 eq = operator.eq 275 for key in ('SeriesInstanceUID', 276 'SeriesNumber', 277 'ImageType', 278 'SequenceName', 279 'EchoNumbers'): 280 signature[key] = (self.get(key), eq) 281 signature['image_shape'] = (self.image_shape, eq) 282 signature['iop'] = (self.image_orient_patient, none_or_close) 283 signature['vox'] = (self.voxel_sizes, none_or_close) 284 return signature 285 286 def __getitem__(self, key): 287 """ Return values from DICOM object""" 288 if key not in self.dcm_data: 289 raise KeyError(f'"{key}" not in self.dcm_data') 290 return self.dcm_data.get(key) 291 292 def get(self, key, default=None): 293 """ Get values from underlying dicom data """ 294 return self.dcm_data.get(key, default) 295 296 @deprecate_with_version('get_affine method is deprecated.\n' 297 'Please use the ``img.affine`` property ' 298 'instead.', 299 '2.5.1', '4.0') 300 def get_affine(self): 301 return self.affine 302 303 @property 304 def affine(self): 305 """ Mapping between voxel and DICOM coordinate system 306 307 (4, 4) affine matrix giving transformation between voxels in data array 308 and mm in the DICOM patient coordinate system. 309 """ 310 # rotation matrix already accounts for the ij transpose in the 311 # DICOM image orientation patient transform. So. column 0 is 312 # direction cosine for changes in row index, column 1 is 313 # direction cosine for changes in column index 314 orient = self.rotation_matrix 315 # therefore, these voxel sizes are in the right order (row, 316 # column, slice) 317 vox = self.voxel_sizes 318 ipp = self.image_position 319 if any(p is None for p in (orient, vox, ipp)): 320 raise WrapperError('Not enough information for affine') 321 aff = np.eye(4) 322 aff[:3, :3] = orient * np.array(vox) 323 aff[:3, 3] = ipp 324 return aff 325 326 def get_pixel_array(self): 327 """ Return unscaled pixel array from DICOM """ 328 data = self.dcm_data.get('pixel_array') 329 if data is None: 330 raise WrapperError('Cannot find data in DICOM') 331 return data 332 333 def get_data(self): 334 """ Get scaled image data from DICOMs 335 336 We return the data as DICOM understands it, first dimension is 337 rows, second dimension is columns 338 339 Returns 340 ------- 341 data : array 342 array with data as scaled from any scaling in the DICOM 343 fields. 344 """ 345 return self._scale_data(self.get_pixel_array()) 346 347 def is_same_series(self, other): 348 """ Return True if `other` appears to be in same series 349 350 Parameters 351 ---------- 352 other : object 353 object with ``series_signature`` attribute that is a 354 mapping. Usually it's a ``Wrapper`` or sub-class instance. 355 356 Returns 357 ------- 358 tf : bool 359 True if `other` might be in the same series as `self`, False 360 otherwise. 361 """ 362 # compare signature dictionaries. The dictionaries each contain 363 # comparison rules, we prefer our own when we have them. If a 364 # key is not present in either dictionary, assume the value is 365 # None. 366 my_sig = self.series_signature 367 your_sig = other.series_signature 368 my_keys = set(my_sig) 369 your_keys = set(your_sig) 370 # we have values in both signatures 371 for key in my_keys.intersection(your_keys): 372 v1, func = my_sig[key] 373 v2, _ = your_sig[key] 374 if not func(v1, v2): 375 return False 376 # values present in one or the other but not both 377 for keys, sig in ((my_keys - your_keys, my_sig), 378 (your_keys - my_keys, your_sig)): 379 for key in keys: 380 v1, func = sig[key] 381 if not func(v1, None): 382 return False 383 return True 384 385 def _scale_data(self, data): 386 # depending on pydicom and dicom files, values might need casting from 387 # Decimal to float 388 scale = float(self.get('RescaleSlope', 1)) 389 offset = float(self.get('RescaleIntercept', 0)) 390 return self._apply_scale_offset(data, scale, offset) 391 392 def _apply_scale_offset(self, data, scale, offset): 393 # a little optimization. If we are applying either the scale or 394 # the offset, we need to allow upcasting to float. 395 if scale != 1: 396 if offset == 0: 397 return data * scale 398 return data * scale + offset 399 if offset != 0: 400 return data + offset 401 return data 402 403 @one_time 404 def b_value(self): 405 """ Return b value for diffusion or None if not available 406 """ 407 q_vec = self.q_vector 408 if q_vec is None: 409 return None 410 return q2bg(q_vec)[0] 411 412 @one_time 413 def b_vector(self): 414 """ Return b vector for diffusion or None if not available 415 """ 416 q_vec = self.q_vector 417 if q_vec is None: 418 return None 419 return q2bg(q_vec)[1] 420 421 422class MultiframeWrapper(Wrapper): 423 """Wrapper for Enhanced MR Storage SOP Class 424 425 Tested with Philips' Enhanced DICOM implementation. 426 427 The specification for the Enhanced MR image IOP / SOP began life as `DICOM 428 supplement 49 <ftp://medical.nema.org/medical/dicom/final/sup49_ft.pdf>`_, 429 but as of 2016 it is part of the standard. In particular see: 430 431 * `A.36 Enhanced MR Information Object Definitions 432 <http://dicom.nema.org/medical/dicom/current/output/pdf/part03.pdf#sect_A.36>`_; 433 * `C.7.6.16 Multi-Frame Functional Groups Module 434 <http://dicom.nema.org/medical/dicom/current/output/pdf/part03.pdf#sect_C.7.6.16>`_; 435 * `C.7.6.17 Multi-Frame Dimension Module 436 <http://dicom.nema.org/medical/dicom/current/output/pdf/part03.pdf#sect_C.7.6.17>`_. 437 438 Attributes 439 ---------- 440 is_multiframe : boolean 441 Identifies `dcmdata` as multi-frame 442 frames : sequence 443 A sequence of ``dicom.dataset.Dataset`` objects populated by the 444 ``dicom.dataset.Dataset.PerFrameFunctionalGroupsSequence`` attribute 445 shared : object 446 The first (and only) ``dicom.dataset.Dataset`` object from a 447 ``dicom.dataset.Dataset.SharedFunctionalgroupSequence``. 448 449 Methods 450 ------- 451 image_shape(self) 452 image_orient_patient(self) 453 voxel_sizes(self) 454 image_position(self) 455 series_signature(self) 456 get_data(self) 457 """ 458 is_multiframe = True 459 460 def __init__(self, dcm_data): 461 """Initializes MultiframeWrapper 462 463 Parameters 464 ---------- 465 dcm_data : object 466 object should allow 'get' and '__getitem__' access. Usually this 467 will be a ``dicom.dataset.Dataset`` object resulting from reading a 468 DICOM file, but a dictionary should also work. 469 """ 470 Wrapper.__init__(self, dcm_data) 471 self.dcm_data = dcm_data 472 self.frames = dcm_data.get('PerFrameFunctionalGroupsSequence') 473 try: 474 self.frames[0] 475 except TypeError: 476 raise WrapperError("PerFrameFunctionalGroupsSequence is empty.") 477 try: 478 self.shared = dcm_data.get('SharedFunctionalGroupsSequence')[0] 479 except TypeError: 480 raise WrapperError("SharedFunctionalGroupsSequence is empty.") 481 self._shape = None 482 483 @one_time 484 def image_shape(self): 485 """The array shape as it will be returned by ``get_data()`` 486 487 The shape is determined by the *Rows* DICOM attribute, *Columns* 488 DICOM attribute, and the set of frame indices given by the 489 *FrameContentSequence[0].DimensionIndexValues* DICOM attribute of each 490 element in the *PerFrameFunctionalGroupsSequence*. The first two 491 axes of the returned shape correspond to the rows, and columns 492 respectively. The remaining axes correspond to those of the frame 493 indices with order preserved. 494 495 What each axis in the frame indices refers to is given by the 496 corresponding entry in the *DimensionIndexSequence* DICOM attribute. 497 **WARNING**: Any axis refering to the *StackID* DICOM attribute will 498 have been removed from the frame indices in determining the shape. This 499 is because only a file containing a single stack is currently allowed by 500 this wrapper. 501 502 References 503 ---------- 504 * C.7.6.16 Multi-Frame Functional Groups Module: 505 http://dicom.nema.org/medical/dicom/current/output/pdf/part03.pdf#sect_C.7.6.16 506 * C.7.6.17 Multi-Frame Dimension Module: 507 http://dicom.nema.org/medical/dicom/current/output/pdf/part03.pdf#sect_C.7.6.17 508 * Diagram of DimensionIndexSequence and DimensionIndexValues: 509 http://dicom.nema.org/medical/dicom/current/output/pdf/part03.pdf#figure_C.7.6.17-1 510 """ 511 rows, cols = self.get('Rows'), self.get('Columns') 512 if None in (rows, cols): 513 raise WrapperError("Rows and/or Columns are empty.") 514 515 # Check number of frames 516 first_frame = self.frames[0] 517 n_frames = self.get('NumberOfFrames') 518 # some Philips may have derived images appended 519 has_derived = False 520 if hasattr(first_frame, 'get') and first_frame.get([0x18, 0x9117]): 521 # DWI image may include derived isotropic, ADC or trace volume 522 try: 523 self.frames = Sequence( 524 frame for frame in self.frames if 525 frame.MRDiffusionSequence[0].DiffusionDirectionality 526 != 'ISOTROPIC' 527 ) 528 except IndexError: 529 # Sequence tag is found but missing items! 530 raise WrapperError("Diffusion file missing information") 531 except AttributeError: 532 # DiffusionDirectionality tag is not required 533 pass 534 else: 535 if n_frames != len(self.frames): 536 warnings.warn("Derived images found and removed") 537 n_frames = len(self.frames) 538 has_derived = True 539 540 assert len(self.frames) == n_frames 541 frame_indices = np.array( 542 [frame.FrameContentSequence[0].DimensionIndexValues 543 for frame in self.frames]) 544 # Check that there is only one multiframe stack index 545 stack_ids = set(frame.FrameContentSequence[0].StackID 546 for frame in self.frames) 547 if len(stack_ids) > 1: 548 raise WrapperError("File contains more than one StackID. " 549 "Cannot handle multi-stack files") 550 # Determine if one of the dimension indices refers to the stack id 551 dim_seq = [dim.DimensionIndexPointer 552 for dim in self.get('DimensionIndexSequence')] 553 stackid_tag = tag_for_keyword('StackID') 554 # remove the stack id axis if present 555 if stackid_tag in dim_seq: 556 stackid_dim_idx = dim_seq.index(stackid_tag) 557 frame_indices = np.delete(frame_indices, stackid_dim_idx, axis=1) 558 dim_seq.pop(stackid_dim_idx) 559 if has_derived: 560 # derived volume is included 561 derived_tag = tag_for_keyword("DiffusionBValue") 562 if derived_tag not in dim_seq: 563 raise WrapperError("Missing information, cannot remove indices " 564 "with confidence.") 565 derived_dim_idx = dim_seq.index(derived_tag) 566 frame_indices = np.delete(frame_indices, derived_dim_idx, axis=1) 567 # account for the 2 additional dimensions (row and column) not included 568 # in the indices 569 n_dim = frame_indices.shape[1] + 2 570 # Store frame indices 571 self._frame_indices = frame_indices 572 if n_dim < 4: # 3D volume 573 return rows, cols, n_frames 574 # More than 3 dimensions 575 ns_unique = [len(np.unique(row)) for row in self._frame_indices.T] 576 shape = (rows, cols) + tuple(ns_unique) 577 n_vols = np.prod(shape[3:]) 578 if n_frames != n_vols * shape[2]: 579 raise WrapperError("Calculated shape does not match number of " 580 "frames.") 581 return tuple(shape) 582 583 @one_time 584 def image_orient_patient(self): 585 """ 586 Note that this is _not_ LR flipped 587 """ 588 try: 589 iop = self.shared.PlaneOrientationSequence[0].ImageOrientationPatient 590 except AttributeError: 591 try: 592 iop = self.frames[0].PlaneOrientationSequence[0].ImageOrientationPatient 593 except AttributeError: 594 raise WrapperError("Not enough information for " 595 "image_orient_patient") 596 if iop is None: 597 return None 598 iop = np.array(list(map(float, iop))) 599 return np.array(iop).reshape(2, 3).T 600 601 @one_time 602 def voxel_sizes(self): 603 """ Get i, j, k voxel sizes """ 604 try: 605 pix_measures = self.shared.PixelMeasuresSequence[0] 606 except AttributeError: 607 try: 608 pix_measures = self.frames[0].PixelMeasuresSequence[0] 609 except AttributeError: 610 raise WrapperError("Not enough data for pixel spacing") 611 pix_space = pix_measures.PixelSpacing 612 try: 613 zs = pix_measures.SliceThickness 614 except AttributeError: 615 zs = self.get('SpacingBetweenSlices') 616 if zs is None: 617 raise WrapperError('Not enough data for slice thickness') 618 # Ensure values are float rather than Decimal 619 return tuple(map(float, list(pix_space) + [zs])) 620 621 @one_time 622 def image_position(self): 623 try: 624 ipp = self.shared.PlanePositionSequence[0].ImagePositionPatient 625 except AttributeError: 626 try: 627 ipp = self.frames[0].PlanePositionSequence[0].ImagePositionPatient 628 except AttributeError: 629 raise WrapperError('Cannot get image position from dicom') 630 if ipp is None: 631 return None 632 return np.array(list(map(float, ipp))) 633 634 @one_time 635 def series_signature(self): 636 signature = {} 637 eq = operator.eq 638 for key in ('SeriesInstanceUID', 639 'SeriesNumber', 640 'ImageType'): 641 signature[key] = (self.get(key), eq) 642 signature['image_shape'] = (self.image_shape, eq) 643 signature['iop'] = (self.image_orient_patient, none_or_close) 644 signature['vox'] = (self.voxel_sizes, none_or_close) 645 return signature 646 647 def get_data(self): 648 shape = self.image_shape 649 if shape is None: 650 raise WrapperError('No valid information for image shape') 651 data = self.get_pixel_array() 652 # Roll frames axis to last 653 data = data.transpose((1, 2, 0)) 654 # Sort frames with first index changing fastest, last slowest 655 sorted_indices = np.lexsort(self._frame_indices.T) 656 data = data[..., sorted_indices] 657 data = data.reshape(shape, order='F') 658 return self._scale_data(data) 659 660 def _scale_data(self, data): 661 pix_trans = getattr( 662 self.frames[0], 'PixelValueTransformationSequence', None) 663 if pix_trans is None: 664 return super(MultiframeWrapper, self)._scale_data(data) 665 scale = float(pix_trans[0].RescaleSlope) 666 offset = float(pix_trans[0].RescaleIntercept) 667 return self._apply_scale_offset(data, scale, offset) 668 669 670class SiemensWrapper(Wrapper): 671 """ Wrapper for Siemens format DICOMs 672 673 Adds attributes: 674 675 * csa_header : mapping 676 * b_matrix : (3,3) array 677 * q_vector : (3,) array 678 """ 679 is_csa = True 680 681 def __init__(self, dcm_data, csa_header=None): 682 """ Initialize Siemens wrapper 683 684 The Siemens-specific information is in the `csa_header`, either 685 passed in here, or read from the input `dcm_data`. 686 687 Parameters 688 ---------- 689 dcm_data : object 690 object should allow 'get' and '__getitem__' access. If `csa_header` 691 is None, it should also be possible to extract a CSA header from 692 `dcm_data`. Usually this will be a ``dicom.dataset.Dataset`` object 693 resulting from reading a DICOM file. A dict should also work. 694 csa_header : None or mapping, optional 695 mapping giving values for Siemens CSA image sub-header. If 696 None, we try and read the CSA information from `dcm_data`. 697 If this fails, we fall back to an empty dict. 698 """ 699 super(SiemensWrapper, self).__init__(dcm_data) 700 if dcm_data is None: 701 dcm_data = {} 702 self.dcm_data = dcm_data 703 if csa_header is None: 704 csa_header = csar.get_csa_header(dcm_data) 705 if csa_header is None: 706 csa_header = {} 707 self.csa_header = csa_header 708 709 @one_time 710 def slice_normal(self): 711 # The std_slice_normal comes from the cross product of the directions 712 # in the ImageOrientationPatient 713 std_slice_normal = super(SiemensWrapper, self).slice_normal 714 csa_slice_normal = csar.get_slice_normal(self.csa_header) 715 if std_slice_normal is None and csa_slice_normal is None: 716 return None 717 elif std_slice_normal is None: 718 return np.array(csa_slice_normal) 719 elif csa_slice_normal is None: 720 return std_slice_normal 721 else: 722 # Make sure the two normals are very close to parallel unit vectors 723 dot_prod = np.dot(csa_slice_normal, std_slice_normal) 724 assert np.allclose(np.fabs(dot_prod), 1.0, atol=1e-5) 725 # Use the slice normal computed with the cross product as it will 726 # always be the most orthogonal, but take the sign from the CSA 727 # slice normal 728 if dot_prod < 0: 729 return -std_slice_normal 730 else: 731 return std_slice_normal 732 733 @one_time 734 def series_signature(self): 735 """ Add ICE dims from CSA header to signature """ 736 signature = super(SiemensWrapper, self).series_signature 737 ice = csar.get_ice_dims(self.csa_header) 738 if ice is not None: 739 ice = ice[:6] + ice[8:9] 740 signature['ICE_Dims'] = (ice, lambda x, y: x == y) 741 return signature 742 743 @one_time 744 def b_matrix(self): 745 """ Get DWI B matrix referring to voxel space 746 747 Parameters 748 ---------- 749 None 750 751 Returns 752 ------- 753 B : (3,3) array or None 754 B matrix in *voxel* orientation space. Returns None if this is 755 not a Siemens header with the required information. We return 756 None if this is a b0 acquisition 757 """ 758 hdr = self.csa_header 759 # read B matrix as recorded in CSA header. This matrix refers to 760 # the space of the DICOM patient coordinate space. 761 B = csar.get_b_matrix(hdr) 762 if B is None: # may be not diffusion or B0 image 763 bval_requested = csar.get_b_value(hdr) 764 if bval_requested is None: 765 return None 766 if bval_requested != 0: 767 raise csar.CSAError('No B matrix and b value != 0') 768 return np.zeros((3, 3)) 769 # rotation from voxels to DICOM PCS, inverted to give the rotation 770 # from DPCS to voxels. Because this is an orthonormal matrix, its 771 # transpose is its inverse 772 R = self.rotation_matrix.T 773 # because B results from V dot V.T, the rotation B is given by R dot 774 # V dot V.T dot R.T == R dot B dot R.T 775 B_vox = np.dot(R, np.dot(B, R.T)) 776 # fix presumed rounding errors in the B matrix by making it positive 777 # semi-definite. 778 return nearest_pos_semi_def(B_vox) 779 780 @one_time 781 def q_vector(self): 782 """ Get DWI q vector referring to voxel space 783 784 Parameters 785 ---------- 786 None 787 788 Returns 789 ------- 790 q: (3,) array 791 Estimated DWI q vector in *voxel* orientation space. Returns 792 None if this is not (detectably) a DWI 793 """ 794 B = self.b_matrix 795 if B is None: 796 return None 797 # We've enforced more or less positive semi definite with the 798 # b_matrix routine 799 return B2q(B, tol=1e-8) 800 801 802class MosaicWrapper(SiemensWrapper): 803 """ Class for Siemens mosaic format data 804 805 Mosaic format is a way of storing a 3D image in a 2D slice - and 806 it's as simple as you'd imagine it would be - just storing the slices 807 in a mosaic similar to a light-box print. 808 809 We need to allow for this when getting the data and (because of an 810 idiosyncrasy in the way Siemens stores the images) calculating the 811 position of the first voxel. 812 813 Adds attributes: 814 815 * n_mosaic : int 816 * mosaic_size : int 817 """ 818 is_mosaic = True 819 820 def __init__(self, dcm_data, csa_header=None, n_mosaic=None): 821 """ Initialize Siemens Mosaic wrapper 822 823 The Siemens-specific information is in the `csa_header`, either 824 passed in here, or read from the input `dcm_data`. 825 826 Parameters 827 ---------- 828 dcm_data : object 829 object should allow 'get' and '__getitem__' access. If `csa_header` 830 is None, it should also be possible for to extract a CSA header from 831 `dcm_data`. Usually this will be a ``dicom.dataset.Dataset`` object 832 resulting from reading a DICOM file. A dict should also work. 833 csa_header : None or mapping, optional 834 mapping giving values for Siemens CSA image sub-header. 835 n_mosaic : None or int, optional 836 number of images in mosaic. If None, try to get this number 837 from `csa_header`. If this fails, raise an error 838 """ 839 SiemensWrapper.__init__(self, dcm_data, csa_header) 840 if n_mosaic is None: 841 try: 842 n_mosaic = csar.get_n_mosaic(self.csa_header) 843 except KeyError: 844 pass 845 if n_mosaic is None or n_mosaic == 0: 846 raise WrapperError('No valid mosaic number in CSA ' 847 'header; is this really ' 848 'Siemens mosiac data?') 849 self.n_mosaic = n_mosaic 850 self.mosaic_size = int(np.ceil(np.sqrt(n_mosaic))) 851 852 @one_time 853 def image_shape(self): 854 """ Return image shape as returned by ``get_data()`` """ 855 # reshape pixel slice array back from mosaic 856 rows = self.get('Rows') 857 cols = self.get('Columns') 858 if None in (rows, cols): 859 return None 860 mosaic_size = self.mosaic_size 861 return (int(rows / mosaic_size), 862 int(cols / mosaic_size), 863 self.n_mosaic) 864 865 @one_time 866 def image_position(self): 867 """ Return position of first voxel in data block 868 869 Adjusts Siemens mosaic position vector for bug in mosaic format 870 position. See ``dicom_mosaic`` in doc/theory for details. 871 872 Parameters 873 ---------- 874 None 875 876 Returns 877 ------- 878 img_pos : (3,) array 879 position in mm of voxel (0,0,0) in Mosaic array 880 """ 881 ipp = super(MosaicWrapper, self).image_position 882 # mosaic image size 883 md_rows, md_cols = (self.get('Rows'), self.get('Columns')) 884 iop = self.image_orient_patient 885 pix_spacing = self.get('PixelSpacing') 886 if any(x is None for x in (ipp, md_rows, md_cols, iop, pix_spacing)): 887 return None 888 # PixelSpacing values are python Decimal in pydicom 0.9.7 889 pix_spacing = np.array(list(map(float, pix_spacing))) 890 # size of mosaic array before rearranging to 3D. 891 md_rc = np.array([md_rows, md_cols]) 892 # size of slice array after reshaping to 3D 893 rd_rc = md_rc / self.mosaic_size 894 # apply algorithm for undoing mosaic translation error - see 895 # ``dicom_mosaic`` doc 896 vox_trans_fixes = (md_rc - rd_rc) / 2 897 # flip IOP field to refer to rows then columns index change - 898 # see dicom_orientation doc 899 Q = np.fliplr(iop) * pix_spacing 900 return ipp + np.dot(Q, vox_trans_fixes[:, None]).ravel() 901 902 def get_data(self): 903 """ Get scaled image data from DICOMs 904 905 Resorts data block from mosaic to 3D 906 907 Returns 908 ------- 909 data : array 910 array with data as scaled from any scaling in the DICOM 911 fields. 912 913 Notes 914 ----- 915 The apparent image in the DICOM file is a 2D array that consists of 916 blocks, that are the output 2D slices. Let's call the original array 917 the *slab*, and the contained slices *slices*. The slices are of 918 pixel dimension ``n_slice_rows`` x ``n_slice_cols``. The slab is of 919 pixel dimension ``n_slab_rows`` x ``n_slab_cols``. Because the 920 arrangement of blocks in the slab is defined as being square, the 921 number of blocks per slab row and slab column is the same. Let 922 ``n_blocks`` be the number of blocks contained in the slab. There is 923 also ``n_slices`` - the number of slices actually collected, some 924 number <= ``n_blocks``. We have the value ``n_slices`` from the 925 'NumberOfImagesInMosaic' field of the Siemens private (CSA) header. 926 ``n_row_blocks`` and ``n_col_blocks`` are therefore given by 927 ``ceil(sqrt(n_slices))``, and ``n_blocks`` is ``n_row_blocks ** 2``. 928 Also ``n_slice_rows == n_slab_rows / n_row_blocks``, etc. Using these 929 numbers we can therefore reconstruct the slices from the 2D DICOM pixel 930 array. 931 """ 932 shape = self.image_shape 933 if shape is None: 934 raise WrapperError('No valid information for image shape') 935 n_slice_rows, n_slice_cols, n_mosaic = shape 936 n_slab_rows = self.mosaic_size 937 n_blocks = n_slab_rows ** 2 938 data = self.get_pixel_array() 939 v4 = data.reshape(n_slab_rows, n_slice_rows, 940 n_slab_rows, n_slice_cols) 941 # move the mosaic dims to the end 942 v4 = v4.transpose((1, 3, 0, 2)) 943 # pool mosaic-generated dims 944 v3 = v4.reshape((n_slice_rows, n_slice_cols, n_blocks)) 945 # delete any padding slices 946 v3 = v3[..., :n_mosaic] 947 return self._scale_data(v3) 948 949 950def none_or_close(val1, val2, rtol=1e-5, atol=1e-6): 951 """ Match if `val1` and `val2` are both None, or are close 952 953 Parameters 954 ---------- 955 val1 : None or array-like 956 val2 : None or array-like 957 rtol : float, optional 958 Relative tolerance; see ``np.allclose`` 959 atol : float, optional 960 Absolute tolerance; see ``np.allclose`` 961 962 Returns 963 ------- 964 tf : bool 965 True iff (both `val1` and `val2` are None) or (`val1` and `val2` 966 are close arrays, as detected by ``np.allclose`` with parameters 967 `rtol` and `atal`). 968 969 Examples 970 -------- 971 >>> none_or_close(None, None) 972 True 973 >>> none_or_close(1, None) 974 False 975 >>> none_or_close(None, 1) 976 False 977 >>> none_or_close([1,2], [1,2]) 978 True 979 >>> none_or_close([0,1], [0,2]) 980 False 981 """ 982 if val1 is None and val2 is None: 983 return True 984 if val1 is None or val2 is None: 985 return False 986 return np.allclose(val1, val2, rtol, atol) 987