innolitics / dicom-numpy Goto Github PK

_sort_by_slice_position in combine_slices.py would be useful to dicom_numpy callers so that instance metadata (e.g. SliceThickness) can be associated with the volume generated from combine_slices.

from dicom_numpy import *

Produces the following error:
TypeError: Item in dicom_numpy.__all__ must be str, not function

Some diffusion MRI series combine several scans into a single series, with no acquisition number to distinguish them. This causes issues for dicom-numpy's generated image transformation.

Because slices are sorted by slice position, a diffusion series can appear to have multiple slices at the same location in space. This means the median space between slices will be 0, causing the image transformation to have all zeros in the slice axis.

To resolve this, slices should first be attempted to be sorted by instance number. If that fails, we can fall back to slice position sorting.

Hi, after combining slices, how could I check the voxel spacing between each slice? Is that same as the pixel spacing in the 2D slice? Thanks!

combine_slices requires HighBit and BitsStored (HighBit + 1 as required by standard) to be the same across slices being combined. However, the standard does not require this, only for the values to be uniform within the same slice.

While some CT scans are generated with a uniform HighBit (either for the maximum pixel value in the series, or just 15 when BitsAllocated is 16), others are generated with HighBit using the slice-specific maximum value. Scans 1, 5, and 8 at coronacases.org are generated this way. All these cases validate with the following:

import pydicom
import sys
import glob
import numpy as np

# https://wiki.python.org/moin/BitManipulation
def bitLen(int_type):
    length = 0
    while (int_type):
        int_type >>= 1
        length += 1
    return(length)

for dcm_file in glob.glob(f'{sys.argv[1]}/**/*.dcm', recursive=True):
    ds = pydicom.dcmread(dcm_file)
    max_pixel = np.max(ds.pixel_array)
    assert ds.HighBit + 1 == ds.BitsStored

    bit_length = bitLen(max_pixel)

    info = f'{dcm_file} {bit_length} {ds.HighBit} {max_pixel}'
    assert bit_length == ds.HighBit or ds.BitsStored == ds.BitsAllocated, info

Hi there, thanks for the very nice software. I have a question. In trying to convert a single dcm to numpy array
I get an output with 2D while it is a 3D image --can you help ?
thanks!
david
IM-0002-0005.dcm.zip
']

import glob
import pydicom
import dicom_numpy
dicomlist = ['/Users/dsb/Downloads/Dicoms/IM-0002-0005
def extract_voxel_data(list_of_dicom_files):
datasets = [pydicom.read_file(f) for f in list_of_dicom_files]
try:
voxel_ndarray, ijk_to_xyz = dicom_numpy.combine_slices(datasets)
except dicom_numpy.DicomImportException as e:
# invalid DICOM data
raise
return voxel_ndarray

outp = extract_voxel_data(dicomlist)

Our documentation linked in the README is out of date. We should:

1. Update the documentation published on read the docs, and
2. Automate the publishing of the documentation as a GitHub action that runs when we tag and push a release

Can you access the affine transform of a dicom file with the dicom-numpy library?

https://pydicom.github.io/pydicom/stable/transition_to_pydicom1.html

_sort_by_slice_position in combine_slices.py would be useful to dicom_numpy callers so that instance metadata (e.g. SliceThickness) can be associated with the volume generated from combine_slices.

Hello authors,

Thanks for the useful toolkit. I just wonder if the library can interpolate between slices so that the resulting numpy array would be even grid in distance unit along all three axes?

thanks

This is a great work, actually, I am looking for DICOM image stitching module that can stitch two or three images into one 2D image like a panorama. Is there any module or library to do so?

abad26d
The changes introduced in this commit were justified by the fact that in DICOM, the first slice is considered the top of the patient, whereas the reverse is true: the Z axis increases from the feet towards the head, see e.g. https://dicom.innolitics.com/ciods/ct-image/image-plane/00200032
I'm trying to combine DICOM slices and convert the resulting stack into a NiFTI file, but when I render it in 3D Slicer, the volume is flipped along the Z axis.

Hi,

combine_slices() works very well for me for multiple 2D DICOM images, which is great.

However, it's not working correctly for Enhanced MR Images, which have a hierarchical structure. If I call subcategories in DICOM, it won't find most variables. The same happens if I call the function for the whole DICOM file (which is inside a list).

AttributeError: 'FileDataset' object has no attribute 'ImageOrientationPatient'

Other times when I run the same function, it complains about the SliceLocation.

I'm aware that this involves plenty of work, so I would like to suggest the Enhanced MR as improvement.

Kind regards,

João Sousa

After inspecting in detail the combine_slices code, I found some optimization possibilities:

• _sort_by_slice_position is called twice:
  First in the voxels (_merge_slice_pixel_arrays),
  Second in the affine (_ijk_to_patient_xyz_transform_matrix).
  Obviously, it can be done once in the outer combine_slices function.
• _slice_positions is called 4 times:
  _validate_slices_form_uniform_grid - 1 time
  _slice_spacing - 1 time
  _sort_by_slice_position - 1 time, but _sort_by_slice_position is called twice --> 2 times
• _extract_cosines is called 6 times:
  _validate_image_orientation - 1 time
  _ijk_to_patient_xyz_transform_matrix - 1 time
  _slice_positions - 1 time, but _slice_positions is called 4 times --> 4 times

All these calls better be done once (and can be).
However,

1. this may incorporate substantial redesign of the code,
2. it will change the functionality of _merge_slice_pixel_arrays and _ijk_to_patient_xyz_transform_matrix, which should not be used alone, but there are some special tests for _merge_slice_pixel_arrays.

Originally posted by @jond01 in #16 (comment)

Hi, I have recently reviewed the contrib-pydicom code. The script input-output/pydicom_series.py there has (more or less) the same purpose as dicom-numpy.
I noticed there two points that we may want to include also here:

1. Read a folder: add an API to read a folder and extract the image array and affine.
   A similar idea is in the example from the docs, which receives a list of files:

   import pydicom
   import dicom_numpy
   
   def extract_voxel_data(list_of_dicom_files):
       datasets = [pydicom.dcmread(f) for f in list_of_dicom_files]
       try:
           voxel_ndarray, ijk_to_xyz = dicom_numpy.combine_slices(datasets)
       except dicom_numpy.DicomImportException as e:
           # invalid DICOM data
           raise
       return voxel_ndarray

   We may go one level up, and receive only the path of the folder containing these files.
   The files within the folder can be filtered to only dicoms with pydicom's built-in is_dicom function, and further split into distinct series (according to the SeriesInstanceUID). combine_slices will be called for each series, and the returned data will be a list of [(voxels0, affine0), (voxels1, affine1), ...].

2. Tighten the dtype of rescaled images: currently, dicom-numpy uses np.float32 every time there is a RescaleSlope or RescaleIntercept:
   

   dicom-numpy/dicom_numpy/combine_slices.py

   Line 108 in 204e955

   voxels_dtype = np.float32 if rescale else slice_dtype

   
   However, many times it is not necessary - for example, if the rescale slope is 1.0 and the rescale intercept is -1024.0, the image can still be of an integer dtype.
   contrib-pydicom seems to have some clever way to determine the proper dtype.

• most of the code can be taken from the autoreg project
• should handle the DICOM timezone tag, when present

See: https://groups.google.com/forum/#!topic/comp.protocols.dicom/rCOq1TUcQ2Q

We should probably reference this, and other relevant links, in our documentation.

Hi,
Combining a single slice dataset (extracting its image and affine) -

from pydicom import dcmread
from dicom_numpy import combine_slices

ds0 = dcmread('dcm_file.dcm')
voxels, transform = combine_slices([ds0])

- produces an error in _check_for_missing_slices, L221:

    if not np.allclose(slice_positions_diffs, slice_positions_diffs[0], atol=0, rtol=1e-5):

slice_positions_diffs is empty because there is 1 slice -> no diff.
I suggest adding the following check - as in _slice_spacing, L231:

    if len(slice_datasets) > 1:

Hi there,
For PET images, it's typical for slices within the same time point to have different values for RescaleSlope and RescaleIntercept because of large differences of uptake from slice to slice. So it might make sense to just rescale the pixel data based on these values and then incorporate everything into the same 3D matrix rather than just giving an exception. This should be reasonably straight forward so if you don't have time I can just fork and file a pull request?

Thanks!