`seg.get_pixels_by_source_instance` -- support search by SOP_INSTANCE_IDS in the case of multiple sources.

[nicolay-r]

Dear highdicom team,

I have the question which relates to searching for particular pixel data (layer) in DICOMSEG file that has Multiple Souces.
In the attached example those sources are: liver, tumor.

...
(0062,0002)  Segment Sequence  2 item(s) ---- 
   (0062,0004) Segment Number                      US: 1
   (0062,0005) Segment Label                       LO: 'Liver'
   ---------
   (0062,0004) Segment Number                      US: 2
   (0062,0005) Segment Label                       LO: 'Tumor-1'
...

Using the following snippet:

seg = hd.seg.segread(buffer)
mask = seg.get_pixels_by_source_instance(
    source_sop_instance_uids=["1.2.826.0.1.3680043.8.274.1.1.2081054905.266506.1768410047.96982"],
    ignore_spatial_locations=True,
    assert_missing_frames_are_empty=True
)

Would result in:

RuntimeError: Indexing via source frames is not permissible since some frames in the image specify multiple source frames.

It works flawlessly in the case when the source is single (seg._single_source_frame_per_frame == 1); but in this case we have two sources.

Question: is there a way to pefrom a search by sop_instance_uids with known source?
(like set liver or tumor in that example)

seg_multisource_example.dicom.zip

[CPBridge]

Hi @nicolay-r

The get_pixels_by_source_instance method matches the requested UIDs against the information contained in the PerFrameFunctionalGroupsSequence in the file. There is one item of the PerFrameFunctionalGroupsSequence per frame in the segmentation image, and optionally this item will contain a DerivationImageSequence listing the UIDs of the image(s) from which that segmentation frame is derived.

The issue with the file you sent is that the DerivationImageSequence is not present, so there is no way to match segmentation frames to source images given UIDs alone, this information is simply not present in the Segmentation file:

>>> seg.PerFrameFunctionalGroupsSequence[0]
(0020,9111)  Frame Content Sequence  1 item(s) ----
   (0020,9157) Dimension Index Values              UL: [1, 1]
   ---------
(0020,9113)  Plane Position Sequence  1 item(s) ----
   (0020,0032) Image Position (Patient)            DS: [187.007034, 109.286133, -52.6307831]
   ---------
(0062,000A)  Segment Identification Sequence  1 item(s) ----
   (0062,000B) Referenced Segment Number           US: 1
   ---------

The error message you sent is a little unfortunate as it misrepresents the problem. It is intended to convey that the requested process is only possible if each segmentation frame specifies exactly one (i.e. not zero or 2+) source images. In the latest release of the library (0.26.0+) you should find that the error message is clearer (and furthermore the constraint is relaxed such that it can still work if 2+ images are referenced by a given segmentation frame):

RuntimeError: Indexing via source frames is not possible because frames do not contain information about source frames.

So in short, what you are asking for is impossible I'm afraid (the information is just not in the segmentation file). You will have to find some other way to match segmentation frames to images, or if you control the creation process, ensure that the necessary information is included in the segmentation. Most likely you will have to match based on spatial position. Unfortunately this is not something highdicom can currently handle, though it is something I plan to add

[CPBridge]

Edit: there is actually a way to do what you want as long as the segmentation exists on the same grid as the source images (i.e. same pixel spacing, same orientation), but it is perhaps not the most elegant at the moment:

# Read the source image
source_image = hd.imread("path/to/source_image.dcm")

# Get the volume geometry from the source image
source_geometry = source_image.get_volume_gemoetry()

# Read the segmentation
seg = hd.seg.segread("/path/to/segmentation.dcm")

# Get a volume from the segmentation
seg_vol = seg.get_volume()

# Use the match_geometry method to crop the segmentation to the geometry of the source image
cropped_seg_vol = seg_vol.match_geometry(source_geometry)

# Now you have a volume containing segmentation pixel data for the same volume that the source image occupies

This is perhaps not ideal because the entire volume of pixel data has to be created from the segmentation and then cropped to a single slice

[CPBridge]

Re-reading your question, maybe you are asking more about how to choose a particular segment (e.g. liver or tumor). In which case yes, you just specify the segment_numbers parameter on any of the get_pixels... methods or the get_volume method. See this and this section in the docs.

[nicolay-r]

@CPBridge , first of all thank you for such a verbose response, including workargound suggestions.

Initially experiment with highdicom==0.25.0.

The issue with the file you sent is that the DerivationImageSequence is not present,

Got it, but coming for the first time to DerivationImageSequence.

The error message you sent is a little unfortunate as it misrepresents the problem.
Agree, there must be something missed on my side in relation to metadata and segments matching during construction.

The attached segmentation file represent a baked DICOMSEG from source image (NIFTI) and numerous annotations (NIFTI). Those were first converted with pyplastimatch and then tied together into DICOMSEG. As I see from your comments there might be various sources of misalignments that causes challenged we see in this thread

Edit: there is actually a way to do what you want as long as the segmentation exists on the same grid as the source images (i.e. same pixel spacing, same orientation), but it is perhaps not the most elegant at the moment:
Re-reading your question, maybe you are asking more about how to choose a particular segment (e.g. liver or tumor)
Since I have
Got an Idea, so that is how I can have the related volume:

seg = hd.seg.segread("data/seg_multisource_example.dicom")
seg_vol = seg.get_volume(segment_numbers=[1])

However, the similar query:

mask = seg.get_pixels_by_source_instance(
    source_sop_instance_uids=["1.2.826.0.1.3680043.8.274.1.1.2081054905.266506.1768410047.96982"],
    segment_numbers=[1],   #  <- pikcing the segment number 1.
    ignore_spatial_locations=True,
    assert_missing_frames_are_empty=True
)

Would fall to the same exception.

Anyway, I believe that eitherway the problem with absence DerivationImageSequence, i.e. it is not supposed to work as several internal information for matching is not there.

Might be easier to follow your suggestion on finding another way for baking DICOMSEG with only a single source

[CPBridge]

However, the similar query, [...] Would fall to the same exception.

Yes absolutely, choosing a specific segment does not solve the issue with the lack of DerivationImageSequence, these are essentially unrelated concerns (your original post seemed to conflate the two). Your snippet would work if the file had DerivationImageSequence information regardless of whether or not you choose to specify a specific segment (with the segment_numbers parameter) or request all segments (by omitting that parameter as you were originally doing). In the latter case it would return an array with the multiple segments stacked down the last array dimension.

Might be easier to follow your suggestion on finding another way for baking DICOMSEG with only a single source

I'm not sure what exactly you mean by source here. Perhaps you mean a single segment (aka class), e.g. just liver or just tumor but not both in the same file. If so, again, the number of segments is a totally separate consideration to the inclusion of the DerivationImageSequence. In other words I would not expect using a single segment per file to solve anything.

It may be that your segmentation mask(s) is (are) simply not pixel-for-pixel aligned with the original file, in which case you might need to transpose/resample the mask(s) before creating the DICOM seg. Or it might just mean that the tool you are using to create them (dcmqi it seems) just isn't adding the per-frame derivation information. @fedorov might be able to help with that

[fedorov]

The attached segmentation file represent a baked DICOMSEG from source image (NIFTI) and numerous annotations (NIFTI).
...
It may be that your segmentation mask(s) is (are) simply not pixel-for-pixel aligned with the original file, in which case you might need to transpose/resample the mask(s) before creating the DICOM seg. Or it might just mean that the tool you are using to create them (dcmqi it seems) just isn't adding the per-frame derivation information. @fedorov might be able to help with that

I completely agree with @CPBridge. Most likely reason you do not have the DerivationImageSequence populated is because the geometry of the image does not mach that of the segmentation (at least according to the heuristic implemented in dcmqi to find matching geometry slices). dcmqi does attempt to match segmentation frames to the segmented image slices and will populate that sequence if matching was successful.

I recommend that you take your NIfTIs and write an ITK script that reorients both the images and the segmentations into the same orientation using for example OrientImageFilter, and also consider resampling segmentations to the geometry of the image (if there are slight differences in orientation vectors/pixel spacing). I would then save the result into NRRD, which does not mess with orientation like NIfTI and reconvert.