From c49dff290f6113327eaa62bbd8aff4da924dd54a Mon Sep 17 00:00:00 2001 From: Brendan Moloney Date: Tue, 18 Jun 2024 16:44:03 -0700 Subject: [PATCH 1/8] BF: Fix for 'split' (concatenated?) multiframe DICOM Can't just use number of frame indices to determine shape of data, as the actual frames could still be split into different files. Also can't assume a multiframe file is more than a single slice. --- nibabel/nicom/dicomwrappers.py | 34 ++++++++++++++++------------------ 1 file changed, 16 insertions(+), 18 deletions(-) diff --git a/nibabel/nicom/dicomwrappers.py b/nibabel/nicom/dicomwrappers.py index 2270ed3f0..894a0ed21 100755 --- a/nibabel/nicom/dicomwrappers.py +++ b/nibabel/nicom/dicomwrappers.py @@ -554,23 +554,20 @@ 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 # Store frame indices self._frame_indices = frame_indices - if n_dim < 4: # 3D volume - return rows, cols, n_frames - # More than 3 dimensions + # Determine size of any extra-spatial 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}.' - ) + shape = (rows, cols) + tuple(x for i, x in enumerate(ns_unique) if i == 0 or x != 1) + n_dim = len(shape) + if n_dim > 3: + 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 @@ -640,10 +637,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) From 4063114c2bde09f34d88c1193a5fd20adc8c1932 Mon Sep 17 00:00:00 2001 From: Brendan Moloney Date: Wed, 24 Jul 2024 15:29:26 -0700 Subject: [PATCH 2/8] BF+TST: Test and fix a bunch of multiframe fixes Corrects issue where order of slice indices was assumed to match the order needed to move along the direction of the slice normal, which resulted in slice orientation flips. Ignores indices that don't evenly divide data, and at the end will try to combine those indices (if needed) into a single tuple index. --- nibabel/nicom/dicomwrappers.py | 124 +++++++++++++++----- nibabel/nicom/tests/test_dicomwrappers.py | 132 ++++++++++++++++++---- 2 files changed, 203 insertions(+), 53 deletions(-) diff --git a/nibabel/nicom/dicomwrappers.py b/nibabel/nicom/dicomwrappers.py index 894a0ed21..c3f484a00 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 = [self.shared.PlanePositionSequence[0].ImagePositionPatient] + except AttributeError: + try: + frame_ipps = [f.PlanePositionSequence[0].ImagePositionPatient for f in self.frames] + 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,20 +575,70 @@ 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) + # Determine the shape and which indices to use + shape = [rows, cols] + curr_parts = n_frames + frames_per_part = 1 + del_indices = {} + for row_idx, row in enumerate(frame_indices.T): + if curr_parts == 1: + break + unique = np.unique(row) + count = len(unique) + if count == 1: + continue + # Replace slice indices with order determined from slice positions along normal + if len(shape) == 2: + 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") + new_parts, leftover = divmod(curr_parts, count) + allowed_val_counts = [new_parts * frames_per_part] + if len(shape) > 2: + # Except for the slice dim, having a unique value for each frame is valid + allowed_val_counts.append(n_frames) + if leftover != 0 or any( + np.count_nonzero(row == val) not in allowed_val_counts for val in unique + ): + if len(shape) == 2: + raise WrapperError('Missing slices from multiframe') + del_indices[row_idx] = count + continue + frames_per_part *= count + shape.append(count) + curr_parts = new_parts + 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 - # Determine size of any extra-spatial dimensions - ns_unique = [len(np.unique(row)) for row in self._frame_indices.T] - shape = (rows, cols) + tuple(x for i, x in enumerate(ns_unique) if i == 0 or x != 1) - n_dim = len(shape) - if n_dim > 3: - 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 @@ -607,18 +678,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): diff --git a/nibabel/nicom/tests/test_dicomwrappers.py b/nibabel/nicom/tests/test_dicomwrappers.py index d14c35dcd..25a58d70e 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 == 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,86 @@ 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) + assert num_of_frames % n_slices == 0 + 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 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 +538,15 @@ 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 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) @@ -541,6 +607,18 @@ 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) + # 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) def test_iop(self): # Test Image orient patient for multiframe @@ -608,22 +686,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 +730,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 +751,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') From 14c24ef7fc156d2a0bb760304e482cfde4694bc3 Mon Sep 17 00:00:00 2001 From: Brendan Moloney Date: Wed, 24 Jul 2024 16:34:11 -0700 Subject: [PATCH 3/8] BF: Trim unneeded trailing indices from _frame_indices --- nibabel/nicom/dicomwrappers.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/nibabel/nicom/dicomwrappers.py b/nibabel/nicom/dicomwrappers.py index c3f484a00..eab0471ec 100755 --- a/nibabel/nicom/dicomwrappers.py +++ b/nibabel/nicom/dicomwrappers.py @@ -581,11 +581,10 @@ def image_shape(self): frames_per_part = 1 del_indices = {} for row_idx, row in enumerate(frame_indices.T): - if curr_parts == 1: - break unique = np.unique(row) count = len(unique) - if count == 1: + if count == 1 or curr_parts == 1: + del_indices[row_idx] = count continue # Replace slice indices with order determined from slice positions along normal if len(shape) == 2: From 019f448c9924e352ed5503aae384b59918bb1d95 Mon Sep 17 00:00:00 2001 From: Brendan Moloney Date: Wed, 24 Jul 2024 17:09:09 -0700 Subject: [PATCH 4/8] BF+TST: Fix 2D plus time case Explicitly use `InStackPositionNumber` to identify the slice dim, produce correct output for 2D + time data. --- nibabel/nicom/dicomwrappers.py | 17 +++++++++++------ nibabel/nicom/tests/test_dicomwrappers.py | 9 ++++++++- 2 files changed, 19 insertions(+), 7 deletions(-) diff --git a/nibabel/nicom/dicomwrappers.py b/nibabel/nicom/dicomwrappers.py index eab0471ec..14041e631 100755 --- a/nibabel/nicom/dicomwrappers.py +++ b/nibabel/nicom/dicomwrappers.py @@ -470,10 +470,10 @@ def __init__(self, dcm_data): # Try to determine slice order and minimal image position patient self._frame_slc_ord = self._ipp = None try: - frame_ipps = [self.shared.PlanePositionSequence[0].ImagePositionPatient] + frame_ipps = [f.PlanePositionSequence[0].ImagePositionPatient for f in self.frames] except AttributeError: try: - frame_ipps = [f.PlanePositionSequence[0].ImagePositionPatient for f in self.frames] + 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): @@ -575,19 +575,24 @@ 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) + 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 count == 1 or curr_parts == 1: + 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 len(shape) == 2: + 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) @@ -595,13 +600,13 @@ def image_shape(self): raise WrapperError("Number of slice indices and positions don't match") new_parts, leftover = divmod(curr_parts, count) allowed_val_counts = [new_parts * frames_per_part] - if len(shape) > 2: + if row_idx != slice_dim_idx: # Except for the slice dim, having a unique value for each frame is valid allowed_val_counts.append(n_frames) if leftover != 0 or any( np.count_nonzero(row == val) not in allowed_val_counts for val in unique ): - if len(shape) == 2: + if row_idx == slice_dim_idx: raise WrapperError('Missing slices from multiframe') del_indices[row_idx] = count continue diff --git a/nibabel/nicom/tests/test_dicomwrappers.py b/nibabel/nicom/tests/test_dicomwrappers.py index 25a58d70e..040242162 100755 --- a/nibabel/nicom/tests/test_dicomwrappers.py +++ b/nibabel/nicom/tests/test_dicomwrappers.py @@ -488,10 +488,13 @@ def __init__(self, div, sid, ipp, iop): 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 = [ @@ -546,6 +549,10 @@ def test_shape(self): 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)) From 0215ce5db008f32e6001335f2b4d4f39d5a0a346 Mon Sep 17 00:00:00 2001 From: Brendan Moloney Date: Wed, 24 Jul 2024 18:33:56 -0700 Subject: [PATCH 5/8] BF+TST: Handle case with extra-spatial index that is unique per frame Not sure if this ever actually happens in real multiframe data, but it does in non-multiframe and I can imagine if a DimensionIndexSequence element refrences a per-frame AcquisitionTime then this could happen. --- nibabel/nicom/dicomwrappers.py | 25 +++++++++++++------ nibabel/nicom/tests/test_dicomwrappers.py | 29 +++++++++++++++++++++++ 2 files changed, 47 insertions(+), 7 deletions(-) diff --git a/nibabel/nicom/dicomwrappers.py b/nibabel/nicom/dicomwrappers.py index 14041e631..374387870 100755 --- a/nibabel/nicom/dicomwrappers.py +++ b/nibabel/nicom/dicomwrappers.py @@ -598,21 +598,32 @@ def image_shape(self): 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) - allowed_val_counts = [new_parts * frames_per_part] - if row_idx != slice_dim_idx: - # Except for the slice dim, having a unique value for each frame is valid - allowed_val_counts.append(n_frames) - if leftover != 0 or any( - np.count_nonzero(row == val) not in allowed_val_counts for val in unique - ): + 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] diff --git a/nibabel/nicom/tests/test_dicomwrappers.py b/nibabel/nicom/tests/test_dicomwrappers.py index 040242162..b50535a4b 100755 --- a/nibabel/nicom/tests/test_dicomwrappers.py +++ b/nibabel/nicom/tests/test_dicomwrappers.py @@ -626,6 +626,35 @@ def test_shape(self): ) 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) def test_iop(self): # Test Image orient patient for multiframe From 259483f1f5412e4e8deb919b800b72abddccd439 Mon Sep 17 00:00:00 2001 From: Brendan Moloney Date: Wed, 24 Jul 2024 23:24:53 -0700 Subject: [PATCH 6/8] TST: Expand test coverage for multiframe dicom shape determination --- nibabel/nicom/tests/test_dicomwrappers.py | 33 ++++++++++++++++++++++- 1 file changed, 32 insertions(+), 1 deletion(-) diff --git a/nibabel/nicom/tests/test_dicomwrappers.py b/nibabel/nicom/tests/test_dicomwrappers.py index b50535a4b..2168476bb 100755 --- a/nibabel/nicom/tests/test_dicomwrappers.py +++ b/nibabel/nicom/tests/test_dicomwrappers.py @@ -473,7 +473,6 @@ def __init__(self, div, sid, ipp, iop): frame_slc_indices = np.array(div_seq)[:, slice_dim] uniq_slc_indices = np.unique(frame_slc_indices) n_slices = len(uniq_slc_indices) - assert num_of_frames % n_slices == 0 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) @@ -579,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) @@ -614,6 +624,11 @@ 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 = ( @@ -655,6 +670,22 @@ def test_shape(self): ) 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 From 52c31052e4f22ff7f0a01883129584c6091e9ac9 Mon Sep 17 00:00:00 2001 From: Brendan Moloney Date: Thu, 25 Jul 2024 09:58:19 -0700 Subject: [PATCH 7/8] TST+CLN: More slice ordering testing, minor cleanup --- nibabel/nicom/tests/test_dicomwrappers.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/nibabel/nicom/tests/test_dicomwrappers.py b/nibabel/nicom/tests/test_dicomwrappers.py index 2168476bb..e01759c86 100755 --- a/nibabel/nicom/tests/test_dicomwrappers.py +++ b/nibabel/nicom/tests/test_dicomwrappers.py @@ -388,7 +388,7 @@ def fake_frames(seq_name, field_name, value_seq, frame_seq=None): class Fake: pass - if frame_seq == None: + 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() @@ -868,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], From 629dbb52e14e813203d1f9c355de95399fd70dda Mon Sep 17 00:00:00 2001 From: Brendan Moloney Date: Thu, 25 Jul 2024 10:05:32 -0700 Subject: [PATCH 8/8] DOC: Add some notes to the changelog --- Changelog | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) 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) ===============================