Permalink
Browse files

First Commit

  • Loading branch information...
adityaapte committed Mar 13, 2012
0 parents commit 14a5dfd8f6e29073659c27f145b5c3fef80c9eb5
Showing 2,078 changed files with 633,704 additions and 0 deletions.
Binary file not shown.
@@ -0,0 +1,9 @@
/StructNameMap.m/1.1/Mon Jan 04 22:31:08 2010/-ko/
/getCERRfiles.m/1.1/Thu Jul 09 22:30:30 2009/-ko/
/getDoseCohort.m/1.1/Mon Feb 22 21:00:58 2010/-ko/
/getMatchingIndex.m/1.1/Wed Jul 28 21:49:19 2010/-ko/
/myMetric.m/1.1/Thu Apr 15 16:46:13 2010/-ko/
/repSpaceHyp.m/1.1/Tue Aug 10 21:19:59 2010/-ko/
/select_outcomes.m/1.1/Sun Dec 27 23:23:07 2009/-ko/
/structureNameMapGUI.m/1.1/Tue Aug 10 21:19:31 2010/-ko/
/sumDose.m/1.1/Wed Dec 23 23:27:44 2009/-ko/
@@ -0,0 +1 @@
CERR_Data_Extraction
@@ -0,0 +1 @@
:ext:aapte@bior-03.ro.wucon.wustl.edu:/home/cvsroot/
@@ -0,0 +1,323 @@
function DVH_out = DVHwithShifts(planC,structNum,doseNum)
%function DVH_out = DVHwithShifts(planC,structNum,doseNum)
%
%APA, 11/09/2010
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by: Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
% This version of CERR is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR. If not, see <http://www.gnu.org/licenses/>.
%Number of Fractions and Trials
numFractions = 30;
numTrials = 100;
% Systematic Shifts in cm (Inter-treatment)
XdispSys_Std = 0;
YdispSys_Std = 0;
ZdispSys_Std = 0;
% Random Shifts in cm (inter-fraction)
XdispRnd_Std = 1;
YdispRnd_Std = 1;
ZdispRnd_Std = 1;
indexS = planC{end};
%Distances are in cm
XpdfMean = 0;
YpdfMean = 0;
ZpdfMean = 0;
%Get scan associated with doseNum
scanNum = getAssociatedScan(planC{indexS.dose}(doseNum).assocScanUID, planC);
if isempty(scanNum) %Assume dose is associated with this scan
scanNum = getStructureAssociatedScan(structNum,planC);
end
%Get reference transformation matrix for doseNum
if ~isempty(scanNum) && isempty(planC{indexS.dose}(doseNum).transM)
referenceTransM = planC{indexS.scan}(scanNum).transM;
else
referenceTransM = planC{indexS.dose}(doseNum).transM;
end
%Store the reference doseArray
referenceDoseArray = planC{indexS.dose}(doseNum).doseArray;
%Try and get a binWidth from stateS. If it doesnt exist, get it from
%the CERROptions file (allows this function to be called outside CERR)
global stateS;
if ~isempty(stateS) && isfield(stateS, 'optS') && isfield(stateS.optS, 'DVHBinWidth') && ~isempty(stateS.optS.DVHBinWidth)
binWidth = stateS.optS.DVHBinWidth;
else
optS = CERROptions;
binWidth = optS.DVHBinWidth;
end
%Compute DVH at given dose
[dosesCurrentV, volsV] = getDVH(structNum, doseNum, planC);
[doseBinsCurrentV, volsHistCurrentV] = doseHist(dosesCurrentV, volsV, binWidth);
%Divide doseArray into fractions
planC{indexS.dose}(doseNum).doseArray = planC{indexS.dose}(doseNum).doseArray / numFractions;
%Get systematic errors
deltaX_systematic = randn(1, numTrials) * XdispSys_Std + XpdfMean;
deltaY_systematic = randn(1, numTrials) * YdispSys_Std + YpdfMean;
deltaZ_systematic = randn(1, numTrials) * ZdispSys_Std + ZpdfMean;
hWait = waitbar(0,'Computing plan robustness...');
if isempty(referenceTransM)
referenceTransM_matrix = eye(4);
else
referenceTransM_matrix = referenceTransM;
end
%Obtain DVH calculation points in blocks
optS = planC{indexS.CERROptions};
ROIImageSize = [planC{indexS.scan}(scanNum).scanInfo(1).sizeOfDimension1 planC{indexS.scan}(scanNum).scanInfo(1).sizeOfDimension2];
deltaY = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
%Get raster segments for structure.
[segmentsM, planC, isError] = getRasterSegments(structNum, planC);
if isempty(segmentsM)
isError = 1;
end
numSegs = size(segmentsM,1);
%Relative sampling of ROI voxels in this place, compared to CT spacing.
%Set when rasterSegments are generated (usually on import).
sampleRate = optS.ROISampleRate;
%Sample the rows
indFullV = 1 : numSegs;
if sampleRate ~= 1
rV = 1 : length(indFullV);
rV([rem(rV+sampleRate-1,sampleRate)~=0]) = [];
indFullV = rV;
end
%Block process to avoid swamping on large structures
if isfield(optS, 'DVHBlockSize') & ~isempty(optS.DVHBlockSize)
DVHBlockSize = optS.DVHBlockSize;
else
DVHBlockSize = 5000;
end
blocks = ceil(length(indFullV)/DVHBlockSize);
start = 1;
for b = 1 : blocks
%Build the interpolation points matrix
dummy = zeros(1,DVHBlockSize * ROIImageSize(1));
x1V = dummy;
y1V = dummy;
z1V = dummy;
volsSectionV = dummy;
if start+DVHBlockSize > length(indFullV)
stop = length(indFullV);
else
stop = start + DVHBlockSize - 1;
end
indV = indFullV(start:stop);
mark = 1;
for i = indV
tmpV = segmentsM(i,1:10);
delta = tmpV(5) * sampleRate;
xV = tmpV(3): delta : tmpV(4);
len = length(xV);
rangeV = ones(1,len);
yV = tmpV(2) * rangeV;
zV = tmpV(1) * rangeV;
sliceThickness = tmpV(10);
%v = delta^2 * sliceThickness;
v = delta * (deltaY*sampleRate) * sliceThickness;
x1V(mark : mark + len - 1) = xV;
y1V(mark : mark + len - 1) = yV;
z1V(mark : mark + len - 1) = zV;
volsSectionV(mark : mark + len - 1) = v;
mark = mark + len;
end
%cut unused matrix elements
x1V = x1V(1:mark-1);
y1V = y1V(1:mark-1);
z1V = z1V(1:mark-1);
volsSectionV = volsSectionV(1:mark-1);
%Get transformation matrices for both dose and structure.
transMDose = getTransM('dose', doseNum, planC);
transMStruct = getTransM('struct', structNum, planC);
%Forward transform the structure's coordinates.
if ~isempty(transMStruct)
[x1V, y1V, z1V] = applyTransM(transMStruct, x1V, y1V, z1V);
end
dvhCalcPtsC{b} = [x1V', y1V', z1V'];
%Interpolate.
% [dosesSectionV] = getDoseAt(doseNum, x1V, y1V, z1V, planC);
% dosesV = [dosesV, dosesSectionV];
% volsV = [volsV, volsSectionV];
start = stop + 1;
end
%Loop over to generate DVH
for iTrial = 1:numTrials
deltaX = deltaX_systematic(iTrial) + randn(1, numFractions) * XdispRnd_Std + XpdfMean;
deltaY = deltaY_systematic(iTrial) + randn(1, numFractions) * YdispRnd_Std + YpdfMean;
deltaZ = deltaZ_systematic(iTrial) + randn(1, numFractions) * ZdispRnd_Std + ZpdfMean;
tmpDoseV = zeros(1,length(dosesCurrentV));
for iFraction = 1:numFractions
waitbar(((iTrial-1)*numFractions + iFraction)/(numTrials*numFractions),hWait)
transM = referenceTransM_matrix;
transM(1:3,4) = transM(1:3,4) + [deltaX(iFraction); deltaY(iFraction); deltaZ(iFraction)];
%Apply the new transM to dose
planC{indexS.dose}(doseNum).transM = transM;
%Get doses and volumes of points in structure.
%[dosesV, volsV] = getDVH(structNum, doseNum, planC);
volsV = [];
dosesV = [];
for b = 1 : blocks
x1V = dvhCalcPtsC{b}(:,1);
y1V = dvhCalcPtsC{b}(:,2);
z1V = dvhCalcPtsC{b}(:,3);
%Back transform the coordinates into the doses' coordinate system.
[x1V, y1V, z1V] = applyTransM(inv(transM), x1V, y1V, z1V);
%Interpolate.
[dosesSectionV] = getDoseAt(doseNum, x1V, y1V, z1V, planC);
dosesV = [dosesV, dosesSectionV];
volsV = [volsV, volsSectionV];
end
tmpDoseV = tmpDoseV + dosesV;
end
DVHm(:,iTrial) = single(tmpDoseV(:));
%Compute Mean and Std for binned dose
%for iTrial=1:numTrials
doseTrialV = DVHm(:,iTrial);
[doseBinsV, volsHistV] = doseHist(doseTrialV, volsV, binWidth);
doseBinsTmpV = doseBinsV;
volsHistTmpV = volsHistV;
if iTrial > 1
AddToCurrent = ~ismember(doseBins{1},doseBinsV);
AddToPrevious = ~ismember(doseBinsV,doseBins{1});
doseBinsV = [doseBinsV doseBins{1}(AddToCurrent)];
volsHistV = [volsHistV volsHist{1}(AddToCurrent)];
for jTrial = 1:iTrial-1
doseBins{jTrial} = [doseBins{jTrial} doseBinsTmpV(AddToPrevious)];
volsHist{jTrial} = [volsHist{jTrial} volsHistTmpV(AddToPrevious)];
[doseBins{jTrial},indSort] = sort(doseBins{jTrial});
volsHist{jTrial} = volsHist{jTrial}(indSort);
end
end
doseBins{iTrial} = doseBinsV;
volsHist{iTrial} = volsHistV;
[doseBins{iTrial},indSort] = sort(doseBins{iTrial});
volsHist{iTrial} = volsHist{iTrial}(indSort);
%end
end
close(hWait)
%Compute Std of binned volumes
blockSize = 100;
for iBlock = 1:ceil(length(volsHist{1})/blockSize)
if iBlock == ceil(length(volsHist{1})/blockSize)
indicesV = (iBlock-1)*blockSize+1:length(volsHist{1});
else
indicesV = (iBlock-1)*blockSize+1:iBlock*blockSize;
end
volsHistStdM = [];
for iTrial = 1:length(volsHist)
volsHistStdM(:,iTrial) = volsHist{iTrial}(indicesV);
end
volsHistStdV(indicesV) = std(volsHistStdM,0,2)';
end
%Reassign reference transM and doseArray to doseNum
planC{indexS.dose}(doseNum).transM = referenceTransM;
planC{indexS.dose}(doseNum).doseArray = referenceDoseArray;
%Compute Mean and Std Dev in blocks
numVoxels = length(volsV);
DVHBlockSize = 5000;
blocks = ceil(numVoxels/DVHBlockSize);
start = 1;
for b = 1 : blocks
if start+DVHBlockSize > numVoxels
stop = numVoxels;
else
stop = start + DVHBlockSize - 1;
end
DVH_block = DVHm(start:stop,:);
meanDoseV(start:stop) = mean(DVH_block,2);
stdDoseV(start:stop) = std(DVH_block,0,2);
start = stop + 1;
end
DVH_out = [meanDoseV(:)'; volsV(:)'];
Oops, something went wrong.

0 comments on commit 14a5dfd

Please sign in to comment.