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')

prescFile = 'presc.img';                                                            %locate prescription and plan filenames
planFile1 = 'plan1.img';
planFile2 = 'plan2.img';
planFile3 = 'plan3.img';

dim1 = 512;
dim2 = dim1;
dim3 = 35;                                                                          %3D matrix dimensions

fid = fopen(prescFile,'rb');                                                        %open plans in Matlab
presc = reshape(fread(fid,'float'),dim1,dim2,dim3);
fclose(fid);

fid = fopen(planFile1,'rb');                                                        %open plans in Matlab
plan1 = reshape(fread(fid,'float'),dim1,dim2,dim3);
fclose(fid);

fid = fopen(planFile2,'rb');                                                        %open plans in Matlab
plan2 = reshape(fread(fid,'float'),dim1,dim2,dim3);
fclose(fid);

fid = fopen(prescFile3,'rb');                                                       %open plans in Matlab
plan3 = reshape(fread(fid,'float'),dim1,dim2,dim3);
fclose(fid);

PTVmask = single(~~(presc));                                                        %define tissue mask

PTVdose1 = plan1.*tumorMask;                                                        %dose distributions in each ROI
PTVdose2 = plan2.*tumorMask;                                                        %dose distributions in each ROI
PTVdose3 = plan3.*tumorMask;                                                        %dose distributions in each ROI

PTVdoseSparse1 = sparse3D(PTVdose1);
PTVdoseSparse2 = sparse3D(PTVdose2);
PTVdoseSparse3 = sparse3D(PTVdose3);

PTVnormDose1 = ((PTVdoseSparse1.non_zero_values)./(PTVpresc.non_zero_values));      %calculate normalized dose
ind1 = find(((PTVnormDose1-1.).^2).^(1/2) >= 0.05);                                 %find voxels whose plans deviate from prescription > 5 percent
vol1 = 100.-double(length(ind1))./(double(PTVnormDose1.Nind))*100;                  %calculate PTV percentage whose plan is within 5 percent prescription

PTVnormDose2 = ((PTVdoseSparse2.non_zero_values)./(PTVpresc.non_zero_values));    
ind2 = find(((PTVnormDose2-1.).^2).^(1/2) >= 0.05);                                
vol2 = 100.-double(length(ind2))./(double(PTVnormDose2.Nind))*100;                  

PTVnormDose3 = ((PTVdoseSparse3.non_zero_values)./(PTVpresc.non_zero_values));    
ind3 = find(((PTVnormDose3-1.).^2).^(1/2) >= 0.05);                                
vol3 = 100.-double(length(ind3))./(double(PTVnormDose3.Nind))*100;                  

Ndvhbins = 1000;                                                                    %number of dvh bins

dMax = 1.1*max(PTVnormDose1);                                                       %define maximum bin for QVH
dvhbins = [0:Ndvhbins-1]*dMax/Ndvhbins;                                             %define bins

dHistPTVnormDose1 = histc(PTVnormDose1,dvhbins);                                    %bin normalized PTV dose 
dHistPTVnormDose1 = 100*dHistPTVnormDose1./(double(PTVnormDose1.Nind));             %calculate frequency percentage of integral

dHistPTVnormDose2 = histc(PTVnormDose2,dvhbins);                                    %bin normalized PTV dose 
dHistPTVnormDose2 = 100*dHistPTVnormDose2./(double(PTVnormDose2.Nind));             %calculate frequency percentage of integral

dHistPTVnormDose3 = histc(PTVnormDose3,dvhbins);                                    %bin normalized PTV dose 
dHistPTVnormDose3 = 100*dHistPTVnormDose3./(double(PTVnormDose3.Nind));             %calculate frequency percentage of integral


QVHptv1 = flipud(cumsum(flipud(dHistPTVnormDose1)));                                %cumulatively sum histogram from max min to min bin
QVHptv2 = flipud(cumsum(flipud(dHistPTVnormDose2)));                                %cumulatively sum histogram from max min to min bin
QVHptv3 = flipud(cumsum(flipud(dHistPTVnormDose3)));                                %cumulatively sum histogram from max min to min bin

plotQVH(dvhbins,QVHptv1, QVHptv2, QVHptv3)                                                         %plot QVH