Error in matRad_generateStf
justingm opened this issue · 9 comments
Hello!
I created a 10x10x12 mm3 water phantom with a 2x2x2 mm3 target located at a depth of 6 mm (following matRad_example1_phantom.m). The voxel size is 0.1 mm.
Since I will use matRad to create proton treatment plans for animals, both the water phantom and target are small. I also created a new proton base data using the Bortfeld's formulation for the IDD (matRad_getDDDfromAnalyCacl.m) and Highland's formula for the sigma (mb_calcMCSsigma.m) for energies that are appropriate for animal irradiation (12.78 MeV to 46.94 MeV w/c is equivalent to a range of 20 mm to 2 mm). Since in reality, small beams are very difficult to create, we will use a collimated beam (1 mm collimator) so I set the machine.data.initFocus.SisFWHMAtIso to 1 mm w/c corresponds to a sigma of 0.43. I put this as a constant for all energies (I hope that is not gonna cause an issue)
Basically in the plan that I have created, I just shoot protons in one direction. I set the bixelWidth and longitudinalSpotSpacing to 0.5. When I try to calculate the plan, I get an error in matRad_generateStf file
(line 276) focusIx(k) = find(machine.data(vEnergyIx(k)).initFocus.SisFWHMAtIso > currentMinimumFWHM,1,'first');
with a message
"Unable to perform assignment because the left and right sides have a different number of elements."
I am not really sure what the problem is.
Best regards,
Justin
PS I am using the master branch
Dear @justingm ,
the the focus assignment can be a very annoying thing. Basically, matRad allows multiple optical beam settings to be selected given the spot spacing (bixelWidth) and FWHM at iso. There's a look-up table to control this in machine.meta (LUT_bxWidthminFWHM) which is used to select a minimum FWHM at Iso (second row) for a given spot spacing (first row). If your minimum FWHM there is too large, this error can occur. I think you can circumvent this if you just set the second row to zeros (not sure) in the machine file.
To check if this is the issue indeed causes the problem, you can also set a breakpoint in that line and just evaluate machine.data(vEnergyIx(k)).initFocus.SisFWHMAtIso > currentMinimumFWHM when you reach it. If this result is empty or has a weird dimension, this is indeed the cause.
@amitantony can we implement a check for this and either throw a meaningful error or, even better, throw a warning and select the next best beam focus?
Yes, I am aware of that so I also made sure that the 2nd row of LUT_bxWidthminFWHM is 1 1.
I tried your suggestion of zeroing it out and now I do not see that error anymore.
However, I am now seeing a new error because of matRad_calcLateralParticleCutOff
(line 172) dose_r = matRad_calcParticleDoseBixel(depthValues(j) + baseData.offset, radialDist_sq, largestSigmaSq4uniqueEnergies(cnt), baseData);
I think it is caused by the fact that I am getting NaN values for the 4th column of energySigmaLUT.
I've encountered this problem before for the beam data base that I made for our cyclotron but I thought I just messed up the beam model so I just shifted to generating the IDDs and sigmas analytically. Unfortunately it still throws this error.
Ah the lateral cut-off... It is unfortunately a messed up thing... but was the solution to keep absoulte dosimetry when cutting of dose...
Your code seems to not find the largest initial beam width. Some things to check:
- in initFocus.dist & sigma, are the dist values the distances from the source in mm?
- Is the source to surface distance actually in the range of these dist values such that matRad can interpolate the beam width from that? Do the values in energySigmaLUT when reaching the loop in line 7 make sense for the first to the third column? (energy, index 1, sensible SSD)
- if you don't find an issue there, can you give me the exact values of the variables passed to the function call in line 82:
matRad/matRad_calcLateralParticleCutOff.m
Lines 82 to 84 in cb6b04a
The dist value seems to be causing the problem. I set it too low.
Thanks for all the help!
Hi! I do have something that I have been wondering about. From what I understood in the code, the depth-dependent sigma (of an infinitesimally small pencil beam) is convolved with the initFocus.sigma to get the actual lateral spread in the object. Now for example, I have modelled our beam using Monte Carlo and scored the 3d dose distribution in water to get the IDDs and depth-dependent sigma. Of course, this is not an infinitesimally small pencil beam. the depth dependent sigmas that I get are the real lateral spread of the beam. How should I set the initFocus.sigma then for this case?
one pragmatic way to do it is to take the fitted width of the first "layer" / depth coordinate and set this as initFocus.sigma. Then undo the convolution on the remaining tabulated depth dependent widths / sigmas to get the respective kernel width/sigma.
problem is, that you now only have one sigma and one position. So you either need more information from your beam model (that is, a few measurments of the width in air at different distannces from the source, which is exactly what's in initFocus.sigma/dist, or something like a phase space parameterization to use, e.g., the CS-Formula to compute the beam width. You could also assume a non-convergent/divergent beam, and just create a initFocus.dist = [0,Inf] & initFocus.sigma = [c,c] (for a constant beam width of c that is independent of the source to surface distance).
We have a complete model of our beamline in BDSIM so extracting Twiss parameters and beam profiles in air shouldn't be a problem. I was actually thinking of just taking the phase space information (x,y position, x,y direction and energy of individual particles) after the collimator and then, I will zero out the positions to create a pencil beam without losing the divergence information of the beam. Calculate the 3d dose distribution in water to extract the IDD and sigmas. Then, I will just set the initFocus.sigma to the initial beam width at the surface of the water phantom for the original beam (not pencil beam).
I am still a bit confused about the purpose of initFocus.sigma and initFocus.SisFWHMAtIso. It says in the wiki that it is for the modelling of the lateral spreading of the beam in air. I think for my case this is unnecessary because when I do MC simulations to get the IDD and sigmas, I start from the exit of the collimator which means that the scattering in air is already taken into account. Or am I confusing the purpose of this table?
if the spread in air information makes sense for your model depends on your collimator is.
Here's some comments regarding your approach / questions:
- initFocus.SisFWHMAtIso & initFocus.sigma are somewhat redundant, as the value could just be added to the sigma list at exactly SAD (spot width at isocenter). Yet we store the beam width at isocenter separately, because it corresponds to the value that a TPS would export to quantify the spot size at isocenter (Dicom Tag Scanning Spot Size, for example), so it is easy to assign and look up.
- The reason you need the initFocus.sigma depending on dist is that in the "normal" spot scanning IMPT plan, your surface is not always at the same distance, so you calculate the source-to-surface distance (SSD) and then look up how wide your beam is when it hits the patient surface, so you can afterwards convolve your kernel with this value. If you mount a collimator with, say, 1mm diameter on the surface, then this does of course not make that much sense, as your beam will always have the same width. If the collimator is at a fixed distance, the beam might still evolve after that until it hits the surface. In any case, then also the convolution of this surface Gaussian with the kernel is not accurate, as the fluence coming out of your collimator will not be Gaussian (it will probably look like a mixture of a truncated Gaussian and a circular uniform distribution or something like that). So you might need to interfer with the dose calculation to make a more accurate convolution with the surface fluence.
If you need more help, we could also set up a meeting to get a better insight and help you out with some implementation details. Maybe this results in some interesting code contribution to matRad. Write us an e-mail if you think this would help.
Ah yes, now it makes sense. I didn't really consider that the source to surface distance is not always the same. However, I don't think it would matter too much because we are dealing with animals and the collimator will be placed at a fixed distance of 20 mm from the surface of the animal. The changes would be very minimal so the effect wouldn't be large.
Anyway, thanks for the help. I am still trying to familiarize myself with how matRad works. We are most likely going to use it as TPS in our small animal irradiation facility (not much commercially available TPS for protons that can be used for the preclinical case and they are very expensive). I will send an email maybe next week when I have things more in order.