JP5620829B2 - Data processing apparatus, image processing apparatus, and methods thereof - Google Patents

Description

  The present invention relates to digital data filtering.

  As a method for speeding up filter processing with a large filter size, for example, a method of limiting a frequency region for performing filter processing is known (for example, Patent Document 1). Patent Document 1 discloses a technique for correcting the defocusing of a bifocal lens. In other words, when convolution is performed between the inverse function of the function representing blur and the captured image, the low frequency component of the captured image is extracted by band division and convolved, considering that the focus blur exists in the low frequency component. Reduce the amount of computation.

  Also, considering the convolution operation by the filter as the interaction between the target pixel to be filtered and its surrounding pixels, it is a multi-body that solves a field where some force interaction occurs between many particles existing in space. It can be regarded as a problem (N-body problem). A fast multipole method is known as a method for speeding up many-body problems (for example, Non-Patent Document 1). The fast multipole method pays attention to the fact that the acting force between distant objects is small in the field where gravity or Coulomb force interacts. This is a method for calculating efficiently together.

  However, the technique disclosed in Patent Document 1 can filter only a specific band, but cannot be applied to the case of filtering the entire frequency band.

  In addition, the technology disclosed in Non-Patent Document 1 reduces the amount of calculation by transforming a force that is inversely proportional to the square of the distance between two points, such as gravity and Coulomb force, but the interaction force Is not applicable if is represented by any more general function.

JP 2007-72558

L. Greengard and V. Rohklin `` A Fast Algorithm for Particle Simulations '' JOURNAL OF COMPUTATIONAL PHYSICS 73, 315-348 (1987)

  An object of the present invention is to reduce the amount of calculation in the filtering process of digital data and to speed up the filtering process.

  The present invention has the following configuration as one means for achieving the above object.

  The data processing according to the present invention divides input data into a plurality of divided data of a predetermined data size, divides a filter into a plurality of divided filters based on the predetermined data size, and sets a filter coefficient for each of the plurality of divided filters. An approximate polynomial is created, and a filter process using the filter is performed on the input data using a plurality of divided data, the plurality of divided filters, and the approximate polynomial.

  According to the present invention, it is possible to reduce the amount of calculation in the filtering process of digital data and to speed up the filtering process.

1 is a block diagram illustrating a configuration example of an image processing apparatus according to an embodiment. 6 is a flowchart for describing filter processing executed by the image processing apparatus. The flowchart explaining the detail of a filter process. The figure explaining division | segmentation of filter size. The figure explaining the relationship between the image data A and filter B 'in a convolution calculation. The figure which shows the enlarged view around division area Rdi and Rdj. The figure explaining expansion of a filter range. The figure which shows an example of the image A used for confirmation of a processing result. The figure which shows an example of the image A used for confirmation of a processing result. The figure which shows the image C acquired by performing the process by filter B ', and the image C acquired by performing the process by the filter B with respect to the image A. The figure explaining the comparison result of the process which uses filter B, and the process which uses filter B '. The figure which shows the image and spatial frequency distribution after filter processing by making the size of a division area into 20x20 pixels. The flowchart explaining a one-dimensional filter process. The flowchart explaining the detail of a filter process.

  Hereinafter, data processing and image processing according to an embodiment of the present invention will be described in detail with reference to the drawings.

[Device configuration]
A configuration example of the image processing apparatus according to the embodiment will be described with reference to the block diagram of FIG.

  The CPU 101 uses the RAM 103 as a work memory, executes a program stored in the ROM 102 or the hard disk drive (HDD) 108, and controls each component to be described later via the system bus 104. Thereby, various processes described later are executed.

  A general-purpose interface (I / F) 105 is a serial bus interface such as USB or IEEE1394 that connects various devices 106 such as a keyboard, mouse, digital camera, scanner, memory card reader, and hard disk drive. The CPU 101 can read various data from the device 106 via the general-purpose I / F 105 and write various data to the device 106.

  The HDD interface (I / F) 107 is an interface such as a serial ATA (SATA) for connecting a secondary storage device such as the HDD 108 or an optical disk drive. The CPU 101 can read data from the HDD 108 and write data to the HDD 108 via the HDD I / F 107. Further, the CPU 101 can expand the data stored in the HDD 108 to the RAM 103, and can similarly store the data expanded in the RAM 103 to the HDD 108. The CPU 101 can execute the data expanded in the RAM 103 as a program.

  A video card (VC) 109 is an interface for the monitor 110. The CPU 101 displays various images including a user interface (UI), which will be described later, an image to be filtered, and an image after filtering, on the monitor 110 via the VC 109.

  The user operates the UI displayed on the monitor 110 with a keyboard or mouse connected to the general-purpose I / F 105 as the device 106, and instructs the CPU 101 to execute, for example, a filter process described later. The CPU 101 instructed to execute the filter process executes a filter process program stored in the ROM 102 or the HDD 108. Then, the CPU 101 reads image data specified by the user via the UI from the HDD 108 or the device 106, and performs a filtering process on the read image data. Further, the image data after the filtering process is stored in the area of the HDD 108 or the device 106 designated by the user via the UI, or the image represented by the image data after the filtering process is displayed on the monitor 110.

[Filter processing]
The filtering process executed by the image processing apparatus will be described with reference to the flowchart of FIG. 2A.

  The CPU 101 inputs image data A designated by the user (S101), and divides the image represented by the input image data A into a plurality of regions of horizontal Nx pixels × vertical Ny pixels (predetermined size) (S102). Note that, when an image is divided in the horizontal direction (X direction), an area having a width less than Nx pixels may occur at the right end of the image, for example. Such an area is also treated as a divided area. Similarly, when an image is divided in the vertical direction (Y direction), an area having a height less than Ny pixels may occur at the lower end of the image, for example. Such an area is also treated as a divided area. In the following description, the number of divided areas is M.

  Next, the CPU 101 inputs data representing the filter B for filter processing from, for example, the HDD 108 (S103). The filter size of filter B is horizontal Nfx pixels × vertical Nfy pixels. The filter B may be expressed as a function of the distance from the target pixel, or may be data representing each component (filter coefficient) of the filter. However, as will be described later, in order to increase the filter size of the filter B to a natural number times the size of the divided region, it is preferable that the component (filter coefficient) of the portion corresponding to the tail of the filter processing is substantially zero.

Next, the CPU 101 divides the filter B so as to correspond to the size of the divided area Nx × Ny (S104). At this time, the CPU 101 enlarges the filter size so that the divided filter is composed of an odd number of divided filters (size Nx × Ny) both vertically and horizontally. The filter size division will be described with reference to FIG. That is, assuming that the number of divided filters is vertical Lx × horizontal Ly, the filter size is expanded so that Lx and Ly satisfy the following expression. Hereinafter, the number of division filters is expressed as Lx × Ly = LL.
Nx × Lx ≧ Nfx
Ny × Ly ≧ Nfy… (1)
Here, Lx and Ly are odd numbers.

The size of the enlarged filter is expressed by the following formula.
Nfx '× Nfy' = (Nx × Lx) × (Ny × Ly)… (2)

  Next, for each divided filter, the CPU 101 approximates the filter coefficient with a polynomial function of the distance from the filter center to obtain a coefficient of each term of the polynomial (S105). Note that a filter in which the original filter B is represented by an approximate polynomial is a filter B ′. Then, as will be described in detail later, the CPU 101 performs filter processing using the filter B ′ on the image data A using a plurality of divided regions, a plurality of divided filters, and an approximate polynomial (S106).

● Approximation polynomial Here, the creation of the approximation polynomial is explained.

  When the filter is expressed by a function capable of second order differentiation, an approximate polynomial may be created by performing Taylor expansion around the center of each divided filter. More generally, a least square method is used in which the coefficients of the approximate polynomial are determined so as to minimize the sum of squared errors. The least squares method can be applied when the filter function includes discontinuous or non-differentiable points, or when the filter coefficient is given in advance for each cell corresponding to a pixel instead of a function of the distance from the center. It is.

Here, an example of obtaining an approximate polynomial using Taylor expansion will be described. For example, it is assumed that the filter coefficient of the filter B is expressed as a function f (x, y) of the distance (x, y) from the filter center. When the function f (x, y) is a function capable of second order differentiation and Taylor expansion is performed up to a quadratic term around the position (a, b), the following equation is obtained.
f (x, y) ≒ f (a, b) + f _x (a, b) (xa) + f _y (a, b) (yb)
_{+ F xx (a, b)} (xa) 2/2 + f yy (a, b) (yb) 2/2
+ f _xy (a, b) (xa) (yb)… (3)
Where f _x = ∂f / ∂x, f _y = ∂f / ∂y,
f _xx = ∂ ² f / ∂x ² , f _yy = ∂ ² f / ∂y ² ,
f _xy = ∂ ² f / ∂x / ∂y

Furthermore, when the equation (3) is modified by defining the following, the equation (4) is obtained.
a _xx = f _xx (a, b) / 2, a _yy = f _yy (a, b) / 2, a _xy = f _xy (a, b),
a _x = f _x (a, b)-2af _xx (a, b)-bf _xy (a, b), a _y = f _y (a, b)-2bf _yy (a, b)-af _xy (a , b),
_{a 0 = f (a, b} ) - af x (a, b) - bf y (a, b) + a 2 f xx (a, b) / 2 + abf xy (a, b) + b 2 f yy (a, b) / 2,
f (x, y) ≒ a _xx x ² + a _xy xy + a _yy y ² + a _x x + a _y y + a ₀ … (4)

In addition, although the example of the 2nd order polynomial was shown above, it is not limited to this, The polynomial which has the order of 1st order or 3rd order or 4th order or more may be used. For example, in the case of an approximate polynomial up to a first order term, equation (5) is used instead of equation (4). If Expression (5) is used, the number of terms is halved, but the approximation accuracy (the error between the original function f (x, y) and the polynomial shown in Expressions (3) and (4)) decreases.
f (x, y) ≒ a _x x + a _y y + a ₀ … (5)

On the other hand, if an approximation polynomial having a degree of 3rd order or 4th order or more is used, the approximation accuracy is improved. However, for example, the number of terms increases to 10 for the third order and 15 for the fourth order, so that the calculation amount increases. In general, when f (x, y) is approximated by an nth order polynomial, it is expressed as equation (6).
f (x, y) ≒ Σ _p Σ _q a _pq x ^p q ^j … (6)
Where Σ _p ranges from p = 0 to n,
Σ _q ranges from q = 0 to n,
n is the degree of the polynomial.

When filter B is divided into LL divided filters, LL approximation polynomials of Equation (6) are also generated, and LL sets of coefficients {a _pq } are also obtained. Note that it is not necessary for all the LL approximation polynomials to have the same order. For example, a division filter close to the filter center may be a second-order approximation polynomial, and a division filter far from the filter center may be a first-order approximation polynomial.

Filter Process Details of the filter process (S106) will be described with reference to the flowchart of FIG. 2B.

  After creation of the approximate polynomial (S105), the CPU 101 performs a convolution operation of each pixel and the filter on the image data A obtained by dividing the region using the approximate polynomial of the filter B ′. The CPU 101 initializes the counter i to 1 (S111), and selects the i-th divided area Rdi (first divided area) of the image data A (S112).

  The relationship between the image data A and the filter B ′ in the convolution operation will be described with reference to FIG. As shown in FIG. 4, the CPU 101 superimposes the filter B ′ on the image data A so that the filter center matches the divided region Rdi, and determines the target region Rf of the convolution operation by the divided region Rdi and the filter B ′ ( S113). Then, the counter j is initialized to 1 (S114), and the jth divided region Rdj (second divided region) included in the target region Rf is selected (S115).

  FIG. 4 shows a state in which the divided region Rdj corresponding to the counter j is selected in the target region Rf centered on the divided region Rdi corresponding to the counter i. The CPU 101 selects a set of coefficients of the approximate polynomial of the filter B ′ based on the position of the divided area Rdj with respect to the position of the divided area Rdi (S116).

When the pixel (x1, y1) in the divided region Rdi is the target pixel, the convolution operation of the target pixel and the filter B ′ is expressed as follows.
h '(x1, y1) = h (x1, y1) * f (x1, y1)
= Σ _{x '} Σ _y' h (x ', y') f (x'-x1, y'-y1)… (7)
Where h (x1, y1) is the pixel value of the target pixel,
h ′ (x1, y1) is the pixel value after filtering,
f (x1, y1) is the filter coefficient,
* Represents a convolution operation.

In Equation (7), the contribution from the pixels included in the divided region Rdj is assumed that the filter B ′ is approximated by a quadratic polynomial (4), and the position of the pixel included in the divided region Rdj is (x2, y2) is expressed by the following equation.
h '(x1, y1) | _Rdj = ΣΣh (x2, y2) f (x2-x1, y2-y1)
≒ ΣΣh (x2, y2) {a _xx (x2-x1) ² + a _xy (x2-x1) (y2-y1)
+ a _yy (y2-y1) ² + a _x (x2-x1) + a _y (y2-y1) + a ₀ }
= ΣΣh (x2, y2) {a _xx x1 ² + a _xy x1y1 + a _yy y1 ²
+ (-2a _xx x2-a _xy y2-a _x ) x1 + (-2a _yy y2-a _xy x2-a _y ) y1
+ a _xx x2 ² + a _xy x2y2 + a _yy y2 ² + a _x x2 + a _y y2 + a ₀ }… (8)
Here, the range of Σ is the range of pixels (x2, y2) included in the divided region Rdj.

Further, when the equation (8) is modified by defining the following, the equation (9) is obtained.
G _xx = ΣΣh (x2, y2) x2 ² , G _xy = ΣΣh (x2, y2) x2y2, G _yy = ΣΣh (x2, y2) y2 ² ,
_{G x = ΣΣh (x2, y2} ) x2, G y = ΣΣh (x2, y2) y2, G 0 = ΣΣh (x2, y2).
h '(x1, y1) | _Rdj ≒ a _xx G ₀ x1 ² + a _xy G ₀ x1y1 + a _yy G ₀ y1 ²
+ (-2a _xx G _x -a _xy G _y -a _x G ₀ ) x1 + (-2a _yy G _y -a _xy G _x -a _y G ₀ ) y1
+ (a _xx G _xx + a _xy G _xy + a _yy G _yy + a _x G _x + a _y G _y + a ₀ G ₀ )
= φ _xx x1 ² + φ _xy x1y1 + φ _yy y1 ² + φ _x x1 + φ _y y1 + φ ₀ … (9)
Where φ _xx = a _xx G ₀ , φ _xy = a _xy G ₀ , φ _yy = a _yy G ₀ ,
φ _x = -2a _xx G _x -a _xy G _y -a _x G ₀ , φ _y = -2a _yy G _y -a _xy G _x -a _y G ₀ ,
φ ₀ = a _xx G _xx + a _xy G _xy + a _yy G _yy + a _x G _x + a _y G _y + a ₀ G ₀

Expression (9) is an expression obtained by approximating the filter B ′ with the above-described second-order polynomial (4). An expression that approximates the filter B ′ with the above-described more general n-order polynomial (6) is as follows.
h '(x1, y1) | _Rdj = Σ _x2 Σ _y2 h (x2, y2) f (x2-x1, y2-y1)
≒ Σ _x2 Σ _y2 h (x2, y2) {Σ _p Σ _q a _pq (x2-x1) ^p (y2-y1) ^q }… (10)

When Expression (10) is rearranged for x1 and y1 as in Expression (9), the following expression is obtained.
h '(x1, y1) | _Rdj ≒ Σ _p Σ _q φ _pq x1 ^p y1 ^q … (11)

{Φ _pq } (0 ≦ p ≦ n, 0 ≦ q ≦ n) in Equation (11) is a coefficient of the x1 ^p y1 ^q term, the coefficient {a _pq } of the approximate polynomial (6), and the divided region Rdj It is obtained only from the position of the pixel in FIG. 5 shows an enlarged view around the divided areas Rdi and Rdj. As shown in FIG. 5, the pixel (x3, y3) different from the pixel (x1, y1) in the divided region Rdi is set as the target pixel, and from the divided region Rdj in the convolution operation of the target pixel (x3, y3) and the filter B ′ Is considered, the following equation is obtained in the same manner as equation (11).
h '(x3, y3) | _Rdj ≒ Σ _p Σ _q φ _pq x3 ^p y3 ^q … (12)

The value of {φ _pq } (0 ≦ p ≦ n, 0 ≦ q ≦ n) in Equation (12) is the value of the pixel (x1, y1) or pixel (x3, y3) of the divided region Rdi as described above. And has the same value as {φ _pq } in equation (11). That is, for each pixel included in the divided region Rdi, the calculation of the contribution from the divided region Rdj in the convolution operation with the filter B ′ may be performed according to the following procedure.
Select a set of coefficients of the approximate polynomial of the filter B ′ from the positional relationship between the divided region Rdi and the divided region Rdj (S116),
-Calculate the coefficient {φ _pq } of the divided region Rdj, save the obtained value (S117),
For each pixel in the divided region Rdi, the contribution from the divided region Rdj is calculated using the coefficient {φ _pq } and Expression (11) (S118).

  Next, the CPU 101 increments the counter j (S119), compares the count value of the counter j with the number LL of divided filters (S120), and if j ≦ LL, returns the process to step S115. If j> LL, the counter i is incremented (S121), the count value of the counter j is compared with the number M of the divided areas Rd (S122), and if i ≦ M, the process returns to step S112.

  When j> M, the CPU 101 adds the result of the convolution operation calculated by Expression (11) for each pixel, and the image of the convolution operation result (filter processing result) of the image data A and the filter B ′ Output as data A '(S123).

[Calculation amount]
When the filter is not divided, the number of calculation times of the filter process for one pixel of interest is Nfx × Nfy times. Therefore, the number of calculation times NR of all the pixels included in one divided region Rd is expressed by the following equation.
NR = Nfx / Nfy / Nx / Ny (13)

  One calculation includes a process of multiplying a pixel value included in the divided region Rd by a filter coefficient and integrating the multiplication result as a new pixel value. Further, the amount of calculation in the filter processing of this embodiment is estimated as follows.

Calculation of coefficient {φ _pq } of divided region Rdj (S117)
The amount of calculation is proportional to the product of the number of pixels in the divided area and the number of divided filters in the filter B ′.
α1 ・ (Nx × Ny) (Nfx '/ Nx) (Nfy' / Ny) = α1 ・ Nfx '・ Nfy'… (14)
Here, α1 is a constant indicating the amount of calculation per time.

● Calculation of contribution from divided area Rdj (S118)
That is, the calculation of the equation (11) using the coefficient {φ _pq }, the calculation of the equation (11) is repeated by the number of divided filters included in the filter B ′ per pixel, and the pixels of the divided region are further calculated. Repeat for a few minutes.
α2 ・ (Nfx '/ Nx) (Nfy' / Ny) ・ Nx ・ Ny = α2 ・ Nfx '・ Nfy'… (15)
Here, α2 is a constant indicating the amount of calculation per expression (11).

Therefore, the result of adding Expressions (14) and (15) is the calculation amount LF required for the convolution operation of the divided region and the filter B ′ in the filter processing of the present embodiment.
LF = α ・ Nfx '・ Nfy'… (16)
Where α = α1 ・ α2

Further, the following expression is obtained from the number of calculations (expression (13)) and expression (16) when the filter is not divided.
α ・ Nfx '・ Nfy' / (Nfx ・ Nfy ・ Nx ・ Ny) = α '/ (Nx ・ Ny)… (17)
Where α '= α ・ (Nfx' / Nfx) ・ (Nfy '/ Nfy)

  Since the filter B ′ is composed of an odd number of divided filters, the size of the filter B ′ is larger than the size of the filter B. Α ′ is taken into account that the increase in the filter size affects the calculation amount. For example, if α ′ is about 9.0 and Nx = Ny = 8, the filter processing of this embodiment can reduce the calculation amount to 9 / (8 × 8) ≈1 / 7.

The reason why the amount of calculation is reduced The reason why the amount of calculation can be reduced as shown in equation (17) is to replace equation (7) or equation (8) of convolution operation with the approximate equation shown in equation (11). There is.

Here, consider a case where a convolution operation is performed on a filter whose filter coefficient is determined by the distance from the filter center without dividing the filter. When the target pixel (x1, y1) to be filtered changes, the distance from the filter center of the peripheral pixels (x2, y2) included in the filter range changes, and the contribution of the peripheral pixels (x2, y2) also changes. Therefore, the contribution C from the peripheral pixel (x2, y2) is expressed by the following equation.
C = h (x2, y2) × f (x2-x1, y2-y1)… (18)
Where h (x2, y2) is the pixel value of the surrounding pixel (x2, y2),
f (x, y) is the filter function,
(x1, y1) is the position of the target pixel.

When the target pixel moves to (x3, y3), the contribution C from the peripheral pixel (x2, y2) is expressed by the following equation.
C = h (x2, y2) × f (x2-x3, y2-y3)… (19)

  That is, even if the contribution C from the same surrounding pixel (x2, y2) is calculated, the value of the filter function f (x, y) changes as the pixel of interest moves, so the contribution C is calculated for each pixel of interest. There is a need to calculate.

On the other hand, in Expression (11), {φ _pq } is determined only by the pixel position and pixel value in the divided region Rdj, and is separated from the position (x1, y1) of the target pixel in the divided region Rdi to be filtered. . Therefore, even when the target pixel moves, {φ _pq } representing the contribution from the pixel in the divided region Rdj is the same (see equations (11) and (12)), and recalculation of {φ _pq } is unnecessary. Thus, {φ _pq } can be regarded as a constant. Also, the expression method shown in equation (11) is replaced with a polynomial function as shown in the modified example of equations (equations (8) and (9)) when the original filter is approximated by a quadratic polynomial (4) This can be easily realized.

  In general, an arbitrary continuous function can be approximated by a polynomial. However, if the error between the original function and the approximate function is reduced, the degree of the polynomial needs to be increased and the amount of calculation increases. Therefore, if the original filter range is divided into regions, and an approximate polynomial in one divided filter range is obtained, an approximation with sufficiently small error within the range of the divided filter can be obtained without increasing the order of the polynomial. A polynomial can be obtained. Furthermore, if the method of determining the shape of the approximate polynomial in advance and obtaining the coefficient of the approximate polynomial using the least square method is used, the original filter does not necessarily have to be a continuous function.

  As for the computational complexity shown in equations (13) to (17), the computational complexity is the sum of “one multiplication” and “one addition” as described above when the filter is not divided. It is. On the other hand, the amount of calculation in this embodiment includes, for example, five additions and multiplications in the equation (9) using a quadratic approximate polynomial, for example. In addition, there is an effect that the size of the filter is expanded from Nfx × Nfy to Nfx ′ × Nfy ′ so that the size of the filter becomes a division filter unit. Therefore, the constant α in the equation (17) is larger than 1, and becomes about 9.0 when a quadratic approximate polynomial is used, for example.

[Expand filter range]
Since the filter processing in the present embodiment is performed in units of an odd number of divided filters, the range of the filter B ′ is expanded as compared with the range of the filter B. The expansion of the filter range will be described with reference to FIG.

  FIG. 6 shows an image A obtained by dividing a 72 × 80 pixel image indicated by a solid line into 8 × 8 pixels, and a filter B having a size indicated by a broken line of 49 × 49 pixels divided into 8 × 8 pixel filters. A state in which filters B ′ having the divided filters are overlapped is shown. In this case, the size of the filter B ′ is 56 × 56 pixels.

  When the upper left pixel (32, 40) in the divided area shown in FIG. 6 is the target pixel 501, the range of the filter B (49 × 49 pixels) is the upper left (8, 16) and lower right (56, 64) indicated by dotted lines Range. When the lower right pixel (39, 47) is the target pixel 502, the range of the filter B is the upper left (15, 23) and lower right (63, 71) indicated by the dotted line. On the other hand, in the present embodiment, the range of the filter B ′ does not change with respect to the target pixels 501 and 502 in the same divided region, and the upper left (8, 16) and lower right (63, 71) ranges indicated by broken lines It is. In this embodiment, for the pixel of interest 501, a pixel to the right of x = 57 and a pixel below y = 65 is added to the reference range, and for the pixel of interest 502, a pixel to the left of x = 14 and a pixel above y = 22 is referenced. Added to the pixel.

  Such a difference in the filter range is an error that occurs due to the execution of the convolution operation for each divided filter. If the filter coefficient of the portion where the filter range is expanded is almost 0, the influence of the error is slight. Therefore, it is most preferable to apply the present embodiment to a filter in which the coefficient at the end of the filter range far from the filter center is almost zero.

[Processing result]
● Gaussian filter Figure 7 shows an example of image A used to check the processing results. An image A shown in FIG. 7 is an 8-bit monochrome image of 2048 × 2048 pixels, and a random square and a circle are superimposed at random positions, and further random noise is added. The filter B is a Gaussian filter expressed by the following equation.
f (x, y) = C1 × exp {-(x ² + y ² ) / 2 / σ ² }… (20)
Where C1 is a constant.

  In the filter B (Expression (20)), σ is 12 pixels, and the filter size is 73 × 73 pixels. The constant C1 is set to 0.00249 so that the sum of the coefficients f (x, y) of the filter B becomes 1.0. The size of the divided area is 8 × 8 pixels. Each divided filter of the filter B ′ is approximated by a second order polynomial expressed by the equation (4), and the size of the filter B is an 11 × 11 divided filter. That is, the filter B is enlarged from 73 × 73 pixels to 88 × 88 pixels to be a filter B ′.

For image A, processing by filter B was performed to obtain image C, and processing by filter B ′ was performed to obtain image D. Then, the difference between the pixel values for each pixel of the image C and the image D was taken, and the maximum and minimum values of the difference values in the entire image area were obtained. Further, the root mean square (RMS) of the difference value of each pixel was calculated by the following equation.
RMS = √ [Σ _x Σ _y {h _C '(x, y)-h _D ' (x, y)} ² } / N _A ]… (21)
Where h _C (x, y) is the pixel value of image C,
h _D (x, y) is the pixel value of image D,
N _A is the number of pixels.

  As a result, the maximum difference value in all image regions was 1, the minimum was -1, the pixel value difference was within ± 1, and the RMS value was 0.196. Furthermore, when the distribution of the spatial frequency component was examined by applying Fast Fourier Transform (FFT) to the distribution of the image D and the difference value, no spatial frequency characteristic related to the size of the divided region was found. . Therefore, it was confirmed that the image C and the image D are almost the same image, and the difference due to the difference in filter processing hardly occurs.

  Further, the processing time by the filter B ′ was 1 / 7.5 of the processing time by the filter B.

● LoG filter Fig. 8 shows an example of image A used to check the processing results. An image A shown in FIG. 8 is an 8-bit monochrome image of 1024 × 1024 pixels, and squares and circles of a random size are superimposed at random positions, and further random noise is added. The filter B is a LoG (Laplacian of Gaussian) filter represented by the following equation.
f (x, y) = ( x 2 + y 2 - 2σ 2) / (2πσ 6) × exp {- (x 2 + y 2) / 2 / σ 2} ... (22)

  In the filter B (Equation (22)), σ is 24 pixels, and the filter size is 145 × 145 pixels. The size of the divided area is 8 × 8 pixels. Each division filter of the filter B ′ is approximated by a second order polynomial expressed by the equation (4), and the size of the filter B is a 19 × 19 division filter. That is, the filter B is enlarged from 145 × 145 pixels to 152 × 152 pixels to be a filter B ′.

  FIG. 9 shows an image C (FIG. 9 (a)) acquired by performing processing using the filter B on the image A, and an image D (FIG. 9 (b)) acquired by performing processing using the filter B ′. In FIG. 8, positive value pixels obtained by the filter processing are represented in black, and negative value pixels are represented in white. Since the σ value in equation (22) is large, the width of the line formed by the positive pixels in FIG. 9 is thick, but it is clear that if it is made thin, it corresponds to the edge of the figure shown in FIG. . Similarly, edges are extracted in both FIGS. 9A and 9B.

  Further, the processing time by the filter B ′ was 1 / 6.2 of the processing time by the filter B.

Changing the size of the divided area In the above, the example where the size of the divided area is 8 x 8 pixels has been explained. However, in the following, when the size of the divided area is changed, the processing time and the filter are not divided A description will be given of how the error changes.

  Image A is an 2048 × 2048 pixel 8-bit monochrome image shown in FIG. The filter B is a Gaussian filter (73 × 73 pixels) shown in Expression (20). Then, the size of the divided area is 8 types of 6 × 6, 8 × 8, 10 × 10, 12 × 12, 14 × 14, 16 × 16, 18 × 18, and 20 × 20 pixels. While changing the size of the divided areas, the accuracy and processing time of the process using the filter B and the process using the filter B ′ were compared.

  A comparison result between the process using the filter B and the process using the filter B ′ will be described with reference to FIG. In FIG. 10, the horizontal axis is “division area side size / filter σ value”, the left vertical axis is the processing time ratio (= processing using filter B ′ / processing using filter B), and right vertical axis. Represents the RMS value shown in Formula (21), respectively.

  As can be seen from FIG. 10, in the process using the filter B ′, the processing time is shortened as the size of the divided region is increased, and the error is reduced as the size of the divided region is decreased.

  FIG. 11 shows an image (FIG. 11 (a)) and a spatial frequency distribution (FIG. 11 (b)) after the filter processing with the size of the divided region being 20 × 20 pixels. Note that the image shown in FIG. 11 (a) is a partially enlarged view so that a striped pattern with 20 pixel intervals in the divided region can be seen.

  Referring to FIG. 11 (a), patterns with pitches corresponding to ± 20 pixels and ± 20 pixels / integer are generated with respect to the XY axes, respectively. This is presumably because the divided area is larger than the size of the original filter, so that the accuracy of the approximate polynomial is lowered and an error occurs.

  As described above, the filter processing of this embodiment for dividing the filter can be expected to shorten the processing time as the size of the divided region is increased. On the other hand, the error is increased and the accuracy is lowered. Therefore, by appropriately determining the size of the divided region, it is possible to shorten the processing time of the filter processing without reducing the accuracy.

  The image processing according to the second embodiment of the present invention will be described below. Note that the same reference numerals in the second embodiment denote the same parts as in the first embodiment, and a detailed description thereof will be omitted.

  In the first embodiment, the two-dimensional filter processing has been described. In the second embodiment, the one-dimensional filter processing will be described with reference to the flowchart of FIG. 12A.

  The CPU 101 inputs one-dimensional data A (for example, sound data) (S201), and generates divided data obtained by dividing the input data A into Nx data (predetermined data size) (S202). When the data A is two-dimensional data such as image data, for example, there is a possibility that divided data less than Nx pixels may be generated at the right end of the image. Such a case is also treated as divided data. In the following description, the number of divided data is M.

  Next, the CPU 101 inputs data representing the filter B for one-dimensional filter processing from, for example, the HDD 108 (S203). Let the filter size of filter B be Nfx data. The filter B may be expressed as a function of the distance from the data of interest, or data representing each component (filter coefficient) of the filter. However, in order to increase the filter size of the filter B to a natural number multiple of the number of pieces of divided data, it is preferable that the component (filter coefficient) of the portion corresponding to the tail of the filter process is substantially zero.

Next, the CPU 101 divides the filter B in correspondence with the number of divided data Nx (S204). At that time, the CPU 101 enlarges the filter size so that the divided filter is composed of an odd number of divided filters (size Nx). That is, if the number of division filters is Lx, the filter size is expanded so that Lx satisfies the following expression.
Nx × Lx ≧ Nfx (23)
Where Lx is an odd number.

The size of the enlarged filter is expressed by the following formula.
Nfx '= Nx × Lx (24)

  Next, for each divided filter, the CPU 101 approximates the filter coefficient with a polynomial function of the distance from the filter center to obtain the coefficient of each term of the polynomial (S205). Note that a filter in which the original filter B is represented by an approximate polynomial is a filter B ′. Then, as will be described in detail later, the CPU 101 performs a filter process using the filter B on the data A using the plurality of divided data, the plurality of divided filters, and the approximate polynomial (S206).

Approximation polynomial When filter B is expressed as a function capable of second order differentiation, for example, Taylor expansion is performed on the center of each divided filter to obtain an approximation polynomial. For example, it is assumed that the filter coefficient of the filter B is expressed as a function f (x) of the distance x from the filter center. When the function f (x) is a function capable of second order differentiation and Taylor expansion is performed up to a second-order term with respect to the position a, the following expression is obtained.
f (x) ≒ f (a) + f _x (a) (xa) + f _xx (a) (xa) ² /2… (25)
Where f _x = df / dx,
f _xx = d ² f / dx ²

Further, the following equation is defined by defining the following and obtaining the equation (26).
a _xx = f _xx (a) / 2,
a _x = f _x (a)-2af _xx (a),
a ₀ = f (a)-af _x (a)-a ² f _xx (a) / 2,
f (x, y) ≒ a _xx x ² + a _x x + a ₀ … (26)

In addition, although the example of the 2nd order polynomial was shown above, it is not limited to this, The polynomial which has the order of 1st order or 3rd order or 4th order or more may be used. In general, when f (x) is approximated by an nth order polynomial, it is expressed as shown in Expression (27).
f (x) ≒ Σ _p a _p x ^p (27)
Where Σ _p ranges from p = 0 to n,
n is the degree of the polynomial.

When filter B is divided into Lx divided filters, Lx approximate polynomials of Expression (27) are also generated, and Lx sets of coefficients {a _p } are also obtained. Note that it is not necessary for all the Lx approximate polynomials to have the same order. For example, a division filter close to the filter center may be a second-order approximation polynomial, and a division filter far from the filter center may be a first-order approximation polynomial.

Filter Process Details of the filter process (S206) will be described with reference to the flowchart of FIG. 12B.

  After creating the approximate polynomial, the CPU 101 performs a convolution operation of each data and filter on the divided data A using the approximate polynomial of the filter B ′. The CPU 101 initializes the counter i to 1 (S211), and selects the i-th divided data Ddi (first divided data) of the data A (S212).

  The CPU 101 superimposes the filter B ′ on the data A so that the divided data Ddi and the filter center coincide with each other, and determines the target data set Sf of the convolution operation using the divided data Ddi and the filter B ′ (S213). Then, the counter j is initialized to 1 (S214), and the j-th divided data Ddj (second divided data) included in the target data set Sf is selected (S215).

The CPU 101 selects a set of coefficients of the approximate polynomial of the filter B ′ based on the position of the divided data Ddj with respect to the position of the divided data Ddi (S216). Assuming that the data x1 in the divided data Ddi is the attention data, the convolution operation of the attention data and the filter B ′ is expressed as follows.
h '(x1) = h (x1) * f (x1)
= Σ _{x '} h (x') f (x'-x1)… (28)
Where h (x1) is the value of the data of interest,
h '(x1) is the value after filtering,
f (x1) is the filter coefficient,
* Represents a convolution operation.

In Equation (28), the contribution from the data included in the divided data Ddj is assumed that the filter B ′ is approximated by a second order polynomial (26), and the position of the data included in the divided data Ddj is x2. It is expressed by the following formula.
h '(x1) | _Ddj = Σh (x2) f (x2-x1)
= Σh (x2) {a _xx (x2-x1) ² + a _x (x2-x1) + a ₀ }
= Σh (x2) {a _xx x1 ² + (-2a _xx x2-a _x ) x1 + a _xx x2 ² + a _x x2 + a ₀ }… (29)
Here, the range of Σ is the range of data x2 included in the divided data Ddj.

Furthermore, the equation (30) is obtained by modifying the equation (29) by defining the following.
G _xx = Σh (x2, y2) x2 ² , G _x = Σh (x2, y2) x2, G ₀ = Σh (x2, y2)
h '(x1) | _Ddj ≒ a _xx G ₀ x1 ² + (-2a _xx G _x -a _x G ₀ ) x1 + (a _xx G _xx + a _x G _x + a ₀ G ₀ )
≒ φ _xx x1 ² + φ _x x1 + φ ₀ … (30)
Where φ _xx = a _xx G ₀ , φ _x = -2a _xx G _x -a _x G ₀ , φ ₀ = a _xx G _xx + a _x G _x + a ₀ G ₀

Expression (30) is an expression obtained by approximating the filter B ′ with the above-described second-order polynomial (26). An expression that approximates the filter B ′ with the above-described more general n-order polynomial (27) is as follows.
h '(x1) | _Ddj = Σ _x2 h (x2) f (x2-x1)
≒ Σ _x2 h (x2) {Σ _p a _p (x2-x1) ^p }… (31)

When Expression (31) is rearranged with respect to x1 similarly to Expression (30), the following expression is obtained.
h '(x1) | _Ddj ≒ Σ _p φ _p x1 ^p … (32)

{Φ _p } (0 ≦ p ≦ n) in equation (32) is the coefficient of the x1 ^p term, and is obtained only from the coefficient {a _p } of the approximate polynomial (26) and the position of the data in the divided data Ddj. It is done. That is, considering the data x3 different from the data x1 in the divided data Ddi as the data of interest, and considering the contribution from the divided data Ddj in the convolution operation of the data of interest x3 and the filter B ′, as in equation (32), The following equation is obtained.
h '(x3) | _Ddj ≒ Σ _p φ _p x3 ^p … (33)

As described above, the value of {φ _p } (0 ≦ p ≦ n) in Equation (33) is independent of the value of data x1 or data x3 of divided data Ddi, and {φ _{p in} Equation (31) } Has the same value. That is, for each pixel included in the divided data Ddi, the calculation of the contribution from the divided data Ddj in the convolution operation with the filter B ′ may be performed according to the following procedure.
Select a set of coefficients of the approximate polynomial of filter B ′ from the positional relationship between the divided data Ddi and the divided data Ddj (S216),
Calculate the coefficient {φ _p } of the divided data Ddj and store the obtained value (S217)
For each data of the divided data Ddi, the contribution from the divided data Ddj is calculated using the coefficient {φ _p } and the equation (32) (S218).

  Next, the CPU 101 increments the counter j (S219), compares the count value of the counter j with the number Lx of divided filters (S220), and if j ≦ Lx, returns the process to step S215. If j> Lx, the counter i is incremented (S221), the count value of the counter j is compared with the number M of the divided data Rd (S222), and if i ≦ M, the process returns to step S212.

  When j> M, the CPU 101 adds the result of the convolution operation calculated by the equation (32) for each data, and the data A of the convolution operation result (filter processing result) of the data A and the filter B ′ Output as' (S223).

[Calculation amount]
When the filter is not divided, the number of times of calculation of the filtering process for the data of interest is Nfx, and the number of calculations NR of all data included in one piece of divided data Dd is expressed by the following equation.
NR = Nfx · Nx (34)

One calculation includes the multiplication of the data included in the divided data Dd and the filter coefficient, and the process of integrating the multiplication results as new data. Further, as described in the first embodiment, the calculation amount LF in the filter processing of the present embodiment is estimated as follows.
LF = α ・ Nfx '… (35)

Further, the following expression is obtained from the number of calculations (expression (34)) and expression (35) when the filter is not divided.
α ・ Nfx '/ (Nfx ・ Nx) = α' / Nx… (36)
Where α '= α ・ (Nfx' / Nfx)

  Since the filter B ′ is composed of an odd number of divided filters, the size of the filter B ′ is larger than the size of the filter B. Α ′ is taken into account that the increase in the filter size affects the calculation amount. For example, if α ′ is about 4 and Nx = 8, the filter processing of this embodiment can reduce the calculation amount to 4/8 = 1/2.

[Other Examples]
The present invention can also be realized by executing the following processing. That is, software (program) that realizes the functions of the above-described embodiments is supplied to a system or apparatus via a network or various storage media, and a computer (or CPU, MPU, etc.) of the system or apparatus reads the program. It is a process to be executed.

Claims (9)

Data dividing means for dividing the input data into a plurality of divided data of a predetermined data size;
Filter dividing means for dividing the filter into a plurality of divided filters based on the predetermined data size;
Creating means for creating an approximate polynomial of filter coefficients for each of the plurality of divided filters;
A data processing apparatus comprising: processing means for performing filter processing on the input data using the plurality of divided data, the plurality of divided filters, and the approximate polynomial. The processing means includes means for selecting first divided data;
Means for determining a reference range of divided data centered on the first divided data;
Means for selecting, from the reference range, second divided data to be referred to when filtering the first divided data;
Means for calculating a coefficient of the approximate polynomial based on a positional relationship between the first divided data and the second divided data in the filtering process of the data of interest included in the first divided data;
2. The data processing apparatus according to claim 1, further comprising means for calculating data contributing to the data of interest from the second divided data using the coefficient and the approximate polynomial.   3. The data according to claim 1, wherein the filter dividing unit divides the filter into an odd number of divided filters having the predetermined data size by enlarging a filter size of the filter. Processing equipment. Image dividing means for dividing the image represented by the image data into a plurality of divided regions of a predetermined size;
Filter dividing means for dividing the two-dimensional filter into a plurality of divided filters based on the predetermined size;
Creating means for creating an approximate polynomial of filter coefficients for each of the plurality of divided filters;
An image processing apparatus comprising: processing means for performing filter processing on the image data using the two-dimensional filter using the plurality of divided regions, the plurality of divided filters, and the approximate polynomial. Means for selecting a first divided region;
Means for determining a reference range of a divided area centered on the first divided area;
Means for selecting, from the reference range, a second divided area to be referred to when filtering the first divided area;
Means for calculating a coefficient of the approximate polynomial based on a positional relationship between the first divided area and the second divided area in the filtering process of the target pixel included in the first divided area;
5. The image processing apparatus according to claim 4, further comprising means for calculating data contributing to the target pixel from the second divided region using the coefficient and the approximate polynomial. A data processing method for a data processing apparatus having data dividing means, filter dividing means, creating means, and processing means,
The data dividing means divides input data into a plurality of divided data of a predetermined data size,
The filter dividing means divides the filter into a plurality of divided filters based on the predetermined data size,
The creating means creates an approximate polynomial of filter coefficients for each of the plurality of divided filters,
The data processing method characterized in that the processing means performs filter processing by the filter on the input data using the plurality of divided data, the plurality of divided filters, and the approximate polynomial. An image processing method of an image processing apparatus having a data dividing means, a filter dividing means, a creating means, and a processing means,
The data dividing means divides the image represented by the image data into a plurality of divided areas of a predetermined size;
The filter dividing means divides the two-dimensional filter into a plurality of divided filters based on the predetermined data size,
The creating means creates an approximate polynomial of filter coefficients for each of the plurality of divided filters,
The image processing method characterized in that the processing means performs filter processing by the two-dimensional filter on the image data using the plurality of divided regions, the plurality of divided filters, and the approximate polynomial.   A program that causes a computer to function as each unit of the data processing device according to any one of claims 1 to 3.   6. A program for causing a computer to function as each unit of the image processing apparatus according to claim 4 or 5.