matRad_siochiLeafSequencing

 multileaf collimator leaf sequencing algorithm 
 for intensity modulated beams with multiple static segments according to 
 Siochi (1999)International Journal of Radiation Oncology * Biology * Physics,
 originally implemented in PLUNC (https://sites.google.com/site/planunc/)

 Implemented in matRad by Eric Christiansen, Emily Heath, and Tong Xu

 call
   resultGUI =
   matRad_siochiLeafSequencing(resultGUI,stf,dij,numOfLevels)
   resultGUI =
   matRad_siochiLeafSequencing(resultGUI,stf,dij,numOfLevels,visBool)

 input
   resultGUI:          resultGUI struct to which the output data will be
                       added, if this field is empty resultGUI struct will
                       be created
   stf:                matRad steering information struct
   dij:                matRad's dij matrix
   numOfLevels:        number of stratification levels
   visBool:            toggle on/off visualization (optional)

 output
   resultGUI:          matRad result struct containing the new dose cube
                       as well as the corresponding weights

 References
   [1] https://www.ncbi.nlm.nih.gov/pubmed/10078655

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 Copyright 2015 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.

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 if visBool not set toogle off visualization

0001 function resultGUI = matRad_siochiLeafSequencing(resultGUI,stf,dij,numOfLevels,visBool)
0002 % multileaf collimator leaf sequencing algorithm
0003 % for intensity modulated beams with multiple static segments according to
0004 % Siochi (1999)International Journal of Radiation Oncology * Biology * Physics,
0005 % originally implemented in PLUNC (https://sites.google.com/site/planunc/)
0006 %
0007 % Implemented in matRad by Eric Christiansen, Emily Heath, and Tong Xu
0008 %
0009 % call
0010 %   resultGUI =
0011 %   matRad_siochiLeafSequencing(resultGUI,stf,dij,numOfLevels)
0012 %   resultGUI =
0013 %   matRad_siochiLeafSequencing(resultGUI,stf,dij,numOfLevels,visBool)
0014 %
0015 % input
0016 %   resultGUI:          resultGUI struct to which the output data will be
0017 %                       added, if this field is empty resultGUI struct will
0018 %                       be created
0019 %   stf:                matRad steering information struct
0020 %   dij:                matRad's dij matrix
0021 %   numOfLevels:        number of stratification levels
0022 %   visBool:            toggle on/off visualization (optional)
0023 %
0024 % output
0025 %   resultGUI:          matRad result struct containing the new dose cube
0026 %                       as well as the corresponding weights
0027 %
0028 % References
0029 %   [1] https://www.ncbi.nlm.nih.gov/pubmed/10078655
0030 %
0031 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0032 %
0033 % Copyright 2015 the matRad development team.
0034 %
0035 % This file is part of the matRad project. It is subject to the license
0036 % terms in the LICENSE file found in the top-level directory of this
0037 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
0038 % of the matRad project, including this file, may be copied, modified,
0039 % propagated, or distributed except according to the terms contained in the
0040 % LICENSE file.
0041 %
0042 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0043 % if visBool not set toogle off visualization
0044 if nargin < 5
0045     visBool = 0;
0046 end
0047 
0048 numOfBeams = numel(stf);
0049 
0050 if visBool
0051     % create the sequencing figure
0052     sz = [800 1000]; % figure size
0053     screensize = get(0,'ScreenSize');
0054     xpos = ceil((screensize(3)-sz(2))/2); % center the figure on the screen horizontally
0055     ypos = ceil((screensize(4)-sz(1))/2); % center the figure on the screen vertically
0056     seqFig = figure('position',[xpos,ypos,sz(2),sz(1)]);     
0057 end
0058 
0059 offset = 0;
0060 
0061 for i = 1:numOfBeams
0062     
0063     numOfRaysPerBeam = stf(i).numOfRays;
0064     
0065     % get relevant weights for current beam
0066     wOfCurrBeams = resultGUI.wUnsequenced(1+offset:numOfRaysPerBeam+offset);%REVIEW OFFSET
0067     
0068     X = ones(numOfRaysPerBeam,1)*NaN;
0069     Z = ones(numOfRaysPerBeam,1)*NaN;
0070     
0071     for j = 1:stf(i).numOfRays
0072         X(j) = stf(i).ray(j).rayPos_bev(:,1);
0073         Z(j) = stf(i).ray(j).rayPos_bev(:,3);
0074     end
0075     
0076     % sort bixels into matrix
0077     minX = min(X);
0078     maxX = max(X);
0079     minZ = min(Z);
0080     maxZ = max(Z);
0081     
0082     dimOfFluenceMxX = (maxX-minX)/stf(i).bixelWidth + 1;
0083     dimOfFluenceMxZ = (maxZ-minZ)/stf(i).bixelWidth + 1;
0084     
0085     %Create the fluence matrix.
0086     fluenceMx = zeros(dimOfFluenceMxZ,dimOfFluenceMxX);
0087     
0088     % Calculate X and Z position of every fluence's matrix spot z axis =
0089     % axis of leaf movement!
0090     xPos = (X-minX)/stf(i).bixelWidth+1;
0091     zPos = (Z-minZ)/stf(i).bixelWidth+1;
0092     
0093     % Make subscripts for fluence matrix
0094     indInFluenceMx = zPos + (xPos-1)*dimOfFluenceMxZ;
0095     
0096     %Save weights in fluence matrix.
0097     fluenceMx(indInFluenceMx) = wOfCurrBeams;
0098     
0099     % Stratification
0100     calFac = max(fluenceMx(:));
0101     D_k = round(fluenceMx/calFac*numOfLevels);
0102     
0103     % Save the stratification in the initial intensity matrix D_0.
0104     D_0 = D_k;
0105     
0106     % container to remember generated shapes; allocate space for 10000
0107     % shapes
0108     shapes = NaN*ones(dimOfFluenceMxZ,dimOfFluenceMxX,10000);
0109     shapesWeight = zeros(10000,1);
0110     k = 0;
0111     
0112     if visBool
0113         clf(seqFig);
0114         colormap(seqFig,'jet');
0115         
0116         seqSubPlots(1) = subplot(2,2,1,'parent',seqFig);
0117         imagesc(D_k,'parent',seqSubPlots(1));
0118         set(seqSubPlots(1),'CLim',[0 numOfLevels],'YDir','normal');
0119         title(seqSubPlots(1),['Beam # ' num2str(i) ': max(D_0) = ' num2str(max(D_0(:))) ' - ' num2str(numel(unique(D_0))) ' intensity levels']);
0120         xlabel(seqSubPlots(1),'x - direction parallel to leaf motion ')
0121         ylabel(seqSubPlots(1),'z - direction perpendicular to leaf motion ')
0122         colorbar;
0123         drawnow
0124     end
0125     
0126     D_k_nonZero = (D_k~=0);
0127     [D_k_Z, D_k_X] = ind2sub([dimOfFluenceMxZ,dimOfFluenceMxX],find(D_k_nonZero));
0128     D_k_MinZ = min(D_k_Z);
0129     D_k_MaxZ = max(D_k_Z);
0130     D_k_MinX = min(D_k_X);
0131     D_k_MaxX = max(D_k_X);
0132     
0133     %Decompose the port, do rod pushing
0134     [tops, bases] = matRad_siochiDecomposePort(D_k,dimOfFluenceMxZ,dimOfFluenceMxX,D_k_MinZ,D_k_MaxZ,D_k_MinX,D_k_MaxX);
0135     %Form segments with and without visualization
0136     if visBool
0137         [shapes,shapesWeight,k,D_k]=matRad_siochiConvertToSegments(shapes,shapesWeight,k,tops,bases,visBool,i,D_k,numOfLevels,seqFig,seqSubPlots);
0138     else
0139         [shapes,shapesWeight,k]=matRad_siochiConvertToSegments(shapes,shapesWeight,k,tops,bases);
0140     end
0141     
0142     sequencing.beam(i).numOfShapes  = k;
0143     sequencing.beam(i).shapes       = shapes(:,:,1:k);
0144     sequencing.beam(i).shapesWeight = shapesWeight(1:k)/numOfLevels*calFac;
0145     sequencing.beam(i).bixelIx      = 1+offset:numOfRaysPerBeam+offset;
0146     sequencing.beam(i).fluence      = D_0;
0147     sequencing.beam(i).sum          = zeros(dimOfFluenceMxZ,dimOfFluenceMxX);
0148     
0149     for j = 1:k
0150         sequencing.beam(i).sum = sequencing.beam(i).sum+sequencing.beam(i).shapes(:,:,j)*sequencing.beam(i).shapesWeight(j);
0151     end
0152     sequencing.w(1+offset:numOfRaysPerBeam+offset,1) = sequencing.beam(i).sum(indInFluenceMx);
0153     
0154     offset = offset + numOfRaysPerBeam;
0155 
0156 end
0157 
0158 resultGUI.w          = sequencing.w;
0159 resultGUI.wSequenced = sequencing.w;
0160 
0161 resultGUI.sequencing   = sequencing;
0162 resultGUI.apertureInfo = matRad_sequencing2ApertureInfo(sequencing,stf);
0163 
0164 doseSequencedDoseGrid = reshape(dij.physicalDose{1} * sequencing.w,dij.doseGrid.dimensions);
0165 % interpolate to ct grid for visualiation & analysis
0166 resultGUI.physicalDose = matRad_interp3(dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z, ...
0167                                         doseSequencedDoseGrid, ...
0168                                         dij.ctGrid.x,dij.ctGrid.y',dij.ctGrid.z);
0169 
0170 % if weights exists from an former DAO remove it
0171 if isfield(resultGUI,'wDao')
0172     resultGUI = rmfield(resultGUI,'wDao');
0173 end
0174 
0175 end
0176 
0177 function [tops, bases] = matRad_siochiDecomposePort(map,dimZ,dimX,minZ,maxZ,minX,maxX)
0178 %Returns tops and bases of a fluence matrix "map" for Siochi leaf
0179 %sequencing algorithm (rod pushing part).  Accounts for collisions and
0180 %tongue and groove (Tng) effects.
0181 
0182 tops = zeros(dimZ, dimX);
0183 bases = zeros(dimZ, dimX);
0184 
0185 for i = minX:maxX
0186     maxTop = -1;
0187     TnG = 1;
0188     for j = minZ:maxZ
0189         if i == minX
0190             bases(j,i) = 1;
0191             tops(j,i) = bases(j,i)+map(j,i)-1;
0192         else %assign trial base positions
0193             if map(j,i) >= map(j,i-1) %current rod >= previous, match the bases
0194                 bases(j,i) = bases(j,i-1);
0195                 tops(j,i) = bases(j,i)+map(j,i)-1;
0196             else %current rod <previous
0197                 if map(j,i) == 0 %rod length=0, put in in next slab after top of previous
0198                     bases(j,i) = tops(j,i-1)+1;
0199                     tops(j,i) = bases(j,i)-1;
0200                 else %rod length~=0, match tops
0201                     tops(j,i) = tops(j,i-1);
0202                     bases(j,i) = tops(j,i)-map(j,i)+1;
0203                 end
0204             end
0205         end
0206         %determine which rod has the highest top in column
0207         if tops(j,i) > maxTop
0208             maxTop = tops(j,i);
0209             maxRow = j;
0210         end
0211     end
0212     
0213     %Correct for collision and tongue and groove error
0214     while(TnG)
0215         %go from maxRow down checking for TnG.  This occurs when a shorter
0216         %rod is "peeking over" a longer one in the direction transverse to
0217         %the leaf motion.  To fix this, match either the tops or bases of
0218         %the rods.
0219         for j = (maxRow-1):-1:minZ
0220             if map(j,i) < map(j+1,i)
0221                 if tops(j,i) > tops(j+1,i)
0222                     tops(j+1,i) = tops(j,i);
0223                     bases(j+1,i) = tops(j+1,i)-map(j+1,i)+1;
0224                 elseif bases(j,i) < bases(j+1,i)
0225                     bases(j,i) = bases(j+1,i);
0226                     tops(j,i) = bases(j,i)+map(j,i)-1;
0227                 end
0228             else
0229                 if tops(j,i) < tops(j+1,i)
0230                     tops(j,i) = tops(j+1,i);
0231                     bases(j,i) = tops(j,i)-map(j,i)+1;
0232                 elseif bases(j,i) > bases(j+1,i)
0233                     bases(j+1,i) = bases(j,i);
0234                     tops(j+1,i) = bases(j+1,i)+map(j+1,i)-1;
0235                 end
0236             end
0237         end
0238         %go from maxRow up checking for TnG
0239         for j = (maxRow+1):maxZ
0240             if map(j,i) < map(j-1,i)
0241                 if tops(j,i) > tops(j-1,i)
0242                     tops(j-1,i) = tops(j,i);
0243                     bases(j-1,i) = tops(j-1,i)-map(j-1,i)+1;
0244                 elseif bases(j,i) < bases(j-1,i)
0245                     bases(j,i) = bases(j-1,i);
0246                     tops(j,i) = bases(j,i)+map(j,i)-1;
0247                 end
0248             else
0249                 if tops(j,i) < tops(j-1,i)
0250                     tops(j,i) = tops(j-1,i);
0251                     bases(j,i) = tops(j,i)-map(j,i)+1;
0252                 elseif bases(j,i) > bases(j-1,i)
0253                     bases(j-1,i) = bases(j,i);
0254                     tops(j-1,i) = bases(j-1,i)+map(j-1,i)-1;
0255                 end
0256             end
0257         end
0258         %now check if all TnG conditions have been removed
0259         TnG = 0;
0260         for j = (minZ+1):maxZ
0261             if map(j,i) < map(j-1,i);
0262                 if tops(j,i) > tops(j-1,i)
0263                     TnG = 1;
0264                 elseif bases(j,i) < bases(j-1,i)
0265                     TnG = 1;
0266                 end
0267             else
0268                 if tops(j,i) < tops(j-1,i)
0269                     TnG = 1;
0270                 elseif bases(j,i) > bases(j-1,i)
0271                     TnG = 1;
0272                 end
0273             end
0274         end
0275     end
0276 end
0277 
0278 end
0279 
0280 function [shapes,shapesWeight,k,D_k] = matRad_siochiConvertToSegments(shapes,shapesWeight,k,tops,bases,visBool,i,D_k,numOfLevels,seqFig,seqSubPlots)
0281 %Convert tops and bases to shape matrices.  These are taken as to be the
0282 %shapes of uniform level/elevation after the rods are pushed.
0283 if nargin < 6
0284     visBool = 0;
0285 end
0286 
0287 
0288 levels = max(tops(:));
0289 
0290 for level = 1:levels
0291     %check if slab is new
0292     if matRad_siochiDifferentSlab(tops,bases,level)
0293         k = k+1; %increment number of unique slabs
0294         shape_k = (bases <= level).*(level <= tops); %shape of current slab
0295         shapes(:,:,k) = shape_k;
0296     end
0297     shapesWeight(k) = shapesWeight(k)+1; %if slab is not unique, this increments weight again
0298     
0299     if visBool
0300         %show the leaf positions
0301         [dimZ,dimX] = size(tops);
0302         seqSubPlots(4) = subplot(2,2,3.5,'parent',seqFig);
0303         imagesc(shape_k,'parent',seqSubPlots(4));
0304         hold(seqSubPlots(4),'on');
0305         set(seqSubPlots(4),'YDir','normal')
0306         xlabel(seqSubPlots(4),'x - direction parallel to leaf motion ')
0307         ylabel(seqSubPlots(4),'z - direction perpendicular to leaf motion ')
0308         title(seqSubPlots(4),['beam # ' num2str(i) ' shape # ' num2str(k) ' d_k = ' num2str(shapesWeight(k))]);
0309         for j = 1:dimZ
0310             leftLeafIx = find(shape_k(j,:)>0,1,'first');
0311             rightLeafIx = find(shape_k(j,:)>0,1,'last');
0312             if leftLeafIx > 1
0313                 plot(seqSubPlots(4),[.5 leftLeafIx-.5],j-[.5 .5] ,'w','LineWidth',2)
0314                 plot(seqSubPlots(4),[.5 leftLeafIx-.5],j+[.5 .5] ,'w','LineWidth',2)
0315                 plot(seqSubPlots(4),[ leftLeafIx-.5 leftLeafIx-.5],j+[.5 -.5] ,'w','LineWidth',2)
0316             end
0317             if rightLeafIx<dimX
0318                 plot(seqSubPlots(4),[dimX+.5 rightLeafIx+.5],j-[.5 .5] ,'w','LineWidth',2)
0319                 plot(seqSubPlots(4),[dimX+.5 rightLeafIx+.5],j+[.5 .5] ,'w','LineWidth',2)
0320                 plot(seqSubPlots(4),[ rightLeafIx+.5 rightLeafIx+.5],j+[.5 -.5] ,'w','LineWidth',2)
0321             end
0322             if isempty(rightLeafIx) && isempty (leftLeafIx)
0323                 plot(seqSubPlots(4),[dimX+.5 .5],j-[.5 .5] ,'w','LineWidth',2)
0324                 plot(seqSubPlots(4),[dimX+.5 .5],j+[.5 .5] ,'w','LineWidth',2)
0325                 plot(seqSubPlots(4),.5*dimX*[1 1]+[0.5],j+[.5 -.5] ,'w','LineWidth',2)
0326             end
0327         end
0328         pause(1);
0329         
0330         %Plot residual intensity matrix.
0331         D_k = D_k-shape_k; %residual intensity matrix for visualization
0332         seqSubPlots(2) = subplot(2,2,2,'parent',seqFig);
0333         imagesc(D_k,'parent',seqSubPlots(2));
0334         set(seqSubPlots(2),'CLim',[0 numOfLevels],'YDir','normal');
0335         title(seqSubPlots(2),['k = ' num2str(k)]);
0336         colorbar
0337         drawnow
0338         
0339         axis tight
0340         drawnow
0341     end
0342 end
0343 
0344 end
0345 
0346 function diffSlab = matRad_siochiDifferentSlab(tops,bases,level)
0347 %Returns 1 if slab level is different than slab level-1 0 otherwise
0348 
0349 if level == 1 %first slab is automatically different
0350     diffSlab = 1;
0351 else
0352     shapeLevel = (bases <= level).*(level <= tops); %shape of slab with current level
0353     shapeLevel_1 = (bases <= level-1).*(level-1 <= tops); %shape of slab with previous level
0354     diffSlab = ~isequal(shapeLevel,shapeLevel_1); %tests if slabs are equal; isequaln was not giving correct results
0355 end
0356 
0357 end
0358