Skip to content

Commit

Permalink
Merge pull request e0404#610 from tobiasbecher/dev_varRBErobOptPhanto…
Browse files Browse the repository at this point in the history 
…mBuilder

Phantom builder added
  • Loading branch information
@wahln
wahln authored Jul 20, 2023
2 parents acd73e7 + 2958dc4 commit 87832ed
Show file tree
Hide file tree
Showing 7 changed files with 390 additions and 148 deletions.
2 changes: 1 addition & 1 deletion .gitignore
View file
Original file line number Diff line number Diff line change
@@ -1,2 +1,2 @@
MCsquare/bin/Materials/
topas/MCrun
topas/MCrun
1 change: 1 addition & 0 deletions AUTHORS.txt
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
List of all matRad developers that contributed code (alphabetical)

* Mark Bangert
* Tobias Becher
* Amit Ben Antony Bennan
* Lucas Burigo
* Gonzalo Cabal
Expand Down
187 changes: 40 additions & 147 deletions examples/matRad_example1_phantom.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2018 the matRad development team.
% Copyright 2023 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
Expand All @@ -19,162 +19,53 @@
% (iii) generate a treatment plan for this phantom

%% set matRad runtime configuration
%clear all; %somewhat needed for the phantom builder
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 = 50;

ct.cubeDim = [yDim xDim zDim]; % second cube dimension represents the x-coordinate
ct.resolution.x = 2;
ct.resolution.y = 2;
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;
ctDim = [200,200,100]; % x,y,z dimensions
ctResolution = [2,2,3]; % x,y,z the same here!

%% 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 = 'spheric'; % either 'cubic' or 'spheric'

% first the OAR
cubeHelper = zeros(ct.cubeDim);
%This uses the phantombuilder class, which helps to easily implement a
%water phantom containing geometrical 3D shapes as targets and organs
builder = matRad_PhantomBuilder(ctDim,ctResolution,1);

switch TYPE

case {'cubic'}

xLowOAR = round(xDim/2 - xDim/4);
xHighOAR = round(xDim/2 + xDim/4);
yLowOAR = round(yDim/2 - yDim/4);
yHighOAR = round(yDim/2 + yDim/4);
zLowOAR = round(zDim/2 - zDim/4);
zHighOAR = round(zDim/2 + zDim/4);

for x = xLowOAR:1:xHighOAR
for y = yLowOAR:1:yHighOAR
for z = zLowOAR:1:zHighOAR
cubeHelper(y,x,z) = 1;
end
end
end

case {'spheric'}

radiusOAR = xDim/4;

for x = 1:xDim
for y = 1:yDim
for z = 1:zDim
currPost = [y x z] - round([ct.cubeDim./2]);
if sqrt(sum(currPost.^2)) < radiusOAR
cubeHelper(y,x,z) = 1;
end
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);
%% Create the VOI data for the phantom
%To add a VOI there are (so far) two different options
%either a Box or a spherical Volume (either OAR or Target)
%NOTE: The order in which the objectives are initialized matters!
% In case of overlaps in the objectives, the firstly created objectives have
% a higher priority! This means that if two VOI have an overlap with
% different HU, then the value of the firstly initialized objective will be
% set in the overlap region

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

case {'spheric'}

radiusPTV = xDim/12;

for x = 1:xDim
for y = 1:yDim
for z = 1:zDim
currPost = [x y z] - round([ct.cubeDim./2]);
if sqrt(sum(currPost.^2)) < radiusPTV
cubeHelper(y,x,z) = 1;
end
end
end
end

end

%define objectives for the VOI

objective1 = struct(DoseObjectives.matRad_SquaredDeviation(800,45));
objective2 = struct(DoseObjectives.matRad_SquaredOverdosing(400,0));
objective3 = struct(DoseObjectives.matRad_SquaredOverdosing(10,0));

% extract the voxel indices and save it in the cst
cst{ixPTV,4}{1} = find(cubeHelper);
builder.addSphericalTarget('Volume1',20,'objectives',objective1,'HU',0);
builder.addBoxOAR('Volume2',[60,30,60],'offset',[0 -15 0],'objectives',objective2);
builder.addBoxOAR('Volume3',[60,30,60],'offset',[0 15 0],'objectives',objective3);

% This adds two Box OAR and one Spherical Target in the middle
% For Box Volumes a name (here Volume2 and Volume3) has to be specified,
% as well as the dimension of the Volume.
% For spherical Volumes a name has to be specified, as well as the radius
% of the sphere
%(Optional) To move the VOI from the center of the ct an offset can be set.
%(Optional) The objectives can already be set here, however this can also
%be done later on
%(Optional) The HU of the VOI can be set (normal value: 0 (water))

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

%% Get the ct and cst (stored as properties of the phantomBuilder class)

%% Assign relative electron densities
vIxOAR = cst{ixOAR,4}{1};
vIxPTV = cst{ixPTV,4}{1};
[ct,cst] = builder.getctcst();

ct.cubeHU{1}(vIxOAR) = 0; % assign HU of water
ct.cubeHU{1}(vIxPTV) = 0; % assign HU of water
%% Treatment Plan
% The next step is to define your treatment plan labeled as 'pln'. This
% structure requires input from the treatment planner and defines the most
Expand Down Expand Up @@ -205,8 +96,8 @@
%%
% The remaining plan parameters are set like in the previous example files
pln.numOfFractions = 30;
pln.propStf.gantryAngles = [0 45];
pln.propStf.couchAngles = [0 0];
pln.propStf.gantryAngles = [0:70:355];
pln.propStf.couchAngles = zeros(size(pln.propStf.gantryAngles));
pln.propStf.bixelWidth = 5;
pln.propStf.numOfBeams = numel(pln.propStf.gantryAngles);
pln.propStf.isoCenter = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0);
Expand All @@ -223,6 +114,8 @@
pln.propDoseCalc.doseGrid.resolution.x = 3; % [mm]
pln.propDoseCalc.doseGrid.resolution.y = 3; % [mm]
pln.propDoseCalc.doseGrid.resolution.z = 3; % [mm]
%%
matRadGUI;

%% Generate Beam Geometry STF
stf = matRad_generateStf(ct,cst,pln);
Expand All @@ -231,14 +124,14 @@
dij = matRad_calcPhotonDose(ct,stf,pln,cst);

%% Export dij matrix
matRad_exportDij('dij.bin',dij,stf);

%matRad_exportDij('dij.bin',dij,stf);
%% Inverse Optimization for intensity-modulated photon therapy
% The goal of the fluence optimization is to find a set of bixel/spot
% weights which yield the best possible dose distribution according to the
% clinical objectives and constraints underlying the radiation treatment.
resultGUI = matRad_fluenceOptimization(dij,cst,pln);

%%
matRadGUI
%% Plot the resulting dose slice
plane = 3;
slice = round(pln.propStf.isoCenter(1,3)./ct.resolution.z);
Expand Down
143 changes: 143 additions & 0 deletions phantoms/builder/matRad_PhantomBuilder.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,143 @@
classdef matRad_PhantomBuilder < handle
% matRad_PhantomBuilder
% Class that helps to create radiotherapy phantoms
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2023 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.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

properties (Access = public)
volumes = {};
end

properties (Access = private)
ct;
cst = {};
end

methods (Access = public)
function obj = matRad_PhantomBuilder(ctDim,ctResolution,numOfCtScen)
obj.ct = struct();
obj.ct.cubeDim = [ctDim(2),ctDim(1),ctDim(3)];
obj.ct.resolution.x = ctResolution(1);
obj.ct.resolution.y = ctResolution(2);
obj.ct.resolution.z = ctResolution(3);
obj.ct.numOfCtScen = numOfCtScen;
obj.ct.cubeHU{1} = ones(obj.ct.cubeDim) * -1000;
end

%functions to create Targets %TODO: Option to extend volumes
function addBoxTarget(obj,name,dimensions,varargin)
% Adds a box target
%
% input:
% name: Name of VOI as string
% dimensions: Dimensions of the box as [x,y,z] array
%
% Name-Value pairs:
% 'offset': The offset of the VOI with respect to the
% center of the geometry as [x,y,z] array
% 'objectives': Either a single objective or a cell array
% of objectives
% 'HU': Houndsfield unit of the volume

obj.volumes(end+1) = {matRad_PhantomVOIBox(name,'TARGET',dimensions,varargin{:})};
obj.updatecst();
end

function addSphericalTarget(obj,name,radius,varargin)
% Adds a spherical target
%
% input:
% name: Name of VOI as string
% radius: Radius of the sphere
%
% Name-Value pairs:
% 'offset': The offset of the VOI with respect to the
% center of the geometry as [x,y,z] array
% 'objectives': Either a single objective or a cell array
% of objectives
% 'HU': Houndsfield unit of the volume

obj.volumes(end+1) = {matRad_PhantomVOISphere(name,'TARGET',radius,varargin{:})};
obj.updatecst();
end


function addBoxOAR(obj,name,dimensions,varargin)
% Adds a box OAR
%
% input:
% name: Name of VOI as string
% dimensions: Dimensions of the box as [x,y,z] array
%
% Name-Value pairs:
% 'offset': The offset of the VOI with respect to the
% center of the geometry as [x,y,z] array
% 'objectives': Either a single objective or a cell array
% of objectives
% 'HU': Houndsfield unit of the volume

obj.volumes(end+1) = {matRad_PhantomVOIBox(name,'OAR',dimensions,varargin{:})};
obj.updatecst();
end

function addSphericalOAR(obj,name,radius,varargin)
% Adds a spherical OAR
%
% input:
% name: Name of VOI as string
% radius: Radius of the sphere
%
% Name-Value pairs:
% 'offset': The offset of the VOI with respect to the
% center of the geometry as [x,y,z] array
% 'objectives': Either a single objective or a cell array
% of objectives
% 'HU': Houndsfield unit of the volume

obj.volumes(end+1) ={matRad_PhantomVOISphere(name,'OAR',radius,varargin{:})};
obj.updatecst();
end


function [ct,cst] = getctcst(obj)
% Returns the ct and struct. The function also initializes
% the HUs in reverse order of defintion
%
% output
% ct: matRad ct struct
% cst: matRad cst struct

%initialize the HU in reverse order of definition (objectives
%defined at the start will have the highest priority in case of
%overlaps)

for i = 1:size(obj.cst,1)
vIxVOI = obj.cst{end-i+1,4}{1};
obj.ct.cubeHU{1}(vIxVOI) = 0; % assign HU
end

ct = obj.ct;
cst = obj.cst;
end

end

methods (Access = private)

function updatecst(obj)
obj.cst = obj.volumes{end}.initializeParameters(obj.ct,obj.cst);
end

end
end
Loading

0 comments on commit 87832ed

Please sign in to comment.