The question about primary fluence used in matRad_calcPhotonDose

[Little-Stars520]

In photon computing, there is a comment stating: "% Toggle custom primary fluence on/off. If set to 0, we assume a homogeneous primary fluence; if set to 1, we use measured radially symmetric data." However, this segment of code does not actually function as intended. I would like to ask: why is a uniform fluence of 1 (unity fluence) used instead of the diagonal fluence measured at dmax ? In reality, shouldn’t the fluence be inhomogeneous? Wouldn’t adopting a uniform fluence of 1 lead to calculation errors? And where is the inhomogeneity of the primary fluence applied?

[wahln]

I would say this switch is somewhat a mix between a half-executed idea and some internal testing.

The general motivation between using a homogeneous fluence is to spare time in the convolution (only executed once for all beamlets).

For LINACs with flattening filter this makes sense because, as you say, you can use a homogeneous primary fluence and then do the slight weighting according to the lateral diagonal distance during or after dose calculation. I think the reweighting was never really implemented as we don't have so many use cases to generate deliverable, device-correct plans with matRad.

For FFF LINACs the diagonal weighting will be too rough. you need to actually convolve the corresponding fluence patch for the beamlet. There's where the switch comes from, to be able to trigger the behavior in those conditions, but you are correct that this currently does not work as intended (The convolution is not correctly performed at the patch coordinates but on the central ray).

In general, the full field convolution centered around the beams central ray is always used in forward calculation from shapes loaded from DICOM.

Would be nice to clean this up:

• The diagonal reweighting should be reliably done when homogeneous fluence is used during convolution. Either directly during dose calculation or the corresponding weight needs to be stored in the dij.
• when the switch is turned the kernels should be convolved at the beamlet coordinates with the corresponding fluence patch at that position
• probably the machine needs a FFF metadata field to toggle the default behavior. Then the switch can be used as intended to run with/without homogeneous fluence reweighting to observe the distance.

Anybody got time for that? :-D

[wahln]

Okay, I still don't full understand where the "d_max" reference comes from. The "primary fluence" is determined at a "suitable" depth. Conceptually, this depth does not make a difference for the algorithm itself (as the energy-fluence is assumed depth-invariant and we are in the fan-line system). However, this assumption (depth dependent kernel convolved with fixed primary fluence) creates an approximation error due to spectrum hardening or softening. I took a quick look at the Bortfeld 1993 paper again, commenting on this:

The error made by this assumption [of lateral invaraiance of the energy spectrum] can be reduced if the energy fluence is defined in a mean treatment depth in tissue instead of air. We found that for 15 MV photons and small or medium sized [...} fields, the error was negligible. However, [...][for wedges and compensators], an error of about 2.5% must be expected.

So one could argue about using the max TPR point, a more deeply measured fluence for deep seated tumors, and so on. Will all lead to different slight approximation errors in different treatment scenarios. But such approximation errors are most certainy overshadowed by the other errors made by a PB algorithm in heterogeneous patients.

[Little-Stars520]

The "d_max" is derived from Martin Siggel's diploma thesis provided by Matrad. Below is the original text from the thesis along with the accompanying figure.


In my view, the maximum point of TPR, i.e., TMR, is attributed to the fact that absolute dose is calibrated at the Dmax point. As a matter of fact, I still do not fully understand the significance of selecting the dmax depth for primary fluence. If we consider the common tumor depth, wouldn't d10cm be a better choice?
Furthermore, I do not quite understand the "lateral invariance of the energy spectrum" you mentioned here. Do you mean that the off-axis ratio of fluence does not change with the variation of depth? Another hypothesis is the lateral invariance of the PB kernel, meaning that the same kernel is used at different lateral positions.
@wahln

[wahln]

yeah,as said, depth at maximum dose is just a practical choice, with the repercussions as explained above. And the maximum field size is the closest practical approximation to infinity.

Basically the concept of this pencil-beam algorithm assumes that the primary (energy) fluence retains its shape and spectrum in the fan-line system, as all scattering contributions are captured in the kernel. Pratically this means that all primary photons travel straight on the depth-coordinate in the fan-line system with constant energy spectrum. So in the fan-line system, I can cut out an x/y plane in BEV at any depth and always see the same shape of primary fluence (and yes, off-axis ratio of fluence does not change).
My choice of depth for determining primary fluence then determines where I have the best approximation accuracy, so to say. One could say that maximum dose makes sense, because there the dose is relatively the largest so my approximation should work best. You could hypothetically also say "I only treat prostate" and take the depth of the average prostate tumor.

The hypothesis of lateral invariance then ensures that from there the computation of the dose is just a simple 2D convolution in beams eye view with the kernel. If the energy spectrum is laterally invariant, then the kernel can be applied at each lateral position in the beams eye view system.

It is not necessary to stick to these assumptions, they are just simplifying the computation a lot. You could also have different kernels for a laterally changing energy spectrum and modify the computation accordingly. It just doesn't seem to be worth the effort.

[Little-Stars520]

Thanks for your explanation — I think I understand your point now. The choice of depth essentially determines where we want to achieve the most accurate dose distribution.

I revisited the matRad code and the Bortfeld 1993 paper. Regarding the primary fluence, Bortfeld describes it as follows:

It mentions an inverse square correction to make it independent of the source distance. I'm puzzled by this correction — is it done by multiplying with the (dref).^2 /r.^2 , r is the distance from the source to the voxel where the fluence is defined? dref may be the dmax? Is this correction intended to match the fanline system's parallel photon beam assumption? This seems consistent with matRad's use of a uniform fluence (i.e., all values are 1). I might have misunderstood. Here, "being independent of the source distance" means that the off-axis ratio does not change with the depth, though I have no idea how this is achieved, it might merely be a conversion from the gantry coordinate system to the fan-line coordinate system.

Additionally, Bortfeld's paper mentions:

I'm unclear on how the inverse square factor is implemented here. In matRad’s matRad_calcPhotonDoseBixel function, the following code is used.

Is SAD applied because the coordinate system origin is fixed at the SAD (i.e., the source-to-origin distance is set to SAD)? In fact, I don't quite understand why the inverse square correction is applied using SAD.
This correction seems to allow the profile at the same depth to be simply adjusted, rather than using primary fluence to correct for off-axis dose.

@wahln

[wahln]

The purpose of the inverse square correction is the transformation to the fan-line system (a fan-beam originating from a point-source loses intensity with 1/r^2, so correcting for that transforms you into a system where the intensity doesn't change with the depth). It is not related directly to the shape of the fluence.

When we do a forward calculation from imported shapes, we do exactly what is described in your first paragraph:
We compute the primary fluence stored in the machine radially and multiply it with the transmission factor (which is basically the MLC shape convolved with the gaussian representing a finite source) to get the collimated fluence:

matRad/matRad/doseCalc/+DoseEngines/matRad_PhotonPencilBeamSVDEngine.m

Lines 461 to 486 in 0994481

	if this.useCustomPrimaryPhotonFluence || this.isFieldBasedDoseCalc
	% overwrite field opening if necessary
	if this.isFieldBasedDoseCalc
	F = ray.shape;
	else
	F = this.Fpre;
	end
	% prepare primary fluence array
	primaryFluence = this.machine.data.primaryFluence;
	r = sqrt( (this.F_X-ray.rayPos(1)).^2 + (this.F_Z-ray.rayPos(3)).^2 );
	Psi = interp1(primaryFluence(:,1)',primaryFluence(:,2)',r,'linear',0);
	% apply the primary fluence to the field
	Fx = F .* Psi;
	% convolve with the gaussian
	Fx = real( ifft2(fft2(Fx,this.gaussConvSize,this.gaussConvSize).* fft2(this.gaussFilter,this.gaussConvSize,this.gaussConvSize)) );
	% Get kernel interpolators
	ray.interpKernels = this.getKernelInterpolators(Fx);
	else
	ray.interpKernels = this.interpKernelCache;
	end

Note that in a beamlet-based calculation for planning, we do not do this. This is the issue with uniform fluence. And looking at the code now, I think, that setting the useCustomPrimaryFluence switch does actually perform the weighting (sorry for my confusion earlier).

Regarding the transformation of the system with the SAD. The isocentric plane is the reference plane for the transformation (i.e. the lateral distances agree between both systems there).