plotting off-axis dose profile

[naninano1]

Hello,
I have generated my own machine data, following your guidance. Currently, I want to benchmark the validity of my model. Therefore, I need to compare the measured and calculated PDDs and Profiles. For PDD, I used the codes which have been posted in other topics.

1-I was wondering if you could please let me know how I can plot dose profiles at different depths?
2- Is there a way to export the 2D dose distribution in the main GUI to ASCII format? So that, the PDD and Profile extraction might be easier.

Thanks.

[wahln]

1. Depending on the geometry, this may be either quite simple. If you use an orthogonal geometry (e.g. the TG119 or BOXPHANTOM) with 0,90,... gantry angle, you could just extract the profile as "lines" from the dose cube, e.g. plot(resultGUI.physicalDose(83,:,65)) would give you a lateral plot along the x axis (note that there is an annoying permutation of 1st & 2nd index vs coordinate - index 1 -> y, index 2 -> x). Needs some tuning for the resolution and exact depth etc, but should give you some results.
2. I think in Matlab you can save 2D arrays in ascii using save('filename.txt','variableName','-ascii');

[naninano1]

Thank you very much for your help Wahln,
Since the BOXPHANTOM doesn't support field sizes larger than ~20 cm, I use a larger phantom.
I am trying to run the following code to obtain PDDs and profiles of 5, 8, 15, 30 cm filed sizes. However, I receive these errors:

Reference to non-existent field 'ray'.

Error in matRad_computeSSD (line 53)
stf(i).ray(j).targetPoint, ...
Reference to non-existent field 'cube'.
Error in matRad_computeSSD (line 54)
{ct.cube{1}});

Error in test (line 195)
stf_tmp = matRad_computeSSD(stf,ct);

I appreciate your help in advance.

matRad_rc; %If this throws an error, run it from the parent directory first to set the paths

%% Create a CT image series
xDim = 200;
yDim = 200;
zDim = 200;

ct.cubeDim      = [yDim xDim zDim]; % second cube dimension represents the x-coordinate
ct.resolution.x = 3;
ct.resolution.y = 3;
ct.resolution.z = 3;
ct.numOfCtScen  = 1;
 
% create an ct image series with zeros - it will be filled later
ct.cubeHU{1} = ones(ct.cubeDim) * -1000; % assign HU of Air

%% Create the VOI data for the phantom
% Now we define structures a contour for the phantom and a target

ixOAR = 1;
ixPTV = 2;

% define general VOI properties
cst{ixOAR,1} = 0;
cst{ixOAR,2} = 'contour';
cst{ixOAR,3} = 'OAR';

cst{ixPTV,1} = 1;
cst{ixPTV,2} = 'target';
cst{ixPTV,3} = 'TARGET';
 
% define optimization parameter for both VOIs
cst{ixOAR,5}.TissueClass  = 1;
cst{ixOAR,5}.alphaX       = 0.1000;
cst{ixOAR,5}.betaX        = 0.0500;
cst{ixOAR,5}.Priority     = 2;
cst{ixOAR,5}.Visible      = 1;
cst{ixOAR,5}.visibleColor = [0 0 0];

% define objective as struct for compatibility with GNU Octave I/O
cst{ixOAR,6}{1} = struct(DoseObjectives.matRad_SquaredOverdosing(10,30));

cst{ixPTV,5}.TissueClass = 1;
cst{ixPTV,5}.alphaX      = 0.1000;
cst{ixPTV,5}.betaX       = 0.0500;
cst{ixPTV,5}.Priority    = 1;
cst{ixPTV,5}.Visible     = 1;
cst{ixPTV,5}.visibleColor = [1 1 1];

% define objective as struct for compatibility with GNU Octave I/O
cst{ixPTV,6}{1} = struct(DoseObjectives.matRad_SquaredDeviation(800,60));

%% Lets create either a cubic or a spheric phantom

TYPE = 'cubic';   % either 'cubic' or 'spheric'

% first the OAR
cubeHelper = zeros(ct.cubeDim);

switch TYPE
   
   case {'cubic'}
      
      xLowOAR  = 20;
      xHighOAR = 180;
      yLowOAR  = 20;
      yHighOAR = 180;
      zLowOAR  = 20;
      zHighOAR = 180;
      
      for x = xLowOAR:1:xHighOAR
         for y = yLowOAR:1:yHighOAR
            for z = zLowOAR:1:zHighOAR
               cubeHelper(y,x,z) = 1;
            end
         end
      end
      
   
      
end

% extract the voxel indices and save it in the cst
cst{ixOAR,4}{1} = find(cubeHelper);


% second the PTV
cubeHelper = zeros(ct.cubeDim);

switch TYPE
   
   case {'cubic'}
      
      xLowPTV  = round(xDim/2 - xDim/8);
      xHighPTV = round(xDim/2 + xDim/8);
      yLowPTV  = round(yDim/2 - yDim/8);
      yHighPTV = round(yDim/2 + yDim/8);
      zLowPTV  = round(zDim/2 - zDim/8);
      zHighPTV = round(zDim/2 + zDim/8);
      
      cubeHelper = zeros(ct.cubeDim);
      
      for x = xLowPTV:1:xHighPTV
         for y = yLowPTV:1:yHighPTV
            for z = zLowPTV:1:zHighPTV
               cubeHelper(y,x,z) = 1;
            end
         end
      end
      
   
      
end



% extract the voxel indices and save it in the cst
cst{ixPTV,4}{1} = find(cubeHelper);


% now we have ct data and cst data for a new phantom
display(ct);
display(cst);

%% Assign relative electron densities
vIxOAR = cst{ixOAR,4}{1};
vIxPTV = cst{ixPTV,4}{1};

ct.cubeHU{1}(vIxOAR) = 0;
ct.cubeHU{1}(vIxPTV) = 0;

pln.radiationMode   = 'photons';     % either photons / protons / carbon
pln.machine         = 'GENERIC';

pln.numOfFractions  = 30;

% beam geometry settings
pln.propStf.bixelWidth      = 250; % [mm] / also corresponds to lateral spot spacing for particles
pln.propStf.gantryAngles    = 0; % [°]
pln.propStf.couchAngles     = 0; % [°]
pln.propStf.numOfBeams      = numel(pln.propStf.gantryAngles);
pln.propStf.isoCenter       = [300,60,300];

% dose calculation settings
pln.propDoseCalc.doseGrid.resolution.x = ct.resolution.x; % [mm]
pln.propDoseCalc.doseGrid.resolution.y = ct.resolution.y; % [mm]
pln.propDoseCalc.doseGrid.resolution.z = ct.resolution.z; % [mm]

% optimization settings
pln.propOpt.optimizer       = 'IPOPT';
pln.propOpt.bioOptimization = 'none'; % none: physical optimization;             const_RBExD; constant RBE of 1.1;
                                      % LEMIV_effect: effect-based optimization; LEMIV_RBExD: optimization of RBE-weighted dose
pln.propOpt.runDAO          = false;  % 1/true: run DAO, 0/false: don't / will be ignored for particles
pln.propOpt.runSequencing   = false;  % 1/true: run sequencing, 0/false: don't / will be ignored for particles and also triggered by runDAO below

%% define steering file by hand

bixelWidths = [50 80 150 300];

load(['basedata/' pln.radiationMode '_' pln.machine '.mat']);
stf.gantryAngle       = pln.propStf.gantryAngles;
stf.couchAngle        = pln.propStf.couchAngles;
stf.bixelWidth        = pln.propStf.bixelWidth;
stf.radiationMode     = pln.radiationMode;
stf.SAD               = 1000;
stf.isoCenter         = pln.propStf.isoCenter;
stf.numOfRays         = 1;
stf.sourcePoint_bev   = [0 -stf.SAD 0];
stf.sourcePoint       = [0 -stf.SAD 0];
stf.numOfBixelsPerRay = 1;
stf.totalNumOfBixels  = 1;


stf_tmp = matRad_computeSSD(stf,ct);
stf.isoCenter(2) = stf.isoCenter(2) - (stf.SAD - stf_tmp.ray.SSD);
pln.propStf.isoCenter = stf.isoCenter;

hF = figure;
hAxAbs = subplot(1,2,1);
hprofile = subplot(1,2,2);

for width = bixelWidths
    pln.propStf.bixelWidth = width;
    stf.bixelWidth = width;   
    
    resultGUI_tmp = matRad_calcDoseDirect(ct,stf,pln,cst,1);
    plot(hAxAbs,ct.resolution.x*(1:ct.cubeDim(2)),resultGUI_tmp.physicalDose(:,100,100),'DisplayName',sprintf('%dx%d',width,width),'LineWidth',1.5); hold(hAxAbs,'on');
	plot(hprofile,ct.resolution.y*(1:ct.cubeDim(2)),resultGUI_tmp.physicalDose(40,:,100),'DisplayName',sprintf('%dy%d',width,width),'LineWidth',1.5); hold(hprofile,'on');
    resultGUI.(['physicalDose_fSize' num2str(width)]) = resultGUI_tmp.physicalDose;
end
%plot(resultGUI.physicalDose(83,:,65))



%% start gui for visualization of result
matRadGUI

[mingersming]

naninano1 what a coincidence. :)

Glad to hear everything is working now, if you need further help with your own base data feel free to contact us.

[naninano1]

naninano1 what a coincidence. :)

Glad to hear everything is working now, if you need further help with your own base data feel free to contact us.

Yes :)
I do appreciate your kind consideration and time.

[wahln]

Wait a minute @mingersming, I think you brought me on the right track... I don't think the kernelLimit here should be changed. this just gives the extent of the scattering kernel, which later on gets convolved correctly with the full sized fluence (i.e., kernelConvLimit is the relevant quantity and correctly dependens on kernelLimit AND fieldSize). Isn't the problem that we don't care for the quadratic shape of the beamlets when computing effectiveLateralCutoff, because raytracing & coordinate computation assume circles around the central ray when using the cutoff? It doesn't show up for small fields, because a cutoff of 50 results in a large enough circle that encompasses the quadratic beamlets, but as soon as you get large enough you run into problems.

Probably it would just be enough to change

matRad/matRad_calcPhotonDose.m

Line 137 in 6fad1b8

effectiveLateralCutoff = lateralCutoff + fieldWidth/2;

to
effectiveLateralCutoff = lateralCutoff + fieldWidth/sqrt(2);

Then you also do not need to increase the cutoff, @naninano1, which keeps computation times low.
I will prepare a hotfix and put you as a reviewer, ok?

We have to thank you @naninano1, we really rely on the input of curious people like you to find such bugs!

[wahln]

Check out #492

[naninano1]

Dear Niklas,
It is my pleasure to discuss the issues and take advantage of the bilateral conversations.
I compared these two recommendations, this:

There is already something like an automatic adjustment:

matRad/matRad_calcPhotonDose.m

Line 137 in 6fad1b8

effectiveLateralCutoff = lateralCutoff + fieldWidth/2;

But I think the code from line 137 has to be placed before line 122 to get the supposed behavior in the end:

matRad/matRad_calcPhotonDose.m

Line 122 in 6fad1b8

kernelLimit = ceil(lateralCutoff/intConvResolution);

And a change to:

kernelLimit = ceil(effectiveLateralCutoff/intConvResolution);

and this:

Wait a minute @mingersming, I think you brought me on the right track... I don't think the kernelLimit here should be changed. this just gives the extent of the scattering kernel, which later on gets convolved correctly with the full sized fluence (i.e., kernelConvLimit is the relevant quantity and correctly dependens on kernelLimit AND fieldSize). Isn't the problem that we don't care for the quadratic shape of the beamlets when computing effectiveLateralCutoff, because raytracing & coordinate computation assume circles around the central ray when using the cutoff? It doesn't show up for small fields, because a cutoff of 50 results in a large enough circle that encompasses the quadratic beamlets, but as soon as you get large enough you run into problems.

Probably it would just be enough to change

matRad/matRad_calcPhotonDose.m

Line 137 in 6fad1b8

effectiveLateralCutoff = lateralCutoff + fieldWidth/2;

to
effectiveLateralCutoff = lateralCutoff + fieldWidth/sqrt(2);

I am going to put the results on the created topic: #492