JP2016122430A - Image filter arithmetic device, gaussian kernel arithmetic device, and program - Google Patents

Abstract

PROBLEM TO BE SOLVED: To provide an image filter arithmetic device in which calculation costs are smaller than conventional, and with which it is possible to control a calculation error as desired.SOLUTION: Provided is an image filter arithmetic device for performing image smoothing on an input image {x(p)} by the convolution arithmetic of a space kernel gand a range kernel g, the image filter arithmetic device comprising: approximation parameter optimization means for determining the set (T, K) of a truncation number K and a significant cycle T in an approximate range kernel g^in which the range kernel gis Fourier-series approximated for a scale parameter σ, a permissible error τ, and a dynamic range R, so that an approximation error becomes τ or less; expansion coefficient arithmetic means for calculating the expansion coefficient {a, and so on, a} of the approximate range kernel g^; and filter arithmetic means for executing the convolution arithmetic using the expansion coefficient {a, and so on, a}.SELECTED DRAWING: Figure 12

Description

  The present invention relates to a Gaussian kernel arithmetic device that performs an operation of Gaussian kernel g (σ, x) on an input sample value x and an input scale parameter σ, and an image filter arithmetic device using the same. About.

  In the field of image processing technology, a bilateral filter is known as a filter for reducing noise while preserving edges. In general, BLF is defined as:

(Definition 1) (Bilateral filter)
Consider an input image (target image) of a D-dimensional grayscale image. Z is a set of whole integers, Ω⊂Z ^D is an image domain, p∈Ω is a pixel position (vector), and N _p ⊂Ω is a neighboring region (filter window region) of the pixel position p. To do. Let R be a set of all real numbers, and let x (p) and z (p) εR be the pixel values of the pixel at the position p (hereinafter “pixel p”) of the input image and output image, respectively. At this time, a filter defined by the following equation is referred to as a “bilateral filter (BLF)”.

Here, g (σ, x) is a Gaussian kernel, σ _s is a scale parameter (standard deviation) that determines the correlation distance in the spatial region, and σ _r is a scale parameter that determines the correlation distance in the pixel value region. (Standard deviation), D is the number of dimensions of the vector q, and ‖q is the norm (l ₂ -norm) of the vector q. In addition, _{g s} a "space kernel (spatial kernel)", the _{g r} referred to as "range-kernel (range kernel)".
(End of definition)

  BLF is a basic tool for edge preserving smoothing in image processing, computer vision, computer graphics, and the like. BLF is designed based on the clear concept of defining filter coefficients from two laterals (ie, pixel position and pixcel intensity). BLF is a typical other edge-preserving smoothing technique such as anisotropic diffusion (non-patent document 6) or weighted least square (non-patent document 7). In contrast, a high smoothing effect can be obtained without using an iterative operation process. Therefore, BLF is widely used in various applications such as noise removal, high dynamic range imaging (HDR), and stereo matching.

  Conventionally, in the BLF, various extensions have been performed in order to improve the smoothing effect and reduce the calculation cost.

  First, as a natural extension of BLF, a cross / joint extension (cross BLF / joint BLF) has been devised (Non-Patent Documents 8 and 9). In filtering an input image with noise, the BLF cancels out the smoothing effect to some extent in order to define filter coefficients from the input image with noise. On the other hand, in the cross / joint BLF, a filter coefficient is defined using a guidance image shot under shooting conditions different from the input image. Thereby, a higher smoothing effect can be obtained.

In the BLF of Expression (1), g _r (•) changes depending on the pixel value x (•) of the input image, and therefore the filter coefficient dynamically changes in each target pixel. The original BLF calculates the filter coefficient from x (•), whereas the cross / joint BLF calculates the filter coefficient from the pixel value y (•) of the guidance image corresponding to x (•).

(Definition 2) (Cross / Joint Bilateral Filter)
Let y (p) εR be the pixel value at position p of the guidance image. At this time, the filter defined by the following equation is called a “cross / joint extension of the BLF”.

(End of definition)

  On the other hand, as an algorithmic extension for measuring a decrease in calculation cost, several high-speed BLFs including fixed-time BLF (O (1) BLF) have been devised so far (Non-Patent Documents 10, 11, 13-). 16). In the original BLF shown in Expression (1a), the calculation cost per pixel depends on the number of pixels M (= | Ω |) in the filter window region Ω (that is, the calculation amount is O (M)). ). Therefore, a large calculation cost is required. Considering the trend toward higher resolution and higher-dimensional image processing in recent years, it is an important technical problem to reduce the cost and speed of BLF. The high-speed BLF was devised to deal with this problem, and has been improved so that it can be executed with a calculation amount of O (log M) or O (1) at the expense of some accuracy. is there.

  In the conventional high-speed BLF, various improvements are made based on the fact that the BLF can be generalized as convolution in a high-dimensional space (bilateral space) measured by the space axis and the range axis. For linear filters, many high-speed algorithms have been devised in the past, but a large memory cost is required to directly handle bilateral space. Therefore, the conventional high-speed BLF algorithm solves this problem by various strategies.

  Non-Patent Document 10 discloses a high-speed BLF that simplifies a bilateral space by limiting a box spatial filter and using a histogram method.

  Non-Patent Document 11 discloses a fixed time O (1) BLF based on an approach similar to the histogram method (Non-Patent Document 12). Non-Patent Document 11 describes a polynomial range kernel using an arbitrary spatial filter and a Taylor series.

  Non-Patent Document 13 describes a high-speed BLF algorithm that approximates quantization and pixel-unit linear interpolation on both the space axis and the range axis, and performs linear filtering by fast Fourier transform (FFT). .

  Non-Patent Document 14 describes a speed-up BLF algorithm (Yang approximation algorithm) that is speeded up by using O (1) linear filtering by quantization and interpolation only in the range axis.

  Non-Patent Document 15 describes a more general downsizing approach of BLF using a three-dimensional FFT.

  Non-Patent Document 16 discloses a method of calculating BLF as a Gaussian transformation in a bilateral space and performing a fast Gaussian transformation (FGT).

  Non-Patent Documents 1-3 disclose O (1) BLF, which is speeded up by approximating a Gaussian range kernel with a triangular function and a polynomial function. According to this method, a translation-variant problem caused by bilateral space quantization or interpolation can be solved. The basis of this idea is the function's shiftability (the property that the function form is invariant even if the function is translated by an arbitrary amount by approximating it with the sum of several basis functions). (See formula (6) below).

  In Non-Patent Document 1-3, the Gaussian range kernel g (σ, x) is approximated by a power of a cosine function or a power of a polynomial function using the following expression (3a) or (3b) ( (Hereinafter referred to as “Chaudhury approximation”).

Equation (3b) can be derived by applying the definition of Napier number as it is. Equation (3a) is derived by expanding the cosine into a Taylor expansion and reducing it to equation (3b) under the condition of N >> 1. The condition of N> N ₀ in the expression (3a) is a condition for guaranteeing that the approximate expression is positive and monotonous in the domain [−T / 2, T / 2] of x.

  Non-Patent Document 2 describes the following approximate expression (4) as an improved expression that further accelerates the convergence of expression (3a).

  Further, in Non-Patent Documents 4 and 5, as another Gaussian kernel approximation formula, Gaussian kernel g (σ, x) is obtained by multiplying the cosine function polynomial and the exponential function by the following formula (5). An approximation method is described (hereinafter referred to as “Deriche approximation”).

Expression (5) is a fitting expression in the domain x∈ [0, ∞], and a ₀ , a ₁ , b ₀ , b ₁ , c ₀ , c ₁ , w ₀ , and w ₁ are fitting constants.

By the way, when the Gaussian kernel approximation calculation method is applied to a bilateral filter, the approximation formula of the Gaussian kernel can be translated in order to use a constant time filter (a filter whose calculation amount does not depend on the size M of the filter window). Must be (shiftable) (Non-Patent Document 3). Here, the function φ (x) (x∈R ^d ) is translatable, so that there is a finite number of functions φ (x−τ) for an arbitrary translational movement τ (∈R ^d ). Linear sum of basis functions φ ₁ (z), ..., φ _N (z)

It is expressed by. If the approximate expression of the Gaussian kernel can be translated, the sum of the numerator and denominator of the bilateral filter (1a) can be reduced to a moving sum (non-standardized moving average) expression. This is because the calculation using fixed time or the O (1) algorithm is possible by calculation using street recursion (Non-patent Document 3).

  Any of the approximate expressions (3a), (4), and (5) can be translated, and can be configured as a constant-time filter when applied to a bilateral filter.

K. N. Chaudhury, D. Sage, and M. Unser, "Fast O (1) bilateral filtering using trigonometric range kernels", IEEE Transactions on Image Processing, vol.20 (12), pp.3376-3382, 2011. K. N. Chaudhury, "Acceleration of the shiftable O (1) algorithm for bilateral filtering and non-local means", IEEE Transactions on Image Processing, vol. 22 (4), pp. 1291-1300, 2013. K. N. Chaudhury, "Constant-time filtering using shiftable kernels", IEEE Signal Processing Letters, Vol.18, no.11, pp.651-654, 2011. Rachid Deriche, "Recursively implementing the Gaussian and its derivatives", Technical report, INRIA, 1993. Tan, Jason L. Dale, and Alan Johnston, “Performance of three recursive algorithms for fast space-variant Gaussian filtering”. Real-Time Imaging 9, pp.251-228, 2003. P. Perona, J. Malik, “Scale-space and edge detection using anisotropic diffusion”, IEEE Trans. On Pattern Recognition and Machine Intelligence, vol. 12, no. 7, pp. 629-639, 1990. R.L.Lagendijk, J. Biemond, and D.E.Boekee, "Regularized iterative image restoration with ringing reduction", IEEE Transactions on Acoustics, Vol.36, no.12, pp.1874-1888, 1998. G. Petschnigg, et al., "Digital Photography with Flash and No-Flash Image Pairs", ACM Trans. On Graph. (Proc. Of ACM SIGGRAPH 2004), Vol.23 no.3, pp.664-672, 2004 . E. Eisemann and F. Durand, "Flash Photography Enhancement via Intrinsic Relighting", ACM Trans. On Graph. (Proc. Of ACM SIGGRAPH 2004), Vol.23, no.3, pp.673-678, 2004. B. Weiss, "Fast median and bilateral filtering", ACM Trans. On Graph. (Proc. Of ACM SIGGRAPH 2006), Vol.25 no.3, pp.519-526, 2006. F. Porikli, "Constant time O (1) bilateral filtering", IEEE Conference on Computer Vision and Pattern Recognition 2008, pp.1-8, 2008. F. Porikli, "Integral histogram: a fast way to extract histograms in Cartesian spaces", IEEE Computer Soc. Conf. On Computer Vision and Pattern Recognition 2005, Vol.1, pp.829-836, 2005. F. Durand and J. Dorsey, "Fast Bilateral Filtering for the Display of High-Dynamic-Range Images", ACM Trans. On Graph. (Proc. Of ACM SIGGRAPH 2002), Vol.21, no.3, pp.257 -266, 2002. Q. Yang, K.-H. Tan, N. Ahuja, "Real-Time O (1) Bilateral Filtering", IEEE Conf. On Comput. Vis. And Pattern Recognit. 2009 (CVPR 2009), pp.557-564 , 2009. S. Paris and F. Durand, "A Fast Approximation of the Bilateral Filter using a Signal Processing Approach", Int. J. of Comput. Vis. (IJCV), Vol.81 no.1, pp.24-52, 2009 . S. Yoshizawa, A. Belyaev, and H. Yokota, "Fast Gauss Bilateral Filtering", Computer Graphics Forum, Vol.29, no.1, pp.60-74, 2010. W. Ray and R. Driver, "Further decomposition of the Karhunen-Loeve series representation of a stationary random process", IEEE Transactions on Information Theory, Vol.16, no.6, pp.663-668, 1970. K. Sugimoto and S. Kamata, "Fast Gaussian filter with second-order shift property of DCT-5", Proc. IEEE Int. Coef. Image Process. (ICIP), pp.514-518, 2013. K. Sugimoto and S. Kamata, "Fast Gaussian filter using DCT-V", Proc. IEEE Int. Coef. Image Process. (ICIP), pp.514-518, 2013. R. Deriche, "Fast algorithms for low-level vision", IEEE 9th Int. Conf. On Pattern Recognition, Vol.1, pp.434-438, 1988. I. T. Young, L. J. Vliet, "Recursive implementation of the Gaussian filter", Signal Processing, Vol.44, no.2, pp.139-151, 1995. L. J. Vliet, I. T. Young and P. W. Verbeek, "Recursive Gaussian Derivative Filters", IEEE Comp. Soc. Press, Vol.1, pp.509-514, 1998. G. Farnemack and C. F. Westin, "Improving Deriche-style Recursive Gaussian Filters", J. Math. Imaging Vol.26, pp.293-299, 2006.

  As shown in Non-Patent Documents 15 and 16, the D-dimensional BLF can be expanded to a ratio of (D + 1) -dimensional convolution as follows by introducing a new intensity axis t.

Here, δ (x) is the Kronecker delta.

It spanned by the coordinates P∈R ^D and strength t∈R (D + 1) a dimensional space called "bilateral space (bilateral space)". First, ξ (p, t), ξ ₀ (p, t) is generated from x (p), and this is regarded as a point in the bilateral space. From ξ (p, t) and ξ ₀ (p, t), their blended point clouds β (p, t), β ₀ (p, t) (see the following equation) are expressed by equation (7a) Is calculated as the numerator and denominator.

These operations can be regarded as (D + 1) -dimensional Gaussian convolution. This is because [p, t] ^T ∈ R ^{D + 1} is the position of the target pixel in the (D + 1) -dimensional space, and [p, x (q)] ^T ∈ R ^{D + 1} is one of the neighboring positions of the target pixel in the (D + 1) -dimensional space. Because it can be regarded. The BLF output is calculated from these ratios as follows.

  In equation (9a), ζ (p, t) is a BLF processed image in the bilateral space, and equation (9b) is projection (degeneration) from the bilateral space to the original coordinate space. FIG. 1 shows a conceptual diagram of a BLF calculation process when the D-dimensional BLF is regarded as a projection of a (D + 1) -dimensional Gaussian convolution ratio onto the D-dimensional space.

As described above, the D-dimensional BLF is calculated as a ratio of two (D + 1) -dimensional Gaussian convolutions. However, this increase in dimension greatly increases the calculation cost. In the conventional BLF acceleration algorithm described above, β (p, t), β ₀ (p, t), and ζ (p, t) are simplified by various methods to reduce the calculation difficulty. Can be considered.

In the high-speed BLF described in Non-Patent Documents 13, 14, and 15, β (p, t), β ₀ (p, t) are obtained by quantization in a bilateral space and piece-wise linear interpolation. , Ζ (p, t) is actively reduced. Quantization significantly reduces the computational cost, but excessive simplification reduces approximation accuracy. Furthermore, its accuracy is not robust to shifts in the spatial or intensity axis direction (translation-variant problem). This is because the accuracy strongly depends on the quantization threshold.

On the other hand, in Non-Patent Document 1-3, β (p, t) and β ₀ (p, t) are approximated by a shiftable kernel. That is, the Gaussian range kernel was replaced with a translatable kernel such as a triangular function kernel (base) or a polynomial kernel (base) (see equations (3a) and (4)). This idea solves the translation change problem and provides a reasonable way to control the approximation accuracy using the tolerance parameter. However, in order to obtain an accuracy of approximation that is acceptable in actual applications, a considerably large number of iterations are still required in the calculation of the kernel.

Actually, when Gaussian kernel calculation is performed by Chaudhury approximation (formulas (3a) and (4)), the convergence speed changes depending on the magnitude of σ (for example, see Non-Patent Document 2, Table I, Table III). . In particular, as σ becomes smaller, the convergence speed becomes slower and the number of calculation truncation N becomes larger. For example, in the case of the expression (3a), if the range of x is [−T / 2, T / 2], in order to ensure that the approximate expression is positive and monotonic in the range, in principle, N Must be at least greater than N ₀ = 4T ² / π ² σ ² = 0.405 (T / σ) ² (Non-patent Document 2) (In Non-Patent Documents 2 and 3, the Gaussian kernel is | Using the fact that x |> 3σ is almost zero, the condition of N ₀ is further relaxed). In general, T = 255 is used in a gray scale image, and at this time, N ₀ = 413, σ = 5 when σ = 8, N ₀ = 1053 when σ = 5, and N ₀ = 2929 when σ = 3. In the approximation of Expression (4), convergence is faster than Expression (3a), but according to Non-Patent Document 2 and Table III, N ₀ = 88 when σ = 8, N ₀ = 92, σ when σ = 5. When N = 3, N ₀ = 96. When σ is small, N is still large, and the BLF calculation for one pixel requires 100 order iterations. Therefore, there is still a problem that a large calculation cost is required. Since a narrow range kernel is required to enhance edge preservation performance, solving this problem is particularly important when applying BLF to online image processing.

  Since the Deriche approximation (Expression (5)) can also be translated, the translational change problem can be overcome even if Expression (5) is used instead of Expressions (3a) and (4). However, the Deriche approximation is a Gaussian kernel fitting equation, and has a relatively large approximation error and poor accuracy. Further, since the approximation error changes depending on σ, there is a problem that the approximation error cannot be controlled.

In the high-speed BLF described in Non-Patent Document 16, ξ (p, t) and ξ ₀ (p, t) are treated as point clouds, and instead of D-dimensional Gaussian convolution, Then, β (p, t) and β ₀ (p, t) are obtained by performing high-speed Gaussian transformation. However, this algorithm tends to increase the computation time when σ _s becomes small.

  In addition, the high-speed BLF described in Non-Patent Documents 10 and 11 has a problem in that the management cost increases because the histogram is used.

  As described above, the conventional high-speed BLF described above is not sufficient in both calculation cost and accuracy, particularly in application to online video processing and application to higher-dimensional / high-resolution image processing. Therefore, in O (1) BLF, there is a strong demand for further optimizing the right between calculation cost and approximation accuracy. In addition, since the required approximation accuracy varies depending on the use of BLF, it is also important that calculation errors can be controlled freely.

  SUMMARY OF THE INVENTION An object of the present invention is to provide a Gaussian kernel arithmetic device that is less computationally expensive than the prior art and that can freely control calculation errors, and an image filter arithmetic device using the same.

  As apparent from the equations (1a) and (2), the BLF and the cross / joint BLF have only a slight difference in the definition equations, so only the BLF will be described as a representative in this specification. The analysis and approach in the present invention described below can be applied to the cross / joint BLF in exactly the same manner. Hereinafter, [1] to [5] will explain the concept and basic principle of the present invention, and [6] will explain the configuration and operation of the present invention. In the following description, assuming that σ is given, Gaussian kernel g (σ, x) is abbreviated as g (x) unless otherwise specified.

[1] Compression Approach In the present invention, from the viewpoint of image compression, the above-described β (p, t), β ₀ (p, t), and ζ (p, t) are efficiently compressed, thereby calculating cost. (Hereinafter referred to as “compression approach”). FIG. 2 shows β (p, t), β ₀ (p, t), ζ (p) sliced at regular intervals along the t-axis in a standard image (“lenna”) developed into a bilateral space. , T). In FIG. 2, the dynamic range of the grayscale original image is normalized to x (p) ε [0, 1]. It can be seen that the layer images of each series change very smoothly with respect to t, and show high correlation between adjacent layer images. In lossy compression and the like, it is known that such redundant data can be compressed by linear transformation approximation at a slight sacrifice in accuracy.

Therefore, consider irreversible compression of β (p, t) as shown in FIG. First, β (p, t) is transformed by a linear transformation F _t along the t-axis, then some insignificant spectral components are removed in the transformation domain, and finally by an inverse transformation F _t ^−1. Reconstruct β ^ (p, t), which is the compressed β (p, t). β ₀ (p, t) and ζ (p, t) can be compressed in the same manner. The BLF compressed in this manner is hereinafter referred to as “compressive bilateral filtering (CBLF)”. Now referred to as the compression process and _{C t.} Since C _t uses only t as a variable, the compression process C _t results in the following Gaussian kernel compression process.

The Gaussian kernel compressed by C _t is given as the sum of the basis functions of the linear transformation F _t as follows.

Here, N is the number of significant spectral components, and w _n (•) and φ _n (•) are the n-th transform coefficient and the basis function of F _t , respectively. Substituting equation (11) into equation (10) yields:

Here, Ω _n (p) is the nth component image.

Ω _n (p) is interpreted as a Gaussian convolution of the product image φ _n (x (q)) x (q). Thus, this compression approach decomposes the BLF into N Gaussian convolutions. N corresponds to the number of iterations of the convolution operation. If N is decreased, the calculation cost is reduced, but the approximation accuracy is degraded.

[2] Compression of Gaussian kernel by Fourier basis It is known that Karhunen-Loeve transform (KLT) may be used as linear transform F _t in order to minimize the l ₂ error in image compression. In KLT, an optimum basis function φ _n (x (q)) is obtained from a set of g _r (x (q) −t) shifted with respect to various t. However, in KLT, generally derived basis functions φ _n (x (q)) and transformation coefficients w _n (t) are mathematically complicated, and the eigenvalue problem is solved numerically, resulting in excessive calculation costs. There is a problem that is necessary.

A suitable alternative to KLT is the family of Fourier transforms. In the field of image compression, it is widely known that KLT can be approximated by discrete cosine transform when the target is a stationary first-order Markov signal with high correlation (Non-patent Document 17). Therefore, in the present invention, it is assumed that the Gaussian kernel set {g _r (x (q) −t)} shifted with respect to t is a stationary first-order Markov signal having high correlation. (See FIG. 2), the Gaussian range kernel is compressed using a Fourier series (FS).

In the following, consider a Gaussian range kernel compressed by FS. When it is clear that the kernel is a Gaussian range kernel, the subscript r of g _r (•) is omitted and simply denoted as g (•).

  The Gaussian kernel g (x) in the domain x∈ [−R, R] can be compressed as follows using a Fourier basis.

Here, K∈N (N is a set of all natural numbers) is the number of low-frequency components used for approximation (hereinafter referred to as “truncation term number”). a _k is a Fourier coefficient. T / 2 represents a range where g (x) takes a significant value, and T is hereinafter referred to as a “significant period”. Since the discrete Gaussian kernel g (x) is an even function, only the cos term needs to be considered. Therefore, Formula (11) becomes as follows.

  From the above equation, this compression approach requires N = 2K + 2 iterations and one sum operation for Gaussian convolution. Therefore, it can be seen that K can be made as small as possible in order to speed up the BLF calculation.

[3] Optimization of Periodic Gaussian Kernel By approximation of equation (13a), Gaussian kernel g (x) in domain x∈ [−T / 2, T / 2] becomes domain x∈ (−∞ , ∞) is approximated by a periodic function g ^ (x; K, T). Hereinafter, g ^ (x; K, T) is referred to as a “periodic Gaussian kernel”. In order to clarify the compression rate of the compression approach by Fourier transform, first, _ak is approximated in the following closed form.

Here, g (x) ≈0 when | x |> T / 2.

Since a _k decreases exponentially, it can be approximated with sufficient accuracy from equation (13a) even if K is reduced to some extent. If σ is given, the attenuation rate of a _k with respect to k is determined. It can be seen that T is the only control factor. The smaller the T, the faster the attenuation of _ak .

  FIG. 4 shows a periodic Gaussian kernel in three different periodic conditions. In FIG. 4, the dynamic range of the input image is [0, R]. In general, R = 255 is widely used in grayscale images, and R = 1 is widely used in binary images. The domain that the range kernel must support for equation (11) is [−R, R]. The double arrows in the center of each figure in FIG. 4 represent the center period of each periodic Gaussian kernel. When the significant period T of the periodic Gaussian kernel is T / 2 = R, the periodic Gaussian kernel has an accurate Gaussian shape throughout the domain. Since Gaussian has a long tail with a value of almost zero, even when the significant period T is T / 2 = 0.7R, the periodic Gaussian kernel still maintains a Gaussian shape throughout the domain. However, when the significant period T is shortened to T / 2 = 0.5R, a Gaussian profile having a period adjacent to the right and left of the domain [−R, R] enters the domain [−R, R]. Therefore, the periodic Gaussian kernel is greatly collapsed from the Gaussian shape. This suggests that the optimum significant period T exists from the viewpoint of approximation error when the number of truncation terms K is given in the equation (13a).

  In order to obtain the optimal combination (K, T) of the number of truncation terms K and the significant period T, first, an approximation error E (K, T), which is an evaluation function, is defined by a “kernel error” of the following equation: .

Given an allowable error τ (0 <τ << 1), the set (K, T) is constrained by the constraint condition E (K, T) ≦ τ ² . Since E (K, T) is a monotonic function with respect to the K axis, the minimum K that satisfies the constraint condition E (K, T) ≦ τ ² is the T function K (T). That is, the periodic Gaussian kernel optimization problem can be formulated as follows.

Here, N represents a set of natural numbers, and (K _τ , T _τ ) represents the optimum value of (K, T).

  Hereinafter, the equations (17a) and (17b) are referred to as “periodic Gaussian kernel optimization”. Equations (17a) and (17b) provide a method for rationally determining (K, T) while freely controlling the approximation error.

  On the other hand, the optimization of the periodic Gaussian kernel is also a new method for O (1) BLF. In the conventional approach using a translatable kernel (Non-Patent Documents 1-3), attention is paid only to the case where T = R, whereas in the approach using the periodic Gaussian kernel of the present invention, FIG. This is because it is generalized even when T ≦ R. This is a new focus in the present invention. In the conventional translatable kernel approach, the effective dynamic range of the target image is actually measured in advance by pre-scanning called MAXFILTER to reduce the dynamic range R itself. However, this method is only effective for low-contrast images and depends on the content of the target image. Due to this content dependency, there is a problem that the worst case performance cannot be improved. On the other hand, in the approach based on the optimization of the periodic Gaussian kernel in the present invention, the worst case performance can be surely improved by fully utilizing the geometric characteristics of the Gaussian kernel.

[4] Evaluation of Approximation Error Next, a specific method for solving the periodic Gaussian kernel optimization problem of equations (17a) and (17b) will be described. When T = R, equation (17a) can be solved analytically using the orthogonality of equation (13a), but not when T <R. Further, if the equation (17a) is solved by an iterative method, there is a problem from the viewpoint of speeding up because a large calculation cost is required for the evaluation of the equation (16). Therefore, in the present invention, the value of Expression (16) is estimated by the following closed form approximate expression under some appropriate assumptions.

Here, erfc (·) is a complementary error function.

  E ^ (K, T) in equation (18) is called "estimated kernel error". In the expression (18), the right side of the expression (16) is divided into the section [−T / 2, T / 2] (the first term on the right side of the expression (18)) and the other part (the second right side of the expression (18)). And approximated under the above assumptions. The detailed proof of the equation (18) is omitted because it is not directly related to the implementation of the present invention.

  In the present invention, the approximate error amount is evaluated using the equation (18) instead of the equation (16) (that is, E (K, T) ≈E ^ (K, T)).

In Expression (18), each term on the right side is non-negative. Accordingly, the necessary conditions for satisfying the constraint condition E (K, T) ≦ τ ² are erfc (πσ (2K + 1) / T) ∈ (0, τ ² ], erfc ((TR) / σ) ∈ [ 0, τ ² ] Therefore, the upper and lower limits of T are determined as follows.

From equation (19a), σζ _τ + R ≦ (πσ (2K + 1)) / ζ _τ , and considering that K is an integer, the minimum satisfying the constraint condition E ² (K, T) ≦ τ ² K is determined as follows.

Since the significant period T is within the range of the equation (19a), T _τ that minimizes the energy error E ^ (K, T) within this range is Tτ. That is,

  From the equation (18), the optimization problem of the objective function E ^ (K, T) of the equation (21) results in the problem of solving the following equation.

The solution _Tτ of the above equation (22) can be numerically solved by the bisection method.

FIG. 5 shows a comparison of the calculated values of the kernel error E (K, T) and the approximate kernel error E ^ (K, T) for several K when R = 1 and σ = 0.1. In FIG. 5, the horizontal axis indicates the period T, and the vertical axis indicates the kernel error E (K, T) or the approximate kernel error E ^ (K, T). FIG. 5 shows that when K is sufficiently large, E (K, T) ≈E ^ (K, T) is satisfied. Under the condition of R = 1 and σ = 0.1, when the allowable error is τ ² = 0.005, the minimum K is a constant value K = 4.

[5] Simplified periodic Gaussian kernel sub-optimization The optimal (K, T) can be determined by using the technique described in [4], If sub-optimal (K, T) is obtained at the expense of calculation cost, a more simplified method described below can be adopted.

From FIG. 4, the periodic Gaussian kernel g ^ (x; K, T) is the Gaussian kernel g (x) in the closed interval [−T / 2, T / 2] when K is sufficiently large. It can be seen that this is expressed as the sum of peak functions g _n (x) translated by 2 nT (n is an integer) in the x-axis direction. Therefore, g ^ (x; K, T) can be expressed as follows using the Heaviside step function H ₀ (x).

Where δ (ξ) is the Dirac delta function.

When K is sufficiently large, when T = R, g ^ (x; K, T) matches g (x) in the domain [−R, R]. When T <R, it is affected by the peak function g _{± 1} (x) of the adjacent period of g ₀ (x) in the vicinity of both ends (x = ± R) of the domain [−R, R]. Error. Therefore, T may be determined so as to compensate for g _{± 1} (± R) ≈0 at x = ± R as follows.

If g (x) ≈0 is compensated for x ≧ ασ, in order to compensate so as not to be affected by the peak function g _{± 1} (x) of the adjacent period, TR = ασ may be used. Therefore, the optimum value _Tτ of the significant period T can be easily determined as follows.

Here, α is a parameter representing a significance level (a ratio at which a function value is considered significant). For example, if the significance level is 99.73%, α = 3.

If it is determined the optimum value T _tau significant period T by the equation (25), the above equation (17b), it is possible to determine the optimum value K _tau of truncation section number K from (19a). For example, the optimum value _Kτ can be calculated by the following equation.

Since this method does not optimize T and K at the same time, generally, the calculated K _τ tends to be larger than the method described in [4], and the calculated (K _τ , T _τ ) is Suboptimal value. However, it is possible to determine the T _tau easily compared to [4], which is one of effective methods.

[6] Configuration of the Present Invention The first configuration of the image filter arithmetic device according to the present invention is based on the input image {x (p)} (p is the pixel position and x (p) is the pixel value at position p). The spatial deviation g _s (qp) defined by the Gaussian function of standard deviation (hereinafter referred to as “scale parameter”) σ _s (p and q are pixel positions), and the Gaussian function of scale parameter σ _r An image filter arithmetic device for smoothing an image by a convolution operation with a defined range kernel g _r (x (q) -x (p)),
For the scale parameter σ _r , the predetermined tolerance τ and the dynamic range R of the input image {x (p)}, the range kernel g _r (x (q) −x (p)) is a finite term. The series truncation term K in the approximate range kernel g ^ _r (x (q) -x (p)) approximated by the Fourier series of the number and the approximate range kernel g ^ _r (x (q) -x ( p)) an approximation that is a set of significant periods T that define an integration region [-T / 2, T / 2] of Fourier integration when calculating the expansion coefficient a _k (k = 0,..., K) of each term in p)) Approximate parameter optimizing means for determining and outputting parameters (T, K) such that the approximate error of the approximate range kernel g ^ _r (x (q) -x (p)) is equal to or less than the allowable error τ; ,
Expansion coefficients from the 0th term to the Kth term of the approximate range kernel g ^ _r (x (q) -x (p)) according to the approximate parameter (T, K) calculated by the approximate parameter optimization means expansion coefficient calculating means for calculating {a ₀ ,..., a _K };
Using the expansion coefficients {a ₀ ,..., A _K } calculated by the expansion coefficient calculating means, the spatial kernel g _s (qp) and the approximate range kernel g are applied to the input image {x (p)}. ^ _r (x (q) -x (p)) and a filter operation means for executing a convolution operation.

  According to this configuration, the approximation parameter optimization means optimizes by changing not only the number of truncation terms K but also the significant period T in accordance with the allowable error τ, so that the calculation cost is lower than conventional, In addition, calculation error can be freely controlled.

According to a second configuration of the image filter arithmetic apparatus according to the present invention, in the first configuration, the approximate parameter optimization means includes a truncation term number calculating means for calculating the truncation term number K,
The truncation term number calculating means calculates the truncation term number K with respect to the scale parameter σ _r , a predetermined allowable error τ, and the dynamic range R of the input image {x (p)}.

It is characterized by calculating by.

According to a third configuration of the image filter arithmetic apparatus according to the present invention, in the first or second configuration, the approximate parameter optimization unit includes a significant cycle optimization unit that calculates the significant cycle T,
The expansion coefficient calculating means is configured to calculate an approximate expansion coefficient obtained by performing Fourier integration of the range kernel g _r (x (q) −x (p)) in an integration region [−T / 2, T / 2], and each expansion coefficient. {a ₀ , ..., a _K }
The significant period optimization means includes an equation of a variable T with respect to the scale parameter σ _r , the truncation term number K, and the dynamic range R.

T is calculated as the significant period T.

According to a fourth configuration of the image filter arithmetic apparatus according to the present invention, in the first configuration, the approximate parameter optimization unit includes:
Based on the scale parameter σ _{r, the} range upper limit value R, and a preset effective width parameter α (ασ _r <R), the significant period T is

A calculation cycle setting means calculated by:
Based on the allowable error τ, the significant period T, and the scale parameter σ _r , the truncation number K is

And a truncation term number determining means for calculating by the following.

In a fifth configuration of the image filter arithmetic device according to the present invention, in any one of the first to fourth configurations, the expansion coefficient calculation means includes the expansion coefficients {a ₁ ,..., A _K }. The

It is a thing calculated by these.

  An image filter calculation program according to the present invention is read and executed by a computer, thereby causing the computer to function as the image filter calculation device having any one of the first to fifth configurations.

The first configuration of the Gaussian kernel arithmetic device according to the present invention calculates an approximate value of a Gaussian kernel g _r (x) defined by a Gaussian function of standard deviation (hereinafter referred to as “scale parameter”) σ _r . A Gaussian kernel computing device that
For the scale parameter σ _r , the predetermined tolerance τ and the dynamic range R of the input image {x (p)}, the range kernel g _r (x (q) −x (p)) is a finite term. The series truncation term K in the approximate range kernel g ^ _r (x (q) -x (p)) approximated by the Fourier series of the number and the approximate range kernel g ^ _r (x (q) -x ( p)) an approximation that is a set of significant periods T that define an integration region [-T / 2, T / 2] of Fourier integration when calculating the expansion coefficient a _k (k = 0,..., K) of each term in p)) Approximate parameter optimizing means for determining and outputting parameters (T, K) such that the approximate error of the approximate range kernel g ^ _r (x (q) -x (p)) is equal to or less than the allowable error τ; ,
The expansion coefficients {a ₀ ,..., A _K } from the 0th term to the Kth term of the approximate Gaussian kernel g _r (x) according to the approximate parameters (T, K) calculated by the approximate parameter optimization means. Expansion coefficient calculating means for calculating
Approximate function calculation means for calculating the approximate Gaussian kernel g _r (x) using the expansion coefficients {a ₀ ,..., A _K } calculated by the expansion coefficient calculation means. .

According to a second configuration of the Gaussian kernel arithmetic device according to the present invention, in the first configuration, the approximate parameter optimization unit includes:
A truncation term number calculating means for calculating the truncation term number K;
The truncation term number calculating means calculates the truncation term number K with respect to the scale parameter σ _r , a predetermined allowable error τ, and the dynamic range R of the input image {x (p)}.

It is characterized by calculating by.

According to a third configuration of the Gaussian kernel arithmetic device according to the present invention, in the first or second configuration, the approximate parameter optimization unit includes a significant cycle optimization unit that calculates the significant cycle T,
The expansion coefficient computing means calculates an approximate expansion coefficient obtained by performing Fourier integration of the approximate Gaussian kernel g _r (x) in the integration region [−T / 2, T / 2], and the expansion coefficients {a ₀ ,. _K }, and
The significant period optimization means includes an equation of a variable T with respect to the scale parameter σ _r , the truncation term number K, and the dynamic range R.

T is calculated as the significant period T.

According to a fourth configuration of the Gaussian kernel arithmetic apparatus according to the present invention, in the first configuration, the approximate parameter optimizing means includes:
Based on the scale parameter σ _{r, the} range upper limit value R, and a preset effective width parameter α (ασ _r <R), the significant period T is

A calculation cycle setting means calculated by:
Based on the allowable error τ, the significant period T, and the scale parameter σ _r , the truncation number K is

And a truncation term number determining means for calculating by the following.

In a fifth configuration of the Gaussian kernel calculation device according to the present invention, in any one of the first to fourth configurations, the expansion coefficient calculation means includes the expansion coefficients {a ₁ ,..., A _K }

It is a thing calculated by these.

  A Gaussian kernel arithmetic program according to the present invention is read and executed by a computer, thereby causing the computer to function as a Gaussian kernel arithmetic device having any one of the first to fifth configurations.

  As described above, according to the image filter arithmetic device and the Gaussian kernel arithmetic device according to the present invention, the calculation cost is lower than the conventional one, and the calculation error can be freely controlled.

  Although the image filter arithmetic device and Gaussian kernel arithmetic device of the present invention are mainly premised on application to BLF, the application range is not limited to BLF, and is applicable to the above-described cross / joint BLF and the like.

It is a conceptual diagram of a BLF calculation process in the case where the D-dimensional BLF is regarded as a projection of a (D + 1) -dimensional Gaussian convolution ratio to the D-dimensional space. In the bilateral space, layer images of β (p, t), β ₀ (p, t), and ζ (p, t) sliced at regular intervals along the t-axis. It is a figure which shows the structure procedure of the compression bilateral filter by the linear transformation approach along t-axis. FIG. 4 shows a periodic Gaussian kernel in three different periodic conditions. It is a comparison figure of the calculated value of kernel error E (K, T) and rough kernel error E ^ (K, T) for some K in the case of R = 1 and σ = 0.1. FIG. 3 is a diagram illustrating a hardware configuration of a Gaussian kernel arithmetic device according to the first embodiment. 3 is a block diagram illustrating a functional configuration of a Gaussian kernel arithmetic device according to Embodiment 1. FIG. 5 is a PAD representing the overall flow of the calculation processing of the optimum approximate parameters (K _τ , T _τ ) and the significant coefficient a _k in the Gaussian kernel arithmetic device 10 of the first embodiment. In the Gaussian kernel arithmetic unit 10 of the first embodiment, a PAD indicating the optimum value calculation processing significant period T _tau. FIG. 6 is a diagram illustrating a relationship between a Gaussian kernel truncation term K approximated by series expansion and an approximation error E; It is a figure which shows the relationship between (sigma) _r and (tau). FIG. 6 is a block diagram illustrating a functional configuration of a filter arithmetic device according to a second embodiment. 10 is a PAD showing a flow of bilateral calculation processing in the filter calculation device 20 of the second embodiment. 5 is a PAD representing a two-dimensional discrete Gaussian filter calculation process by the GF calculator 29. It is the figure which compared the approximation accuracy [dB] and calculation time in the Chaudhury approximation algorithm ((epsilon) = 0.8) and CBLF ((tau) = 0.1) of a nonpatent literature 2. FIG. It is the figure which compared the calculation accuracy with the approximation accuracy [dB] in Yang approximation algorithm (B = 8) of nonpatent literature 14, and CBLF. An original image, a smoothed image generated by the original BLF, Chaudhury approximation algorithm, Yang approximation algorithm, and CBLF algorithm, and an error image amplified by 10 times. It is a block diagram which shows the whole structure of the filter calculating apparatus which concerns on Example 3. FIG. It is a block diagram showing the structure of the calculation parameter memory 51 and the approximate parameter optimization part 52 of FIG. FIG. 19 is a block diagram illustrating a configuration of a filter operation circuit 54 in FIG. 18. FIG. 21 is a block diagram illustrating a configuration of a B circuit 78 (B circuits 78-1 to 78-K) in FIG. 20.

  Hereinafter, embodiments for carrying out the present invention will be described with reference to the drawings.

  FIG. 6 is a diagram showing a hardware configuration of the Gaussian kernel arithmetic device according to the present invention.

  In FIG. 6, the Gaussian kernel arithmetic device of the present invention includes a computer 1, an input device 2, an output device 3, and an external storage device 4. The input device 2 is a general computer input device such as a keyboard, a mouse, or a CD-ROM. The output device 3 is a general computer output device such as a printer or a display. The external storage device 4 is a general external storage device such as a magnetic disk. The Gaussian kernel arithmetic device according to the present embodiment is realized by reading the provided Gaussian kernel arithmetic program from the input device 2 into the computer 1 and executing it.

FIG. 7 is a block diagram showing a functional configuration of the Gaussian kernel arithmetic device 10 according to the present invention. The functional configuration of FIG. 7 is realized by causing the computer 1 to execute the Gaussian kernel arithmetic program.

  The operation parameter memory 11 stores operation parameters σ, τ, R, etc. of Gaussian kernels input in advance from the input device 2. Here, σ is a Gaussian kernel scale parameter, τ is an allowable approximation error, and R is a maximum value (range upper limit value) that can be taken by the sample value.

The approximate parameter optimization means 12 determines the number K of truncation terms and the significant period T of the periodic Gaussian kernel シ ア ン (x; K, T) described above. The approximate parameter optimization means 12 includes a truncation number calculator 13, a significant period optimizer 14, and a KT register 14a. The truncation term number calculator 13 calculates the truncation term number K _τ that is the minimum with respect to the allowable error τ, using Equation (20). The significant period optimizer 14 calculates the value of the significant period T _τ that minimizes the approximation error of Expression (18) with respect to the number of truncation terms K _τ . The KT register 14 a holds the optimized approximate parameters (K _τ , T _τ ) calculated by the truncation term number calculator 13 and the significant period optimizer 14.

The significant coefficient calculator 15 uses the truncation term number K and the optimum value (K _τ , T _τ ) of the significant period T calculated by the approximate parameter optimization unit 12 and held in the KT register 14a according to the equation (13b). A Fourier coefficient (hereinafter referred to as “significant coefficient”) a _k (k = 0, 1,..., K) of the periodic Gaussian kernel g ^ (x; K _τ , T _τ ) is calculated. The significant coefficient memory 16 stores the significant coefficient a _k (k = 0, 1,..., K _τ ) calculated by the significant coefficient calculator 15.

The approximate function calculator 17 uses the sample value x input from the input device 2 and the significant coefficient a _k (k = 0, 1,..., K) stored in the significant coefficient memory 16 according to Expression (13a). The value of the periodic Gaussian kernel g ^ (x; K, T) is calculated and output to the output device 3.

The operation of the Gaussian kernel arithmetic apparatus 10 according to this embodiment configured as described above will be described below. In the Gaussian kernel arithmetic device 10 according to the present embodiment, when arithmetic parameters (R, σ, τ) are first input from the input device 2 and these are set in the arithmetic parameter memory 11, first, approximate parameter optimization is performed. The means 12 determines the number K of truncation terms and the optimum value (K _τ , T _τ ) of the significant period T, and then calculates the significant coefficient a _k (k = 0, 1,..., K) by the significant coefficient calculator 15. These significant coefficients _ak are stored in the significant coefficient memory 16. Thereafter, each time the sample value x is input from the input device 2, the approximate function calculator 17 uses the censored term number K and the significant coefficient _ak calculated in advance, and the periodic Gaussian kernel g ^ (x; The value of K _τ , T _τ ) is calculated by equation (13a), and this is output to the output device 3.

[1] Calculation of the number of truncation terms K and the significant coefficient a _k FIGS. 8 and 9 show the calculation processing of the number K of truncation terms, the significant period T, and the significant coefficient a _k in the Gaussian kernel computation device 10 of the first embodiment. It is a flowchart. Here, it is assumed that calculation parameters (R, σ, τ) are input from the input device 2 in advance and are set in the calculation parameter memory 11.

[1.1] Overall Flow FIG. 8 is a PAD (Problem Analysis Diagram) (JIS X 0128: 1988) showing the overall flow of the calculation processing of the optimum approximate parameters (K _τ , T _τ ) and the significant coefficient a _k. / ISO 8631). Here, first, (K _τ , T _τ ) is set so that the approximation error E ^ in the equation (18) becomes the allowable error τ and K is minimized, and then the significant coefficient (Fourier coefficient) a _k is calculated. I do.

(Step S101)
First, the truncation term number calculator 13 reads each parameter (R, σ, τ) set in the computation parameter memory 11, and calculates the optimum number of truncation terms K _τ for the allowable error τ by Expression (20). To the KT register 14a.

(Step S102)
Next, the significant cycle optimizer 14 reads each parameter (R, σ, τ) set in the calculation parameter memory 11 and solves the equation (22) of T in the closed interval of the equation (19a). calculates the optimal significant period T _tau for tolerance tau, saving to KT register 14a. Since equation (22) cannot be solved analytically, it must be solved numerically by an iterative method. However, since the range that T can take is restricted to a narrow range by the equation (19a), it is preferable to calculate by using the bisection method. A specific example of the calculation process of the significant period _Tτ will be described later.

(Step S103)
Next, the significant coefficient calculator 15 calculates the zeroth Fourier coefficient a ₀ (= 1 / T _τ ) and stores it in the significant coefficient memory 16.

(Steps S103-S104)
Next, the significant coefficient calculator 15 calculates the first to K _τ- th Fourier coefficients a _k (k = 1,..., K) from the optimum values (K _τ , T _τ ) previously held in the KT register 14a. The approximate value of _{τ 2} ) is calculated by equation (15) and stored in the significant coefficient memory 16.

Through the above steps S101 to S104, the significant coefficient memory 16 is in a state where the Fourier coefficient a _k and the number of truncation terms K _τ necessary for the calculation of the Gaussian kernel are stored.

[1.2] calculation of significant period T _tau Next, an example of calculation of significant period T _tau in step S102 in FIG. 8. Here, since the search range of T is restricted to a narrow range by the equation (19a), an example of searching for the optimum value _Tτ using a dichotomy algorithm is shown. Figure 9 is a PAD indicating the calculation of significant period T _tau in step S102.

(Step S201)
First, the significant period optimizer 14 sets the initial value [T ₁ , T ₂ ] of the search range of T, the maximum number of iterations i _max , and the accuracy ε _x of the root of equation (22). The initial value [T ₁ , T ₂ ] of the search range is calculated by the equation (19a). The maximum number of iterations i _max may be set as appropriate in consideration of convergence, but it is sufficient to set iter = 12.

(Step S202)
Next, the significant period optimizer 14 performs an initial evaluation of the objective function f (T; K _τ , R, σ) at T = T ₁ and T = T ₂ , and z ₁ = f (T ₁ ; K _τ , R, σ), z _mid = f (T ₂ ; K _τ , R, σ). Here, the objective function f (T) is the left side of the equation (22).

(Steps S203 to S205)
Next, the significant cycle optimizer 14 sets a search direction. That is, when z ₁ <0 (S203), the search direction dT is set to the positive direction (T ₂ −T ₁ ), and the search reference point T _rtb is set to T ₁ (S204). When z ₁ ≧ 0 (S203), the search direction dT is set to the negative direction (T ₁ −T ₂ ), and the search reference point T _rtb is set to T ₂ (S205).

(Steps S206 to S211)
Next, the significant period optimizer 14 calculates a solution of the equation (22) by binary search. That is, the processes of S207 to S211 are repeatedly executed.

First, dT is halved, an intermediate point T _mid of the current search range is calculated by T _rtb + dT, and an objective function value f (T _mid ; K _τ , R, σ) at the intermediate point T _mid is calculated. This is set as _zmid (S207). Next, when z _mid is 0 or less, the search reference point T _rtb is shifted to the intermediate point T _mid (S208, S209). Then, | dT | is when it becomes an objective function value z _mid are 0 or intermediate point T _mid reached following accuracy epsilon _x outputs the search reference point T _rtb as the solution (the optimum value T _tau) (S210, S211), otherwise return to step S207.

By performing the above processing, the optimum significant period _Tτ that is the root of the equation (22) can be easily calculated.

In the present embodiment, the example in which the approximate parameter optimization unit 12 determines (K _τ , T _τ ) according to the equations (18) and (19a) has been described. However, as a more simplified method, Then, the approximate parameter optimization means 12 calculates T _{τ according} to the equation (22), and using this T _τ , the truncation term number calculator 13 calculates the equation (19a) (specifically, the K (T) calculation in FIG. 8). It may be configured to calculate the optimal truncation claim number K _tau by treatment).

Finally, the calculation cost of the Gaussian kernel arithmetic device 10 according to the present embodiment will be described in comparison with a conventional Gaussian kernel arithmetic device and the like. FIG. 10 is a diagram showing the relationship between the Gaussian kernel truncation number K approximated by series expansion and the approximation error E. In FIG. 10, (a), (b), and (c) show cases where the scale parameter σ is 0.05, 0.1, and 0.2, respectively. The dynamic range R is normalized to 1. The horizontal axis represents the number of truncation terms of spectral components (the number of spectral components considered significant) N, and the vertical axis represents the square root √E of approximation error (Expression (16)). In (a) to (c), “Chaudhury'13” represents the approximation error of the approximate Gaussian kernel calculated by the Chaudhury approximation formula (4) mentioned in the background art (in this case, the horizontal axis is truncated). The number of terms K is the number of terms N of the sum of the right side of the equation (4)). “KLT” represents a case where the Gaussian kernel is approximated by the Karoonen-Leve transform (KLT), and is the approximation that is theoretically most optimal. “Ours (T = 2R)” represents an approximation error in the case where the significant period T is fixed to the dynamic range R in the Gaussian kernel arithmetic device of the first embodiment (that is, normal Fourier series approximation). “Ours (T ≦ 2R)” is determined by the processing according to the flowchart of FIG. 9 in the Gaussian kernel arithmetic unit of the first embodiment, and the optimal value (K _τ , T _τ ) of the truncation term K and the significant period T is determined. Represents the approximation error in the case. In each of the drawings (a) to (c), the line E = 0.05 (5%) is the allowable error τ in this calculation example.

Comparing the value of K at E = 0.5, it is clear that the Gaussian kernel arithmetic unit according to the present embodiment can obtain a desired approximation accuracy with a smaller number of truncation terms K than the conventional Chaudhury approximation. I understand. The tendency becomes more prominent as the scale parameter σ becomes smaller. In addition, in the Chaudhury approximation and the approximation of the Gaussian kernel arithmetic unit of this embodiment, the approximation error E asymptotically approaches a certain value when K increases, and from a certain constant value (hereinafter referred to as “asymptotic approximation error value E _∞ ”). Does not get smaller. asymptotic approximation error value E _∞ as σ is large tends to increase. However, as is clear from FIG. 10 (c), the in approximation of Gaussian kernel computation device of this embodiment, it is possible to reduce the asymptotic approximation error value E _∞ than Chaudhury approximation, good accuracy over a large σ An approximate value can be obtained.

  In addition, comparing “Ours (T = 2R)” and “Ours (T ≦ 2R)”, the significant period T is optimized to a value smaller than the dynamic range R in the Gaussian kernel arithmetic unit of this embodiment. Thus, it can be seen that the desired approximate accuracy can be obtained with a small number of truncation terms K. This is because the shape of the periodic Gaussian kernel approaches the shape of the cosine function by appropriately reducing the significant period T from R with respect to the spread of the Gaussian kernel defined by the scale parameter σ. This is because the high frequency component is reduced. That is, when the dynamic range R is very large with respect to the scale parameter σ, the shape of the periodic Gaussian kernel is a series of sharp peaks with a constant period. Although necessary, the significant period T is made smaller than the dynamic range R, the interval between the peaks is narrowed, and the shape of the periodic Gaussian kernel is brought close to the shape of a smooth wave with a constant period, so that high frequency components necessary for approximation are required. Can be reduced. Therefore, in the Gaussian kernel arithmetic device of the present embodiment, since the truncation term number K is small, the calculation amount (the number of iterations) is reduced as compared with the conventional case when the approximation function calculator 17 performs the calculation of Expression (13a). be able to.

Furthermore, FIG. 10 shows that the smaller the range scale parameter σ _r , the greater the number N (ie, K) of spectral components required. This is because, as σ _r becomes smaller, the glue point groups β (p, t) and β ₀ (p, t) with smaller redundancy are generated from the viewpoint of compressibility. In order to observe this behavior in more detail, FIG. 11 shows the relationship between _σr and τ. Since the approximation effect of the Gaussian range kernel depends only on σ _r / R and τ, the number of spatial filters M = 4K + 1 is known from the contour map of FIG.

[1] Arithmetic Algorithm In this embodiment, a filter arithmetic device in which the Gaussian kernel arithmetic device according to the present invention is applied to a bilateral filter will be described. From equation (1a), the bilateral filter is expressed as follows.

Here, x (p) is an input image, z (p) is an output image, and β (p) and β ₀ (p) are a numerator image and a denominator image, respectively.

By substituting the periodic Gaussian kernel (formula (13a)) according to the present invention into the formulas (33b) and (33c), the numerator image β (p) and the denominator image β ₀ (p) are respectively It is approximated as follows.

z ^ (p) is referred to as a “compressive bilateral filter (CBLF)”, β ^ (p) is referred to as an “approximate numerator image”, and β ^ ₀ (p) is referred to as an “approximate denominator image”.

In the approximate numerator image β ^ (p) and the approximate denominator image β ^ ₀ (p), the sum of q (portion enclosed by {•}) is the function cos (πkx (p) / L) x ( p), sin (πkx (p) / L) x (p), cos (πkx (p) / L), and sin (πkx (p) / L) are Gaussian filter processed values. Therefore, when the Gaussian filter function GF (σ, p, {f}) (σ is a scale factor, p is a pixel position vector, and {f} is a set of function values of the scalar function f), an approximate molecular image β ^ (P), the approximate denominator image β ^ ₀ (p) can be expressed as follows.

  Regarding the calculation of the Gaussian filter function GF (σ, p, {f}), since many O (1) Gaussian filter algorithms are already known (see, for example, Non-Patent Documents 18-23), Can be used. As a result, the bilateral filter operation algorithm using the Gaussian kernel operation method according to the present invention is as follows.

  Here, cos (•) and sin (•) are defined as vector operations in element units. Also,

Represents elemental summation, integration, and division. The function Calc_Significant_Period (·) represents solving Equation (22) by an appropriate solution. Specifically, the process of FIG. 9 of the first embodiment can be used. The function GF (•) represents the O (1) algorithm of the spatial Gaussian filter. Here, the efficient O (1) algorithm described in Non-Patent Documents 18 and 19 is used. This algorithm operates approximately twice as fast as a recursive Gaussian filter (Non-Patent Documents 20-23) widely known as a conventional O (1) Gaussian filter, and is stable over a wide range of σ _s. Works. The calculation cost of the CBLF in this embodiment is mainly the number of (4K + 1) spatial Gaussian filters ((2K + 1) β (p, t) and 2K β ₀ (p, t)). Dependent.

[2] Device Configuration and Operation In this embodiment, the overall hardware configuration of the filter operation device is the same as that shown in FIG. 6, and the filter operation device according to this embodiment uses a provided filter operation program. It is assumed that it is realized by being read from the input device 2 to the computer 1 and executed.

In the following description, for the sake of simplicity, unless otherwise specified, the approximate molecular image β ^ (p) and the approximate denominator image β ^ ₀ (p) are omitted, and the molecular image β (p) is simply omitted. , Denominator image β ₀ (p).

FIG. 12 is a block diagram illustrating a functional configuration of the filter arithmetic apparatus according to the second embodiment. In FIG. 12, the same parts as those in FIG. 7 of the first embodiment are denoted by the same reference numerals. The filter operation device according to the present embodiment may be configured as a hardware as an embedded board using a microcomputer or the like, but the provided filter operation program is read from the input device 2 to the computer 1 and executed. This may be realized by doing so.

  In addition to the calculation parameter memory 11, the approximate parameter optimization unit 12, the significant coefficient calculator 15, the significant coefficient memory 16, and the approximate function calculator 17 described in the first embodiment, the filter calculation device 20 according to the second embodiment includes: An input image memory 21, an output image memory 22, and a term counter 23 are provided. The approximate parameter optimizing means 12, the significant coefficient calculator 15, and the significant coefficient memory 16 are the same as those in the first embodiment, and thus description thereof is omitted.

The operation parameter memory 11 stores operation parameters σ _r , σ _s , τ, R of Gaussian kernels input in advance from the input device 2. Here, σ _r and σ _s are spatial and range dimension scale parameters (standard deviation), τ is an allowable approximation error, and R is a maximum value (range upper limit value) that a pixel value can take.

The input image memory 21 stores an input image x = (x _ij ) input from the input device 2. The input image x is image data having a size of N × M pixels, and x _ij represents the pixel value of the pixel at the position (i, j). The pixel value of each pixel takes a value from 0 to R.

The output image memory 22 stores an output image z = (z _ij ) obtained by performing bilateral arithmetic processing on the input image x.

The term number counter 23 counts the Fourier series term number k for which approximate calculation is performed. As in the first embodiment, the approximate parameter optimization unit 12 calculates and outputs the optimum number of truncation terms K _τ, so the term number counter 23 counts the term number k from 0 to K _τ .

The approximate function calculator 17 is calculated by the approximate parameter optimization means 12 (K _τ , T _τ ), and the significant coefficient a _k (k = 0, 1,..., K _τ ) calculated by the significant coefficient calculator 15. ) To perform bilateral filter function operations. Here, as described above, the numerator image β (x) and the denominator image β ₀ (x) of equations (31a) and (31b) are calculated separately, and finally the output image z is calculated by equation (30c). To do.

  The approximate function calculator 17 includes a basis function calculator 25, a basis function table 26, a denominator image calculator 27, a numerator image calculator 28, a GF calculator 29, a denominator image memory 30, an adder 30a, a numerator image memory 31, and an addition. A device 31a and an image divider 32 are provided.

  The basis function calculator 25 calculates the basis function of the Fourier expansion equation shown in equations (34a) and (34b) for each pixel p (i, j) of the input image x at the term number k output from the term counter 23.

Perform the operation.

  The basis function table 26 stores the value of each basis function calculated by the basis function calculator 25 as a lookup table (LUT).

The denominator image calculator 27 calculates the pixel value β ^(k) ₀ (p) = β ^(k) _0ij of the k-th term of the denominator image β ₀ with respect to each pixel p by Expression (35b). The molecular image calculator 28 calculates the pixel value β ^(k) (p) = β ^(k) _ij of the k-th term of the molecular image β for each pixel p according to the equation (35a). The GF calculator 29 performs the Gaussian convolution calculation GF (•) shown in the equations (34a) and (34b).

The denominator image memory 30 stores the denominator image β ₀ output from the denominator image calculator 27. The molecular image memory 31 stores the molecular image β output from the molecular image calculator 28. The adders 30a and 31a read the pixel values β _0ij and β _ij of the respective pixels of the denominator image β ₀ and the numerator image β from the denominator image memory 30 and the image memory 31, respectively, and the denominator image calculator 27 and the numerator image calculator. The pixel values β _0ij and β _ij of the denominator image memory 30 and the image memory 31 are updated by adding the values β ^(k) _0ij and β ^(k) _ij of the k-th term output by the No. 28.

  The image divider 32 calculates the output image z by the equation (33a) after the calculations of the equations (34b) and (34c) by the denominator image calculator 27 and the numerator image calculator 28 are completed, and saves them in the output image memory 22. To do. The output image z stored in the output image memory 22 is output to the output device 3.

The operation of the filter arithmetic device 20 according to the present embodiment configured as described above will be described below. FIG. 13 is a PAD showing the flow of bilateral calculation processing in the filter calculation device 20 of the second embodiment. As in the first embodiment, first, (K _τ , T _τ ) is set so that the approximate error E ^ in the equation (18) becomes the allowable error τ and K is minimized, and then the significant coefficient (Fourier coefficient) a Calculate _k and filter function.

(Step S401)
First, the truncation term number calculator 13 reads each parameter (R, σ, τ) set in the computation parameter memory 11, and calculates the optimum number of truncation terms K _τ for the allowable error τ by Expression (20). To the KT register 14a.

(Step S402)
Next, the significant cycle optimizer 14 reads each parameter (R, σ, τ) set in the calculation parameter memory 11 and solves the equation (22) of T in the closed interval of the equation (19a). calculates the optimal significant period T _tau for tolerance tau, saving to KT register 14a. The details of the (K, T) optimum value search process are as described in the first embodiment.

(Step S403)
Next, the term counter 23 sets the current term number k to 0, and the significant coefficient calculator 15 calculates the zeroth Fourier coefficient a ₀ (= 1 / 2T _τ ) and stores it in the significant coefficient memory 16. To do.

(Step S404)
Next, the k = 0th term (DC term) of the denominator image β ₀ and the numerator image β is calculated. The denominator image calculator 27 sets all pixel values β _0ij of the denominator image β ₀ to a ₀ and stores them in the denominator image memory 30. The molecular image calculator 28 calculates a pixel value β _ij for each pixel (i, j) of the molecular image β by a ₀ GF (σ _s , (i, j), {x _kl }) Save to the image memory 31. Where Gaussian convolution operation

Is executed by the GF calculator 29.

(Step S405)
Next, while the term number counter 23 increments the term number k by 1, the calculation of the k (> 0) -th term (AC term) of the denominator image β ₀ and the numerator image β is performed by the following steps S406 to S411. Execute by repeating sequentially.

(Step S406)
First, the significant coefficient calculator 15 calculates the k-th Fourier coefficient a _k by the equation (15) using T _{τ of the} KT register 14a and the scale parameter σ _r of the calculation parameter memory 11, and the significant coefficient memory 16 Save to

(Step S407)
Next, the base function calculator 25 uses the pixel values of T _tau and the input image memory 21 of the KT register 14a the (x _ij), each pixel (i, j) base function cos term number k for (kπx _ij The value of / T _τ ), sin (kπx _ij / T _τ ) is calculated and stored in the basis function table 26.

(Steps S408-S409)
Next, the k-th term of the denominator image is calculated. First, for each pixel (i, j), the denominator image calculator 27 uses the scale parameter σ _s of the calculation parameter memory 11 and the cosine basis c _ij and sine basis s _ij of the basis function table 26 to convolution Ω. _Cij = GF (σ _s , (i, j), {c _kl }) and Ω _sij = GF (σ _s , (i, j), {s _kl }) are calculated (S408). Where Gaussian convolution operation

Is executed by the GF calculator 29. Next, the denominator image calculator 27 calculates the value a _k (c _ij Ω _cij + s _ij Ω _sij ) of the denominator image using the Fourier coefficient a _k of the significant coefficient memory 16, and the adder 30a This is added to the pixel value β _0ij in the denominator image memory 30 to update the pixel value β _0ij in the denominator image memory 30. This is executed for all pixels.

(Steps S410-S411)
Next, the k-th term of the molecular image is calculated. First, for each pixel (i, j), the molecular image calculator 28 uses the scale parameter σ _s of the calculation parameter memory 11 and the cosine basis c _ij and sine basis s _ij of the basis function table 26 to convolution Ω. _cij = GF performed _{(σ s, (i, j} ), {c kl x kl}), Ω sij = GF calculation of _{(σ s, (i, j} ), {s kl x kl}) (S410). Where Gaussian convolution operation

Is executed by the GF calculator 29. Next, the molecular image calculator 28 uses the Fourier coefficient a _k of the significant coefficient memory 16 to calculate the value a _k (c _ij Ω _cij + s _ij Ω _sij ) of the molecular image, and the adder 31a This is added to the pixel value β _ij of the molecular image memory 31 to update the pixel value β _ij of the molecular image memory 31. This is executed for all pixels.

(Step S407)
After the processing of steps S405~S411 is completed, the image divider 32, for each pixel (i, j), the pixel value beta _ij molecules image memory 31 by dividing the pixel value beta _0Ij denominator image memory 30 The pixel value z _ij of the output image z is calculated and stored in the output image memory 22. When the pixel values of the output image z are calculated for all the pixels, the output image z is output to the output device 3, and the bilateral calculation process is terminated.

  With the above operation, an output image z obtained by performing bilateral arithmetic processing on the input image x can be obtained. Next, the Gaussian convolution calculation GF (σ, p, {f}) by the GF calculator 29 will be supplementarily described. The Gaussian convolution calculation is performed using a known algorithm here, but in this embodiment, an example configured with reference to those described in Non-Patent Documents 18 and 19 is shown as an example.

[2.1] GF Operation Unit The GF operation unit 29 performs Gaussian convolution calculation (Gaussian filter calculation) using a constant time algorithm per pixel. Since the actual Gaussian convolution operation is performed using discrete values, the discrete value of the continuous function f (p) of _p is expressed as f _p below, and the Gaussian kernel g (σ, The discrete value of p) is expressed as g _p (σ is omitted unless otherwise required, and p is a scalar x or a vector (x, y)). Also, the Gaussian convolution operation GF (σ, p, {f (q)}) at the position p is abbreviated as GF [f _p ]. The GF calculator 29 performs a two-dimensional Gaussian filter operation. Since the Gaussian kernel g (σ, q) can be separated into variables with respect to the two-dimensional coordinates (x, y), the Gaussian filter operation is performed in the x direction. It can be performed independently in the y direction. That is,

Here, GF _x [•] and GF _y [•] represent Gaussian filter operations for the x-coordinate and y-coordinate, respectively, and W represents the size of the window region where the Gaussian filter operation is performed. Usually, the window size W is set in a range where g _x is significant (for example, W = 3σ when the significance level is 99.73%).

Therefore, it is sufficient to consider constant time per pixel for one-dimensional Gaussian filter calculation. First, when Gaussian kernel g _x is expanded by DCT-5, it becomes as follows.

Here, _{G k} is the k-th DCT coefficients of _{g x, K 1} is the discontinuation number of terms. K ₁ may be set to K ₁ = 3 / (φσ) if the DCT coefficient G _k is in a significant range, for example, if the significance level is 99.73%. Thus, convolution of a function f _x and the Gaussian kernel g _x is expressed as follows.

_Here, ^{F k (x)} is the k-th short DCT coefficients of the input sequence _{f x} at position _x, ^{C k (u)} is a cosine kernel. The short-time DCT coefficient F _k ^(x) can be expressed by the following initial value and recurrence formula, and by calculating using this recurrence formula, the calculation amount of formula (35) per pixel is O (1).

  In particular, when k = 0, the recurrence formula is simplified as follows.

In the actual calculation, rounding errors in C _k ⁽¹⁾ and C _k ^(R) propagate in the equation (36b). In order to suppress this error propagation, at the right side of equation (35a), out factoring the maximum weighting factor G _0.

At this time, the recurrence formula for Z _k ^(x) is as follows.

As described above, the one-dimensional Gaussian filter operation can be executed by the following algorithm.
(1) First, G ₀ , {γ _k C _k ^(u) | k = 0,..., K ₁ ; u = −W,..., W}, {2C _k ⁽¹⁾ | k = 0,. K ₁ } is calculated in advance and prepared as a lookup table (LUT).
(2) First, F ₀ ⁽⁰⁾ , Z _k ⁽⁰⁾ , Z _k ⁽¹⁾ (k = 1,..., K ₁ ) are calculated as initial values according to the equations (36a) and (37b), First, (f * g) ₀ is calculated by the equation (37a).
(3) Next, F ₀ ^(x) is updated by the recurrence formula (36b).
(4) Next, (f * g) _x is calculated by the equation (37a), F ₀ ^(x) is updated by the recurrence equation (36b), and Z _k ^(x−1 ) is calculated by the recurrence equation (38). ⁾ , Z _k ^(x) (k = 0,..., K ₁ ) is updated from x = 1 to N−1 (N is an input sequence (f ₀ , f ₁ ,..., F _N−1 )). (F * g) _x for all x.

Incidentally, f _x in x <0 and x ≧ N, for example, may be supplemented by the following extrapolation rule.

From the above, the two-dimensional discrete Gaussian filter calculation GF [f _xy ] by the GF calculator 29 can be executed by the process shown in the PAD of FIG.

(1) Flow of overall processing (two-dimensional discrete Gaussian filter processing) (FIG. 14A)
First, the GF computing unit 29 performs an input function {f _xy } (given as an N × M matrix) that is a function of two-dimensional coordinates (x, y) at all coordinate positions (N × M points). A one-dimensional Gaussian filter process is executed for the y coordinate, and the function after the filter process is set to {f ^(y) _xy } (S501). The function {f ^(y) _xy } is an N × M matrix. Next, the GF computing unit 29 executes a one-dimensional Gaussian filter process on the y coordinate at all coordinate positions (N × M points) for the function {f ^(y) _xy } (S502). The function {f ^(xy) _xy } (N × M matrix) is output as a function after the two-dimensional discrete Gaussian filter processing (S503). Here, the one-dimensional Gaussian filter processing for the x-coordinate and y-coordinate in steps S501 and S502 is performed as shown in FIG.

(2) Flow of one-dimensional discrete Gaussian filter processing (FIG. 14B)
FIG. 14B shows details of the processing of GF _y [•], GF _x [•] in steps S501 and S502 (GF _x is shown as a representative). The input function {f _x } is a function of the coordinate x (x = 0, 1,..., N−1). Also, in FIG. 14B, the value indicated by {} is calculated in advance by the GF calculator 29 and stored as an LUT in the internal memory.

(Step S601) First, the GF calculator 29 calculates the initial value F ₀ ⁽⁰⁾ of the recurrence formula (44b ⁾ using the formula (44a), and stores the result in the variable F ₀ . Here, C ₀ ^(u) = 1.

(Steps S602 to S603) Next, the GF computing unit 29 converts the initial values Z _k ⁽⁰⁾ , Z _k ⁽¹⁾ (k = 1,..., K ₁ ) of the recurrence formula (47) into the formula (46b). , (44a) and store the results in variables Z _k ⁻ and Z _k , respectively.

(Step S604) Next, GF calculator 29 is now variable _{F 0} and ^{_{{Z k -; k = 1}} , ..., K 1} with the value stored in the, by the equation (46a), first x Calculate (f * g) ₀ at = 0.

(Step S605) Next, the GF calculator 29 calculates F ₀ ⁽¹⁾ by the recurrence formula (45) using the variable F ₀ and the input function {f _x }, and sets the result to the variable F ₀ . Store again. Here, x function value f _x for <0 are complemented by extrapolation rule of formula (48) (in the case of the function value f _x of x ≧ N required after the same).

(Steps S606 to S610) Next, for x = 1,..., N−1, the following steps S606 to S610 are repeatedly executed while incrementing _x, thereby sequentially calculating (f * g) _x .

First, GF calculator 29 sets the value of x (S606), the current, variable _{F 0} and ^{_{{Z k -; k = 1}} , ..., K 1} with the value stored in the formula (46a ), (F * g) _x in the current _x is calculated (S607). Next, the GF calculator 29 calculates F ₀ ^{(x + 1)} by the recurrence formula (45) using the current variable F ₀ and the input function {f _x }, and stores the result in the variable F ₀ again. At the same time, Δ ^(x) is calculated by the equation (44c) using the input function {f _x } (S608). Finally, the GF calculator 29 calculates {Z _k ⁻ ; k = 1,..., K ₁ }, {Z _k ; k = 1,... For each k (k = 1,..., K ₁ ). using the value stored in the _{K 1},} it calculates the _Z ^{k (x + 1)} by a recursion formula (47), _Z ^k currently stored in _{Z k} a ^(x) _{Z k} ^- stored again in At the same time, the newly calculated Z _k ^{(x + 1)} is stored again in Z _k (S609, S610).

(Step S611) {(f * g) _x ; x = 0,..., N−1} calculated by the above processing is output as a function after the one-dimensional discrete Gaussian filter processing.

[2.2] Example Using Box Filter for Spatial Kernel In the above embodiment, a Gaussian filter (GF) calculator is used as the calculator for the spatial kernel. However, in the present invention, Instead of the GF calculator, a box filter calculator (BF calculator) can be used.

  A box spatial kernel can be defined as follows.

Here, “‖ · ‖ _∞ ” represents l _∞ −norm, and takes the value of the largest of the elements of the vector q.

  Using the above box space kernel, the BF calculator is defined as a convolution of an arbitrary function f (p) in the space domain and the box space kernel as follows.

  When a BF calculator is used, GF (•) in formulas (35a) and (35b) may be replaced with BF (•).

[2.3] Experimental Example (1) Basic Performance in Filtering Finally, the basic performance in filtering of the filter arithmetic device 20 according to the second embodiment will be described in comparison with a conventional bilateral filter arithmetic algorithm. .

In this experiment, the performance trade-off between approximation accuracy and execution time was evaluated in two-dimensional filtering. In the experiment, a Kodak Photo CD including a 24RGB image having a size of 768 × 512 pixels or 512 × 768 pixels was used, and average approximation accuracy and average execution time were measured for 24 images. The execution time includes memory allocation time and calculation time. The approximation accuracy of the target algorithm was quantified by the peak signal-to-noise ratio (PSNR) [dB] between the original image and the smoothed image generated from the original image by the target algorithm and the original BLF. It is also assumed that a Gaussian space kernel or Gaussian range kernel that supports ± 4σ _s or ± 4σ _r achieves accurate accuracy in the original BLF. In general, it is evaluated that a sufficient approximation accuracy can be realized at about 50 dB.

First, we examine the superiority of CBLF over the Chaudhury approximation algorithm using Gaussian space / range kernels. FIG. 15 is a diagram comparing approximation accuracy [dB] and calculation time in the Chaudhury approximation algorithm (ε = 0.8) and CBLF (τ = 0.1) of Non-Patent Document 2. FIG. 15A shows the relationship of the approximation accuracy (PSNR) [dB] with respect to the range scale parameter σ _r when the spatial scale parameter σ _s is 2 and 6. FIG. 15B shows the relationship of the calculation time [s / Mpixels] per 1M pixel with respect to the range scale parameter σ _r when the spatial scale parameter σ _s is 2 and 6. In each figure, “Chaudhury '13” indicates the result of the Chaudhury approximation algorithm, and “Ours” indicates the result of the filter arithmetic unit 20 (CBLF) according to the present invention. In the Chaudhury approximation algorithm, the allowable error ε is 0.8, and in the filter arithmetic unit 20 (CBLF) according to the present invention, the allowable error τ is 0.1. As a result of the experiment, the filter operation device 20 according to the present invention has a calculation speed about 5 to 8 times faster than that in the case where the Chaudhury approximation algorithm is used (FIG. 8A) with the same approximation accuracy (about 50 dB). It can be seen (FIG. 8B). Here, in the calculation of the spatial filter in the filter arithmetic unit 20 according to the present invention, the constant time GF arithmetic unit described in [2.1] is used. These results also demonstrate that the CBLF execution time does not depend on σ _s . That is, the calculation amount per pixel is O (1). The inventor conducted thorough experiments for various σ _{s and} confirmed this tendency. Therefore, it has been proved that the filter operation device 20 (CBLF) according to the present invention is considerably superior to the conventional Chaudhury approximation algorithm in theory and experiment.

Next, the superiority of the filter arithmetic unit 20 (CBLF) according to the present invention over the Yang approximation algorithm of Non-Patent Document 14 using the Gaussian range kernel and the box space kernel (formula (49a)) will be confirmed.
According to Non-Patent Document 14, the Yang approximation algorithm can achieve sufficient accuracy by setting the size B of the box space kernel to B = 8 in many cases.
FIG. 16 is a diagram comparing the calculation accuracy with the approximation accuracy [dB] in Yang approximation algorithm (B = 8) and CBLF of Non-Patent Document 14. FIG. 16A shows the relationship of the approximation accuracy (PSNR) [dB] with respect to the range scale parameter σ _r when the spatial scale parameter σ _s is 2 and 6. FIG. 16B shows the relationship of the calculation time [s / Mpixels] per 1M pixel with respect to the range scale parameter σ _r when the spatial scale parameter σ _s is 2 and 6. In each figure, "Yang + '09 (B = 8)" and "Yang + '09 (B = 4)" show the results of the Yang approximation algorithm for B = 8, 4, respectively, and "Ours" is the present invention. The result of the filter arithmetic unit 20 (CBLF) concerning is shown.
From FIG. 16A, in order to obtain an allowable error comparable to the Yang approximation algorithm (B = 8), in the filter arithmetic unit 20 (CBLF) according to the present invention, τ = 2 with respect to σ _s = 2. It can be seen that it is appropriate to set τ = 0.07 for σ _s = 6. As is clear from FIG. 16, in CBLF, a zigzag change in approximation accuracy and a stepwise change in execution time are observed. This is a phenomenon caused by a change in the number K of significant spectral terms. The number K of significant spectral terms for a given allowable error τ is determined from the intersection of a contour line of K and a straight line of the given allowable error τ (dotted line in FIG. 11) in FIG. Note that the number M of spatial filters is M = 2B in the Yang approximation algorithm and M = 4K + 1 in CBLF. Further, as shown in FIG. 16A, in the Yang approximation algorithm, it is possible to obtain sufficient approximation accuracy with B = 8 for a wide range of σ _r , but with B = 4, sufficient approximation accuracy can be obtained. I can't understand. In quantization-based algorithms such as the Yang approximation algorithm, the approximate lower bound of M is caused by insufficient approximation accuracy. As a result, redundancy cannot be eliminated as σ _r increases. On the other hand, in CBLF, M can be reduced to about half when 0.15 ≦ σ _r . Therefore, CBLF has a better performance trade-off than the Yang approximation algorithm, particularly when σ _r is large (FIG. 16B).

(2) Visual Evaluation of Noise Finally, the result of visual evaluation of the approximation error of the original BLF, Chaudhury approximation algorithm, Yang approximation algorithm, and CBLF will be described. As the test images, three RGB images of 256 × 256 pixels in the standard image Data-BAse3 were used.

  FIG. 17 shows an original image, a smoothed image generated by the original BLF, Chaudhury approximation algorithm, Yang approximation algorithm, and CBLF algorithm, and an error image amplified 10 times. All O (1) BLF algorithms appear to produce a smoothed image that closely approximates the original BLF, and in human vision these differences appear to be small. Obviously, the Chaudhury approximation algorithm and CBLF produce almost the same error image. This is because both are based on a similar framework based on cosine term expansion. On the other hand, the Yang approximation algorithm appears to show some offset error in each segment. In all these O (1) algorithms, errors occur as high frequency noise mainly at sharp edges. Since these edges are visually amplified, the actual error visually felt is sufficiently small and fine.

  Next, an embodiment in which the filter arithmetic device according to the present invention is configured by a parallel arithmetic circuit such as a coprocessor or FPGA will be described.

  FIG. 18 is a block diagram illustrating the overall configuration of the filter arithmetic apparatus according to the third embodiment. 18A is an overall view, and FIG. 18B is a configuration diagram of an input line buffer and an output line buffer. The filter operation device 50 of this embodiment includes an operation parameter memory 51, an approximate parameter optimization unit 52, an input line buffer array 53, a filter operation circuit 54, and an output line buffer array 55.

  The calculation parameter memory 51 is a memory for storing calculation parameters necessary for the filter calculation. The approximate parameter optimization unit 52 is a circuit block that calculates each parameter necessary for the approximate calculation of the Gaussian kernel based on the calculation parameter stored in the calculation parameter memory 51. The input line buffer array 53 is a memory in which a partial image cut out from the original image is input and temporarily stored. Let the size of the original image be W × H pixels (horizontal W, vertical H). The input line buffer array 53 can store partial images of (2M + 1) rows (where 2M + 1 is a filter window size and is set by a calculation parameter). As shown in FIG. 18B, the input line buffer array 53 is configured by arranging 2M + 1 rows of W register arrays, and the register column of each column is a shift register, and a register is provided for each row. It is comprised so that the value of can be read. Current image data is input to the input line buffer array 53 line by line in order from the upper end line of the original image. The filter operation circuit 54 is a circuit that performs a filter operation on the original image held in the input line buffer array 53 and outputs a smoothed image. The output line buffer array 55 is a memory that temporarily holds the smoothed image output from the filter arithmetic circuit 54, and has the same configuration as the input line buffer array 53.

  FIG. 19 is a block diagram showing the configuration of the calculation parameter memory 51 and the approximate parameter optimization unit 52 of FIG.

The calculation parameter memory 51 stores values of the allowable error τ, the dynamic range R of the original image, the scale parameter σ _r of the range region, the scale parameter σ _s of the spatial region, and the filter size parameter M as calculation parameters. . The values of these calculation parameters are input from the outside.

The approximate parameter optimization unit 52 includes a ζ _τ operation circuit 61, a ζ _τ register 62, a K operation circuit 63, a K register 64, a D operation circuit 65, a D register 66, a T operation circuit 67, a T register 68, and an _ak operation circuit. 69, an _ak register 70, an m computing unit 71, and an m register 72.

zeta _tau arithmetic circuit 61 on the basis of the tolerance tau stored in the calculation parameter memory 51, calculates the zeta _tau by performing the calculation of the inverse complementary error function of Equation (19b), stores the zeta _tau register 62. The ζ _τ operation circuit 61 is a circuit that performs a square operation and an inverse complementary error function operation, and can be configured by a look-up table (LUT).

K arithmetic circuit 63, the dynamic range R stored in the operation parameter memory 51, based on the scale parameter sigma _r, and zeta _tau stored in zeta _tau register 62, by performing the calculation of Equation (20), discontinuation The number of terms K is calculated and stored in the K register 64.

The D arithmetic circuit 65 has a lower limit value (σ _r ζ _τ + R) and an upper limit value of the range D = [σ _r ζ _τ + R, πσ _r (2K + 1) / ζ _τ ] of the significant period T shown in Expression (19a). (Πσ _r (2K + 1) / ζ _τ ) is calculated, and these values are stored in the D register 66.

  The T arithmetic circuit 67 is a circuit that calculates the significant period T by solving the equation (22) by the bisection method in the range D of the significant period T stored in the D register 66. The significant period T calculated by the T arithmetic circuit 67 is stored in the T register 68.

The a _k arithmetic circuit 69 reads the number of truncation terms K stored in the K register 64, the scale parameter σ _r stored in the arithmetic parameter memory 51, and the significant period T stored in the T register 68, and k = Each coefficient a _k is stored in the a _k register 70 by calculating Expression (15) from 1 to K. Here, since K is variable, the _ak register 70 needs a sufficient number of registers corresponding to the change width of K. From FIG. 11, the upper limit of K is about 4 to 6 at most. .

  The m calculator 71 calculates the filter window size m = 2M + 1 from the filter size parameter M stored in the calculation parameter memory 51 and stores it in the m register 72.

  FIG. 20 is a block diagram showing the configuration of the filter arithmetic circuit 54 of FIG. The filter operation circuit 54 includes a spatial filter operation unit 75, a DC component b register array 76, a DC component b0 register array 77, a B circuit 78-k (k = 1 to K), and a b register array 79-k (k = 1 to 1). K), b0 register array 80-k (k = 1 to K), adders 81 and 82, and a divider 83.

  The spatial filter calculator 75 is the GF calculator described in the second embodiment, and reads out pixel data (x [−M] to x [M]) for m (= 2M + 1) rows from the input line buffer array 53. The spatial filter operation of the DC component shown in the equation (38) is sequentially performed on each pixel in the center row (0th row (see FIG. 18B)), and W approximated molecular images obtained as a result thereof. The value of the DC component (k = 0) of β ^ (formula (35a)) is stored in the DC component b register array 76. The DC component b register array 76 is a register that holds the value of the CD component of one row (W) of approximate molecular images β ^ (formula (35a)).

The DC component b ₀ register array 77 is a register that holds the value of the DC component (k = 0) of one row (W) of approximate denominator images β ₀ ^ (formula (35b)). Has been.

The B circuit 78-k (k = 1 to K) is a circuit that calculates the k-th frequency component of the approximate numerator image β ^ and the approximate denominator image β ₀ ^ in the expressions (35a) and (35b). The b register array 79-k (k = 1 to K) is a register that holds the k-th frequency component of one row (W pieces) of approximate molecular images β ^ calculated by the B circuit 78-k. The b ₀ register array 80-k (k = 1 to K) is a register that holds the k-th frequency component of one row (W) of approximate denominator images β ₀ ^ calculated by the B circuit 78-k. .

  The adder 81 adds the frequency components of the approximate molecular image β ^ of each pixel (W pieces) held in the b register arrays 76 and 79-1 to 79-K to each pixel, and adds the expression (35a ) To calculate the value of the approximate molecular image β ^ for each pixel of one row.

The adder 82 adds the frequency components of the approximate denominator image β ₀ ^ of each pixel (W pixels) held in the b ₀ register arrays 77 and 80-1 to 80-K to each pixel, and The sum operation of the right side of (35b) is performed, and the value of the approximate denominator image β ₀ ^ for each pixel of one row is calculated.

The divider 83 divides the value of the approximate numerator image β ^ output from the adder 81 by the value of the approximate denominator image β ₀ ^ output from the adder 82 for each pixel, and is expressed by Expression (33a). The final smoothed image is calculated and output to the output line buffer array 55.

  FIG. 21 is a block diagram showing a configuration of B circuit 78 (B circuits 78-1 to 78-K) in FIG.

The B circuit 78 includes a pixel value selector 91, a term number register 92, a phase calculator 93, a _ck calculator 94, a _kk register array 95, a multiplier 96, a cf _k register array 97, an _sk calculator 98, and _sk. A register array 99, a multiplier 100, an sf _k register array 101, spatial filter arithmetic units 102 to 105, multipliers 106, 107, 110, 111, 109, 113, and adders 108, 112 are provided.

  The pixel value selector 91 sequentially selects the pixels q from the original image data of the pixels held in the input line buffer array 53, and outputs the pixel value x (q) of the pixel q. Here, the pixel values of the original image are additionally input to the input line buffer array 53 for each row. However, the pixel value selector 91 selects x [− in the newly added row (FIG. 18B). It is assumed that each pixel (W pixels) in row M] is sequentially selected.

  The term number register 92 is a register that holds the value of the term number k (k = 1,..., K) of the B circuit 78.

The phase calculator 93 calculates and outputs the phase value φ _k (q) = 2πkx (q) / T for the pixel q selected by the pixel value selector 91.

The c _k calculator 94 calculates and outputs a cosine value c _k (q) = cos (φ _k (q)) of the phase value φ _k (q) output from the phase calculator 93. The c _k register array 95 is a register that stores the cosine value c _k (q) output from the c _k arithmetic unit 94. The _ck register array 95 is a shift register of (2M + 1) rows and W columns having the same structure as that of the line buffer array 53 shown in FIG. 18B, and the cosine value c _k (q ) Is added to the -Mth line.

The multiplier 96 multiplies the cosine value c _k (q) output from the _ck arithmetic unit 94 by the pixel value x (q) output from the pixel value selector 91 to obtain a multiplied value c _k (q) x (q ) Is output. The cf _k register array 97 is a register that stores the multiplication value c _k (q) x (q) output from the multiplier 96. The cf _k register array 97 is also a shift register of (2M + 1) rows and W columns having the same structure as the line buffer array 53 shown in FIG. 18B, and the multiplication value c _k (q ) x (q) is added to the −Mth row.

The s _k calculator 98 calculates and outputs a sine value s _k (q) = sin (φ _k (q)) of the phase value φ _k (q) output from the phase calculator 93. The s _k register array 99 is a register that stores the sine value s _k (q) output from the s _k calculator 98. The s _k register array 99 is also a shift register of (2M + 1) rows and W columns having the same structure as that of the line buffer array 53 shown in FIG. 18B, and the sign value s _k (q ) Is added to the -Mth line.

The multiplier 100 multiplies the sine value s _k (q) output from the s _k calculator 98 by the pixel value x (q) output from the pixel value selector 91, and the multiplied value s _k (q) x (q ) Is output. The sf _k register array 101 is a register that stores a multiplication value s _k (q) x (q) output from the multiplier 100. The sf _k register array 101 is also a shift register of (2M + 1) rows and W columns having a structure similar to that of the line buffer array 53 shown in FIG. 18B, and the multiplication value s _k (q ) x (q) is added to the −Mth row.

The spatial filter computing unit 102 calculates the Gaussian filter function GF (σ _s , p, {c _k (q)}) (p is the first row in the center row (FIG. 18B). This is a spatial filter that performs an operation on W pixels) in the 0th row). σ _s is input from the operation parameter memory 51, and {c _k (q)} is input from the _ck register array 95.

The spatial filter calculator 103 calculates the Gaussian filter function GF (σ _s , p, {s _k (q)}) (p is the first row in the center row (FIG. 18B). This is a spatial filter that performs an operation on W pixels) in the 0th row). σ _s is input from the operation parameter memory 51, and {s _k (q)} is input from the _sk register array 99.

The multiplier 106 multiplies the value GF (σ _s , p, {c _k (q)}) output from the spatial filter calculator 102 by the cosine value c _k (p) read from the _ck register array 95. And the multiplication value c _k (p) GF (σ _s , p, {c _k (q)}) of the first term on the right side of the equation (35b) is output.

Multiplier 107, the value _{GF (σ s, p, {} s k (q)}) for outputting the spatial filter operation unit 103 and the multiplication of the sign value s _k read from s _k register array 99 (p) The multiplication value s _k (p) GF (σ _s , p, {s _k (q)}) of the second term on the right side of the equation (35b) is output.

The adder 108 adds the output value of the multiplier 106 and the output value of the multiplier 107, and calculates the sum in the right parenthesis of Expression (35b). The multiplier 109 multiplies the output value of the adder 108 by the coefficient value a _k read from the a _k register 70, and calculates the approximate denominator image β ^ ₀ (p) of Expression (35b) for each pixel p. and, save to _{b 0} register array 80-k.

The spatial filter computing unit 104 calculates the Gaussian filter function GF (σ _s , p, {c _k (q) x (q)}) (p in the middle row (FIG. 18 ( It is a spatial filter that performs the calculation of W pixels) in the 0th row) of b). σ _s is input from the operation parameter memory 51, and {c _k (q) x (q)} is input from the cf _k register array 97.

The spatial filter computing unit 105 calculates the Gaussian filter function GF (σ _s , p, {s _k (q) x (q)}) (p in the middle row (FIG. 18 ( It is a spatial filter that performs the calculation of W pixels) in the 0th row) of b). σ _s is input from the operation parameter memory 51, and {s _k (q) x (q)} is input from the sf _k register array 101.

The multiplier 110 outputs the value GF (σ _s , p, {c _k (q) x (q)}) output from the spatial filter calculator 104 and the cosine value c _k (p) read from the _ck register array 95. And the multiplication value c _k (p) GF (σ _s , p, {c _k (q) x (q)}) of the first term on the right side of the equation (35a) is output.

The multiplier 111 outputs the value GF (σ _s , p, {s _k (q) x (q)}) output from the spatial filter calculator 105 and the sine value s _k (p) read from the s _k register array 99. And the multiplication value s _k (p) GF (σ _s , p, {s _k (q) x (q)}) of the second term on the right side of the equation (35a) is output.

The adder 112 adds the output value of the multiplier 110 and the output value of the multiplier 111, and calculates the sum in the right parenthesis of Expression (35a). The multiplier 113 multiplies the output value of the adder 112 by the coefficient value a _k read from the a _k register 70, and calculates the approximate molecular image β ^ (p) of the equation (35a) for each pixel p. , B are stored in the register array 79-k.

  With the configuration as described above, the bilateral filter can be operated in the same manner as in the second embodiment.

DESCRIPTION OF SYMBOLS 1 Computer 2 Input device 3 Output device 4 External storage device 10 Gaussian kernel operation device 11 Operation parameter memory 12 Approximate parameter optimization means 13 Censored term number calculator 14 Significant period optimizer 14a KT register 15 Significant coefficient calculator 16 Significant Coefficient memory 17 Approximate function calculator 20 Filter calculator 21 Input image memory 22 Output image memory 23 Term counter 25 Base function calculator 26 Base function table 27 Denominator image calculator 28 Numerator image calculator 29 GF calculator 30 Denominator image memory 30a adder 31 molecular image memory 31a adder 32 image divider 50 filter operation device 51 operation parameter memory 52 approximate parameter optimization unit 53 input line buffer array 54 filter operation circuit 55 output line buffer array 61 ζ _τ operation circuit 62 ζ _τ register 63 K arithmetic circuit 64 K register 65 D arithmetic circuit 66 D register 67 T arithmetic circuit 68 T register 69 a _k arithmetic circuit 70 a _k register 71 m arithmetic unit 72 m register 75 Spatial filter arithmetic unit 76 DC component b register array 77 DC component b ₀ register array 78, 78-k (k = 1 to K) B circuit 9-k (k = 1 to K) b register array 7
80-k (k = 1 to K) b ₀ register array 81, 82 adder 83 divider 91 pixel value selector 92 term number register 93 phase calculator 94 _ck calculator 95 _ck register array 96 multiplier 97 cf _k Register array 98 s _k arithmetic unit 99 s _k register array 100 multiplier 101 sf _k register array 102 to 105 spatial filter arithmetic units 106, 107, 110, 111, 109, 113 multipliers 108, 112 adders

Claims (12)

The input image {x (p)} (where p is the pixel position and x (p) is the pixel value at position p) is defined by a Gaussian function with a standard deviation (hereinafter referred to as “scale parameter”) σ _s The spatial kernel g _s (qp) (p and q are pixel positions) and the range kernel g _r (x (q) -x (p)) defined by the Gaussian function of the scale parameter σ _r An image filter arithmetic device that smoothes an image by a volume calculation,
For the scale parameter σ _r , the predetermined tolerance τ and the dynamic range R of the input image {x (p)}, the range kernel g _r (x (q) −x (p)) is a finite term. The series truncation term K in the approximate range kernel g ^ _r (x (q) -x (p)) approximated by the Fourier series of the number and the approximate range kernel g ^ _r (x (q) -x ( p)) an approximation that is a set of significant periods T that define an integration region [-T / 2, T / 2] of Fourier integration when calculating the expansion coefficient a _k (k = 0,..., K) of each term in p)) Approximate parameter optimizing means for determining and outputting parameters (T, K) such that the approximate error of the approximate range kernel g ^ _r (x (q) -x (p)) is equal to or less than the allowable error τ; ,
Expansion coefficients from the 0th term to the Kth term of the approximate range kernel g ^ _r (x (q) -x (p)) according to the approximate parameter (T, K) calculated by the approximate parameter optimization means expansion coefficient calculating means for calculating {a ₀ ,..., a _K };
Using the expansion coefficients {a ₀ ,..., A _K } calculated by the expansion coefficient calculating means, the spatial kernel g _s (qp) and the approximate range kernel g are applied to the input image {x (p)}. An image filter arithmetic device comprising: filter arithmetic means for executing a convolution operation with ^ _r (x (q) -x (p)). The approximate parameter optimization means includes truncation number calculation means for calculating the truncation number K.
The truncation term number calculating means calculates the truncation term number K with respect to the scale parameter σ _r , a predetermined allowable error τ, and the dynamic range R of the input image {x (p)}.
The image filter arithmetic device according to claim 1, wherein the image filter arithmetic device is calculated by: The approximate parameter optimization means includes significant period optimization means for calculating the significant period T,
The expansion coefficient calculating means is configured to calculate an approximate expansion coefficient obtained by performing Fourier integration of the range kernel g _r (x (q) −x (p)) in an integration region [−T / 2, T / 2], and each expansion coefficient. {a ₀ , ..., a _K }
The significant period optimization means includes an equation of a variable T with respect to the scale parameter σ _r , the truncation term number K, and the dynamic range R.
3. The image filter arithmetic device according to claim 1, wherein T satisfying the condition is calculated as the significant period T. 4. The approximate parameter optimization means includes:
Based on the scale parameter σ _{r, the} range upper limit value R, and a preset effective width parameter α (ασ _r <R), the significant period T is
A calculation cycle setting means calculated by:
Based on the allowable error τ, the significant period T, and the scale parameter σ _r , the truncation number K is
The image filter arithmetic apparatus according to claim 1, further comprising: a truncation term number determining means for calculating by the following. The expansion coefficient computing means calculates the expansion coefficients {a ₁ ,..., A _K },
The image filter arithmetic device according to claim 1, wherein the image filter arithmetic device is calculated by:   A program that causes a computer to function as the image filter arithmetic device according to any one of claims 1 to 5 by being read and executed by a computer. A Gaussian kernel arithmetic unit that calculates an approximate value of a Gaussian kernel g _r (x) defined by a Gaussian function of standard deviation (hereinafter referred to as “scale parameter”) σ _r ,
For the scale parameter σ _r , the predetermined tolerance τ and the dynamic range R of the input image {x (p)}, the range kernel g _r (x (q) −x (p)) is a finite term. The series truncation term K in the approximate range kernel g ^ _r (x (q) -x (p)) approximated by the Fourier series of the number and the approximate range kernel g ^ _r (x (q) -x ( p)) an approximation that is a set of significant periods T that define an integration region [-T / 2, T / 2] of Fourier integration when calculating the expansion coefficient a _k (k = 0,..., K) of each term in p)) Approximate parameter optimizing means for determining and outputting parameters (T, K) such that the approximate error of the approximate range kernel g ^ _r (x (q) -x (p)) is equal to or less than the allowable error τ; ,
The expansion coefficients {a ₀ ,..., A _K } from the 0th term to the Kth term of the approximate Gaussian kernel g _r (x) according to the approximate parameters (T, K) calculated by the approximate parameter optimization means. Expansion coefficient calculating means for calculating
Approximation function calculation means for calculating the approximate Gaussian kernel g _r (x) using the expansion coefficients {a ₀ ,..., A _K } calculated by the expansion coefficient calculation means, and a Gaussian kernel calculation device comprising: . The approximate parameter optimization means includes truncation number calculation means for calculating the truncation number K.
The truncation term number calculating means calculates the truncation term number K with respect to the scale parameter σ _r , a predetermined allowable error τ, and the dynamic range R of the input image {x (p)}.
8. The image filter arithmetic device according to claim 7, wherein the image filter arithmetic device is calculated by: The approximate parameter optimization means includes significant period optimization means for calculating the significant period T,
The expansion coefficient computing means calculates an approximate expansion coefficient obtained by performing Fourier integration of the approximate Gaussian kernel g _r (x) in the integration region [−T / 2, T / 2], and the expansion coefficients {a ₀ ,. _K }, and
The significant period optimization means includes an equation of a variable T with respect to the scale parameter σ _r , the truncation term number K, and the dynamic range R.
9. The Gaussian kernel arithmetic apparatus according to claim 8, wherein T satisfying the condition is calculated as the significant period T. The approximate parameter optimization means includes:
Based on the scale parameter σ _{r, the} range upper limit value R, and a preset effective width parameter α (ασ _r <R), the significant period T is
A calculation cycle setting means calculated by:
Based on the allowable error τ, the significant period T, and the scale parameter σ _r , the truncation number K is
The Gaussian kernel arithmetic device according to claim 7, further comprising: a truncation term number determining means for calculating by the following. The expansion coefficient computing means calculates the expansion coefficients {a ₁ ,..., A _K },
The Gaussian kernel arithmetic device according to claim 7, wherein the Gaussian kernel arithmetic device is calculated by:   A program that causes a computer to function as the Gaussian kernel arithmetic device according to any one of claims 7 to 11 by being read and executed by a computer.