From 87fca8f664486858410527e934343e1fede706ad Mon Sep 17 00:00:00 2001 From: JenHardt Date: Wed, 6 Dec 2023 15:35:27 +0100 Subject: [PATCH 1/5] Update: dose Accumulation x and y direction were switched in the dose accumulation and in the dvf generation, fixed that. Added dvfType and doseeAccumulation Type as varriable --- 4D/matRad_addMovement.m | 22 ++++++++++++++++------ 4D/matRad_calc4dDose.m | 19 ++++++++++++------- 4D/matRad_doseAcc.m | 20 ++++++++++---------- tools/matRad_compareDose.m | 2 +- 4 files changed, 39 insertions(+), 24 deletions(-) diff --git a/4D/matRad_addMovement.m b/4D/matRad_addMovement.m index 9f075deeb..46a29a421 100644 --- a/4D/matRad_addMovement.m +++ b/4D/matRad_addMovement.m @@ -1,4 +1,4 @@ -function [ct, cst] = matRad_addMovement(ct, cst, motionPeriod, numOfCtScen, amp,visBool) +function [ct, cst] = matRad_addMovement(ct, cst, motionPeriod, numOfCtScen, amp, dvfType, visBool) % adds artificial sinosodal patient motion by creating a deformation vector % field and applying it to the ct.cube by geometric transformation % @@ -11,6 +11,7 @@ % motionPeriod: the length of a whole breathing cycle (in seconds) % numOfCtScen: number of ct phases % amp: amplitude of the sinosoidal movement (in pixels) +% dvfType: push or pull dvf % visBool boolean flag for visualization % % note: 1st dim --> x LPS coordinate system @@ -43,6 +44,9 @@ if ~exist('visBool','var') visBool = false; end +if ~exist('dvfType','var') + dvfType = 'pull'; +end matRad_cfg = MatRad_Config.instance(); @@ -50,8 +54,15 @@ ct.motionPeriod = motionPeriod; ct.numOfCtScen = numOfCtScen; +if ~(strcmp(dvfType, 'pull') || strcmp(dvfType,'push')) + matRad_cfg.disperror('Chosse Push or Pull DVF type') +end + % set type -ct.dvfType = 'pull'; % push or pull +ct.dvfMetadata.dvfType = dvfType; +if strcmp(dvfType,'push') + amp = -amp; %dvf_pull = -dvf_push; +end env = matRad_getEnvironment(); @@ -110,10 +121,9 @@ end % convert dvfs to [mm] - tmp = ct.dvf{i}(:,:,:,1); - ct.dvf{i}(:,:,:,1) = -ct.dvf{i}(:,:,:,2) * ct.resolution.x; - ct.dvf{i}(:,:,:,2) = -tmp * ct.resolution.y; - ct.dvf{i}(:,:,:,3) = -ct.dvf{i}(:,:,:,3) * ct.resolution.z; + ct.dvf{i}(:,:,:,1) = ct.dvf{i}(:,:,:,1).* ct.resolution.x; + ct.dvf{i}(:,:,:,2) = ct.dvf{i}(:,:,:,2).*ct.resolution.y; + ct.dvf{i}(:,:,:,3) = ct.dvf{i}(:,:,:,3).* ct.resolution.z; ct.dvf{i} = permute(ct.dvf{i}, [4,1,2,3]); end diff --git a/4D/matRad_calc4dDose.m b/4D/matRad_calc4dDose.m index 1f783e492..c54723d5a 100644 --- a/4D/matRad_calc4dDose.m +++ b/4D/matRad_calc4dDose.m @@ -1,4 +1,4 @@ -function [resultGUI, timeSequence] = matRad_calc4dDose(ct, pln, dij, stf, cst, resultGUI, totalPhaseMatrix) +function [resultGUI, timeSequence] = matRad_calc4dDose(ct, pln, dij, stf, cst, resultGUI, totalPhaseMatrix,accType) % wrapper for the whole 4D dose calculation pipeline and calculated dose % accumulation % @@ -13,6 +13,7 @@ % cst: matRad cst struct % resultGUI: struct containing optimized fluence vector % totalPhaseMatrix optional intput for totalPhaseMatrix +% accType: witch algorithim for dose accumulation % output % resultGUI: structure containing phase dose, RBE weighted dose, etc % timeSequence: timing information about the irradiation @@ -33,6 +34,10 @@ % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +if ~exist('accType','var') + accType = 'DDM'; +end + if ~exist('totalPhaseMatrix','var') % make a time sequence for when each bixel is irradiated, the sequence % follows the backforth spot scanning @@ -62,7 +67,7 @@ elseif strcmp(pln.bioParam.model,'constRBE') resultGUI.phaseRBExD{i} = tmpResultGUI.RBExD; % compute all fields - elseif any(strcmp(pln.bioParam.model,{'MCN','LEM','WED'})) + elseif any(strcmp(pln.bioParam.model,{'MCN','LEM','WED','HEL'})) resultGUI.phaseAlphaDose{i} = tmpResultGUI.alpha .* tmpResultGUI.physicalDose; resultGUI.phaseSqrtBetaDose{i} = sqrt(tmpResultGUI.beta) .* tmpResultGUI.physicalDose; end @@ -72,16 +77,16 @@ % accumulation if strcmp(pln.bioParam.model,'none') - resultGUI.accPhysicalDose = matRad_doseAcc(ct,resultGUI.phaseDose, cst, 'DDM'); + resultGUI.accPhysicalDose = matRad_doseAcc(ct,resultGUI.phaseDose, cst, accType); elseif strcmp(pln.bioParam.model,'constRBE') - resultGUI.accRBExD = matRad_doseAcc(ct,resultGUI.phaseRBExD, cst, 'DDM'); + resultGUI.accRBExD = matRad_doseAcc(ct,resultGUI.phaseRBExD, cst, accType); -elseif any(strcmp(pln.bioParam.model,{'MCN','LEM','WED'})) +elseif any(strcmp(pln.bioParam.model,{'MCN','LEM','WED','HEL'})) - resultGUI.accAlphaDose = matRad_doseAcc(ct,resultGUI.phaseAlphaDose, cst, 'DDM'); - resultGUI.accSqrtBetaDose = matRad_doseAcc(ct,resultGUI.phaseSqrtBetaDose, cst, 'DDM'); + resultGUI.accAlphaDose = matRad_doseAcc(ct,resultGUI.phaseAlphaDose, cst,accType); + resultGUI.accSqrtBetaDose = matRad_doseAcc(ct,resultGUI.phaseSqrtBetaDose, cst, accType); % only compute where we have biologically defined tissue ix = dij.ax(:,1)~=0; diff --git a/4D/matRad_doseAcc.m b/4D/matRad_doseAcc.m index b1119a603..6e4eca84a 100644 --- a/4D/matRad_doseAcc.m +++ b/4D/matRad_doseAcc.m @@ -57,7 +57,7 @@ if strcmp(accMethod,'DDM') - if ~strcmp(ct.dvfType,'pull') + if ~strcmp(ct.dvfMetadata.dvfType,'pull') error('dose accumulation via direct dose mapping (DDM) requires pull dvfs'); end @@ -75,19 +75,16 @@ dvf_y_i = squeeze(ct.dvf{1,i}(2,:,:,:))/ct.resolution.y; dvf_z_i = squeeze(ct.dvf{1,i}(3,:,:,:))/ct.resolution.z; - d_ref = matRad_interp3(yGridVec,xGridVec,zGridVec,permute(phaseCubes{i},[2 1 3]), ... - Y(ix) + dvf_y_i(ix), ... - X(ix) + dvf_x_i(ix), ... - Z(ix) + dvf_z_i(ix), ... - 'linear',0); + d_ref = matRad_interp3(X, Y, Z,... + phaseCubes{i}, X-dvf_x_i,Y-dvf_y_i,Z-dvf_z_i,'linear',0); - dAcc(ix) = dAcc(ix) + d_ref; + dAcc(ix) = dAcc(ix) + d_ref(ix); end elseif strcmp(accMethod,'EMT') % funktioniert nicht wenn Dosis in einer Phase = 0 ist... - if ~strcmp(ct.dvfType,'push') + if ~strcmp(ct.dvfMetadata.dvfType,'push') error('dose accumulation via interpolation requires push dvfs'); end @@ -99,8 +96,10 @@ dvf_y_i = squeeze(ct.dvf{1,i}(2,:,:,:))/ct.resolution.y; dvf_z_i = squeeze(ct.dvf{1,i}(3,:,:,:))/ct.resolution.z; - m_i = ct.cube{i}; - e_i = phaseCubes{i}.*m_i; + m_i = ct.cube{i}; + m_i = permute(m_i,[2 1 3]); + e_i = phaseCubes{i}.*ct.cube{i}; + e_i = permute(e_i,[2 1 3]); ix = e_i>0; @@ -163,6 +162,7 @@ k = find(m_ref); dAcc(k) = dAcc(k) + e_ref(k)./m_ref(k); + dAcc = permute(dAcc,[2 1 3]); end % elseif strcmp(accMethod,'DDMM') diff --git a/tools/matRad_compareDose.m b/tools/matRad_compareDose.m index 876f1d3a1..bdb6ccf8c 100644 --- a/tools/matRad_compareDose.m +++ b/tools/matRad_compareDose.m @@ -142,7 +142,7 @@ % Calculate absolute difference cube and dose windows for plots differenceCube = cube1-cube2; - doseDiffWindow = [-max(differenceCube(:)) max(differenceCube(:))]; + doseDiffWindow = [-max(abs(differenceCube(:))) max(abs(differenceCube(:)))]; %doseGammaWindow = [0 max(gammaCube(:))]; doseGammaWindow = [0 2]; %We choose 2 as maximum value since the gamma colormap has a sharp cut in the middle From 57495ad3558e108bd4d462fecca3f10d5e3acefe Mon Sep 17 00:00:00 2001 From: JenHardt Date: Thu, 7 Dec 2023 08:49:07 +0100 Subject: [PATCH 2/5] update to pull request updated example and error message --- 4D/matRad_addMovement.m | 2 +- examples/matRad_example9_4DDoseCalcMinimal.m | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/4D/matRad_addMovement.m b/4D/matRad_addMovement.m index 46a29a421..c1a30ffb3 100644 --- a/4D/matRad_addMovement.m +++ b/4D/matRad_addMovement.m @@ -55,7 +55,7 @@ ct.numOfCtScen = numOfCtScen; if ~(strcmp(dvfType, 'pull') || strcmp(dvfType,'push')) - matRad_cfg.disperror('Chosse Push or Pull DVF type') + matRad_cfg.dispError('Chosse Push or Pull DVF type') end % set type diff --git a/examples/matRad_example9_4DDoseCalcMinimal.m b/examples/matRad_example9_4DDoseCalcMinimal.m index 816ea1cbf..f2e596822 100644 --- a/examples/matRad_example9_4DDoseCalcMinimal.m +++ b/examples/matRad_example9_4DDoseCalcMinimal.m @@ -29,7 +29,7 @@ numOfCtScen = 5; motionPeriod = 2.5; % [s] -[ct,cst] = matRad_addMovement(ct, cst,motionPeriod, numOfCtScen, amplitude); +[ct,cst] = matRad_addMovement(ct, cst,motionPeriod, numOfCtScen, amplitude,'pull'); % Set up a plan, compute dose influence on all phases, conventional optimization % meta information for treatment plan pln.numOfFractions = 30; From f15149c593e4d5e402c383e582ba024cc3b813f0 Mon Sep 17 00:00:00 2001 From: JenHardt Date: Thu, 7 Dec 2023 09:49:06 +0100 Subject: [PATCH 3/5] Update matRad_example10_4DphotonRobust.m --- examples/matRad_example10_4DphotonRobust.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/matRad_example10_4DphotonRobust.m b/examples/matRad_example10_4DphotonRobust.m index 13211ab22..ad1949fb1 100644 --- a/examples/matRad_example10_4DphotonRobust.m +++ b/examples/matRad_example10_4DphotonRobust.m @@ -122,7 +122,7 @@ numOfCtScen = 5; motionPeriod = 2.5; % [s] -[ct,cst] = matRad_addMovement(ct, cst,motionPeriod, numOfCtScen, amplitude,1); +[ct,cst] = matRad_addMovement(ct, cst,motionPeriod, numOfCtScen, amplitude,'pull',1); % show the deformation vector field slice = 25; % select a specific slice and to plot the vector field From 31a773dc8f278fcd874a882f7690656c6ef75f15 Mon Sep 17 00:00:00 2001 From: JenHardt Date: Mon, 11 Dec 2023 15:35:11 +0100 Subject: [PATCH 4/5] typo fix --- 4D/matRad_addMovement.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/4D/matRad_addMovement.m b/4D/matRad_addMovement.m index c1a30ffb3..35d40161a 100644 --- a/4D/matRad_addMovement.m +++ b/4D/matRad_addMovement.m @@ -55,7 +55,7 @@ ct.numOfCtScen = numOfCtScen; if ~(strcmp(dvfType, 'pull') || strcmp(dvfType,'push')) - matRad_cfg.dispError('Chosse Push or Pull DVF type') + matRad_cfg.dispError('Choose Push or Pull DVF type') end % set type From 16e064df84c60e903de657d2000986893623c773 Mon Sep 17 00:00:00 2001 From: JenHardt Date: Tue, 12 Dec 2023 08:55:12 +0100 Subject: [PATCH 5/5] input parser update --- 4D/matRad_addMovement.m | 28 +++++++++----------- examples/matRad_example10_4DphotonRobust.m | 2 +- examples/matRad_example9_4DDoseCalcMinimal.m | 2 +- 3 files changed, 14 insertions(+), 18 deletions(-) diff --git a/4D/matRad_addMovement.m b/4D/matRad_addMovement.m index 35d40161a..93fda54b9 100644 --- a/4D/matRad_addMovement.m +++ b/4D/matRad_addMovement.m @@ -1,4 +1,4 @@ -function [ct, cst] = matRad_addMovement(ct, cst, motionPeriod, numOfCtScen, amp, dvfType, visBool) +function [ct, cst] = matRad_addMovement(ct, cst, motionPeriod, numOfCtScen, amp, varargin) % adds artificial sinosodal patient motion by creating a deformation vector % field and applying it to the ct.cube by geometric transformation % @@ -11,8 +11,8 @@ % motionPeriod: the length of a whole breathing cycle (in seconds) % numOfCtScen: number of ct phases % amp: amplitude of the sinosoidal movement (in pixels) -% dvfType: push or pull dvf -% visBool boolean flag for visualization +% varargin: dvfType: push or pull dvf +% visBool boolean flag for visualization % % note: 1st dim --> x LPS coordinate system % 2nd dim --> y LPS coordinate system @@ -41,12 +41,12 @@ % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -if ~exist('visBool','var') - visBool = false; -end -if ~exist('dvfType','var') - dvfType = 'pull'; -end +expectedDVF = {'pull', 'push'}; + +p = inputParser; +addParameter(p,'dvfType','pull', @(x) any(validatestring(x,expectedDVF))) +addParameter(p,'visBool',false, @islogical); +parse(p,varargin{:}); matRad_cfg = MatRad_Config.instance(); @@ -54,13 +54,9 @@ ct.motionPeriod = motionPeriod; ct.numOfCtScen = numOfCtScen; -if ~(strcmp(dvfType, 'pull') || strcmp(dvfType,'push')) - matRad_cfg.dispError('Choose Push or Pull DVF type') -end - % set type -ct.dvfMetadata.dvfType = dvfType; -if strcmp(dvfType,'push') +ct.dvfMetadata.dvfType = p.Results.dvfType; +if strcmp(p.Results.dvfType,'push') amp = -amp; %dvf_pull = -dvf_push; end @@ -130,7 +126,7 @@ -if visBool +if p.Results.visBool slice = round(ct.cubeDim(3)/2); figure, for i = 1:numOfCtScen diff --git a/examples/matRad_example10_4DphotonRobust.m b/examples/matRad_example10_4DphotonRobust.m index ad1949fb1..d611843ac 100644 --- a/examples/matRad_example10_4DphotonRobust.m +++ b/examples/matRad_example10_4DphotonRobust.m @@ -122,7 +122,7 @@ numOfCtScen = 5; motionPeriod = 2.5; % [s] -[ct,cst] = matRad_addMovement(ct, cst,motionPeriod, numOfCtScen, amplitude,'pull',1); +[ct,cst] = matRad_addMovement(ct, cst,motionPeriod, numOfCtScen, amplitude,'visBool',true); % show the deformation vector field slice = 25; % select a specific slice and to plot the vector field diff --git a/examples/matRad_example9_4DDoseCalcMinimal.m b/examples/matRad_example9_4DDoseCalcMinimal.m index f2e596822..063bcdfff 100644 --- a/examples/matRad_example9_4DDoseCalcMinimal.m +++ b/examples/matRad_example9_4DDoseCalcMinimal.m @@ -29,7 +29,7 @@ numOfCtScen = 5; motionPeriod = 2.5; % [s] -[ct,cst] = matRad_addMovement(ct, cst,motionPeriod, numOfCtScen, amplitude,'pull'); +[ct,cst] = matRad_addMovement(ct, cst,motionPeriod, numOfCtScen, amplitude,'dvfType','pull'); % Set up a plan, compute dose influence on all phases, conventional optimization % meta information for treatment plan pln.numOfFractions = 30;