Open access

Analytical Design of Two-Dimensional Filters and Applications in Biomedical Image Processing

Written By

Radu Matei and Daniela Matei

Submitted: 06 April 2012 Published: 16 January 2013

DOI: 10.5772/52195

From the Edited Volume

Digital Filters and Signal Processing

Edited by Fausto Pedro García Márquez and Noor Zaman

Chapter metrics overview

3,337 Chapter Downloads

View Full Metrics

1. Introduction

The field of two-dimensional filters and their design methods has known a large development due to its importance in image processing (Lim, 1990; Lu & Antoniou, 1992). There are methods based on numerical optimization and also analytical methods relying on 1D prototypes. A commonly-used design technique for 2D filters is to start from a specified 1D prototype filter and transform its transfer function using various frequency mappings in order to obtain a 2D filter with a desired frequency response. These are essentially spectral transformations from s to z plane, followed by z to (z1,z2) mappings, approached in early papers (Chakrabarti & Mitra, 1977; Hirano & Aggarwal, 1978; Harn & Shenoi, 1986; Nie & Unbehauen, 1989). Generally these transformations conserve stability, so from 1D prototypes various stable recursive 2D filters can be obtained. The most common types are directional, fan-shaped, diamond-shaped and circular filters. Diamond filters are commonly used as anti-aliasing filters in the conversion between signals sampled on the rectangular sampling grid and the quincunx sampling grid. Various design methods for diamond-shaped filters were studied in (Tosic, 1997; Lim & Low, 1997; Low & Lim, 1998; Ito, 2010; Matei, 2010).

There are several classes of filters with orientation-selective frequency response, useful in tasks like edge detection, motion analysis, texture segmentation etc. Some relevant papers on directional filters and their applications are (Danielsson, 1980; Paplinski, 1998; Austvoll, 2000). An important class of orientation-selective filters are steerable filters, synthesized as a linear combination of a set of basis filters (Freeman & Adelson, 1991) and steerable wedge filters (Simoncelli & Farid, 1996). A directional filter bank (DFB) for image decomposition in the frequency domain was proposed in (Bamberger, 1992). In (Qunshan & Swamy, 1994) various 2D recursive filters are approached. Fan-shaped, also known as wedge-shaped filters find interesting applications. Design methods for IIR and FIR fan filters are presented in some early papers (Kayran & King, 1983; Ansari, 1987). An efficient design method for recursive fan filters is presented in (Zhu & Zhenya, 1990). An implementation of recursive fan filters using all-pass sections is given in (Zhu & Nakamura, 1996). In (Mollova, 1997), an analytical least-squares technique for FIR filters, in particular fan-type, is proposed. Design methods for efficient 2D FIR filters were treated in papers like (Zhu et al., 1999; Zhu et al., 2006). Zero-phase filters were studied as well (Psarakis, 1990). Different types of 2D filters derived from 1D prototypes through spectral transformations were treated in (Matei, 2011a).

We propose in this chapter some new design procedures for particular classes of 2D filters; the described methods are mainly analytical but also include numerical approximations. Various types of 2D filters will be approached, both recursive (IIR) and non-recursive (FIR). The design methods will focus however on recursive filters, since they are the most efficient.

The proposed design methods start from either digital or analog 1D prototypes with a desired characteristic. In this chapter we will mainly use analog prototypes, since the design turns out to be simpler and the 2D filters result of lower complexity. This analog prototype filter is described by a transfer function in the complex variable s, which can be factorized as a product of elementary functions of first or second order. The prototype transfer function results from an usual approximation (Butterworth, Chebyshev, elliptic) and the shape of the frequency response corresponds to the desired characteristic of the 2D filter.

The next design stage consists in finding the specific complex frequency transformation from the axis s to the complex plane (z1,z2), of the general form F:2, sF(z1,z2). This mapping will be determined for each type of 2D filter separately, starting from the geometrical specification of its shape in the frequency plane. Once found this particular mapping, the 2D filter function results directly by applying this transformation to each factor function of the prototype. Thus, the 2D filter transfer function H(z1,z2) results directly factorized, which is a major advantage in its implementation. The proposed design method applies the bilinear transform as an intermediate step in determining the 1D to 2D frequency mapping. All the proposed design techniques are mainly analytical but also involve numerical optimization, in particular rational approximations (e.g. Chebyshev-Padé). Some of the designed 2D filters result with complex coefficients. This should not be a serious shortcoming, since such IIR filters are also used (Nikolova et al., 2011).

In this chapter we will approach two main classes of 2D filters. The first one comprises three types of orientation-selective filters, as follows: square-shaped (diamond-type) IIR filters, with arbitrary orientation in the frequency plane; fan-type IIR filters with specified orientation and aperture angles; and very selective IIR multi-directional filters (in particular two-directional and three-directional), which are useful in detecting and extracting simultaneously lines with different orientations from an image.

The other class discussed here refers to FIR filters. From this category we will approach zero-phase filters with circular frequency response. Zero-phase filters, with real transfer functions, are often used in image processing since they do not introduce any phase distortions. All these types of 2D filters are analyzed in detail in the following sections.

Stability of the two-dimensional recursive filters is also an important issue and is much more complicated than for 1D filters. For 2D filters, in general, it is quite difficult to take stability constraints into account during approximation stage (O’Connor & Huang, 1978). Therefore, various techniques were developed to separate stability from approximation. If the designed filter becomes unstable, some stabilization procedures are needed (Jury et al., 1977). Various stability conditions for 2D filters have been found (Mastorakis, 2000).

The medical image processing field has known a rapid development due to imaging value in assisting and assessing clinical diagnosis (Semmlow, 2004; Berry, 2007; Dougherty, 2011). In particular, the currently used vascular imaging technique is x-ray angiography, mainly in diagnosing cardio-vascular pathologies. A frequent application of cardiac imaging is the localization of narrowed or blocked coronary arteries. Fluorescein angiography is the best technique to view the retinal circulation and is useful for diagnosing retinal or optic nerve condition and assessing disorders like diabetic retinopathy, macular degeneration, retinal vein occlusions etc. There are many papers approaching various methods and techniques aiming at improving angiogram images. In papers like (Frangi et al., 1998) the multiscale analysis is used, with the purpose of vessel enhancement and detection. Usual approaches include Hessian-based filtering, based on the multiscale local structure of an image and directional features of vessels (Truc et al., 2007). In cardio-vascular imaging, an essential pre-processing task is the enhancement of coronary arterial tree, commonly using gradient or other local operators. In (Khan et al., 2004) a decimation-free directional filter bank is used. An adaptive vessel detection scheme is proposed in (Wu et al., 2006) based on Gabor filter response. Filtering is an elementary operation in low level computer vision and a pre-processing stage in many biomedical image processing applications. Some edge-preserving filtering techniques for biomedical image smoothing have been proposed (Rydell et al., 2008; Wong et al., 2004). At the end of this chapter some simulation results are given for biomedical image filtering using some of the proposed 2D filters, namely the directional narrow fan-filter with specified orientation and the zero-phase circular filter.

Advertisement

2. Analog and digital 1D prototype filters used in 2D filter design

This section presents the types of analog and digital 1D recursive prototype filters which will be further used to derive the desired 2D filter characteristics. An analog IIR prototype filter of order N has a transfer function in variable s of the general form:

HP(s)=P(s)Q(s)=i=0Mpisi/j=0NqjsjE1

This general transfer function can be factorized into simpler rational functions of first and second order. Such a second-order rational function (biquad) can be written:

Hb(s)=k(s2+b1s+b0)/(s2+a1s+a0)E2

where generally the second-order polynomials at numerator and denominator have complex-conjugated roots, and k is a constant. For typical approximations – Chebyshev or elliptic – usually b1=0, therefore the numerator has imaginary zeros. For odd-order filters, the denominator contains at least a first-order factor (s+α). An elliptic approximation with very low ripple can be used for an almost maximally-flat low-order filter. Next we consider two such low-pass (LP) prototypes with imposed specifications. The first is an elliptic LP analog filter of order N=6, cutoff frequency ωc=0.4π, peak-to-peak ripple Rp=0.04dB, stop-band attenuation Rs=38dB. Its transfer function can be factorized into three biquad functions like (2): HP(s)=kHb1(s)Hb2(s)Hb3(s) where k=2.375 and:

Hb1(s)=(s2+39.195)/(s2+0.2221s+2.8797)E3
Hb2(s)=(s2+6.5057)/(s2+0.9172s+2.4291)E4
Hb3(s)=(s2+4.2217)/(s2+2.0448s+1.5454)E5

The frequency response magnitude of this LP filter for ω[π,π] is shown in Fig. 1(a).

The second prototype is an elliptic LP analog filter with parameters: N=4, ωc=0.4π, Rp=0.05db,Rs=36db. Its transfer function is written as a product of two biquad functions like (2): HP(s)=kHb1(s)Hb2(s), where k0.01 and

Hb1(s)=(s2+33.385)/(s2+0.5894s+2.2398) E6
Hb2(s)=(s2+6.4226)/(s2+1.9691s+1.5266)E7

The frequency response magnitude of this LP filter for ω[π,π] is shown in Fig.1 (b). The simplest analog LP filter has a transfer function Hj(s)=α/(s+α), where the value α gives the selectivity (Fig.1(c)). If the filter characteristic is shifted to a given frequency ω01[π,π], the transfer function becomes:

HjS(s)=α/(s+α+jω01)E8

In Fig.1 (d) the shifted filter response magnitude for ω01=0.416π is shown. Another useful analog prototype is the selective second-order (resonant) filter with central frequency ω0:

Hr(s)=αs/(s2+αs+ω02)E9

The transfer function magnitude for such a filter with α=0.1 and ω0=1.3 is shown in Fig. 1 (e). This will be further used as a prototype for two-directional filters.

Figure 1.

Frequency response magnitudes of: (a) LP elliptic prototype of order 6; (b) LP elliptic prototype of order 4; very selective first-order filter with central frequencies ω0=0 (c) and ω0=0.416π (d); (e) selective band-pass filter with ω0=1.3

A useful zero-phase prototype can be obtained from the general function (1) by preserving only the magnitude characteristics of the 1D filter; this prototype will be further used to design 2D zero-phase FIR filters of different types, specifically circular filters, with real-valued transfer functions. In order to obtain a zero-phase filter, we consider the magnitude characteristics of HP(jω), defined by the absolute value |HP(jω)|=|P(jω)|/|Q(jω)|. We look for a series expansion of the magnitude |HP(jω)| that has to be an approximation as accurate as possible on the frequency domain [π,π]. The most convenient for our purpose is the Chebyshev series expansion, because it yields an efficient approximation of a given function, which is uniform along the desired interval. The Chebyshev series in powers of the frequency variable ω for a given function on a specified interval can be easily found using a symbolic computation software like MAPLE. However, we will finally need a trigonometric expansion of |HP(jω)|, namely in cos(nω), rather than a polynomial expansion in powers of ω. Therefore, prior to Chebyshev series calculation, we apply the change of variable:

ω=arccos(x)x=cos(ω)E10

and so we get the polynomial expansion in variable x:

|HP(arccos(x))|n=0Nanxn=a0+a0x+a2x2+a3x3+...+aNxNE11

where the number of terms N is chosen large enough to ensure the desired precision. The next step is to substitute back x=cosω in the polynomial expression (11), therefore we obtain the factorized function in cosω, with n+2m=N:

|HP(ω)|n=0Nancosn(ω)=ki=1n(cosω+ai)j=1m(cos2ω+a1jcosω+a2j)E12

Next let us consider a recursive digital filter of order N with the transfer function:

HP(z)=P(z)Q(z)=i=0Mpizi/j=0NqjzjE13

This general transfer function with M=N can be factorized into first and second order rational functions. For an odd order filter, HP(z) has at least one first-order factor:

H1(z)=(b1z+b0)/(z+a0)E14

The transfer function also contains second-order (biquad) functions, where in general the numerator and denominator polynomials have complex-conjugated roots:

H2(z)=(b2z2+b1z+b0)/(z2+a1z+a0)E15

We will further use the term template, common in the field of cellular neural networks, for the coefficient matrices of the numerator and denominator of a 2D transfer function H(z1,z2).

Advertisement

3. Diamond-type recursive filters

In this section a design method is proposed for 2D square-shaped (diamond-type) IIR filters. The design relies on an analog 1D maximally-flat low-pass prototype filter. To this filter a frequency transformation is applied, which yields a 2D filter with the desired square shape in the frequency plane. The proposed method combines the analytical approach with numerical approximations.

3.1. Specification of diamond-type filters in the frequency plane

The standard diamond filter has the shape in the frequency plane as shown in Fig.2 (a). It is a square with a side length of π2, while its axis is tilted by an angle of φ=π/4 radians about the two frequency axes. Next we will consider the orientation angle φ about the ω2 (vertical) axis. In this chapter a more general case is approached, i.e. a 2D diamond-type filter with a square shape in the frequency plane, but with arbitrary axis orientation angle, as shown in Fig. 2(e). Next we refer to them as diamond-type filters, since they are more general than the diamond filter from Fig. 2 (a).

The diamond-type filter in Fig.2 (e) is derived as the intersection of two oriented low-pass filters whose axes are perpendicular to each other, for which the shape in the frequency plane is given in Fig.2 (c), (d). Correspondingly, the diamond-type filter transfer function HD(z1,z2) results as a product of two partial transfer functions:

HD(z1,z2)=H1(z1,z2)H2(z1,z2)E16

The frequency characteristic of H2(z1,z2) is ideally identical to the frequency characteristic of H1(z1,z2) rotated by an angle of φ=π/2. Since this rotation of axes implies the frequency variable change: ω1ω2, ω2ω1, the transfer function H2(z1,z2) can be derived from H1(z1,z2) as H2(z1,z2)=H1(z2,z11). A more general filter belonging to this class is a rhomboidal filter, as shown in Fig.2 (f). In this case the two oriented LP filters may have different bandwidths and their axes are no longer perpendicular to each other.

Figure 2.

(a) diamond filter; (b) wide-band oriented filter; (c), (d) wide-band oriented filters with orientations forming an angle φ=π/2; (e) square-shaped filter resulted as product of the above oriented filters; (f) rhomboidal filter

3.2. Design method for diamond-type filters

The issue of this section is to find the transfer function H2D(z1,z2) of the desired 2D filter using a complex frequency transformation sF(z1,z2). From a prototype HP(s)=HP(jω) (which varies on one axis only), a 2D oriented filter is obtained by rotating the axes of the plane (ω1,ω2) by an angle φ. The rotation is defined by the following linear transformation, where ω1,ω2 are the original frequency variables and ω¯1,ω¯2 the rotated ones:

[ω1ω2]=[cosφsinφsinφcosφ][ω¯1ω¯2]E17

The spatial orientation is specified by an angle φ with respect to ω1axis, defined by the 1D to 2D frequency mapping ωω1cosφ+ω2sinφ. By substitution, we obtain the oriented filter transfer function Hφ(ω1,ω2)=HP(ω1cosφ+ω2sinφ). In the complex plane (s1,s2) the above frequency transformation becomes:

ss1cosφ+s2sinφE18

The oriented filter Hφ(ω1,ω2) has the frequency response magnitude section along the line ω1cosφ+ω2sinφ=0, identical with prototype HP(ω), and constant along the perpendicular line (filter longitudinal axis) ω1sinφω2cosφ=0. The usual method to obtain a discrete filter from an analog prototype is the bilinear transform. If the sample interval takes the value T=1, the bilinear transform for s1 and s2 in the complex plane (s1,s2) has the form:

s1=2(z11)/(z1+1)s2=2(z21)/(z2+1)E19

This method is straightforward, still the resulted 2D filter will present linearity distortions in its shape, which increase towards the limits of the frequency plane as compared to the ideal frequency response. This is mainly due to the so-called frequency warping effect of the bilinear transform, expressed by the continuous to discrete frequency mapping:

ω=(2/T)arctg(ωaT/2)E20

where ω is a frequency of the discrete filter and ωa is the corresponding frequency of the analog filter. This error can be corrected by applying a pre-warping. Taking T=1 in (20), we substitute the mappings:

ω12arctg(ω1/2)         ω22arctg(ω2/2)E21

In order to include the nonlinear mappings (21) into the frequency transformation, a rational approximation is needed. One of the most efficient is Chebyshev-Padé, which gives uniform approximation over a specified range. We get the accurate approximation for ω[π,π]:

arctg(ω/2)0.4751ω/(1+0.05ω2)E22

Substituting the nonlinear mappings (21) with approximate expression (22) into (18) we get the 1D to 2D mapping which includes the pre-warping along both frequency axes:

sFφ(s1,s2)=0.95[s1cosφ10.05s12+s2sinφ10.05s22]E23

Applying the bilinear transform (19) along the two axes we obtain the mapping sFφ(z1,z2) in matrix form, where z1=[1   z1   z12] and z2=[1   z2   z22]:

sFφ(z1,z2)=kMφ(z1,z2)/Nφ(z1,z2)=k(z1×Mφ×z2T)/(z1×Nφ×z2T)E24

Here k=1.5233 and the matrices Mφ and Nφ of size 3×3 are given by:

Mφ=cosφ[131000131]+sinφ[101303101]Nφ=[131393131]E25

Substituting the mapping (24) into the expression (2) of the biquad transfer function Hb(s) with b1=0, we get the 2D transfer function HB(z1,z2) in the matrix form:

HB(z1,z2)=(z1×B1×z2T)/(z1×A1×z2T)E26

where z1 and z2 are the vectors:

z1=[1z1z12z13z14]z2=[1z2z22z23z24]E27

and the 5×5 templates B1, A1 are given by the expressions:

B1=k2MφMφ+b0NφNφ ;    A1=k2MφMφ+a1kMφNφ+a0NφNφE28

For instance, corresponding to the third biquad function Hb3(s) given by (5), the following 5×5 templates result according to the expressions (28):

B1=[0.24640.94071.14180.40270.06710.94073.22333.89171.60920.40271.14183.89176.14843.89171.14180.40271.60923.89173.22330.94070.06710.40271.14180.94070.2464] A1=[0.09470.19410.07320.01630.02450.19410.27380.71810.07740.31120.07320.71811.00002.88511.27430.01630.07742.88513.65701.17680.02450.31121.27431.17680.3131]

Figure 3.

(a) LP correction filter characteristic; frequency response magnitudes and contour plots of: (b), (c) uncorrected diamond-type filter; (d), (e) corrected diamond-type filter

The characteristics of a diamond-type filter with orientation angle φ=π/4 and based on the prototype filter of order 6 given by the factors (3)-(5) is shown in Fig.3 (b), (c). As can be noticed, the filter characteristic corrected by pre-warping has a good linearity, however it still twists towards the margins of the frequency plane. These marginal linearity distortions can be corrected using an additional LP filter. For instance, we can choose as prototype an 1D elliptic digital filter of order N=3, pass-band ripple Rp=0.1dB, stop-band attenuation Rs=40 dB and cutoff frequency ωc=0.6, which has the coefficients given by the vectors:

BC1=[0.35131.011.010.3513]           AC1=[10.96440.67010.088]E29

The 2D low-pass filter is separable and results by applying successively the 1D filter along the two frequency axes; the 4×4 matrices of the correction filter result as: BC=BC1TBC1, AC=AC1TAC1, where the symbol denotes outer product of vectors. The correction filter has the following transfer function, where z1=[1   z1   z12   z13], z2=[1   z2   z22   z23]:

HC(z1,z2)=(z1×BC×z2T)/(z1×AC×z2T)E30

The resulted 2D square-shaped correction filter characteristic is shown in Fig.3 (a) and is almost maximally-flat, as required. The corrected version of the diamond-type filter from Fig. 3 (b), (c) has the magnitude and the contour plot shown in Fig. 3 (d), (e). It can be easily noticed that the initial distortions have been eliminated. Another two diamond-type filters with orientation angles φ=π/12 and φ=π/6 are shown in Fig. 4 (a)-(d).

Figure 4.

Diamond-type filters with orientation angle: (a), (b) φ=π/12; (c), (d) φ=π/6

Advertisement

4. Fan-type recursive filters

In this section an analytical design method in the frequency domain for 2D fan-type filters is proposed, starting from an 1D analog prototype filter, with a transfer function decomposed as a product of elementary functions. Since we envisage designing efficient 2D filters, of minimum order, recursive filters are used as prototypes, and the 2D fan-type filters will result recursive as well.

In Fig.5 (a) a general fan-type filter is shown, with an aperture angle BOD=θ, oriented along an axis CC' and its longitudinal axis forming an angle AOC=φ with frequency axis Oω2. A particular case is the two-quadrant fan filter, shown in Fig.5 (b). Fig.5 (c) shows a DFB with 8-band frequency partition (Bamberger, 1992), an angularly-oriented image decomposition which splits the frequency plane into fan-shaped sub-bands (channels).

Figure 5.

(a) Ideal fan filter with given aperture, oriented at an angle φ; (b) Ideal two-quadrant fan filter (c) 8-band partitions of the frequency plane

Figure 6.

Frequency response magnitudes and contour plots for: (a) fan-type filter with aperture θ=0.1π and orientation φ=π/7; corrected filters with θ=0.1π, φ=π/7 (b) and φ=0(c)

The 1D analog filter discussed in section 2 is used as prototype. The general fan-type filter can be derived from a LP prototype using the frequency mapping (Matei & Matei, 2012):

ωfφ(ω1,ω2)=a(ω1cosφω2sinφ)/(ω1sinφ+ω2cosφ)E31

In (31), a=1/tg(θ/2) is the aperture coefficient, where θ is the aperture angle of the fan-type filter. This frequency mapping in the complex variables s1=jω1, s2=jω2 is:

sfφ(s1,s2)=ja(s1cosφs2sinφ)/(s1sinφ+s2cosφ)E32

Applying the same steps as in Section 3.2 in order to obtain a discrete form of the above frequency mapping, using relations (21), (22) and (32) we obtain the 1D to 2D mapping which includes pre-warping along both axes of the frequency plane:

sFφ(s1,s2)=ja(s1(10.05s22)cosφs2(10.05s12)sinφ)(s1(10.05s22)sinφ+s2(10.05s12)cosφ)E33

We now apply the bilinear transform (19) along the two axes and obtain the mapping sFφ(z1,z2) in matrix form, where z1=[1   z1   z12] and z2=[1   z2   z22]:

sFφ(z1,z2)=jaPφ(z1,z2)/Qφ(z1,z2)=ja(z1×Pφ×z2T)/(z1×Qφ×z2T)E34

and the 3×3 matrices Pφ and Qφ are given by:

Pφ=cosφ[131000131]sinφ[101303101];Qφ=sinφ[131000131]+cosφ[101303101]E35

Substituting the mapping (34) into the biquad expression (2) with b1=0, we get the 2D transfer function in matrix form HW1(z1,z2)=(z1×B2×z2T)/(z1×A2×z2T), similar to (26), where the vectors z1, z2 are given by (27). The 5×5 templates B2 and A2 are given by:

B2=b0QφQφa2PφPφ ;     A2=a0QφQφa2PφPφ+jaa1PφQφE36

The 2D transfer function for each biquad is complex. The characteristics of a fan-type filter designed with this method and based on the prototype filter of order 4 given by (6)-(7) is shown in Fig.6 (a), for the indicated parameters. As with the diamond-type filter analyzed in the previous section, the fan-type filter characteristic features marginal linearity distortions which can be corrected using a LP filter, similar with the correction filter used in Section 3.2 and having the frequency characteristic shown in Fig. 3 (a).

Figure 7.

(a) Frequency response magnitude and (b) contour plot for the 2-quadrant fan filter

Two corrected fan-type filters with specified parameters have the magnitudes and contour plots shown in Fig.6 (b), (c). The initial distortions have been eliminated. With the same correction filter, we obtain the two-quadrant fan filter, shown in Fig. 7, by setting the aperture angle θ=π/2 and orientation angle φ=π/4.

Advertisement

5. Very selective multidirectional IIR Filters

In this section a design method based on spectral transformations is proposed for another class of 2D IIR filters, namely multi-directional filters. The design starts from an analog prototype with specified parameters. Applying an appropriate frequency transformation to the 1D transfer function, the desired 2D filter is directly obtained in a factorized form, like the filters designed in the previous sections. For two-directional filters, an example is given of extracting lines with two different orientations from a test image. The spectral transformation used in the case of multi-directional filters is similar to the one presented in the previous section, derived for fan-type filters and given by (34), (35). In this section the design of two-directional and three-directional filters with specified orientation is detailed. The method can be easily generalized to arbitrary multi-directional filters.

Figure 8.

Ideal shapes of some directional filters in the frequency plane: (a) two-directional filter; (b) three-directional filter; (c) three-band selective filter with ω01=0.59π, ω02=0.36π

5.1. Two-directional fan-type filters

A two-directional 2D filter is orientation-selective along two directions in the frequency plane. It is based on a selective resonant IIR prototype as given in section 2. Applying the same frequency transformation sFφ(z1,z2) derived for fan-type filters and given by (33) to the prototype filter (9) we get the 2D two-directional transfer function H2(z1,z2) in matrix form H2(z1,z2)=(z1×B3×z2T)/(z1×A3×z2T), similar to (26), but the 5×5 matrices B3, A3 now have the form:

B3=aαPφQφ          A3=aαPφQφ+j(a2PφPφω02QφQφ)E37

Figure 9.

(a) Contour plot of a two-directional filter with θ=π/6, φ=π/5; frequency response magnitudes and contour plots of two corrected two-directional filters with parameters: θ=π/6, φ=π/5 (b), (c) and θ=π/2, φ=π/4 (d), (e)

The denominator matrix A3 has complex elements. In Fig.9 (a), the contour plot of the frequency response magnitude is shown for a two-directional filter with aperture θ=π/6 and orientation φ=π/5. As with the previous types of filters, the marginal linearity distortions can be corrected using an additional LP square-shaped filter.

The templates B3C and A3C of the corrected filter result by convolution: B3C=B3BC and A3C=A3AC. In Fig.9 (b), (c) the frequency response magnitudes and contour plots are displayed for the corrected two-directional filter with specified aperture and orientation. The initial distortions have been eliminated.

The second two-directional filter in Fig.9 (d), (e) is a particular case, being oriented along the two frequency axes (θ=π/2, φ=π/4), therefore can be used to detect simultaneously horizontal and vertical lines from an image, as shown in the next section.

5.2. Three-directional fan-type filters

In order to design a three-directional filter like the one depicted in Fig. 8 (b), we must start from an analog three-band selective filter, like the one with frequency response shown in Fig. 8 (c). For a three directional filter, the middle peak frequency can always be taken ω0=0, and the other two on each side at specified values. The prototype transfer function HP(s) in variable s will be in this case the sum of three elementary functions:

HP(s)=BP(s)AP(s)=αs+α+αs+α+jω01+αs+α+jω02E38

The frequency response of a filter of this kind with parameter values α=0.03, ω01=0.59π and ω02=0.36π is shown in Fig. 8 (c). Substituting the mapping (34) into the expression (8) of the elementary function HjS(s), we get the 2D transfer function H1(z1,z2) in matrix form: H1(z1,z2)=(z1×Bb×z2T)/(z1×Ab×z2T), where z1=[1   z1   z12], z2=[1   z2   z22] and the 3×3 templates Bb,Ab are given by:

Bb=αQφ;Ab=αQφ+j(aPφ+ω0Qφ)E39

Each of the three elementary terms in (38) corresponds to a pair of 3×3 templates Bb and Ab given by (39). If the three elementary filters are given by the pairs of templates (Bb1, Ab1), (Bb2, Ab2) and (Bb3, Ab3), the templates of size 7×7 of the entire three-directional filter will result by summing up the convolutions of elementary templates:

B3=Bb1Ab2Ab3+Ab1Bb2Ab3+Ab1Ab2Bb3 ;     A3=Ab1Ab2Ab3E40

The numerator BP(s) of HP(s) from (38) has the general form:

BP(s)=α(a2s2+a1s+a0)E41

where a2=3, a1=6α+j2(ω01+ω02) and a0=3α2ω01ω02+j2α(ω01+ω02).

We see that a2 is real and a0, a1 are generally complex. The coefficients a0, a1 are real only when ω02=ω01, i.e. for symmetric frequency values around the origin. Finally, for any specified set of values α, ω01, ω02 the denominator factorizes as Bj(s)=3α(s+r1)(s+r2), where r1 and r2 are complex roots. Therefore the factorized prototype transfer function is:

Hj(s)=Bj(s)Aj(s)=3α(s+r1)(s+r2)(s+α)(s+p1)(s+p2)E42

Figure 10.

Frequency response magnitude (a) and contour plot (b) of a corrected three-directional filter with parameters: θ=0.23π, φ=0.27π, ω01=0.59π, ω02=0.36π

At the denominator, we denoted p1=α+jω01, p2=α+jω02. Applying to each factor the frequency transformation sFφ(z1,z2) given by (34), after some algebraic manipulations, we finally obtain the templates of the three-directional filter as discrete convolutions:

B3=3αQφ(r1Qφ+jaPφ)(r2Qφ+jaPφ)E43
A3=(αQφ+jaPφ)(p1Qφ+jaPφ)(p2Qφ+jaPφ)E44

This implies the fact that the transfer function H3(z1,z2) of the 2D three-directional filter with templates B3 and A3 of size 7×7 results directly in a factorized form, which is an important advantage in implementation. As a general remark on the method, using an analog prototype instead of a digital one, as is currently done, simplifies the design in this case, as the frequency mapping results simpler and leads to a 2D filter of lower complexity. The designed filters result with complex coefficients, however such IIR filters can also be implemented (Nikolova et al., 2011).

Advertisement

6. Directional IIR filters designed in polar coordinates

We approach here a particular class of 2D filters, namely filters whose frequency response is symmetric about the origin and has at the same time an angular periodicity. The contour plots of their frequency response, resulted as sections with planes parallel with the frequency plane, can be defined as closed curves which can be described in terms of a variable radius which is a periodic function of the current angle formed with one of the axes.

It can be described in polar coordinates by ρ=ρ(φ), where φ is the angle formed by the radius op with ω1-axis, as shown in Fig.8(a) for a four-lobe filter. Therefore ρ(φ) is a periodic function of the angle φ in the range φ[0,2π].

6.1. Spectral transformation for filters designed in polar coordinates

The main issue approached here is to find the transfer function of the desired 2D filter H2D(z1,z2) using appropriate frequency transformations of the form ωF(ω1,ω2). The elementary transfer functions (14) and (15) have the complex frequency responses:

H1(jω)=(b0+b1cosω+jb1sinω)/(a0+cosω+jsinω)E45
H2(jω)=b1+(b2+b0)cosω+j(b2b0)sinωa1+(1+a0)cosω+j(1a0)sinω=P(ω)Q(ω)E46

The proposed design method for these 2D filters is based on the frequency transformation:

F:2,ω2F(z1,z2)=Bf(z1,z2)/Af(z1,z2)E47

which maps the real frequency axis ω onto the complex plane (z1,z2), defined by the real frequency mapping:

F1:2,ω2F(ω1,ω2)=(ω12+ω22)/ρ(ω1,ω2)E48

In (48) ρ(ω1,ω2) is initially determined in the angle variable φ as ρ(φ) and can be referred to as a radial compressing function. In the frequency plane (ω1,ω2) we have:

cosφ=ω1/ω12+ω22E49

Figure 11.

(a) contour plot of a four-lobe filter; (b) periodic function ρ(φ); (c) LP prototype

If the radial function ρ(φ) can be expressed in the variable cosφ, using (49) we obtain by substitution the function ρ(ω1,ω2). The function ρ(φ) will result as a polynomial or a ratio of polynomials in cosφ. For instance, the four-lobe filter with the contour plot given in Fig.11 (a) corresponds to a function:

ρ(φ)=a+bcos4φ=a+b8bcos2φ+8bcos4φE50

plotted in Fig.11 (b) on the range φ[0,2π]. More generally, the 2D filter can be rotated in the frequency plane with a specified angle φ0 about one of the frequency axes, e.g. Oω2. For instance, in a four-lobe filter, two opposite lobes are oriented along a direction at an angle φ0, and the other two at φ0+π/2, as shown in Fig.12 (b). It can be shown that the cosine of the current angle φ with initial phase φ0 can be expressed as:

cos2(φ+φ0)=(cos2φ0ω12+sin2φ0ω22+0.5sin2φ0ω1ω2)/(ω12+ω22)E51

For filters with an even number of lobes, the radial function ρ(φ) is expressed in even powers of cosφ or cos(φ+φ0). The frequency transformation (48) can be also expressed as:

ω(ω12+ω12)/ρ(ω1,ω2)=F1(ω1,ω2)E52

In order to obtain a rational expression for the frequency response of the 2D filter from an elementary 1D prototype of the form (45) or (46) by applying the frequency mapping (52), we need to derive rational expressions for the functions cosω and sinω. Using the Chebyshev-Padé method in a symbolic computation software, the following second-order rational approximations were found:

cosω(1.05590.086514ω0.1304ω2)/(1+0.75ω0.110583ω2)=CS(ω)/AS(ω)E53
sinω(0.167+1.46287ω0.259815ω2)/(1+0.75ω0.110583ω2)=SS(ω)/AS(ω)E54

which are sufficiently accurate on the range ω[0,π]. Since these functions are developed on the range [0,π], their approximations result neither odd nor even. However, using the above approximations will lead to relatively complex 2D filters, described by templates of size at least 9×9. For the type of filters approached, namely selective two-directional (four-lobe) filters, the approximations for cosω and sinω need not necessarily hold throughout the range [0,π], but only on a smaller range near the origin, corresponding to filter pass-band. Using now the Padé method we get the first-order approximations:

sinω(s0+s1ω)/(1+rω)               cosω(c0+c1ω)/(1+rω)E55

with s0=0.0928, s1=2.5218, c0=1.0104, c1=1.2193, r=1.979, which hold only on a

narrower range around zero of the interval ω[0,π]. Using (55) instead of (53), (54) will result in much more efficient 2D filters, which fully satisfy the imposed specifications.

We will use here a Chebyshev low-pass second-order filter of the general form (15). For this type of filter we have the coefficient symmetry b2=b0. According to (46) we can write:

H2(jω)=b1+2b0cosωa1+(1+a0)cosω+j(1a0)sinω=P(ω)Q(ω)E56

The numerator results real because the imaginary part is cancelled. Substituting the expressions (55) into this complex frequency response we get the rational approximation:

H2(jω)=b1(1+rω)+2b0(c0+c1ω)a1(1+rω)+(1+a0)(c0+c1ω)+j(1a0)(s0+s1ω)E57

which can be also written as:

H2(jω)=b1(1+rω2)+2b0(c0+c1ω2)a1(1+rω2)+(1+a0)(c0+c1ω2)+j(1a0)(s0+s1ω2)E58

The function (58) has even parity, since it is expressed as a rational function in ω2.

6.2. Two-directional filter design

We approach now the design of a particular filter type designed in polar coordinates, namely two-directional (selective four-lobe) filters along the two plane axes or with a specified orientation angle. Let us consider the radial function given by:

Hr(φ)=1/(pB˜(φ)p+1)E59

where B˜(φ) is a periodic function; let B˜(φ)=cos(4φ). We use it to design a 2D filter with four narrow lobes in the plane (ω1,ω2). Using trigonometric identities, (59) becomes:

Hr(φ)=1/(1+8p(cosφ)28p(cosφ)4)E60

and is plotted for φ[π,π] in Fig.12(a). This is a periodic function with period Φ=π/4 and has the shape of a multi-band (“comb”) filter. In order to control the amplitude of this function, we introduce another parameter k, such that the radial function ρ(φ) takes the form ρ(φ)=kHr(φ). We get using (49):

ω2F(ω1,ω2)=(ω14+(2+8p)ω12ω22+ω24)/(k(ω12+ω22))E61

and the function F2(s1,s2) of the form:

F2(s1,s2)=(s14+(2+8p)s12s22+s24)/(k(s12+s22))E62

Finally we derive a transfer function of the 2D filter H(z1,z2) in the complex plane (z1,z2). This can be achieved if we find a discrete counterpart of the function ρ(ω1,ω2), denoted R(z1,z2). A possible method is to express the function ρ(ω1,ω2) in the complex plane (s1,s2) and then find the appropriate mapping to (z1,z2) using the bilinear transform or the Euler approximation. Even if generally the bilinear transform is more often used, being a more accurate mapping, especially if a frequency pre-warping is applied to compensate for distortions, for the particular type of 2D filter approached here the Euler formula leads to more efficient filters and with better characteristics, at the same filter order. We will use the backward Euler method, which approximates the spatial derivative X/x by X[n]X[n1], replacing s by s=1z1. On the two directions of the plane we have: s1=1z11, s2=1z21. The operators s12, s22 and s1s2 correspond to second-order partial derivatives: 2/x2s12=z1+z112, 2/y2s22=z2+z212, 2/xys1s2. For the mixed operator s1s2, using repeatedly the Euler formula, we get the expression (Matei, 2011 a): 2s1s2=z1+z11+z2+z212z1z21z11z2. Substituting the above relations into (62) we obtain a frequency mapping similar to (47), with the templates:

Bf=[0010002+8p816p2+8p01816p20+32p816p102+8p816p2+8p000100]Af=k[010141010][010141010]=kA1A1E63

The template Af results as a convolution of two 3×3 matrices. The last step in the design of this 2D filter is to apply the frequency transformation (61) to the frequency response (58) and we find the filter templates B and A as linear combinations of Bf and Af:

B=(b1+2b0c0)Af+(rb1+2b0c1)BfE64
A=(a1+(1+a0)c0+j(1a0)s0)Af+(a1r+(1+a0)c1+j(1a0)s1)BfE65

where b2, b1, b0, a1, a0 are the coefficients of prototype (15). Finally the 2D filter transfer function in z1 and z2 has the following expression, with z1 and z2 given by (27):

H2D(z1,z2)=B(z1,z2)/A(z1,z2)=(z1×B×z2T)/(z1×A×z2T)E66

Let us design a two-directional filter following this procedure. As 1D prototype let us consider a type-2 low-pass Chebyshev digital filter with the parameter values: orderN=2, stop-band attenuation Rs=40db and passband-edge frequency ωp=0.5 (1.0 is half the sampling frequency). The transfer function in z is:

Hp(z)=(0.012277z20.012525z+0.012277)/(z21.850147z+0.862316)E67

For a good directional selectivity we also choose p=30 and k=10. The 2D filter frequency response magnitude is displayed in Fig.12 (c) and shows a very good linearity along the two directions and practically no distortions in the stop band. The constant level contour in the plane (ω1,ω2) is given in Fig.12 (d). Calculating the singular values of filter templates for the above parameters we find the vectors SA, SB for the templates A, B:

SA=[2.093470.002250.000500]      SB=[2.142970.019390.0038700]E68

Taking into account the fact that the first singular value of the templates A and B is much larger than the other four, the filter designed above can be approximated by a separable filter.

The singular value decomposition of a matrix M is written as M=U×S×V where U and V are unitary matrices and S is a diagonal matrix containing the singular values. Thus we can

write for the filter templates A and B:

A=UA×SA×VAB=UB×SB×VBE69

If UA1 and VΑ1 are the first columns of the matrices UA, VA, then A can be approximated by a matrix A1=sA1UA1VA1T, where sA1 is the largest singular value of A, UA1 and VΑ1 are the corresponding columns of UA and VA, stands for outer product, T for transposition. Similarly for B we find B1=sB1UB1VB1T. For the specified filter parameters we obtain sB1=2.14297 and for template B the column vectors UB1 and VB1 result identical: UB1=VB1=[-0.004240.3971-0.827430.3971-0.00424]T.

For the template A we get sA1=2.09347 and the vectors UA1,VΑ1 have complex elements. The frequency response of the resulted filter is given in Fig.12 (e). As can be noticed, the effect of the above approximation is an “overshoot” at zero frequency. This should not affect the filter functionality in detecting lines parallel with the two axes. Moreover, since the marginal elements of the 5×1 vectors UA1, VA1, UB1, VB1 have negligible values, by discarding them we obtain the 3×1 vectors: UB2=VB2=[0.3971-0.82740.3971]T

UA2=[-0.3150.6334-0.315]T+j[0.2573-0.51750.2573]T

VA2=[0.4067-0.8180.4067]Tj[0.0024-0.00480.0024]T

We finally obtain a very selective two-directional 2D filter implemented with two minimum size (3×3) templates. The template B is real while A is complex. The frequency response magnitude of this filter is shown in Fig.12 (f) and is practically similar to the one in Fig.12 (e). Similarly we can design a two-directional (four-lobe) filter with a specified orientation angle. Using the previously described method and based on the Euler approximation, the expression (10) of cos2(φ+φ0) corresponds to a frequency transformation in the complex variables z1 and z2, written in matrix form as:

Figure 12.

(a) Periodic radial function; (b) contour plot of a two-directional filter with orientation φ0=π/12; (c) frequency response and (d) contour plot of the two-directional filter with 5×5 templates; frequency response of the filter with separable 5×5 (e) and 3×3 templates (f)

cos2(φ+φ0)F(z1,z2)=Bφ0(z1,z2)/Aφ0(z1,z2)=(Z1×Bφ0×Z2T)/(Z1×Aφ0×Z2T)E70
Bφ0=cos2φ0[010020010]+sin2φ0[000121000]+0.25sin(2φ0)[011111110]E71

and Aφ0 is identical to A1 from (63). The 5×5 templates of the mapping (47) are given by:

Bf=Aφ0Aφ0Af=Aφ0Aφ0+8pBφ0Aφ08pBφ0Bφ0E72

The final filter templates result according to relations (64) and (65).

Regarding the proposed method, the frequency responses of this class of 2D filters can be viewed as derived through a radial distortion from a generic maximally-flat circular filter. Indeed, referring to (48), the circular filter is the trivial case for which ρ(ω1,ω2)=1 and the mapping (48) reduces to ω2ω12+ω12, as expected. This method allows one to design 2D filters with a non-convex shape in the frequency plane. The proposed design method does not involve global numerical optimization techniques, but only a few numerical approximations. The method is more general and can be applied as well to design fan filters, diamond filters, multi-directional filters etc. (Matei, 2011 a).

Advertisement

7. Zero-phase FIR circular filters

Filters with circular symmetry are very useful in image processing. We propose an efficient design technique for 2D circularly-symmetric filters, based on the previous 1D filters, considered as prototypes. Given a 1D prototype HP(ω), the corresponding 2D circular filter function HC(ω1,ω2) results using the frequency mapping ωω12+ω22:

HC(ω1,ω2)=HP(ω12+ω22)E73

The currently-used approximation of the circular function cosω12+ω22 is given by:

cosω12+ω22C(ω1,ω2)=0.5+0.5(cosω1+cosω2)+0.5cosω1cosω2E74

which corresponds to the 3×3 array:

C=[0.1250.250.1250.250.50.250.1250.250.125]E75

Let us consider as prototype a LP analog elliptic filter of order N=4, pass-band peak-to-peak ripple RP=0.04 dB, stop-band attenuation RS=40 dB and passband-edge frequency ΩP=π/2. Its transfer function in variable s is:

HP(s)=0.1037(s4+19.864s2+84.041)/(s4+3.2041s3+8.4315s2+13.126s+14.082)E76

Using MAPLE or another symbolic computation program and following the design steps described in section 2, we obtain a a polynomial approximation of the magnitude |HP(jω)| through Chebyshev expansion, which has the following factorized form, with x=cosω:

|HP(ω)|48.6(x+0.8491)(x+0.7717)(x1.087)(x2+1.9934x+0.994)(x2+1.0797x+0.318)(x20.3849x+0.1766)(x21.2882x+0.5314)(x21.9338x+0.9726)E77

In order to obtain a filter with circular symmetry from the factorized 1D prototype function, we replace in (12) cosω with the circular cosine function (74). For instance, corresponding to (12), the filter template A results in general as the discrete convolution:

A=kA11A12A1nA21A22A2mE78

where A1i (i=1n) are 3×3 templates and A2j (j=1m) are 5×5 templates, given by: A1i=C+aiA01 and A2j=CC+a1jC0+a2jA02, where A01 is a 3×3 zero template and A02 a 5×5 zero template with the central element equal to one; C0 is a 5×5 template

obtained by bordering C with zeros. The above expressions correspond to the factors in (12).

The frequency response HC(ω1,ω2) of the 2D circular filter results in a factorized form by substituting x=C(ω1,ω2) in (77). Even if the filter results of high order, with very large templates, next we show that using the Singular Value Decomposition (SVD), the resulted 2D filter can be approximated with a negligible error. For the filter template B we can write B=UB×SB×VB. The vector of singular values SB of size 1×27 has 14 non-zero elements:

SB1=[0.50536   0.086111   0.032794   0.013627   0.00521   0.002937   0.001935              0.001061  0.000639   0.000451   0.000418   0.0000385   0.0000196   0.00000144]

Figure 13.

Frequency response magnitude (a) and contour plot (b) of a circular FIR filter

Figure 14.

Frequency response magnitudes and contour plots for the circular filter resulted by taking into account the first largest: (a), (b) 8 singular values; (c), (d) 5 singular values

Let us denote the vector above as SB1=[sk], with k=1...14 in our case. The exact filter matrix B can be written as: B=UB1×SB1×VB1, where UB1 and VB1 are made up of the first 14 columns of the unitary matrices UB and VB. If we consider the first largest M values of the vector SB1=[sk], the matrix B can be approximated as:

BBM=k=1MskUBkVBkTE79

Here BM is the approximation of matrix B taking into account the first M singular values (in our case M14), while UBk, VBk are the k-th columns of the matrices UB1 and VB1; stands for outer product and the superscript T for transposition.

Fig.14 shows the frequency response magnitudes of the designed circular filter approximated by taking into account the first largest 8 singular values and 5 singular values. It can be noticed that even retaining only the first 5 singular values, the 2D filter preserves its circular shape without large distortions. In this case the filter template B is approximated by BM from (79), for M=5. Therefore, the template B can be written as a sum of only 5 separable matrices according to (79). This is an important aspect in the filter implementation.

Advertisement

8. Applications and simulation results

An example of image filtering with a two-directional filter is given. We use the filter shown in Fig.3(e), (f). This type of filter can be used in simultaneously detecting perpendicular lines from an image. The binary test image in Fig.15 (a) contains straight lines with different orientations and lengths, and a few curves. It is known that the spectrum of a straight line is oriented in the plane (ω1,ω2) at an angle of π/2 with respect to the line direction. Depending on filter selectivity, only the lines with the spectrum oriented more or less along the filter pass-bands will remain in the filtered image. In the output image in Fig.15 (b), the lines roughly oriented horizontally and vertically are preserved, while the others are filtered out or appear very blurred, due to directional low-pass filtering. The joints of detected lines appear as darker pixels and can be detected, if after filtering a proper threshold is applied.

Let us apply the designed fan-type filters, which can be regarded as components of a DFB, in filtering a typical retinal vascular image. Clinicians usually search in angiograms relevant features like number and position of vessels (arteries, capillaries). An angular-oriented filter bank may be used in analyzing angiography images by detecting vessels with a given orientation. Let us consider the retinal fluorescein angiogram from Fig.16(a), featuring some pathological elements which indicate a diabetic retinopathy. This image is filtered using 5 oriented wedge filters with narrow aperture (θ=π/24), designed using the method described in section 4. Fig.16 (b)-(f) show the directionally filtered angiography images.

Figure 15.

(a) Test image; (b) filtered image

Figure 16.

(a) Retinal fluorescein angiogram; (b)-(f) images resulted as output of five component filters of the fan-type filter bank

Figure 17.

(a) Retinal fluorescein angiogram; (b), (c) results of image filtering with a FIR LP circular filter with cutoff frequency ΩC1=0.16π and ΩC2=0.08π

As can be easily noticed, the vessels for which the frequency spectrum overlaps more or less with the filter characteristic remain visible, while the others are blurred, an effect of the directional low-pass filtering (Matei & Matei, 2012). The directional resolution depends on the filter angular selectivity given by θ.

The designed zero-phase circularly-symmetric FIR filters may be useful as well in pre-processing tasks on biomedical images, having a blurring effect on the image which depends on its selectivity given by the circular filter bandwidth. The effect is somewhat similar to the Gaussian smoothing, which is used as a pre-processing stage in computer vision tasks to enhance image structures at different scales. Applying the presented design procedure, a circularly-symmetric filter bank can be derived, with components having desired bandwidths. Let us consider another retinal fluorescein angiogram, displayed in Fig.17(a). In the simulation result shown in Fig.17 (b) and (c), the two circular filters introduce gradual blurring which is visible on the fine image details, like small vessels and capillaries. In the image in Fig.17 (c) all the finer details have been almost completely smoothed out.

Advertisement

9. Conclusion

The design methods presented in this chapter combine the analytical approach based on 1D prototype filters and frequency transformations with numerical optimization techniques. For the classes of 2D filters designed here, we have used mainly analog filters as prototypes, which turn out to make simpler the expressions of the derived frequency mappings, and therefore the complexity of the designed 2D filters is lower in the analyzed cases. The prototypes used here were both maximally-flat or very selective, either low-pass or band-pass. For each type of 2D filter, a particular spectral transformation is derived. An important advantage is that these spectral transformations include some parameters which depend on the 2D filter specifications, like bandwidth, orientation, aperture etc. Once found the specific frequency mapping, the 2D filter results from its factorized prototype function by a simple substitution in each factor. The designed filters are versatile in the sense that the prototype parameters (bandwidth, selectivity) can be adjusted and the 2D filter will inherit these properties.

An advantage of the analytical approach over the completely numerical optimization techniques is the possibility to control the 2D filter parameters by adjusting the prototype. Another novelty is the proposed analytical design method in polar coordinates, which can yield selective two-directional and even multi-directional filters, and also fan and diamond filters. In polar coordinates more general filters with a specified rotation angle can be synthesized.

The design methods approached here are rather simple, efficient and flexible, since by starting from different specifications, the matrices of a new 2D filter result directly by applying the determined frequency mapping, and there is no need to resume every time the whole design procedure.

Stability of the designed filters is also an important problem and will be studied in detail in future work on this topic. In principle the spectral transformations used preserve the stability of the 1D prototype. The derived 2D filter could become unstable only if the numerical approximations introduce large errors. In this case the precision of approximation has to be increased by considering higher order terms, which would increase in turn the filter complexity; however, this is the price paid for obtaining efficient and stable 2D filters. Further research will focus on an efficient implementation of the designed filters and also on their applications in real-life image processing.

References

  1. 1. AnsariR1987Efficient IIR and FIR fan filters, IEEE Transactions on Circuits and Systems, 34Aug. 1987, 941945
  2. 2. AustvollI2000Directional filters and a new structure for estimation of optical flow. Proc. of Int. Conf. on Image Processing ICIP 2000, Vancouver, Canada, 2574577
  3. 3. BambergerRSmithM1992A filter bank for the directional decomposition of images: theory and design, IEEE Trans. Signal Processing, 40Apr. 1992, 882893
  4. 4. BerryE2007A Practical Approach to Medical Image Processing. Taylor & Francis, 2007
  5. 5. ChakrabartiSMitraS. K1977Design of two-dimensional digital filters via spectral transformations, Proceedings of the IEEE, 65June 1977, 905914
  6. 6. DanielssonP. E1980Rotation-invariant linear operators with directional response. Proc. of 5th International Conf. on Pattern Recognition, Miami, USA, Dec. 1980
  7. 7. DoughertyGeditor). (2011Medical Image Processing: Techniques and Applications. Springer, 2011
  8. 8. FrangiA. Fet al1998Multiscale vessel enhancement filtering, Intl. Conf. on Medical Image Computing Computer-Assisted Intervention, 1496Berlin, 1998, 130137
  9. 9. FreemanW. TAdelsonE. H1991The design and use of steerable filters. IEEE Trans. on Pattern Analysis and Machine Intelligence, 13Sept. 1991, 891906
  10. 10. HarnLShenoiB1986Design of stable two-dimensional IIR filters using digital spectral transformations. IEEE Transactions on Circuits and Systems, 33May 1986, 483490
  11. 11. HiranoKAggarwalJ. K1978Design of two-dimensional recursive digital filters. IEEE Trans. on Circuits and Systems, 25Dec. 1978, 10661076
  12. 12. ItoN2010Efficient design of two-dimensional diamond-shaped filters. Proceedings of Int. Symposium ISPACS 2010, Chengdu, China, 6-8 Dec. 2010, 14
  13. 13. JuryE. IKolavennuV. RAndersonB. D1977Stabilization of certain two-dimensional recursive digital filters. Proceedings of the IEEE, 656887892
  14. 14. KayranAKingR1983Design of recursive and nonrecursive fan filters with complex transformations, IEEE Trans. on Circuits and Systems, CAS-30(12), 1983, 849857
  15. 15. KhanM. A. Uet al2004Coronary angiogram image enhancement using decimation-free directional filter banks, IEEE ICASSP, Montreal, May 17-21, 2004, 5441444
  16. 16. LimJ. S1990Two-Dimensional Signal and Image Processing. Prentice-Hall 1990
  17. 17. LimY. CLowS. H1997The synthesis of sharp diamond-shaped filters using the frequency response masking approach. Proc. of IEEE Int. Conf. on Acoustics, Speech & Signal Processing ICASSP-97, 21812184Munich, Germany, Apr. 21-24, 1997
  18. 18. LowS. HLimY. C1998A new approach to design sharp diamond-shaped filters. Signal Processing, 67May 1998, 3548
  19. 19. LuW. SAntoniouA1992Two-Dimensional Digital Filters, CRC Press, 1992
  20. 20. MastorakisN. E2000New necessary stability conditions for 2D systems, IEEE Trans. on Circuits and Systems, Part I, 4711031105July 2000
  21. 21. MateiR2010A new design method for IIR diamond-shaped filters, Proc. of the 18th European Signal Processing Conference EUSIPCO 2010, Aalborg, Denmark, 6569
  22. 22. MateiR2011aNew Design Methods for Two-Dimensional Filters Based on 1D Prototypes and Spectral Transformations, In: "Digital Filters", Fausto Pedro García Márquez (Ed.), 91121IN-TECH Open Access Publisher, Vienna, 2011
  23. 23. MateiR2011bA class of 2D recursive filters with two-directional selectivity, Proc. of the WSEAS International Conference on Applied, Numerical and Computational Mathematics (ICANCM’11), 212217Barcelona, Spain, Sept. 15-17, 2011
  24. 24. MateiRMateiD2012Vascular image processing using recursive directional filters, World Congress on Medical Physics and Biomedical Engineering, Beijing, China, May 26-31, 2012, IFMBE Proceedings 39947950
  25. 25. MollovaG. S1997Analytical least squares design of 2-D fan type FIR filter, International Conference on Digital Signal Processing DSP’97, 2625628July 1997
  26. 26. NieXUnbehauenR1989D IIR filter design using the extended McClellan transformation, Proc. of International Conference ICASSP’89, 23-26 May 1989, 315721574
  27. 27. NikolovaZet al2011Complex coefficient IIR digital filters, In: “Digital Filters“, InTech, April 2011, 209239
  28. 28. OConnorB. THuangT.S. (1978Stability of general two-dimensional recursive digital filters, IEEE Trans. Acoustics, Speech & Signal Processing, 26550560
  29. 29. PaplinskiA. P1998Directional filtering in edge detection. IEEE Transactions on Image Processing, Apr. 1998, 7611615
  30. 30. PsarakisE. Zet al1990Design of two-dimensional zero phase FIR fan filters via the McClellan transform. IEEE Trans. Circuits & Systems, Jan. 1990, 371016
  31. 31. QunshanGSwamyM. N. S1994On the design of a broad class of 2D recursive digital filters with fan, diamond and elliptically-symmetric responses, IEEE Trans. on Circuits and Systems II, 41Sep.1994, 603614
  32. 32. RydellJet al2008Bilateral filtering of fMRI data. IEEE Journal of Selected Topics in Signal Processing, Dec 2008, 2891896
  33. 33. SemmlowJ. L2004Biosignal and Biomedical Image Processing. Marcel Dekker, 2004
  34. 34. SimoncelliE. PFaridH1996Steerable wedge filters for local orientation analysis. IEEE Transactions on Image Processing, 5Sep 1996, 13771382
  35. 35. TosicD. Vet al1997Symbolic approach to 2D biorthogonal diamond-shaped filter design, 21st Int. Conf. on Microelectronics, 1997, Nis, Yugoslavia, 2709712
  36. 36. TrucPet al2007A new approach to vessel enhancement in angiography images, Int. Conf. on Complex Medical Engineering CME 2007, Beijing, 23-27 May 2007, 878884
  37. 37. WongW. C. Ket al2004Trilateral filtering for biomedical images. IEEE Int. Symposium on Biomedical Imaging, 15-18 Apr. 2004, 1820823
  38. 38. WuDet al2006On the adaptive detection of blood vessels in retinal images, IEEE Trans. Biomedical Engineering, 53Feb 2006, 341343
  39. 39. ZhuWPZhenyaH. (1990A design method for complementary recursive fan filters, IEEE ISCAS 1990, 1-3 May 1990, 321532156
  40. 40. ZhuWPNakamuraS. (1996An efficient approach for the synthesis of 2-D recursive fan filters using 1-D prototypes, IEEE Transactions on Signal Processing, 44979983Apr. 1996
  41. 41. ZhuWP. et al1999A least-square design approach for 2D FIR filters with arbitrary frequency response. IEEE Transactions on Circuits and Systems II, 4610271034Aug. 1999
  42. 42. ZhuWP. et al2006Realization of 2D FIR filters using generalized polyphase structure combined with singular-value decomposition, Proc. of ISCAS 2006, Kos, Greece

Written By

Radu Matei and Daniela Matei

Submitted: 06 April 2012 Published: 16 January 2013