diff --git a/matRad/gui/widgets/matRad_WorkflowWidget.m b/matRad/gui/widgets/matRad_WorkflowWidget.m index e66759b76..25f915940 100644 --- a/matRad/gui/widgets/matRad_WorkflowWidget.m +++ b/matRad/gui/widgets/matRad_WorkflowWidget.m @@ -424,9 +424,11 @@ function btnCalcDose_Callback(this, hObject, eventdata) % prepare dij for 3d conformal if isfield(pln,'propOpt') && isfield(pln.propOpt,'conf3D') && pln.propOpt.conf3D dij = matRad_collapseDij(dij); + stf = matRad_collapseStf(stf); end % assign results to base worksapce assignin('base','dij',dij); + assignin('base','stf',stf); catch ME diff --git a/matRad/optimization/matRad_collapseDij.m b/matRad/optimization/matRad_collapseDij.m index 03c66895e..945a8b639 100644 --- a/matRad/optimization/matRad_collapseDij.m +++ b/matRad/optimization/matRad_collapseDij.m @@ -28,14 +28,14 @@ % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -dijNew.totalNumOfBixels = 1; -dijNew.totalNumOfRays = 1; -dijNew.numOfBeams = 1; -dijNew.numOfRaysPerBeam = 1; +dijNew.totalNumOfBixels = dij.numOfBeams; +dijNew.totalNumOfRays = dij.numOfBeams; +dijNew.numOfBeams = dij.numOfBeams; +dijNew.numOfRaysPerBeam = ones(dij.numOfBeams,1); -dijNew.beamNum = 1; -dijNew.bixelNum = 1; -dijNew.rayNum = 1; +dijNew.beamNum = [1:dij.numOfBeams]'; +dijNew.bixelNum = [1:dij.numOfBeams]'; +dijNew.rayNum = [1:dij.numOfBeams]'; dijNew.doseGrid = dij.doseGrid; dijNew.ctGrid = dij.ctGrid; @@ -43,6 +43,11 @@ dijNew.numOfScenarios = dij.numOfScenarios; for i = 1:dij.numOfScenarios - dijNew.physicalDose{i} = sum(dij.physicalDose{i},2); + tmp = sparse(dij.doseGrid.numOfVoxels,dij.numOfBeams); % initialize sparse matrix + for j = 1:dij.numOfBeams + % Sum only the columns corresponding to beam j + tmp(:, j) = sum(dij.physicalDose{i}(:, dij.beamNum == j), 2); + end + dijNew.physicalDose{i} = tmp; end diff --git a/matRad/optimization/matRad_collapseStf.m b/matRad/optimization/matRad_collapseStf.m new file mode 100644 index 000000000..7213d6bdf --- /dev/null +++ b/matRad/optimization/matRad_collapseStf.m @@ -0,0 +1,36 @@ +function stf = matRad_collapseStf(stf) +% matRad collapse dij function for simulation of 3D conformal treatments. +% Function to supress intensity-modulation for photons in order to simulate +% 3D conformal treatments. +% +% call +% dijNew = matRad_collapseDij(dij) +% +% input +% dij: dose influence matrix +% +% output +% dijNew: collapsed dose influence matrix +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% 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/LICENSE.md. 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% dummy collapse so that it works with sequencing + +for i = 1:size(stf,2) + stf(i).numOfRays = 1; + stf(i).totalNumOfBixels = 1; +end \ No newline at end of file diff --git a/matRad/sequencing/matRad_engelLeafSequencing.m b/matRad/sequencing/matRad_engelLeafSequencing.m index a1d944374..bc1436574 100644 --- a/matRad/sequencing/matRad_engelLeafSequencing.m +++ b/matRad/sequencing/matRad_engelLeafSequencing.m @@ -57,12 +57,12 @@ numOfRaysPerBeam = stf(i).numOfRays; % get relevant weights for current beam - wOfCurrBeams = resultGUI.w(1+offset:numOfRaysPerBeam+offset); + wOfCurrBeams = resultGUI.w(1+offset:numOfRaysPerBeam+offset).* ones(size(stf(i).ray,2),1); - X = ones(numOfRaysPerBeam,1)*NaN; - Z = ones(numOfRaysPerBeam,1)*NaN; + X = ones(size(stf(i).ray,2),1)*NaN; %this way it also works with3dconformal + Z = ones(size(stf(i).ray,2),1)*NaN; - for j=1:stf(i).numOfRays + for j=1:size(stf(i).ray,2) X(j) = stf(i).ray(j).rayPos_bev(:,1); Z(j) = stf(i).ray(j).rayPos_bev(:,3); end @@ -352,15 +352,24 @@ clear rightIntLimit; end - - sequencing.beam(i).numOfShapes = k; - sequencing.beam(i).shapes = shapes(:,:,1:k); - sequencing.beam(i).shapesWeight = shapesWeight(1:k)/numOfLevels*calFac; - sequencing.beam(i).bixelIx = 1+offset:numOfRaysPerBeam+offset; - sequencing.beam(i).fluence = D_0; + if sum(wOfCurrBeams)>0 + + sequencing.beam(i).numOfShapes = k; + sequencing.beam(i).shapes = shapes(:,:,1:k); + sequencing.beam(i).shapesWeight = shapesWeight(1:k)/numOfLevels*calFac; + sequencing.beam(i).bixelIx = 1+offset:numOfRaysPerBeam+offset; + sequencing.beam(i).fluence = D_0; + + else + sequencing.beam(i).numOfShapes = 1; + sequencing.beam(i).shapes = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); + sequencing.beam(i).shapesWeight = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); + sequencing.beam(i).bixelIx = 1+offset:numOfRaysPerBeam+offset; + sequencing.beam(i).fluence = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); + end sequencing.w(1+offset:numOfRaysPerBeam+offset,1) = D_0(indInFluenceMx)/numOfLevels*calFac; - + offset = offset + numOfRaysPerBeam; end diff --git a/matRad/sequencing/matRad_siochiLeafSequencing.m b/matRad/sequencing/matRad_siochiLeafSequencing.m index 4240a6c9a..02841d5b0 100644 --- a/matRad/sequencing/matRad_siochiLeafSequencing.m +++ b/matRad/sequencing/matRad_siochiLeafSequencing.m @@ -70,12 +70,12 @@ numOfRaysPerBeam = stf(i).numOfRays; % get relevant weights for current beam - wOfCurrBeams = wUnsequenced(1+offset:numOfRaysPerBeam+offset);%REVIEW OFFSET + wOfCurrBeams = wUnsequenced(1+offset:numOfRaysPerBeam+offset).* ones(size(stf(i).ray,2),1);%REVIEW OFFSET - X = ones(numOfRaysPerBeam,1)*NaN; - Z = ones(numOfRaysPerBeam,1)*NaN; + X = ones(size(stf(i).ray,2),1)*NaN; %this way it also works with3dconformal + Z = ones(size(stf(i).ray,2),1)*NaN; - for j = 1:stf(i).numOfRays + for j = 1:size(stf(i).ray,2) X(j) = stf(i).ray(j).rayPos_bev(:,1); Z(j) = stf(i).ray(j).rayPos_bev(:,3); end @@ -137,27 +137,41 @@ D_k_MinX = min(D_k_X); D_k_MaxX = max(D_k_X); - %Decompose the port, do rod pushing - [tops, bases] = matRad_siochiDecomposePort(D_k,dimOfFluenceMxZ,dimOfFluenceMxX,D_k_MinZ,D_k_MaxZ,D_k_MinX,D_k_MaxX); - %Form segments with and without visualization - if visBool - [shapes,shapesWeight,k,D_k]=matRad_siochiConvertToSegments(shapes,shapesWeight,k,tops,bases,visBool,i,D_k,numOfLevels,seqFig,seqSubPlots); + if sum(wOfCurrBeams)>0 + %Decompose the port, do rod pushing + [tops, bases] = matRad_siochiDecomposePort(D_k,dimOfFluenceMxZ,dimOfFluenceMxX,D_k_MinZ,D_k_MaxZ,D_k_MinX,D_k_MaxX); + %Form segments with and without visualization + if visBool + [shapes,shapesWeight,k,D_k]=matRad_siochiConvertToSegments(shapes,shapesWeight,k,tops,bases,visBool,i,D_k,numOfLevels,seqFig,seqSubPlots); + else + [shapes,shapesWeight,k]=matRad_siochiConvertToSegments(shapes,shapesWeight,k,tops,bases); + end + + sequencing.beam(i).numOfShapes = k; + sequencing.beam(i).shapes = shapes(:,:,1:k); + sequencing.beam(i).shapesWeight = shapesWeight(1:k)/numOfLevels*calFac; + sequencing.beam(i).bixelIx = 1+offset:numOfRaysPerBeam+offset; + sequencing.beam(i).fluence = D_0; + sequencing.beam(i).sum = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); + + for j = 1:k + sequencing.beam(i).sum = sequencing.beam(i).sum+sequencing.beam(i).shapes(:,:,j)*sequencing.beam(i).shapesWeight(j); + end + else - [shapes,shapesWeight,k]=matRad_siochiConvertToSegments(shapes,shapesWeight,k,tops,bases); + sequencing.beam(i).numOfShapes = 1; + sequencing.beam(i).shapes = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); + sequencing.beam(i).shapesWeight = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); + sequencing.beam(i).bixelIx = 1+offset:numOfRaysPerBeam+offset; + sequencing.beam(i).fluence = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); + sequencing.beam(i).sum = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); end - sequencing.beam(i).numOfShapes = k; - sequencing.beam(i).shapes = shapes(:,:,1:k); - sequencing.beam(i).shapesWeight = shapesWeight(1:k)/numOfLevels*calFac; - sequencing.beam(i).bixelIx = 1+offset:numOfRaysPerBeam+offset; - sequencing.beam(i).fluence = D_0; - sequencing.beam(i).sum = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); - - for j = 1:k - sequencing.beam(i).sum = sequencing.beam(i).sum+sequencing.beam(i).shapes(:,:,j)*sequencing.beam(i).shapesWeight(j); + if numOfRaysPerBeam >1 + sequencing.w(1+offset:numOfRaysPerBeam+offset,1) = sequencing.beam(i).sum(indInFluenceMx); + else + sequencing.w(1+offset:numOfRaysPerBeam+offset,1) = wOfCurrBeams(1); end - sequencing.w(1+offset:numOfRaysPerBeam+offset,1) = sequencing.beam(i).sum(indInFluenceMx); - offset = offset + numOfRaysPerBeam; end diff --git a/matRad/sequencing/matRad_xiaLeafSequencing.m b/matRad/sequencing/matRad_xiaLeafSequencing.m index aaa893324..052c918c6 100644 --- a/matRad/sequencing/matRad_xiaLeafSequencing.m +++ b/matRad/sequencing/matRad_xiaLeafSequencing.m @@ -60,12 +60,12 @@ numOfRaysPerBeam = stf(i).numOfRays; % get relevant weights for current beam - wOfCurrBeams = resultGUI.w(1+offset:numOfRaysPerBeam+offset);%REVIEW OFFSET + wOfCurrBeams = resultGUI.w(1+offset:numOfRaysPerBeam+offset).* ones(size(stf(i).ray,2),1);%;%REVIEW OFFSET - X = ones(numOfRaysPerBeam,1)*NaN; - Z = ones(numOfRaysPerBeam,1)*NaN; - - for j=1:stf(i).numOfRays + X = ones(size(stf(i).ray,2),1)*NaN; %this way it also works with3dconformal + Z = ones(size(stf(i).ray,2),1)*NaN; + + for j=1:size(stf(i).ray,2) X(j) = stf(i).ray(j).rayPos_bev(:,1); Z(j) = stf(i).ray(j).rayPos_bev(:,3); end @@ -239,15 +239,26 @@ L_k = max(D_k(:)); % eq 5 end + + if sum(wOfCurrBeams)>0 + + sequencing.beam(i).numOfShapes = k; + sequencing.beam(i).shapes = shapes(:,:,1:k); + sequencing.beam(i).shapesWeight = shapesWeight(1:k)/numOfLevels*calFac; + sequencing.beam(i).bixelIx = 1+offset:numOfRaysPerBeam+offset; + sequencing.beam(i).fluence = D_0; + + else + sequencing.beam(i).numOfShapes = 1; + sequencing.beam(i).shapes = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); + sequencing.beam(i).shapesWeight = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); + sequencing.beam(i).bixelIx = 1+offset:numOfRaysPerBeam+offset; + sequencing.beam(i).fluence = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); + end - sequencing.beam(i).numOfShapes = k; - sequencing.beam(i).shapes = shapes(:,:,1:k); - sequencing.beam(i).shapesWeight = shapesWeight(1:k)/numOfLevels*calFac; - sequencing.beam(i).bixelIx = 1+offset:numOfRaysPerBeam+offset; - sequencing.beam(i).fluence = D_0; sequencing.w(1+offset:numOfRaysPerBeam+offset,1) = D_0(indInFluenceMx)/numOfLevels*calFac; - + offset = offset + numOfRaysPerBeam; end