Closed ferdymercury closed 9 months ago
wahln
commented
9 months ago 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!
ferdymercury
commented
9 months ago 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
tobiasbecher
commented
9 months ago 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 :)
tobiasbecher
commented
9 months ago Two more comments:
ferdymercury
commented
9 months ago 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?
wahln
commented
9 months ago I will just enter this conversation and say stupid things that confuse everybody due to lack of in-depth knowledge about the TROTS format.
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.beamNum will, 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.rayNum will 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.bixelNum will 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));
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
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.
ferdymercury
commented
9 months ago
- 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!
tobiasbecher
commented
9 months ago 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
wahln
commented
9 months ago dij.beamNum is needed for the calculation of the beam-wise cubes in the end.
If there is undersampling, I have twosuggestions.
ferdymercury
commented
9 months ago 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.mIf 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.
Primzahl123
commented
9 months ago 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.
ferdymercury
commented
9 months ago 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.
tobiasbecher
commented
9 months ago First 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:
However, this should not lead to issues with an actual plan calculation. The bigger problem is
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.
wahln
commented
9 months ago 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.
tobiasbecher
commented
9 months ago @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
commented
9 months ago @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.
ferdymercury
commented
9 months ago Are there really also targets with undersampled voxels? These would be the only ones where minDose dose make sense.
To some extent, see:
ferdymercury
commented
9 months ago 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?
And the resulting DVHs:
tobiasbecher
commented
9 months ago This shouldnt actually matter since you use only the sampled voxels
ferdymercury
commented
9 months ago Thanks for the help, for the time being I have a workable version of the script now that gets reasonable results!
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:
```Matlab % 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; end ```I 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.