matRad_engelLeafSequencing

 multileaf collimator leaf sequencing algorithm 
 for intensity modulated beams with multiple static segments accroding 
 to Engel et al. 2005 Discrete Applied Mathematics
 
 call
   resultGUI = matRad_engelLeafSequencing(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] http://www.sciencedirect.com/science/article/pii/S0166218X05001411

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

 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.

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

0001 function resultGUI = matRad_engelLeafSequencing(resultGUI,stf,dij,numOfLevels,visBool)
0002 % multileaf collimator leaf sequencing algorithm
0003 % for intensity modulated beams with multiple static segments accroding
0004 % to Engel et al. 2005 Discrete Applied Mathematics
0005 %
0006 % call
0007 %   resultGUI = matRad_engelLeafSequencing(resultGUI,stf,dij,numOfLevels,visBool)
0008 %
0009 % input
0010 %   resultGUI:          resultGUI struct to which the output data will be added, if
0011 %                       this field is empty resultGUI struct will be created
0012 %   stf:                matRad steering information struct
0013 %   dij:                matRad's dij matrix
0014 %   numOfLevels:        number of stratification levels
0015 %   visBool:            toggle on/off visualization (optional)
0016 %
0017 % output
0018 %   resultGUI:          matRad result struct containing the new dose cube
0019 %                       as well as the corresponding weights
0020 %
0021 % References
0022 %   [1] http://www.sciencedirect.com/science/article/pii/S0166218X05001411
0023 %
0024 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0025 %
0026 % Copyright 2015 the matRad development team.
0027 %
0028 % This file is part of the matRad project. It is subject to the license
0029 % terms in the LICENSE file found in the top-level directory of this
0030 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
0031 % of the matRad project, including this file, may be copied, modified,
0032 % propagated, or distributed except according to the terms contained in the
0033 % LICENSE file.
0034 %
0035 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0036 
0037 % if visBool not set toogle off visualization
0038 if nargin < 5
0039     visBool = 0;
0040 end
0041 
0042 numOfBeams = numel(stf);
0043 
0044 if visBool
0045     % create the sequencing figure
0046     sz = [800 1000]; % figure size
0047     screensize = get(0,'ScreenSize');
0048     xpos = ceil((screensize(3)-sz(2))/2); % center the figure on the screen horizontally
0049     ypos = ceil((screensize(4)-sz(1))/2); % center the figure on the screen vertically
0050     seqFig = figure('position',[xpos,ypos,sz(2),sz(1)]);     
0051 end
0052 
0053 offset = 0;
0054 
0055 for i = 1:numOfBeams
0056     
0057     numOfRaysPerBeam = stf(i).numOfRays; 
0058     
0059     % get relevant weights for current beam
0060     wOfCurrBeams = resultGUI.w(1+offset:numOfRaysPerBeam+offset);
0061     
0062     X = ones(numOfRaysPerBeam,1)*NaN;
0063     Z = ones(numOfRaysPerBeam,1)*NaN;
0064         
0065     for j=1:stf(i).numOfRays
0066       X(j) = stf(i).ray(j).rayPos_bev(:,1);
0067       Z(j) = stf(i).ray(j).rayPos_bev(:,3);
0068     end
0069         
0070     % sort bixels into matrix
0071     minX = min(X);
0072     maxX = max(X);
0073     minZ = min(Z);
0074     maxZ = max(Z);
0075     
0076     dimOfFluenceMxX = (maxX-minX)/stf(i).bixelWidth + 1;
0077     dimOfFluenceMxZ = (maxZ-minZ)/stf(i).bixelWidth + 1;
0078     
0079     %Create the fluence matrix.
0080     fluenceMx = zeros(dimOfFluenceMxZ,dimOfFluenceMxX);
0081     
0082     % Calculate X and Z position of every fluence's matrix spot
0083     % z axis = axis of leaf movement!
0084     xPos = (X-minX)/stf(i).bixelWidth+1;
0085     zPos = (Z-minZ)/stf(i).bixelWidth+1;
0086     
0087     % Make subscripts for fluence matrix
0088     indInFluenceMx = zPos + (xPos-1)*dimOfFluenceMxZ;
0089     
0090     %Save weights in fluence matrix.
0091     fluenceMx(indInFluenceMx) = wOfCurrBeams;
0092     
0093     % Stratification
0094     calFac = max(fluenceMx(:));
0095     D_k = round(fluenceMx/calFac*numOfLevels); 
0096     
0097     % Save the stratification in the initial intensity matrix D_0.
0098     D_0 = D_k;
0099     
0100     % container to remember generated shapes; allocate space for 10000 shapes
0101     shapes = NaN*ones(dimOfFluenceMxZ,dimOfFluenceMxX,10000);
0102   
0103     k = 0;
0104     
0105     if visBool
0106         clf(seqFig);
0107         colormap(seqFig,'jet');
0108         
0109         seqSubPlots(1) = subplot(2,2,1,'parent',seqFig);
0110         imagesc(D_k,'parent',seqSubPlots(1));
0111         set(seqSubPlots(1),'CLim',[0 numOfLevels],'YDir','normal');
0112         title(seqSubPlots(1),['Beam # ' num2str(i) ': max(D_0) = ' num2str(max(D_0(:))) ' - ' num2str(numel(unique(D_0))) ' intensity levels']);
0113         xlabel(seqSubPlots(1),'x - direction parallel to leaf motion ')
0114         ylabel(seqSubPlots(1),'z - direction perpendicular to leaf motion ')
0115         colorbar;
0116         drawnow
0117     end
0118     
0119     % start sequencer
0120     while max(D_k(:) > 0)
0121         
0122         %calculate the difference matrix diffMat
0123         diffMat = diff([zeros(size(D_k,1),1) D_k zeros(size(D_k,1),1)],[],2);
0124         
0125         %calculate complexities
0126         c = sum(max(0,diffMat),2); %TNMU-row-complexity
0127         com = max(c); %TNMU complexity
0128         g = com - c; %row complexity gap
0129         
0130         %initialize segment
0131         segment = zeros(size(D_k));
0132         
0133         k = k + 1;
0134         
0135         %Plot residual intensity matrix.
0136         if visBool
0137             seqSubPlots(2) = subplot(2,2,2,'parent',seqFig);
0138             imagesc(D_k,'parent',seqSubPlots(2));
0139             set(seqSubPlots(2),'CLim',[0 numOfLevels],'YDir','normal');
0140             title(seqSubPlots(2),['k = ' num2str(k)]);
0141             colorbar
0142             drawnow
0143         end
0144        
0145         
0146         %loop over all rows
0147         for j=1:size(D_0,1)
0148             
0149             %determine essential intervals
0150             data(j).left(1) = 0; %left interval limit, actual for an empty interval
0151             data(j).right(1) = 0; %right interal limit, actual for an empty interval
0152             data(j).v(1) = g(j);  %greatest number such that the inequalities (6) resp. (7) is satisfied with u=v
0153             data(j).w(1) = inf; %smallest number in the interval
0154             data(j).u(1) = data(j).v(1); %min(v,w)
0155             
0156             [~, pos, ~] = find(diffMat(j,:) > 0); % indices of all positive elements in the j. row of diffmat
0157             [~, neg, ~] = find(diffMat(j,:) < 0); % indices of all negative elements in the j. row of diffMat
0158             
0159             n=2;
0160             
0161             %loop over the positive elements in the j. row of diffmat ->
0162             %possible left interval limits
0163             for m=1:size(pos,2)
0164                 
0165                 %loop over the negative elements in the j. row of diffMat ->
0166                 %possible right interval limit
0167                 for l=1:size(neg,2)
0168                     
0169                     %take only intervals I=[l,r] with l<=r
0170                     if pos(m) <= neg(l)-1
0171                         
0172                         %set interval limits
0173                         data(j).left(n) = pos(m);
0174                         data(j).right(n) = neg(l)-1;
0175                         
0176                         %calculate v according to Lemma 8
0177                         if g(j) <= abs( diffMat(j,pos(m)) + diffMat(j,neg(l)) )
0178                             data(j).v(n) = min( diffMat(j,pos(m)), -diffMat(j,neg(l)) ) + g(j);
0179                         else
0180                             data(j).v(n) = ( diffMat(j, pos(m)) - diffMat(j, neg(l)) + g(j)) / 2;
0181                         end
0182                         
0183                         %calculate w and u according to equality (11) and
0184                         %(12)
0185                         data(j).w(n) = min(D_k(j,pos(m):(neg(l)-1)));
0186                         data(j).u(n) = min(data(j).v(n), data(j).w(n));
0187                         
0188                         n = n+1;
0189                     end
0190                 end
0191             end
0192             
0193             u(j) = max(data(j).u);
0194             
0195         end
0196         
0197         %calculate u_max from theorem 9
0198         d_k = min(u);
0199         
0200         %loop over all rows
0201         for j=1:size(D_0,1)
0202  
0203             %find all possible (and essential) intervals
0204             candidate = find(data(j).u >= d_k); 
0205             
0206             %calculate the potential of the possible intervals
0207             
0208             %initialize p as -Inf
0209             data(j).p(1:length(data(j).left)) = -Inf;
0210             
0211             %loop over all possible intervals
0212             for s=1:size(candidate,2)
0213                 
0214                 if (s==1 && data(j).left(candidate(s)) == 0)
0215                     data(j).p(candidate(1)) = 0;
0216                     
0217                     
0218                 else
0219                     %calculate p1 according to equality (17)
0220                     if (d_k == diffMat(j, data(j).left(candidate(s))) && d_k ~= D_k(j, data(j).left(candidate(s))))
0221                         p1 = 1;
0222                
0223                     else
0224                         p1 = 0;
0225                     
0226                     end
0227                     
0228                     %calculate p2 according to equalitiy (18)
0229                    % if data(j).right(candidate(s)) < size(D_0, 2)
0230                         
0231                         if (d_k == -diffMat(j, data(j).right(candidate(s))+1) && d_k ~= D_k(j, data(j).right(candidate(s))))
0232                             p2 = 1;
0233                         else
0234                             p2 = 0;
0235                         end
0236                         
0237 %                     else
0238 %
0239 %                         if d_k == -diffMat(j, data(j).right(candidate(s))+1)
0240 %                             p2 = 1;
0241 %                         else
0242 %                             p2 = 0;
0243 %                         end
0244 %
0245 %                     end
0246                     
0247                     %calculate p3 according to equality (19)
0248                     p3 = size(find(D_k(j, data(j).left(candidate(s)):data(j).right(candidate(s))) == d_k),2);
0249                     
0250                     data(j).p(candidate(s)) = p1 + p2+ p3;
0251                     
0252                 end
0253                 
0254             end
0255                 
0256                 %determinate intervals with maximum potential
0257                 maxPot = find(data(j).p == max(data(j).p));
0258                 
0259                 %if several intervals have maximum potential, select
0260                 %the interval which has maximum length
0261                 if size(maxPot,2) > 1
0262 
0263                     for t=1:size(maxPot,2)
0264                         if t==1 && data(j).left(maxPot(t)) == 0
0265                             data(j).l(1) = 0;
0266                         else
0267                             data(j).l(maxPot(t)) = data(j).right(maxPot(t)) - data(j).left(maxPot(t)) + 1;
0268                         end
0269                     end
0270                     
0271                     %data(j).l(maxPot) = data(j).right(maxPot) - data(j).left(maxPot) + 1;
0272                     
0273                     maxLength = find(data(j).l == max(data(j).l));
0274                      
0275                     %left and right interval limits of the selected
0276                     %interval
0277                     leftIntLimit(j) = data(j).left(maxLength(1));
0278                     rightIntLimit(j) = data(j).right(maxLength(1));
0279                     
0280  
0281                 else
0282                     
0283                     %left and right interval limits of the selected
0284                     %interval
0285                     leftIntLimit(j) = data(j).left(maxPot);
0286                     rightIntLimit(j) = data(j).right(maxPot);
0287                     
0288                     
0289                 end
0290                 
0291                 %create segment associated by the selected interval
0292                 if leftIntLimit(j) ~= 0 
0293                     
0294                     segment(j,leftIntLimit(j):rightIntLimit(j)) = 1;
0295  
0296                 end
0297 
0298         end
0299         
0300         %write the segment in shape_k
0301         shape_k = segment;
0302 
0303         %show the leaf positions
0304         if visBool
0305             seqSubPlots(4) = subplot(2,2,3.5,'parent',seqFig);
0306             imagesc(shape_k,'parent',seqSubPlots(4));
0307             hold(seqSubPlots(4),'on');
0308             set(seqSubPlots(4),'YDir','normal')
0309             xlabel(seqSubPlots(4),'x - direction parallel to leaf motion ')
0310             ylabel(seqSubPlots(4),'z - direction perpendicular to leaf motion ')
0311             title(seqSubPlots(4),['beam # ' num2str(i) ' shape # ' num2str(k) ' d_k = ' num2str(d_k)]);
0312             for j = 1:dimOfFluenceMxZ
0313                 leftLeafIx = find(shape_k(j,:)>0,1,'first');
0314                 rightLeafIx = find(shape_k(j,:)>0,1,'last');
0315                 if leftLeafIx > 1
0316                     plot(seqSubPlots(4),[.5 leftLeafIx-.5],j-[.5 .5] ,'w','LineWidth',2)
0317                     plot(seqSubPlots(4),[.5 leftLeafIx-.5],j+[.5 .5] ,'w','LineWidth',2)
0318                     plot(seqSubPlots(4),[ leftLeafIx-.5 leftLeafIx-.5],j+[.5 -.5] ,'w','LineWidth',2)
0319                 end
0320                 if rightLeafIx<dimOfFluenceMxX
0321                     plot(seqSubPlots(4),[dimOfFluenceMxX+.5 rightLeafIx+.5],j-[.5 .5] ,'w','LineWidth',2)
0322                     plot(seqSubPlots(4),[dimOfFluenceMxX+.5 rightLeafIx+.5],j+[.5 .5] ,'w','LineWidth',2)
0323                     plot(seqSubPlots(4),[ rightLeafIx+.5 rightLeafIx+.5],j+[.5 -.5] ,'w','LineWidth',2)
0324                 end
0325                 if isempty(rightLeafIx) && isempty (leftLeafIx)
0326                     plot(seqSubPlots(4),[dimOfFluenceMxX+.5 .5],j-[.5 .5] ,'w','LineWidth',2)
0327                     plot(seqSubPlots(4),[dimOfFluenceMxX+.5 .5],j+[.5 .5] ,'w','LineWidth',2)
0328                     plot(seqSubPlots(4),.5*dimOfFluenceMxX*[1 1]+[0.5],j+[.5 -.5] ,'w','LineWidth',2)
0329                 end
0330             end
0331             
0332             axis tight
0333             drawnow
0334             pause(1);
0335         end
0336     
0337 
0338         %save shape_k in container
0339         shapes(:,:,k) =  shape_k;
0340         
0341         %save the calculated MU
0342         shapesWeight(k) = d_k; 
0343         
0344         %calculate  new matrix, the  diference matrix and complexities
0345         D_k = D_k - d_k*shape_k; 
0346         
0347         %delete variables
0348         clear data;
0349         clear segment;
0350         clear u;
0351         clear leftIntLimit;
0352         clear rightIntLimit;
0353        
0354     end
0355     
0356     sequencing.beam(i).numOfShapes  = k;
0357     sequencing.beam(i).shapes       = shapes(:,:,1:k);
0358     sequencing.beam(i).shapesWeight = shapesWeight(1:k)/numOfLevels*calFac;
0359     sequencing.beam(i).bixelIx      = 1+offset:numOfRaysPerBeam+offset;
0360     sequencing.beam(i).fluence      = D_0;
0361     
0362     sequencing.w(1+offset:numOfRaysPerBeam+offset,1) = D_0(indInFluenceMx)/numOfLevels*calFac;
0363 
0364     offset = offset + numOfRaysPerBeam;
0365 
0366 end
0367 
0368 resultGUI.w          = sequencing.w;
0369 resultGUI.wSequenced = sequencing.w;
0370 
0371 resultGUI.sequencing   = sequencing;
0372 resultGUI.apertureInfo = matRad_sequencing2ApertureInfo(sequencing,stf);
0373 
0374 doseSequencedDoseGrid = reshape(dij.physicalDose{1} * sequencing.w,dij.doseGrid.dimensions);
0375 % interpolate to ct grid for visualiation & analysis
0376 resultGUI.physicalDose = matRad_interp3(dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z, ...
0377                                         doseSequencedDoseGrid, ...
0378                                         dij.ctGrid.x,dij.ctGrid.y',dij.ctGrid.z);
0379 
0380 % if weights exists from an former DAO remove it
0381 if isfield(resultGUI,'wDao')
0382     resultGUI = rmfield(resultGUI,'wDao');
0383 end
0384 
0385 end
0386