matRad_example1_phantom

% Example: Generate your own phantom geometry

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 Copyright 2018 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: Generate your own phantom geometry
0002 %
0003 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 % Copyright 2018 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 create arbitrary ct data (resolution, ct numbers)
0019 % (ii) how to create a cst structure containing the volume of interests of the phantom
0020 % (iii) generate a treatment plan for this phantom
0021 
0022 matRad_rc; %If this throws an error, run it from the parent directory first to set the paths
0023 
0024 %% Create a CT image series
0025 xDim = 200;
0026 yDim = 200;
0027 zDim = 50;
0028 
0029 ct.cubeDim      = [yDim xDim zDim]; % second cube dimension represents the x-coordinate
0030 ct.resolution.x = 2;
0031 ct.resolution.y = 2;
0032 ct.resolution.z = 3;
0033 ct.numOfCtScen  = 1;
0034  
0035 % create an ct image series with zeros - it will be filled later
0036 ct.cubeHU{1} = ones(ct.cubeDim) * -1000; % assign HU of Air
0037 
0038 %% Create the VOI data for the phantom
0039 % Now we define structures a contour for the phantom and a target
0040 
0041 ixOAR = 1;
0042 ixPTV = 2;
0043 
0044 % define general VOI properties
0045 cst{ixOAR,1} = 0;
0046 cst{ixOAR,2} = 'contour';
0047 cst{ixOAR,3} = 'OAR';
0048 
0049 cst{ixPTV,1} = 1;
0050 cst{ixPTV,2} = 'target';
0051 cst{ixPTV,3} = 'TARGET';
0052  
0053 % define optimization parameter for both VOIs
0054 cst{ixOAR,5}.TissueClass  = 1;
0055 cst{ixOAR,5}.alphaX       = 0.1000;
0056 cst{ixOAR,5}.betaX        = 0.0500;
0057 cst{ixOAR,5}.Priority     = 2;
0058 cst{ixOAR,5}.Visible      = 1;
0059 cst{ixOAR,5}.visibleColor = [0 0 0];
0060 
0061 % define objective as struct for compatibility with GNU Octave I/O
0062 cst{ixOAR,6}{1} = struct(DoseObjectives.matRad_SquaredOverdosing(10,30));
0063 
0064 cst{ixPTV,5}.TissueClass = 1;
0065 cst{ixPTV,5}.alphaX      = 0.1000;
0066 cst{ixPTV,5}.betaX       = 0.0500;
0067 cst{ixPTV,5}.Priority    = 1;
0068 cst{ixPTV,5}.Visible     = 1;
0069 cst{ixPTV,5}.visibleColor = [1 1 1];
0070 
0071 % define objective as struct for compatibility with GNU Octave I/O
0072 cst{ixPTV,6}{1} = struct(DoseObjectives.matRad_SquaredDeviation(800,60));
0073 
0074 %% Lets create either a cubic or a spheric phantom
0075 
0076 TYPE = 'spheric';   % either 'cubic' or 'spheric'
0077 
0078 % first the OAR
0079 cubeHelper = zeros(ct.cubeDim);
0080 
0081 switch TYPE
0082    
0083    case {'cubic'}
0084       
0085       xLowOAR  = round(xDim/2 - xDim/4);
0086       xHighOAR = round(xDim/2 + xDim/4);
0087       yLowOAR  = round(yDim/2 - yDim/4);
0088       yHighOAR = round(yDim/2 + yDim/4);
0089       zLowOAR  = round(zDim/2 - zDim/4);
0090       zHighOAR = round(zDim/2 + zDim/4);
0091       
0092       for x = xLowOAR:1:xHighOAR
0093          for y = yLowOAR:1:yHighOAR
0094             for z = zLowOAR:1:zHighOAR
0095                cubeHelper(y,x,z) = 1;
0096             end
0097          end
0098       end
0099       
0100    case {'spheric'}
0101       
0102       radiusOAR = xDim/4;
0103       
0104       for x = 1:xDim
0105          for y = 1:yDim
0106             for z = 1:zDim
0107                currPost = [y x z] - round([ct.cubeDim./2]);
0108                if  sqrt(sum(currPost.^2)) < radiusOAR
0109                   cubeHelper(y,x,z) = 1;
0110                end
0111             end
0112          end
0113       end
0114       
0115 end
0116 
0117 % extract the voxel indices and save it in the cst
0118 cst{ixOAR,4}{1} = find(cubeHelper);
0119 
0120 
0121 % second the PTV
0122 cubeHelper = zeros(ct.cubeDim);
0123 
0124 switch TYPE
0125    
0126    case {'cubic'}
0127       
0128       xLowPTV  = round(xDim/2 - xDim/8);
0129       xHighPTV = round(xDim/2 + xDim/8);
0130       yLowPTV  = round(yDim/2 - yDim/8);
0131       yHighPTV = round(yDim/2 + yDim/8);
0132       zLowPTV  = round(zDim/2 - zDim/8);
0133       zHighPTV = round(zDim/2 + zDim/8);
0134       
0135       cubeHelper = zeros(ct.cubeDim);
0136       
0137       for x = xLowPTV:1:xHighPTV
0138          for y = yLowPTV:1:yHighPTV
0139             for z = zLowPTV:1:zHighPTV
0140                cubeHelper(y,x,z) = 1;
0141             end
0142          end
0143       end
0144       
0145    case {'spheric'}
0146       
0147       radiusPTV = xDim/12;
0148       
0149       for x = 1:xDim
0150          for y = 1:yDim
0151             for z = 1:zDim
0152                currPost = [x y z] - round([ct.cubeDim./2]);
0153                if  sqrt(sum(currPost.^2)) < radiusPTV
0154                   cubeHelper(y,x,z) = 1;
0155                end
0156             end
0157          end
0158       end
0159       
0160 end
0161 
0162 
0163 
0164 % extract the voxel indices and save it in the cst
0165 cst{ixPTV,4}{1} = find(cubeHelper);
0166 
0167 
0168 % now we have ct data and cst data for a new phantom
0169 display(ct);
0170 display(cst);
0171 
0172 %% Assign relative electron densities
0173 vIxOAR = cst{ixOAR,4}{1};
0174 vIxPTV = cst{ixPTV,4}{1};
0175 
0176 ct.cubeHU{1}(vIxOAR) = 0;
0177 ct.cubeHU{1}(vIxPTV) = 0;
0178 
0179 %% Treatment Plan
0180 % The next step is to define your treatment plan labeled as 'pln'. This
0181 % structure requires input from the treatment planner and defines the most
0182 % important cornerstones of your treatment plan.
0183 %%
0184 % First of all, we need to define what kind of radiation modality we would
0185 % like to use. Possible values are photons, protons or carbon. In this
0186 % example we would like to use photons for treatment planning. Next, we
0187 % need to define a treatment machine to correctly load the corresponding
0188 % base data. matRad features generic base data in the file
0189 % 'photons_Generic.mat'; consequently the machine has to be set to 'Generic'
0190 pln.radiationMode = 'photons';            
0191 pln.machine       = 'Generic';
0192 
0193 %%
0194 % Define the flavor of biological optimization for treatment planning along
0195 % with the quantity that should be used for optimization. Possible values
0196 % are (none: physical dose based optimization; const_RBExD: constant RBE of 1.1;
0197 % LEMIV_effect: effect-based optimization; LEMIV_RBExD: optimization of
0198 % RBE-weighted dose. As we use photons, we select 'none' as we want to optimize the
0199 % physical dose.
0200 pln.propOpt.bioOptimization = 'none';                                              
0201 
0202 %%
0203 % The remaining plan parameters are set like in the previous example files
0204 pln.numOfFractions        = 30;
0205 pln.propStf.gantryAngles  = [0 45];
0206 pln.propStf.couchAngles   = [0 0];
0207 pln.propStf.bixelWidth    = 5;
0208 pln.propStf.numOfBeams    = numel(pln.propStf.gantryAngles);
0209 pln.propStf.isoCenter     = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0);
0210 pln.propOpt.runDAO        = 0;
0211 pln.propOpt.runSequencing = 0;
0212 
0213 % dose calculation settings
0214 pln.propDoseCalc.doseGrid.resolution.x = 3; % [mm]
0215 pln.propDoseCalc.doseGrid.resolution.y = 3; % [mm]
0216 pln.propDoseCalc.doseGrid.resolution.z = 3; % [mm]
0217 
0218 %% Generate Beam Geometry STF
0219 stf = matRad_generateStf(ct,cst,pln);
0220 
0221 %% Dose Calculation
0222 dij = matRad_calcPhotonDose(ct,stf,pln,cst);
0223 
0224 %% Inverse Optimization for intensity-modulated photon therapy
0225 % The goal of the fluence optimization is to find a set of bixel/spot
0226 % weights which yield the best possible dose distribution according to the
0227 % clinical objectives and constraints underlying the radiation treatment.
0228 resultGUI = matRad_fluenceOptimization(dij,cst,pln);
0229 
0230 %% Plot the resulting dose slice
0231 plane      = 3;
0232 slice      = round(pln.propStf.isoCenter(1,3)./ct.resolution.z);
0233 doseWindow = [0 max([resultGUI.physicalDose(:)])];
0234 
0235 figure,title('phantom plan')
0236 matRad_plotSliceWrapper(gca,ct,cst,1,resultGUI.physicalDose,plane,slice,[],[],colorcube,[],doseWindow,[]);
0237 
0238 %%
0239 % We export the the created phantom & dose as dicom. This is handled by the
0240 % class matRad_DicomExporter. When no arguments are given, the exporter searches
0241 % the workspace itself for matRad-structures. The output directory can be set by
0242 % the property dicomDir. While the different DICOM datasets (ct, RTStruct, etc)
0243 % can be exported individually, we call the wrapper to do all possible exports.
0244 dcmExport = matRad_DicomExporter();
0245 dcmExport.dicomDir = [pwd filesep 'dicomExport'];
0246 dcmExport.matRad_exportDicom();
0247