close all; clear all;

k = 50;     % image size is k x k

m = k*k;    % total number of pixels

% define regions
tumor   = false(k);
crit    = false(k);

tumor(21:30, 20:32) = true;    % tumor
crit(33:40, 16:25) = true;     % critical structure 1
crit(23:32, 38:45) = true;     % critical structure 2
other = ~(tumor | crit);       % other region

tumor_min = 0.9; tumor_max = 1.5;
other_min = 0.0; other_max = 0.5;  
crit_min  = 0.0; crit_max  = 0.1;  

Dmin = zeros(k);
Dmin(other) = other_min;    % regular cells
Dmin(tumor) = tumor_min;    % tumor
Dmin(crit)  = crit_min;     % critical regions

Dmax = zeros(k);
Dmax(other) = other_max;    % regular cells
Dmax(tumor) = tumor_max;    % tumor
Dmax(crit)  = crit_max;     % critical regions


Ntheta = 20;    % number of angles (equally spaced from 0 to pi)
Nd = 20;        % number of parallel lines for each angle
n = Nd*Ntheta;  % total number of measurements

% build the matrix A
i = 1;
for itheta = 1:Ntheta,
    for id = 1:Nd,
        lines_d = 0.2*k*(id - Nd/2-0.5)/(Nd/2);
        lines_theta = pi*itheta/Ntheta;
        
        L = line_pixel_length(lines_d, lines_theta, k);
        
        A(:,i) = L(:);
        
        i = i+1;
    end
end

% dmin and dmax contain the vectorized version of the images Dmin and Dmax
dmin = Dmin(:);
dmax = Dmax(:);

% ----------------------------------------------------------
% Plot area geometry
% ----------------------------------------------------------
% Sets the colormap so that:
%   tumor regions are in red
%   critical regions are in blue
%   other regions are in green
%
figure
imagesc(Dmax)
C = [0 0 0.5;0 0.5 0; 0 0.5 0; 0 0.5 0; 0.5 0 0];
colormap(C)
axis image; axis off
title('Area geometry')
% print -depsc -zbuffer Dmax.eps

% ----------------------------------------------------------
% Plot dosage level for uniform beam levels
% ----------------------------------------------------------
b = 0.05*ones(n,1);     % uniform level of 0.05 for all beams
d = A*b;

figure
imagesc(reshape(d,k,k))
axis image; axis off
title('Dosage (uniform beam pattern)')
% print -depsc -zbuffer unif_D.eps

% ----------------------------------------------------------
% Plot histogram of pixel dosages for the three regions
% ----------------------------------------------------------
% You may simply copy and paste the following code to plot 
% the histogram for your solution (assuming you store it in
% the variable d).
%
figure
subplot(3,1,1)
hist(d(tumor),0:0.05:1.5); hold on; v = axis;
plot([tumor_min tumor_min]-0.025, [0 v(4)], 'k--', ...
     [tumor_max tumor_max]+0.025, [0 v(4)], 'k--', 'LineWidth', 2)
hold off;
xlim([-0.1 1.6]); ylabel('tumor');
h = findobj(gca,'Type','patch');
set(h, 'FaceColor',[0.5 0 0])

title('Dosage histograms (uniform beam pattern)') % set title of figure

subplot(3,1,2)
hist(d(crit),0:0.05:1.5);  hold on; v = axis;
plot([crit_min crit_min]-0.025, [0 v(4)], 'k--', ...
     [crit_max crit_max]+0.025, [0 v(4)], 'k--', 'LineWidth', 2)
hold off;
xlim([-0.1 1.6]); ylabel('critical')
h = findobj(gca,'Type','patch');
set(h, 'FaceColor',[0 0 0.5])

subplot(3,1,3)
hist(d(other),0:0.05:1.5); hold on; v = axis;
plot([other_min other_min]-0.025, [0 v(4)], 'k--', ...
     [other_max other_max]+0.025, [0 v(4)], 'k--', 'LineWidth', 2)
hold off;
xlim([-0.1 1.6]); ylabel('other')
h = findobj(gca,'Type','patch');
set(h, 'FaceColor',[0 0.5 0])
% print -depsc unif_hist.eps
