From 473fd3b81a10d9de89aba76cfc44f80e4892388e Mon Sep 17 00:00:00 2001
From: Stanic <g906c@ad.dkfz-heidelberg.de>
Date: Mon, 21 Nov 2022 10:21:03 +0100
Subject: [PATCH 1/4] Update matRad_saveStructs.m

---
 matRad_saveStructs.m | 13 ++++++++++---
 1 file changed, 10 insertions(+), 3 deletions(-)

diff --git a/matRad_saveStructs.m b/matRad_saveStructs.m
index bfb3681b9..1e7695b6c 100644
--- a/matRad_saveStructs.m
+++ b/matRad_saveStructs.m
@@ -1,4 +1,4 @@
-function status = matRad_saveStructs(path)
+function status = matRad_saveStructs(path, engine)
 % matRad_saveStructs Mat file transfer for python interface
 %
 % input
@@ -20,8 +20,15 @@
 
 load(path);
 
-save('ct.mat', 'ct');
-save('cst.mat', 'cst');
+if engine=='matlab'
+    save('ct.mat', 'ct');
+    save('cst.mat', 'cst');
+else
+    save('ct.mat', '-mat7-binary', 'ct');
+    save('cst.mat', '-mat7-binary', 'cst');
+end
+
+%Choosing the engine is necessary because Octave has trouble reading .mat files. Doesn't recognize them as binary.
 
 status = 'Files written';
 

From 72f9e5accbfd1e7780798610830533c55c00f73f Mon Sep 17 00:00:00 2001
From: g-stanic <goran.stanic@dkfz-heidelberg.de>
Date: Fri, 24 Feb 2023 16:06:57 +0100
Subject: [PATCH 2/4] IO_stf added

This functions helps load and save stf that comes from the StfGenerator_temporary.py
---
 IO_stf.asv             |  21 ++
 IO_stf.m               |  23 ++
 matRad_generateStf.asv | 516 +++++++++++++++++++++++++++++++++++++++++
 matRad_generateStf.m   |  10 +-
 4 files changed, 569 insertions(+), 1 deletion(-)
 create mode 100644 IO_stf.asv
 create mode 100644 IO_stf.m
 create mode 100644 matRad_generateStf.asv

diff --git a/IO_stf.asv b/IO_stf.asv
new file mode 100644
index 000000000..5593dfe55
--- /dev/null
+++ b/IO_stf.asv
@@ -0,0 +1,21 @@
+clear all;
+matRad_rc;
+
+%%
+stf_init = load('C:\Users\g906c\Projects\pyRadPlan\stf_with_separate_rays.mat', 'stf');
+ray_init = load('C:\Users\g906c\Projects\pyRadPlan\stf_with_separate_rays.mat', 'rays');
+
+%%
+stf = [stf_init.stf{:}];
+ray = cell(size(ray_init.rays, 2), 1);
+
+for i=1:size(ray_init.rays, 2)
+    ray{i} = [ray_init.rays{i}{:}];
+end
+
+for i=1:size(ray_init.rays, 2)
+    stf(i).ray = ray{i};
+end
+
+%%
+save('stf.m')
\ No newline at end of file
diff --git a/IO_stf.m b/IO_stf.m
new file mode 100644
index 000000000..665e7a25f
--- /dev/null
+++ b/IO_stf.m
@@ -0,0 +1,23 @@
+function status = IO_stf()
+
+stf_init = load('C:\Users\g906c\Projects\pyRadPlan\stf_with_separate_rays.mat', 'stf');
+ray_init = load('C:\Users\g906c\Projects\pyRadPlan\stf_with_separate_rays.mat', 'rays');
+
+%%
+stf = [stf_init.stf{:}];
+ray = cell(size(ray_init.rays, 2), 1);
+
+for i=1:size(ray_init.rays, 2)
+    ray{i} = [ray_init.rays{i}{:}];
+end
+
+for i=1:size(ray_init.rays, 2)
+    stf(i).ray = ray{i};
+end
+
+%%
+save('C:\Users\g906c\Projects\pyRadPlan\stf.mat', 'stf')
+
+status = 'STF written';
+
+end
\ No newline at end of file
diff --git a/matRad_generateStf.asv b/matRad_generateStf.asv
new file mode 100644
index 000000000..e4ce32668
--- /dev/null
+++ b/matRad_generateStf.asv
@@ -0,0 +1,516 @@
+function stf = matRad_generateStf(ct,cst,pln,visMode)
+% matRad steering information generation
+% 
+% call
+%   stf = matRad_generateStf(ct,cst,pln,visMode)
+%
+% input
+%   ct:         ct cube
+%   cst:        matRad cst struct
+%   pln:        matRad plan meta information struct
+%   visMode:    toggle on/off different visualizations by setting this value to 1,2,3 (optional)
+%
+% output
+%   stf:        matRad steering information struct
+%
+% References
+%   -
+%
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Copyright 2015 the matRad development team. 
+% 
+% This file is part of the matRad project. It is subject to the license 
+% terms in the LICENSE file found in the top-level directory of this 
+% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part 
+% of the matRad project, including this file, may be copied, modified, 
+% propagated, or distributed except according to the terms contained in the 
+% LICENSE file.
+%
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+matRad_cfg = MatRad_Config.instance();
+
+matRad_cfg.dispInfo('matRad: Generating stf struct... ');
+
+if nargin < 4
+    visMode = 0;
+end
+
+if numel(pln.propStf.gantryAngles) ~= numel(pln.propStf.couchAngles)
+    matRad_cfg.dispError('Inconsistent number of gantry and couch angles.');
+end
+
+if ~isnumeric(pln.propStf.bixelWidth) || pln.propStf.bixelWidth < 0 || ~isfinite(pln.propStf.bixelWidth)
+   matRad_cfg.dispError('bixel width (spot distance) needs to be a real number [mm] larger than zero.');
+end
+
+% find all target voxels from cst cell array
+V = [];
+for i=1:size(cst,1)
+    if isequal(cst{i,3},'TARGET') && ~isempty(cst{i,6})
+        V = [V;vertcat(cst{i,4}{:})];
+    end
+end
+
+disp(length(V));
+
+% Remove double voxels
+V = unique(V);
+% generate voi cube for targets
+voiTarget    = zeros(ct.cubeDim);
+voiTarget(V) = 1;
+    
+% add margin
+addmarginBool = matRad_cfg.propStf.defaultAddMargin;
+if isfield(pln,'propStf') && isfield(pln.propStf,'addMargin')
+   addmarginBool = pln.propStf.addMargin; 
+end
+
+if addmarginBool
+    voiTarget = matRad_addMargin(voiTarget,cst,ct.resolution,ct.resolution,true);
+    V   = find(voiTarget>0);
+end
+
+% throw error message if no target is found
+if isempty(V)
+    matRad_cfg.dispError('Could not find target.');
+end
+
+% Convert linear indices to 3D voxel coordinates
+[coordsY_vox, coordsX_vox, coordsZ_vox] = ind2sub(ct.cubeDim,V);
+
+% prepare structures necessary for particles
+fileName = [pln.radiationMode '_' pln.machine];
+try
+   load([fileparts(mfilename('fullpath')) filesep 'basedata' filesep fileName]);
+   SAD = machine.meta.SAD;
+catch
+   matRad_cfg.dispError('Could not find the following machine file: %s',fileName); 
+end
+
+if strcmp(pln.radiationMode,'protons') || strcmp(pln.radiationMode,'carbon')
+      
+    availableEnergies = [machine.data.energy];
+    availablePeakPos  = [machine.data.peakPos] + [machine.data.offset];
+    
+    if sum(availablePeakPos<0)>0
+       matRad_cfg.dispError('at least one available peak position is negative - inconsistent machine file') 
+    end
+    %clear machine;
+end
+
+% calculate rED or rSP from HU
+ct = matRad_calcWaterEqD(ct, pln);
+
+% take only voxels inside patient
+V = [cst{:,4}];
+V = unique(vertcat(V{:}));
+
+% ignore densities outside of contours
+eraseCtDensMask = ones(prod(ct.cubeDim),1);
+eraseCtDensMask(V) = 0;
+for i = 1:ct.numOfCtScen
+    ct.cube{i}(eraseCtDensMask == 1) = 0;
+end
+
+% Define steering file like struct. Prellocating for speed.
+stf = struct;
+
+save('coordsX_vox', 'coordsX_vox');
+save('coordsY_vox', 'coordsY_vox');
+save('coordsZ_vox', 'coordsZ_vox');
+
+% loop over all angles
+for i = 1:length(pln.propStf.gantryAngles)
+    
+    % Correct for iso center position. Whit this correction Isocenter is
+    % (0,0,0) [mm]
+    coordsX = coordsX_vox*ct.resolution.x - pln.propStf.isoCenter(i,1);
+    coordsY = coordsY_vox*ct.resolution.y - pln.propStf.isoCenter(i,2);
+    coordsZ = coordsZ_vox*ct.resolution.z - pln.propStf.isoCenter(i,3);
+    
+    save('coordsX', 'coordsX');
+
+    % Save meta information for treatment plan
+    stf(i).gantryAngle   = pln.propStf.gantryAngles(i);
+    stf(i).couchAngle    = pln.propStf.couchAngles(i);
+    stf(i).bixelWidth    = pln.propStf.bixelWidth;
+    stf(i).radiationMode = pln.radiationMode;
+    stf(i).SAD           = SAD;
+    stf(i).isoCenter     = pln.propStf.isoCenter(i,:);
+        
+    % Get the (active) rotation matrix. We perform a passive/system 
+    % rotation with row vector coordinates, which would introduce two 
+    % inversions / transpositions of the matrix, thus no changes to the
+    % rotation matrix are necessary
+    rotMat_system_T = matRad_getRotationMatrix(pln.propStf.gantryAngles(i),pln.propStf.couchAngles(i));
+    
+    rot_coords = [coordsX coordsY coordsZ]*rotMat_system_T;
+    
+    % project x and z coordinates to isocenter
+    coordsAtIsoCenterPlane(:,1) = (rot_coords(:,1)*SAD)./(SAD + rot_coords(:,2));
+    coordsAtIsoCenterPlane(:,2) = (rot_coords(:,3)*SAD)./(SAD + rot_coords(:,2));
+    
+    % Take unique rows values for beamlets positions. Calculate position of
+    % central ray for every bixel    
+    rayPos = unique(pln.propStf.bixelWidth*round([           coordsAtIsoCenterPlane(:,1) ... 
+                                                  zeros(size(coordsAtIsoCenterPlane,1),1) ...
+                                                             coordsAtIsoCenterPlane(:,2)]/pln.propStf.bixelWidth),'rows');
+                                                         
+    disp(rayPos);
+                                                  
+    % pad ray position array if resolution of target voxel grid not sufficient
+    maxCtResolution = max([ct.resolution.x ct.resolution.y ct.resolution.z]);
+    if pln.propStf.bixelWidth < maxCtResolution
+        origRayPos = rayPos;
+        for j = -floor(maxCtResolution/pln.propStf.bixelWidth):floor(maxCtResolution/pln.propStf.bixelWidth)
+            for k = -floor(maxCtResolution/pln.propStf.bixelWidth):floor(maxCtResolution/pln.propStf.bixelWidth)
+                if abs(j)+abs(k)==0
+                    continue;
+                end                
+                rayPos = [rayPos; origRayPos(:,1)+j*pln.propStf.bixelWidth origRayPos(:,2) origRayPos(:,3)+k*pln.propStf.bixelWidth];
+            end
+        end
+     end
+
+     % remove spaces within rows of bixels for DAO
+     if pln.propOpt.runDAO
+         % create single x,y,z vectors
+         x = rayPos(:,1);
+         y = rayPos(:,2);
+         z = rayPos(:,3);
+         uniZ = unique(z);
+         for j = 1:numel(uniZ)
+             x_loc = x(z == uniZ(j));
+             x_min = min(x_loc);
+             x_max = max(x_loc);
+             x = [x; [x_min:pln.propStf.bixelWidth:x_max]'];
+             y = [y; zeros((x_max-x_min)/pln.propStf.bixelWidth+1,1)];
+             z = [z; uniZ(j)*ones((x_max-x_min)/pln.propStf.bixelWidth+1,1)];             
+         end
+         
+         rayPos = [x,y,z];
+     end
+     
+    % remove double rays
+    rayPos = unique(rayPos,'rows');
+    
+    % Save the number of rays
+    stf(i).numOfRays = size(rayPos,1);
+    
+    % Save ray and target position in beam eye's view (bev)
+    for j = 1:stf(i).numOfRays
+        stf(i).ray(j).rayPos_bev = rayPos(j,:);
+        stf(i).ray(j).targetPoint_bev = [2*stf(i).ray(j).rayPos_bev(1) ...
+                                                               SAD ...
+                                         2*stf(i).ray(j).rayPos_bev(3)];
+    end
+    
+    % source position in bev
+    stf(i).sourcePoint_bev = [0 -SAD 0];
+    
+    % get (active) rotation matrix 
+    % transpose matrix because we are working with row vectors
+    rotMat_vectors_T = transpose(matRad_getRotationMatrix(pln.propStf.gantryAngles(i),pln.propStf.couchAngles(i)));
+    
+    
+    stf(i).sourcePoint = stf(i).sourcePoint_bev*rotMat_vectors_T;
+    
+    % Save ray and target position in lps system.
+    for j = 1:stf(i).numOfRays
+        stf(i).ray(j).rayPos      = stf(i).ray(j).rayPos_bev*rotMat_vectors_T;
+        stf(i).ray(j).targetPoint = stf(i).ray(j).targetPoint_bev*rotMat_vectors_T;
+        if strcmp(pln.radiationMode,'photons') 
+            stf(i).ray(j).beamletCornersAtIso = [rayPos(j,:) + [+stf(i).bixelWidth/2,0,+stf(i).bixelWidth/2];...
+                                                 rayPos(j,:) + [-stf(i).bixelWidth/2,0,+stf(i).bixelWidth/2];...
+                                                 rayPos(j,:) + [-stf(i).bixelWidth/2,0,-stf(i).bixelWidth/2];...
+                                                 rayPos(j,:) + [+stf(i).bixelWidth/2,0,-stf(i).bixelWidth/2]]*rotMat_vectors_T;
+            stf(i).ray(j).rayCorners_SCD = (repmat([0, machine.meta.SCD - SAD, 0],4,1)+ (machine.meta.SCD/SAD) * ...
+                                                             [rayPos(j,:) + [+stf(i).bixelWidth/2,0,+stf(i).bixelWidth/2];...
+                                                              rayPos(j,:) + [-stf(i).bixelWidth/2,0,+stf(i).bixelWidth/2];...
+                                                              rayPos(j,:) + [-stf(i).bixelWidth/2,0,-stf(i).bixelWidth/2];...
+                                                              rayPos(j,:) + [+stf(i).bixelWidth/2,0,-stf(i).bixelWidth/2]])*rotMat_vectors_T;
+        end
+    end
+    
+    % loop over all rays to determine meta information for each ray    
+    stf(i).numOfBixelsPerRay = ones(1,stf(i).numOfRays);
+    
+    for j = stf(i).numOfRays:-1:1
+
+        % ray tracing necessary to determine depth of the target
+        [~,l,rho,~,~] = matRad_siddonRayTracer(stf(i).isoCenter, ...
+                             ct.resolution, ...
+                             stf(i).sourcePoint, ...
+                             stf(i).ray(j).targetPoint, ...
+                             [{ct.cube{1}} {voiTarget}]);
+
+        % find appropriate energies for particles
+       if strcmp(stf(i).radiationMode,'protons') || strcmp(stf(i).radiationMode,'carbon')
+
+           % target hit
+           if sum(rho{2}) > 0 
+
+                % compute radiological depths
+                % http://www.ncbi.nlm.nih.gov/pubmed/4000088, eq 14
+                radDepths = cumsum(l .* rho{1}); 
+
+                % find target entry & exit
+                diff_voi    = diff([rho{2}]);
+                targetEntry = radDepths(diff_voi == 1);
+                targetExit  = radDepths(diff_voi == -1);
+
+                if numel(targetEntry) ~= numel(targetExit)
+                    matRad_cfg.dispError('Inconsistency during ray tracing. Please check correct assignment and overlap priorities of structure types OAR & TARGET.');
+                end
+
+                stf(i).ray(j).energy = [];
+
+                % Save energies in stf struct
+                for k = 1:numel(targetEntry)
+                    stf(i).ray(j).energy = [stf(i).ray(j).energy availableEnergies(availablePeakPos>=targetEntry(k)&availablePeakPos<=targetExit(k))];
+                end
+  
+                % book keeping & calculate focus index
+                stf(i).numOfBixelsPerRay(j) = numel([stf(i).ray(j).energy]);
+                currentMinimumFWHM = matRad_interp1(machine.meta.LUT_bxWidthminFWHM(1,:)',...
+                                             machine.meta.LUT_bxWidthminFWHM(2,:)',...
+                                             pln.propStf.bixelWidth);
+                focusIx  =  ones(stf(i).numOfBixelsPerRay(j),1);
+                [~, vEnergyIx] = min(abs(bsxfun(@minus,[machine.data.energy]',...
+                                repmat(stf(i).ray(j).energy,length([machine.data]),1))));
+
+                % get for each spot the focus index
+                for k = 1:stf(i).numOfBixelsPerRay(j)                    
+                    focusIx(k) = find(machine.data(vEnergyIx(k)).initFocus.SisFWHMAtIso > currentMinimumFWHM,1,'first');
+                end
+
+                stf(i).ray(j).focusIx = focusIx';
+                 
+            else % target not hit
+                stf(i).ray(j)               = [];
+                stf(i).numOfBixelsPerRay(j) = [];
+           end
+           
+       elseif strcmp(stf(i).radiationMode,'photons')
+           
+         % book keeping for photons
+         stf(i).ray(j).energy = machine.data.energy;
+         
+       else
+          matRad_cfg.dispError('Error generating stf struct: invalid radiation modality.');
+       end
+       
+    end
+    
+    % store total number of rays for beam-i
+    stf(i).numOfRays = size(stf(i).ray,2);
+     
+    % post processing for particle remove energy slices
+    if strcmp(stf(i).radiationMode,'protons') || strcmp(stf(i).radiationMode,'carbon')
+        
+        % get minimum energy per field
+        minEnergy = min([stf(i).ray.energy]);
+        maxEnergy = max([stf(i).ray.energy]);
+        
+        % get corresponding peak position
+        availableEnergies = [machine.data.energy];
+        minPeakPos  = machine.data(minEnergy == availableEnergies).peakPos;
+        maxPeakPos  = machine.data(maxEnergy == availableEnergies).peakPos;
+        
+        % find set of energyies with adequate spacing
+        if ~isfield(pln.propStf, 'longitudinalSpotSpacing')
+            longitudinalSpotSpacing = matRad_cfg.propStf.defaultLongitudinalSpotSpacing;
+        else
+            longitudinalSpotSpacing = pln.propStf.longitudinalSpotSpacing;
+        end
+        
+        stf(i).longitudinalSpotSpacing = longitudinalSpotSpacing;
+        
+        tolerance              = longitudinalSpotSpacing/10;
+        availablePeakPos       = [machine.data.peakPos];
+        
+        useEnergyBool = availablePeakPos >= minPeakPos & availablePeakPos <= maxPeakPos;
+        
+        ixCurr = find(useEnergyBool,1,'first');
+        ixRun  = ixCurr + 1;
+        ixEnd  = find(useEnergyBool,1,'last');
+
+        while ixRun <= ixEnd
+            if abs(availablePeakPos(ixRun)-availablePeakPos(ixCurr)) < ...
+                                    longitudinalSpotSpacing - tolerance
+                useEnergyBool(ixRun) = 0;
+            else
+                ixCurr = ixRun;
+            end
+            ixRun = ixRun + 1;
+        end
+        
+        for j = stf(i).numOfRays:-1:1
+            for k = stf(i).numOfBixelsPerRay(j):-1:1
+                maskEnergy = stf(i).ray(j).energy(k) == availableEnergies;
+                if ~useEnergyBool(maskEnergy)
+                    stf(i).ray(j).energy(k)     = [];
+                    stf(i).ray(j).focusIx(k)    = [];
+                    stf(i).numOfBixelsPerRay(j) = stf(i).numOfBixelsPerRay(j) - 1;
+                end
+            end
+            if isempty(stf(i).ray(j).energy)
+                stf(i).ray(j) = [];
+                stf(i).numOfBixelsPerRay(j) = [];
+                stf(i).numOfRays = stf(i).numOfRays - 1;
+            end
+        end
+        
+    end
+    
+    % save total number of bixels
+    stf(i).totalNumOfBixels = sum(stf(i).numOfBixelsPerRay);
+    
+    % Show progress
+    matRad_progress(i,length(pln.propStf.gantryAngles));
+
+    %% visualization
+    if visMode > 0
+        
+        clf;
+        % first subplot: visualization in bev
+        subplot(1,2,1)
+        hold on
+        
+        % plot rotated target coordinates
+        plot3(rot_coords(:,1),rot_coords(:,2),rot_coords(:,3),'r.')
+        
+        % surface rendering
+        if visMode == 2
+            
+            % generate a 3D rectangular grid centered at isocenter in
+            % voxel coordinates
+            [X,Y,Z] = meshgrid((1:ct.cubeDim(2))-stf(i).isoCenter(1)/ct.resolution.x, ...
+                               (1:ct.cubeDim(1))-stf(i).isoCenter(2)/ct.resolution.y, ...
+                               (1:ct.cubeDim(3))-stf(i).isoCenter(3)/ct.resolution.z);
+            
+            % computes surface
+            patSurfCube      = 0*ct.cube{1};
+            idx              = [cst{:,4}];
+            idx              = unique(vertcat(idx{:}));
+            patSurfCube(idx) = 1;
+            
+            [f,v] = isosurface(X,Y,Z,patSurfCube,.5);
+            
+            % convert isosurface from voxel to [mm]
+            v(:,1) = v(:,1)*ct.resolution.x;
+            v(:,2) = v(:,2)*ct.resolution.y;
+            v(:,3) = v(:,3)*ct.resolution.z;
+            
+            % rotate surface
+            rotated_surface = v*rotMat_system_T;
+            
+            % surface rendering
+            surface = patch('Faces',f,'Vertices',rotated_surface);
+            set(surface,'FaceColor',[0 0 1],'EdgeColor','none','FaceAlpha',.4);
+            lighting gouraud;
+        
+        end
+        
+        % plot projection matrix: coordinates at isocenter
+        plot3(rayPos(:,1),rayPos(:,2),rayPos(:,3),'k.');
+        
+        % Plot matrix border of matrix at isocenter
+        for j = 1:stf(i).numOfRays
+            
+            % Compute border for every bixels
+            targetPoint_vox_X_1 = stf(i).ray(j).targetPoint_bev(:,1) + pln.propStf.bixelWidth;
+            targetPoint_vox_Y_1 = stf(i).ray(j).targetPoint_bev(:,2);
+            targetPoint_vox_Z_1 = stf(i).ray(j).targetPoint_bev(:,3) + pln.propStf.bixelWidth;
+            
+            targetPoint_vox_X_2 = stf(i).ray(j).targetPoint_bev(:,1) + pln.propStf.bixelWidth;
+            targetPoint_vox_Y_2 = stf(i).ray(j).targetPoint_bev(:,2);
+            targetPoint_vox_Z_2 = stf(i).ray(j).targetPoint_bev(:,3) - pln.propStf.bixelWidth;
+            
+            targetPoint_vox_X_3 = stf(i).ray(j).targetPoint_bev(:,1) - pln.propStf.bixelWidth;
+            targetPoint_vox_Y_3 = stf(i).ray(j).targetPoint_bev(:,2);
+            targetPoint_vox_Z_3 = stf(i).ray(j).targetPoint_bev(:,3) - pln.propStf.bixelWidth;
+            
+            targetPoint_vox_X_4 = stf(i).ray(j).targetPoint_bev(:,1) - pln.propStf.bixelWidth;
+            targetPoint_vox_Y_4 = stf(i).ray(j).targetPoint_bev(:,2);
+            targetPoint_vox_Z_4 = stf(i).ray(j).targetPoint_bev(:,3) + pln.propStf.bixelWidth;
+            
+            % plot
+            plot3([stf(i).sourcePoint_bev(1) targetPoint_vox_X_1],[stf(i).sourcePoint_bev(2) targetPoint_vox_Y_1],[stf(i).sourcePoint_bev(3) targetPoint_vox_Z_1],'g')
+            plot3([stf(i).sourcePoint_bev(1) targetPoint_vox_X_2],[stf(i).sourcePoint_bev(2) targetPoint_vox_Y_2],[stf(i).sourcePoint_bev(3) targetPoint_vox_Z_2],'g')
+            plot3([stf(i).sourcePoint_bev(1) targetPoint_vox_X_3],[stf(i).sourcePoint_bev(2) targetPoint_vox_Y_3],[stf(i).sourcePoint_bev(3) targetPoint_vox_Z_3],'g')
+            plot3([stf(i).sourcePoint_bev(1) targetPoint_vox_X_4],[stf(i).sourcePoint_bev(2) targetPoint_vox_Y_4],[stf(i).sourcePoint_bev(3) targetPoint_vox_Z_4],'g')
+            
+        end
+        
+        % Plot properties
+        daspect([1 1 1]);
+        view(0,-90);
+        xlabel 'X [mm]'
+        ylabel 'Y [mm]'
+        zlabel 'Z [mm]'
+        title ('Beam''s eye view')
+        axis([-300 300 -300 300 -300 300]);
+        
+        % second subplot: visualization in lps coordinate system
+        subplot(1,2,2)
+        
+        % Plot target coordinates whitout any rotation
+        plot3(coordsX,coordsY,coordsZ,'r.')
+        hold on;
+        
+        % Rotated projection matrix at isocenter
+        isocenter_plane_coor = rayPos*rotMat_vectors_T;
+        
+        % Plot isocenter plane
+        plot3(isocenter_plane_coor(:,1),isocenter_plane_coor(:,2),isocenter_plane_coor(:,3),'y.');
+        
+        % Plot rotated bixels border.
+        for j = 1:stf(i).numOfRays
+            % Generate rotated projection target points.
+            targetPoint_vox_1_rotated = [stf(i).ray(j).targetPoint_bev(:,1) + pln.propStf.bixelWidth,stf(i).ray(j).targetPoint_bev(:,2),stf(i).ray(j).targetPoint_bev(:,3) + pln.propStf.bixelWidth]*rotMat_vectors_T;
+            targetPoint_vox_2_rotated = [stf(i).ray(j).targetPoint_bev(:,1) + pln.propStf.bixelWidth,stf(i).ray(j).targetPoint_bev(:,2),stf(i).ray(j).targetPoint_bev(:,3) - pln.propStf.bixelWidth]*rotMat_vectors_T;
+            targetPoint_vox_3_rotated = [stf(i).ray(j).targetPoint_bev(:,1) - pln.propStf.bixelWidth,stf(i).ray(j).targetPoint_bev(:,2),stf(i).ray(j).targetPoint_bev(:,3) - pln.propStf.bixelWidth]*rotMat_vectors_T;
+            targetPoint_vox_4_rotated = [stf(i).ray(j).targetPoint_bev(:,1) - pln.propStf.bixelWidth,stf(i).ray(j).targetPoint_bev(:,2),stf(i).ray(j).targetPoint_bev(:,3) + pln.propStf.bixelWidth]*rotMat_vectors_T;
+            
+            % Plot rotated target points.
+            plot3([stf(i).sourcePoint(1) targetPoint_vox_1_rotated(:,1)],[stf(i).sourcePoint(2) targetPoint_vox_1_rotated(:,2)],[stf(i).sourcePoint(3) targetPoint_vox_1_rotated(:,3)],'g')
+            plot3([stf(i).sourcePoint(1) targetPoint_vox_2_rotated(:,1)],[stf(i).sourcePoint(2) targetPoint_vox_2_rotated(:,2)],[stf(i).sourcePoint(3) targetPoint_vox_2_rotated(:,3)],'g')
+            plot3([stf(i).sourcePoint(1) targetPoint_vox_3_rotated(:,1)],[stf(i).sourcePoint(2) targetPoint_vox_3_rotated(:,2)],[stf(i).sourcePoint(3) targetPoint_vox_3_rotated(:,3)],'g')
+            plot3([stf(i).sourcePoint(1) targetPoint_vox_4_rotated(:,1)],[stf(i).sourcePoint(2) targetPoint_vox_4_rotated(:,2)],[stf(i).sourcePoint(3) targetPoint_vox_4_rotated(:,3)],'g')
+        end
+        
+        % surface rendering
+        if visMode == 2
+            surface = patch('Faces',f,'Vertices',v);
+            set(surface,'FaceColor',[0 0 1],'EdgeColor','none','FaceAlpha',.4);
+            lighting gouraud;
+        end
+        
+        % labels etc.
+        daspect([1 1 1]);
+        view(0,-90);
+        xlabel 'X [mm]'
+        ylabel 'Y [mm]'
+        zlabel 'Z [mm]'
+        title 'lps coordinate system'
+        axis([-300 300 -300 300 -300 300]);
+        %pause(1);
+    end
+    
+    % include rangeshifter data if not yet available 
+    if strcmp(pln.radiationMode, 'protons') || strcmp(pln.radiationMode, 'carbon')
+        for j = 1:stf(i).numOfRays
+            for k = 1:numel(stf(i).ray(j).energy)
+                stf(i).ray(j).rangeShifter(k).ID = 0;
+                stf(i).ray(j).rangeShifter(k).eqThickness = 0;
+                stf(i).ray(j).rangeShifter(k).sourceRashiDistance = 0;
+            end
+        end
+    end
+        
+end    
+
+end
diff --git a/matRad_generateStf.m b/matRad_generateStf.m
index f128f0db0..d3059e1f5 100644
--- a/matRad_generateStf.m
+++ b/matRad_generateStf.m
@@ -53,6 +53,8 @@
     end
 end
 
+disp(length(V));
+
 % Remove double voxels
 V = unique(V);
 % generate voi cube for targets
@@ -115,6 +117,10 @@
 % Define steering file like struct. Prellocating for speed.
 stf = struct;
 
+save('coordsX_vox', 'coordsX_vox');
+save('coordsY_vox', 'coordsY_vox');
+save('coordsZ_vox', 'coordsZ_vox');
+
 % loop over all angles
 for i = 1:length(pln.propStf.gantryAngles)
     
@@ -123,6 +129,8 @@
     coordsX = coordsX_vox*ct.resolution.x - pln.propStf.isoCenter(i,1);
     coordsY = coordsY_vox*ct.resolution.y - pln.propStf.isoCenter(i,2);
     coordsZ = coordsZ_vox*ct.resolution.z - pln.propStf.isoCenter(i,3);
+    
+    save('coordsX', 'coordsX');
 
     % Save meta information for treatment plan
     stf(i).gantryAngle   = pln.propStf.gantryAngles(i);
@@ -182,7 +190,7 @@
          
          rayPos = [x,y,z];
      end
-    
+     
     % remove double rays
     rayPos = unique(rayPos,'rows');
     

From 92c34d36e7930cac3e9ae62c691817675a0fc9ce Mon Sep 17 00:00:00 2001
From: g-stanic <goran.stanic@dkfz-heidelberg.de>
Date: Thu, 31 Aug 2023 10:00:53 +0200
Subject: [PATCH 3/4] Update of IO scripts

---
 IO_stf.asv                           |  21 --
 IO_stf.m                             |   8 +-
 examples/matRad_example2_photons.m   |   6 +-
 examples/matRad_example4_photonsMC.m |   2 +-
 matRad_callFromPython.m              |  14 +-
 matRad_generateStf.asv               | 516 ---------------------------
 matRad_saveStructs.m                 |  12 +-
 7 files changed, 23 insertions(+), 556 deletions(-)
 delete mode 100644 IO_stf.asv
 delete mode 100644 matRad_generateStf.asv

diff --git a/IO_stf.asv b/IO_stf.asv
deleted file mode 100644
index 5593dfe55..000000000
--- a/IO_stf.asv
+++ /dev/null
@@ -1,21 +0,0 @@
-clear all;
-matRad_rc;
-
-%%
-stf_init = load('C:\Users\g906c\Projects\pyRadPlan\stf_with_separate_rays.mat', 'stf');
-ray_init = load('C:\Users\g906c\Projects\pyRadPlan\stf_with_separate_rays.mat', 'rays');
-
-%%
-stf = [stf_init.stf{:}];
-ray = cell(size(ray_init.rays, 2), 1);
-
-for i=1:size(ray_init.rays, 2)
-    ray{i} = [ray_init.rays{i}{:}];
-end
-
-for i=1:size(ray_init.rays, 2)
-    stf(i).ray = ray{i};
-end
-
-%%
-save('stf.m')
\ No newline at end of file
diff --git a/IO_stf.m b/IO_stf.m
index 665e7a25f..9a8c1c8be 100644
--- a/IO_stf.m
+++ b/IO_stf.m
@@ -1,7 +1,7 @@
-function status = IO_stf()
+function status = IO_stf(path)
 
-stf_init = load('C:\Users\g906c\Projects\pyRadPlan\stf_with_separate_rays.mat', 'stf');
-ray_init = load('C:\Users\g906c\Projects\pyRadPlan\stf_with_separate_rays.mat', 'rays');
+stf_init = load(strcat(path,'stf_with_separate_rays.mat'), 'stf');
+ray_init = load(strcat(path,'stf_with_separate_rays.mat'), 'rays');
 
 %%
 stf = [stf_init.stf{:}];
@@ -16,7 +16,7 @@
 end
 
 %%
-save('C:\Users\g906c\Projects\pyRadPlan\stf.mat', 'stf')
+save(strcat(path,'stf.mat'), 'stf')
 
 status = 'STF written';
 
diff --git a/examples/matRad_example2_photons.m b/examples/matRad_example2_photons.m
index 328276735..5c637f1a1 100644
--- a/examples/matRad_example2_photons.m
+++ b/examples/matRad_example2_photons.m
@@ -106,7 +106,7 @@
 % to 30. Internally, matRad considers the fraction dose for optimization, 
 % however, objetives and constraints are defined for the entire treatment.
 pln.numOfFractions         = 30;
-pln.propStf.gantryAngles   = [0:40:359];
+pln.propStf.gantryAngles   = [0, 90, 180, 270];
 pln.propStf.couchAngles    = zeros(1,numel(pln.propStf.gantryAngles));
 pln.propStf.bixelWidth     = 5;
 
@@ -146,7 +146,7 @@
 % Let's generate dosimetric information by pre-computing dose influence 
 % matrices for unit beamlet intensities. Having dose influences available 
 % allows subsequent inverse optimization.
-dij = matRad_calcPhotonDose(ct,stf,pln,cst);
+dij = matRad_calcPhotonDoseMC(ct,stf,pln,cst,1000);
 
 %% Inverse Optimization for IMRT
 % The goal of the fluence optimization is to find a set of beamlet/pencil 
@@ -155,7 +155,7 @@
 % treatment. Once the optimization has finished, trigger once the GUI to 
 % visualize the optimized dose cubes.
 resultGUI = matRad_fluenceOptimization(dij,cst,pln);
-matRadGUI;
+%matRadGUI;
 
 %% Plot the Resulting Dose Slice
 % Let's plot the transversal iso-center dose slice
diff --git a/examples/matRad_example4_photonsMC.m b/examples/matRad_example4_photonsMC.m
index e6e34fbe3..aa5afa7ab 100644
--- a/examples/matRad_example4_photonsMC.m
+++ b/examples/matRad_example4_photonsMC.m
@@ -56,7 +56,7 @@
 %% Dose Calculation
 % Calculate dose influence matrix for unit pencil beam intensities using 
 % a Monte Carlo algorithm
-dij = matRad_calcPhotonDoseMC(ct,stf,pln,cst);
+dij = matRad_calcPhotonDoseMC(ct,stf,pln,cst,1000);
 
 %% Inverse Optimization for IMRT
 resultGUI = matRad_fluenceOptimization(dij,cst,pln);
diff --git a/matRad_callFromPython.m b/matRad_callFromPython.m
index 9bd8b5f65..f2f37f75c 100644
--- a/matRad_callFromPython.m
+++ b/matRad_callFromPython.m
@@ -1,10 +1,14 @@
-function matRad_callFromPython(functionName, outputName, varargin)
+function matRad_callFromPython(functionName, outputName, inputPath, outputPath, varargin)
 %matRad_callFromPython Function that uses temporary mat file to call any function from within python
 
 for i=1:length(varargin)
-    load(varargin{i});
-    var=varargin{i}(1:end-4);
-    functionVars{i}=var;
+    if contains(string(varargin{i}), string('.mat'))
+        load(strcat(inputPath, varargin{i}));
+        [path, var, ext]=fileparts(varargin{i});
+        functionVars{i}=var;
+    else
+        functionVars{i} = num2str(varargin{i});
+    end
 end
 
 execFunc = sprintf('%s = %s(%s);', outputName, functionName, strjoin(functionVars,','));
@@ -14,6 +18,6 @@ function matRad_callFromPython(functionName, outputName, varargin)
 %end
 
 eval(execFunc);
-save(strcat(outputName,'.mat'), outputName);
+save(strcat(outputPath, outputName,'.mat'), outputName);
 
 end
diff --git a/matRad_generateStf.asv b/matRad_generateStf.asv
deleted file mode 100644
index e4ce32668..000000000
--- a/matRad_generateStf.asv
+++ /dev/null
@@ -1,516 +0,0 @@
-function stf = matRad_generateStf(ct,cst,pln,visMode)
-% matRad steering information generation
-% 
-% call
-%   stf = matRad_generateStf(ct,cst,pln,visMode)
-%
-% input
-%   ct:         ct cube
-%   cst:        matRad cst struct
-%   pln:        matRad plan meta information struct
-%   visMode:    toggle on/off different visualizations by setting this value to 1,2,3 (optional)
-%
-% output
-%   stf:        matRad steering information struct
-%
-% References
-%   -
-%
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%
-% Copyright 2015 the matRad development team. 
-% 
-% This file is part of the matRad project. It is subject to the license 
-% terms in the LICENSE file found in the top-level directory of this 
-% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part 
-% of the matRad project, including this file, may be copied, modified, 
-% propagated, or distributed except according to the terms contained in the 
-% LICENSE file.
-%
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-matRad_cfg = MatRad_Config.instance();
-
-matRad_cfg.dispInfo('matRad: Generating stf struct... ');
-
-if nargin < 4
-    visMode = 0;
-end
-
-if numel(pln.propStf.gantryAngles) ~= numel(pln.propStf.couchAngles)
-    matRad_cfg.dispError('Inconsistent number of gantry and couch angles.');
-end
-
-if ~isnumeric(pln.propStf.bixelWidth) || pln.propStf.bixelWidth < 0 || ~isfinite(pln.propStf.bixelWidth)
-   matRad_cfg.dispError('bixel width (spot distance) needs to be a real number [mm] larger than zero.');
-end
-
-% find all target voxels from cst cell array
-V = [];
-for i=1:size(cst,1)
-    if isequal(cst{i,3},'TARGET') && ~isempty(cst{i,6})
-        V = [V;vertcat(cst{i,4}{:})];
-    end
-end
-
-disp(length(V));
-
-% Remove double voxels
-V = unique(V);
-% generate voi cube for targets
-voiTarget    = zeros(ct.cubeDim);
-voiTarget(V) = 1;
-    
-% add margin
-addmarginBool = matRad_cfg.propStf.defaultAddMargin;
-if isfield(pln,'propStf') && isfield(pln.propStf,'addMargin')
-   addmarginBool = pln.propStf.addMargin; 
-end
-
-if addmarginBool
-    voiTarget = matRad_addMargin(voiTarget,cst,ct.resolution,ct.resolution,true);
-    V   = find(voiTarget>0);
-end
-
-% throw error message if no target is found
-if isempty(V)
-    matRad_cfg.dispError('Could not find target.');
-end
-
-% Convert linear indices to 3D voxel coordinates
-[coordsY_vox, coordsX_vox, coordsZ_vox] = ind2sub(ct.cubeDim,V);
-
-% prepare structures necessary for particles
-fileName = [pln.radiationMode '_' pln.machine];
-try
-   load([fileparts(mfilename('fullpath')) filesep 'basedata' filesep fileName]);
-   SAD = machine.meta.SAD;
-catch
-   matRad_cfg.dispError('Could not find the following machine file: %s',fileName); 
-end
-
-if strcmp(pln.radiationMode,'protons') || strcmp(pln.radiationMode,'carbon')
-      
-    availableEnergies = [machine.data.energy];
-    availablePeakPos  = [machine.data.peakPos] + [machine.data.offset];
-    
-    if sum(availablePeakPos<0)>0
-       matRad_cfg.dispError('at least one available peak position is negative - inconsistent machine file') 
-    end
-    %clear machine;
-end
-
-% calculate rED or rSP from HU
-ct = matRad_calcWaterEqD(ct, pln);
-
-% take only voxels inside patient
-V = [cst{:,4}];
-V = unique(vertcat(V{:}));
-
-% ignore densities outside of contours
-eraseCtDensMask = ones(prod(ct.cubeDim),1);
-eraseCtDensMask(V) = 0;
-for i = 1:ct.numOfCtScen
-    ct.cube{i}(eraseCtDensMask == 1) = 0;
-end
-
-% Define steering file like struct. Prellocating for speed.
-stf = struct;
-
-save('coordsX_vox', 'coordsX_vox');
-save('coordsY_vox', 'coordsY_vox');
-save('coordsZ_vox', 'coordsZ_vox');
-
-% loop over all angles
-for i = 1:length(pln.propStf.gantryAngles)
-    
-    % Correct for iso center position. Whit this correction Isocenter is
-    % (0,0,0) [mm]
-    coordsX = coordsX_vox*ct.resolution.x - pln.propStf.isoCenter(i,1);
-    coordsY = coordsY_vox*ct.resolution.y - pln.propStf.isoCenter(i,2);
-    coordsZ = coordsZ_vox*ct.resolution.z - pln.propStf.isoCenter(i,3);
-    
-    save('coordsX', 'coordsX');
-
-    % Save meta information for treatment plan
-    stf(i).gantryAngle   = pln.propStf.gantryAngles(i);
-    stf(i).couchAngle    = pln.propStf.couchAngles(i);
-    stf(i).bixelWidth    = pln.propStf.bixelWidth;
-    stf(i).radiationMode = pln.radiationMode;
-    stf(i).SAD           = SAD;
-    stf(i).isoCenter     = pln.propStf.isoCenter(i,:);
-        
-    % Get the (active) rotation matrix. We perform a passive/system 
-    % rotation with row vector coordinates, which would introduce two 
-    % inversions / transpositions of the matrix, thus no changes to the
-    % rotation matrix are necessary
-    rotMat_system_T = matRad_getRotationMatrix(pln.propStf.gantryAngles(i),pln.propStf.couchAngles(i));
-    
-    rot_coords = [coordsX coordsY coordsZ]*rotMat_system_T;
-    
-    % project x and z coordinates to isocenter
-    coordsAtIsoCenterPlane(:,1) = (rot_coords(:,1)*SAD)./(SAD + rot_coords(:,2));
-    coordsAtIsoCenterPlane(:,2) = (rot_coords(:,3)*SAD)./(SAD + rot_coords(:,2));
-    
-    % Take unique rows values for beamlets positions. Calculate position of
-    % central ray for every bixel    
-    rayPos = unique(pln.propStf.bixelWidth*round([           coordsAtIsoCenterPlane(:,1) ... 
-                                                  zeros(size(coordsAtIsoCenterPlane,1),1) ...
-                                                             coordsAtIsoCenterPlane(:,2)]/pln.propStf.bixelWidth),'rows');
-                                                         
-    disp(rayPos);
-                                                  
-    % pad ray position array if resolution of target voxel grid not sufficient
-    maxCtResolution = max([ct.resolution.x ct.resolution.y ct.resolution.z]);
-    if pln.propStf.bixelWidth < maxCtResolution
-        origRayPos = rayPos;
-        for j = -floor(maxCtResolution/pln.propStf.bixelWidth):floor(maxCtResolution/pln.propStf.bixelWidth)
-            for k = -floor(maxCtResolution/pln.propStf.bixelWidth):floor(maxCtResolution/pln.propStf.bixelWidth)
-                if abs(j)+abs(k)==0
-                    continue;
-                end                
-                rayPos = [rayPos; origRayPos(:,1)+j*pln.propStf.bixelWidth origRayPos(:,2) origRayPos(:,3)+k*pln.propStf.bixelWidth];
-            end
-        end
-     end
-
-     % remove spaces within rows of bixels for DAO
-     if pln.propOpt.runDAO
-         % create single x,y,z vectors
-         x = rayPos(:,1);
-         y = rayPos(:,2);
-         z = rayPos(:,3);
-         uniZ = unique(z);
-         for j = 1:numel(uniZ)
-             x_loc = x(z == uniZ(j));
-             x_min = min(x_loc);
-             x_max = max(x_loc);
-             x = [x; [x_min:pln.propStf.bixelWidth:x_max]'];
-             y = [y; zeros((x_max-x_min)/pln.propStf.bixelWidth+1,1)];
-             z = [z; uniZ(j)*ones((x_max-x_min)/pln.propStf.bixelWidth+1,1)];             
-         end
-         
-         rayPos = [x,y,z];
-     end
-     
-    % remove double rays
-    rayPos = unique(rayPos,'rows');
-    
-    % Save the number of rays
-    stf(i).numOfRays = size(rayPos,1);
-    
-    % Save ray and target position in beam eye's view (bev)
-    for j = 1:stf(i).numOfRays
-        stf(i).ray(j).rayPos_bev = rayPos(j,:);
-        stf(i).ray(j).targetPoint_bev = [2*stf(i).ray(j).rayPos_bev(1) ...
-                                                               SAD ...
-                                         2*stf(i).ray(j).rayPos_bev(3)];
-    end
-    
-    % source position in bev
-    stf(i).sourcePoint_bev = [0 -SAD 0];
-    
-    % get (active) rotation matrix 
-    % transpose matrix because we are working with row vectors
-    rotMat_vectors_T = transpose(matRad_getRotationMatrix(pln.propStf.gantryAngles(i),pln.propStf.couchAngles(i)));
-    
-    
-    stf(i).sourcePoint = stf(i).sourcePoint_bev*rotMat_vectors_T;
-    
-    % Save ray and target position in lps system.
-    for j = 1:stf(i).numOfRays
-        stf(i).ray(j).rayPos      = stf(i).ray(j).rayPos_bev*rotMat_vectors_T;
-        stf(i).ray(j).targetPoint = stf(i).ray(j).targetPoint_bev*rotMat_vectors_T;
-        if strcmp(pln.radiationMode,'photons') 
-            stf(i).ray(j).beamletCornersAtIso = [rayPos(j,:) + [+stf(i).bixelWidth/2,0,+stf(i).bixelWidth/2];...
-                                                 rayPos(j,:) + [-stf(i).bixelWidth/2,0,+stf(i).bixelWidth/2];...
-                                                 rayPos(j,:) + [-stf(i).bixelWidth/2,0,-stf(i).bixelWidth/2];...
-                                                 rayPos(j,:) + [+stf(i).bixelWidth/2,0,-stf(i).bixelWidth/2]]*rotMat_vectors_T;
-            stf(i).ray(j).rayCorners_SCD = (repmat([0, machine.meta.SCD - SAD, 0],4,1)+ (machine.meta.SCD/SAD) * ...
-                                                             [rayPos(j,:) + [+stf(i).bixelWidth/2,0,+stf(i).bixelWidth/2];...
-                                                              rayPos(j,:) + [-stf(i).bixelWidth/2,0,+stf(i).bixelWidth/2];...
-                                                              rayPos(j,:) + [-stf(i).bixelWidth/2,0,-stf(i).bixelWidth/2];...
-                                                              rayPos(j,:) + [+stf(i).bixelWidth/2,0,-stf(i).bixelWidth/2]])*rotMat_vectors_T;
-        end
-    end
-    
-    % loop over all rays to determine meta information for each ray    
-    stf(i).numOfBixelsPerRay = ones(1,stf(i).numOfRays);
-    
-    for j = stf(i).numOfRays:-1:1
-
-        % ray tracing necessary to determine depth of the target
-        [~,l,rho,~,~] = matRad_siddonRayTracer(stf(i).isoCenter, ...
-                             ct.resolution, ...
-                             stf(i).sourcePoint, ...
-                             stf(i).ray(j).targetPoint, ...
-                             [{ct.cube{1}} {voiTarget}]);
-
-        % find appropriate energies for particles
-       if strcmp(stf(i).radiationMode,'protons') || strcmp(stf(i).radiationMode,'carbon')
-
-           % target hit
-           if sum(rho{2}) > 0 
-
-                % compute radiological depths
-                % http://www.ncbi.nlm.nih.gov/pubmed/4000088, eq 14
-                radDepths = cumsum(l .* rho{1}); 
-
-                % find target entry & exit
-                diff_voi    = diff([rho{2}]);
-                targetEntry = radDepths(diff_voi == 1);
-                targetExit  = radDepths(diff_voi == -1);
-
-                if numel(targetEntry) ~= numel(targetExit)
-                    matRad_cfg.dispError('Inconsistency during ray tracing. Please check correct assignment and overlap priorities of structure types OAR & TARGET.');
-                end
-
-                stf(i).ray(j).energy = [];
-
-                % Save energies in stf struct
-                for k = 1:numel(targetEntry)
-                    stf(i).ray(j).energy = [stf(i).ray(j).energy availableEnergies(availablePeakPos>=targetEntry(k)&availablePeakPos<=targetExit(k))];
-                end
-  
-                % book keeping & calculate focus index
-                stf(i).numOfBixelsPerRay(j) = numel([stf(i).ray(j).energy]);
-                currentMinimumFWHM = matRad_interp1(machine.meta.LUT_bxWidthminFWHM(1,:)',...
-                                             machine.meta.LUT_bxWidthminFWHM(2,:)',...
-                                             pln.propStf.bixelWidth);
-                focusIx  =  ones(stf(i).numOfBixelsPerRay(j),1);
-                [~, vEnergyIx] = min(abs(bsxfun(@minus,[machine.data.energy]',...
-                                repmat(stf(i).ray(j).energy,length([machine.data]),1))));
-
-                % get for each spot the focus index
-                for k = 1:stf(i).numOfBixelsPerRay(j)                    
-                    focusIx(k) = find(machine.data(vEnergyIx(k)).initFocus.SisFWHMAtIso > currentMinimumFWHM,1,'first');
-                end
-
-                stf(i).ray(j).focusIx = focusIx';
-                 
-            else % target not hit
-                stf(i).ray(j)               = [];
-                stf(i).numOfBixelsPerRay(j) = [];
-           end
-           
-       elseif strcmp(stf(i).radiationMode,'photons')
-           
-         % book keeping for photons
-         stf(i).ray(j).energy = machine.data.energy;
-         
-       else
-          matRad_cfg.dispError('Error generating stf struct: invalid radiation modality.');
-       end
-       
-    end
-    
-    % store total number of rays for beam-i
-    stf(i).numOfRays = size(stf(i).ray,2);
-     
-    % post processing for particle remove energy slices
-    if strcmp(stf(i).radiationMode,'protons') || strcmp(stf(i).radiationMode,'carbon')
-        
-        % get minimum energy per field
-        minEnergy = min([stf(i).ray.energy]);
-        maxEnergy = max([stf(i).ray.energy]);
-        
-        % get corresponding peak position
-        availableEnergies = [machine.data.energy];
-        minPeakPos  = machine.data(minEnergy == availableEnergies).peakPos;
-        maxPeakPos  = machine.data(maxEnergy == availableEnergies).peakPos;
-        
-        % find set of energyies with adequate spacing
-        if ~isfield(pln.propStf, 'longitudinalSpotSpacing')
-            longitudinalSpotSpacing = matRad_cfg.propStf.defaultLongitudinalSpotSpacing;
-        else
-            longitudinalSpotSpacing = pln.propStf.longitudinalSpotSpacing;
-        end
-        
-        stf(i).longitudinalSpotSpacing = longitudinalSpotSpacing;
-        
-        tolerance              = longitudinalSpotSpacing/10;
-        availablePeakPos       = [machine.data.peakPos];
-        
-        useEnergyBool = availablePeakPos >= minPeakPos & availablePeakPos <= maxPeakPos;
-        
-        ixCurr = find(useEnergyBool,1,'first');
-        ixRun  = ixCurr + 1;
-        ixEnd  = find(useEnergyBool,1,'last');
-
-        while ixRun <= ixEnd
-            if abs(availablePeakPos(ixRun)-availablePeakPos(ixCurr)) < ...
-                                    longitudinalSpotSpacing - tolerance
-                useEnergyBool(ixRun) = 0;
-            else
-                ixCurr = ixRun;
-            end
-            ixRun = ixRun + 1;
-        end
-        
-        for j = stf(i).numOfRays:-1:1
-            for k = stf(i).numOfBixelsPerRay(j):-1:1
-                maskEnergy = stf(i).ray(j).energy(k) == availableEnergies;
-                if ~useEnergyBool(maskEnergy)
-                    stf(i).ray(j).energy(k)     = [];
-                    stf(i).ray(j).focusIx(k)    = [];
-                    stf(i).numOfBixelsPerRay(j) = stf(i).numOfBixelsPerRay(j) - 1;
-                end
-            end
-            if isempty(stf(i).ray(j).energy)
-                stf(i).ray(j) = [];
-                stf(i).numOfBixelsPerRay(j) = [];
-                stf(i).numOfRays = stf(i).numOfRays - 1;
-            end
-        end
-        
-    end
-    
-    % save total number of bixels
-    stf(i).totalNumOfBixels = sum(stf(i).numOfBixelsPerRay);
-    
-    % Show progress
-    matRad_progress(i,length(pln.propStf.gantryAngles));
-
-    %% visualization
-    if visMode > 0
-        
-        clf;
-        % first subplot: visualization in bev
-        subplot(1,2,1)
-        hold on
-        
-        % plot rotated target coordinates
-        plot3(rot_coords(:,1),rot_coords(:,2),rot_coords(:,3),'r.')
-        
-        % surface rendering
-        if visMode == 2
-            
-            % generate a 3D rectangular grid centered at isocenter in
-            % voxel coordinates
-            [X,Y,Z] = meshgrid((1:ct.cubeDim(2))-stf(i).isoCenter(1)/ct.resolution.x, ...
-                               (1:ct.cubeDim(1))-stf(i).isoCenter(2)/ct.resolution.y, ...
-                               (1:ct.cubeDim(3))-stf(i).isoCenter(3)/ct.resolution.z);
-            
-            % computes surface
-            patSurfCube      = 0*ct.cube{1};
-            idx              = [cst{:,4}];
-            idx              = unique(vertcat(idx{:}));
-            patSurfCube(idx) = 1;
-            
-            [f,v] = isosurface(X,Y,Z,patSurfCube,.5);
-            
-            % convert isosurface from voxel to [mm]
-            v(:,1) = v(:,1)*ct.resolution.x;
-            v(:,2) = v(:,2)*ct.resolution.y;
-            v(:,3) = v(:,3)*ct.resolution.z;
-            
-            % rotate surface
-            rotated_surface = v*rotMat_system_T;
-            
-            % surface rendering
-            surface = patch('Faces',f,'Vertices',rotated_surface);
-            set(surface,'FaceColor',[0 0 1],'EdgeColor','none','FaceAlpha',.4);
-            lighting gouraud;
-        
-        end
-        
-        % plot projection matrix: coordinates at isocenter
-        plot3(rayPos(:,1),rayPos(:,2),rayPos(:,3),'k.');
-        
-        % Plot matrix border of matrix at isocenter
-        for j = 1:stf(i).numOfRays
-            
-            % Compute border for every bixels
-            targetPoint_vox_X_1 = stf(i).ray(j).targetPoint_bev(:,1) + pln.propStf.bixelWidth;
-            targetPoint_vox_Y_1 = stf(i).ray(j).targetPoint_bev(:,2);
-            targetPoint_vox_Z_1 = stf(i).ray(j).targetPoint_bev(:,3) + pln.propStf.bixelWidth;
-            
-            targetPoint_vox_X_2 = stf(i).ray(j).targetPoint_bev(:,1) + pln.propStf.bixelWidth;
-            targetPoint_vox_Y_2 = stf(i).ray(j).targetPoint_bev(:,2);
-            targetPoint_vox_Z_2 = stf(i).ray(j).targetPoint_bev(:,3) - pln.propStf.bixelWidth;
-            
-            targetPoint_vox_X_3 = stf(i).ray(j).targetPoint_bev(:,1) - pln.propStf.bixelWidth;
-            targetPoint_vox_Y_3 = stf(i).ray(j).targetPoint_bev(:,2);
-            targetPoint_vox_Z_3 = stf(i).ray(j).targetPoint_bev(:,3) - pln.propStf.bixelWidth;
-            
-            targetPoint_vox_X_4 = stf(i).ray(j).targetPoint_bev(:,1) - pln.propStf.bixelWidth;
-            targetPoint_vox_Y_4 = stf(i).ray(j).targetPoint_bev(:,2);
-            targetPoint_vox_Z_4 = stf(i).ray(j).targetPoint_bev(:,3) + pln.propStf.bixelWidth;
-            
-            % plot
-            plot3([stf(i).sourcePoint_bev(1) targetPoint_vox_X_1],[stf(i).sourcePoint_bev(2) targetPoint_vox_Y_1],[stf(i).sourcePoint_bev(3) targetPoint_vox_Z_1],'g')
-            plot3([stf(i).sourcePoint_bev(1) targetPoint_vox_X_2],[stf(i).sourcePoint_bev(2) targetPoint_vox_Y_2],[stf(i).sourcePoint_bev(3) targetPoint_vox_Z_2],'g')
-            plot3([stf(i).sourcePoint_bev(1) targetPoint_vox_X_3],[stf(i).sourcePoint_bev(2) targetPoint_vox_Y_3],[stf(i).sourcePoint_bev(3) targetPoint_vox_Z_3],'g')
-            plot3([stf(i).sourcePoint_bev(1) targetPoint_vox_X_4],[stf(i).sourcePoint_bev(2) targetPoint_vox_Y_4],[stf(i).sourcePoint_bev(3) targetPoint_vox_Z_4],'g')
-            
-        end
-        
-        % Plot properties
-        daspect([1 1 1]);
-        view(0,-90);
-        xlabel 'X [mm]'
-        ylabel 'Y [mm]'
-        zlabel 'Z [mm]'
-        title ('Beam''s eye view')
-        axis([-300 300 -300 300 -300 300]);
-        
-        % second subplot: visualization in lps coordinate system
-        subplot(1,2,2)
-        
-        % Plot target coordinates whitout any rotation
-        plot3(coordsX,coordsY,coordsZ,'r.')
-        hold on;
-        
-        % Rotated projection matrix at isocenter
-        isocenter_plane_coor = rayPos*rotMat_vectors_T;
-        
-        % Plot isocenter plane
-        plot3(isocenter_plane_coor(:,1),isocenter_plane_coor(:,2),isocenter_plane_coor(:,3),'y.');
-        
-        % Plot rotated bixels border.
-        for j = 1:stf(i).numOfRays
-            % Generate rotated projection target points.
-            targetPoint_vox_1_rotated = [stf(i).ray(j).targetPoint_bev(:,1) + pln.propStf.bixelWidth,stf(i).ray(j).targetPoint_bev(:,2),stf(i).ray(j).targetPoint_bev(:,3) + pln.propStf.bixelWidth]*rotMat_vectors_T;
-            targetPoint_vox_2_rotated = [stf(i).ray(j).targetPoint_bev(:,1) + pln.propStf.bixelWidth,stf(i).ray(j).targetPoint_bev(:,2),stf(i).ray(j).targetPoint_bev(:,3) - pln.propStf.bixelWidth]*rotMat_vectors_T;
-            targetPoint_vox_3_rotated = [stf(i).ray(j).targetPoint_bev(:,1) - pln.propStf.bixelWidth,stf(i).ray(j).targetPoint_bev(:,2),stf(i).ray(j).targetPoint_bev(:,3) - pln.propStf.bixelWidth]*rotMat_vectors_T;
-            targetPoint_vox_4_rotated = [stf(i).ray(j).targetPoint_bev(:,1) - pln.propStf.bixelWidth,stf(i).ray(j).targetPoint_bev(:,2),stf(i).ray(j).targetPoint_bev(:,3) + pln.propStf.bixelWidth]*rotMat_vectors_T;
-            
-            % Plot rotated target points.
-            plot3([stf(i).sourcePoint(1) targetPoint_vox_1_rotated(:,1)],[stf(i).sourcePoint(2) targetPoint_vox_1_rotated(:,2)],[stf(i).sourcePoint(3) targetPoint_vox_1_rotated(:,3)],'g')
-            plot3([stf(i).sourcePoint(1) targetPoint_vox_2_rotated(:,1)],[stf(i).sourcePoint(2) targetPoint_vox_2_rotated(:,2)],[stf(i).sourcePoint(3) targetPoint_vox_2_rotated(:,3)],'g')
-            plot3([stf(i).sourcePoint(1) targetPoint_vox_3_rotated(:,1)],[stf(i).sourcePoint(2) targetPoint_vox_3_rotated(:,2)],[stf(i).sourcePoint(3) targetPoint_vox_3_rotated(:,3)],'g')
-            plot3([stf(i).sourcePoint(1) targetPoint_vox_4_rotated(:,1)],[stf(i).sourcePoint(2) targetPoint_vox_4_rotated(:,2)],[stf(i).sourcePoint(3) targetPoint_vox_4_rotated(:,3)],'g')
-        end
-        
-        % surface rendering
-        if visMode == 2
-            surface = patch('Faces',f,'Vertices',v);
-            set(surface,'FaceColor',[0 0 1],'EdgeColor','none','FaceAlpha',.4);
-            lighting gouraud;
-        end
-        
-        % labels etc.
-        daspect([1 1 1]);
-        view(0,-90);
-        xlabel 'X [mm]'
-        ylabel 'Y [mm]'
-        zlabel 'Z [mm]'
-        title 'lps coordinate system'
-        axis([-300 300 -300 300 -300 300]);
-        %pause(1);
-    end
-    
-    % include rangeshifter data if not yet available 
-    if strcmp(pln.radiationMode, 'protons') || strcmp(pln.radiationMode, 'carbon')
-        for j = 1:stf(i).numOfRays
-            for k = 1:numel(stf(i).ray(j).energy)
-                stf(i).ray(j).rangeShifter(k).ID = 0;
-                stf(i).ray(j).rangeShifter(k).eqThickness = 0;
-                stf(i).ray(j).rangeShifter(k).sourceRashiDistance = 0;
-            end
-        end
-    end
-        
-end    
-
-end
diff --git a/matRad_saveStructs.m b/matRad_saveStructs.m
index 1e7695b6c..a756b04ce 100644
--- a/matRad_saveStructs.m
+++ b/matRad_saveStructs.m
@@ -1,4 +1,4 @@
-function status = matRad_saveStructs(path, engine)
+function status = matRad_saveStructs(load_path, save_path, engine)
 % matRad_saveStructs Mat file transfer for python interface
 %
 % input
@@ -18,14 +18,14 @@
 %
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-load(path);
+load(load_path);
 
 if engine=='matlab'
-    save('ct.mat', 'ct');
-    save('cst.mat', 'cst');
+    save(append(save_path, 'ct.mat'), 'ct');
+    save(append(save_path, 'cst.mat'), 'cst');
 else
-    save('ct.mat', '-mat7-binary', 'ct');
-    save('cst.mat', '-mat7-binary', 'cst');
+    save(append(save_path, 'ct.mat'), '-mat7-binary', 'ct');
+    save(append(save_path, 'ct.mat'), '-mat7-binary', 'cst');
 end
 
 %Choosing the engine is necessary because Octave has trouble reading .mat files. Doesn't recognize them as binary.

From bbfd333faad7f1cfbd5e30c4c779764cae8734f5 Mon Sep 17 00:00:00 2001
From: g-stanic <goran.stanic@dkfz-heidelberg.de>
Date: Thu, 31 Aug 2023 10:16:39 +0200
Subject: [PATCH 4/4] Corrections of example scripts which caused the tests to
 error out

---
 examples/matRad_example2_photons.m   | 2 +-
 examples/matRad_example4_photonsMC.m | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/examples/matRad_example2_photons.m b/examples/matRad_example2_photons.m
index 5c637f1a1..4bbc84477 100644
--- a/examples/matRad_example2_photons.m
+++ b/examples/matRad_example2_photons.m
@@ -106,7 +106,7 @@
 % to 30. Internally, matRad considers the fraction dose for optimization, 
 % however, objetives and constraints are defined for the entire treatment.
 pln.numOfFractions         = 30;
-pln.propStf.gantryAngles   = [0, 90, 180, 270];
+pln.propStf.gantryAngles   = [0:40:359];
 pln.propStf.couchAngles    = zeros(1,numel(pln.propStf.gantryAngles));
 pln.propStf.bixelWidth     = 5;
 
diff --git a/examples/matRad_example4_photonsMC.m b/examples/matRad_example4_photonsMC.m
index aa5afa7ab..e6e34fbe3 100644
--- a/examples/matRad_example4_photonsMC.m
+++ b/examples/matRad_example4_photonsMC.m
@@ -56,7 +56,7 @@
 %% Dose Calculation
 % Calculate dose influence matrix for unit pencil beam intensities using 
 % a Monte Carlo algorithm
-dij = matRad_calcPhotonDoseMC(ct,stf,pln,cst,1000);
+dij = matRad_calcPhotonDoseMC(ct,stf,pln,cst);
 
 %% Inverse Optimization for IMRT
 resultGUI = matRad_fluenceOptimization(dij,cst,pln);