WiscPlan_olderMV/opt/linlsqUtil/normalDVH.asv

function normalDVH(uni245cm,lin245cm,root245cm,quad245cm,gomp245cm)

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

prescTumorUni = open_presc('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input\HN003uniPrescVolNew\PTV70dosePlus.bin');
prescTumorLin = open_presc('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input\HN003powerLaw\HN003linPrescVol\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\HN003quadPrescVol\PTV70dosePlus.bin');
prescTumorGomp = open_presc('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input\sigmoidLaw\HN003gompPrescVol\PTV70dosePlus.bin');
prescNormal = open_presc('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input\HN003uniPrescVolNew\NormalDosePlus.bin');
prescLarynx = open_presc('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input\HN003uniPrescVolNew\larynxDosePlus.bin');
prescRt = open_presc('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input\HN003uniPrescVolNew\RtparotidDosePlus.bin');
prescLt = open_presc('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input\HN003uniPrescVolNew\LtparotidDosePlus.bin');
prescCord = open_presc('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input\HN003uniPrescVolNew\CordDosePlus.bin');

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

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

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(gomp245cm,'rb');                                           %open plans in Matlab
planGomp245cm = reshape(fread(fid,'float'),72,72,dim);
fclose(fid);

for i=1:prescNormal.Nind;
    prescNormal.non_zero_values(i)=1;
end

for i=1:prescLarynx.Nind;
    prescLarynx.non_zero_values(i)=1;
end

for i=1:prescRt.Nind;
    prescRt.non_zero_values(i)=1;
end

for i=1:prescLt.Nind;
    prescLt.non_zero_values(i)=1;
end

for i=1:prescCord.Nind;
    prescCord.non_zero_values(i)=1;
end

tumorMaskUni = single(~~(full3D(prescTumorUni)));
tumorMaskLin = single(~~(full3D(prescTumorLin)));                                 %define tissue masks
tumorMaskRoot = single(~~(full3D(prescTumorRoot)));
tumorMaskQuad = single(~~(full3D(prescTumorQuad)));
tumorMaskGomp = single(~~(full3D(prescTumorGomp)));
larynxMask = single(~~(full3D(prescLarynx)));
rtMask = single(~~(full3D(prescRt)));
cordMask = single(~~(full3D(prescCord)));
normalMask = single(~~(full3D(prescNormal)));

tumorDoseUni245cm = planUni245cm.*tumorMaskUni;
tumorDoseLin245cm = planLin245cm.*tumorMaskLin;                                                %dose distributions in each ROI
tumorDoseRoot245cm = planRoot245cm.*tumorMaskRoot;
tumorDoseQuad245cm = planQuad245cm.*tumorMaskQuad;
tumorDoseGomp245cm = planGomp245cm.*tumorMaskGomp;


larynxDoseUni245cm = planUni245cm.*larynxMask;
larynxDoseLin245cm = planLin245cm.*larynxMask;
larynxDoseRoot245cm = planRoot245cm.*larynxMask;
larynxDoseQuad245cm = planQuad245cm.*larynxMask;
larynxDoseGomp245cm = planGomp245cm.*larynxMask;

rtDoseUni245cm = planUni245cm.*rtMask;
rtDoseLin245cm = planLin245cm.*rtMask;
rtDoseRoot245cm = planRoot245cm.*rtMask;
rtDoseQuad245cm = planQuad245cm.*rtMask;
rtDoseGomp245cm = planGomp245cm.*rtMask;

cordDoseUni245cm = planUni245cm.*cordMask;
cordDoseLin245cm = planLin245cm.*cordMask;
cordDoseRoot245cm = planRoot245cm.*cordMask;
cordDoseQuad245cm = planQuad245cm.*cordMask;
cordDoseGomp245cm = planGomp245cm.*cordMask;

normalDoseUni245cm = planUni245cm.*normalMask;
normalDoseLin245cm = planLin245cm.*normalMask;
normalDoseRoot245cm = planRoot245cm.*normalMask;
normalDoseQuad245cm = planQuad245cm.*normalMask;
normalDoseGomp245cm = planGomp245cm.*normalMask;


tumorDoseSparseUni245cm = sparse3D(tumorDoseUni245cm);                                      %convert 3D matrices to sparse matrices
tumorDoseSparseLin245cm = sparse3D(tumorDoseLin245cm);
tumorDoseSparseRoot245cm = sparse3D(tumorDoseRoot245cm);
tumorDoseSparseQuad245cm = sparse3D(tumorDoseQuad245cm);
tumorDoseSparseGomp245cm = sparse3D(tumorDoseGomp245cm);

larynxDoseSparseUni245cm = sparse3D(larynxDoseUni245cm);
larynxDoseSparseLin245cm = sparse3D(larynxDoseLin245cm);
larynxDoseSparseRoot245cm = sparse3D(larynxDoseRoot245cm);
larynxDoseSparseQuad245cm = sparse3D(larynxDoseQuad245cm);
larynxDoseSparseGomp245cm = sparse3D(larynxDoseGomp245cm);

rtDoseSparseUni245cm = sparse3D(rtDoseUni245cm);
rtDoseSparseLin245cm = sparse3D(rtDoseLin245cm);
rtDoseSparseRoot245cm = sparse3D(rtDoseRoot245cm);
rtDoseSparseQuad245cm = sparse3D(rtDoseQuad245cm);
rtDoseSparseGomp245cm = sparse3D(rtDoseGomp245cm);

cordDoseSparseUni245cm = sparse3D(cordDoseUni245cm);
cordDoseSparseLin245cm = sparse3D(cordDoseLin245cm);
cordDoseSparseRoot245cm = sparse3D(cordDoseRoot245cm);
cordDoseSparseQuad245cm = sparse3D(cordDoseQuad245cm);
cordDoseSparseGomp245cm = sparse3D(cordDoseGomp245cm);

normalDoseSparseUni245cm = sparse3D(normalDoseUni245cm);
normalDoseSparseLin245cm = sparse3D(normalDoseLin245cm);
normalDoseSparseRoot245cm = sparse3D(normalDoseRoot245cm);
normalDoseSparseQuad245cm = sparse3D(normalDoseQuad245cm);
normalDoseSparseGomp245cm = sparse3D(normalDoseGomp245cm);

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

tumorNormDoseGomp245cm = ((tumorDoseSparseGomp245cm.non_zero_values)./(prescTumorGomp.non_zero_values));  %calculate equivalent uniform dose
tumorAverageDoseGomp245cm = mean(prescTumorGomp.non_zero_values);
tumorEdoseGomp245cm = tumorNormDoseGomp245cm.*tumorAverageDoseGomp245cm;

dMax = 1.1*max(tumorEdoseLin245cm);

dvhbins = [0:Ndvhbins-1]*dMax./Ndvhbins;

dHistTumorLin245cm = histc(tumorEdoseLin245cm,dvhbins);
dHistTumorRoot245cm = histc(tumorEdoseRoot245cm,dvhbins);
dHistTumorQuad245cm = histc(tumorEdoseQuad245cm,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);

dHistLarynxUni245cm = histc(larynxDoseSparseUni245cm.non_zero_values,dvhbins);
dHistLarynxLin245cm = histc(larynxDoseSparseLin245cm.non_zero_values,dvhbins);
dHistLarynxRoot245cm = histc(larynxDoseSparseRoot245cm.non_zero_values,dvhbins);
dHistLarynxQuad245cm = histc(larynxDoseSparseQuad245cm.non_zero_values,dvhbins);
dHistLarynxGomp245cm = histc(larynxDoseSparseGomp245cm.non_zero_values,dvhbins);

dHistRtUni245cm = histc(rtDoseSparseUni245cm.non_zero_values,dvhbins);
dHistRtLin245cm = histc(rtDoseSparseLin245cm.non_zero_values,dvhbins);
dHistRtRoot245cm = histc(rtDoseSparseRoot245cm.non_zero_values,dvhbins);
dHistRtQuad245cm = histc(rtDoseSparseQuad245cm.non_zero_values,dvhbins);
dHistRtGomp245cm = histc(rtDoseSparseGomp245cm.non_zero_values,dvhbins);

dHistCordUni245cm = histc(cordDoseSparseUni245cm.non_zero_values,dvhbins);
dHistCordLin245cm = histc(cordDoseSparseLin245cm.non_zero_values,dvhbins);
dHistCordRoot245cm = histc(cordDoseSparseRoot245cm.non_zero_values,dvhbins);
dHistCordQuad245cm = histc(cordDoseSparseQuad245cm.non_zero_values,dvhbins);
dHistCordGomp245cm = histc(cordDoseSparseGomp245cm.non_zero_values,dvhbins);

dHistNormalUni245cm = histc(normalDoseSparseUni245cm.non_zero_values,dvhbins);
dHistNormalLin245cm = histc(normalDoseSparseLin245cm.non_zero_values,dvhbins);
dHistNormalRoot245cm = histc(normalDoseSparseRoot245cm.non_zero_values,dvhbins);
dHistNormalQuad245cm = histc(normalDoseSparseQuad245cm.non_zero_values,dvhbins);
dHistNormalGomp245cm = histc(normalDoseSparseGomp245cm.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));
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);

dHistLarynxNormUni245cm = 100*dHistLarynxUni245cm./(double(larynxDoseSparseUni245cm.Nind));
dHistLarynxNormLin245cm = 100*dHistLarynxLin245cm./(double(larynxDoseSparseLin245cm.Nind));
dHistLarynxNormRoot245cm = 100*dHistLarynxRoot245cm./(double(larynxDoseSparseRoot245cm.Nind));
dHistLarynxNormQuad245cm = 100*dHistLarynxQuad245cm./(double(larynxDoseSparseQuad245cm.Nind));
dHistLarynxNormGomp245cm = 100*dHistLarynxGomp245cm./(double(larynxDoseSparseGomp245cm.Nind));

dHistRtNormUni245cm = 100*dHistRtUni245cm./(double(rtDoseSparseUni245cm.Nind));
dHistRtNormLin245cm = 100*dHistRtLin245cm./(double(rtDoseSparseLin245cm.Nind));
dHistRtNormRoot245cm = 100*dHistRtRoot245cm./(double(rtDoseSparseRoot245cm.Nind));
dHistRtNormQuad245cm = 100*dHistRtQuad245cm./(double(rtDoseSparseQuad245cm.Nind));
dHistRtNormGomp245cm = 100*dHistRtGomp245cm./(double(rtDoseSparseGomp245cm.Nind));

dHistCordNormUni245cm = 100*dHistCordUni245cm./(double(cordDoseSparseUni245cm.Nind));
dHistCordNormLin245cm = 100*dHistCordLin245cm./(double(cordDoseSparseLin245cm.Nind));
dHistCordNormRoot245cm = 100*dHistCordRoot245cm./(double(cordDoseSparseRoot245cm.Nind));
dHistCordNormQuad245cm = 100*dHistCordQuad245cm./(double(cordDoseSparseQuad245cm.Nind));
dHistCordNormGomp245cm = 100*dHistCordGomp245cm./(double(cordDoseSparseGomp245cm.Nind));

dHistNormalNormUni245cm = 100*dHistNormalUni245cm./(double(normalDoseSparseUni245cm.Nind));
dHistNormalNormLin245cm = 100*dHistNormalLin245cm./(double(normalDoseSparseLin245cm.Nind));
dHistNormalNormRoot245cm = 100*dHistNormalRoot245cm./(double(normalDoseSparseRoot245cm.Nind));
dHistNormalNormQuad245cm = 100*dHistNormalQuad245cm./(double(normalDoseSparseQuad245cm.Nind));
dHistNormalNormGomp245cm = 100*dHistNormalGomp245cm./(double(normalDoseSparseGomp245cm.Nind));


DVHtumorUni245cm = flipud(cumsum(flipud(dHistTumorNormUni245cm)));
DVHtumorLin245cm = flipud(cumsum(flipud(dHistTumorNormLin245cm)));
DVHtumorRoot245cm = flipud(cumsum(flipud(dHistTumorNormRoot245cm)));
DVHtumorQuad245cm = flipud(cumsum(flipud(dHistTumorNormQuad245cm)));
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)));

DVHlarynxUni245cm = flipud(cumsum(flipud(dHistLarynxNormUni245cm)));
DVHlarynxLin245cm = flipud(cumsum(flipud(dHistLarynxNormLin245cm)));
DVHlarynxRoot245cm = flipud(cumsum(flipud(dHistLarynxNormRoot245cm)));
DVHlarynxQuad245cm = flipud(cumsum(flipud(dHistLarynxNormQuad245cm)));
DVHlarynxGomp245cm = flipud(cumsum(flipud(dHistLarynxNormGomp245cm)));


DVHrtUni245cm = flipud(cumsum(flipud(dHistRtNormUni245cm)));
DVHrtLin245cm = flipud(cumsum(flipud(dHistRtNormLin245cm)));
DVHrtRoot245cm = flipud(cumsum(flipud(dHistRtNormRoot245cm)));
DVHrtQuad245cm = flipud(cumsum(flipud(dHistRtNormQuad245cm)));
DVHrtGomp245cm = flipud(cumsum(flipud(dHistRtNormGomp245cm)));

DVHcordUni245cm = flipud(cumsum(flipud(dHistCordNormUni245cm)));
DVHcordLin245cm = flipud(cumsum(flipud(dHistCordNormLin245cm)));
DVHcordRoot245cm = flipud(cumsum(flipud(dHistCordNormRoot245cm)));
DVHcordQuad245cm = flipud(cumsum(flipud(dHistCordNormQuad245cm)));
DVHcordGomp245cm = flipud(cumsum(flipud(dHistCordNormGomp245cm)));

DVHnormalUni245cm = flipud(cumsum(flipud(dHistNormalNormUni245cm)));
DVHnormalLin245cm = flipud(cumsum(flipud(dHistNormalNormLin245cm)));
DVHnormalRoot245cm = flipud(cumsum(flipud(dHistNormalNormRoot245cm)));
DVHnormalQuad245cm = flipud(cumsum(flipud(dHistNormalNormQuad245cm)));
DVHnormalGomp245cm = flipud(cumsum(flipud(dHistNormalNormGomp245cm)));


figNormal(dvhbins,DVHrtLin245cm,DVHrtRoot245cm,DVHrtQuad245cm,...
     DVHcordLin245cm,DVHcordRoot245cm,DVHcordQuad245cm,...
     DVHlarynxLin245cm,DVHlarynxRoot245cm,DVHlarynxQuad245cm,...
     DVHnormalLin245cm,DVHnormalRoot245cm,DVHnormalQuad245cm)


end