matRad_generateStf

Purpose ^

matRad steering information generation

Synopsis ^

function stf = matRad_generateStf(ct,cst,pln,visMode)

Description ^

 matRad steering information generation
 
 call
   stf = matRad_generateStf(ct,cst,pln,visMode)

 input
   ct:         ct cube
   cst:        matRad cst struct
   pln:        matRad plan meta information struct
   visMode:    toggle on/off different visualizations by setting this value to 1,2,3 (optional)

 output
   stf:        matRad steering information struct

 References
   -

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

 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 stf = matRad_generateStf(ct,cst,pln,visMode)
0002 % matRad steering information generation
0003 %
0004 % call
0005 %   stf = matRad_generateStf(ct,cst,pln,visMode)
0006 %
0007 % input
0008 %   ct:         ct cube
0009 %   cst:        matRad cst struct
0010 %   pln:        matRad plan meta information struct
0011 %   visMode:    toggle on/off different visualizations by setting this value to 1,2,3 (optional)
0012 %
0013 % output
0014 %   stf:        matRad steering information struct
0015 %
0016 % References
0017 %   -
0018 %
0019 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0020 %
0021 % Copyright 2015 the matRad development team.
0022 %
0023 % This file is part of the matRad project. It is subject to the license
0024 % terms in the LICENSE file found in the top-level directory of this
0025 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
0026 % of the matRad project, including this file, may be copied, modified,
0027 % propagated, or distributed except according to the terms contained in the
0028 % LICENSE file.
0029 %
0030 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0031 
0032 matRad_cfg = MatRad_Config.instance();
0033 
0034 matRad_cfg.dispInfo('matRad: Generating stf struct... ');
0035 
0036 if nargin < 4
0037     visMode = 0;
0038 end
0039 
0040 if numel(pln.propStf.gantryAngles) ~= numel(pln.propStf.couchAngles)
0041     matRad_cfg.dispError('Inconsistent number of gantry and couch angles.');
0042 end
0043 
0044 if pln.propStf.bixelWidth < 0 || ~isfinite(pln.propStf.bixelWidth)
0045    matRad_cfg.dispError('bixel width (spot distance) needs to be a real number [mm] larger than zero.');
0046 end
0047 
0048 % find all target voxels from cst cell array
0049 V = [];
0050 for i=1:size(cst,1)
0051     if isequal(cst{i,3},'TARGET') && ~isempty(cst{i,6})
0052         V = [V;vertcat(cst{i,4}{:})];
0053     end
0054 end
0055 
0056 % Remove double voxels
0057 V = unique(V);
0058 % generate voi cube for targets
0059 voiTarget    = zeros(ct.cubeDim);
0060 voiTarget(V) = 1;
0061     
0062 % add margin
0063 addmarginBool = matRad_cfg.propStf.defaultAddMargin;
0064 if isfield(pln,'propStf') && isfield(pln.propStf,'addMargin')
0065    addmarginBool = pln.propStf.addMargin; 
0066 end
0067 
0068 if addmarginBool
0069     voiTarget = matRad_addMargin(voiTarget,cst,ct.resolution,ct.resolution,true);
0070     V   = find(voiTarget>0);
0071 end
0072 
0073 % throw error message if no target is found
0074 if isempty(V)
0075     matRad_cfg.dispError('Could not find target.');
0076 end
0077 
0078 % Convert linear indices to 3D voxel coordinates
0079 [coordsY_vox, coordsX_vox, coordsZ_vox] = ind2sub(ct.cubeDim,V);
0080 
0081 % prepare structures necessary for particles
0082 fileName = [pln.radiationMode '_' pln.machine];
0083 try
0084    load([fileparts(mfilename('fullpath')) filesep 'basedata' filesep fileName]);
0085    SAD = machine.meta.SAD;
0086 catch
0087    matRad_cfg.dispError('Could not find the following machine file: %s',fileName); 
0088 end
0089 
0090 if strcmp(pln.radiationMode,'protons') || strcmp(pln.radiationMode,'carbon')
0091       
0092     availableEnergies = [machine.data.energy];
0093     availablePeakPos  = [machine.data.peakPos] + [machine.data.offset];
0094     
0095     if sum(availablePeakPos<0)>0
0096        matRad_cfg.dispError('at least one available peak position is negative - inconsistent machine file') 
0097     end
0098     %clear machine;
0099 end
0100 
0101 % calculate rED or rSP from HU
0102 ct = matRad_calcWaterEqD(ct, pln);
0103 
0104 % take only voxels inside patient
0105 V = [cst{:,4}];
0106 V = unique(vertcat(V{:}));
0107 
0108 % ignore densities outside of contours
0109 eraseCtDensMask = ones(prod(ct.cubeDim),1);
0110 eraseCtDensMask(V) = 0;
0111 for i = 1:ct.numOfCtScen
0112     ct.cube{i}(eraseCtDensMask == 1) = 0;
0113 end
0114 
0115 % Define steering file like struct. Prellocating for speed.
0116 stf = struct;
0117 
0118 % loop over all angles
0119 for i = 1:length(pln.propStf.gantryAngles)
0120     
0121     % Correct for iso center position. Whit this correction Isocenter is
0122     % (0,0,0) [mm]
0123     coordsX = coordsX_vox*ct.resolution.x - pln.propStf.isoCenter(i,1);
0124     coordsY = coordsY_vox*ct.resolution.y - pln.propStf.isoCenter(i,2);
0125     coordsZ = coordsZ_vox*ct.resolution.z - pln.propStf.isoCenter(i,3);
0126 
0127     % Save meta information for treatment plan
0128     stf(i).gantryAngle   = pln.propStf.gantryAngles(i);
0129     stf(i).couchAngle    = pln.propStf.couchAngles(i);
0130     stf(i).bixelWidth    = pln.propStf.bixelWidth;
0131     stf(i).radiationMode = pln.radiationMode;
0132     stf(i).SAD           = SAD;
0133     stf(i).isoCenter     = pln.propStf.isoCenter(i,:);
0134         
0135     % Get the (active) rotation matrix. We perform a passive/system
0136     % rotation with row vector coordinates, which would introduce two
0137     % inversions / transpositions of the matrix, thus no changes to the
0138     % rotation matrix are necessary
0139     rotMat_system_T = matRad_getRotationMatrix(pln.propStf.gantryAngles(i),pln.propStf.couchAngles(i));
0140     
0141     rot_coords = [coordsX coordsY coordsZ]*rotMat_system_T;
0142     
0143     % project x and z coordinates to isocenter
0144     coordsAtIsoCenterPlane(:,1) = (rot_coords(:,1)*SAD)./(SAD + rot_coords(:,2));
0145     coordsAtIsoCenterPlane(:,2) = (rot_coords(:,3)*SAD)./(SAD + rot_coords(:,2));
0146     
0147     % Take unique rows values for beamlets positions. Calculate position of
0148     % central ray for every bixel
0149     rayPos = unique(pln.propStf.bixelWidth*round([           coordsAtIsoCenterPlane(:,1) ... 
0150                                                   zeros(size(coordsAtIsoCenterPlane,1),1) ...
0151                                                              coordsAtIsoCenterPlane(:,2)]/pln.propStf.bixelWidth),'rows');
0152                                                   
0153     % pad ray position array if resolution of target voxel grid not sufficient
0154     maxCtResolution = max([ct.resolution.x ct.resolution.y ct.resolution.z]);
0155     if pln.propStf.bixelWidth < maxCtResolution
0156         origRayPos = rayPos;
0157         for j = -floor(maxCtResolution/pln.propStf.bixelWidth):floor(maxCtResolution/pln.propStf.bixelWidth)
0158             for k = -floor(maxCtResolution/pln.propStf.bixelWidth):floor(maxCtResolution/pln.propStf.bixelWidth)
0159                 if abs(j)+abs(k)==0
0160                     continue;
0161                 end                
0162                 rayPos = [rayPos; origRayPos(:,1)+j*pln.propStf.bixelWidth origRayPos(:,2) origRayPos(:,3)+k*pln.propStf.bixelWidth];
0163             end
0164         end
0165      end
0166 
0167      % remove spaces within rows of bixels for DAO
0168      if pln.propOpt.runDAO
0169          % create single x,y,z vectors
0170          x = rayPos(:,1);
0171          y = rayPos(:,2);
0172          z = rayPos(:,3);
0173          uniZ = unique(z);
0174          for j = 1:numel(uniZ)
0175              x_loc = x(z == uniZ(j));
0176              x_min = min(x_loc);
0177              x_max = max(x_loc);
0178              x = [x; [x_min:pln.propStf.bixelWidth:x_max]'];
0179              y = [y; zeros((x_max-x_min)/pln.propStf.bixelWidth+1,1)];
0180              z = [z; uniZ(j)*ones((x_max-x_min)/pln.propStf.bixelWidth+1,1)];             
0181          end
0182          
0183          rayPos = [x,y,z];
0184      end
0185     
0186     % remove double rays
0187     rayPos = unique(rayPos,'rows');
0188     
0189     % Save the number of rays
0190     stf(i).numOfRays = size(rayPos,1);
0191     
0192     % Save ray and target position in beam eye's view (bev)
0193     for j = 1:stf(i).numOfRays
0194         stf(i).ray(j).rayPos_bev = rayPos(j,:);
0195         stf(i).ray(j).targetPoint_bev = [2*stf(i).ray(j).rayPos_bev(1) ...
0196                                                                SAD ...
0197                                          2*stf(i).ray(j).rayPos_bev(3)];
0198     end
0199     
0200     % source position in bev
0201     stf(i).sourcePoint_bev = [0 -SAD 0];
0202     
0203     % get (active) rotation matrix
0204     % transpose matrix because we are working with row vectors
0205     rotMat_vectors_T = transpose(matRad_getRotationMatrix(pln.propStf.gantryAngles(i),pln.propStf.couchAngles(i)));
0206     
0207     
0208     stf(i).sourcePoint = stf(i).sourcePoint_bev*rotMat_vectors_T;
0209     
0210     % Save ray and target position in lps system.
0211     for j = 1:stf(i).numOfRays
0212         stf(i).ray(j).rayPos      = stf(i).ray(j).rayPos_bev*rotMat_vectors_T;
0213         stf(i).ray(j).targetPoint = stf(i).ray(j).targetPoint_bev*rotMat_vectors_T;
0214         if strcmp(pln.radiationMode,'photons') 
0215             stf(i).ray(j).beamletCornersAtIso = [rayPos(j,:) + [+stf(i).bixelWidth/2,0,+stf(i).bixelWidth/2];...
0216                                                  rayPos(j,:) + [-stf(i).bixelWidth/2,0,+stf(i).bixelWidth/2];...
0217                                                  rayPos(j,:) + [-stf(i).bixelWidth/2,0,-stf(i).bixelWidth/2];...
0218                                                  rayPos(j,:) + [+stf(i).bixelWidth/2,0,-stf(i).bixelWidth/2]]*rotMat_vectors_T;
0219             stf(i).ray(j).rayCorners_SCD = (repmat([0, machine.meta.SCD - SAD, 0],4,1)+ (machine.meta.SCD/SAD) * ...
0220                                                              [rayPos(j,:) + [+stf(i).bixelWidth/2,0,+stf(i).bixelWidth/2];...
0221                                                               rayPos(j,:) + [-stf(i).bixelWidth/2,0,+stf(i).bixelWidth/2];...
0222                                                               rayPos(j,:) + [-stf(i).bixelWidth/2,0,-stf(i).bixelWidth/2];...
0223                                                               rayPos(j,:) + [+stf(i).bixelWidth/2,0,-stf(i).bixelWidth/2]])*rotMat_vectors_T;
0224         end
0225     end
0226     
0227     % loop over all rays to determine meta information for each ray
0228     stf(i).numOfBixelsPerRay = ones(1,stf(i).numOfRays);
0229     
0230     for j = stf(i).numOfRays:-1:1
0231 
0232         % ray tracing necessary to determine depth of the target
0233         [~,l,rho,~,~] = matRad_siddonRayTracer(stf(i).isoCenter, ...
0234                              ct.resolution, ...
0235                              stf(i).sourcePoint, ...
0236                              stf(i).ray(j).targetPoint, ...
0237                              [{ct.cube{1}} {voiTarget}]);
0238 
0239         % find appropriate energies for particles
0240        if strcmp(stf(i).radiationMode,'protons') || strcmp(stf(i).radiationMode,'carbon')
0241 
0242            % target hit
0243            if sum(rho{2}) > 0 
0244 
0245                 % compute radiological depths
0246                 % http://www.ncbi.nlm.nih.gov/pubmed/4000088, eq 14
0247                 radDepths = cumsum(l .* rho{1}); 
0248 
0249                 % find target entry & exit
0250                 diff_voi    = diff([rho{2}]);
0251                 targetEntry = radDepths(diff_voi == 1);
0252                 targetExit  = radDepths(diff_voi == -1);
0253 
0254                 if numel(targetEntry) ~= numel(targetExit)
0255                     matRad_cfg.dispError('Inconsistency during ray tracing. Please check correct assignment and overlap priorities of structure types OAR & TARGET.');
0256                 end
0257 
0258                 stf(i).ray(j).energy = [];
0259 
0260                 % Save energies in stf struct
0261                 for k = 1:numel(targetEntry)
0262                     stf(i).ray(j).energy = [stf(i).ray(j).energy availableEnergies(availablePeakPos>=targetEntry(k)&availablePeakPos<=targetExit(k))];
0263                 end
0264   
0265                 % book keeping & calculate focus index
0266                 stf(i).numOfBixelsPerRay(j) = numel([stf(i).ray(j).energy]);
0267                 currentMinimumFWHM = matRad_interp1(machine.meta.LUT_bxWidthminFWHM(1,:)',...
0268                                              machine.meta.LUT_bxWidthminFWHM(2,:)',...
0269                                              pln.propStf.bixelWidth);
0270                 focusIx  =  ones(stf(i).numOfBixelsPerRay(j),1);
0271                 [~, vEnergyIx] = min(abs(bsxfun(@minus,[machine.data.energy]',...
0272                                 repmat(stf(i).ray(j).energy,length([machine.data]),1))));
0273 
0274                 % get for each spot the focus index
0275                 for k = 1:stf(i).numOfBixelsPerRay(j)                    
0276                     focusIx(k) = find(machine.data(vEnergyIx(k)).initFocus.SisFWHMAtIso > currentMinimumFWHM,1,'first');
0277                 end
0278 
0279                 stf(i).ray(j).focusIx = focusIx';
0280                  
0281             else % target not hit
0282                 stf(i).ray(j)               = [];
0283                 stf(i).numOfBixelsPerRay(j) = [];
0284            end
0285            
0286        elseif strcmp(stf(i).radiationMode,'photons')
0287            
0288          % book keeping for photons
0289          stf(i).ray(j).energy = machine.data.energy;
0290          
0291        else
0292           matRad_cfg.dispError('Error generating stf struct: invalid radiation modality.');
0293        end
0294        
0295     end
0296     
0297     % store total number of rays for beam-i
0298     stf(i).numOfRays = size(stf(i).ray,2);
0299      
0300     % post processing for particle remove energy slices
0301     if strcmp(stf(i).radiationMode,'protons') || strcmp(stf(i).radiationMode,'carbon')
0302         
0303         % get minimum energy per field
0304         minEnergy = min([stf(i).ray.energy]);
0305         maxEnergy = max([stf(i).ray.energy]);
0306         
0307         % get corresponding peak position
0308         availableEnergies = [machine.data.energy];
0309         minPeakPos  = machine.data(minEnergy == availableEnergies).peakPos;
0310         maxPeakPos  = machine.data(maxEnergy == availableEnergies).peakPos;
0311         
0312         % find set of energyies with adequate spacing
0313         if ~isfield(pln.propStf, 'longitudinalSpotSpacing')
0314             longitudinalSpotSpacing = matRad_cfg.propStf.defaultLongitudinalSpotSpacing;
0315         else
0316             longitudinalSpotSpacing = pln.propStf.longitudinalSpotSpacing;
0317         end
0318         
0319         stf(i).longitudinalSpotSpacing = longitudinalSpotSpacing;
0320         
0321         tolerance              = longitudinalSpotSpacing/10;
0322         availablePeakPos       = [machine.data.peakPos];
0323         
0324         useEnergyBool = availablePeakPos >= minPeakPos & availablePeakPos <= maxPeakPos;
0325         
0326         ixCurr = find(useEnergyBool,1,'first');
0327         ixRun  = ixCurr + 1;
0328         ixEnd  = find(useEnergyBool,1,'last');
0329 
0330         while ixRun <= ixEnd
0331             if abs(availablePeakPos(ixRun)-availablePeakPos(ixCurr)) < ...
0332                                     longitudinalSpotSpacing - tolerance
0333                 useEnergyBool(ixRun) = 0;
0334             else
0335                 ixCurr = ixRun;
0336             end
0337             ixRun = ixRun + 1;
0338         end
0339         
0340         for j = stf(i).numOfRays:-1:1
0341             for k = stf(i).numOfBixelsPerRay(j):-1:1
0342                 maskEnergy = stf(i).ray(j).energy(k) == availableEnergies;
0343                 if ~useEnergyBool(maskEnergy)
0344                     stf(i).ray(j).energy(k)     = [];
0345                     stf(i).ray(j).focusIx(k)    = [];
0346                     stf(i).numOfBixelsPerRay(j) = stf(i).numOfBixelsPerRay(j) - 1;
0347                 end
0348             end
0349             if isempty(stf(i).ray(j).energy)
0350                 stf(i).ray(j) = [];
0351                 stf(i).numOfBixelsPerRay(j) = [];
0352                 stf(i).numOfRays = stf(i).numOfRays - 1;
0353             end
0354         end
0355         
0356     end
0357     
0358     % save total number of bixels
0359     stf(i).totalNumOfBixels = sum(stf(i).numOfBixelsPerRay);
0360     
0361     % Show progress
0362     matRad_progress(i,length(pln.propStf.gantryAngles));
0363 
0364     %% visualization
0365     if visMode > 0
0366         
0367         clf;
0368         % first subplot: visualization in bev
0369         subplot(1,2,1)
0370         hold on
0371         
0372         % plot rotated target coordinates
0373         plot3(rot_coords(:,1),rot_coords(:,2),rot_coords(:,3),'r.')
0374         
0375         % surface rendering
0376         if visMode == 2
0377             
0378             % generate a 3D rectangular grid centered at isocenter in
0379             % voxel coordinates
0380             [X,Y,Z] = meshgrid((1:ct.cubeDim(2))-stf(i).isoCenter(1)/ct.resolution.x, ...
0381                                (1:ct.cubeDim(1))-stf(i).isoCenter(2)/ct.resolution.y, ...
0382                                (1:ct.cubeDim(3))-stf(i).isoCenter(3)/ct.resolution.z);
0383             
0384             % computes surface
0385             patSurfCube      = 0*ct.cube{1};
0386             idx              = [cst{:,4}];
0387             idx              = unique(vertcat(idx{:}));
0388             patSurfCube(idx) = 1;
0389             
0390             [f,v] = isosurface(X,Y,Z,patSurfCube,.5);
0391             
0392             % convert isosurface from voxel to [mm]
0393             v(:,1) = v(:,1)*ct.resolution.x;
0394             v(:,2) = v(:,2)*ct.resolution.y;
0395             v(:,3) = v(:,3)*ct.resolution.z;
0396             
0397             % rotate surface
0398             rotated_surface = v*rotMat_system_T;
0399             
0400             % surface rendering
0401             surface = patch('Faces',f,'Vertices',rotated_surface);
0402             set(surface,'FaceColor',[0 0 1],'EdgeColor','none','FaceAlpha',.4);
0403             lighting gouraud;
0404         
0405         end
0406         
0407         % plot projection matrix: coordinates at isocenter
0408         plot3(rayPos(:,1),rayPos(:,2),rayPos(:,3),'k.');
0409         
0410         % Plot matrix border of matrix at isocenter
0411         for j = 1:stf(i).numOfRays
0412             
0413             % Compute border for every bixels
0414             targetPoint_vox_X_1 = stf(i).ray(j).targetPoint_bev(:,1) + pln.propStf.bixelWidth;
0415             targetPoint_vox_Y_1 = stf(i).ray(j).targetPoint_bev(:,2);
0416             targetPoint_vox_Z_1 = stf(i).ray(j).targetPoint_bev(:,3) + pln.propStf.bixelWidth;
0417             
0418             targetPoint_vox_X_2 = stf(i).ray(j).targetPoint_bev(:,1) + pln.propStf.bixelWidth;
0419             targetPoint_vox_Y_2 = stf(i).ray(j).targetPoint_bev(:,2);
0420             targetPoint_vox_Z_2 = stf(i).ray(j).targetPoint_bev(:,3) - pln.propStf.bixelWidth;
0421             
0422             targetPoint_vox_X_3 = stf(i).ray(j).targetPoint_bev(:,1) - pln.propStf.bixelWidth;
0423             targetPoint_vox_Y_3 = stf(i).ray(j).targetPoint_bev(:,2);
0424             targetPoint_vox_Z_3 = stf(i).ray(j).targetPoint_bev(:,3) - pln.propStf.bixelWidth;
0425             
0426             targetPoint_vox_X_4 = stf(i).ray(j).targetPoint_bev(:,1) - pln.propStf.bixelWidth;
0427             targetPoint_vox_Y_4 = stf(i).ray(j).targetPoint_bev(:,2);
0428             targetPoint_vox_Z_4 = stf(i).ray(j).targetPoint_bev(:,3) + pln.propStf.bixelWidth;
0429             
0430             % plot
0431             plot3([stf(i).sourcePoint_bev(1) targetPoint_vox_X_1],[stf(i).sourcePoint_bev(2) targetPoint_vox_Y_1],[stf(i).sourcePoint_bev(3) targetPoint_vox_Z_1],'g')
0432             plot3([stf(i).sourcePoint_bev(1) targetPoint_vox_X_2],[stf(i).sourcePoint_bev(2) targetPoint_vox_Y_2],[stf(i).sourcePoint_bev(3) targetPoint_vox_Z_2],'g')
0433             plot3([stf(i).sourcePoint_bev(1) targetPoint_vox_X_3],[stf(i).sourcePoint_bev(2) targetPoint_vox_Y_3],[stf(i).sourcePoint_bev(3) targetPoint_vox_Z_3],'g')
0434             plot3([stf(i).sourcePoint_bev(1) targetPoint_vox_X_4],[stf(i).sourcePoint_bev(2) targetPoint_vox_Y_4],[stf(i).sourcePoint_bev(3) targetPoint_vox_Z_4],'g')
0435             
0436         end
0437         
0438         % Plot properties
0439         daspect([1 1 1]);
0440         view(0,-90);
0441         xlabel 'X [mm]'
0442         ylabel 'Y [mm]'
0443         zlabel 'Z [mm]'
0444         title ('Beam''s eye view')
0445         axis([-300 300 -300 300 -300 300]);
0446         
0447         % second subplot: visualization in lps coordinate system
0448         subplot(1,2,2)
0449         
0450         % Plot target coordinates whitout any rotation
0451         plot3(coordsX,coordsY,coordsZ,'r.')
0452         hold on;
0453         
0454         % Rotated projection matrix at isocenter
0455         isocenter_plane_coor = rayPos*rotMat_vectors_T;
0456         
0457         % Plot isocenter plane
0458         plot3(isocenter_plane_coor(:,1),isocenter_plane_coor(:,2),isocenter_plane_coor(:,3),'y.');
0459         
0460         % Plot rotated bixels border.
0461         for j = 1:stf(i).numOfRays
0462             % Generate rotated projection target points.
0463             targetPoint_vox_1_rotated = [stf(i).ray(j).targetPoint_bev(:,1) + pln.propStf.bixelWidth,stf(i).ray(j).targetPoint_bev(:,2),stf(i).ray(j).targetPoint_bev(:,3) + pln.propStf.bixelWidth]*rotMat_vectors_T;
0464             targetPoint_vox_2_rotated = [stf(i).ray(j).targetPoint_bev(:,1) + pln.propStf.bixelWidth,stf(i).ray(j).targetPoint_bev(:,2),stf(i).ray(j).targetPoint_bev(:,3) - pln.propStf.bixelWidth]*rotMat_vectors_T;
0465             targetPoint_vox_3_rotated = [stf(i).ray(j).targetPoint_bev(:,1) - pln.propStf.bixelWidth,stf(i).ray(j).targetPoint_bev(:,2),stf(i).ray(j).targetPoint_bev(:,3) - pln.propStf.bixelWidth]*rotMat_vectors_T;
0466             targetPoint_vox_4_rotated = [stf(i).ray(j).targetPoint_bev(:,1) - pln.propStf.bixelWidth,stf(i).ray(j).targetPoint_bev(:,2),stf(i).ray(j).targetPoint_bev(:,3) + pln.propStf.bixelWidth]*rotMat_vectors_T;
0467             
0468             % Plot rotated target points.
0469             plot3([stf(i).sourcePoint(1) targetPoint_vox_1_rotated(:,1)],[stf(i).sourcePoint(2) targetPoint_vox_1_rotated(:,2)],[stf(i).sourcePoint(3) targetPoint_vox_1_rotated(:,3)],'g')
0470             plot3([stf(i).sourcePoint(1) targetPoint_vox_2_rotated(:,1)],[stf(i).sourcePoint(2) targetPoint_vox_2_rotated(:,2)],[stf(i).sourcePoint(3) targetPoint_vox_2_rotated(:,3)],'g')
0471             plot3([stf(i).sourcePoint(1) targetPoint_vox_3_rotated(:,1)],[stf(i).sourcePoint(2) targetPoint_vox_3_rotated(:,2)],[stf(i).sourcePoint(3) targetPoint_vox_3_rotated(:,3)],'g')
0472             plot3([stf(i).sourcePoint(1) targetPoint_vox_4_rotated(:,1)],[stf(i).sourcePoint(2) targetPoint_vox_4_rotated(:,2)],[stf(i).sourcePoint(3) targetPoint_vox_4_rotated(:,3)],'g')
0473         end
0474         
0475         % surface rendering
0476         if visMode == 2
0477             surface = patch('Faces',f,'Vertices',v);
0478             set(surface,'FaceColor',[0 0 1],'EdgeColor','none','FaceAlpha',.4);
0479             lighting gouraud;
0480         end
0481         
0482         % labels etc.
0483         daspect([1 1 1]);
0484         view(0,-90);
0485         xlabel 'X [mm]'
0486         ylabel 'Y [mm]'
0487         zlabel 'Z [mm]'
0488         title 'lps coordinate system'
0489         axis([-300 300 -300 300 -300 300]);
0490         %pause(1);
0491     end
0492     
0493     % include rangeshifter data if not yet available
0494     if strcmp(pln.radiationMode, 'protons') || strcmp(pln.radiationMode, 'carbon')
0495         for j = 1:stf(i).numOfRays
0496             for k = 1:numel(stf(i).ray(j).energy)
0497                 stf(i).ray(j).rangeShifter(k).ID = 0;
0498                 stf(i).ray(j).rangeShifter(k).eqThickness = 0;
0499                 stf(i).ray(j).rangeShifter(k).sourceRashiDistance = 0;
0500             end
0501         end
0502     end
0503         
0504 end    
0505 
0506 end
| Generated by m2html © 2005