% Example: Proton Treatment Plan with subsequent Isocenter shift %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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 subsequent Isocenter shift 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 simulate a lateral patient displacement by shifting the iso-center 0022 % (v) how to recalculated the dose considering the shifted geometry and the previously optimized pencil beam intensities 0023 % (vi) how to compare the two results 0024 0025 %% Patient Data Import 0026 % Let's begin with a clear Matlab environment and import the prostate 0027 % patient into your workspace 0028 0029 matRad_rc; %If this throws an error, run it from the parent directory first to set the paths 0030 0031 load('PROSTATE.mat'); 0032 0033 %% Treatment Plan 0034 % The next step is to define your treatment plan labeled as 'pln'. This 0035 % structure requires input from the treatment planner and defines the most 0036 % important cornerstones of your treatment plan. 0037 0038 %% 0039 % First of all, we need to define what kind of radiation modality we would 0040 % like to use. Possible values are photons, protons or carbon. In this 0041 % example we would like to use protons for treatment planning. Next, we 0042 % need to define a treatment machine to correctly load the corresponding 0043 % base data. matRad features generic base data in the file 0044 % 'proton_Generic.mat'; consequently the machine has to be set accordingly 0045 pln.radiationMode = 'protons'; 0046 pln.machine = 'Generic'; 0047 0048 %% 0049 % Define the flavor of biological optimization for treatment planning along 0050 % with the quantity that should be used for optimization. Possible values 0051 % are (none: physical optimization; const_RBExD: constant RBE of 1.1; 0052 % LEMIV_effect: effect-based optimization; LEMIV_RBExD: optimization of 0053 % RBE-weighted dose. As we use protons, we follow here the clinical 0054 % standard and use a constant relative biological effectiveness of 1.1. 0055 % Therefore we set bioOptimization to const_RBExD 0056 pln.propOpt.bioOptimization = 'const_RBExD'; 0057 0058 %% 0059 % for particles it is possible to also calculate the LET disutribution 0060 % alongside the physical dose. Therefore you need to activate the 0061 % corresponding option during dose calculcation 0062 pln.propDoseCalc.calcLET = 1; 0063 0064 %% 0065 % Now we have to set the remaining plan parameters. 0066 pln.numOfFractions = 30; 0067 pln.propStf.gantryAngles = [90 270]; 0068 pln.propStf.couchAngles = [0 0]; 0069 pln.propStf.bixelWidth = 3; 0070 pln.propStf.numOfBeams = numel(pln.propStf.gantryAngles); 0071 pln.propStf.isoCenter = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0); 0072 pln.propOpt.runDAO = 0; 0073 pln.propOpt.runSequencing = 0; 0074 0075 % dose calculation settings 0076 pln.propDoseCalc.doseGrid.resolution.x = 3; % [mm] 0077 pln.propDoseCalc.doseGrid.resolution.y = 3; % [mm] 0078 pln.propDoseCalc.doseGrid.resolution.z = 3; % [mm] 0079 0080 %% Generate Beam Geometry STF 0081 stf = matRad_generateStf(ct,cst,pln); 0082 0083 %% Dose Calculation 0084 % Lets generate dosimetric information by pre-computing dose influence 0085 % matrices for unit beamlet intensities. Having dose influences available 0086 % allows for subsequent inverse optimization. 0087 dij = matRad_calcParticleDose(ct,stf,pln,cst); 0088 0089 %% Inverse Optimization for IMPT 0090 % The goal of the fluence optimization is to find a set of bixel/spot 0091 % weights which yield the best possible dose distribution according to the 0092 % clinical objectives and constraints underlying the radiation treatment 0093 resultGUI = matRad_fluenceOptimization(dij,cst,pln); 0094 0095 %% Plot the Resulting Dose Slice 0096 % Let's plot the transversal iso-center dose slice 0097 slice = round(pln.propStf.isoCenter(1,3)./ct.resolution.z); 0098 figure 0099 imagesc(resultGUI.RBExDose(:,:,slice)),colorbar,colormap(jet) 0100 0101 %% Plot the Resulting Beam Dose Slice 0102 % Let's plot the transversal iso-center dose slice of beam 1 and beam 2 0103 % separately 0104 figure 0105 subplot(121),imagesc(resultGUI.RBExDose_beam1(:,:,slice)),colorbar,colormap(jet),title('dose of beam 1') 0106 subplot(122),imagesc(resultGUI.RBExDose_beam2(:,:,slice)),colorbar,colormap(jet),title('dose of beam 2') 0107 0108 %% and the corresponding LET distribution 0109 % Transversal iso-center slice 0110 figure 0111 imagesc(resultGUI.LET(:,:,slice)),colormap(jet),colorbar,title('LET [keV/µm]') 0112 0113 %% 0114 % Now let's simulate a patient shift in y direction for both beams 0115 stf(1).isoCenter(2) = stf(1).isoCenter(2) - 4; 0116 stf(2).isoCenter(2) = stf(2).isoCenter(2) - 4; 0117 pln.propStf.isoCenter = reshape([stf.isoCenter],[3 pln.propStf.numOfBeams])'; 0118 0119 %% Recalculate Plan 0120 % Let's use the existing optimized pencil beam weights and recalculate the RBE weighted dose 0121 resultGUI_isoShift = matRad_calcDoseDirect(ct,stf,pln,cst,resultGUI.w); 0122 0123 %% Visual Comparison of results 0124 % Let's compare the new recalculation against the optimization result. 0125 plane = 3; 0126 doseWindow = [0 max([resultGUI.RBExDose(:); resultGUI_isoShift.RBExDose(:)])]; 0127 0128 figure,title('original plan') 0129 matRad_plotSliceWrapper(gca,ct,cst,1,resultGUI.RBExDose,plane,slice,[],0.75,colorcube,[],doseWindow,[]); 0130 figure,title('shifted plan') 0131 matRad_plotSliceWrapper(gca,ct,cst,1,resultGUI_isoShift.RBExDose,plane,slice,[],0.75,colorcube,[],doseWindow,[]); 0132 0133 absDiffCube = resultGUI.RBExDose-resultGUI_isoShift.RBExDose; 0134 figure,title('absolute difference') 0135 matRad_plotSliceWrapper(gca,ct,cst,1,absDiffCube,plane,slice,[],[],colorcube); 0136 0137 % Let's plot single profiles that are perpendicular to the beam direction 0138 ixProfileY = round(pln.propStf.isoCenter(1,2)./ct.resolution.y); 0139 0140 profileOrginal = resultGUI.RBExDose(:,ixProfileY,slice); 0141 profileShifted = resultGUI_isoShift.RBExDose(:,ixProfileY,slice); 0142 0143 figure,plot(profileOrginal,'LineWidth',2),grid on,hold on, 0144 plot(profileShifted,'LineWidth',2),legend({'original profile','shifted profile'}), 0145 xlabel('mm'),ylabel('Gy(RBE)'),title('profile plot') 0146 0147 %% Quantitative Comparison of results 0148 % Compare the two dose cubes using a gamma-index analysis. The gamma index 0149 % is a composite quality distribution equally taking into account a dose 0150 % difference and a distance to agreement criterion in order to quantify 0151 % differences between two dose cubes. A gamma-index value of smaller than 1 0152 % indicates a successful test and a value greater than 1 illustrates a 0153 % failed test. 0154 0155 doseDifference = 2; 0156 distToAgreement = 2; 0157 n = 1; 0158 0159 [gammaCube,gammaPassRateCell] = matRad_gammaIndex(... 0160 resultGUI_isoShift.RBExDose,resultGUI.RBExDose,... 0161 [ct.resolution.x, ct.resolution.y, ct.resolution.z],... 0162 [doseDifference distToAgreement],slice,n,'global',cst); 0163 0164 0165 [env, ~] = matRad_getEnvironment(); 0166 % Let's plot the gamma index histogram 0167 figure 0168 hist(gammaCube(gammaCube>0),100) 0169 title('gamma index histogram') 0170 0171 0172