matRad_example6_protonsNoise

Purpose ^

% Example: Proton Treatment Plan with Manipulated CT values

Synopsis ^

This is a script file.

Description ^

% 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.

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

Cross-reference information ^

This function calls: This function is called by:

Source code ^

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

| Generated by m2html © 2005