Merge branch 'dev_MixedModality' of https://github.com/e0404/matRad i…

…nto dev_MixedModality

.gitattributes

@@ -1,5 +0,0 @@
-phantoms/BOXPHANTOM.mat filter=lfs diff=lfs merge=lfs -text
-phantoms/TG119.mat filter=lfs diff=lfs merge=lfs -text
-phantoms/HEAD_AND_NECK.mat filter=lfs diff=lfs merge=lfs -text
-phantoms/PROSTATE.mat filter=lfs diff=lfs merge=lfs -text
-phantoms/LIVER.mat filter=lfs diff=lfs merge=lfs -text

.travis/before_install_linux.sh → .github/before_install_linux.sh

@@ -1,6 +1,6 @@
 #!/usr/bin/env bash
 
-sudo chmod +x .travis/runtests.sh
+sudo chmod +x .github/runtests.sh
 sudo chmod +x MCsquare/bin/MCsquare_linux
 
 mv optimization/optimizer/ipopt.m optimization/optimizer/ipopt.m.bak

.travis/runtests.sh → .github/runtests.sh

@@ -4,7 +4,8 @@
 ## Make sure some failures are detected by the CI runners
 function exitIfError {
 	# pass "$?" as argument: i.e. the exit status of the last call
-	if [ "$1" -ne 0 ]; then
+	# currently octave 6 can finish with a segfault when the program is closed due to some bug, for now we try to ignore it
+	if [ "$1" -ne 0 ] && [ "$1" -ne 139 ]; then
 		exit $1;
 	fi
 }

.github/workflows/tests.yml

@@ -3,36 +3,53 @@ name: Tests
 # Controls when the action will run. 
 on: [push, pull_request, workflow_dispatch]  
 jobs:
-  test-matlab: #Matlab test Job
+  test-matlab-stable: #Matlab test Job for supported Release
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2 # Checks-out repository under $GITHUB_WORKSPACE
+      - uses: actions/checkout@v3 # Checks-out repository under $GITHUB_WORKSPACE
       # Install MATLAB
       - name: Install MATLAB
-        uses: matlab-actions/setup-matlab@v1        
+        uses: matlab-actions/setup-matlab@v1
+        with:
+          release: R2022b
+      # Runs test command
+      - name: Run Tests
+        uses: matlab-actions/run-command@v1
+        with:
+          command: cd test; matRad_runTests
+  test-matlab-latest: #Matlab test Job for latest Matlab release
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v3 # Checks-out repository under $GITHUB_WORKSPACE
+      # Install MATLAB
+      - name: Install MATLAB
+        uses: matlab-actions/setup-matlab@v1
+        with:
+          release: latest
       # Runs test command
       - name: Run Tests
         uses: matlab-actions/run-command@v1
         with:
           command: cd test; matRad_runTests
-  test-octave: #Octave test Job
-    runs-on: ubuntu-20.04 # We use Ubuntu-20.04 because it has Octave 5.2
+  test-octave-6: #Octave test Job
+    runs-on: ubuntu-22.04 # We use Ubuntu-22.04 because it has Octave 6.4
     steps:
-      - uses: actions/checkout@v2 # Checks-out repository under $GITHUB_WORKSPACE
+      - uses: actions/checkout@v3 # Checks-out repository under $GITHUB_WORKSPACE
       - name: Install OCTAVE
-        run: |
-          sudo apt-get update
+        run: | 
+          sudo apt update
           sudo apt-get install -y gdb gfortran fonts-freefont-otf gnuplot-x11 libgdcm-dev octave liboctave-dev
       - name: Prepare Test Environment
         run: |
-          sudo chmod +x .travis/before_install_linux.sh
-          sudo .travis/before_install_linux.sh
-      - name: Run Tests
-        uses: GabrielBB/xvfb-action@v1 #For Headless tests
-        with:
-          run: .travis/runtests.sh octave-cli
-      - name: Upload Logs if failure
-        uses: actions/upload-artifact@v2
+          sudo chmod +x .github/before_install_linux.sh
+          sudo .github/before_install_linux.sh
+      - name: Run Tests 
+        run: xvfb-run -a .github/runtests.sh octave-cli
+        # uses: GabrielBB/xvfb-action@v1 #For Headless tests
+        # with:
+        #   run: .github/runtests.sh octave-cli
+      - name: Upload logs if test fails
+        uses: actions/upload-artifact@v3
         if: failure()
         with:
           name: Test Log

.gitignore

@@ -1 +1,2 @@
 MCsquare/bin/Materials/
+topas/MCrun

4D/matRad_addMovement.m

@@ -1,4 +1,4 @@
-function [ct, cst] = matRad_addMovement(ct, cst, motionPeriod, numOfCtScen, amp,visBool)
+function [ct, cst] = matRad_addMovement(ct, cst, motionPeriod, numOfCtScen, amp, varargin)
 % adds artificial sinosodal patient motion by creating a deformation vector
 % field and applying it to the ct.cube by geometric transformation
 %
@@ -11,7 +11,8 @@
 %   motionPeriod:   the length of a whole breathing cycle (in seconds)
 %   numOfCtScen:    number of ct phases
 %   amp:            amplitude of the sinosoidal movement (in pixels)
-%   visBool         boolean flag for visualization
+%   varargin:   dvfType:        push or pull dvf
+%               visBool         boolean flag for visualization
 %
 %   note:           1st dim --> x LPS coordinate system
 %                   2nd dim --> y LPS coordinate system
@@ -40,9 +41,12 @@
 %
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-if ~exist('visBool','var')
-    visBool = false;
-end
+expectedDVF = {'pull', 'push'};
+
+p = inputParser; 
+addParameter(p,'dvfType','pull', @(x) any(validatestring(x,expectedDVF)))
+addParameter(p,'visBool',false, @islogical);
+parse(p,varargin{:});
 
 matRad_cfg = MatRad_Config.instance();
 
@@ -51,7 +55,10 @@
 ct.numOfCtScen = numOfCtScen;
 
 % set type
-ct.dvfType = 'pull'; % push or pull
+ct.dvfMetadata.dvfType = p.Results.dvfType;
+if strcmp(p.Results.dvfType,'push')
+    amp = -amp; %dvf_pull = -dvf_push;
+end
 
 env = matRad_getEnvironment();
 
@@ -110,17 +117,16 @@
     end
 
     % convert dvfs to [mm]
-    tmp = ct.dvf{i}(:,:,:,1);
-    ct.dvf{i}(:,:,:,1) = -ct.dvf{i}(:,:,:,2) * ct.resolution.x;
-    ct.dvf{i}(:,:,:,2) = -tmp * ct.resolution.y;
-    ct.dvf{i}(:,:,:,3) = -ct.dvf{i}(:,:,:,3) * ct.resolution.z;
+    ct.dvf{i}(:,:,:,1) = ct.dvf{i}(:,:,:,1).* ct.resolution.x;
+    ct.dvf{i}(:,:,:,2) = ct.dvf{i}(:,:,:,2).*ct.resolution.y;
+    ct.dvf{i}(:,:,:,3) = ct.dvf{i}(:,:,:,3).* ct.resolution.z;
 
     ct.dvf{i} = permute(ct.dvf{i}, [4,1,2,3]);    
 end
 
 
 
-if visBool
+if p.Results.visBool
     slice = round(ct.cubeDim(3)/2);
     figure,
     for i = 1:numOfCtScen

4D/matRad_calc4dDose.m

@@ -1,4 +1,4 @@
-function [resultGUI, timeSequence] = matRad_calc4dDose(ct, pln, dij, stf, cst, resultGUI, totalPhaseMatrix)
+function [resultGUI, timeSequence] = matRad_calc4dDose(ct, pln, dij, stf, cst, resultGUI, totalPhaseMatrix,accType)
 % wrapper for the whole 4D dose calculation pipeline and calculated dose
 % accumulation
 %
@@ -13,6 +13,7 @@
 %   cst:             matRad cst struct
 %   resultGUI:       struct containing optimized fluence vector
 %   totalPhaseMatrix optional intput for totalPhaseMatrix
+%   accType:         witch algorithim for dose accumulation
 % output
 %   resultGUI:      structure containing phase dose, RBE weighted dose, etc
 %   timeSequence:   timing information about the irradiation
@@ -33,6 +34,10 @@
 %
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+if ~exist('accType','var')
+    accType = 'DDM';
+end
+
 if ~exist('totalPhaseMatrix','var')
    % make a time sequence for when each bixel is irradiated, the sequence
    % follows the backforth spot scanning
@@ -62,7 +67,7 @@
     elseif strcmp(pln.bioParam.model,'constRBE')
         resultGUI.phaseRBExD{i} = tmpResultGUI.RBExD;
     % compute all fields
-    elseif any(strcmp(pln.bioParam.model,{'MCN','LEM','WED'}))
+    elseif any(strcmp(pln.bioParam.model,{'MCN','LEM','WED','HEL'}))
         resultGUI.phaseAlphaDose{i}    = tmpResultGUI.alpha .* tmpResultGUI.physicalDose;
         resultGUI.phaseSqrtBetaDose{i} = sqrt(tmpResultGUI.beta) .* tmpResultGUI.physicalDose;
     end
@@ -72,24 +77,24 @@
 % accumulation
 if strcmp(pln.bioParam.model,'none')       
 
-    resultGUI.accPhysicalDose = matRad_doseAcc(ct,resultGUI.phaseDose, cst, 'DDM');
+    resultGUI.accPhysicalDose = matRad_doseAcc(ct,resultGUI.phaseDose, cst, accType);
 
 elseif strcmp(pln.bioParam.model,'constRBE')
 
-    resultGUI.accRBExD = matRad_doseAcc(ct,resultGUI.phaseRBExD, cst, 'DDM');
+    resultGUI.accRBExD = matRad_doseAcc(ct,resultGUI.phaseRBExD, cst, accType);
 
-elseif any(strcmp(pln.bioParam.model,{'MCN','LEM','WED'}))
+elseif any(strcmp(pln.bioParam.model,{'MCN','LEM','WED','HEL'}))
 
-    resultGUI.accAlphaDose    = matRad_doseAcc(ct,resultGUI.phaseAlphaDose, cst, 'DDM');
-    resultGUI.accSqrtBetaDose = matRad_doseAcc(ct,resultGUI.phaseSqrtBetaDose, cst, 'DDM');
+    resultGUI.accAlphaDose    = matRad_doseAcc(ct,resultGUI.phaseAlphaDose, cst,accType);
+    resultGUI.accSqrtBetaDose = matRad_doseAcc(ct,resultGUI.phaseSqrtBetaDose, cst, accType);
 
     % only compute where we have biologically defined tissue
-    ix = dij.alphaX~=0; 
+    ix = dij.ax(:,1)~=0; 
 
     resultGUI.accEffect = resultGUI.accAlphaDose + resultGUI.accSqrtBetaDose.^2;
 
     resultGUI.accRBExD     = zeros(ct.cubeDim);
-    resultGUI.accRBExD(ix) = ((sqrt(dij.alphaX(ix).^2 + 4 .* dij.betaX(ix) .* resultGUI.accEffect(ix)) - dij.alphaX(ix))./(2.*dij.betaX(ix)));
+    resultGUI.accRBExD(ix) = ((sqrt(dij.ax(ix).^2 + 4 .* dij.bx(ix) .* resultGUI.accEffect(ix)) - dij.ax(ix))./(2.*dij.bx(ix)));
 
 end
 

4D/matRad_doseAcc.m

@@ -57,14 +57,14 @@
 
 if strcmp(accMethod,'DDM')
 
-    if ~strcmp(ct.dvfType,'pull')
+    if ~strcmp(ct.dvfMetadata.dvfType,'pull')
         error('dose accumulation via direct dose mapping (DDM) requires pull dvfs');
     end    
 
-    [Y,X,Z] = meshgrid(yGridVec,xGridVec,zGridVec);
+    [X,Y,Z] = meshgrid(xGridVec,yGridVec,zGridVec);
 
     % TODODODODODOD %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    ix = [];1:prod(ct.cubeDim);%[];
+    ix = [];
     for i = 1:size(cst,1)
         ix = unique([ix; cst{i,4}{1}]);
     end
@@ -74,20 +74,17 @@
         dvf_x_i = squeeze(ct.dvf{1,i}(1,:,:,:))/ct.resolution.x;
         dvf_y_i = squeeze(ct.dvf{1,i}(2,:,:,:))/ct.resolution.y;
         dvf_z_i = squeeze(ct.dvf{1,i}(3,:,:,:))/ct.resolution.z;
-
-        d_ref = matRad_interp3(yGridVec,xGridVec',zGridVec,phaseCubes{i}, ...
-                         Y(ix) + dvf_y_i(ix), ...     
-                         X(ix) + dvf_x_i(ix), ... 
-                         Z(ix) + dvf_z_i(ix), ...
-                         'linear',0);  
 
-        dAcc(ix) = dAcc(ix) + d_ref;
+        d_ref = matRad_interp3(X, Y, Z,...
+                            phaseCubes{i}, X-dvf_x_i,Y-dvf_y_i,Z-dvf_z_i,'linear',0);   
+
+        dAcc(ix) = dAcc(ix) + d_ref(ix);
 
     end
 
 elseif strcmp(accMethod,'EMT')   % funktioniert nicht wenn Dosis in einer Phase = 0 ist...
 
-    if ~strcmp(ct.dvfType,'push')
+    if ~strcmp(ct.dvfMetadata.dvfType,'push')
         error('dose accumulation via interpolation requires push dvfs');
     end
 
@@ -99,8 +96,10 @@
         dvf_y_i = squeeze(ct.dvf{1,i}(2,:,:,:))/ct.resolution.y;
         dvf_z_i = squeeze(ct.dvf{1,i}(3,:,:,:))/ct.resolution.z;
 
-        m_i     = ct.cube{i};
-        e_i     = phaseCubes{i}.*m_i;
+        m_i = ct.cube{i};
+        m_i = permute(m_i,[2 1 3]);
+        e_i = phaseCubes{i}.*ct.cube{i};
+        e_i     = permute(e_i,[2 1 3]);
 
         ix = e_i>0;
 
@@ -163,6 +162,7 @@
          k = find(m_ref);
          dAcc(k) = dAcc(k) + e_ref(k)./m_ref(k);
 
+          dAcc = permute(dAcc,[2 1 3]);
     end
 
 % elseif strcmp(accMethod,'DDMM')

AUTHORS.txt

@@ -1,6 +1,9 @@
 List of all matRad developers that contributed code (alphabetical)
 
+* Nelly Abbani
+* Nabe Al-Hasnawi
 * Mark Bangert
+* Tobias Becher
 * Amit Ben Antony Bennan
 * Lucas Burigo
 * Gonzalo Cabal

ChangeLog.txt

@@ -1,3 +1,37 @@
+Monte-Carlo Update
+    Workflow of the Monte Carlo pipeline including MCsquare and TOPAS have been completely overhauled
+    - Restructured the MCEmittanceBaseData class, fit and calculation pipeline
+        - Added example 13 for generating analytical data file by fitting to given machine emittance
+            - added function to plot particleBaseDataEntry
+            - edited function to fit base data
+            - added function to generate a single pencil beam
+        - MCemmittanceBaseData can calculate meanEnergy and spread also for carbon and helium ions
+    - Added support for 4D calculation
+        - added 4D accumulation in calcDoseDirectMC
+        - Fixed number of errors in 4D calculation workflow
+    - TOPAS Updates:
+        - Added comments to whole pipeline
+        - Implemented dij calculation
+        - Restructured resampling in calcParticleDoseMCtopas in separate function
+        - calcParticleDoseMCtopas now generates dij in matRad format for calcDoseDirect
+        - Merged support functions for TOPAS into topasConfig class
+        - modular Schneider Converter: Converter is generated on demand rather than read from file which allows a variety of different options
+        - modular TOPAS scorer which can be individually turned on and off
+        - Export feature for TOPAS to run externally (includes functions to read from external folders) 
+    - MCsquare Updates:
+        - Merged support functions for MCsquare into MCsquareConfig class
+        - added variable in pln to contain an existing BDL file if available
+        - renamed MCsquare property "Num_Primaries" -> "numHistories" to be in line with other Monte Carlo (this is written to the BDL file in original format)
+    - Implemented calculation of std of physicalDose
+
+    Changes to matRad workflow:
+    - Added a class constructor for pln in MatRad_config, which loads requested classes in pln and writes default values that were not set
+    - Changed matRad_calcCubes to accept a variety of different fields for Monte Carlo, without changing the current usage
+    - Enabled helium in calcdoseDirect
+    - Bugfix for coordinate system in resizeCstToGrid function 
+    - Added optional initial weights to fluence Optimization
+    - Added flag in stf to catch specific error where no energies could be found
+
 Development Changes
 	- Bugfixes:
 		- Fix SFUD wrapper for new pln & cst format