I couldn't find a good all-in-one function to calculate DVH parameters from dose maps + structure masks (GTV, OAR, ...), so I made one myself. Any feedback is very welcome.
The primary purpose of this function is to calculate DVH parameters, like D99%, V40Gy, D0.5cc and the like. In my experience, the actual DVH itself is desired less often, but since it needs to be calculated anyway before parameters can be extracted, the function can also return that for free.
This function is supposed to be very "Matlab-native", meaning it needs and returns typical Matlab objects such as logical arrays and tables. Therefore, it is expected you have already imported the dose files (which typically begin as RTdose dicoms) as arrays, converted structure contours (RTstruct dicoms) to binary arrays, and so forth.
A limitation of this function is that it cannot calculate multiple DVHs for multiple structures at the same time. However, it is fast enough that it should be no problem calling it multiple times for every structure. As long as the underlying dose map remains the same, you can plot the resulting DVHs on top of each other; the dose bins (DVH x-axis) is calculated using the dose map maximum.
Name-value pair arguments are used to obtain the DVH parameters, e.g.
'Dperc', 99 -> returns the minimal dose to 99% of the volume
'Dvol', 0.5 -> returns the minimal dose to 0.5 cc
'Vperc', 45 -> returns the volume percentage receiving at least 45 Gy
'Vvol', 30 -> returns the volume in cc receiving at least 30 Gy
The rest of the syntax and input arguments are fairly straightforward, please see the code. Most important things to remember:
- Structure mask must be a logical and same size as the dose map
- Voxel sizes are required and are in cm
>> [params, dvh_values] = dvh(dose_map, gtv_mask, [0.3; 0.3; 0.3], 'Dperc', 99, 'Vvol', 45);
>> plot(dvh_values(:, 1), dvh_values(:, 2))
>> params = dvh(dose_map, gtv_mask, [0.3; 0.3; 0.3], 'Vperc', 45)
>> [~, dvh_values] = dvh(dose_map, gtv_mask, [0.3; 0.3; 0.3], 'Normalize', false);