% Example: Proton Treatment Plan with Manipulated CT values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Copyright 2017 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. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0001 %% Example: Proton Treatment Plan with Manipulated CT values 0002 % 0003 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0004 % 0005 % Copyright 2017 the matRad development team. 0006 % 0007 % This file is part of the matRad project. It is subject to the license 0008 % terms in the LICENSE file found in the top-level directory of this 0009 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part 0010 % of the matRad project, including this file, may be copied, modified, 0011 % propagated, or distributed except according to the terms contained in the 0012 % LICENSE file. 0013 % 0014 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0015 0016 %% 0017 % In this example we will show 0018 % (i) how to load patient data into matRad 0019 % (ii) how to setup a proton dose calculation 0020 % (iii) how to inversely optimize the pencil beam intensities directly from command window in MATLAB. 0021 % (iv) how to re-optimize a treatment plan 0022 % (v) how to manipulate the CT cube by adding noise to the cube 0023 % (vi) how to recalculate the dose considering the manipulated CT cube and the previously optimized pencil beam intensities 0024 % (vii) how to compare the two results 0025 0026 %% Patient Data Import 0027 % Let's begin with a clear Matlab environment and import the prostate 0028 % patient into your workspace. 0029 0030 matRad_rc; %If this throws an error, run it from the parent directory first to set the paths 0031 0032 load('PROSTATE.mat'); 0033 0034 %% Treatment Plan 0035 % The next step is to define your treatment plan labeled as 'pln'. This 0036 % structure requires input from the treatment planner and defines 0037 % the most important cornerstones of your treatment plan. 0038 0039 pln.radiationMode = 'protons'; 0040 pln.machine = 'generic_MCsquare'; 0041 pln.numOfFractions = 30; 0042 pln.propOpt.bioOptimization = 'const_RBExD'; 0043 pln.propStf.gantryAngles = [90 270]; 0044 pln.propStf.couchAngles = [0 0]; 0045 pln.propStf.bixelWidth = 3; 0046 pln.propStf.numOfBeams = numel(pln.propStf.gantryAngles); 0047 pln.propStf.isoCenter = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0); 0048 pln.propOpt.runDAO = 0; 0049 pln.propOpt.runSequencing = 0; 0050 0051 % dose calculation settings 0052 pln.propDoseCalc.doseGrid.resolution.x = 3; % [mm] 0053 pln.propDoseCalc.doseGrid.resolution.y = 3; % [mm] 0054 pln.propDoseCalc.doseGrid.resolution.z = 3; % [mm] 0055 0056 %% Generate Beam Geometry STF 0057 stf = matRad_generateStf(ct,cst,pln); 0058 0059 %% Dose Calculation 0060 dij = matRad_calcParticleDoseMC(ct,stf,pln,cst); 0061 0062 %% Inverse Optimization for IMPT 0063 resultGUI = matRad_fluenceOptimization(dij,cst,pln); 0064 0065 %% Calculate quality indicators 0066 [dvh,qi] = matRad_indicatorWrapper(cst,pln,resultGUI); 0067 ixRectum = 1; 0068 display(qi(ixRectum).D_5); 0069 0070 %% 0071 % Let's change the optimization parameter of the rectum in such a way that it 0072 % will be better spared. We increase the penalty and lower the threshold 0073 % of the squared overdose objective function. Afterwards we re-optimize 0074 % the treatment plan and evaluate dose statistics one more time. 0075 0076 objective = cst{ixRectum,6}{1}; %This gives a struct 0077 objective = matRad_DoseOptimizationFunction.createInstanceFromStruct(objective); %Now we turn it into a class 0078 objective = objective.setDoseParameters(40); %We can simply call this function to change the/all dose parameter(s) 0079 cst{ixRectum,6}{1} = struct(objective); % We put it back as struct 0080 0081 cst{ixRectum,6}{1}.parameters{1} = 40; 0082 cst{ixRectum,6}{1}.penalty = 500; 0083 resultGUI = matRad_fluenceOptimization(dij,cst,pln); 0084 [dvh2,qi2] = matRad_indicatorWrapper(cst,pln,resultGUI); 0085 display(qi2(ixRectum).D_5); 0086 0087 %% Plot the Resulting Dose Slice 0088 % Let's plot the transversal iso-center dose slice 0089 slice = round(pln.propStf.isoCenter(1,3)./ct.resolution.z); 0090 figure 0091 imagesc(resultGUI.RBExDose(:,:,slice)),colorbar, colormap(jet) 0092 0093 %% 0094 % Now let's simulate a range undershoot by scaling the relative stopping power cube by 3.5% percent 0095 ct_manip = ct; 0096 ct_manip.cubeHU{1} = 1.035*ct_manip.cubeHU{1}; 0097 0098 %% Recalculate Plan with MC square 0099 % Let's use the existing optimized pencil beam weights and recalculate the RBE weighted dose 0100 resultGUI_noise = matRad_calcDoseDirectMC(ct_manip,stf,pln,cst,resultGUI.w); 0101 0102 %% Visual Comparison of results 0103 % Let's compare the new recalculation against the optimization result. 0104 plane = 3; 0105 doseWindow = [0 max([resultGUI.RBExDose(:); resultGUI_noise.RBExDose(:)])]; 0106 0107 figure,title('original plan') 0108 matRad_plotSliceWrapper(gca,ct,cst,1,resultGUI.RBExDose,plane,slice,[],0.75,colorcube,[],doseWindow,[]); 0109 figure,title('manipulated plan') 0110 matRad_plotSliceWrapper(gca,ct_manip,cst,1,resultGUI_noise.RBExDose,plane,slice,[],0.75,colorcube,[],doseWindow,[]); 0111 0112 % Let's plot single profiles along the beam direction 0113 ixProfileY = round(pln.propStf.isoCenter(1,1)./ct.resolution.x); 0114 0115 profileOrginal = resultGUI.RBExDose(:,ixProfileY,slice); 0116 profileNoise = resultGUI_noise.RBExDose(:,ixProfileY,slice); 0117 0118 figure,plot(profileOrginal,'LineWidth',2),grid on,hold on, 0119 plot(profileNoise,'LineWidth',2),legend({'original profile','noise profile'}), 0120 xlabel('mm'),ylabel('Gy(RBE)'),title('profile plot') 0121 0122 %% Quantitative Comparison of results 0123 % Compare the two dose cubes using a gamma-index analysis. 0124 0125 doseDifference = 2; 0126 distToAgreement = 2; 0127 n = 1; 0128 0129 [gammaCube,gammaPassRateCell] = matRad_gammaIndex(... 0130 resultGUI_noise.RBExDose,resultGUI.RBExDose,... 0131 [ct.resolution.x, ct.resolution.y, ct.resolution.z],... 0132 [doseDifference distToAgreement],slice,n,'global',cst); 0133