FPGA Implementation of 2D Signals Encoder Using QMF Based Dyadic DWT: Application to Neutron Tomography Projections

INTRODUCTION

---

Wavelet Transform has been successfully applied to different fields, ranging from pure mathematics to applied science (Nibouche et al., 2000). Numerous studies, carried out on Wavelet Transform, have proven their advantages in image processing and data compression and have made it a basic encoding technique in recent data compression standards. Purely software implementations of the Discrete Wavelet Transform, however, appear to present a problem when real time systems are required in terms of performance. Therefore, hardware acceleration of the DWT has become a topic of recent research.

The basic motivation of our application to radio isotopic images and neutron tomography projections, (Fig. 1) is to reduce the data volume and to achieve a low bit rate in the digital representation images without apparent loss of image quality, in order to reach good 3D reconstruction on a distant Computer or storage for future retrieval.

Radiological image is a digital 2D image acquired by a certain radiological procedure which can be: X-rays, neutron radiography, gamma camera image or a nuclear magnetic reso nance image. For tomography, the gray level value represents the relative linear attenuation coefficient of the object.

Fig. 1:	Neutron imaging system (Kharfi et al., 2005)

Higher quality images show finer structural or functional information of body organs and support more reliable diagnostic outcomes (Wong et al., 1995).

Present approach for image compression helps a lot for 2D image storage and transmission. It is based on projection images collection, encoding/decoding, TCP/IP transmission and therefore 3D reconstruction. The design is based on the MALLAT’s Fast Wavelet Transform Algorithm, which is a fast implementation of the Discrete Wavelet Transform. The design utilizes various digital techniques and specific features of the Xilinx Virtex-IIV2MB1000 FPGA (Appendix A and B) to accelerate the transform.

Characteristics of digital radiological images: Digital radiological image is a digital image acquired by a certain radiological procedure which can be X-rays, neutron radiography, gamma camera image or a nuclear magnetic resonance image. It is a two-dimensional MxN array of non-negative integers (gray levels) (Fig. 1).

The performance of radiological image compression depends not only on the compression ratio but also the image quality of reconstructed images. Higher quality images show finer structural or functional information of body organs and support more reliable diagnostic outcomes.

Image compression framework: The general framework for radiological image compression based on wavelet transform (Fig. 2) is similar to other digital compression fields, the framework includes three major stages: Image transformation, quantization (irreversible compression only) and entropy encoding. The relative importance of each stage varies from one technique to another. All reversible compression techniques do not involve the stage of quantization (Wong et al., 1995). It is useful to consider entropy encoding as a two-step process. The first step uses a statistical model to convert coefficients into an intermediate sequence of symbols. The second step converts the symbol to a data stream in which the symbols no longer have externally identifiable boundaries. The entropy coding uses a variable-length code to achieve higher compression rates, such as Huffman coding. The same code tables used to compress an image are needed to decompress it. A promising approach of image compression is progressive transmission that transmits image data in stages and at each stage, an approximation to the original image is reconstructed at the receiver.

A radiological image compression algorithm that provides multi-resolution image representation offers many advantages.

Wavelet transform theory: Formally, the wavelet transform is defined by many authors as a mathematical technique in which a particular signal is analyzed (or synthesized) in the time domain by using different versions of a dilated (or contracted) and translated (or shifted) basis function called the wavelet prototype or the mother wavelet.

The CWT as defined by Morlet and Grossman is (Nibouche et al., 2002):

(1)

Fig. 2:	Wavelet-based encoding scheme (Nibouche et al., 2002)

Where f(t) belongs to the square integrable functions space, L²(R). In the same way, the inverse CWT can be defined as:

(2)

Where a and b are, the scaling and the shifting factors and to ease the implementation of a wavelet system, these factors have been adopted to be a factor of two. i.e., (Nibouche et al., 2002):

a = 2-j, b = 2-j.k

The reconstruction is only possible if the factor C_ψ is defined (<∞), thus ψ(t) (mother wavelet) must have zero mean and finite energy (admissibility condition). The digital realization of CWT is DWT, given as the inner product of the signal x(t) being transformed with each of the discrete orthogonal/orthonormal set of basis wavelets.

Based on Mallat’s work on the relationships between the Quadrature Mirror Filters (QMF), pyramid algorithms and orthonormal wavelet bases, Daubechies succeeded to construct a set of wavelet orthonormal basis functions which have become the cornerstone of many applications such as JPEG2000 compression standard (Mallat, 1989).

Multiresolution, filter banks and DWT implementation: This concept has been introduced first by Mallat (1989) and then Nibouche (2002). It defines clearly the relationships between the QMF, pyramid algorithms and orthonormal wavelet bases. The strength of multiresolution lies in its ability to decompose a signal in finer and finer details. Most importantly, it allows the description of a signal in terms of time-frequency or time-scale analysis.

Scaling function: Let the scaling function be defined as:

φk(t) = φ(t-k) K∈Z φ∈L2(R)

By scaling and translating, a 2-dimensional scaling function is generated:

To ease the implementation:

Thus: φ_j,k(t) = 2^j/2 φ(2^j t-k), a = 2^-j, b = 2^-j.k

By spanning the whole space:

(3)

Where h(k) are the scaling function coefficients. Equation 3 is called the multiresolution analysis equation

By the same way the wavelet function:

ψk(t)

=

ψ(t-k) K∈Z ψ∈L2(R)

Two-dimensional scaled and translated wavelet function is:

ψj,k(t) = 2j/2 ψ(2jt – n); j,k∈Ζ such that - ∞<j,k<∞

where Z is set of integers.

In terms of the scaling function:

(4)

Where g(k) are the wavelet coefficients. Equation 4 leads to a dyadic decomposition.

For practical and computational reasons, discrete time filter banks are required. Such structures decompose a signal into a coarse representation along with added details (Nibouche et al., 2002). We define the relationship between the expansion coefficients at lower and higher scale levels, by using a scaled and shifted version of Eq. 3 and 4:

(5)

(6)

For finite length signals (which is the case for digital images), the use of a Finite Impulse Response filter (FIR) is the most appropriate choice. Since Eq. 5 and 6 compute one output for each two consecutive inputs, a basic operation required here consists of using a down-sampler or decimator by a factor of two (Fliege, 1994) (Fig. 3).

Fig. 3:	One dimensional discrete wavelet transform

Fig. 4:	Two-dimensional Mallat’s analysis and synthesis algorithm

The filter bank, used to implement Mallat’s Algorithm, is defined as a combination of a low pass filter and high pass filter, both followed by a factor of two decimation. To allow further level of decomposition, identical stages are cascaded leading to a multiresolution analysis. Figure 4 shows analysis/synthesis:

For two-dimensional (Image) analysis and synthesis. The filter bank used is the (2,2) Cohen-Debuchies-Feaveu wavelet filter because of its good compression characteristics (Cohen et al., 1992). It is also the fastest and the simplest one to implement in hardware. Polyphase decomposition is adopted to remedy to the problem of resources lost (Fliege, 1994).

The range of coefficients for each sub-band is divided into 16 quantization levels. The coefficients are quantized into one of the 16 possible levels and those of small magnitude can be neglected without considerable distortion to the image. Huffman trees are an efficient way of coming up with a variable length encoding for a set of characters, given the relative occurrence probabilities.

FPGA implementation of DWT: Configuring FPGA chips is not a straightforward operation. Normally, it involves a number of different steps (Fig. 5) in which different tools can be used:

•	Choosing a design capture tool (VHDL or schematic).
•	Converting the file into a netlist format.
•	Generating the bit stream necessary for configuring the FPGA, the netlist file is passed through implementation tools (e.g., map, place and route tools).

The input image to the encoder is raw gray scale frames of 512 by 512 pixels (at maximum). Each pixel is represented by 256 gray scale levels (8 bits).

Fig. 5:	Design Flow (Nibouche et al., 2002)

Input frames are loaded to the 16 bits word memory by a computer and DWT results are read back, once the FPGAv has processed it. Dynamic vector quantization and entropy encoding on the wavelet coefficients are then made.