8th WSEAS International Conference on SIGNAL PROCESSING, COMPUTATIONAL GEOMETRY and ARTIFICIAL VISION (ISCGAV’08) Rhodes, Greece, August 20-22, 2008

Arbitrarily-Oriented Anisotropic 3D Gaussian Filtering Computed with 1D Convolutions without Interpolation Vladim´ır Ulman Centre for Biomedical Image Analysis Faculty of Informatics, Masaryk University Botanick´a 68a, 602 00 Brno Czech Republic [email protected] Abstract: The paper presents a procedure that allows for good approximation of anisotropic arbitrarily-oriented arbitrarily-dimensional Gaussian filter. Even though the method is generally nD, it is demonstrated and discussed in 3D. It, essentially, substitutes a given 3D Gaussian by six 1D oriented recursive convolutions. Since many such orientations may be available, a set of fast constraints is presented to narrow the number of them and to select the best ones. Firstly, only integer elements are allowed in the directional vector of each orientation not to break translation-invariance of the method. Secondly, validations due to requirements on separability of the approximated filter are applied. Finally, the best directions are selected in the convolution test. We compare the method to a similar recent approach of [Lampert and Wirjadi, 2007] in 3D. A possibility how to improve the performance of bank filtering by using this method is outlined. Key–Words: Gaussian separability, spatial image filtering.

1

Introduction and Motivation

tion is required whenever the 1D convolution gets off the voxel grid. The interpolation pulls in the drawback of the Lampert’s method: the translation-variance — impulse responses vary as the their spatial positions change, see Lampert in Figure 1. Following the demonstration of this for 2D filtering by [6], where it is called positional variability, the similar behaviour can be observed for 3D filtering in Figure 3. Unfortunately, the average error magnitude changes with filter’s orientation. Such unwanted filtering noise can be completely suppressed when the filtering becomes translation-invariant. Relaxing the need for interpolation, by using only directions given by vectors with only integer elements, is the way to achieve that [6]. The key idea of our approach, derived from Lam and Shi [6], is to replace the nD naive-convolution with a cascade of m ≥ n 1D convolutions that can’t run off the voxel grid. This trades the interpolation for increased number of convolutions. However, as we shall see later in the Results section, this does not increase the computational time enormously while it improves accuracy. Moreover, as the number of convolutions per filter increases, the possibility to share partial results when computing many similar filters increases as well. Total number of convolutions may be decreased in this way when computing bank of filters. In this paper, we would like to extend the recent

Gaussian filtering is one of the most fundamental concepts in signal processing, without any doubt. Although many rivaling procedures have been developed, Gaussian filtering still remained among the most popular — owing in great part to its pleasant properties. Most notably, it was the separability that attracted when moving from 1D signal processing to nD allowing for relatively fast convolutions. Anyway, Gaussians started to gain attraction as early as in the 80’s [2]. Since 90’s, 1D recursive filtering joined the separability to form considerably fast yet reasonably accurate framework [3, 10, 5]. This still hadn’t permitted for effective computation of arbitrarily-oriented nD Gaussian, until recently [7]. The Lampert’s method [7] replaces the nD naiveconvolution, i.e. the O(k n mn ) algorithm for convolving image of size m with kernel of size k, with a cascade gn ∗v~n (

...

(g2 ∗v~2 (g1 ∗v~1 I)) )

(1)

of n carefully oriented convolutions ∗v~i with 1D filters gi . Input image is denoted as I. Each convolution runs along the direction v~i . Lampert’s method is shown to be optimal in the terms of minimum number of oriented convolutions together with the minimum number of interpolations necessary [7]. An interpola-

ISSN: 1790-5109

56

ISBN: 978-960-6766-95-4

8th WSEAS International Conference on SIGNAL PROCESSING, COMPUTATIONAL GEOMETRY and ARTIFICIAL VISION (ISCGAV’08) Rhodes, Greece, August 20-22, 2008

Lam and Shi

naive−convolution Lampert and Wirjadi

-10 -15

y

(1;0)

1.

1. x

y

Figure 1

y

-20

(1;0)

dB

1.

-20 -40 -60 -80 -100 -120 -140 0

x

x y

Figure 2 y

2.

(0.5;1)

2.

(0;1)

x z

y β

y

x

2.’

(0.5;1)

x z

3.

(1;1)

x

-25 -30 -35 -40 5 15 25 35 45 55 65 75 85

beta (deg)

10

20

30 beta 40 50 60 (deg) 70

x y

α

dB

y

alpha (deg) 80

90

Figure 3: Positional variability of responses computed according to Lampert and Wirjadi [7]. For every kernel configuration (α, β, 5, 2.5, 2.5): input images with differently positioned impulses were filtered and, afterwards, shifted so that original impulse positions aligned; we aggregated squared differences between voxel and mean values at every spatial coordinate. The sum was divided by the sum of squared means (total energy of the mean impulse) and 10 log() was computed from it to express in dB (deciBell). The left-hand-side plot displays the variability as a function of orientation in polar coordinates (β cos(α), β sin(α)). The right-hand-side plot shows variability profile along β by fixing α = 25 (the crosses in both plots).

x

Figure 1: Illustrations of the kernels used in the three methods for 2D and their applications on single impulses positioned at x or – x . The naive-convolution requires only one convolution step with the “full” kernel. The Lampert’s method separates the task into two 1D convolutions: the first runs along xaxis, the second is oriented. The last method uses three 1D convolutions, each stays within the voxel grid. Notice the different responses of the same kernel for single impulses if being positioned at x and – x owing to the interpolation in the situation 2’. Figure 2: A certain isoline of a quarter of Gaussian kernel. The upper plot shows it in the initial position while the other shows it after the rotations.

2

translation-invariant approach of [6] from 2D to nD, much like [7] extended [4], for comparison refer to Figure 1. Especially, we aim to describe a procedure for selecting the appropriate directional vectors and associated variances in order to be able to convolve given image with any 3D Gaussian. The Gaussian kernel given as (α, β, Σ, σ1 , σ2 ) is understood as follows:

Method

• filter’s mean is zero,

Our filtering method consists of m = (n(n + 1)/2) pairs, each of which associates a variance σi2 with a certain nD vector v~i . We call the nD vector the base vector and a m-tuple of these the base. Bases are drawn from the initial base set. We denote m pairs 2 ) as the admissible base when(~ vi , σi2 ), . . . , (v~m , σm ever they satisfy the separability constraint (explained in the next section) and σi2 ≥ 0, ∀i ∈ h1, mi. A cascade of 1D Gaussian convolutions then runs each along its base vector v~i with σi2 . The filtering setup of our method for an arbitrary input nD Gaussian is such admissible base that yields filtering result closest to the nD naive-convolution with this particular kernel.

• variance in the main direction is Σ2 , variances in the two orthogonal directions are σ12 , σ22 ≤ Σ2 ,

2.1

• initially, filter’s axes are identical to the axes in the Euclidean space; without loss of generality, filter’s main direction is along the x-axis,

Gaussian decomposition

Assume, we have parameters of arbitrary input nD Gaussian compactly stored in the covariance matrix C. The input nD Gaussian kernel g n (~x, C) has the form

• filter is firstly rotated around the z-axis by α, • filter is then rotated around newly positioned yaxis by β; refer to Figure 2.

g n (~x, C) = ISSN: 1790-5109

57

1 n 2

1 T −1 C ~ x

1/2

(2π) |C|

· e− 2 ~x

ISBN: 978-960-6766-95-4

(2)

8th WSEAS International Conference on SIGNAL PROCESSING, COMPUTATIONAL GEOMETRY and ARTIFICIAL VISION (ISCGAV’08) Rhodes, Greece, August 20-22, 2008

where ~x = (x1 , . . . , xn )T is a vector with spatial position within the coordinate system, typically Euclidean, in which convolutions shall take place. Let us select the following m base vectors v~i = (a(i,1) , . . . , a(i,n) )T , i = 1, . . . , m, a(i,k) ∈ Z, Z be the domain of integer numbers. We can establish the following transition matrix between ~x and ~u = (u1 , . . . , um )T : ~x = A · ~u = [v~1

...

v~m ] · ~u.

the input filter’s main axis. Unfortunately for the 3D case, we were not able to find any reasonable criterion to proceed similarly. Note that, when moving into 3D, the number of base vectors as well as Gaussian parameters increases. Hence, we adopted the approach of searching the entire initial base set for the most accurate admissible base. This involves, firstly, to decide whether a base is admissible by examining a solution of the equation (5) and then, eventually, to compute the error, which will be explained in the next subsection. The full search ensures that the best available base is used. For the 3D case, we use the following 91 bases for the initial base set. Each base consists of base vectors from the classes base, A, B, C and D:

(3)

2 ) be the diagonal coLet E = diag(σ12 , . . . , σm variance matrix of some Gaussian kernel in the mD separable along v~i . Consider the development of the argument of the Gaussian’s exponential:

~uT E −1 ~u = (A−1 ~x)T E −1 (A−1 ~x) = ~xT (AEAT )−1 ~x (4) The right-hand-side of the equation explains the shape of the separable mD Gaussian in the context of coordinate system bond with the computing platform. We observe that C = AEAT , which can be further developed, following [6], into the separability constraint: X

C=

h

base : A: B: C: D:

Each base is composed using the rules: ~v1,2,3 ∈ base, v~4 ∈ A ∪ D, v~5 ∈ B ∪ D, v~6 ∈ C ∪ D, no vector is repeated, up to one vector from D is used, maximum number of vectors containing element 2 is limited to 2.

i

σi2 v~k · v~k T .

(5)

k=1,...,m

Since C is the Gaussian’s covariance matrix, by the definition, it is symmetric as such. It is easy to show that the right-hand-side of the eq. (5) is symmetric as well. For the eq. (5) to hold, it is enough to compare diagonals and both upper triangles, which constitutes 2 un(n(n + 1)/2) linear equations with σ12 , . . . , σm knowns — C was given and v~1 , . . . , v~m were fixed. Note that we have chosen m = (n(n + 1)/2) previously. The system of linear equations based on the equation (5), basically, completes the base m-tuple by defining appropriate σi2 to v~i so that the base mathematically approximates input Gaussian given by the C. However, the base is technically realizable only if it is admissible, i.e. σi2 ≥ 0. Hence, cascaded convolutions based on an admissible base result in something close to the nD Gaussian. Every i-th 1D convolution computes 1D Gaussian g 1 (~u, [σi ]) that iterates only ui in the ~u; the corresponding ~x position is given by the eq. (3). We skip i-th convolution ∗v~i if σi2 = 0.

2.2

2.3

Filtering setup determination

To determine the base which explains the input filter the best, the admissible bases from bases in the initial base set need to be computed and tested afterwards. In the test, an impulse response of the examined admissible base is computed and compared to a reference impulse response of the input Gaussian kernel. Thus, every admissible base requires one run of the cascaded convolutions. The comparison is quantified with the NSE error measure, defined in eq. (8). The final filtering setup is the admissible base showing the smallest deviation. In order to minimize the computational load of the tests, the size of a test image with sample impulse is changing with parameters of the input Gaussian. Namely, the image size of 8Σ × 8σ1 × 8σ2 is used for (0, 0, Σ, σ1 , σ2 ). For the oriented kernel, the image size approximates the size of an axis-aligned bounding box that wraps the kernel in the same manner. Minimum image size is 9 × 9 × 9. Note that the test image size is constant during the tests because it depends solely on the input parameters. Also note that an impulse response resulting from a single-voxel impulse naive-convolved with sampled kernel results in exactly the same sampled kernel multiplied by the magnitude of the impulse. Hence, the reference test image can be precomputed fast.

Computation of admissible bases

It is apparent from the previous section that the goal is to choose such base vectors that form an admissible base, ideally directly the filtering setup. In the 2D case, where m = 3, Lam and Shi [6] seemed to manually select appropriate base vectors by considering firstly the ratio σ1 /Σ and secondly the orientation of

ISSN: 1790-5109

{(1, 0, 0), (0, 1, 0), (0, 0, 1)}, {(1, 1, 0), (2, 1, 0), (1, 2, 0)}, {(1, 0, 1), (2, 0, 1), (1, 0, 2)}, {(0, 1, 1), (0, 2, 1), (0, 1, 2)}, {(1, 1, 1), (2, 1, 1), (1, 2, 1), (1, 1, 2)}. (6)

58

ISBN: 978-960-6766-95-4

8th WSEAS International Conference on SIGNAL PROCESSING, COMPUTATIONAL GEOMETRY and ARTIFICIAL VISION (ISCGAV’08) Rhodes, Greece, August 20-22, 2008

3

Results and Discussion

set size 27132 560 91

In the rest of this paper, we considered only 3D Gaussians with the following restrictions. Both variances for the axes orthogonal to the main direction were always equal. Every Gaussian was defined by some j = 1, . . . , 3610 as (αj , βj , Σj , σj , σj ). Both αj and βj were allowed to span from 0◦ to 90◦ in steps of 5◦ . For every such (αj , βj ) pair, we had used all Σj from {2, 3, 4, 5} and consequently all σj from {1, 2, . . . , Σj − 1}. In this section, we will show results on the determination of filtering setup before results on the performance of cascaded convolutions.

3.1

RMSE 2.87 ± 4.50 3.31 ± 5.81 3.36 ± 5.88

Table 1: The averaged errors of filtering setups, over all configurations (αj , βj , Σj , σj , σj )j=1,...,3610 , achieved for different initial base sets. Note that reducing the size of the initial base has not increased the error considerably. compute the common parts of similar filters only once. Clearly, the specification of filters should be known apriori. Originally, we used 19 base vectors whose elements were from {0, 1, 2} and up to one element with number 2 was allowed in every base vector. This constituted what we called the full initial base set of 27132 bases. The evaluation of linear equations, eq. (5), was relatively fast. It didn’t incur any substantial slowdown as the search over the initial base set required less than 0.085sec for a given Gaussian in our experiments. When considering all 3610 tested configurations, the total number of admissible bases was 1646871, which gave only ef = 1646871/(27132 × 3610) = 1.68% (eq. (7)). The convolution tests, however, demanded 0.563sec per filter on average. After an analysis of the admissible bases resulting from the full initial base set, 3 vectors in every base were fixed, namely (1, 0, 0), (0, 1, 0) and (0, 0, 1). The size of new initial base set dropped to 560, similarly the total number of admissible bases dropped from 1646871 to 109507. The efficiency ratio increased to ef = 5.42%. The vast majority of discarded bases would not constitute filtering setups anyway as the overall accuracy hadn’t worsen dramatically, see Table 1. The final initial base set of 91 bases allowed for 100383 admissible bases in total and for the efficiency of 30.56%. The number of all admissible bases as well as the accuracy have been roughly kept, Table 1. The effect of the initial set’s reduction on the distribution of errors is displayed in Figure 4. Determination of all admissible bases for a single Gaussian has reduced to less than 0.0003sec, convolution tests to 0.037sec.

Narrowing the initial base set

As it was already explained in subsection 2.3, convolution test was conducted for every admissible base in order to determine the best one — the best configuration for the cascaded convolutions. This testing proved itself to be critical step of the method in terms of time consumption, especially when many admissible bases emerged. Narrowing the initial base set, to decrease the number of admissible bases, while preserving the level of accuracy, by not discarding good bases from it, had significant impact on the usability of the entire method. We had established the efficiency ratio ef for every initial base set as a ratio between an average number of admissible bases for every of 3610 tested Gaussians and the number of all possible bases: total number of admissible bases × 100%. size of initial base set × 3610 (7) The ratio balances between the two counter-effects. The smaller the ratio is, the more unnecessary bases are present in the initial base set — the computational load is needlessly higher. The greater the number is, the worse selectivity of the Gaussian decomposition is — many admissible bases exist and are passed to the convolution test. The best practice seems to raise the efficiency ratio while to lower the size of the initial base set. By decreasing the size of the initial base set also decreases the number of total admissible bases because there is simply less bases to choose from. This, in turn, decreases the number of convolution tests, which results in a faster determination of the filtering setup. If the numerator decreases slower than the denumerator of eq. (7), the efficiency raises. Such reasoning has also positive influence on number of convolutions with a bank of filters. Small set of initial bases forces for selecting setups that may share some pairs (~ vi , σi2 ). Since there is no restriction on the order of 1D convolutions in the method, we may ef =

ISSN: 1790-5109

NSE (dB) -42.78 ± 11.69 -41.90 ± 13.91 -41.46 ± 13.48

3.2

Results of filtering

As emphasized in the introductory section, we see the proposed method as the competitor to the method of Lampert and Wirjadi [7]. However, implementation of both methods shared essential parts since both utilized cascade of 1D convolutions, each along arbitrarily oriented axis. We implemented the 1D Young’s recur-

59

ISBN: 978-960-6766-95-4

8th WSEAS International Conference on SIGNAL PROCESSING, COMPUTATIONAL GEOMETRY and ARTIFICIAL VISION (ISCGAV’08) Rhodes, Greece, August 20-22, 2008

Σj -σj 2-1 3-1 3-2 4-1 4-2 4-3 5-1 5-2 5-3 5-4 2-1 3-2 5-2

sive filters [11] with improved boundary handling [9]. The main distinction was that the proposed method used more convolutions traded for no interpolations at all while the Lampert’s method required linear and bilinear interpolations for 3D. In this subsection, we compared our method only to their and to the naiveconvolution. The first study one typically starts with is to preview the result itself. We computed 3D impulse responses for particular Gaussian, thresholded them and plotted the result in Figure 5. The immediate observation follows the one made by Lam and Shi [6]: the isolines of the Lampert’s method tend to appear more rectangular when compared to the true Gaussian. On the other hand, the proposed method seems to oversmooth the filtering as the outer contour is broader and further to the centre. Quality of a filtering method was described by an error measure which, in this case, compared a cascade-convolved image with a naive-convolved one. Absolute relative error was chosen in [8] whereas maximum absolute error or root of squared differences (RMSE) were chosen in [10, 11, 4, 7]. We opted for RMSE and the normalized squared error (NSE). The NSE, essentially the error-to-signal ratio as in [1, 6], gives the average squared magnitude of errors in the examined image T with respect to the total energy of the reference image R: P

NSE (T, R) = 10 log10

− R(~x))2 . x))2 ~ x (R (~

x) ~ x (T (~ P

this paper (dB) -25.87 ± 9.39 -26.50 ± 9.79 -45.88 ± 3.95 -25.53 ± 10.54 -46.49 ± 2.38 -51.66 ± 1.51 -25.87 ± 11.07 -47.05 ± 2.51 -50.59 ± 2.04 -52.88 ± 1.63 -103.81 ± 6.84 -100.46 ± 3.06 -82.88 ± 4.96

Table 2: The averaged errors achieved for different settings of filter’s smoothing. The numbers are averaged, over all orientations, NSE errors computed against 3D naive-convolutions. The first 10 lines were evaluated on the impulse image whereas the last 3 lines on a sample image. Note that the lower the value is, the better the result is. true for sample image, which is also evident both from the last 3 rows of the table and from the figure. On the other hand, the proposed method at the same smoothing ratio σj /Σj = 0.5 displayed the same feature of decreased accuracy at α ∈ h5, 85i, β ∈ {80, 85} both in Figures 6 and 4, i.e. on both types of convolved images. This could signal that the smoothing-related part of the error dominates the orientation-related part. The interpolation in Lampert’s method adds a constant performance penalty attached to every input voxel. Nevertheless, in 3D we have achieved 76ms and 110ms on average for the Lampert’s and the proposed method, respectively. Only convolutions were measured to see the difference between interpolations and additional convolutions. The average was over all configurations (αj , βj , Σj , σj , σj )j=1,...,3610 for a float-valued sample image of 130 × 130 × 90px (5.8MB uncompressed). Note that the Σj and σj didn’t affect the speed of the filtering owing to the employed recursive filters. All the times in the paper were measured on Linux 2.6 on a PC-type computer with Intel(R) Core(TM)2 Quad 2.4GHz cpu with 4GB DDR2 memory. We see a potential of the proposed method in a decrease of total number of 1D convolutions when computing bank of similar 3D Gaussians. This could alleviate the total computational demand of the bank. If we allowed to use more of the tested admissible bases that passed over reasonable minimum NSE quality

(8)

The NSE, unlike the RMSE, readily gives the magnitude of the error regardless the average intensity in the filtered image. The accuracy is summarized in Table 2. The proposed method seems to outperform very decently the competing approach except for filters with small variances. In these cases, the method even doesn’t keep filtering quality at any stable level, which is demonstrated with high standard deviations (and also in Figure 6). However, the opposite seems to be true when applied on a sample image. Distinctions in the smoothing parameters seemed to be available when studying errors in Table 2 and in Figure 6. Notice that two data layers emerged in both plots in the figure, which, together with the table data, suggests that different Σj -σj combinations may provide different levels of accuracy. Indeed, the core of our implementations of both Lampert’s and the proposed method’s are Young’s et al. 1D recursive filters. These were evaluated on impulse responses and reported to give better accuracy for greater sigmas [10, 11]. This is acknowledged in the first 10 rows in Table 2. However, the opposite tendency seems to be

ISSN: 1790-5109

Lampert (dB) -26.82 ± 3.76 -27.87 ± 3.97 -41.71 ± 2.27 -28.47 ± 4.13 -42.50 ± 2.44 -48.88 ± 1.21 -28.52 ± 4.24 -42.95 ± 2.50 -49.24 ± 1.17 -50.81 ± 0.92 -95.09 ± 4.92 -94.45 ± 2.95 -82.09 ± 4.63

60

ISBN: 978-960-6766-95-4

-30

90 85 80 75 70 65 60 55 50 45 40 35 30 25 20 15 10 5 0

-35

beta (deg)

beta (deg)

8th WSEAS International Conference on SIGNAL PROCESSING, COMPUTATIONAL GEOMETRY and ARTIFICIAL VISION (ISCGAV’08) Rhodes, Greece, August 20-22, 2008

-40

-45

-50

-55

-30

90 85 80 75 70 65 60 55 50 45 40 35 30 25 20 15 10 5 0

-35

-40

-45

-50

-55

0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90

0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90

alpha (deg)

alpha (deg)

Figure 4: The NSE error, eq. (8), distribution for the proposed method over arbitrarily oriented 3D Gaussians (αj , βj , 4, 2, 2) computed on an impulse image. The mapping of a color to an error, expressed in dB, is described in the vertical bars. The left image shows errors when the full initial base set of 27132 bases was used while the right image is with only 91 bases considered. The error has increased predominantly in the central, brighter, region.

Figure 5: 3D impulse responses, shown as 2D cuts, of Gaussian (30, 20, 4, 2, 2) were computed by using, from left to right, Lampert and Wirjadi [7], naive-convolution and the proposed method. Each filtered image was thresholded with two fixed narrow intervals to display two isolines of the same response.

-70 -75 -80 -85 -90 -95 -100 -105 -110

-70 -75 -80 -85 -90 -95 -100 -105 -110

dB

dB

Lampert for (alpha,beta,2,1) Lampert for (alpha,beta,5,2) 0 30 60 alpha 90 (deg)

0

10

20

30

40 50 60 beta (deg)

70

proposed method for (alpha,beta,2,1) proposed method for (alpha,beta,5,2) 80

0 30 60 alpha 90 (deg)

90

0

10

20

30

40 50 60 beta (deg)

70

80

90

Figure 6: NSE error distributions computed on a sample image with oriented 3D Gaussians (αj , βj , 2, 1, 1) and (αj , βj , 5, 2, 2). Errors for the former and latter filters are displayed with down- and up-triangles, respectively.

ISSN: 1790-5109

61

ISBN: 978-960-6766-95-4

8th WSEAS International Conference on SIGNAL PROCESSING, COMPUTATIONAL GEOMETRY and ARTIFICIAL VISION (ISCGAV’08) Rhodes, Greece, August 20-22, 2008

-25

in the processing pipeline, we also indicated a possible application of a bank filtering whose performance may be improved with the method. We have partly outlined the way to achieve that. Last but not least, C++ routines are freely available at the address: http://cbia.fi.muni.cz/projects/ fast-gaussian-filtering.html.

(30,20,3,2,2) (60,20,3,2,2) accuracy at -40dB

-30 -35 -40 -45 -50 0

5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90

Acknowledgments: I would like to thank my colleagues Michal Kozubek and, especially, Pavel Matula. This work has been supported by the Ministry of Education of the Czech Republic (Grants No. MSM0021622419, LC535 and 2B06052).

Figure 7: The diagram of admissible bases and accuracies they provide. The abscissa shows ID of every base in the initial base set. Different filters may use the same base but with some different σi2 . threshold, instead of allowing the use of a single best one, we would usually have a few of these to choose from for every Gaussian. For example, for filters (30, 20, 3, 2, 2) and (60, 20, 3, 2, 2) we may find many such admissible bases, Figure 7. The goal then would be to find one such base for each of the two filters that would share the most of the pairs (~ vi , σi2 ). In this example, we found v~3 = (0, 0, 1) and v~4 = (1, 1, 0) with σ3 = 1.86743 and σ4 = 1.38268 in both admissible bases. The corresponding convolutions are adepts for being placed at the beginning of both cascades so that the result can be shared between filters with minimum decrease in accuracy. The computation of two 1D convolutions can be saved in this way.

References: [1] A. Bernardino and J. Santos-Victor. Fast IIR isotropic 2-d complex Gabor filters with boundary initialization. IEEE Transactions on Image Processing, 15(11), 2006, pp. 3338–3348. [2] P.J. Burt. Fast algorithms for estimating local image properties. CVGIP, 21(3), March 1983, pp. 368–382. [3] R. Deriche. Recursively implementing the Gaussian and its derivatives. Technical Report 1893, INRIA, May 1993. [4] Jan-Mark Geusebroek, Arnold W. M. Smeulders, and J. van de Weijer. Fast anisotropic Gauss filtering. In Proceedings of the 7th ECCV, 2002, pp. 99–112. [5] J. S. Jin and Y. Gao. Recursive implementation of log filtering. Real-Time Imaging, 3(1), 1997, pp. 59–65. [6] S.Y.M. Lam and B.E. Shi. Recursive anisotropic 2-d Gaussian filtering based on a triple-axis decomposition. IEEE Trans. on Image Processing, 16(7), 2007, pp. 1925–30. [7] C.H. Lampert and O. Wirjadi. An optimal nonorthogonal separation of the anisotropic Gaussian convolution filter. IEEE Transactions on Image Processing, 15(11), 2006, pp. 3501– 3513. [8] S. Tan, J.L. Dale, and A. Johnston. Performance of three recursive algorithms for fast spacevariant Gaussian filtering. Real-Time Imaging, 9(3), 2003, pp. 215–228. [9] B. Triggs and M. Sdika. Boundary conditions for Young - van Vliet recursive filtering. IEEE Transactions on Signal Processing, 54(5), 2006. [10] I. T. Young and L. J. van Vliet. Recursive implementation of the Gaussian filter. Signal processing, 44(2), 1995, pp. 139–151. [11] I.T. Young, L.J. van Vliet, and M. van Ginkel. Recursive Gabor filtering. Signal processing, 50(11), 2002, pp. 2798–2805.

4 Conclusion The purpose of this paper was to give an insight into the problem of separately convolving anisotropic arbitrarily-oriented 3D Gaussian filter in the spatial domain. We have given brief overview of the situation in the field. The most importantly, we described the procedure of turning input 3D Gaussian parameters into parameters for a translation-invariant cascade of 1D oriented convolutions. We consider this to be our contribution to the field. We also attempted to discuss stumbling aspects of the procedure. We compared the cascade-convolved filtering results to the naive-convolved ones. For larger variances, the proposed method allowed for slightly improved accuracy over the Lampert’s method. On the other hand, the former worked up to 1.5 times slower than its competitor. For the demonstration and evaluation purposes, we had posed certain restrictions and limited the generality of input Gaussian to some 3610 filters. However, the proposed method really has the potential to compute truly every Gaussian filter. It has also the ability to scale into arbitrary nD Gaussian filtering. Since Gaussian filtering is rarely the final task

ISSN: 1790-5109

62

ISBN: 978-960-6766-95-4