Design Technique and Efficient Polyphase Implementation for 2D Elliptically Shaped FIR Filters

TL;DR

This paper introduces a new method for designing 2D FIR filters with elliptical frequency responses, enabling efficient image processing and accurate reconstruction.

Contribution

A novel analytical approach for designing elliptically shaped 2D FIR filters with directional selectivity and orientation control.

Findings

The proposed filters achieve accurate shapes with negligible distortions.

The polyphase implementation reduces computational complexity compared to direct filtering.

The filters efficiently extract straight lines from images with specified orientations.

Abstract

This paper presents a novel analytical approach for the efficient design of a particular class of 2D FIR filters, having a frequency response with an elliptically shaped support in the frequency plane. The filter design is based on a Gaussian shaped prototype filter, which is frequently used in signal and image processing. In order to express the Gaussian prototype frequency response as a trigonometric polynomial, we developed it into a Fourier series up to a specified order, given by the imposed approximation precision. We determined analytically a 1D to 2D frequency transformation, which was applied to the factored frequency response of the prototype, yielding directly the factored frequency response of a directional, elliptically shaped 2D filter, with specified selectivity and an orientation angle. The designed filters have accurate shapes and negligible distortions. We also designed a 2D uniform filter bank of elliptical filters, which was then applied in decomposing a test image into sub-band images, thus proving its usefulness as an analysis filter bank. Then, the original image was accurately reconstructed from its sub-band images. Very selective directional elliptical filters can be used in efficiently extracting straight lines with specified orientations from images, as shown in simulation examples. A computationally efficient implementation at the system level was also discussed, based on a polyphase and block filtering approach. The proposed implementation is illustrated for a smaller size of the filter kernel and input image and is shown to have reduced computational complexity due to its parallel structure, being much more arithmetically efficient compared not only to the direct filtering approach but also with the most recent similar implementations.

Figures

Equations

Full text

1. Introduction

In recent years, the technology, architecture and dedicated image processing algorithms for modern image sensors have known a spectacular development. For example, high-resolution image sensors for remote sensing applications must yield clear, low-noise images, mostly in visible but also in infrared or microwave domains. In order to provide the required high quality, relevant and precise information, images acquired by such sensors must undergo pre-processing stages using various enhancement and restoration operations. In this wide range of image processing tasks, digital filters and filter banks play a fundamental role, for instance, in analyzing images and decomposing them in a set of sub-band images, extracting relevant features and details. They are also extensively used in automotive, where high-speed image processing and feature extraction is highly demanded in real-time applications, such as driver assistance systems, etc.

Two-dimensional filters have been extensively investigated by many researchers in the vast field of digital signal and image processing due to their important applications in image processing, and a wide range of design techniques have been developed [1]. Unlike global numerical optimization algorithms, analytical design methods are based on 1D prototypes with imposed shapes and parameters, to which frequency transformations are applied, leading directly to the desired 2D filters. Perhaps the major benefit of the analytical approach is that the factored frequency response results, and the designed 2D filter is parametric and therefore easily adjustable.

There exists in the literature a wide variety of 2D filters, of the FIR or IIR type, with various shapes, each of them having applications in image processing tasks. One of the most currently used methods used in the analytical design of FIR filters is the McClellan transform [2,3]. Some important works in the field propose various computationally efficient 2D FIR filter structures also using frequency transformations, like multiplier-less circularly symmetric FIR filters using harmony search algorithms [4], sampling-kernel based interpolation [5], a design based on Gaussian approximation [6], and a single-parameter tunable 2-D Farrow structure [7]. Various transformations for the design of FIR and IIR filters are given in [8], and the efficient synthesis of 2D FIR filters with improved selectivity is approached in [9]. Two-dimensional filter banks of various types and shapes are widely used in essential applications, like texture analysis and segmentation, or feature extraction tasks. They decompose the frequency spectrum of the image into several sub-bands, which is required in coding applications, etc.

Different types of two-dimensional filter banks have been investigated [10], in particular directional [11] and nonuniform frequency plane partitioning [12]. An application of directional filter banks in fingerprints analysis is given in [13]. Elliptically shaped 2D filters, including directional ones, were approached in early papers such as [14,15]. A 2D wavelet filter bank with elliptical support was applied in iris recognition [16] and directional elliptical filters were used in robust human pose detection [17]. The design of 2D FIR filters with a circular and elliptical shape based on a swarm optimization algorithm is proposed in [18].

Gaussian separable filtering can also be combined with other methods, such as integral image techniques to construct a Gaussian pyramid following an SIFT algorithm; it may be followed, for instance, by Canny edge detection and dilation operations in order to obtain an optimized edge detection, with superior performances and reduced complexity, also suitable for real-time processing [19]. Such analytical design methods could also be extended to 3D image processing, with interesting applications in physics. For instance, paper [20] presents a 3D digital image correlation approach for accurately measuring strain and deformations in vinyl chloride-coated metal multilayer sheets.

Some computationally efficient structures of 2D FIR filters using block filtering, systolic arrays, etc., are elaborated on in [21,22].

Several analytical design techniques for various 2D filters have also been proposed in previous works. Thus, we studied a class of 2D FIR circular filter banks and their efficient implementation using a polyphase approach in [23], the design and applications of adjustable elliptical and circular 2D IIR filters in [24] and another type of directional elliptically shaped IIR filter, with applications in the detection of oriented straight lines and other details in images, in [25].

In this paper, an analytical design procedure is elaborated on for a particular class of 2D filter, specifically 2D FIR Gaussian elliptical filters. A uniform elliptical 2D filter bank is also designed, with a specified number of component filters, based on a corresponding 1D prototype with a Gaussian frequency response. This is easily approximated using Fourier series expansion. By simply shifting the low-pass filter to equally spaced peak frequencies, we also obtain band-pass filters, thus generating a uniform filter bank. We have determined a particular 1D to 2D frequency transformation, which is applied to the 1D prototype and yields directly the desired 2D filter with elliptically shaped support. The 2D filters result in accurate shapes and negligible distortions. The designed elliptical FB was tested on a greyscale image, which was decomposed into sub-band components, and then the original image was reconstructed very accurately from these components. In a particular case, very selective low-pass directional filters can detect and extract straight lines with specified orientation from images, as shown by the simulation results.

A novel and more efficient implementation solution is also proposed for the designed 2D FIR filter. It has a significantly lower complexity than a recent implementation in [21] while maintaining the same throughput. The proposed solution is based on a new algorithm for 2D FIR filtering with a smaller kernel of 3 × 3, which has been extended at a larger kernel size of 30 × 30 using a polyphase decomposition of the given 2D filtering operation and block filtering. In order to facilitate the reformulation of the algorithm for 2D FIR filtering for a polyphase implementation, we have introduced some auxiliary vectors. Thus, it was possible to extend the proposed 2D filtering algorithm for a significantly larger size.

This paper is organized as follows. Section 2 presents the analytical design technique, first approximating a Gaussian prototype using Fourier series, then generating an uniform Gaussian filter bank with a specified number of components, peak frequencies and bandwidths. Next, a specific 1D to 2D frequency transformation is determined, which will be applied to the given prototype and will lead directly to the desired elliptically shaped 2D filter with a specified orientation in the frequency plane. Next, several design examples are given, showing the 2D filter characteristics (the frequency response and contour plot), for the components of an elliptical filter bank, and then for very selective filters with various orientations. In Section 3, we present an example of image analysis and reconstruction using the elliptical filter bank by decomposing it into sub-band images. A simulation of straight line extraction from a test image, using narrow directional elliptical filters, is also given. The proposed filtering implementation at the system level, based on the polyphase and block structure, is described in Section 4. Some comparative discussions regarding the filter design and computational complexity of the proposed implementation are provided in Section 5. Finally, conclusions are given in the last section.

2. Analytical Design Technique for 2D Elliptical FIR Filters

In this main section, a novel analytical design procedure for a class of 2D FIR filters with elliptically shaped frequency response is presented. The proposed technique is based on a suitable 1D prototype with imposed specifications (peak frequency and bandwidth), to which 1D to 2D frequency mapping is applied, which directly leads to the desired 2D filter. We have chosen as a prototype a Gaussian-shaped filter, which has the useful property of scalability on the frequency axis. A uniform filter bank prototype is first derived, which yields a 2D elliptically shaped filter bank through the elaborated frequency mapping.

2.1. Design of Gaussian FIR Filter Prototype Based on Fourier Series Approximation

As is well known, the Gaussian filter function in the frequency domain can be expressed as $[eqn]$ , where parameter $[eqn]$ represents variance. In order to derive a simpler, more practical form, we will apply the substitution $[eqn]$ , otherwise written as $[eqn]$ . In this way, the low-pass filter with a Gaussian shape will be expressed in the more convenient form $[eqn]$ . We consider that the Gaussian filter is defined by the range $[eqn]$ as its support by extending it on the entire frequency axis as a periodic function with period $[eqn]$ . We derive the following approximation through a trigonometric polynomial, by developing the Gaussian $[eqn]$ in Fourier series, up to a given order, N:

$$[eqn]$$

From this expression of the Gaussian low-pass prototype, we also obtain a band-pass (BP) prototype by shifting this low-pass filter to the symmetric frequencies $[eqn]$ :

$$[eqn]$$

Expressions (1) and (2) above can easily be calculated in a MATLAB routine. Next, the low-pass, band-pass and high-pass components of the specified FIR filter bank prototype will be derived.

2.2. The Design of a Gaussian Uniform FIR Filter Bank Prototype

In this section, a uniform, Gaussian filter bank (FB) prototype with seven-component filters will be designed, consisting in one low-pass (LP) filter, five band-pass (BP) filters and one high-pass (HP) filter.

As in every uniform filter bank, the peak frequencies are equally spaced on the frequency axis. For the 5 BP filters, we impose a bandwidth equal to $[eqn]$ , while the LP and HP filters have each half of this bandwidth, namely $[eqn]$ . The k-th ideal Gaussian band-pass filter is derived by shifting the low-pass prototype to the frequency $[eqn]$ , and its frequency response has the expression

$$[eqn]$$

Next, scaling parameter p has to be calculated for the specified bandwidth value. For our purpose, the filter bandwidth is defined at 0.5 of the peak value (at 6 dB). The magnitude characteristics of any two consecutive, adjacent filters partially overlap and intersect at the value 0.5. Considering the LP filter given by $[eqn]$ , we apply the condition $[eqn]$ , also written $[eqn]$ ; thus, the value of scaling parameter p is $[eqn]$ . Since the bandwidth of each component filter is $[eqn]$ , we obtain the scaling parameter value $[eqn]$ . The ideal uniform Gaussian filter bank is plotted in Figure 1a. For a very accurate approximation and a small ripple level, the Fourier series given by (1) and (2) was truncated at a number of N = 12 terms. The approximation order must be a trade-off between a small distortion level and a reasonable size of the filter matrices for a convenient implementation.

In the next step, the components of the uniform prototype filter bank are determined. Following the steps of the design method described above, knowing the bandwidths and peak frequencies of component filters, using Equations (1) and (2), the frequency responses of the FB components are calculated.

For instance, the LP Gaussian prototype, approximated by a truncated Fourier series using (1), results in

$$[eqn]$$

Using trigonometric identities, this is further written as follows:

$$[eqn]$$

Next, the frequency responses of the BP filter components of the bank are given as

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

As can be noticed from the expressions of frequency responses given by (5)–(11), the filters of the prototype filter bank with symmetric peak frequencies with respect to the value $[eqn]$ (at the middle of the interval) have symmetric zeros, which is according to digital filter theory. For instance, the zeros of the fifth BP filter are the zeros of the first BP filter with a changed sign, the second and forth BP filter also have opposite zeros, etc. As the FB is made up of an odd number of filters, the middle (third) BP filter, with peak frequency $[eqn]$ , has no pair, and its transfer function has pairs of complementary zeros. The uniform Gaussian prototype filter bank designed using the above procedure is plotted in Figure 1b and it looks very similar to its ideal counterpart in Figure 1a, with a very low level of ripple. In Table 1, the maximum ripple in the stop-band for all the filters of the filter bank are given. Due to the Gaussian shape, these filters have no ripple in the pass-band.

2.3. Design of 2D FIR Filters with Elliptically Shaped Frequency Response Using Frequency Transformation

In this work, we approached the analytical design of a particular class of anisotropic (directional) filters, with elliptically shaped support in the frequency plane. The shape of their characteristics is determined by specifying the values of the major and minor semi-axes of the ellipse (E and F). Their orientation is given by angle $[eqn]$ , formed by its longitudinal axis with the $[eqn]$ -axis of the frequency plane. In the following paragraphs, we determine the expression of the 1D to 2D frequency transformation that leads to 2D elliptically shaped filters.

According to the proposed analytical design technique, a 2D filter with an elliptically shaped frequency response is obtained by applying to the specified 1D prototype the following 1D to 2D frequency transformation $[eqn]$ , with the expression

$$[eqn]$$

In the above mapping, E and F are the values of major and minor semi-axes of the ellipse, respectively. The expressions of the coefficients $[eqn]$ , $[eqn]$ and $[eqn]$ are found by direct identification. In our case, applying this frequency mapping to a Gaussian filter with the specified peak frequency and bandwidth, we derive a 2D filter with elliptically shaped support in the frequency plane, depending on the parameters E, F and $[eqn]$ . These define the aspect of the ellipse and its orientation through the mapping (12).

A more convenient form of the frequency mapping (12) results using the simple identity $[eqn]$ ; making this substitution in (12) and re-arranging, we obtain

$$[eqn]$$

Next, introducing the notations $[eqn]$ , $[eqn]$ , we derive the following expressions of the coefficients a, b and c:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Applying the general 1D to 2D mapping (13) to the ideal Gaussian prototype $[eqn]$ , we obtain the ideal 2D frequency response:

$$[eqn]$$

The Gaussian 1D prototype, with specified selectivity, was approximated as a factored trigonometric polynomial in variable $[eqn]$ as in Relations (5)–(11). We will next determine the following frequency mapping, which is derived from (13):

$$[eqn]$$

The final purpose is to find a 1D to 2D mapping of the form $[eqn]$ , which can lead directly to the transfer function of a 2D digital filter with the desired shape.

In this step, the efficient Chebyshev series approximation is used, since it has the major advantage of yielding a precise and uniform approximation of a given function on a specified range of values. Using the symbolic calculation software MAPLE 6.02 and an intermediate change in variable, we derive the following accurate approximation (for $[eqn]$ , specifically $[eqn]$ ):

$$[eqn]$$

This approximation is plotted in Figure 2a in blue and is practically superposed on the exact function $[eqn]$ , which shows that the approximation is very precise.

Considering Equation (18) as well, we obtain

$$[eqn]$$

Using again the Chebyshev series method, the following low-order trigonometric approximation for the function $[eqn]$ results in

$$[eqn]$$

This approximation is plotted as the blue curve in Figure 2b, while the exact function $[eqn]$ is plotted in red. It can be noticed that the approximation is very accurate at least on the interval $[eqn]$ , and it tends to diverge visibly towards the margins of the interval.

Regarding this marginal error, even if it is large towards the limits of the interval, it does not affect the designed filters much since the component filters of the bank are derived from a narrow-band low-pass filter, and the marginal errors cancel out. Therefore, it is essential that the approximation is very precise on a sufficiently large interval around zero, as the plot in Figure 2b shows.

The above approximation is then written separately for the frequency variables $[eqn]$ and $[eqn]$ and their sum $[eqn]$ as follows:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

The next design step is to substitute Relations (22)–(24) into (13), obtaining an expression of $[eqn]$ which is then introduced in (20); finally, a trigonometric expression of the function $[eqn]$ in the variables $[eqn]$ and $[eqn]$ is determined. Considering this two-variable function as a 2D discrete Fourier transform, by directly identifying the corresponding coefficients, we finally reach a matrix of size $[eqn]$ , with the general form, symmetric with respect to the central element:

$$[eqn]$$

where the matrix elements have the following expressions, which depend on the ellipse semi-axes parameters E and F and on the angle of orientation, $[eqn]$ :

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Because the parametric matrix, $[eqn]$ , in (25) corresponds to the function $[eqn]$ , then taking into account Equation (20), the function $[eqn]$ will in turn correspond to matrix $[eqn]$ of size $[eqn]$ , according to the following relation:

$$[eqn]$$

where by the symbol $[eqn]$ , we denote the matrix convolution operation. In the above relation, $[eqn]$ is a $[eqn]$ null matrix, with a central element of value one. The matrix $[eqn]$ is obtained by bordering $[eqn]$ ( $[eqn]$ ) with zero values up to size $[eqn]$ .

Taking into account Expressions (5)–(11) for the factored frequency responses of the 1D prototype filter bank components, we can express them in the following form, which is, in general, a product of first-order and second-order factors:

$$[eqn]$$

In (34), the coefficients $[eqn]$ , $[eqn]$ and $[eqn]$ are generic coefficients which result in each factor in the prototype frequency response $[eqn]$ , while k is the constant resulting in front.

In the next step, simply by substituting in each factor of the general frequency response (34) the variable $[eqn]$ with the function $[eqn]$ derived as before, we find directly the factored frequency response of any 2D elliptically filter, generally denoted as $[eqn]$ and containing both first-order and second-order factors:

$$[eqn]$$

where S is a concise notation for the two-variable function $[eqn]$ . Thus, the large kernel H corresponding to the filter function $[eqn]$ can be expressed directly as a convolution product of smaller matrices, corresponding to each factor from (35):

$$[eqn]$$

Matrix Expression (36) corresponds to the factored frequency response (35). Using matrix $[eqn]$ of size $[eqn]$ given in (25) and taking into account Equations (22)–(24), each of the matrices $[eqn]$ from (36) is derived by adding coefficient $[eqn]$ occurring in first-order factors in (34) to the center element in matrix $[eqn]$ . Matrix $[eqn]$ , corresponding to second-order factors, is

$$[eqn]$$

where $[eqn]$ is a null matrix of size $[eqn]$ with a central element of value one. Matrix $[eqn]$ results from bordering matrix $[eqn]$ (size $[eqn]$ ) with zeros, symmetrically on all sides; here, the symbol * denotes matrix convolution.

Another advantage of the proposed method is that it yields the 2D filter frequency response directly factored, which may be a major advantage in a cascade structure (sequential) implementation of these filters. Some design examples of 2D elliptical filters are presented as follows.

2.4. Design Examples of 2D FIR Filters with Elliptical Shape

As a first design example using the proposed analytical technique, we have applied the frequency transformation determined above to each component of the prototype filter bank of a Gaussian shape, with the frequency responses given by Expressions (6)–(11). For the aspect ratio (between the semi-axes of the ellipse), we considered the value $[eqn]$ and the orientation angle of value $[eqn]$ . It can be noticed from Figure 3 that the resulting elliptical filters have a very accurate elliptical shape, steep transition and negligible ripple in the stop-band region.

If the aspect ratio of the ellipse (the ratio of semi-axes) $[eqn]$ is given a larger value, a very stretched ellipse is obtained, and the corresponding 2D elliptically shaped filters will be very narrow, with a very selective directional characteristic. In Figure 4, the frequency responses and contour plots were displayed for four low-pass very selective directional filters. They have the semi-axes values $[eqn]$ and $[eqn]$ and four orientation angles, namely $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ . As will be shown, such filters can be applied in detecting straight lines and other oriented objects in images.

The next section presents two important applications of the designed 2D elliptical filters, namely the analysis of a test image by the filter bank and the detection of oriented straight lines from an image, using very selective directional filters.

3. Applications and Simulation Results

This section presents a simulation example of the decomposition and reconstruction of a real-life test image using the designed elliptical uniform filter bank and a quantitative evaluation of the reconstruction quality, using specific measures. Then, a simulation example of directional filtering using very selective low-pass elliptical filters, with various orientation angles in the frequency plane, is given. This shows the capability of these narrow filters in extracting oriented straight lines from images.

3.1. Image Analysis Using the Designed Elliptical Uniform Filter Bank

In this section, we show a useful application of the uniform elliptically shaped FIR filter bank designed through the proposed analytical technique. The following simulation will be performed on a grayscale test image of size $[eqn]$ shown in Figure 5a, representing some branches of trees without foliage in winter time. This particular test image was chosen due to its structure of fine details, represented by the tree ramifications. The image was filtered by passing it through the six components of the filter bank (one low-pass filter and five band-pass filters). In this case, the elliptical uniform filter bank actually decomposes the test image into sub-band images, thus behaving as a usual analysis filter bank. The original image shown in Figure 5a is first low-pass filtered by the filter $[eqn]$ , which has a relatively narrow bandwidth, both transversally and longitudinally, for the specified values of the ellipse semi-axes. Figure 5b showing the output of the LP filter is very blurred, and the fine details, like thin twigs, can be no longer distinguished. Next, Figure 5c–f showing the outputs of the first four band-pass filters contain specific details, with higher and higher frequencies, corresponding to the selected region from the frequency plane.

As in every synthesis filter bank, the original test image that was decomposed by the analysis FB can be reconstructed by adding the sub-band component images. For instance, (g) is created by summing the first two component (b) and (c). Next, (h) is obtained through the summation of the first three components, namely (b–d), and so on. It can be observed that even the partially reconstructed (h) has regained in some measure its finer details. By summing all the six-component images, the reconstructed (i) will be obtained. It can be noticed that (i) is visually, from a subjective point of view, quite close to the original image, as it displays even the finest details quite clearly. This proves experimentally that the designed elliptically shaped uniform filter bank practically behaves and can be used as analysis filter bank for decomposing images into sub-band components.

Next, a quantitative analysis of the quality of image decomposition and reconstruction is performed using specific measures.

3.2. Quantitative Evaluation of Filtering Performance Using Relative Energy, RMSE and PSNR

In order to evaluate the relevance of image decomposition using the designed elliptical uniform filter bank, we can also evaluate the energy of each component sub-band image using the well-known formula $[eqn]$ , where the image is of size M × N. We will use the more useful expression of the relative energy as a percentage:

$$[eqn]$$

Calculating the energies of the seven filtered images created as the output of the designed filter bank, the values given in Table 2 are obtained; summing these values, we find out that they add up to 93.13%. It can also be noticed that about 53% of the image energy is contained in the low-pass component (in the frequency domain around zero), while almost 82% is contained in the first four components at the output of LP filter and first three BP filters. As can be seen from the bar plot in Figure 6, the relative energies of the sub-band images decrease monotonically.

In this way, the decomposition of the test image into sub-band images can also be described in terms of image energy. For the designed elliptical filter bank, the monotonical decrease in the energies of the sub-band images, the highest energy corresponding to the LP filter, explains the fact that the original image can be gradually reconstructed by successively summing up the sub-band images from low to high frequencies and that the higher frequency components have the least effect on the recovered image.

The simulation results of the image reconstruction for its component sub-band images show that the original test image is recovered with more and more precision by adding up the components, from the low-pass component to the high-pass one. Even if visually the reconstructed images seem more and more similar to the original one, it would be useful to evaluate their performance in a quantitative manner as well, using some objective metrics commonly applied in image processing.

The most appropriate measures for our purpose would be the root mean squared error (RMSE) and, derived from it, the peak signal-to-noise ratio (PSNR), which is the ratio between the maximum value or power of a signal and the power of noise that affects the signal. In our case, the images are not actually affected by noise. Therefore, we consider as “noise” the difference or error between the intermediate reconstructed image and the original image. The larger the value of the PSNR, the more similar the reconstructed image is to the original image, and the smaller the reconstruction error is.

Because, usually, signals or images have a wide dynamic range, the PSNR is most commonly expressed as a logarithmic quantity in decibels (dB).

Given a grayscale image, I, unaffected by noise, and its noisy or distorted version, J, the RMSE is defined mathematically as

$$[eqn]$$

where by $[eqn]$ and $[eqn]$ we denote the value of the current pixel (i,j) from the two images, respectively. Based on the root mean squared error, the PSNR (in decibels) is defined by the expression

$$[eqn]$$

In the above formula, $[eqn]$ is the maximum possible pixel value of the image. In the most common situation, when pixels are represented on 8 bits, this value is 255. Thus, in calculating the RMSE and PSNR, noisy image $[eqn]$ in (39) will be assimilated to the intermediate reconstructed image, with a given degree of distortion (in our case, blurring) compared to the ideal image, $[eqn]$ .

The gradually reconstructed images are denoted as follows, Ir1_2, Ir1_3, Ir1_4, Ir1_5, Ir1_6, and Ir1_7, where Ir1_2 is the image obtained by summing the first two components (LP and BP1), Ir1_3 is the image obtained by summing the first three components, and so on, while Ir1_7 is obtained by summing all seven components. In Table 3, the resulting values for the RMSE and PSNR, using Expressions (39) and (40), are displayed. Therefore, the PSNR obtained by adding up all component sub-band images is 21.79 dB. In future work, we will try to improve this figure by a better design and optimization of the 2D elliptical filter bank.

3.3. Application of Selective Directional Filtering in Oriented Straight Line Detection

The following set of simulations highlight the effect of very selective directional filtering on images containing oriented details. As is well known in image processing, the spectrum of a straight line from an image also has the shape of a straight line in the frequency plane $[eqn]$ , passing through the origin and orthogonal to the line direction.

Let us consider the image of size 690 × 690 pixels shown in Figure 7a, representing a structure of metallic beams, seen in perspective. This image is relevant for our directional filtering example, as it contains many straight lines oriented in various directions.

The logarithmic FFT spectrum of this image is displayed in Figure 7b. It has an interesting structure, consisting of fine straight lines of various orientations, corresponding to the various straight lines in the image. Taking into account the spectrum structure, the directional filtering can be easily understood. Thus, for a specified orientation angle of the directional filter, only the lines whose spectra overlap along the axis with the narrow filter characteristic will be preserved unaltered in the output image. The lines with other orientation will be more or less blurred. Lines with the spectrum perpendicular to the filter axis will be practically deleted due to blurring.

Applying a very selective Gaussian elliptical filter, like the ones shown in Figure 4, with various orientation angles (as specified in the figure), the filtered images in Figure 7c–g are obtained. For a specified orientation of the filter, different straight lines (various edges of metallic beams and other linear details) are outlined, while others are more or less affected by blurring, depending on the angle of orientation.

In the following section, an efficient implementation at the system level is proposed.

4. Polyphase Implementation of Designed 2D Elliptical FIR Filters

In the following, an efficient implementation is proposed at the system level for the designed 2D FIR elliptically shaped filters, relying on a polyphase decomposition of a 2D filtering operation with a convolution kernel of a relatively large size. In order to handle convolution with such a large kernel, block processing [21,22] and a polyphase decomposition approach are used in this work.

Firstly, using sub-expression sharing techniques, we obtain a fast algorithm that computes four outputs with a three-tap FIR filter for a 3 × 4 filtering operation, with a minimum number of multiplications. It has been shown that the minimal filtering algorithm for computing m outputs with an r-tap FIR filter, which we call F(m, r), requires µ(F(m, r)) = m + r − 1 multiplications [24]. Then, we can nest minimal 1D algorithms to form minimal 2D algorithms for computing 4 × 4 outputs with an 3 × 3 kernel filter [25]. These require the design of an efficient 2D filtering algorithm with a 3 × 3 kernel, but which can compute 4 × 4 output samples in parallel. In this section, a set of matrix equations describing the proposed novel filtering algorithm is presented. The obtained algorithm has been verified in MATLAB, version R2017a.

Thus, two partial output images, denoted as $[eqn]$ and $[eqn]$ , given by Equations (41) and (42), are obtained:

$$[eqn]$$ $$[eqn]$$

In the above equations, $[eqn]$ means the reversed version of vector $[eqn]$ , in which its elements are flipped (mirrored) horizontally.

Using the previous partial results, $[eqn]$ and $[eqn]$ , we next compute the output vector containing 16 samples of the 2D filtering operation using the following equation:

$$[eqn]$$

where the output samples are given by

$$[eqn]$$

and the block matrices $[eqn]$ , $[eqn]$ and $[eqn]$ have the following structure:

$$[eqn]$$

and $[eqn]$ , $[eqn]$ , $[eqn]$ are zero matrices of size $[eqn]$ , $[eqn]$ , and $[eqn]$ , respectively.

Also, we have vector H of the form

$$[eqn]$$

while input vector $[eqn]$ has the following structure:

$$[eqn]$$

The above algorithm for a 2D FIR filtering operation has been derived in order to reduce the number of operations, especially multiplications, due to the fact that there are redundant operations in the direct implementation of a 2D FIR filter. Using a new 2D algorithm, we significantly reduced the number of the multiplications not only compared with the direct implementation of a 2D FIR filter but also as compared with a recently proposed algorithm [23]. Using an efficient subexpression sharing technique, it was possible to obtain an increased reduction in the arithmetic complexity especially of the number of multiplications, as can be seen in the above algorithm for a 2D FIR filtering operation.

Using a block polyphase decomposition and the previous fast algorithm, we can obtain an efficient algorithm for the computation of the 2D FIR filter with a large kernel. In order to extend 2D filtering operations at a larger kernel, a block processing approach using polyphase decomposition will be used. In order to obtain 2D polyphase implementation using the proposed 2D algorithm from Equations (41)–(43), a decimation factor of 3 was used. With this value, a 12 × 12 input matrix was obtained.

For an easier understanding of our method, our discussion was restricted to a particular case where the kernel matrix was 9 × 9 and the input matrix was 12 × 12. However, the obtained results can be easily extended for our designed filter matrix. In order to use our previous algorithm for a 2D filter, a decimation factor of 3 was applied, both for the kernel matrix and for the input matrix.

Equations (41)–(43) can now be reformulated for a polyphase block implementation as will be shown below. Vector h is replaced with a new vector, defined as

$$[eqn]$$

Vectors $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ for the case where the kernel matrix is 9 × 9 and the input matrix is 12 × 12 are defined below.

First, vector $[eqn]$ has the following generic, parameterized form:

$$[eqn]$$

from which we derive nine vectors obtained for i, j = 0, 1, 2; for instance, the following vectors are derived for the indicated indices i and j:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Applying decimation with 3 instead of the kernel matrix of size 9 × 9, we obtain nine matrices of size 3 × 3. For example, for $[eqn]$ using decimation with factor 3, the following 3 × 3 block matrix was obtained:

$$[eqn]$$

We take the rows of matrix $[eqn]$ and concatenate them, and then by reversing the order of the elements, vector $[eqn]$ is obtained:

$$[eqn]$$

Vector $[eqn]$ is also replaced with vector $[eqn]$ , defined as

$$[eqn]$$

Vectors $[eqn]$ , have the following generic, parameterized form:

$$[eqn]$$

From this general form, by giving successive values (0, 1, 2,…,5) to indices i and j, all 36 vectors are derived. For instance, the following vectors are obtained for the indicated index values:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Applying decimation with a factor 3, instead of the input matrix of 12 × 12, we obtain nine matrices of size 3 × 3. As an example, for vector X21, using decimation with a factor of 3, the following 3 × 3 block matrix results:

$$[eqn]$$

At this step, we take the rows of matrix $[eqn]$ , concatenate them, and then by reversing the order of the elements, the following matrix $[eqn]$ is obtained:

$$[eqn]$$

Thus, we decomposed a 2D FIR filtering operation with a 9 × 9 kernel into 49 inner products, as shown below.

In the preprocessing stage, it is first necessary to obtain, from a 12 × 12 input matrix, matrix $[eqn]$ of size 3 × 3, using a decimation operation. Then, by taking the rows of matrix $[eqn]$ , vectors $[eqn]$ will be obtained, as has been shown previously for matrix $[eqn]$ .

Using polyphase decomposition and the obtained equations for the 2D FIR filter (41)–(43), the following equations are used for the polyphase implementation of our 2D FIR filter with a 9 × 9 kernel:

$$[eqn]$$

where $[eqn]$ is a null matrix of size 4 × 7, $[eqn]$ , and also, we have the following vectors:

$$[eqn]$$ $$[eqn]$$

In order to obtain the final output samples, a postprocessing stage is required, in which the partial results $[eqn]$ and $[eqn]$ are added.

To determine the elements of vector $[eqn]$ , we first compute

$$[eqn]$$ $$[eqn]$$

Then, vectors $[eqn]$ and $[eqn]$ are split into smaller sections of nine elements, and the elements of auxiliary vector $[eqn]$ are computed using an inner product operation, as follows:

$$[eqn]$$

where “ $[eqn]$ ” denotes the inner product operation.

In a similar way, the elements of auxiliary vector $[eqn]$ are computed by splitting the two vectors, $[eqn]$ and $[eqn]$ , into smaller sections of nine elements, using an inner product operation. Vectors $[eqn]$ and $[eqn]$ result from the following equations:

$$[eqn]$$ $$[eqn]$$

with $[eqn]$ and $[eqn]$ , where matrix $[eqn]$ is the $[eqn]$ identity matrix (with ones on the main diagonal and zeros elsewhere), with the structure

$$[eqn]$$

and $[eqn]$ and $[eqn]$ are zero matrices of sizes 63 × 54 and 63 × 27, respectively. Then, we can compute, as before, the intermediate values:

$$[eqn]$$

Also, as can be seen from Equations (68) and (72), in order to introduce polyphase decomposition, it was necessary to introduce a post-processing block that performs a segmentation of vectors $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ in segments of nine elements that are used in 49 inner product operations.

The two auxilliary vectors, $[eqn]$ and $[eqn]$ , were introduced to facilitate the reformulation of the algorithm for 2D FIR filtering, with polyphase implementation.

This polyphase algorithm was developed for a 2D FIR filtering operation, with the aim to reduce computational complexity, since it is well known that there are redundant operations in the direct implementation of 2D convolution. In the direct approach, there are, for instance, overlapping input data blocks. For this reason, eliminating redundancy as much as possible is essential for reducing the number of arithmetic operations. The novel approach proposed here, using a subexpression sharing technique, allows for an efficient implementation structure with a high degree of parallelism and therefore an increased speed of operation.

5. Discussion

The proposed design method for 2D elliptically shaped filters is entirely analytical, without using any global numerical optimization techniques. Analytical design techniques lead to closed-form and parametric filters, with adjustable, tunable frequency responses. To the best of the authors’ knowledge, analytical designs for 2D FIR elliptical filter banks have not been previously studied. For this reason, a rigorous comparison in terms of performance with other elliptical filters found in the literature, like [14,15,16,17,18], is quite difficult to make since the design approaches, filter characteristics and purposes are very different.

For instance, paper [18] approaches the design and optimization of the performance of 2D FIR filters with circular and elliptical symmetry. However, it differs significantly from our analytical design method since it uses a global optimization approach, based on a particular, very efficient meta-heuristic algorithm, namely the Salp-Swarm optimization algorithm. As in every numerical optimization algorithm, it starts from an initial set of random parameters, regarded as a population, which is updated over a number of iterations until some fitness condition is met, specifically an imposed maximum and average ripple in the pass-band and stop-band. The 2D FIR filter kernel resulting from the optimization process cannot be, in general, decomposed and the frequency response cannot be factored, as in our analytical method.

For the elliptically shaped filter, the authors have chosen some fixed semi-axes values (the aspect ratio) and the ellipse axes are parallel with the frequency plane axes. Also, in [18], only maximally flat, low-pass circular and elliptical filters with a wide band have been designed. By comparison, in our paper, we have designed narrow-band elliptically shaped filters based on a Gaussian prototype (of the LP, BP and HP type, forming an uniform filter bank) and they are oriented at a specified angle in the frequency plane. In particular, we have designed an elliptical Gaussian filter bank based on a uniform 1D prototype with specified peak frequencies and bandwidths.

However, the filters can be roughly compared regarding the maximum ripple in the stop-band. Due to the specific Gaussian shape, our filters have no ripple in the pass-band. The stop-band ripple of the derived 2D filters would be rather difficult to measure exactly, but in principle, its maximum values are approximately the ones of their prototypes, given in Table 1 from Section 2. It can be noticed that these maximum stop-band ripple values are lower than the ones of the elliptical filters designed in [18].

The elliptically shaped filter bank designed in this work can be regarded as a continuation of the work [23], in which circular 2D FIR filter banks were designed. However, the design technique is more complex for elliptical filters, the expression of the frequency transformation being more complicated than in the circular case. As with circular filter banks, we have analyzed a test image, decomposing and then reconstructing it relatively accurate from its sub-band components. This suggests the eventual possibility of using such filters in image coding. Referring to existing works on related topics, other analytical techniques for designing 2D elliptically shaped filters of the IIR type were previously proposed in [24,25].

For this application, we chose the Gaussian filter as a prototype for the 2D filter bank, due to its well-known advantages. Among others, it can be easily and efficiently approximated and can also be scaled on the frequency axis to modify its selectivity. Because its frequency response is of a zero phase, the frequency components of an image will not be phase-shifted, and thus, the filter will not produce image distortions. The resulting filters have accurate shapes, with negligible distortions.

The proposed solution was based on a new efficient algorithm for a 2D FIR filtering operation with lower arithmetic complexity for a smaller kernel of 3 × 3, which was then extended at a larger kernel size of 30 × 30 using a polyphase decomposition of the given 2D filtering operation and block filtering.

In order to facilitate the reformulation of the algorithm for 2D FIR filtering for a polyphase implementation, we introduced some auxiliary vectors. So, it was possible to extend the proposed algorithm 2D algorithm for a significantly larger size.

Using the proposed algorithm we have significantly reduced the complexity of the algorithm while maintaining the same throughput of 16 output samples at a time, as in [23]. Thus, instead of 100 inner products as in [23], we have only 49 inner products in the proposed algorithm. In the simpler case that was used, where the kernel size was 9 × 9, each inner product involved nine multiplications. Thus, instead of 900 multiplications as in [23], we now have only 49 × 9 = 441 multiplications. This advantage is even more important in the case of hardware implementation.

6. Conclusions

The proposed design technique for 2D elliptically shaped filter banks is fully analytical, without using any complicated numerical optimization algorithms. It is mainly based on a 1D to 2D frequency transformation which is especially developed using various approximations and is applied to the 1D prototype, yielding directly a factored 2D frequency response, which may be useful in a sequential implementation. Equivalently, the overall filter kernel results were directly decomposed as a convolution product of smaller matrices. The major advantage of this novel design approach is that the resulting filter is parametric. The filter specifications (the selectivity and orientation angle) appear explicitly in the filter kernel and the filter is adjustable. Thus, for any given specifications, the filter is determined directly, and there is no need to resume every time the entire design process. This approach leads to filters with accurate shapes and low distortions. As shown through simulations, the uniform elliptical filter bank decomposes an image into sub-band components, with increasing frequency content and decreasing energy.

Then, by adding these sub-band images one by one, starting from the low-frequency component, the initial image is gradually reconstructed. By summing all the sub-band images, the reconstructed image is visually quite close to the original test image, the finest details being quite clearly visible. The reconstruction quality was evaluated numerically using measures such as the RMSE and PSNR.

Regarding the implementation part, the proposed solution based on the polyphase and block filtering approach is very efficient and leads to a filter structure with reduced arithmetical complexity. This is true compared not only to the direct filtering approach (which can be very complex for large filter kernels and large images) but also to the most recent similar solution presented in [23].

In a future continuation of this work, we will study more rigorously the design of 2D elliptical filter banks and their applications in image analysis, in particular, focusing on the improvement of image reconstruction quality.