diff --git a/src/highdicom/image.py b/src/highdicom/image.py index aaeda870..a4b86662 100644 --- a/src/highdicom/image.py +++ b/src/highdicom/image.py @@ -63,6 +63,7 @@ from highdicom.spatial import ( _are_orientations_coplanar, _are_orientations_equal, + _check_orientation_consistency, _get_spatial_information, ImageToReferenceTransformer, compute_tile_positions_per_frame, @@ -7014,6 +7015,7 @@ def get_volume_from_series( atol: float | None = None, rtol: float | None = None, perpendicular_tol: float | None = None, + orientation_tol: float | None = None, ) -> Volume: """Create volume from a series of single frame images. @@ -7162,11 +7164,22 @@ def get_volume_from_series( ): raise ValueError('Images do not share a frame of reference.') - if not all( - ds.ImageOrientationPatient == image_orientation - for ds in series_datasets - ): - raise ValueError('Images do not have the same orientation.') + if orientation_tol is None: + image_orientation = series_datasets[0].ImageOrientationPatient + if not all( + ds.ImageOrientationPatient == image_orientation + for ds in series_datasets + ): + raise ValueError('Images do not have the same orientation.') + else: + image_orientation = _check_orientation_consistency( + series_datasets, orientation_tol + ) + if image_orientation is None: + raise ValueError( + 'Orientations are not consistent within the specified ' + 'tolerance.' + ) if not all( ds.PixelSpacing == pixel_spacing @@ -7183,6 +7196,7 @@ def get_volume_from_series( atol=atol, rtol=rtol, perpendicular_tol=perpendicular_tol, + orientation_tol=orientation_tol, ) if slice_spacing is None: raise ValueError('Series is not a regularly-spaced volume.') diff --git a/src/highdicom/spatial.py b/src/highdicom/spatial.py index 26917de5..763ec64e 100644 --- a/src/highdicom/spatial.py +++ b/src/highdicom/spatial.py @@ -1,4 +1,5 @@ import itertools +from collections import Counter from collections.abc import Generator, Iterator, Sequence import logging from typing_extensions import Self @@ -3403,12 +3404,65 @@ def are_points_coplanar( return max_dev <= tol +def _check_orientation_consistency( + series_datasets: Sequence[Dataset], + orientation_tol: float, +) -> Sequence[float] | None: + """Check orientations are consistent within a tolerance. + + First establishes the most commonly occurring orientation value within the + series, then checks that all other orientation values found are within the + given tolerance. + + Parameters + ---------- + series_datasets: Sequence[pydicom.Dataset] + Series of legacy datasets whose orientations should be checked. + orientation_tol: float + Tolerance value used to determine whether two orientation values are + different. The checks fails if the cosine similarity of either of the + two direction unit vectors within an orientation has a cosine + similarity of less than 1 minus this value from the corresponding unit + vector in the most common orientation. + + Returns + ------- + list[float] | None: + If orientations are the similar within the given tolerance, returns the + most commonly occurring orientation in the ImageOrientationPatient + direction cosine format of a list of six floats. Otherwise returns None. + + """ + orientations = [ + tuple(ds.ImageOrientationPatient) + for ds in series_datasets + ] + counter = Counter(orientations) + + orientation, _ = counter.most_common(1)[0] + v1 = np.array(orientation[:3]) + v2 = np.array(orientation[3:]) + + for ori, _ in counter.items(): + if ori == orientation: + continue + + if ( + np.array(ori[:3]) @ v1 < (1.0 - orientation_tol) or + np.array(ori[3:]) @ v2 < (1.0 - orientation_tol) + ): + return None + + return list(orientation) + + def get_series_volume_positions( datasets: Sequence[pydicom.Dataset], *, rtol: float | None = None, atol: float | None = None, perpendicular_tol: float | None = None, + orientation_tol: float | None = None, sort: bool = True, allow_missing_positions: bool = False, allow_duplicate_positions: bool = False, @@ -3522,9 +3576,17 @@ def get_series_volume_positions( ) # Check image orientations are consistent - image_orientation = datasets[0].ImageOrientationPatient - for ds in datasets[1:]: - if ds.ImageOrientationPatient != image_orientation: + if orientation_tol is None: + image_orientation = datasets[0].ImageOrientationPatient + for ds in datasets[1:]: + if ds.ImageOrientationPatient != image_orientation: + return None, None + else: + image_orientation = _check_orientation_consistency( + datasets, + orientation_tol + ) + if image_orientation is None: return None, None positions = [ds.ImagePositionPatient for ds in datasets] diff --git a/tests/test_image.py b/tests/test_image.py index a963c4c0..47205154 100644 --- a/tests/test_image.py +++ b/tests/test_image.py @@ -22,6 +22,7 @@ from highdicom.content import VOILUTTransformation from highdicom.image import ( _CombinedPixelTransform, + get_volume_from_series, ) from highdicom.pixels import ( apply_voi_window, @@ -1042,6 +1043,39 @@ def test_get_volume_multiframe_ct(): assert volume.array.dtype == dtype +def test_get_volume_from_series_with_tolerance(): + ct_files = [ + get_testdata_file('dicomdirtests/77654033/CT2/17136'), + get_testdata_file('dicomdirtests/77654033/CT2/17196'), + get_testdata_file('dicomdirtests/77654033/CT2/17166'), + ] + ct_series = [pydicom.dcmread(f) for f in ct_files] + + theta = 0.05 + ct_series[0].ImageOrientationPatient = [ + pydicom.valuerep.format_number_as_ds(x) for x in [ + np.cos(theta), np.sin(theta), 0.0, + 0.0, 1.0, 0.0 + ] + ] + + # With no tolerance (default), this should fail + msg = "Images do not have the same orientation." + with pytest.raises(ValueError, match=msg): + get_volume_from_series(ct_series) + + # With a tolerance that is too low, we should see a different error + # msg = "Images do not have the same orientation." + msg = 'Orientations are not consistent within the specified tolerance.' + with pytest.raises(ValueError, match=msg): + get_volume_from_series(ct_series, orientation_tol=0.0000001) + + vol = get_volume_from_series(ct_series, orientation_tol=0.01) + + # The majority orientation should have been used + assert vol.direction_cosines == tuple(ct_series[1].ImageOrientationPatient) + + def test_tiled_full_no_dimension_index(): # The dimension index sequence is optional with TILED_FULL images # Check that the image is read correctly and the same total pixel matrix is