matRad_importFieldShapes

Purpose ^

function to import collimator shapes from a DICOM RT plan

Synopsis ^

function collimation = matRad_importFieldShapes(beamSequence, fractionSequence)

Description ^

 function to import collimator shapes from a DICOM RT plan
 
 call
   collimation = matRad_importFieldShapes(beamSequence, fractionSequence)

 input
   beamSequence: struct containing the beamSequence elements from the RT plan    
   fractionSequence: struct containing the fractionGroupSequence elements from the RT plan    

 output
   collimation: struct with all meta information about the collimators and
   all field shape matrices 

 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 collimation = matRad_importFieldShapes(beamSequence, fractionSequence)
0002 % function to import collimator shapes from a DICOM RT plan
0003 %
0004 % call
0005 %   collimation = matRad_importFieldShapes(beamSequence, fractionSequence)
0006 %
0007 % input
0008 %   beamSequence: struct containing the beamSequence elements from the RT plan
0009 %   fractionSequence: struct containing the fractionGroupSequence elements from the RT plan
0010 %
0011 % output
0012 %   collimation: struct with all meta information about the collimators and
0013 %   all field shape matrices
0014 %
0015 % References
0016 %   -
0017 %
0018 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0019 %
0020 % Copyright 2015 the matRad development team.
0021 %
0022 % This file is part of the matRad project. It is subject to the license
0023 % terms in the LICENSE file found in the top-level directory of this
0024 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
0025 % of the matRad project, including this file, may be copied, modified,
0026 % propagated, or distributed except according to the terms contained in the
0027 % LICENSE file.
0028 %
0029 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0030 
0031 counter = 0;
0032 maximumExtent = 0;
0033 tmpCollimation.Fields = struct;
0034 
0035 % check for following meta data in every control point sequence
0036 % Format: 'DICOM Name Tag' 'Name in struct'; ...
0037 meta =  {'NominalBeamEnergy' 'Energy';'GantryAngle' 'GantryAngle';...
0038         'PatientSupportAngle' 'CouchAngle';'SourceToSurfaceDistance' 'SSD'};
0039 
0040 % extract field information
0041 beamSeqNames = fields(beamSequence);
0042 for i = 1:length(beamSeqNames)
0043     
0044     currBeamSeq = beamSequence.(beamSeqNames{i});
0045     cumWeight = 0;
0046     
0047     % get total MU applied by beam i
0048     tmpCollimation.beamMeterset(i) = fractionSequence.ReferencedBeamSequence.(beamSeqNames{i}).BeamMeterset;
0049     
0050     % get collimator device types
0051     currDeviceSeq = beamSequence.(beamSeqNames{i}).BeamLimitingDeviceSequence;
0052     currDeviceSeqNames = fieldnames(currDeviceSeq);
0053     
0054     % set device specific parameters
0055     device = struct;
0056     for j = 1:length(currDeviceSeqNames)
0057         currLimitsSeq = currDeviceSeq.(currDeviceSeqNames{j});
0058         device(j).DeviceType = currLimitsSeq.RTBeamLimitingDeviceType;
0059         device(j).NumOfLeafs = currLimitsSeq.NumberOfLeafJawPairs;
0060         device(j).Direction = device(j).DeviceType(end); 
0061         if strncmpi(device(j).DeviceType,'MLC',3)
0062            device(j).Limits = currLimitsSeq.LeafPositionBoundaries; 
0063         end
0064     end
0065     tmpCollimation.Devices{i} = device;
0066     
0067     currControlPointSeqNames = fieldnames(currBeamSeq.ControlPointSequence);
0068     % all meta informations must be defined in the first control point sequence
0069     % see DICOM Doc Sec. C.8.8.14.5
0070     % http://dicom.nema.org/MEDICAL/Dicom/2015c/output/chtml/part03/sect_C.8.8.14.5.html
0071     for j = 1:length(meta)
0072         try
0073             FieldMeta.(meta{j,2}) = currBeamSeq.ControlPointSequence.(currControlPointSeqNames{1}).(meta{j,1});
0074         catch
0075             warning(['Field ' meta{j,1} ' not found on beam sequence ' beamSeqNames{i} ...
0076                      '. No field shape import performed!']);
0077             return;
0078         end
0079     end
0080     
0081     for j = 1:length(currControlPointSeqNames)
0082        counter = counter + 1;
0083        currControlPointElement = currBeamSeq.ControlPointSequence.(currControlPointSeqNames{j});      
0084         
0085        if isfield(currControlPointElement, 'BeamLimitingDevicePositionSequence')
0086            % get the leaf position for every device
0087            tmpCollimation.Fields(counter).LeafPos{length(currDeviceSeqNames),1} = [];
0088 
0089            % beam limiting device position sequence has to be defined on
0090            % the first control point and has to be defined on following
0091            % points only if it changes -> default initilation if counter > 1
0092            if counter > 1
0093                for k = 1:length(currDeviceSeqNames)
0094                    tmpCollimation.Fields(counter).LeafPos{k} = tmpCollimation.Fields(counter-1).LeafPos{k};
0095                end
0096            end
0097 
0098            for k = 1:length(currDeviceSeqNames)
0099 
0100                if isfield(currControlPointElement.BeamLimitingDevicePositionSequence,currDeviceSeqNames{k})
0101                    currLeafPos = currControlPointElement.BeamLimitingDevicePositionSequence.(currDeviceSeqNames{k}).LeafJawPositions;          
0102 
0103                    deviceIx = find(strcmp({device(:).DeviceType}, ...
0104                        currControlPointElement.BeamLimitingDevicePositionSequence.(currDeviceSeqNames{k}).RTBeamLimitingDeviceType));
0105 
0106                    if (length(currLeafPos) ~= 2 * device(deviceIx).NumOfLeafs)
0107                        warning(['Number of leafs/jaws does not match given number of leaf/jaw positions in control point sequence ' ...
0108                                 currControlPointSeqNames{j} ' on beam sequence ' beamSeqNames{i} ' for device ' ...
0109                                 device(deviceIx).DeviceType '. No field shape import performed!']);
0110                        return;
0111                    end
0112                    
0113                    % set left and right leaf positions
0114                    tmpCollimation.Fields(counter).LeafPos{deviceIx}(:,1) = currLeafPos(1:device(deviceIx).NumOfLeafs);
0115                    tmpCollimation.Fields(counter).LeafPos{deviceIx}(:,2) = currLeafPos(device(deviceIx).NumOfLeafs+1:end);
0116                    % find the total maximum extent of one beam (in any direction)
0117                    maximumExtent = max(maximumExtent, max(abs(currLeafPos))); % check opening direction
0118                    % check direction perpendicular to the openening for MLC
0119                    if strncmpi(device(k).DeviceType,'MLC',3)
0120                        maximumExtent = max(maximumExtent,max(abs(device(deviceIx).Limits)));
0121                    end
0122                end
0123            end
0124        else
0125            tmpCollimation.Fields(counter) = tmpCollimation.Fields(counter - 1);
0126        end
0127        
0128        % get field meta information
0129        if isfield(currControlPointElement, 'CumulativeMetersetWeight')      
0130            newCumWeight = currControlPointElement.CumulativeMetersetWeight;
0131            tmpCollimation.Fields(counter).Weight = (newCumWeight - cumWeight) / ...
0132                                     currBeamSeq.FinalCumulativeMetersetWeight * ...
0133                                     tmpCollimation.beamMeterset(i)/100;
0134            cumWeight = newCumWeight;
0135        else
0136            warning(['No CumulativeMetersetWeight found in control point sequence ' currControlPointSeqNames{j} ...
0137                     ' on beam ' beamSeqNames{i} '. No field shape import performed!']);
0138            return;
0139        end
0140        tmpCollimation.Fields(counter).SAD = currBeamSeq.SourceAxisDistance;
0141         
0142        % other meta information is only included in all control point
0143        % sequences if it changes during treatment, otherwise use FieldMeta
0144        for k = 1:length(meta)
0145            if isfield(currControlPointElement,meta{k,1})
0146                tmpCollimation.Fields(counter).(meta{k,2}) = currControlPointElement.(meta{k,1});
0147            else
0148                tmpCollimation.Fields(counter).(meta{k,2}) = FieldMeta.(meta{k,2});
0149            end
0150        end
0151        % save information which control point sequence belongs to which beam sequence
0152        tmpCollimation.Fields(counter).BeamIndex = i;
0153     end
0154 end
0155 tmpCollimation.numOfFields = counter;
0156 
0157 % field import works only if the leaf width is a multiple of the conv
0158 % resolution
0159 convResolution = .5; % [mm]
0160 tmpCollimation.convResolution = convResolution;
0161 
0162 % get temporary shape limits to calculate the shapes
0163 shapeLimit = ceil(maximumExtent / convResolution);
0164 
0165 % calculate field shapes from leaf positions
0166 maximumVoxelExtent = 0;
0167 [X,Y] = meshgrid(-shapeLimit:shapeLimit-1);
0168 for i = 1:length(tmpCollimation.Fields)
0169     shape = ones(2*shapeLimit); 
0170     beamIndex = tmpCollimation.Fields(i).BeamIndex;
0171     for j = 1:length(tmpCollimation.Devices{beamIndex})
0172         % check for ASYM and SYM jaws == type 1
0173         if strncmpi(tmpCollimation.Devices{beamIndex}(j).DeviceType,'ASYM',4)
0174             type = 1;
0175         elseif (strcmpi(tmpCollimation.Devices{beamIndex}(j).DeviceType,'X') || ...
0176                 strcmpi(tmpCollimation.Devices{beamIndex}(j).DeviceType,'Y'))
0177             type = 1;
0178         % MLC == type 2
0179         elseif strncmpi(tmpCollimation.Devices{beamIndex}(j).DeviceType,'MLC',3)
0180             type = 2;
0181         else
0182             warning(['Device type ' tmpCollimation.Devices{beamIndex}(j).DeviceType ...
0183                     ' not supported. Field shapes could not be imported!']);
0184             return;
0185         end
0186         for k = 1:tmpCollimation.Devices{beamIndex}(j).NumOfLeafs
0187             % determine corner points of the open area
0188             p1 = ceil(tmpCollimation.Fields(i).LeafPos{j}(k,1)/convResolution)+shapeLimit;
0189             p2 = ceil(tmpCollimation.Fields(i).LeafPos{j}(k,2)/convResolution)+shapeLimit+1;
0190             if type == 2
0191                 p3 = ceil(tmpCollimation.Devices{beamIndex}(j).Limits(k)/convResolution)+shapeLimit+1;
0192                 p4 = ceil(tmpCollimation.Devices{beamIndex}(j).Limits(k+1)/convResolution)+shapeLimit;
0193             else % for one dimensional collimation (ASMX/Y) other direction is fully open
0194                 p3 = 1;
0195                 p4 = 2*shapeLimit;
0196             end
0197 
0198             % set elements covered by the collimator to 0
0199             % differentiate between x and y direction
0200             if (p1 > 0) && (p1 <= 2*shapeLimit) && ...
0201                (p2 > 0) && (p2 <= 2*shapeLimit) && ...
0202                (p3 > 0) && (p3 <= 2*shapeLimit) && ...
0203                (p4 > 0) && (p4 <= 2*shapeLimit)
0204                 try
0205                     if strcmpi(tmpCollimation.Devices{beamIndex}(j).Direction, 'X')
0206                         shape(p3:p4,1:p1) = 0;
0207                         shape(p3:p4,p2:end) = 0;
0208                     elseif strcmpi(tmpCollimation.Devices{beamIndex}(j).Direction, 'Y')
0209                         shape(1:p1,p3:p4) = 0;
0210                         shape(p2:end,p3:p4) = 0;
0211                     else
0212                         warning(['Wrong collimation direction ' tmpCollimation.Devices{beamIndex}(j).Direction ...
0213                                  ' given for device ' tmpCollimation.Devices{beamIndex}(j).DeviceType ...
0214                                  '. Fields could not be imported.']);
0215                         return;
0216                     end
0217                 catch
0218                     warning('Error in setting field shapes. Field shapes could not be imported.')
0219                     return;
0220                 end
0221             end
0222         end
0223     end
0224     openVoxelDistance = [X(shape == 1); Y(shape == 1)];
0225     maximumVoxelExtent = max(maximumVoxelExtent, max(abs(openVoxelDistance)));
0226     tmpCollimation.Fields(i).Shape = shape;
0227 end
0228 
0229 tmpCollimation.fieldWidth = 2 * maximumVoxelExtent * convResolution;
0230 % truncate field shapes to a symmetrical field with limits maximumVoxelExtent
0231 voxelRange = (-maximumVoxelExtent+1:maximumVoxelExtent) + shapeLimit;
0232 for i = 1:length(tmpCollimation.Fields) 
0233      tmpCollimation.Fields(i).Shape = tmpCollimation.Fields(i).Shape(voxelRange, voxelRange);
0234 end
0235 collimation = tmpCollimation;
0236 end

| Generated by m2html © 2005