%
% 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);
[m,n,NumSlices,q] = size(prescription);

% Generate Images for middle slice
Slice = round(NumSlices/2);

for i=1:n
   for j=1:m
        DisplayMatrix(i,j)=Prob.RadLevel((Slice-1)*m*n + m*(i-1) + j);
   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;
contour(prescription(:,:,Slice,2)');
saveas(gcf, ['../images/',dirname,'/',probname,'-PatientContour'], SaveType);

figure;
contour(prescription(:,:,Slice,2)');
hold on;
contour(DisplayMatrix)
hold off;
saveas(gcf, ['../images/',dirname,'/',probname,'-TreatmentContour'], SaveType);

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


% JOSH'S NEW STUFF 9/14/05

% Lets See if we can some DVHs

% Find the Critical Structures
CritIndices = unique(prescription(:,:,:,2));
CritIndices = CritIndices(find(CritIndices >= 3));
ind = 1;
for k = 1:length(CritIndices)
   RwInd = find(prescription(:,:,:,2) == CritIndices(k));
   LengthIndex(k) = length(RwInd);
   CritRadLevel(ind : ind + LengthIndex(k) - 1) = ...
      (1 / min(Prob.TLB))*(Prob.DoseMatrix(RwInd,:)*Prob.ThePlan);
   ind = ind + LengthIndex(k);
end

TumorRadLevel = (1 / min(Prob.TLB))*(Prob.TumorDoseMatrix*Prob.x);
RegularRadLevel = (1 / min(Prob.TLB))*(Prob.RegularDoseMatrix*Prob.x);

k = 1;
TumorVol = zeros(length(0:.01:1.5),1);
RegVol = TumorVol;
for i = 0:.01:1.5
   TumorVol(k) = length(find(TumorRadLevel >= i)) / length(TumorRadLevel);
   RegVol(k) = length(find(RegularRadLevel >= i)) / length(RegularRadLevel);
   ind = 1;
   for j = 1:length(CritIndices)
      CritVol(k,j) = ...
         length(find(CritRadLevel(ind : ind + LengthIndex(j) - 1) >= i)) ...
            / length(CritRadLevel(ind : ind + LengthIndex(j) - 1));
      ind = ind + LengthIndex(j);
   end
   k = k + 1;
end
figure
plot(0:.01:1.5, TumorVol, 'b*-')
hold on
plot(0:.01:1.5, RegVol, 'gd-')
for k = 1:length(CritIndices)
   plot(0:.01:1.5, CritVol(:,k), 'rs-')
end
axis([0, 1.5, 0, 1.1])
saveas(gcf, ['../images/',dirname,'/',probname,'-DVH'], SaveType);


