matRad_importDicomSteeringParticles

 matRad function to import a matRad stf struct from dicom RTPLAN data
 
 call
   [stf, pln] = matRad_importDicomSteeringParticles(ct, pln, rtPlanFile)

 input
   ct:             ct imported by the matRad_importDicomCt function
   pln:            matRad pln struct with meta information
   rtPlanFile:       name of RTPLAN DICOM file

 output
   stf             matRad stf struct
   pln:            matRad pln struct. 
                   Note: pln is input and output since pln.bixelWidth is 
                   determined here.

 References
   -
 Note
 not implemented - compensator. Fixed SAD.
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 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 [stf, pln] = matRad_importDicomSteeringParticles(ct, pln, rtPlanFile)
0002 % matRad function to import a matRad stf struct from dicom RTPLAN data
0003 %
0004 % call
0005 %   [stf, pln] = matRad_importDicomSteeringParticles(ct, pln, rtPlanFile)
0006 %
0007 % input
0008 %   ct:             ct imported by the matRad_importDicomCt function
0009 %   pln:            matRad pln struct with meta information
0010 %   rtPlanFile:       name of RTPLAN DICOM file
0011 %
0012 % output
0013 %   stf             matRad stf struct
0014 %   pln:            matRad pln struct.
0015 %                   Note: pln is input and output since pln.bixelWidth is
0016 %                   determined here.
0017 %
0018 % References
0019 %   -
0020 % Note
0021 % not implemented - compensator. Fixed SAD.
0022 %
0023 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0024 %
0025 % Copyright 2015 the matRad development team.
0026 %
0027 % This file is part of the matRad project. It is subject to the license
0028 % terms in the LICENSE file found in the top-level directory of this
0029 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
0030 % of the matRad project, including this file, may be copied, modified,
0031 % propagated, or distributed except according to the terms contained in the
0032 % LICENSE file.
0033 %
0034 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0035 
0036 %% load plan file
0037 % load machine data
0038 
0039 dlgBaseDataText = ['Import steering information from DICOM Plan.','Choose corresponding matRad base data for ', ...
0040         pln.radiationMode, '.'];
0041 % messagebox only necessary for non windows users
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 % RT Plan consists only on meta information
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 % Number of Beams from plan
0060 numOfBeamsPlan = length(pln.propStf.gantryAngles);
0061 
0062 % use only the treatment beams
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 % reinitialize the BeamSeqNames and length, as the Seq itself is reduced.
0076 BeamSeqNames = fieldnames(BeamSeq);
0077 
0078 % remove empty ControlPointSequences
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 % check if number of beams correspond
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 %% generate stf struct
0099 % surfaceEntry = BeamSeq.Item_1.IonControlPointSequence.Item_1.SurfaceEntryPoint;
0100 
0101 % Preallocate stf
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     % there might be several SAD's, e.g. compensator?
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         % now loop over ControlPointSequences
0126     ControlPointSeqNames = fieldnames(ControlPointSeq);
0127     numOfContrPointSeq = length(ControlPointSeqNames);
0128     % create empty helper matrix
0129     temporarySteering = zeros(0,8);
0130     for currContr = 1:numOfContrPointSeq
0131         currContrSeq = ControlPointSeq.(ControlPointSeqNames{currContr});
0132         % get energy, equal for all coming elements in the next loop
0133         currEnergy = currContrSeq.NominalBeamEnergy;
0134         % get focusValue
0135         currFocus = unique(currContrSeq.ScanningSpotSize);
0136         % get the Spotpositions
0137         numOfScanSpots = currContrSeq.NumberOfScanSpotPositions;
0138         % x is 1, 3, 5 ...; y 2, 4, 6,
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             % rangeshifter identification
0144             rashiID = currContrSeq.RangeShifterSettingsSequence.Item_1.ReferencedRangeShifterNumber;
0145             % rangeshifter waterequivalent thickness
0146             rashiWeThickness = currContrSeq.RangeShifterSettingsSequence.Item_1.RangeShifterWaterEquivalentThickness;
0147             % rangeshifter isocenter to range shifter distance
0148             rashiIsoRangeDist = currContrSeq.RangeShifterSettingsSequence.Item_1.IsocenterToRangeShifterDistance;
0149         elseif currContr == 1
0150             rashiID = 0;
0151             rashiWeThickness = 0;
0152             rashiIsoRangeDist = 0;            
0153         else
0154             % in this case range shifter settings has not changed between this
0155             % and previous control sequence, so reuse values.
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     % finds all unique rays and saves them in to the stf
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     % saves all energies and weights to their corresponding ray
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         % helpers to construct something like a(:).b = c.b(:) after this
0184         % loop
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     % reassign to preserve data structure
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     % getting some information of the rays
0201     % clean up energies, so they appear only one time per energy
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     % save total number of bixels
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 %         w = [w stf(currBeam).ray(j).weight];
0217     end
0218     
0219     stf(i).totalNumOfBixels = numOfBixels;
0220     
0221     % get bixelwidth
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     % coordinate transformation with rotation matrix.
0239     % use transpose matrix because we are working with row vectors
0240     rotMat_vectors_T = transpose(matRad_getRotationMatrix(pln.propStf.gantryAngles(i),pln.propStf.couchAngles(i)));
0241 
0242     % Rotated Source point (1st gantry, 2nd couch)
0243     stf(i).sourcePoint = stf(i).sourcePoint_bev*rotMat_vectors_T;
0244     
0245     % Save ray and target position in lps system.
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     % book keeping & calculate focus index
0252     for j = 1:stf(i).numOfRays
0253             stf(i).numOfBixelsPerRay(j) = numel([stf(i).ray(j).energy]);
0254     end
0255     
0256     % use the original machine energies
0257     for j = 1:stf(i).numOfRays
0258         % loop over all energies
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     % get focusIx instead of focusFWHM
0271     for j = 1:stf(i).numOfRays
0272         % loop over all energies
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