Simulink HDL Coder 1.6

Image Reconstruction Using Cosimulation

This demo shows how to use Simulink® HDL Coder™ to check, generate and verify HDL for a fixed-point image reconstruction model using the Filtered Back-Projection algorithm.

Contents

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® 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 demo Reconstructing an Image from Projection Data.

Additional Requirements:

Open the Demo Model

modelname = 'hdlcoderrecon';
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, dpram_model.vhd models a dual-port RAM with a read and a write port. In our test, the RAM models were successfully mapped onto Block RAMs in a Xilinx® Virtex-4 device during synthesis using Synplify® Pro®. The RAM models were also mapped onto M4K blocks in an Altera® 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);

The function hdlcoderreconcmds creates and returns the necessary Tcl commands to cosimulate the VHDL RAM model in ModelSim®.

Before simulation starts, ModelSim must be completely initialized. If ModelSim initialization does not complete within the wait time, the function pingHdlSim returns -1, and the cosimulation will not run properly.

wait_time = 5;   % in minutes
vsim('tclstart', hdlcoderreconcmds);
disp('Waiting for the HDL simulator to initialize...');
if pingHdlSim(wait_time*60) == -1
    error('hdlcoder:hdlcoderdemos:hdlcoderrecon_m:nohdlsimulator',...
          ['HDL simulator not ready for ' num2str(wait_time) ...
           ' minutes -- aborting']);
else
    disp('HDL simulator initialized.');
end

Waiting for the HDL simulator to initialize...
HDL simulator initialized.

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.');

Simulating the model...
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 Toolbox™. 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

checkhdl([modelname '/image reconstruction'],...
         'TargetDirectory', workingdir);
makehdl([modelname '/image reconstruction'],...
        'TargetDirectory', workingdir);

### Starting HDL Check.
### HDL Check Complete with 0 errors, 2 warnings and 0 messages.

### Generating HDL for 'hdlcoderrecon/image reconstruction'

### Applying HDL Code Generation Control Statements
### Starting HDL Check.
### HDL Check Complete with 0 errors, 2 warnings and 0 messages.

Warning: The full input range (2^input_wordlen) of Lookup Tables blocks
should be covered for HDL code generation.
Warning: The full input range (2^input_wordlen) of Lookup Tables blocks
should be covered for HDL code generation.


### Begin VHDL Code Generation
### Working on hdlcoderrecon/image reconstruction/back projection as C:\TEMP
\R2009bd_340_3120\tp878bcf78_4fd2_4e09_9830_16aecb322df5\back_projection.vhd
### Working on hdlcoderrecon/image reconstruction/counters as C:\TEMP\R2009b
d_340_3120\tp878bcf78_4fd2_4e09_9830_16aecb322df5\counters.vhd
### Working on hdlcoderrecon/image reconstruction/filtering as C:\TEMP\R2009
bd_340_3120\tp878bcf78_4fd2_4e09_9830_16aecb322df5\filtering.vhd
### Working on hdlcoderrecon/image reconstruction/filtering/FIR as C:\TEMP\R
2009bd_340_3120\tp878bcf78_4fd2_4e09_9830_16aecb322df5\FIR.vhd
### Starting VHDL code generation process for filter: FIR
### Starting VHDL code generation process for filter: FIR
### Generating: C:\TEMP\R2009bd_340_3120\tp878bcf78_4fd2_4e09_9830_16aecb322
df5\FIR.vhd
### Starting generation of FIR VHDL entity
### Starting generation of FIR VHDL architecture
### HDL latency is 0 samples
### Successful completion of VHDL code generation process for filter: FIR

### Working on hdlcoderrecon/image reconstruction as C:\TEMP\R2009bd_340_312
0\tp878bcf78_4fd2_4e09_9830_16aecb322df5\image_reconstruction.vhd
### HDL Code Generation Complete.

The HDL test bench is not generated in this demo, due to the large size of the test bench data. The following display shows results displayed by the ModelSim HDL simulator when a test bench was generated. The display shows simulation waveforms near the end of the test bench.

If you would like to simulate or synthesize code generated from this model, you will need to copy the RAM model file dpram_model.vhd to your target directory. Then, modify the generated test bench compilation script and synthesis script to include dpram_model.vhd. Finally, edit the generated file filtering.vhd and remove the duplicated dpram_model1 instance.

This concludes the image reconstruction demo. ModelSim will be closed.

tclHdlSim('after 1000 quit -f;');
cd(currentdir);