%Save this file in Matlab as eudmodel.m
%EUDMODEL(DVH), where DVH is a 2 column matrix corresponding to the cumulative, not
%differential, dose volume histogram. The 1st column corresponds to increasing absolute dose or 
%percentage dose values, and the 2nd column to the corresponding absolute or relative volume value.
%The matrix must have a minimum of two rows, and both columns must be of equal length.
%by Hiram A. Gay, MD
%Revised July 8 2007

function probability = eudmodel(dvh)
%user input section
clc; disp('Welcome to the Equivalent Uniform Dose (EUD)-Based Model Program'); disp(' ');
disp('Please note that: 1) the variable dvh should be a CUMULATIVE, not differential, DVH');
disp('                  2) the program assumes that all treatment fractions are equal');
disp(' '); disp(' ');
%end of user input section

%verifying that the cumulative DVH has at least 2 rows and columns
 [nb,N]=size(dvh);
if (nb < 2)
    disp('Error: Cumulative dvh must have at least 2 rows.'); return;
end
if (N < 2)
    disp('Error: Cumulative dvh must have at least 2 columns.'); return;
end
%verifying that the cumulative DVH has no negative numbers in the dose or volume columns
for i=1:nb
    if (dvh(i,1) < 0)
        message = sprintf('Error: Dose data error. dvh column 1, row %g is negative',i);
        disp(message); return;
    end
    if (dvh(i,2) < 0)
        message = sprintf('Error: Volume data error. dvh column 2, row %g is negative',i);
        disp(message); return;
    end
end
% Converting cumulative DVH to differential DVH, and checking for DVH errors
for i=2:nb
    dvh(i-1,1)=dvh(i-1,1)+(dvh(i,1)-dvh(i-1,1))/2;
    if (dvh(i,1)-dvh(i-1,1) <= 0)
        message = sprintf('Error: Dose data error. dvh column 1, row %g <= dvh column 1, row %g',i,i-1);
        disp(message); return;
    end
    dvh(i-1,2)=(dvh(i-1,2)-dvh(i,2));
    if (dvh(i-1,2) < 0)
        message = sprintf('Error: Volume bin < 0. Verify dvh column 2, rows %g and %g',i-1,i);
        disp(message); return;
    end
end
dvh(nb,:)=[];
[nb,N]=size(dvh);
nf=input('Enter the number of treatment fractions: '); disp(' '); disp(' ');
disp('Is the DVH dose data in: ');
disp('     1. percentage dose format');
disp('     2. absolute dose format');
dose_type=input('Enter 1 or 2: '); disp(' '); disp(' ');
%if DVH dose data is in percentage dose format
if (dose_type==1)
    normalized_fraction=input('Enter the dose in Gy (not cGy) corresponding to the 100% dose for ONE fraction: ');
    disp(' '); disp(' ');
   %converting percentage dose bins into absolute dose bins
    for i=1:nb
        dvh(i,1)=dvh(i,1)*nf*normalized_fraction/100;
    end
    message = sprintf('The maximum dose was %g Gy. Is this number reasonable?',dvh(nb,1));
    disp(message);
    disp('     1. yes');
    disp('     2. no');
    answer=input('Enter 1 or 2: '); disp(' '); disp(' ');
%if DVH dose data is in absolute dose format
    if (answer == 2 )
        disp('Error: if the maximum dose was too high: ');
        disp('         1) the dose data could be in ABSOLUTE, not percentage, dose format.');
        disp('         2) the 100% dose entered was for more than 1 fraction.'); return;
    elseif (answer == 1)
    else
        disp('Error: Invalid choice. Exiting program.'); return;
    end
elseif (dose_type==2)
    disp('Is the DVH absolute dose data in: ');
    disp('     1. Gy');
    disp('     2. cGy');
    answer2 = input('Enter 1 or 2: ');
%if DVH dose data is in cGy it is converted to Gy
    if (answer2 == 2)
        for i=1:nb
            dvh(i,1)=dvh(i,1)/100;
        end
    elseif (answer2 == 1)
    else
        disp('Error: Invalid choice. Exiting program.'); return;
    end
else
    disp('Error: Invalid choice. Exiting program.'); return;
end
%EUD mathematical model parameters input section
clc; disp('Does the DVH correspond to:');
disp('     1. tumor target');
disp('     2. normal tissue')
tissue_type=input('Enter 1 or 2: '); disp(' ');
if (tissue_type==1)
    clc
    disp('Structure (Source)                  End-point            a*'); disp(' ');
    disp('Breast (Brenner28)                    Local control        -7.2');
    disp('Melanoma (Brenner28)                  Local control       -10');
    disp('Squamous cc (Brenner28)               Local control       -13');
    disp('* = Niemierko'); disp(' ');
    a=input('Enter the value of parameter a: ' );
    gamma50=input('Enter the value of parameter gamma50 (recommend 2 if unknown): ' );
    tcd50=input('Enter the TCD50 (Gy): ');
    standard_fractionation=input('Enter the source data''s dose per fraction (Gy): ');
    ab=input('Enter the tumor alpha/beta ratio (Gy): ');
elseif (tissue_type==2)
    clc
    disp('Normal tissue EUD Parameters:'); disp(' ');
    disp('Structure (Source)    End-point                    a* / a**  g50**  TD50***  DPF****');
    disp(' ');
    disp('BRAIN (Emami)         Necrosis                        /  5    3      60      1.8 - 2');
    disp('Brainstem (Emami)     Necrosis                        /  7    3      65      1.8 - 2');
    disp('Optic chiasm (Emami)       Blindness                  / 25    3      65      1.8 - 2');
    disp('Colon (Emami)         Obstruction/perforation         /  6    4      55      1.8 - 2');
    disp('Ear(mid/ext) (Emami)  Acute serous otitis             / 31    3      40      1.8 - 2');
    disp('Ear(mid/ext) (Emami)  Chronic serous otitis           / 31    4      65      1.8 - 2');
    disp('Esophagus (Emami)     Perforation                     / 19    4      68      1.8 - 2');
    disp('Heart (Emami)         Pericarditis                    /  3    3      50      1.8 - 2');
    disp('Kidney (Emami)        Nephritis                       /  1    3      28      1.8 - 2');
    disp('Lens (Emami)          Cataract                        /  3    1      18      1.8 - 2');
    disp('Liver (Dawson29)        Liver failure               0.9 /');
    disp('Liver (Emami)         Liver failure                   /  3    3      40      1.8 - 2');
    disp('Liver (Lawrence30)      Liver failure               0.6 /');
    disp('Lung (Emami)          Pneumonitis                     /  1    2    24.5      1.8 - 2');
    disp('Lung (Kwa31)            Pneumonitis                 1.0 /');
    disp('Optic nerve (Emami)   Blindness                       / 25    3      65      1.8 - 2');
    disp('Parotids (Chao32)       Salivary function (<25%)    0.5 /');
    disp('Parotids (Eisbruch24)   Salivary function (<25%)   <0.5 /');
    disp('Retina (Emami)        Blindness                       / 15    2      65      1.8 - 2');
    disp('Spinal cord (Powers33)  White matter necrosis        13 /');
    disp('* = Niemierko / ** = Gay  / ***= Emami / **** = dose per fraction'); disp(' ');
    a=input('Enter the value of parameter a: ' );
    gamma50=input('Enter the value of parameter gamma50 (recommend 4 if unknown): ' );
    td50=input('Enter the TD50 (Gy): ');
    standard_fractionation=input('Enter the source data''s dose per fraction (Gy): ');
    ab=input('Enter the normal tissue alpha/beta ratio (Gy): ');
else
    disp('Error: Invalid choice. Exiting program.'); return;
end
%end of EUD mathematical model parameters input section
total_volume=0;
%calculating the biologically equivalent dose and the total volume
for i=1:nb
    bndvh(i,1)=dvh(i,1)*((ab+dvh(i,1)/nf))/(ab+standard_fractionation);
    total_volume=dvh(i,2)+total_volume;
end
%normalizing volume data to 1 (therefore, total volume corresponds to 1)
for i=1:nb
    dvh(i,2)=dvh(i,2)/total_volume;
    bndvh(i,2)=dvh(i,2);
end
eud=0;
%calculating the EUD
for i=1:nb
    eud=eud+(bndvh(i,2))*(bndvh(i,1))^a;
end
eud=eud^(1/a);
disp(' '); disp(' ');
message = sprintf('The equivalent uniform dose = %g Gy',eud);
disp(message); disp(' ');
%Results section
if (tissue_type==1)
    %calculating tumor contol probability
    tcp=1/(1+((tcd50/eud)^(4*gamma50)));
    tcp=tcp*100;
    message = sprintf('The tumor control probability = %10.10f %%',tcp);
    disp(message); disp(' ');
elseif (tissue_type==2)
    %calculating normal tissue complication probability
    ntcp=1/(1+((td50/eud)^(4*gamma50)));
    ntcp=ntcp*100;
    message = sprintf('The normal tissue complication probability = %10.10f %%',ntcp);
    disp(message); disp(' ');
end
%end of results section