% Example: Photon Treatment Plan %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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: Photon Treatment Plan 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 photon dose calculation and 0020 % (iii) how to inversely optimize beamlet intensities 0021 % (iv) how to visually and quantitatively evaluate the result 0022 0023 %% Patient Data Import 0024 % Let's begin with a clear Matlab environment. Then, import the TG119 0025 % phantom into your workspace. The phantom is comprised of a 'ct' and 'cst' 0026 % structure defining the CT images and the structure set. Make sure the 0027 % matRad root directory with all its subdirectories is added to the Matlab 0028 % search path. 0029 0030 matRad_rc; %If this throws an error, run it from the parent directory first to set the paths 0031 0032 load('TG119.mat'); 0033 0034 %% 0035 % The file TG119.mat contains two Matlab variables. Let's check what we 0036 % have just imported. First, the 'ct' variable comprises the ct cube along 0037 %with some meta information describing properties of the ct cube (cube 0038 % dimensions, resolution, number of CT scenarios). Please note that 0039 %multiple ct cubes (e.g. 4D CT) can be stored in the cell array ct.cube{} 0040 display(ct); 0041 0042 %% 0043 % The 'cst' cell array defines volumes of interests along with information 0044 % required for optimization. Each row belongs to one certain volume of 0045 % interest (VOI), whereas each column defines different properties. 0046 % Specifically, the second and third column show the name and the type of 0047 % the structure. The type can be set to OAR, TARGET or IGNORED. The fourth 0048 % column contains a linear index vector that lists all voxels belonging to 0049 % a certain VOI. 0050 display(cst); 0051 0052 %% 0053 % The fifth column represents meta parameters for optimization. The overlap 0054 % priority is used to resolve ambiguities of overlapping structures (voxels 0055 % belonging to multiple structures will only be assigned to the VOI(s) with 0056 % the highest overlap priority, i.e.. the lowest value). The parameters 0057 % alphaX and betaX correspond to the tissue's photon-radiosensitivity 0058 % parameter of the linear quadratic model. These parameter are required for 0059 % biological treatment planning using a variable RBE. Let's output the meta 0060 % optimization parameter of the target, which is stored in the thrid row: 0061 ixTarget = 3; 0062 display(cst{ixTarget,5}); 0063 0064 %% 0065 % The sixth column contains optimization information such as objectives and 0066 % constraints which are required to calculate the objective function value. 0067 % Please note, that multiple objectives/constraints can be defined for 0068 % individual structures. Here, we have defined a squared deviation 0069 % objective making it 'expensive/costly' for the optimizer to over- and 0070 % underdose the target structure (both are equally important). 0071 display(cst{ixTarget,6}); 0072 0073 %% Treatment Plan 0074 % The next step is to define your treatment plan labeled as 'pln'. This 0075 % matlab structure requires input from the treatment planner and defines 0076 % the most important cornerstones of your treatment plan. 0077 0078 %% 0079 % First of all, we need to define what kind of radiation modality we would 0080 % like to use. Possible values are photons, protons or carbon. In this case 0081 % we want to use photons. Then, we need to define a treatment machine to 0082 % correctly load the corresponding base data. matRad includes base data for 0083 % generic photon linear accelerator called 'Generic'. By this means matRad 0084 % will look for 'photons_Generic.mat' in our root directory and will use 0085 % the data provided in there for dose calculation 0086 pln.radiationMode = 'photons'; 0087 pln.machine = 'Generic'; 0088 0089 %% 0090 % Define the flavor of optimization along with the quantity that should be 0091 % used for optimization. Possible values are (none: physical optimization; 0092 % const_RBExD: constant RBE of 1.1; LEMIV_effect: effect-based 0093 % optimization; LEMIV_RBExD: optimization of RBE-weighted dose. As we are 0094 % using photons, we simply set the parameter to 'none' thereby indicating 0095 % the physical dose should be optimized. 0096 pln.propOpt.bioOptimization = 'none'; 0097 0098 %% 0099 % Now we have to set some beam parameters. We can define multiple beam 0100 % angles for the treatment and pass these to the plan as a vector. matRad 0101 % will then interpret the vector as multiple beams. In this case, we define 0102 % linear spaced beams from 0 degree to 359 degree in 40 degree steps. This 0103 % results in 9 beams. All corresponding couch angles are set to 0 at this 0104 % point. Moreover, we set the bixelWidth to 5, which results in a beamlet 0105 % size of 5 x 5 mm in the isocenter plane. The number of fractions is set 0106 % to 30. Internally, matRad considers the fraction dose for optimization, 0107 % however, objetives and constraints are defined for the entire treatment. 0108 pln.numOfFractions = 30; 0109 pln.propStf.gantryAngles = [0:40:359]; 0110 pln.propStf.couchAngles = zeros(1,numel(pln.propStf.gantryAngles)); 0111 pln.propStf.bixelWidth = 5; 0112 0113 %% 0114 % Obtain the number of beams and voxels from the existing variables and 0115 % calculate the iso-center which is per default the center of gravity of 0116 % all target voxels. 0117 pln.propStf.numOfBeams = numel(pln.propStf.gantryAngles); 0118 pln.propStf.isoCenter = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0); 0119 0120 %% dose calculation settings 0121 % set resolution of dose calculation and optimization 0122 pln.propDoseCalc.doseGrid.resolution.x = 3; % [mm] 0123 pln.propDoseCalc.doseGrid.resolution.y = 3; % [mm] 0124 pln.propDoseCalc.doseGrid.resolution.z = 3; % [mm] 0125 0126 %% 0127 % Enable sequencing and disable direct aperture optimization (DAO) for now. 0128 % A DAO optimization is shown in a seperate example. 0129 pln.propOpt.runSequencing = 1; 0130 pln.propOpt.runDAO = 0; 0131 0132 %% 0133 % and et voila our treatment plan structure is ready. Lets have a look: 0134 display(pln); 0135 0136 %% Generate Beam Geometry STF 0137 % The steering file struct comprises the complete beam geometry along with 0138 % ray position, pencil beam positions and energies, source to axis distance (SAD) etc. 0139 stf = matRad_generateStf(ct,cst,pln); 0140 0141 %% 0142 % Let's display the beam geometry information of the 6th beam 0143 display(stf(6)); 0144 0145 %% Dose Calculation 0146 % Let's generate dosimetric information by pre-computing dose influence 0147 % matrices for unit beamlet intensities. Having dose influences available 0148 % allows subsequent inverse optimization. 0149 dij = matRad_calcPhotonDose(ct,stf,pln,cst); 0150 0151 %% Inverse Optimization for IMRT 0152 % The goal of the fluence optimization is to find a set of beamlet/pencil 0153 % beam weights which yield the best possible dose distribution according to 0154 % the clinical objectives and constraints underlying the radiation 0155 % treatment. Once the optimization has finished, trigger once the GUI to 0156 % visualize the optimized dose cubes. 0157 resultGUI = matRad_fluenceOptimization(dij,cst,pln); 0158 matRadGUI; 0159 0160 %% Plot the Resulting Dose Slice 0161 % Let's plot the transversal iso-center dose slice 0162 slice = round(pln.propStf.isoCenter(1,3)./ct.resolution.z); 0163 figure 0164 imagesc(resultGUI.physicalDose(:,:,slice)),colorbar, colormap(jet); 0165 0166 %% Now let's create another treatment plan but this time use a coarser beam spacing. 0167 % Instead of 40 degree spacing use a 50 degree geantry beam spacing 0168 pln.propStf.gantryAngles = [0:50:359]; 0169 pln.propStf.couchAngles = zeros(1,numel(pln.propStf.gantryAngles)); 0170 pln.propStf.numOfBeams = numel(pln.propStf.gantryAngles); 0171 stf = matRad_generateStf(ct,cst,pln); 0172 pln.propStf.isoCenter = vertcat(stf.isoCenter); 0173 dij = matRad_calcPhotonDose(ct,stf,pln,cst); 0174 resultGUI_coarse = matRad_fluenceOptimization(dij,cst,pln); 0175 0176 %% Visual Comparison of results 0177 % Let's compare the new recalculation against the optimization result. 0178 % Check if you have added all subdirectories to the Matlab search path, 0179 % otherwise it will not find the plotting function 0180 plane = 3; 0181 doseWindow = [0 max([resultGUI.physicalDose(:); resultGUI_coarse.physicalDose(:)])]; 0182 0183 figure,title('original plan - fine beam spacing') 0184 matRad_plotSliceWrapper(gca,ct,cst,1,resultGUI.physicalDose,plane,slice,[],0.75,colorcube,[],doseWindow,[]); 0185 figure,title('modified plan - coarse beam spacing') 0186 matRad_plotSliceWrapper(gca,ct,cst,1,resultGUI_coarse.physicalDose,plane,slice,[],0.75,colorcube,[],doseWindow,[]); 0187 0188 %% 0189 % At this point we would like to see the absolute difference of the first 0190 % optimization (finer beam spacing) and the second optimization (coarser 0191 % beam spacing) 0192 absDiffCube = resultGUI.physicalDose-resultGUI_coarse.physicalDose; 0193 figure,title( 'fine beam spacing plan - coarse beam spacing plan') 0194 matRad_plotSliceWrapper(gca,ct,cst,1,absDiffCube,plane,slice,[],[],colorcube); 0195 0196 %% Obtain dose statistics 0197 % Two more columns will be added to the cst structure depicting the DVH and 0198 % standard dose statistics such as D95,D98, mean dose, max dose etc. 0199 [dvh,qi] = matRad_indicatorWrapper(cst,pln,resultGUI); 0200 [dvh_coarse,qi_coarse] = matRad_indicatorWrapper(cst,pln,resultGUI_coarse); 0201 0202 %% 0203 % The treatment plan using more beams should in principle result in a 0204 % better OAR sparing. Therefore lets have a look at the D95 of the OAR of 0205 % both plans 0206 ixOAR = 2; 0207 display(qi(ixOAR).D_95); 0208 display(qi_coarse(ixOAR).D_95); 0209 0210 0211 %% 0212 % Export the dose in binary format. matRad can write multiple formats, here we 0213 % choose to export in nrrd format using the matRad_writeCube function, which 0214 % chooses the appropriate subroutine from the extension. The metadata struct can 0215 % also store more optional parameters, but requires only the resolution to be s 0216 % set. 0217 metadata = struct('resolution',[ct.resolution.x ct.resolution.y ct.resolution.z]); 0218 matRad_writeCube([pwd filesep 'photonDose_example2.nrrd'],resultGUI.physicalDose,'double',metadata); 0219 0220