diff --git a/matRad/basedata/matRad_MCemittanceBaseData.m b/matRad/basedata/matRad_MCemittanceBaseData.m index e081dc28c..0f20655f3 100644 --- a/matRad/basedata/matRad_MCemittanceBaseData.m +++ b/matRad/basedata/matRad_MCemittanceBaseData.m @@ -613,7 +613,7 @@ end end - methods (Access = protected) + methods %(Access = protected) function obj = getRangeShiftersFromStf(obj,stf) allRays = [stf.ray]; raShis = [allRays.rangeShifter]; diff --git a/matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m b/matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m index 5d653adcc..8b117dde6 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m +++ b/matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m @@ -31,6 +31,8 @@ %Other Dose Calculation Properties calcLET = true; + + externalCalculation = 'off'; end properties (SetAccess = protected, GetAccess = public) @@ -140,6 +142,13 @@ function setDefaults(this) matRad_cfg = MatRad_Config.instance(); + if isfolder(this.externalCalculation) + %dij = this.readFiles(this.externalCalculation); + matRad_cfg.dispError('MCsquare not yet configured for loading external simulaiton results.'); + return; + else + end + %Now we can run initDoseCalc as usual dij = this.initDoseCalc(ct,cst,stf); @@ -162,6 +171,7 @@ function setDefaults(this) % Calculate MCsquare base data % Argument stf is optional, if given, calculation only for energies given in stf MCsquareBDL = matRad_MCsquareBaseData(this.machine); + MCsquareBDL = MCsquareBDL.getRangeShiftersFromStf(stf); %matRad_createMCsquareBaseDataFile(bdFile,machine,1); bdlFolder = fullfile(this.workingDir,'BDL'); @@ -275,15 +285,22 @@ function setDefaults(this) counter = 0; + stfMCsquare = []; + for i = 1:length(stf) - %Create new stf for MCsquare with energy layer ordering and - %shifted scenario isocenter - stfMCsquare(i).isoCenter = matRad_world2cubeCoords(stf(i).isoCenter, ct) + isoCenterShift; %MCsquare uses the isoCenter in cubeCoords - stfMCsquare(i).gantryAngle = mod(180-stf(i).gantryAngle,360); %Different MCsquare geometry - stfMCsquare(i).couchAngle = stf(i).couchAngle; - stfMCsquare(i).energies = unique([stf(i).ray.energy]); - stfMCsquare(i).SAD = stf(i).SAD; - + + stfFieldMCsquare = []; + + stfFieldMCsquare.isoCenter = matRad_world2cubeCoords(stf(i).isoCenter, ct) + isoCenterShift; %MCsquare uses the isoCenter in cubeCoords + stfFieldMCsquare.gantryAngle = mod(180-stf(i).gantryAngle,360); %Different MCsquare geometry + stfFieldMCsquare.couchAngle = stf(i).couchAngle; + stfFieldMCsquare.energies = unique([stf(i).ray.energy]); + stfFieldMCsquare.SAD = stf(i).SAD; + stfFieldMCsquare.originalStfFieldIndex = i; % Required for ordering later + + stfFieldMCsquare.rangeShifterID = 0; + stfFieldMCsquare.rangeShifterType = 'binary'; + %Let's check if we have a unique or no range shifter, because MCsquare %only allows one range shifter type per field which can be IN or OUT %per spot @@ -302,23 +319,68 @@ function setDefaults(this) matRad_cfg.dispError('MCsquare does not support different range shifter IDs per field! Aborting.\n'); end + + %Create new stf for MCsquare with energy layer ordering and + %shifted scenario isocenter + % Need to split the current stf field into two separate fields for MCsquare, with and without RaSh + stfFieldMCsquareRaShi = []; if ~isempty(raShiField) - stfMCsquare(i).rangeShifterID = raShiField; - stfMCsquare(i).rangeShifterType = 'binary'; - else - stfMCsquare(i).rangeShifterID = 0; - stfMCsquare(i).rangeShifterType = 'binary'; + % Copy the field information + stfFieldMCsquareRaShi = stfFieldMCsquare; + + stfFieldMCsquareRaShi.rangeShifterID = raShiField; + stfFieldMCsquareRaShi.rangeShifterType = 'binary'; + + % Select the energies that have a RaShi for + % this stf field + raShiLayers = []; + for j = 1:stf(i).numOfRays + currentRay = stf(i).ray(j); + raShiLayers = [raShiLayers, currentRay.energy([currentRay.rangeShifter.ID] == stfFieldMCsquareRaShi.rangeShifterID)]; + end + stfFieldMCsquareRaShi.energies = unique(raShiLayers); + + % Need to delete an energy layer from non rashi + % field if the layer is only delivered with + % rashi + + % Extract all single bixel RaShiIDs + allRangeShifterIDs = arrayfun(@(rashi) rashi.ID, [stf(i).ray.rangeShifter]); + allEnergies = [stf(i).ray.energy]; + for layerEnergy = stfFieldMCsquareRaShi.energies + raShiIDSinLayer = unique(allRangeShifterIDs(allEnergies == layerEnergy)); + + % If there is only one raShiId for this + % layer, and its the current range shifter, + % we can eliminate the layer from the + % non-RaShi field + if isscalar(raShiIDSinLayer) && raShiIDSinLayer == stfFieldMCsquareRaShi.rangeShifterID + stfFieldMCsquare(stfFieldMCsquare.energies == layerEnergy) = []; + end + end + + % allocate empty target point container + for j = 1:numel(stfFieldMCsquareRaShi.energies) + stfFieldMCsquareRaShi.energyLayer(j).targetPoints = []; + stfFieldMCsquareRaShi.energyLayer(j).numOfPrimaries = []; + stfFieldMCsquareRaShi.energyLayer(j).MU = []; + stfFieldMCsquareRaShi.energyLayer(j).rayNum = []; + stfFieldMCsquareRaShi.energyLayer(j).bixelNum = []; + end + + end % allocate empty target point container - for j = 1:numel(stfMCsquare(i).energies) - stfMCsquare(i).energyLayer(j).targetPoints = []; - stfMCsquare(i).energyLayer(j).numOfPrimaries = []; - stfMCsquare(i).energyLayer(j).MU = []; - stfMCsquare(i).energyLayer(j).rayNum = []; - stfMCsquare(i).energyLayer(j).bixelNum = []; + for j = 1:numel(stfFieldMCsquare.energies) + stfFieldMCsquare.energyLayer(j).targetPoints = []; + stfFieldMCsquare.energyLayer(j).numOfPrimaries = []; + stfFieldMCsquare.energyLayer(j).MU = []; + stfFieldMCsquare.energyLayer(j).rayNum = []; + stfFieldMCsquare.energyLayer(j).bixelNum = []; end + for j = 1:stf(i).numOfRays for k = 1:stf(i).numOfBixelsPerRay(j) counter = counter + 1; @@ -327,59 +389,111 @@ function setDefaults(this) dij.bixelNum(counter) = k; end - for k = 1:numel(stfMCsquare(i).energies) + for k = 1:numel(stfFieldMCsquare.energies) %Check if ray has a spot in the current energy layer - if any(stf(i).ray(j).energy == stfMCsquare(i).energies(k)) - energyIx = find(stf(i).ray(j).energy == stfMCsquare(i).energies(k)); - stfMCsquare(i).energyLayer(k).rayNum = [stfMCsquare(i).energyLayer(k).rayNum j]; - stfMCsquare(i).energyLayer(k).bixelNum = [stfMCsquare(i).energyLayer(k).bixelNum energyIx]; - stfMCsquare(i).energyLayer(k).targetPoints = [stfMCsquare(i).energyLayer(k).targetPoints; ... + if any(stf(i).ray(j).energy == stfFieldMCsquare.energies(k)) + energyIx = find(stf(i).ray(j).energy == stfFieldMCsquare.energies(k)); + + % If more than one energy layer is + % found, one of them is for the + % RaShiField + energyIx = energyIx([stf(i).ray(j).rangeShifter(energyIx).ID] == 0); % Select the one with no RaShi; + + if isempty(energyIx) + continue; + end + + stfFieldMCsquare.energyLayer(k).rayNum = [stfFieldMCsquare.energyLayer(k).rayNum j]; + stfFieldMCsquare.energyLayer(k).bixelNum = [stfFieldMCsquare.energyLayer(k).bixelNum energyIx]; + stfFieldMCsquare.energyLayer(k).targetPoints = [stfFieldMCsquare.energyLayer(k).targetPoints; ... -stf(i).ray(j).rayPos_bev(1) stf(i).ray(j).rayPos_bev(3)]; %Number of primaries depending on beamlet-wise or field-based compuation (direct dose calculation) if this.calcDoseDirect - stfMCsquare(i).energyLayer(k).numOfPrimaries = [stfMCsquare(i).energyLayer(k).numOfPrimaries ... - round(stf(i).ray(j).weight(stf(i).ray(j).energy == stfMCsquare(i).energies(k))*this.numHistoriesDirect)]; + stfFieldMCsquare.energyLayer(k).numOfPrimaries = [stfFieldMCsquare.energyLayer(k).numOfPrimaries ... + round(stf(i).ray(j).weight(energyIx)*this.numHistoriesDirect)]; - stfMCsquare(i).energyLayer(k).MU = [stfMCsquare(i).energyLayer(k).MU ... - round(stf(i).ray(j).weight(stf(i).ray(j).energy == stfMCsquare(i).energies(k))*this.numHistoriesDirect)]; + stfFieldMCsquare.energyLayer(k).MU = [stfFieldMCsquare.energyLayer(k).MU ... + round(stf(i).ray(j).weight(energyIx)*this.numHistoriesDirect)]; - totalWeights = totalWeights + stf(i).ray(j).weight(stf(i).ray(j).energy == stfMCsquare(i).energies(k)); + totalWeights = totalWeights + stf(i).ray(j).weight(energyIx); else - stfMCsquare(i).energyLayer(k).numOfPrimaries = [stfMCsquare(i).energyLayer(k).numOfPrimaries ... + stfFieldMCsquare.energyLayer(k).numOfPrimaries = [stfFieldMCsquare.energyLayer(k).numOfPrimaries ... this.numHistoriesPerBeamlet]; - stfMCsquare(i).energyLayer(k).MU = [stfMCsquare(i).energyLayer(k).MU ... + stfFieldMCsquare.energyLayer(k).MU = [stfFieldMCsquare.energyLayer(k).MU ... this.numHistoriesPerBeamlet]; end - %Now add the range shifter - raShis = stf(i).ray(j).rangeShifter(energyIx); - %sanity check range shifters - raShiIDs = unique([raShis.ID]); - %raShiIDs = raShiIDs(raShiIDs ~= 0); + end + end - if ~isscalar(raShiIDs) - matRad_cfg.dispError('MCsquare only supports one range shifter setting (on or off) per energy! Aborting.\n'); + % Add the bixels to the RaShi field if any + if ~isempty(raShiField) + for k = 1:numel(stfFieldMCsquareRaShi.energies) + if any(stf(i).ray(j).energy == stfFieldMCsquareRaShi.energies(k)) + + energyIx = find(stf(i).ray(j).energy == stfFieldMCsquareRaShi.energies(k)); + + % If more than one energy layer is + % found, one of them is for the + % RaShiField + energyIx = energyIx([stf(i).ray(j).rangeShifter(energyIx).ID] == stfFieldMCsquareRaShi.rangeShifterID); % Select the one with no RaShi; + + stfFieldMCsquareRaShi.energyLayer(k).rayNum = [stfFieldMCsquareRaShi.energyLayer(k).rayNum j]; + stfFieldMCsquareRaShi.energyLayer(k).bixelNum = [stfFieldMCsquareRaShi.energyLayer(k).bixelNum energyIx]; + stfFieldMCsquareRaShi.energyLayer(k).targetPoints = [stfFieldMCsquareRaShi.energyLayer(k).targetPoints; ... + -stf(i).ray(j).rayPos_bev(1) stf(i).ray(j).rayPos_bev(3)]; + + %Number of primaries depending on beamlet-wise or field-based computation (direct dose calculation) + if this.calcDoseDirect + stfFieldMCsquareRaShi.energyLayer(k).numOfPrimaries = [stfFieldMCsquareRaShi.energyLayer(k).numOfPrimaries ... + round(stf(i).ray(j).weight(energyIx)*this.numHistoriesDirect)]; + + stfFieldMCsquareRaShi.energyLayer(k).MU = [stfFieldMCsquareRaShi.energyLayer(k).MU ... + round(stf(i).ray(j).weight(energyIx)*this.numHistoriesDirect)]; + + totalWeights = totalWeights + stf(i).ray(j).weight(energyIx); + else + stfFieldMCsquareRaShi.energyLayer(k).numOfPrimaries = [stfFieldMCsquareRaShi.energyLayer(k).numOfPrimaries ... + this.numHistoriesPerBeamlet]; + + stfFieldMCsquareRaShi.energyLayer(k).MU = [stfFieldMCsquareRaShi.energyLayer(k).MU ... + this.numHistoriesPerBeamlet]; + end + + %Now add the range shifter + raShis = stf(i).ray(j).rangeShifter(energyIx); + + %sanity check range shifters + raShiIDs = unique([raShis.ID]); + %raShiIDs = raShiIDs(raShiIDs ~= 0); + + if ~isscalar(raShiIDs) + matRad_cfg.dispError('MCsquare only supports one range shifter setting (on or off) per energy! Aborting.\n'); + end + + stfFieldMCsquareRaShi.energyLayer(k).rangeShifter = raShis(1); end - - stfMCsquare(i).energyLayer(k).rangeShifter = raShis(1); end end end + + % Fill the stfMCsquare + stfMCsquare = [stfMCsquare,stfFieldMCsquare, stfFieldMCsquareRaShi]; end % remember order counterMCsquare = 0; MCsquareOrder = NaN * ones(dij.totalNumOfBixels,1); - for i = 1:length(stf) + for i = 1:length(stfMCsquare) for j = 1:numel(stfMCsquare(i).energies) for k = 1:numel(stfMCsquare(i).energyLayer(j).numOfPrimaries) counterMCsquare = counterMCsquare + 1; - ix = find(i == dij.beamNum & ... - stfMCsquare(i).energyLayer(j).rayNum(k) == dij.rayNum & ... - stfMCsquare(i).energyLayer(j).bixelNum(k) == dij.bixelNum); + ix = find(stfMCsquare(i).originalStfFieldIndex == dij.beamNum & ... + stfMCsquare(i).energyLayer(j).rayNum(k) == dij.rayNum & ... + stfMCsquare(i).energyLayer(j).bixelNum(k) == dij.bixelNum); MCsquareOrder(ix) = counterMCsquare; end @@ -403,104 +517,108 @@ function setDefaults(this) %% MC computation and dij filling - % run MCsquare - mcSquareCall = [this.mcSquareBinary ' ' MCsquareConfigFile]; - matRad_cfg.dispInfo(['Calling Monte Carlo Engine: ' mcSquareCall]); - if matRad_cfg.logLevel >= 3 - [status,cmdout] = system(mcSquareCall,'-echo'); - else - [status,cmdout] = system(mcSquareCall); - matRad_cfg.dispInfo(cmdout); - end - if status == 0 - matRad_cfg.dispInfo('MCsquare exited successfully with status %d!',status); + if strcmp(this.externalCalculation,'write') + matRad_cfg.dispInfo(['MCsquare simulation skipped for external calculation\nFiles have been written to: "',strrep(this.workingDir,'\','\\'),'"']); else - matRad_cfg.dispInfo('MCsquare did not exit successfully with status %d! Results might be compromised!',status); - end - - mask = false(dij.doseGrid.numOfVoxels,1); - mask(this.VdoseGrid) = true; - - if this.calcDoseDirect - if abs(totalWeights-sum(this.directWeights)) > 1e-2 - matRad_cfg.dispWarning('Sum of provided weights and weights used in MCsquare inconsistent!'); + % run MCsquare + mcSquareCall = [this.mcSquareBinary ' ' sprintf('"%s"', MCsquareConfigFile)]; + matRad_cfg.dispInfo(['Calling Monte Carlo Engine: ' mcSquareCall]); + if matRad_cfg.logLevel >= 3 + [status,cmdout] = system(mcSquareCall,'-echo'); + else + [status,cmdout] = system(mcSquareCall); + matRad_cfg.dispInfo(cmdout); end - finalResultWeight = absCalibrationFactorMC2 * totalWeights; - else - finalResultWeight = absCalibrationFactorMC2; - end - - % read sparse matrix - if ~this.calcDoseDirect - dij.physicalDose{ctScen,shiftScen,rangeShiftScen} = finalResultWeight * matRad_sparseBeamletsReaderMCsquare ( ... - [this.config.Output_Directory filesep 'Sparse_Dose.bin'], ... - dij.doseGrid.dimensions, ... - dij.totalNumOfBixels, ... - mask); - - %Read sparse LET - if this.calcLET - dij.mLETDose{ctScen,shiftScen,rangeShiftScen} = dij.physicalDose{ctScen,shiftScen,rangeShiftScen} .* matRad_sparseBeamletsReaderMCsquare ( ... - [this.config.Output_Directory filesep 'Sparse_LET.bin'], ... + if status == 0 + matRad_cfg.dispInfo('MCsquare exited successfully with status %d!',status); + else + matRad_cfg.dispInfo('MCsquare did not exit successfully with status %d! Results might be compromised!',status); + end + + mask = false(dij.doseGrid.numOfVoxels,1); + mask(this.VdoseGrid) = true; + + if this.calcDoseDirect + if abs(totalWeights-sum(this.directWeights)) > 1e-2 + matRad_cfg.dispWarning('Sum of provided weights and weights used in MCsquare inconsistent!'); + end + finalResultWeight = absCalibrationFactorMC2 * totalWeights; + else + finalResultWeight = absCalibrationFactorMC2; + end + + % read sparse matrix + if ~this.calcDoseDirect + dij.physicalDose{ctScen,shiftScen,rangeShiftScen} = finalResultWeight * matRad_sparseBeamletsReaderMCsquare ( ... + [this.config.Output_Directory filesep 'Sparse_Dose.bin'], ... dij.doseGrid.dimensions, ... dij.totalNumOfBixels, ... mask); + + %Read sparse LET + if this.calcLET + dij.mLETDose{ctScen,shiftScen,rangeShiftScen} = dij.physicalDose{ctScen,shiftScen,rangeShiftScen} .* matRad_sparseBeamletsReaderMCsquare ( ... + [this.config.Output_Directory filesep 'Sparse_LET.bin'], ... + dij.doseGrid.dimensions, ... + dij.totalNumOfBixels, ... + mask); + end + + % reorder influence matrix to comply with matRad default ordering + dij.physicalDose = cellfun(@(mx) mx(:,MCsquareOrder),dij.physicalDose,'UniformOutput',false); + if this.calcLET + dij.mLETDose = cellfun(@(mx) mx(:,MCsquareOrder),dij.mLETDose,'UniformOutput',false); + end + else + cube = this.readMhd('Dose.mhd'); + dij.physicalDose{ctScen,shiftScen,rangeShiftScen} = sparse(this.VdoseGrid,ones(numel(this.VdoseGrid),1), ... + finalResultWeight * cube(this.VdoseGrid), ... + dij.doseGrid.numOfVoxels,1); + + %Read LET cube + if this.calcLET + cube = this.readMhd('LET.mhd'); + dij.mLETDose{ctScen,shiftScen,rangeShiftScen} = dij.physicalDose{ctScen,shiftScen,rangeShiftScen} .* sparse(this.VdoseGrid,ones(numel(this.VdoseGrid),1), ... + cube(this.VdoseGrid), ... + dij.doseGrid.numOfVoxels,1); + end + + % Postprocessing for dij: + % This is already the combined dose over all bixels, so all parameters are 1 in this case + dij = rmfield(dij,'MCsquareCalcOrder'); + + dij.numOfBeams = 1; + dij.beamNum = 1; + dij.bixelNum = 1; + dij.rayNum = 1; + dij.totalNumOfBixels = 1; + dij.totalNumOfRays = 1; + dij.numOfRaysPerBeam = 1; end - - % reorder influence matrix to comply with matRad default ordering - dij.physicalDose = cellfun(@(mx) mx(:,MCsquareOrder),dij.physicalDose,'UniformOutput',false); - if this.calcLET - dij.mLETDose = cellfun(@(mx) mx(:,MCsquareOrder),dij.mLETDose,'UniformOutput',false); + + + if this.config.Beamlet_Mode + end - else - cube = this.readMhd('Dose.mhd'); - dij.physicalDose{ctScen,shiftScen,rangeShiftScen} = sparse(this.VdoseGrid,ones(numel(this.VdoseGrid),1), ... - finalResultWeight * cube(this.VdoseGrid), ... - dij.doseGrid.numOfVoxels,1); - - %Read LET cube - if this.calcLET - cube = this.readMhd('LET.mhd'); - dij.mLETDose{ctScen,shiftScen,rangeShiftScen} = dij.physicalDose{ctScen,shiftScen,rangeShiftScen} .* sparse(this.VdoseGrid,ones(numel(this.VdoseGrid),1), ... - cube(this.VdoseGrid), ... - dij.doseGrid.numOfVoxels,1); + + matRad_cfg.dispInfo('Scenario %d of %d finished!\n',scenarioIx,this.multScen.totNumScen); + + %% clear all data + %could also be moved to the "finalize" function + delete([this.config.CT_File(1:end-4) '.*']); + fullfile(this.workingDir,'currBixels.txt'); + fullfile(this.workingDir,'MCsquareConfig.txt'); + + %For Octave temporarily disable confirmation for recursive rmdir + if strcmp(matRad_cfg.env,'OCTAVE') + rmdirConfirmState = confirm_recursive_rmdir(0); + end + rmdir(this.config.Output_Directory,'s'); + + %Reset to old confirmatoin state + if strcmp(matRad_cfg.env,'OCTAVE') + confirm_recursive_rmdir(rmdirConfirmState); end - - % Postprocessing for dij: - % This is already the combined dose over all bixels, so all parameters are 1 in this case - dij = rmfield(dij,'MCsquareCalcOrder'); - - dij.numOfBeams = 1; - dij.beamNum = 1; - dij.bixelNum = 1; - dij.rayNum = 1; - dij.totalNumOfBixels = 1; - dij.totalNumOfRays = 1; - dij.numOfRaysPerBeam = 1; - end - - - if this.config.Beamlet_Mode - - end - - matRad_cfg.dispInfo('Scenario %d of %d finished!\n',scenarioIx,this.multScen.totNumScen); - - %% clear all data - %could also be moved to the "finalize" function - delete([this.config.CT_File(1:end-4) '.*']); - fullfile(this.workingDir,'currBixels.txt'); - fullfile(this.workingDir,'MCsquareConfig.txt'); - - %For Octave temporarily disable confirmation for recursive rmdir - if strcmp(matRad_cfg.env,'OCTAVE') - rmdirConfirmState = confirm_recursive_rmdir(0); - end - rmdir(this.config.Output_Directory,'s'); - - %Reset to old confirmatoin state - if strcmp(matRad_cfg.env,'OCTAVE') - confirm_recursive_rmdir(rmdirConfirmState); end % cd back @@ -509,6 +627,30 @@ function setDefaults(this) end matRad_cfg.dispInfo('matRad: Simulation finished!\n'); + + if strcmp(this.externalCalculation, 'write') + dij.beamNum = 1; + dij.bixelNum = 1; + dij.doseGrid = this.doseGrid; + dij.numOfBeams = 1; + dij.numOfRaysPerBeam = 1; + dij.numOfScenarios = this.multScen.totNumScen; + for i = 1:this.multScen.numOfCtScen + for j = 1:this.multScen.totNumShiftScen + for k = 1:this.multScen.totNumRangeScen + if this.multScen.scenMask(i,j,k) + %TODO: loop over all expected output quantities + dij.physicalDose{i,j,k} = zeros(dij.ctGrid.numOfVoxels,1); + dij.physicalDose_std{i,j,k} = zeros(dij.ctGrid.numOfVoxels,1); + end + + end + end + end + dij.rayNum = 1; + dij.totalNumOfBixels = 1; + dij.totalNumOfRays = 1; + end %Finalize dose calculation dij = this.finalizeDose(dij); @@ -710,8 +852,8 @@ function writeInputFiles(obj,filename,stf) fprintf(fileHandle,'####RangeShifterSetting\n%s\n','IN'); pmma_rsp = 1.165; %TODO: hardcoded for now rsWidth = rangeShifter.eqThickness / pmma_rsp; - isoToRaShi = stf(i).SAD - rangeShifter.sourceRashiDistance + rsWidth; - fprintf(fileHandle,'####IsocenterToRangeShifterDistance\n%f\n',-isoToRaShi/10); %in cm + isoToRaShi = stf(i).SAD - rangeShifter.sourceRashiDistance - rsWidth; + fprintf(fileHandle,'####IsocenterToRangeShifterDistance\n%f\n',isoToRaShi/10); %in cm fprintf(fileHandle,'####RangeShifterWaterEquivalentThickness\n%f\n',rangeShifter.eqThickness); else fprintf(fileHandle,'####RangeShifterSetting\n%s\n','OUT'); diff --git a/matRad/doseCalc/+DoseEngines/matRad_ParticlePencilBeamEngineAbstract.m b/matRad/doseCalc/+DoseEngines/matRad_ParticlePencilBeamEngineAbstract.m index 64749b5f2..ae35ad73e 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_ParticlePencilBeamEngineAbstract.m +++ b/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; diff --git a/matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m b/matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m index 58ab20995..618b08c30 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m +++ b/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 diff --git a/matRad/geometry/matRad_getIsoCenter.m b/matRad/geometry/matRad_getIsoCenter.m index 16a0670b7..7c66048b3 100644 --- a/matRad/geometry/matRad_getIsoCenter.m +++ b/matRad/geometry/matRad_getIsoCenter.m @@ -69,7 +69,7 @@ % Calculated isocenter. -isoCenter = mean(coord); +isoCenter = mean(coord, 1); % Visualization diff --git a/matRad/steering/matRad_StfGeneratorParticleIMPT.m b/matRad/steering/matRad_StfGeneratorParticleIMPT.m index 7e8e2845b..b8c8b3e35 100644 --- a/matRad/steering/matRad_StfGeneratorParticleIMPT.m +++ b/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 + diff_voi = [diff([rho{shiftScen}{end}])]; + entryIx = find(diff_voi == 1); + end + + if rho{shiftScen}{end}(end)~=0 + matRad_cfg.dispWarning('Target exit on the last voxel'); + exitIx = numel(rho{shiftScen}{1}); + else + diff_voi = [diff([rho{shiftScen}{end}])]; + exitIx = find(diff_voi == -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))]; diff --git a/matRad/steering/matRad_StfGeneratorParticleRayBixelAbstract.m b/matRad/steering/matRad_StfGeneratorParticleRayBixelAbstract.m index 7862cfecb..0cb196841 100644 --- a/matRad/steering/matRad_StfGeneratorParticleRayBixelAbstract.m +++ b/matRad/steering/matRad_StfGeneratorParticleRayBixelAbstract.m @@ -18,6 +18,7 @@ properties useRangeShifter = false; + rangeShifterEqD end properties (Access = protected) @@ -66,10 +67,18 @@ 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); + + if isempty(this.rangeShifterEqD) + this.rangeShifterEqD = round(min(this.availablePeakPos)* 1.25); + end + + this.availablePeakPosRaShi = this.availablePeakPos - this.rangeShifterEqD; + + % Available PeakPositionRaShi has to have same size() as + % availablePeaPos for indexing + this.availablePeakPosRaShi(this.availablePeakPosRaShi<0) = 0; + + 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.',this.rangeShifterEqD); end if sum(this.availablePeakPos<0)>0 diff --git a/matRad/steering/matRad_StfGeneratorParticleSingleBeamlet.m b/matRad/steering/matRad_StfGeneratorParticleSingleBeamlet.m index 775617a6f..3c99d8560 100644 --- a/matRad/steering/matRad_StfGeneratorParticleSingleBeamlet.m +++ b/matRad/steering/matRad_StfGeneratorParticleSingleBeamlet.m @@ -18,6 +18,7 @@ properties energy; raShiThickness = 50; %Range shifter to be used if useRangeShifter = true; + visBool = 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.visBool); 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.visBool); end if size(this.isoCenter,1) == 1 diff --git a/test/doseCalc/test_HongPB.m b/test/doseCalc/test_HongPB.m index 40f5300b4..4688ccd12 100644 --- a/test/doseCalc/test_HongPB.m +++ b/test/doseCalc/test_HongPB.m @@ -97,6 +97,23 @@ assertExceptionThrown(@() DoseEngines.matRad_ParticleHongPencilBeamEngine.isAvailable(testData.pln)); assertFalse(DoseEngines.matRad_ParticleHongPencilBeamEngine.isAvailable(testData.pln,[])); - - - \ No newline at end of file + 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))); + + \ No newline at end of file diff --git a/test/doseCalc/test_MCsquareEngine.m b/test/doseCalc/test_MCsquareEngine.m new file mode 100644 index 000000000..37c565471 --- /dev/null +++ b/test/doseCalc/test_MCsquareEngine.m @@ -0,0 +1,56 @@ +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 + % linesConfigFile = readlines(fullfile(matRad_cfg.primaryUserFolder, 'MCsquare', 'MCsquareConfig.txt')); + 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 + % linesBixelFile = readlines(fullfile(matRad_cfg.primaryUserFolder, 'MCsquare', 'currBixels.txt')); + 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 + + diff --git a/test/doseCalc/test_sigmaRashi.m b/test/doseCalc/test_sigmaRashi.m new file mode 100644 index 000000000..4b00db312 --- /dev/null +++ b/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); diff --git a/test/steering/test_stfGEneratorParticleIMPT.m b/test/steering/test_stfGEneratorParticleIMPT.m index 5ad451e76..440764a55 100644 --- a/test/steering/test_stfGEneratorParticleIMPT.m +++ b/test/steering/test_stfGEneratorParticleIMPT.m @@ -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); \ No newline at end of file