diff --git a/matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m b/matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m index 58ab20995..efd950bd5 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m +++ b/matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m @@ -246,14 +246,27 @@ function writeAllFiles(obj,ct,cst,stf,machine,w) end end - % Get alpha beta parameters from bioParam struct - if isfield(obj.bioParameters, 'tissuseAlphaX') - obj.bioParameters.AlphaX = obj.bioModel.tissueAlphaX(1); - obj.bioParameters.BetaX = obj.bioModel.tissueBetaX(1); + tmpAlphaX = []; + tmpBetaX = []; + for idx = 1:size(cst, 1) + if ~isempty(cst{idx,5}) && isfield(cst{idx,5}, 'alphaX') + tmpAlphaX = [tmpAlphaX cst{idx,5}.alphaX]; + tmpBetaX = [tmpBetaX cst{idx,5}.betaX]; + end end - if numel(obj.bioParameters.AlphaX)>1 - matRad_cfg.dispWarning('!!! Only a unique alpha/beta ratio supported at the moment. Found multiple, only the first one will be used !!!!'); + abX = [tmpAlphaX(:) tmpBetaX(:)]; + unique_abX = unique(abX, 'rows'); + obj.bioParameters.AlphaX = unique_abX(:, 1); + obj.bioParameters.BetaX = unique_abX(:, 2); + + % Get alpha beta parameters from bioParam struct + if isfield(obj.bioParameters, 'tissueAlphaX') + obj.bioParameters.AlphaX = obj.bioModel.tissueAlphaX; + obj.bioParameters.BetaX = obj.bioModel.tissueBetaX; end + %if numel(obj.bioParameters.AlphaX)>1 + % matRad_cfg.dispWarning('!!! Only a unique alpha/beta ratio supported at the moment. Found multiple, only the first one will be used !!!!'); + %end end if obj.scorer.LET @@ -508,10 +521,15 @@ function writeAllFiles(obj,ct,cst,stf,machine,w) end % Get photon parameters for RBExDose calculation - if this.calcBioDose + if this.calcBioDose || this.scorer.RBE this.scorer.RBE = true; - [dij.ax,dij.bx] = matRad_getPhotonLQMParameters(cst,dij.doseGrid.numOfVoxels,1,VdoseGrid); - dij.abx(dij.bx>0) = dij.ax(dij.bx>0)./dij.bx(dij.bx>0); + this.calcBioDose = true; + [dij.ax,dij.bx] = matRad_getPhotonLQMParameters(cst,dij.doseGrid.numOfVoxels,this.VdoseGrid); + dij.abx = zeros(size(dij.ax{1})); + ax = dij.ax{1}; + bx = dij.bx{1}; + dij.abx(bx>0) = ax(bx>0)./bx(bx>0); + dij.abx = {dij.abx}; end % save current directory to revert back to later @@ -1171,16 +1189,39 @@ function writeAllFiles(obj,ct,cst,stf,machine,w) end % Handle RBE-related quantities (not multiplied by sum(w)!) elseif ~isempty(strfind(lower(topasCubesTallies{j}),'alpha')) - modelName = strsplit(topasCubesTallies{j},'_'); - modelName = modelName{end}; + talliesFlags = strsplit(topasCubesTallies{j},'_'); + modelName = talliesFlags{end}; if isfield(topasCubes,[topasCubesTallies{j} '_beam' num2str(d)]) && iscell(topasCubes.([topasCubesTallies{j} '_beam' num2str(d)])) - dij.(['mAlphaDose_' modelName]){ctScen}(:,d) = reshape(topasCubes.([topasCubesTallies{j} '_beam',num2str(d)]){ctScen},[],1) .* dij.physicalDose{ctScen}(:,d); + + if contains(topasCubesTallies{j}, 'CellType') + ab_idx = str2num(talliesFlags{2}); + organAlpha = obj.bioParameters.AlphaX(ab_idx); + organBeta = obj.bioParameters.BetaX(ab_idx); + mask = find( (dij.ax{1} == organAlpha) & (dij.bx{1} == organBeta)); + topasCube_values = reshape(topasCubes.([topasCubesTallies{j} '_beam',num2str(d)]){ctScen},[],1); + topasCube_values = topasCube_values(mask,d); + dij.(['mAlphaDose_' modelName]){ctScen}(mask,d) = topasCube_values .* dij.physicalDose{ctScen}(mask,d); + else + + dij.(['mAlphaDose_' modelName]){ctScen}(:,d) = reshape(topasCubes.([topasCubesTallies{j} '_beam',num2str(d)]){ctScen},[],1) .* dij.physicalDose{ctScen}(:,d); + end end elseif ~isempty(strfind(lower(topasCubesTallies{j}),'beta')) modelName = strsplit(topasCubesTallies{j},'_'); modelName = modelName{end}; if isfield(topasCubes,[topasCubesTallies{j} '_beam' num2str(d)]) && iscell(topasCubes.([topasCubesTallies{j} '_beam' num2str(d)])) + if contains(topasCubesTallies{j}, 'CellType') + ab_idx = str2num(talliesFlags{2}); + organAlpha = obj.bioParameters.AlphaX(ab_idx); + organBeta = obj.bioParameters.BetaX(ab_idx); + mask = find( (dij.ax{1} == organAlpha) & (dij.bx{1} == organBeta)); + topasCube_values = reshape(topasCubes.([topasCubesTallies{j} '_beam',num2str(d)]){ctScen},[],1); + topasCube_values = topasCube_values(mask,d); + dij.(['mBetaDose_' modelName]){ctScen}(mask,d) = topasCube_values .* dij.physicalDose{ctScen}(mask,d); + else + dij.(['mSqrtBetaDose_' modelName]){ctScen}(:,d) = sqrt(reshape(topasCubes.([topasCubesTallies{j} '_beam',num2str(d)]){ctScen},[],1)) .* dij.physicalDose{ctScen}(:,d); + end end elseif ~isempty(strfind(topasCubesTallies{j},'LET')) if isfield(topasCubes,[topasCubesTallies{j} '_beam' num2str(d)]) && iscell(topasCubes.([topasCubesTallies{j} '_beam' num2str(d)])) @@ -1355,7 +1396,31 @@ function writeScorers(obj,fID,beamIx) % Read appropriate scorer from file and write to config file matRad_cfg.dispDebug('Reading RBE Scorer from %s\n',fname); scorerName = fileread(fname); - fprintf(fID,'\n%s\n\n',scorerName); + generalScorer = scorerName; + + if length(obj.bioParameters.AlphaX) ==1 + fprintf(fID,'\n%s\n\n',scorerName); + else + for idxCell = 1:length(obj.bioParameters.AlphaX) + scorerName = generalScorer; + insertText = ['_CellType_' num2str(idxCell)]; + scorerName = strrep(scorerName, ... + 'Alpha/', ['Alpha' insertText '/']); + scorerName = strrep(scorerName, ... + 'Beta/', ['Beta' insertText '/']); + scorerName = strrep(scorerName, ... + 'Sc/PrescribedDose Gy', ['Sc/PrescribedDose' insertText ' Gy']); + scorerName = strrep(scorerName, ... + 'Sc/CellLines', ['Sc/CellLines' insertText]); + scorerName = strrep(scorerName, ... + 'Sc/SimultaneousExposure', ['Sc/SimultaneousExposure' insertText]); + scorerName = strrep(scorerName, ... + 'Sim/ScoreLabel + "', ['Sim/ScoreLabel + "' insertText]); + + fprintf(fID,'\n%s\n\n',scorerName); + end + + end if obj.calc4DInterplay for PhaseNum = obj.MCparam.Phases{beamIx}' @@ -1402,28 +1467,60 @@ function writeScorers(obj,fID,beamIx) % Begin writing biological scorer components: cell lines switch obj.radiationMode case 'protons' - fprintf(fID,'\n### Biological Parameters ###\n'); - fprintf(fID,'sv:Sc/CellLines = 1 "CellLineGeneric"\n'); - fprintf(fID,'d:Sc/CellLineGeneric/Alphax = Sc/AlphaX /Gy\n'); - fprintf(fID,'d:Sc/CellLineGeneric/Betax = Sc/BetaX /Gy2\n'); - fprintf(fID,'d:Sc/CellLineGeneric/AlphaBetaRatiox = Sc/AlphaBetaX Gy\n\n'); + if length(obj.bioParameters.AlphaX) ==1 + fprintf(fID,'\n### Biological Parameters ###\n'); + fprintf(fID,'sv:Sc/CellLines = 1 "CellLineGeneric"\n'); + fprintf(fID,'d:Sc/CellLineGeneric/Alphax = Sc/AlphaX /Gy\n'); + fprintf(fID,'d:Sc/CellLineGeneric/Betax = Sc/BetaX /Gy2\n'); + fprintf(fID,'d:Sc/CellLineGeneric/AlphaBetaRatiox = Sc/AlphaBetaX Gy\n\n'); + else + for idxCell = 1:length(obj.bioParameters.AlphaX) + insertText = ['_CellType_' num2str(idxCell)]; + fprintf(fID,'\n### Biological Parameters ###\n'); + fprintf(fID, ['sv:Sc/CellLines' insertText ' = 1 "CellLineGeneric' insertText '"\n']); + fprintf(fID, ['d:Sc/CellLineGeneric' insertText '/Alphax = Sc/AlphaX' insertText ' /Gy\n']); + fprintf(fID, ['d:Sc/CellLineGeneric' insertText '/Betax = Sc/BetaX' insertText ' /Gy2\n']); + fprintf(fID, ['d:Sc/CellLineGeneric' insertText '/AlphaBetaRatiox = Sc/AlphaBetaX' insertText ' Gy\n\n']); + end + end case {'carbon','helium'} - fprintf(fID,'\n### Biological Parameters ###\n'); - fprintf(fID,'sv:Sc/CellLines = 1 "CellGeneric_abR2"\n'); - fprintf(fID,'d:Sc/CellGeneric_abR2/Alphax = Sc/AlphaX /Gy\n'); - fprintf(fID,'d:Sc/CellGeneric_abR2/Betax = Sc/BetaX /Gy2\n\n'); - % fprintf(fID,'d:Sc/CellGeneric_abR2/AlphaBetaRatiox = Sc/AlphaBetaX Gy\n'); + if length(obj.bioParameters.AlphaX) ==1 + fprintf(fID,'\n### Biological Parameters ###\n'); + fprintf(fID,'sv:Sc/CellLines = 1 "CellGeneric_abR2"\n'); + fprintf(fID,'d:Sc/CellGeneric_abR2/Alphax = Sc/AlphaX /Gy\n'); + fprintf(fID,'d:Sc/CellGeneric_abR2/Betax = Sc/BetaX /Gy2\n\n'); + % fprintf(fID,'d:Sc/CellGeneric_abR2/AlphaBetaRatiox = Sc/AlphaBetaX Gy\n'); + else + for idxCell = 1:length(obj.bioParameters.AlphaX) + insertText = ['_CellType_' num2str(idxCell)]; + fprintf(fID,'\n### Biological Parameters ###\n'); + fprintf(fID, ['sv:Sc/CellLines' insertText ' = 1 "CellLineGeneric_abR2' insertText '"\n']); + fprintf(fID, ['d:Sc/CellLineGeneric_abR2' insertText '/Alphax = Sc/AlphaX' inserText ' /Gy\n']); + fprintf(fID, ['d:Sc/CellLineGeneric_abR2' insertText '/Betax = Sc/BetaX' insertText ' /Gy2\n\n']); + end + end otherwise matRad_cfg.dispError([obj.radiationMode ' not implemented']); end % write biological scorer components: dose parameters matRad_cfg.dispDebug('Writing Biologial Scorer components.\n'); - fprintf(fID,'d:Sc/PrescribedDose = %.4f Gy\n',obj.bioParameters.PrescribedDose); - fprintf(fID,'b:Sc/SimultaneousExposure = %s\n',obj.bioParameters.SimultaneousExposure); - fprintf(fID,'d:Sc/AlphaX = %.4f /Gy\n',obj.bioParameters.AlphaX); - fprintf(fID,'d:Sc/BetaX = %.4f /Gy2\n',obj.bioParameters.BetaX); - fprintf(fID,'d:Sc/AlphaBetaX = %.4f Gy\n',obj.bioParameters.AlphaX/obj.bioParameters.BetaX); + if length(obj.bioParameters.AlphaX) ==1 + fprintf(fID,'d:Sc/PrescribedDose = %.4f Gy\n',obj.bioParameters.PrescribedDose); + fprintf(fID,'b:Sc/SimultaneousExposure = %s\n',obj.bioParameters.SimultaneousExposure); + fprintf(fID,'d:Sc/AlphaX = %.4f /Gy\n',obj.bioParameters.AlphaX); + fprintf(fID,'d:Sc/BetaX = %.4f /Gy2\n',obj.bioParameters.BetaX); + fprintf(fID,'d:Sc/AlphaBetaX = %.4f Gy\n',obj.bioParameters.AlphaX/obj.bioParameters.BetaX); + else + for idxCell = 1:length(obj.bioParameters.AlphaX) + insertText = ['_CellType_' num2str(idxCell)]; + fprintf(fID, ['d:Sc/PrescribedDose' insertText ' = %.4f Gy\n'],obj.bioParameters.PrescribedDose); + fprintf(fID, ['b:Sc/SimultaneousExposure' insertText ' = %s\n'],obj.bioParameters.SimultaneousExposure); + fprintf(fID, ['d:Sc/AlphaX' insertText ' = %.4f /Gy\n'],obj.bioParameters.AlphaX(idxCell)); + fprintf(fID, ['d:Sc/BetaX' insertText ' = %.4f /Gy2\n'],obj.bioParameters.BetaX(idxCell)); + fprintf(fID, ['d:Sc/AlphaBetaX' insertText ' = %.4f Gy\n\n'],obj.bioParameters.AlphaX(idxCell)/obj.bioParameters.BetaX(idxCell)); + end + end % Update MCparam.tallies with processed scorer for i = 1:length(obj.scorer.RBE_model) @@ -1451,9 +1548,23 @@ function writeScorers(obj,fID,beamIx) % Write subscorer to config files for s = 1:length(scorerNames) if strcmp(obj.radiationMode,'protons') - fprintf(fID,'s:Sc/%s%s/ReferencedSubScorer_LET = "ProtonLET"\n',scorerPrefix,scorerNames{s}); + if length(obj.bioParameters.AlphaX) ==1 + fprintf(fID,'s:Sc/%s%s/ReferencedSubScorer_LET = "ProtonLET"\n',scorerPrefix,scorerNames{s}); + else + for idxCell = 1:length(obj.bioParameters.AlphaX) + insertText = ['_CellType_' num2str(idxCell)]; + fprintf(fID,['s:Sc/%s%s' insertText '/ReferencedSubScorer_LET = "ProtonLET"\n'],scorerPrefix,scorerNames{s}); + end + end end + if length(obj.bioParameters.AlphaX) ==1 fprintf(fID,'s:Sc/%s%s/ReferencedSubScorer_Dose = "Tally_DoseToWater"\n',scorerPrefix,scorerNames{s}); + else + for idxCell = 1:length(obj.bioParameters.AlphaX) + insertText = ['_CellType_' num2str(idxCell)]; + fprintf(fID,['s:Sc/%s%s' insertText '/ReferencedSubScorer_Dose = "Tally_DoseToWater"\n'],scorerPrefix,scorerNames{s}); + end + end if obj.calc4DInterplay for PhaseNum = obj.MCparam.Phases{beamIx}' if strcmp(obj.radiationMode,'protons') diff --git a/matRad/planAnalysis/matRad_compareDose.m b/matRad/planAnalysis/matRad_compareDose.m index f79feb831..3adcd050b 100644 --- a/matRad/planAnalysis/matRad_compareDose.m +++ b/matRad/planAnalysis/matRad_compareDose.m @@ -148,6 +148,9 @@ % Calculate absolute difference cube and dose windows for plots differenceCube = cube1-cube2; doseDiffWindow = [-max(abs(differenceCube(:))) max(abs(differenceCube(:)))]; + if doseDiffWindow(1)==0 && doseDiffWindow(2)==0 + doseDiffWindow(2) = 1; + end %doseGammaWindow = [0 max(gammaCube(:))]; doseGammaWindow = [0 2]; %We choose 2 as maximum value since the gamma colormap has a sharp cut in the middle