Abstract: In this study, we present an implementation on FPGA of 2D signals Encoder/Decoder using dyadic Discrete Wavelet Transform based on quadrature mirror filters, by applying fast wavelet Mallat`s algorithm. The wavelet coefficients will be encoded by Huffman code in order to be transmitted progressively through an Ethernet TCP/IP based connection. The proposed study is implemented and synthesized in VHDL for Xilinx Virtex-IIV2MB1000 FPGA device using ISE 8.1 and simulated on Modelsim PE 6.0d. The synthesis results are presented in detail. The proposed design can substantially accelerate the DWT and the possible reconfiguration can be exploited to reach a higher performance in the future. The system is designed to be integrated as an extension to the nuclear imaging system implemented around our nuclear research reactor. Assuming a Pentium4 processor with clock frequency of 3.3 GHz for the Matlab software implementation, a speed up of over 5 times for a picture size of 256 x 256 was achieved.
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, L2(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) = 2j/2 φ(2j 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.
RESULTS
The VHDL implementation (Sanchez, 2003) of DWT/IDWT is simulated on ModelSim PE 6.0d, for correct functionality, before being synthesized using XILINX ISE8.1. Radiological images are converted to memory bit files and vice versa, initial data to be processed is written inside the memory (8 bits for each pixel) and the results (8 bits for integer part of coefficients and 10 bits for fractional part) are stored back in the same memory. Some processed images are shown in Appendix C-F. A comparison based on compression ratio (CR), Root mean square error (RMSE) and compression speed is made between VHDL implementations and Matlab7 (Gonzalez, 2004) programs, the CR and the RMSE are varying with image type and size as in Table 1. Performance results for 1-levet DWT are shown in ISE reports in Table 2.
| Table 1: | Compression ratios and RMS errors for different images |
| Table 2: | ISE synthesis reports |
CONCLUSION
We have implemented a Wavelet transform based image encoder/decoder on re-programmable hardware FPGA. The system has multiple configurations which support different compression levels, depending on the user input choices. Lessons learned from this study will help us enhance similar implementations in the future for further improvements using the same reconfigurable IC, especially when the whole tomography system will be achieved.
APPENDIX
| Fig. A: | Virtex-IIV2MB1000 System Board |
| Fig. B: | Download setup |
| Fig. C: | Cameraman image and its level 3 DWT |
| Fig. D: | Neutron radiography image and its level 3 DWT |
| Fig. E: | Radio isotopic image and its level 3 DWT |
| Fig. F: | Computed tomography 3D image (a), its level 3 DWT (b) and the reconstructed image from DWT coefficients (c) |