WiscPlan_olderMV/opt/linlsqUtil/QVH.asv

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

prescFile = 'presc.img';
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(~~(full3D(presc)));                                            %define tissue masks

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 equivalent uniform 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));    %calculate equivalent uniform dose
ind2 = find(((PTVnormDose2-1.).^2).^(1/2) >= 0.05);                                %find voxels whose plans deviate from prescription > 5 percent
vol2 = 100.-double(length(ind2))./(double(PTVnormDose2.Nind))*100;                  %calculate PTV percentage whose plan is within 5 percent prescription

PTVnormDose3 = ((PTVdoseSparse3.non_zero_values)./(PTVpresc.non_zero_values));    %calculate equivalent uniform dose
ind3 = find(((PTVnormDose3-1.).^2).^(1/2) >= 0.05);                                %find voxels whose plans deviate from prescription > 5 percent
vol3 = 100.-double(length(ind3))./(double(PTVnormDose3.Nind))*100;                  %calculate PTV percentage whose plan is within 5 percent prescription



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*dHistPTVnormDose1./(double(PTVnormDose1.Nind));            %calculate frequency percentage of integral


QVHptv = flipud(cumsum(flipud(dHistPTVnormDose)));                              %cumulatively sum histogram from max min to min bin

plotQVH(dvhbins,QVHptv)                                                         %plot QVH