diff --git a/docs/itk_lib.rst b/docs/itk_lib.rst index 00cbce94..86d65446 100644 --- a/docs/itk_lib.rst +++ b/docs/itk_lib.rst @@ -3,16 +3,23 @@ ITK === -`ITK`_ is a widely-used library for volumetric image processing. In order to -use ITK with highdicom, the ``itk`` python package must be installed separately. -Version 5.4.0 or later is required. +`ITK`_ is a widely-used library for volumetric image processing. While many +simple processing operations can be performed within highdicom itself, ITK +provides a much larger set of tools and many existing image processing +pipelines are built in ITK. By integrating highdicom with ITK, you can +therefore benefit from both ITK's large set of processing features and +highdicom's support for DICOM reading and writing with minimal "glue code". In +particular, correctly transferring spatial metadata between highdicom and ITK +representations is straightforwardly handled for you. In order to use ITK with +highdicom, the ``itk`` python package must be installed separately. Version +5.4.0 or later is required. .. _itk_vol: Volume Conversions ------------------ -Highdicom supports conversions with the ``itk.Image`` class through the +Highdicom supports conversions to/from the ``itk.Image`` class through the :meth:`highdicom.Volume.to_itk` and :meth:`highdicom.Volume.from_itk` methods. Like highdicom, ITK uses the "LPS" convention. However, when converting to and from NumPy arrays, ITK reverses the order of dimensions. This permutation is @@ -42,4 +49,136 @@ Creating a volume from an ITK Image: vol = hd.Volume.from_itk(itk_image) +Use Cases +--------- + +**Load a DICOM Segmentation into ITK:** + +This allows you to work primarily in ITK, while still benefitting from +highdicom's full feature set for reading segmentations, such as correctly +combining multiple segments, filtering segments, filling in missing slices, +lazily retrieving frames, etc. This behavior is not limited to Segmentations, +and would work equally well with any DICOM image the :class:`highdicom.Image` +class supports. The resulting ITK image will have its spatial affine matrix +correctly populated from the source DICOM file. + +.. code-block:: python + + import highdicom as hd + + # Here we load in an example DICOM segmentation from the highdicom repo test + # data that contains two segments. Parameters of the get_volume() method + # control the volume that is extracted. For example, here we choose one of the + # two segments and cast to uint32 before converting to ITK + itk_image = ( + hd.seg.segread("data/test_files/seg_image_ct_binary_overlap.dcm") + .get_volume( + segment_numbers=[2], + relabel=True, + combine_segments=True, + dtype='uint32', + ) + .to_itk() + ) + +**Create a DICOM Segmentation from an ITK image:** + +This allows you to build processing pipelines primarily in ITK, but then +benefit from highdicom's full versatile support for creating DICOM +segmentations, including support for multiple segmentation types, multiple +compression methods (transfer syntaxes) etc. The spatial metadata is correctly +carried throughout this whole snippet and stored in the segmentation without +requiring any user code. + + +.. code-block:: python + + import itk + import highdicom as hd + from pydicom.sr.codedict import codes + from pydicom.uid import JPEGLSLossless + + + # An example CT file in highdicom test data repo + ct_file = "data/test_files/ct_image.dcm" + + # Read the image file using itk + itk_image = itk.imread(ct_file, pixel_type=itk.SS) + + # We do still need to load the source CT image file with highdicom (or + # pydicom) to get the metadata, but we can skip the pixel data + source_image = hd.imread(ct_file, lazy_frame_retrieval=True) + + + def itk_segmentation_method(image: itk.Image) -> itk.Image: + """Toy example of a processing method built with ITK. + + Here we just apply a simple intensity threshold at 300HU as an example. + + """ + thresholder = itk.BinaryThresholdImageFilter[ + type(image), + itk.Image[itk.UC, 3] # produce unsigned 8 bit output image + ].New() + thresholder.SetInput(image) + thresholder.SetLowerThreshold(300) + thresholder.SetInsideValue(1) + thresholder.Update() + return thresholder.GetOutput() + + + # Run the ITK processing + seg_image = itk_segmentation_method(itk_image) + + # Create a highdicom Volume from the ITK image to pass to the + # highdicom.Segmentation constructor + seg_volume = hd.Volume.from_itk(seg_image) + + # In this case, the Volume has shape (128, 128, 1). highdicom.Segmentation + # stores volumes split into frames split down the first dimension, so to + # get the most sensible output we will transpose the volume. We could do + # this explicitly, or, more straightforwardly, just transpose to match to + # the orientation of the source image + seg_volume = seg_volume.match_orientation(source_image.get_volume_geometry()) + + # Now to start creating the segmentation. First we have to describe the segment. + # Here, we describe a bone segment produced by an automatic algorithm) + bone_description = hd.seg.SegmentDescription( + segment_number=1, + segment_label='Bone', + segmented_property_category=codes.SCT.AnatomicalStructure, + segmented_property_type=codes.SCT.Bone, + algorithm_type=hd.seg.SegmentAlgorithmTypeValues.AUTOMATIC, + algorithm_identification=hd.AlgorithmIdentificationSequence( + name='Thresholder 5000X', + version='v1.0', + family=codes.cid7162.MorphologicalOperations + ) + ) + + # Now construct the DICOM Segmentation Image from the ITK segmentation mask. + # Highdicom will take the spatial information from the itk image (via the + # highdicom Volume) and compare it to the source image metadata to establish + # the spatial relationship between the segmentation mask and the source image + # automatically + seg = hd.seg.Segmentation( + source_images=[source_image], + pixel_array=seg_volume, + segmentation_type=hd.seg.SegmentationTypeValues.LABELMAP, + segment_descriptions=[bone_description], + series_instance_uid=hd.UID(), + series_number=1, + sop_instance_uid=hd.UID(), + instance_number=1, + manufacturer='Foo Corp.', + manufacturer_model_name='Bone Segmentation Algorithm', + software_versions='0.0.1', + device_serial_number='1234567890', + transfer_syntax_uid=JPEGLSLossless, + series_description='Bone Threshold Segmentation', + ) + + # Save output file + seg.save_as("bone_segmentation.dcm") + .. _`ITK`: https://itk.org/ diff --git a/docs/sitk_lib.rst b/docs/sitk_lib.rst index 91752149..25822b38 100644 --- a/docs/sitk_lib.rst +++ b/docs/sitk_lib.rst @@ -3,9 +3,16 @@ SimpleITK ========= -`SimpleITK`_ provides a simplified interface for the algorithms and data structures -of `ITK`_. In order to use ITK with highdicom, the ``SimpleITK`` python package must -be installed separately. Version 2.2.1 or later is required. +`SimpleITK`_ is a widely-used library for volumetric image processing. While +many simple processing operations can be performed within highdicom itself, +SimpleITK provides a much larger set of tools and many existing image +processing pipelines are built in SimpleITK. By integrating highdicom with +SimpleITK you can therefore benefit from both ITK's large set of processing +features and highdicom's support for DICOM reading and writing with minimal +"glue code". In particular, correctly transferring spatial metadata between +highdicom and SimpleITK representations is straightforwardly handled for you. +In order to use ITK with highdicom, the ``SimpleITK`` python package must be +installed separately. Version 2.2.1 or later is required. .. _sitk_vol: @@ -46,3 +53,132 @@ Creating a volume from a SimpleITK Image: .. _`ITK`: https://itk.org/ .. _`SimpleITK`: https://simpleitk.org/ + + +Use Cases +--------- + +**Load a DICOM Segmentation into SimpleITK:** + +This allows you to work primarily in SimpleITK, while still benefitting from +highdicom's full feature set for reading segmentations, such as correctly +combining multiple segments, filtering segments, filling in missing slices, +lazily retrieving frames, etc. This behavior is not limited to Segmentations, +and would work equally well with any DICOM image the :class:`highdicom.Image` +class supports. The resulting SimpleITK image will have its spatial affine +matrix correctly populated from the source DICOM file. + +.. code-block:: python + + import highdicom as hd + + # Here we load in an example DICOM segmentation from the highdicom repo test + # data that contains two segments. Parameters of the get_volume() method + # control the volume that is extracted. For example, here we choose one of the + # two segments and cast to uint32 before converting to ITK + sitk_image = ( + hd.seg.segread("data/test_files/seg_image_ct_binary_overlap.dcm") + .get_volume( + segment_numbers=[2], + relabel=True, + combine_segments=True, + dtype='uint32', + ) + .to_simpleitk() + ) + +**Create a DICOM Segmentation from a SimpleITK image:** + +This allows you to build processing pipelines primarily in SimpleITK, but then +benefit from highdicom's full versatile support for creating DICOM +segmentations, including support for multiple segmentation types, multiple +compression methods (transfer syntaxes) etc. The spatial metadata is correctly +carried throughout this whole snippet and stored in the segmentation without +requiring any user code. + + +.. code-block:: python + + import SimpleITK as sitk + import highdicom as hd + from pydicom.sr.codedict import codes + from pydicom.uid import JPEGLSLossless + + + # An example CT file in highdicom test data repo + ct_file = "data/test_files/ct_image.dcm" + + # Read the image file using itk + sitk_image = sitk.ReadImage(ct_file, outputPixelType=sitk.sitkInt16) + + # We do still need to load the source CT image file with highdicom (or + # pydicom) to get the metadata, but we can skip the pixel data + source_image = hd.imread(ct_file, lazy_frame_retrieval=True) + + + def sitk_segmentation_method(image: sitk.Image) -> sitk.Image: + """Toy example of a processing method built with SimpleITK. + + Here we just apply a simple intensity threshold at 300HU as an example. + + """ + thresholder = sitk.BinaryThresholdImageFilter() + thresholder.SetLowerThreshold(300) + thresholder.SetUpperThreshold(10000) + thresholder.SetInsideValue(1) + return thresholder.Execute(image) + + + # Run the ITK processing + seg_image = sitk_segmentation_method(sitk_image) + + # Create a highdicom Volume from the ITK image to pass to the + # highdicom.Segmentation constructor + seg_volume = hd.Volume.from_simpleitk(seg_image) + + # In this case, the Volume has shape (128, 128, 1). highdicom.Segmentation + # stores volumes split into frames split down the first dimension, so to + # get the most sensible output we will transpose the volume. We could do + # this explicitly, or, more straightforwardly, just transpose to match to + # the orientation of the source image + seg_volume = seg_volume.match_orientation(source_image.get_volume_geometry()) + + # Now to start creating the segmentation. First we have to describe the segment + # Here, we describe a bone segment produced by an automatic algorithm) + bone_description = hd.seg.SegmentDescription( + segment_number=1, + segment_label='Bone', + segmented_property_category=codes.SCT.AnatomicalStructure, + segmented_property_type=codes.SCT.Bone, + algorithm_type=hd.seg.SegmentAlgorithmTypeValues.AUTOMATIC, + algorithm_identification=hd.AlgorithmIdentificationSequence( + name='Thresholder 5000X', + version='v1.0', + family=codes.cid7162.MorphologicalOperations + ) + ) + + # Now construct the DICOM Segmentation Image from the ITK segmentation mask. + # Highdicom will take the spatial information from the itk image (via the + # highdicom Volume) and compare it to the source image metadata to establish + # the spatial relationship between the segmentation mask and the source image + # automatically + seg = hd.seg.Segmentation( + source_images=[source_image], + pixel_array=seg_volume, + segmentation_type=hd.seg.SegmentationTypeValues.LABELMAP, + segment_descriptions=[bone_description], + series_instance_uid=hd.UID(), + series_number=1, + sop_instance_uid=hd.UID(), + instance_number=1, + manufacturer='Foo Corp.', + manufacturer_model_name='Bone Segmentation Algorithm', + software_versions='0.0.1', + device_serial_number='1234567890', + transfer_syntax_uid=JPEGLSLossless, + series_description='Bone Threshold Segmentation', + ) + + # Save output file + seg.save_as("bone_segmentation.dcm")