%% Image Reconstruction Using the MATLAB Function Block
% This example shows how to use HDL Coder(TM) to check, generate
% and verify HDL for a fixed-point image reconstruction model using the
% Filtered Back-Projection algorithm.
%
% To run this example, the following additional products are required:
% Image Processing Toolbox(TM), DSP System Toolbox(TM).

% Copyright 2007-2012 The MathWorks, Inc.


%% Introduction
% The algorithm used in this model is based on the function |iradon| from
% the Image Processing Toolbox. The model takes serialized parallel-beam
% projection data, performs reconstruction with the specified image size
% and angle resolution, and outputs a partial pixel data stream. The
% partial pixels are stored to the MATLAB(R) workspace, and then summed up
% in fixed-point arithmetic to form the final image.
%
% For more information on the Inverse Radon transform and the Filtered
% Back-Projection algorithm, see the Image Processing Toolbox documentation
% on |iradon| and the example |Reconstructing an Image from Projection Data|.
%
% *Additional Requirements:*
%
% * Signal Processing Toolbox
% * Image Processing Toolbox
% * DSP System Toolbox

%% Open the Example Model
modelname = 'hdlcoderrecon2'; 
open_system(modelname);


%% Model Setup
% The test image is the Shepp-Logan head phantom generated from the
% function |phantom|. The image is converted to synthetic projections as
% the input data to the model, using the function |radon|. Both functions
% are from the Image Processing Toolbox.

% Define image size (max = 175)
img_siz = 120;
P = phantom(img_siz);
imshow(P, 'InitialMagnification', 180)
title('Head phantom')

%%

% Define angle resolution (in degree)
ang_res = 4;
theta = 0:ang_res:180-ang_res;
[R,xp] = radon(P,theta);     
figure('Color', 'w'), imagesc(theta,xp,R)
colormap(hot), colorbar
title('Projections of the head phantom')
xlabel('Parallel Rotation Angle - \theta (degrees)')
ylabel('Parallel Sensor Position - x\prime (pixels)')

%%
% The filtering is performed in the spatial domain. The highpass ramp
% filter is first created in the frequency domain, and then transformed 
% into a 121-tap symmetric FIR filter. This approach saves the need to
% implement FFT and IFFT in hardware.

% FIR Filter definition
order = 256;
w = linspace(0, 2*pi, order*2);
H1 = linspace(0, 1, order+1);
H1 = [H1 H1(end-1:-1:2)];

fil_siz = 60;           % FIR 1/2 order
h1 = ifftshift(ifft(H1));
h1 = h1(length(h1)/2-fil_siz+1:length(h1)/2+fil_siz+1);
h1 = h1/max(h1);

subplot(211), plot(w,fftshift(H1)), axis tight
set(gca, 'XTick', 0:pi/2:2*pi)
set(gca, 'XTickLabel', {'-pi', '-pi/2', '0', 'pi/2', 'pi'})
title('Ramp filter in frequency domain')
subplot(212), plot(h1), axis tight
title('Ramp filter in spatial domain')

%%
% The following code sets up model parameters and fixed-point precisions.

% Pixel parameters - center pixel, pixel range, bit width
x_center = floor((img_siz+1)/2);
x_range = [1-x_center img_siz-x_center];
y_top = x_center - 1;
y_range = [1+y_top-img_siz y_top];
pix_bits = nextpow2(img_siz+1);

% Projection parameters - projection center, length, bit width
prj_center = ceil(size(R,1)/2);
prj_length = size(R,1)-1;
prj_bits = nextpow2(size(R,1));

% Simulation parameters
latency = 2;
pipeline = 4;
sim_time = img_siz^2*(length(theta)+latency) + pipeline - 1;
simout_len = img_siz^2*length(theta);

% Model fixed-point precisions
fpw_cos = 12;                                        % cos, sin width
fpw_prj = 18;                                        % projection width 
fps_prj = fpw_prj - 1 - nextpow2(max(max(abs(R))));  % projection scaling
fpw_fil = 16;                                        % filter width
fps_fil = fpw_fil - 1 - nextpow2(max(abs(h1)));      % filter scaling
fpw_out = 18;                                        % output width
fps_out = fps_prj - (fpw_prj - fpw_out);             % output scaling


%% Modeling FPGA Memory Blocks
% In the _filtering_ subsystem, the output of the filter is stored in a
% ping-pong memory buffer. While the filtered projection of the current
% angle is written into one bank of the buffer, the filtered projection of
% the previous angle is read out from the other bank and used in the
% following Back-Projection step.
%
% This kind of data buffering scheme can be effectively implemented using
% the dedicated RAM resources in an FPGA. In this design, the Dual Port RAM
% block behaviorally simulates a RAM with a write port and a read port, and
% generates HDL code that infers RAM in an FPGA. In our test, the generated
% code was successfully mapped onto Block RAMs in a Xilinx(R) Virtex-4
% device during synthesis using Synplify(R) Pro(R). The generated code was
% also mapped onto M4K blocks in an Altera(R) Stratix II device.

open_system([modelname '/image reconstruction/filtering']);


%% Simulating the Design
% On each clock cycle, the model outputs a new 'partial' pixel containing
% the contribution from one angle, after the initial frame and pipeline
% latency. To form the final image, the contribution from one angle must be
% added to that of the next angle, until the last angle is calculated.
% Therefore it takes a total of (180/angle resolution)*(image size ^2)
% cycles to reconstruct one frame of the image.

currentdir = pwd;
workingdir = tempname;
mkdir(tempdir,strrep(workingdir,tempdir,''));
cd(workingdir);

%%
% Depending on the image size and angle resolution setting and the speed of
% your computer, the simulation may take a few minutes.

disp('Simulating the model...');
sim(modelname);
disp('Model simulation complete.');

%% 
% To simplify the design, the summation of the partial pixels is not
% implemented in the model. However, the effect of fixed-point arithmetic
% is modeled using the Fixed-Point Designer(TM). In hardware, the summation
% requires a large amount of memory storage, and is best implemented using
% RAM blocks in FPGA or external memory such as SRAM or DRAM.

% Set fixed-point math properties
F = fimath('RoundMode', 'floor', 'OverflowMode', 'wrap', ...
    'SumMode', 'SpecifyPrecision', 'SumWordLength', fpw_out, ...
    'SumFractionLength', fps_out-nextpow2(length(theta)),...
    'CastBeforeSum', 0);

% Sum up the partial pixels from all theta
recon_img = simout;
recon_img.fimath = F;
recon_img = reshape(recon_img,img_siz^2,length(theta));
recon_img = sum(recon_img,2);

% Transform back into original dimension
recon_img = reshape(recon_img,img_siz,img_siz);
recon_img = rot90(recon_img);

%%
% Compared to the original head phantom, the reconstructed image has
% considerable noise, due to the relatively large angle resolution setting
% (ang_res = 4). A smaller ang_res will result in a much clearer image, at
% the cost of longer simulation time.

% Plot reconstructed image
m1 = min(min(double(recon_img)));
m2 = max(max(double(recon_img)));
subplot(121), imshow(double(recon_img),[m1 m2],'InitialMagnification',180)
title(['Reconstructed head phantom at ' num2str(ang_res) '\circ resolution'])
subplot(122), imshow(P,'InitialMagnification',180)
title('Original head phantom')


%% Check, Generate and Verify HDL

makehdl([modelname '/image reconstruction'],...
        'TargetDirectory', workingdir);
    
%% 
% The HDL test bench is not generated in this example, due to the large size
% of the test bench data. The following display shows results displayed by
% the ModelSim(R) HDL simulator when a test bench was generated. The display
% shows simulation waveforms near the end of the test bench.
% 
% <<recon_msim_screen.jpg>>


%%
% This concludes the image reconstruction example.

cd(currentdir);


displayEndOfDemoMessage(mfilename)
