Simulink HDL Coder 1.6

Image Reconstruction Using the Embedded MATLAB™ Function Block

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.

To run this demo, the following additional products are required: Image Processing Toolbox™, Signal Processing Blockset™.

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 = '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® Virtex-4 device during synthesis using Synplify® Pro®. The generated code was 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);

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...
Embedded MATLAB parsing for model "hdlcoderrecon2"...Done
Embedded MATLAB code generation for model "hdlcoderrecon2"....Done
Embedded MATLAB compilation for model "hdlcoderrecon2"...        1 file(s) c
opied.
Done
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 'hdlcoderrecon2/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 hdlcoderrecon2/image reconstruction/back projection as C:\TEM
P\R2009bd_340_3120\tp7bcdcde0_0a8c_41bf_af59_40353fc686ab\back_projection.vh
d
### Working on hdlcoderrecon2/image reconstruction/counters as C:\TEMP\R2009
bd_340_3120\tp7bcdcde0_0a8c_41bf_af59_40353fc686ab\counters.vhd
### Working on hdlcoderrecon2/image reconstruction/filtering/Dual Port RAM a
s C:\TEMP\R2009bd_340_3120\tp7bcdcde0_0a8c_41bf_af59_40353fc686ab\Dual_Port_
RAM.vhd
### Working on hdlcoderrecon2/image reconstruction/filtering/Dual Port RAM/D
ualPortRAM_Inst0 as C:\TEMP\R2009bd_340_3120\tp7bcdcde0_0a8c_41bf_af59_40353
fc686ab\DualPortRAM_Inst0.vhd
### Working on hdlcoderrecon2/image reconstruction/filtering/Dual Port RAM1 
as C:\TEMP\R2009bd_340_3120\tp7bcdcde0_0a8c_41bf_af59_40353fc686ab\Dual_Port
_RAM1.vhd
### Working on hdlcoderrecon2/image reconstruction/filtering/Dual Port RAM1/
DualPortRAM_Inst1 as C:\TEMP\R2009bd_340_3120\tp7bcdcde0_0a8c_41bf_af59_4035
3fc686ab\DualPortRAM_Inst1.vhd
### Working on hdlcoderrecon2/image reconstruction/filtering as C:\TEMP\R200
9bd_340_3120\tp7bcdcde0_0a8c_41bf_af59_40353fc686ab\filtering.vhd
### Working on hdlcoderrecon2/image reconstruction/filtering/FIR as C:\TEMP\
R2009bd_340_3120\tp7bcdcde0_0a8c_41bf_af59_40353fc686ab\FIR.vhd
### Starting VHDL code generation process for filter: FIR
### Starting VHDL code generation process for filter: FIR
### Generating: C:\TEMP\R2009bd_340_3120\tp7bcdcde0_0a8c_41bf_af59_40353fc68
6ab\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 hdlcoderrecon2/image reconstruction/filtering/filtering cntrl
 as C:\TEMP\R2009bd_340_3120\tp7bcdcde0_0a8c_41bf_af59_40353fc686ab\filterin
g_cntrl.vhd
Embedded MATLAB parsing for model "hdlcoderrecon2"...Done
Embedded MATLAB code generation for model "hdlcoderrecon2"....Done
### Working on hdlcoderrecon2/image reconstruction as C:\TEMP\R2009bd_340_31
20\tp7bcdcde0_0a8c_41bf_af59_40353fc686ab\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.

This concludes the image reconstruction demo.

cd(currentdir);