WiscPlan_olderMV/opt/linlsqUtil/tumorDVH.m

function tumorDVH(lin0625cm,root0625cm,quad0625cm)

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\HN003rootPrescVol\PTV70DosePlus.bin');
prescTumorQuad = open_presc('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input\HN003powerLaw\HN003quadPrescVol\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(lin0625cm,'rb');                                           %open plans in Matlab
planLin0625cm = reshape(fread(fid,'float'),72,72,dim);
fclose(fid);

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

fid = fopen(quad0625cm,'rb');                                           %open plans in Matlab
planQuad0625cm = 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)));


% tumorDoseUni245cm = planUni245cm.*tumorMaskUni;
tumorDoseLin0625cm = planLin0625cm.*tumorMaskLin;                                                %dose distributions in each ROI
tumorDoseRoot0625cm = planRoot0625cm.*tumorMaskRoot;
tumorDoseQuad0625cm = planQuad0625cm.*tumorMaskQuad;


% tumorDoseSparseUni245cm = sparse3D(tumorDoseUni245cm);                                      %convert 3D matrices to sparse matrices
tumorDoseSparseLin0625cm = sparse3D(tumorDoseLin0625cm);
tumorDoseSparseRoot0625cm = sparse3D(tumorDoseRoot0625cm);
tumorDoseSparseQuad0625cm = sparse3D(tumorDoseQuad0625cm);


tumorNormDoseLin0625cm = ((tumorDoseSparseLin0625cm.non_zero_values)./(prescTumorLin.non_zero_values));  %calculate equivalent uniform dose
tumorAverageDoseLin0625cm = mean(prescTumorLin.non_zero_values);
tumorEdoseLin0625cm = tumorNormDoseLin0625cm.*tumorAverageDoseLin0625cm;

tumorNormDoseRoot0625cm = ((tumorDoseSparseRoot0625cm.non_zero_values)./(prescTumorRoot.non_zero_values));  %calculate equivalent uniform dose
tumorAverageDoseRoot0625cm = mean(prescTumorRoot.non_zero_values);
tumorEdoseRoot0625cm = tumorNormDoseRoot0625cm.*tumorAverageDoseRoot0625cm;

tumorNormDoseQuad0625cm = ((tumorDoseSparseQuad0625cm.non_zero_values)./(prescTumorQuad.non_zero_values));  %calculate equivalent uniform dose
tumorAverageDoseQuad0625cm = mean(prescTumorQuad.non_zero_values);
tumorEdoseQuad0625cm = tumorNormDoseQuad0625cm.*tumorAverageDoseQuad0625cm;

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

dHistTumorLin0625cm = histc(tumorEdoseLin0625cm,dvhbins);
dHistTumorRoot0625cm = histc(tumorEdoseRoot0625cm,dvhbins);
dHistTumorQuad0625cm = histc(tumorEdoseQuad0625cm,dvhbins);


% dHistTumorNormUni245cm = 100*dHistTumorUni245cm./(double(tumorDoseSparseUni245cm.Nind));
dHistTumorNormLin0625cm = 100*dHistTumorLin0625cm./(double(tumorDoseSparseLin0625cm.Nind));
dHistTumorNormRoot0625cm = 100*dHistTumorRoot0625cm./(double(tumorDoseSparseRoot0625cm.Nind));
dHistTumorNormQuad0625cm = 100*dHistTumorQuad0625cm./(double(tumorDoseSparseQuad0625cm.Nind));


% DVHtumorUni245cm = flipud(cumsum(flipud(dHistTumorNormUni245cm)));
DVHtumorLin0625cm = flipud(cumsum(flipud(dHistTumorNormLin0625cm)));
DVHtumorRoot0625cm = flipud(cumsum(flipud(dHistTumorNormRoot0625cm)));
DVHtumorQuad0625cm = flipud(cumsum(flipud(dHistTumorNormQuad0625cm)));



tumorfigure(dvhbins,DVHtumorLin0625cm,DVHtumorRoot0625cm,DVHtumorQuad0625cm)

end