TROTS dataset in matRad
ferdymercury opened this issue · 20 comments
I want to implement and contribute the following feature to matRad:
The TROTS dataset is a very interesting set of CTs and treatment plans for benchmarking optimizers. I am working on importing this dataset into matRad.
Current state of the contribution:
I have already correctly imported the CT, all contours, as well as treatment plan constraints and objectives.
TROTS screenshot:
matRad screenshot:
Script below:
% The code expects patientFolder to contain a TROTS *.mat file
% See e.g. https://zenodo.org/records/2708302/files/Protons.zip?download=1
% To open the result, if you did not close matRad, just open matRadGUI and
% click Refresh. No need to touch nothing else.
% If you already closed MATLAB, open it again, start matRadGUI, and then
% click on the LoadMat button, and select mini.mat, not the original TROTS
% Authored by:
% S. Tattenberg - TRIUMF and NOSM
% F. Hueso-González - IFIC (CSIC/UV)
clear, clc
patientFolder = '/tmp/'; % with TROTS mat file
TrotsMatFile = patientFolder + "Protons_01.mat";
load(TrotsMatFile);
%TROTS has: data, patient, problem, problem_lex, solutionX
%matRad needs: cst, ct, patientFolder, pln, stf, dij
%% Define CT structure.
ct.resolution.x = patient.Resolution(1);
ct.resolution.y = patient.Resolution(2);
ct.resolution.z = patient.Resolution(3);
nRows = size(patient.CT, 2); % DICOM Rows, goes with y, is index 2 because of how matrix is stored
nColumns = size(patient.CT, 1);% DICOM Columns, goes with x, is index 1 because of how matrix is stored
nSlices = size(patient.CT, 3);% DICOM slices, goes with z
ct.x = linspace(patient.Offset(1), patient.Offset(1)+(nColumns-1)*patient.Resolution(1), nColumns);
ct.y = linspace(patient.Offset(2), patient.Offset(2)+(nRows-1)*patient.Resolution(2), nRows);
ct.z = linspace(patient.Offset(3), patient.Offset(3)+(nSlices-1)*patient.Resolution(3), nSlices);
ct.cubeDim = [nRows nColumns nSlices];
ct.numOfCtScen = 1;
ct.timeStamp = string(datetime("now"));
ct.cubeHU = {double(permute(patient.CT,[2 1 3]))};
%% Create cst structure
% see https://github.com/e0404/matRad/blob/master/dicom/matRad_createCst.m
% and https://github.com/e0404/matRad/blob/master/dicom/matRad_dummyCst.m
% and https://github.com/e0404/matRad/wiki/The-cst-cell
nStructures = length(patient.StructureNames);
defaultColors = colorcube(nStructures);
cst = cell(nStructures,6);
[grx,gry] = ndgrid(ct.x,ct.y);
disp('Calculating contours')
for i = 1:nStructures
cst{i,1} = i-1;
cst{i,2} = patient.StructureNames{i};
if contains(cst{i,2}, 'ctv', 'IgnoreCase', true) ...
|| contains(cst{i,2}, 'ptv', 'IgnoreCase', true)
cst{i,3} = 'TARGET';
else
cst{i,3} = 'OAR';
end
linvoxs = [];
disp(cst{i,2})
for s = 1:nSlices
csXY = patient.Contours{1,1}{s,i};
for sc = 1:length(csXY)
contoursXY = csXY{sc};
in = inpolygon(grx,gry,contoursXY(:,1),contoursXY(:,2)).';
ind = find(in) + (s-1)*nRows*nColumns;
linvoxs = cat(1,linvoxs, ind);
end
if ~isempty(csXY)
disp(s)
end
end
cst{i,4}{1} = linvoxs;
if strcmp(cst{i,3}, 'OAR')
cst{i,5}.Priority = 2;
else
cst{i,5}.Priority = 1;
end
cst{i,5}.alphaX = 0.1;
cst{i,5}.betaX = 0.05;
cst{i,5}.Visible = 1;
cst{i,5}.visibleColor = defaultColors(i,:);
cst{i,6} = [];
end
%% Define constraints and objectives in cst
for OARIndex=1:size(cst,1)
totalIndices = size(problem, 2);
objectiveIndex = 1;
for index=1:totalIndices
if contains(cst(OARIndex,2),problem(index).Name)
cst{OARIndex,6}{objectiveIndex} = struct();
if problem(index).IsConstraint == 1
cst{OARIndex,6}{objectiveIndex}.className = 'DoseConstraints.matRad_MinMaxDose';
if problem(index).Minimise == 1
cst{OARIndex,6}{objectiveIndex}.parameters = cell({0,problem(index).Objective,1}); % 1 is approx, 2 is voxel-wise
else
cst{OARIndex,6}{objectiveIndex}.parameters = cell({problem(index).Objective,100,1});
end
cst{OARIndex,6}{objectiveIndex}.epsilon = 1.0000e-03;
else
if problem(index).Minimise == 1
cst{OARIndex,6}{objectiveIndex}.className = 'DoseObjectives.matRad_SquaredOverdosing';
else
cst{OARIndex,6}{objectiveIndex}.className = 'DoseObjectives.matRad_SquaredUnderdosing';
end
cst{OARIndex,6}{objectiveIndex}.parameters = cell({problem(index).Objective});
cst{OARIndex,6}{objectiveIndex}.penalty = problem(index).Weight;
end
objectiveIndex = objectiveIndex + 1;
end
end
end
%% clear and save
clearvars -except cst ct patientFolder
save([patientFolder 'mini.mat'], '-v7.3')
%% Definition of pln
pln = struct;
pln.propStf = struct;
pln.propStf.bixelWidth = 5;
numBeams = size(patient.Beams.BeamConfig,2);
pln.propStf.gantryAngles = [];
for beamIndex=1:numBeams
pln.propStf.gantryAngles(beamIndex) = patient.Beams.BeamConfig(beamIndex).Gantry;
end
numCouches= size(patient.Beams.BeamConfig,2);
pln.propStf.couchAngles = [];
for couchIndex=1:numCouches
pln.propStf.couchAngles(couchIndex) = patient.Beams.BeamConfig(couchIndex).Couch;
end
pln.propStf.numOfBeams = numBeams;
for i=1:numBeams
img_isox = patient.Isocentre(1) - ct.x(1);% matrad subtracts offset
img_isoy = patient.Isocentre(2) - ct.y(1);% matrad subtracts offset
img_isoz = patient.Isocentre(3) - ct.z(1);% matrad subtracts offset
pln.propStf.isoCenter(i,:) = [img_isox img_isoy img_isoz]; % not sure what this is
end
pln.voxelDimensions = ct.cubeDim;
pln.numOfVoxels = prod(pln.voxelDimensions);
pln.numOfFractions = 1;
if contains(patientFolder,'Proton', 'IgnoreCase', true)
pln.radiationMode = 'protons';
pln.propOpt.bioOptimization = 'const_RBExD';
else
pln.radiationMode = 'photons';
pln.propOpt.bioOptimization = 'none';
end
pln.machine = 'Generic';
pln.propOpt.runDAO = 0;
pln.propOpt.runSequencing = 0;
stf = matRad_generateStf(ct,cst,pln);
%% Pass patient TROTS Dij to MATLAB
patientStructIdx = find(strcmpi(patient.StructureNames,'Patient'));
sampVox = cell2mat(patient.SampledVoxels(patientStructIdx));
svoxels = size(sampVox, 2);
doseNames = struct2cell(data.matrix);
patientDoseIdx = find(strcmpi(doseNames(1,:),'Patient'));
pat_dij = data.matrix(patientDoseIdx).A;
voxels = size(pat_dij,1);
spots = size(pat_dij,2);
if data.matrix(patientDoseIdx).b ~= 0
fprintf('Plan recalculation not yet supported')
return
end
if ~isempty(data.matrix(patientDoseIdx).c)
fprintf('Quadratic cost functions not yet supported')
return
end
if mod(voxels,9) == 0 && voxels/9 == svoxels
fprintf('Multiple scenarios not yet supported')
return
end
if voxels ~= svoxels
fprintf('Wrong number of bixels')
return
end
dij.doseGrid.resolution.x = ct.resolution.x;
dij.doseGrid.resolution.y = ct.resolution.y;
dij.doseGrid.resolution.z = ct.resolution.z;
dij.doseGrid.x = ct.x;
dij.doseGrid.y = ct.y;
dij.doseGrid.z = ct.z;
dij.doseGrid.dimensions = ct.cubeDim;
dij.doseGrid.numOfVoxels = prod(dij.doseGrid.dimensions);
dij.ctGrid.resolution.x = ct.resolution.x;
dij.ctGrid.resolution.y = ct.resolution.y;
dij.ctGrid.resolution.z = ct.resolution.z;
dij.ctGrid.x = ct.x;
dij.ctGrid.y = ct.y;
dij.ctGrid.z = ct.z;
dij.ctGrid.dimensions = ct.cubeDim;
dij.ctGrid.numOfVoxels = prod(dij.ctGrid.dimensions);
% % adjust isocenter internally for different dose grid
% offset = [dij.doseGrid.resolution.x - dij.ctGrid.resolution.x ...
% dij.doseGrid.resolution.y - dij.ctGrid.resolution.y ...
% dij.doseGrid.resolution.z - dij.ctGrid.resolution.z];
%
% for i = 1:numel(stf)
% stf(i).isoCenter = stf(i).isoCenter + offset;
% end
% %set up HU to rED or rSP conversion
% if ~isfield(pln,'propDoseCalc') || ~isfield(pln.propDoseCalc,'useGivenEqDensityCube')
% disableHUconversion = matRad_cfg.propDoseCalc.defaultUseGivenEqDensityCube;
% else
% disableHUconversion = pln.propDoseCalc.useGivenEqDensityCube;
% end
%
% %If we want to omit HU conversion check if we have a ct.cube ready
% if disableHUconversion && ~isfield(ct,'cube')
% matRad_cfg.dispWarning('HU Conversion requested to be omitted but no ct.cube exists! Will override and do the conversion anyway!');
% disableHUconversion = false;
% end
%
% % calculate rED or rSP from HU
% if disableHUconversion
% matRad_cfg.dispInfo('Omitting HU to rED/rSP conversion and using existing ct.cube!\n');
% else
% ct = matRad_calcWaterEqD(ct, pln);
% end
%
% meta information for dij
dij.numOfBeams = pln.propStf.numOfBeams;
dij.numOfScenarios = 1; % for the moment we exclude the 9 scenarios TROTS
dij.numOfRaysPerBeam = patient.Beams.ElementIndex;
dij.totalNumOfBixels = voxels*spots;% sum([stf(:).totalNumOfBixels]);
dij.totalNumOfRays = sum(dij.numOfRaysPerBeam);
if contains(patientFolder,'Proton', 'IgnoreCase', true)
dij.RBE = 1.1;
endI have the following issues:
I am now stuck on how to pass the Dij to matRad, so any help in this regard is welcome.
pat_dij is 24760x1080 sparse double, where 1080 is the sum of all spots (rays).
sampVox is 3x24760 uint16, which tells you the x,y,z index of each voxel of the dij. I can easily linearize those indices into a 1x24760 array.
I think @tobiasbecher is your best person to talk to here. He started with running the TROTS dataset in matRad as well. Maybe you can work together so you don't double your workload!
Thanks for the info!! I'd be glad to join efforts with @tobiasbecher, please let me know how far you went into the TROTS conversion to matRad
Hey @ferdymercury
maybe first on: I have not looked too much into an actual conversion of TROTS data to matRad, but I think I found a way this should work.
But to clarify one or two things first:
dij.numOfRaysPerBeam = patient.Beams.ElementIndex;
dij.totalNumOfBixels = voxels*spots;% sum([stf(:).totalNumOfBixels]);
dij.totalNumOfRays = sum(dij.numOfRaysPerBeam);
This does not fully makes sense with respect to the ray bixel concept in matRad: https://github.com/e0404/matRad/wiki/Dose-influence-matrix-calculation
As far as I can tell there is no information provided by TROTS on the rays - but you do not really need them for the calculation (so just leave them blank). And totalNumOfBixels is just be the number of spots (1080 for the example).
Regarding the actual dij: Just to make sure that we are on the same page: You want to get the dij matrix to match
voxel * bixel? (so (ct_x * ct_y * ct_z) x bixel) (the latter should be 1080 here - the size of solutionX)
As far as I can tell, the main issue there is the interpolation from sampled voxels to full voxels.
A dirty way to solve this would be to use https://github.com/SebastiaanBreedveld/TROTS/blob/master/Scripts%20Matlab/TROTSComputeDose.m
and just run this with solutionX as a vector of just 0's and a single 1. Depending on the position of the 1 it should return a different column of the Dij matrix, which can be then stitched together.
Does this make sense to you?
There are probably some nicer ways to do this, but I guess since the dij has to be only converted once it would work.
(Other methods should work similarly to the griddata approach in the TROTSComputeDose method)
I have to see how much more time I can use for the problem this week, but can look into it again next week :)
Two more comments:
- I have tried to run an optimization with a matRad generated dij based on your ct and ct which works nicely
- I would however recommend to update the minMaxDose approximation. There is an issue with the epsilon approximation on the current version which is fixed here: https://github.com/tobiasbecher/matRad/blob/dev_ParetoNav/optimization/%2BDoseConstraints/matRad_MinMaxDose.m
Thanks a lot Tobias!!
As far as I can tell there is no information provided by TROTS on the rays - but you do not really need them for the calculation (so just leave them blank).
patient.Beams.ElementIndex tells you how many spots per beam angle you have.
And totalNumOfBixels is just be the number of spots (1080 for the example).
ohh, I see, so we are imposing here that 1 ray == 1 spot == 1 bixel ?
You want to get the dij matrix to match
voxel * bixel? (so (ct_x * ct_y * ct_z) x bixel) (the latter should be 1080 here - the size of solutionX)
What I really want is to use the Dij from TROTS, that includes all the information on the machine beam parameters and doses, rather than using a default generic machine with their own dose calculation algorithm.
So, I thought the best way would be to use the dose matrix from TROTS, which is sampled in the CT grid, and convert it to matRad format. Ideally keeping the same 'sparseness'. But I am open to other ideas :)
Depending on the position of the 1 it should return a different column of the Dij matrix, which can be then stitched together.
Does this make sense to you?
Ok, yes, I understand, thanks for the idea.
So, to sum up:
dij.bixelNum wil go from 1 to 1080
dij.beamNum will go from 1 to 3
dij.rayNum will go e.g. from 1 to 200, then from 1 to 300, then from 1 to 250 (exact numbers are in patient.ElementIndex)
and then dij.physicalDose{1} = spalloc(numOfVoxels, 1080), right?
I will just enter this conversation and say stupid things that confuse everybody due to lack of in-depth knowledge about the TROTS format.
dij parameters for bixel, ray, and beam number
Check this out: https://github.com/e0404/matRad/wiki/The-dij-struct
dij.bixelNum, dij.beamNum and dij.rayNum all have number of beamlets entries. They are not important for optimization, but to identify which part of the dij belongs to what when accessing the stf.
dij.beamNumwill, for each bixel, tell to which beam it belongs. So if your first beam has 100 beamlets/bixels, the first 100 entries will be 1 and the beamlets will be stored in the first 100 columns of the dij.physicalDose{1}.dij.rayNumwill tell which ray in the beam each bixel belongs to. So for above beam, the first 100 entries will run from 1 to 100, and start at 1 again at the next beam.dij.bixelNumwill tell you which bixel on the ray of the beam it is. For photons, one ray has one bixel (for particles, a ray can have multiple bixels / spots with different energies). So it is just a vector of ones.
How does this fit together? For an arbitrary bixel j You can find, for example, a specific bixels energy level in a ray of the stf by:
bixelEnergy = stf(dij.beamNum(j)).ray(dij.rayNum(j)).energy(dij.bixelNum(j));
Resolution
Since TROTS is apparantly running on ctGrid, make sure dij.ctGrid and dij.doseGrid contain correct/the same resolutions such that matRad does not want to resample stuff
importing the sparse dij
As far as I remember, TROTS stores the dijs per structure, so you need to piece them together. One particular problem you will have is that there are also these smoothing matrices of bixel x bixel dimension for quadratic form objectives (i.e. w'*A*w), which we don't have an official implementation in matRad for. For now, I would advise to create a new field in the dij, like quadForms, and store them in there with some index to which structure they belong. We can figure out a way to implement such an objective nicely together.
- So it is just a vector of ones
Ahh ok got it. TROTS does not say which energies or positions correspond to each spot. It may be that there are spots on same xy and different energy, but we don't know, so we just put all in separate bixels (all 1s as you said).
TROTS stores the dijs per structure, so you need to piece them together
Yes, but there is also one Dij for the whole patient contour, so that one should be enough?
Thanks for the help! I'm a newbie on matRad so your comments are warmly appreciated!
Regarding the bixel/rays: You can adopt this information if you want, but it is not strictly necessary for the matRad optimization. Once you have the dij matrix the only thing you need to do is call fluenceOptimization and while this does require the dij, it never uses the ray information in there (I dont know where this is used exactly in matRad, but maybe Niklas can tell)
Yes, but there is also one Dij for the whole patient contour, so that one should be enough?
So the thing with their dij is, that it only uses a fraction of voxels in each structure and in the end uses some grid interpolation to get the final dose
% WARNING!!! CAUTION!!! WARNING!!! CAUTION!!! WARNING!!! CAUTION!!!
%
% Since this dose is computed based on undersampled data
% (especially outside the delineated structures), this dose
% is ONLY an approximation to get an idea of what the dose
% looks like.
%
% There are deviations in target coverage (which is generally
% lower in this reconstruction than in reality).
%
% Nevertheless, the correspondence with the real dose is
% still pretty awesome, just not at the volumes with
% high dose gradients.
%
% WARNING!!! CAUTION!!! WARNING!!! CAUTION!!! WARNING!!! CAUTION!!!
So lets say your structure consists of voxels [1,2,3,4,5]. TROTS now uses e.g. only voxels [1,3,5] for the actual "problem" and then finds the dose for voxels 2 and 4 through grid interpolation.
If you want to use this you also need to adapt the cst, meaning that you only store the sample voxels (here[1,3,5]) in there. Right now your cst stores all the voxels for the structure.
So either you change the cst to only store the sampled voxels (which would however no longer allow to show the visualization of the volume) or you have to scale the dij to full resolution.
And regarding the patient structure: It should be possible to only use the patient, but this would reduce the precision of the interpolation.
E.g lets say you have voxels [1,2,3,4,5,6,7,8,9,10]. The body samples voxels [1,3,8,10] and Organ1 samples [5,6]. You could use only the body for interpolation to find the remaining information, but it should be more accurate when also including the information from Organ1.
I checked the difference of sampled voxels with
f = setdiff(patient.SampledVoxels{1}',patient.SampledVoxels{2}',"rows") ) and out of the 24760 sampled voxels in the patient structure only 22 are also sampled for the spinal cord
dij.beamNum is needed for the calculation of the beam-wise cubes in the end.
If there is undersampling, I have twosuggestions.
- Store two versions, one with all voxels and one with only undersampled voxels.
- Do nothing and treat the dose in the other voxels which are not sampled as zero? Or am I missing something here?
Thanks for the help!!
This is what I have advanced so far. What do you think?
dij.bixelNum = ones(spots,1);
dij.rayNum = [];
for r=1:length(dij.numOfRaysPerBeam)
dij.rayNum = horzcat(dij.rayNum, 1:1:dij.numOfRaysPerBeam(r));
end
dij.beamNum = [];
for b=1:dij.numOfBeams
dij.beamNum = horzcat(dij.beamNum, ones(1, patient.Beams.ElementIndex(b))*b);
end
% Allocate space for dij.physicalDose sparse matrix
for i = 1:dij.numOfScenarios
dij.physicalDose{i} = spalloc(dij.doseGrid.numOfVoxels,spots,1);
end
sInd = sub2ind(ct.cubeDim, sampVox(2,:), sampVox(1,:), sampVox(3,:));
dij.physicalDose{1}(sInd,:) = pat_dij;
and then:
%% Required to avoid issues with overlap
for OARIndex=1:size(cst,1)
cst{OARIndex,5}.Priority = 1;
end
optimized_data = matRad_fluenceOptimization(dij, cst, pln);
% I had to comment line 226 in matRad_fluenceOptimization.m
If you think this is correct, then I will loop over the high-res Dij of all the other TROTS structures, and start overwriting the physicalDose matrix with those more-oversampled ones.
As a last step, I could also take a look into the interpolation Tobias mentioned.
Hi Fernando,
dij.rayNum & dij.beamNum need to be transposed. Otherwise, the beam weights vector becomes quadratic matrix when using logical indexing with the bixels.
for r=1:length(dij.numOfRaysPerBeam)
%dij.rayNum = horzcat(dij.rayNum, 1:1:dij.numOfRaysPerBeam(r));
dij.rayNum = [dij.rayNum; (1:1:dij.numOfRaysPerBeam(r))'];
end
dij.beamNum = [];
for b=1:dij.numOfBeams
%dij.beamNum = horzcat(dij.beamNum; ones(patient.Beams.ElementIndex(b))*b, 1);
dij.beamNum = [dij.beamNum; ones(patient.Beams.ElementIndex(b),1)*b];
end
then enable function 'matRad_calcCubes' in line 226.
Thanks for the hints!
Below a new version using all Dijs, not just the patient one. Optimization results still look weird. So maybe it's the interpolation issue. Right now, we have: nnz(dij.physicalDose{1})/pln.numOfVoxels/spots = 0.013%, or 13% if we do not divide by nSpots.
Script below:
% The code expects patientFolder to contain a TROTS *.mat file
% See e.g. https://zenodo.org/records/2708302/files/Protons.zip?download=1
% To open the result, if you did not close matRad, just open matRadGUI and
% click Refresh. No need to touch nothing else.
% If you already closed MATLAB, open it again, start matRadGUI, and then
% click on the LoadMat button, and select mini.mat, not the original TROTS
% Authored by:
% S. Tattenberg - TRIUMF and NOSM
% F. Hueso-González - IFIC (CSIC/UV)
% Assistance:
% N. Wahl, T. Becher, A. Hammi
% See https://github.com/e0404/matRad/issues/695
clear, clc
patientFolder = '/tmp/'
TrotsMatFile = patientFolder + "Protons_01.mat";
load(TrotsMatFile);
%TROTS has: data, patient, problem, problem_lex, solutionX
%matRad needs: cst, ct, patientFolder, pln, stf, dij
%% Define CT structure.
ct.resolution.x = patient.Resolution(1);
ct.resolution.y = patient.Resolution(2);
ct.resolution.z = patient.Resolution(3);
nRows = size(patient.CT, 2); % DICOM Rows, goes with y, is index 2 because of how matrix is stored
nColumns = size(patient.CT, 1);% DICOM Columns, goes with x, is index 1 because of how matrix is stored
nSlices = size(patient.CT, 3);% DICOM slices, goes with z
ct.x = linspace(patient.Offset(1), patient.Offset(1)+(nColumns-1)*patient.Resolution(1), nColumns);
ct.y = linspace(patient.Offset(2), patient.Offset(2)+(nRows-1)*patient.Resolution(2), nRows);
ct.z = linspace(patient.Offset(3), patient.Offset(3)+(nSlices-1)*patient.Resolution(3), nSlices);
ct.cubeDim = [nRows nColumns nSlices];
ct.numOfCtScen = 1;
ct.timeStamp = string(datetime("now"));
ct.cubeHU = {double(permute(patient.CT,[2 1 3]))};
%% Create cst structure
% see https://github.com/e0404/matRad/blob/master/dicom/matRad_createCst.m
% and https://github.com/e0404/matRad/blob/master/dicom/matRad_dummyCst.m
% and https://github.com/e0404/matRad/wiki/The-cst-cell
nStructures = length(patient.StructureNames);
defaultColors = colorcube(nStructures);
cst = cell(nStructures,6);
[grx,gry] = ndgrid(ct.x,ct.y);
disp('Calculating contours')
for i = 1:nStructures
cst{i,1} = i-1;
cst{i,2} = patient.StructureNames{i};
if contains(cst{i,2}, 'ctv', 'IgnoreCase', true) ...
|| contains(cst{i,2}, 'ptv', 'IgnoreCase', true)
cst{i,3} = 'TARGET';
else
cst{i,3} = 'OAR';
end
linvoxs = [];
disp(cst{i,2})
for s = 1:nSlices
csXY = patient.Contours{1,1}{s,i};
for sc = 1:length(csXY)
contoursXY = csXY{sc};
in = inpolygon(grx,gry,contoursXY(:,1),contoursXY(:,2)).';
ind = find(in) + (s-1)*nRows*nColumns;
linvoxs = cat(1, linvoxs, ind);
end
if ~isempty(csXY)
disp(s)
end
end
cst{i,4}{1} = linvoxs;
if strcmp(cst{i,3}, 'OAR')
cst{i,5}.Priority = 2;
else
cst{i,5}.Priority = 1;
end
cst{i,5}.alphaX = 0.1;
cst{i,5}.betaX = 0.05;
cst{i,5}.Visible = 1;
cst{i,5}.visibleColor = defaultColors(i,:);
cst{i,6} = [];
end
%% Define constraints and objectives in cst
for OARIndex=1:size(cst,1)
totalIndices = size(problem, 2);
objectiveIndex = 1;
for index=1:totalIndices
if contains(cst(OARIndex,2),problem(index).Name)
cst{OARIndex,6}{objectiveIndex} = struct();
if problem(index).IsConstraint == 1
cst{OARIndex,6}{objectiveIndex}.className = 'DoseConstraints.matRad_MinMaxDose';
if problem(index).Minimise == 1
cst{OARIndex,6}{objectiveIndex}.parameters = cell({0,problem(index).Objective,1}); % 1 is approx, 2 is voxel-wise
else
cst{OARIndex,6}{objectiveIndex}.parameters = cell({problem(index).Objective,100,1});
end
cst{OARIndex,6}{objectiveIndex}.epsilon = 1.0000e-03;
else
if problem(index).Minimise == 1
cst{OARIndex,6}{objectiveIndex}.className = 'DoseObjectives.matRad_SquaredOverdosing';
else
cst{OARIndex,6}{objectiveIndex}.className = 'DoseObjectives.matRad_SquaredUnderdosing';
end
cst{OARIndex,6}{objectiveIndex}.parameters = cell({problem(index).Objective});
cst{OARIndex,6}{objectiveIndex}.penalty = problem(index).Weight;
end
objectiveIndex = objectiveIndex + 1;
end
end
end
%% clear and save
%clearvars -except cst ct patientFolder
%save([patientFolder 'mini.mat'], '-v7.3')
%% Definition of pln
pln = struct;
pln.propStf = struct;
pln.propStf.bixelWidth = ct.resolution.x;
numBeams = size(patient.Beams.BeamConfig,2);
pln.propStf.gantryAngles = [];
for beamIndex=1:numBeams
pln.propStf.gantryAngles(beamIndex) = patient.Beams.BeamConfig(beamIndex).Gantry;
end
numCouches= size(patient.Beams.BeamConfig,2);
pln.propStf.couchAngles = [];
for couchIndex=1:numCouches
pln.propStf.couchAngles(couchIndex) = patient.Beams.BeamConfig(couchIndex).Couch;
end
pln.propStf.numOfBeams = numBeams;
for i=1:numBeams
img_isox = patient.Isocentre(1) - ct.x(1);% matrad subtracts offset
img_isoy = patient.Isocentre(2) - ct.y(1);% matrad subtracts offset
img_isoz = patient.Isocentre(3) - ct.z(1);% matrad subtracts offset
pln.propStf.isoCenter(i,:) = [img_isox img_isoy img_isoz]; % not sure what this is
end
pln.voxelDimensions = ct.cubeDim;
pln.numOfVoxels = prod(pln.voxelDimensions);
pln.numOfFractions = 1;
if contains(patientFolder,'Proton', 'IgnoreCase', true)
pln.radiationMode = 'protons';
pln.propOpt.bioOptimization = 'const_RBExD';
else
pln.radiationMode = 'photons';
pln.propOpt.bioOptimization = 'none';
end
pln.machine = 'Generic';
pln.propOpt.runDAO = 0;
pln.propOpt.runSequencing = 0;
%stf = matRad_generateStf(ct,cst,pln);
%% Pass patient TROTS Dij to MATLAB
spots = size(solutionX, 1);
dij.doseGrid.resolution.x = ct.resolution.x;
dij.doseGrid.resolution.y = ct.resolution.y;
dij.doseGrid.resolution.z = ct.resolution.z;
dij.doseGrid.x = ct.x;
dij.doseGrid.y = ct.y;
dij.doseGrid.z = ct.z;
dij.doseGrid.dimensions = ct.cubeDim;
dij.doseGrid.numOfVoxels = prod(dij.doseGrid.dimensions);
dij.ctGrid.resolution.x = ct.resolution.x;
dij.ctGrid.resolution.y = ct.resolution.y;
dij.ctGrid.resolution.z = ct.resolution.z;
dij.ctGrid.x = ct.x;
dij.ctGrid.y = ct.y;
dij.ctGrid.z = ct.z;
dij.ctGrid.dimensions = ct.cubeDim;
dij.ctGrid.numOfVoxels = prod(dij.ctGrid.dimensions);
% % adjust isocenter internally for different dose grid
% offset = [dij.doseGrid.resolution.x - dij.ctGrid.resolution.x ...
% dij.doseGrid.resolution.y - dij.ctGrid.resolution.y ...
% dij.doseGrid.resolution.z - dij.ctGrid.resolution.z];
%
% for i = 1:numel(stf)
% stf(i).isoCenter = stf(i).isoCenter + offset;
% end
% %set up HU to rED or rSP conversion
% if ~isfield(pln,'propDoseCalc') || ~isfield(pln.propDoseCalc,'useGivenEqDensityCube')
% disableHUconversion = matRad_cfg.propDoseCalc.defaultUseGivenEqDensityCube;
% else
% disableHUconversion = pln.propDoseCalc.useGivenEqDensityCube;
% end
%
% %If we want to omit HU conversion check if we have a ct.cube ready
% if disableHUconversion && ~isfield(ct,'cube')
% matRad_cfg.dispWarning('HU Conversion requested to be omitted but no ct.cube exists! Will override and do the conversion anyway!');
% disableHUconversion = false;
% end
%
% % calculate rED or rSP from HU
% if disableHUconversion
% matRad_cfg.dispInfo('Omitting HU to rED/rSP conversion and using existing ct.cube!\n');
% else
% ct = matRad_calcWaterEqD(ct, pln);
% end
%
% meta information for dij
dij.numOfBeams = pln.propStf.numOfBeams;
dij.numOfScenarios = 1; % for the moment we exclude the 9 scenarios TROTS
dij.numOfRaysPerBeam = patient.Beams.ElementIndex;
dij.totalNumOfBixels = spots;% sum([stf(:).totalNumOfBixels]);
dij.totalNumOfRays = sum(dij.numOfRaysPerBeam);
if contains(patientFolder,'Proton', 'IgnoreCase', true)
dij.RBE = 1.1;
end
dij.bixelNum = ones(spots,1);
dij.rayNum = [];
for r=1:length(dij.numOfRaysPerBeam)
dij.rayNum = [dij.rayNum; (1:1:dij.numOfRaysPerBeam(r))'];
end
dij.beamNum = [];
for b=1:dij.numOfBeams
dij.beamNum = [dij.beamNum; ones(patient.Beams.ElementIndex(b),1)*b];
end
% Allocate space for dij.physicalDose sparse matrix
doseNames = struct2cell(data.matrix);
for i = 1:dij.numOfScenarios
dij.physicalDose{i} = spalloc(dij.doseGrid.numOfVoxels,spots,1);
for structIdx = 1:nStructures
disp(patient.StructureNames(structIdx));
sampVox = cell2mat(patient.SampledVoxels(structIdx));
svoxels = size(sampVox, 2);
doseIdx = find(strcmpi(doseNames(1,:),patient.StructureNames(structIdx)), 1);
struct_dij = data.matrix(doseIdx).A;
voxels = size(struct_dij,1);
if data.matrix(doseIdx).b ~= 0
disp('Plan recalculation not yet supported')
return
end
if ~isempty(data.matrix(doseIdx).c)
disp('Quadratic cost functions not yet supported')
return
end
if mod(voxels,9) == 0 && voxels/9 == svoxels
disp('Warning: Multiple scenarios not yet supported')
struct_dij = struct_dij(1:9:end,:);
voxels = size(struct_dij,1);
end
if voxels ~= svoxels
disp('Wrong number of voxels')
return
end
sInd = sub2ind(ct.cubeDim, sampVox(2,:), sampVox(1,:), sampVox(3,:));
dij.physicalDose{1}(sInd,:) = struct_dij;
end
end
%
% % Allocate memory for dose_temp cell array
% doseTmpContainer = cell(numOfBixelsContainer,dij.numOfScenarios);
%
% % take only voxels inside patient
% VctGrid = [cst{:,4}];
% VctGrid = unique(vertcat(VctGrid{:}));
%
% % ignore densities outside of contours
% if ~isfield(pln,'propDoseCalc') || ~isfield(pln.propDoseCalc,'ignoreOutsideDensities')
% ignoreOutsideDensities = matRad_cfg.propDoseCalc.defaultIgnoreOutsideDensities;
% else
% ignoreOutsideDensities = pln.propDoseCalc.ignoreOutsideDensities;
% end
%
% if ignoreOutsideDensities
% eraseCtDensMask = ones(prod(ct.cubeDim),1);
% eraseCtDensMask(VctGrid) = 0;
% for i = 1:ct.numOfCtScen
% ct.cube{i}(eraseCtDensMask == 1) = 0;
% end
% end
%
% % Convert CT subscripts to linear indices.
% [yCoordsV_vox, xCoordsV_vox, zCoordsV_vox] = ind2sub(ct.cubeDim,VctGrid);
%
% % receive linear indices and grid locations from the dose grid
% tmpCube = zeros(ct.cubeDim);
% tmpCube(VctGrid) = 1;
% % interpolate cube
% VdoseGrid = find(matRad_interp3(dij.ctGrid.x, dij.ctGrid.y, dij.ctGrid.z,tmpCube, ...
% dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z,'nearest'));
%
% % Convert CT subscripts to coarse linear indices.
% [yCoordsV_voxDoseGrid, xCoordsV_voxDoseGrid, zCoordsV_voxDoseGrid] = ind2sub(dij.doseGrid.dimensions,VdoseGrid);
%
% % load base data% load machine file
% fileName = [pln.radiationMode '_' pln.machine];
% try
% load([fileparts(mfilename('fullpath')) filesep 'basedata' filesep fileName]);
% catch
% matRad_cfg.dispError('Could not find the following machine file: %s\n',fileName);
% end
%
% % compute SSDs
% stf = matRad_computeSSD(stf,ct);
%% Required to avoid issues with overlap
for OARIndex=1:size(cst,1)
cst{OARIndex,5}.Priority = 1;
end
optimized_data = matRad_fluenceOptimization(dij, cst, pln);
%%
% ct.dicomInfo.PixelSpacing = [ct.resolution.x ct.resolution.y];
% ct.dicomInfo.SlicePositions = ct.z;
% ct.dicomInfo.SliceThickness = repmat(patient.Resolution(3),1,nSlices);
% ct.dicomInfo.ImagePositionPatient = [patient.Offset(2) patient.Offset(1) patient.Offset(3)];
% if patient.PatientPosition ~= "HFS"
% return
% end
% ct.dicomInfo.PatientPosition = patient.PatientPosition;
% ct.dicomInfo.ImagePositionPatient = [ct.x(1) ct.y(1) ct.z(1)]; % this works only as HFS
% ct.dicomInfo.ImageOrientationPatient = [1;0;0;0;1;0]; % this works only as HFS
% ct.dicomInfo.Width = nColumns;
% ct.dicomInfo.Height = nRows;
% ct.dicomInfo.RescaleSlope = 1;
% ct.dicomInfo.RescaleIntercept = 0;
% ct.dicomInfo.Manufacturer = 'TROTS';
% ct.dicomInfo.PatientName = patient.Identifier;
% ct.dicomMeta = {}; % not neededFirst on sorry that the response took this long.
I think the idea of the approach is nice and if you just calculate the dose with solutionX and load it in the matRadGUI the voxel structure looks reasonable - so the dij should be correct.
However there are two issues right now, that show up when using the full cst:
- Objectives are often normalized with respect to the number of voxels in the structure. This causes issues with the penalties, since they are set based on the undersampled number of voxels.
However, this should not lead to issues with an actual plan calculation. The bigger problem is
2. Because a lot of voxels are always 0, minDose constraints will never work - and they are used for target structures.
So the way to solve this is to store only the sampled voxels in cst{i,4}. Unfortunately you cannot store both, full voxel and sampled voxel information in the same cst - you either have to overwrite your initial cst or create a new one (e.g cst2) that stores only the sampled voxels and is then passed to fluenceOptimization.
Maybe not the nicest option if you want full support, but it should (hopefully) give you a correct optimization.
Might be an interesting option though, to add sampledVoxel support to the cst, but that would be a big change and I dont know if @wahln wants to do this at the moment.
I think with this undersampling the bigger issue is that there is a whole different data structure idea, i.e., that dij's are stored per structure. Getting this undersampling nicely into a full-cube-dij as we have it would require, probably, a lot of inefficient indexing on the matrix, so I think the better approach would be to resample these parts to full resolution.
Are there really also targets with undersampled voxels? These would be the only ones where minDose dose make sense.
@ferdymercury Do you mainly want to use the dataset for work of your own, or do you want to implement it sustainably for further usage?
@ferdymercury Do you mainly want to use the dataset for work of your own, or do you want to implement it sustainably for further usage?
Main motivation is own work, but I'd be happy if it serves the community in improving their optimizers by benchmarking against a standard dataset.
Are there really also targets with undersampled voxels? These would be the only ones where minDose dose make sense.
To some extent, see:
So the way to solve this is to store only the sampled voxels in cst{i,4}. Unfortunately you cannot store both, full voxel and sampled voxel information in the same cst - you either have to overwrite your initial cst or create a new one (e.g cst2) that stores only the sampled voxels and is then passed to fluenceOptimization.
Thanks for the idea! I now did the trick to duplicate the structures, ones for visualization, the undersampled ones just for optimization.
So now the optimization result is looking better :), wUnsequenced has values in the right order of magnitude.
Next step: take a look at what @tobiasbecher mentioned about the penalties. You propose that instead of taking problem.Weight, we take the number of voxels as penalty?
Below the current version of the script.
% The code expects patientFolder to contain a TROTS *.mat file
% See e.g. https://zenodo.org/records/2708302/files/Protons.zip?download=1
% To open the result, if you did not close matRad, just open matRadGUI and
% click Refresh. No need to touch nothing else.
% If you already closed MATLAB, open it again, start matRadGUI, and then
% click on the LoadMat button, and select mini.mat, not the original TROTS
% Authored by:
% S. Tattenberg - TRIUMF and NOSM
% F. Hueso-González - IFIC (CSIC/UV)
% Assistance:
% N. Wahl, T. Becher, A. Hammi
% See https://github.com/e0404/matRad/issues/695
clear, clc, close all
patientFolder = '/tmp/'; % with TROTS mat file
TrotsMatFile = patientFolder + "Protons_01.mat";
load(TrotsMatFile);
%TROTS has: data, patient, problem, problem_lex, solutionX
%matRad needs: cst, ct, patientFolder, pln, stf, dij
%% Define CT structure.
ct.resolution.x = patient.Resolution(1);
ct.resolution.y = patient.Resolution(2);
ct.resolution.z = patient.Resolution(3);
nRows = size(patient.CT, 2); % DICOM Rows, goes with y, is index 2 because of how matrix is stored
nColumns = size(patient.CT, 1);% DICOM Columns, goes with x, is index 1 because of how matrix is stored
nSlices = size(patient.CT, 3);% DICOM slices, goes with z
ct.x = linspace(patient.Offset(1), patient.Offset(1)+(nColumns-1)*patient.Resolution(1), nColumns);
ct.y = linspace(patient.Offset(2), patient.Offset(2)+(nRows-1)*patient.Resolution(2), nRows);
ct.z = linspace(patient.Offset(3), patient.Offset(3)+(nSlices-1)*patient.Resolution(3), nSlices);
ct.cubeDim = [nRows nColumns nSlices];
ct.numOfCtScen = 1;
ct.timeStamp = string(datetime("now"));
ct.cubeHU = {double(permute(patient.CT,[2 1 3]))};
%% Create cst structure
% see https://github.com/e0404/matRad/blob/master/dicom/matRad_createCst.m
% and https://github.com/e0404/matRad/blob/master/dicom/matRad_dummyCst.m
% and https://github.com/e0404/matRad/wiki/The-cst-cell
nStructures = length(patient.StructureNames);
defaultColors = colorcube(nStructures);
cst = cell(nStructures*2,6);
[grx,gry] = ndgrid(ct.x,ct.y);
disp('Calculating contours')
for i = 1:nStructures
cst{i,1} = i-1;
cst{i,2} = patient.StructureNames{i};
if contains(cst{i,2}, 'ctv', 'IgnoreCase', true) ...
|| contains(cst{i,2}, 'ptv', 'IgnoreCase', true)
cst{i,3} = 'TARGET';
else
cst{i,3} = 'OAR';
end
linvoxs = [];
disp(cst{i,2})
for s = 1:nSlices
csXY = patient.Contours{1,1}{s,i};
for sc = 1:length(csXY)
contoursXY = csXY{sc};
in = inpolygon(grx,gry,contoursXY(:,1),contoursXY(:,2)).';
ind = find(in) + (s-1)*nRows*nColumns;
linvoxs = cat(1, linvoxs, ind);
end
if ~isempty(csXY)
disp(s)
end
end
cst{i,4}{1} = linvoxs;
if strcmp(cst{i,3}, 'OAR')
cst{i,5}.Priority = 2;
else
cst{i,5}.Priority = 1;
end
cst{i,5}.alphaX = 0.1;
cst{i,5}.betaX = 0.05;
cst{i,5}.Visible = 1;
cst{i,5}.visibleColor = defaultColors(i,:);
cst{i,6} = [];
end
% Define now sparse structures for optimization
for i = 1:nStructures
cst{i+nStructures,1} = i+nStructures-1;
cst{i+nStructures,2} = ['sp', patient.StructureNames{i}];
if contains(cst{i+nStructures,2}, 'ctv', 'IgnoreCase', true) ...
|| contains(cst{i+nStructures,2}, 'ptv', 'IgnoreCase', true)
cst{i+nStructures,3} = 'TARGET';
else
cst{i+nStructures,3} = 'OAR';
end
sampVox = cell2mat(patient.SampledVoxels(i));
linvoxs = sub2ind(ct.cubeDim, sampVox(2,:), sampVox(1,:), sampVox(3,:));
disp(cst{i+nStructures,2})
cst{i+nStructures,4}{1} = linvoxs;
if strcmp(cst{i+nStructures,3}, 'OAR')
cst{i+nStructures,5}.Priority = 2;
else
cst{i+nStructures,5}.Priority = 1;
end
cst{i+nStructures,5}.alphaX = 0.1;
cst{i+nStructures,5}.betaX = 0.05;
cst{i+nStructures,5}.Visible = 1;
cst{i+nStructures,5}.visibleColor = defaultColors(i,:);
cst{i+nStructures,6} = [];
end
%% Define constraints and objectives in cst
for i=1:nStructures
OARIndex = i + nStructures;
totalIndices = size(problem, 2);
objectiveIndex = 1;
for index=1:totalIndices
if contains(cst(OARIndex,2),problem(index).Name)
cst{OARIndex,6}{objectiveIndex} = struct();
if problem(index).IsConstraint == 1
cst{OARIndex,6}{objectiveIndex}.className = 'DoseConstraints.matRad_MinMaxDose';
if problem(index).Minimise == 1
cst{OARIndex,6}{objectiveIndex}.parameters = cell({0,problem(index).Objective,1}); % 1 is approx, 2 is voxel-wise
else
cst{OARIndex,6}{objectiveIndex}.parameters = cell({problem(index).Objective,100,1});
end
cst{OARIndex,6}{objectiveIndex}.epsilon = 1.0000e-03;
else
if problem(index).Minimise == 1
cst{OARIndex,6}{objectiveIndex}.className = 'DoseObjectives.matRad_SquaredOverdosing';
else
cst{OARIndex,6}{objectiveIndex}.className = 'DoseObjectives.matRad_SquaredUnderdosing';
end
cst{OARIndex,6}{objectiveIndex}.parameters = cell({problem(index).Objective});
cst{OARIndex,6}{objectiveIndex}.penalty = problem(index).Weight;
end
objectiveIndex = objectiveIndex + 1;
end
end
end
%% clear and save
%clearvars -except cst ct patientFolder
%save([patientFolder 'mini.mat'], '-v7.3')
%% Definition of pln
pln = struct;
pln.propStf = struct;
pln.propStf.bixelWidth = ct.resolution.x;
numBeams = size(patient.Beams.BeamConfig,2);
pln.propStf.gantryAngles = [];
for beamIndex=1:numBeams
pln.propStf.gantryAngles(beamIndex) = patient.Beams.BeamConfig(beamIndex).Gantry;
end
numCouches= size(patient.Beams.BeamConfig,2);
pln.propStf.couchAngles = [];
for couchIndex=1:numCouches
pln.propStf.couchAngles(couchIndex) = patient.Beams.BeamConfig(couchIndex).Couch;
end
pln.propStf.numOfBeams = numBeams;
for i=1:numBeams
img_isox = patient.Isocentre(1) - ct.x(1);% matrad subtracts offset
img_isoy = patient.Isocentre(2) - ct.y(1);% matrad subtracts offset
img_isoz = patient.Isocentre(3) - ct.z(1);% matrad subtracts offset
pln.propStf.isoCenter(i,:) = [img_isox img_isoy img_isoz]; % not sure what this is
end
pln.voxelDimensions = ct.cubeDim;
pln.numOfVoxels = prod(pln.voxelDimensions);
pln.numOfFractions = 1;
if contains(TrotsMatFile, 'Proton', 'IgnoreCase', true)
pln.radiationMode = 'protons';
pln.propOpt.bioOptimization = 'const_RBExD';
else
pln.radiationMode = 'photons';
pln.propOpt.bioOptimization = 'none';
end
pln.machine = 'Generic';
pln.propOpt.runDAO = 0;
pln.propOpt.runSequencing = 0;
%% Pass patient TROTS Dij to MATLAB
spots = size(solutionX, 1);
dij.doseGrid.resolution.x = ct.resolution.x;
dij.doseGrid.resolution.y = ct.resolution.y;
dij.doseGrid.resolution.z = ct.resolution.z;
dij.doseGrid.x = ct.x;
dij.doseGrid.y = ct.y;
dij.doseGrid.z = ct.z;
dij.doseGrid.dimensions = ct.cubeDim;
dij.doseGrid.numOfVoxels = prod(dij.doseGrid.dimensions);
dij.ctGrid.resolution.x = ct.resolution.x;
dij.ctGrid.resolution.y = ct.resolution.y;
dij.ctGrid.resolution.z = ct.resolution.z;
dij.ctGrid.x = ct.x;
dij.ctGrid.y = ct.y;
dij.ctGrid.z = ct.z;
dij.ctGrid.dimensions = ct.cubeDim;
dij.ctGrid.numOfVoxels = prod(dij.ctGrid.dimensions);
dij.numOfBeams = pln.propStf.numOfBeams;
dij.numOfScenarios = 1; % for the moment we exclude the 9 scenarios TROTS
dij.numOfRaysPerBeam = patient.Beams.ElementIndex;
dij.totalNumOfBixels = spots;% sum([stf(:).totalNumOfBixels]);
dij.totalNumOfRays = sum(dij.numOfRaysPerBeam);
if contains(patientFolder,'Proton', 'IgnoreCase', true)
dij.RBE = 1.1;
end
dij.bixelNum = ones(spots,1);
dij.rayNum = [];
for r=1:length(dij.numOfRaysPerBeam)
dij.rayNum = [dij.rayNum; (1:1:dij.numOfRaysPerBeam(r))'];
end
dij.beamNum = [];
for b=1:dij.numOfBeams
dij.beamNum = [dij.beamNum; ones(patient.Beams.ElementIndex(b),1)*b];
end
% Allocate space for dij.physicalDose sparse matrix
doseNames = struct2cell(data.matrix);
for i = 1:dij.numOfScenarios
dij.physicalDose{i} = spalloc(dij.doseGrid.numOfVoxels,spots,1);
for structIdx = 1:nStructures
disp(patient.StructureNames(structIdx));
sampVox = cell2mat(patient.SampledVoxels(structIdx));
svoxels = size(sampVox, 2);
if length(find(strcmpi(doseNames(1,:),patient.StructureNames(structIdx)))) > 1
disp('Warning: Multiple Dij for one struct not yet supported, taking first one')
end
doseIdx = find(strcmpi(doseNames(1,:),patient.StructureNames(structIdx)), 1);%TODO handle double ones!
struct_dij = data.matrix(doseIdx).A;
voxels = size(struct_dij,1);
if data.matrix(doseIdx).b ~= 0
disp('Plan recalculation not yet supported')
return
end
if ~isempty(data.matrix(doseIdx).c)
disp('Quadratic cost functions not yet supported')
return
end
if mod(voxels,9) == 0 && voxels/9 == svoxels
disp('Warning: Multiple scenarios not yet supported')
struct_dij = struct_dij(1:9:end,:);
voxels = size(struct_dij,1);
end
if voxels ~= svoxels
disp('Wrong number of voxels')
return
end
sInd = sub2ind(ct.cubeDim, sampVox(2,:), sampVox(1,:), sampVox(3,:));
dij.physicalDose{1}(sInd,:) = struct_dij;
end
end
%% Optimize
%pln.propOpt.optimizer = 'fmincon';
[resultGUI, optimizer] = matRad_fluenceOptimization(dij, cst, pln);
solutionM = resultGUI.wUnsequenced;
dvh = matRad_calcDVH(cst(nStructures+1:end,:), resultGUI.physicalDose);
%qi = matRad_calcQualityIndicators(cst, pln, resultGUI.physicalDose);
figure;
for i = 1:length(dvh)
plot(dvh(i).doseGrid.', dvh(i).volumePoints.', 'color', defaultColors(i,:));
hold on;
end
hold off;
matRadGUI;
%clearvars -except cst ct dij pln resultGUI patientFolder
%save([patientFolder 'full.mat'], '-v7.3')And the resulting DVHs:
This shouldnt actually matter since you use only the sampled voxels
Thanks for the help, for the time being I have a workable version of the script now that gets reasonable results!