FFF Kernel Creation & Clarification on Input Data

[edwardwang1]

Hello,

I am attempting to use this repository to create an 6MV FFF kernel for use in matRad. I have obtained the relevant input data, and was able to run the ppbkc.m script to generate the machine data. However, the calculated kernels seem incorrect (see Figure), and dose calculation in matRad fails. I suspect that there is an error in my input data.

To be specific, here are the steps to generate the input data.

1. TPR: We calculated TPR from PDD data by simply normalizing the data to the value at 5cm (the standard at our institution). Our measured data only goes down to 30x30mm field size. I understand that we need to perform an extrapolation to obtain the TPR at 0 field size. What was the type of model that was fitted to perform the extrapolation?

2. OF: We have the OF for a variety of field sizes, however it seems that this repository only takes in one field size. Based on the of.dat files in the referenceData and synergy folders, the field size appears to be 10cm. Is this correct always, or should we be picking a different field size based on our other measurements?

Below, please find the other images generated using the ppbkc.m script.

Thank you

[edwardwang1]

Just wanted to follow up with my progress thus far. Let me know if this issue is more suited for the main matRad repo.

It appears that dose calculation was failing because in matRad_calcPhotonDoseBixel, some of the dose values are less than 0. I see that matRad already performs a correction factor for numerical instability.

dose(dose < 0 & dose > -1e-14) = 0

I have changed it so that all values below 0 are set to 0, and was able to perform dose calculation without issue.

dose(dose < 0) = 0

Obviously this is a non-ideal solution, but I'm not sure why I'm getting negative dose values as large as -3e-5.

My colleague has built a script to export plan data from matRad into Eclipse for dose calculation. It appears that even with this error, our FFF matRad machine file is more accurate than the generic one.

One last note, I am not setting the pln.propDoseCalc.useCustomPrimaryPhotonFluence filed to True (as described here: #6). I found that if I did set this field to true, the Dij calculation takes much longer.

I would love to hear about others' experience with generating an FFF machine file!

[edwardwang1]

Hi @markbangert @wahln @mingersming, apologies for tagging, just doing so in case notifications for this repo are disabled.

I have some more updates on this project.

1. I have corrected my output factor, since I realized it is a function of r (which is a function of x and y). However making this correction did not substantially change my results.

2. My colleague and I have created a cube water phantom in matRad (based on example 1). We exported the phantom as an NRRD using matRad, converted it to a DICOM, and imported it into Eclipse. We then performed a dose calculation in matRad using our custom kernel, and compared it to the dose calculation in Eclipse. To generate the dose in matRad, we skipped optimization and used calcDoseDirect

w = ones(stf.numOfRays,1);
resultGUI = matRad_calcDoseDirect(ct,stf,pln,cst,w)

We are seeing a large discrepancy in centerline PDD between the two cases. The PDD fall off is much steeper in matRad than in Eclipse. Please see the figures below. Does anyone have any insight on what could be causing this? We are happy to send our kernel/and or base data for analysis in case there is something weird with our data.

[wahln]

Hi @edwardwang1,
sorry for not responding - I didn't get notifications for this repo, as you suspected, and thus didn't see your progress.
A few months ago, we used the code to do some preliminary tests on creating machine data for an FFF device (Ethos - are you trying, by chance, Ethos/Halcyon as well?). I don't recall the exact process we went through, but I remember that we had some issues in getting everything to run (and I even remember some annoying stuff with negative dose values), but eventually we figured it out and got a somewhat good enough dataset. I can look into this end of next week or beginning of the week after that.

So far I have the following comments:

• pln.propDoseCalc.useCustomPrimaryPhotonFluence = true; is, in this case, important. There is a shortcut in pencil-beam dose calculation that you can take if you assume that the fluence is homogeneous and you have one beamlet size -> you only need to do one convolution of the pencil beam kernels with the primary fluence by beam. If the fluence is dependent on the lateral position in the beams eye view (as it is substantially for FFF), you need to do the convolution for each individual beamlet. So for dosimetric verification, this parameter needs to be on (even if the dose calculation takes longer). You could also play the beamlet weights afterwards to approximate the effect of the inhomogeneous fluence, but we didn't implement this yet.
• It is weird that useCustomPrimaryPhotonFluence does have basically no effect in your case. For it to work, make sure that machine.data.primaryFluence is correctly populated in the machine file.
• The kernel's you show above don't look unreasonable. The reference kernel's from matRad look a bit different because they used a less stable radial integration (which was, I think, to mimick some previous code from a commissioned clinical system).
• In general, your data looks reasonable. I have to figure out if we did some substantial changes to the code as soon as I get a chance to look at it

[wahln]

OF: We have the OF for a variety of field sizes, however it seems that this repository only takes in one field size. Based on the of.dat files in the referenceData and synergy folders, the field size appears to be 10cm. Is this correct always, or should we be picking a different field size based on our other measurements?

The output factor should be normalized to the field 10x10 cm^2, but the OF depends on the field size (even in your graph). I don't get the question.

[wahln]

TPR: We calculated TPR from PDD data by simply normalizing the data to the value at 5cm (the standard at our institution). Our measured data only goes down to 30x30mm field size. I understand that we need to perform an extrapolation to obtain the TPR at 0 field size. What was the type of model that was fitted to perform the extrapolation?

In my stuff it seems to be normalized to 10cm depth at 990 SSD (so normalized in the isocenter), but as long as the correct SSD and depth are set in the script it should work anyways. But just to make sure - are you confident the TPR is correct? Since it is not only a normalization, since PDD is usually measured at fixed SSD with moving detector, while TPR is defined for fixed SAD) and consequently if corrections need to be applied (e.g. inverse square).

I think the extrapolation to 0 is done within the script.

[edwardwang1]

Hi @wahln

Thank you for the reply! We are currently trying to build a machine file for Varian Truebeam. However, our center is getting an Ethos soon, so we would like to move towards that eventually as well.

Regarding the output factor, our physicist collaborator gave us this table below. However, the the OF.dat in the base data only has 2 columns. Originally, to create OF.dat we took the values from the table at a single column (there by fixing the field size in the X-direction). I'm pretty sure this is incorrect. What we are doing now is taking the diagonal values of the table, and calculating r=sqrt(x^2 + y^2). If this is also incorrect please let me know.

I will give the TPR a more thorough look. We are using the data that was collected when the machine was first commissioned, so I assume the necessary adjustments were made, but I will double check with the physics team.

[wahln]

Ah yes, the output factor is not perfectly symmetric in x/y, which we assume in the script (the script only considers square fields).
So this indeed introduces a small approximation error.

Regarding your PDDs between Eclipse and matRad, Something is fishy there, maybe at import / depth plot. If you look at the intensity plot, the dose reduced to around 30% of the max value at 200 mm depth in water within matRad, but in your line-plot, it some gets close to zero after 160mm.

[wahln]

I think your TPR would make sense compared to the also shallow falloff of your depth dose. But something in the comparison is fishy, as said above.

[wahln]

By the way if you want to have faster calculation time, you can also try to create a very large single bixel. Something like pln.bixelWidth = 50 (depending on your field size you want to validate).
You can also set up the stf for that manually (not to difficult if you come from 0 degree).

[edwardwang1]

Hi Niklas,

We have fixed our issue with the water phantom PDDs not matching. Embarrassingly, we were using a different sized cube between matRad and Eclipse. After matching the cube sizes, the PDDs became very similar. We still have negative doses during matRad_calcPhotonDoseBixel, but hopefully flooring them to 0 will not affect our downstream tasks. Thanks for your help.

[edwardwang1]

Hi all,

@wahln @markbangert @mingersming (tagging in case notifications are off)

Reopening this issue because I have some further questions. All experiments below are run in the master branch.

I have modified my protocol such that I directly importing an Eclipse plan into matRad, to avoid any possibility of setup error. I am using a 80x80x80cm phantom with a 40x40x40cm PTV inside it.

After importing the CT, RTPlan, RTStruct and RTDose from Eclipse, I perform the dose calculation in matRad using:
pln.propDoseCalc.useCustomPrimaryPhotonFluence = false;
pln.machine = 'TBFFF_CustomFinal';
w = ones(stf.numOfRays,1);
resultGUI = matRad_calcDoseDirect(ct,stf,pln,cst,w);

I have been able to reproduce my results that my custom machine in matRad matches the exported Eclipse dose in both cross beam profile and the PDD at 10x10cm field size.

At 40x40cm field size, there are some slight deviations in the crossbeam profiles, but not the PDD.

Lastly, I created a IMRT plan for lung in Eclipse, and imported it into matRad. Then, I used the matRad GUI to recalculate the dose (with the "Recalc" button). The animation below is the line profile in the y-direction. There are some differences that are visible, which I believe is due to the heterogeneity in the lung. I see that matRad has a heterogenietyCorrection branch - I will look into it to see if it improves the results.

Following this issue,, I tried setting `pln.propDoseCalc.useCustomPrimaryPhotonFluence' to True, but this did not have any effect in the matRad calculated dose (the doses are identical) for the water phantom. This is the case for my custom machine, as well as the generic machine. I also tried changing the flag when recalculating the lung IMRT case - it also had no effect.

I am happy to send my waterphantom mat file in case someone wants to try and reproduce!

[edwardwang1]

I realized I should probably add a single beam lung case.

Crossbeam Profile:

PDD:
The PDD is taken at the axial slice that passes through the isocenter. The isocenter is at approximately 200mm in the x (left-right) direction, and the animation below shows the PDD f rom 150mm to 250mm.

[wahln]

Hi @edwardwang1,
the heterogeneity correction branch does not really have anything to do with photon dose calculation in Lung. The pencil-beam algorithm will not yield good results within lung, especially after the tissue/lung interfaces. One could introduce some kernel correction, but within matRad we are not focusing so much on photon dose calculation development, unfortunately.

I'm still a bit worried about the lateral water profile differences. First of all, the dose calculation seems to cut the profile way to early. The penumbra also seems too sharp in the low-dose region.
Could you try recalculating with the current dev branch (there are some changes for the coordinate system there which will probable need readjustment of the isocenter).

Regarding the parameter pln.propDoseCalc.useCustomPrimaryPhotonFluence having no effect: If you import a plan, afterwards the bixelWidth might be set to 'field'. When this is the case, matRad will always use the primary photon fluence in the machine. The idea is that in this case we directly import shape information and shape weights, which only make sense combined with the fluence. For planning dose calculation, you might want to do the beamlet dose calculation with a homogenous fluence (save time on the convolutions) and later just reweight according to the fluence at the beamlet position - however, this may not be such a good idea for FFF beams. Probably we should implement some explanation for this / rethink the default settings / make this dependent on the machine...

[edwardwang1]

Hi @wahln,

Thanks for the response. Good eye! I didn't even notice the profile being cut off. I have repeated the water phantom crossbeam profile analysis using the dev branch, and it seems that the profile is still being cut off (animations below). I used matRadGUI in the dev branch to re-import my DICOM to take care of any differences in coordinates.

The code that I used was:

stf.machine       = 'TBFFF_CustomFinal'; 
w = ones(stf.numOfRays,1);
resultGUI = matRad_calcDoseDirect(ct,stf,pln,cst,w);

I believe I came across a small bug. Line 158 of matRad_PhotonPencilBeamSVDEngine.m this.intConvResolution = this.collimation.convResolution; throws an error because the 'this' object doesn't have the 'collimation' field.

To fix this, I simply added the line below to the matRad_PhotonPencilBeamSVDEngine(pln) function in the same file.
this.collimation = pln.propStf.collimation;

Additionally, running the pencil beam calculation gives the following warning. Perhaps it's related to the profile being cut off? Warning: Kernel Cut-Off 'Inf mm' larger than machine data range of '179.500000 mm'. Using '179.500000 mm'!

Regarding 'pln.propDoseCalc.useCustomPrimaryPhotonFluence', you are correct that the bixelWidth is indeed being set to 'field' when I import the plan. For my future dose calculations, I will set the flag to true (and bite the bullet when it comes to computation time.

Lastly, how is field size specified/calculated in matRad? I know for a single bixel, field size is equal to the bixel width. In Eclipse, the user can specify a field size. How would I do this in matRad for a beam composed of multiple bixels?

10x10:

40x40:

[wahln]

It might indeed be related to the kernel cut off, but I am not sure. Could also be the

geometricLateralCutOff;     %lateral geometric cut-off in mm, used for raytracing and geometry
dosimetricLateralCutOff;    %relative dosimetric cut-off (in fraction of values calculated)

Those can also be defined in pln.propDoseCalc, are defined in the abstract base class matRad_PencilBeamEngineAbstract and have default values of 50 mm and 99.5%, respectively. You can check their default values by constructing an engine object:
DoseEngines.matRad_PhotonPencilBeamSVDEngine() and look at the properties.

I would suggest trying to increase the geometric cut-off before thinking about wider kernels.

When you use bixels / beamlets for IMRT, matRad does not really consider available field size (or field size limits set with jaws, for example). More generally, matRad does not consider machine parameters related to beam limiting devices (if you look in the machine file, there is also no information regarding this). We rely on the user then to make sure what they compute lines up with their machines.
I would really like to expand on this functionality in matRad at some point, but currently we don't really have the resources.

And thanks about the error report, I will fix it on the dev branch (check e0404/matRad@1e107b3).