Range shifter positioning for MC and minor bug fixing

matRad/basedata/matRad_MCemittanceBaseData.m

@@ -613,7 +613,7 @@
         end
     end
 
-    methods (Access = protected)
+    methods
         function obj = getRangeShiftersFromStf(obj,stf)
             allRays = [stf.ray];
             raShis = [allRays.rangeShifter];

matRad/doseCalc/+DoseEngines/matRad_ParticlePencilBeamEngineAbstract.m

@@ -374,7 +374,7 @@ function chooseLateralModel(this)
                 end
 
                 % compute!
-                sigmaRashi = matRad_calcSigmaRashi(currBixel.baseData.energy, ...
+                sigmaRashi = matRad_calcSigmaRashi(currBixel.baseData, ...
                     currBixel.rangeShifter, ...
                     currBixel.SSD);
 
@@ -629,7 +629,7 @@ function calcLateralParticleCutOff(this,cutOffLevel,stfElement)
                 if  strcmp(this.machine.meta.radiationMode,'protons') && rangeShifterLUT(i).eqThickness > 0 
 
                     %get max range shift
-                    sigmaRashi = matRad_calcSigmaRashi(this.machine.data(energyIx).energy, ...
+                    sigmaRashi = matRad_calcSigmaRashi(this.machine.data(energyIx), ...
                         rangeShifterLUT(i), ...
                         energySigmaLUT(i,3));
 
@@ -818,7 +818,7 @@ function calcLateralParticleCutOff(this,cutOffLevel,stfElement)
                     if rangeShifter.eqThickness > 0 && strcmp(pln.radiationMode,'protons')
 
                         % compute!
-                        sigmaRashi = matRad_calcSigmaRashi(this.machine.data(energyIx).energy,rangeShifter,maxSSD);
+                        sigmaRashi = matRad_calcSigmaRashi(this.machine.data(energyIx),rangeShifter,maxSSD);
 
                         % add to initial sigma in quadrature
                         sigmaIni_sq = sigmaIni_sq +  sigmaRashi^2;

matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m

@@ -2684,7 +2684,7 @@ function writeRangeShifter(~,fID,rangeShifter,sourceToNozzleDistance)
             fprintf(fID,'d:Ge/%s/HLZ      = %f  mm\n',rangeShifter.topasID,rsWidth/2);
             fprintf(fID,'d:Ge/%s/TransX   = 500 mm * Tf/Beam/%sOut/Value\n',rangeShifter.topasID,rangeShifter.topasID);
             fprintf(fID,'d:Ge/%s/TransY   = 0   mm\n',rangeShifter.topasID);
-            fprintf(fID,'d:Ge/%s/TransZ   = %f mm\n',rangeShifter.topasID,rangeShifter.sourceRashiDistance - sourceToNozzleDistance);
+            fprintf(fID,'d:Ge/%s/TransZ   = %f mm\n',rangeShifter.topasID,rangeShifter.sourceRashiDistance - sourceToNozzleDistance + rsWidth/2);
 
         end
 

matRad/geometry/matRad_getIsoCenter.m

@@ -69,7 +69,7 @@
 
 
 % Calculated isocenter.
-isoCenter = mean(coord);
+isoCenter = mean(coord, 1);
 
 
 % Visualization

matRad/steering/matRad_StfGeneratorParticleIMPT.m

@@ -125,9 +125,21 @@
                                     end
 
                                     % find target entry & exit
-                                    diff_voi    = [diff([rho{shiftScen}{end}])];
-                                    entryIx = find(diff_voi == 1);
-                                    exitIx = find(diff_voi == -1);
+                                    if rho{shiftScen}{end}(1)~=0
+                                        matRad_cfg.dispWarning('Target entry on the first voxel');
+                                        entryIx = 1;
+                                    else
+                                        diffVoi    = [diff([rho{shiftScen}{end}])];
+                                        entryIx = find(diffVoi == 1);
+                                    end
+
+                                    if rho{shiftScen}{end}(end)~=0
+                                        matRad_cfg.dispWarning('Target exit on the last voxel');
+                                        exitIx = numel(rho{shiftScen}{end});
+                                    else
+                                        diffVoi    = [diff([rho{shiftScen}{end}])];
+                                        exitIx = find(diffVoi == -1);
+                                    end
 
                                     %We approximate the interface using the rad depth between the last voxel before and the first voxel after the interface 
                                     % This captures the case that the first relevant voxel is a target voxel
@@ -176,9 +188,13 @@
                             %non-reachable low-range spots
                             raShiEnergies = this.availableEnergies(this.availablePeakPosRaShi >= targetEntry(k) & min(this.availablePeakPos) > this.availablePeakPosRaShi);
 
+                            if isempty(raShiEnergies)
+                                matRad_cfg.dispWarning('No energies available for range shifting, please change the range shifter thickness');
+                            end
+
                             raShi.ID = 1;
-                            raShi.eqThickness = rangeShifterEqD;
-                            raShi.sourceRashiDistance = round(min(ctEntryPoint) - 2*rangeShifterEqD,-1); %place a little away from entry, round to cms to reduce number of unique settings
+                            raShi.eqThickness = this.rangeShifterEqD;
+                           raShi.sourceRashiDistance = 10 * round((min(ctEntryPoint) - 2*this.rangeShifterEqD) / 10); %place a little away from entry, round to cms to reduce number of unique settings
 
                             beam.ray(j).energy = [beam.ray(j).energy raShiEnergies];
                             beam.ray(j).rangeShifter = [beam.ray(j).rangeShifter repmat(raShi,1,length(raShiEnergies))];

matRad/steering/matRad_StfGeneratorParticleRayBixelAbstract.m

@@ -18,6 +18,7 @@
 
     properties
         useRangeShifter = false;
+        rangeShifterEqD
     end
 
     properties (Access = protected)
@@ -66,14 +67,25 @@ function initialize(this)
             if this.useRangeShifter
                 %For now only a generic range shifter is used whose thickness is
                 %determined by the minimum peak width to play with
-                rangeShifterEqD = round(min(this.availablePeakPos)* 1.25);
-                this.availablePeakPosRaShi = this.availablePeakPos - rangeShifterEqD;
-
-                matRad_cfg.dispWarning('Use of range shifter enabled. matRad will generate a generic range shifter with WEPL %f to enable ranges below the shortest base data entry.',rangeShifterEqD);
+
+                matRad_cfg.dispInfo('\nUse of range shifter active');
+
+                if ~isempty(this.rangeShifterEqD)
+                      matRad_cfg.dispInfo('\nUsing provided range shifter thickness of %f mm\n', this.rangeShifterEqD);                  
+                else
+                    this.rangeShifterEqD = round(min(this.availablePeakPos)* 1.25);
+                    matRad_cfg.dispInfo('\nUsing generic range shifter thickness of %f mm determined to allow ranges below the shortest base data entry\n', this.rangeShifterEqD);
+                end
+
+                this.availablePeakPosRaShi = this.availablePeakPos - this.rangeShifterEqD;
+
+                % Available PeakPositionRaShi has to have same size() as
+                % availablePeakPos for indexing
+                this.availablePeakPosRaShi(this.availablePeakPosRaShi<0) = 0;
             end
 
             if sum(this.availablePeakPos<0)>0
-                matRad_cfg.dispError('at least one available peak position is negative - inconsistent machine file')
+                matRad_cfg.dispError('At least one available peak position is negative - inconsistent machine file')
             end
 
             %Create Water equivalent cube in ct

matRad/steering/matRad_StfGeneratorParticleSingleBeamlet.m

@@ -18,6 +18,7 @@
     properties 
         energy;
         raShiThickness = 50; %Range shifter to be used if useRangeShifter = true;
+        visualize = false;
     end
 
     properties (Constant)
@@ -50,12 +51,12 @@ function createPatientGeometry(this)
             matRad_cfg = MatRad_Config.instance();
 
             if isempty(this.isoCenter)
-                this.isoCenter = matRad_getIsoCenter(this.cst,this.ct,visBool);
+                this.isoCenter = matRad_getIsoCenter(this.cst,this.ct,this.visualize);
             end
 
             if ~isequal(size(this.isoCenter),[this.numOfBeams,3]) && ~size(this.isoCenter,1) ~= 1
                 matRad_cfg.dispWarning('IsoCenter invalid, creating new one automatically!');
-                this.isoCenter = matRad_getIsoCenter(this.cst,this.ct,visBool);
+                this.isoCenter = matRad_getIsoCenter(this.cst,this.ct,this.visualize);
             end
 
             if size(this.isoCenter,1) == 1          
@@ -197,7 +198,10 @@ function createPatientGeometry(this)
 
                 %Place range shifter 2 times the range away from isocenter, but
                 %at least 10 cm
-                sourceRaShi = round(ctEntryPoint - 2*this.raShiThickness,-1); %place a little away from entry, round to cms to reduce number of unique settings;
+
+                roundToDigit = @(x,n) round(x * 10^n)/10^n;
+
+                sourceRaShi = roundToDigit((ctEntryPoint - 2*this.raShiThickness),-1); %place a little away from entry, round to cms to reduce number of unique settings;
                 beam.ray.rangeShifter.sourceRashiDistance = sourceRaShi;
             else
                 beam.ray.rangeShifter.ID = 0;

test/doseCalc/test_HongPB.m

@@ -89,6 +89,23 @@
         assertExceptionThrown(@() DoseEngines.matRad_ParticleHongPencilBeamEngine.isAvailable(testData.pln));
         assertFalse(DoseEngines.matRad_ParticleHongPencilBeamEngine.isAvailable(testData.pln,[]));
 
-
-
-
+    function test_doseCalcWithRashi
+        testData = load('protons_testData.mat');
+        engine   = DoseEngines.matRad_ParticleHongPencilBeamEngine(testData.pln);
+
+        stf = testData.stf;
+
+        % Add rangeShifter
+        stf(1).ray(1).rangeShifter.ID =1;
+        stf(1).ray(1).rangeShifter.eqThickness =1;
+        stf(1).ray(1).rangeShifter.sourceRashiDistance = -(stf(1).sourcePoint(2) + 100);
+
+        % Add rangeShifter
+        stf(2).ray(2).rangeShifter.ID =1;
+        stf(2).ray(2).rangeShifter.eqThickness =1;
+        stf(2).ray(2).rangeShifter.sourceRashiDistance = -(stf(2).sourcePoint(2) + 100);
+
+        resultGUI = engine.calcDoseForward(testData.ct,testData.cst,stf,ones(sum([stf.totalNumOfBixels]),1));
+        assertTrue(isequal(fieldnames(resultGUI),fieldnames(testData.resultGUI)));
+
+

test/doseCalc/test_MCsquareEngine.m

@@ -0,0 +1,89 @@
+function test_suite = test_MCsquareEngine
+
+test_functions=localfunctions();
+
+initTestSuite;
+
+
+function test_MCsquareDoseCalcBasic
+
+    matRad_cfg = MatRad_Config.instance();
+    radModes = DoseEngines.matRad_ParticleMCsquareEngine.possibleRadiationModes;
+
+    for i = 1:numel(radModes)
+        load([radModes{i} '_testData.mat']);
+        pln.bioModel = matRad_bioModel(radModes{i},'none');
+
+        w = ones(1,sum([stf(:).totalNumOfBixels]));
+
+        pln.propDoseCalc.engine = 'MCsquare';
+        pln.propDoseCalc.externalCalculation = 'write';
+        pln.propDoseCalc.numHistoriesDirect = 42;
+        resultGUI = matRad_calcDoseForward(ct,cst,stf,pln, w);
+
+        assertTrue(exist(fullfile(matRad_cfg.primaryUserFolder, 'MCsquare'), 'dir')==7); % Check it exists and its a folder
+        assertTrue(exist(fullfile(matRad_cfg.primaryUserFolder, 'MCsquare', 'MCsquareConfig.txt'), 'file')==2);
+        assertTrue(exist(fullfile(matRad_cfg.primaryUserFolder, 'MCsquare', 'currBixels.txt'), 'file')==2);
+
+
+        % Check parameters
+        % Read config file
+        fid = fopen(fullfile(matRad_cfg.primaryUserFolder, 'MCsquare', 'MCsquareConfig.txt'),'r');
+            linesConfigFile = {};
+            while ~feof(fid)
+                linesConfigFile{end+1,1} = fgetl(fid);
+            end
+        fclose(fid);
+
+
+        assertTrue(any(strcmp(linesConfigFile, "Num_Primaries 42")));
+
+        % Read currBixel file
+        fid = fopen(fullfile(matRad_cfg.primaryUserFolder, 'MCsquare', 'currBixels.txt'),'r');
+            linesBixelFile = {};
+            while ~feof(fid)
+                linesBixelFile{end+1,1} = fgetl(fid);
+            end
+        fclose(fid);
+
+        assertTrue(any(strcmp(linesBixelFile, "##NumberOfFields")));
+        assertTrue(str2double(linesBixelFile(find(strcmp(linesBixelFile, "##NumberOfFields"))+1)) == numel(stf));
+
+    end
+
+function test_MCsquareRaShi
+    matRad_cfg = MatRad_Config.instance();
+    radModes = DoseEngines.matRad_ParticleMCsquareEngine.possibleRadiationModes;
+
+    for i = 1:numel(radModes)
+        load([radModes{i} '_testData.mat']);
+        pln.bioModel = matRad_bioModel(radModes{i},'none');
+
+        pln.propStf.gantryAngles = 0;
+        pln.propStf.couchAngles = 0;
+        stfGenerator = matRad_StfGeneratorParticleSingleBeamlet(pln);
+        stfGenerator.useRangeShifter = true;
+        stf = stfGenerator.generate(ct,cst);
+
+
+        w = 1;
+
+        pln.propDoseCalc.engine = 'MCsquare';
+        pln.propDoseCalc.externalCalculation = 'write';
+        pln.propDoseCalc.numHistoriesDirect = 42;
+        resultGUI = matRad_calcDoseForward(ct,cst,stf,pln, w);
+
+        assertTrue(exist(fullfile(matRad_cfg.primaryUserFolder, 'MCsquare'), 'dir')==7); % Check it exists and its a folder
+        assertTrue(exist(fullfile(matRad_cfg.primaryUserFolder, 'MCsquare', 'MCsquareConfig.txt'), 'file')==2);
+        assertTrue(exist(fullfile(matRad_cfg.primaryUserFolder, 'MCsquare', 'currBixels.txt'), 'file')==2);
+
+        fid = fopen(fullfile(matRad_cfg.primaryUserFolder, 'MCsquare', 'currBixels.txt'),'r');
+            linesBixelFile = {};
+            while ~feof(fid)
+                linesBixelFile{end+1,1} = fgetl(fid);
+            end
+        fclose(fid);
+
+        assertTrue(any(strcmp(linesBixelFile, "####RangeShifterWaterEquivalentThickness")));
+        assertTrue(str2double(linesBixelFile(find(strcmp(linesBixelFile, "####RangeShifterWaterEquivalentThickness"))+1)) == stf.ray.rangeShifter.eqThickness);
+    end

test/doseCalc/test_sigmaRashi.m

@@ -0,0 +1,19 @@
+function test_suite= test_sigmaRashi
+
+test_functions=localfunctions();
+
+initTestSuite;
+
+function test_calcSigmaRashi
+
+    baseDataEntry.range = 100;
+
+    rangeShifter.ID = 1;
+    rangeShifter.eqThickness = 1;
+    rangeShifter.sourceRashiDistance = 9000;
+
+    SSD = 10000;
+
+    sigma = matRad_calcSigmaRashi(baseDataEntry, rangeShifter, SSD);
+
+    assertElementsAlmostEqual(sigma,1.1,'relative',1e-2,1e-2);

test/steering/test_stfGEneratorParticleIMPT.m

@@ -1,4 +1,4 @@
-function test_suite = test_stfGeneratorPhotonIMRT
+function test_suite = test_stfGEneratorParticleIMPT
 
     test_functions=localfunctions();
 
@@ -83,3 +83,24 @@ function test_generate_multibeams()
 
             %assertTrue(isscalar(stf2(i).ray.energy));
         end
+
+function test_generateRangeShifterStf()
+        % geometry settings
+        load protons_testData.mat
+
+        % Move Target shallower so that range shifter calculation
+        ct.resolution.y=5;
+        VolHelper = false(ct.cubeDim);
+        VolHelper(2:3,5:6,5:6) = true;
+        ixTarget = find(VolHelper);
+
+        cst{1,4}{1} = ixTarget;
+
+        stfGen = matRad_StfGeneratorParticleIMPT(pln);
+        stfGen.useRangeShifter = true;
+        stfGen.rangeShifterEqD = 2;
+
+        stf = stfGen.generate(ct,cst);
+
+        assertTrue(stf(1).ray(1).rangeShifter(1).ID==1);
+        assertTrue(stf(1).ray(1).rangeShifter(1).eqThickness==2);