Trying to understand how to get kernel for matrad

[Erhushenshou]

Hi Mark,
I am new to TPS calculation algorithms. I have trouble in reading your code in ppbkc.m as to generate kernel.

%%
[tprMax,tprMaxIx] = max(tpr);
meanMaxPos_mm     = ceil(mean(tprDepths(tprMaxIx)));

tpr_0 = tpr(:,1)/tprMax(1);

% compute mu for TPR0 with exponential fit, only use data points behind max
% (neglects build-up) 
[~,ix] = min(abs(tprDepths-meanMaxPos_mm)); %% the depth nearest to the meanMaxPos in tprDepth

fSx  = sum(tprDepths(ix+1:end));
fSxx = sum(tprDepths(ix+1:end).^2);
ftmp = -log(tpr_0(ix+1:end));
fSy  = sum(ftmp);
fSxy = sum(ftmp.*tprDepths(ix+1:end));

% mu = 0.005066; % reference value from literature
mu = (fSxy-((fSx*fSy)/length(tpr_0(ix+1:end)))) / ...
            (fSxx-(fSx^2/length(tpr_0(ix+1:end))));
        

%% compute betas

maxPos_fun = @(x) (log(mu)-log(x))/(mu-x);

beta(1) = fmincon(@(x) (maxPos_fun(x) - meanMaxPos_mm).^2,1,[],[],[],[],0,1000);
beta(2) = fmincon(@(x) (maxPos_fun(x) - (meanMaxPos_mm+1/mu)/2).^2,1,[],[],[],[],0,1000);
beta(3) = fmincon(@(x) (maxPos_fun(x) - 1/mu).^2,1,[],[],[],[],0,1000);

%% get weights for indivdual compoents of all TPRs
D_1 = (beta(1)/(beta(1)-mu))*(exp(-mu*tprDepths(ix+1:end))-exp(-beta(1)*tprDepths(ix+1:end)));
D_2 = (beta(2)/(beta(2)-mu))*(exp(-mu*tprDepths(ix+1:end))-exp(-beta(2)*tprDepths(ix+1:end)));
D_3 = (beta(3)/(beta(3)-mu))*(exp(-mu*tprDepths(ix+1:end))-exp(-beta(3)*tprDepths(ix+1:end)));

mx1 = [D_1 D_2 D_3]'*[D_1 D_2 D_3];
mx2 = [D_1 D_2 D_3]'*scaledTpr(ix+1:end,:);

W_ri = (mx1\mx2)';

%% compute normalization for kernel
kernelExtension = 512; % pixel
kernelResolution = 0.5; % mm
fRNorm = ppbkc_calcKernelNorm(kernelExtension,kernelResolution,primaryFluence);

%% computation of correction factors for output factor

FWHM  = params('fwhm_gauss'); % mm
sigma = FWHM/(sqrt(8*log(2))); % mm
sigmaVox = sigma/kernelResolution; % vox
iCenter = kernelExtension/2;

correctionFactors = ones(size(outputFactor,1),1);

[X,Y] = meshgrid(-kernelExtension/2+1:kernelExtension/2);
Z = sqrt(X.^2 + Y.^2)*0.5;
primflu = interp1(primaryFluence(:,1),primaryFluence(:,2),Z,'linear',0);

for i = 1:numel(outputFactor(:,1))
        
    if outputFactor(i,1) < 5.4 * sqrt(2)*sigma
        
        lowerLimit = round(iCenter - outputFactor(i,1)/2/kernelResolution + 1);
        upperLimit = floor(iCenter + outputFactor(i,1)/2/kernelResolution);

        fieldShape = 0*primflu;
        fieldShape(lowerLimit:upperLimit,lowerLimit:upperLimit) = 1;

         gaussFilter = 1/(sqrt(2*pi)*sigmaVox) .* exp( -X.^2 / (2*sigmaVox^2) ) ...
                    .* 1/(sqrt(2*pi)*sigmaVox) .* exp( -Y.^2 / (2*sigmaVox^2) );
               

        convRes = fftshift( ifft2(    fft2(fieldShape .*primflu,kernelExtension,kernelExtension) ...
                                   .* fft2(gaussFilter,kernelExtension,kernelExtension) ) );
        correctionFactors(i) = 1/convRes(iCenter-1,iCenter-1);
                     
    end
    
end

I don't quite understand this algorithm. I can't find any detailed formulas in Bortfeld et al.(1993) article. Can you help me get through this that I can understand the whole algorithm better. If I learn more, I have more capabilities to join the open source team for TPS algorithms. Thanks.

Tonio

[wahln]

I have not fully familiar with this code, since it was written by @markbangert and @mingersming, but here's my guesswork at the first part of the code:

photonPencilBeamKernelCalc/ppbkc.m

Lines 42 to 51 in c2296dd

	maxPos_fun = @(x) (log(mu)-log(x))/(mu-x);
	beta(1) = fmincon(@(x) (maxPos_fun(x) - meanMaxPos_mm).^2,1,[],[],[],[],0,1000);
	beta(2) = fmincon(@(x) (maxPos_fun(x) - (meanMaxPos_mm+1/mu)/2).^2,1,[],[],[],[],0,1000);
	beta(3) = fmincon(@(x) (maxPos_fun(x) - 1/mu).^2,1,[],[],[],[],0,1000);
	%% get weights for indivdual compoents of all TPRs
	D_1 = (beta(1)/(beta(1)-mu))*(exp(-mu*tprDepths(ix+1:end))-exp(-beta(1)*tprDepths(ix+1:end)));
	D_2 = (beta(2)/(beta(2)-mu))*(exp(-mu*tprDepths(ix+1:end))-exp(-beta(2)*tprDepths(ix+1:end)));
	D_3 = (beta(3)/(beta(3)-mu))*(exp(-mu*tprDepths(ix+1:end))-exp(-beta(3)*tprDepths(ix+1:end)));

Note that in the paper it say that the determination of the beta values is not critical. The fmincon performs a fit, and I think these fits are somewhat heuristic to get to similar magnitudes in the paper, but I could be wrong.

photonPencilBeamKernelCalc/ppbkc.m

Lines 53 to 56 in c2296dd

	mx1 = [D_1 D_2 D_3]'*[D_1 D_2 D_3];
	mx2 = [D_1 D_2 D_3]'*scaledTpr(ix+1:end,:);
	W_ri = (mx1\mx2)';

The weights are obtained not directly by a least squares fit, but by building a linear system. A matrix of all pairwise multiplications of D_i is built (mx1) together with a matrix containing the products of D_i with the measured TPRs. Solving this with Matlab’s backslash operator fits a set of weights for each TPR so it is closely reproduced by multiplication with the weights.

[Erhushenshou]

I still find it hard to understand why beta could be calculated through optimizing (log(mu)-log(x))/(mu-x) = dmax,dmax-1/mu,1/,u. And why is it that Di.T * D/(D.T*Tpr)=Wi, while the formula in the paper should be
. If it is Singular value decomposion in matrix, why is that? Can Mark's email be shared so that confusio could be explained. Great thanks.

[wahln]

The equation you showed, which is Eq. 18 in the paper, is handled later directly together with the kernel part from Eq. 19. https://github.com/e0404/photonPencilBeamKernelCalc/blob/master/ppbkc.m#L105-116
The derivative from (18) is numerically approximated from the linearly fitted W_ri for all field sizes. From this, the lateral kernels are then directly computed.

Regarding the determination of the beta, it is really purely heuristic. The maxPos_fun describes the maximum position of a depth dose component. The betas are chosen to create depth dose components with (1) a maximum close to the one averaged of all TPR measurments, (3) a deep maximum at 1/mu, and (2) a maximum in-between.
One could surely find other heuristics or find the "best" beta values / depth dose components, but there does not seem to be the need for it as long as you are choosing somewhat reasonable betas given the physical motivation behind the curves.

[Erhushenshou]

Thanks for reply. So what are the physical considerations and aspects of choosing to create depth dose components with (1) a maximum close to the one averaged of all TPR measurments, (3) a deep maximum at 1/mu, and (2) a maximum in-between?

[wahln]

For the exact values, there is none afaik. These values are chosen to "mimic" what isolation of primary and small/large field scatter would give according to the paper. Sorry that I can't give a more satisfying answer here. But note that the decomposition does not need a physical motivation, you could also choose a purely mathematical one, or, as in our case, a heuristic. If you want to choose different values, do it.