From 4871c51e78db530dfbe0ba3aeef6a0056123b27c Mon Sep 17 00:00:00 2001 From: Yiheng Wang Date: Tue, 21 Jun 2022 21:21:13 +0800 Subject: [PATCH 1/8] implement pydicomreader Signed-off-by: Yiheng Wang --- monai/data/__init__.py | 2 +- monai/data/image_reader.py | 300 ++++++++++++++++++++++++++++++++++++- 2 files changed, 299 insertions(+), 3 deletions(-) diff --git a/monai/data/__init__.py b/monai/data/__init__.py index 293f058acf..6f834294fe 100644 --- a/monai/data/__init__.py +++ b/monai/data/__init__.py @@ -47,7 +47,7 @@ from .folder_layout import FolderLayout from .grid_dataset import GridPatchDataset, PatchDataset, PatchIter, PatchIterd from .image_dataset import ImageDataset -from .image_reader import ImageReader, ITKReader, NibabelReader, NrrdReader, NumpyReader, PILReader +from .image_reader import ImageReader, ITKReader, NibabelReader, NrrdReader, NumpyReader, PILReader, PydicomReader from .image_writer import ( SUPPORTED_WRITERS, ImageWriter, diff --git a/monai/data/image_reader.py b/monai/data/image_reader.py index 13872d5c61..98ba52e537 100644 --- a/monai/data/image_reader.py +++ b/monai/data/image_reader.py @@ -9,6 +9,8 @@ # See the License for the specific language governing permissions and # limitations under the License. +import glob +import os import warnings from abc import ABC, abstractmethod from dataclasses import dataclass @@ -27,22 +29,33 @@ import itk import nibabel as nib import nrrd + import pydicom from nibabel.nifti1 import Nifti1Image from PIL import Image as PILImage - has_nrrd = has_itk = has_nib = has_pil = True + has_nrrd = has_itk = has_nib = has_pil = has_pydicom = True else: itk, has_itk = optional_import("itk", allow_namespace_pkg=True) nib, has_nib = optional_import("nibabel") Nifti1Image, _ = optional_import("nibabel.nifti1", name="Nifti1Image") PILImage, has_pil = optional_import("PIL.Image") + pydicom, has_pydicom = optional_import("pydicom") nrrd, has_nrrd = optional_import("nrrd", allow_namespace_pkg=True) OpenSlide, _ = optional_import("openslide", name="OpenSlide") CuImage, _ = optional_import("cucim", name="CuImage") TiffFile, _ = optional_import("tifffile", name="TiffFile") -__all__ = ["ImageReader", "ITKReader", "NibabelReader", "NumpyReader", "PILReader", "WSIReader", "NrrdReader"] +__all__ = [ + "ImageReader", + "ITKReader", + "NibabelReader", + "NumpyReader", + "PILReader", + "PydicomReader", + "WSIReader", + "NrrdReader", +] class ImageReader(ABC): @@ -354,6 +367,289 @@ def _get_array_data(self, img): return np_img if self.reverse_indexing else np.moveaxis(np_img.T, 0, -1) +@require_pkg(pkg_name="pydicom") +class PydicomReader(ImageReader): + """ + Load medical images based on Pydicom library. + All the supported image formats can be found at: + https://dicom.nema.org/medical/dicom/current/output/chtml/part10/chapter_7.html + The loaded data array will be in C order, for example, a 3D image NumPy + array index order will be `CDWH`. + + This class refers to: + https://nipy.org/nibabel/dicom/dicom_orientation.html#dicom-affine-formula + https://github.com/pydicom/contrib-pydicom/blob/master/input-output/pydicom_series.py + + Args: + channel_dim: the channel dimension of the input image, default is None. + This is used to set original_channel_dim in the metadata, EnsureChannelFirstD reads this field. + If None, `original_channel_dim` will be either `no_channel` or `-1`. + affine_lps_to_ras: whether to convert the affine matrix from "LPS" to "RAS". Defaults to ``True``. + Set to ``True`` to be consistent with ``NibabelReader``, otherwise the affine matrix remains in the ITK convention. + kwargs: additional args for `pydicom.dcmread` API. more details about available args: + https://pydicom.github.io/pydicom/stable/reference/generated/pydicom.filereader.dcmread.html#pydicom.filereader.dcmread + + """ + + def __init__( + self, + channel_dim: Optional[int] = None, + affine_lps_to_ras: bool = True, + **kwargs, + ): + super().__init__() + self.kwargs = kwargs + self.channel_dim = channel_dim + self.affine_lps_to_ras = affine_lps_to_ras + + def verify_suffix(self, filename: Union[Sequence[PathLike], PathLike]) -> bool: + """ + Verify whether the specified file or files format is supported by Pydicom reader. + + Args: + filename: file name or a list of file names to read. + if a list of files, verify all the suffixes. + + """ + return has_pydicom + + def read(self, data: Union[Sequence[PathLike], PathLike], **kwargs): + """ + Read image data from specified file or files, it can read a list of images + and stack them together as multi-channel data in `get_data()`. + If passing directory path instead of file path, will treat it as DICOM images series and read. + + Args: + data: file name or a list of file names to read, + kwargs: additional args for `pydicom.dcmread` API, will override `self.kwargs` for existing keys. + + Returns: + If `data` represents a filename or a directory: return a tuple. + If `data` represents a list of filenames or a list of directories, return a list of tuples. + The tuple is consisted with (data array, metadata). + + """ + img_ = [] + + filenames: Sequence[PathLike] = ensure_tuple(data) + kwargs_ = self.kwargs.copy() + kwargs_.update(kwargs) + for name in filenames: + name = f"{name}" + if Path(name).is_dir(): + # read DICOM series + img_.append(self._read_dicom_slices(name, **kwargs_)) + else: + ds = self._read_single_file(filename=name, **kwargs_) + ds_array = self._get_array_data(ds) + ds_metadata = self._get_meta_dict(ds) + ds_metadata["spatial_shape"] = ds_array.shape + ds_metadata["spacing"] = np.asarray(ds_metadata["00280030"]["Value"]) + img_.append((ds_array, ds_metadata)) + return img_ if len(filenames) > 1 else img_[0] + + def _read_single_file(self, filename: PathLike, **kwargs): + """ + Read image data from specified filename. + + Args: + filename: file name to read. It should also include the path. + kwargs: additional args for `pydicom.dcmread` API. + + """ + ds = pydicom.dcmread(filename, **kwargs) + if hasattr(ds, "pixel_array"): + return ds + elif hasattr(ds, "waveform_array"): + raise NotImplementedError("Waveform data is not implemented.") + elif hasattr(ds, "overlay_array"): + raise NotImplementedError("Overlay data is not implemented.") + else: + raise ValueError("PydicomReader can only read pixel data, waveform data and overlay data.") + + def _read_dicom_slices(self, name: PathLike, **kwargs): + """ + Read a list of dicom slices. Their data arrays will be stacked together at a new dimension as the + last dimension. The stack order depends on Instance Number. The metadata will be produced by the + first slice's metadata, as well as the new spacing. + + Args: + name: the directory that contains DICOM images series. + kwargs: additional args for `pydicom.dcmread` API. + + Returns: + a tuple that consisted with data array and metadata. + + """ + slices = [] + for slc in glob.glob(os.path.join(name, "**")): + slc_ds = self._read_single_file(filename=slc, **kwargs) + if hasattr(slc_ds, "InstanceNumber"): + slices.append(slc_ds) + else: + warnings.warn(f"slice: {slc} in directory {name} does not have InstanceNumber tag, skip it.") + slices = sorted(slices, key=lambda s: s.InstanceNumber) # type: ignore + + if len(slices) == 0: + raise ValueError(f"the directory: {name} does not have valid slices.") + if len(slices) == 1: + return slices[0] + first_slice = slices[0] + average_distance = 0.0 + first_array = self._get_array_data(first_slice) + shape = first_array.shape + spacing, pos = first_slice.PixelSpacing, first_slice.ImagePositionPatient[2] + stack_array = [] + stack_array.append(first_array) + for idx in range(1, len(slices)): + slc = slices[idx] + slc_array = self._get_array_data(slc) + slc_shape = slc_array.shape + slc_spacing, slc_pos = slc.PixelSpacing, slc.ImagePositionPatient[2] # type: ignore + if spacing != slc_spacing: + raise ValueError(f"the directory: {name} contains slices that have different spacings.") + if shape != slc_shape: + raise ValueError(f"the directory: {name} contains slices that have different shapes.") + average_distance += abs(pos - slc_pos) + pos = slc_pos + stack_array.append(slc_array) + average_distance /= len(slices) - 1 + spacing.append(average_distance) + + stack_metadata = self._get_meta_dict(first_slice) + stack_metadata["spacing"] = np.asarray(spacing) + stack_metadata["lastImagePositionPatient"] = slices[-1].ImagePositionPatient + stack_metadata["spatial_shape"] = shape + (len(slices),) + + stack_array = np.stack(stack_array, axis=-1) + + return stack_array, stack_metadata + + def get_data(self, dicom_data): + """ + Extract data array and metadata from loaded image and return them. + This function returns two objects, first is numpy array of image data, second is dict of metadata. + It constructs `affine`, `original_affine`, and `spatial_shape` and stores them in meta dict. + When loading a list of files, they are stacked together at a new dimension as the first dimension, + and the metadata of the first image is used to represent the output metadata. + + Args: + dicom_data: a tuple or a list of tuple. The tuple should consist with data array and metadata. + + """ + img_array: List[np.ndarray] = [] + compatible_meta: Dict = {} + if not isinstance(dicom_data, List): + # for single image, put (data array, metadata) into a list + dicom_data = [dicom_data] + for (data, header) in ensure_tuple(dicom_data): + img_array.append(data) + header["original_affine"] = self._get_affine(header, self.affine_lps_to_ras) + header["affine"] = header["original_affine"].copy() + if self.channel_dim is None: # default to "no_channel" or -1 + header["original_channel_dim"] = "no_channel" if len(data.shape) == len(header["spatial_shape"]) else -1 + else: + header["original_channel_dim"] = self.channel_dim + _copy_compatible_dict(header, compatible_meta) + + return _stack_images(img_array, compatible_meta), compatible_meta + + def _get_meta_dict(self, img) -> Dict: + """ + Get all the metadata of the image and convert to dict type. + + Args: + img: a Pydicom dataset object. + + """ + meta_dict = img.to_json_dict() + # remove Pixel Data "7FE00008" or "7FE00009" or "7FE00010" + # remove Data Set Trailing Padding "FFFCFFFC" + for key in ["7FE00008", "7FE00009", "7FE00010", "FFFCFFFC"]: + if key in meta_dict.keys(): + meta_dict.pop(key) + + return meta_dict # type: ignore + + def _get_affine(self, metadata: Dict, lps_to_ras: bool = True): + """ + Get or construct the affine matrix of the image, it can be used to correct + spacing, orientation or execute spatial transforms. + + Args: + metadata: metadata with dict type. + lps_to_ras: whether to convert the affine matrix from "LPS" to "RAS". Defaults to True. + + """ + affine: np.ndarray = np.eye(4) + rx, ry, rz, cx, cy, cz = metadata["00200037"]["Value"] + sx, sy, sz = metadata["00200032"]["Value"] + dr, dc = metadata["spacing"][:2] + affine[0, 0] = rx * dr + affine[0, 1] = cx * dc + affine[0, 3] = sx + affine[1, 0] = ry * dr + affine[1, 1] = cy * dc + affine[1, 3] = sy + affine[2, 0] = rz * dr + affine[2, 1] = cz * dc + affine[2, 2] = 0 + affine[2, 3] = sz + + # 3d + if "lastImagePositionPatient" in metadata.keys(): + t1n, t2n, t3n = metadata["lastImagePositionPatient"] + n = metadata["spatial_shape"][-1] + k1, k2, k3 = (t1n - sx) / (n - 1), (t2n - sy) / (n - 1), (t3n - sz) / (n - 1) + affine[0, 2] = k1 + affine[1, 2] = k2 + affine[2, 2] = k3 + + if lps_to_ras: + affine = orientation_ras_lps(affine) + return affine + + def _get_spatial_shape(self, img): + """ + Get the spatial shape of `img`. + + Args: + img: a Pydicom dataset object loaded from an image file. + + """ + sr = itk.array_from_matrix(img.GetDirection()).shape[0] + sr = max(min(sr, 3), 1) + _size = list(itk.size(img)) + if self.channel_dim is not None: + _size.pop(self.channel_dim) + return np.asarray(_size[:sr]) + + def _get_array_data(self, img): + """ + Get the array data of the image. If `RescaleSlope` and `RescaleIntercept` are available, the raw array data + will be rescaled. The output data has the type: np.float32 + + Args: + img: a Pydicom dataset object. + + """ + # process Dicom series + data = img.pixel_array.astype(np.float32) + + slope, offset = 1, 0 + rescale_flag = False + if hasattr(img, "RescaleSlope"): + slope = img.RescaleSlope + rescale_flag = True + if hasattr(img, "RescaleIntercept"): + offset = img.RescaleIntercept + rescale_flag = True + if rescale_flag is True: + data = data * slope + offset + + return np.transpose(data, (1, 0)) + + @require_pkg(pkg_name="nibabel") class NibabelReader(ImageReader): """ From 9c76e72b59419b48d3ef2456a56dbe8daa2efbf7 Mon Sep 17 00:00:00 2001 From: monai-bot Date: Tue, 21 Jun 2022 13:43:14 +0000 Subject: [PATCH 2/8] [MONAI] code formatting Signed-off-by: monai-bot --- monai/data/image_reader.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/monai/data/image_reader.py b/monai/data/image_reader.py index 98ba52e537..b5cdd0bd11 100644 --- a/monai/data/image_reader.py +++ b/monai/data/image_reader.py @@ -391,12 +391,7 @@ class PydicomReader(ImageReader): """ - def __init__( - self, - channel_dim: Optional[int] = None, - affine_lps_to_ras: bool = True, - **kwargs, - ): + def __init__(self, channel_dim: Optional[int] = None, affine_lps_to_ras: bool = True, **kwargs): super().__init__() self.kwargs = kwargs self.channel_dim = channel_dim From 068af3b262addf1112a4bc76202b3a292699fdef Mon Sep 17 00:00:00 2001 From: Yiheng Wang Date: Wed, 22 Jun 2022 23:42:53 +0800 Subject: [PATCH 3/8] modify read function to return obj Signed-off-by: Yiheng Wang --- monai/data/image_reader.py | 122 ++++++++++++++++++----------------- monai/transforms/io/array.py | 11 +++- 2 files changed, 72 insertions(+), 61 deletions(-) diff --git a/monai/data/image_reader.py b/monai/data/image_reader.py index b5cdd0bd11..80c30f4a3b 100644 --- a/monai/data/image_reader.py +++ b/monai/data/image_reader.py @@ -385,7 +385,7 @@ class PydicomReader(ImageReader): This is used to set original_channel_dim in the metadata, EnsureChannelFirstD reads this field. If None, `original_channel_dim` will be either `no_channel` or `-1`. affine_lps_to_ras: whether to convert the affine matrix from "LPS" to "RAS". Defaults to ``True``. - Set to ``True`` to be consistent with ``NibabelReader``, otherwise the affine matrix remains in the ITK convention. + Set to ``True`` to be consistent with ``NibabelReader``, otherwise the affine matrix remains in the Dicom convention. kwargs: additional args for `pydicom.dcmread` API. more details about available args: https://pydicom.github.io/pydicom/stable/reference/generated/pydicom.filereader.dcmread.html#pydicom.filereader.dcmread @@ -419,9 +419,9 @@ def read(self, data: Union[Sequence[PathLike], PathLike], **kwargs): kwargs: additional args for `pydicom.dcmread` API, will override `self.kwargs` for existing keys. Returns: - If `data` represents a filename or a directory: return a tuple. - If `data` represents a list of filenames or a list of directories, return a list of tuples. - The tuple is consisted with (data array, metadata). + If `data` represents a filename: return a pydicom dataset object. + If `data` represents a list of filenames or a directory: return a list of pydicom dataset object. + If `data` represents a list of directories: return a list of list of pydicom dataset object. """ img_ = [] @@ -429,66 +429,45 @@ def read(self, data: Union[Sequence[PathLike], PathLike], **kwargs): filenames: Sequence[PathLike] = ensure_tuple(data) kwargs_ = self.kwargs.copy() kwargs_.update(kwargs) + + self.has_series = False + for name in filenames: name = f"{name}" if Path(name).is_dir(): # read DICOM series - img_.append(self._read_dicom_slices(name, **kwargs_)) + slices = [] + for slc in glob.glob(os.path.join(name, "**")): + slices.append(pydicom.dcmread(fp=slc, **kwargs_)) + img_.append(slices if len(slices) > 1 else slices[0]) + self.has_series = True else: - ds = self._read_single_file(filename=name, **kwargs_) - ds_array = self._get_array_data(ds) - ds_metadata = self._get_meta_dict(ds) - ds_metadata["spatial_shape"] = ds_array.shape - ds_metadata["spacing"] = np.asarray(ds_metadata["00280030"]["Value"]) - img_.append((ds_array, ds_metadata)) + ds = pydicom.dcmread(fp=name, **kwargs_) + img_.append(ds) return img_ if len(filenames) > 1 else img_[0] - def _read_single_file(self, filename: PathLike, **kwargs): + def _combine_dicom_series(self, data): """ - Read image data from specified filename. - - Args: - filename: file name to read. It should also include the path. - kwargs: additional args for `pydicom.dcmread` API. - - """ - ds = pydicom.dcmread(filename, **kwargs) - if hasattr(ds, "pixel_array"): - return ds - elif hasattr(ds, "waveform_array"): - raise NotImplementedError("Waveform data is not implemented.") - elif hasattr(ds, "overlay_array"): - raise NotImplementedError("Overlay data is not implemented.") - else: - raise ValueError("PydicomReader can only read pixel data, waveform data and overlay data.") - - def _read_dicom_slices(self, name: PathLike, **kwargs): - """ - Read a list of dicom slices. Their data arrays will be stacked together at a new dimension as the + Combine dicom series. Their data arrays will be stacked together at a new dimension as the last dimension. The stack order depends on Instance Number. The metadata will be produced by the first slice's metadata, as well as the new spacing. - Args: - name: the directory that contains DICOM images series. - kwargs: additional args for `pydicom.dcmread` API. - + data: a list of pydicom dataset objects. Returns: a tuple that consisted with data array and metadata. - """ slices = [] - for slc in glob.glob(os.path.join(name, "**")): - slc_ds = self._read_single_file(filename=slc, **kwargs) + # for a dicom s + for slc_ds in data: if hasattr(slc_ds, "InstanceNumber"): slices.append(slc_ds) else: - warnings.warn(f"slice: {slc} in directory {name} does not have InstanceNumber tag, skip it.") - slices = sorted(slices, key=lambda s: s.InstanceNumber) # type: ignore + warnings.warn(f"slice: {slc_ds[0x0008, 0x0018]} does not have InstanceNumber tag, skip it.") + slices = sorted(slices, key=lambda s: s.InstanceNumber) if len(slices) == 0: - raise ValueError(f"the directory: {name} does not have valid slices.") - if len(slices) == 1: - return slices[0] + raise ValueError("the input does not have valid slices.") + first_slice = slices[0] average_distance = 0.0 first_array = self._get_array_data(first_slice) @@ -500,27 +479,32 @@ def _read_dicom_slices(self, name: PathLike, **kwargs): slc = slices[idx] slc_array = self._get_array_data(slc) slc_shape = slc_array.shape - slc_spacing, slc_pos = slc.PixelSpacing, slc.ImagePositionPatient[2] # type: ignore + slc_spacing, slc_pos = slc.PixelSpacing, slc.ImagePositionPatient[2] if spacing != slc_spacing: - raise ValueError(f"the directory: {name} contains slices that have different spacings.") + raise ValueError("the list contains slices that have different spacings.") if shape != slc_shape: - raise ValueError(f"the directory: {name} contains slices that have different shapes.") + raise ValueError("the list contains slices that have different shapes.") average_distance += abs(pos - slc_pos) pos = slc_pos stack_array.append(slc_array) - average_distance /= len(slices) - 1 - spacing.append(average_distance) - stack_metadata = self._get_meta_dict(first_slice) - stack_metadata["spacing"] = np.asarray(spacing) - stack_metadata["lastImagePositionPatient"] = slices[-1].ImagePositionPatient - stack_metadata["spatial_shape"] = shape + (len(slices),) - - stack_array = np.stack(stack_array, axis=-1) + if len(slices) > 1: + average_distance /= len(slices) - 1 + spacing.append(average_distance) + stack_array = np.stack(stack_array, axis=-1) + stack_metadata = self._get_meta_dict(first_slice) + stack_metadata["spacing"] = np.asarray(spacing) + stack_metadata["lastImagePositionPatient"] = np.asarray(slices[-1].ImagePositionPatient) + stack_metadata["spatial_shape"] = shape + (len(slices),) + else: + stack_array = stack_array[0] + stack_metadata = self._get_meta_dict(first_slice) + stack_metadata["spacing"] = np.asarray(spacing) + stack_metadata["spatial_shape"] = shape return stack_array, stack_metadata - def get_data(self, dicom_data): + def get_data(self, data): """ Extract data array and metadata from loaded image and return them. This function returns two objects, first is numpy array of image data, second is dict of metadata. @@ -529,14 +513,32 @@ def get_data(self, dicom_data): and the metadata of the first image is used to represent the output metadata. Args: - dicom_data: a tuple or a list of tuple. The tuple should consist with data array and metadata. + data: a pydicom dataset object, or a list of pydicom dataset objects, or a list of list of + pydicom dataset objects. """ + + dicom_data = [] + if self.has_series is True: + # a list, or a list of list + if not isinstance(data[0], List): + dicom_data.append(self._combine_dicom_series(data)) + else: + for series in data: + dicom_data.append(series) + else: + if not isinstance(data, List): + data = [data] + for d in data: + data_array = self._get_array_data(d) + metadata = self._get_meta_dict(d) + metadata["spatial_shape"] = data_array.shape + metadata["spacing"] = np.asarray(metadata["00280030"]["Value"]) + dicom_data.append((data_array, metadata)) + img_array: List[np.ndarray] = [] compatible_meta: Dict = {} - if not isinstance(dicom_data, List): - # for single image, put (data array, metadata) into a list - dicom_data = [dicom_data] + for (data, header) in ensure_tuple(dicom_data): img_array.append(data) header["original_affine"] = self._get_affine(header, self.affine_lps_to_ras) diff --git a/monai/transforms/io/array.py b/monai/transforms/io/array.py index 87f9c6676d..6e2f32d156 100644 --- a/monai/transforms/io/array.py +++ b/monai/transforms/io/array.py @@ -28,7 +28,15 @@ from monai.config import DtypeLike, NdarrayOrTensor, PathLike from monai.data import image_writer from monai.data.folder_layout import FolderLayout -from monai.data.image_reader import ImageReader, ITKReader, NibabelReader, NrrdReader, NumpyReader, PILReader +from monai.data.image_reader import ( + ImageReader, + ITKReader, + NibabelReader, + NrrdReader, + NumpyReader, + PILReader, + PydicomReader, +) from monai.transforms.transform import Transform from monai.transforms.utility.array import EnsureChannelFirst from monai.utils import GridSampleMode, GridSamplePadMode @@ -47,6 +55,7 @@ "numpyreader": NumpyReader, "pilreader": PILReader, "nibabelreader": NibabelReader, + "pydicomreader": PydicomReader, } From 8b5c473a2ba7bc6fc48ca26cb1a7712b4db05414 Mon Sep 17 00:00:00 2001 From: Yiheng Wang Date: Thu, 23 Jun 2022 23:33:11 +0800 Subject: [PATCH 4/8] enhance docstrings and fix errors Signed-off-by: Yiheng Wang --- monai/data/image_reader.py | 92 ++++++++++++++++++++++---------------- tests/test_init_reader.py | 8 +++- 2 files changed, 60 insertions(+), 40 deletions(-) diff --git a/monai/data/image_reader.py b/monai/data/image_reader.py index 80c30f4a3b..dfdb23593d 100644 --- a/monai/data/image_reader.py +++ b/monai/data/image_reader.py @@ -373,8 +373,6 @@ class PydicomReader(ImageReader): Load medical images based on Pydicom library. All the supported image formats can be found at: https://dicom.nema.org/medical/dicom/current/output/chtml/part10/chapter_7.html - The loaded data array will be in C order, for example, a 3D image NumPy - array index order will be `CDWH`. This class refers to: https://nipy.org/nibabel/dicom/dicom_orientation.html#dicom-affine-formula @@ -417,6 +415,10 @@ def read(self, data: Union[Sequence[PathLike], PathLike], **kwargs): Args: data: file name or a list of file names to read, kwargs: additional args for `pydicom.dcmread` API, will override `self.kwargs` for existing keys. + If the `get_data` function will be called later, please ensure that the argument `stop_before_pixels` + is `True`, and `specific_tags` covers all necessary tags, + such as `PixelSpacing`, `ImagePositionPatient`, `ImageOrientationPatient` and all `pixel_array` + related tags. Returns: If `data` represents a filename: return a pydicom dataset object. @@ -448,21 +450,33 @@ def read(self, data: Union[Sequence[PathLike], PathLike], **kwargs): def _combine_dicom_series(self, data): """ - Combine dicom series. Their data arrays will be stacked together at a new dimension as the - last dimension. The stack order depends on Instance Number. The metadata will be produced by the - first slice's metadata, as well as the new spacing. + Combine dicom series (a list of pydicom dataset objects). Their data arrays will be stacked together at a new + dimension as the last dimension. + + The stack order depends on Instance Number. The metadata will be based on the + first slice's metadata, and some new items will be added: + + "spacing": the new spacing of the stacked volume. + "lastImagePositionPatient": `ImagePositionPatient` for the last slice, it will be used to achieve the affine + matrix. + "spatial_shape": the spatial shape of the stacked volume. + Args: data: a list of pydicom dataset objects. Returns: a tuple that consisted with data array and metadata. """ slices = [] - # for a dicom s + # for a dicom series for slc_ds in data: + if not hasattr(slc_ds, "PixelSpacing"): + raise ValueError(f"slice: {slc_ds.filename} does not have PixelSpacing tag.") + if not hasattr(slc_ds, "ImagePositionPatient"): + raise ValueError(f"slice: {slc_ds.filename} does not have ImagePositionPatient tag.") if hasattr(slc_ds, "InstanceNumber"): slices.append(slc_ds) else: - warnings.warn(f"slice: {slc_ds[0x0008, 0x0018]} does not have InstanceNumber tag, skip it.") + warnings.warn(f"slice: {slc_ds.filename} does not have InstanceNumber tag, skip it.") slices = sorted(slices, key=lambda s: s.InstanceNumber) if len(slices) == 0: @@ -512,6 +526,9 @@ def get_data(self, data): When loading a list of files, they are stacked together at a new dimension as the first dimension, and the metadata of the first image is used to represent the output metadata. + To use this function, all pydicom dataset objects should contain: `pixel_array`, `PixelSpacing`, + `ImagePositionPatient` and `ImageOrientationPatient`. + Args: data: a pydicom dataset object, or a list of pydicom dataset objects, or a list of list of pydicom dataset objects. @@ -519,14 +536,17 @@ def get_data(self, data): """ dicom_data = [] + # dicom series if self.has_series is True: - # a list, or a list of list + # a list, all objects within a list belong to one dicom series if not isinstance(data[0], List): dicom_data.append(self._combine_dicom_series(data)) + # a list of list, each inner list represents a dicom series else: for series in data: - dicom_data.append(series) + dicom_data.append(self._combine_dicom_series(series)) else: + # a single pydicom dataset object if not isinstance(data, List): data = [data] for d in data: @@ -539,15 +559,17 @@ def get_data(self, data): img_array: List[np.ndarray] = [] compatible_meta: Dict = {} - for (data, header) in ensure_tuple(dicom_data): - img_array.append(data) - header["original_affine"] = self._get_affine(header, self.affine_lps_to_ras) - header["affine"] = header["original_affine"].copy() + for (data_array, metadata) in ensure_tuple(dicom_data): + img_array.append(data_array) + metadata["original_affine"] = self._get_affine(metadata, self.affine_lps_to_ras) + metadata["affine"] = metadata["original_affine"].copy() if self.channel_dim is None: # default to "no_channel" or -1 - header["original_channel_dim"] = "no_channel" if len(data.shape) == len(header["spatial_shape"]) else -1 + metadata["original_channel_dim"] = ( + "no_channel" if len(data_array.shape) == len(metadata["spatial_shape"]) else -1 + ) else: - header["original_channel_dim"] = self.channel_dim - _copy_compatible_dict(header, compatible_meta) + metadata["original_channel_dim"] = self.channel_dim + _copy_compatible_dict(metadata, compatible_meta) return _stack_images(img_array, compatible_meta), compatible_meta @@ -559,6 +581,11 @@ def _get_meta_dict(self, img) -> Dict: img: a Pydicom dataset object. """ + if not hasattr(img, "ImagePositionPatient"): + raise ValueError(f"dicom data: {img.filename} does not have ImagePositionPatient.") + if not hasattr(img, "ImageOrientationPatient"): + raise ValueError(f"dicom data: {img.filename} does not have ImageOrientationPatient.") + meta_dict = img.to_json_dict() # remove Pixel Data "7FE00008" or "7FE00009" or "7FE00010" # remove Data Set Trailing Padding "FFFCFFFC" @@ -579,17 +606,19 @@ def _get_affine(self, metadata: Dict, lps_to_ras: bool = True): """ affine: np.ndarray = np.eye(4) + # "00200037" is the tag of `ImageOrientationPatient` rx, ry, rz, cx, cy, cz = metadata["00200037"]["Value"] + # "00200032" is the tag of `ImagePositionPatient` sx, sy, sz = metadata["00200032"]["Value"] dr, dc = metadata["spacing"][:2] - affine[0, 0] = rx * dr - affine[0, 1] = cx * dc + affine[0, 0] = cx * dr + affine[0, 1] = rx * dc affine[0, 3] = sx - affine[1, 0] = ry * dr - affine[1, 1] = cy * dc + affine[1, 0] = cy * dr + affine[1, 1] = ry * dc affine[1, 3] = sy - affine[2, 0] = rz * dr - affine[2, 1] = cz * dc + affine[2, 0] = cz * dr + affine[2, 1] = rz * dc affine[2, 2] = 0 affine[2, 3] = sz @@ -606,21 +635,6 @@ def _get_affine(self, metadata: Dict, lps_to_ras: bool = True): affine = orientation_ras_lps(affine) return affine - def _get_spatial_shape(self, img): - """ - Get the spatial shape of `img`. - - Args: - img: a Pydicom dataset object loaded from an image file. - - """ - sr = itk.array_from_matrix(img.GetDirection()).shape[0] - sr = max(min(sr, 3), 1) - _size = list(itk.size(img)) - if self.channel_dim is not None: - _size.pop(self.channel_dim) - return np.asarray(_size[:sr]) - def _get_array_data(self, img): """ Get the array data of the image. If `RescaleSlope` and `RescaleIntercept` are available, the raw array data @@ -631,6 +645,8 @@ def _get_array_data(self, img): """ # process Dicom series + if not hasattr(img, "pixel_array"): + raise ValueError(f"dicom data: {img.filename} does not have pixel_array.") data = img.pixel_array.astype(np.float32) slope, offset = 1, 0 @@ -644,7 +660,7 @@ def _get_array_data(self, img): if rescale_flag is True: data = data * slope + offset - return np.transpose(data, (1, 0)) + return data @require_pkg(pkg_name="nibabel") diff --git a/tests/test_init_reader.py b/tests/test_init_reader.py index df055e571c..e9a8ec96f7 100644 --- a/tests/test_init_reader.py +++ b/tests/test_init_reader.py @@ -11,7 +11,7 @@ import unittest -from monai.data import ITKReader, NibabelReader, NrrdReader, NumpyReader, PILReader +from monai.data import ITKReader, NibabelReader, NrrdReader, NumpyReader, PILReader, PydicomReader from monai.transforms import LoadImage, LoadImaged from tests.utils import SkipIfNoModule @@ -23,7 +23,7 @@ def test_load_image(self): self.assertIsInstance(instance1, LoadImage) self.assertIsInstance(instance2, LoadImage) - for r in ["NibabelReader", "PILReader", "ITKReader", "NumpyReader", "NrrdReader", None]: + for r in ["NibabelReader", "PILReader", "ITKReader", "NumpyReader", "NrrdReader", "PydicomReader", None]: inst = LoadImaged("image", reader=r) self.assertIsInstance(inst, LoadImaged) @@ -31,6 +31,7 @@ def test_load_image(self): @SkipIfNoModule("nibabel") @SkipIfNoModule("PIL") @SkipIfNoModule("nrrd") + @SkipIfNoModule("Pydicom") def test_readers(self): inst = ITKReader() self.assertIsInstance(inst, ITKReader) @@ -40,6 +41,9 @@ def test_readers(self): inst = NibabelReader(as_closest_canonical=True) self.assertIsInstance(inst, NibabelReader) + inst = PydicomReader() + self.assertIsInstance(inst, PydicomReader) + inst = NumpyReader() self.assertIsInstance(inst, NumpyReader) inst = NumpyReader(npz_keys="test") From f590247c524bc89208c41804dbaf8ff9241288d4 Mon Sep 17 00:00:00 2001 From: Yiheng Wang Date: Fri, 24 Jun 2022 13:34:56 +0800 Subject: [PATCH 5/8] add unittest Signed-off-by: Yiheng Wang --- monai/data/image_reader.py | 22 ++++++------ monai/transforms/io/array.py | 4 +-- requirements-dev.txt | 1 + tests/test_load_image.py | 65 +++++++++++++++++++++++++++++++++++- 4 files changed, 79 insertions(+), 13 deletions(-) diff --git a/monai/data/image_reader.py b/monai/data/image_reader.py index dfdb23593d..f578b896b6 100644 --- a/monai/data/image_reader.py +++ b/monai/data/image_reader.py @@ -383,9 +383,14 @@ class PydicomReader(ImageReader): This is used to set original_channel_dim in the metadata, EnsureChannelFirstD reads this field. If None, `original_channel_dim` will be either `no_channel` or `-1`. affine_lps_to_ras: whether to convert the affine matrix from "LPS" to "RAS". Defaults to ``True``. - Set to ``True`` to be consistent with ``NibabelReader``, otherwise the affine matrix remains in the Dicom convention. + Set to ``True`` to be consistent with ``NibabelReader``, + otherwise the affine matrix remains in the Dicom convention. kwargs: additional args for `pydicom.dcmread` API. more details about available args: https://pydicom.github.io/pydicom/stable/reference/generated/pydicom.filereader.dcmread.html#pydicom.filereader.dcmread + If the `get_data` function will be called + (for example, when using this reader with `monai.transforms.LoadImage`), please ensure that the argument + `stop_before_pixels` is `True`, and `specific_tags` covers all necessary tags, such as `PixelSpacing`, + `ImagePositionPatient`, `ImageOrientationPatient` and all `pixel_array` related tags. """ @@ -415,10 +420,6 @@ def read(self, data: Union[Sequence[PathLike], PathLike], **kwargs): Args: data: file name or a list of file names to read, kwargs: additional args for `pydicom.dcmread` API, will override `self.kwargs` for existing keys. - If the `get_data` function will be called later, please ensure that the argument `stop_before_pixels` - is `True`, and `specific_tags` covers all necessary tags, - such as `PixelSpacing`, `ImagePositionPatient`, `ImageOrientationPatient` and all `pixel_array` - related tags. Returns: If `data` represents a filename: return a pydicom dataset object. @@ -456,10 +457,10 @@ def _combine_dicom_series(self, data): The stack order depends on Instance Number. The metadata will be based on the first slice's metadata, and some new items will be added: - "spacing": the new spacing of the stacked volume. + "spacing": the new spacing of the stacked slices. "lastImagePositionPatient": `ImagePositionPatient` for the last slice, it will be used to achieve the affine matrix. - "spatial_shape": the spatial shape of the stacked volume. + "spatial_shape": the spatial shape of the stacked slices. Args: data: a list of pydicom dataset objects. @@ -523,8 +524,9 @@ def get_data(self, data): Extract data array and metadata from loaded image and return them. This function returns two objects, first is numpy array of image data, second is dict of metadata. It constructs `affine`, `original_affine`, and `spatial_shape` and stores them in meta dict. - When loading a list of files, they are stacked together at a new dimension as the first dimension, - and the metadata of the first image is used to represent the output metadata. + For dicom series within the input, all slices will be stacked first, + When loading a list of files (dicom file, or stacked dicom series), they are stacked together at a new + dimension as the first dimension, and the metadata of the first image is used to represent the output metadata. To use this function, all pydicom dataset objects should contain: `pixel_array`, `PixelSpacing`, `ImagePositionPatient` and `ImageOrientationPatient`. @@ -536,7 +538,7 @@ def get_data(self, data): """ dicom_data = [] - # dicom series + # combine dicom series if exists if self.has_series is True: # a list, all objects within a list belong to one dicom series if not isinstance(data[0], List): diff --git a/monai/transforms/io/array.py b/monai/transforms/io/array.py index 6e2f32d156..811c01b88e 100644 --- a/monai/transforms/io/array.py +++ b/monai/transforms/io/array.py @@ -50,12 +50,12 @@ __all__ = ["LoadImage", "SaveImage", "SUPPORTED_READERS"] SUPPORTED_READERS = { + "pydicomreader": PydicomReader, "itkreader": ITKReader, "nrrdreader": NrrdReader, "numpyreader": NumpyReader, "pilreader": PILReader, "nibabelreader": NibabelReader, - "pydicomreader": PydicomReader, } @@ -119,7 +119,7 @@ def __init__( - if `reader` is None, a default set of `SUPPORTED_READERS` will be used. - if `reader` is a string, it's treated as a class name or dotted path (such as ``"monai.data.ITKReader"``), the supported built-in reader classes are - ``"ITKReader"``, ``"NibabelReader"``, ``"NumpyReader"``. + ``"ITKReader"``, ``"NibabelReader"``, ``"NumpyReader"``, ``"PydicomReader"``. a reader instance will be constructed with the `*args` and `**kwargs` parameters. - if `reader` is a reader class/instance, it will be registered to this loader accordingly. image_only: if True return only the image volume, otherwise return image data array and header dict. diff --git a/requirements-dev.txt b/requirements-dev.txt index df1fe67f5e..994c9c67be 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -49,3 +49,4 @@ fire jsonschema pynrrd pre-commit +pydicom==2.3.0 diff --git a/tests/test_load_image.py b/tests/test_load_image.py index 201fe2fd5b..ea1adba40c 100644 --- a/tests/test_load_image.py +++ b/tests/test_load_image.py @@ -20,7 +20,7 @@ from parameterized import parameterized from PIL import Image -from monai.data import ITKReader, NibabelReader +from monai.data import ITKReader, NibabelReader, PydicomReader from monai.transforms import LoadImage @@ -134,6 +134,31 @@ def get_data(self, _obj): (128, 128, 3, 128), ] +# test same dicom data with PydicomReader +TEST_CASE_19 = [ + {"image_only": False, "reader": PydicomReader()}, + "tests/testing_data/CT_DICOM", + (16, 16, 4), + (16, 16, 4), +] + +TEST_CASE_20 = [ + {"image_only": False, "reader": "PydicomReader", "ensure_channel_first": True}, + "tests/testing_data/CT_DICOM", + (16, 16, 4), + (1, 16, 16, 4), +] + +TEST_CASE_21 = [ + {"image_only": False, "reader": "PydicomReader", "affine_lps_to_ras": True, "defer_size": "2 MB"}, + "tests/testing_data/CT_DICOM", + (16, 16, 4), + (16, 16, 4), +] + +# test reader consistency between PydicomReader and ITKReader on dicom data +TEST_CASE_22 = ["tests/testing_data/CT_DICOM"] + class TestLoadImage(unittest.TestCase): @parameterized.expand( @@ -208,6 +233,44 @@ def test_itk_reader_multichannel(self): np.testing.assert_allclose(result[:, :, 1], test_image[:, :, 1]) np.testing.assert_allclose(result[:, :, 2], test_image[:, :, 2]) + @parameterized.expand([TEST_CASE_19, TEST_CASE_20, TEST_CASE_21]) + def test_pydicom_dicom_series_reader(self, input_param, filenames, expected_shape, expected_np_shape): + result, header = LoadImage(**input_param)(filenames) + self.assertTrue("affine" in header) + self.assertEqual(header["filename_or_obj"], f"{Path(filenames)}") + np.testing.assert_allclose( + header["affine"], + np.array( + [ + [0.0, -0.488281, 0.0, 125.0], + [-0.488281, 0.0, 0.0, 128.100006], + [0.0, 0.0, 68.33333333, -99.480003], + [0.0, 0.0, 0.0, 1.0], + ] + ), + ) + self.assertTupleEqual(tuple(header["spatial_shape"]), expected_shape) + self.assertTupleEqual(result.shape, expected_np_shape) + + @parameterized.expand([TEST_CASE_22]) + def test_dicom_reader_consistency(self, filenames): + itk_param = {"reader": "ITKReader"} + pydicom_param = {"reader": "PydicomReader"} + for affine_flag in [True, False]: + itk_param["affine_lps_to_ras"] = affine_flag + pydicom_param["affine_lps_to_ras"] = affine_flag + itk_result, itk_header = LoadImage(**itk_param)(filenames) + pydicom_result, pydicom_header = LoadImage(**pydicom_param)(filenames) + + # the row and column order is different + # see also: + # https://nipy.org/nibabel/dicom/dicom_orientation.html#i-j-columns-rows-in-dicom + + self.assertTrue(np.all(pydicom_result == np.transpose(itk_result, [1, 0, 2]))) + + r_c_transform = np.array([[0, 1, 0, 0], [1, 0, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]) + np.testing.assert_allclose(itk_header["affine"], pydicom_header["affine"] @ r_c_transform) + def test_load_nifti_multichannel(self): test_image = np.random.randint(0, 256, size=(31, 64, 16, 2)).astype(np.float32) with tempfile.TemporaryDirectory() as tempdir: From 0564e1ed026078a4eee4c2e1f0e88e48bc9ad9f0 Mon Sep 17 00:00:00 2001 From: Wenqi Li Date: Sun, 26 Jun 2022 16:20:56 +0100 Subject: [PATCH 6/8] adds deps Signed-off-by: Wenqi Li --- docs/requirements.txt | 1 + docs/source/installation.md | 4 ++-- environment-dev.yml | 1 + requirements-dev.txt | 2 +- setup.cfg | 3 +++ 5 files changed, 8 insertions(+), 3 deletions(-) diff --git a/docs/requirements.txt b/docs/requirements.txt index b7edff27fa..50c3f07b5d 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -29,3 +29,4 @@ pyyaml fire jsonschema pynrrd +pydicom diff --git a/docs/source/installation.md b/docs/source/installation.md index 393205cd1b..8623faee21 100644 --- a/docs/source/installation.md +++ b/docs/source/installation.md @@ -190,9 +190,9 @@ Since MONAI v0.2.0, the extras syntax such as `pip install 'monai[nibabel]'` is - The options are ``` -[nibabel, skimage, pillow, tensorboard, gdown, ignite, torchvision, itk, tqdm, lmdb, psutil, cucim, openslide, pandas, einops, transformers, mlflow, matplotlib, tensorboardX, tifffile, imagecodecs, pyyaml, fire, jsonschema, pynrrd] +[nibabel, skimage, pillow, tensorboard, gdown, ignite, torchvision, itk, tqdm, lmdb, psutil, cucim, openslide, pandas, einops, transformers, mlflow, matplotlib, tensorboardX, tifffile, imagecodecs, pyyaml, fire, jsonschema, pynrrd, pydicom] ``` which correspond to `nibabel`, `scikit-image`, `pillow`, `tensorboard`, -`gdown`, `pytorch-ignite`, `torchvision`, `itk`, `tqdm`, `lmdb`, `psutil`, `cucim`, `openslide-python`, `pandas`, `einops`, `transformers`, `mlflow`, `matplotlib`, `tensorboardX`, `tifffile`, `imagecodecs`, `pyyaml`, `fire`, `jsonschema`, `pynrrd`, respectively. +`gdown`, `pytorch-ignite`, `torchvision`, `itk`, `tqdm`, `lmdb`, `psutil`, `cucim`, `openslide-python`, `pandas`, `einops`, `transformers`, `mlflow`, `matplotlib`, `tensorboardX`, `tifffile`, `imagecodecs`, `pyyaml`, `fire`, `jsonschema`, `pynrrd`, `pydicom`, respectively. - `pip install 'monai[all]'` installs all the optional dependencies. diff --git a/environment-dev.yml b/environment-dev.yml index dc9a0970cd..e71f434dc4 100644 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -46,6 +46,7 @@ dependencies: - fire - jsonschema - pynrrd + - pydicom - pip - pip: # pip for itk as conda-forge version only up to v5.1 diff --git a/requirements-dev.txt b/requirements-dev.txt index 994c9c67be..0c043ec420 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -49,4 +49,4 @@ fire jsonschema pynrrd pre-commit -pydicom==2.3.0 +pydicom diff --git a/setup.cfg b/setup.cfg index a25a64797c..d10ca8ff49 100644 --- a/setup.cfg +++ b/setup.cfg @@ -54,6 +54,7 @@ all = fire jsonschema pynrrd + pydicom nibabel = nibabel skimage = @@ -104,6 +105,8 @@ jsonschema = jsonschema pynrrd = pynrrd +pydicom = + pydicom [flake8] select = B,C,E,F,N,P,T4,W,B9 From f01a7cf0340ff2538ed4c3b23436e39995dbd25c Mon Sep 17 00:00:00 2001 From: Wenqi Li Date: Sun, 26 Jun 2022 18:13:33 +0100 Subject: [PATCH 7/8] optional meta Signed-off-by: Wenqi Li --- monai/data/image_reader.py | 45 ++++++++++++++++++-------------------- 1 file changed, 21 insertions(+), 24 deletions(-) diff --git a/monai/data/image_reader.py b/monai/data/image_reader.py index f578b896b6..d6dd0188bc 100644 --- a/monai/data/image_reader.py +++ b/monai/data/image_reader.py @@ -439,9 +439,7 @@ def read(self, data: Union[Sequence[PathLike], PathLike], **kwargs): name = f"{name}" if Path(name).is_dir(): # read DICOM series - slices = [] - for slc in glob.glob(os.path.join(name, "**")): - slices.append(pydicom.dcmread(fp=slc, **kwargs_)) + slices = [pydicom.dcmread(fp=slc, **kwargs_) for slc in glob.glob(os.path.join(name, "**"))] img_.append(slices if len(slices) > 1 else slices[0]) self.has_series = True else: @@ -470,10 +468,6 @@ def _combine_dicom_series(self, data): slices = [] # for a dicom series for slc_ds in data: - if not hasattr(slc_ds, "PixelSpacing"): - raise ValueError(f"slice: {slc_ds.filename} does not have PixelSpacing tag.") - if not hasattr(slc_ds, "ImagePositionPatient"): - raise ValueError(f"slice: {slc_ds.filename} does not have ImagePositionPatient tag.") if hasattr(slc_ds, "InstanceNumber"): slices.append(slc_ds) else: @@ -487,18 +481,18 @@ def _combine_dicom_series(self, data): average_distance = 0.0 first_array = self._get_array_data(first_slice) shape = first_array.shape - spacing, pos = first_slice.PixelSpacing, first_slice.ImagePositionPatient[2] - stack_array = [] - stack_array.append(first_array) + spacing = getattr(first_slice, "PixelSpacing", (1.0, 1.0, 1.0)) + pos = getattr(first_slice, "ImagePositionPatient", (0.0, 0.0, 0.0))[2] + stack_array = [first_array] for idx in range(1, len(slices)): - slc = slices[idx] - slc_array = self._get_array_data(slc) + slc_array = self._get_array_data(slices[idx]) slc_shape = slc_array.shape - slc_spacing, slc_pos = slc.PixelSpacing, slc.ImagePositionPatient[2] + slc_spacing = getattr(first_slice, "PixelSpacing", (1.0, 1.0, 1.0)) + slc_pos = getattr(first_slice, "ImagePositionPatient", (0.0, 0.0, float(idx)))[2] if spacing != slc_spacing: - raise ValueError("the list contains slices that have different spacings.") + warnings.warn(f"the list contains slices that have different spacings {spacing} and {slc_spacing}.") if shape != slc_shape: - raise ValueError("the list contains slices that have different shapes.") + warnings.warn(f"the list contains slices that have different shapes {shape} and {slc_shape}.") average_distance += abs(pos - slc_pos) pos = slc_pos stack_array.append(slc_array) @@ -509,7 +503,8 @@ def _combine_dicom_series(self, data): stack_array = np.stack(stack_array, axis=-1) stack_metadata = self._get_meta_dict(first_slice) stack_metadata["spacing"] = np.asarray(spacing) - stack_metadata["lastImagePositionPatient"] = np.asarray(slices[-1].ImagePositionPatient) + if hasattr(slices[-1], "ImagePositionPatient"): + stack_metadata["lastImagePositionPatient"] = np.asarray(slices[-1].ImagePositionPatient) stack_metadata["spatial_shape"] = shape + (len(slices),) else: stack_array = stack_array[0] @@ -555,7 +550,7 @@ def get_data(self, data): data_array = self._get_array_data(d) metadata = self._get_meta_dict(d) metadata["spatial_shape"] = data_array.shape - metadata["spacing"] = np.asarray(metadata["00280030"]["Value"]) + metadata["spacing"] = np.asarray(metadata.get("00280030", {}).get("Value", (1.0, 1.0, 1.0))) dicom_data.append((data_array, metadata)) img_array: List[np.ndarray] = [] @@ -608,11 +603,13 @@ def _get_affine(self, metadata: Dict, lps_to_ras: bool = True): """ affine: np.ndarray = np.eye(4) + if not ("00200037" in metadata and "00200032" in metadata): + return affine # "00200037" is the tag of `ImageOrientationPatient` rx, ry, rz, cx, cy, cz = metadata["00200037"]["Value"] # "00200032" is the tag of `ImagePositionPatient` sx, sy, sz = metadata["00200032"]["Value"] - dr, dc = metadata["spacing"][:2] + dr, dc = metadata.get("spacing", (1.0, 1.0))[:2] affine[0, 0] = cx * dr affine[0, 1] = rx * dc affine[0, 3] = sx @@ -625,7 +622,7 @@ def _get_affine(self, metadata: Dict, lps_to_ras: bool = True): affine[2, 3] = sz # 3d - if "lastImagePositionPatient" in metadata.keys(): + if "lastImagePositionPatient" in metadata: t1n, t2n, t3n = metadata["lastImagePositionPatient"] n = metadata["spatial_shape"][-1] k1, k2, k3 = (t1n - sx) / (n - 1), (t2n - sy) / (n - 1), (t3n - sz) / (n - 1) @@ -640,7 +637,7 @@ def _get_affine(self, metadata: Dict, lps_to_ras: bool = True): def _get_array_data(self, img): """ Get the array data of the image. If `RescaleSlope` and `RescaleIntercept` are available, the raw array data - will be rescaled. The output data has the type: np.float32 + will be rescaled. The output data has the dtype np.float32 if the rescaling is applied. Args: img: a Pydicom dataset object. @@ -649,9 +646,9 @@ def _get_array_data(self, img): # process Dicom series if not hasattr(img, "pixel_array"): raise ValueError(f"dicom data: {img.filename} does not have pixel_array.") - data = img.pixel_array.astype(np.float32) + data = img.pixel_array - slope, offset = 1, 0 + slope, offset = 1.0, 0.0 rescale_flag = False if hasattr(img, "RescaleSlope"): slope = img.RescaleSlope @@ -659,8 +656,8 @@ def _get_array_data(self, img): if hasattr(img, "RescaleIntercept"): offset = img.RescaleIntercept rescale_flag = True - if rescale_flag is True: - data = data * slope + offset + if rescale_flag: + data = data.astype(np.float32) * slope + offset return data From 16ea235e3d32e481aacecf3a616879c62dff46ff Mon Sep 17 00:00:00 2001 From: Wenqi Li Date: Sun, 26 Jun 2022 18:55:33 +0100 Subject: [PATCH 8/8] adds consistency Signed-off-by: Wenqi Li --- monai/data/image_reader.py | 31 ++++++++++++++++++++++++++----- tests/test_load_image.py | 32 +++----------------------------- 2 files changed, 29 insertions(+), 34 deletions(-) diff --git a/monai/data/image_reader.py b/monai/data/image_reader.py index d6dd0188bc..bccd689a37 100644 --- a/monai/data/image_reader.py +++ b/monai/data/image_reader.py @@ -21,7 +21,12 @@ from torch.utils.data._utils.collate import np_str_obj_array_pattern from monai.config import DtypeLike, KeysCollection, PathLike -from monai.data.utils import correct_nifti_header_if_necessary, is_supported_format, orientation_ras_lps +from monai.data.utils import ( + affine_to_spacing, + correct_nifti_header_if_necessary, + is_supported_format, + orientation_ras_lps, +) from monai.transforms.utility.array import EnsureChannelFirst from monai.utils import ensure_tuple, ensure_tuple_rep, optional_import, require_pkg @@ -385,6 +390,8 @@ class PydicomReader(ImageReader): affine_lps_to_ras: whether to convert the affine matrix from "LPS" to "RAS". Defaults to ``True``. Set to ``True`` to be consistent with ``NibabelReader``, otherwise the affine matrix remains in the Dicom convention. + swap_ij: whether to swap the first two spatial axes. Default to ``True``, so that the outputs + are consistent with the other readers. kwargs: additional args for `pydicom.dcmread` API. more details about available args: https://pydicom.github.io/pydicom/stable/reference/generated/pydicom.filereader.dcmread.html#pydicom.filereader.dcmread If the `get_data` function will be called @@ -392,13 +399,20 @@ class PydicomReader(ImageReader): `stop_before_pixels` is `True`, and `specific_tags` covers all necessary tags, such as `PixelSpacing`, `ImagePositionPatient`, `ImageOrientationPatient` and all `pixel_array` related tags. + Note:: + + the current + """ - def __init__(self, channel_dim: Optional[int] = None, affine_lps_to_ras: bool = True, **kwargs): + def __init__( + self, channel_dim: Optional[int] = None, affine_lps_to_ras: bool = True, swap_ij: bool = True, **kwargs + ): super().__init__() self.kwargs = kwargs self.channel_dim = channel_dim self.affine_lps_to_ras = affine_lps_to_ras + self.swap_ij = swap_ij def verify_suffix(self, filename: Union[Sequence[PathLike], PathLike]) -> bool: """ @@ -557,15 +571,22 @@ def get_data(self, data): compatible_meta: Dict = {} for (data_array, metadata) in ensure_tuple(dicom_data): - img_array.append(data_array) - metadata["original_affine"] = self._get_affine(metadata, self.affine_lps_to_ras) - metadata["affine"] = metadata["original_affine"].copy() + img_array.append(np.ascontiguousarray(np.swapaxes(data_array, 0, 1) if self.swap_ij else data_array)) + affine = self._get_affine(metadata, self.affine_lps_to_ras) + if self.swap_ij: + affine = affine @ np.array([[0, 1, 0, 0], [1, 0, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]) + sp_size = list(metadata["spatial_shape"]) + sp_size[0], sp_size[1] = sp_size[1], sp_size[0] + metadata["spatial_shape"] = ensure_tuple(sp_size) + metadata["original_affine"] = affine + metadata["affine"] = affine.copy() if self.channel_dim is None: # default to "no_channel" or -1 metadata["original_channel_dim"] = ( "no_channel" if len(data_array.shape) == len(metadata["spatial_shape"]) else -1 ) else: metadata["original_channel_dim"] = self.channel_dim + metadata["spacing"] = affine_to_spacing(metadata["original_affine"], r=len(metadata["spatial_shape"])) _copy_compatible_dict(metadata, compatible_meta) return _stack_images(img_array, compatible_meta), compatible_meta diff --git a/tests/test_load_image.py b/tests/test_load_image.py index ea1adba40c..7b5572f5fc 100644 --- a/tests/test_load_image.py +++ b/tests/test_load_image.py @@ -199,7 +199,7 @@ def test_itk_reader(self, input_param, filenames, expected_shape): np.testing.assert_allclose(header["original_affine"], np_diag) self.assertTupleEqual(result.shape, expected_shape) - @parameterized.expand([TEST_CASE_10, TEST_CASE_11, TEST_CASE_12]) + @parameterized.expand([TEST_CASE_10, TEST_CASE_11, TEST_CASE_12, TEST_CASE_19, TEST_CASE_20, TEST_CASE_21]) def test_itk_dicom_series_reader(self, input_param, filenames, expected_shape, expected_np_shape): result, header = LoadImage(**input_param)(filenames) self.assertTrue("affine" in header) @@ -233,25 +233,6 @@ def test_itk_reader_multichannel(self): np.testing.assert_allclose(result[:, :, 1], test_image[:, :, 1]) np.testing.assert_allclose(result[:, :, 2], test_image[:, :, 2]) - @parameterized.expand([TEST_CASE_19, TEST_CASE_20, TEST_CASE_21]) - def test_pydicom_dicom_series_reader(self, input_param, filenames, expected_shape, expected_np_shape): - result, header = LoadImage(**input_param)(filenames) - self.assertTrue("affine" in header) - self.assertEqual(header["filename_or_obj"], f"{Path(filenames)}") - np.testing.assert_allclose( - header["affine"], - np.array( - [ - [0.0, -0.488281, 0.0, 125.0], - [-0.488281, 0.0, 0.0, 128.100006], - [0.0, 0.0, 68.33333333, -99.480003], - [0.0, 0.0, 0.0, 1.0], - ] - ), - ) - self.assertTupleEqual(tuple(header["spatial_shape"]), expected_shape) - self.assertTupleEqual(result.shape, expected_np_shape) - @parameterized.expand([TEST_CASE_22]) def test_dicom_reader_consistency(self, filenames): itk_param = {"reader": "ITKReader"} @@ -261,15 +242,8 @@ def test_dicom_reader_consistency(self, filenames): pydicom_param["affine_lps_to_ras"] = affine_flag itk_result, itk_header = LoadImage(**itk_param)(filenames) pydicom_result, pydicom_header = LoadImage(**pydicom_param)(filenames) - - # the row and column order is different - # see also: - # https://nipy.org/nibabel/dicom/dicom_orientation.html#i-j-columns-rows-in-dicom - - self.assertTrue(np.all(pydicom_result == np.transpose(itk_result, [1, 0, 2]))) - - r_c_transform = np.array([[0, 1, 0, 0], [1, 0, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]) - np.testing.assert_allclose(itk_header["affine"], pydicom_header["affine"] @ r_c_transform) + np.testing.assert_allclose(pydicom_result, itk_result) + np.testing.assert_allclose(itk_header["affine"], pydicom_header["affine"]) def test_load_nifti_multichannel(self): test_image = np.random.randint(0, 256, size=(31, 64, 16, 2)).astype(np.float32)