matRad_calcDoseInit

0001 
0002 matRad_cfg =  MatRad_Config.instance();
0003 
0004 % default: dose influence matrix computation
0005 if ~exist('calcDoseDirect','var')
0006     calcDoseDirect = false;
0007 end
0008 
0009 % to guarantee downwards compatibility with data that does not have
0010 % ct.x/y/z
0011 if ~any(isfield(ct,{'x','y','z'}))
0012     ct.x = ct.resolution.x*[0:ct.cubeDim(2)-1]-ct.resolution.x/2;
0013     ct.y = ct.resolution.y*[0:ct.cubeDim(1)-1]-ct.resolution.y/2;
0014     ct.z = ct.resolution.z*[0:ct.cubeDim(3)-1]-ct.resolution.z/2;
0015 end
0016 
0017 % set grids
0018 if ~isfield(pln,'propDoseCalc') || ...
0019    ~isfield(pln.propDoseCalc,'doseGrid') || ...
0020    ~isfield(pln.propDoseCalc.doseGrid,'resolution')
0021     % default values
0022     dij.doseGrid.resolution = matRad_cfg.propDoseCalc.defaultResolution;
0023 else
0024     % take values from pln strcut
0025     dij.doseGrid.resolution.x = pln.propDoseCalc.doseGrid.resolution.x;
0026     dij.doseGrid.resolution.y = pln.propDoseCalc.doseGrid.resolution.y;
0027     dij.doseGrid.resolution.z = pln.propDoseCalc.doseGrid.resolution.z;
0028 end
0029 
0030 dij.doseGrid.x = ct.x(1):dij.doseGrid.resolution.x:ct.x(end);
0031 dij.doseGrid.y = ct.y(1):dij.doseGrid.resolution.y:ct.y(end);
0032 dij.doseGrid.z = ct.z(1):dij.doseGrid.resolution.z:ct.z(end);
0033 
0034 dij.doseGrid.dimensions  = [numel(dij.doseGrid.y) numel(dij.doseGrid.x) numel(dij.doseGrid.z)];
0035 dij.doseGrid.numOfVoxels = prod(dij.doseGrid.dimensions);
0036 
0037 dij.ctGrid.resolution.x = ct.resolution.x;
0038 dij.ctGrid.resolution.y = ct.resolution.y;
0039 dij.ctGrid.resolution.z = ct.resolution.z;
0040 
0041 dij.ctGrid.x = ct.x;
0042 dij.ctGrid.y = ct.y;
0043 dij.ctGrid.z = ct.z;
0044 
0045 dij.ctGrid.dimensions  = [numel(dij.ctGrid.y) numel(dij.ctGrid.x) numel(dij.ctGrid.z)];
0046 dij.ctGrid.numOfVoxels = prod(dij.ctGrid.dimensions);
0047 
0048 % adjust isocenter internally for different dose grid
0049 offset = [dij.doseGrid.resolution.x - dij.ctGrid.resolution.x ...
0050           dij.doseGrid.resolution.y - dij.ctGrid.resolution.y ...
0051           dij.doseGrid.resolution.z - dij.ctGrid.resolution.z];
0052     
0053 for i = 1:numel(stf)
0054     stf(i).isoCenter = stf(i).isoCenter + offset;
0055 end
0056 
0057 %set up HU to rED or rSP conversion
0058 if ~isfield(pln,'propDoseCalc') || ~isfield(pln.propDoseCalc,'useGivenEqDensityCube')
0059     disableHUconversion = matRad_cfg.propDoseCalc.defaultUseGivenEqDensityCube;
0060 else
0061     disableHUconversion = pln.propDoseCalc.useGivenEqDensityCube;
0062 end
0063 
0064 %If we want to omit HU conversion check if we have a ct.cube ready
0065 if disableHUconversion && ~isfield(ct,'cube')
0066     matRad_cfg.dispWarning('HU Conversion requested to be omitted but no ct.cube exists! Will override and do the conversion anyway!');
0067     disableHUconversion = false;
0068 end
0069     
0070 % calculate rED or rSP from HU
0071 if disableHUconversion
0072     matRad_cfg.dispInfo('Omitting HU to rED/rSP conversion and using existing ct.cube!\n');
0073 else
0074     ct = matRad_calcWaterEqD(ct, pln);
0075 end
0076 
0077 % meta information for dij
0078 dij.numOfBeams         = pln.propStf.numOfBeams;
0079 dij.numOfScenarios     = 1;
0080 dij.numOfRaysPerBeam   = [stf(:).numOfRays];
0081 dij.totalNumOfBixels   = sum([stf(:).totalNumOfBixels]);
0082 dij.totalNumOfRays     = sum(dij.numOfRaysPerBeam);
0083 
0084 % check if full dose influence data is required
0085 if calcDoseDirect 
0086     numOfColumnsDij      = length(stf);
0087     numOfBixelsContainer = 1;
0088 else
0089     numOfColumnsDij      = dij.totalNumOfBixels;
0090     numOfBixelsContainer = ceil(dij.totalNumOfBixels/10);
0091 end
0092 
0093 % set up arrays for book keeping
0094 dij.bixelNum = NaN*ones(numOfColumnsDij,1);
0095 dij.rayNum   = NaN*ones(numOfColumnsDij,1);
0096 dij.beamNum  = NaN*ones(numOfColumnsDij,1);
0097 
0098 
0099 % Allocate space for dij.physicalDose sparse matrix
0100 for i = 1:dij.numOfScenarios
0101     dij.physicalDose{i} = spalloc(dij.doseGrid.numOfVoxels,numOfColumnsDij,1);
0102 end
0103 
0104 % Allocate memory for dose_temp cell array
0105 doseTmpContainer     = cell(numOfBixelsContainer,dij.numOfScenarios);
0106 
0107 % take only voxels inside patient
0108 VctGrid = [cst{:,4}];
0109 VctGrid = unique(vertcat(VctGrid{:}));
0110 
0111 % ignore densities outside of contours
0112 if ~isfield(pln,'propDoseCalc') || ~isfield(pln.propDoseCalc,'ignoreOutsideDensities')
0113     ignoreOutsideDensities = matRad_cfg.propDoseCalc.defaultIgnoreOutsideDensities;
0114 else
0115     ignoreOutsideDensities = pln.propDoseCalc.ignoreOutsideDensities;
0116 end
0117 
0118 if ignoreOutsideDensities
0119     eraseCtDensMask = ones(prod(ct.cubeDim),1);
0120     eraseCtDensMask(VctGrid) = 0;
0121     for i = 1:ct.numOfCtScen
0122         ct.cube{i}(eraseCtDensMask == 1) = 0;
0123     end
0124 end
0125 
0126 % Convert CT subscripts to linear indices.
0127 [yCoordsV_vox, xCoordsV_vox, zCoordsV_vox] = ind2sub(ct.cubeDim,VctGrid);
0128 
0129 % receive linear indices and grid locations from the dose grid
0130 tmpCube    = zeros(ct.cubeDim);
0131 tmpCube(VctGrid) = 1;
0132 % interpolate cube
0133 VdoseGrid = find(matRad_interp3(dij.ctGrid.x,  dij.ctGrid.y,   dij.ctGrid.z,tmpCube, ...
0134                                 dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z,'nearest'));
0135 
0136 % Convert CT subscripts to coarse linear indices.
0137 [yCoordsV_voxDoseGrid, xCoordsV_voxDoseGrid, zCoordsV_voxDoseGrid] = ind2sub(dij.doseGrid.dimensions,VdoseGrid);
0138 
0139 % load base data% load machine file
0140 fileName = [pln.radiationMode '_' pln.machine];
0141 try
0142    load([fileparts(mfilename('fullpath')) filesep 'basedata' filesep fileName]);
0143 catch
0144    matRad_cfg.dispError('Could not find the following machine file: %s\n',fileName); 
0145 end
0146 
0147 % compute SSDs
0148 stf = matRad_computeSSD(stf,ct);