[BUG] Recalculating a plan in matRad imported from Eclipse yields dosimetric differences when MLCs are present.

[edwardwang1]

Describe the bug
When importing an RTPlan (with MLCS) and RTDose from Eclipse, recalculating the RTPlan within matRad will lead to differences in the dose compared to the imported one. The recalculation appears to be correct for open fields, but not for fields with MLCs. We have performed our experiments for a FFF 6MV kernel in using a cubic water phantom.

To Reproduce
Steps to reproduce the behavior:

1. Generate a machine specific pencil beam kernel using https://github.com/e0404/photonPencilBeamKernelCalc
2. In Eclipse (or presumably another TPS), create an IMRT plan (single or multiple beams). Add MLCs (easiest to just fit the MLCs to the target).
3. Calculate the dose in Eclipse
4. Export the dose and plan from Eclipse and import into matRad.
5. Recalculate the dose in matRad using the custom kernel and compare to the imported plan.

Expected behavior
For a water phantom, we expected the recalculated doses to be very similar to the imported doses.

Screenshots
1-beam IMRT, 10x10cm field, open field
DVH from imported dose

DVH from recalculated dose

1-beam IMRT, 10x10cm field, MLC fit to target
Gamma pass fraction at 2%/2mm over the phantom = 96.79%
DVH from imported dose

DVH from recalculated dose

2-beam IMRT, 10x10cm field, MLC fit to target
Gamma pass fraction at 2%/2mm over the phantom = 94.69%
DVH from imported Dose

DVH from recalculated dose

Desktop (please complete the following information):

• OS: [e.g. Windows]
• Environment & Version [e.g. Matlab 2023a]
• matRad master branch

Additional context
When recalculating the dose, we are using pln.propDoseCalc.useCustomPrimaryPhotonFluence = true;, but this does not seem to have any effect since the dij calculation is very fast.

[wahln]

Hi @edwardwang1
First of all - nice to see that you figured out the kernel generation and it works on the water phantom with an open field. And thanks for putting so much work into this, as at the time when we validated the dose calculation we didn't have FFF data to work with.
It would be good to see the dose distributions to see where the issues lie, especially for the single MLC-collimated beam. Since this is only phantom data, if there's no trouble sharing it, I could also take a look at it myself.
And just to make sure again - this is not some weird double-MLC like in the Halcyon, right?

[edwardwang1]

Hi @wahln

I very much appreciate all of your help along the way! There is one caveat with using the kernel that I've created. I've had to change line 62 of matRad_calcPhotonDoseBixel.m from
dose(dose < 0 & dose > -1e-14) = 0; to
dose(dose < 0) = 0;
Please find attached a .zip file that contains the .mat files of the examples in this issue, as well as our custom kernel. Let me know also if you prefer me to email it.

ExportedEclipsePlansAndDoses.zip

Also, for your convenience here is a script to compare dose profiles

function visualizeCrossbeamProfilePair(doseDistribution1, doseDistribution2, plane)
    % Scrollable comparison of the crossbeam profiles of two 3D dose distributions on the same plot
    %
    % Inputs:
    %   doseDistribution1 - 3D matrix representing the first dose distribution
    %   doseDistribution2 - 3D matrix representing the second dose distribution
    %   plane - The plane to extract the profile ('xy', 'xz', or 'yz')

    % Validate inputs
    if nargin < 3
        error('Three input arguments are required: doseDistribution1, doseDistribution2, and plane');
    end
    
    if ~ismember(plane, {'xy', 'xz', 'yz'})
        error('Invalid plane. Choose from ''xy'', ''xz'', or ''yz''.');
    end
    
    % Determine the number of slices in the specified plane
    switch plane
        case 'xy'
            numSlices = size(doseDistribution1, 3);
            xLabel = 'X-axis (cm)';
        case 'xz'
            numSlices = size(doseDistribution1, 2);
            xLabel = 'X-axis (cm)';
        case 'yz'
            numSlices = size(doseDistribution1, 1);
            xLabel = 'Y-axis (cm)';
    end
    
    % Initial plot for the first slice
    currentSlice = round(numSlices / 2);
    figureHandle = figure('Name', 'Scrollable Crossbeam Profile Comparison');
    
    plotHandle1 = plotProfile(doseDistribution1, plane, currentSlice, xLabel, 'Dose Distribution 1');
    hold on;
    plotHandle2 = plotProfile(doseDistribution2, plane, currentSlice, xLabel, 'Dose Distribution 2');
    hold off;
    
    legend('Dose Distribution 1', 'Dose Distribution 2');
    title(sprintf('Crossbeam Profile in the %s plane at slice %d', plane, currentSlice));
    
    % Add a slider to the figure for scrolling through slices
    uicontrol('Style', 'slider', ...
              'Min', 1, 'Max', numSlices, ...
              'Value', currentSlice, ...
              'Units', 'normalized', ...
              'Position', [0.15 0.02 0.7 0.05], ...
              'Callback', @(src, event) updatePlot(doseDistribution1, doseDistribution2, plane, round(get(src, 'Value')), plotHandle1, plotHandle2, xLabel));
    
    % Add a text box to display the current slice number
    sliceText = uicontrol('Style', 'text', ...
                          'Units', 'normalized', ...
                          'Position', [0.45 0.1 0.1 0.05], ...
                          'String', sprintf('Slice: %d', currentSlice));
                      
    % Nested function to update the plots when the slider is moved
    function updatePlot(doseDistribution1, doseDistribution2, plane, sliceIndex, plotHandle1, plotHandle2, xLabel)
        % Update the first plot
        set(plotHandle1, 'YData', extractProfile(doseDistribution1, plane, sliceIndex));
        
        % Update the second plot
        set(plotHandle2, 'YData', extractProfile(doseDistribution2, plane, sliceIndex));
        
        % Update title and slice text
        title(sprintf('Crossbeam Profile in the %s plane at slice %d', plane, sliceIndex));
        xlabel(xLabel);
        ylabel('Dose (Gy)');
        
        set(sliceText, 'String', sprintf('Slice: %d', sliceIndex));
    end

    % Nested function to extract and plot the profile
    function plotHandle = plotProfile(doseDistribution, plane, sliceIndex, xLabel, plotTitle)
        profile = extractProfile(doseDistribution, plane, sliceIndex);
        plotHandle = plot(profile, 'LineWidth', 2);
        hold on;
        xlabel(xLabel);
        ylabel('Dose (Gy)');
        grid on;
    end

    % Nested function to extract the profile based on the plane and slice index
    function profile = extractProfile(doseDistribution, plane, sliceIndex)
        switch plane
            case 'xy'
                % disp(round(end/2));
                % disp(sliceIndex)
                profile = doseDistribution(round(end/2), :, sliceIndex);
            case 'xz'
                profile = doseDistribution(:, sliceIndex, round(end/2));
            case 'yz'
                profile = doseDistribution(sliceIndex, :, round(end/2));
        end
    end
end

You can use it like so after loading in the mat file

resultGUI2 = resultGUI;
% Recalc dose using matRadGUI to referesh resultGUI with the new calc
visualizeCrossbeamProfilePair(resultGUI.physicalDose, resultGUI2.physicalDose, 'xz');
visualizeCrossbeamProfilePair(resultGUI.physicalDose, resultGUI2.physicalDose, 'yz');

[wahln]

Thanks.
Just to report for now - I can load and visualize your data.
Taking a first look, I have to say I expected worse. The things that strike me are:

• First look it still looks to me as there is some issue with the overall penumbra of the matRad kernels. have you tried playing with the penumbra-setting during kernel fitting?
• There seem to be some weird things happening at some boundaries. I need to investigate this in more detail.
  I guess these are the issues annoying you as well?

By the way, the usePrimaryPhotonFluence setting has no effect because when fields are imported and recalculated, matRad is no longer in "beamlet mode", but calculates the field as a whole (indicated by bixelWidth = 'field'). Here, the primary fluence is always used, because the impact of the fluence convolution is negliglibe for field computation as you only need to do one convolution with the primary fluence anyway.

[edwardwang1]

Thanks for taking a look so quickly!

I remember adjusting the FWHM of the kernel, but I recall it didn't have a significant effect. However, I didn't do this rigorously, so I can definitely go back and look at it again. I am hoping to create a simple test case that I can use to test the DICOM exporter you sent me previously, however in order to do so I need to ensure that a plan exported from Eclipse and into matRad matches.

I have a question about "field" vs "beamlet" mode. I understand that when I import an open field into matRad, it simply performs a field calculation. However, in the case of a collimated beam (please see attached image), the MLCs form a non-rectangular shape, so I don't understand how a "field" calculation is appropriate, as certain bixels would be blocked by the MLCs. Since we're using a cube as a target, this is not the best demonstration, but if you look at the corners, you can see that the opening by the MLCs is not fully square.

Lastly, I forgot to mention in my previous comment that we are using a Varian Truebeam equipped with HDMLCs. There are 60 leaf pairs. The first 14 and last 14 leaf pairs have 5mm width, and the middle 32 have 2.5mm width.

[wahln]

Okay, regarding the "field" setting, don't take the name too literal. matRad will take the MLC and other device limiters from the DICOM, and define the field as a mask image of the collimated field / shapes. This mask will then be used on the primary fluence, which will then be convolved with the kernel. The "field" setting is thus the right choice here.

I am not completely sure about the different leaf sizes, but I that matRad defines the mask on a much finer grid than the actual leaf sizes, so it should not be an issue for recalculation.

Now to the penumbra. I am quite certain that this is a major part of the differences.
There's two important values here for the PB dose calulaction, the source-to-collimator distance in machine.meta.SCD and the machine.meta.data.penumbraFWHMatIso defining the penumbra of a field at isocenter (which needs careful definition with FFF fields I guess). During dose calculation, these values are used, among others, to determine the correct Gaussian filter applied to the collimated fluence. In your case there's SCD = 0 and penumbraFWHMatIso = 10 stored, and it seems like this is off (changing these values directly in the machine changed the calculation results a bit). With a little bit of playing around I could get the DVHs closer to each other and better gamma pass rate.

Here's an example for the 1 Beam MLC after recalculation with a penumbra of 5 mm and an SCD of 500

Also note that I did this on the most recent "dev" branch. In the commits just pushed, I also added the option to provide pln.propDoseCalc.ignoreInvalidValues = true, which then allows you to use your kernel without editing the code. Apart from that, it also fixes some GUI display bugs such that you can load your data in the GUI without issues.
To work with the new branch, there's a coordinate transformation necessary, as matRad now interprets the "isoCenter" to be defined in the world coordinate system:

%Transform into correct coordinate system on the new branch
pln.propStf.isoCenter = matRad_cubeCoords2worldCoords(pln.propStf.isoCenter,ct);
for i = 1:numel(stf)
    stf(i).isoCenter = matRad_cubeCoords2worldCoords(stf(i).isoCenter,ct);
end

%Now run your dose calculation and ignore nan / negative values
pln.propDoseCalc.ignoreInvalidValues = true;
resultGUI2 = matRad_calcDoseForward(ct,cst,stf,pln);

%Run matRad's compare dose
matRad_compareDose(resultGUI.physicalDose,resultGUI2.physicalDose,ct,cst);

%Put the recalculated dose in the resultGUI for GUI display
resultGUI.physicalDose_recalc = resultGUI2.physicalDose;

%Run the GUI
matRadGUI

Probably these values are also not correctly stored by the kernel fitting code.

[edwardwang1]

Hi Niklas,

Thanks for pointing this out to me. The SCD should definitely not be 0, that's a mistake on my part. I asked our physics associate and apparently for the Varian TrueBeam it's 550mm. I will note that the results are much less sensitive to the SCD than the FWHM.

To determine the FWHM, I iterated through a list of values at 0.1 increment for the single collimated beam case and found it to be 5.4. The results are much better for both the single and double collimated beam cases (pretty similar to what you have). For the next test, I imported another collimated 2-beam case, but this time the collimation is much more complex. Here is how I obtained it:

1. Use matRad to create an IMRT plan
2. Export the IMRT plan from matRad (using the code you sent me previously - this is also something I want to talk about, but probably best in a separate discussion).
3. Import the plan into Eclipse
4. Calculate the dose in Eclipse
5. Import the plan and dose into matRad

As you can see first, this added collimation worsens the results quite drastically. Do you have any thoughts as to why this would be the case? I have attached the .mat file (in zip format) as well as a zip of the raw DICOM.
exportedmatRadPlanCalcInEclipse.zip

CUBE_DICOM.zip

NEW_photons_TB_E.zip

Thanks!

[wahln]

I need some time to look into it.
I'm a bit confused what exactly you compare at the moment, can you also make a native IMRT plan in Eclipse and try to recalculate it with matRad with your data set?
I also started a shared discussion thread on external machine integration for photons: #772

[edwardwang1]

Hi Niklas,

I have added my experience to the discussion. Also, it was a good idea to natively calculate an IMRT plan in Eclipse, the results are pretty good. Now I just need to figure out why my other plan doesn't work.

To clarify, in the other plan, I obtained MLC leaf positions through loading an RTPlan generated with matRad. Then, I calculate the dose in Eclipse using loaded plan. From here, I export the dose and plan generated from Eclipse and import back into matRad. Except for how the MLCs were determined, there is nothing else different between the native Eclipse plan and the matRad plan, so I'm not sure why I'm seeing dose discrepancies.

NativeEclipse2BeamIMRTSlidingWindow.zip

[edwardwang1]

I've developed a theory as to what could be causing the discrepancy. I noticed that in the native Eclipse plan, each beam had 166 control points, whereas in the plan optimized within matRad, there are only 9 control points. It may be possible that Eclipse and matRad handle interpolating between control points differently.

To test this, I modified the matRad RTPlan exporter to interpolate between the control points. Now, the exported plan has 161 control points when importing into Eclipse. When repeating the analysis now the results look much better (although there is still a larger discrepancy than compared to the native Eclipse plan).

[wahln]

Just bumping this Issue. Is there anything weird still going on / are you at the same state as a month ago?

[edwardwang1]

I haven't done any additional experiments yet with the water phantom. I did experiment with some lung cases, but there is pretty large discrepancy so I am waiting for the VMC++ license at the moment.