function powerDVH(lin245cm,root245cm,quad245cm,power025,power4)

addpath('C:\Documents and Settings\Steve\My Documents\Research\matlab')
addpath('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqUtil')
addpath('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input')

Ndvhbins = 1000;                                                            %number of dvh bins

prescTumorLin = open_presc('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input\HN003powerLaw\normPower1\PTV70dosePlus.bin');                            %open prescription in Matlab
prescTumorRoot = open_presc('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input\HN003powerLaw\normPower2\PTV70dosePlus.bin');
prescTumorQuad = open_presc('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input\HN003powerLaw\normPower3\PTV70dosePlus.bin');
prescTumorPower025 = open_presc('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input\HN003powerLaw\normPower4\PTV70dosePlus.bin');
prescTumorPower4 = open_presc('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input\HN003powerLaw\normPower5\PTV70dosePlus.bin');

%dim = 1;                                                                    %2D
dim = 30;                                                                   %3D

fid = fopen(lin245cm,'rb');                                           %open plans in Matlab
planLin245cm = reshape(fread(fid,'float'),72,72,dim);
fclose(fid);

fid = fopen(root245cm,'rb');                                           %open plans in Matlab
planRoot245cm = reshape(fread(fid,'float'),72,72,dim);
fclose(fid);

fid = fopen(quad245cm,'rb');                                           %open plans in Matlab
planQuad245cm = reshape(fread(fid,'float'),72,72,dim);
fclose(fid);

fid = fopen(power025,'rb');                                           %open plans in Matlab
planPower025 = reshape(fread(fid,'float'),72,72,dim);
fclose(fid);

fid = fopen(power4,'rb');                                           %open plans in Matlab
planPower4 = reshape(fread(fid,'float'),72,72,dim);
fclose(fid);


tumorMaskLin = single(~~(full3D(prescTumorLin)));                                 %define tissue masks
tumorMaskRoot = single(~~(full3D(prescTumorRoot)));
tumorMaskQuad = single(~~(full3D(prescTumorQuad)));
tumorMaskPower025 = single(~~(full3D(prescTumorPower025)));
tumorMaskPower4 = single(~~(full3D(prescTumorPower4)));


tumorDoseLin245cm = planLin245cm.*tumorMaskLin;                                                %dose distributions in each ROI
tumorDoseRoot245cm = planRoot245cm.*tumorMaskRoot;
tumorDoseQuad245cm = planQuad245cm.*tumorMaskQuad;
tumorDosePower025 = planPower025.*tumorMaskPower025;
tumorDosePower4 = planPower4.*tumorMaskPower4;


%convert 3D matrices to sparse matrices
tumorDoseSparseLin245cm = sparse3D(tumorDoseLin245cm);
tumorDoseSparseRoot245cm = sparse3D(tumorDoseRoot245cm);
tumorDoseSparseQuad245cm = sparse3D(tumorDoseQuad245cm);
tumorDoseSparsePower025 = sparse3D(tumorDosePower025);
tumorDoseSparsePower4 = sparse3D(tumorDosePower4);


tumorNormDoseLin245cm = ((tumorDoseSparseLin245cm.non_zero_values)./(prescTumorLin.non_zero_values));  %calculate equivalent uniform dose
tumorAverageDoseLin245cm = mean(prescTumorLin.non_zero_values);
tumorEdoseLin245cm = tumorNormDoseLin245cm.*tumorAverageDoseLin245cm;
tumorEdoseLin245cm = tumorNormDoseLin245cm;
indLin = find(((tumorNormDoseLin245cm-1.).^2).^(1/2) >= 0.05);
volLin = 100.-double(length(indLin))./(double(tumorDoseSparseLin245cm.Nind))*100

tumorNormDoseRoot245cm = ((tumorDoseSparseRoot245cm.non_zero_values)./(prescTumorRoot.non_zero_values));  %calculate equivalent uniform dose
tumorAverageDoseRoot245cm = mean(prescTumorRoot.non_zero_values);
tumorEdoseRoot245cm = tumorNormDoseRoot245cm.*tumorAverageDoseRoot245cm;
tumorEdoseRoot245cm = tumorNormDoseRoot245cm;
indRoot = find(((tumorNormDoseRoot245cm-1.).^2).^(1/2) >= 0.05);
volRoot = 100.-double(length(indRoot))./(double(tumorDoseSparseRoot245cm.Nind))*100

tumorNormDoseQuad245cm = ((tumorDoseSparseQuad245cm.non_zero_values)./(prescTumorQuad.non_zero_values));  %calculate equivalent uniform dose
tumorAverageDoseQuad245cm = mean(prescTumorQuad.non_zero_values);
tumorEdoseQuad245cm = tumorNormDoseQuad245cm.*tumorAverageDoseQuad245cm;
tumorEdoseQuad245cm = tumorNormDoseQuad245cm;
indQuad = find(((tumorNormDoseQuad245cm-1.).^2).^(1/2) >= 0.05);
volQuad = 100.-double(length(indQuad))./(double(tumorDoseSparseQuad245cm.Nind))*100

tumorNormDosePower025 = ((tumorDoseSparsePower025.non_zero_values)./(prescTumorPower025.non_zero_values));  %calculate equivalent uniform dose
tumorAverageDosePower025 = mean(prescTumorPower025.non_zero_values);
tumorEdosePower025 = tumorNormDosePower025.*tumorAverageDosePower025;
tumorEdosePower025 = tumorNormDosePower025;
indPower025 = find(((tumorNormDosePower025-1.).^2).^(1/2) >= 0.15);
min(((tumorNormDosePower025-1.).^2).^(1/2))
volPower025 = 100.-double(length(indPower025))./(double(tumorDoseSparsePower025.Nind))*100

tumorNormDosePower4 = ((tumorDoseSparsePower4.non_zero_values)./(prescTumorPower4.non_zero_values));  %calculate equivalent uniform dose
tumorAverageDosePower4 = mean(prescTumorPower4.non_zero_values);
tumorEdosePower4 = tumorNormDosePower4.*tumorAverageDosePower4;
tumorEdosePower4 = tumorNormDosePower4;
indPower4 = find(((tumorNormDosePower4-1.).^2).^(1/2) >= 0.15);
volPower4 = 100.-double(length(indPower4))./(double(tumorDoseSparsePower4.Nind))*100

% 
dMax = 1.1*max(tumorEdoseLin245cm);
dvhbins = [0:Ndvhbins-1]*dMax/Ndvhbins;

dHistTumorLin245cm = histc(tumorEdoseLin245cm,dvhbins);
dHistTumorRoot245cm = histc(tumorEdoseRoot245cm,dvhbins);
dHistTumorQuad245cm = histc(tumorEdoseQuad245cm,dvhbins);
dHistTumorPower025 = histc(tumorEdosePower025,dvhbins);
dHistTumorPower4 = histc(tumorEdosePower4,dvhbins);

% dHistTumorGomp245cm = histc(tumorEdoseGomp245cm,dvhbins);
% 
% 
% dHistTumorUni245cm = histc(tumorDoseSparseUni245cm.non_zero_values,dvhbins);
% dHistTumorLin = histc(tumorDoseSparseLin245cm.non_zero_values,dvhbins);
% dHistTumorRoot = histc(tumorDoseSparseRoot245cm.non_zero_values,dvhbins);
% dHistTumorQuad = histc(tumorDoseSparseQuad245cm.non_zero_values,dvhbins);
% dHistTumorGomp = histc(tumorDoseSparseGomp245cm.non_zero_values,dvhbins);
% 
% 
% 
% dHistTumorNormUni245cm = 100*dHistTumorUni245cm./(double(tumorDoseSparseUni245cm.Nind));
dHistTumorNormLin245cm = 100*dHistTumorLin245cm./(double(tumorDoseSparseLin245cm.Nind));
dHistTumorNormRoot245cm = 100*dHistTumorRoot245cm./(double(tumorDoseSparseRoot245cm.Nind));
dHistTumorNormQuad245cm = 100*dHistTumorQuad245cm./(double(tumorDoseSparseQuad245cm.Nind));
dHistTumorNormPower025 = 100*dHistTumorPower025./(double(tumorDoseSparsePower025.Nind));
dHistTumorNormPower4 = 100*dHistTumorPower4./(double(tumorDoseSparsePower4.Nind));

% dHistTumorNormGomp245cm = 100*dHistTumorGomp245cm./(double(tumorDoseSparseGomp245cm.Nind));
% 
% dHistTumorNormLin = 100*dHistTumorLin./(double(tumorDoseSparseLin245cm.Nind));
% dHistTumorNormRoot = 100*dHistTumorRoot./(double(tumorDoseSparseRoot245cm.Nind));
% dHistTumorNormQuad = 100*dHistTumorQuad./(double(tumorDoseSparseQuad245cm.Nind));
% dHistTumorNormGomp = 100*dHistTumorGomp/(prescTumorGomp.Nind);
% 
% 
% DVHtumorUni245cm = flipud(cumsum(flipud(dHistTumorNormUni245cm)));
DVHtumorLin245cm = flipud(cumsum(flipud(dHistTumorNormLin245cm)));
DVHtumorRoot245cm = flipud(cumsum(flipud(dHistTumorNormRoot245cm)));
DVHtumorQuad245cm = flipud(cumsum(flipud(dHistTumorNormQuad245cm)));
DVHtumorPower025 = flipud(cumsum(flipud(dHistTumorNormPower025)));
DVHtumorPower4 = flipud(cumsum(flipud(dHistTumorNormPower4)));

% DVHtumorGomp245cm = flipud(cumsum(flipud(dHistTumorNormGomp245cm)));
% 
% DVHtumorLin = flipud(cumsum(flipud(dHistTumorNormLin)));
% DVHtumorRoot = flipud(cumsum(flipud(dHistTumorNormRoot)));
% DVHtumorQuad = flipud(cumsum(flipud(dHistTumorNormQuad)));
% DVHtumorGomp = flipud(cumsum(flipud(dHistTumorNormGomp)));



normalfigure(dvhbins,DVHtumorLin245cm,DVHtumorRoot245cm,DVHtumorQuad245cm,DVHtumorPower025,DVHtumorPower4)

% DVHrtLin245cm,DVHrtRoot245cm,DVHrtQuad245cm,...
%     DVHcordLin245cm,DVHcordRoot245cm,DVHcordQuad245cm,...
%     DVHlarynxLin245cm,DVHlarynxRoot245cm,DVHlarynxQuad245cm)
end