0001 function collimation = matRad_importFieldShapes(beamSequence, fractionSequence)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031 counter = 0;
0032 maximumExtent = 0;
0033 tmpCollimation.Fields = struct;
0034
0035
0036
0037 meta = {'NominalBeamEnergy' 'Energy';'GantryAngle' 'GantryAngle';...
0038 'PatientSupportAngle' 'CouchAngle';'SourceToSurfaceDistance' 'SSD'};
0039
0040
0041 beamSeqNames = fields(beamSequence);
0042 for i = 1:length(beamSeqNames)
0043
0044 currBeamSeq = beamSequence.(beamSeqNames{i});
0045 cumWeight = 0;
0046
0047
0048 tmpCollimation.beamMeterset(i) = fractionSequence.ReferencedBeamSequence.(beamSeqNames{i}).BeamMeterset;
0049
0050
0051 currDeviceSeq = beamSequence.(beamSeqNames{i}).BeamLimitingDeviceSequence;
0052 currDeviceSeqNames = fieldnames(currDeviceSeq);
0053
0054
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
0069
0070
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
0087 tmpCollimation.Fields(counter).LeafPos{length(currDeviceSeqNames),1} = [];
0088
0089
0090
0091
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
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
0117 maximumExtent = max(maximumExtent, max(abs(currLeafPos)));
0118
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
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
0143
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
0152 tmpCollimation.Fields(counter).BeamIndex = i;
0153 end
0154 end
0155 tmpCollimation.numOfFields = counter;
0156
0157
0158
0159 convResolution = .5;
0160 tmpCollimation.convResolution = convResolution;
0161
0162
0163 shapeLimit = ceil(maximumExtent / convResolution);
0164
0165
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
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
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
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
0194 p3 = 1;
0195 p4 = 2*shapeLimit;
0196 end
0197
0198
0199
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
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