matRad_xiaLeafSequencing

Purpose ^

multileaf collimator leaf sequencing algorithm

Synopsis ^

function resultGUI = matRad_xiaLeafSequencing(resultGUI,stf,dij,numOfLevels,visBool)

Description ^

 multileaf collimator leaf sequencing algorithm 
 for intensity modulated beams with multiple static segments according to 
 Xia et al. (1998) Medical Physics
 
 call
   resultGUI = matRad_xiaLeafSequencing(resultGUI,stf,dij,numOfLevels)
   resultGUI = matRad_xiaLeafSequencing(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://online.medphys.org/resource/1/mphya6/v25/i8/p1424_s1

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

 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.

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

Cross-reference information ^

This function calls: This function is called by:

Source code ^

0001 function resultGUI = matRad_xiaLeafSequencing(resultGUI,stf,dij,numOfLevels,visBool)
0002 % multileaf collimator leaf sequencing algorithm
0003 % for intensity modulated beams with multiple static segments according to
0004 % Xia et al. (1998) Medical Physics
0005 %
0006 % call
0007 %   resultGUI = matRad_xiaLeafSequencing(resultGUI,stf,dij,numOfLevels)
0008 %   resultGUI = matRad_xiaLeafSequencing(resultGUI,stf,dij,numOfLevels,visBool)
0009 %
0010 % input
0011 %   resultGUI:          resultGUI struct to which the output data will be added, if
0012 %                       this field is empty resultGUI struct will be created
0013 %   stf:                matRad steering information struct
0014 %   dij:                matRad's dij matrix
0015 %   numOfLevels:        number of stratification levels
0016 %   visBool:            toggle on/off visualization (optional)
0017 %
0018 % output
0019 %   resultGUI:          matRad result struct containing the new dose cube as well as
0020 %                       the corresponding weights
0021 %
0022 % References
0023 %   [1] http://online.medphys.org/resource/1/mphya6/v25/i8/p1424_s1
0024 %
0025 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0026 %
0027 % Copyright 2015 the matRad development team.
0028 %
0029 % This file is part of the matRad project. It is subject to the license
0030 % terms in the LICENSE file found in the top-level directory of this
0031 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
0032 % of the matRad project, including this file, may be copied, modified,
0033 % propagated, or distributed except according to the terms contained in the
0034 % LICENSE file.
0035 %
0036 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0037 
0038 % if visBool not set toogle off visualization
0039 if nargin < 5
0040     visBool = 0;
0041 end
0042 
0043 mode = 'rl'; % sliding window (sw) or reducing level (rl)
0044 
0045 numOfBeams = numel(stf);
0046 
0047 if visBool
0048     % create the sequencing figure
0049     sz = [800 1000]; % figure size
0050     screensize = get(0,'ScreenSize');
0051     xpos = ceil((screensize(3)-sz(2))/2); % center the figure on the screen horizontally
0052     ypos = ceil((screensize(4)-sz(1))/2); % center the figure on the screen vertically
0053     seqFig = figure('position',[xpos,ypos,sz(2),sz(1)]);     
0054 end
0055 
0056 offset = 0;
0057 
0058 for i = 1:numOfBeams
0059     
0060     numOfRaysPerBeam = stf(i).numOfRays;
0061     
0062     % get relevant weights for current beam
0063     wOfCurrBeams = resultGUI.w(1+offset:numOfRaysPerBeam+offset);%REVIEW OFFSET
0064     
0065     X = ones(numOfRaysPerBeam,1)*NaN;
0066     Z = ones(numOfRaysPerBeam,1)*NaN;
0067         
0068     for j=1:stf(i).numOfRays
0069       X(j) = stf(i).ray(j).rayPos_bev(:,1);
0070       Z(j) = stf(i).ray(j).rayPos_bev(:,3);
0071     end
0072         
0073     % sort bixels into matrix
0074     minX = min(X);
0075     maxX = max(X);
0076     minZ = min(Z);
0077     maxZ = max(Z);
0078     
0079     dimOfFluenceMxX = (maxX-minX)/stf(i).bixelWidth + 1;
0080     dimOfFluenceMxZ = (maxZ-minZ)/stf(i).bixelWidth + 1;
0081     
0082     %Create the fluence matrix.
0083     fluenceMx = zeros(dimOfFluenceMxZ,dimOfFluenceMxX);
0084     
0085     % Calculate X and Z position of every fluence's matrix spot
0086     % z axis = axis of leaf movement!
0087     xPos = (X-minX)/stf(i).bixelWidth+1;
0088     zPos = (Z-minZ)/stf(i).bixelWidth+1;
0089     
0090     % Make subscripts for fluence matrix
0091     indInFluenceMx = zPos + (xPos-1)*dimOfFluenceMxZ;
0092     
0093     %Save weights in fluence matrix.
0094     fluenceMx(indInFluenceMx) = wOfCurrBeams;
0095     
0096     % Stratification
0097     calFac = max(fluenceMx(:));
0098     D_k = round(fluenceMx/calFac*numOfLevels); 
0099     
0100     % Save the stratification in the initial intensity matrix D_0.
0101     D_0 = D_k;
0102     
0103     % Save the maximun intensity (Equation 5)
0104     L_k = max(D_k(:));
0105     
0106     % Save the maximun initial intensity matrix value in L_0.
0107     L_0 = L_k;
0108     
0109     % Set k=0, this variable is used for residuals intensity matrices D_k.
0110     k = 0;
0111     
0112     % container to remember generated shapes; allocate space for 10000 shapes
0113     shapes = NaN*ones(dimOfFluenceMxZ,dimOfFluenceMxX,10000);
0114     
0115     if visBool
0116         clf(seqFig);
0117         colormap(seqFig,'jet');
0118         
0119         seqSubPlots(1) = subplot(2,2,1,'parent',seqFig);
0120         imagesc(D_k,'parent',seqSubPlots(1));
0121         set(seqSubPlots(1),'CLim',[0 L_0],'YDir','normal');
0122         title(seqSubPlots(1),['Beam # ' num2str(i) ': L_0 = ' num2str(L_0) ' - ' num2str(numel(unique(D_0))) ' intensity levels'])
0123         xlabel(seqSubPlots(1),'x - direction parallel to leaf motion ')
0124         ylabel(seqSubPlots(1),'z - direction perpendicular to leaf motion ')
0125         colorbar;
0126         drawnow
0127     end
0128     
0129     % start sequencer
0130     while L_k > 0
0131         
0132         k = k + 1;
0133         
0134         %Plot residual intensity matrix.
0135         if visBool
0136             seqSubPlots(2) = subplot(2,2,2,'parent',seqFig);
0137             imagesc(D_k,'parent',seqSubPlots(2));
0138             set(seqSubPlots(2),'CLim',[0 L_0],'YDir','normal');
0139             title(seqSubPlots(2),['k = ' num2str(k) ' - ' num2str(numel(unique(D_k))) ' intensity levels remaining...']);
0140             xlabel(seqSubPlots(2),'x - direction parallel to leaf motion ');
0141             ylabel(seqSubPlots(2),'z - direction perpendicular to leaf motion ');
0142             colorbar
0143             drawnow
0144         end
0145         
0146         %Rounded off integer. Equation 7.
0147         m = floor(log2(L_k));
0148         
0149         % Convert m=1 if is less than 1. This happens when L_k belong to ]0,2[
0150         if m < 1
0151             m = 1;
0152         end
0153         
0154         %Calculate the delivery intensity unit. Equation 6.
0155         d_k = floor(2^(m-1));
0156         
0157         % Opening matrix.
0158         openingMx = D_k >= d_k;
0159         
0160         % Plot opening matrix.
0161         if visBool
0162             seqSubPlots(3) = subplot(2,2,3,'parent',seqFig);
0163             imagesc(openingMx,'parent',seqSubPlots(3));
0164             set(seqSubPlots(3),'YDir','normal')
0165             xlabel(seqSubPlots(3),'x - direction parallel to leaf motion ')
0166             ylabel(seqSubPlots(3),'z - direction perpendicular to leaf motion ')
0167             title(seqSubPlots(3),'Opening matrix');
0168             drawnow
0169         end
0170         
0171         if strcmp(mode,'sw') % sliding window technique!
0172             for j = 1:dimOfFluenceMxZ
0173                 openIx = find(openingMx(j,:) == 1,1,'first');
0174                 if ~isempty(openIx)
0175                     closeIx = find(openingMx(j,openIx+1:end) == 0,1,'first');
0176                     if ~isempty(closeIx)
0177                         openingMx(j,openIx+closeIx:end) = 0;
0178                     end
0179                 end
0180                 
0181             end
0182         elseif strcmp(mode,'rl') % reducing levels technique!
0183             for j = 1:dimOfFluenceMxZ
0184                 [maxVal,maxIx] = max(openingMx(j,:) .* D_k(j,:));
0185                 if maxVal > 0
0186                     closeIx = maxIx + find(openingMx(j,maxIx+1:end) == 0,1,'first');
0187                     if ~isempty(closeIx)
0188                         openingMx(j,closeIx:end) = 0;
0189                     end
0190                     openIx = find(openingMx(j,1:maxIx-1) == 0,1,'last');
0191                     if ~isempty(openIx)
0192                         openingMx(j,1:openIx) = 0;
0193                     end
0194                 end
0195                                 
0196             end
0197             
0198         end
0199         
0200         shape_k       = openingMx * d_k;
0201                 
0202                 if visBool
0203                     seqSubPlots(4) = subplot(2,2,4,'parent',seqFig);
0204                     imagesc(shape_k,'parent',seqSubPlots(4));
0205                     set(seqSubPlots(4),'YDir','normal')
0206                     hold(seqSubPlots(4),'on');
0207                     xlabel(seqSubPlots(4),'x - direction parallel to leaf motion ')
0208                     ylabel(seqSubPlots(4),'z - direction perpendicular to leaf motion ')
0209                     title(seqSubPlots(4),['beam # ' num2str(i) ' shape # ' num2str(k) ' d_k = ' num2str(d_k)]);
0210                     for j = 1:dimOfFluenceMxZ
0211                        leftLeafIx = find(shape_k(j,:)>0,1,'first');
0212                        rightLeafIx = find(shape_k(j,:)>0,1,'last');
0213                        if leftLeafIx > 1
0214                            plot(seqSubPlots(4),[.5 leftLeafIx-.5],j-[.5 .5] ,'w','LineWidth',2)
0215                            plot(seqSubPlots(4),[.5 leftLeafIx-.5],j+[.5 .5] ,'w','LineWidth',2)
0216                            plot(seqSubPlots(4),[ leftLeafIx-.5 leftLeafIx-.5],j+[.5 -.5] ,'w','LineWidth',2)
0217                        end
0218                        if rightLeafIx<dimOfFluenceMxX
0219                            plot(seqSubPlots(4),[dimOfFluenceMxX+.5 rightLeafIx+.5],j-[.5 .5] ,'w','LineWidth',2)
0220                            plot(seqSubPlots(4),[dimOfFluenceMxX+.5 rightLeafIx+.5],j+[.5 .5] ,'w','LineWidth',2)
0221                            plot(seqSubPlots(4),[ rightLeafIx+.5 rightLeafIx+.5],j+[.5 -.5] ,'w','LineWidth',2)
0222                        end
0223                        if isempty(rightLeafIx) && isempty (leftLeafIx)
0224                            plot(seqSubPlots(4),[dimOfFluenceMxX+.5 .5],j-[.5 .5] ,'w','LineWidth',2)
0225                            plot(seqSubPlots(4),[dimOfFluenceMxX+.5 .5],j+[.5 .5] ,'w','LineWidth',2)
0226                            plot(seqSubPlots(4),0.5*dimOfFluenceMxX*[1 1]+[0.5],j+[.5 -.5] ,'w','LineWidth',2)
0227                        end
0228                     end
0229                     
0230                     axis(seqSubPlots(4),'tight');
0231                     drawnow
0232                     pause(1);
0233                 end 
0234 
0235         shapes(:,:,k) = shape_k;
0236         shapesWeight(k) = d_k;
0237         D_k = D_k - shape_k;
0238         
0239         L_k = max(D_k(:)); % eq 5
0240         
0241     end
0242     
0243     sequencing.beam(i).numOfShapes  = k;
0244     sequencing.beam(i).shapes       = shapes(:,:,1:k);
0245     sequencing.beam(i).shapesWeight = shapesWeight(1:k)/numOfLevels*calFac;
0246     sequencing.beam(i).bixelIx      = 1+offset:numOfRaysPerBeam+offset;
0247     sequencing.beam(i).fluence      = D_0;
0248     
0249     sequencing.w(1+offset:numOfRaysPerBeam+offset,1) = D_0(indInFluenceMx)/numOfLevels*calFac;
0250 
0251     offset = offset + numOfRaysPerBeam;
0252 
0253 end
0254 
0255 resultGUI.w          = sequencing.w;
0256 resultGUI.wSequenced = sequencing.w;
0257 
0258 resultGUI.sequencing   = sequencing;
0259 resultGUI.apertureInfo = matRad_sequencing2ApertureInfo(sequencing,stf);
0260 
0261 doseSequencedDoseGrid = reshape(dij.physicalDose{1} * sequencing.w,dij.doseGrid.dimensions);
0262 % interpolate to ct grid for visualiation & analysis
0263 resultGUI.physicalDose = matRad_interp3(dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z, ...
0264                                         doseSequencedDoseGrid, ...
0265                                         dij.ctGrid.x,dij.ctGrid.y',dij.ctGrid.z);
0266 
0267 % if weights exists from an former DAO remove it
0268 if isfield(resultGUI,'wDao')
0269     resultGUI = rmfield(resultGUI,'wDao');
0270 end
0271 
0272 end
0273
| Generated by m2html © 2005