diff --git a/.github/workflows/run_unit_tests.yml b/.github/workflows/run_unit_tests.yml index 13b92d31..1ee0cc75 100644 --- a/.github/workflows/run_unit_tests.yml +++ b/.github/workflows/run_unit_tests.yml @@ -15,7 +15,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: ["3.7", "3.8", "3.9", "3.10", "3.11", "3.12"] + python-version: ["3.10", "3.11", "3.12"] dependencies: [".", "'.[libjpeg]'"] steps: diff --git a/data/test_files/seg_image_sm_control.dcm b/data/test_files/seg_image_sm_control.dcm index bf695236..2140bda6 100644 Binary files a/data/test_files/seg_image_sm_control.dcm and b/data/test_files/seg_image_sm_control.dcm differ diff --git a/data/test_files/seg_image_sm_dots.dcm b/data/test_files/seg_image_sm_dots.dcm index 21572272..62d2b5c2 100644 Binary files a/data/test_files/seg_image_sm_dots.dcm and b/data/test_files/seg_image_sm_dots.dcm differ diff --git a/data/test_files/seg_image_sm_dots_tiled_full.dcm b/data/test_files/seg_image_sm_dots_tiled_full.dcm index 8421dddc..efaddfae 100644 Binary files a/data/test_files/seg_image_sm_dots_tiled_full.dcm and b/data/test_files/seg_image_sm_dots_tiled_full.dcm differ diff --git a/data/test_files/seg_image_sm_numbers.dcm b/data/test_files/seg_image_sm_numbers.dcm index 374672b3..7cdf052a 100644 Binary files a/data/test_files/seg_image_sm_numbers.dcm and b/data/test_files/seg_image_sm_numbers.dcm differ diff --git a/docs/seg.rst b/docs/seg.rst index 6152f65c..2abdda37 100644 --- a/docs/seg.rst +++ b/docs/seg.rst @@ -807,22 +807,31 @@ We recommend that if you do this, you specify ``max_fractional_value=1`` to clearly communicate that the segmentation is inherently binary in nature. Why would you want to make this seemingly rather strange choice? Well, -``"FRACTIONAL"`` SEGs tend to compress *much* better than ``"BINARY"`` ones -(see next section). Note however, that this is arguably an misuse of the intent -of the standard, so *caveat emptor*. +``"FRACTIONAL"`` SEGs tend to compress better than ``"BINARY"`` ones (see next +section). Note however, that this is arguably an misuse of the intent of the +standard, so *caveat emptor*. Also note that while this used to be a more +serious issue it is less serious now that ``"JPEG2000Lossless"`` compression is +now supported for ``"BINARY"`` segmentations as of highdicom v0.23.0. Compression ----------- The types of pixel compression available in segmentation images depends on the -segmentation type. Pixels in a ``"BINARY"`` segmentation image are "bit-packed" -such that 8 pixels are grouped into 1 byte in the stored array. If a given frame -contains a number of pixels that is not divisible by 8 exactly, a single byte +segmentation type. + +Pixels in an uncompressed ``"BINARY"`` segmentation image are "bit-packed" such +that 8 pixels are grouped into 1 byte in the stored array. If a given frame +contains a number of pixels that is not divisible by 8 exactly, a single byte will straddle a frame boundary into the next frame if there is one, or the byte will be padded with zeroes of there are no further frames. This means that -retrieving individual frames from segmentation images in which each frame -size is not divisible by 8 becomes problematic. No further compression may be -applied to frames of ``"BINARY"`` segmentation images. +retrieving individual frames from segmentation images in which each frame size +is not divisible by 8 becomes problematic. For this reason, as well as for +space efficiency (sparse segmentations tend to compress very well), we +therefore strongly recommend using ``"JPEG2000Lossless"`` compression with +``"BINARY"`` segmentations. This is the only compression method currently +supported for ``"BINARY"`` segmentations. However, beware that reading these +single-bit JPEG 2000 images may not be supported by all other tools and +viewers. Pixels in ``"FRACTIONAL"`` segmentation images may be compressed using one of the lossless compression methods available within DICOM. Currently *highdicom* @@ -830,16 +839,6 @@ supports the following compressed transfer syntaxes when creating ``"FRACTIONAL"`` segmentation images: ``"RLELossless"``, ``"JPEG2000Lossless"``, and ``"JPEGLSLossless"``. -Note that there may be advantages to using ``"FRACTIONAL"`` segmentations to -store segmentation images that are binary in nature (i.e. only taking values 0 -and 1): - -- If the segmentation is very simple or sparse, the lossless compression methods - available in ``"FRACTIONAL"`` images may be more effective than the - "bit-packing" method required by ``"BINARY"`` segmentations. -- The clear frame boundaries make retrieving individual frames from - ``"FRACTIONAL"`` image files possible. - Multiprocessing --------------- diff --git a/pyproject.toml b/pyproject.toml index 185e2f26..e8898135 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -7,7 +7,7 @@ name = "highdicom" dynamic = ["version"] description = "High-level DICOM abstractions." readme = "README.md" -requires-python = ">=3.6" +requires-python = ">=3.10" authors = [ { name = "Markus D. Herrmann" }, ] @@ -15,6 +15,7 @@ maintainers = [ { name = "Markus D. Herrmann" }, { name = "Christopher P. Bridge" }, ] +license = { text = "LICENSE" } classifiers = [ "Development Status :: 4 - Beta", "Intended Audience :: Science/Research", @@ -23,10 +24,6 @@ classifiers = [ "Operating System :: Microsoft :: Windows", "Operating System :: POSIX :: Linux", "Programming Language :: Python :: 3", - "Programming Language :: Python :: 3.6", - "Programming Language :: Python :: 3.7", - "Programming Language :: Python :: 3.8", - "Programming Language :: Python :: 3.9", "Programming Language :: Python :: 3.10", "Programming Language :: Python :: 3.11", "Programming Language :: Python :: 3.12", @@ -35,16 +32,16 @@ classifiers = [ ] dependencies = [ "numpy>=1.19", - "pillow-jpls>=1.0", "pillow>=8.3", - "pydicom>=2.3.0,!=2.4.0", + "pydicom>=3.0.1", + "pyjpegls>=1.0.0", ] [project.optional-dependencies] libjpeg = [ - "pylibjpeg-libjpeg>=1.3", - "pylibjpeg-openjpeg>=1.2", - "pylibjpeg>=1.4", + "pylibjpeg-libjpeg>=2.1", + "pylibjpeg-openjpeg>=2.0.0", + "pylibjpeg>=2.0", ] test = [ "mypy==0.971", diff --git a/src/highdicom/_module_utils.py b/src/highdicom/_module_utils.py index eccd5d90..5ebc2673 100644 --- a/src/highdicom/_module_utils.py +++ b/src/highdicom/_module_utils.py @@ -281,3 +281,25 @@ def does_iod_have_pixel_data(sop_class_uid: str) -> bool: return any( is_attribute_in_iod(attr, sop_class_uid) for attr in pixel_attrs ) + + +def is_multiframe_image(dataset: Dataset): + """Determine whether an image is a multiframe image. + The definition used is whether the IOD allows for multiple frames, not + whether this particular instance has more than one frame. + + Parameters + ---------- + dataset: pydicom.Dataset + A dataset to check. + + Returns + ------- + bool: + Whether the image belongs to a multiframe IOD. + + """ + return is_attribute_in_iod( + 'PerFrameFunctionalGroupsSequence', + dataset.SOPClassUID, + ) diff --git a/src/highdicom/base.py b/src/highdicom/base.py index 9a80673a..4583537b 100644 --- a/src/highdicom/base.py +++ b/src/highdicom/base.py @@ -139,8 +139,6 @@ def __init__( "Big Endian transfer syntaxes are retired and no longer " "supported by highdicom." ) - self.is_little_endian = True # backwards compatibility - self.is_implicit_VR = transfer_syntax_uid.is_implicit_VR # Include all File Meta Information required for writing SOP instance # to a file in PS3.10 format. @@ -154,7 +152,6 @@ def __init__( '1.2.826.0.1.3680043.9.7433.1.1' ) self.file_meta.ImplementationVersionName = f'highdicom{__version__}' - self.fix_meta_info(enforce_standard=True) with BytesIO() as fp: write_file_meta_info(fp, self.file_meta, enforce_standard=True) self.file_meta.FileMetaInformationGroupLength = len(fp.getvalue()) diff --git a/src/highdicom/content.py b/src/highdicom/content.py index 2f3f287a..87fea27e 100644 --- a/src/highdicom/content.py +++ b/src/highdicom/content.py @@ -6,11 +6,12 @@ import numpy as np from pydicom.dataset import Dataset +from pydicom import DataElement from pydicom.sequence import Sequence as DataElementSequence from pydicom.sr.coding import Code from pydicom.sr.codedict import codes from pydicom.valuerep import DS, format_number_as_ds -from pydicom._storage_sopclass_uids import SegmentationStorage +from pydicom.uid import SegmentationStorage from highdicom.enum import ( CoordinateSystemNames, @@ -469,18 +470,27 @@ def __init__( 'Position in Pixel Matrix must be specified for ' 'slide coordinate system.' ) - col_position, row_position = pixel_matrix_position x, y, z = image_position - item.XOffsetInSlideCoordinateSystem = DS(x, auto_format=True) - item.YOffsetInSlideCoordinateSystem = DS(y, auto_format=True) - item.ZOffsetInSlideCoordinateSystem = DS(z, auto_format=True) + col_position, row_position = pixel_matrix_position if row_position < 0 or col_position < 0: raise ValueError( 'Both items in "pixel_matrix_position" must be positive ' 'integers.' ) - item.RowPositionInTotalImagePixelMatrix = row_position - item.ColumnPositionInTotalImagePixelMatrix = col_position + + # Use hard-coded tags to avoid the keyword dictionary lookup + # (this constructor is called a large number of times in large + # multiframe images, so some optimization makes sense) + x_tag = 0x0040072a # XOffsetInSlideCoordinateSystem + y_tag = 0x0040073a # YOffsetInSlideCoordinateSystem + z_tag = 0x0040074a # ZOffsetInSlideCoordinateSystem + row_tag = 0x0048021f # RowPositionInTotalImagePixelMatrix + column_tag = 0x0048021e # ColumnPositionInTotalImagePixelMatrix + item.add(DataElement(x_tag, 'DS', DS(x, auto_format=True))) + item.add(DataElement(y_tag, 'DS', DS(y, auto_format=True))) + item.add(DataElement(z_tag, 'DS', DS(z, auto_format=True))) + item.add(DataElement(row_tag, 'SL', int(row_position))) + item.add(DataElement(column_tag, 'SL', int(col_position))) elif coordinate_system == CoordinateSystemNames.PATIENT: item.ImagePositionPatient = [ DS(ip, auto_format=True) for ip in image_position diff --git a/src/highdicom/frame.py b/src/highdicom/frame.py index 6a98fd8b..47b959b7 100644 --- a/src/highdicom/frame.py +++ b/src/highdicom/frame.py @@ -6,14 +6,16 @@ from PIL import Image from pydicom.dataset import Dataset, FileMetaDataset from pydicom.encaps import encapsulate -from pydicom.pixel_data_handlers.numpy_handler import pack_bits -from pydicom.pixel_data_handlers.rle_handler import rle_encode_frame +from pydicom.pixels.utils import pack_bits +from pydicom.pixels.encoders.base import get_encoder from pydicom.uid import ( ExplicitVRLittleEndian, ImplicitVRLittleEndian, JPEG2000Lossless, + JPEG2000, JPEGBaseline8Bit, JPEGLSLossless, + JPEGLSNearLossless, UID, RLELossless, ) @@ -51,7 +53,15 @@ def encode_frame( bits_stored: int Number of bits that are required to store a pixel sample photometric_interpretation: Union[PhotometricInterpretationValues, str] - Photometric interpretation + Photometric interpretation that will be used to store data. Usually, + this will match the photometric interpretation of the input pixel + array, however for ``"JPEGBaseline8Bit"``, ``"JPEG2000"``, and + ``"JPEG2000Lossless"`` transfer syntaxes with color images, the pixel + data must be passed in in RGB format and will be converted and stored + as ``"YBR_FULL_422"`` (``"JPEGBaseline8Bit"``), ``"YBR_ICT"`` + (``"JPEG2000"``), or `"YBR_RCT"`` (``"JPEG2000Lossless"``). In these + cases the values of photometric metric passed must match those given + above. pixel_representation: Union[highdicom.PixelRepresentationValues, int, None], optional Whether pixel samples are represented as unsigned integers or 2's complements @@ -111,8 +121,10 @@ def encode_frame( } compressed_transfer_syntaxes = { JPEGBaseline8Bit, + JPEG2000, JPEG2000Lossless, JPEGLSLossless, + JPEGLSNearLossless, RLELossless, } supported_transfer_syntaxes = uncompressed_transfer_syntaxes.union( @@ -132,6 +144,16 @@ def encode_frame( 'Planar configuration must be 0 for color image frames ' 'with native encoding.' ) + allowable_pis = { + 1: ['MONOCHROME1', 'MONOCHROME2', 'PALETTE_COLOR'], + 3: ['RGB', 'YBR_FULL'], + }[samples_per_pixel] + if photometric_interpretation not in allowable_pis: + raise ValueError( + 'Photometric_interpretation of ' + f"'{photometric_interpretation}' " + f'not supported for samples_per_pixel={samples_per_pixel}.' + ) if bits_allocated == 1: if (rows * cols * samples_per_pixel) % 8 != 0: raise ValueError( @@ -142,129 +164,86 @@ def encode_frame( else: return array.flatten().tobytes() - else: - compression_lut = { - JPEGBaseline8Bit: ( - 'jpeg', - { - 'quality': 95 - }, - ), - JPEG2000Lossless: ( - 'jpeg2000', - { - 'tile_size': None, - 'num_resolutions': 1, - 'irreversible': False, - 'no_jp2': True, - }, - ), - JPEGLSLossless: ( - 'JPEG-LS', - { - 'near_lossless': 0, - } - ) - } - - if transfer_syntax_uid == JPEGBaseline8Bit: - if samples_per_pixel == 1: - if planar_configuration is not None: - raise ValueError( - 'Planar configuration must be absent for encoding of ' - 'monochrome image frames with JPEG Baseline codec.' - ) - if photometric_interpretation not in ( - 'MONOCHROME1', 'MONOCHROME2' - ): - raise ValueError( - 'Photometric intpretation must be either "MONOCHROME1" ' - 'or "MONOCHROME2" for encoding of monochrome image ' - 'frames with JPEG Baseline codec.' - ) - elif samples_per_pixel == 3: - if photometric_interpretation != 'YBR_FULL_422': - raise ValueError( - 'Photometric intpretation must be "YBR_FULL_422" for ' - 'encoding of color image frames with ' - 'JPEG Baseline codec.' - ) - if planar_configuration != 0: - raise ValueError( - 'Planar configuration must be 0 for encoding of ' - 'color image frames with JPEG Baseline codec.' - ) - else: + elif transfer_syntax_uid == JPEGBaseline8Bit: + if samples_per_pixel == 1: + if planar_configuration is not None: raise ValueError( - 'Samples per pixel must be 1 or 3 for ' - 'encoding of image frames with JPEG Baseline codec.' + 'Planar configuration must be absent for encoding of ' + 'monochrome image frames with JPEG Baseline codec.' ) - if bits_allocated != 8 or bits_stored != 8: + if photometric_interpretation not in ( + 'MONOCHROME1', 'MONOCHROME2' + ): raise ValueError( - 'Bits allocated and bits stored must be 8 for ' - 'encoding of image frames with JPEG Baseline codec.' + 'Photometric intpretation must be either "MONOCHROME1" ' + 'or "MONOCHROME2" for encoding of monochrome image ' + 'frames with JPEG Baseline codec.' ) - if pixel_representation != 0: + elif samples_per_pixel == 3: + if photometric_interpretation != 'YBR_FULL_422': raise ValueError( - 'Pixel representation must be 0 for ' - 'encoding of image frames with JPEG Baseline codec.' + 'Photometric intpretation must be "YBR_FULL_422" for ' + 'encoding of color image frames with ' + 'JPEG Baseline codec.' ) - - elif transfer_syntax_uid == JPEG2000Lossless: - if samples_per_pixel == 1: - if planar_configuration is not None: - raise ValueError( - 'Planar configuration must be absent for encoding of ' - 'monochrome image frames with Lossless JPEG 2000 codec.' - ) - if photometric_interpretation not in ( - 'MONOCHROME1', 'MONOCHROME2' - ): - raise ValueError( - 'Photometric intpretation must be either "MONOCHROME1" ' - 'or "MONOCHROME2" for encoding of monochrome image ' - 'frames with Lossless JPEG 2000 codec.' - ) - if bits_allocated not in (8, 16): - raise ValueError( - 'Bits Allocated must be 8 or 16 for encoding of ' - 'monochrome image frames with Lossless JPEG 2000 codec.' - ) - elif samples_per_pixel == 3: - if photometric_interpretation != 'YBR_FULL': - raise ValueError( - 'Photometric interpretation must be "YBR_FULL" for ' - 'encoding of color image frames with ' - 'Lossless JPEG 2000 codec.' - ) - if planar_configuration != 0: - raise ValueError( - 'Planar configuration must be 0 for encoding of ' - 'color image frames with Lossless JPEG 2000 codec.' - ) - if bits_allocated != 8: - raise ValueError( - 'Bits Allocated must be 8 for encoding of ' - 'color image frames with Lossless JPEG 2000 codec.' - ) - else: + if planar_configuration != 0: raise ValueError( - 'Samples per pixel must be 1 or 3 for ' - 'encoding of image frames with Lossless JPEG 2000 codec.' + 'Planar configuration must be 0 for encoding of ' + 'color image frames with JPEG Baseline codec.' ) + else: + raise ValueError( + 'Samples per pixel must be 1 or 3 for ' + 'encoding of image frames with JPEG Baseline codec.' + ) + if bits_allocated != 8 or bits_stored != 8: + raise ValueError( + 'Bits allocated and bits stored must be 8 for ' + 'encoding of image frames with JPEG Baseline codec.' + ) + if pixel_representation != 0: + raise ValueError( + 'Pixel representation must be 0 for ' + 'encoding of image frames with JPEG Baseline codec.' + ) + + # Pydicom does not have an encoder for JPEGBaseline8Bit so + # we do this manually + if samples_per_pixel == 3: + image = Image.fromarray(array, mode='RGB') + else: + image = Image.fromarray(array) + with BytesIO() as buf: + image.save(buf, format='jpeg', quality=95) + data = buf.getvalue() + else: + name = { + JPEG2000: "JPEG 2000", + JPEG2000Lossless: "Lossless JPEG 2000", + JPEGLSLossless: "Lossless JPEG-LS", + JPEGLSNearLossless: "Near-Lossless JPEG-LS", + RLELossless: "RLE Lossless", + }[transfer_syntax_uid] + + kwargs = {} + + if samples_per_pixel not in (1, 3): + raise ValueError( + 'Samples per pixel must be 1 or 3 for ' + f'encoding of image frames with {name} codec.' + ) + + if transfer_syntax_uid != RLELossless: if pixel_representation != 0: raise ValueError( 'Pixel representation must be 0 for ' - 'encoding of image frames with Lossless JPEG 2000 codec.' + f'encoding of image frames with {name} codec.' ) - - elif transfer_syntax_uid == JPEGLSLossless: - import pillow_jpls # noqa if samples_per_pixel == 1: if planar_configuration is not None: raise ValueError( 'Planar configuration must be absent for encoding of ' - 'monochrome image frames with Lossless JPEG-LS codec.' + f'monochrome image frames with {name} codec.' ) if photometric_interpretation not in ( 'MONOCHROME1', 'MONOCHROME2' @@ -272,55 +251,98 @@ def encode_frame( raise ValueError( 'Photometric intpretation must be either "MONOCHROME1" ' 'or "MONOCHROME2" for encoding of monochrome image ' - 'frames with Lossless JPEG-LS codec.' - ) - if bits_allocated not in (8, 16): - raise ValueError( - 'Bits Allocated must be 8 or 16 for encoding of ' - 'monochrome image frames with Lossless JPEG-LS codec.' + f'frames with with {name} codec.' ) + if transfer_syntax_uid == JPEG2000Lossless: + if bits_allocated not in (1, 8, 16): + raise ValueError( + 'Bits Allocated must be 1, 8, or 16 for encoding ' + f'of monochrome image frames with with {name} ' + 'codec.' + ) + else: + if bits_allocated not in (8, 16): + raise ValueError( + 'Bits Allocated must be 8 or 16 for encoding of ' + f'monochrome image frames with with {name} codec.' + ) elif samples_per_pixel == 3: - if photometric_interpretation != 'YBR_FULL': - raise ValueError( - 'Photometric interpretation must be "YBR_FULL" for ' - 'encoding of color image frames with ' - 'Lossless JPEG-LS codec.' - ) if planar_configuration != 0: raise ValueError( 'Planar configuration must be 0 for encoding of ' - 'color image frames with Lossless JPEG-LS codec.' + f'color image frames with {name} codec.' + ) + if bits_allocated not in (8, 16): + raise ValueError( + 'Bits Allocated must be 8 or 16 for encoding of ' + f'color image frames with {name} codec.' ) - if bits_allocated != 8: + + required_pi = { + JPEG2000: PhotometricInterpretationValues.YBR_ICT, + JPEG2000Lossless: PhotometricInterpretationValues.YBR_RCT, + JPEGLSLossless: PhotometricInterpretationValues.RGB, + JPEGLSNearLossless: PhotometricInterpretationValues.RGB, + }[transfer_syntax_uid] + + if photometric_interpretation != required_pi.value: raise ValueError( - 'Bits Allocated must be 8 for encoding of ' - 'color image frames with Lossless JPEG-LS codec.' + f'Photometric interpretation must be ' + f'"{required_pi.value}" for encoding of color image ' + f'frames with {name} codec.' ) - else: + + if transfer_syntax_uid == JPEG2000: + kwargs = {'j2k_psnr': [100]} + + if transfer_syntax_uid in (JPEG2000, JPEG2000Lossless): + # This seems to be an openjpeg limitation + if array.shape[0] < 32 or array.shape[1] < 32: raise ValueError( - 'Samples per pixel must be 1 or 3 for ' - 'encoding of image frames with Lossless JPEG-LS codec.' + 'Images smaller than 32 pixels along both dimensions ' + f'cannot be encoded with {name} codec.' ) - if pixel_representation != 0: - raise ValueError( - 'Pixel representation must be 0 for ' - 'encoding of image frames with Lossless JPEG-LS codec.' + + if transfer_syntax_uid == JPEG2000Lossless and bits_allocated == 1: + # Single bit JPEG2000 compression. Pydicom doesn't (yet) support + # this case + if array.dtype != bool: + if array.max() > 1: + raise ValueError( + 'Array must contain only 0 and 1 for bits_allocated = 1' + ) + array = array.astype(bool) + + try: + from openjpeg.utils import encode_array + except ModuleNotFoundError: + raise ModuleNotFoundError( + "Highdicom requires the pylibjpeg-openjpeg package to " + "compress frames using the JPEG2000Lossless transfer " + "syntax." ) - if transfer_syntax_uid in compression_lut: - image_format, kwargs = compression_lut[transfer_syntax_uid] - if samples_per_pixel == 3: - image = Image.fromarray(array, mode='RGB') - else: - image = Image.fromarray(array) - with BytesIO() as buf: - image.save(buf, format=image_format, **kwargs) - data = buf.getvalue() - elif transfer_syntax_uid == RLELossless: - data = rle_encode_frame(array) + data = encode_array( + array, + bits_stored=1, + photometric_interpretation=2, + use_mct=False, + ) else: - raise ValueError( - f'Transfer Syntax "{transfer_syntax_uid}" is not supported.' + encoder = get_encoder(transfer_syntax_uid) + + data = encoder.encode( + array, + rows=array.shape[0], + columns=array.shape[1], + samples_per_pixel=samples_per_pixel, + number_of_frames=1, + bits_allocated=bits_allocated, + bits_stored=bits_stored, + photometric_interpretation=photometric_interpretation, + pixel_representation=pixel_representation, + planar_configuration=planar_configuration, + **kwargs, ) return data @@ -437,16 +459,4 @@ def decode_frame( array = ds.pixel_array - # In case of the JPEG baseline transfer syntax, the pixel_array property - # does not convert the pixel data into the correct (or let's say expected) - # color space after decompression. - if ( - 'YBR' in ds.PhotometricInterpretation and - ds.SamplesPerPixel == 3 and - transfer_syntax_uid == JPEGBaseline8Bit - ): - image = Image.fromarray(array, mode='YCbCr') - image = image.convert(mode='RGB') - array = np.asarray(image) - return array diff --git a/src/highdicom/io.py b/src/highdicom/io.py index 67417f8c..1bccfdef 100644 --- a/src/highdicom/io.py +++ b/src/highdicom/io.py @@ -8,7 +8,7 @@ import numpy as np import pydicom from pydicom.dataset import Dataset -from pydicom.encaps import get_frame_offsets +from pydicom.encaps import parse_basic_offsets from pydicom.filebase import DicomFile, DicomFileLike, DicomBytesIO from pydicom.filereader import ( data_element_offset_to_value, @@ -16,7 +16,7 @@ read_file_meta_info, read_partial ) -from pydicom.pixel_data_handlers.numpy_handler import unpack_bits +from pydicom.pixels.utils import unpack_bits from pydicom.tag import TupleTag, ItemTag, SequenceDelimiterTag from pydicom.uid import UID @@ -120,7 +120,7 @@ def _read_bot(fp: DicomFileLike) -> List[int]: fp.is_implicit_VR, 'OB' ) fp.seek(pixel_data_element_value_offset - 4, 1) - is_empty, offsets = get_frame_offsets(fp) + offsets = parse_basic_offsets(fp) return offsets @@ -249,7 +249,7 @@ def __init__(self, filename: Union[str, Path, DicomFileLike]): DICOM Part10 file containing a dataset of an image SOP Instance """ - if isinstance(filename, DicomFileLike): + if isinstance(filename, (DicomFileLike, DicomBytesIO)): fp = filename self._fp = fp if isinstance(filename, DicomBytesIO): @@ -261,8 +261,8 @@ def __init__(self, filename: Union[str, Path, DicomFileLike]): self._fp = None else: raise TypeError( - 'Argument "filename" must either an open DICOM file object or ' - 'the path to a DICOM file stored on disk.' + 'Argument "filename" must be either an open DICOM file object ' + 'or the path to a DICOM file stored on disk.' ) self._metadata = None diff --git a/src/highdicom/pr/content.py b/src/highdicom/pr/content.py index 701e6cba..ea061dac 100644 --- a/src/highdicom/pr/content.py +++ b/src/highdicom/pr/content.py @@ -34,7 +34,7 @@ ) from highdicom.sr.coding import CodedConcept from highdicom.uid import UID -from highdicom.utils import is_tiled_image +from highdicom.spatial import is_tiled_image from highdicom.valuerep import ( check_person_name, _check_code_string, diff --git a/src/highdicom/pr/sop.py b/src/highdicom/pr/sop.py index d3b7afca..49b86861 100644 --- a/src/highdicom/pr/sop.py +++ b/src/highdicom/pr/sop.py @@ -8,7 +8,7 @@ from pydicom import Dataset from pydicom.sr.coding import Code from pydicom.uid import ExplicitVRLittleEndian -from pydicom._storage_sopclass_uids import ( +from pydicom.uid import ( AdvancedBlendingPresentationStateStorage, ColorSoftcopyPresentationStateStorage, GrayscaleSoftcopyPresentationStateStorage, diff --git a/src/highdicom/sc/sop.py b/src/highdicom/sc/sop.py index df8affce..0646d4ae 100644 --- a/src/highdicom/sc/sop.py +++ b/src/highdicom/sc/sop.py @@ -16,8 +16,10 @@ ExplicitVRLittleEndian, RLELossless, JPEGBaseline8Bit, + JPEG2000, JPEG2000Lossless, JPEGLSLossless, + JPEGLSNearLossless, ) from highdicom.base import SOPClass @@ -108,9 +110,17 @@ def __init__( image; either a 2D grayscale image or a 3D color image (RGB color space) photometric_interpretation: Union[str, highdicom.PhotometricInterpretationValues] - Interpretation of pixel data; either ``"MONOCHROME1"`` or - ``"MONOCHROME2"`` for 2D grayscale images or ``"RGB"`` or - ``"YBR_FULL"`` for 3D color images + Interpretation with which to store pixel data in the dataset; + either ``"MONOCHROME1"`` or ``"MONOCHROME2"`` for 2D grayscale + images or ``"RGB"`` or ``"YBR_FULL"`` for 3D color images. Note + that this should match the photometric interpretation of the input + pixel array, except in the following cases: if + ``transfer_syntax_uid`` is ``"JPEGBaseline8Bit"``, + ``photometric_interpretation must be ``"YBR_FULL_422"``, if + ``transfer_syntax_uid`` is ``"JPEG2000"``, + ``photometric_interpretation must be ``"YBR_ICT"``, if + ``transfer_syntax_uid`` is ``"JPEG2000Lossless"``, + ``photometric_interpretation must be ``"YBR_RCT"``. bits_allocated: int Number of bits that should be allocated per pixel value coordinate_system: Union[str, highdicom.CoordinateSystemNames] @@ -188,8 +198,10 @@ def __init__( ExplicitVRLittleEndian, RLELossless, JPEGBaseline8Bit, + JPEG2000, JPEG2000Lossless, JPEGLSLossless, + JPEGLSNearLossless, } if transfer_syntax_uid not in supported_transfer_syntaxes: raise ValueError( @@ -326,12 +338,23 @@ def __init__( photometric_interpretation ) if pixel_array.ndim == 3: - accepted_interpretations = { - PhotometricInterpretationValues.RGB.value, - PhotometricInterpretationValues.YBR_FULL.value, - PhotometricInterpretationValues.YBR_FULL_422.value, - PhotometricInterpretationValues.YBR_PARTIAL_420.value, - } + if transfer_syntax_uid == JPEGBaseline8Bit: + accepted_interpretations = { + PhotometricInterpretationValues.YBR_FULL_422.value, + } + elif transfer_syntax_uid == JPEG2000: + accepted_interpretations = { + PhotometricInterpretationValues.YBR_ICT.value, + } + elif transfer_syntax_uid == JPEG2000Lossless: + accepted_interpretations = { + PhotometricInterpretationValues.YBR_RCT.value, + } + else: + accepted_interpretations = { + PhotometricInterpretationValues.RGB.value, + PhotometricInterpretationValues.YBR_FULL.value, + } if photometric_interpretation.value not in accepted_interpretations: raise ValueError( 'Pixel array has an unexpected photometric interpretation.' diff --git a/src/highdicom/seg/pyramid.py b/src/highdicom/seg/pyramid.py index af7f7ba0..e7e64fa4 100644 --- a/src/highdicom/seg/pyramid.py +++ b/src/highdicom/seg/pyramid.py @@ -77,7 +77,10 @@ def create_segmentation_pyramid( List of source images. If there are multiple source images, they should be from the same series and pyramid. pixel_arrays: Sequence[numpy.ndarray] - List of segmentation pixel arrays. Each should be a total pixel matrix. + List of segmentation pixel arrays. Each should be a total pixel matrix, + i.e. have shape (rows, columns), (1, rows, columns), or (1, rows, + columns, segments). Otherwise all options supported by the constructor + of :class:`highdicom.seg.Segmentation` are permitted. segmentation_type: Union[str, highdicom.seg.SegmentationTypeValues] Type of segmentation, either ``"BINARY"`` or ``"FRACTIONAL"`` segment_descriptions: Sequence[highdicom.seg.SegmentDescription] @@ -123,16 +126,17 @@ def create_segmentation_pyramid( A human readable label for the output pyramid. **kwargs: Any Any further parameters are passed directly to the constructor of the - :class:highdicom.seg.Segmentation object. However the following + :class:`highdicom.seg.Segmentation` object. However the following parameters are disallowed: ``instance_number``, ``sop_instance_uid``, ``plane_orientation``, ``plane_positions``, ``pixel_measures``, ``pixel_array``, ``tile_pixel_array``. Note ---- - Downsampling is performed via simple nearest neighbor interpolation. If - more control is needed over the downsampling process (for example - anti-aliasing), explicitly pass the downsampled arrays. + Downsampling is performed via simple nearest neighbor interpolation (for + ``BINARY`` segmentations) or bi-linear interpolation (for ``FRACTIONAL`` + segmentations). If more control is needed over the downsampling process + (for example anti-aliasing), explicitly pass the downsampled arrays. """ # Disallow duplicate items in kwargs @@ -149,9 +153,11 @@ def create_segmentation_pyramid( if len(error_keys) > 0: raise TypeError( f'kwargs supplied to the create_segmentation_pyramid function ' - f'should not contain a value for parameter {error_keys[0]}.' + f'should not contain a value for parameter {list(error_keys)[0]}.' ) + segmentation_type = SegmentationTypeValues(segmentation_type) + if pyramid_uid is None: pyramid_uid = UID() if series_instance_uid is None: @@ -263,6 +269,30 @@ def create_segmentation_pyramid( ) # Check that pixel arrays have an appropriate shape + if len(set(p.ndim for p in pixel_arrays)) != 1: + raise ValueError( + 'Each item of argument "pixel_arrays" must have the same number of ' + 'dimensions.' + ) + if pixel_arrays[0].ndim == 4: + n_segment_channels = pixel_arrays[0].shape[3] + else: + n_segment_channels = None + + dtype = pixel_arrays[0].dtype + if dtype in (np.bool_, np.uint8, np.uint16): + resampler = Image.Resampling.NEAREST + elif dtype in (np.float32, np.float64): + if segmentation_type == SegmentationTypeValues.FRACTIONAL: + resampler = Image.Resampling.BILINEAR + else: + # This is a floating point image that will ultimately be treated as + # binary + resampler = Image.Resampling.NEAREST + else: + raise TypeError('Pixel array has an invalid data type.') + + # Checks on consistency of the pixel arrays for pixel_array in pixel_arrays: if pixel_array.ndim not in (2, 3, 4): raise ValueError( @@ -274,6 +304,17 @@ def create_segmentation_pyramid( 'Each item of argument "pixel_arrays" must contain a single ' 'frame, with a size of 1 along dimension 0.' ) + if pixel_array.dtype != dtype: + raise TypeError( + 'Each item of argument "pixel_arrays" must have ' + 'the same datatype.' + ) + if pixel_array.ndim == 4: + if pixel_array.shape[3] != n_segment_channels: + raise ValueError( + 'Each item of argument "pixel_arrays" must have ' + 'the same shape down axis 3.' + ) # Check the pixel arrays are appropriately ordered for index in range(1, len(pixel_arrays)): @@ -322,7 +363,17 @@ def create_segmentation_pyramid( if n_pix_arrays == 1: # Create a pillow image for use later with resizing - mask_image = Image.fromarray(pixel_arrays[0]) + if pixel_arrays[0].ndim == 2: + mask_images = [Image.fromarray(pixel_arrays[0])] + elif pixel_arrays[0].ndim == 3: + # Remove frame dimension before casting + mask_images = [Image.fromarray(pixel_arrays[0][0])] + else: # ndim = 4 + # One "Image" for each segment + mask_images = [ + Image.fromarray(pixel_arrays[0][0, :, :, i]) + for i in range(pixel_arrays[0].shape[3]) + ] all_segs = [] @@ -336,16 +387,14 @@ def create_segmentation_pyramid( if n_pix_arrays > 1: pixel_array = pixel_arrays[output_level] else: - need_resize = True - if n_sources > 1: - output_size = ( - source_image.TotalPixelMatrixColumns, - source_image.TotalPixelMatrixRows - ) + if output_level == 0: + pixel_array = pixel_arrays[0] else: - if output_level == 0: - pixel_array = pixel_arrays[0] - need_resize = False + if n_sources > 1: + output_size = ( + source_image.TotalPixelMatrixColumns, + source_image.TotalPixelMatrixRows + ) else: f = downsample_factors[output_level - 1] output_size = ( @@ -353,10 +402,14 @@ def create_segmentation_pyramid( int(source_images[0].TotalPixelMatrixRows / f) ) - if need_resize: - pixel_array = np.array( - mask_image.resize(output_size, Image.Resampling.NEAREST) - ) + resized_masks = [ + np.array(im.resize(output_size, resampler)) + for im in mask_images + ] + if len(resized_masks) > 1: + pixel_array = np.stack(resized_masks, axis=-1)[None] + else: + pixel_array = resized_masks[0][None] if n_sources == 1: source_pixel_measures = ( diff --git a/src/highdicom/seg/sop.py b/src/highdicom/seg/sop.py index f1bcd50a..00945d63 100644 --- a/src/highdicom/seg/sop.py +++ b/src/highdicom/seg/sop.py @@ -4,6 +4,7 @@ from concurrent.futures import Executor, Future, ProcessPoolExecutor from contextlib import contextmanager from copy import deepcopy +from itertools import chain from os import PathLike import sqlite3 from typing import ( @@ -24,11 +25,12 @@ import warnings import numpy as np +from pydicom.dataelem import DataElement from pydicom.dataset import Dataset from pydicom.datadict import get_entry, keyword_for_tag, tag_for_keyword from pydicom.encaps import encapsulate from pydicom.multival import MultiValue -from pydicom.pixel_data_handlers.numpy_handler import pack_bits +from pydicom.pixels.utils import pack_bits from pydicom.tag import BaseTag, Tag from pydicom.uid import ( ExplicitVRLittleEndian, @@ -43,7 +45,11 @@ from pydicom.sr.coding import Code from pydicom.filereader import dcmread -from highdicom._module_utils import ModuleUsageValues, get_module_usage +from highdicom._module_utils import ( + ModuleUsageValues, + get_module_usage, + does_iod_have_pixel_data, +) from highdicom.base import SOPClass, _check_little_endian from highdicom.content import ( ContentCreatorIdentificationCodeSequence, @@ -58,11 +64,6 @@ from highdicom.frame import encode_frame from highdicom.utils import ( are_plane_positions_tiled_full, - compute_plane_position_tiled_full, - is_tiled_image, - get_tile_array, - iter_tiled_full_frame_data, - tile_pixel_matrix, ) from highdicom.seg.content import ( DimensionIndexSequence, @@ -76,7 +77,14 @@ SegmentAlgorithmTypeValues, ) from highdicom.seg.utils import iter_segments -from highdicom.spatial import ImageToReferenceTransformer +from highdicom.spatial import ( + ImageToReferenceTransformer, + compute_tile_positions_per_frame, + get_image_coordinate_system, + get_tile_array, + is_tiled_image, + iter_tiled_full_frame_data, +) from highdicom.sr.coding import CodedConcept from highdicom.valuerep import ( check_person_name, @@ -91,6 +99,14 @@ _NO_FRAME_REF_VALUE = -1 +# These codes are needed many times in loops so we precompute them +_DERIVATION_CODE = CodedConcept.from_code( + codes.cid7203.SegmentationImageDerivation +) +_PURPOSE_CODE = CodedConcept.from_code( + codes.cid7202.SourceImageForImageProcessingOperation +) + def _get_unsigned_dtype(max_val: Union[int, np.integer]) -> type: """Get the smallest unsigned NumPy datatype to accommodate a value. @@ -1180,6 +1196,7 @@ def __init__( tile_size: Union[Sequence[int], None] = None, pyramid_uid: Optional[str] = None, pyramid_label: Optional[str] = None, + further_source_images: Optional[Sequence[Dataset]] = None, **kwargs: Any ) -> None: """ @@ -1389,6 +1406,13 @@ def __init__( Human readable label for the pyramid containing this segmentation. Should only be used if this segmentation is part of a multi-resolution pyramid. + further_source_images: Optional[Sequence[pydicom.Dataset]], optional + Additional images to record as source images in the segmentation. + Unlike the main ``source_images`` parameter, these images will + *not* be used to infer the position and orientation of the + ``pixel_array`` in the case that no plane positions are supplied. + Images from multiple series may be passed, however they must all + belong to the same study. **kwargs: Any, optional Additional keyword arguments that will be passed to the constructor of `highdicom.base.SOPClass` @@ -1515,6 +1539,7 @@ def __init__( CoordinateSystemNames.SLIDE else: self._coordinate_system = CoordinateSystemNames.PATIENT + self._coordinate_system = get_image_coordinate_system(src_img) else: # Only allow missing FrameOfReferenceUID if it is not required # for this IOD @@ -1551,14 +1576,42 @@ def __init__( ) self._coordinate_system = None + # Remember whether these values were provided by the user, or inferred + # from the source image. If inferred, we can skip some checks + user_provided_orientation = plane_orientation is not None + user_provided_measures = pixel_measures is not None + # General Reference + if further_source_images is not None: + # We make no requirement here that images should be from the same + # series etc, but they should belong to the same study and be image + # objects + for s_img in further_source_images: + if not isinstance(s_img, Dataset): + raise TypeError( + "All items in 'further_source_images' should be " + "of type 'pydicom.Dataset'." + ) + if s_img.StudyInstanceUID != self.StudyInstanceUID: + raise ValueError( + "All items in 'further_source_images' should belong " + "to the same study as 'source_images'." + ) + if not does_iod_have_pixel_data(s_img.SOPClassUID): + raise ValueError( + "All items in 'further_source_images' should be " + "image objects." + ) + else: + further_source_images = [] + # Note that appending directly to the SourceImageSequence is typically # slow so it's more efficient to build as a Python list then convert # later. We save conversion for after the main loop source_image_seq: List[Dataset] = [] referenced_series: Dict[str, List[Dataset]] = defaultdict(list) - for s_img in source_images: + for s_img in chain(source_images, further_source_images): ref = Dataset() ref.ReferencedSOPClassUID = s_img.SOPClassUID ref.ReferencedSOPInstanceUID = s_img.SOPInstanceUID @@ -1617,7 +1670,10 @@ def __init__( if self.SegmentationType == SegmentationTypeValues.BINARY.value: self.BitsAllocated = 1 self.HighBit = 0 - if self.file_meta.TransferSyntaxUID.is_encapsulated: + if ( + self.file_meta.TransferSyntaxUID != JPEG2000Lossless and + self.file_meta.TransferSyntaxUID.is_encapsulated + ): raise ValueError( 'The chosen transfer syntax ' f'{self.file_meta.TransferSyntaxUID} ' @@ -1717,18 +1773,6 @@ def __init__( self.DimensionIndexSequence[0].DimensionOrganizationUID self.DimensionOrganizationSequence = [dimension_organization] - if has_ref_frame_uid: - if is_multiframe: - source_plane_positions = \ - self.DimensionIndexSequence.get_plane_positions_of_image( - src_img - ) - else: - source_plane_positions = \ - self.DimensionIndexSequence.get_plane_positions_of_series( - source_images - ) - if pixel_measures is not None: sffg_item.PixelMeasuresSequence = pixel_measures if ( @@ -1758,12 +1802,29 @@ def __init__( if has_ref_frame_uid: if tile_pixel_array: + src_origin_seq = src_img.TotalPixelMatrixOriginSequence[0] + src_x_offset = src_origin_seq.XOffsetInSlideCoordinateSystem + src_y_offset = src_origin_seq.YOffsetInSlideCoordinateSystem + src_z_offset = src_origin_seq.get( + 'ZOffsetInSlideCoordinateSystem', + 0.0, + ) + if plane_positions is None: # Use the origin of the source image - origin_seq = src_img.TotalPixelMatrixOriginSequence[0] - x_offset = origin_seq.XOffsetInSlideCoordinateSystem - y_offset = origin_seq.YOffsetInSlideCoordinateSystem + x_offset = src_x_offset + y_offset = src_y_offset + z_offset = src_z_offset + origin_preserved = True else: + if len(plane_positions) != 1: + raise ValueError( + "If specifying plane_positions when the " + '"tile_pixel_array" argument is True, a ' + "single plane position should be provided " + "representing the position of the top " + "left corner of the total pixel matrix." + ) # Use the provided image origin pp = plane_positions[0][0] rp = pp.RowPositionInTotalImagePixelMatrix @@ -1779,28 +1840,129 @@ def __init__( ) x_offset = pp.XOffsetInSlideCoordinateSystem y_offset = pp.YOffsetInSlideCoordinateSystem + z_offset = pp.get( + 'ZOffsetInSlideCoordinateSystem', + 0.0, + ) + origin_preserved = ( + x_offset == src_x_offset and + y_offset == src_y_offset and + z_offset == src_z_offset + ) + orientation = plane_orientation[0].ImageOrientationSlide + image_position = [x_offset, y_offset, z_offset] - plane_positions = [ - compute_plane_position_tiled_full( - row_index=r, - column_index=c, - x_offset=x_offset, - y_offset=y_offset, - rows=self.Rows, - columns=self.Columns, - image_orientation=orientation, - pixel_spacing=pixel_measures[0].PixelSpacing, + are_total_pixel_matrix_locations_preserved = ( + origin_preserved and + ( + not user_provided_orientation or + plane_orientation == source_plane_orientation + ) and + ( + not user_provided_measures or + pixel_measures == source_pixel_measures ) - for c, r in tile_pixel_matrix( - total_pixel_matrix_rows=pixel_array.shape[1], - total_pixel_matrix_columns=pixel_array.shape[2], - rows=self.Rows, - columns=self.Columns, + ) + + if are_total_pixel_matrix_locations_preserved: + if ( + pixel_array.shape[1:3] != + ( + src_img.TotalPixelMatrixRows, + src_img.TotalPixelMatrixColumns + ) + ): + raise ValueError( + "Shape of input pixel_array does not match shape " + "of the total pixel matrix of the source image." + ) + + # The overall total pixel matrix can match the source + # image's but if the image is tiled differently, spatial + # locations within each frame are not preserved + are_spatial_locations_preserved = ( + tile_size == (src_img.Rows, src_img.Columns) ) + else: + are_spatial_locations_preserved = False + + raw_plane_positions = compute_tile_positions_per_frame( + rows=self.Rows, + columns=self.Columns, + total_pixel_matrix_rows=pixel_array.shape[1], + total_pixel_matrix_columns=pixel_array.shape[2], + total_pixel_matrix_image_position=image_position, + image_orientation=orientation, + pixel_spacing=pixel_measures[0].PixelSpacing, + ) + plane_sort_index = np.arange(len(raw_plane_positions)) + + # Only need to create the plane position DICOM objects if + # they will be placed into the object. Otherwise skip this + # as it is really inefficient + if ( + dimension_organization_type != + DimensionOrganizationTypeValues.TILED_FULL + ): + plane_positions = [ + PlanePositionSequence( + CoordinateSystemNames.SLIDE, + image_position=coords, + pixel_matrix_position=offsets, + ) + for offsets, coords in raw_plane_positions + ] + else: + # Unneeded + plane_positions = [None] + + # Match the format used elsewhere + plane_position_values = np.array( + [ + [*offsets, *coords] + for offsets, coords in raw_plane_positions + ] + ) + + # compute_tile_positions_per_frame returns + # (c, r, x, y, z) but the dimension index sequence + # requires (r, c, x, y z). Swap here to correct for + # this + plane_position_values = plane_position_values[ + :, [1, 0, 2, 3, 4] ] else: + are_measures_and_orientation_preserved = ( + ( + not user_provided_orientation or + plane_orientation == source_plane_orientation + ) and + ( + not user_provided_measures or + pixel_measures == source_pixel_measures + ) + ) + + if ( + plane_positions is None or + are_measures_and_orientation_preserved + ): + # Calculating source positions can be slow, so avoid unless + # necessary + dim_ind = self.DimensionIndexSequence + if is_multiframe: + source_plane_positions = \ + dim_ind.get_plane_positions_of_image( + src_img + ) + else: + source_plane_positions = \ + dim_ind.get_plane_positions_of_series( + source_images + ) + if plane_positions is None: if pixel_array.shape[0] != len(source_plane_positions): raise ValueError( @@ -1809,6 +1971,8 @@ def __init__( '"pixel_array" argument.' ) plane_positions = source_plane_positions + are_spatial_locations_preserved = \ + are_measures_and_orientation_preserved else: if pixel_array.shape[0] != len(plane_positions): raise ValueError( @@ -1816,28 +1980,26 @@ def __init__( 'via "plane_positions" argument does not match ' 'size of first dimension of "pixel_array" argument.' ) - - # plane_position_values is an array giving, for each plane of - # the input array, the raw values of all attributes that - # describe its position. The first dimension is sorted the same - # way as the input pixel array and the second is sorted the - # same way as the dimension index sequence (without segment - # number) plane_sort_index is a list of indices into the input - # planes giving the order in which they should be arranged to - # correctly sort them for inclusion into the segmentation - plane_position_values, plane_sort_index = \ - self.DimensionIndexSequence.get_index_values( - plane_positions - ) - - are_spatial_locations_preserved = ( - all( - plane_positions[i] == source_plane_positions[i] - for i in range(len(plane_positions)) - ) and - plane_orientation == source_plane_orientation and - pixel_measures == source_pixel_measures - ) + if are_measures_and_orientation_preserved: + are_spatial_locations_preserved = all( + plane_positions[i] == source_plane_positions[i] + for i in range(len(plane_positions)) + ) + else: + are_spatial_locations_preserved = False + + # plane_position_values is an array giving, for each plane of + # the input array, the raw values of all attributes that + # describe its position. The first dimension is sorted the same + # way as the input pixel array and the second is sorted the + # same way as the dimension index sequence (without segment + # number) plane_sort_index is a list of indices into the input + # planes giving the order in which they should be arranged to + # correctly sort them for inclusion into the segmentation + plane_position_values, plane_sort_index = \ + self.DimensionIndexSequence.get_index_values( + plane_positions + ) else: # Only one spatial location supported @@ -1846,34 +2008,20 @@ def __init__( plane_sort_index = np.array([0]) are_spatial_locations_preserved = True - if are_spatial_locations_preserved: - if tile_pixel_array: - if ( - pixel_array.shape[1:3] != - ( - src_img.TotalPixelMatrixRows, - src_img.TotalPixelMatrixColumns - ) - ): - raise ValueError( - "Shape of input pixel_array does not match shape of " - "the total pixel matrix of the source image." - ) - else: - if pixel_array.shape[1:3] != (src_img.Rows, src_img.Columns): - raise ValueError( - "Shape of input pixel_array does not match shape of " - "the source image." - ) + if are_spatial_locations_preserved and not tile_pixel_array: + if pixel_array.shape[1:3] != (src_img.Rows, src_img.Columns): + raise ValueError( + "Shape of input pixel_array does not match shape of " + "the source image." + ) # Dimension Organization Type dimension_organization_type = self._check_dimension_organization_type( dimension_organization_type=dimension_organization_type, is_tiled=is_tiled, - are_spatial_locations_preserved=are_spatial_locations_preserved, omit_empty_frames=omit_empty_frames, - source_image=src_img, plane_positions=plane_positions, + tile_pixel_array=tile_pixel_array, rows=self.Rows, columns=self.Columns, ) @@ -1908,7 +2056,11 @@ def __init__( else: included_plane_indices = list(range(len(plane_positions))) - if has_ref_frame_uid: + if ( + has_ref_frame_uid and + dimension_organization_type != + DimensionOrganizationTypeValues.TILED_FULL + ): # Get unique values of attributes in the Plane Position Sequence or # Plane Position Slide Sequence, which define the position of the # plane with respect to the three dimensional patient or slide @@ -1942,6 +2094,16 @@ def __init__( total_pixel_matrix_size=total_pixel_matrix_size, ) + plane_position_names = ( + self.DimensionIndexSequence.get_index_keywords() + ) + row_dim_index = plane_position_names.index( + 'RowPositionInTotalImagePixelMatrix' + ) + col_dim_index = plane_position_names.index( + 'ColumnPositionInTotalImagePixelMatrix' + ) + is_encaps = self.file_meta.TransferSyntaxUID.is_encapsulated process_pool: Optional[Executor] = None @@ -1961,6 +2123,13 @@ def __init__( # at the end frames: Union[List[bytes], List[np.ndarray]] = [] + # In the case of native encoding when the number pixels in a frame is + # not a multiple of 8. This array carries "leftover" pixels that + # couldn't be encoded in previous iterations, to future iterations This + # saves having to keep the entire un-endoded array in memory, which can + # get extremely heavy on memory in the case of very large arrays + remainder_pixels = np.empty((0, ), dtype=np.uint8) + if is_encaps: if using_multiprocessing: # In the case of encapsulated transfer syntaxes with multiple @@ -2008,11 +2177,27 @@ def __init__( for plane_index in plane_sort_index: if tile_pixel_array: - pos = plane_positions[plane_index][0] + if ( + dimension_organization_type == + DimensionOrganizationTypeValues.TILED_FULL + ): + row_offset = int( + plane_position_values[plane_index, row_dim_index] + ) + column_offset = int( + plane_position_values[plane_index, col_dim_index] + ) + else: + pos = plane_positions[plane_index][0] + row_offset = pos.RowPositionInTotalImagePixelMatrix + column_offset = ( + pos.ColumnPositionInTotalImagePixelMatrix + ) + plane_array = get_tile_array( pixel_array[0], - row_offset=pos.RowPositionInTotalImagePixelMatrix, - column_offset=pos.ColumnPositionInTotalImagePixelMatrix, + row_offset=row_offset, + column_offset=column_offset, tile_rows=self.Rows, tile_columns=self.Columns, ) @@ -2045,32 +2230,34 @@ def __init__( f'add plane #{plane_index} for segment #{segment_number}' ) - # Get the item of the PerFrameFunctionalGroupsSequence for this - # segmentation frame - if has_ref_frame_uid: - plane_pos_val = plane_position_values[plane_index] - try: - dimension_index_values = ( - self._get_dimension_index_values( - unique_dimension_values=unique_dimension_values, - plane_position_value=plane_pos_val, - coordinate_system=self._coordinate_system, - ) - ) - except IndexError as error: - raise IndexError( - 'Could not determine position of plane ' - f'#{plane_index} in three dimensional coordinate ' - f'system based on dimension index values: {error}' - ) from error - else: - dimension_index_values = [] - if ( dimension_organization_type != DimensionOrganizationTypeValues.TILED_FULL ): # No per-frame functional group for TILED FULL + + # Get the item of the PerFrameFunctionalGroupsSequence for + # this segmentation frame + if has_ref_frame_uid: + plane_pos_val = plane_position_values[plane_index] + try: + dimension_index_values = ( + self._get_dimension_index_values( + unique_dimension_values=unique_dimension_values, # noqa: E501 + plane_position_value=plane_pos_val, + coordinate_system=self._coordinate_system, + ) + ) + except IndexError as error: + raise IndexError( + 'Could not determine position of plane ' + f'#{plane_index} in three dimensional ' + 'coordinate system based on dimension index ' + f'values: {error}' + ) from error + else: + dimension_index_values = [] + pffg_item = self._get_pffg_item( segment_number=segment_number, dimension_index_values=dimension_index_values, @@ -2104,8 +2291,25 @@ def __init__( ) frame_futures.append(future) else: - # Concatenate the 1D array for encoding at the end - frames.append(segment_array.flatten()) + flat_array = segment_array.flatten() + if ( + self.SegmentationType == + SegmentationTypeValues.BINARY.value and + (self.Rows * self.Columns) // 8 != 0 + ): + # Need to encode a multiple of 8 pixels at a time + full_array = np.concatenate( + [remainder_pixels, flat_array] + ) + # Round down to closest multiple of 8 + n_pixels_to_take = 8 * (len(full_array) // 8) + to_encode = full_array[:n_pixels_to_take] + remainder_pixels = full_array[n_pixels_to_take:] + else: + # Simple - each frame can be individually encoded + to_encode = flat_array + + frames.append(self._encode_pixels_native(to_encode)) if ( dimension_organization_type != @@ -2128,13 +2332,13 @@ def __init__( self.NumberOfFrames = len(frames) self.PixelData = encapsulate(frames) else: - # Encode the whole pixel array at once - # This allows for correct bit-packing in cases where - # number of pixels per frame is not a multiple of 8 self.NumberOfFrames = len(frames) - self.PixelData = self._encode_pixels_native( - np.concatenate(frames) - ) + + # May need to add in a final set of pixels + if len(remainder_pixels) > 0: + frames.append(self._encode_pixels_native(remainder_pixels)) + + self.PixelData = b''.join(frames) # Add a null trailing byte if required if len(self.PixelData) % 2 == 1: @@ -2338,6 +2542,8 @@ def _add_slide_coordinate_metadata( format_number_as_ds(x_origin) origin_item.YOffsetInSlideCoordinateSystem = \ format_number_as_ds(y_origin) + origin_item.ZOffsetInSlideCoordinateSystem = \ + format_number_as_ds(z_origin) self.TotalPixelMatrixOriginSequence = [origin_item] self.TotalPixelMatrixFocalPlanes = 1 if total_pixel_matrix_size is None: @@ -2355,7 +2561,9 @@ def _add_slide_coordinate_metadata( else: transform = ImageToReferenceTransformer( image_position=(x_origin, y_origin, z_origin), - image_orientation=plane_orientation, + image_orientation=( + plane_orientation[0].ImageOrientationSlide + ), pixel_spacing=pixel_measures[0].PixelSpacing ) center_image_coordinates = np.array( @@ -2385,10 +2593,9 @@ def _check_dimension_organization_type( None, ], is_tiled: bool, - are_spatial_locations_preserved: bool, omit_empty_frames: bool, - source_image: Dataset, plane_positions: Sequence[PlanePositionSequence], + tile_pixel_array: bool, rows: int, columns: int, ) -> Optional[DimensionOrganizationTypeValues]: @@ -2400,13 +2607,10 @@ def _check_dimension_organization_type( The specified DimensionOrganizationType for the output Segmentation. is_tiled: bool Whether the source image is a tiled image. - are_spatial_locations_preserved: bool - Whether spatial locations are preserved between the source image - and the segmentation pixel array. omit_empty_frames: bool Whether it was specified to omit empty frames. - source_image: pydicom.Dataset - Representative dataset of the source images. + tile_pixel_array: bool + Whether the total pixel matrix was passed. plane_positions: Sequence[highdicom.PlanePositionSequence] Plane positions of all frames. rows: int @@ -2448,10 +2652,15 @@ def _check_dimension_organization_type( dimension_organization_type == DimensionOrganizationTypeValues.TILED_FULL ): - if not are_plane_positions_tiled_full( - plane_positions, - rows, - columns, + # Need to check positions if they were not generated by us + # when using tile_pixel_array + if ( + not tile_pixel_array and + not are_plane_positions_tiled_full( + plane_positions, + rows, + columns, + ) ): raise ValueError( 'A value of "TILED_FULL" for parameter ' @@ -2897,66 +3106,140 @@ def _get_pffg_item( Per Frame Functional Groups Sequence for this segmentation frame. """ + # NB this function is called many times in a loop when there are a + # large number of frames, and has been observed to dominate the + # creation time of some segmentations. Therefore we use low-level + # pydicom primitives to improve performance as much as possible pffg_item = Dataset() frame_content_item = Dataset() - frame_content_item.DimensionIndexValues = ( - [int(segment_number)] + dimension_index_values + frame_content_item.add( + DataElement( + 0x00209157, # DimensionIndexValues + 'UL', + [int(segment_number)] + dimension_index_values + ) + ) + pffg_item.add( + DataElement( + 0x00209111, # FrameContentSequence + 'SQ', + [frame_content_item] + ) ) - pffg_item.FrameContentSequence = [frame_content_item] if has_ref_frame_uid: if coordinate_system == CoordinateSystemNames.SLIDE: - pffg_item.PlanePositionSlideSequence = plane_position + pffg_item.add( + DataElement( + 0x0048021a, # PlanePositionSlideSequence + 'SQ', + plane_position + ) + ) else: - pffg_item.PlanePositionSequence = plane_position - - # Determining the source images that map to the frame is not - # always trivial. Since DerivationImageSequence is a type 2 - # attribute, we leave its value empty. - pffg_item.DerivationImageSequence = [] + pffg_item.add( + DataElement( + 0x00209113, # PlanePositionSequence + 'SQ', + plane_position + ) + ) if are_spatial_locations_preserved: derivation_image_item = Dataset() - derivation_code = codes.cid7203.Segmentation - derivation_image_item.DerivationCodeSequence = [ - CodedConcept.from_code(derivation_code) - ] + derivation_image_item.add( + DataElement( + 0x00089215, # DerivationCodeSequence + 'SQ', + [_DERIVATION_CODE] + ) + ) derivation_src_img_item = Dataset() - if hasattr(source_images[0], 'NumberOfFrames'): + if 0x00280008 in source_images[0]: # NumberOfFrames # A single multi-frame source image src_img_item = source_images[0] # Frame numbers are one-based - derivation_src_img_item.ReferencedFrameNumber = ( - source_image_index + 1 + derivation_src_img_item.add( + DataElement( + 0x00081160, # ReferencedFrameNumber + 'IS', + source_image_index + 1 + ) ) else: # Multiple single-frame source images src_img_item = source_images[source_image_index] - derivation_src_img_item.ReferencedSOPClassUID = \ - src_img_item.SOPClassUID - derivation_src_img_item.ReferencedSOPInstanceUID = \ - src_img_item.SOPInstanceUID - purpose_code = \ - codes.cid7202.SourceImageForImageProcessingOperation - derivation_src_img_item.PurposeOfReferenceCodeSequence = [ - CodedConcept.from_code(purpose_code) - ] - derivation_src_img_item.SpatialLocationsPreserved = 'YES' - derivation_image_item.SourceImageSequence = [ - derivation_src_img_item, - ] - pffg_item.DerivationImageSequence.append( - derivation_image_item + derivation_src_img_item.add( + DataElement( + 0x00081150, # ReferencedSOPClassUID + 'UI', + src_img_item[0x00080016].value # SOPClassUID + ) + ) + derivation_src_img_item.add( + DataElement( + 0x00081155, # ReferencedSOPInstanceUID + 'UI', + src_img_item[0x00080018].value # SOPInstanceUID + ) + ) + derivation_src_img_item.add( + DataElement( + 0x0040a170, # PurposeOfReferenceCodeSequence + 'SQ', + [_PURPOSE_CODE] + ) + ) + derivation_src_img_item.add( + DataElement( + 0x0028135a, # SpatialLocationsPreserved + 'CS', + 'YES' + ) + ) + derivation_image_item.add( + DataElement( + 0x00082112, # SourceImageSequence + 'SQ', + [derivation_src_img_item] + ) + ) + pffg_item.add( + DataElement( + 0x00089124, # DerivationImageSequence + 'SQ', + [derivation_image_item] + ) ) else: + # Determining the source images that map to the frame is not + # always trivial. Since DerivationImageSequence is a type 2 + # attribute, we leave its value empty. + pffg_item.add( + DataElement( + 0x00089124, # DerivationImageSequence + 'SQ', + [] + ) + ) logger.debug('spatial locations not preserved') identification = Dataset() - identification.ReferencedSegmentNumber = int(segment_number) - pffg_item.SegmentIdentificationSequence = [ - identification, - ] + identification.add( + DataElement( + 0x0062000b, # ReferencedSegmentNumber + 'US', + int(segment_number) + ) + ) + pffg_item.add( + DataElement( + 0x0062000a, # SegmentIdentificationSequence + 'SQ', + [identification] + ) + ) return pffg_item @@ -2977,7 +3260,7 @@ def _encode_pixels_native(self, planes: np.ndarray) -> bytes: """ if self.SegmentationType == SegmentationTypeValues.BINARY.value: - return pack_bits(planes) + return pack_bits(planes, pad=False) else: return planes.tobytes() diff --git a/src/highdicom/spatial.py b/src/highdicom/spatial.py index d5fab854..de4a7557 100644 --- a/src/highdicom/spatial.py +++ b/src/highdicom/spatial.py @@ -1,7 +1,648 @@ -from typing import Sequence, Tuple +import itertools +from typing import Generator, Iterator, List, Optional, Sequence, Tuple +from pydicom import Dataset import numpy as np +from highdicom.enum import CoordinateSystemNames +from highdicom._module_utils import is_multiframe_image + + +# Tolerance value used by default in tests for equality +_DEFAULT_TOLERANCE = 1e-5 + + +def is_tiled_image(dataset: Dataset) -> bool: + """Determine whether a dataset represents a tiled image. + + Returns + ------- + bool: + True if the dataset is a tiled image. False otherwise. + + """ + if ( + hasattr(dataset, 'TotalPixelMatrixRows') and + hasattr(dataset, 'TotalPixelMatrixColumns') and + hasattr(dataset, 'NumberOfFrames') + ): + return True + return False + + +def tile_pixel_matrix( + total_pixel_matrix_rows: int, + total_pixel_matrix_columns: int, + rows: int, + columns: int, +) -> Iterator[Tuple[int, int]]: + """Tiles an image into smaller frames (rectangular regions). + + Follows the convention used in image with Dimension Organization Type + "TILED_FULL" images. + + Parameters + ---------- + total_pixel_matrix_rows: int + Number of rows in the Total Pixel Matrix + total_pixel_matrix_columns: int + Number of columns in the Total Pixel Matrix + rows: int + Number of rows per Frame (tile) + columns: int + Number of columns per Frame (tile) + + Returns + ------- + Iterator + One-based (Column, Row) index of each Frame (tile) + + """ + tiles_per_col = int(np.ceil(total_pixel_matrix_rows / rows)) + tiles_per_row = int(np.ceil(total_pixel_matrix_columns / columns)) + tile_row_indices = iter(range(1, tiles_per_col + 1)) + tile_col_indices = iter(range(1, tiles_per_row + 1)) + return ( + (c, r) for (r, c) in itertools.product( + tile_row_indices, + tile_col_indices + ) + ) + + +def get_tile_array( + pixel_array: np.ndarray, + row_offset: int, + column_offset: int, + tile_rows: int, + tile_columns: int, + pad: bool = True, +) -> np.ndarray: + """Extract a tile from a total pixel matrix array. + + Parameters + ---------- + pixel_array: np.ndarray + Array representing a total pixel matrix. The first two dimensions + are treated as the rows and columns, respectively, of the total pixel + matrix. Any subsequent dimensions are not used but are retained in the + output array. + row_offset: int + Offset of the first row of the requested tile from the top of the total + pixel matrix (1-based index). + column_offset: int + Offset of the first column of the requested tile from the left of the + total pixel matrix (1-based index). + tile_rows: int + Number of rows per tile. + tile_columns: + Number of columns per tile. + pad: bool + Whether to pad the returned array with zeros at the right and/or bottom + to ensure that it matches the correct tile size. Otherwise, the returned + array is not padded and may be smaller than the full tile size. + + Returns + ------- + np.ndarray: + Returned pixel array for the requested tile. + + """ + if row_offset < 1 or row_offset > pixel_array.shape[0]: + raise ValueError( + "Row offset must be between 1 and the size of dimension 0 of the " + "pixel array." + ) + if column_offset < 1 or column_offset > pixel_array.shape[1]: + raise ValueError( + "Column offset must be between 1 and the size of dimension 1 of " + "the pixel array." + ) + # Move to pythonic 1-based indexing + row_offset -= 1 + column_offset -= 1 + row_end = row_offset + tile_rows + if row_end > pixel_array.shape[0]: + pad_rows = row_end - pixel_array.shape[0] + row_end = pixel_array.shape[0] + else: + pad_rows = 0 + column_end = column_offset + tile_columns + if column_end > pixel_array.shape[1]: + pad_columns = column_end - pixel_array.shape[1] + column_end = pixel_array.shape[1] + else: + pad_columns = 0 + # Account for 1-based to 0-based index conversion + tile_array = pixel_array[row_offset:row_end, column_offset:column_end] + if pad and (pad_rows > 0 or pad_columns > 0): + extra_dims = pixel_array.ndim - 2 + padding = [(0, pad_rows), (0, pad_columns)] + [(0, 0)] * extra_dims + tile_array = np.pad(tile_array, padding) + + return tile_array + + +def compute_tile_positions_per_frame( + rows: int, + columns: int, + total_pixel_matrix_rows: int, + total_pixel_matrix_columns: int, + total_pixel_matrix_image_position: Sequence[float], + image_orientation: Sequence[float], + pixel_spacing: Sequence[float], +) -> List[Tuple[List[int], List[float]]]: + """Get positions of each tile in a TILED_FULL image. + + A TILED_FULL image is one with DimensionOrganizationType of "TILED_FULL". + + Parameters + ---------- + rows: int + Number of rows per tile. + columns: int + Number of columns per tile. + total_pixel_matrix_rows: int + Number of rows in the total pixel matrix. + total_pixel_matrix_columns: int + Number of columns in the total pixel matrix. + total_pixel_matrix_image_position: Sequence[float] + Position of the top left pixel of the total pixel matrix in the frame + of reference. Sequence of length 3. + image_orientation: Sequence[float] + Orientation cosines of the total pixel matrix. Sequence of length 6. + pixel_spacing: Sequence[float] + Pixel spacing between the (row, columns) in mm. Sequence of length 2. + + Returns + ------- + List[Tuple[List[int], List[float]]]: + List with positions for each of the tiles in the tiled image. The + first tuple contains the (column offset, row offset) values, which + are one-based offsets of the tile in pixel units from the top left + of the total pixel matrix. The second tuple contains the image + position in the frame of reference for the tile. + + """ + if len(total_pixel_matrix_image_position) != 3: + raise ValueError( + "Argument 'total_pixel_matrix_image_position' must have length 3." + ) + if len(image_orientation) != 6: + raise ValueError( + "Argument 'image_orientation' must have length 6." + ) + if len(pixel_spacing) != 2: + raise ValueError( + "Argument 'pixel_spacing' must have length 2." + ) + + tiles_per_column = ( + (total_pixel_matrix_columns - 1) // columns + 1 + ) + tiles_per_row = (total_pixel_matrix_rows - 1) // rows + 1 + + # N x 2 array of (c, r) tile indices + tile_indices = np.stack( + np.meshgrid( + range(tiles_per_column), + range(tiles_per_row), + indexing='xy', + ) + ).reshape(2, -1).T + + # N x 2 array of (c, r) pixel indices + pixel_indices = tile_indices * [columns, rows] + + transformer = PixelToReferenceTransformer( + image_position=total_pixel_matrix_image_position, + image_orientation=image_orientation, + pixel_spacing=pixel_spacing, + ) + image_positions = transformer(pixel_indices) + + # Convert 0-based to 1-based indexing for the output + pixel_indices += 1 + + return list( + zip( + pixel_indices.tolist(), + image_positions.tolist() + ) + ) + + +def iter_tiled_full_frame_data( + dataset: Dataset, +) -> Generator[Tuple[int, int, int, int, float, float, float], None, None]: + """Get data on the position of each tile in a TILED_FULL image. + + This works only with images with Dimension Organization Type of + "TILED_FULL". + + Unlike :func:`highdicom.utils.compute_plane_position_slide_per_frame`, + this functions returns the data in their basic Python types rather than + wrapping as :class:`highdicom.PlanePositionSequence` + + Parameters + ---------- + dataset: pydicom.dataset.Dataset + VL Whole Slide Microscopy Image or Segmentation Image using the + "TILED_FULL" DimensionOrganizationType. + + Returns + ------- + channel: int + 1-based integer index of the "channel". The meaning of "channel" + depends on the image type. For segmentation images, the channel is the + segment number. For other images, it is the optical path number. + focal_plane_index: int + 1-based integer index of the focal plane. + column_position: int + 1-based column position of the tile (measured left from the left side + of the total pixel matrix). + row_position: int + 1-based row position of the tile (measured down from the top of the + total pixel matrix). + x: float + X coordinate in the frame-of-reference coordinate system in millimeter + units. + y: float + Y coordinate in the frame-of-reference coordinate system in millimeter + units. + z: float + Z coordinate in the frame-of-reference coordinate system in millimeter + units. + + """ + allowed_sop_class_uids = { + '1.2.840.10008.5.1.4.1.1.77.1.6', # VL Whole Slide Microscopy Image + '1.2.840.10008.5.1.4.1.1.66.4', # Segmentation Image + } + if dataset.SOPClassUID not in allowed_sop_class_uids: + raise ValueError( + 'Expected a VL Whole Slide Microscopy Image or Segmentation Image.' + ) + if ( + not hasattr(dataset, "DimensionOrganizationType") or + dataset.DimensionOrganizationType != "TILED_FULL" + ): + raise ValueError( + 'Expected an image with "TILED_FULL" dimension organization type.' + ) + + image_origin = dataset.TotalPixelMatrixOriginSequence[0] + image_orientation = ( + float(dataset.ImageOrientationSlide[0]), + float(dataset.ImageOrientationSlide[1]), + float(dataset.ImageOrientationSlide[2]), + float(dataset.ImageOrientationSlide[3]), + float(dataset.ImageOrientationSlide[4]), + float(dataset.ImageOrientationSlide[5]), + ) + num_focal_planes = getattr( + dataset, + 'TotalPixelMatrixFocalPlanes', + 1 + ) + + is_segmentation = dataset.SOPClassUID == '1.2.840.10008.5.1.4.1.1.66.4' + + # The "channels" output is either segment for segmentations, or optical + # path for other images + if is_segmentation: + num_channels = len(dataset.SegmentSequence) + else: + num_channels = getattr( + dataset, + 'NumberOfOpticalPaths', + len(dataset.OpticalPathSequence) + ) + + shared_fg = dataset.SharedFunctionalGroupsSequence[0] + pixel_measures = shared_fg.PixelMeasuresSequence[0] + pixel_spacing = ( + float(pixel_measures.PixelSpacing[0]), + float(pixel_measures.PixelSpacing[1]), + ) + spacing_between_slices = float( + getattr( + pixel_measures, + 'SpacingBetweenSlices', + 1.0 + ) + ) + x_offset = image_origin.XOffsetInSlideCoordinateSystem + y_offset = image_origin.YOffsetInSlideCoordinateSystem + + for channel in range(1, num_channels + 1): + for slice_index in range(1, num_focal_planes + 1): + z_offset = float(slice_index - 1) * spacing_between_slices + + for offsets, coords in compute_tile_positions_per_frame( + rows=dataset.Rows, + columns=dataset.Columns, + total_pixel_matrix_rows=dataset.TotalPixelMatrixRows, + total_pixel_matrix_columns=dataset.TotalPixelMatrixColumns, + total_pixel_matrix_image_position=( + x_offset, y_offset, z_offset + ), + image_orientation=image_orientation, + pixel_spacing=pixel_spacing + ): + yield ( + channel, + slice_index, + int(offsets[0]), + int(offsets[1]), + float(coords[0]), + float(coords[1]), + float(coords[2]), + ) + + +def get_image_coordinate_system( + dataset: Dataset +) -> Optional[CoordinateSystemNames]: + """Get the coordinate system used by an image. + + Parameters + ---------- + dataset: pydicom.Dataset + Dataset representing an image. + + Returns + -------- + Union[highdicom.enum.CoordinateSystemNames, None] + Coordinate system used by the image, if any. + + """ + if not hasattr(dataset, "FrameOfReferenceUID"): + return None + + # Using Container Type Code Sequence attribute would be more + # elegant, but unfortunately it is a type 2 attribute. + if ( + hasattr(dataset, 'ImageOrientationSlide') or + hasattr(dataset, 'ImageCenterPointCoordinatesSequence') + ): + return CoordinateSystemNames.SLIDE + else: + return CoordinateSystemNames.PATIENT + + +def _get_spatial_information( + dataset: Dataset, + frame_number: Optional[int] = None, + for_total_pixel_matrix: bool = False, +) -> Tuple[ + List[float], + List[float], + List[float], + Optional[float], +]: + """Get spatial information from an image dataset. + + Parameters + ---------- + dataset: pydicom.Dataset + Dataset representing an image. + frame_number: Union[int, None], optional + Specific 1-based frame number. Required if dataset is a multi-frame + image. Should be None otherwise. + for_total_pixel_matrix: bool, optional + If True, get spatial information for the total pixel matrix of a tiled + image. This should only be True if the image is a tiled image and is + incompatible with specifying a frame number. + + Returns + ------- + image_position: List[float] + Image position (3 elements) of the dataset or frame. + image_orientation: List[float] + Image orientation (6 elements) of the dataset or frame. + pixel_spacing: List[float] + Pixel spacing (2 elements) of the dataset or frame. + spacing_between_slices: Union[float, None] + Spacing between adjacent slices, if found in the dataset. Note that + there is no guarantee in general that slices will be consistently + spaced or even parallel. This function does not attempt to calculate + spacing if it is not found in the dataset. + + """ + coordinate_system = get_image_coordinate_system(dataset) + + if coordinate_system is None: + raise ValueError( + 'The input "dataset" has no spatial information ' + 'as it has no frame of reference.' + ) + + if for_total_pixel_matrix: + if not hasattr(dataset, 'TotalPixelMatrixOriginSequence'): + raise ValueError('Image is not a tiled image.') + origin_seq = dataset.TotalPixelMatrixOriginSequence[0] + position = ( + origin_seq.XOffsetInSlideCoordinateSystem, + origin_seq.YOffsetInSlideCoordinateSystem, + getattr(origin_seq, 'ZOffsetInSlideCoordinateSystem', 0.0) + ) + shared_seq = dataset.SharedFunctionalGroupsSequence[0] + if hasattr(shared_seq, 'PixelMeasuresSequence'): + pixel_spacing = shared_seq.PixelMeasuresSequence[0].PixelSpacing + spacing_between_slices = getattr( + shared_seq.PixelMeasuresSequence[0], + "SpacingBetweenSlices", + None, + ) + else: + raise ValueError( + "PixelMeasuresSequence not found in the " + "SharedFunctionalGroupsSequence." + ) + + orientation = dataset.ImageOrientationSlide + return position, orientation, pixel_spacing, spacing_between_slices + + if is_multiframe_image(dataset): + if frame_number is None: + raise TypeError( + 'Argument "frame_number" must be specified for a multi-frame ' + 'image.' + ) + shared_seq = dataset.SharedFunctionalGroupsSequence[0] + is_tiled_full = ( + dataset.get("DimensionOrganizationType", "") == "TILED_FULL" + ) + if is_tiled_full: + frame_seq = None + else: + frame_seq = dataset.PerFrameFunctionalGroupsSequence[ + frame_number - 1 + ] + + # Find spacing in either shared or per-frame sequences (this logic is + # the same for patient or slide coordinate system) + if hasattr(shared_seq, 'PixelMeasuresSequence'): + pixel_measures = shared_seq.PixelMeasuresSequence[0] + elif ( + frame_seq is not None and + hasattr(frame_seq, 'PixelMeasuresSequence') + ): + pixel_measures = frame_seq.PixelMeasuresSequence[0] + else: + raise ValueError('No pixel measures information found.') + pixel_spacing = pixel_measures.PixelSpacing + spacing_between_slices = getattr( + pixel_measures, + 'SpacingBetweenSlices', + None, + ) + + if coordinate_system == CoordinateSystemNames.SLIDE: + if is_tiled_full: + # TODO this iteration is probably rather inefficient + _, _, _, _, *position = next( + itertools.islice( + iter_tiled_full_frame_data(dataset), + frame_number - 1, + frame_number, + ) + ) + else: + # Find position in either shared or per-frame sequences + if hasattr(shared_seq, 'PlanePositionSlideSequence'): + pos_seq = shared_seq.PlanePositionSlideSequence[0] + elif ( + frame_seq is not None and + hasattr(frame_seq, 'PlanePositionSlideSequence') + ): + pos_seq = frame_seq.PlanePositionSlideSequence[0] + else: + raise ValueError('No frame position information found.') + + position = [ + pos_seq.XOffsetInSlideCoordinateSystem, + pos_seq.YOffsetInSlideCoordinateSystem, + pos_seq.ZOffsetInSlideCoordinateSystem, + ] + + orientation = dataset.ImageOrientationSlide + + else: # PATIENT coordinate system (multiframe) + # Find position in either shared or per-frame sequences + if hasattr(shared_seq, 'PlanePositionSequence'): + pos_seq = shared_seq.PlanePositionSequence[0] + elif ( + frame_seq is not None and + hasattr(frame_seq, 'PlanePositionSequence') + ): + pos_seq = frame_seq.PlanePositionSequence[0] + else: + raise ValueError('No frame position information found.') + + position = pos_seq.ImagePositionPatient + + # Find orientation in either shared or per-frame sequences + if hasattr(shared_seq, 'PlaneOrientationSequence'): + pos_seq = shared_seq.PlaneOrientationSequence[0] + elif ( + frame_seq is not None and + hasattr(frame_seq, 'PlaneOrientationSequence') + ): + pos_seq = frame_seq.PlaneOrientationSequence[0] + else: + raise ValueError('No frame orientation information found.') + + orientation = pos_seq.ImageOrientationPatient + + else: # Single-frame image + if frame_number is not None and frame_number != 1: + raise TypeError( + 'Argument "frame_number" must be None or 1 for a single-frame ' + "image." + ) + position = dataset.ImagePositionPatient + orientation = dataset.ImageOrientationPatient + pixel_spacing = dataset.PixelSpacing + spacing_between_slices = getattr( + dataset, + "SpacingBetweenSlices", + None, + ) + + return position, orientation, pixel_spacing, spacing_between_slices + + +def _get_normal_vector(image_orientation: Sequence[float]) -> np.ndarray: + """Get normal vector given image cosines. + + Parameters + ---------- + image_orientation: Sequence[float] + Row and column cosines (6 element list) giving the orientation of the + image. + + Returns + ------- + np.ndarray + Array of shape (3, ) giving the normal vector to the image plane. + + """ + row_cosines = np.array(image_orientation[:3], dtype=float) + column_cosines = np.array(image_orientation[3:], dtype=float) + n = np.cross(row_cosines.T, column_cosines.T) + return n + + +def _are_images_coplanar( + image_position_a: Sequence[float], + image_orientation_a: Sequence[float], + image_position_b: Sequence[float], + image_orientation_b: Sequence[float], + tol: float = _DEFAULT_TOLERANCE, +) -> bool: + """Determine whether two images or image frames are coplanar. + + Two images are coplanar in the frame of reference coordinate system if and + only if their vectors have the same (or opposite direction) and the + shortest distance from the plane to the coordinate system origin is + the same for both planes. + + Parameters + ---------- + image_position_a: Sequence[float] + Image position (3 element list) giving the position of the center of + the top left pixel of the first image. + image_orientation_a: Sequence[float] + Row and column cosines (6 element list) giving the orientation of the + first image. + image_position_b: Sequence[float] + Image position (3 element list) giving the position of the center of + the top left pixel of the second image. + image_orientation_b: Sequence[float] + Row and column cosines (6 element list) giving the orientation of the + second image. + tol: float + Tolerance to use to determine equality. + + Returns + ------- + bool + True if the two images are coplanar. False otherwise. + + """ + n_a = _get_normal_vector(image_orientation_a) + n_b = _get_normal_vector(image_orientation_b) + if 1.0 - np.abs(n_a @ n_b) > tol: + return False + + # Find distances of both planes along n_a + dis_a = np.array(image_position_a, dtype=float) @ n_a + dis_b = np.array(image_position_b, dtype=float) @ n_a + + return abs(dis_a - dis_b) < tol + def create_rotation_matrix( image_orientation: Sequence[float], @@ -101,7 +742,7 @@ def _create_affine_transformation_matrix( rotation[:, 1] *= column_spacing # 4x4 transformation matrix - return np.row_stack( + return np.vstack( [ np.column_stack([ rotation, @@ -187,7 +828,7 @@ def _create_inv_affine_transformation_matrix( inv_rotation = np.linalg.inv(rotation) # 4x4 transformation matrix - return np.row_stack( + return np.vstack( [ np.column_stack([ inv_rotation, @@ -202,8 +843,8 @@ class PixelToReferenceTransformer: """Class for transforming pixel indices to reference coordinates. - This class facilitates the mapping of pixel indices to the pixel matrix of - an image or an image frame (tile or plane) into the patient or slide + This class facilitates the mapping of pixel indices of the pixel matrix + of an image or an image frame (tile or plane) into the patient or slide coordinate system defined by the frame of reference. Pixel indices are (column, row) pairs of zero-based integer values, where @@ -329,7 +970,7 @@ def __call__(self, indices: np.ndarray) -> np.ndarray: 'Argument "indices" must be a two-dimensional array ' 'of integers.' ) - pixel_matrix_coordinates = np.row_stack([ + pixel_matrix_coordinates = np.vstack([ indices.T.astype(float), np.zeros((indices.shape[0], ), dtype=float), np.ones((indices.shape[0], ), dtype=float), @@ -337,6 +978,47 @@ def __call__(self, indices: np.ndarray) -> np.ndarray: reference_coordinates = np.dot(self._affine, pixel_matrix_coordinates) return reference_coordinates[:3, :].T + @classmethod + def for_image( + cls, + dataset: Dataset, + frame_number: Optional[int] = None, + for_total_pixel_matrix: bool = False, + ) -> 'PixelToReferenceTransformer': + """Construct a transformer for a given image or image frame. + + Parameters + ---------- + dataset: pydicom.Dataset + Dataset representing an image. + frame_number: Union[int, None], optional + Frame number (using 1-based indexing) of the frame for which to get + the transformer. This should be provided if and only if the dataset + is a multi-frame image. + for_total_pixel_matrix: bool, optional + If True, use the spatial information for the total pixel matrix of + a tiled image. The result will be a transformer that maps pixel + indices of the total pixel matrix to frame of reference + coordinates. This should only be True if the image is a tiled image + and is incompatible with specifying a frame number. + + Returns + ------- + highdicom.spatial.PixelToReferenceTransformer: + Transformer object for the given image, or image frame. + + """ + position, orientation, spacing, _ = _get_spatial_information( + dataset, + frame_number=frame_number, + for_total_pixel_matrix=for_total_pixel_matrix, + ) + return cls( + image_position=position, + image_orientation=orientation, + pixel_spacing=spacing, + ) + class ReferenceToPixelTransformer: @@ -352,7 +1034,9 @@ class ReferenceToPixelTransformer: Pixel indices are (column, row) pairs of zero-based integer values, where the (0, 0) index is located at the **center** of the top left hand corner - pixel of the pixel matrix. + pixel of the pixel matrix. The result of the transform also contains a + third coordinate, giving position along the normal vector of the imaging + plane. Examples -------- @@ -382,7 +1066,9 @@ def __init__( image_position: Sequence[float], image_orientation: Sequence[float], pixel_spacing: Sequence[float], - spacing_between_slices: float = 1.0 + spacing_between_slices: float = 1.0, + round_output: bool = True, + drop_slice_index: bool = False, ): """Construct transformation object. @@ -412,6 +1098,15 @@ def __init__( spacing_between_slices: float, optional Distance (in the coordinate defined by the frame of reference) between neighboring slices. Default: 1 + round_output: bool, optional + If True, outputs are rounded to the nearest integer. Otherwise, + they are returned as float. + drop_slice_index: bool, optional + Whether to remove the 3rd column of the output array + (representing the out-of-plane coordinate) and return a 2D output + array. If this option is taken, and the resulting coordinates + do not lie in the range -0.5 to 0.5, a ``RuntimeError`` will be + triggered. Raises ------ @@ -423,6 +1118,8 @@ def __init__( an incorrect length. """ + self._round_output = round_output + self._drop_slice_index = drop_slice_index self._affine = _create_inv_affine_transformation_matrix( image_position=image_position, image_orientation=image_orientation, @@ -450,12 +1147,17 @@ def __call__(self, coordinates: np.ndarray) -> np.ndarray: Returns ------- numpy.ndarray - Array of (column, row) zero-based indices at pixel resolution. - Array of integer values with shape ``(n, 2)``, where *n* is - the number of indices, the first column represents the `column` - index and the second column represents the `row` index. - The ``(0, 0)`` coordinate is located at the **center** of the top - left pixel in the total pixel matrix. + Array of (column, row, slice) zero-based indices at pixel + resolution. Array of integer or floating point values with shape + ``(n, 3)``, where *n* is the number of indices, the first column + represents the `column` index, the second column represents the + `row` index, and the third column represents the `slice` coordinate + in the direction normal to the image plane (with scale given by the + ``spacing_between_slices_to`` parameter). The ``(0, 0, 0)`` + coordinate is located at the **center** of the top left pixel in + the total pixel matrix. The datatype of the array will be integer + if ``round_output`` is True (the default), or float if + ``round_output`` is False. Note ---- @@ -473,12 +1175,337 @@ def __call__(self, coordinates: np.ndarray) -> np.ndarray: 'Argument "coordinates" must be a two-dimensional array ' 'with shape [n, 3].' ) - reference_coordinates = np.row_stack([ + reference_coordinates = np.vstack([ coordinates.T.astype(float), np.ones((coordinates.shape[0], ), dtype=float) ]) pixel_matrix_coordinates = np.dot(self._affine, reference_coordinates) - return np.around(pixel_matrix_coordinates[:3, :].T).astype(int) + pixel_matrix_coordinates = pixel_matrix_coordinates[:3, :].T + if self._drop_slice_index: + if np.abs(pixel_matrix_coordinates[:, 2]).max() > 0.5: + raise RuntimeError( + "Output indices do not lie within the given image " + "plane." + ) + pixel_matrix_coordinates = pixel_matrix_coordinates[:, :2] + if self._round_output: + return np.around(pixel_matrix_coordinates).astype(int) + else: + return pixel_matrix_coordinates + + @classmethod + def for_image( + cls, + dataset: Dataset, + frame_number: Optional[int] = None, + for_total_pixel_matrix: bool = False, + round_output: bool = True, + drop_slice_index: bool = False, + ) -> 'ReferenceToPixelTransformer': + """Construct a transformer for a given image or image frame. + + Parameters + ---------- + dataset: pydicom.Dataset + Dataset representing an image. + frame_number: Union[int, None], optional + Frame number (using 1-based indexing) of the frame for which to get + the transformer. This should be provided if and only if the dataset + is a multi-frame image. + for_total_pixel_matrix: bool, optional + If True, use the spatial information for the total pixel matrix of + a tiled image. The result will be a transformer that maps frame of + reference coordinates to indices of the total pixel matrix. This + should only be True if the image is a tiled image and is + incompatible with specifying a frame number. + round_output: bool, optional + If True, outputs are rounded to the nearest integer. Otherwise, + they are returned as float. + drop_slice_index: bool, optional + Whether to remove the 3rd element of the output array + (representing the out-of-plane coordinate) and return a 2D output + array. If this option is taken, and the resulting coordinates + do not lie in the range -0.5 to 0.5, a ``RuntimeError`` will be + triggered. + + Returns + ------- + highdicom.spatial.ReferenceToPixelTransformer: + Transformer object for the given image, or image frame. + + """ + ( + position, + orientation, + spacing, + slice_spacing, + ) = _get_spatial_information( + dataset, + frame_number=frame_number, + for_total_pixel_matrix=for_total_pixel_matrix, + ) + if slice_spacing is None: + slice_spacing = 1.0 + return cls( + image_position=position, + image_orientation=orientation, + pixel_spacing=spacing, + spacing_between_slices=slice_spacing, + round_output=round_output, + drop_slice_index=drop_slice_index, + ) + + +class PixelToPixelTransformer: + + """Class for transforming pixel indices between two images. + + This class facilitates the mapping of pixel indices of the pixel matrix of + an image or an image frame (tile or plane) into those of another image or + image frame in the same frame of reference. This can include (but is not + limited) to mapping between different frames of the same image, or + different images within the same series (e.g. two levels of a spatial + pyramid). However, it is required that the two images be coplanar + within the frame-of-reference coordinate system. + + Pixel indices are (column, row) pairs of zero-based integer values, where + the (0, 0) index is located at the **center** of the top left hand corner + pixel of the pixel matrix. + + Examples + -------- + + Create a transformer for two images, where the second image has an axis + flipped relative to the first. + + >>> transformer = PixelToPixelTransformer( + ... image_position_from=[0.0, 0.0, 0.0], + ... image_orientation_from=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + ... pixel_spacing_from=[1.0, 1.0], + ... image_position_to=[0.0, 100.0, 0.0], + ... image_orientation_to=[1.0, 0.0, 0.0, 0.0, -1.0, 0.0], + ... pixel_spacing_to=[1.0, 1.0], + ... ) + + >>> indices_in = np.array([[0, 0], [50, 50]]) + >>> indices_out = transformer(indices_in) + >>> print(indices_out) + [[ 0 100] + [ 50 50]] + + Warning + ------- + This class shall not be used to map spatial coordinates (SCOORD) + between images. Use the + :class:`highdicom.spatial.ImageToImageTransformer` class instead. + + """ + + def __init__( + self, + image_position_from: Sequence[float], + image_orientation_from: Sequence[float], + pixel_spacing_from: Sequence[float], + image_position_to: Sequence[float], + image_orientation_to: Sequence[float], + pixel_spacing_to: Sequence[float], + round_output: bool = True, + ): + """Construct transformation object. + + The resulting object will map pixel indices of the "from" image to + pixel indices of the "to" image. + + Parameters + ---------- + image_position_from: Sequence[float] + Position of the "from" image in the frame of reference, + i.e., the offset of the top left hand corner pixel in the pixel + matrix from the origin of the reference coordinate system along the + X, Y, and Z axis + image_orientation_from: Sequence[float] + Cosines of the row direction (first triplet: horizontal, left to + right, increasing column index) and the column direction (second + triplet: vertical, top to bottom, increasing row index) direction + of the "from" image expressed in the three-dimensional patient or + slide coordinate system defined by the frame of reference + pixel_spacing_from: Sequence[float] + Spacing between pixels of the "from" imagem in millimeter unit + along the column direction (first value: spacing between rows, + vertical, top to bottom, increasing row index) and the rows + direction (second value: spacing between columns: horizontal, left + to right, increasing column index) + image_position_to: Sequence[float] + Position of the "to" image using the same definition as the "from" + image. + image_orientation_to: Sequence[float] + Orientation cosines of the "to" image using the same definition as + the "from" image. + pixel_spacing_to: Sequence[float] + Pixel spacing of the "to" image using the same definition as + the "from" image. + round_output: bool, optional + If True, outputs are rounded to the nearest integer. Otherwise, + they are returned as float. + + Raises + ------ + TypeError + When any of the arguments is not a sequence. + ValueError + When any of the arguments has an incorrect length, or if the two + images are not coplanar in the frame of reference coordinate + system. + + """ + self._round_output = round_output + if not _are_images_coplanar( + image_position_a=image_position_from, + image_orientation_a=image_orientation_from, + image_position_b=image_position_to, + image_orientation_b=image_orientation_to, + ): + raise ValueError( + "To two images do not exist within the same plane " + "in the frame of reference. and therefore pixel-to-pixel " + "transformation is not possible." + ) + pix_to_ref = _create_affine_transformation_matrix( + image_position=image_position_from, + image_orientation=image_orientation_from, + pixel_spacing=pixel_spacing_from, + ) + ref_to_pix = _create_inv_affine_transformation_matrix( + image_position=image_position_to, + image_orientation=image_orientation_to, + pixel_spacing=pixel_spacing_to, + ) + self._affine = np.dot(ref_to_pix, pix_to_ref) + + @property + def affine(self) -> np.ndarray: + """numpy.ndarray: 4x4 affine transformation matrix""" + return self._affine + + def __call__(self, indices: np.ndarray) -> np.ndarray: + """Transform pixel indices between two images. + + Parameters + ---------- + indices: numpy.ndarray + Array of (column, row) zero-based pixel indices of the "from" image + in the range [0, Columns - 1] and [0, Rows - 1], respectively. + Array of integer values with shape ``(n, 2)``, where *n* is the + number of indices, the first column represents the `column` index + and the second column represents the `row` index. The ``(0, 0)`` + coordinate is located at the **center** of the top left pixel in + the total pixel matrix. + + Returns + ------- + numpy.ndarray + Array of (column, row) zero-based pixel indices of the "to" image. + + Raises + ------ + ValueError + When `indices` has incorrect shape. + TypeError + When `indices` don't have integer data type. + + """ + if indices.shape[1] != 2: + raise ValueError( + 'Argument "indices" must be a two-dimensional array ' + 'with shape [n, 2].' + ) + if indices.dtype.kind not in ('u', 'i'): + raise TypeError( + 'Argument "indices" must be a two-dimensional array ' + 'of integers.' + ) + pixel_matrix_coordinates = np.vstack([ + indices.T.astype(float), + np.zeros((indices.shape[0], ), dtype=float), + np.ones((indices.shape[0], ), dtype=float), + ]) + output_coordinates = np.dot(self._affine, pixel_matrix_coordinates) + output_coordinates = output_coordinates[:2, :].T + if self._round_output: + return np.around(output_coordinates).astype(int) + else: + return output_coordinates + + @classmethod + def for_images( + cls, + dataset_from: Dataset, + dataset_to: Dataset, + frame_number_from: Optional[int] = None, + frame_number_to: Optional[int] = None, + for_total_pixel_matrix_from: bool = False, + for_total_pixel_matrix_to: bool = False, + round_output: bool = True, + ) -> 'PixelToPixelTransformer': + """Construct a transformer for two given images or image frames. + + Parameters + ---------- + dataset: pydicom.Dataset + Dataset representing an image. + frame_number: Union[int, None], optional + Frame number (using 1-based indexing) of the frame for which to get + the transformer. This should be provided if and only if the dataset + is a multi-frame image. + for_total_pixel_matrix: bool, optional + If True, use the spatial information for the total pixel matrix of + a tiled image. The result will be a transformer that maps pixel + indices of the total pixel matrix to frame of reference + coordinates. This should only be True if the image is a tiled image + and is incompatible with specifying a frame number. + round_output: bool, optional + If True, outputs are rounded to the nearest integer. Otherwise, + they are returned as float. + + Returns + ------- + highdicom.spatial.PixelToPixelTransformer: + Transformer object for the given image, or image frame. + + """ + if ( + not hasattr(dataset_from, 'FrameOfReferenceUID') or + not hasattr(dataset_to, 'FrameOfReferenceUID') + ): + raise ValueError( + 'Cannot determine spatial relationship because datasets ' + 'lack a frame of reference UID.' + ) + if dataset_from.FrameOfReferenceUID != dataset_to.FrameOfReferenceUID: + raise ValueError( + 'Datasets do not share a frame of reference, so the spatial ' + 'relationship between them is not defined.' + ) + + pos_f, ori_f, spa_f, _ = _get_spatial_information( + dataset_from, + frame_number=frame_number_from, + for_total_pixel_matrix=for_total_pixel_matrix_from, + ) + pos_t, ori_t, spa_t, _ = _get_spatial_information( + dataset_to, + frame_number=frame_number_to, + for_total_pixel_matrix=for_total_pixel_matrix_to, + ) + return cls( + image_position_from=pos_f, + image_orientation_from=ori_f, + pixel_spacing_from=spa_f, + image_position_to=pos_t, + image_orientation_to=ori_t, + pixel_spacing_to=spa_t, + round_output=round_output, + ) class ImageToReferenceTransformer: @@ -610,7 +1637,7 @@ def __call__(self, coordinates: np.ndarray) -> np.ndarray: 'Argument "coordinates" must be a two-dimensional array ' 'with shape [n, 2].' ) - image_coordinates = np.row_stack([ + image_coordinates = np.vstack([ coordinates.T.astype(float), np.zeros((coordinates.shape[0], ), dtype=float), np.ones((coordinates.shape[0], ), dtype=float), @@ -618,6 +1645,47 @@ def __call__(self, coordinates: np.ndarray) -> np.ndarray: reference_coordinates = np.dot(self._affine, image_coordinates) return reference_coordinates[:3, :].T + @classmethod + def for_image( + cls, + dataset: Dataset, + frame_number: Optional[int] = None, + for_total_pixel_matrix: bool = False, + ) -> 'ImageToReferenceTransformer': + """Construct a transformer for a given image or image frame. + + Parameters + ---------- + dataset: pydicom.Dataset + Dataset representing an image. + frame_number: Union[int, None], optional + Frame number (using 1-based indexing) of the frame for which to get + the transformer. This should be provided if and only if the dataset + is a multi-frame image. + for_total_pixel_matrix: bool, optional + If True, use the spatial information for the total pixel matrix of + a tiled image. The result will be a transformer that maps image + coordinates of the total pixel matrix to frame of reference + coordinates. This should only be True if the image is a tiled image + and is incompatible with specifying a frame number. + + Returns + ------- + highdicom.spatial.ImageToReferenceTransformer: + Transformer object for the given image, or image frame. + + """ + position, orientation, spacing, _ = _get_spatial_information( + dataset, + frame_number=frame_number, + for_total_pixel_matrix=for_total_pixel_matrix, + ) + return cls( + image_position=position, + image_orientation=orientation, + pixel_spacing=spacing, + ) + class ReferenceToImageTransformer: @@ -668,7 +1736,8 @@ def __init__( image_position: Sequence[float], image_orientation: Sequence[float], pixel_spacing: Sequence[float], - spacing_between_slices: float = 1.0 + spacing_between_slices: float = 1.0, + drop_slice_coord: bool = False, ): """Construct transformation object. @@ -698,6 +1767,12 @@ def __init__( spacing_between_slices: float, optional Distance (in the coordinate defined by the frame of reference) between neighboring slices. Default: 1 + drop_slice_coord: bool, optional + Whether to remove the 3rd column of the output array + (representing the out-of-plane coordinate) and return a 2D output + array. If this option is taken, and the resulting coordinates + do not lie in the range -0.5 to 0.5, a ``RuntimeError`` will be + triggered. Raises ------ @@ -709,6 +1784,7 @@ def __init__( an incorrect length. """ + self._drop_slice_coord = drop_slice_coord # Image coordinates are shifted relative to pixel matrix indices by # 0.5 pixels and we thus have to correct for this shift. correction_affine = np.array([ @@ -748,13 +1824,13 @@ def __call__(self, coordinates: np.ndarray) -> np.ndarray: ------- numpy.ndarray Array of (column, row, slice) indices, where `column` and `row` are - zero-based indices to the total pixel matrix and the `slice` index - represents the signed distance of the input coordinate in the - direction normal to the plane of the total pixel matrix. - The `row` and `column` indices are constrained by the dimension of - the total pixel matrix. Note, however, that in general, the - resulting coordinate may not lie within the imaging plane, and - consequently the `slice` offset may be non-zero. + zero-based coordinates in the total pixel matrix and the `slice` + index represents the signed distance of the input coordinate in the + direction normal to the plane of the total pixel matrix. The `row` + and `column` coordinates are constrained by the dimension of the + total pixel matrix. Note, however, that in general, the resulting + coordinate may not lie within the imaging plane, and consequently + the `slice` offset may be non-zero. Raises ------ @@ -767,12 +1843,322 @@ def __call__(self, coordinates: np.ndarray) -> np.ndarray: 'Argument "coordinates" must be a two-dimensional array ' 'with shape [n, 3].' ) - reference_coordinates = np.row_stack([ + reference_coordinates = np.vstack([ coordinates.T.astype(float), np.ones((coordinates.shape[0], ), dtype=float) ]) image_coordinates = np.dot(self._affine, reference_coordinates) - return image_coordinates[:3, :].T + image_coordinates = image_coordinates[:3, :].T + if self._drop_slice_coord: + if np.abs(image_coordinates[:, 2]).max() > 0.5: + raise RuntimeError( + "Output coordinates do not lie within the given image " + "plane." + ) + image_coordinates = image_coordinates[:, :2] + return image_coordinates + + @classmethod + def for_image( + cls, + dataset: Dataset, + frame_number: Optional[int] = None, + for_total_pixel_matrix: bool = False, + drop_slice_coord: bool = False, + ) -> 'ReferenceToImageTransformer': + """Construct a transformer for a given image or image frame. + + Parameters + ---------- + dataset: pydicom.Dataset + Dataset representing an image. + frame_number: Union[int, None], optional + Frame number (using 1-based indexing) of the frame for which to get + the transformer. This should be provided if and only if the dataset + is a multi-frame image. + for_total_pixel_matrix: bool, optional + If True, use the spatial information for the total pixel matrix of + a tiled image. The result will be a transformer that maps frame of + reference coordinates to indices of the total pixel matrix. This + should only be True if the image is a tiled image and is + incompatible with specifying a frame number. + drop_slice_coord: bool, optional + Whether to remove the 3rd column of the output array + (representing the out-of-plane coordinate) and return a 2D output + array. If this option is taken, and the resulting coordinates + do not lie in the range -0.5 to 0.5, a ``RuntimeError`` will be + triggered. + + Returns + ------- + highdicom.spatial.ReferenceToImageTransformer: + Transformer object for the given image, or image frame. + + """ + ( + position, + orientation, + spacing, + slice_spacing, + ) = _get_spatial_information( + dataset, + frame_number=frame_number, + for_total_pixel_matrix=for_total_pixel_matrix, + ) + if slice_spacing is None: + slice_spacing = 1.0 + return cls( + image_position=position, + image_orientation=orientation, + pixel_spacing=spacing, + spacing_between_slices=slice_spacing, + drop_slice_coord=drop_slice_coord, + ) + + +class ImageToImageTransformer: + + """Class for transforming image coordinates between two images. + + This class facilitates the mapping of image coordinates of + an image or an image frame (tile or plane) into those of another image or + image frame in the same frame of reference. This can include (but is not + limited) to mapping between different frames of the same image, or + different images within the same series (e.g. two levels of a spatial + pyramid). However, it is required that the two images be coplanar + within the frame-of-reference coordinate system. + + Image coordinates are (column, row) pairs of floating-point values, where + the (0.0, 0.0) point is located at the top left corner of the top left hand + corner pixel of the pixel matrix. Image coordinates have pixel units at + sub-pixel resolution. + + Examples + -------- + + Create a transformer for two images, where the second image has an axis + flipped relative to the first. + + >>> transformer = ImageToImageTransformer( + ... image_position_from=[0.0, 0.0, 0.0], + ... image_orientation_from=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + ... pixel_spacing_from=[1.0, 1.0], + ... image_position_to=[0.0, 100.0, 0.0], + ... image_orientation_to=[1.0, 0.0, 0.0, 0.0, -1.0, 0.0], + ... pixel_spacing_to=[1.0, 1.0], + ... ) + + >>> coords_in = np.array([[0, 0], [50, 50]]) + >>> coords_out = transformer(coords_in) + >>> print(coords_out) + [[ 0. 101.] + [ 50. 51.]] + + Warning + ------- + This class shall not be used to pixel indices between images. Use the + :class:`highdicom.spatial.PixelToPixelTransformer` class instead. + + """ + + def __init__( + self, + image_position_from: Sequence[float], + image_orientation_from: Sequence[float], + pixel_spacing_from: Sequence[float], + image_position_to: Sequence[float], + image_orientation_to: Sequence[float], + pixel_spacing_to: Sequence[float], + ): + """Construct transformation object. + + The resulting object will map image coordinates of the "from" image to + image coordinates of the "to" image. + + Parameters + ---------- + image_position_from: Sequence[float] + Position of the "from" image in the frame of reference, + i.e., the offset of the top left hand corner pixel in the pixel + matrix from the origin of the reference coordinate system along the + X, Y, and Z axis + image_orientation_from: Sequence[float] + Cosines of the row direction (first triplet: horizontal, left to + right, increasing column index) and the column direction (second + triplet: vertical, top to bottom, increasing row index) direction + of the "from" image expressed in the three-dimensional patient or + slide coordinate system defined by the frame of reference + pixel_spacing_from: Sequence[float] + Spacing between pixels of the "from" imagem in millimeter unit + along the column direction (first value: spacing between rows, + vertical, top to bottom, increasing row index) and the rows + direction (second value: spacing between columns: horizontal, left + to right, increasing column index) + image_position_to: Sequence[float] + Position of the "to" image using the same definition as the "from" + image. + image_orientation_to: Sequence[float] + Orientation cosines of the "to" image using the same definition as + the "from" image. + pixel_spacing_to: Sequence[float] + Pixel spacing of the "to" image using the same definition as + the "from" image. + + Raises + ------ + TypeError + When any of the arguments is not a sequence. + ValueError + When any of the arguments has an incorrect length, or if the two + images are not coplanar in the frame of reference coordinate + system. + + """ + if not _are_images_coplanar( + image_position_a=image_position_from, + image_orientation_a=image_orientation_from, + image_position_b=image_position_to, + image_orientation_b=image_orientation_to, + ): + raise ValueError( + "To two images do not exist within the same plane " + "in the frame of reference. and therefore pixel-to-pixel " + "transformation is not possible." + ) + # Image coordinates are shifted relative to pixel matrix indices by + # 0.5 pixels and we thus have to correct for this shift. + pix_to_im = np.array([ + [1.0, 0.0, 0.0, 0.5], + [0.0, 1.0, 0.0, 0.5], + [0.0, 0.0, 1.0, 0.0], + [0.0, 0.0, 0.0, 1.0], + ]) + ref_to_pix = _create_inv_affine_transformation_matrix( + image_position=image_position_to, + image_orientation=image_orientation_to, + pixel_spacing=pixel_spacing_to, + ) + pix_to_ref = _create_affine_transformation_matrix( + image_position=image_position_from, + image_orientation=image_orientation_from, + pixel_spacing=pixel_spacing_from, + ) + im_to_pix = np.array([ + [1.0, 0.0, 0.0, -0.5], + [0.0, 1.0, 0.0, -0.5], + [0.0, 0.0, 1.0, 0.0], + [0.0, 0.0, 0.0, 1.0], + ]) + self._affine = pix_to_im @ ref_to_pix @ pix_to_ref @ im_to_pix + + @property + def affine(self) -> np.ndarray: + """numpy.ndarray: 4x4 affine transformation matrix""" + return self._affine + + def __call__(self, coordinates: np.ndarray) -> np.ndarray: + """Transform pixel indices between two images. + + Parameters + ---------- + indices: numpy.ndarray + Array of (column, row) coordinates at sub-pixel resolution in the + range [0, Columns] and [0, Rows], respectively. + Array of floating-point values with shape ``(n, 2)``, where *n* is + the number of coordinates, the first column represents the `column` + values and the second column represents the `row` values. + The ``(0.0, 0.0)`` coordinate is located at the top left corner + of the top left hand corner pixel in the total pixel matrix. + + Returns + ------- + numpy.ndarray + Array of (column, row) image coordinates in the "to" image. + + Raises + ------ + ValueError + When `coordinates` has incorrect shape. + + """ + if coordinates.shape[1] != 2: + raise ValueError( + 'Argument "coordinates" must be a two-dimensional array ' + 'with shape [n, 2].' + ) + image_coordinates = np.vstack([ + coordinates.T.astype(float), + np.zeros((coordinates.shape[0], ), dtype=float), + np.ones((coordinates.shape[0], ), dtype=float), + ]) + output_coordinates = np.dot(self._affine, image_coordinates) + return output_coordinates[:2, :].T + + @classmethod + def for_images( + cls, + dataset_from: Dataset, + dataset_to: Dataset, + frame_number_from: Optional[int] = None, + frame_number_to: Optional[int] = None, + for_total_pixel_matrix_from: bool = False, + for_total_pixel_matrix_to: bool = False, + ) -> 'ImageToImageTransformer': + """Construct a transformer for two given images or image frames. + + Parameters + ---------- + dataset: pydicom.Dataset + Dataset representing an image. + frame_number: Union[int, None], optional + Frame number (using 1-based indexing) of the frame for which to get + the transformer. This should be provided if and only if the dataset + is a multi-frame image. + for_total_pixel_matrix: bool, optional + If True, use the spatial information for the total pixel matrix of + a tiled image. The result will be a transformer that maps image + coordinates of the total pixel matrix to frame of reference + coordinates. This should only be True if the image is a tiled image + and is incompatible with specifying a frame number. + + Returns + ------- + highdicom.spatial.ImageToImageTransformer: + Transformer object for the given image, or image frame. + + """ + if ( + not hasattr(dataset_from, 'FrameOfReferenceUID') or + not hasattr(dataset_to, 'FrameOfReferenceUID') + ): + raise ValueError( + 'Cannot determine spatial relationship because datasets ' + 'lack a frame of reference UID.' + ) + if dataset_from.FrameOfReferenceUID != dataset_to.FrameOfReferenceUID: + raise ValueError( + 'Datasets do not share a frame of reference, so the spatial ' + 'relationship between them is not defined.' + ) + + pos_f, ori_f, spa_f, _ = _get_spatial_information( + dataset_from, + frame_number=frame_number_from, + for_total_pixel_matrix=for_total_pixel_matrix_from, + ) + pos_t, ori_t, spa_t, _ = _get_spatial_information( + dataset_to, + frame_number=frame_number_to, + for_total_pixel_matrix=for_total_pixel_matrix_to, + ) + return cls( + image_position_from=pos_f, + image_orientation_from=ori_f, + pixel_spacing_from=spa_f, + image_position_to=pos_t, + image_orientation_to=ori_t, + pixel_spacing_to=spa_t, + ) def map_pixel_into_coordinate_system( diff --git a/src/highdicom/sr/sop.py b/src/highdicom/sr/sop.py index 44ed97a2..f60eaeeb 100644 --- a/src/highdicom/sr/sop.py +++ b/src/highdicom/sr/sop.py @@ -1,11 +1,13 @@ """Module for SOP Classes of Structured Report (SR) IODs.""" import datetime +from itertools import chain import logging from collections import defaultdict from copy import deepcopy from os import PathLike from typing import ( Any, + Generator, cast, Mapping, List, @@ -351,6 +353,116 @@ def content(self) -> ContentSequence: """highdicom.sr.value_types.ContentSequence: SR document content""" return self._content + def get_evidence( + self, + current_procedure_only: bool = False, + ) -> List[Tuple[UID, UID, UID, UID]]: + """Get a list of all SOP Instances listed as evidence in this SR. + + Parameters + ---------- + current_procedure_only: bool, optional + If True, return only evidence created in order to satisfy the + current requested procedure (found in the + *CurrentRequestedProcedureEvidenceSequence*). If False, also + include other evidence (found in the + *PertinentOtherEvidenceSequence*). + + Returns + ------- + List[Tuple[highdicom.UID, highdicom.UID, highdicom.UID, highdicom.UID]]: + List of tuples of UIDs, each representing a single instance. Each + tuple consists of (StudyInstanceUID, SeriesInstanceUID, + SOPInstanceUID, SOPClassUID). + + """ + def extract_evidence( + sequence: DataElementSequence, + ) -> Generator[Tuple[UID, UID, UID, UID], None, None]: + for item in sequence: + for series_ds in item.ReferencedSeriesSequence: + for instance_ds in series_ds.ReferencedSOPSequence: + yield ( + UID(item.StudyInstanceUID), + UID(series_ds.SeriesInstanceUID), + UID(instance_ds.ReferencedSOPInstanceUID), + UID(instance_ds.ReferencedSOPClassUID), + ) + + current_evidence_seq = self.get( + 'CurrentRequestedProcedureEvidenceSequence' + ) + if current_evidence_seq is not None: + current_evidence = extract_evidence(current_evidence_seq) + else: + current_evidence = [] + + other_evidence_seq = self.get('PertinentOtherEvidenceSequence') + if other_evidence_seq is not None and not current_procedure_only: + other_evidence = extract_evidence(other_evidence_seq) + else: + other_evidence = [] + + evidence = list(chain(current_evidence, other_evidence)) + + # Deduplicate the list while preserving order + evidence = list(dict.fromkeys(evidence)) + + return evidence + + def get_evidence_series( + self, + current_procedure_only: bool = False, + ) -> List[Tuple[UID, UID]]: + """Get a list of all series listed as evidence in this SR. + + Parameters + ---------- + current_procedure_only: bool, optional + If True, return only evidence created in order to satisfy the + current requested procedure (found in the + *CurrentRequestedProcedureEvidenceSequence*). If False, also + include other evidence (found in the + *PertinentOtherEvidenceSequence*). + + Returns + ------- + List[Tuple[highdicom.UID, highdicom.UID]]: + List of tuples of UIDs, each representing a single series. Each + tuple consists of (StudyInstanceUID, SeriesInstanceUID). + + """ + def extract_evidence_series( + sequence: DataElementSequence, + ) -> Generator[Tuple[UID, UID], None, None]: + for item in sequence: + for series_ds in item.ReferencedSeriesSequence: + yield ( + UID(item.StudyInstanceUID), + UID(series_ds.SeriesInstanceUID), + ) + + current_evidence_seq = self.get( + 'CurrentRequestedProcedureEvidenceSequence' + ) + if current_evidence_seq is not None: + current_evidence = extract_evidence_series(current_evidence_seq) + else: + current_evidence = [] + + other_evidence_seq = self.get('PertinentOtherEvidenceSequence') + if other_evidence_seq is not None and not current_procedure_only: + other_evidence = extract_evidence_series(other_evidence_seq) + else: + other_evidence = [] + + evidence = list(chain(current_evidence, other_evidence)) + + # Deduplicate the list while preserving order + evidence = list(dict.fromkeys(evidence)) + + return evidence + class EnhancedSR(_SR): diff --git a/src/highdicom/sr/templates.py b/src/highdicom/sr/templates.py index f102d357..0eb54607 100644 --- a/src/highdicom/sr/templates.py +++ b/src/highdicom/sr/templates.py @@ -562,7 +562,7 @@ def _get_coded_modality(sop_class_uid: str) -> Code: '1.2.840.10008.5.1.4.1.1.14.2': codes.cid29.IntravascularOpticalCoherenceTomography, # noqa: E501 '1.2.840.10008.5.1.4.1.1.20': codes.cid29.NuclearMedicine, '1.2.840.10008.5.1.4.1.1.66.1': codes.cid32.Registration, - '1.2.840.10008.5.1.4.1.1.66.2': codes.cid32.SpatialFiducials, + '1.2.840.10008.5.1.4.1.1.66.2': codes.cid32.SpatialFiducialsProducer, '1.2.840.10008.5.1.4.1.1.66.3': codes.cid32.Registration, '1.2.840.10008.5.1.4.1.1.66.4': codes.cid32.Segmentation, '1.2.840.10008.5.1.4.1.1.67': codes.cid32.RealWorldValueMap, @@ -1703,7 +1703,7 @@ def __init__( Identifier of the observed specimen (may have limited scope, e.g., only relevant with respect to the corresponding container) container_identifier: Union[str, None], optional - Identifier of the container holding the speciment (e.g., a glass + Identifier of the container holding the specimen (e.g., a glass slide) specimen_type: Union[pydicom.sr.coding.Code, highdicom.sr.CodedConcept, None], optional Type of the specimen (see @@ -1809,6 +1809,49 @@ def specimen_type(self) -> Union[CodedConcept, None]: return matches[0].value return None + @classmethod + def from_image( + cls, + image: Dataset, + ) -> 'SubjectContextSpecimen': + """Deduce specimen information from an existing image. + + This is appropriate, for example, when copying the specimen information + from a source image into a derived SR or similar object. + + Parameters + ---------- + image: pydicom.Dataset + An image from which to infer specimen information. There is no + limitation on the type of image, however it must have the Specimen + module included. + + Raises + ------ + ValueError: + If the input image does not contain specimen information. + + """ + if not hasattr(image, 'ContainerIdentifier'): + raise ValueError("Image does not contain specimen information.") + + description = image.SpecimenDescriptionSequence[0] + + # Specimen type code sequence is optional + if hasattr(description, 'SpecimenTypeCodeSequence'): + specimen_type: Optional[CodedConcept] = CodedConcept.from_dataset( + description.SpecimenTypeCodeSequence[0] + ) + else: + specimen_type = None + + return cls( + container_identifier=image.ContainerIdentifier, + identifier=description.SpecimenIdentifier, + uid=description.SpecimenUID, + specimen_type=specimen_type, + ) + @classmethod def from_sequence( cls, @@ -2124,6 +2167,41 @@ def __init__( raise TypeError('Unexpected subject class specific context.') self.extend(subject_class_specific_context) + @classmethod + def from_image(cls, image: Dataset) -> 'Optional[SubjectContext]': + """Get a subject context inferred from an existing image. + + Currently this is only supported for subjects that are specimens. + + Parameters + ---------- + image: pydicom.Dataset + Dataset of an existing DICOM image object + containing metadata on the imaging subject. Highdicom will attempt + to infer the subject context from this image. If successful, it + will be returned as a ``SubjectContext``, otherwise ``None``. + + Returns + ------- + Optional[highdicom.sr.SubjectContext]: + SubjectContext, if it can be inferred from the image. Otherwise, + ``None``. + + """ + try: + subject_context_specimen = SubjectContextSpecimen.from_image( + image + ) + except ValueError: + pass + else: + return cls( + subject_class=codes.DCM.Specimen, + subject_class_specific_context=subject_context_specimen, + ) + + return None + @property def subject_class(self) -> CodedConcept: """highdicom.sr.CodedConcept: type of subject""" diff --git a/src/highdicom/utils.py b/src/highdicom/utils.py index 7ef56207..99230701 100644 --- a/src/highdicom/utils.py +++ b/src/highdicom/utils.py @@ -1,128 +1,32 @@ import itertools -from typing import Iterator, Generator, List, Optional, Sequence, Tuple +from typing import List, Optional, Sequence +import warnings -import numpy as np from pydicom.dataset import Dataset from highdicom.content import PlanePositionSequence from highdicom.enum import CoordinateSystemNames from highdicom.spatial import ( map_pixel_into_coordinate_system, - PixelToReferenceTransformer, + tile_pixel_matrix, + get_tile_array, + iter_tiled_full_frame_data, + is_tiled_image, ) -def tile_pixel_matrix( - total_pixel_matrix_rows: int, - total_pixel_matrix_columns: int, - rows: int, - columns: int, -) -> Iterator[Tuple[int, int]]: - """Tiles an image into smaller frames (rectangular regions). - - Follows the convention used in image with Dimension Organization Type - "TILED_FULL" images. - - Parameters - ---------- - total_pixel_matrix_rows: int - Number of rows in the Total Pixel Matrix - total_pixel_matrix_columns: int - Number of columns in the Total Pixel Matrix - rows: int - Number of rows per Frame (tile) - columns: int - Number of columns per Frame (tile) - - Returns - ------- - Iterator - One-based (Column, Row) index of each Frame (tile) - - """ - tiles_per_col = int(np.ceil(total_pixel_matrix_rows / rows)) - tiles_per_row = int(np.ceil(total_pixel_matrix_columns / columns)) - tile_row_indices = iter(range(1, tiles_per_col + 1)) - tile_col_indices = iter(range(1, tiles_per_row + 1)) - return ( - (c, r) for (r, c) in itertools.product( - tile_row_indices, - tile_col_indices - ) - ) - - -def get_tile_array( - pixel_array: np.ndarray, - row_offset: int, - column_offset: int, - tile_rows: int, - tile_columns: int, - pad: bool = True, -) -> np.ndarray: - """Extract a tile from a total pixel matrix array. - - Parameters - ---------- - pixel_array: np.ndarray - Array representing a total pixel matrix. The first two dimensions - are treated as the rows and columns, respectively, of the total pixel - matrix. Any subsequent dimensions are not used but are retained in the - output array. - row_offset: int - Offset of the first row of the requested tile from the top of the total - pixel matrix (1-based index). - column_offset: int - Offset of the first column of the requested tile from the left of the - total pixel matrix (1-based index). - tile_rows: int - Number of rows per tile. - tile_columns: - Number of columns per tile. - pad: bool - Whether to pad the returned array with zeros at the right and/or bottom - to ensure that it matches the correct tile size. Otherwise, the returned - array is not padded and may be smaller than the full tile size. - - Returns - ------- - np.ndarray: - Returned pixel array for the requested tile. - - """ - if row_offset < 1 or row_offset > pixel_array.shape[0]: - raise ValueError( - "Row offset must be between 1 and the size of dimension 0 of the " - "pixel array." - ) - if column_offset < 1 or column_offset > pixel_array.shape[1]: - raise ValueError( - "Column offset must be between 1 and the size of dimension 1 of " - "the pixel array." - ) - # Move to pythonic 1-based indexing - row_offset -= 1 - column_offset -= 1 - row_end = row_offset + tile_rows - if row_end > pixel_array.shape[0]: - pad_rows = row_end - pixel_array.shape[0] - row_end = pixel_array.shape[0] - else: - pad_rows = 0 - column_end = column_offset + tile_columns - if column_end > pixel_array.shape[1]: - pad_columns = column_end - pixel_array.shape[1] - column_end = pixel_array.shape[1] - else: - pad_columns = 0 - # Account for 1-based to 0-based index conversion - tile_array = pixel_array[row_offset:row_end, column_offset:column_end] - if pad_rows > 0 or pad_columns > 0: - extra_dims = pixel_array.ndim - 2 - padding = [(0, pad_rows), (0, pad_columns)] + [(0, 0)] * extra_dims - tile_array = np.pad(tile_array, padding) - - return tile_array +# Several functions that were initially defined in this module were moved to +# highdicom.spatial to consolidate similar functionality and prevent circular +# dependencies. Therefore they are re-exported here for backwards compatibility +__all__ = [ + "tile_pixel_matrix", # backwards compatibility + "get_tile_array", # backwards compatibility + "iter_tiled_full_frame_data", # backwards compatibility + "is_tiled_image", # backwards compatibility + "compute_plane_position_slide_per_frame", + "compute_plane_position_tiled_full", + "are_plane_positions_tiled_full", +] def compute_plane_position_tiled_full( @@ -134,7 +38,7 @@ def compute_plane_position_tiled_full( columns: int, image_orientation: Sequence[float], pixel_spacing: Sequence[float], - slice_thickness: Optional[float] = None, + slice_thickness: Optional[float] = None, # unused (deprecated) spacing_between_slices: Optional[float] = None, slice_index: Optional[int] = None ) -> PlanePositionSequence: @@ -177,7 +81,10 @@ def compute_plane_position_tiled_full( increasing Row index) and the row direction (second value: spacing between columns, horizontal, left to right, increasing Column index) slice_thickness: Union[float, None], optional - Thickness of a focal plane in micrometers + This parameter is unused and passing anything other than None will + cause a warning to be issued. Use spacing_between_slices to specify the + spacing between neighboring slices. This parameter will be removed in a + future version of the library. spacing_between_slices: Union[float, None], optional Distance between neighboring focal planes in micrometers slice_index: Union[int, None], optional @@ -195,6 +102,12 @@ def compute_plane_position_tiled_full( When only one of `slice_index` and `spacing_between_slices` is provided """ + if slice_thickness is not None: + warnings.warn( + "Passing a slice_thickness other than None has no effect and " + "will be deprecated in a future version of the library.", + UserWarning + ) if row_index < 1 or column_index < 1: raise ValueError("Row and column indices must be positive integers.") row_offset_frame = ((row_index - 1) * rows) @@ -232,158 +145,6 @@ def compute_plane_position_tiled_full( ) -def iter_tiled_full_frame_data( - dataset: Dataset, -) -> Generator[Tuple[int, int, int, int, float, float, float], None, None]: - """Get data on the position of each tile in a TILED_FULL image. - - This works only with images with Dimension Organization Type of - "TILED_FULL". - - Unlike :func:`highdicom.utils.compute_plane_position_slide_per_frame`, - this functions returns the data in their basic Python types rather than - wrapping as :class:`highdicom.PlanePositionSequence` - - Parameters - ---------- - dataset: pydicom.dataset.Dataset - VL Whole Slide Microscopy Image or Segmentation Image using the - "TILED_FULL" DimensionOrganizationType. - - Returns - ------- - channel: int - 1-based integer index of the "channel". The meaning of "channel" - depends on the image type. For segmentation images, the channel is the - segment number. For other images, it is the optical path number. - focal_plane_index: int - 1-based integer index of the focal plane. - column_position: int - 1-based column position of the tile (measured left from the left side - of the total pixel matrix). - row_position: int - 1-based row position of the tile (measured down from the top of the - total pixel matrix). - x: float - X coordinate in the frame-of-reference coordinate system in millimeter - units. - y: float - Y coordinate in the frame-of-reference coordinate system in millimeter - units. - z: float - Z coordinate in the frame-of-reference coordinate system in millimeter - units. - - """ - allowed_sop_class_uids = { - '1.2.840.10008.5.1.4.1.1.77.1.6', # VL Whole Slide Microscopy Image - '1.2.840.10008.5.1.4.1.1.66.4', # Segmentation Image - } - if dataset.SOPClassUID not in allowed_sop_class_uids: - raise ValueError( - 'Expected a VL Whole Slide Microscopy Image or Segmentation Image.' - ) - if ( - not hasattr(dataset, "DimensionOrganizationType") or - dataset.DimensionOrganizationType != "TILED_FULL" - ): - raise ValueError( - 'Expected an image with "TILED_FULL" dimension organization type.' - ) - - image_origin = dataset.TotalPixelMatrixOriginSequence[0] - image_orientation = ( - float(dataset.ImageOrientationSlide[0]), - float(dataset.ImageOrientationSlide[1]), - float(dataset.ImageOrientationSlide[2]), - float(dataset.ImageOrientationSlide[3]), - float(dataset.ImageOrientationSlide[4]), - float(dataset.ImageOrientationSlide[5]), - ) - tiles_per_column = int( - np.ceil(dataset.TotalPixelMatrixRows / dataset.Rows) - ) - tiles_per_row = int( - np.ceil(dataset.TotalPixelMatrixColumns / dataset.Columns) - ) - num_focal_planes = getattr( - dataset, - 'TotalPixelMatrixFocalPlanes', - 1 - ) - - is_segmentation = dataset.SOPClassUID == '1.2.840.10008.5.1.4.1.1.66.4' - - # The "channels" output is either segment for segmentations, or optical - # path for other images - if is_segmentation: - num_channels = len(dataset.SegmentSequence) - else: - num_channels = getattr( - dataset, - 'NumberOfOpticalPaths', - len(dataset.OpticalPathSequence) - ) - - shared_fg = dataset.SharedFunctionalGroupsSequence[0] - pixel_measures = shared_fg.PixelMeasuresSequence[0] - pixel_spacing = ( - float(pixel_measures.PixelSpacing[0]), - float(pixel_measures.PixelSpacing[1]), - ) - spacing_between_slices = float( - getattr( - pixel_measures, - 'SpacingBetweenSlices', - 1.0 - ) - ) - x_offset = image_origin.XOffsetInSlideCoordinateSystem - y_offset = image_origin.YOffsetInSlideCoordinateSystem - - # Array of tile indices (col_index, row_index) - tile_indices = np.array( - [ - (c, r) for (r, c) in - itertools.product( - range(1, tiles_per_column + 1), - range(1, tiles_per_row + 1) - ) - ] - ) - - # Pixel offsets of each in the total pixel matrix - frame_pixel_offsets = ( - (tile_indices - 1) * np.array([dataset.Columns, dataset.Rows]) - ) - - for channel in range(1, num_channels + 1): - for slice_index in range(1, num_focal_planes + 1): - # These checks are needed for mypy to determine the correct type - z_offset = float(slice_index - 1) * spacing_between_slices - transformer = PixelToReferenceTransformer( - image_position=(x_offset, y_offset, z_offset), - image_orientation=image_orientation, - pixel_spacing=pixel_spacing - ) - - reference_coordinates = transformer(frame_pixel_offsets) - - for offsets, coords in zip( - frame_pixel_offsets, - reference_coordinates - ): - yield ( - channel, - slice_index, - int(offsets[0] + 1), - int(offsets[1] + 1), - float(coords[0]), - float(coords[1]), - float(coords[2]), - ) - - def compute_plane_position_slide_per_frame( dataset: Dataset ) -> List[PlanePositionSequence]: @@ -420,24 +181,6 @@ def compute_plane_position_slide_per_frame( ] -def is_tiled_image(dataset: Dataset) -> bool: - """Determine whether a dataset represents a tiled image. - - Returns - ------- - bool: - True if the dataset is a tiled image. False otherwise. - - """ - if ( - hasattr(dataset, 'TotalPixelMatrixRows') and - hasattr(dataset, 'TotalPixelMatrixColumns') and - hasattr(dataset, 'NumberOfFrames') - ): - return True - return False - - def are_plane_positions_tiled_full( plane_positions: Sequence[PlanePositionSequence], rows: int, diff --git a/src/highdicom/version.py b/src/highdicom/version.py index 81edede8..08a9dbff 100644 --- a/src/highdicom/version.py +++ b/src/highdicom/version.py @@ -1 +1 @@ -__version__ = '0.22.0' +__version__ = '0.23.0' diff --git a/tests/test_base.py b/tests/test_base.py index 01198392..ce547eb3 100644 --- a/tests/test_base.py +++ b/tests/test_base.py @@ -74,7 +74,7 @@ def test_big_endian(self): ) def test_explicit_vr(self): - sop_class = SOPClass( + _ = SOPClass( study_instance_uid=UID(), series_instance_uid=UID(), series_number=1, @@ -85,11 +85,9 @@ def test_explicit_vr(self): manufacturer='highdicom', transfer_syntax_uid=ExplicitVRLittleEndian, ) - assert not sop_class.is_implicit_VR - assert sop_class.is_little_endian def test_implicit_vr(self): - sop_class = SOPClass( + _ = SOPClass( study_instance_uid=UID(), series_instance_uid=UID(), series_number=1, @@ -100,8 +98,6 @@ def test_implicit_vr(self): manufacturer='highdicom', transfer_syntax_uid=ImplicitVRLittleEndian, ) - assert sop_class.is_implicit_VR - assert sop_class.is_little_endian class TestEndianCheck(unittest.TestCase): diff --git a/tests/test_frame.py b/tests/test_frame.py index 35477244..d99e5131 100644 --- a/tests/test_frame.py +++ b/tests/test_frame.py @@ -4,8 +4,10 @@ import numpy as np import pytest from pydicom.uid import ( + JPEG2000, JPEG2000Lossless, JPEGLSLossless, + JPEGLSNearLossless, JPEGBaseline8Bit, ) @@ -23,6 +25,7 @@ def setUp(self): ) def test_jpeg_rgb(self): + pytest.importorskip("libjpeg") filepath = str(self._test_files_dir.joinpath('frame_rgb.jpeg')) with open(filepath, 'br') as fp: compressed_frame = fp.read() @@ -68,6 +71,7 @@ def test_jpeg_rgb(self): ) def test_jpeg_rgb_empty(self): + pytest.importorskip("libjpeg") filepath = str(self._test_files_dir.joinpath('frame_rgb_empty.jpeg')) with open(filepath, 'br') as fp: compressed_frame = fp.read() @@ -159,15 +163,75 @@ def test_jpeg_monochrome(self): assert compressed_frame.endswith(b'\xFF\xD9') def test_jpeg2000_rgb(self): + pytest.importorskip("openjpeg") bits_allocated = 8 - frame = np.ones((16, 32, 3), dtype=np.dtype(f'uint{bits_allocated}')) + frame = np.ones((48, 32, 3), dtype=np.dtype(f'uint{bits_allocated}')) + frame[2:4, 5:30, 0] = 7 + compressed_frame = encode_frame( + frame, + transfer_syntax_uid=JPEG2000, + bits_allocated=bits_allocated, + bits_stored=bits_allocated, + photometric_interpretation='YBR_ICT', + pixel_representation=0, + planar_configuration=0 + ) + assert compressed_frame.startswith(b"\xFF\x4F\xFF\x51") + assert compressed_frame.endswith(b'\xFF\xD9') + decoded_frame = decode_frame( + value=compressed_frame, + transfer_syntax_uid=JPEG2000, + rows=frame.shape[0], + columns=frame.shape[1], + samples_per_pixel=frame.shape[2], + bits_allocated=bits_allocated, + bits_stored=bits_allocated, + photometric_interpretation='YBR_ICT', + pixel_representation=0, + planar_configuration=0 + ) + np.testing.assert_allclose(frame, decoded_frame, atol=2) + + def test_jpeg2000_monochrome(self): + pytest.importorskip("openjpeg") + bits_allocated = 8 + frame = np.zeros((48, 32), dtype=np.dtype(f'uint{bits_allocated}')) + frame[2:4, 5:30] = 7 + compressed_frame = encode_frame( + frame, + transfer_syntax_uid=JPEG2000, + bits_allocated=bits_allocated, + bits_stored=bits_allocated, + photometric_interpretation='MONOCHROME2', + pixel_representation=0, + ) + assert compressed_frame.startswith(b"\xFF\x4F\xFF\x51") + assert compressed_frame.endswith(b'\xFF\xD9') + decoded_frame = decode_frame( + value=compressed_frame, + transfer_syntax_uid=JPEG2000, + rows=frame.shape[0], + columns=frame.shape[1], + samples_per_pixel=1, + bits_allocated=bits_allocated, + bits_stored=bits_allocated, + photometric_interpretation='MONOCHROME2', + pixel_representation=0, + planar_configuration=0 + ) + np.testing.assert_allclose(frame, decoded_frame, atol=2) + + def test_jpeg2000lossless_rgb(self): + pytest.importorskip("openjpeg") + bits_allocated = 8 + frame = np.ones((48, 32, 3), dtype=np.dtype(f'uint{bits_allocated}')) frame *= 255 compressed_frame = encode_frame( frame, transfer_syntax_uid=JPEG2000Lossless, bits_allocated=bits_allocated, bits_stored=bits_allocated, - photometric_interpretation='YBR_FULL', + photometric_interpretation='YBR_RCT', pixel_representation=0, planar_configuration=0 ) @@ -181,15 +245,45 @@ def test_jpeg2000_rgb(self): samples_per_pixel=frame.shape[2], bits_allocated=bits_allocated, bits_stored=bits_allocated, - photometric_interpretation='YBR_FULL', + photometric_interpretation='YBR_RCT', pixel_representation=0, planar_configuration=0 ) np.testing.assert_array_equal(frame, decoded_frame) - def test_jpeg2000_monochrome(self): + def test_jpeg2000lossless_monochrome(self): + pytest.importorskip("openjpeg") bits_allocated = 16 - frame = np.zeros((16, 32), dtype=np.dtype(f'uint{bits_allocated}')) + frame = np.zeros((48, 32), dtype=np.dtype(f'uint{bits_allocated}')) + compressed_frame = encode_frame( + frame, + transfer_syntax_uid=JPEG2000Lossless, + bits_allocated=bits_allocated, + bits_stored=bits_allocated, + photometric_interpretation='MONOCHROME2', + pixel_representation=0, + ) + assert compressed_frame.startswith(b"\xFF\x4F\xFF\x51") + assert compressed_frame.endswith(b'\xFF\xD9') + decoded_frame = decode_frame( + value=compressed_frame, + transfer_syntax_uid=JPEG2000Lossless, + rows=frame.shape[0], + columns=frame.shape[1], + samples_per_pixel=1, + bits_allocated=bits_allocated, + bits_stored=bits_allocated, + photometric_interpretation='MONOCHROME2', + pixel_representation=0, + planar_configuration=0 + ) + np.testing.assert_array_equal(frame, decoded_frame) + + def test_jpeg2000lossless_single_bit(self): + pytest.importorskip("openjpeg") + bits_allocated = 1 + frame = np.zeros((48, 32), dtype=np.dtype('uint8')) + frame[12:45, 3:6] = 1 compressed_frame = encode_frame( frame, transfer_syntax_uid=JPEG2000Lossless, @@ -224,7 +318,7 @@ def test_jpegls_rgb(self): transfer_syntax_uid=JPEGLSLossless, bits_allocated=bits_allocated, bits_stored=bits_allocated, - photometric_interpretation='YBR_FULL', + photometric_interpretation='RGB', pixel_representation=0, planar_configuration=0 ) @@ -238,7 +332,7 @@ def test_jpegls_rgb(self): samples_per_pixel=frame.shape[2], bits_allocated=bits_allocated, bits_stored=bits_allocated, - photometric_interpretation='YBR_FULL', + photometric_interpretation='RGB', pixel_representation=0, planar_configuration=0 ) @@ -272,6 +366,64 @@ def test_jpegls_monochrome(self): ) np.testing.assert_array_equal(frame, decoded_frame) + def test_jpeglsnearlossless_rgb(self): + pytest.importorskip("libjpeg") + bits_allocated = 8 + frame = np.ones((16, 32, 3), dtype=np.dtype(f'uint{bits_allocated}')) + frame *= 255 + compressed_frame = encode_frame( + frame, + transfer_syntax_uid=JPEGLSNearLossless, + bits_allocated=bits_allocated, + bits_stored=bits_allocated, + photometric_interpretation='RGB', + pixel_representation=0, + planar_configuration=0 + ) + assert compressed_frame.startswith(b'\xFF\xD8') + assert compressed_frame.endswith(b'\xFF\xD9') + decoded_frame = decode_frame( + value=compressed_frame, + transfer_syntax_uid=JPEGLSNearLossless, + rows=frame.shape[0], + columns=frame.shape[1], + samples_per_pixel=frame.shape[2], + bits_allocated=bits_allocated, + bits_stored=bits_allocated, + photometric_interpretation='RGB', + pixel_representation=0, + planar_configuration=0 + ) + np.testing.assert_allclose(frame, decoded_frame) + + def test_jpeglsnearlossless_monochrome(self): + pytest.importorskip("libjpeg") + bits_allocated = 16 + frame = np.zeros((16, 32), dtype=np.dtype(f'uint{bits_allocated}')) + compressed_frame = encode_frame( + frame, + transfer_syntax_uid=JPEGLSNearLossless, + bits_allocated=bits_allocated, + bits_stored=bits_allocated, + photometric_interpretation='MONOCHROME2', + pixel_representation=0, + ) + assert compressed_frame.startswith(b'\xFF\xD8') + assert compressed_frame.endswith(b'\xFF\xD9') + decoded_frame = decode_frame( + value=compressed_frame, + transfer_syntax_uid=JPEGLSNearLossless, + rows=frame.shape[0], + columns=frame.shape[1], + samples_per_pixel=1, + bits_allocated=bits_allocated, + bits_stored=bits_allocated, + photometric_interpretation='MONOCHROME2', + pixel_representation=0, + planar_configuration=0 + ) + np.testing.assert_allclose(frame, decoded_frame) + def test_jpeg_rgb_wrong_photometric_interpretation(self): frame = np.ones((16, 32, 3), dtype=np.uint8) with pytest.raises(ValueError): diff --git a/tests/test_pm.py b/tests/test_pm.py index ae18d5ce..3538f300 100644 --- a/tests/test_pm.py +++ b/tests/test_pm.py @@ -6,7 +6,7 @@ import numpy as np import pytest from pydicom import dcmread -from pydicom.data import get_testdata_files +from pydicom.data import get_testdata_file, get_testdata_files from pydicom.sr.codedict import codes from pydicom.sr.coding import Code from pydicom.uid import ( @@ -232,6 +232,10 @@ def setUp(self): str(data_dir.joinpath('test_files', 'ct_image.dcm')) ) + self._ct_multiframe_image = dcmread( + get_testdata_file('eCT_Supplemental.dcm') + ) + self._sm_image = dcmread( str(data_dir.joinpath('test_files', 'sm_image.dcm')) ) @@ -337,6 +341,11 @@ def test_multi_frame_sm_image_single_native(self): window_width=window_width, content_label=content_label ) + + # Work around pydicom 3 decoding issue (should be able to remove this + # soon) + pmap.pixel_array_options(use_v2_backend=True) + assert pmap.SOPClassUID == '1.2.840.10008.5.1.4.1.1.30' assert pmap.SOPInstanceUID == self._sop_instance_uid assert pmap.SeriesInstanceUID == self._series_instance_uid @@ -546,10 +555,11 @@ def test_multi_frame_sm_image_ushort_palette_lut(self): assert instance.PixelPresentation == 'MONOCHROME' def test_multi_frame_sm_image_ushort_encapsulated_jpeg2000(self): + pytest.importorskip("openjpeg") pixel_array = np.random.randint( low=0, high=2**8, - size=self._sm_image.pixel_array.shape[:3], + size=self._ct_multiframe_image.pixel_array.shape[:3], dtype=np.uint8 ) window_center = 128 @@ -564,7 +574,7 @@ def test_multi_frame_sm_image_ushort_encapsulated_jpeg2000(self): slope=1 ) pmap = ParametricMap( - [self._sm_image], + [self._ct_multiframe_image], pixel_array, self._series_instance_uid, self._series_number, @@ -651,6 +661,11 @@ def test_single_frame_ct_image_double(self): window_width=window_width, ) assert pmap.BitsAllocated == 64 + + # Work around pydicom 3 decoding issue (should be able to remove this + # soon) + pmap.pixel_array_options(use_v2_backend=True) + assert np.array_equal(pmap.pixel_array, pixel_array) def test_single_frame_ct_image_ushort_native(self): @@ -766,6 +781,10 @@ def test_series_single_frame_ct_image_single(self): window_width=window_width, ) + # Work around pydicom 3 decoding issue (should be able to remove this + # soon) + pmap.pixel_array_options(use_v2_backend=True) + assert np.array_equal(pmap.pixel_array, pixel_array) def test_multi_frame_sm_image_with_spatial_positions_not_preserved(self): diff --git a/tests/test_sc.py b/tests/test_sc.py index b2c63b60..6a6e7c82 100644 --- a/tests/test_sc.py +++ b/tests/test_sc.py @@ -4,13 +4,15 @@ import numpy as np import pytest -from pydicom.encaps import generate_pixel_data_frame +from pydicom.encaps import generate_fragmented_frames from pydicom.filereader import dcmread from pydicom.uid import ( RLELossless, JPEGBaseline8Bit, + JPEG2000, JPEG2000Lossless, JPEGLSLossless, + JPEGLSNearLossless, ) from pydicom.valuerep import DA, TM @@ -65,7 +67,7 @@ def get_array_after_writing(instance): pixel_representation=ds.PixelRepresentation, planar_configuration=getattr(ds, 'PlanarConfiguration', None) ) - for value in generate_pixel_data_frame(instance.PixelData) + for (value, ) in generate_fragmented_frames(instance.PixelData) ] if len(decoded_frame_arrays) > 1: return np.stack(decoded_frame_arrays) @@ -318,6 +320,7 @@ def test_rgb_rle(self): ) def test_monochrome_jpeg_baseline(self): + pytest.importorskip("libjpeg") bits_allocated = 8 photometric_interpretation = 'MONOCHROME2' coordinate_system = 'PATIENT' @@ -347,6 +350,7 @@ def test_monochrome_jpeg_baseline(self): ) def test_rgb_jpeg_baseline(self): + pytest.importorskip("libjpeg") bits_allocated = 8 photometric_interpretation = 'YBR_FULL_422' coordinate_system = 'PATIENT' @@ -373,7 +377,8 @@ def test_rgb_jpeg_baseline(self): reread_frame = self.get_array_after_writing(instance) np.testing.assert_allclose(frame, reread_frame, rtol=1.2) - def test_monochrome_jpeg2000(self): + def test_monochrome_jpeg2000lossless(self): + pytest.importorskip("openjpeg") bits_allocated = 8 photometric_interpretation = 'MONOCHROME2' coordinate_system = 'PATIENT' @@ -400,9 +405,39 @@ def test_monochrome_jpeg2000(self): frame ) + def test_monochrome_jpeg2000(self): + pytest.importorskip("openjpeg") + bits_allocated = 8 + photometric_interpretation = 'MONOCHROME2' + coordinate_system = 'PATIENT' + frame = np.random.randint(0, 256, size=(256, 256), dtype=np.uint8) + instance = SCImage( + pixel_array=frame, + photometric_interpretation=photometric_interpretation, + bits_allocated=bits_allocated, + coordinate_system=coordinate_system, + study_instance_uid=self._study_instance_uid, + series_instance_uid=self._series_instance_uid, + sop_instance_uid=self._sop_instance_uid, + series_number=self._series_number, + instance_number=self._instance_number, + manufacturer=self._manufacturer, + patient_orientation=self._patient_orientation, + transfer_syntax_uid=JPEG2000 + ) + + assert instance.file_meta.TransferSyntaxUID == JPEG2000 + + assert np.allclose( + self.get_array_after_writing(instance), + frame, + atol=2 + ) + def test_rgb_jpeg2000(self): + pytest.importorskip("openjpeg") bits_allocated = 8 - photometric_interpretation = 'YBR_FULL' + photometric_interpretation = 'YBR_RCT' coordinate_system = 'PATIENT' frame = np.random.randint(0, 256, size=(256, 256, 3), dtype=np.uint8) instance = SCImage( @@ -455,10 +490,39 @@ def test_monochrome_jpegls(self): frame ) + def test_monochrome_jpegls_near_lossless(self): + pytest.importorskip("libjpeg") + bits_allocated = 16 + photometric_interpretation = 'MONOCHROME2' + coordinate_system = 'PATIENT' + frame = np.random.randint(0, 2**16, size=(256, 256), dtype=np.uint16) + instance = SCImage( + pixel_array=frame, + photometric_interpretation=photometric_interpretation, + bits_allocated=bits_allocated, + coordinate_system=coordinate_system, + study_instance_uid=self._study_instance_uid, + series_instance_uid=self._series_instance_uid, + sop_instance_uid=self._sop_instance_uid, + series_number=self._series_number, + instance_number=self._instance_number, + manufacturer=self._manufacturer, + patient_orientation=self._patient_orientation, + transfer_syntax_uid=JPEGLSNearLossless + ) + + assert instance.file_meta.TransferSyntaxUID == JPEGLSNearLossless + + assert np.allclose( + self.get_array_after_writing(instance), + frame, + atol=2 + ) + def test_rgb_jpegls(self): pytest.importorskip("libjpeg") bits_allocated = 8 - photometric_interpretation = 'YBR_FULL' + photometric_interpretation = 'RGB' coordinate_system = 'PATIENT' frame = np.random.randint(0, 256, size=(256, 256, 3), dtype=np.uint8) instance = SCImage( diff --git a/tests/test_seg.py b/tests/test_seg.py index cbcd5d8c..a25a40ae 100644 --- a/tests/test_seg.py +++ b/tests/test_seg.py @@ -657,6 +657,22 @@ def setUp(self): self._sm_image = dcmread( str(data_dir.joinpath('test_files', 'sm_image.dcm')) ) + + # Hack around the openjpeg bug encoding images smaller than 32 pixels + frame = self._sm_image.pixel_array + frame = np.pad(frame, ((0, 0), (12, 12), (12, 12), (0, 0))) + self._sm_image.PixelData = frame.flatten().tobytes() + self._sm_image.TotalPixelMatrixRows = ( + frame.shape[1] * + int(self._sm_image.TotalPixelMatrixRows / self._sm_image.Rows) + ) + self._sm_image.TotalPixelMatrixColumns = ( + frame.shape[2] * + int(self._sm_image.TotalPixelMatrixColumns / self._sm_image.Columns) + ) + self._sm_image.Rows = frame.shape[1] + self._sm_image.Columns = frame.shape[2] + # Override te existing ImageOrientationSlide to make the frame ordering # simpler for the tests self._sm_pixel_array = np.zeros( @@ -697,6 +713,16 @@ def setUp(self): ct_series, key=lambda x: x.ImagePositionPatient[2] ) + + # Hack around the fact that the images are too small to be encoded by + # openjpeg + for im in self._ct_series: + frame = im.pixel_array + frame = np.pad(frame, ((8, 8), (8, 8))) + im.Rows = frame.shape[0] + im.Columns = frame.shape[1] + im.PixelData = frame.flatten().tobytes() + self._ct_series_mask_array = np.zeros( (len(self._ct_series), ) + self._ct_series[0].pixel_array.shape, dtype=bool @@ -725,7 +751,13 @@ def setUp(self): # Fixtures to use to parametrize segmentation creation # Using this fixture mechanism, we can parametrize class methods @staticmethod - @pytest.fixture(params=[ExplicitVRLittleEndian, ImplicitVRLittleEndian]) + @pytest.fixture( + params=[ + ExplicitVRLittleEndian, + ImplicitVRLittleEndian, + JPEG2000Lossless, + ] + ) def binary_transfer_syntax_uid(request): return request.param @@ -1497,6 +1529,36 @@ def test_construction_tiled_full(self): assert instance.DimensionOrganizationType == "TILED_FULL" assert not hasattr(instance, "PerFrameFunctionalGroupsSequence") + def test_construction_further_source_images(self): + further_source_image = deepcopy(self._ct_image) + series_uid = UID() + sop_uid = UID() + further_source_image.SeriesInstanceUID = series_uid + further_source_image.SOPInstanceUID = sop_uid + instance = Segmentation( + [self._ct_image], + self._ct_pixel_array, + SegmentationTypeValues.FRACTIONAL.value, + self._segment_descriptions, + self._series_instance_uid, + self._series_number, + self._sop_instance_uid, + self._instance_number, + self._manufacturer, + self._manufacturer_model_name, + self._software_versions, + self._device_serial_number, + content_label=self._content_label, + further_source_images=[further_source_image], + ) + assert len(instance.SourceImageSequence) == 2 + further_item = instance.SourceImageSequence[1] + assert further_item.ReferencedSOPInstanceUID == sop_uid + + assert len(instance.ReferencedSeriesSequence) == 2 + further_item = instance.ReferencedSeriesSequence[1] + assert further_item.SeriesInstanceUID == series_uid + @staticmethod @pytest.fixture( params=[ @@ -1519,10 +1581,8 @@ def segmentation_type(request): @pytest.fixture( params=[ None, - (10, 10), - (10, 25), - (25, 25), - (30, 30), + (32, 32), + (48, 48), ]) def tile_size(request): return request.param @@ -1579,7 +1639,17 @@ def test_construction_autotile( else: omit_empty_frames_values = [False, True] - transfer_syntax_uids = [ExplicitVRLittleEndian] + transfer_syntax_uids = [ + ExplicitVRLittleEndian, + ] + try: + import openjpeg # noqa: F401 + except ModuleNotFoundError: + pass + else: + transfer_syntax_uids += [ + JPEG2000Lossless, + ] if segmentation_type.value == 'FRACTIONAL': try: import libjpeg # noqa: F401 @@ -1587,7 +1657,6 @@ def test_construction_autotile( pass else: transfer_syntax_uids += [ - JPEG2000Lossless, JPEGLSLossless, ] @@ -1687,6 +1756,8 @@ def test_pixel_types_fractional( pix_type, test_data, ): + if fractional_transfer_syntax_uid == JPEG2000Lossless: + pytest.importorskip("openjpeg") if fractional_transfer_syntax_uid == JPEGLSLossless: pytest.importorskip("libjpeg") @@ -1845,6 +1916,8 @@ def test_pixel_types_binary( pix_type, test_data, ): + if binary_transfer_syntax_uid == JPEG2000Lossless: + pytest.importorskip("openjpeg") sources, mask = self._tests[test_data] # Two segments, overlapping @@ -3761,9 +3834,31 @@ def setUp(self): ) self._seg_pix[5:15, 3:8] = 1 + self._seg_pix_multisegment = np.zeros( + (1, *tpm_size, 3), + dtype=np.uint8, + ) + self._seg_pix_multisegment[0, 5:15, 3:8, 0] = 1 + self._seg_pix_multisegment[0, 23:42, 13:18, 1] = 1 + self._seg_pix_multisegment[0, 45:48, 25:30, 1] = 1 + + self._seg_pix_fractional = np.zeros( + tpm_size, + dtype=np.float32, + ) + self._seg_pix_fractional[5:15, 3:8] = 0.5 + self._seg_pix_fractional[3:5, 10:12] = 1.0 + self._n_downsamples = 3 self._downsampled_pix_arrays = [self._seg_pix] + self._downsampled_pix_arrays_multisegment = [self._seg_pix_multisegment] + self._downsampled_pix_arrays_fractional = [self._seg_pix_fractional] seg_pil = Image.fromarray(self._seg_pix) + seg_pil_fractional = Image.fromarray(self._seg_pix_fractional) + seg_pil_multisegment = [ + Image.fromarray(self._seg_pix_multisegment[0, :, :, i]) + for i in range(self._seg_pix_multisegment.shape[3]) + ] pyramid_uid = UID() self._source_pyramid = [deepcopy(self._sm_image)] self._source_pyramid[0].PyramidUID = pyramid_uid @@ -3780,6 +3875,22 @@ def setUp(self): ) self._downsampled_pix_arrays.append(resized) + resized_fractional = np.array( + seg_pil_fractional.resize(out_size, Image.Resampling.BILINEAR) + ) + self._downsampled_pix_arrays_fractional.append(resized_fractional) + + resized_multisegment = np.stack( + [ + im.resize(out_size, Image.Resampling.NEAREST) + for im in seg_pil_multisegment + ], + axis=-1 + )[None] + self._downsampled_pix_arrays_multisegment.append( + resized_multisegment + ) + # Mock lower-resolution source images. No need to have their pixel # data correctly set as it isn't used. Just update the relevant # metadata @@ -3820,6 +3931,21 @@ def setUp(self): ) ), ] + self._segment_descriptions_multi = [ + SegmentDescription( + segment_number=i, + segment_label=f'Segment #{i}', + segmented_property_category=self._segmented_property_category, + segmented_property_type=self._segmented_property_type, + algorithm_type=SegmentAlgorithmTypeValues.AUTOMATIC.value, + algorithm_identification=AlgorithmIdentificationSequence( + name='bla', + family=codes.DCM.ArtificialIntelligence, + version='v1' + ) + ) + for i in range(1, 4) + ] def test_pyramid_factors(self): downsample_factors = [2.0, 5.0] @@ -3925,6 +4051,83 @@ def test_multiple_source_single_pixel_array(self): pix ) + def test_multiple_source_single_pixel_array_multisegment(self): + # Test construction when given multiple source images and a single + # segmentation image + segs = create_segmentation_pyramid( + source_images=self._source_pyramid, + pixel_arrays=[self._seg_pix_multisegment], + segmentation_type=SegmentationTypeValues.BINARY, + segment_descriptions=self._segment_descriptions_multi, + series_instance_uid=UID(), + series_number=1, + manufacturer='Foo', + manufacturer_model_name='Bar', + software_versions='1', + device_serial_number='123', + ) + + assert len(segs) == len(self._source_pyramid) + for pix, seg in zip(self._downsampled_pix_arrays_multisegment, segs): + assert hasattr(seg, 'PyramidUID') + seg_pix = seg.get_total_pixel_matrix() + print("pixel array", pix.shape, seg_pix.shape) + print("pixel array", pix.max(), seg_pix.max()) + assert np.array_equal( + seg.get_total_pixel_matrix(), + pix[0] + ) + + def test_multiple_source_single_pixel_array_float(self): + # Test construction when given multiple source images and a single + # segmentation image + segs = create_segmentation_pyramid( + source_images=self._source_pyramid, + pixel_arrays=[self._seg_pix.astype(np.float32)], + segmentation_type=SegmentationTypeValues.BINARY, + segment_descriptions=self._segment_descriptions, + series_instance_uid=UID(), + series_number=1, + manufacturer='Foo', + manufacturer_model_name='Bar', + software_versions='1', + device_serial_number='123', + ) + + assert len(segs) == len(self._source_pyramid) + for pix, seg in zip(self._downsampled_pix_arrays, segs): + assert hasattr(seg, 'PyramidUID') + assert np.array_equal( + seg.get_total_pixel_matrix(combine_segments=True), + pix + ) + + def test_multiple_source_single_pixel_array_fractional(self): + # Test construction when given multiple source images and a single + # segmentation image + segs = create_segmentation_pyramid( + source_images=self._source_pyramid, + pixel_arrays=[self._seg_pix_fractional], + segmentation_type=SegmentationTypeValues.FRACTIONAL, + segment_descriptions=self._segment_descriptions, + series_instance_uid=UID(), + series_number=1, + manufacturer='Foo', + manufacturer_model_name='Bar', + software_versions='1', + device_serial_number='123', + max_fractional_value=200, + ) + + assert len(segs) == len(self._source_pyramid) + for pix, seg in zip(self._downsampled_pix_arrays_fractional, segs): + assert hasattr(seg, 'PyramidUID') + assert np.allclose( + seg.get_total_pixel_matrix(combine_segments=False)[:, :, 0], + pix, + atol=0.005, # 1/200 + ) + def test_multiple_source_multiple_pixel_arrays(self): # Test construction when given multiple source images and multiple # segmentation images @@ -3948,3 +4151,53 @@ def test_multiple_source_multiple_pixel_arrays(self): seg.get_total_pixel_matrix(combine_segments=True), pix ) + + def test_multiple_source_multiple_pixel_arrays_multisegment(self): + # Test construction when given multiple source images and multiple + # segmentation images + segs = create_segmentation_pyramid( + source_images=self._source_pyramid, + pixel_arrays=self._downsampled_pix_arrays_multisegment, + segmentation_type=SegmentationTypeValues.BINARY, + segment_descriptions=self._segment_descriptions_multi, + series_instance_uid=UID(), + series_number=1, + manufacturer='Foo', + manufacturer_model_name='Bar', + software_versions='1', + device_serial_number='123', + ) + + assert len(segs) == len(self._source_pyramid) + for pix, seg in zip(self._downsampled_pix_arrays_multisegment, segs): + assert hasattr(seg, 'PyramidUID') + assert np.array_equal( + seg.get_total_pixel_matrix(), + pix[0] + ) + + def test_multiple_source_multiple_pixel_arrays_multisegment_labelmap(self): + # Test construction when given multiple source images and multiple + # segmentation images + mask = np.argmax(self._seg_pix_multisegment, axis=3).astype(np.uint8) + segs = create_segmentation_pyramid( + source_images=self._source_pyramid, + pixel_arrays=mask, + segmentation_type=SegmentationTypeValues.BINARY, + segment_descriptions=self._segment_descriptions_multi, + series_instance_uid=UID(), + series_number=1, + manufacturer='Foo', + manufacturer_model_name='Bar', + software_versions='1', + device_serial_number='123', + ) + + assert len(segs) == len(self._source_pyramid) + for pix, seg in zip(self._downsampled_pix_arrays_multisegment, segs): + mask = np.argmax(pix, axis=3).astype(np.uint8) + assert hasattr(seg, 'PyramidUID') + assert np.array_equal( + seg.get_total_pixel_matrix(combine_segments=True), + mask[0] + ) diff --git a/tests/test_spatial.py b/tests/test_spatial.py index bce8d1f5..1d53cdc6 100644 --- a/tests/test_spatial.py +++ b/tests/test_spatial.py @@ -1,14 +1,51 @@ +from pathlib import Path import numpy as np +import pydicom +from pydicom.data import get_testdata_file import pytest from highdicom.spatial import ( + ImageToImageTransformer, ImageToReferenceTransformer, + PixelToPixelTransformer, PixelToReferenceTransformer, ReferenceToImageTransformer, ReferenceToPixelTransformer, + is_tiled_image, + _are_images_coplanar, ) +@pytest.mark.parametrize( + 'filepath,expected_output', + [ + ( + Path(__file__).parents[1].joinpath('data/test_files/ct_image.dcm'), + False + ), + ( + Path(__file__).parents[1].joinpath('data/test_files/sm_image.dcm'), + True + ), + ( + Path(__file__).parents[1].joinpath( + 'data/test_files/seg_image_ct_binary.dcm' + ), + False + ), + ( + Path(__file__).parents[1].joinpath( + 'data/test_files/seg_image_sm_control.dcm' + ), + True + ), + ] +) +def test_is_tiled_image(filepath, expected_output): + dcm = pydicom.dcmread(filepath) + assert is_tiled_image(dcm) == expected_output + + params_pixel_to_reference = [ # Slide pytest.param( @@ -251,7 +288,7 @@ (89.0, 47.0, -10.0), ]), np.array([ - (158, -6, -10), + (158, -6.5, -10), ]) ), # Patient @@ -413,6 +450,143 @@ ] +params_pixel_to_pixel = [ + # Pure translation + pytest.param( + dict( + image_position_from=(56.0, 34.2, 1.0), + image_orientation_from=(1.0, 0.0, 0.0, 0.0, 1.0, 0.0), + pixel_spacing_from=(1.0, 1.0), + image_position_to=(66.0, 32.2, 1.0), + image_orientation_to=(1.0, 0.0, 0.0, 0.0, 1.0, 0.0), + pixel_spacing_to=(1.0, 1.0), + ), + np.array([ + (0, 0), + (50, 50), + ]), + np.array([ + (-10, 2), + (40, 52), + ]) + ), + # Two images with different spacings (e.g. different pyramid levels) + pytest.param( + dict( + image_position_from=(56.0, 34.2, 1.0), + image_orientation_from=(1.0, 0.0, 0.0, 0.0, 1.0, 0.0), + pixel_spacing_from=(1.0, 1.0), + image_position_to=(56.0, 34.2, 1.0), + image_orientation_to=(1.0, 0.0, 0.0, 0.0, 1.0, 0.0), + pixel_spacing_to=(0.5, 0.5), + ), + np.array([ + (0, 0), + (50, 50), + ]), + np.array([ + (0, 0), + (100, 100), + ]) + ), + # 30 degree rotation, anisotropic scale + pytest.param( + dict( + image_position_from=(56.0, 34.2, 1.0), + image_orientation_from=(1.0, 0.0, 0.0, 0.0, 1.0, 0.0), + pixel_spacing_from=(1.0, 1.0), + image_position_to=(56.0, 34.2, 1.0), + image_orientation_to=( + np.cos(np.radians(30)), -np.sin(np.radians(30)), 0.0, + np.sin(np.radians(30)), np.cos(np.radians(30)), 0.0, + ), + pixel_spacing_to=(0.25, 0.5), + ), + np.array([ + (0, 0), + (50, 50), + ]), + np.array([ + (0, 0), + ( + 2.0 * ( + 50 * np.cos(np.radians(30)) - + 50 * np.sin(np.radians(30)) + ), + 4.0 * ( + 50 * np.sin(np.radians(30)) + + 50 * np.cos(np.radians(30)) + ), + ), + ]) + ), +] + + +params_image_to_image = [ + # Pure translation + pytest.param( + dict( + image_position_from=(56.0, 34.2, 1.0), + image_orientation_from=(1.0, 0.0, 0.0, 0.0, 1.0, 0.0), + pixel_spacing_from=(1.0, 1.0), + image_position_to=(66.0, 32.2, 1.0), + image_orientation_to=(1.0, 0.0, 0.0, 0.0, 1.0, 0.0), + pixel_spacing_to=(1.0, 1.0), + ), + np.array([ + (0.0, 0.0), + (50.0, 50.0), + ]), + np.array([ + (-10.0, 2.0), + (40.0, 52.0), + ]) + ), + # Two images with different spacings (e.g. different pyramid levels) + pytest.param( + dict( + image_position_from=(56.0, 34.2, 1.0), + image_orientation_from=(1.0, 0.0, 0.0, 0.0, 1.0, 0.0), + pixel_spacing_from=(1.0, 1.0), + image_position_to=(56.0, 34.2, 1.0), + image_orientation_to=(1.0, 0.0, 0.0, 0.0, 1.0, 0.0), + pixel_spacing_to=(0.5, 0.5), + ), + np.array([ + (0.5, 0.5), + (25, 50), + ]), + np.array([ # y = 2x - 0.5 + (0.5, 0.5), + (49.5, 99.5), + ]) + ), + # 30 degree rotation, anisotropic scale + pytest.param( + dict( + image_position_from=(56.0, 34.2, 1.0), + image_orientation_from=(1.0, 0.0, 0.0, 0.0, 1.0, 0.0), + pixel_spacing_from=(1.0, 1.0), + image_position_to=(56.0, 34.2, 1.0), + image_orientation_to=( + np.cos(np.radians(30)), -np.sin(np.radians(30)), 0.0, + np.sin(np.radians(30)), np.cos(np.radians(30)), 0.0, + ), + pixel_spacing_to=(0.25, 0.5), + ), + np.array([ + (0, 0), + (50, 50), + ]), + np.array([ + (0.133974596, -2.23205081), + (36.7365150, 270.973030), + ]) + ), +] + + @pytest.mark.parametrize( 'params,inputs,expected_outputs', params_pixel_to_reference @@ -423,13 +597,46 @@ def test_map_pixel_into_coordinate_system(params, inputs, expected_outputs): np.testing.assert_array_almost_equal(outputs, expected_outputs) +@pytest.mark.parametrize( + 'round_output', + [False, True], +) +@pytest.mark.parametrize( + 'drop_slice_index', + [False, True], +) @pytest.mark.parametrize( 'params,inputs,expected_outputs', params_reference_to_pixel ) -def test_map_coordinate_into_pixel_matrix(params, inputs, expected_outputs): - transform = ReferenceToPixelTransformer(**params) +def test_map_coordinate_into_pixel_matrix( + params, + inputs, + expected_outputs, + round_output, + drop_slice_index, +): + transform = ReferenceToPixelTransformer( + round_output=round_output, + drop_slice_index=drop_slice_index, + **params, + ) + if round_output: + expected_outputs = np.around(expected_outputs) + expected_dtype = np.int64 + else: + expected_dtype = np.float64 + if drop_slice_index: + if np.abs(expected_outputs[:, 2]).max() >= 0.5: + # In these cases, the transform should fail + # because the dropped slice index is no close + # to zero + with pytest.raises(RuntimeError): + transform(inputs) + return + expected_outputs = expected_outputs[:, :2] outputs = transform(inputs) + assert outputs.dtype == expected_dtype np.testing.assert_array_almost_equal(outputs, expected_outputs) @@ -443,11 +650,218 @@ def test_map_image_to_reference_coordinate(params, inputs, expected_outputs): np.testing.assert_array_almost_equal(outputs, expected_outputs) +@pytest.mark.parametrize( + 'drop_slice_coord', + [False, True], +) @pytest.mark.parametrize( 'params,inputs,expected_outputs', params_reference_to_image ) -def test_map_reference_to_image_coordinate(params, inputs, expected_outputs): - transform = ReferenceToImageTransformer(**params) +def test_map_reference_to_image_coordinate( + params, + inputs, + expected_outputs, + drop_slice_coord, +): + transform = ReferenceToImageTransformer( + drop_slice_coord=drop_slice_coord, + **params, + ) + if drop_slice_coord: + if np.abs(expected_outputs[:, 2]).max() >= 0.5: + # In these cases, the transform should fail + # because the dropped slice index is no close + # to zero + with pytest.raises(RuntimeError): + transform(inputs) + return + expected_outputs = expected_outputs[:, :2] + outputs = transform(inputs) + np.testing.assert_array_almost_equal(outputs, expected_outputs) + + +@pytest.mark.parametrize( + 'round_output', + [False, True], +) +@pytest.mark.parametrize( + 'params,inputs,expected_outputs', + params_pixel_to_pixel +) +def test_map_indices_between_images( + params, + inputs, + expected_outputs, + round_output, +): + transform = PixelToPixelTransformer( + round_output=round_output, + **params, + ) outputs = transform(inputs) + if round_output: + expected_outputs = np.around(expected_outputs) + assert outputs.dtype == np.int64 + else: + assert outputs.dtype == np.float64 + np.testing.assert_array_almost_equal(outputs, expected_outputs) + + +@pytest.mark.parametrize( + 'params,inputs,expected_outputs', + params_image_to_image +) +def test_map_coordinates_between_images(params, inputs, expected_outputs): + transform = ImageToImageTransformer(**params) + outputs = transform(inputs) + np.testing.assert_array_almost_equal(outputs, expected_outputs) + + +all_single_image_transformer_classes = [ + ImageToReferenceTransformer, + PixelToReferenceTransformer, + ReferenceToPixelTransformer, + ReferenceToImageTransformer, +] + + +all_image_pair_transformer_classes = [ + ImageToImageTransformer, + PixelToPixelTransformer, +] + + +@pytest.mark.parametrize( + 'transformer_cls', + all_single_image_transformer_classes +) +def test_transformers_ct_image(transformer_cls): + file_path = Path(__file__) + data_dir = file_path.parent.parent.joinpath('data/test_files') + dcm = pydicom.dcmread(data_dir / 'ct_image.dcm') + transformer_cls.for_image(dcm) + + +@pytest.mark.parametrize( + 'transformer_cls', + all_single_image_transformer_classes +) +def test_transformers_sm_image(transformer_cls): + file_path = Path(__file__) + data_dir = file_path.parent.parent.joinpath('data/test_files') + dcm = pydicom.dcmread(data_dir / 'sm_image.dcm') + transformer_cls.for_image(dcm, frame_number=3) + + +@pytest.mark.parametrize( + 'transformer_cls', + all_single_image_transformer_classes +) +def test_transformers_sm_image_total_pixel_matrix(transformer_cls): + file_path = Path(__file__) + data_dir = file_path.parent.parent.joinpath('data/test_files') + dcm = pydicom.dcmread(data_dir / 'sm_image.dcm') + transformer_cls.for_image(dcm, for_total_pixel_matrix=True) + + +@pytest.mark.parametrize( + 'transformer_cls', + all_single_image_transformer_classes +) +def test_transformers_seg_sm_image_total_pixel_matrix(transformer_cls): + file_path = Path(__file__) + data_dir = file_path.parent.parent.joinpath('data/test_files') + dcm = pydicom.dcmread(data_dir / 'seg_image_sm_dots.dcm') + transformer_cls.for_image(dcm, for_total_pixel_matrix=True) + + +@pytest.mark.parametrize( + 'transformer_cls', + all_single_image_transformer_classes +) +def test_transformers_enhanced_ct_image(transformer_cls): + dcm = pydicom.dcmread(get_testdata_file('eCT_Supplemental.dcm')) + transformer_cls.for_image(dcm, frame_number=2) + + +@pytest.mark.parametrize( + 'transformer_cls', + all_image_pair_transformer_classes +) +def test_pair_transformers_sm_image_frame_to_frame(transformer_cls): + file_path = Path(__file__) + data_dir = file_path.parent.parent.joinpath('data/test_files') + dcm = pydicom.dcmread(data_dir / 'sm_image.dcm') + transformer_cls.for_images( + dataset_from=dcm, + dataset_to=dcm, + frame_number_from=1, + frame_number_to=3, + ) + + +@pytest.mark.parametrize( + 'transformer_cls', + all_image_pair_transformer_classes +) +def test_pair_transformers_sm_image_tpm_to_frame(transformer_cls): + file_path = Path(__file__) + data_dir = file_path.parent.parent.joinpath('data/test_files') + dcm = pydicom.dcmread(data_dir / 'sm_image.dcm') + transformer_cls.for_images( + dataset_from=dcm, + dataset_to=dcm, + for_total_pixel_matrix_from=True, + frame_number_to=3, + ) + + +@pytest.mark.parametrize( + 'pos_a,ori_a,pos_b,ori_b,result', + [ + pytest.param( + [1.0, 1.0, 1.0], + [1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + [1.0, 1.0, 1.0], + [1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + True, + ), + pytest.param( + [1.0, 1.0, 1.0], + [1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + [1.0, 1.0, 1.0], + [1.0, 0.0, 0.0, 0.0, -1.0, 0.0], # flipped + True, + ), + pytest.param( + [211.0, 11.0, 1.0], # shifted in plane + [1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + [1.0, 1.0, 1.0], + [1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + True, + ), + pytest.param( + [1.0, 1.0, 11.0], # shifted out of plane + [1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + [1.0, 1.0, 1.0], + [1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + False, + ), + pytest.param( + [1.0, 1.0, 1.0], + [1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + [1.0, 1.0, 1.0], + [1.0, 0.0, 0.0, 0.0, 0.0, 1.0], # different orientation + False, + ), + ] +) +def test_are_images_coplanar(pos_a, ori_a, pos_b, ori_b, result): + assert _are_images_coplanar( + image_position_a=pos_a, + image_orientation_a=ori_a, + image_position_b=pos_b, + image_orientation_b=ori_b, + ) == result diff --git a/tests/test_sr.py b/tests/test_sr.py index 4f05cabb..e774d4a7 100644 --- a/tests/test_sr.py +++ b/tests/test_sr.py @@ -1398,6 +1398,114 @@ def test_construction_all(self): assert context[5].TextValue == self._physical_location +class TestSubjectContextSpecimen(unittest.TestCase): + + def setUp(self): + super().setUp() + file_path = Path(__file__) + data_dir = file_path.parent.parent.joinpath('data') + self._sm_image = dcmread( + str(data_dir.joinpath('test_files', 'sm_image.dcm')) + ) + + def test_from_image(self): + specimen_context = SubjectContextSpecimen.from_image( + self._sm_image + ) + + assert ( + specimen_context.container_identifier == + self._sm_image.ContainerIdentifier + ) + description = self._sm_image.SpecimenDescriptionSequence[0] + assert ( + specimen_context.specimen_identifier == + description.SpecimenIdentifier + ) + assert ( + specimen_context.specimen_uid == description.SpecimenUID + ) + assert specimen_context.specimen_type == codes.SCT.TissueSection + + def test_from_image_no_specimen_type(self): + # Specimen type is optional, for_image method should cope correctly if + # it is missing + image = deepcopy(self._sm_image) + delattr( + image.SpecimenDescriptionSequence[0], + 'SpecimenTypeCodeSequence' + ) + specimen_context = SubjectContextSpecimen.from_image(image) + + assert ( + specimen_context.container_identifier == + self._sm_image.ContainerIdentifier + ) + description = self._sm_image.SpecimenDescriptionSequence[0] + assert ( + specimen_context.specimen_identifier == + description.SpecimenIdentifier + ) + assert ( + specimen_context.specimen_uid == description.SpecimenUID + ) + assert specimen_context.specimen_type is None + + +class TestSubjectContext(unittest.TestCase): + + def setUp(self): + super().setUp() + file_path = Path(__file__) + data_dir = file_path.parent.parent.joinpath('data') + self._sm_image = dcmread( + str(data_dir.joinpath('test_files', 'sm_image.dcm')) + ) + self._ct_image = dcmread( + str(data_dir.joinpath('test_files', 'ct_image.dcm')) + ) + + def test_from_image(self): + subject_context = SubjectContext.from_image( + self._sm_image + ) + + assert subject_context is not None + + has_specimen_uid = False + has_specimen_id = False + has_container_id = False + + for item in subject_context: + + # SpecimenUID + specimen_uid = '2.25.281821656492584880365678271074145532563' + if item.ConceptNameCodeSequence[0].CodeValue == '121039': + assert item.UID == specimen_uid + has_specimen_uid = True + + # Specimen Identifier + elif item.ConceptNameCodeSequence[0].CodeValue == '121041': + assert item.TextValue == 'S19-1_A_1_1' + has_specimen_id = True + + # Specimen Container Identifier + elif item.ConceptNameCodeSequence[0].CodeValue == '111700': + assert item.TextValue == 'S19-1_A_1_1' + has_container_id = True + + assert has_specimen_uid + assert has_specimen_id + assert has_container_id + + def test_from_image_no_subject_info(self): + subject_context = SubjectContext.from_image( + self._ct_image + ) + + assert subject_context is None + + class TestObservationContext(unittest.TestCase): def setUp(self): @@ -3772,6 +3880,20 @@ def test_construction(self): performed_procedure_codes=self._performed_procedures ) assert report.SOPClassUID == '1.2.840.10008.5.1.4.1.1.88.22' + evidence = report.get_evidence() + assert len(evidence) == 1 + assert evidence[0] == ( + self._ref_dataset.StudyInstanceUID, + self._ref_dataset.SeriesInstanceUID, + self._ref_dataset.SOPInstanceUID, + self._ref_dataset.SOPClassUID, + ) + evidence_series = report.get_evidence_series() + assert len(evidence_series) == 1 + assert evidence_series[0] == ( + self._ref_dataset.StudyInstanceUID, + self._ref_dataset.SeriesInstanceUID, + ) def test_construction_content_is_sequence(self): report = EnhancedSR( @@ -3913,6 +4035,21 @@ def test_construction(self): with pytest.raises(AttributeError): assert report.PertinentOtherEvidenceSequence + evidence = report.get_evidence() + assert len(evidence) == 1 + assert evidence[0] == ( + self._ref_dataset.StudyInstanceUID, + self._ref_dataset.SeriesInstanceUID, + self._ref_dataset.SOPInstanceUID, + self._ref_dataset.SOPClassUID, + ) + evidence_series = report.get_evidence_series() + assert len(evidence_series) == 1 + assert evidence_series[0] == ( + self._ref_dataset.StudyInstanceUID, + self._ref_dataset.SeriesInstanceUID, + ) + def test_construction_content_is_sequence(self): report = ComprehensiveSR( evidence=[self._ref_dataset], @@ -4029,6 +4166,96 @@ def test_evidence_multiple_studies(self): with pytest.raises(AttributeError): assert report.PertinentOtherEvidenceSequence + evidence = report.get_evidence() + assert len(evidence) == 2 + assert evidence[0] == ( + self._ref_dataset.StudyInstanceUID, + self._ref_dataset.SeriesInstanceUID, + self._ref_dataset.SOPInstanceUID, + self._ref_dataset.SOPClassUID, + ) + assert evidence[1] == ( + ref_dataset.StudyInstanceUID, + ref_dataset.SeriesInstanceUID, + ref_dataset.SOPInstanceUID, + ref_dataset.SOPClassUID, + ) + evidence_series = report.get_evidence_series() + assert len(evidence_series) == 2 + assert evidence_series[0] == ( + self._ref_dataset.StudyInstanceUID, + self._ref_dataset.SeriesInstanceUID, + ) + assert evidence_series[1] == ( + ref_dataset.StudyInstanceUID, + ref_dataset.SeriesInstanceUID, + ) + + def test_current_and_other_evidence(self): + ref_dataset2 = deepcopy(self._ref_dataset) + ref_dataset2.SeriesInstanceUID = '1.2.3' + ref_dataset2.SOPInstanceUID = '1.2.3' + + report = Comprehensive3DSR( + evidence=[self._ref_dataset, ref_dataset2], + content=self._content, + series_instance_uid=self._series_instance_uid, + series_number=self._series_number, + sop_instance_uid=self._sop_instance_uid, + instance_number=self._instance_number, + institution_name=self._institution_name, + institutional_department_name=self._department_name, + manufacturer=self._manufacturer + ) + ref_evd_items = report.CurrentRequestedProcedureEvidenceSequence + assert len(ref_evd_items) == 1 + unref_evd_items = report.PertinentOtherEvidenceSequence + assert len(unref_evd_items) == 1 + + evidence = report.get_evidence() + assert len(evidence) == 2 + assert evidence[0] == ( + self._ref_dataset.StudyInstanceUID, + self._ref_dataset.SeriesInstanceUID, + self._ref_dataset.SOPInstanceUID, + self._ref_dataset.SOPClassUID, + ) + assert evidence[1] == ( + ref_dataset2.StudyInstanceUID, + ref_dataset2.SeriesInstanceUID, + ref_dataset2.SOPInstanceUID, + ref_dataset2.SOPClassUID, + ) + evidence_series = report.get_evidence_series() + assert len(evidence_series) == 2 + assert evidence_series[0] == ( + self._ref_dataset.StudyInstanceUID, + self._ref_dataset.SeriesInstanceUID, + ) + assert evidence_series[1] == ( + ref_dataset2.StudyInstanceUID, + ref_dataset2.SeriesInstanceUID, + ) + + current_evidence = report.get_evidence( + current_procedure_only=True + ) + assert len(current_evidence) == 1 + assert current_evidence[0] == ( + self._ref_dataset.StudyInstanceUID, + self._ref_dataset.SeriesInstanceUID, + self._ref_dataset.SOPInstanceUID, + self._ref_dataset.SOPClassUID, + ) + current_evidence_series = report.get_evidence_series( + current_procedure_only=True + ) + assert len(current_evidence_series) == 1 + assert current_evidence_series[0] == ( + self._ref_dataset.StudyInstanceUID, + self._ref_dataset.SeriesInstanceUID, + ) + def test_srread(self): report = ComprehensiveSR( evidence=[self._ref_dataset], @@ -5178,7 +5405,7 @@ def test_ct_construction(self): value_item = group[4].MeasuredValueSequence[0] unit_code_item = value_item.MeasurementUnitsCodeSequence[0] assert unit_code_item.CodeValue == 'mm' - assert unit_code_item.CodeMeaning == 'millimeter' + assert unit_code_item.CodeMeaning == 'mm' assert unit_code_item.CodingSchemeDesignator == 'UCUM' assert isinstance(group[5], NumContentItem) assert group[5].name == codes.DCM.VerticalPixelSpacing @@ -5189,7 +5416,7 @@ def test_ct_construction(self): value_item = group[6].MeasuredValueSequence[0] unit_code_item = value_item.MeasurementUnitsCodeSequence[0] assert unit_code_item.CodeValue == 'mm' - assert unit_code_item.CodeMeaning == 'millimeter' + assert unit_code_item.CodeMeaning == 'mm' assert unit_code_item.CodingSchemeDesignator == 'UCUM' assert isinstance(group[7], NumContentItem) assert group[7].name == codes.DCM.SliceThickness diff --git a/tests/test_utils.py b/tests/test_utils.py index e036f889..f13d2f3d 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -14,7 +14,6 @@ from highdicom.utils import ( compute_plane_position_tiled_full, compute_plane_position_slide_per_frame, - is_tiled_image, are_plane_positions_tiled_full, ) @@ -144,36 +143,6 @@ def test_compute_plane_position_tiled_full_with_missing_parameters(): ) -@pytest.mark.parametrize( - 'filepath,expected_output', - [ - ( - Path(__file__).parents[1].joinpath('data/test_files/ct_image.dcm'), - False - ), - ( - Path(__file__).parents[1].joinpath('data/test_files/sm_image.dcm'), - True - ), - ( - Path(__file__).parents[1].joinpath( - 'data/test_files/seg_image_ct_binary.dcm' - ), - False - ), - ( - Path(__file__).parents[1].joinpath( - 'data/test_files/seg_image_sm_control.dcm' - ), - True - ), - ] -) -def test_is_tiled_image(filepath, expected_output): - dcm = dcmread(filepath) - assert is_tiled_image(dcm) == expected_output - - def test_compute_plane_position_slide_per_frame(): iterator = itertools.product(range(1, 4), range(1, 3)) for num_optical_paths, num_focal_planes in iterator: diff --git a/tests/utils.py b/tests/utils.py index f1a18c16..f70233a8 100644 --- a/tests/utils.py +++ b/tests/utils.py @@ -7,15 +7,17 @@ def write_and_read_dataset(dataset: Dataset): """Write DICOM dataset to buffer and read it back from buffer.""" clone = Dataset(dataset) - clone.is_little_endian = True if hasattr(dataset, 'file_meta'): clone.file_meta = FileMetaDataset(dataset.file_meta) - if dataset.file_meta.TransferSyntaxUID == '1.2.840.10008.1.2': - clone.is_implicit_VR = True - else: - clone.is_implicit_VR = False + little_endian = None + implicit_vr = None else: - clone.is_implicit_VR = False + little_endian = True + implicit_vr = True with BytesIO() as fp: - clone.save_as(fp) + clone.save_as( + fp, + implicit_vr=implicit_vr, + little_endian=little_endian, + ) return dcmread(fp, force=True)