qurit / rt-utils

A minimal Python library to facilitate the creation and manipulation of DICOM RTStructs.
MIT License
181 stars 56 forks source link

cv.fillPoly() exception when running rtstruct.get_roi_mask_by_name #87

Open tomaroberts opened 1 year ago

tomaroberts commented 1 year ago

Hi there,

For my project, I've got multiple patients and I've generated an individual RTStructBuilder object per patient. Within the RTStructBuilder objects I have multiple ROIs (>50).

I am trying to use the rtstruct.get_roi_mask_by_name function to generate binary masks of each region, however I am getting an exception within image_helper.py line 278.

The exception occurs with some – but not all – of the ROIs within the RTStructBuilder object.

For instance, it will work with region1, e.g.:

mask = rtstruct.get_roi_mask_by_name('region1')

whereas say region2 will fail with the following error message:

Traceback (most recent call last):
  File "/Applications/PyCharm.app/Contents/plugins/python/helpers/pydev/pydevd.py", line 1483, in _exec
    pydev_imports.execfile(file, globals, locals)  # execute the script
  File "/Applications/PyCharm.app/Contents/plugins/python/helpers/pydev/_pydev_imps/_pydev_execfile.py", line 18, in execfile
    exec(compile(contents+"\n", file, 'exec'), glob, loc)
  File "/Users/tr17/Library/Application Support/JetBrains/PyCharm2021.2/scratches/load_rt.py", line 123, in <module>
    mask = rtstruct_TS.get_roi_mask_by_name('lung_upper_lobe_right')
  File "/Users/tr17/venv/prague/lib/python3.9/site-packages/rt_utils/rtstruct.py", line 113, in get_roi_mask_by_name
    return image_helper.create_series_mask_from_contour_sequence(
  File "/Users/tr17/venv/prague/lib/python3.9/site-packages/rt_utils/image_helper.py", line 247, in create_series_mask_from_contour_sequence
    mask[:, :, i] = get_slice_mask_from_slice_contour_data(
  File "/Users/tr17/venv/prague/lib/python3.9/site-packages/rt_utils/image_helper.py", line 281, in get_slice_mask_from_slice_contour_data
    cv.fillPoly(img=slice_mask, pts = polygons, color = 1)
cv2.error: OpenCV(4.7.0) /Users/xperience/GHA-OCV-Python/_work/opencv-python/opencv-python/opencv/modules/imgproc/src/drawing.cpp:2402: error: (-215:Assertion failed) p.checkVector(2, CV_32S) >= 0 in function 'fillPoly'

I've looked up the fillPoly error online and there are various posts about casting as int32, but I think you are already doing that properly.

I'm a bit stuck... much appreciated if you can help. Thanks in advance.

tomaroberts commented 1 year ago

To add to this, I've just noticed in Debug mode that the polygons variable within the function get_slice_mask_from_slice_contour_data becomes a list of length 2 on the iterations where cv.fillPoly fails.

e.g.

image

EDIT: Ohhhhh, is it because in this case polygons[0] is only a 2x1-sized array, thus there are not enough points for fillPoly? If that's true... have you seen this behaviour before? Is there a workaround you would suggest? I guess some logic or a check of some sort could be implemented to exclude small polygons from the slice_mask? But then I'm not sure how reflective the mask is of the original contours... Hmm. Be good to get your thoughts.

tomaroberts commented 1 year ago

Ok, I've actually added some code to the get_slice_mask_from_slice_contour_data() function which effectively ignores any contours consisting of <3 points:

def get_slice_mask_from_slice_contour_data(
    series_slice: Dataset, slice_contour_data, transformation_matrix: np.ndarray
):
    # Go through all contours in a slice, create polygons in correct space and with a correct format 
    # and append to polygons array (appropriate for fillPoly) 
    polygons = []
    for contour_coords in slice_contour_data:
        reshaped_contour_data = np.reshape(contour_coords, [len(contour_coords) // 3, 3])
        translated_contour_data = apply_transformation_to_3d_points(reshaped_contour_data, transformation_matrix)
        polygon = [np.around([translated_contour_data[:, :2]]).astype(np.int32)]
        polygon = np.array(polygon).squeeze()
        polygons.append(polygon)

    # <---INSERTED
    # Find and remove any polygons which consist of less than three points, otherwise cv.fillPoly() generates exception
    polygons_to_remove = []
    for idx, poly in enumerate(polygons[:]):
        if len(poly) < 3:
            polygons_to_remove.append(idx)
    polygons = [p for idx, p in enumerate(polygons) if idx not in polygons_to_remove]
    # end INSERTED--->

    slice_mask = create_empty_slice_mask(series_slice).astype(np.uint8)
    cv.fillPoly(img=slice_mask, pts = polygons, color = 1)
    return slice_mask

I'm now able to run rtstruct.get_roi_mask_by_name() on previously failing regions. Not sure if this is a potentially unstable approach? I'm not super familiar with RTStructs and contours.

plesqui commented 11 months ago

Hello @tomaroberts,

Thanks for catching this up, and apologies for the delay in responding. The original maintainers that created the tool have become very busy with other endeavors and our group is now finding new contributors (like myself) to help moving forward this tool.

Regarding your analysis, I think you found the reason for the exception (need at least 3 contour points to fill) and your approach to mitigate it sounds reasonable to me. Although perhaps the bigger question is... where do these 2-point contours come from and were they supposed to be like that? Maybe, if we catch them at the creation step we can avoid writing them so that fillPoly does not complain downstream.

Either way, we can place a temporary patch until we understand the root cause. I'll work on a PR based on your sample code to patch this issue for now.

Would you be able to share some sample data that caused this issue with me so that I can test the fix?

Thanks again!

Pedro

tomaroberts commented 11 months ago

Hi @plesqui

Thanks for the reply and update. Good luck in your new contributor role!

Regarding where the 2-point contours come from... I'm not super familiar with RTStructs and ContourData, but I can explain the workflow I have. I am creating an app based on TotalSegmentator. The outputs from TotalSegmentator are multi-slice 2D voxel masks saved to NIFTI files. I import the masks into Python as numpy arrays which I feed into rt-utils (relevant code if for reference).

I assume at the edges of the segmentation regions it's possible to have slices which only contain a couple of pixels, so when rt-utils performs its mask to contour conversion, you are left with 2 contour points in those peripheral regions.

Hope that makes some sense. Whether it's the true source of the issue, I'm not 100% sure.

tomaroberts commented 11 months ago

PS: will look into sharing some data, sorry forgot to mention.

plesqui commented 11 months ago

Thank you for the explanation!

In that case, it makes sense that 2-point contours would be generated from such masks. I was worried that the 2-point contours were generated from valid masks due to a internal bug in rt-utils (which could still be possible).

I am debating if this post-processing step where invalid contours are removed should fall within rt-utils (like in your proposed approach), or if the user should make sure that the input masks are valid (i.e., won't have isolated pixels or other kind of noise that would cause these problems downstream). I am just worried that if we (rt-utils) discard some of the input data like these odd contours, we might loose relevant information unintentionally.

Let me think a bit more about it. Maybe we can throw a warning so that the user knows that contour data was lost and can go back and review why. What would you think?

tomaroberts commented 11 months ago

@plesqui

I get your quandary about whether to do this on the rt-utils side – makes sense – definitely don't want to potentially erase information.

It's not clear to me how to fix it entirely on the input mask side, hence it I think it would be a useful feature to erase 2-point contours within rt-utils.

Suggestion: how about putting a function into rt-utils such that the user can optionally erase 2-point contours? Perhaps a bit like the code I have in post above? I guess it would need to be accessible from rtstruct.add_roi()?

Perhaps something like:

rtstruct.add_roi(
  mask=MASK_FROM_ML_MODEL, 
  color=[255, 0, 255], 
  name="RT-Utils ROI!",
  erase_small_contours=True
)

Or however you think best to integrate with existing codebase.

MMdeB commented 11 months ago

Hi, I am writing an analysis script also based on totalsegmentator and have the same exact problem. I also export the segmentations from total segmentator to an RT struct using rt-utils. If I then want to load these RT structs I get an error on some, but not all. I think it would be great if this feature could be added to rt-utils!

jrenslo commented 8 months ago

Hello! I encountered the same problem with contour slices containing only a single point. There is an example in the open data set on cancerimagingarchive.net.

Collection: CPTAC-PDA Subject: C3N-00512 RTSTRUCT: Pre-Dose, Pancreas -1 (there are two, I believe both have the issue)

I have attached the manifest file for downloading this data (they provide a downloader application) manifest-1698234614986.tcia.txt (just remove the .txt extension to use the file)

As for suggested changes, I think a warning + omitting the faulty contours when building the mask would be appropriate, in the interest of robustness to input, but as a novice coder I defer to your judgment. At the very least a more informative error message would be helpful.

I would suggest editing the following function to check for faulty contours and raise an error. https://github.com/qurit/rt-utils/blob/8ba5cf64bd0ea5119523da0bafd10174115eebe4/rt_utils/rtstruct_builder.py#L42

mffield commented 8 months ago

Hi all, I encountered the same issue with reading contours from RTSTRUCTs and I was going to create another issue until I saw this. It arose with clinical contours where on some slices the region shrank to a single point.

I am glad others have found solutions. Here is code (showing whole function and inserted section) for how I dealt with the offending contours. The problem to me seemed that the function had an issue with the array changing shape (probably the squeeze operation just before doing this?). I just forced it to remain a consistent shape. It still means a polygon cannot be filled as it is less than 3 points like others have pointed out. But it at least passes through fillpoly and would put in the point that was in the RTSTRUCT rather than fully ignoring it. Ideally we'd also raise a warning when this situation has occurred.

def get_slice_mask_from_slice_contour_data(series_slice: Dataset, slice_contour_data, transformation_matrix: np.ndarray):
    # Go through all contours in a slice, create polygons in correct space and with a correct format 
    # and append to polygons array (appropriate for fillPoly) 
    polygons = []
    for contour_coords in slice_contour_data:
        reshaped_contour_data = np.reshape(contour_coords, [len(contour_coords) // 3, 3])
        translated_contour_data = apply_transformation_to_3d_points(reshaped_contour_data, transformation_matrix)
        polygon = [np.around([translated_contour_data[:, :2]]).astype(np.int32)]
        polygon = np.array(polygon).squeeze()
# Inserted ---------------------
        if len(polygon.shape)<2:
            temp = np.zeros([1,2])
            np.put(temp,[0,1],[polygon[0],polygon[1]])
            polygon = temp.astype(np.int32)
# End Inserted -------------------
        polygons.append(polygon)
    slice_mask = create_empty_slice_mask(series_slice).astype(np.uint8)
    cv.fillPoly(img=slice_mask, pts = polygons, color = 1)
    return slice_mask

I hope something like this can be included in the main package that can smoothly handle these situations because it will be more common as use of autosegmentation increases.

Thanks, Matt