Moving a Contour and then assigning to RT-Struct about dvh-analytics

[cutright]

I haven't written any code in DVHA for that, unfortunately you'll need to dive into DICOM land. This is also really a pydicom question, they might be able to help you out better than I can since I haven't done this before. But I do recommend downloading the dicompyler desktop application to peruse the DICOM structure of your structure file. This might be more complicated than it seems.

Your ROI lives in ds.ROIContourSequence[roi_index]. This contains an array of contours in ContourSequence, where each contour has a ContourData tag containing the points.

If you're only looking to apply an axial shift, this is relatively trivial to simply edit the CountourData tag directly. Involving the z-direction may take more trial and error for you. I'm pretty sure each slice has to live on a plane coincident with imaging data? But again, new territory for me.

[cutright]

For sure in DVHA, editing Z like that will fail. DICOM is really only 2.5D. As far as I know, all commercial TPS's require contours to be coincident with imaging data (but if they can, fairly certain they have to conform to imaging slices on export). This is why their sup/inf expansions look so terrible haha

It might help if you could provide some details about what you're trying to do ultimately? Do these edits really need to go back in DICOM? It's likely a TPS won't interpret the way you hope (if at all). If you plan to go to completely 3D, these Shapely methods will likely not be appropriate.

[cutright]

I don't see why not. Are you converting your structure data to a DICOM compliant format? You should be able to edit the pydicom dataset directly (dicompyler-core stores this in its ds property), and then call save_as on this pydicom dataset.

[sama2689]

Sorry, to clarify, I'm unsure of how the assignment would happen (i.e. how I would update the struct file itself rather than just extracting the values to the gtv_planes dict and then altering the dict?) .

I have extracted gtv_planes and made the necessary alterations by looping through, but am not sure how to now update all relevant properties of the dicom file to reflect this new translated contour. Is there some function for this I've missed?

[sama2689]

Thanks for the guidance, If a slice has to live on a plaine coincident with imaging data, does that mean that applying translation via editing the output of get_planes_from_string isn't valid?

My translation function is

def Contour_Translate(contour_planes,shift_vector):
    """
    Applies translational shift to all points in a contour plane

    Parameters
    ----------
    contour_planes : Dict
        a "sets of points" formatted dictionary of input contour.
    shift_vector : List
        A list of size 3, specifying the contour shift vector as [x,y,z].

    Returns
    -------
    contour_planes_translated : Dict
        Output contour planes after translation, formatted as a 
        "sets of points" dictionary.

    """
    contour_planes_translated={}
    for key,plane in list(contour_planes.items()): #iterate over dictionary, k are keys representing each plane
        key=str(float(key)+shift_vector[2]) #the new key should match the z location
        #print(key)
        plane_copy=copy.deepcopy(plane)
        for point_index in range(0,len(plane[0])): #loop over points in this plane
            plane_copy[0][point_index]=list(np.around(np.array(plane[0][point_index])
                                       +np.array(shift_vector),decimals=2)) #add shift vector to all points    
        contour_planes_translated[key]=plane_copy
    
    return contour_planes_translated

I am calling this function in a simple example

from dicompylercore import dicomparser, dvh, dvhcalc
from dvha.tools.roi_geometry import overlap_volume
from dvha.tools.roi_formatter import dicompyler_roi_coord_to_db_string, get_planes_from_string

dp = dicomparser.DicomParser("RT-STRUCT.dcm") #parse RTStruct dicom

gtv_coords=dp.GetStructureCoordinates(16)
gtv_coords_str=dicompyler_roi_coord_to_db_string(gtv_coords)
gtv_planes=get_planes_from_string(gtv_coords_str)

gtv_translated=Contour_Translate(gtv_planes, [-1,-1,-1])


overlap=overlap_volume(gtv_planes, gtv_translated)

When I apply this function to apply nonzero translation in z, I notice that I get 0 overlap_volume compared with the pre-translation contour. Could this be because the z direction has been changed which is somehow setting the if len(thicknesses) > 0: in overlap_volume to false and so is returning 0?

[sama2689]

Sure, I suppose the edits don't really need to go back into the DICOM. My long term objective is to determine a translation vector that will maximize volume overlap between two GTV's (which i will call GTV_base and GTV_ada) from 2 different struct files. GTV_base stays fixed and I search for a translation vector to apply to GTV_ada such that volume overlap between the two is maximized. I was doing this step by step (as you may have noticed from the various issues you have helped me with in the past week 😊) as follows

1. Apply several different translations to GTV_ada
2. Compute Overlap of all translated versions of GTV_ada with GTV_base and return the translation vector of the maximal overlap

So coding-wise I've written a translation function which I thought could do this, but I suppose the z-editing is proving difficult.

[cutright]

That's actually a fairly complicated problem. If you want to stick with Shapely and DVHA code, you'll need to write a slice interpolator. I've written something for that before, but it's slow and technically incorrect.

I iterated through relative perimeter space along two z-adjacent contours/polygons, then did a 3D linear interpolation between those two points (of same relative perimeter space distance) to get the next point on the interpolated slice. This is approximately true if the two adjacent contours are not wildly different. Shapely offers methods to sample your polygon in this manner (you can also do it in numpy like I do with DVHA's DTH calculation in v0.9.7). But you also have to be careful about where you start your polygon for this slice interpolation. I wrote this up in a manuscript, but ultimately chose not to publish it... saving its application for a new upcoming project, stay tuned =)

Ultimately, if you want to do it the "right" way, you need to pixelize your contours into a binary mask, then use some imaging interpolation (e.g., with scipy). Alternatively, you could voxelize your contours. DVHA has code for this (it's how I calculate OVH). But I expect it to be insanely slow for the resolution you'd want. You'd also have to write code to calculate overlap between voxelized structures (I have no code for that).

Hope that helps...

[sama2689]

Ah, by slice interpolator you mean to say that I'd have to make something so that the planes can be kept the same even in the case of a translation with z, and then call overlap volume? That is indeed a fairly complicated problem, but I'll give it a try. Fortunately in our case speed isn't really a huge issue as we are not dealing with too many patients at once.

I had tried to go the route of extracting the binary mask, but had some trouble extracting it =/. I can write code for overlap of voxelized structures I think. Fortunately I am not restricted by speed or memory so I think binary mask may be the way to go.

[sama2689]

I have decided to keep any z shifts in-plane for now and that seems to fix the issue, So long as the increment of the shift vector is in units of the slice thickness no interpolation is needed and overlap_volume works out.

[sama2689]

Could the dicompyler-core's interpolate_between_planes function do this? As time is not a factor I thought I could generate an interpolated dict of the translated GTV and then use the planes in that as inputs to the interpolate_between_planes function.

[bastula]

@sama2689 Actually it won't, that method is poorly named and doesn't interpolate contours but makes copies and inserts between slices. See the last line of that function 😳

[sama2689]

By the way, I have found a workaround to this problem. I extract the point cloud itself using dicmpyler core and then compute the convex hull created by the two point clouds using scipy's convex hull function. This convex hull function can then be combined with a translation function which takes in some input shift vector to create a function that translates a contour and computes the convex hull volume of the two contour point clouds.

This function can be sent to scipy's minimize and then becomes an optimization problem as a function of a shift vector applied to a point cloud. The shift vector which minimizes the convex of the point cloud consisting of two structures is also the shift vector which will create maximal overlap between the two. Essentially, this produces an automated way to extract the couch shift needed to maximize dose cube GTV overlap with fraction CT. Technical note on this is currently being drafted, thank you so much @cutright and @bastula for your help on this.