Skip to content

Commit

Permalink
4D update to dose engines. Also changes optimization
Browse files Browse the repository at this point in the history
  • Loading branch information
@wahln
wahln committed Apr 18, 2024
1 parent aa5bf3e commit fd85a05
Show file tree
Hide file tree
Showing 11 changed files with 234 additions and 118 deletions.
2 changes: 1 addition & 1 deletion 4D/matRad_calc4dDose.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,7 +89,7 @@
resultGUI.accSqrtBetaDose = matRad_doseAcc(ct,resultGUI.phaseSqrtBetaDose, cst, accType);

% only compute where we have biologically defined tissue
ix = dij.ax(:,1)~=0;
ix = any(cell2mat(cellfun(@(ax) ax ~= 0,dij.ax','UniformOutput',false)),2);

resultGUI.accEffect = resultGUI.accAlphaDose + resultGUI.accSqrtBetaDose.^2;

Expand Down
44 changes: 31 additions & 13 deletions doseCalc/+DoseEngines/@matRad_DoseEngineBase/initDoseCalc.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -136,22 +136,40 @@

if isempty(this.voxelSubIx)
% take only voxels inside patient
VctGrid = [cst{:,4}];
VctGrid = unique(vertcat(VctGrid{:}));
tmpVctGridScen = cell(1,ct.numOfCtScen);
for s = 1:ct.numOfCtScen
tmpScen = cellfun(@(c) c{s},cst(:,4),'UniformOutput',false);
tmpVctGridScen{s} = unique(vertcat(tmpScen{:}));
end
else
VctGrid = this.voxelSubIx;
if iscell(this.voxelSubIx)
tmpVctGridScen = cell(1,ct.numOfCtScen);
for s = 1:ct.numOfCtScen
tmpVctGridScen{s} = this.voxelSubIx;
end
else
tmpVctGridScen = this.voxelSubIx;
end
end
this.VctGrid = unique(vertcat(tmpVctGridScen{:}));
% No we find the subindexes for the indivdual scenarios. This helps us
% doing a subselection later on.
this.VctGridScenIx = cellfun(@(c) ismember(this.VctGrid,c),tmpVctGridScen,'UniformOutput',false);


tmpVdoseGridScen = cell(1,ct.numOfCtScen);
for s = 1:ct.numOfCtScen
% receive linear indices and grid locations from the dose grid
tmpCube = zeros(ct.cubeDim);
tmpCube(tmpVctGridScen{s}) = 1;

% interpolate cube
tmpVdoseGridScen{s} = find(matRad_interp3(dij.ctGrid.x, dij.ctGrid.y, dij.ctGrid.z,tmpCube, ...
dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z,'nearest'));
end
this.VdoseGrid = unique(vertcat(tmpVdoseGridScen{:}));
this.VdoseGridScenIx = cellfun(@(c) ismember(this.VdoseGrid,c), tmpVdoseGridScen,'UniformOutput',false);

% receive linear indices and grid locations from the dose grid
tmpCube = zeros(ct.cubeDim);
tmpCube(VctGrid) = 1;
% interpolate cube
this.VdoseGrid = find(matRad_interp3(dij.ctGrid.x, dij.ctGrid.y, dij.ctGrid.z,tmpCube, ...
dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z,'nearest'));

% save vct grid as own property in order to allow sub-classes
% to access it
this.VctGrid = VctGrid;

% Convert CT subscripts to linear indices.
[this.yCoordsV_vox, this.xCoordsV_vox, this.zCoordsV_vox] = ind2sub(ct.cubeDim,this.VctGrid);
Expand Down
4 changes: 3 additions & 1 deletion doseCalc/+DoseEngines/@matRad_DoseEngineBase/matRad_DoseEngineBase.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,9 @@
%offset; % offset adjustment for isocenter

VctGrid; % voxel grid inside patient
VdoseGrid; % voxel dose grid
VdoseGrid; % voxel dose grid
VctGridScenIx; %logical subindexes of scenarios in ct grid
VdoseGridScenIx; %logical subindexes of scenarios in dose grid

VctGridMask; % voxel grid inside patient as logical mask
VdoseGridMask; % voxel dose grid inside patient as logical mask
Expand Down
10 changes: 8 additions & 2 deletions doseCalc/+DoseEngines/matRad_ParticleAnalyticalBortfeldEngine.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,7 +96,8 @@ function chooseLateralModel(this)
elseif ~strcmp(this.lateralModel,'single')
matRad_cfg.dispWarning('Engine only supports analytically computed singleGaussian lateral Model!');
this.lateralModel = 'single';
end
end
matRad_cfg = MatRad_Config.instance();

matRad_cfg.dispInfo('Using an analytically computed %s Gaussian pencil-beam kernel model!\n');
end
Expand Down Expand Up @@ -319,11 +320,16 @@ function chooseLateralModel(this)
end

function calcLateralParticleCutOff(this,cutOffLevel,~)
calcRange = false;
if ~isfield(this.machine.data,'range')
calcRange = true;
end

for i = 1:numel(this.machine.data)
this.machine.data(i).LatCutOff.CompFac = 1-cutOffLevel;
this.machine.data(i).LatCutOff.numSig = sqrt(2) * sqrt(gammaincinv(0.995,1)); %For a 2D symmetric gaussian we need the inverse of the incomplete Gamma function for defining the CutOff
this.machine.data(i).LatCutOff.maxSigmaIni = max([this.machine.data(i).initFocus(:).SisFWHMAtIso]) ./ 2.3548;
if ~isfield(this.machine.data(i),'range')
if calcRange
this.machine.data(i).range = 10 * this.alpha*this.machine.data(i).energy.^this.p;
end
this.machine.data(i).LatCutOff.CutOff = this.machine.data(i).LatCutOff.numSig * sqrt(this.machine.data(i).LatCutOff.maxSigmaIni^2 + this.calcSigmaLatMCS(this.machine.data(i).range,this.machine.data(i).energy)^2);
Expand Down
20 changes: 14 additions & 6 deletions doseCalc/+DoseEngines/matRad_ParticleFineSamplingPencilBeamEngine.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,17 +71,25 @@ function setDefaults(this)
lateralRayCutOff = this.getLateralDistanceFromDoseCutOffOnRay(ray);

% Ray tracing for beam i and ray j
[ix,ray.radialDist_sq,ray.latDists] = this.calcGeoDists(currBeam.bevCoords, ...
[ix,radialDist_sq,latDists] = this.calcGeoDists(currBeam.bevCoords, ...
ray.sourcePoint_bev, ...
ray.targetPoint_bev, ...
ray.SAD, ...
currBeam.ixRadDepths, ...
currBeam.validCoordsAll, ...
lateralRayCutOff);

%ray.geoDepths = cellfun(@(rD) rD(ix),currBeam.geoDepths,'UniformOutput',false); %Not needed for this engine
ray.radDepths = cellfun(@(rD) rD(ix),currBeam.radDepths,'UniformOutput',false);
ray.ix = currBeam.ixRadDepths(ix);
ray.subIxVdoseGrid = currBeam.subIxVdoseGrid(ix);
ray.validCoords = cellfun(@(beamIx) beamIx & ix,currBeam.validCoords,'UniformOutput',false);
ray.ix = cellfun(@(ixInGrid) this.VdoseGrid(ixInGrid),ray.validCoords,'UniformOutput',false);

%subCoords = cellfun(@(beamIx) beamIx(ix),currBeam.validCoords,'UniformOutput',false);
%ray.radialDist_sq = cellfun(@(subix) radialDist_sq(subix),radialDist_sq,subCoords);
ray.radialDist_sq = cellfun(@(beamIx) radialDist_sq(beamIx(ix)),currBeam.validCoords,'UniformOutput',false);
ray.latDists = cellfun(@(beamIx) latDists(beamIx(ix),:),currBeam.validCoords,'UniformOutput',false);

ray.validCoordsAll = any(cell2mat(ray.validCoords),2);

ray.geoDepths = cellfun(@(rD,ix) rD(ix),currBeam.geoDepths,ray.validCoords,'UniformOutput',false); %usually not needed for particle beams
ray.radDepths = cellfun(@(rD,ix) rD(ix),currBeam.radDepths,ray.validCoords,'UniformOutput',false);
end

function bixel = initBixel(this,currRay,k)
Expand Down
85 changes: 56 additions & 29 deletions doseCalc/+DoseEngines/matRad_ParticlePencilBeamEngineAbstract.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -235,17 +235,26 @@ function chooseLateralModel(this)
lateralRayCutOff = this.getLateralDistanceFromDoseCutOffOnRay(ray);

% Ray tracing for beam i and ray j
[ix,ray.radialDist_sq] = this.calcGeoDists(currBeam.bevCoords, ...
[ix,radialDist_sq] = this.calcGeoDists(currBeam.bevCoords, ...
ray.sourcePoint_bev, ...
ray.targetPoint_bev, ...
ray.SAD, ...
currBeam.ixRadDepths, ...
currBeam.validCoordsAll, ...
lateralRayCutOff);

ray.validCoords = cellfun(@(beamIx) beamIx & ix,currBeam.validCoords,'UniformOutput',false);
ray.ix = cellfun(@(ixInGrid) this.VdoseGrid(ixInGrid),ray.validCoords,'UniformOutput',false);

%subCoords = cellfun(@(beamIx) beamIx(ix),currBeam.validCoords,'UniformOutput',false);
%ray.radialDist_sq = cellfun(@(subix) radialDist_sq(subix),radialDist_sq,subCoords);
ray.radialDist_sq = cellfun(@(beamIx) radialDist_sq(beamIx(ix)),currBeam.validCoords,'UniformOutput',false);

ray.validCoordsAll = any(cell2mat(ray.validCoords),2);

%ray.geoDepths = cellfun(@(rD) rD(ix),currBeam.geoDepths,'UniformOutput',false); %usually not needed for particle beams
ray.radDepths = cellfun(@(rD) rD(ix),currBeam.radDepths,'UniformOutput',false);
ray.ix = currBeam.ixRadDepths(ix);
ray.subIxVdoseGrid = currBeam.subIxVdoseGrid(ix);
ray.geoDepths = cellfun(@(rD,ix) rD(ix),currBeam.geoDepths,ray.validCoords,'UniformOutput',false); %usually not needed for particle beams
ray.radDepths = cellfun(@(rD,ix) rD(ix),currBeam.radDepths,ray.validCoords,'UniformOutput',false);
%ray.ix = currBeam.ixRadDepths(ix);
%ray.subIxVdoseGrid = currBeam.subIxVdoseGrid(ix);
end

function [currBixel] = getBixelIndicesOnRay(this,currBixel,currRay)
Expand Down Expand Up @@ -391,11 +400,16 @@ function chooseLateralModel(this)
radDepthOffset = this.machine.data(maxEnergyIx).offset + minRaShi;

% apply limit in depth
subSelectIx = currBeam.radDepths{1} < (this.machine.data(maxEnergyIx).depths(end) - radDepthOffset);
currBeam.ixRadDepths = currBeam.ixRadDepths(subSelectIx);
currBeam.subIxVdoseGrid = currBeam.subIxVdoseGrid(subSelectIx);
currBeam.radDepths = cellfun(@(rd) rd(subSelectIx),currBeam.radDepths,'UniformOutput',false);
currBeam.bevCoords = currBeam.bevCoords(subSelectIx,:);
%subSelectIx = currBeam.radDepths{1} < (this.machine.data(maxEnergyIx).depths(end) - radDepthOffset);

subSelectIx = cellfun(@(rD) rD < (this.machine.data(maxEnergyIx).depths(end) - radDepthOffset),currBeam.radDepths,'UniformOutput',false);
currBeam.validCoords = cellfun(@and,subSelectIx,currBeam.validCoords,'UniformOutput',false);
currBeam.validCoordsAll = any(cell2mat(currBeam.validCoords),2);

%currBeam.ixRadDepths = currBeam.ixRadDepths(subSelectIx);
%currBeam.subIxVdoseGrid = currBeam.subIxVdoseGrid(subSelectIx);
%currBeam.radDepths = cellfun(@(rd) rd(subSelectIx),currBeam.radDepths,'UniformOutput',false);
%currBeam.bevCoords = currBeam.bevCoords(subSelectIx,:);

%Precompute CutOff
this.calcLateralParticleCutOff(this.dosimetricLateralCutOff,currBeam);
Expand All @@ -406,26 +420,31 @@ function chooseLateralModel(this)

matRad_cfg.dispInfo('Initializing biological dose calculation...\n');

dij.ax = zeros(dij.doseGrid.numOfVoxels,1);
dij.bx = zeros(dij.doseGrid.numOfVoxels,1);
numOfCtScen = numel(this.VdoseGridScenIx);


cstDownsampled = matRad_setOverlapPriorities(cst);

% resizing cst to dose cube resolution
cstDownsampled = matRad_resizeCstToGrid(cstDownsampled,dij.ctGrid.x,dij.ctGrid.y,dij.ctGrid.z,...
dij.doseGrid.x,dij.doseGrid.y,dij.doseGrid.z);
% retrieve photon LQM parameter for the current dose grid voxels
[dij.ax,dij.bx] = matRad_getPhotonLQMParameters(cstDownsampled,dij.doseGrid.numOfVoxels,1,this.VdoseGrid);

% vAlphaX and vBetaX for parameters in VdoseGrid
this.vAlphaX = dij.ax(this.VdoseGridMask);
this.vBetaX = dij.bx(this.VdoseGridMask);
this.vTissueIndex = zeros(size(this.VdoseGrid,1),1);
tmpScenVdoseGrid = cell(numOfCtScen,1);

if strcmp(this.bioParam.model,'LEM')
matRad_cfg.dispInfo('\tUsing LEM model with precomputed kernels\n');
[dij.ax,dij.bx] = matRad_getPhotonLQMParameters(cstDownsampled,dij.doseGrid.numOfVoxels,this.VdoseGrid);


for s = 1:numOfCtScen
tmpScenVdoseGrid{s} = this.VdoseGrid(this.VdoseGridScenIx{s});
% retrieve photon LQM parameter for the current dose grid voxels

% vAlphaX and vBetaX for parameters in VdoseGrid
this.vAlphaX{s} = dij.ax{s}(tmpScenVdoseGrid{s});
this.vBetaX{s} = dij.bx{s}(tmpScenVdoseGrid{s});
this.vTissueIndex{s} = zeros(size(tmpScenVdoseGrid{s},1),1);
end

if strcmp(this.bioParam.model,'LEM')
matRad_cfg.dispInfo('\tUsing LEM model with precomputed kernels\n');

if isfield(this.machine.data,'alphaX') && isfield(this.machine.data,'betaX')
for i = 1:size(cstDownsampled,1)
Expand All @@ -439,15 +458,20 @@ function chooseLateralModel(this)

% check consitency of biological baseData and cst settings
if ~isempty(IdxTissue)
isInVdoseGrid = ismember(this.VdoseGrid,cstDownsampled{i,4}{1});
this.vTissueIndex(isInVdoseGrid) = IdxTissue;
for s = 1:numOfCtScen
tmpScenVdoseGrid = this.VdoseGrid(this.VdoseGridScenIx{s});
isInVdoseGrid = ismember(tmpScenVdoseGrid,cstDownsampled{i,4}{s});
this.vTissueIndex{s}(isInVdoseGrid) = IdxTissue;
end
else
matRad_cfg.dispError('Biological base data and cst are inconsistent!');
end

else
this.vTissueIndex(row) = 1;
matRad_cfg.dispInfo('\tTissue type of %s was set to 1\n',cstDownsampled{i,2});
for s = 1:numOfCtScen
this.vTissueIndex{s}(:) = 1;
end
matRad_cfg.dispWarning('\tTissue type of %s was set to 1\n',cstDownsampled{i,2});
end
end
dij.vTissueIndex = this.vTissueIndex;
Expand All @@ -464,6 +488,7 @@ function chooseLateralModel(this)
else
matRad_cfg.dispError('Unknown Biological Model!');
end

end

function dij = allocateBioDoseContainer(this,dij)
Expand Down Expand Up @@ -896,9 +921,11 @@ function calcLateralParticleCutOff(this,cutOffLevel,stfElement)

% just use tissue classes of voxels found by ray tracer
if this.calcBioDose
ray.vTissueIndex = this.vTissueIndex(ray.subIxVdoseGrid,:);
ray.vAlphaX = this.vAlphaX(ray.subIxVdoseGrid);
ray.vBetaX = this.vBetaX(ray.subIxVdoseGrid);
for s = 1:numel(this.vTissueIndex)
ray.vTissueIndex{s} = this.vTissueIndex{s}(ray.validCoords{s},:);
ray.vAlphaX{s} = this.vAlphaX{s}(ray.validCoords{s});
ray.vBetaX{s} = this.vBetaX{s}(ray.validCoords{s});
end
end
end

Expand Down
Loading

0 comments on commit fd85a05

Please sign in to comment.