Hi there,
I had dosimetric indices that were 0 so I dug up on the code and I think there is a problem. On the below dictionnary you will see that the dosimetric indices of the heart have a problem. We cannot have a D1cc at 0 and D0.1cc > 0.
{
"V100%(PTV)": 38.336,
"V150%(PTV)": 13.946,
"V200%(PTV)": 5.655,
"CI(PTV)": 0.8126274903013004,
"HI(PTV)": 0.6102973991047177,
"D90%(CTV)": 4.567464192708334,
"CI(CTV)": 0.011427092990664281,
"HI(CTV)": 0.5674990396734081,
"D2cc(Skin)": 1.9542788984763677,
"D1cc(Skin)": 2.0903842139790076,
"D0.1cc(Skin)": 2.357177734375,
"D2cc(Chestwall)": 2.380239835349462,
"D1cc(Chestwall)": 2.5163130326704546,
"D0.1cc(Chestwall)": 2.718912760416666,
"D2cc(Lungs)": 1.4137550636574077,
"D1cc(Lungs)": 1.4971816748903508,
"D0.1cc(Lungs)": 1.629638671875,
"D2cc(Heart)": 0.22303549091467698,
"D1cc(Heart)": 0,
"D0.1cc(Heart)": 0.23793003015350878
}
The problem is that you are identifying the first index which is bigger than your threshold but after that you are taking an interval that does not include this threshold to interpolate your dosimetric index, resulting in potential negative dosimetric indices, which you then turn into 0 if they are negative. But a volume or a dose, should never be negative. Here is an example with prints to highlight the problem in a D2cc computation, we can see the interval we are interpolation into is [1.693, 1.913] but 2cc is not included in this interval…
Here is my proposed fix in computeDcc() and computeDx() functions, that ensures there is never a negative volume or dose interpolated:
def computeDx(self, percentile:float, return_percentage:bool=False) -> float:
"""
Compute Dx metric (e.g. D95% if x=95, dose that is reveived in at least 95% of the volume)
Parameters
----------
percentile: float
Percentage of volume
return_percentage: bool
Whether to return the dose in Gy on % of the prescription
Return
------
Dx: float
Dose received in at least x % of the volume contour
"""
index = np.searchsorted(-self._volume, -percentile)
if (index > len(self._volume) - 2): index = len(self._volume) - 2
volume = self._volume[index - 1]
volume2 = self._volume[index]
if (volume == volume2):
Dx = self._dose[index]
else:
assert volume >= percentile and volume2 <= percentile, (f"Volume interpolation error: volume1: {volume}, volume2: {volume2}, requested percentile: {percentile}")
w2 = (volume - percentile) / (volume - volume2)
w1 = (percentile - volume2) / (volume - volume2)
Dx = w1 * self._dose[index] + w2 * self._dose[index + 1]
if Dx < 0:
raise ValueError("Dx computation resulted in 0 dose. ")
if return_percentage:
assert self._prescription is not None
return (Dx / self._prescription) * 100
return Dx
def computeDcc(self, x:float, return_percentage:bool=False) -> float:
"""
Compute Dcc metric (e.g. D200cc if x=200 for minimal dose that is received the most irradiated 200cm^3)
Parameters
----------
x: float
Volume in cm^3
return_percentage: bool
Whether to return the dose in Gy on % of the prescription
Return
------
Dcc: float
Dose received
"""
index = np.searchsorted(-self._volume_absolute, -x)
if (index > len(self._volume) - 2): index = len(self._volume) - 2
volume = self._volume_absolute[index - 1]
volume2 = self._volume_absolute[index]
if (volume == volume2):
Dcc = self._dose[index]
else:
assert volume >= x and volume2 <= x, (f"Volume interpolation error: volume1: {volume}, volume2: {volume2}, requested volume: {x}")
w2 = (volume - x) / (volume - volume2)
w1 = (x - volume2) / (volume - volume2)
Dcc = w1 * self._dose[index] + w2 * self._dose[index + 1]
if Dcc < 0:
raise ValueError("Dcc computation resulted in 0 dose. ")
if return_percentage:
assert self._prescription is not None
return (Dcc / self._prescription) * 100
return Dcc
I am not sure if computeVg() and computeVx() functions are also affected since they don’t use interpolation, I have not tested. My group and I have our own fork of your repository which we would happily discuss merging.
Best,
SĂ©bastien
), you might see that I already merged it in the merge request above (after some discussions and help from your colleague Hossein), I am just awaiting for review from my colleagues.