Dev mix mod effect prescription bug fix

.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_calc4dDose.m

@@ -84,12 +84,12 @@
     resultGUI.accSqrtBetaDose = matRad_doseAcc(ct,resultGUI.phaseSqrtBetaDose, cst, 'DDM');
 
     % 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

@@ -61,10 +61,10 @@
         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,8 +74,8 @@
         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}, ...
+
+        d_ref = matRad_interp3(yGridVec,xGridVec,zGridVec,permute(phaseCubes{i},[2 1 3]), ...
                          Y(ix) + dvf_y_i(ix), ...     
                          X(ix) + dvf_x_i(ix), ... 
                          Z(ix) + dvf_z_i(ix), ...

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

IO/matRad_readNRRD.m

@@ -27,23 +27,25 @@
 %
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+matRad_cfg = MatRad_Config.instance();
+
 % Open file.
 hFile = fopen(filename, 'r');
 if hFile <= 0
-    error('Could not open NRRD file!');
+    matRad_cfg.dispError('Could not open NRRD file!');
 end
 cleaner = onCleanup(@() fclose(hFile));
 
 %% Determine NRRD Version
 nrrdLine = fgetl(hFile);
 regTokens = regexp(nrrdLine,'NRRD00(?:\.)?0([1-5])','tokens');
 if isempty(regTokens)
-    error('Invalid Header line! Could not identify NRRD version!');
+    matRad_cfg.dispError('Invalid Header line! Could not identify NRRD version!');
 end
 nrrdVersion = str2num(regTokens{1}{1});
 
 if nrrdVersion > 5
-    error('NRRD version > 5 not supported!');
+    matRad_cfg.dispError('NRRD version > 5 not supported!');
 end
 
 %% Read header
@@ -60,15 +62,15 @@
         %Parse the line
         lineContent = regexp(currentLine, '(.+):(=|\s)(.+)', 'tokens');
         if isempty(lineContent)
-            warning(['Could not parse line: "' lineContent '"']);
+            matRad_cfg.dispWarning(['Could not parse line: "' currentLine '"']);
         elseif isequal(lineContent{1}{2},' ') %space after colon refers to "field"
             nrrdMetaData.fields{end+1,1} = lineContent{1}{1}; %Fieldname
             nrrdMetaData.fields{end,2} = lineContent{1}{3}; %Information
         elseif isequal(lineContent{1}{2},'=') %= after colon refers to key-value pair
             nrrdMetaData.keys{end+1,1} = lineContent{1}{1}; %Key
             nrrdMetaData.keys{end,2} = lineContent{1}{3}; %Value
         else
-            warning(['Could not parse line: "' lineContent '"']);
+            matRad_cfg.dispWarning(['Could not parse line: "' currentLine '"']);
         end        
     end
     currentLine = fgetl(hFile);
@@ -85,7 +87,7 @@
         endianFieldIx = find(ismember(nrrdMetaData.fields(:,1),'endian'));
         if ~isempty(endianFieldIx)
             if ~isequal(nrrdMetaData.fields{endianFieldIx,2},'little') && ~isequal(nrrdMetaData.fields{endianFieldIx,2},'big')
-                error(['Datatype is ' datatype ', thus endian information is required but could not be interpreted!']);
+                matRad_cfg.dispError(['Datatype is ' datatype ', thus endian information is required but could not be interpreted!']);
             end;
             %Now we compare the file endian to the system endian
             %First acquire system endian
@@ -101,36 +103,36 @@
                 doSwapBytes = true;
             end
         else
-            error(['Datatype is ' datatype ', thus endian information is required but could not be found!']);
+            matRad_cfg.dispError(['Datatype is ' datatype ', thus endian information is required but could not be found!']);
         end
     end
     metadata.datatype = datatype;
 
 else
-    error('Could not find required "type" field!');
+    matRad_cfg.dispError('Could not find required "type" field!');
 end
 
 %Check for the always required image dimension
 dimFieldIx = find(ismember(nrrdMetaData.fields(:,1), 'dimension')); 
 if ~isempty(dimFieldIx)
     [metadata.dimension,success] = str2num(nrrdMetaData.fields{dimFieldIx,2});
     if ~success
-        error('Could not read required dimension field');
+        matRad_cfg.dispError('Could not read required dimension field');
     end
 else
-    error('Could not find required "dimension" field!');
+    matRad_cfg.dispError('Could not find required "dimension" field!');
 end
 
 %Check for size / dim length
 sizeFieldIx = find(ismember(nrrdMetaData.fields(:,1), 'sizes')); 
 if ~isempty(sizeFieldIx)
     sizes = textscan(nrrdMetaData.fields{sizeFieldIx,2},'%d');
     if numel(sizes{1}) ~= metadata.dimension || ~all(sizes{1} > 0) 
-        error('Incorrect size definition!');
+        matRad_cfg.dispError('Incorrect size definition!');
     end
     metadata.cubeDim = sizes{1}';
 else
-    error('Could not find required "dimension" field!');
+    matRad_cfg.dispError('Could not find required "dimension" field!');
 end
 
 %Check for resolution
@@ -140,7 +142,7 @@
 if ~isempty(spacingFieldIx)
     resolutions = textscan(nrrdMetaData.fields{spacingFieldIx,2},'%f');
     if numel(resolutions{1}) ~= metadata.dimension
-        error('Incorrect spacings definition');
+        matRad_cfg.dispError('Incorrect spacings definition');
     end
     metadata.resolution = resolutions{1}';    
 
@@ -167,14 +169,14 @@
         currentAxis = find(vectors{c});
 
         if numel(find(vectors{c})) ~= 1
-            error('Sorry! We currently only support spaces with cartesian basis!'); 
+            matRad_cfg.dispError('Sorry! We currently only support spaces with cartesian basis!'); 
         end        
         metadata.axisPermutation(c) = currentAxis*sign(vectors{c}(currentAxis));
         metadata.resolution(c) = vectors{c}(currentAxis);       
     end
 
 else
-    warning('No Resolution Information available');
+    matRad_cfg.dispWarning('No Resolution Information available');
 end
 
 %find the origin if we have one
@@ -202,7 +204,7 @@
 data_fileFieldIx = find(ismember(nrrdMetaData.fields(:,1), 'data file'));
 datafileFieldIx = find(ismember(nrrdMetaData.fields(:,1), 'datafile'));
 if ~isempty(data_fileFieldIx) || ~isempty(datafileFieldIx)
-    error('Sorry! We currently do not support detached data files!');
+    matRad_cfg.dispError('Sorry! We currently do not support detached data files!');
     %Proposed workflow:
     %check for data file    
     %close file
@@ -215,15 +217,15 @@
 %Check for encoding
 encodingFieldIx = find(ismember(nrrdMetaData.fields(:,1), 'encoding'));
 if isempty(encodingFieldIx)
-    error('Could not find required "encoding" field!');
+    matRad_cfg.dispError('Could not find required "encoding" field!');
 end
 switch nrrdMetaData.fields{encodingFieldIx,2}
     case 'raw'
         cube = fread(hFile,prod(metadata.cubeDim), metadata.datatype);
     case {'txt','text','ascii'}
         cube = cast(fscanf(hFile,'%f'),metadata.datatype);
     case 'hex'
-        error('Sorry: NRRD hex file not yet supported!');
+        matRad_cfg.dispError('Sorry: NRRD hex file not yet supported!');
     case {'gz','gzip'}
         compressedByteArray = fread(hFile, inf, 'uint8');
 
@@ -242,7 +244,7 @@
                 javaByteOutputStream.close();
                 javaUnpackSuccessful = true;
             catch
-                warning('Java unpacking failed... using temporary files!');
+                matRad_cfg.dispWarning('Java unpacking failed... using temporary files!');
             end
         end
         if ~javaUnpackSuccessful
@@ -252,7 +254,7 @@
             hFileTmp = fopen(tmpFile, 'wb');
 
             if hFileTmp <= 0
-                error('Could not open temporary file for GZIP!');
+                matRad_cfg.dispError('Could not open temporary file for GZIP!');
             end
 
             fwrite(hFileTmp, compressedByteArray, 'uint8');
@@ -264,16 +266,16 @@
             %Read the uncompressed file
             hFileTmp = fopen(tmpName, 'rb');
             if hFileTmp <= 0
-                error('Could not open unpacked file!');
+                matRad_cfg.dispError('Could not open unpacked file!');
             end
             cleanTmpFile = onCleanup(@() fclose(hFileTmp));
 
             cube = fread(hFileTmp, prod(metadata.cubeDim), metadata.datatype);
         end
     case {'bz2','bzip2'}
-        error('Sorry: bzip compression not yet supported!');
+        matRad_cfg.dispError('Sorry: bzip compression not yet supported!');
     otherwise 
-        error(['Undefined NRRD encoding scheme: ' nrrdMetaData.encoding]);
+        matRad_cfg.dispError(['Undefined NRRD encoding scheme: ' nrrdMetaData.encoding]);
 end
 
 %maybe we need to correct the byte ordering (endian)
@@ -327,7 +329,7 @@
     case 'double'
         datatype = 'double';        
     otherwise
-        error('Could not identify datatype for NRRD data');
+        matRad_cfg.dispError('Could not identify datatype for NRRD data');
 end
 
 end