diff --git a/4D/matRad_calc4dDose.m b/4D/matRad_calc4dDose.m index c54723d5a..1249d54e8 100644 --- a/4D/matRad_calc4dDose.m +++ b/4D/matRad_calc4dDose.m @@ -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; diff --git a/doseCalc/+DoseEngines/@matRad_DoseEngineBase/initDoseCalc.m b/doseCalc/+DoseEngines/@matRad_DoseEngineBase/initDoseCalc.m index d1c905e0d..681f6ddb4 100644 --- a/doseCalc/+DoseEngines/@matRad_DoseEngineBase/initDoseCalc.m +++ b/doseCalc/+DoseEngines/@matRad_DoseEngineBase/initDoseCalc.m @@ -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); diff --git a/doseCalc/+DoseEngines/@matRad_DoseEngineBase/matRad_DoseEngineBase.m b/doseCalc/+DoseEngines/@matRad_DoseEngineBase/matRad_DoseEngineBase.m index 1917cc9f9..f81f56fbd 100644 --- a/doseCalc/+DoseEngines/@matRad_DoseEngineBase/matRad_DoseEngineBase.m +++ b/doseCalc/+DoseEngines/@matRad_DoseEngineBase/matRad_DoseEngineBase.m @@ -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 diff --git a/doseCalc/+DoseEngines/matRad_ParticleAnalyticalBortfeldEngine.m b/doseCalc/+DoseEngines/matRad_ParticleAnalyticalBortfeldEngine.m index 5f07583d0..052c4d569 100644 --- a/doseCalc/+DoseEngines/matRad_ParticleAnalyticalBortfeldEngine.m +++ b/doseCalc/+DoseEngines/matRad_ParticleAnalyticalBortfeldEngine.m @@ -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 @@ -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); diff --git a/doseCalc/+DoseEngines/matRad_ParticleFineSamplingPencilBeamEngine.m b/doseCalc/+DoseEngines/matRad_ParticleFineSamplingPencilBeamEngine.m index 34e0148c7..f958652e1 100644 --- a/doseCalc/+DoseEngines/matRad_ParticleFineSamplingPencilBeamEngine.m +++ b/doseCalc/+DoseEngines/matRad_ParticleFineSamplingPencilBeamEngine.m @@ -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) diff --git a/doseCalc/+DoseEngines/matRad_ParticlePencilBeamEngineAbstract.m b/doseCalc/+DoseEngines/matRad_ParticlePencilBeamEngineAbstract.m index ade8e5619..26b58e188 100644 --- a/doseCalc/+DoseEngines/matRad_ParticlePencilBeamEngineAbstract.m +++ b/doseCalc/+DoseEngines/matRad_ParticlePencilBeamEngineAbstract.m @@ -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) @@ -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); @@ -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) @@ -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; @@ -464,6 +488,7 @@ function chooseLateralModel(this) else matRad_cfg.dispError('Unknown Biological Model!'); end + end function dij = allocateBioDoseContainer(this,dij) @@ -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 diff --git a/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m b/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m index 8496d35ac..05318d032 100644 --- a/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m +++ b/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m @@ -134,6 +134,9 @@ function setDefaults(this) scenRay = currRay; scenRay.radDepths = scenRay.radDepths{ctScen}; scenRay.radDepths = (1+this.multScen.relRangeShift(scenNum))*scenRay.radDepths + this.multScen.absRangeShift(scenNum); + scenRay.radialDist_sq = scenRay.radialDist_sq{ctScen}; + scenRay.ix = scenRay.ix{ctScen}; + if this.multScen.absRangeShift(scenNum) < 0 %TODO: better way to handle this? scenRay.radDepths(scenRay.radDepths < 0) = 0; @@ -143,6 +146,20 @@ function setDefaults(this) scenRay.geoDepths = scenRay.geoDepths{ctScen}; end + if isfield(scenRay,'latDists') + scenRay.latDists = scenRay.latDists{ctScen}; + end + + if isfield(scenRay,'isoLatDists') + scenRay.isoLatDists = scenRay.isoLatDists{ctScen}; + end + + if isfield(scenRay,'vTissueIndex') + scenRay.vTissueIndex = scenRay.vTissueIndex{ctScen}; + scenRay.vAlphaX = scenRay.vAlphaX{ctScen}; + scenRay.vBetaX = scenRay.vBetaX{ctScen}; + end + for k = 1:currRay.numOfBixels %Bixel Computation currBixel = this.computeBixel(scenRay,k); @@ -310,14 +327,13 @@ function setDefaults(this) %radDepthsMat = cellfun(@(radDepthCube) matRad_interp3(dij.ctGrid.x, dij.ctGrid.y, dij.ctGrid.z,radDepthCube,dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z,'nearest'),radDepthsMat,'UniformOutput',false); %Find valid coordinates - %TODO: currently we are taking coordinates valid in the first - %cube. Might cause issues in other cubes - coordIsValid = ~isnan(radDepthVdoseGrid{1}); - currBeam.subIxVdoseGrid = find(coordIsValid); - currBeam.ixRadDepths = this.VdoseGrid(coordIsValid); - currBeam.radDepths = cellfun(@(rd) rd(coordIsValid),radDepthVdoseGrid,'UniformOutput',false); - currBeam.geoDepths = cellfun(@(gd) gd(coordIsValid),geoDistVdoseGrid,'UniformOutput',false); - currBeam.bevCoords = rot_coordsVdoseGrid(coordIsValid,:); + coordIsValid = cellfun(@isfinite, radDepthVdoseGrid,'UniformOutput',false); %Reduce coordinates for finite values + currBeam.validCoords = cellfun(@and,coordIsValid,this.VdoseGridScenIx,'UniformOutput',false); %Reduce coordinates according to scenario + currBeam.validCoordsAll = any(cell2mat(coordIsValid),2); + + currBeam.radDepths = radDepthVdoseGrid; + currBeam.geoDepths = geoDistVdoseGrid; + currBeam.bevCoords = rot_coordsVdoseGrid; % compute SSDs currBeam = matRad_computeSSD(currBeam,ct,'densityThreshold',this.ssdDensityThreshold); @@ -368,18 +384,30 @@ function setDefaults(this) lateralRayCutOff = this.getLateralDistanceFromDoseCutOffOnRay(ray); % Ray tracing for beam i and ray j - [ix,ray.radialDist_sq,ray.latDists,ray.isoLatDists] = this.calcGeoDists(currBeam.bevCoords, ... + [ix,radialDist_sq,latDists,isoLatDists] = this.calcGeoDists(currBeam.bevCoords, ... ray.sourcePoint_bev, ... ray.targetPoint_bev, ... ray.SAD, ... - currBeam.ixRadDepths, ... + currBeam.validCoordsAll, ... lateralRayCutOff); %Subindex given the relevant indices from the geometric %distance calculation - ray.geoDepths = cellfun(@(rD) rD(ix),currBeam.geoDepths,'UniformOutput',false); - ray.radDepths = cellfun(@(rD) rD(ix),currBeam.radDepths,'UniformOutput',false); - ray.ix = currBeam.ixRadDepths(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.isoLatDists = cellfun(@(beamIx) isoLatDists(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); + %ray.ix = currBeam.ixRadDepths(ix); + %ray.subIxVdoseGrid = currBeam.subIxVdoseGrid(ix); end function lateralRayCutOff = getLateralDistanceFromDoseCutOffOnRay(this,ray) @@ -587,14 +615,14 @@ function setDefaults(this) % Define function for obtain rotation matrix. if all(a==b) % rotation matrix corresponds to eye matrix if the vectors are the same - rot_coords_temp = rot_coords_bev; + rot_coords_temp = rot_coords_bev(radDepthIx,:); else % Define rotation matrix ssc = @(v) [0 -v(3) v(2); v(3) 0 -v(1); -v(2) v(1) 0]; R = eye(3) + ssc(cross(a,b)) + ssc(cross(a,b))^2*(1-dot(a,b))/(norm(cross(a,b))^2); % Rotate every CT voxel - rot_coords_temp = rot_coords_bev*R; + rot_coords_temp = rot_coords_bev(radDepthIx,:)*R; end % Put [0 0 0] position CT in center of the beamlet. @@ -613,7 +641,8 @@ function setDefaults(this) % index list within considered voxels %ix = radDepthIx(subsetMask); - ix = subsetMask; + ix = radDepthIx; + ix(ix) = subsetMask; % return radial distances squared if nargout > 1 diff --git a/matRad_calcCubes.m b/matRad_calcCubes.m index 8336c402d..6a1511347 100644 --- a/matRad_calcCubes.m +++ b/matRad_calcCubes.m @@ -112,7 +112,7 @@ wBeam = (resultGUI.w .* beamInfo(i).logIx); % consider biological optimization - ix = dij.bx(:,scenNum)~=0 & resultGUI.(['physicalDose', beamInfo(i).suffix])(:) > 0; + ix = dij.bx{scenNum}~=0 & resultGUI.(['physicalDose', beamInfo(i).suffix])(:) > 0; % Calculate effect from alpha- and sqrtBetaDose resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix]) = full(dij.(['mAlphaDose' RBE_model{j}]){scenNum} * wBeam + (dij.(['mSqrtBetaDose' RBE_model{j}]){scenNum} * wBeam).^2); @@ -120,7 +120,7 @@ % Calculate RBExD from the effect resultGUI.(['RBExD', RBE_model{j}, beamInfo(i).suffix]) = zeros(size(resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix]))); - resultGUI.(['RBExD', RBE_model{j}, beamInfo(i).suffix])(ix) = (sqrt(dij.ax(ix).^2 + 4 .* dij.bx(ix) .* resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix])(ix)) - dij.ax(ix))./(2.*dij.bx(ix)); + resultGUI.(['RBExD', RBE_model{j}, beamInfo(i).suffix])(ix) = (sqrt(dij.ax{scenNum}(ix).^2 + 4 .* dij.bx{scenNum}(ix) .* resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix])(ix)) - dij.ax{scenNum}(ix))./(2.*dij.bx{scenNum}(ix)); % Divide RBExD with the physicalDose to get the plain RBE cube resultGUI.(['RBE', RBE_model{j}, beamInfo(i).suffix]) = resultGUI.(['RBExD', RBE_model{j}, beamInfo(i).suffix])./resultGUI.(['physicalDose', beamInfo(i).suffix]); diff --git a/matRad_fluenceOptimization.m b/matRad_fluenceOptimization.m index 402078f7f..84f8fc1a3 100755 --- a/matRad_fluenceOptimization.m +++ b/matRad_fluenceOptimization.m @@ -134,14 +134,15 @@ matRad_cfg.dispInfo('chosen uniform weight of %f!\n',bixelWeight); elseif pln.bioParam.bioOpt - % retrieve photon LQM parameter - [ax,bx] = matRad_getPhotonLQMParameters(cst,dij.doseGrid.numOfVoxels,dij.numOfScenarios); - - if ~isequal(dij.ax(dij.ax~=0),ax(dij.ax~=0)) || ... - ~isequal(dij.bx(dij.bx~=0),bx(dij.bx~=0)) - matRad_cfg.dispError('Inconsistent biological parameter - please recalculate dose influence matrix!\n'); + % retrieve photon LQM parameter + [ax,bx] = matRad_getPhotonLQMParameters(cst,dij.doseGrid.numOfVoxels); + checkAxBx = cellfun(@(ax1,bx1,ax2,bx2) isequal(ax1(ax1~=0),ax2(ax1~=0)) && isequal(bx1(bx1~=0),bx2(bx1~=0)),dij.ax,dij.bx,ax,bx); + if ~all(checkAxBx) + matRad_cfg.dispError('Inconsistent biological parameters in dij.ax and/or dij.bx - please recalculate dose influence matrix before optimization!\n'); end + + for i = 1:size(cst,1) for j = 1:size(cst{i,6},2) @@ -152,8 +153,10 @@ end end - - dij.ixDose = dij.bx~=0; + + for s = 1:numel(dij.bx) + dij.ixDose{s} = dij.bx{s}~=0; + end if isequal(pln.bioParam.quantityOpt,'effect') @@ -168,8 +171,11 @@ elseif isequal(pln.bioParam.quantityOpt,'RBExD') %pre-calculations - dij.gamma = zeros(dij.doseGrid.numOfVoxels,dij.numOfScenarios); - dij.gamma(dij.ixDose) = dij.ax(dij.ixDose)./(2*dij.bx(dij.ixDose)); + for s = 1:numel(dij.ixDose) + dij.gamma{s} = zeros(dij.doseGrid.numOfVoxels,dij.numOfScenarios); + dij.gamma{s}(dij.ixDose{s}) = dij.ax{s}(dij.ixDose{s})./(2*dij.bx{s}(dij.ixDose{s})); + end + % calculate current effect in target aTmp = dij.mAlphaDose{1}*wOnes; @@ -262,6 +268,7 @@ backProjection.scenarios = ixForOpt; backProjection.scenarioProb = pln.multScen.scenProb; backProjection.nominalCtScenarios = linIxDIJ_nominalCT; +%backProjection.scenDim = pln.multScen optiProb = matRad_OptimizationProblem(backProjection); optiProb.quantityOpt = pln.bioParam.quantityOpt; diff --git a/matRad_getPhotonLQMParameters.m b/matRad_getPhotonLQMParameters.m index df32788c6..11ed01f54 100644 --- a/matRad_getPhotonLQMParameters.m +++ b/matRad_getPhotonLQMParameters.m @@ -1,15 +1,14 @@ -function [ax,bx] = matRad_getPhotonLQMParameters(cst,numVoxel,ctScen,VdoseGrid) -% matRad function to receive the photon LQM reference parameter -% +function [ax,bx] = matRad_getPhotonLQMParameters(cst,numVoxel,VdoseGrid) +% matRad function to receive the photon LQM reference parameter +% % call % [ax,bx] = matRad_getPhotonLQMParameters(cst,numVoxel,ctScen,VdoseGrid) % % input -% cst: matRad cst struct -% numVoxel: number of voxels of the dose cube -% ctScen: CT scenarios for alpha_x and beta_x should be calculated -% VdoseGrid: optional linear index vector that allows to specify subindices -% for which ax and bx will be computed +% cst: matRad cst struct +% numVoxel: number of voxels of the dose cube +% VdoseGrid: optional linear index vector that allows to specify subindices +% for which ax and bx will be computed % % output % ax: vector containing for each linear voxel index alpha_x @@ -20,37 +19,47 @@ % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % -% Copyright 2018 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 +% Copyright 2018 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. % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ax = zeros(numVoxel,ctScen); -bx = zeros(numVoxel,ctScen); +numOfCtScen = unique(cellfun(@numel,cst(:,4))); -for i = 1:size(cst,1) - if isequal(cst{i,3},'OAR') || isequal(cst{i,3},'TARGET') - for ctScen = 1:numel(cst{i,4}) - - if exist('VdoseGrid','var') - isInVdoseGrid = ismember(VdoseGrid,cst{i,4}{ctScen}); - ax(VdoseGrid(isInVdoseGrid)) = cst{i,5}.alphaX; - bx(VdoseGrid(isInVdoseGrid)) = cst{i,5}.betaX; - else - ax(cst{i,4}{ctScen},ctScen) = cst{i,5}.alphaX; - bx(cst{i,4}{ctScen},ctScen) = cst{i,5}.betaX; - end - - end - end +if numel(numOfCtScen) > 1 + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Inconsinstent number of ct scnearios in cst!'); end +ax = cell(numOfCtScen,1); +bx = cell(numOfCtScen,1); +[ax{:}] = deal(zeros(numVoxel,1)); +[bx{:}] = deal(zeros(numVoxel,1)); + +for i = 1:size(cst,1) + if isequal(cst{i,3},'OAR') || isequal(cst{i,3},'TARGET') + for s = 1:numOfCtScen + if exist('VdoseGrid','var') + if iscell(VdoseGrid) + isInVdoseGrid = ismember(VdoseGrid{s},cst{i,4}{s}); + else + isInVdoseGrid = ismember(VdoseGrid,cst{i,4}{s}); + end + ax{s}(VdoseGrid(isInVdoseGrid)) = cst{i,5}.alphaX; + bx{s}(VdoseGrid(isInVdoseGrid)) = cst{i,5}.betaX; + else + ax{s}(cst{i,4}{s}) = cst{i,5}.alphaX; + bx{s}(cst{i,4}{s}) = cst{i,5}.betaX; + end + end + end +end \ No newline at end of file diff --git a/optimization/projections/matRad_VariableRBEProjection.m b/optimization/projections/matRad_VariableRBEProjection.m index 29bf45919..36dc53177 100644 --- a/optimization/projections/matRad_VariableRBEProjection.m +++ b/optimization/projections/matRad_VariableRBEProjection.m @@ -21,7 +21,8 @@ function RBExD = computeSingleScenario(obj,dij,scen,w) effect = computeSingleScenario@matRad_EffectProjection(obj,dij,scen,w); %First compute effect RBExD = zeros(dij.doseGrid.numOfVoxels,1); - RBExD(dij.ixDose(:,scen)) = sqrt((effect(dij.ixDose(:,scen))./dij.bx(dij.ixDose(:,scen)))+(dij.gamma(dij.ixDose(:,scen)).^2)) - dij.gamma(dij.ixDose(:,scen)); + ctScen = ind2sub(size(dij.physicalDose),scen); + RBExD(dij.ixDose{ctScen}) = sqrt((effect(dij.ixDose{ctScen})./dij.bx{ctScen}(dij.ixDose{ctScen}))+(dij.gamma{ctScen}(dij.ixDose{ctScen}).^2)) - dij.gamma{ctScen}(dij.ixDose{ctScen}); end function wGrad = projectSingleScenarioGradient(obj,dij,doseGrad,scen,w) @@ -34,10 +35,13 @@ %a computation (will skip if weights are equal to cache) obj = obj.compute(dij,w); + %Get corresponding ct scenario + ctScen = ind2sub(size(dij.physicalDose),scen); + %Scaling vor variable RBExD - scaledEffect = obj.d{scen} + dij.gamma(:,scen); + scaledEffect = obj.d{scen} + dij.gamma{ctScen}; doseGradTmp = zeros(dij.doseGrid.numOfVoxels,1); - doseGradTmp(dij.ixDose(:,scen)) = doseGrad{scen}(dij.ixDose(:,scen)) ./ (2*dij.bx(dij.ixDose(:,scen)).*scaledEffect(dij.ixDose(:,scen))); + doseGradTmp(dij.ixDose{ctScen}) = doseGrad{scen}(dij.ixDose{ctScen}) ./ (2*dij.bx{ctScen}(dij.ixDose{ctScen}).*scaledEffect(dij.ixDose{ctScen})); %Now modify the effect computation vBias = (doseGradTmp' * dij.mAlphaDose{scen})'; @@ -60,8 +64,11 @@ eExpSqTerm = dij.mSqrtBetaDose{scen}*w; eExp = eExpLinTerm + eExpSqTerm.^2; + %Get corresponding ct scenario + ctScen = ind2sub(size(dij.physicalDose),scen); + RBExDexp = zeros(dij.doseGrid.numOfVoxels,1); - RBExDexp(dij.ixDose(:,scen)) = sqrt((eExp(dij.ixDose(:,scen))./dij.bx(dij.ixDose(:,scen)))+(dij.gamma(dij.ixDose(:,scen)).^2)) - dij.gamma(dij.ixDose(:,scen)); + RBExDexp(dij.ixDose{ctScen}) = sqrt((eExp(dij.ixDose{ctScen})./dij.bx{ctScen}(dij.ixDose{ctScen}))+(dij.gamma{ctScen}(dij.ixDose{ctScen}).^2)) - dij.gamma{ctScen}(dij.ixDose{ctScen}); for i = 1:size(dij.physicalDoseOmega,2) dOmegaV{scen,i} = dij.mAlphaDoseOmega{scen,i} * w; @@ -76,11 +83,14 @@ %While the dose cache should be up to date here, we ask for %a computation (will skip if weights are equal to cache) obj = obj.compute(dij,w); + + %Get corresponding ct scenario + ctScen = ind2sub(size(dij.physicalDose),scen); %Scaling vor variable RBExD - scaledEffect = obj.dExp{scen} + dij.gamma; + scaledEffect = obj.dExp{scen} + dij.gamma{ctScen}; doseGradTmp = zeros(dij.doseGrid.numOfVoxels,1); - doseGradTmp(dij.ixDose(:,scen)) = dExpGrad{scen}(dij.ixDose(:,scen)) ./ (2*dij.bx(dij.ixDose(:,scen)).*scaledEffect(dij.ixDose(:,scen))); + doseGradTmp(dij.ixDose{ctScen}) = dExpGrad{scen}(dij.ixDose{ctScen}) ./ (2*dij.bx{ctScen}(dij.ixDose{ctScen}).*scaledEffect(dij.ixDose{ctScen})); %Now modify the effect computation vBias = (doseGradTmp' * dij.mAlphaDoseExp{scen})';