Skip to content

Commit

Permalink
Added option to specify dose agreement as the percent of max dose or …
Browse files Browse the repository at this point in the history 
…dose at each voxel in the reference dose.
  • Loading branch information
@adityaapte
adityaapte committed May 20, 2019
1 parent be3b918 commit 1643aec
Show file tree
Hide file tree
Showing 4 changed files with 128 additions and 144 deletions.
83 changes: 62 additions & 21 deletions CERR_core/PlanAnalysis/gamma/CERRGammafnc.m
View file
Expand Up @@ -112,46 +112,79 @@ function CERRGammafnc(gamaCommand, varargin)

% Get User Input for 2D dose difference and DTA
gammaGUIFig = figure('numbertitle','off','name','3-D Gamma Calculation','tag','CERRgammaInputGUI','position',...
[ScreenSize(3)/2 ScreenSize(4)/2 300 350],'MenuBar','none','Resize','off','CloseRequestFcn','CERRGammafnc(''GAMMACANCLE'')');
[ScreenSize(3)/2 ScreenSize(4)/2 320 350],'MenuBar','none','Resize','off','CloseRequestFcn','CERRGammafnc(''GAMMACANCLE'')');

bgColor = get(gammaGUIFig,'color');

% Text dose difference
uicontrol('parent',gammaGUIFig,'style','text','backgroundcolor',bgColor,'position',[10 310 140 30],'String', 'Dose Difference (% of max ref dose)');
uicontrol('parent',gammaGUIFig,'style','text','backgroundcolor',...
bgColor,'position',[5 310 70 30],'String', 'Dose Difference');
% Input dose difference
uicontrol('parent',gammaGUIFig,'style','Edit','backgroundcolor',[1 1 1],'position',[20 290 90 20],'String', '3','tag','InputDoseDiff' );

uicontrol('parent',gammaGUIFig,'style','Edit','backgroundcolor',...
[1 1 1],'position',[15 290 50 20],'String', '3','tag','InputDoseDiff' );

% Text dose difference method
uicontrol('parent',gammaGUIFig,'style','text','backgroundcolor',...
bgColor,'position',[80 310 140 30],'String', 'Method');
% Input dose difference method
uicontrol('parent',gammaGUIFig,'style','popup','backgroundcolor',...
[1 1 1],'position',[80 290 140 20],'String',...
{'% of max ref dose','% ref dose at each voxel'},'value',1,...
'tag','InputDoseDiffMethod' );

% Text DTA
uicontrol('parent',gammaGUIFig,'style','text','backgroundcolor',bgColor,'position',[160 310 140 30],'String', 'DTA (mm)');
uicontrol('parent',gammaGUIFig,'style','text','backgroundcolor',...
bgColor,'position',[200 310 120 30],'String', 'DTA (mm)');
% Input DTA
uicontrol('parent',gammaGUIFig,'style','Edit','backgroundcolor',[1 1 1],'position',[170 290 90 20],'String', '3','tag','InputDTA');
uicontrol('parent',gammaGUIFig,'style','Edit','backgroundcolor',...
[1 1 1],'position',[240 290 50 20],'String', '3','tag','InputDTA');

%Text Base Dose
uicontrol('parent',gammaGUIFig,'style','text','backgroundcolor',bgColor,'position',[10 250 150 20],'HorizontalAlignment','Left','String', 'Reference Dose');
uicontrol('Style', 'popup','String', doseFraction,'Position', [8 205 290 50],'tag','baseDoseGamma','backgroundcolor', [1 1 1]);
uicontrol('parent',gammaGUIFig,'style','text','backgroundcolor',bgColor,...
'position',[10 250 150 20],'HorizontalAlignment','Left','String', 'Reference Dose');
uicontrol('Style', 'popup','String', doseFraction,'Position', ...
[8 205 290 50],'tag','baseDoseGamma','backgroundcolor', [1 1 1]);

% Text Ref Dose
uicontrol('parent',gammaGUIFig,'style','text','backgroundcolor',bgColor,'position',[10 205 150 20],'HorizontalAlignment','Left','String', 'Evaluation Dose');
uicontrol('Style', 'popup','String', doseFraction','Position', [8 160 290 50], 'tag','refDoseGamma','backgroundcolor',[1 1 1]);
uicontrol('parent',gammaGUIFig,'style','text','backgroundcolor',bgColor,...
'position',[10 205 150 20],'HorizontalAlignment','Left','String', 'Evaluation Dose');
uicontrol('Style', 'popup','String', doseFraction','Position', ...
[8 160 290 50], 'tag','refDoseGamma','backgroundcolor',[1 1 1]);

% Text Threshold
uicontrol('parent',gammaGUIFig,'style','text','backgroundcolor',bgColor,'position',[10 140 140 30],'String', 'Threshold low dose (% max ref dose)');
uicontrol('parent',gammaGUIFig,'style','text','backgroundcolor',...
bgColor,'position',[10 140 100 30],'String', 'Threshold low dose (% max ref dose)');
% Input Threshold
uicontrol('parent',gammaGUIFig,'style','Edit','backgroundcolor',[1 1 1],'position',[160 150 50 20],'String', '5','tag','InputThreshold');

uicontrol('parent',gammaGUIFig,'style','Edit','backgroundcolor',...
[1 1 1],'position',[110 140 40 30],'String', '0.001','tag','InputThreshold');

% Text MaxDistance
uicontrol('parent',gammaGUIFig,'style','text','backgroundcolor',...
bgColor,'position',[170 140 110 30],'String', 'Max search distance factor (times DTA)');
% Input MaxDistance
uicontrol('parent',gammaGUIFig,'style','Edit','backgroundcolor',...
[1 1 1],'position',[280 140 30 30],'String', '3','tag','MaxDistanceFactor');

%Text Structure
uicontrol('parent',gammaGUIFig,'style','text','backgroundcolor',bgColor,'position',[10 110 150 20],'HorizontalAlignment','Left','String', 'Structure');
uicontrol('Style', 'popup','String', structNamC,'Position', [8 65 290 50],'tag','structure','backgroundcolor', [1 1 1]);
uicontrol('parent',gammaGUIFig,'style','text','backgroundcolor',...
bgColor,'position',[10 110 150 20],'HorizontalAlignment','Left','String', 'Structure');
uicontrol('Style', 'popup','String', structNamC,'Position', ...
[8 65 290 50],'tag','structure','backgroundcolor', [1 1 1]);

% Button GO
uicontrol('parent',gammaGUIFig,'style','pushbutton','position',[50 50 70 20],'String', 'Calculate','Callback', ['CERRGammafnc(''CALCGAMMA3D'')' ] );
uicontrol('parent',gammaGUIFig,'style','pushbutton','position',...
[50 50 70 20],'String', 'Calculate','Callback', 'CERRGammafnc(''CALCGAMMA3D'')');

% Button Cancel
uicontrol('parent',gammaGUIFig,'style','pushbutton','position',[180 50 70 20],'String', 'Cancel','Callback', 'CERRGammafnc(''GAMMACANCLE'')');
uicontrol('parent',gammaGUIFig,'style','pushbutton','position',...
[180 50 70 20],'String', 'Cancel','Callback', 'CERRGammafnc(''GAMMACANCLE'')');

% Waitbar
axes('units','pixels','Position', [8 15 280 15], 'ytick',[],'xtick',[], 'box', 'on', 'parent', gammaGUIFig, 'color', bgColor);
uicontrol(gammaGUIFig, 'style', 'text', 'units', 'pixels', 'position', [8 35 40 12], 'string', 'Status', 'fontweight', 'bold');
axes('units','pixels','Position', [8 15 280 15], 'ytick',[],...
'xtick',[], 'box', 'on', 'parent', gammaGUIFig, 'color', bgColor);
uicontrol(gammaGUIFig, 'style', 'text', 'units', 'pixels', ...
'position', [8 35 40 12], 'string', 'Status', 'fontweight', 'bold');

%Waitbar axis, part of wb.
ud.wb.wbAxis = axes('units', 'pixels', 'Position', [8 15 280 15], 'color', [.9 .9 .9], 'ytick',[],'xtick',[], 'box', 'on', 'xlimmode', 'manual', 'ylimmode', 'manual', 'parent', gammaGUIFig);
ud.wb.patch = patch([0 0 0 0], [0 1 1 0], 'red', 'parent', ud.wb.wbAxis);
Expand All @@ -178,7 +211,6 @@ function CERRGammafnc(gamaCommand, varargin)
end

doseDiffIN = str2num(get(findobj('tag','InputDoseDiff'),'string'));

if isempty(doseDiffIN)
warndlg('Please Enter Dose Difference');
return;
Expand All @@ -189,11 +221,19 @@ function CERRGammafnc(gamaCommand, varargin)
return;
end
end

doseDiffMethod = get(findobj('tag','InputDoseDiffMethod'),'value');


DTA = str2num(get(findobj('tag','InputDTA'),'string'))/10;
if isempty(DTA)
warndlg('Please Enter Distance to Agreement (DTA)');
end

maxDistanceFactor = str2num(get(findobj('tag','MaxDistanceFactor'),'string'));
if isempty(DTA)
warndlg('Please Enter MaxDistance factor.');
end


threshold = str2num(get(findobj('tag','InputThreshold'),'string'));
Expand All @@ -219,7 +259,8 @@ function CERRGammafnc(gamaCommand, varargin)
% return;
% end

createGammaDose(baseDose,refDose,structNum,doseDiffIN,DTA,threshold);
createGammaDose(baseDose,refDose,structNum,doseDiffIN,DTA,...
threshold,maxDistanceFactor,doseDiffMethod);

CERRGammafnc('DISPCERRGAMMA', length(planC{indexS.dose}),structNum);

Expand Down
29 changes: 17 additions & 12 deletions CERR_core/Utilities/createGammaDose.m
View file
@@ -1,5 +1,7 @@
function gammaM = createGammaDose(doseNum1,doseNum2,strNum,dosePercent,distAgreement,thresholdPercentMax, planC)
% function gammaM = createGammaDose(doseNum1,doseNum2,dosePercent,distAgreement,thresholdPercentMax, planC)
function gammaM = createGammaDose(doseNum1,doseNum2,strNum,doseAgreement,...
distAgreement,thresholdPercentMax,maxDistanceFactor,doseDiffMethod,planC)
% function gammaM = createGammaDose(doseNum1,doseNum2,strNum,...
% doseAgreement,distAgreement,thresholdPercentMax,maxDistanceFactor,doseDiffMethod,planC);
%
% This function creates and adds a gamma dose to planC.
% 3D Gamma is calculated between doseNum1 and doseNum2.
Expand All @@ -16,21 +18,28 @@
indexS = planC{end};

% Prepare inputs for gamma calculation
%[newXgrid, newYgrid, newZgrid, doseArray1, doseArray2] = prepareDosesForGamma(doseNum1,doseNum2,1, planC);
assocScan = 1;
assocScan = getDoseAssociatedScan(doseNum1,planC);
if isempty(assocScan)
assocScan = 1;
end
[newXgrid, newYgrid, newZgrid, doseArray1, doseArray2, interpMask3M] = ...
prepareDosesForGamma(doseNum1,doseNum2,strNum,assocScan,planC);

deltaX = abs(newXgrid(2) - newXgrid(1));
deltaY = abs(newYgrid(2) - newYgrid(1));
deltaZ = abs(newZgrid(2) - newZgrid(1));

doseAgreement = dosePercent*max(planC{indexS.dose}(doseNum1).doseArray(:))/100;
if doseDiffMethod == 1
doseAgreement = doseAgreement*max(planC{indexS.dose}(doseNum1).doseArray(:))/100;
end

thresholdAbsolute = thresholdPercentMax*max(planC{indexS.dose}(doseNum1).doseArray(:))/100;

gammaM = gammaDose3d(doseArray1, doseArray2, interpMask3M, [deltaX deltaY deltaZ], doseAgreement, distAgreement, [], thresholdAbsolute);
%gammaM = gammaDose3d_new(doseArray1, doseArray2, [deltaX deltaY deltaZ], doseAgreement, distAgreement, [], thresholdAbsolute);
maxCalcDistance = distAgreement * maxDistanceFactor;

gammaM = gammaDose3d(doseArray1, doseArray2, interpMask3M, ...
[deltaX deltaY deltaZ], doseAgreement, distAgreement, ...
maxCalcDistance, thresholdAbsolute, doseDiffMethod);

% Assume doses within the filter threshold pass gamma
%gammaM(isnan(gammaM)) = 0;
Expand Down Expand Up @@ -58,12 +67,8 @@
planC{indexS.dose}(newDoseNum).doseUID = createUID('dose');
strName = planC{indexS.structures}(strNum).structureName;
planC{indexS.dose}(newDoseNum).fractionGroupID = ['Gamma_',...
num2str(dosePercent),'%_',num2str(distAgreement*10),'mm','_',strName];
num2str(doseAgreement),'%_',num2str(distAgreement*10),'mm','_',strName];

%Switch to new dose
sliceCallBack('selectDose', num2str(newDoseNum));


% stateS.doseSet = newDoseNum;
% stateS.doseChanged = 1;
% CERRRefresh
158 changes: 48 additions & 110 deletions CERR_core/Utilities/gammaDose3d.m
View file
@@ -1,7 +1,39 @@
function gammaM = gammaDose3d(doseArray1, doseArray2, strMask3M, deltaXYZv, doseAgreement, distAgreement, maxDistance, thresholdAbsolute)
% function gammaM = gammaDose3d(doseArray1, doseArray2, strMask3M, deltaXYZv, doseAgreement, distAgreement, maxDistance, thresholdAbsolute)
function gammaM = gammaDose3d(doseArray1, doseArray2, strMask3M, deltaXYZv,...
doseAgreement, distAgreement, maxSearchDistance, thresholdAbsolute, doseDiffMethod)
% function gammaM = gammaDose3d(doseArray1, doseArray2, strMask3M, deltaXYZv,...
% doseAgreement, distAgreement, maxSearchDistance, thresholdAbsolute, doseDiffMethod)
%
% Block-wise comparison for each voxel.
% This function returns the Gamma calculation for input dose distributions.
%
% Computation uses the following algorithm for speedup:
% "A fast algorithm for gamma evaluation in 3D”,
% Wendling et al, Medical physics, 34 (5), pp. 1647, 2007.
%
% which achieves speed and computer memory efficiency by searching around
% each reference point with increasing distance in a sphere,
% which has a radius of a chosen maximum search distance.
%
% doseArray1: reference dose
%
% doseArray2: evaluation dose
%
% strMask3M: binary mask. Gamma value is computed only for voxels that
% have a strMask3M value of 1
%
% deltaXYZv: Voxel size (dx,dy,dz) in cm
%
% doseDiffMethod: method to account for dose difference.
% 1: single dose value, usually percentage of max reference dose.
% 2: percentage of dose at each voxel.
%
% doseAgreement: when doseDiffMethod = 1, this is specified in Gy.
% when doseDiffMethod = 2, this is specified as a percentage.
%
% distAgreement: distance agreement in cm.
%
% maxSearchDistance: maximum radius of sphere as described in Wendling et al.
%
% thresholdAbsolute: Dose in Gy. for which gamma calculation is ignored.
%
% APA, 08/28/2015

Expand All @@ -10,13 +42,14 @@
deltaZ = deltaXYZv(3);

% Get Block size to process
if ~exist('maxDistance', 'var') || (exist('maxDistance', 'var') && isempty(maxDistance))
maxDistance = distAgreement*2;
if ~exist('maxSearchDistance', 'var') || ...
(exist('maxSearchDistance', 'var') && isempty(maxSearchDistance))
maxSearchDistance = distAgreement*2;
end
%maxDistance = 1; % cm
slcWindow = floor(2*maxDistance/deltaZ);
rowWindow = floor(2*maxDistance/deltaY);
colWindow = floor(2*maxDistance/deltaX);
%maxSearchDistance = 1; % cm
slcWindow = floor(2*maxSearchDistance/deltaZ);
rowWindow = floor(2*maxSearchDistance/deltaY);
colWindow = floor(2*maxSearchDistance/deltaX);

% Make sure that the window is of odd size
if mod(slcWindow,2) == 0
Expand Down Expand Up @@ -101,7 +134,12 @@
for slc = slcNum:slcWindow+slcNum-1
slc2M = doseArray2(:,:,slc);
tmpGammaV = bsxfun(@minus,slc2M(indSlcM),slc1M(calcSlcIndV)');
tmpGammaV = bsxfun(@plus,tmpGammaV.^2/doseAgreement^2 , (xysq + zV(slcCount)^2) / distAgreement^2);
if doseDiffMethod == 1
tmpGammaV = bsxfun(@plus,tmpGammaV.^2/doseAgreement^2 , (xysq + zV(slcCount)^2) / distAgreement^2);
elseif doseDiffMethod == 2
tmpGammaV = bsxfun(@rdivide,tmpGammaV.^2,(doseAgreement/100*slc1M(calcSlcIndV)').^2);
tmpGammaV = bsxfun(@plus, tmpGammaV, (xysq + zV(slcCount)^2) / distAgreement^2);
end
gammaV(calcSlcIndV) = min(gammaV(calcSlcIndV),min(tmpGammaV));
gammaM(:,:,slcNum) = reshape(gammaV,siz);
slcCount = slcCount + 1;
Expand All @@ -118,103 +156,3 @@
end




% % APA, 04/27/2012
%
% deltaX = deltaXYZv(1);
% deltaY = deltaXYZv(2);
% deltaZ = deltaXYZv(3);
% incrementRadius = min([deltaX deltaY deltaZ]);
%
% % Initial gamma
% gammaM = ((doseArray1-doseArray2).^2).^0.5/doseAgreement;
% %convergedM = false(size(gammaM));
%
% % Find regions of zero dose and exclude from calculation
% if ~exist('thresholdAbsolute', 'var')
% thresholdAbsolute = 0;
% end
% convergedM = doseArray1 <= thresholdAbsolute;
% %convergedM = doseArray1 >= thresholdAbsolute; % for MR scan
% gammaM(convergedM) = NaN;
%
% % Calculate until 4 times the permissible distance to agreement.
% if ~exist('maxDistance', 'var') || (exist('maxDistance', 'var') && isempty(maxDistance))
% maxDistance = distAgreement*4;
% end
% outerRadiusV = incrementRadius:incrementRadius:maxDistance;
% numIters = length(outerRadiusV);
%
% siz = size(gammaM);
% minDoseDiffM = zeros(siz,'single');
% maxDoseDiffM = minDoseDiffM;
%
% % convergenceCountM = zeros(siz,'uint8');
%
% % Update waitbar on gamma GUI
% gammaGUIFig = findobj('tag','CERRgammaInputGUI');
% if ~isempty(gammaGUIFig)
% ud = get(gammaGUIFig,'userdata');
% set(ud.wb.patch,'xData',[0 0 0 0])
% end
%
% for radNum = 1:numIters
%
% disp(['--- Gamma Calculation Iteraion ', num2str(radNum), ' ----'])
%
% % Create an ellipsoid ring neighborhood
% outerRadius = outerRadiusV(radNum);
%
% rowOutV = floor(outerRadius/deltaY);
% colOutV = floor(outerRadius/deltaX);
% slcOutV = floor(outerRadius/deltaZ);
%
% NHOOD_outer = createEllipsoidNHOOD(1:rowOutV,1:colOutV,1:slcOutV);
%
% % Compute Min and Max for the ellipsoid neighborhood
% [minLocalM, maxLocalM] = getMinMaxIM(doseArray2,NHOOD_outer);
%
% % Compute difference between dose1 and (min,max) for dose2
% minDoseDiffM(~convergedM) = doseArray1(~convergedM) - minLocalM(~convergedM);
% maxDoseDiffM(~convergedM) = maxLocalM(~convergedM) - doseArray1(~convergedM);
%
% % If dose1 is contained within (min,max) for dose2, then it converged!
% newConvergedM = ~convergedM & minDoseDiffM >= 0 & maxDoseDiffM >= 0;
%
% % Compute gamma for voxels that converged
% gammaM(newConvergedM) = min(gammaM(newConvergedM), outerRadius/distAgreement);
%
% % Add newly converged voxels to the list
% convergedM = convergedM | newConvergedM;
%
% % Compute gamma for voxels that have not yet converged
% gammaForNotConvergedM = (min((minDoseDiffM(~convergedM)-doseAgreement).^2, (maxDoseDiffM(~convergedM)-doseAgreement).^2)./doseAgreement^2 + (outerRadius/distAgreement)^2).^0.5;
% gammaM(~convergedM) = min(gammaM(~convergedM), gammaForNotConvergedM);
%
% % If gamma is less or equal to outerRadius/distAgreement, then the voxel conveged!
% convergedM(~convergedM) = gammaM(~convergedM) <= outerRadius/distAgreement;
%
% % % Count number of consecutive gamma increments for each uncoverged voxel
% % indConvergeM = false(size(gammaM));
% % indConvergeM(~convergedM) = gammaM(~convergedM) <= gammaForNotConvergedM;
% % convergenceCountM(indConvergeM) = convergenceCountM(indConvergeM) + 1;
% %
% % % If gamma increases 4 consecutive times for a voxel, then assume it has converged
% % convergedM(convergenceCountM > 3) = 1;
%
% if ~isempty(gammaGUIFig)
% set(ud.wb.patch,'xData',[0 0 radNum/numIters radNum/numIters])
% drawnow
% end
%
% if all(convergedM(:))
% disp('All converged!!!')
% if ~isempty(gammaGUIFig)
% set(ud.wb.patch,'xData',[0 0 1 1])
% end
% break
% end
%
% end

2 changes: 1 addition & 1 deletion CERR_core/Utilities/interpStructMaskToDose.m
View file
Expand Up @@ -41,7 +41,7 @@
kMax = min(find(newZgrid > zV(max(kV))));
inputTM = eye(4);
for slcNum=kMin:kMax %1:length(newZgrid)
disp(['Interpolating slice ', num2str(kMin), '/', num2str(kMin)])
disp(['Interpolating slice ', num2str(slcNum), '/', num2str(kMax)])
strTmpM = slice3DVol(mask3M, xV, yV, zV, newZgrid(slcNum), 3, 'linear', inputTM, [], newXgrid, newYgrid);
if ~isempty(strTmpM)
interpMask3M(:,:,slcNum) = strTmpM > 0.5;
Expand Down

0 comments on commit 1643aec

Please sign in to comment.