0001 function [stf, pln] = matRad_importDicomSteeringParticles(ct, pln, rtPlanFile)
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
0032
0033
0034
0035
0036
0037
0038
0039 dlgBaseDataText = ['Import steering information from DICOM Plan.','Choose corresponding matRad base data for ', ...
0040 pln.radiationMode, '.'];
0041
0042 if ~ispc
0043 uiwait(helpdlg(dlgBaseDataText,['DICOM import - ', pln.radiationMode, ' base data' ]));
0044 end
0045 [fileName,pathName] = uigetfile('*.mat', dlgBaseDataText);
0046 load([pathName filesep fileName]);
0047
0048 ix = find(fileName == '_');
0049 pln.machine = fileName(ix(1)+1:end-4);
0050
0051
0052 if verLessThan('matlab','9')
0053 rtPlanInfo = dicominfo(rtPlanFile{1});
0054 else
0055 rtPlanInfo = dicominfo(rtPlanFile{1},'UseDictionaryVR',true);
0056 end
0057 BeamSeq = rtPlanInfo.IonBeamSequence;
0058 BeamSeqNames = fieldnames(BeamSeq);
0059
0060 numOfBeamsPlan = length(pln.propStf.gantryAngles);
0061
0062
0063 for i = 1:length(BeamSeqNames)
0064 currBeamSeq = BeamSeq.(BeamSeqNames{i});
0065 try
0066 treatDelType = currBeamSeq.TreatmentDeliveryType;
0067 if ~strcmpi(treatDelType,'TREATMENT')
0068 BeamSeq = rmfield(BeamSeq,BeamSeqNames{i});
0069 end
0070 catch
0071 warning('Something went wrong while determining the type of the beam.');
0072 end
0073 end
0074
0075
0076 BeamSeqNames = fieldnames(BeamSeq);
0077
0078
0079 for i = 1:length(BeamSeqNames)
0080 currBeamSeq = BeamSeq.(BeamSeqNames{i});
0081 ControlPointSeq = currBeamSeq.IonControlPointSequence;
0082 ControlPointSeqNames = fieldnames(ControlPointSeq);
0083 numOfContrPointSeq = length(ControlPointSeqNames);
0084 for currContr = 1:numOfContrPointSeq
0085 currContrSeq = ControlPointSeq.(ControlPointSeqNames{currContr});
0086 if sum(currContrSeq.ScanSpotMetersetWeights) == 0
0087 ControlPointSeq = rmfield(ControlPointSeq,ControlPointSeqNames{currContr});
0088 end
0089 end
0090 BeamSeq.(BeamSeqNames{i}).IonControlPointSequence = ControlPointSeq;
0091 end
0092
0093
0094 if ~isequal(length(BeamSeqNames),numOfBeamsPlan)
0095 warning('Number of beams from beamsequences do not correspond to number of Gantry Angles');
0096 end
0097
0098
0099
0100
0101
0102 stf(length(BeamSeqNames)).gantryAngle = [];
0103 stf(length(BeamSeqNames)).couchAngle = [];
0104 stf(length(BeamSeqNames)).bixelWidth = [];
0105 stf(length(BeamSeqNames)).radiationMode = [];
0106 stf(length(BeamSeqNames)).SAD = [];
0107 stf(length(BeamSeqNames)).isoCenter = [];
0108 stf(length(BeamSeqNames)).sourcePoint_bev = [];
0109 stf(length(BeamSeqNames)).numOfRays = [];
0110 stf(length(BeamSeqNames)).numOfBixelsPerRay = [];
0111 stf(length(BeamSeqNames)).totalNumOfBixels = [];
0112 stf(length(BeamSeqNames)).ray = [];
0113
0114 for i = 1:length(BeamSeqNames)
0115 currBeamSeq = BeamSeq.(BeamSeqNames{i});
0116 ControlPointSeq = currBeamSeq.IonControlPointSequence;
0117 stf(i).gantryAngle = pln.propStf.gantryAngles(i);
0118 stf(i).couchAngle = pln.propStf.couchAngles(i);
0119 stf(i).bixelWidth = pln.propStf.bixelWidth;
0120 stf(i).radiationMode = pln.radiationMode;
0121
0122 stf(i).SAD = machine.meta.SAD;
0123 stf(i).isoCenter = pln.propStf.isoCenter(i,:);
0124 stf(i).sourcePoint_bev = [0 -stf(i).SAD 0];
0125
0126 ControlPointSeqNames = fieldnames(ControlPointSeq);
0127 numOfContrPointSeq = length(ControlPointSeqNames);
0128
0129 temporarySteering = zeros(0,8);
0130 for currContr = 1:numOfContrPointSeq
0131 currContrSeq = ControlPointSeq.(ControlPointSeqNames{currContr});
0132
0133 currEnergy = currContrSeq.NominalBeamEnergy;
0134
0135 currFocus = unique(currContrSeq.ScanningSpotSize);
0136
0137 numOfScanSpots = currContrSeq.NumberOfScanSpotPositions;
0138
0139 c1_help = currContrSeq.ScanSpotPositionMap(1:2:(2 * numOfScanSpots));
0140 c2_help = currContrSeq.ScanSpotPositionMap(2:2:(2 * numOfScanSpots));
0141 weight_help = currContrSeq.ScanSpotMetersetWeights;
0142 if isfield(currContrSeq, 'RangeShifterSettingsSequence')
0143
0144 rashiID = currContrSeq.RangeShifterSettingsSequence.Item_1.ReferencedRangeShifterNumber;
0145
0146 rashiWeThickness = currContrSeq.RangeShifterSettingsSequence.Item_1.RangeShifterWaterEquivalentThickness;
0147
0148 rashiIsoRangeDist = currContrSeq.RangeShifterSettingsSequence.Item_1.IsocenterToRangeShifterDistance;
0149 elseif currContr == 1
0150 rashiID = 0;
0151 rashiWeThickness = 0;
0152 rashiIsoRangeDist = 0;
0153 else
0154
0155
0156 end
0157 temporarySteering = [temporarySteering; c1_help c2_help ...
0158 (currEnergy * ones(numOfScanSpots,1)) weight_help (currFocus * ones(numOfScanSpots,1)) ...
0159 (rashiID * ones(numOfScanSpots,1)) (rashiWeThickness * ones(numOfScanSpots,1)) (rashiIsoRangeDist * ones(numOfScanSpots,1))];
0160 end
0161
0162
0163 [RayPosTmp, ~, ic] = unique(temporarySteering(:,1:2), 'rows');
0164 clear ray;
0165 for j = 1:size(RayPosTmp,1)
0166 stf(i).ray(j).rayPos_bev = double([RayPosTmp(j,1) 0 RayPosTmp(j,2)]);
0167 stf(i).ray(j).energy = [];
0168 stf(i).ray(j).focusFWHM = [];
0169 stf(i).ray(j).focusIx = [];
0170 stf(i).ray(j).weight = [];
0171 stf(i).ray(j).rangeShifter = struct();
0172 ray(j).ID = [];
0173 ray(j).eqThickness = [];
0174 ray(j).sourceRashiDistance = [];
0175 end
0176
0177
0178 for j = 1:size(temporarySteering,1)
0179 k = ic(j);
0180 stf(i).ray(k).energy = [stf(i).ray(k).energy double(temporarySteering(j,3))];
0181 stf(i).ray(k).focusFWHM = [stf(i).ray(k).focusFWHM double(temporarySteering(j,5))];
0182 stf(i).ray(k).weight = [stf(i).ray(k).weight double(temporarySteering(j,4)) / 1e6];
0183
0184
0185 ray(k).ID = [ray(k).ID double(temporarySteering(j,6))];
0186 ray(k).eqThickness = [ray(k).eqThickness double(temporarySteering(j,7))];
0187 ray(k).sourceRashiDistance = [ray(k).sourceRashiDistance double(temporarySteering(j,8))];
0188 end
0189
0190
0191 for j = 1:numel(ray)
0192 for k = 1:numel(ray(j).ID)
0193 stf(i).ray(j).rangeShifter(k).ID = ray(j).ID(k);
0194 stf(i).ray(j).rangeShifter(k).eqThickness = ray(j).eqThickness(k);
0195 stf(i).ray(j).rangeShifter(k).sourceRashiDistance = stf(i).SAD - ray(j).sourceRashiDistance(k);
0196 end
0197 end
0198
0199
0200
0201
0202 numOfRays = size(stf(i).ray,2);
0203 for l = 1:numOfRays
0204 stf(i).ray(l).energy = unique(stf(i).ray(l).energy);
0205 stf(i).ray(l).targetPoint_bev = [2*stf(i).ray(l).rayPos_bev(1) ...
0206 machine.meta.SAD ...
0207 2*stf(i).ray(l).rayPos_bev(3)];
0208 end
0209 stf(i).numOfRays = numel(stf(i).ray);
0210
0211
0212 numOfBixels = 0;
0213 for j = 1:numel(stf(i).ray)
0214 numOfBixels = numOfBixels + numel(stf(i).ray(j).energy);
0215 stf(i).numOfBixelsPerRay(j) = numel(stf(i).ray(j).energy);
0216
0217 end
0218
0219 stf(i).totalNumOfBixels = numOfBixels;
0220
0221
0222 bixelWidth_help = zeros(size(stf(i).ray,2),2);
0223 for j = 1:stf(i).numOfRays
0224 bixelWidth_help(j,1) = stf(i).ray(j).rayPos_bev(1);
0225 bixelWidth_help(j,2) = stf(i).ray(j).rayPos_bev(3);
0226 end
0227 bixelWidth_help1 = unique(round(1e3*bixelWidth_help(:,1))/1e3,'sorted');
0228 bixelWidth_help2 = unique(round(1e3*bixelWidth_help(:,2))/1e3,'sorted');
0229
0230 bixelWidth = unique([unique(diff(bixelWidth_help1))' unique(diff(bixelWidth_help2))']);
0231
0232 if numel(bixelWidth) == 1
0233 stf(i).bixelWidth = bixelWidth;
0234 else
0235 stf(i).bixelWidth = NaN;
0236 end
0237
0238
0239
0240 rotMat_vectors_T = transpose(matRad_getRotationMatrix(pln.propStf.gantryAngles(i),pln.propStf.couchAngles(i)));
0241
0242
0243 stf(i).sourcePoint = stf(i).sourcePoint_bev*rotMat_vectors_T;
0244
0245
0246 for j = 1:stf(i).numOfRays
0247 stf(i).ray(j).rayPos = stf(i).ray(j).rayPos_bev*rotMat_vectors_T;
0248 stf(i).ray(j).targetPoint = stf(i).ray(j).targetPoint_bev*rotMat_vectors_T;
0249 end
0250
0251
0252 for j = 1:stf(i).numOfRays
0253 stf(i).numOfBixelsPerRay(j) = numel([stf(i).ray(j).energy]);
0254 end
0255
0256
0257 for j = 1:stf(i).numOfRays
0258
0259 numOfEnergy = length(stf(i).ray(j).energy);
0260 for k = 1:numOfEnergy
0261 energyIndex = find(abs([machine.data(:).energy]-stf(i).ray(j).energy(k))<10^-2);
0262 if ~isempty(energyIndex)
0263 stf(i).ray(j).energy(k) = machine.data(energyIndex).energy;
0264 else
0265 error('No match between imported and machine data. Maybe wrong machine loaded.');
0266 end
0267 end
0268 end
0269
0270
0271 for j = 1:stf(i).numOfRays
0272
0273 numOfEnergy = length(stf(i).ray(j).energy);
0274 for k = 1:numOfEnergy
0275 energyTemp = stf(i).ray(j).energy(k);
0276 focusFWHM = stf(i).ray(j).focusFWHM(k);
0277 energyIxTemp = find([machine.data.energy] == energyTemp);
0278 focusIxTemp = find(abs([machine.data(energyIxTemp).initFocus.SisFWHMAtIso] - focusFWHM )< 10^-3);
0279 stf(i).ray(j).focusIx(k) = focusIxTemp;
0280 stf(i).ray(j).focusFWHM(k) = machine.data(energyIxTemp).initFocus.SisFWHMAtIso(stf(i).ray(j).focusIx(k));
0281 end
0282 end
0283
0284 stf(i).timeStamp = datestr(clock);
0285
0286 end
0287
0288 if any(isnan([stf(:).bixelWidth])) || numel(unique([stf(:).bixelWidth])) > 1
0289 pln.propStf.bixelWidth = NaN;
0290 else
0291 pln.propStf.bixelWidth = stf(1).bixelWidth;
0292 end
0293
0294 end