I am actively working on the commissioning of an in-house planning software for HDR brachytherapy. Recently, we encountered surprising observations regarding the gamma index calculations. I am wondering if you can help. I will explain the issue as best as I can.
The idea here is that we want to ensure that the isodose lines displayed in our software match with the ones displayed in the TPS (in this case Oncentra Prostate (OcP) from Elekta). The method we want to use is the gamma index analysis with the 3D doses, because it tells us how far in dose and distance we are from the ground truth (the TPS). When conducting testings with the pyMedPhys gamma module, we obtained poor results for some cases (60% of voxels with gamma <= with 2%/2mm thresholds). Therefore, we visually investigated these cases to understand this result.
For one example case, iso_ocp.png (joined file) shows how the isodose lines look like in the TPS. We added a reference point (lets call it ‘P1’, a point inside source) to have a baseline. Next, we tried to replicate this image based on the exported RTDOSE file of OcP. To read the RTDOSE file, we used pymedphys.dicom.zyx_and_dose_from_dataset as explained in the tutorial and we get ocp_3D_dose_iso.png, where the axes are in mm and the colorbar in Gy. The image was plotted with matplotlib.pyplot.imshow, by using the extremum coordinates (zyx) of the axes as the extent. In this image, the star represents ‘P1’ in the coordinate system extracted from the RTDOSE file by zyx_and_dose_from_dataset. Therefore, we can conclude that there is a shift of a few mm. This is problematic, since this dose distribution is actually the ‘ground truth’ or the reference when conducting the gamma analysis. When recalculating the 3D of this case based on the RTPLAN with our software, we get the image iso_gmco, which seems to be correctly placed. To verify whether there is a bug in the RTDOSE file or in the way pymedphys is reading the file and handling the coordinates, we repeated the process with SlicerRT, which was validated in a paper published in 2012. The image slicerRT_iso shows the imported case (RTDOSE and RTSTRUCT exported from OcP) and the calculated isodose lines with ‘P1’. As seen, slicer correctly placed the 3D dose so that the isodose lines displayed in slicerRT visually matched with the ones displayed in OcP. Therefore, we conclude that there are some unkown issues in pymedphys.dicom.zyx_and_dose_from_dataset. This explains why we can obtain bad results when conducting the gamma analysis. Can you verify on your side if you observe such shifts as well? I wrote to one of the slicerRT maintainer to hope to have some insights on this.
changing imshow to pcolormesh did not do the trick. Regarding the 3D dose, in fact, my software computes it (iso_gmco) since it basically optimize a plan, calculate DVH, 3D dose (evaluation 3D dose), and then I send the RTPLAN with dwell positions and dwell times to Oncentra Prostate for the recalculation of the 3D dose in the TPS (reference 3D dose). I export the recalculated 3D dose from Oncentra, and then I read the RTDOSE with pymedphys for gamma analysis. The mismatch in coordinates seem to be present whether using pcolormesh or imshow. Here is basically what I plot (excluding the anatomy)
reference = pydicom.read_file(RTDOSE_file, force=True)
tps_coordinates, tps_3D_dose = pymedphys.dicom.zyx_and_dose_from_dataset(reference)
fig, ax = plt.subplots()
fig.gca().invert_yaxis() # the ax is inverted in Oncentra
extent_tps = [np.min(tps_coordinates),
im = plt.pcolormesh(tps_coordinates, tps_coordinates, tps_3D_dose[index], shading="nearest")
ax.contour(tps_3D_dose[index], levels=list(isodose_levels.values()), colors=list(isodose_levels.keys()), extent=extent_tps)
index just gives me the slice I want to show. I use contour of matplotlib to display the isodose, which are the classic ones (75%, 100%, 125%, 150% and 200%). We need to define extent for the contour in order to have the coordinates of the isosurfaces.
Also, just to let you know, I have done the gamma analysis for 40 cases with Slicer (1%/1mm, 2%/2mm, 3%/3mm). We will be able to compare afterwards the gamma results obtained with the module of pymedphys.
I encountered this bug since the voxels of the 3D dose with Oncentra Prostate are not exactly cubes. The dimension is different in the x axis and y axis. This is why I observed this problem with Oncentra Prostate, but not Oncentra Brachy, where the voxels are cubic. You can do further testings on your own by using voxels with different x and y size.
Also, just to let you know, in Oncentra Prostate, the GridFrameOffsetVector is like :  0, -1.03125, -2.0625, …, -82.5. I have to somewhat flip the z axis and the 3D dose (axis 0) to be in the right coordinates.
I will make the changes as requested by the end of next week. I am pretty busy right now. I am also investigating dicompyler-core because some of my colleagues encountered weird results regarding dvhcalc with OcP, so this project might also have this serious issue. I will fix this bug in both projects and let you know.