Skip to content
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.

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

BED implementation for optimization #534

Open
wants to merge 16 commits into
base: dev
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 10 commits
Commits
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
38 changes: 30 additions & 8 deletions matRadGUI.m
Original file line number Diff line number Diff line change
Expand Up @@ -182,6 +182,10 @@
end
set(handles.popUpMachine,'String',handles.Machines);

%Add Bed as Biogical Optimization method
stringBioOpt = get(handles.popMenuBioOpt, 'String');
stringBioOpt{5} = 'BED';
set(handles.popMenuBioOpt, 'String', stringBioOpt);

vChar = get(handles.editGantryAngle,'String');
if strcmp(vChar(1,1),'0') && length(vChar)==6
Expand Down Expand Up @@ -716,10 +720,10 @@ function popupRadMode_Callback(hObject, eventdata, handles)
switch RadIdentifier
case 'photons'

set(handles.popMenuBioOpt,'Enable','off');
set(handles.popMenuBioOpt,'Enable','on');
ix = find(strcmp(contentPopUp,'none'));
set(handles.popMenuBioOpt,'Value',ix);
set(handles.btnSetTissue,'Enable','off');
set(handles.btnSetTissue,'Enable','on');

set(handles.btnRunSequencing,'Enable','on');
set(handles.btnRunDAO,'Enable','on');
Expand Down Expand Up @@ -1872,10 +1876,10 @@ function UpdateState(handles)
set(handles.btnSetTissue,'Enable','on');
elseif strcmp(pln.radiationMode,'protons')
set(handles.popMenuBioOpt,'Enable','on');
set(handles.btnSetTissue,'Enable','off');
set(handles.btnSetTissue,'Enable','on');
else
set(handles.popMenuBioOpt,'Enable','off');
set(handles.btnSetTissue,'Enable','off');
set(handles.popMenuBioOpt,'Enable','on');
set(handles.btnSetTissue,'Enable','on');
end

cMapControls = allchild(handles.uipanel_colormapOptions);
Expand Down Expand Up @@ -2027,16 +2031,18 @@ function setPln(handles)
set(handles.popupRadMode,'Value',find(strcmp(get(handles.popupRadMode,'String'),pln.radiationMode)));
set(handles.popUpMachine,'Value',find(strcmp(get(handles.popUpMachine,'String'),pln.machine)));


if ~strcmp(pln.propOpt.bioOptimization,'none')
set(handles.popMenuBioOpt,'Enable','on');
contentPopUp = get(handles.popMenuBioOpt,'String');
ix = find(strcmp(pln.propOpt.bioOptimization,contentPopUp));
set(handles.popMenuBioOpt,'Value',ix);
set(handles.btnSetTissue,'Enable','on');
else
set(handles.popMenuBioOpt,'Enable','off');
set(handles.btnSetTissue,'Enable','off');
set(handles.popMenuBioOpt,'Enable','on');
set(handles.btnSetTissue,'Enable','on');
end

%% enable sequencing and DAO button if radiation mode is set to photons
if strcmp(pln.radiationMode,'photons') && pln.propOpt.runSequencing
set(handles.btnRunSequencing,'Enable','on');
Expand Down Expand Up @@ -2175,12 +2181,19 @@ function getPlnFromGUI(handles)
contents = get(handles.popUpMachine,'String');
pln.machine = contents{get(handles.popUpMachine,'Value')};

%{
if (~strcmp(pln.radiationMode,'photons'))
contentBioOpt = get(handles.popMenuBioOpt,'String');
pln.propOpt.bioOptimization = contentBioOpt{get(handles.popMenuBioOpt,'Value'),:};
elseif isfield(pln, 'propOpt') && strcmp(pln.propOpt.bioOptimization, 'BED')
pln.propOpt.bioOptimization = 'BED';
else
pln.propOpt.bioOptimization = 'none';
end
%}

contentBioOpt = get(handles.popMenuBioOpt,'String');
pln.propOpt.bioOptimization = contentBioOpt{get(handles.popMenuBioOpt,'Value'),:};

pln.propOpt.runSequencing = logical(get(handles.btnRunSequencing,'Value'));
pln.propOpt.runDAO = logical(get(handles.btnRunDAO,'Value'));
Expand Down Expand Up @@ -2322,7 +2335,12 @@ function btnDVH_Callback(~, ~, handles)
%check if one of the default fields is selected
if sum(strcmp(SelectedCube,{'physicalDose','effect','RBE,','RBExDose','alpha','beta'})) > 0
resultGUI_SelectedCube.physicalDose = resultGUI.physicalDose;
resultGUI_SelectedCube.RBExDose = resultGUI.RBExDose;
if isfield(resultGUI,'RBExDose')
resultGUI_SelectedCube.RBExDose = resultGUI.RBExDose;
end
elseif strcmp(SelectedCube,'BED')
resultGUI_SelectedCube.physicalDose = resultGUI.physicalDose;
resultGUI_SelectedCube.BED = resultGUI.BED;
else
Idx = find(SelectedCube == '_');
SelectedSuffix = SelectedCube(Idx(1):end);
Expand Down Expand Up @@ -2876,6 +2894,10 @@ function SaveResultToGUI(~, ~, ~)
if isfield(resultGUI,'RBExDose')
resultGUI.(['RBExDose' Suffix]) = resultGUI.RBExDose;
end
if isfield(resultGUI,'BED')
resultGUI.(['BED' Suffix]) = resultGUI.BED;
end


if strcmp(pln.radiationMode,'carbon') == 1
if isfield(resultGUI,'effect')
Expand Down
35 changes: 33 additions & 2 deletions matRad_calcCubes.m
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
function resultGUI = matRad_calcCubes(w,dij,scenNum)
function resultGUI = matRad_calcCubes(w,dij,pln,scenNum)
% matRad computation of all cubes for the resultGUI struct
% which is used as result container and for visualization in matRad's GUI
%
Expand Down Expand Up @@ -30,7 +30,7 @@
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 3
if ~exist('sceneNum', 'var')
scenNum = 1;
end

Expand Down Expand Up @@ -100,6 +100,37 @@
end
end

% CACLCULATE BED (biological effective dose) BED
if exist('pln','var')
% When depth Dependent alpha beta values are calculated in dij calculation
alphax = reshape(dij.ax, dij.doseGrid.dimensions);
ix = ~(alphax == 0);
if isfield(dij,'mAlphaDose') && isfield(dij,'mSqrtBetaDose')
for i = 1:length(beamInfo)
% photon equivaluent BED = n * effect / alphax
resultGUI.(['BED', beamInfo(i).suffix]) = zeros(dij.doseGrid.dimensions);
resultGUI.(['BED', beamInfo(i).suffix])(ix) = full(pln.numOfFractions.*resultGUI.(['effect', beamInfo(i).suffix])(ix) ./alphax(ix));
resultGUI.(['BED', beamInfo(i).suffix]) = reshape(resultGUI.(['BED', beamInfo(i).suffix]), dij.doseGrid.dimensions);
end

else
% Get Alpha and Beta Values form dij.ax and dij.bx
alphax = reshape(dij.ax, dij.doseGrid.dimensions);
betax = reshape(dij.bx, dij.doseGrid.dimensions);
for i = 1:length(beamInfo)
ix = ~isnan(dij.ax./dij.bx);
if isfield(resultGUI, 'RBExDose')
Dose = resultGUI.(['RBExDose', beamInfo(i).suffix]);
else
Dose = resultGUI.(['physicalDose', beamInfo(i).suffix]);
end
effect = alphax.* Dose + betax.*Dose.^2;
resultGUI.(['BED', beamInfo(i).suffix]) = zeros(dij.doseGrid.dimensions);
resultGUI.(['BED', beamInfo(i).suffix])(ix) = full(pln.numOfFractions.*effect(ix)./alphax(ix));
amitantony marked this conversation as resolved.
Show resolved Hide resolved
resultGUI.(['BED', beamInfo(i).suffix]) = reshape(resultGUI.(['BED', beamInfo(i).suffix]), dij.doseGrid.dimensions);
end
end
end
% group similar fields together
resultGUI = orderfields(resultGUI);

Expand Down
8 changes: 4 additions & 4 deletions matRad_calcDoseFillDij.m
Original file line number Diff line number Diff line change
Expand Up @@ -12,8 +12,8 @@
dij.mLETDose{1}(:,i) = dij.mLETDose{1}(:,i) + stf(i).ray(j).weight(k) * letDoseTmpContainer{1,1};
end

if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')) ...
&& strcmp(pln.radiationMode,'carbon')
if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD') ...
||isequal(pln.propOpt.bioOptimization,'BED') ) && strcmp(pln.radiationMode,'carbon')
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should this not be LEMIV_BED? for consistency


% score alpha and beta matrices
dij.mAlphaDose{1}(:,i) = dij.mAlphaDose{1}(:,i) + stf(i).ray(j).weight(k) * alphaDoseTmpContainer{1,1};
Expand All @@ -33,8 +33,8 @@
dij.mLETDose{1}(:,(ceil(counter/numOfBixelsContainer)-1)*numOfBixelsContainer+1:counter) = [letDoseTmpContainer{1:mod(counter-1,numOfBixelsContainer)+1,1}];
end

if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')) ...
&& strcmp(pln.radiationMode,'carbon')
if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD') ...
||isequal(pln.propOpt.bioOptimization,'BED') )&& strcmp(pln.radiationMode,'carbon')
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

see the previous comment


dij.mAlphaDose{1}(:,(ceil(counter/numOfBixelsContainer)-1)*numOfBixelsContainer+1:counter) = [alphaDoseTmpContainer{1:mod(counter-1,numOfBixelsContainer)+1,1}];
dij.mSqrtBetaDose{1}(:,(ceil(counter/numOfBixelsContainer)-1)*numOfBixelsContainer+1:counter) = [betaDoseTmpContainer{1:mod(counter-1,numOfBixelsContainer)+1,1}];
Expand Down
27 changes: 15 additions & 12 deletions matRad_calcParticleDose.m
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,8 @@
% helper function for energy selection
round2 = @(a,b)round(a*10^b)/10^b;

if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')) ...
if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD') ...
||isequal(pln.propOpt.bioOptimization,'BED')) ...
&& strcmp(pln.radiationMode,'carbon')

alphaDoseTmpContainer = cell(numOfBixelsContainer,dij.numOfScenarios);
Expand All @@ -56,7 +57,8 @@
dij.mSqrtBetaDose{i} = spalloc(dij.doseGrid.numOfVoxels,numOfColumnsDij,1);
end

elseif isequal(pln.propOpt.bioOptimization,'const_RBExD') && strcmp(pln.radiationMode,'protons')
elseif (isequal(pln.propOpt.bioOptimization,'const_RBExD') ||isequal(pln.propOpt.bioOptimization,'BED') )...
&& strcmp(pln.radiationMode,'protons')
dij.RBE = 1.1;
matRad_cfg.dispInfo('matRad: Using a constant RBE of %g\n',dij.RBE);
end
Expand All @@ -76,12 +78,8 @@
end

% generates tissue class matrix for biological optimization
if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')) ...
&& strcmp(pln.radiationMode,'carbon')

if isfield(machine.data,'alphaX') && isfield(machine.data,'betaX')

matRad_cfg.dispInfo('matRad: loading biological base data... ');

matRad_cfg.dispInfo('matRad: loading biological base data... ');
vTissueIndex = zeros(size(VdoseGrid,1),1);
dij.ax = zeros(dij.doseGrid.numOfVoxels,1);
dij.bx = zeros(dij.doseGrid.numOfVoxels,1);
Expand All @@ -93,6 +91,11 @@
dij.doseGrid.x,dij.doseGrid.y,dij.doseGrid.z);
% retrieve photon LQM parameter for the current dose grid voxels
[dij.ax,dij.bx] = matRad_getPhotonLQMParameters(cst,dij.doseGrid.numOfVoxels,1,VdoseGrid);

if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')...
||isequal(pln.propOpt.bioOptimization,'BED') ) && strcmp(pln.radiationMode,'carbon')

if isfield(machine.data,'alphaX') && isfield(machine.data,'betaX')

for i = 1:size(cst,1)

Expand Down Expand Up @@ -200,8 +203,8 @@
end

% just use tissue classes of voxels found by ray tracer
if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')) ...
&& strcmp(pln.radiationMode,'carbon')
if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')...
||isequal(pln.propOpt.bioOptimization,'BED') ) && strcmp(pln.radiationMode,'carbon')
vTissueIndex_j = vTissueIndex(ix,:);
end

Expand Down Expand Up @@ -362,8 +365,8 @@



if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')) ...
&& strcmp(pln.radiationMode,'carbon')
if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD') ...
||isequal(pln.propOpt.bioOptimization,'BED') ) && strcmp(pln.radiationMode,'carbon')
% calculate alpha and beta values for bixel k on ray j of
[bixelAlpha, bixelBeta] = matRad_calcLQParameter(...
currRadDepths,...
Expand Down
11 changes: 11 additions & 0 deletions matRad_calcParticleDoseMC.m
Original file line number Diff line number Diff line change
Expand Up @@ -288,6 +288,17 @@
dij.physicalDose{1} = dij.physicalDose{1}(:,MCsquareOrder);
end

% calculate Alpha and Beta Values for BED projection
dij.ax = zeros(dij.doseGrid.numOfVoxels,1);
dij.bx = zeros(dij.doseGrid.numOfVoxels,1);
cst = matRad_setOverlapPriorities(cst);
% resizing cst to dose cube resolution
cst = matRad_resizeCstToGrid(cst,dij.ctGrid.x,dij.ctGrid.y,dij.ctGrid.z,...
dij.doseGrid.x,dij.doseGrid.y,dij.doseGrid.z);
% retrieve photon LQM parameter for the current dose grid voxels
[dij.ax,dij.bx] = matRad_getPhotonLQMParameters(cst,dij.doseGrid.numOfVoxels,1,VdoseGrid);


matRad_cfg.dispInfo('matRad: done!\n');

try
Expand Down
13 changes: 13 additions & 0 deletions matRad_calcPhotonDose.m
Original file line number Diff line number Diff line change
Expand Up @@ -314,6 +314,19 @@
end
end

%% calculate Alpha and Beta Values for BED projection

dij.ax = zeros(dij.doseGrid.numOfVoxels,1);
dij.bx = zeros(dij.doseGrid.numOfVoxels,1);

cst = matRad_setOverlapPriorities(cst);

% resizing cst to dose cube resolution
cst = matRad_resizeCstToGrid(cst,dij.ctGrid.x,dij.ctGrid.y,dij.ctGrid.z,...
dij.doseGrid.x,dij.doseGrid.y,dij.doseGrid.z);
% retrieve photon LQM parameter for the current dose grid voxels
[dij.ax,dij.bx] = matRad_getPhotonLQMParameters(cst,dij.doseGrid.numOfVoxels,1,VdoseGrid);

%Close Waitbar
if ishandle(figureWait)
delete(figureWait);
Expand Down
11 changes: 11 additions & 0 deletions matRad_calcPhotonDoseMC.m
Original file line number Diff line number Diff line change
Expand Up @@ -328,6 +328,17 @@
end
end

% calculate Alpha and Beta Values for BED projection
dij.ax = zeros(dij.doseGrid.numOfVoxels,1);
dij.bx = zeros(dij.doseGrid.numOfVoxels,1);
cst = matRad_setOverlapPriorities(cst);
% resizing cst to dose cube resolution
cst = matRad_resizeCstToGrid(cst,dij.ctGrid.x,dij.ctGrid.y,dij.ctGrid.z,...
dij.doseGrid.x,dij.doseGrid.y,dij.doseGrid.z);
% retrieve photon LQM parameter for the current dose grid voxels
[dij.ax,dij.bx] = matRad_getPhotonLQMParameters(cst,dij.doseGrid.numOfVoxels,1,VdoseGrid);


matRad_cfg.dispInfo('matRad: MC photon dose calculation done!\n');
matRad_cfg.dispInfo(evalc('toc'));

Expand Down
24 changes: 21 additions & 3 deletions matRad_fluenceOptimization.m
Original file line number Diff line number Diff line change
Expand Up @@ -66,6 +66,7 @@
end

obj = obj.setDoseParameters(obj.getDoseParameters()/pln.numOfFractions);
obj.numOfFractions = pln.numOfFractions;

cst{i,6}{j} = obj;
end
Expand Down Expand Up @@ -168,7 +169,22 @@
4*cst{ixTarget,5}.betaX.*CurrEffectTarget)./(2*cst{ixTarget,5}.betaX*(dij.physicalDose{1}(V,:)*wOnes)));
wInit = ((doseTarget)/(TolEstBio*maxCurrRBE*max(dij.physicalDose{1}(V,:)*wOnes)))* wOnes;
end

elseif strcmp(pln.propOpt.bioOptimization, 'BED')
if isfield(dij, 'mAlphaDose') && isfield(dij, 'mSqrtBetaDose')
abr = cst{ixTarget,5}.alphaX./cst{ixTarget,5}.betaX;
meanBED = mean(pln.numOfFractions.*(dij.mAlphaDose{1}(V,:)*wOnes + (dij.mSqrtBetaDose{1}(V,:)*wOnes).^2)./cst{ixTarget,5}.alphaX);
BEDTarget = pln.numOfFractions.*doseTarget.*(1 + doseTarget./abr);
elseif isfield(dij, 'RBE')
abr = cst{ixTarget,5}.alphaX./cst{ixTarget,5}.betaX;
meanBED = mean(pln.numOfFractions.*dij.RBE.*dij.physicalDose{1}(V,:)*wOnes.*(1+dij.RBE.*dij.physicalDose{1}(V,:)*wOnes./abr));
BEDTarget = pln.numOfFractions.*dij.RBE.*doseTarget.*(1 + dij.RBE.*doseTarget./abr);
else
abr = cst{ixTarget,5}.alphaX./cst{ixTarget,5}.betaX;
meanBED = mean(pln.numOfFractions.*dij.physicalDose{1}(V,:)*wOnes.*(1+dij.physicalDose{1}(V,:)*wOnes./abr));
BEDTarget = pln.numOfFractions.*doseTarget.*(1 + doseTarget./abr);
end
bixelWeight = BEDTarget/meanBED;
wInit = wOnes * bixelWeight;
else
bixelWeight = (doseTarget)/(mean(dij.physicalDose{1}(V,:)*wOnes));
wInit = wOnes * bixelWeight;
Expand All @@ -190,14 +206,16 @@
backProjection = matRad_ConstantRBEProjection;
case 'LEMIV_RBExD'
backProjection = matRad_VariableRBEProjection;
case 'BED'
backProjection = matRad_BEDProjection;
case 'none'
backProjection = matRad_DoseProjection;
otherwise
warning(['Did not recognize bioloigcal setting ''' pln.probOpt.bioOptimization '''!\nUsing physical dose optimization!']);
backProjection = matRad_DoseProjection;
end


backProjection.numOfFractions = pln.numOfFractions;
%backProjection = matRad_DoseProjection();

optiProb = matRad_OptimizationProblem(backProjection);
Expand Down Expand Up @@ -225,7 +243,7 @@
wOpt = optimizer.wResult;
info = optimizer.resultInfo;

resultGUI = matRad_calcCubes(wOpt,dij);
resultGUI = matRad_calcCubes(wOpt,dij,pln);
resultGUI.wUnsequenced = wOpt;
resultGUI.usedOptimizer = optimizer;
resultGUI.info = info;
Expand Down
2 changes: 2 additions & 0 deletions matRad_indicatorWrapper.m
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,8 @@

if isfield(resultGUI,'RBExDose')
doseCube = resultGUI.RBExDose;
elseif isfield(resultGUI,'BED')
doseCube = resultGUI.BED;
else
doseCube = resultGUI.physicalDose;
end
Expand Down
2 changes: 2 additions & 0 deletions matRad_showDVH.m
Original file line number Diff line number Diff line change
Expand Up @@ -96,6 +96,8 @@ function matRad_showDVH(dvh,cst,pln,lineStyleIndicator)
if exist('pln','var') && ~isempty(pln)
if strcmp(pln.propOpt.bioOptimization,'none')
xlabel('Dose [Gy]','FontSize',fontSizeValue);
elseif strcmp(pln.propOpt.bioOptimization,'BED')
xlabel('Dose Gy]','FontSize',fontSizeValue);
else
xlabel('RBE x Dose [Gy(RBE)]','FontSize',fontSizeValue);
end
Expand Down
1 change: 1 addition & 0 deletions optimization/+DoseConstraints/matRad_MinMaxDVH.m
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@
voxelScalingRatio = 1;
referenceScalingVal = 0.01;
parameters = {30,0,100};
numOfFractions = NaN;
end

methods
Expand Down
1 change: 1 addition & 0 deletions optimization/+DoseConstraints/matRad_MinMaxDose.m
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@
properties
parameters = {0,30,1};
epsilon = 1e-3; %slack parameter for the logistic approximation
numOfFractions = NaN;
end

methods
Expand Down
Loading