matRad_importDicomCt

 matRad function to import dicom ct data
 
 call
   ct = matRad_importDicomCt(ctList, resolution, dicomMetaBool)
   ct = matRad_importDicomCt(ctList, resolution, dicomMetaBool, grid)
   ct = matRad_importDicomCt(ctList, resolution, dicomMetaBool, visBool)
   ct = matRad_importDicomCt(ctList, resolution, dicomMetaBool, grid, visBool)

 input
   ctList:         list of dicom ct files
   resolution:       resolution of the imported ct cube, i.e. this function
                   will interpolate to a different resolution if desired
   dicomMetaBool:  store complete dicom information if true
   grid:           optional: a priori grid specified for interpolation
   visBool:        optional: turn on/off visualization

 output
   ct:             matRad ct struct. Note that this 3D matlab array 
                   contains water euqivalent electron denisities.
                   Hounsfield units are converted using a standard lookup
                   table in matRad_calcWaterEqD

 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 ct = matRad_importDicomCt(ctList, resolution, dicomMetaBool, grid, visBool)
0002 % matRad function to import dicom ct data
0003 %
0004 % call
0005 %   ct = matRad_importDicomCt(ctList, resolution, dicomMetaBool)
0006 %   ct = matRad_importDicomCt(ctList, resolution, dicomMetaBool, grid)
0007 %   ct = matRad_importDicomCt(ctList, resolution, dicomMetaBool, visBool)
0008 %   ct = matRad_importDicomCt(ctList, resolution, dicomMetaBool, grid, visBool)
0009 %
0010 % input
0011 %   ctList:         list of dicom ct files
0012 %   resolution:       resolution of the imported ct cube, i.e. this function
0013 %                   will interpolate to a different resolution if desired
0014 %   dicomMetaBool:  store complete dicom information if true
0015 %   grid:           optional: a priori grid specified for interpolation
0016 %   visBool:        optional: turn on/off visualization
0017 %
0018 % output
0019 %   ct:             matRad ct struct. Note that this 3D matlab array
0020 %                   contains water euqivalent electron denisities.
0021 %                   Hounsfield units are converted using a standard lookup
0022 %                   table in matRad_calcWaterEqD
0023 %
0024 % References
0025 %   -
0026 %
0027 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0028 %
0029 % Copyright 2015 the matRad development team.
0030 %
0031 % This file is part of the matRad project. It is subject to the license
0032 % terms in the LICENSE file found in the top-level directory of this
0033 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
0034 % of the matRad project, including this file, may be copied, modified,
0035 % propagated, or distributed except according to the terms contained in the
0036 % LICENSE file.
0037 %
0038 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0039 
0040 matRad_cfg = MatRad_Config.instance();
0041 
0042 matRad_cfg.dispInfo('\nimporting ct-cube...');
0043 
0044 %% processing input variables
0045 if ~exist('visBool','var')
0046   visBool = 0;
0047 end
0048 
0049 % creation of ctInfo list
0050 numOfSlices = size(ctList,1);
0051 matRad_cfg.dispInfo('\ncreating info...')
0052 
0053 sliceThicknessStandard = true;
0054 for i = 1:numOfSlices
0055 
0056     if verLessThan('matlab','9')
0057         tmpDicomInfo = dicominfo(ctList{i,1});
0058     else
0059         tmpDicomInfo = dicominfo(ctList{i,1},'UseDictionaryVR',true);
0060     end
0061     
0062     % remember relevant dicom info - do not record everything as some tags
0063     % might not been defined for individual files
0064     ctInfo(i).PixelSpacing            = tmpDicomInfo.PixelSpacing;
0065     ctInfo(i).ImagePositionPatient    = tmpDicomInfo.ImagePositionPatient;
0066     ctInfo(i).SliceThickness          = tmpDicomInfo.SliceThickness;
0067     ctInfo(i).ImageOrientationPatient = tmpDicomInfo.ImageOrientationPatient;
0068     ctInfo(i).PatientPosition         = tmpDicomInfo.PatientPosition;
0069     ctInfo(i).Rows                    = tmpDicomInfo.Rows;
0070     ctInfo(i).Columns                 = tmpDicomInfo.Columns;
0071     ctInfo(i).Width                   = tmpDicomInfo.Width;
0072     ctInfo(i).Height                  = tmpDicomInfo.Height;
0073     ctInfo(i).RescaleSlope            = tmpDicomInfo.RescaleSlope;
0074     ctInfo(i).RescaleIntercept        = tmpDicomInfo.RescaleIntercept;
0075     
0076     %Problem due to some CT files using non-standard SpacingBetweenSlices
0077     
0078     if isempty(ctInfo(i).SliceThickness)
0079         %Print warning ocne
0080         if sliceThicknessStandard
0081             matRad_cfg.dispWarning('Non-standard use of SliceThickness Attribute (empty), trying to overwrite with SpacingBetweenSlices');
0082             sliceThicknessStandard = false;
0083         end
0084         ctInfo(i).SliceThickness = tmpDicomInfo.SpacingBetweenSlices;
0085     end
0086     
0087     if i == 1
0088         completeDicom = tmpDicomInfo;
0089     end
0090     
0091     matRad_progress(i,numOfSlices);
0092 end
0093 
0094 % adjusting sequence of slices (filenames may not be ordered propperly....
0095 % e.g. CT1.dcm, CT10.dcm, CT100zCoordList = [ctInfo.ImagePositionPatient(1,3)]';.dcm, CT101.dcm,...
0096 CoordList = [ctInfo.ImagePositionPatient]';
0097 [~, indexing] = sort(CoordList(:,3)); % get sortation from z-coordinates
0098 
0099 ctList = ctList(indexing);
0100 ctInfo = ctInfo(indexing);
0101 
0102 %% check data set for consistency
0103 if size(unique([ctInfo.PixelSpacing]','rows'),1) > 1
0104     matRad_cfg.dispError('Different pixel size in different CT slices');
0105 end
0106 
0107 coordsOfFirstPixel = [ctInfo.ImagePositionPatient];
0108 if numel(unique(coordsOfFirstPixel(1,:))) > 1 || numel(unique(coordsOfFirstPixel(2,:))) > 1
0109     matRad_cfg.dispError('Ct slices are not aligned');
0110 end
0111 if sum(diff(coordsOfFirstPixel(3,:))<=0) > 0
0112     matRad_cfg.dispError('Ct slices not monotonically increasing');
0113 end
0114 if numel(unique([ctInfo.Rows])) > 1 || numel(unique([ctInfo.Columns])) > 1
0115     matRad_cfg.dispError('Ct slice sizes inconsistent');
0116 end
0117 
0118 
0119 %% checking the patient position
0120 % As of now, the matRad treatment planning system is only valid for
0121 % patients in a supine position. Other orientations (e.g. prone, decubitus
0122 % left/right) are not supported.
0123 % Defined Terms:
0124 % HFP     Head First-Prone                  (not supported)
0125 % HFS     Head First-Supine                 (supported)
0126 % HFDR    Head First-Decubitus Right        (not supported)
0127 % HFDL    Head First-Decubitus Left         (not supported)
0128 % FFDR    Feet First-Decubitus Right        (not supported)
0129 % FFDL    Feet First-Decubitus Left         (not supported)
0130 % FFP     Feet First-Prone                  (not supported)
0131 % FFS     Feet First-Supine                 (supported)
0132 
0133 if isempty(regexp(ctInfo(1).PatientPosition,{'S','P'}, 'once'))
0134     matRad_cfg.dispError(['This Patient Position is not supported by matRad.'...
0135         ' As of now only ''HFS'' (Head First-Supine), ''FFS'''...
0136         ' (Feet First-Supine), '...    
0137         '''HFP'' (Head First-Prone), and ''FFP'''...
0138         ' (Feet First-Prone) can be processed.'])    
0139 end
0140 
0141 %% creation of ct-cube
0142 matRad_cfg.dispInfo('reading slices...')
0143 origCt = zeros(ctInfo(1).Height, ctInfo(1).Width, numOfSlices);
0144 for i = 1:numOfSlices
0145     currentFilename = ctList{i};
0146     [currentImage, map] = dicomread(currentFilename);
0147     origCt(:,:,i) = currentImage(:,:); % creation of the ct cube
0148     
0149     % draw current ct-slice
0150     if visBool
0151         if ~isempty(map)
0152             image(ind2rgb(uint8(63*currentImage/max(currentImage(:))),map));
0153             xlabel('x [voxelnumber]')
0154             ylabel('y [voxelnumber]')
0155             title(['Slice # ' int2str(i) ' of ' int2str(numOfSlices)])
0156         else
0157             image(ind2rgb(uint8(63*currentImage/max(currentImage(:))),bone));
0158             xlabel('x [voxelnumber]')
0159             ylabel('y [voxelnumber]')
0160             title(['Slice # ' int2str(i) ' of ' int2str(numOfSlices)])
0161         end
0162         axis equal tight;
0163         pause(0.1);
0164     end
0165     matRad_progress(i,numOfSlices);
0166 end
0167 
0168 %% correction if not lps-coordinate-system
0169 % when using the physical coordinates (ctInfo.ImagePositionPatient) to
0170 % arrange the  slices in z-direction, there is no more need for mirroring
0171 % in the z-direction
0172 matRad_cfg.dispInfo('\nz-coordinates taken from ImagePositionPatient\n')
0173 
0174 % The x- & y-direction in lps-coordinates are specified in:
0175 % ImageOrientationPatient
0176 xDir = ctInfo(1).ImageOrientationPatient(1:3); % lps: [1;0;0]
0177 yDir = ctInfo(1).ImageOrientationPatient(4:6); % lps: [0;1;0]
0178 nonStandardDirection = false;
0179 
0180 % correct x- & y-direction
0181 %
0182 % if xDir(1) == 1 && xDir(2) == 0 && xDir(3) == 0
0183 %     matRad_cfg.dispInfo('x-direction OK\n')
0184 % elseif xDir(1) == -1 && xDir(2) == 0 && xDir(3) == 0
0185 %     matRad_cfg.dispInfo('\nMirroring x-direction...')
0186 %     origCt = flip(origCt,1);
0187 %     matRad_cfg.dispInfo('finished!\n')
0188 % else
0189 %     nonStandardDirection = true;
0190 % end
0191 %
0192 % if yDir(1) == 0 && yDir(2) == 1 && yDir(3) == 0
0193 %     matRad_cfg.dispInfo('y-direction OK\n')
0194 % elseif yDir(1) == 0 && yDir(2) == -1 && yDir(3) == 0
0195 %     matRad_cfg.dispInfo('\nMirroring y-direction...')
0196 %     origCt = flip(origCt,2);
0197 %     matRad_cfg.dispInfo('finished!\n')
0198 % else
0199 %     nonStandardDirection = true;
0200 % end
0201 
0202 if nonStandardDirection
0203     matRad_cfg.dispInfo(['Non-standard patient orientation.\n'...
0204         'CT might not fit to contoured structures\n'])
0205 end
0206 
0207 %% interpolate cube
0208 matRad_cfg.dispInfo('\nInterpolating CT cube...');
0209 if exist('grid','var')
0210     ct = matRad_interpDicomCtCube(origCt, ctInfo, resolution, grid);
0211 else
0212     ct = matRad_interpDicomCtCube(origCt, ctInfo, resolution);
0213 end
0214 matRad_cfg.dispInfo('finished!\n');
0215 
0216 %% remember some parameters of original dicom
0217 ct.dicomInfo.PixelSpacing            = ctInfo(1).PixelSpacing;
0218                                        tmp = [ctInfo.ImagePositionPatient];
0219 ct.dicomInfo.SlicePositions          = tmp(3,:);
0220 ct.dicomInfo.SliceThickness          = [ctInfo.SliceThickness];
0221 ct.dicomInfo.ImagePositionPatient    = ctInfo(1).ImagePositionPatient;
0222 ct.dicomInfo.ImageOrientationPatient = ctInfo(1).ImageOrientationPatient;
0223 ct.dicomInfo.PatientPosition         = ctInfo(1).PatientPosition;
0224 ct.dicomInfo.Width                   = ctInfo(1).Width;
0225 ct.dicomInfo.Height                  = ctInfo(1).Height;
0226 ct.dicomInfo.RescaleSlope            = ctInfo(1).RescaleSlope;
0227 ct.dicomInfo.RescaleIntercept        = ctInfo(1).RescaleIntercept;
0228 if isfield(completeDicom, 'Manufacturer')
0229 ct.dicomInfo.Manufacturer            = completeDicom.Manufacturer;
0230 end
0231 if isfield(completeDicom, 'ManufacturerModelName')
0232 ct.dicomInfo.ManufacturerModelName   = completeDicom.ManufacturerModelName;
0233 end
0234 if isfield(completeDicom, 'ConvolutionKernel')
0235 ct.dicomInfo.ConvolutionKernel       = completeDicom.ConvolutionKernel;
0236 end
0237 
0238 % store patientName only if user wants to
0239 if isfield(completeDicom,'PatientName') && dicomMetaBool == true
0240     ct.dicomInfo.PatientName         = completeDicom.PatientName;
0241 end
0242 if dicomMetaBool == true
0243     ct.dicomMeta                     = completeDicom;
0244 end
0245 
0246 ct.timeStamp = datestr(clock);
0247 
0248 % convert to Hounsfield units
0249 matRad_cfg.dispInfo('\nconversion of ct-Cube to Hounsfield units...');
0250 ct = matRad_calcHU(ct);
0251 matRad_cfg.dispInfo('finished!\n');
0252 
0253 end