diff --git a/Changelog b/Changelog index 689295125..24e89095f 100644 --- a/Changelog +++ b/Changelog @@ -25,6 +25,33 @@ Eric Larson (EL), Demian Wassermann, Stephan Gerhard and Ross Markello (RM). References like "pr/298" refer to github pull request numbers. +Upcoming release (To be determined) +=================================== + +New features +------------ + +Enhancements +------------ + * Ability to read data from many multiframe DICOM files that previously generated errors + +Bug fixes +--------- + * Fixed multiframe DICOM issue where data could be flipped along slice dimension relative to the + affine + * Fixed multiframe DICOM issue where ``image_position`` and the translation component in the + ``affine`` could be incorrect + +Documentation +------------- + +Maintenance +----------- + +API changes and deprecations +---------------------------- + + 5.2.1 (Monday 26 February 2024) =============================== diff --git a/nibabel/nicom/dicomwrappers.py b/nibabel/nicom/dicomwrappers.py index 2270ed3f0..374387870 100755 --- a/nibabel/nicom/dicomwrappers.py +++ b/nibabel/nicom/dicomwrappers.py @@ -467,6 +467,25 @@ def __init__(self, dcm_data): self.shared = dcm_data.get('SharedFunctionalGroupsSequence')[0] except TypeError: raise WrapperError('SharedFunctionalGroupsSequence is empty.') + # Try to determine slice order and minimal image position patient + self._frame_slc_ord = self._ipp = None + try: + frame_ipps = [f.PlanePositionSequence[0].ImagePositionPatient for f in self.frames] + except AttributeError: + try: + frame_ipps = [self.shared.PlanePositionSequence[0].ImagePositionPatient] + except AttributeError: + frame_ipps = None + if frame_ipps is not None and all(ipp is not None for ipp in frame_ipps): + frame_ipps = [np.array(list(map(float, ipp))) for ipp in frame_ipps] + frame_slc_pos = [np.inner(ipp, self.slice_normal) for ipp in frame_ipps] + rnd_slc_pos = np.round(frame_slc_pos, 4) + uniq_slc_pos = np.unique(rnd_slc_pos) + pos_ord_map = { + val: order for val, order in zip(uniq_slc_pos, np.argsort(uniq_slc_pos)) + } + self._frame_slc_ord = [pos_ord_map[pos] for pos in rnd_slc_pos] + self._ipp = frame_ipps[np.argmin(frame_slc_pos)] self._shape = None @cached_property @@ -509,14 +528,16 @@ def image_shape(self): if hasattr(first_frame, 'get') and first_frame.get([0x18, 0x9117]): # DWI image may include derived isotropic, ADC or trace volume try: - anisotropic = pydicom.Sequence( - frame - for frame in self.frames - if frame.MRDiffusionSequence[0].DiffusionDirectionality != 'ISOTROPIC' - ) + aniso_frames = pydicom.Sequence() + aniso_slc_ord = [] + for slc_ord, frame in zip(self._frame_slc_ord, self.frames): + if frame.MRDiffusionSequence[0].DiffusionDirectionality != 'ISOTROPIC': + aniso_frames.append(frame) + aniso_slc_ord.append(slc_ord) # Image contains DWI volumes followed by derived images; remove derived images - if len(anisotropic) != 0: - self.frames = anisotropic + if len(aniso_frames) != 0: + self.frames = aniso_frames + self._frame_slc_ord = aniso_slc_ord except IndexError: # Sequence tag is found but missing items! raise WrapperError('Diffusion file missing information') @@ -554,23 +575,85 @@ def image_shape(self): raise WrapperError('Missing information, cannot remove indices with confidence.') derived_dim_idx = dim_seq.index(derived_tag) frame_indices = np.delete(frame_indices, derived_dim_idx, axis=1) - # account for the 2 additional dimensions (row and column) not included - # in the indices - n_dim = frame_indices.shape[1] + 2 + dim_seq.pop(derived_dim_idx) + # Determine the shape and which indices to use + shape = [rows, cols] + curr_parts = n_frames + frames_per_part = 1 + del_indices = {} + stackpos_tag = pydicom.datadict.tag_for_keyword('InStackPositionNumber') + slice_dim_idx = dim_seq.index(stackpos_tag) + for row_idx, row in enumerate(frame_indices.T): + unique = np.unique(row) + count = len(unique) + if curr_parts == 1 or (count == 1 and row_idx != slice_dim_idx): + del_indices[row_idx] = count + continue + # Replace slice indices with order determined from slice positions along normal + if row_idx == slice_dim_idx: + if len(shape) > 2: + raise WrapperError('Non-singular index precedes the slice index') + row = self._frame_slc_ord + frame_indices.T[row_idx, :] = row + unique = np.unique(row) + if len(unique) != count: + raise WrapperError("Number of slice indices and positions don't match") + elif count == n_frames: + if shape[-1] == 'remaining': + raise WrapperError('At most one index have ambiguous size') + shape.append('remaining') + continue + new_parts, leftover = divmod(curr_parts, count) + expected = new_parts * frames_per_part + if leftover != 0 or any(np.count_nonzero(row == val) != expected for val in unique): + if row_idx == slice_dim_idx: + raise WrapperError('Missing slices from multiframe') + del_indices[row_idx] = count + continue + if shape[-1] == 'remaining': + shape[-1] = new_parts + frames_per_part *= shape[-1] + new_parts = 1 + frames_per_part *= count + shape.append(count) + curr_parts = new_parts + if shape[-1] == 'remaining': + if curr_parts > 1: + shape[-1] = curr_parts + curr_parts = 1 + else: + del_indices[len(shape)] = 1 + shape = shape[:-1] + if del_indices: + if curr_parts > 1: + ns_failed = [k for k, v in del_indices.items() if v != 1] + if len(ns_failed) > 1: + # If some indices weren't used yet but we still have unaccounted for + # partitions, try combining indices into single tuple and using that + tup_dtype = np.dtype(','.join(['I'] * len(ns_failed))) + row = [tuple(x for x in vals) for vals in frame_indices[:, ns_failed]] + row = np.array(row, dtype=tup_dtype) + frame_indices = np.delete(frame_indices, np.array(list(del_indices.keys())), axis=1) + if curr_parts > 1 and len(ns_failed) > 1: + unique = np.unique(row, axis=0) + count = len(unique) + new_parts, rem = divmod(curr_parts, count) + allowed_val_counts = [new_parts * frames_per_part, n_frames] + if rem == 0 and all( + np.count_nonzero(row == val) in allowed_val_counts for val in unique + ): + shape.append(count) + curr_parts = new_parts + ord_vals = np.argsort(unique) + order = {tuple(unique[i]): ord_vals[i] for i in range(count)} + ord_row = np.array([order[tuple(v)] for v in row]) + frame_indices = np.hstack( + [frame_indices, np.array(ord_row).reshape((n_frames, 1))] + ) + if curr_parts > 1: + raise WrapperError('Unable to determine sorting of final dimension(s)') # Store frame indices self._frame_indices = frame_indices - if n_dim < 4: # 3D volume - return rows, cols, n_frames - # More than 3 dimensions - ns_unique = [len(np.unique(row)) for row in self._frame_indices.T] - shape = (rows, cols) + tuple(ns_unique) - n_vols = np.prod(shape[3:]) - n_frames_calc = n_vols * shape[2] - if n_frames != n_frames_calc: - raise WrapperError( - f'Calculated # of frames ({n_frames_calc}={n_vols}*{shape[2]}) ' - f'of shape {shape} does not match NumberOfFrames {n_frames}.' - ) return tuple(shape) @cached_property @@ -610,18 +693,11 @@ def voxel_sizes(self): # Ensure values are float rather than Decimal return tuple(map(float, list(pix_space) + [zs])) - @cached_property + @property def image_position(self): - try: - ipp = self.shared.PlanePositionSequence[0].ImagePositionPatient - except AttributeError: - try: - ipp = self.frames[0].PlanePositionSequence[0].ImagePositionPatient - except AttributeError: - raise WrapperError('Cannot get image position from dicom') - if ipp is None: - return None - return np.array(list(map(float, ipp))) + if self._ipp is None: + raise WrapperError('Not enough information for image_position_patient') + return self._ipp @cached_property def series_signature(self): @@ -640,10 +716,11 @@ def get_data(self): raise WrapperError('No valid information for image shape') data = self.get_pixel_array() # Roll frames axis to last - data = data.transpose((1, 2, 0)) - # Sort frames with first index changing fastest, last slowest - sorted_indices = np.lexsort(self._frame_indices.T) - data = data[..., sorted_indices] + if len(data.shape) > 2: + data = data.transpose((1, 2, 0)) + # Sort frames with first index changing fastest, last slowest + sorted_indices = np.lexsort(self._frame_indices.T) + data = data[..., sorted_indices] data = data.reshape(shape, order='F') return self._scale_data(data) diff --git a/nibabel/nicom/tests/test_dicomwrappers.py b/nibabel/nicom/tests/test_dicomwrappers.py index d14c35dcd..e01759c86 100755 --- a/nibabel/nicom/tests/test_dicomwrappers.py +++ b/nibabel/nicom/tests/test_dicomwrappers.py @@ -364,7 +364,7 @@ def test_decimal_rescale(): assert dw.get_data().dtype != np.dtype(object) -def fake_frames(seq_name, field_name, value_seq): +def fake_frames(seq_name, field_name, value_seq, frame_seq=None): """Make fake frames for multiframe testing Parameters @@ -375,6 +375,8 @@ def fake_frames(seq_name, field_name, value_seq): name of field within sequence value_seq : length N sequence sequence of values + frame_seq : length N list + previous result from this function to update Returns ------- @@ -386,19 +388,28 @@ def fake_frames(seq_name, field_name, value_seq): class Fake: pass - frames = [] - for value in value_seq: - fake_frame = Fake() + if frame_seq is None: + frame_seq = [Fake() for _ in range(len(value_seq))] + for value, fake_frame in zip(value_seq, frame_seq): fake_element = Fake() setattr(fake_element, field_name, value) setattr(fake_frame, seq_name, [fake_element]) - frames.append(fake_frame) - return frames + return frame_seq -def fake_shape_dependents(div_seq, sid_seq=None, sid_dim=None): +def fake_shape_dependents( + div_seq, + sid_seq=None, + sid_dim=None, + ipp_seq=None, + slice_dim=None, + flip_ipp_idx_corr=False, +): """Make a fake dictionary of data that ``image_shape`` is dependent on. + If you are providing the ``ipp_seq`` argument, they should be generated using + a slice normal aligned with the z-axis (i.e. iop == (0, 1, 0, 1, 0, 0)). + Parameters ---------- div_seq : list of tuples @@ -407,39 +418,88 @@ def fake_shape_dependents(div_seq, sid_seq=None, sid_dim=None): list of values to use for the `StackID` of each frame. sid_dim : int the index of the column in 'div_seq' to use as 'sid_seq' + ipp_seq : list of tuples + list of values to use for `ImagePositionPatient` for each frame + slice_dim : int + the index of the column in 'div_seq' corresponding to slices + flip_ipp_idx_corr : bool + generate ipp values so slice location is negatively correlated with slice index """ - class DimIdxSeqElem: + class PrintBase: + def __repr__(self): + attr_strs = [] + for attr in dir(self): + if attr[0].isupper(): + attr_strs.append(f'{attr}={getattr(self, attr)}') + return f"{self.__class__.__name__}({', '.join(attr_strs)})" + + class DimIdxSeqElem(PrintBase): def __init__(self, dip=(0, 0), fgp=None): self.DimensionIndexPointer = dip if fgp is not None: self.FunctionalGroupPointer = fgp - class FrmContSeqElem: + class FrmContSeqElem(PrintBase): def __init__(self, div, sid): self.DimensionIndexValues = div self.StackID = sid - class PerFrmFuncGrpSeqElem: - def __init__(self, div, sid): + class PlnPosSeqElem(PrintBase): + def __init__(self, ipp): + self.ImagePositionPatient = ipp + + class PlnOrientSeqElem(PrintBase): + def __init__(self, iop): + self.ImageOrientationPatient = iop + + class PerFrmFuncGrpSeqElem(PrintBase): + def __init__(self, div, sid, ipp, iop): self.FrameContentSequence = [FrmContSeqElem(div, sid)] + self.PlanePositionSequence = [PlnPosSeqElem(ipp)] + self.PlaneOrientationSequence = [PlnOrientSeqElem(iop)] # if no StackID values passed in then use the values at index 'sid_dim' in # the value for DimensionIndexValues for it + n_indices = len(div_seq[0]) if sid_seq is None: if sid_dim is None: sid_dim = 0 sid_seq = [div[sid_dim] for div in div_seq] - # create the DimensionIndexSequence + # Determine slice_dim and create per-slice ipp information + if slice_dim is None: + slice_dim = 1 if sid_dim == 0 else 0 num_of_frames = len(div_seq) - dim_idx_seq = [DimIdxSeqElem()] * num_of_frames + frame_slc_indices = np.array(div_seq)[:, slice_dim] + uniq_slc_indices = np.unique(frame_slc_indices) + n_slices = len(uniq_slc_indices) + iop_seq = [(0.0, 1.0, 0.0, 1.0, 0.0, 0.0) for _ in range(num_of_frames)] + if ipp_seq is None: + slc_locs = np.linspace(-1.0, 1.0, n_slices) + if flip_ipp_idx_corr: + slc_locs = slc_locs[::-1] + slc_idx_loc = { + div_idx: slc_locs[arr_idx] for arr_idx, div_idx in enumerate(np.sort(uniq_slc_indices)) + } + ipp_seq = [(-1.0, -1.0, slc_idx_loc[idx]) for idx in frame_slc_indices] + else: + assert flip_ipp_idx_corr is False # caller can flip it themselves + assert len(ipp_seq) == num_of_frames + # create the DimensionIndexSequence + dim_idx_seq = [DimIdxSeqElem()] * n_indices + # Add entry for InStackPositionNumber to DimensionIndexSequence + fcs_tag = pydicom.datadict.tag_for_keyword('FrameContentSequence') + isp_tag = pydicom.datadict.tag_for_keyword('InStackPositionNumber') + dim_idx_seq[slice_dim] = DimIdxSeqElem(isp_tag, fcs_tag) # add an entry for StackID into the DimensionIndexSequence if sid_dim is not None: sid_tag = pydicom.datadict.tag_for_keyword('StackID') - fcs_tag = pydicom.datadict.tag_for_keyword('FrameContentSequence') dim_idx_seq[sid_dim] = DimIdxSeqElem(sid_tag, fcs_tag) # create the PerFrameFunctionalGroupsSequence - frames = [PerFrmFuncGrpSeqElem(div, sid) for div, sid in zip(div_seq, sid_seq)] + frames = [ + PerFrmFuncGrpSeqElem(div, sid, ipp, iop) + for div, sid, ipp, iop in zip(div_seq, sid_seq, ipp_seq, iop_seq) + ] return { 'NumberOfFrames': num_of_frames, 'DimensionIndexSequence': dim_idx_seq, @@ -480,7 +540,19 @@ def test_shape(self): # PerFrameFunctionalGroupsSequence does not match NumberOfFrames with pytest.raises(AssertionError): dw.image_shape - # check 3D shape when StackID index is 0 + # check 2D shape with StackID index is 0 + div_seq = ((1, 1),) + fake_mf.update(fake_shape_dependents(div_seq, sid_dim=0)) + assert MFW(fake_mf).image_shape == (32, 64) + # Check 2D shape with extraneous extra indices + div_seq = ((1, 1, 2),) + fake_mf.update(fake_shape_dependents(div_seq, sid_dim=0)) + assert MFW(fake_mf).image_shape == (32, 64) + # Check 2D plus time + div_seq = ((1, 1, 1), (1, 1, 2), (1, 1, 3)) + fake_mf.update(fake_shape_dependents(div_seq, sid_dim=0)) + assert MFW(fake_mf).image_shape == (32, 64, 1, 3) + # Check 3D shape when StackID index is 0 div_seq = ((1, 1), (1, 2), (1, 3), (1, 4)) fake_mf.update(fake_shape_dependents(div_seq, sid_dim=0)) assert MFW(fake_mf).image_shape == (32, 64, 4) @@ -506,6 +578,17 @@ def test_shape(self): div_seq = ((1, 1, 0), (1, 2, 0), (1, 1, 3), (1, 2, 3)) fake_mf.update(fake_shape_dependents(div_seq, sid_dim=0)) assert MFW(fake_mf).image_shape == (32, 64, 2, 2) + # Check number of IPP vals match the number of slices or we raise + frames = fake_mf['PerFrameFunctionalGroupsSequence'] + for frame in frames[1:]: + frame.PlanePositionSequence = frames[0].PlanePositionSequence[:] + with pytest.raises(didw.WrapperError): + MFW(fake_mf).image_shape + # Check we raise on missing slices + div_seq = ((1, 1, 0), (1, 2, 0), (1, 1, 1)) + fake_mf.update(fake_shape_dependents(div_seq, sid_dim=0)) + with pytest.raises(didw.WrapperError): + MFW(fake_mf).image_shape # check 3D shape when there is no StackID index div_seq = ((1,), (2,), (3,), (4,)) sid_seq = (1, 1, 1, 1) @@ -541,6 +624,68 @@ def test_shape(self): div_seq = ((1, 1, 1), (2, 1, 1), (1, 1, 2), (2, 1, 2), (1, 1, 3), (2, 1, 3)) fake_mf.update(fake_shape_dependents(div_seq, sid_dim=1)) assert MFW(fake_mf).image_shape == (32, 64, 2, 3) + # Check non-singular dimension preceding slice dim raises + div_seq = ((1, 1, 1), (1, 2, 1), (1, 1, 2), (1, 2, 2), (1, 1, 3), (1, 2, 3)) + fake_mf.update(fake_shape_dependents(div_seq, sid_dim=0, slice_dim=2)) + with pytest.raises(didw.WrapperError): + MFW(fake_mf).image_shape + # Test with combo indices, here with the last two needing to be combined into + # a single index corresponding to [(1, 1), (1, 1), (2, 1), (2, 1), (2, 2), (2, 2)] + div_seq = ( + (1, 1, 1, 1), + (1, 2, 1, 1), + (1, 1, 2, 1), + (1, 2, 2, 1), + (1, 1, 2, 2), + (1, 2, 2, 2), + ) + fake_mf.update(fake_shape_dependents(div_seq, sid_dim=0)) + assert MFW(fake_mf).image_shape == (32, 64, 2, 3) + # Test invalid 4D indices + div_seq = ((1, 1, 1), (1, 2, 1), (1, 1, 2), (1, 2, 2), (1, 1, 3), (1, 2, 4)) + fake_mf.update(fake_shape_dependents(div_seq, sid_dim=0)) + with pytest.raises(didw.WrapperError): + MFW(fake_mf).image_shape + div_seq = ((1, 1, 1), (1, 2, 1), (1, 1, 2), (1, 2, 2), (1, 1, 3), (1, 2, 2)) + fake_mf.update(fake_shape_dependents(div_seq, sid_dim=0)) + with pytest.raises(didw.WrapperError): + MFW(fake_mf).image_shape + # Time index that is unique to each frame + div_seq = ((1, 1, 1), (1, 2, 2), (1, 1, 3), (1, 2, 4), (1, 1, 5), (1, 2, 6)) + fake_mf.update(fake_shape_dependents(div_seq, sid_dim=0)) + assert MFW(fake_mf).image_shape == (32, 64, 2, 3) + div_seq = ( + (1, 1, 1, 1), + (1, 2, 2, 1), + (1, 1, 3, 1), + (1, 2, 4, 1), + (1, 1, 5, 1), + (1, 2, 6, 1), + (1, 1, 7, 2), + (1, 2, 8, 2), + (1, 1, 9, 2), + (1, 2, 10, 2), + (1, 1, 11, 2), + (1, 2, 12, 2), + ) + fake_mf.update(fake_shape_dependents(div_seq, sid_dim=0)) + assert MFW(fake_mf).image_shape == (32, 64, 2, 3, 2) + # Check we only allow one extra spatial dimension with unique val per frame + div_seq = ( + (1, 1, 1, 6), + (1, 2, 2, 5), + (1, 1, 3, 4), + (1, 2, 4, 3), + (1, 1, 5, 2), + (1, 2, 6, 1), + ) + fake_mf.update(fake_shape_dependents(div_seq, sid_dim=0)) + with pytest.raises(didw.WrapperError): + MFW(fake_mf).image_shape + # Check that having unique value per frame works with single volume + div_seq = ((1, 1, 1), (1, 2, 2), (1, 3, 3)) + fake_mf.update(fake_shape_dependents(div_seq, sid_dim=0)) + assert MFW(fake_mf).image_shape == (32, 64, 3) def test_iop(self): # Test Image orient patient for multiframe @@ -608,22 +753,30 @@ def test_image_position(self): with pytest.raises(didw.WrapperError): dw.image_position # Make a fake frame - fake_frame = fake_frames( - 'PlanePositionSequence', 'ImagePositionPatient', [[-2.0, 3.0, 7]] - )[0] - fake_mf['SharedFunctionalGroupsSequence'] = [fake_frame] + iop = [0, 1, 0, 1, 0, 0] + frames = fake_frames('PlaneOrientationSequence', 'ImageOrientationPatient', [iop]) + frames = fake_frames( + 'PlanePositionSequence', 'ImagePositionPatient', [[-2.0, 3.0, 7]], frames + ) + fake_mf['SharedFunctionalGroupsSequence'] = frames assert_array_equal(MFW(fake_mf).image_position, [-2, 3, 7]) fake_mf['SharedFunctionalGroupsSequence'] = [None] with pytest.raises(didw.WrapperError): MFW(fake_mf).image_position - fake_mf['PerFrameFunctionalGroupsSequence'] = [fake_frame] + fake_mf['PerFrameFunctionalGroupsSequence'] = frames assert_array_equal(MFW(fake_mf).image_position, [-2, 3, 7]) # Check lists of Decimals work - fake_frame.PlanePositionSequence[0].ImagePositionPatient = [ + frames[0].PlanePositionSequence[0].ImagePositionPatient = [ Decimal(str(v)) for v in [-2, 3, 7] ] assert_array_equal(MFW(fake_mf).image_position, [-2, 3, 7]) assert MFW(fake_mf).image_position.dtype == float + # We should get minimum along slice normal with multiple frames + frames = fake_frames('PlaneOrientationSequence', 'ImageOrientationPatient', [iop] * 2) + ipps = [[-2.0, 3.0, 7], [-2.0, 3.0, 6]] + frames = fake_frames('PlanePositionSequence', 'ImagePositionPatient', ipps, frames) + fake_mf['PerFrameFunctionalGroupsSequence'] = frames + assert_array_equal(MFW(fake_mf).image_position, [-2, 3, 6]) @dicom_test @pytest.mark.xfail(reason='Not packaged in install', raises=FileNotFoundError) @@ -644,7 +797,7 @@ def test_data_real(self): if endian_codes[data.dtype.byteorder] == '>': data = data.byteswap() dat_str = data.tobytes() - assert sha1(dat_str).hexdigest() == '149323269b0af92baa7508e19ca315240f77fa8c' + assert sha1(dat_str).hexdigest() == 'dc011bb49682fb78f3cebacf965cb65cc9daba7d' @dicom_test def test_slicethickness_fallback(self): @@ -665,7 +818,7 @@ def test_data_derived_shape(self): def test_data_trace(self): # Test that a standalone trace volume is found and not dropped dw = didw.wrapper_from_file(DATA_FILE_SIEMENS_TRACE) - assert dw.image_shape == (72, 72, 39, 1) + assert dw.image_shape == (72, 72, 39) @dicom_test @needs_nibabel_data('nitest-dicom') @@ -715,6 +868,11 @@ def test_data_fake(self): sorted_data = data[..., [3, 1, 2, 0]] fake_mf['pixel_array'] = np.rollaxis(sorted_data, 2) assert_array_equal(MFW(fake_mf).get_data(), data * 2.0 - 1) + # Check slice sorting with negative index / IPP correlation + fake_mf.update(fake_shape_dependents(dim_idxs, sid_dim=0, flip_ipp_idx_corr=True)) + sorted_data = data[..., [0, 2, 1, 3]] + fake_mf['pixel_array'] = np.rollaxis(sorted_data, 2) + assert_array_equal(MFW(fake_mf).get_data(), data * 2.0 - 1) # 5D! dim_idxs = [ [1, 4, 2, 1],