Merge dev_VARBE_MC_Update into dev_VARBE_MC_experimental

commit 3610bf5 Author: Noa Homolka <n.homolka@dkfz.de> Date: Wed Oct 12 13:57:57 2022 +0200 Update TopasConfig - added output of numOfBixels that did not meet the minimum numOfParticles required - removed unnecessary statement regarding alpha beta parameters - added comment to externalCalculation commit 834bb1f Merge: 19e22df 1b13404 Author: Niklas Wahl <n.wahl@dkfz.de> Date: Wed Oct 5 17:21:33 2022 +0200 Merge branch 'master' into pr/536 # Conflicts: # README.md commit 19e22df Author: Niklas Wahl <n.wahl@dkfz.de> Date: Wed Oct 5 16:36:42 2022 +0200 reset to old bixelwidth-dependent dij sampling commit 1ee132d Author: Noa Homolka <n.homolka@dkfz.de> Date: Tue Oct 4 20:31:30 2022 +0200 Merge dev_varRBErobOpt_experimental Important branch Bug Fix: - fixed 2 major bugs in matRad_calcCubes! - fixed bug in TopasConfig, where the empty preallocated dij field was used to calculate tallies instead of the actual dose - made beamSetup more concise with shorter spaces between numbers - added fix for selection of focusIndex for generated phaseSpace data - changed handling of 0 fields - use isprop instead of isfield now that the class system is used - removed several statements that are no longer used or required commit 1b13404 Author: Niklas Wahl <n.wahl@dkfz.de> Date: Sat Sep 24 00:39:04 2022 +0200 Update tests.yml commit d24dc5f Author: Niklas Wahl <n.wahl@dkfz.de> Date: Sat Sep 24 00:14:14 2022 +0200 Update README.md commit 4a3e73a Author: Noa Homolka <n.homolka@dkfz.de> Date: Thu Sep 22 21:35:00 2022 +0200 Update to emittanceBaseData energyspread handling commit 87bf5fe Author: Niklas Wahl <n.wahl@dkfz.de> Date: Fri Sep 23 15:25:24 2022 +0200 Fix indent commit 37d3c58 Author: Noa Homolka <n.homolka@dkfz.de> Date: Thu Sep 22 20:53:11 2022 +0200 Update to matRad_calcCubes.m - removed section that added "non-processed MC tallies", not needed anymore - removed 2 lines from MCsquare calculation that still referenced the outdated MC_tallies system commit dc9552e Author: Noa Homolka <n.homolka@dkfz.de> Date: Thu Sep 22 13:45:13 2022 +0200 MC dij update calcDoseDirect - renamed topasCubes variable into dij to make it clearer (used to be a relict from deprecated workflow) - Brought dij generation of MCsquare and TOPAS in line for calcDoseDirect (all beam/ray/bixel fields are 1) commit be032cc Author: Noa Homolka <n.homolka@dkfz.de> Date: Thu Sep 22 13:42:04 2022 +0200 bug fix for recent resampling update - fixed resampleCTtoGrid function being called incorrectly - fixed input commit 45b6517 Author: Noa Homolka <n.homolka@dkfz.de> Date: Wed Sep 21 16:44:07 2022 +0200 Merge dev_varRBErobOpt_experimental -> dev_varRBErobOpt_Update Update PR: - removed TOPAS scripts - added compatibility patch - Added patch for fluence Optimization and input of precalculated weigths - Added input of LET into calcLQParameter function - Moved resampling function from TOPASconfig in general tools - made externalCalculation parameter clearer to use and saved folder is output to the console - commit 942ea44 Merge: 75201fa f69a007 Author: Niklas Wahl <n.wahl@dkfz.de> Date: Sat Sep 17 01:33:39 2022 +0200 Merge pull request e0404#579 from wahln/master Create CITATION.cff commit f69a007 Author: Niklas Wahl <n.wahl@dkfz.de> Date: Sat Sep 17 01:28:31 2022 +0200 Create CITATION.cff commit 75201fa Merge: 3d29b76 e469b8f Author: Niklas Wahl <n.wahl@dkfz.de> Date: Mon Aug 22 15:37:07 2022 +0200 Merge pull request e0404#574 from e0404/hotfix/dicomDoseImport Fix RTDose import when no RTPlan is imported commit e469b8f Author: Niklas Wahl <n.wahl@dkfz.de> Date: Thu Aug 11 14:11:28 2022 +0200 Fix RTDose import when no RTPlan is imported
HomolkaN · Oct 12, 2022 · da02eae · da02eae
1 parent 38f8bc0
commit da02eae
Show file tree

Hide file tree

Showing 20 changed files with 268 additions and 385 deletions.
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
@@ -1,5 +1,5 @@
 # This is a basic workflow to help you get started with Actions
-name: Testing
+name: Tests
 # Controls when the action will run. 
 on: [push, pull_request, workflow_dispatch]  
 jobs:

diff --git a/CITATION.cff b/CITATION.cff
@@ -0,0 +1,123 @@
+cff-version: 1.2.0
+message: "If you use matRad, consider citing it and the corresponding research paper"
+authors:
+- family-names: "Ackermann"
+  given-names: "Benjamin"
+- family-names: "Bangert"
+  given-names: "Mark"
+- family-names: "Bennan"
+  given-names: "Amit Ben Antony"
+- family-names: "Burigo"
+  given-names: "Lucas"
+- family-names: "Cabal"
+  given-names: "Gonzalo"
+- family-names: "Cisternas"
+  given-names: "Eduardo"
+- family-names: "Charton"
+  given-names: "Louis"
+- family-names: "Christiansen"
+  given-names: "Eric"
+- family-names: "Ecker"
+  given-names: "Swantje"
+- family-names: "Ellerbrock"
+  given-names: "Malte"
+- family-names: "Gabryś"
+  given-names: "Hubert"
+- family-names: "Handrack"
+  given-names: "Josefine"
+- family-names: "Heath"
+  given-names: "Emily"
+- family-names: "Hermann"
+  given-names: "Cindy"
+- family-names: "Homolka"
+  given-names: "Noa"
+- family-names: "Jäkel"
+  given-names: "Oliver"
+- family-names: "Klinge"
+  given-names: "Thomas"
+- family-names: "Mairani"
+  given-names: "Andrea"
+- family-names: "Mescher"
+  given-names: "Hennig"
+- family-names: "Müller"
+  given-names: "Lucas-Raphael"
+- family-names: "Neishabouri"
+  given-names: "Ahmad"
+- family-names: "Palkowitsch"
+  given-names: "Martina"
+- family-names: "Parodi"
+  given-names: "Katia"
+- family-names: "Pezzano"
+  given-names: "Giuseppe"
+- family-names: "Ramirez"
+  given-names: "Daniel"
+- family-names: "Stadler"
+  given-names: "Alexander"
+- family-names: "Ulrich"
+  given-names: "Silke"
+- family-names: "Wahl"
+  given-names: "Niklas"
+  orcid: "https://orcid.org/0000-0002-1451-223X"
+- family-names: "Welsch"
+  given-names: "Jona"
+- family-names: "Wieser"
+  given-names: "Hans-Peter"
+- family-names: "Winter"
+  given-names: "Johanna"
+- family-names: "Xu"
+  given-names: "Tong" 
+
+title: "matRad"
+version: 2.10.1
+doi: 10.5281/zenodo.3879616
+date-released: 2020-06-05
+url: "http://www.matRad.org"
+location: "Heidelberg"
+organization: "Deutsches Krebsforschungszentrum"
+preferred-citation:
+  type: article
+  authors:
+  - family-names: "Wieser"
+    given-names: "Hans-Peter"
+  - family-names: "Cisternas"
+    given-names: "Eduardo"
+  - family-names: "Wahl"
+    given-names: "Niklas"
+    orcid: "https://orcid.org/0000-0002-1451-223X"
+  - family-names: "Ulrich"
+    given-names: "Silke"
+  - family-names: "Stadler"
+    given-names: "Alexander"
+  - family-names: "Mescher"
+    given-names: "Hennig"
+  - family-names: "Müller"
+    given-names: "Lucas-Raphael"
+  - family-names: "Klinge"
+    given-names: "Thomas"
+  - family-names: "Gabryś"
+    given-names: "Hubert"
+  - family-names: "Burigo"
+    given-names: "Lucas"
+  - family-names: "Mairani"
+    given-names: "Andrea"
+  - family-names: "Ecker"
+    given-names: "Swantje"
+  - family-names: "Ackermann"
+    given-names: "Benjamin"
+  - family-names: "Ellerbrock"
+    given-names: "Malte"
+  - family-names: "Parodi"
+    given-names: "Katia"
+  - family-names: "Jäkel"
+    given-names: "Oliver"
+  - family-names: "Bangert"
+    given-names: "Mark"
+
+  doi: 10.1002/mp.12251
+  journal: "Medical Physics"
+  start: 2556 
+  end: 2568
+  title: "Development of the open-source dose calculation and optimization toolkit matRad"
+  issue: 6
+  volume: 46
+  year: 2017
diff --git a/MatRad_Config.m b/MatRad_Config.m
@@ -179,8 +179,6 @@ function setDefaultProperties(obj)
             obj.propOpt.defaultRunDAO = 0;
             obj.propOpt.defaultRunSequencing = 0;
 
-            obj.propMC.defaultCarbonEnergySpread = 0; %[%]
-
             obj.propMC.ompMC_defaultHistories = 1e6;
             obj.propMC.ompMC_defaultOutputVariance = false;
 

diff --git a/MatRad_MCemittanceBaseData.m b/MatRad_MCemittanceBaseData.m
@@ -23,16 +23,16 @@
     % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
     properties
-        machine         %matRad base data machine struct
-        bdl_path = ''   %stores path to generated file
-        nozzleToIso     %Nozzle to Isocenter Distance
-        smx             %Scanning magnet X to isocenter Distance
-        smy             %Scanning magnet y to isocenter Distance
-        monteCarloData  %MC Phase space data struct
-        selectedFocus   %array containing selected focus indices per energy
-        energyspread    %custom energy spread
-        matRad_cfg      %matRad config
-        rangeShifters   %Stores range shifters
+        machine                             %matRad base data machine struct
+        bdl_path = ''                       %stores path to generated file
+        nozzleToIso                         %Nozzle to Isocenter Distance
+        smx                                 %Scanning magnet X to isocenter Distance
+        smy                                 %Scanning magnet y to isocenter Distance
+        monteCarloData                      %MC Phase space data struct
+        selectedFocus                       %array containing selected focus indices per energy
+        defaultRelativeEnergySpread = 0;    %default energy spread
+        matRad_cfg                          %matRad config
+        rangeShifters                       %Stores range shifters
     end
 
     properties (SetAccess = private)
@@ -58,7 +58,6 @@
             end
 
             obj.matRad_cfg = MatRad_Config.instance();
-            obj.energyspread = obj.matRad_cfg.propMC.defaultCarbonEnergySpread;
 
             obj.machine = machine;
             obj.problemSigma = false;
@@ -258,15 +257,17 @@
                     % https://www.nist.gov/system/files/documents/2017/04/26/newstar.pdf
                     meanEnergy = @(x) 11.39 * x^0.628 + 11.24;
                     mcDataEnergy.MeanEnergy = meanEnergy(r80);
-                    mcDataEnergy.EnergySpread = obj.energyspread;
+                    % reading in a potential given energyspread could go here directly. How would you parse the energyspread
+                    % into the function? Through a field in the machine?
+                    mcDataEnergy.EnergySpread = obj.defaultRelativeEnergySpread;
                 case 'helium'
                     % Fit to Range-Energy relationship
                     % Data from "Update to ESTAR, PSTAR, and ASTAR Databases" - ICRU Report 90, 2014
                     % Normalized energy before fit (MeV/u)! Only used ranges [10 350] mm for fit
                     % https://www.nist.gov/system/files/documents/2017/04/26/newstar.pdf
                     meanEnergy = @(x) 7.57* x.^0.5848 + 3.063;
                     mcDataEnergy.MeanEnergy = meanEnergy(r80);
-                    mcDataEnergy.EnergySpread = obj.energyspread;
+                    mcDataEnergy.EnergySpread = obj.defaultRelativeEnergySpread;
                 otherwise
                     error('not implemented')
             end

diff --git a/README.md b/README.md
@@ -3,8 +3,8 @@
 [![Contributors](https://img.shields.io/github/contributors/e0404/matRad)](https://github.com/e0404/matRad/graphs/contributors)
 [![GitPitch](https://gitpitch.com/assets/badge.svg)](https://gitpitch.com/e0404/matRad/master) 
 
-[![TravisCI Build Status](https://travis-ci.org/e0404/matRad.svg?branch=dev_varRBErobOpt)](https://travis-ci.org/e0404/matRad)
-[![Azure Pipelines Build Status](https://dev.azure.com/e0404/matRad/_apis/build/status/e0404.matRad?branchName=dev_varRBErobOpt)](https://dev.azure.com/e0404/matRad)
+[![GitHub Build Status](https://github.com/e0404/matRad/actions/workflows/tests.yml/badge.svg)](https://github.com/e0404/matRad/actions/workflows/tests.yml)
+[![Azure Pipelines Build Status](https://dev.azure.com/e0404/matRad/_apis/build/status/e0404.matRad)](https://dev.azure.com/e0404/matRad)
 
 DOIs:
 - General DOI: [![DOI](https://zenodo.org/badge/doi/10.5281/zenodo.3879615.svg)](https://doi.org/10.5281/zenodo.3879615)
@@ -111,7 +111,7 @@ Copyright 2020 the matRad development team.
 matrad@dkfz.de
 
 All the elements of the compilation of matRad and Ipopt are free software. You can redistribute and/or modify matRad's source code version provided as files with .m and .mat extension under the terms of the GNU GENERAL PUBLIC LICENSE Version 3 (GPL v3). You can also add to matRad the Ipopt functionality by using the precompiled mex files of the Ipopt optimizer in object code version which are licensed under the Eclipse Public License Version 1.0 (EPL v1.0), also made available for download via https://projects.coin-or.org/Ipopt.
-matRad also contains interfaces to an open-source photon Monte Carlo dose calculation engine developed by Edgardo Doerner hosted on GitHub (http://github.com/edoerner/ompMC) and to the open-source proton Monte Carlo project MCsquare (www.openmcsquare.org) from UCLouvain, Louvain-la-Neuve, Belgium. Both interfaces are integrated into matRad as submodules.
+matRad also contains interfaces to an open-source photon Monte Carlo dose calculation engine developed by Edgardo Dörner hosted on GitHub (http://github.com/edoerner/ompMC) and to the open-source proton Monte Carlo project MCsquare (www.openmcsquare.org) from UCLouvain, Louvain-la-Neuve, Belgium. Both interfaces are integrated into matRad as submodules.
 
 In addition, we provide a matlab standalone version of the compilation of matRad and Ipopt, where the files of matRad and Ipopt are licensed under GPL v3 and EPL v1.0 respectively. The matlab standalone version is meant to be used by students for learning and practicing scientific programming and does not yet contain the interfaces to the aforementioned Monte Carlo dose calculation engines.
 

diff --git a/dicom/matRad_importDicom.m b/dicom/matRad_importDicom.m
@@ -141,10 +141,8 @@
     if ~(cellfun(@isempty,files.rtdose(1,1:2))) 
         fprintf('loading Dose files \n', structures(i).structName);
         % parse plan in order to scale dose cubes to a fraction based dose
-        if exist('pln','var')
-            if isfield(pln,'numOfFractions')
-                resultGUI = matRad_importDicomRTDose(ct, files.rtdose, pln);
-            end
+        if exist('pln','var') && ~isempty(pln) && isfield(pln,'numOfFractions')
+            resultGUI = matRad_importDicomRTDose(ct, files.rtdose, pln);
         else
             resultGUI = matRad_importDicomRTDose(ct, files.rtdose);
         end

diff --git a/matRad_calcCubes.m b/matRad_calcCubes.m
@@ -146,24 +146,6 @@
 end
 
 %% Final processing
-% Add non-processed MC tallies
-% Note that the tallies are already computed per beam and altogether
-if isfield(dij,'MC_tallies')
-    for f = 1:numel(dij.MC_tallies)
-        tally = dij.MC_tallies{f};
-        % skip tallies processed above
-        if ~isfield(resultGUI,tally) && isempty(strfind(lower(tally),'std'))
-            tallyCut = strsplit(tally,'_');
-            if size(tallyCut,2) > 1
-                beamNum = str2num(cell2mat(regexp(tally,'\d','Match')));
-                resultGUI.(tally) = reshape(full(dij.(tallyCut{1}){scenNum}(:,beamNum)),dij.doseGrid.dimensions);
-            else
-                resultGUI.(tally) = reshape(full(dij.(tally){scenNum} * wBeam),dij.doseGrid.dimensions);
-            end
-        end
-    end
-end
-
 % Remove suffix for RBExD if there's only one available
 if any(cellfun(@(teststr) ~isempty(strfind(lower(teststr),'alpha')), fieldnames(dij))) && isfield(dij,'RBE_models') && length(dij.RBE_models) == 1
     % Get fieldnames that include the specified RBE model

diff --git a/matRad_calcDoseDirectMC.m b/matRad_calcDoseDirectMC.m
@@ -87,7 +87,7 @@
 
 % calc resulting dose
 if ~isprop(pln.propMC,'externalCalculation') || ~pln.propMC.externalCalculation
-    if pln.multScen.totNumScen == 1
+    if pln.multScen.numOfCtScen == 1
         % calculate cubes; use uniform weights here, weighting with actual fluence
         % already performed in dij construction
         if size(dij.physicalDose{1},2) ~= dij.numOfBeams || size(dij.physicalDose{1},2) ~= numel(dij.beamNum)
@@ -114,9 +114,9 @@
 
     end
 
-    if pln.multScen.totNumScen ~= 1
+    if pln.multScen.numOfCtScen ~= 1
         resultGUI.accPhysicalDose = zeros(size(resultGUI.phaseDose{1}));
-        for i = 1:pln.multScen.totNumScen
+        for i = 1:pln.multScen.numOfCtScen
             resultGUI.accPhysicalDose = resultGUI.accPhysicalDose + resultGUI.phaseDose{i};
         end
     end

diff --git a/matRad_calcParticleDoseMCsquare.m b/matRad_calcParticleDoseMCsquare.m
@@ -398,8 +398,6 @@
                             dij.doseGrid.dimensions, ...
                             dij.totalNumOfBixels, ...
                             mask);
-
-                        dij.MC_tallies{1} = 'LET';
                     end
                 else
                     %Read dose cube
@@ -416,20 +414,17 @@
                             sparse(VdoseGrid,ones(numel(VdoseGrid),1), ...
                                 cube(VdoseGrid), ...
                                 dij.doseGrid.numOfVoxels,1);
-
-                        dij.MC_tallies{1} = 'LET';
                     end
 
-                    % Postprocessing for dij: Remove fields since this is already the combined dose over all bixels
-                    fields = {'MCsquareCalcOrder','bixelNum','rayNum','totalNumOfRays','totalNumOfBixels','numOfRaysPerBeam','numOfBeams'};
-                    for f = 1:length(fields)
-                        if isfield(dij,fields{f})
-                            dij = rmfield(dij,fields{f});
-                        end
-                    end
+                    % Postprocessing for dij:
+                    % This is already the combined dose over all bixels, so all parameters are 1 in this case
+                    dij = rmfield(dij,'MCsquareCalcOrder');
+
                     dij.numOfBeams = 1;
                     dij.beamNum = 1;
                     dij.bixelNum = 1;
+                    dij.rayNum = 1;
+                    dij.totalNumOfBixels = 1;
                     dij.totalNumOfRays = 1;
                     dij.numOfRaysPerBeam = 1;
                 end
@@ -470,6 +465,9 @@
     end   
 end
 
+% Order fields for easier comparison between different dijs
+dij = orderfields(dij);
+
 %% cd back
 cd(currFolder);
 

diff --git a/matRad_calcParticleDoseMCtopas.m b/matRad_calcParticleDoseMCtopas.m
@@ -1,4 +1,4 @@
-function topasCubes = matRad_calcParticleDoseMCtopas(ct,stf,pln,cst,calcDoseDirect)
+function dij = matRad_calcParticleDoseMCtopas(ct,stf,pln,cst,calcDoseDirect)
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % matRad TOPAS Monte Carlo proton dose calculation wrapper
 %   This calls a TOPAS installation (not included in matRad due to
@@ -91,7 +91,7 @@
 matRad_calcDoseInit;
 
 % for TOPAS we explicitly downsample the ct to the dose grid (might not be necessary in future versions with separated grids)
-[ctR,~,~] = pln.propMC.resampleGrid(ct,cst,pln,stf);
+[ctR,~,~] = matRad_resampleCTtoGrid(ct,cst,pln,stf);
 
 % overwrite CT grid in dij in case of modulation.
 if isfield(ctR,'ctGrid')
@@ -167,9 +167,11 @@
     % change director back to original directory
     cd(pln.propMC.workingDir);
 
-    % save dij and weights, they are needed for later reading the data back in
+    % Skip local calculation and data readout with this parameter. All necessary parameters to read the data back in
+    % later are stored in the MCparam file that is stored in the folder. The folder is generated in the working
+    % directory and the matRad_plan*.txt file can be manually called with TOPAS.
     if pln.propMC.externalCalculation
-        matRad_cfg.dispInfo('TOPAS simulation skipped for external calculation\n');
+        matRad_cfg.dispInfo(['TOPAS simulation skipped for external calculation\nFiles have been written to: "',replace(pln.propMC.workingDir,'\','\\'),'"']);
     else
         for ctScen = 1:ct.numOfCtScen
             for beamIx = 1:numel(stf)
@@ -233,11 +235,14 @@
 %% Simulation(s) finished - read out volume scorers from topas simulation
 % Skip readout if external files were generated
 if ~pln.propMC.externalCalculation
-    topasCubes = pln.propMC.readFiles(pln.propMC.workingDir);
+    dij = pln.propMC.readFiles(pln.propMC.workingDir);
 else
-    topasCubes = [];
+    dij = [];
 end
 
+% Order fields for easier comparison between different dijs
+dij = orderfields(dij);
+
 % manipulate isocenter back
 for k = 1:length(stf)
     stf(k).isoCenter = stf(k).isoCenter - pln.multScen.isoShift(shiftScen,:);

diff --git a/matRad_calcPhotonDose.m b/matRad_calcPhotonDose.m
@@ -320,7 +320,7 @@
 
                         % sample dose only for bixel based dose calculation
                         if ~isFieldBasedDoseCalc
-                            r0   = 25;   % [mm] sample beyond the inner core
+                            r0   = 20 + stf(i).bixelWidth;   % [mm] sample beyond the inner core
                             Type = 'radius';
                             [ix,bixelDose] = matRad_DijSampling(ix,bixelDose,manipulatedRadDepthCube,rad_distancesSq,Type,r0);
                         end