Boundary Treatment for Young-van Vliet Recursive Zero-Mean Gabor Filtering

. This paper deals with convolution setting at boundary regions for 1D convolution computed during recursive Gaussian and Gabor ﬁltering as well as staged Gabor ﬁltering computed more e ﬃ ciently as modulation, recursive Gaussian, and demodulation. These are established fast approximations to their ﬁlters. Until recently, all the three applications of recursive ﬁlters su ﬀ ered from distortions near boundary as a result of inappropriate boundary treatment. The extension of input data with constant border value is presumed. We review a recently suggested setting for recursive Gaussian and Gabor ﬁltering. Then, a new convolution setting for the more e ﬃ cient staged Gabor ﬁltering is presented. We also o ﬀ er a formula to compute the scale coe ﬃ cient, using which a zero-mean Gabor ﬁlter can be obtained from either recursive or staged Gabor ﬁlter.


Introduction
Shortly after publication, the paper by Young et al. [1] (we will use the acronym YvV based on their paper [2] where the framework was first introduced) has become a standard reference for the implementation of Gaussian filtering in image processing although the paper was primarily focused on Gabor filtering.The shift towards Gaussian filtering may be in part because the authors originally developed filtering coefficients for Gaussian from which coefficients for Gabor filter were subsequently derived and in part because the Gabor filtering can be implemented with advantage as modulation, Gaussian filtering, and demodulation [3,4]we refer to this as to the staged Gabor filtering.Thus, the Gaussian filtering can become central even to Gabor filtering.
The original YvV's paper utilizes bidirectional recursive infinite-impulse-response (IIR) filtering where a "forward" pass processes input data in a causal direction and a "backward" pass processes results of the forward pass in an anticausal direction The i t represents input value at position t from 1D data of length n, for example, a scanning line of a digital image.Similarly, u t and v t hold results of forward and backward passes at position t, respectively.The a j and b j are the forward and backward pass coefficients, respectively.They are set based on Gaussian sigma parameter [1].The resulting v t is an overscaled filtering result because we, intentionally, adopted both equations and labeling from Triggs and Sdika [5].To comply with YvV's notation introduced in [1], for Gaussian filtering set  10) of [1].To get properly scaled result, multiply v t afterwards by B from (18) of [1].As a matter of fact, the forward pass, (1), and the backward pass, (2), for Gaussian are the same except for the direction.
Even though the above describes only 1D convolution, there exist approaches to 2D, 3D, or nD Gaussian and Gabor filtering techniques utilizing cascades of several such 1D convolutions and the separability of the filters; for example, isotropic Gaussian can be computed with convolutions computed along the main coordinate-system axes.If 1D convolution along general arbitrarily oriented real axis is available, a general anisotropic 2D [6] or nD [7] Gaussian filtering can be computed.A slightly more stable and a bit more computationally demanding arbitrary 2D [8] or nD [9] Gaussian filtering utilizes 1D convolutions only along integer axes (axes where orientation vector consists only of integer components).As a consequence, any arbitrarily oriented anisotropic Gabor filtering can be carried out then by using the staged filtering technique but, as we shall see in Section 4, a special care must be taken at boundaries.
It is very important to handle boundary initialization correctly in the YvV's recursive filter because the filter tends to reflect relatively large portion of its response history, see Figures 1 and 2. The span is estimated up to 3σ [5].Hence, a care must be taken of what values of u −2 , u −1 , u 0 to initiate with the forward pass and what values of v n+1 , v n+2 , v n+3 to initiate with the backward pass (the forwardbackward transition) in order to avoid cumulation of errors in the border regions of image due to the cascading of 1D convolutions.Depending on one's preference, the input data may be prefixed with a constant sequence of i t = i 1 , t = −∞, . . ., 0 and, similarly, postfixed with i t = i n , t = n + 1, . . ., ∞ to facilitate filtering.A padding with zeros or mirroring input data near border are also plausible choices.
In the two figures, we illustrate how different techniques for a Gaussian and Gabor filtering perform near boundary on some randomly generated data.The direct filtering, both Gaussian and Gabor, with YvV boundary initialization assumes that backward pass is initiated with a steady-state response of the filter.But recursive filters, typically, decay to the steady state on constant input only after several iterations.The constant sequence would have to begin already a way ahead the image boundary, which cannot be always guaranteed, so that the forward pass, if continued after the image boundary, would enter the boundary region already with constant steady-state responses.Then, the YvV forward-backward initialization would perform well.In the staged techniques, the use of Gaussian filtering with any of the established boundary treatment may be tempting, but it is incorrect as demonstrated in Figure 2. Again, the assumption of constant boundaries in these cases of staged Gabor techniques is violated because of the reasons to be discussed in detail in Section 4. Note that all of the currently available boundary treatment techniques for YvV recursive filtering are presented in both figures, except for those we are going to discuss in this paper.
In this paper, we deal with approaches to correct initialization of YvV filters in applications of direct Gaussian in Section 2, direct Gabor in Section 3, and staged Gabor filtering in Section 4. The first and the second of the three cases were dealt with in [5] explicitly and implicitly, respectively.We had only replayed their derivation for the second case, a complex direct Gabor filtering, to arrive to some elegant solution presented in this paper.The last case forms the main contribution of this paper.We always consider extending the input 1D data with i 1 and i n in the anti-causal and causal directions, respectively, that is, extending input data line with constant value of the pixel at line's end.Finally, in Section 5, we offer a formula to compute the scale coefficient, using which a zero-mean Gabor filter can be obtained from either direct or staged Gabor filter.This paper is supposed to summarize all necessary information for efficient and correct 1D recursive implementation of any of the three filtering applications.Such an implementation would only require reading the paper by Young et al. [1], in addition to this one, to learn how to find coefficients for the forward and backward passes.

Direct Gaussian Filtering
In this section, we will remind of the results of the recent paper by Triggs and Sdika [5] who derived closed-form equations for computing correct forward and backward initialization values based solely on filter parameters.Indeed, their solution does not require any prior knowledge of the filter, for example, Gaussian sigma, dimensionality, isotropy, the way filter was derived, and so forth.It only requires that the filtering should be realized by means of forward and consecutive backward passes, (1) and (2), even with arbitrary depth of recursion.Based on their assumption that linear filtering can be written using matrix multiplication, they managed to first describe the filtering itself and then, the most importantly, also the forward-backward transition.
Triggs and Sdika [5] confirmed that for 3rd order recursive filters according to (1) and (2), the correct initialization of the forward pass with the steady-state response to an infinite stream of i 1 is For the staged cases, which is always a modulation, Gaussian filtering, and demodulation, we show that any forward-backward transition for direct Gaussian filtering (discussed in Section 2) shall not be used, despite that this transition may be correct for regular Gaussian filtering.The results on the forward-backward transition are compared to the FIR implementation with constantly extended border as well.The FIR serves as a ground-truth result.The error spans 7 data positions.Only real parts of the complex filtering results are shown here.
However, their main contribution was the derivation of an explicit formula for the forward-backward transition.We reproduce their result, contained in their ( 14), (15), and (5) from [5], for k = l = 3 and i where M is a 3 × 3 matrix with elements [m i, j ], i, j = 1, 2, 3, Note that for Gaussian, it holds a j = b j , j = 1, 2, 3 in ( 1) and (2). Figure 3 confirms the desired effect of this forwardbackward transition.This initialization uses the last input value and three last values from the forward pass to initiate the backward pass.

Direct Gabor Filtering
A 1D complex Gabor filter consists of a 1D Gaussian filter associated with some sigma and a complex exponential e kxω where x is position within the filter, ω is the filter frequency, and k 2 = −1 is the imaginary unit.YvV filter coefficients from ( 1) and (2) then change to complex with b YvV j being b j from (10) of [1].Both forward and backward passes happen in the complex domain.
The correct initialization for the YvV's forward pass is (4) with complex coefficients a j , Similarly for the YvV's forward-backward transition, following the derivation from [5] right from the beginning ((3) of [5]) in the complex domain, we arrived to ( 5) and ( 6), and we found that correct initialization matrix M is the Gaussian's matrix M element-wise multiplied with complex exponentials: Note that since in the case of Gabor filtering the forward and backward coefficients are not the same, a j / = b j , the direct substitution into the Gaussian's matrix M does not lead to the valid result.The matrix for Gabor filter must be evaluated with different a j and b j right from the start leading into an overly complicated matrix from which, after realizing that a j and b j differ only in the complex exponential (17), the result of ( 19) is derived.
A more elegant and short derivation exists.During their derivation, Triggs and Sdika [5] arrived to the implicit definition of general M for YvV recursive filtering, (18) of [5], in which I 1 is a matrix with a 1 in the top-left corner and zeros elsewhere, A and B "encodes" the forward and backward passes, respectively, for example, when k = l = 3: For details, refer to their original paper.We now consider again the Gaussian filtering as defined in Section 2, that is, with a j = b j , a j ∈ R, so that we can define A for Gabor filtering explicitly as A = e kω D −1 AD and similarly B = e −kω DBD −1 with Using (20) and generality of their solution, we have M = I 1 + B M A also for a Gabor filter.After substitution to the right hand side of it, the terms e −kω and e kω cancel themselves, and we multiply appropriately from both sides, Since the I 1 has zeros everywhere except for the top-left corner, when multiplying with any matrix, only the value at this corner is not zeroed in the result.Hence, Finally, the forward-backward transition for Gabor filter is expressed by means of Gaussian filtering, from which we realize that M = DMD.
Figure 4 shows an application of the presented forwardbackward transition.

Staged Gabor Filtering
The staged Gabor filtering makes use of the equality e q+w = e q • e w in the 1D convolution expression to yield the convolution The Gabor filtering is then split into three stages: Modulation of input data by e −klω , followed by a Gaussian convolution, and demodulation of the convolution result by e ktω in the end.This is a generally favored approach to a convolution with complex Gabor filter (Section 3) because it requires less real multiplications and additions per pixel [4] and is easier to implement compared to the direct Gabor filtering.
Both are a consequence of replacing truly complex Gabor convolution with two real Gaussian convolutions.A real constant input data i t = c, t = 1, . . ., n is no longer constant after a modulation which gives complex i mod t = ce −ktω .This is an important observation for the subsequent Gaussian filtering in the second stage since the modulated data i mod t is actually an input to it.In particular, the initialization for Gaussian filtering presented in Section 2 cannot be used because it expects that input data is extended with a real constant.In the following, we focus on the two real Gaussian convolutions using YvV's forward and backward passes as defined in (1)-( 3), a j = b j , but with initialization based on the presumption that input data line is extended with complex periodic function of ω.

Application to Higher
Dimensions.Before we get to an initialization of individual 1D convolution, we will first discuss its applicability in a higher dimensional space to ease the comprehension of the rest of this section.This is motivated by the fact that any Gaussian filtering can be realized as a cascade of several 1D Gaussian filtering [7,9].Note that the filtering happens after the modulation stage.Also note that an input to the next filtering in a cascade is actually an output of the preceding one.Consider, for example, a 2D Gabor filtering task with Gaussian envelope separable along the x and y axes and with frequency tuning of |(dx, dy)| in the direction of vector (dx, dy).One would typically compute the product i x,y • e −k(x•dx+y•d y)ω for a (28) But we may easily interleave the modulation with filtering The point of this change was only to show that it is correct to split the 2D, or generally nD, Gabor filtering to a cascade of pairs, each pair consisting of "1D" modulation in the direction of the subsequent 1D convolution and the 1D convolution itself.Hence, in the following subsections, we can assume that the 1D Gaussian filtering with complex period boundary extension was always preceded with an appropriate modulation.This is correct for any 1D filtering in a cascade even when modulation is only conducted prior the nD Gaussian filtering, as in (28).(cross marks).Only the amplitude and phase of the signal are changed.After the backward pass (circle marks), the phase is synchronized again with the modulated input as the Gaussian filtering is fully computed now.The period of signal is kept unchanged after each of the two passes.
the forward pass applied on complex periodic data i mod t = i 1 e −ktω , t = −∞, . . ., 0. Suppose that, based on observation depicted in Figure 5, the response u t has the shape of Re −ktω+kϕ where R is a complex gain and ϕ is some phase shift.Rewriting (1), we obtain which yields (after division by e −ktω+kϕ ) Note that R depends on ϕ, but the value of u t is independent of it.The expression in (32) requires the knowledge of the original value (prior to the modulation) of input data i 1 .Value of i 1 may be also computed by reverting the effect of modulation i 1 = i mod 1 e kω .Alternatively, we realize that the forward initialization values u −2 , u −1 , u 0 are adjacent in position to the first convolution position t = 1.In general, we seek initialization values u t− j , j = 1, 2, 3, for some 1D EURASIP Journal on Advances in Signal Processing convolution starting at position t with a value i mod t (recall that u t = R • e −ktω+kϕ ): We observe that (de)modulation terms including t and ϕ actually cancel out.The only term with t that remains is the i mod t which is the first value the 1D convolution receives on its input.The convolution then uses this value to compute the initialization values and starts the forward pass of the staged Gabor filtering with Note that the right-hand side of (18) for the direct Gabor filtering and the right-hand side of (32) for R are the same except for the phase shift e −kϕ .This evocates impression that we actually perform Gaussian filtering with initialization for the direct Gabor filtering.The handling of forward-backward transition shall support this observation to some extent as well.

Backward Initialization.
For the forward-backward transition, we first demodulate input data, so that it appears constant, and we may reuse previous result for the direct Gabor filtering from Section 3. We examine the original YvV filtering equations ( 1), (2) for Gaussian, that is, a j = b j , on a modulated input data i mod t = i n e −ktω , t = n + 1, . . ., ∞, after both equations are multiplied with demodulation term e ktω : a j e −k jω • v t+ j e ktω e k jω , (36) from which This pair of equations is the YvV Gabor filtering (compare with the first paragraph of Section 3) on constant input data i n .We may use the result for the direct Gabor forwardbackward transition where matrix M from ( 19) is used (39) We aim to rewrite (38) back for the modulated input data.The following equations are useful: In a similar way, (u n , u n−1 , u n−2 ) T and (v + , v + , v + ) T can be expressed.These equations enable us to express (38) only with u n , u + , v n , and v + .We multiply the whole equation with the matrix E from the left, This leaves the vector (v n , v n+1 , v n+2 ) T alone in the left-hand side of the equation.The right-hand side of it can be regarded as a subtraction of the two terms A − B: Regarding the term A, the multiplication of the three 3 × 3 matrices results in where M is the forward-backward transition matrix for Gaussian, that is, the real matrix from (5).
Regarding the term B, it is easy to show that M multiplied from right by a diagonal matrix is equal to M multiplied from left by the same diagonal matrix provided the diagonal is constant.Hence, we may swap the multiplication with M in (43).Furthermore, we let D = 1 − a 1 e kω − a 2 e 2kω − a 3 e 3kω and D be its complex conjugate, so that After some manipulations, we end up with in which the term in outer square brackets can be precomputed prior to the computation of convolution and M is the forward-backward transition matrix for Gabor filtering, that is, the complex matrix from (19).Moreover, since the vector (i mod ) T is constant, the precomputation can be simplified.For instance, the second row of B evaluates to Finally, to summarize the forward-backward transition M and m l, j are defined in Section 2. Figure 6 confirms the desired effect of this forward-backward transition.Note that the initialization for the backward pass works with i mod n , which is an (already modulated) input to the filtering (to the second stage), instead of working with i n , which is an (original, user's) input to the modulation (to the first stage).Thus, the original data i t , after being processed in the first stage, is not required any more during the filtering or in the initialization of it.

Theoretical and Experimental Comparison.
Let us make a comment on the solution by Bernardino and Santos-Victor [4], which was, to our best knowledge, so far the only published solution for the forward-backward transition during the staged Gabor filtering based on the YvV recursive filters.Though they approached it with the Z-transform, we found that our solution is equal to theirs both in terms of accuracy and speed (measured as number of multiplyadd operations per pixel, ops/px, required to process the transition).The major difference, however, is that while we establish initial values v n , v n+1 , and v n+2 for the backward pass only from the filter coefficients a 1 , a 2 , and a 3 , they require poles of Z-transform of the filter in addition to its coefficients.Their formulae for the transition simply require more input parameters.This does not limit applicability of their approach.Only for filters for which we know only their coefficients, for example, the "latest" YvV filters from 2002 [1], the poles must be computed additionally.Note that the poles are roots of a 3rd order polynomial (the 3rd order refers to (1) and ( 2) in which we fixed the order of filter recursivity) that is found in denominator of the filter's transfer function [10]."Our approach is closer to the spirit of Triggs and Sdika [5] since we are able to treat any general recursive filter that is according to the YvV scheme regardless of how the filter coefficients were derived and what its poles are."Note that Bernardino and Santos-Victor [4] used the "older" YvV filter from 1998 [11], where filter poles were explicitly computed to derive filter coefficients, so that all parameters for their forward-backward transition were easily available.Figure 7: Comparison of results of the same Gabor filter computed with different filtering without the zero-mean correction, denoted as with DC.Two input images were used: One was randomly generated, denoted as Set1, while the second is a copy of the first one with a constant value added to every image position.We may also notice how differently respective filtering variants respond to the change of mean value of input data, which is a result of different approximations involved in the respective variants.Only real parts of the complex filtering results are shown here.
We have also consulted and tested their C++ implementation, which the authors advertise in their original paper.We tested it over several filters on some randomly generated data.In each case, we inserted filter coefficients a i established by their implementation to our implementation (it cannot be done the other way round because we do not have and we don't even expect to have the poles at hand) to ensure the same filter is used and compared the results.For 1D tests and with respect to errors due to the floating-point arithmetics, both filtering implementations returned with the same values, not only within the border regions.Minor difference is that they actually compute values of v n+1 , v n+2 , and v n+3 needing, therefore, to conduct one iteration more of the backward pass, which is another 14 ops/px more for every transition yielding 94 ops/px in total in contrast to 54 ops/px required by our solution for every forwardbackward transition.
We argue that our solution appears to be easier to use as it is independent of filter poles.

Towards Zero-Mean Gabor Filter
"When implementing a Gabor filter, one should not forget to consider one of the classical results on Gabor filtering, since the real part of a filtering result is biased, see Figure 7." A generally acknowledged solution to remove the bias is to subtract a scaled Gaussian filtered input image from the real part [4,12,13].A Gaussian filter with exactly the same parameters as those of the Gaussian envelope of the Gabor filter is used.
In the staged filtering, enforcing the zero-mean property is nothing else but yet another exactly the same filtering with the same Gaussian this time on original nonmodulated input image [4].Indeed, the staged 1D zero-mean complex Gabor filtering involves three 1D convolutions with the same Gaussian filters (the same coefficients a j and b j ).Twice the (real) filtering operates on complex-modulated input with boundaries treated according to Section 4 and once the filtering operates on real input image with boundaries treated according to Section 2. The latter result is then scaled and subtracted from the real-part result of the first two filtering.
For the scale coefficient, one typically uses the value of Fourier transform of the Gaussian envelope at the Gabor filter's frequency [4,13].Actually, Fourier transform of the YvV recursive filter is used [4].Note that since the filter is an approximation to Gaussian filter, the explicit formula for Fourier transform of Gaussian can't be used for the computation of the scale coefficient.
Because we deal with complex Gabor filtering, we may alternatively obtain the scale coefficient from Fourier transform of the Gabor filter at the "zero" (the DC) frequency [12].We offer here a simple and fast to compute formula that is alternative and equal to the solution by Bernardino and Santos-Victor [4].We make use of the fact that transfer function of a YvV filter, given with the system in (1) and (2), is If the coefficients for some Gaussian are used, that is, a j = b j given in (3), then the expression BH(z) will become transfer function of the Gaussian filter.Note that it is a property of YvV Gaussian filters from 2002 [1] that B = (1+a 1 +a 2 +a 3 ) 2 , a j ∈ R and that 1+a 1 +a 2 +a 3 equals to B, (10) from [2], and to α, (12) from [11].Following Young et al. [1], we change

Figure 1 :
Figure 1: Example of an error after a direct Gaussian filtering on sample input data: n = 149, Sigma = 3, input extended constantly with the border value.A result of the incorrect forward-backward transition according to YvV is compared to the FIR implementation with constantly extended border as well.The FIR serves as a groundtruth result.The error spans 5 data positions.

Figure 2 :
Figure2: Example of errors after a direct (square marks) and staged (circle marks) Gabor filtering on sample input data: n = 149, Sigma = 3, period = 4 px, input extended constantly with the border value.For the staged cases, which is always a modulation, Gaussian filtering, and demodulation, we show that any forward-backward transition for direct Gaussian filtering (discussed in Section 2) shall not be used, despite that this transition may be correct for regular Gaussian filtering.The results on the forward-backward transition are compared to the FIR implementation with constantly extended border as well.The FIR serves as a ground-truth result.The error spans 7 data positions.Only real parts of the complex filtering results are shown here.

Figure 3 :
Figure 3: The forward-backward transition initiates correctly when compared to an FIR implementation.Compare it with Figure 1.

Figure 4 :
Figure 4: The forward-backward transition for the direct Gabor filtering initiates correctly when compared to an FIR implementation.Compare it with Figure 2.Only real parts of the complex filtering results are shown here.

4. 2 .Figure 5 :
Figure 5: Illustration of a result of the forward pass u t (square marks) during a Gaussian filtering on a modulated input i mod t

Figure 6 :
Figure 6: The forward-backward transition for the staged Gabor filtering initiates correctly when compared to an FIR implementation.Compare it with Figure 2.Only real parts of the complex filtering results are shown here.