Skip to content

Extension to other dose distributions #4

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
fgesualdi wants to merge 21 commits into
base: master
Choose a base branch
from
Open

Extension to other dose distributions #4

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
21 commits
Select commit Hold shift + click to select a range
a66362a
Calculate the cdf of the dose for several distributions
fgesualdi Apr 1, 2022
7b10e3c
Compute different Copulas
fgesualdi Apr 1, 2022
a5d61f3
Add the computation of the moments of the probabilistic DVH using a l…
fgesualdi Apr 4, 2022
af310e1
Fix bug
fgesualdi Apr 4, 2022
1f12904
Return paremters k and theta
fgesualdi Apr 4, 2022
5b85e39
Change name
fgesualdi Apr 4, 2022
7bcf0bd
Add function explanation
fgesualdi Apr 4, 2022
4ca8fac
No change
fgesualdi Apr 4, 2022
a815484
Write as a module to be reachable from any other module
fgesualdi Apr 4, 2022
eb1009b
Change names of transform calls
fgesualdi Apr 4, 2022
a48d85e
Get support of a shifted beta distribution of the dose in a voxel
fgesualdi Apr 4, 2022
b7fdfaf
Get the mean and std of a beta distribution from that of a shifted beta
fgesualdi Apr 4, 2022
eae6abc
Get random multivariate gaussian samples
fgesualdi Apr 4, 2022
30c2044
Change names of transform calls
fgesualdi Apr 4, 2022
b95090b
Modularize
fgesualdi Apr 4, 2022
7c2f528
Modularize the sample validation and add optional plots
fgesualdi Apr 4, 2022
d6d8487
Modularize and extend functionality for other dose distribution assum…
fgesualdi Apr 4, 2022
5b4daea
Modularize
fgesualdi Apr 4, 2022
c2cdb9b
Modularize and extend functionality for other dose distribution assum…
fgesualdi Apr 4, 2022
21f14f4
Correct plotting
fgesualdi Apr 4, 2022
a6f3696
Fix bug
fgesualdi Apr 4, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
203 changes: 197 additions & 6 deletions cProbOpt/apm_DVHprob.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
function [expDvh,stdDvh,covDvh] = apm_DVHprob(expDose,covDose,nBins,dMax,method)
function [expDvh,stdDvh,covDvh] = apm_DVHprob(expDose,covDose,nBins,dMax,method,copula_model,copula_marginals,dmin_shiftedbeta,dmax_shiftedbeta)

if nargin < 3 || isempty(nBins)
nBins = 100;
Expand All @@ -12,6 +12,25 @@
method = 'int_gauss';
end

% Define parameters of the Copula
if nargin < 6
copula_model = 'gauss';
end

if nargin < 7
copula_marginals = cell(1,numel(expDose));
copula_marginals(:) = {'gauss'};
end

% Define support of the shifted beta
if nargin < 8
dmin_shiftedbeta= zeros(1,numel(expDose));
end

if nargin < 9
dmax_shiftedbeta = ones(1, numel(expDose));
end

dAxis = linspace(0,dMax,nBins);
expDvh(1,:) = dAxis;
stdDvh(1,:) = dAxis;
Expand All @@ -35,17 +54,28 @@
reducedVoxNum = numel(expDose);

switch method
case 'copula'
[expDvh(2,:),stdDvh(2,:)] = probDVHCopulaExpStd(expDose,covDose,dAxis,copula_model,copula_marginals,dmin_shiftedbeta,dmax_shiftedbeta);
if nargout > 2
covDvh = calcCovDvhCopula(expDose,covDose,dAxis,expDvh(2,:),copula_model,copula_marginals,dmin_shiftedbeta,dmax_shiftedbeta);
end
case 'int_gauss'
[expDvh(2,:),stdDvh(2,:)] = probDVHGaussExpStd(expDose,covDose,dAxis);
if nargout > 2
covDvh = calcCovDvhGauss(expDose,covDose,dAxis,expDvh(2,:));
end

case 'int_loggauss'
[expDvh(2,:),stdDvh(2,:)] = probDVHLnGaussExpStd(expDose,covDose,dAxis);
if nargout > 2
covDvh = calcCovDvhLnGauss(expDose,covDose,dAxis,expDvh(2,:));
end

case 'int_gauss_cuda'
disp('Setting up kernel');
tic




%single precision
%{
cuKernel = parallel.gpu.CUDAKernel('C:\dev\APM\CUDA\QI\DVH\x64\Release\dvh.ptx','float* , float*, const float*, const float*, const float*, int, int','dvhRawMoments');
Expand Down Expand Up @@ -145,6 +175,85 @@

end

function [expDvh,stdDvh] = probDVHCopulaExpStd(expDose,covDose,dAxis,copula_model,copula_marginals,dmin_shiftedbeta,dmax_shiftedbeta)
nBins = numel(dAxis);
biCum = zeros(1,nBins);
diagCum = zeros(1,nBins);

for i = 1:numel(expDose)
mean_i = expDose(i);
sig_i = sqrt(covDose(i,i));
if ~isreal(sig_i) || isnan(sig_i)
continue;
end
P1 = apm_calcCdfDose(dAxis, mean_i*ones(numel(1,nBins)), sig_i*ones(numel(1,nBins)), copula_marginals(i), 1,dmin_shiftedbeta(i),dmax_shiftedbeta(i));


if any(isnan(P1))
warning('CDF evaluates to NaN');
P1(isnan(P1)) = 0;
end
diagCum = diagCum + P1;

for l = i+1:numel(expDose)
mean_l = expDose(l);
sig_l = sqrt(covDose(l,l));
if ~isreal(sig_l) || isnan(sig_l)
continue;
end
r = (covDose(i,l) / (sig_i*sig_l)) * ones(1,nBins);
u_i = apm_calcCdfDose(dAxis, mean_i, sig_i, copula_marginals(i),0,dmin_shiftedbeta(i),dmax_shiftedbeta(i));
u_l = apm_calcCdfDose(dAxis, mean_l, sig_l, copula_marginals(l),0,dmin_shiftedbeta(l),dmax_shiftedbeta(l));
P2D = apm_calcCopula(u_i, u_l, r, copula_model) + 1 - u_i - u_l;

if any(isnan(P2D))
warning('CDF (2D) evaluates to NaN');
P2D(isnan(P2D)) = 0;
end

biCum = biCum + P2D;
end
matRad_progress(i,numel(expDose));
end

% biCum has only the terms from l=i+1 to numel(expDose),
% then we need to add the terms l=i (diagCum) and l from 0 to i-1 (biCum).
biCum = diagCum + 2*biCum;

expDvh = diagCum ./ numel(expDose);

stdDvh = sqrt( 1/numel(expDose)^2 * (biCum) - expDvh.^2);
end

function covDvh = calcCovDvhCopula(expDose,covDose,dAxis,expDvh,copula_model,copula_marginals,dmin_shiftedbeta,dmax_shiftedbeta)
nBins = numel(dAxis);
covDvh = zeros(nBins);

for p = 1:nBins
for q = 1:nBins
s = 0;
for i = 1:numel(expDose)
mean_i = expDose(i);
sig_i = sqrt(covDose(i,i));
for l = 1:numel(expDose)
mean_l = expDose(l);
sig_l = sqrt(covDose(l,l));
r = covDose(i,l) / (sig_i*sig_l);
u_i = apm_calcCdfDose(dAxis(p), mean_i, sig_i, copula_marginals(i),0,dmin_shiftedbeta(i),dmax_shiftedbeta(i));
u_l = apm_calcCdfDose(dAxis(q), mean_l, sig_l, copula_marginals(l),0,dmin_shiftedbeta(l),dmax_shiftedbeta(l));
s = s + apm_calcCopula(u_i, u_l, r, copula_model) + 1 - u_i - u_l;
end
end

covDvh(p,q) = s;
end
matRad_progress(p,nBins);
end
covDvh = 1/numel(expDose)^2 * covDvh;
covDvh = covDvh - transpose(expDvh)*expDvh;
end


function [expDvh,stdDvh] = probDVHGaussExpStd(expDose,covDose,dAxis)
nBins = numel(dAxis);
biCum = zeros(1,nBins);
Expand Down Expand Up @@ -175,7 +284,6 @@

r = covDose(i,l) / (sig_i*sig_l);
P2D = arrayfun(@(dParam) bvn((dParam - expDose(i)) / sig_i,Inf,(dParam - expDose(l)) / sig_l,Inf,r),dAxis);
%P2D = arrayfun(@(dParam) mexBVNcdf(-[dParam dParam],-[expDose(i) expDose(l)],[covDose(i,i) covDose(i,l); covDose(l,i) covDose(l,l)]),dAxis);

if any(isnan(P2D))
warn('CDF (2D) evaluates to NaN');
Expand All @@ -186,11 +294,14 @@
end
matRad_progress(i,numel(expDose));
end

% biCum has only the terms from l=i+1 to numel(expDose),
% then we need to add the terms l=i (diagCum) and l from 0 to i-1 (biCum).
biCum = diagCum + 2*biCum;

expDvh = diagCum ./ numel(expDose);
stdDvh = 1/numel(expDose)^2 * (diagCum + 2*biCum);

stdDvh = sqrt(stdDvh - expDvh.^2);
stdDvh = sqrt( 1/numel(expDose)^2 * (biCum) - expDvh.^2);
end

function covDvh = calcCovDvhGauss(expDose,covDose,dAxis,expDvh)
Expand All @@ -216,3 +327,83 @@
covDvh = 1/numel(expDose)^2 * covDvh;
covDvh = covDvh - transpose(expDvh)*expDvh;
end


function [expDvh,stdDvh] = probDVHLnGaussExpStd(expDose,covDose,dAxis)
nBins = numel(dAxis);
biCum = zeros(1,nBins);
diagCum = zeros(1,nBins);

%Compute parameters of the lognormal
[lnExpDose, lnCovDose] = apm_transformMeanCovarianceToLogNormalParameters(expDose,covDose);

for i = 1:numel(lnExpDose) %for each voxel
lnSig_i = sqrt(lnCovDose(i,i));

if ~isreal(lnSig_i) || isnan(lnSig_i)
continue;
end

P1 = 1 - apm_normcdf(log(dAxis),lnExpDose(i)*ones(numel(1,nBins)),lnSig_i*ones(numel(1,nBins)));

if any(isnan(P1))
warning('CDF evaluates to NaN');
P1(isnan(P1)) = 0;
end
diagCum = diagCum + P1;

for l = i+1:numel(lnExpDose)
lnSig_l = sqrt(lnCovDose(l,l));

if ~isreal(lnSig_l) || isnan(lnSig_l)
continue;
end

r = lnCovDose(i,l) / (lnSig_i*lnSig_l);
P2D = arrayfun(@(dParam) bvn((log(dParam) - lnExpDose(i)) / lnSig_i,Inf,(log(dParam) - lnExpDose(l)) / lnSig_l,Inf,r), dAxis);

if any(isnan(P2D))
warn('CDF (2D) evaluates to NaN');
P2D(isnan(P2D)) = 0;
end

biCum = biCum + P2D;
end
matRad_progress(i,numel(expDose));
end

% biCum has only the terms from l=i+1 to numel(expDose),
% then we need to add the terms l=i (diagCum) and l from 0 to i-1 (biCum).
biCum = diagCum + 2*biCum;

expDvh = diagCum ./ numel(lnExpDose);

stdDvh = sqrt( 1/numel(lnExpDose)^2 * (biCum) - expDvh.^2);
end

function covDvh = calcCovDvhLnGauss(expDose,covDose,dAxis,expDvh)
nBins = numel(dAxis);
covDvh = zeros(nBins);

%Compute parameters of the lognormal
[lnExpDose, lnCovDose] = apm_transformMeanCovarianceToLogNormalParameters(expDose,covDose);

for p = 1:nBins
for q = 1:nBins
s = 0;
for i = 1:numel(lnExpDose)
lnSig_i = sqrt(lnCovDose(i,i));
for l = 1:numel(lnExpDose)
lnSig_l = sqrt(lnCovDose(l,l));
r = lnCovDose(i,l) / (lnSig_i*lnSig_l);
s = s + bvn((log(dAxis(p)) - lnExpDose(i)) / lnSig_i,Inf,(log(dAxis(q)) - lnExpDose(l)) / lnSig_l,inferiorto,r);
end
end

covDvh(p,q) = s;
end
matRad_progress(p,nBins);
end
covDvh = 1/numel(lnExpDose)^2 * covDvh;
covDvh = covDvh - transpose(expDvh)*expDvh;
end
Loading