matRad_plotPlan3D

 matRad function to visualize a plan in 3D 
 Stf is optional for plotting more detailed field contours in 
 visualization of the impinging beams.
 
 call
  rotMat = matRad_plotPlan3D(axesHandle,pln)
  rotMat = matRad_plotPlan3D(axesHandle,pln,stf)

 input
   axesHandle: handle to the axes the plan should be visualized in.
   pln:        matRad plan meta information struct
   stf:        optional steering information struct. if stf is passed and 
               not empty, the function will use the ray position  
               information to plot more detailed field contours than with 
               pln only

 output
   -

 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.

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

0001 function matRad_plotPlan3D(axesHandle,pln,stf)
0002 % matRad function to visualize a plan in 3D
0003 % Stf is optional for plotting more detailed field contours in
0004 % visualization of the impinging beams.
0005 %
0006 % call
0007 %  rotMat = matRad_plotPlan3D(axesHandle,pln)
0008 %  rotMat = matRad_plotPlan3D(axesHandle,pln,stf)
0009 %
0010 % input
0011 %   axesHandle: handle to the axes the plan should be visualized in.
0012 %   pln:        matRad plan meta information struct
0013 %   stf:        optional steering information struct. if stf is passed and
0014 %               not empty, the function will use the ray position
0015 %               information to plot more detailed field contours than with
0016 %               pln only
0017 %
0018 % output
0019 %   -
0020 %
0021 % References
0022 %   -
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 ishold(axesHandle)
0038     wasHold = true;
0039 else
0040     wasHold = false;
0041     hold(axesHandle,'on');
0042 end
0043 
0044 %nice pink ;)
0045 beamColor = [255 20 147]/255;
0046 
0047 %We perform a rudimentary visualization of the beam angles if there is no stf
0048 if nargin < 3 || isempty(stf)
0049     visFieldExtent = 50;
0050     visFieldPoints_bev = [-visFieldExtent 0 -visFieldExtent; ...
0051         visFieldExtent  0 -visFieldExtent; ...
0052         visFieldExtent  0  visFieldExtent; ...
0053         -visFieldExtent 0  visFieldExtent;
0054         -visFieldExtent 0 -visFieldExtent]; %repeat the first value for closed contour
0055     
0056     visFieldPoints_bev = visFieldPoints_bev ./ 2;
0057     visFieldPoints_bev = visFieldPoints_bev';
0058     
0059     fileName = [pln.radiationMode '_' pln.machine];
0060     %Get a SAD
0061     try
0062         load([fileparts(mfilename('fullpath')) filesep 'basedata' filesep fileName]);
0063         SAD = machine.meta.SAD;
0064     catch
0065         if strcmp(pln.radiationMode,'protons') || strcmp(pln.radiationMode,'carbon')
0066             SAD = 10000;
0067         else
0068             SAD = 1500;
0069         end
0070     end
0071     
0072     % default beam vector
0073     beamVector = [0 SAD 0];
0074    
0075     
0076     for beamIx = 1:pln.propStf.numOfBeams
0077         rotMat = matRad_getRotationMatrix(pln.propStf.gantryAngles(beamIx),pln.propStf.couchAngles(beamIx));
0078         beamIsoCenter = pln.propStf.isoCenter(beamIx,:);
0079         currBeamVector = rotMat*beamVector';        
0080         currBeamSource = beamIsoCenter - currBeamVector';
0081         currBeamOuterTarget = beamIsoCenter + currBeamVector';
0082         
0083         %Central ray
0084         line('XData',[currBeamSource(1) beamIsoCenter(1)],'YData',[currBeamSource(2) beamIsoCenter(2)],'ZData',[currBeamSource(3) beamIsoCenter(3)],'Parent',axesHandle,'LineWidth',2,'Color',beamColor);
0085         line('XData',[currBeamOuterTarget(1) beamIsoCenter(1)],'YData',[currBeamOuterTarget(2) beamIsoCenter(2)],'ZData',[currBeamOuterTarget(3) beamIsoCenter(3)],'Parent',axesHandle,'LineWidth',2,'Color',beamColor,'LineStyle',':');
0086         
0087         %Field rays
0088         for v = 1:size(visFieldPoints_bev,2)
0089             visFieldPoints_world(1:3,v) = rotMat*visFieldPoints_bev(1:3,v) + beamIsoCenter';
0090             
0091             %skipt the drawing of the first vertex since it is there twice
0092             if v==1 
0093                 continue;
0094             end
0095             
0096             
0097             
0098             %Draw the lines from the source point to the contour point
0099             line('XData',[currBeamSource(1) visFieldPoints_world(1,v)],'YData',[currBeamSource(2) visFieldPoints_world(2,v)],'ZData',[currBeamSource(3) visFieldPoints_world(3,v)],'Parent',axesHandle,'LineWidth',1,'Color',beamColor);
0100             %extend further
0101             currPointVector = visFieldPoints_world(1:3,v) - currBeamSource';
0102             currPointOuterTarget = currPointVector + visFieldPoints_world(1:3,v);
0103             line('XData',[currPointOuterTarget(1) visFieldPoints_world(1,v)],'YData',[currPointOuterTarget(2) visFieldPoints_world(2,v)],'ZData',[currPointOuterTarget(3) visFieldPoints_world(3,v)],'Parent',axesHandle,'LineWidth',1,'Color',beamColor,'LineStyle',':');
0104             
0105             %contour
0106             line('XData',[visFieldPoints_world(1,v) visFieldPoints_world(1,v-1)],...
0107                 'YData',[visFieldPoints_world(2,v) visFieldPoints_world(2,v-1)],...
0108                 'ZData',[visFieldPoints_world(3,v) visFieldPoints_world(3,v-1)],...
0109                 'Parent',axesHandle,'LineWidth',2,'Color',beamColor);
0110         end
0111     end
0112 else %We use the steering information to visualize the field contour
0113     for fieldIx = 1:numel(stf)
0114         beamTarget = stf(fieldIx).isoCenter;
0115         beamSource = stf(fieldIx).sourcePoint + stf(fieldIx).isoCenter;
0116         
0117         rotMat = matRad_getRotationMatrix(pln.propStf.gantryAngles(fieldIx),pln.propStf.couchAngles(fieldIx));
0118         
0119         bixelWidth = stf(fieldIx).bixelWidth;
0120         %Accumulate ray positions in matrix
0121         rayPos = [stf(fieldIx).ray(:).rayPos_bev];
0122         rayPos = reshape(rayPos,[3 numel(rayPos)/3])';
0123         
0124         %Compute a ray matrix with ones where a ray is
0125         %Maximum absolute values give extent in one direction
0126         symmMaxExtent = max([abs(min(rayPos)); abs(max(rayPos))]);
0127         %we want to have indices and not mm
0128         symmMaxExtent = symmMaxExtent ./ bixelWidth;
0129         %now make it symmetric
0130         symmMaxExtent = 2*symmMaxExtent + [1 0 1];
0131         %extra padding of one element in each direction to handle contours right
0132         symmMaxExtent = symmMaxExtent + [2 0 2];
0133         
0134         %Fill the ray matrix with ones at positions of rays
0135         rayMat = zeros(symmMaxExtent(1),symmMaxExtent(3));
0136         for rayIx = 1:size(rayPos,1)
0137             centerMat = (size(rayMat)-1) ./ 2 + 1;
0138             el2D = rayPos(rayIx,[1 3]);
0139             el2D = el2D ./ bixelWidth;
0140             el2DIx = centerMat + el2D;
0141             rayMat(el2DIx(1),el2DIx(2)) = 1;
0142         end
0143         %Create contour of the field
0144         fieldContour2D = contourc(rayMat,1);
0145         
0146         %Column in the contour matrix
0147         cColumn = 1;
0148         contourIx = 0;
0149         
0150         %Orientate the contour in 3D world space
0151         while cColumn <= size(fieldContour2D,2)
0152             contourIx = contourIx + 1;
0153             
0154             %Get contour data
0155             fieldContour3D(contourIx).level = fieldContour2D(1,cColumn);
0156             fieldContour3D(contourIx).numVertices = fieldContour2D(2,cColumn);
0157             fieldContour3D(contourIx).coord_bev = [ fieldContour2D(2,cColumn+1:cColumn+fieldContour3D(contourIx).numVertices); ... x
0158                 zeros(1,fieldContour3D(contourIx).numVertices);
0159                 fieldContour2D(1,cColumn+1:cColumn+fieldContour3D(contourIx).numVertices)]; ... y];
0160                 
0161             fieldContour3D(contourIx).coord_bev = bsxfun(@minus,fieldContour3D(contourIx).coord_bev,symmMaxExtent'./2) * bixelWidth;
0162             
0163             % Transform to world space
0164             % compute coordinates in lps coordinate system, i.e. rotate beam
0165             % geometry around fixed patient;
0166             
0167             for v = 1:fieldContour3D(contourIx).numVertices
0168                 fieldContour3D(contourIx).coord(1:3,v) = rotMat * fieldContour3D(contourIx).coord_bev(1:3,v) + beamTarget';
0169                 %Draw the lines from the source point to the contour point
0170                 line('XData',[beamSource(1) fieldContour3D(contourIx).coord(1,v)],'YData',[beamSource(2) fieldContour3D(contourIx).coord(2,v)],'ZData',[beamSource(3) fieldContour3D(contourIx).coord(3,v)],'Parent',axesHandle,'LineWidth',1,'Color',beamColor);
0171             end
0172             
0173             %Draw the contour in the isocenter (shows the field)
0174             line('XData',fieldContour3D(contourIx).coord(1,:),'YData',fieldContour3D(contourIx).coord(2,:),'ZData',fieldContour3D(contourIx).coord(3,:),'Parent',axesHandle,'LineWidth',2,'Color',beamColor);
0175             
0176             cColumn = cColumn + fieldContour3D(contourIx).numVertices + 1;
0177         end
0178         
0179     end
0180 end
0181 
0182 if ~wasHold
0183     hold(axesHandle,'off');
0184 end
0185 
0186 end