pyplati / platipy

Processing Library and Analysis Toolkit for Medical Imaging in Python
https://pyplati.github.io/platipy/
Apache License 2.0
103 stars 24 forks source link

Structure volumes calculated by platipy are systemetically higher compared to clinical output #218

Open YYAN97 opened 1 year ago

YYAN97 commented 1 year ago

Hi @pchlap,

I have ran some dosimetric analysis on clinical RTdose (with RTstruct) using your package. The code was nice and smooth. I was able to get the DVH and corresponding metrics. When I compared the results (11 contours, 44 DVH metrics) calculated by the code with clinical treatment planning system output (MIM, Ohio), I found that while the DVH metrics were very similar (differences were mostly within 2%), the volumes of structures calculated by the code were systematically higher compared to MIM output (a few percents, up to 30% for small structures). This might suggest that your code deals with the RTstruct in a different way from the clinical system, which might potentially affect the dosimetric analysis despite the DVH metrics are very similar in my experiment. This would be less ideal since we want to match everything with clinical practice. Could you provide some insight into this? Does your code process RTstruct in a consistent way with the clinical systems? What might cause this deviation and how to fix this? Is it a common observation that python packages (yours and maybe others) can get a systemetically higher structre volumes compared to clinincal systems?

Thank you and looking forward to your response.

pchlap commented 1 year ago

Hi @YYAN97,

thanks a lot for posting this really great question. Hopefully, I can shed a little bit of light on the topic.

Firstly, it can be difficult, or impossible, to know exactly how commercial systems are computing things like DVH metrics. Since their source code isn't available we would need to rely on reading their documentation or communicating with their support teams to gain more insight.

However, I believe the main cause of discrepancies here is the voxel size used for computation. PlatiPy converts the contours represented in points from the RTSTRUCT DICOM into a binary voxel-wise mask. We use the same voxel resolution as the image for this, however another approach would be to use a finer grid for better accuracy. I know from working on MIM extensions that they use a grid size 2x the resolution of the image.

Using a coarser grid size would explain why volumes are systematically larger when using platipy. So one solution could be to experiment using a finer grid to confirm that this does indeed reduce this effect.

There may however be other approaches to computing volumes from contour points as well as DVHs and their metrics. I know that the dicompyler Python tool does this, would be interesting to compare this as well. This tool has been archived and is no longer maintained but the authors are working on an alternative. https://github.com/bastula/dicompyler/blob/ffba3a83dbbc56ee3b10ca43abcf959fca4c962f/dicompyler/dvhdata.py#L45

I'd also be curious to hear @rnfinnegan's thoughts on this. We have discussed this in the past but I suppose the upshot for us was that this wasn't a major concern for our research work. However as you point out the effect does increase for smaller structures so it would be interesting to work towards solving this.

rnfinnegan commented 1 year ago

Hello,

This is a great question @YYAN97!

Different grid size could definitely be a culprit, as Phil mentioned. Another likely cause is how structures are represented.

As @pchlap said, the format of RTSTRUCT "contour" structures is a list of points that correspond to vertices on a polygon. This has the benefit that the structure resolution is essentially decoupled from the image space, and a contour can go through any part of a voxel. In many clinical tools. images are interpolated to the screen resolution to make the images look nicer, and hence make contouring a bit easier. When the volume is calculated the software can then just compute the volume enclosed by the 3D polyhedron defined by the vertices. Dose is a bit more complicated, since a contour might overlap some arbitrary fraction of the dose grid voxel. Some systems could calculate dose by scaling by the fraction of the dose voxel enclosed, e.g. if 40% of a dose grid voxel with value of 10 Gy is enclosed, then that would contribute 4 Gy (40% x 10 Gy).

In platipy, everything is an voxelised array - images, dose grids, structures, deformation vector fields, etc. This has the benefit that image operations (masking, cropping, applying transformations) is the same for every time of object. But it does unfortunately mean that information can be lost when converting from other representations, with RTSTRUCTs being the most obvious example.

I don't think there is a perfect solution. For example, consider the case of automatic segmentation. Methods like CNNs operate directly on the voxelised array, and the output is also a voxelised array. It wouldn't make much sense to convert this to an list of vertices (like RTSTRUCT format) because we would potentially lose information.

Some ideas that might work:

I'd like to continue this conversation, it's a good challenge and one that probably has different optimal solutions for different use cases 🚀

YYAN97 commented 1 year ago

Hi @pchlap and @rnfinnegan,

Thank you both so much for the response! I agree that MIM should have some algorithms that smooth the contours and partially account for the volumes of the voxels on the boundary when calculating volumes/doses. I like @rnfinnegan 's ideas especially the second one which seems like a balance between accuracy, flexibility and computation complexity/power. In terms of dose calculation, I think the current apporach will not lead to a substantial deviation except for small structures with steep dose fall-off. Looking forward to your update if you decided to improve on this issue.