From 97df3fff41e7e3c41d217b1fc1884bead90c3806 Mon Sep 17 00:00:00 2001 From: s742o Date: Wed, 2 Oct 2024 17:11:21 +0200 Subject: [PATCH 1/4] Cutoff method and Analytical Engine test --- .../matRad_ParticleAnalyticalBortfeldEngine.m | 10 ++- .../matRad_ParticlePencilBeamEngineAbstract.m | 17 +++- test/doseCalc/test_Analytical.m | 80 +++++++++++++++++++ 3 files changed, 104 insertions(+), 3 deletions(-) create mode 100644 test/doseCalc/test_Analytical.m diff --git a/matRad/doseCalc/+DoseEngines/matRad_ParticleAnalyticalBortfeldEngine.m b/matRad/doseCalc/+DoseEngines/matRad_ParticleAnalyticalBortfeldEngine.m index a8f4fd46f..5fdb0be01 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_ParticleAnalyticalBortfeldEngine.m +++ b/matRad/doseCalc/+DoseEngines/matRad_ParticleAnalyticalBortfeldEngine.m @@ -74,7 +74,9 @@ methods (Access = protected) function dij = initDoseCalc(this,ct,cst,stf) - + + matRad_cfg = MatRad_Config.instance(); + if this.calcLET == true matRad_cfg.dispWarning('Engine does not support LET calculation! Disabling!'); this.calcLET = false; @@ -324,6 +326,11 @@ function calcLateralParticleCutOff(this,cutOffLevel,~) calcRange = true; end + if strcmp(this.cutOffMethod, 'relative') + warning('Relative cutoff not yet implemented for Analytical Engine. Default integral method used.') + this.cutOffMethod = 'integral'; + end + for i = 1:numel(this.machine.data) this.machine.data(i).LatCutOff.CompFac = 1/cutOffLevel; this.machine.data(i).LatCutOff.numSig = sqrt(2) * sqrt(gammaincinv(0.995,1)); %For a 2D symmetric gaussian we need the inverse of the incomplete Gamma function for defining the CutOff @@ -362,6 +369,7 @@ function calcLateralParticleCutOff(this,cutOffLevel,~) end catch msg = 'Your machine file is invalid and does not contain the basic field (meta/data/radiationMode)!'; + error(msg); return; end diff --git a/matRad/doseCalc/+DoseEngines/matRad_ParticlePencilBeamEngineAbstract.m b/matRad/doseCalc/+DoseEngines/matRad_ParticlePencilBeamEngineAbstract.m index 1e2f4b7e9..9811aeb6f 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_ParticlePencilBeamEngineAbstract.m +++ b/matRad/doseCalc/+DoseEngines/matRad_ParticlePencilBeamEngineAbstract.m @@ -26,6 +26,8 @@ airOffsetCorrection = true; % Corrects WEPL for SSD difference to kernel database lateralModel = 'auto'; % Lateral Model used. 'auto' uses the most accurate model available (i.e. multiple Gaussians). 'single','double','multi' try to force a singleGaussian or doubleGaussian model, if available + cutOffMethod = 'integral'; % or 'relative' + visBoolLateralCutOff = false; % Boolean switch for visualization during+ LeteralCutOff calculation end @@ -710,8 +712,19 @@ function calcLateralParticleCutOff(this,cutOffLevel,stfElement) matRad_cfg.dispWarning('LateralParticleCutOff: shell integration is wrong !') end - IX = find(cumArea >= idd(j) * cutOffLevel,1, 'first'); - this.machine.data(energyIx).LatCutOff.CompFac = cutOffLevel^-1; + % Find radius at which integrated dose becomes + % bigger than cutoff * IDD + if strcmp(this.cutOffMethod, 'integral') + IX = find(cumArea >= idd(j) * cutOffLevel,1, 'first'); + this.machine.data(energyIx).LatCutOff.CompFac = cutOffLevel^-1; + elseif strcmp(this.cutOffMethod, 'relative') + IX = find(dose_r <= (1-cutOffLevel) * max(dose_r), 1, 'first'); + relFac = cumArea(IX)./cumArea(end); % (or idd(j)) to find the appropriate integral of dose + this.machine.data(energyIx).LatCutOff.CompFac = relFac^-1; + else + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Invalid Cutoff Method'); + end if isempty(IX) depthDoseCutOff = Inf; diff --git a/test/doseCalc/test_Analytical.m b/test/doseCalc/test_Analytical.m new file mode 100644 index 000000000..fa86f7a51 --- /dev/null +++ b/test/doseCalc/test_Analytical.m @@ -0,0 +1,80 @@ +function test_suite = test_Analytical + +test_functions=localfunctions(); + +initTestSuite; + +function test_getAnalyticalEngineFromPln + % Single gaussian lateral model + protonDummyPln = struct('radiationMode','protons','machine','Generic'); + protonDummyPln.propDoseCalc.engine = 'AnalyticalPB'; + engine = DoseEngines.matRad_ParticleAnalyticalBortfeldEngine.getEngineFromPln(protonDummyPln); + assertTrue(isa(engine,'DoseEngines.matRad_ParticleAnalyticalBortfeldEngine')); + + % Double Gaussian lateral model + % If you don't have my clusterDose basedata you cannot try this :P + %{ + protonDummyPln = struct('radiationMode','protons','machine','Generic_clusterDose'); + protonDummyPln.propDoseCalc.engine = 'AnalyticalPB'; + engine = DoseEngines.matRad_ParticleAnalyticalBortfeldEngine.getEngineFromPln(protonDummyPln); + assertTrue(isa(engine,'DoseEngines.matRad_ParticleAnalyticalBortfeldEngine')); + %} + +function test_loadMachineForAnalytical + possibleRadModes = DoseEngines.matRad_ParticleAnalyticalBortfeldEngine.possibleRadiationModes; + for i = 1:numel(possibleRadModes) + machine = DoseEngines.matRad_ParticleAnalyticalBortfeldEngine.loadMachine(possibleRadModes{i},'Generic'); + assertTrue(isstruct(machine)); + assertTrue(isfield(machine, 'meta')); + assertTrue(isfield(machine.meta, 'radiationMode')); + assertTrue(strcmp(machine.meta.radiationMode, 'protons')); + end + +function test_calcDoseAnalytical + matRad_cfg = MatRad_Config.instance(); + + protonDummyPln = struct('radiationMode','protons','machine','Generic'); + protonDummyPln.propDoseCalc.engine = 'AnalyticalPB'; + + load([protonDummyPln.radiationMode '_' protonDummyPln.machine]); + + load BOXPHANTOM.mat + + stf = matRad_generateStf(ct, cst, protonDummyPln); + + resultGUI = matRad_calcDoseForward(ct, cst, stf, protonDummyPln, ones([1, stf(:).totalNumOfBixels])); + + assertTrue(isfield(resultGUI, 'physicalDose')); + assertTrue(isfield(resultGUI, 'w')); + assertTrue(isequal(size(ct.cube{1}), size(resultGUI.physicalDose))) + +function test_nonSupportedSettings + % Radiation mode other than protons not implemented + carbonDummyPln = struct('radiationMode','carbon','machine','Generic'); + carbonDummyPln.propDoseCalc.engine = 'AnalyticalPB'; + engine = DoseEngines.matRad_ParticleAnalyticalBortfeldEngine.getEngineFromPln(carbonDummyPln); + assertTrue(~isa(engine,'DoseEngines.matRad_ParticleAnalyticalBortfeldEngine')); + + % Biological models, LET, other lateral models not implemented + protonDummyPln = struct('radiationMode','protons','machine','Generic'); + protonDummyPln.propDoseCalc.engine = 'AnalyticalPB'; + protonDummyPln.propDoseCalc.calcLET = true; + protonDummyPln.propDoseCalc.calcBioDose = true; + protonDummyPln.propDoseCalc.lateralModel = 'double'; + engine = DoseEngines.matRad_ParticleAnalyticalBortfeldEngine.getEngineFromPln(protonDummyPln); + assertTrue(isa(engine,'DoseEngines.matRad_ParticleAnalyticalBortfeldEngine')); + load BOXPHANTOM.mat + stf = matRad_generateStf(ct, cst, protonDummyPln); + resultGUI = matRad_calcDoseForward(ct, cst, stf, protonDummyPln, ones([1, stf(:).totalNumOfBixels])); + assertTrue(~engine.calcLET) + %assertTrue(~engine.calcBioDose) % Access protected property + + % Invalid machine without radiation mode field + matRad_cfg = MatRad_Config.instance(); + protonDummyPln = struct('radiationMode','protons','machine','Empty'); + protonDummyPln.propDoseCalc.engine = 'AnalyticalPB'; + machine = []; + assertExceptionThrown(@() DoseEngines.matRad_ParticleAnalyticalBortfeldEngine.isAvailable(protonDummyPln, machine)); + + + From 1e454c6eabe9aa4800f3a5a80d490253aae78fc9 Mon Sep 17 00:00:00 2001 From: s742o Date: Fri, 4 Oct 2024 11:39:45 +0200 Subject: [PATCH 2/4] Analytical engine small modification --- .../matRad_ParticleAnalyticalBortfeldEngine.m | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/matRad/doseCalc/+DoseEngines/matRad_ParticleAnalyticalBortfeldEngine.m b/matRad/doseCalc/+DoseEngines/matRad_ParticleAnalyticalBortfeldEngine.m index 5fdb0be01..1a9c5dca9 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_ParticleAnalyticalBortfeldEngine.m +++ b/matRad/doseCalc/+DoseEngines/matRad_ParticleAnalyticalBortfeldEngine.m @@ -326,14 +326,9 @@ function calcLateralParticleCutOff(this,cutOffLevel,~) calcRange = true; end - if strcmp(this.cutOffMethod, 'relative') - warning('Relative cutoff not yet implemented for Analytical Engine. Default integral method used.') - this.cutOffMethod = 'integral'; - end - for i = 1:numel(this.machine.data) this.machine.data(i).LatCutOff.CompFac = 1/cutOffLevel; - this.machine.data(i).LatCutOff.numSig = sqrt(2) * sqrt(gammaincinv(0.995,1)); %For a 2D symmetric gaussian we need the inverse of the incomplete Gamma function for defining the CutOff + this.machine.data(i).LatCutOff.numSig = sqrt(-2*log(1-cutOffLevel)); %For a 2D symmetric gaussian we need the inverse of the incomplete Gamma function for defining the CutOff this.machine.data(i).LatCutOff.maxSigmaIni = max([this.machine.data(i).initFocus(:).SisFWHMAtIso]) ./ 2.3548; if calcRange this.machine.data(i).range = 10 * this.alpha*this.machine.data(i).energy.^this.p; @@ -369,7 +364,6 @@ function calcLateralParticleCutOff(this,cutOffLevel,~) end catch msg = 'Your machine file is invalid and does not contain the basic field (meta/data/radiationMode)!'; - error(msg); return; end From a311bba07eb36ada55c4927a9ccd81ea0742d01a Mon Sep 17 00:00:00 2001 From: Niklas Wahl Date: Fri, 4 Oct 2024 13:57:14 +0200 Subject: [PATCH 3/4] small improvement of code readability and quality --- .../matRad_ParticleAnalyticalBortfeldEngine.m | 7 +++++- .../matRad_ParticlePencilBeamEngineAbstract.m | 25 +++++++++---------- 2 files changed, 18 insertions(+), 14 deletions(-) diff --git a/matRad/doseCalc/+DoseEngines/matRad_ParticleAnalyticalBortfeldEngine.m b/matRad/doseCalc/+DoseEngines/matRad_ParticleAnalyticalBortfeldEngine.m index 1a9c5dca9..a3c53873c 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_ParticleAnalyticalBortfeldEngine.m +++ b/matRad/doseCalc/+DoseEngines/matRad_ParticleAnalyticalBortfeldEngine.m @@ -326,9 +326,14 @@ function calcLateralParticleCutOff(this,cutOffLevel,~) calcRange = true; end + if ~any(strcmp(this.cutOffMethod,{'integral','relative'})) + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('LateralParticleCutOff: Invalid Cutoff Method. Must be ''integral'' or ''relative''!'); + end + for i = 1:numel(this.machine.data) this.machine.data(i).LatCutOff.CompFac = 1/cutOffLevel; - this.machine.data(i).LatCutOff.numSig = sqrt(-2*log(1-cutOffLevel)); %For a 2D symmetric gaussian we need the inverse of the incomplete Gamma function for defining the CutOff + this.machine.data(i).LatCutOff.numSig = sqrt(-2*log(1-cutOffLevel)); this.machine.data(i).LatCutOff.maxSigmaIni = max([this.machine.data(i).initFocus(:).SisFWHMAtIso]) ./ 2.3548; if calcRange this.machine.data(i).range = 10 * this.alpha*this.machine.data(i).energy.^this.p; diff --git a/matRad/doseCalc/+DoseEngines/matRad_ParticlePencilBeamEngineAbstract.m b/matRad/doseCalc/+DoseEngines/matRad_ParticlePencilBeamEngineAbstract.m index 9811aeb6f..908316670 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_ParticlePencilBeamEngineAbstract.m +++ b/matRad/doseCalc/+DoseEngines/matRad_ParticlePencilBeamEngineAbstract.m @@ -26,7 +26,7 @@ airOffsetCorrection = true; % Corrects WEPL for SSD difference to kernel database lateralModel = 'auto'; % Lateral Model used. 'auto' uses the most accurate model available (i.e. multiple Gaussians). 'single','double','multi' try to force a singleGaussian or doubleGaussian model, if available - cutOffMethod = 'integral'; % or 'relative' + cutOffMethod = 'integral'; % or 'relative' - describes how to calculate the lateral dosimetric cutoff visBoolLateralCutOff = false; % Boolean switch for visualization during+ LeteralCutOff calculation end @@ -714,27 +714,26 @@ function calcLateralParticleCutOff(this,cutOffLevel,stfElement) % Find radius at which integrated dose becomes % bigger than cutoff * IDD - if strcmp(this.cutOffMethod, 'integral') - IX = find(cumArea >= idd(j) * cutOffLevel,1, 'first'); - this.machine.data(energyIx).LatCutOff.CompFac = cutOffLevel^-1; - elseif strcmp(this.cutOffMethod, 'relative') - IX = find(dose_r <= (1-cutOffLevel) * max(dose_r), 1, 'first'); - relFac = cumArea(IX)./cumArea(end); % (or idd(j)) to find the appropriate integral of dose - this.machine.data(energyIx).LatCutOff.CompFac = relFac^-1; - else - matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Invalid Cutoff Method'); + switch this.cutOffMethod + case 'integral' + IX = find(cumArea >= idd(j) * cutOffLevel,1, 'first'); + this.machine.data(energyIx).LatCutOff.CompFac = cutOffLevel^-1; + case 'relative' + IX = find(dose_r <= (1-cutOffLevel) * max(dose_r), 1, 'first'); + relFac = cumArea(IX)./cumArea(end); % (or idd(j)) to find the appropriate integral of dose + this.machine.data(energyIx).LatCutOff.CompFac = relFac^-1; + otherwise + matRad_cfg.dispError('LateralParticleCutOff: Invalid Cutoff Method. Must be ''integral'' or ''relative''!'); end if isempty(IX) depthDoseCutOff = Inf; - matRad_cfg.dispWarning('LateralParticleCutOff: Couldnt find lateral cut off !') + matRad_cfg.dispWarning('LateralParticleCutOff: Couldnt find lateral cut off!') elseif isnumeric(IX) depthDoseCutOff = r_mid(IX); end this.machine.data(energyIx).LatCutOff.CutOff(j) = depthDoseCutOff; - end end end From 0d7beb152d8f71fdb036ec7e13fbde9e2d45a51b Mon Sep 17 00:00:00 2001 From: Niklas Wahl Date: Fri, 4 Oct 2024 14:14:53 +0200 Subject: [PATCH 4/4] fix tests --- test/doseCalc/test_Analytical.m | 65 ++++++++++++++------------------- 1 file changed, 27 insertions(+), 38 deletions(-) diff --git a/test/doseCalc/test_Analytical.m b/test/doseCalc/test_Analytical.m index fa86f7a51..a04f53fdb 100644 --- a/test/doseCalc/test_Analytical.m +++ b/test/doseCalc/test_Analytical.m @@ -6,17 +6,17 @@ function test_getAnalyticalEngineFromPln % Single gaussian lateral model - protonDummyPln = struct('radiationMode','protons','machine','Generic'); - protonDummyPln.propDoseCalc.engine = 'AnalyticalPB'; - engine = DoseEngines.matRad_ParticleAnalyticalBortfeldEngine.getEngineFromPln(protonDummyPln); + testData.pln = struct('radiationMode','protons','machine','Generic'); + testData.pln.propDoseCalc.engine = 'AnalyticalPB'; + engine = DoseEngines.matRad_ParticleAnalyticalBortfeldEngine.getEngineFromPln(testData.pln); assertTrue(isa(engine,'DoseEngines.matRad_ParticleAnalyticalBortfeldEngine')); % Double Gaussian lateral model % If you don't have my clusterDose basedata you cannot try this :P %{ - protonDummyPln = struct('radiationMode','protons','machine','Generic_clusterDose'); - protonDummyPln.propDoseCalc.engine = 'AnalyticalPB'; - engine = DoseEngines.matRad_ParticleAnalyticalBortfeldEngine.getEngineFromPln(protonDummyPln); + testData.pln = struct('radiationMode','protons','machine','Generic_clusterDose'); + testData.pln.propDoseCalc.engine = 'AnalyticalPB'; + engine = DoseEngines.matRad_ParticleAnalyticalBortfeldEngine.getEngineFromPln(testData.pln); assertTrue(isa(engine,'DoseEngines.matRad_ParticleAnalyticalBortfeldEngine')); %} @@ -31,50 +31,39 @@ end function test_calcDoseAnalytical - matRad_cfg = MatRad_Config.instance(); - - protonDummyPln = struct('radiationMode','protons','machine','Generic'); - protonDummyPln.propDoseCalc.engine = 'AnalyticalPB'; - - load([protonDummyPln.radiationMode '_' protonDummyPln.machine]); - - load BOXPHANTOM.mat - - stf = matRad_generateStf(ct, cst, protonDummyPln); - - resultGUI = matRad_calcDoseForward(ct, cst, stf, protonDummyPln, ones([1, stf(:).totalNumOfBixels])); + testData = load('protons_testData.mat'); + assertTrue(DoseEngines.matRad_ParticleAnalyticalBortfeldEngine.isAvailable(testData.pln)) + testData.pln.propDoseCalc.engine = 'AnalyticalPB'; + resultGUI = matRad_calcDoseForward(testData.ct, testData.cst, testData.stf, testData.pln, ones(sum([testData.stf(:).totalNumOfBixels]),1)); assertTrue(isfield(resultGUI, 'physicalDose')); assertTrue(isfield(resultGUI, 'w')); - assertTrue(isequal(size(ct.cube{1}), size(resultGUI.physicalDose))) + assertTrue(isequal(testData.ct.cubeDim, size(resultGUI.physicalDose))); function test_nonSupportedSettings % Radiation mode other than protons not implemented - carbonDummyPln = struct('radiationMode','carbon','machine','Generic'); - carbonDummyPln.propDoseCalc.engine = 'AnalyticalPB'; - engine = DoseEngines.matRad_ParticleAnalyticalBortfeldEngine.getEngineFromPln(carbonDummyPln); - assertTrue(~isa(engine,'DoseEngines.matRad_ParticleAnalyticalBortfeldEngine')); - + testData = load('carbon_testData.mat'); + testData.pln.propDoseCalc.engine = 'AnalyticalPB'; + assertFalse(DoseEngines.matRad_ParticleAnalyticalBortfeldEngine.isAvailable(testData.pln)); + % Biological models, LET, other lateral models not implemented - protonDummyPln = struct('radiationMode','protons','machine','Generic'); - protonDummyPln.propDoseCalc.engine = 'AnalyticalPB'; - protonDummyPln.propDoseCalc.calcLET = true; - protonDummyPln.propDoseCalc.calcBioDose = true; - protonDummyPln.propDoseCalc.lateralModel = 'double'; - engine = DoseEngines.matRad_ParticleAnalyticalBortfeldEngine.getEngineFromPln(protonDummyPln); + testData = load('protons_testData.mat'); + testData.pln.propDoseCalc.engine = 'AnalyticalPB'; + testData.pln.propDoseCalc.calcLET = true; + testData.pln.propDoseCalc.calcBioDose = true; + testData.pln.propDoseCalc.lateralModel = 'double'; + engine = DoseEngines.matRad_ParticleAnalyticalBortfeldEngine.getEngineFromPln(testData.pln); assertTrue(isa(engine,'DoseEngines.matRad_ParticleAnalyticalBortfeldEngine')); - load BOXPHANTOM.mat - stf = matRad_generateStf(ct, cst, protonDummyPln); - resultGUI = matRad_calcDoseForward(ct, cst, stf, protonDummyPln, ones([1, stf(:).totalNumOfBixels])); + + resultGUI = matRad_calcDoseForward(testData.ct, testData.cst, testData.stf, testData.pln, ones(sum([testData.stf(:).totalNumOfBixels]),1)); assertTrue(~engine.calcLET) %assertTrue(~engine.calcBioDose) % Access protected property % Invalid machine without radiation mode field - matRad_cfg = MatRad_Config.instance(); - protonDummyPln = struct('radiationMode','protons','machine','Empty'); - protonDummyPln.propDoseCalc.engine = 'AnalyticalPB'; - machine = []; - assertExceptionThrown(@() DoseEngines.matRad_ParticleAnalyticalBortfeldEngine.isAvailable(protonDummyPln, machine)); + testData.pln.machine = 'Empty'; + testData.pln.propDoseCalc.engine = 'AnalyticalPB'; + assertExceptionThrown(@() DoseEngines.matRad_ParticleAnalyticalBortfeldEngine.isAvailable(testData.pln)); + assertFalse(DoseEngines.matRad_ParticleAnalyticalBortfeldEngine.isAvailable(testData.pln,[]));