%
% Figure - A. Holder and J. Reese, 6/4/03
%
% ver. 1.0
%
% This function outputs some useful pictures based upon the
% solution data. 
%

function Figure(Prob, prescription, probname, dirname, SaveType)

%
% This lines up the Radiation Levels into the corresponding
% matrix for displaying the image.
%

n=length(prescription);
DisplayMatrix = zeros(n,n);

for j=1:n
   for i=1:n
        DisplayMatrix(i,j)=Prob.RadLevel(n*(j-1)+i);
   end
end

%
% Pop up some nifty pictures. ViewPrescription is used so that the
% regular restricted regions do not show up on the contour.
%

%
% makesure the directory is there
%
if ~isdir(['../images/',dirname])
   eval(['!mkdir ../images/',dirname])
end

figure;
ViewPrescription = prescription(:,:,1);
ViewPrescription(Prob.RegularRowIndex) = 0;
contour(ViewPrescription)
saveas(gcf, ['../images/',dirname,'/',probname,'-PatientContour'], SaveType);

figure;
ViewPrescription = prescription(:,:,1);
ViewPrescription(Prob.RegularRowIndex) = 0;
contour(ViewPrescription)
hold on;
contour(DisplayMatrix)
hold off;
saveas(gcf, ['../images/',dirname,'/',probname,'-TreatmentContour'], SaveType);

figure;
mesh(DisplayMatrix)
saveas(gcf, ['../images/',dirname,'/',probname,'-TreatmentMesh'], SaveType);

% Lets See if we can some DVHs

TumorRadLevel = (1 / min(Prob.TLB)) * (Prob.TumorDoseMatrix * Prob.x);
CritRadLevel = (1 / min(Prob.TLB)) * (Prob.CriticalDoseMatrix * Prob.x);
k = 1;
TumorVol = zeros(length(0:.01:1.5),1);
CritVol = TumorVol;
for i = 0:.01:1.5
   TumorVol(k) = length(find(TumorRadLevel >= i)) / length(TumorRadLevel);
   CritVol(k) = length(find(CritRadLevel >= i)) / length(CritRadLevel);
   k = k + 1;
end
figure
plot(0:.01:1.5, TumorVol, 'b*-')
hold on
plot(0:.01:1.5, CritVol, 'rs-')
axis([0, 1.5, 0, 1.1])
saveas(gcf, ['../images/',dirname,'/',probname,'-DVH'], SaveType);


