Mismatch between DICOM segmentation volume and original image

[mingxin-zheng]

Describe the bug
Could be related to #6928

Background:

This issue appears when I am working with DICOM segmentations generated by OHIF viewer.
There are dicom images and segmentations (labels) in OHIF.
I can read the dicom images succesfully with MONAI, but the labels in DICOM segmentation "SE00000.dcm" are organized by invidual slices. So if I read it directly with LoadImaged, the volume sizes of image and label will be mismatch.

In my example below, the s0004 (total segmentator) has a image size of [255, 177, 440], but the "label" generated by the DICOM segmentation is (7987, 177, 255)

I have a script to fix this:

import os
import pydicom
import numpy as np
from monai.transforms import LoadImage, SaveImage
from monai.data.meta_tensor import MetaTensor

print("Reading")

# Read the DICOM directory
dicom_dir = "s0004 ANONYMOUS/Unknown Study/CT"
seg_file_path = "s0004 ANONYMOUS/Unknown Study/SEG s0004_seg/SE000000.dcm"
output_dir = "."
# CT/CT_trans.nii.gz contains in the seg in `output_dir`

dicom_img = LoadImage()(dicom_dir)
seg_vol = np.zeros(dicom_img.shape)
# print(seg_vol.shape)

# Read the segmentation DICOM file
seg_dcm = pydicom.dcmread(seg_file_path)
# Extract the seg_record from the segmentation DICOM file
seg_record = seg_dcm.pixel_array  # assuming that seg_dcm.pixel_array gives the segmentation record

print("Shape of the dicom seg: ", seg_record.shape)
print("Shaping")
# # Iterate over the PerFrameFunctionalGroupsSequence and update the seg_vol array
for i, frame_item in enumerate(seg_dcm.PerFrameFunctionalGroupsSequence):
    frame_content_seq = frame_item.FrameContentSequence[0]
    class_id, slice_index = frame_content_seq.DimensionIndexValues
    
    # Assuming seg_record[i, y, x] gives the i-th segmentation record for seg_vol[x, y, slice_index - 1]
    seg_slice = seg_record[i, :, :].T
    seg_vol[:, :, slice_index - 1][seg_slice > 0] = class_id

print("Writing")

seg = MetaTensor(seg_vol, affine=dicom_img.affine, meta=dicom_img.meta)
SaveImage(output_dir=output_dir, output_dtype="uint16")(seg)

print(f"Shape of DICOM read from {dicom_dir}")
print(dicom_img.shape)
seg_output = LoadImage()(os.path.join(output_dir, "CT", "CT_trans.nii.gz"))
print(f"Shape of nii.gz written to {output_dir}")
print(seg_output.shape)

To Reproduce
Steps to reproduce the behavior:

1. Go to '...'
2. Install '....'
3. Run commands '....'

Below is the data and script to reproduce:
data_and_script.zip

Expected behavior
A clear and concise description of what you expected to happen.

Screenshots
If applicable, add screenshots to help explain your problem.

Environment

Ensuring you use the relevant python executable, please paste the output of:

python -c 'import monai; monai.config.print_debug_info()'

Additional context
Add any other context about the problem here.

[wyli]

Loading DICOM seg object is not supported.. @ericspod @vikashg would be great to understand the current best practice, any recommendations for offline preprocessing?

[wyli]

one possibility is to use segimage2itkimage from dcmqi https://qiicr.gitbook.io/dcmqi-guide/opening/cmd_tools/seg/segimage2itkimage converting the segmentation into nifti

then use monai transforms to convert the seg.nii into the same space as the corresponding image:

import monai.transforms as mt

img = mt.LoadImage(ensure_channel_first=True)('image.nii.gz')
seg = mt.LoadImage(ensure_channel_first=True)('seg-1.nii.gz')
out_seg = mt.ResampleToMatch(mode="nearest", padding_mode="zeros")(seg, img)
saved = mt.SaveImage("/Users/Documents/temp")(out_seg)

[mingxin-zheng]

Thanks @wyli for the investigation. This is close to what I have been doing now with the workaround script.

From what I have read from the segmentation object meta, it doesn't contain info of the original image. So it needs to info image (w, h, d) plus the affine to correctly align the label and the image.

[ericspod]

I'm afraid for now this conversion is probably the best idea. We have been discussing the need for DICOM seg support with MONAI Deploy but haven't implementing support in core. @vikashg perhaps would have a better idea of seg software support.

[KumoLiu]

An alternative way using pydicom_seg.

import itk
import pydicom
import pydicom_seg
import numpy as np
import SimpleITK as sitk
import monai.transforms as mt
from monai.data import ITKReader
from monai.data import MetaTensor

img_path = "/workspace/data_and_script/s0004 ANONYMOUS/Unknown Study/CT"
seg_path = "/workspace/data_and_script/s0004 ANONYMOUS/Unknown Study/SEG s0004_seg/SE000000.dcm"

seg_dcm = pydicom.dcmread(seg_path)
reader = pydicom_seg.MultiClassReader()
result = reader.read(seg_dcm)
seg_image = result.image

# sitk image -> itk image
image_dimension = 3
itk_image = itk.GetImageFromArray(sitk.GetArrayFromImage(seg_image))
itk_image.SetOrigin(seg_image.GetOrigin())
itk_image.SetSpacing(seg_image.GetSpacing())   
itk_image.SetDirection(itk.GetMatrixFromArray(np.reshape(np.array(seg_image.GetDirection()), [image_dimension]*2)))

itkreader = ITKReader()
_seg_data = itkreader.get_data(itk_image)
seg_data = MetaTensor(_seg_data[0], meta=_seg_data[1])
img = mt.LoadImage(ensure_channel_first=True)(img_path)
out = mt.ResampleToMatch(mode="nearest", padding_mode="zeros")(seg_data.unsqueeze(0), img)

[vikashg]

This one works.

[vikashg]

Hello @wyli and @ericspod I will post some code tonight and link it here. Sorry I was tied up with few things so I wasn't able to respond earlier

[vikashg]

Hello all,
The code proposed by @KumoLiu works. However, it is based on pydicom_seg which is fine unless we want to worry about one more dependency. I didn't know about pydicom_seg. So I wrote some example code for reading dicomseg. However, this code only supports a dicom seg containing a single label. I am working on adding some support for multi label data (as is the case of the example shared). Here is the example code

https://github.com/Project-MONAI/monai-dicom-utils/tree/vikash/dicomseg/app
Here is the sample image