Is there any reason that the GetIsodosePoints would not work in 3D?

[sama2689]

I am trying to find the 3D centroid of an isodose given an rtdose file. So what I have done for this is loop over all the planeZs, extracted the isodose I am interested in, and computed.

#Extract z coordinates of all planes by using body plane
body_roi=plu.get_structure_id(rtss.GetStructures(), 'BODY') #find the roi number of the structure called 'BODY'
body_coords=rtss.GetStructureCoordinates(body_roi)
plane_Zs=np.array(list(body_coords.keys()),dtype=np.float32) #extract and cast to a numpy array

The code above extracts the values of planeZs in the body structure. I thought I could then simply extract whatever isodose I was interested in and then compute centroid by just averaging out the points.

# Extract this section finds the centroid of a given isodose level
isodoselevelGy=1.0 #Dose in Gy that is desired
threshold=round(isodoselevelGy/float(rtdose.ds.DoseGridScaling)) #scale threshold to be in the dicom file's desired units.

for z in plane_Zs:
isodoseZ=rtdose.GetIsodosePoints(z,threshold,5) # this will give list of (x,y) points
isodoseZ=[point+(z,) for point in isodoseZ] # this adds z values to the (x,y) tuples to put them in (z,y,z) format

#if isodoseZ is not empty, add it to the global list of all isodose points found in previous planeZs
if isodoseZ:
isodosepoints[i:i+len(isodoseZ)]=isodoseZ
i=i+len(isodoseZ)

Is there any reason this should not work? Do I have the syntax correct for extracting the 1Gy isodose from a parts rtdose file? I tried this method and for some reason the centroid of the extract isodose has enormous disagreement with the centroid as determined by Eclipse. I was just wondering if this was due to some fundamental issue or if there is likely to be a bug somewhere in my code. I found the centroid just using the npmean of the isodose points in each direction,

centroid=np.around([np.mean(isodosepoints[:,0]), np.mean(isodosepoints[:,1]), np.mean(isodosepoints[:,2]) ],2 )

[cutright]

Your centroid calculation has an assumption of uniformly distributed points built-in, but I don't believe this is guaranteed in the isodose calculation.

For example, if you had a square defined by four points, one at each corner, your equation would give you the centroid as you expect. However, if you have the same square with 100 more points all on the right edge of the square, the shape has not changed but your centroid calculation would place the centroid very heavily to the right.

I'd recommend using shapely's centroid calculation per slice and then average those. This should be OK since slices are on regular intervals.