Reworked 3d conformal and sequencing

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

matRad/optimization/matRad_collapseDij.m

@@ -28,21 +28,26 @@
 %
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-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;
 
 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
 

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

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

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

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