ImagingDataCommons / highdicom

High-level DICOM abstractions for the Python programming language
https://highdicom.readthedocs.io
MIT License
168 stars 35 forks source link

DICOM SEG for CMR cine data (2D + t) #174

Open apint0-media opened 2 years ago

apint0-media commented 2 years ago

Trying to generate a dicom seg file based on a multiple unsigned np.array, which is of shape (#_frames, #rows, #cols). The way I've approached the task was by generating dicom seg with the same series uid across all frames of the segmentation array. However, when displaying the output the 3 different delineated structures do not appear on a single frame (spatially they do not overlap).

Bellow is an example of the code that I've used to generate the dicom outputs.

Is it possible to show all the 3 structures in a single frame?

# img_dsets (list): list of pydicom.Datasets with lenght equal to the number of short-axis cine data
# seg_vol (np.array): Segmentation volume array with shape (#_frames, #rows, #cols)
import pydicom as pyd
import highdicom as hd

manufacturer_name = "dummy_manufacturer_name"
manufacturer_model_name = "dummy_model_name"
software_version = "0.1"

seg_labels = {
    "myocardium": {"class_number": 1, "sct_code": pyd.sr.codedict.codes.SCT.LeftVentricleMyocardium},
    "left_ventricle": {"class_number": 3, "sct_code": pyd.sr.codedict.codes.SCT.LeftVentricle},
    "right_ventricle": {"class_number": 2, "sct_code": pyd.sr.codedict.codes.SCT.RightVentricle},
}

algo_details = hd.AlgorithmIdentificationSequence(
    name=manufacturer_model_name,
    version=software_version,
    family=pyd.sr.codedict.codes.cid7162.ArtificialIntelligence,
    source=manufacturer_name,
)

segment_descriptions = []
for label_name, label_info in seg_labels.items():

    seg_details = hd.seg.SegmentDescription(
        segment_number=label_info["class_number"],
        segment_label=label_name,
        segmented_property_category=pyd.sr.codedict.codes.cid7150.AnatomicalStructure,
        segmented_property_type=label_info["sct_code"],
        algorithm_type=hd.seg.SegmentAlgorithmTypeValues.AUTOMATIC,
        algorithm_identification=algo_details,
        tracking_uid=hd.UID(),
        tracking_id=f"Cardiac structure #{class_number}",
    )

    segment_descriptions.append(seg_details)

series_uid = hd.UID()

for frame_idx in range(seg_vol.shape[0]):

    seg_dset = hd.seg.Segmentation(
        source_images=[img_dsets[frame_idx]],
        pixel_array=seg_vol[frame_idx],
        segmentation_type=hd.seg.enum.SegmentationTypeValues.BINARY,
        segment_descriptions=segment_descriptions,
        series_description="Segmentation-Test",
        series_number=5,
        series_instance_uid=series_uid,
        sop_instance_uid=hd.UID(),
        instance_number=int(frame_idx) + 1,
        manufacturer=manufacturer,
        manufacturer_model_name=manufacturer_model_name,
        software_versions=software_version,
        device_serial_number=str(img_dsets[frame_idx].DeviceSerialNumber),
    )

    seg_dsets.append(seg_dset)

Also, I've tried to generate a single dataset file for all frames at once, and with that approach the output data only shows a total number of frames equal to the # of classes, not iterating over all the frames.

hackermd commented 2 years ago

Is it possible to show all the 3 structures in a single frame?

I am not sure I understand your question correctly. By "structures" do you mean "segments"?

The highdicom.seg.Segmentation class provides methods for accessing the pixel data, e.g., get_pixels_by_dimension_index_values(). These methods provide a combine_segments parameter. When the argument is True, the methods will combine multiple segments into a single pixel array to form a label mask.

CPBridge commented 2 years ago

Hi @apint0-media , thanks for trying out highdicom.

First of all, please note that the correct way to create the DICOM SEG in a situation like this where you have a series of input files that are segmented is to place the full segmentation in a single DICOM SEG file. This is because DICOM SEG is a newer part of the standard and expects things to be multiframe.

Secondly, also note that, due to the way the standard was written, it is not a simple exercise to figure out which frames in the stored DICOM SEG correspond to which frames in your source images. This is because a single frame of the SEG file is stored for each frame of the source images and each segment ("structure"), but empty frames are not stored. If you are assuming that the .pixel_array of the resulting seg files correspond in a simple way to your input frames, this may make things appear like they are not working correctly when in fact they are. As @hackermd mentioned, we recommend that you use the methods we built into highdicom to reconstruct segmentation masks from the stored file, as we have figured out all this logic for you.

Could you please provide the following information and we'll see if we can figure out what is going on in your case:

Thanks!

apint0-media commented 2 years ago

I am not sure I understand your question correctly. By "structures" do you mean "segments"?

The highdicom.seg.Segmentation class provides methods for accessing the pixel data, e.g., get_pixels_by_dimension_index_values(). These methods provide a combine_segments parameter. When the argument is True, the methods will combine multiple segments into a single pixel array to form a label mask.

Thanks for the help! Think an image might be better to describe what was I referring to. The goal was to obtain something like figure bellow, for each frame (viewer: 3D Slicer)

image-20220411-085648

Currently, with the above code, the output does not show the 3 regions in the same plane and does not contain all the information of the patient along the frames dimension.

I will explore the get_pixels_by_dimension_index_values(), to verify if the label map is generated accordingly.

apint0-media commented 2 years ago

Hi @apint0-media , thanks for trying out highdicom.

First of all, please note that the correct way to create the DICOM SEG in a situation like this where you have a series of input files that are segmented is to place the full segmentation in a single DICOM SEG file. This is because DICOM SEG is a newer part of the standard and expects things to be multiframe.

Secondly, also note that, due to the way the standard was written, it is not a simple exercise to figure out which frames in the stored DICOM SEG correspond to which frames in your source images. This is because a single frame of the SEG file is stored for each frame of the source images and each segment ("structure"), but empty frames are not stored. If you are assuming that the .pixel_array of the resulting seg files correspond in a simple way to your input frames, this may make things appear like they are not working correctly when in fact they are. As @hackermd mentioned, we recommend that you use the methods we built into highdicom to reconstruct segmentation masks from the stored file, as we have figured out all this logic for you.

Could you please provide the following information and we'll see if we can figure out what is going on in your case:

  • What values are in seg_vol (is each pixel either a 0, 1, 2, 3)?
  • How are you displaying the resulting segmentation? Are you using an existing viewer, if so what is it? Or have you written your own python code to display it? If you could you post that code? There are very few pieces of software that will correctly display a DICOM seg file at the moment. One of them is the OHIF viewer, you might want to look into that.
  • Are you able to provide a screenshot of the display so we can better understand the issue?

Thanks!

Hi @CPBridge. Great feedback, thanks!

I will definitely check what is happening to the .pixel_array and its dimensions. Regarding the questions asked:

  1. Yes, each pixel is either one of four values belonging to the set [0, 1, 2, 3].
  2. I am using 3D Slicer, but also MicroDicom. Curiously, I was checking OHIF to discard also if it was on the side of the viewer.
  3. My previous comment shows one of the views of the "desired output". As soon as possible I will provide a screenshot of the results
CPBridge commented 2 years ago

Thanks @apint0-media. I am wondering whether this may have something to do with the time dimension. Unfortunately the way that segmentation frames are indexed in segmentations is extremely flexible and we have only implemented a couple of specific situations in highdicom so far, they may not be dealing with the time dimension correctly. We have not tried with data with a time dimension before. Or the viewers may not be treating the time direction correctly. Or both. It's going to be tricky to tease apart which. Any data that you are able to share is going to help a lot. To start, would you be able to post the contents of the SharedFunctionalGroupsSequence, the PerFrameFunctionalGroupsSequence, and the DimensionIndexSequence of the segmentation instance(s) that you are creating (e.g. as printed by pydicom or dcmdump, or similar)?

the output [...] does not contain all the information of the patient along the frames dimension

Could you clarify what you mean here?

apint0-media commented 2 years ago

the output [...] does not contain all the information of the patient along the frames dimension

Could you clarify what you mean here?

What I'm viewing with MicroDicom is only the information of a single input reference, i.e. 3 object of seg_arraythat concern the first frame. Since I recon that my vocabulary might not relate precisely to what I'm trying to detail, I've attached a small video sample of the current output.

https://user-images.githubusercontent.com/104420558/167875665-efd5c3e6-760b-494a-b9a8-d3d59d16bf1e.mp4

Made some changes to the code based on our interactions (the video sample already reflects them).

With that in mind, here is the output of the contents SharedFunctionalGroupsSequence, PerFrameFunctionalGroupsSequence, and DimensionIndexSequence.

(5200, 9229)  Shared Functional Groups Sequence  1 item(s) ----
   (0020, 9116)  Plane Orientation Sequence  1 item(s) ----
      (0020, 0037) Image Orientation (Patient)         DS: [0.59040129158056, 0.7679697850161, -0.2482916112202, -0.3876612065461, -1.0605152e-008, -0.921801925003]
      ---------
   (0028, 9110)  Pixel Measures Sequence  1 item(s) ----
      (0018, 0050) Slice Thickness                     DS: '8.0'
      (0018, 0088) Spacing Between Slices              DS: '10.0'
      (0028, 0030) Pixel Spacing                       DS: [1.7708333730698, 1.7708333730698]
      ---------
   ---------
(5200, 9230)  Per-frame Functional Groups Sequence  1 item(s) ----
   (0008, 9124)  Derivation Image Sequence  1 item(s) ----
      (0008, 2112)  Source Image Sequence  1 item(s) ----
         (0008, 1150) Referenced SOP Class UID            UI: MR Image Storage
         (0008, 1155) Referenced SOP Instance UID         UI: 1.3.12.2.1107.5.2.30.26258.2018042708245526542901674
         (0028, 135a) Spatial Locations Preserved         CS: 'YES'
         (0040, a170)  Purpose of Reference Code Sequence  1 item(s) ----
            (0008, 0100) Code Value                          SH: '121322'
            (0008, 0102) Coding Scheme Designator            SH: 'DCM'
            (0008, 0104) Code Meaning                        LO: 'Source image for image processing operation'
            ---------
         ---------
      (0008, 9215)  Derivation Code Sequence  1 item(s) ----
         (0008, 0100) Code Value                          SH: '113076'
         (0008, 0102) Coding Scheme Designator            SH: 'DCM'
         (0008, 0104) Code Meaning                        LO: 'Segmentation'
         ---------
      ---------
   (0020, 9111)  Frame Content Sequence  1 item(s) ----
      (0020, 9157) Dimension Index Values              UL: [2, 1]
      ---------
   (0020, 9113)  Plane Position Sequence  1 item(s) ----
      (0020, 0032) Image Position (Patient)            DS: [15.063542066509, -136.90274478107, 165.82297569271]
      ---------
   (0062, 000a)  Segment Identification Sequence  1 item(s) ----
      (0062, 000b) Referenced Segment Number           US: 2
      ---------
   ---------

DimensionIndexSequence

(0020, 9221)  Dimension Organization Sequence  1 item(s) ----
   (0020, 9164) Dimension Organization UID          UI: 1.2.826.0.1.3680043.10.511.3.13262044635184244845644933344388480
   ---------
(0020, 9222)  Dimension Index Sequence  2 item(s) ----
   (0020, 9164) Dimension Organization UID          UI: 1.2.826.0.1.3680043.10.511.3.13262044635184244845644933344388480
   (0020, 9165) Dimension Index Pointer             AT: (0062, 000b)
   (0020, 9167) Functional Group Pointer            AT: (0062, 000a)
   (0020, 9421) Dimension Description Label         LO: 'Segment Number'
   ---------
   (0020, 9164) Dimension Organization UID          UI: 1.2.826.0.1.3680

043.10.511.3.13262044635184244845644933344388480
   (0020, 9165) Dimension Index Pointer             AT: (0020, 0032)
   (0020, 9167) Functional Group Pointer            AT: (0020, 9113)
   (0020, 9421) Dimension Description Label         LO: 'Image Position Patient'
   ---------
apint0-media commented 2 years ago

Also, I've tried to use get_pixels_by_dimension_index_values but I'm having some troubles with the argument dimension_index_values. Dont' know how to define it properly.