Gamma index calculation - possible issues

cedricbelanger · 29 October 2021 16:21

Hi Smon,

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.

Best regards,

Cédric Bélanger

iso_gmco :

slicerRT_iso :

ocp_3D_dose_iso :

iso_ocp :

SimonBiggs · 29 October 2021 20:55

Hi @cedricbelanger,

Thanks for your bug report, and thank you for all the detail you’ve gone into to help reproduce it. It’s massively appreciated .

I have a request regarding the following:

Just to cover our bases can you change from using imshow to instead using the following:

import matplotlib.pyplot as plt
plt.pcolormesh(x, y, dose, shading="nearest")

and then let’s go from there.

Cheers,
Simon

SimonBiggs · 29 October 2021 21:00

Also, you might be interested, a while back (a long while back) I wrote some code to read brachy DICOM files, pull out the seeds, run a TG43 calc and compare the results to the original DICOM dose:

cedricbelanger · 30 October 2021 13:03

Hi @SimonBiggs,

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[2]), 
		       np.max(tps_coordinates[2]), 
		       np.min(tps_coordinates[1]), 
		      np.max(tps_coordinates[1])]
im = plt.pcolormesh(tps_coordinates[2], tps_coordinates[1], 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.

Thank you for looking into this,

Cédric

SimonBiggs · 30 October 2021 22:15

Hi Cédric,

Can you make one more coordinates related change in your plotting just to make sure everything is as expected, can you use the following for your contour plotting instead of extent:

z, y, x = tps_coordinates
dose = tps_3D_dose[index]
ax.contour(
    x,
    y,
    dose,
    levels=list(isodose_levels.values()),
    colors=list(isodose_levels.keys()),
)

cedricbelanger · 1 November 2021 15:06

Extent was adding a shift, but there is still a shift to resolve !

ocp_3D_dose_iso_2|690x433

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.

Cédric

SimonBiggs · 2 November 2021 03:15

Hi Cédric,

Good to hear we resolved that component .

Under the hood zyx_and_dose_from_dataset gets its coords from the xyz_axes_from_dataset function:

github.com

pymedphys/pymedphys/blob/17b99b22bc5868666648fd23fb0dfcd3e39de761/lib/pymedphys/_dicom/dose.py#L33-L38

    
      
          def zyx_and_dose_from_dataset(dataset):
              x, y, z = xyz_axes_from_dataset(dataset)
              coords = (z, y, x)
              dose = dose_from_dataset(dataset)
          
          
    return coords, dose

Here is that coords function:

github.com

pymedphys/pymedphys/blob/17b99b22bc5868666648fd23fb0dfcd3e39de761/lib/pymedphys/_dicom/coords.py#L55-L187

    
      
          def xyz_axes_from_dataset(
              ds: "pydicom.dataset.Dataset", coord_system: str = "DICOM"
          ) -> Tuple["np.ndarray", "np.ndarray", "np.ndarray"]:
              r"""Returns the x, y and z axes of a DICOM dataset's
              pixel array in the specified coordinate system.
          
          
    For DICOM RT Dose datasets, these are the x, y, z axes of the
              dose grid.
          
          
    Parameters
              ----------
              ds : pydicom.dataset.Dataset
                  A DICOM dataset that contains pixel data. Supported modalities
                  include 'CT' and 'RTDOSE'.
          
          
    coord_system : str, optional
                  The coordinate system in which to return the `x`, `y` and `z`
                  axes of the DICOM dataset. The accepted, case-insensitive
                  values of `coord_system` are:

This file has been truncated. show original

Would you be able to directly copy that function into your code and see if you can simplify it right down to just the component that appears to be causing the issue?

Cheers,
Simon

SimonBiggs · 2 November 2021 03:16

@Matthew-Jennings, if you’re available might you be able to be of assistance here?

cedricbelanger · 2 November 2021 12:32

Hi @SimonBiggs,

I found the bug in xyz_axes_from_dataset :

di = float(ds.PixelSpacing[0])
dj = float(ds.PixelSpacing[1])

should be :

di = float(ds.PixelSpacing[1])
dj = float(ds.PixelSpacing[0])

This is because PixelSpacing[0] is associated with the rows (i.e. y) and PixelSpacing[1] is associated with the columns (i.e. x), see Pixel Spacing Attribute – DICOM Standard Browser (innolitics.com) :

col_range = np.arange(0, ds.Columns * di, di)
row_range = np.arange(0, ds.Rows * dj, dj)

Here how it looks like with the change :

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.

Best regards,

Cédric

cedricbelanger · 2 November 2021 12:42

Also, just to let you know, in Oncentra Prostate, the GridFrameOffsetVector is like : [81] 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.

cedricbelanger · 2 November 2021 12:48

New problem, I just ran with the copied pasted function with the resolved bug, and now I get an error for one case :

reference = pydicom.read_file(RTDOSE_file, force=True)
	tps_coordinates, tps_3D_dose = pymedphys.dicom.zyx_and_dose_from_dataset(reference)
	print(tps_coordinates[0].shape, 
		tps_coordinates[1].shape,
		tps_coordinates[2].shape)
	tps_x, tps_y, tps_z = xyz_axes_from_dataset(reference) # -> function I copied and pasted from your post
	tps_coordinates	= (tps_z, tps_y, tps_x)
	print(tps_coordinates[0].shape, 
		tps_coordinates[1].shape,
		tps_coordinates[2].shape)

Which gives me the following output :

(73,) (65,) (81,)
(73,) (65,) (82,)

Why is there a difference? I would guess this is a problem with the np.arange function with voxel size with a lot of decimal numbers? In this case,

print(di, dj)

gives 0.9825 1.015625

cedricbelanger · 2 November 2021 13:09

Changing

col_range = np.arange(0, ds.Columns * di, di)
row_range = np.arange(0, ds.Rows * dj, dj)

with

col_range = np.array([i * di for i in range(ds.Columns)])
row_range = np.array([i * dj for i in range(ds.Rows)])

did the trick!

Best regards,

Cédric

SimonBiggs · 3 November 2021 21:41

cedricbelanger:

Changing

col_range = np.arange(0, ds.Columns * di, di)
row_range = np.arange(0, ds.Rows * dj, dj)

with

col_range = np.array([i * di for i in range(ds.Columns)])
row_range = np.array([i * dj for i in range(ds.Rows)])

did the trick!

Yup, arange has it’s own problems and I’ve found it is best simply never used.

That is a pretty serious bug. Thank you for undergoing the investigation.

Might you be in a position to commit this change as a pull request (PR) to the GitHub repo? It would involve three steps:

creating a dummy plan that has a known dose point value for a given x, y, z value;
then use that point and the PyMedPhys coords and dose extraction to write a failing test demonstrating the issue;
and lastly then apply your above fixes to make that test pass.

I would be more than happy to guide you through any issues you might have in getting this set up. @sjswerdloff might you be in a position to help me onboard @cedricbelanger to create their first PR?

Cheers,
Simon

sjswerdloff · 4 November 2021 04:22

Sure, I can try to help with onboarding @cedricbelanger
I will walk through a fresh setup for myself this weekend.

SimonBiggs · 4 November 2021 05:05

Awesome thanks Stuart

cedricbelanger · 4 November 2021 16:14

Hi @SimonBiggs @sjswerdloff

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.

Best regards,

Cédric

SimonBiggs · 4 November 2021 19:39

Awesome, thanks @cedricbelanger.

If you’re interested, after submitting the PR, you can feel free to have your name and your hospital/uni’s logo up on the active contributors section of the docs:

That’s up to you, just let me know