Image Reconstruction Using the MATLAB Function Block

This example shows how to use HDL Coder™ 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™, DSP System Toolbox™.


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 example Reconstructing an Image from Projection DataReconstructing an Image from Projection Data.

Additional Requirements:

  • Signal Processing Toolbox

  • Image Processing Toolbox

  • DSP System Toolbox

Open the Example Model

modelname = 'hdlcoderrecon2';

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;

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...');
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 Designer™. 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);
### Generating HDL for 'hdlcoderrecon2/image reconstruction'.
### Starting HDL check.
### Starting VHDL code generation process for filter: filter
### Begin VHDL Code Generation for 'hdlcoderrecon2'.
### Working on hdlcoderrecon2/image reconstruction/back projection as /tmp/BR2015ad_189898_78091/tp34aa29b1_824c_4d53_845b_257f5ae7a34a/hdlcoderrecon2/back_projection.vhd.
### Working on hdlcoderrecon2/image reconstruction/counters as /tmp/BR2015ad_189898_78091/tp34aa29b1_824c_4d53_845b_257f5ae7a34a/hdlcoderrecon2/counters.vhd.
### Working on hdlcoderrecon2/image reconstruction/filtering/filtering cntrl as /tmp/BR2015ad_189898_78091/tp34aa29b1_824c_4d53_845b_257f5ae7a34a/hdlcoderrecon2/filtering_cntrl.vhd.
### Working on hdlcoderrecon2/image reconstruction/filtering/DualPortRAM_Wrapper_256x18b/DualPortRAM_256x18b as /tmp/BR2015ad_189898_78091/tp34aa29b1_824c_4d53_845b_257f5ae7a34a/hdlcoderrecon2/DualPortRAM_256x18b.vhd.
### Working on hdlcoderrecon2/image reconstruction/filtering/DualPortRAM_Wrapper_256x18b as /tmp/BR2015ad_189898_78091/tp34aa29b1_824c_4d53_845b_257f5ae7a34a/hdlcoderrecon2/DualPortRAM_Wrapper_256x18b.vhd.
### Working on hdlcoderrecon2/image reconstruction/filtering/FIR as /tmp/BR2015ad_189898_78091/tp34aa29b1_824c_4d53_845b_257f5ae7a34a/hdlcoderrecon2/FIR.vhd.
### Working on hdlcoderrecon2/image reconstruction/filtering as /tmp/BR2015ad_189898_78091/tp34aa29b1_824c_4d53_845b_257f5ae7a34a/hdlcoderrecon2/filtering.vhd.
### Working on hdlcoderrecon2/image reconstruction as /tmp/BR2015ad_189898_78091/tp34aa29b1_824c_4d53_845b_257f5ae7a34a/hdlcoderrecon2/image_reconstruction.vhd.
### Generating package file /tmp/BR2015ad_189898_78091/tp34aa29b1_824c_4d53_845b_257f5ae7a34a/hdlcoderrecon2/image_reconstruction_pkg.vhd.
### Creating HDL Code Generation Check Report file:////tmp/BR2015ad_189898_78091/tp34aa29b1_824c_4d53_845b_257f5ae7a34a/hdlcoderrecon2/image_reconstruction_report.html
### HDL check for 'hdlcoderrecon2' complete with 0 errors, 0 warnings, and 0 messages.
### HDL code generation complete.

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® 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 example.

Was this topic helpful?