Calculate DVH and DVH parameters (D99%, V40Gy, D0.5cc, ...)

A function to calculate the cumulative dose-volume histogram (DVH) and DVH parameters from a dose map and binary structure mask.

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.

Use

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

Examples:

        
>> [params, dvh_values] = dvh(dose_map, gtv_mask, [0.3; 0.3; 0.3], 'Dperc', 99, 'Vvol', 45);
>> params
   params =
     1×2 table
       D99% (Gy)    V45.0cc (cc)
       _________    ____________
        40.711         29.889   
>> plot(dvh_values(:, 1), dvh_values(:, 2)) % To plot the DVH

>> params = dvh(dose_map, gtv_mask, [0.3; 0.3; 0.3], 'Vperc', 45)
    params =
      table
        V45Gy (%)
        ________
         45.782 

>> [~, dvh_values] = dvh(dose_map, gtv_mask, [0.3; 0.3; 0.3], 'Normalize', false); % get an absolute DVH with cc values