Methods of computer graphics which had already become common in design of new optical systems and materials find currently new applications such as stomatology and ophthalmology. Some modern imaging systems are now designed in conjunction with the human vision system which is at their end. As a result simulation of the effects of human vision becomes necessary. These include partial defocusing and resulting “blur” of image, scattering and halo/corona and so on. Such effects are usually simulated convolving the original, “ideal” image with the pixel spread function. The latter frequently has size about that of the source image, so straightforward calculation of convolution would take a giant number of operations. Therefore in case of high resolution a decent speed is usually achieved by using the Fast Fourier Transform (FFT) for convolution, but since FFT operates periodic functions on a lattice with resolution being an integer power of a prime numbers, the required working resolution may considerably increase that of the original image and required memory becomes inadmissible. This paper presents an alternative method that allows calculations in much smaller memory avoiding overheads introduced by FFT requirements.
Methods of computer graphics which had already become common in design of new
optical systems and materials, are now applied in stomathology and ophthalmology.
Some imaging systems are now designed in conjunction with the human vision system
which is at their end. Simulation of many effects of human vision, as well as defocusing
and aberration of lenses etc use 2D convolution of the original, “ideal” image with a
kernel for the pixel spread function [
Copyright © 2020 for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). ? Supported by RFBR, grants 18-01-00569, 18-31-20032 and 19-01-00435.
In simulation of such human vision effects as halo and glare [
Since the kernel depends only on the difference between the source and target
coordinates, efficiency can be drastically improved using Fast Fourier Transform (FFT) [
However FFT is not good regarding memory. First, FFT operates periodic functions (both input, kernel and the result) so that the boundaries of the image are glued. As a result, the FFT convolution to the black source image with the only white pixel at (0,0) produces filtered image which is gray not only near (0,0) but also near the opposite end (x = Nx; y = Ny). To avoid this artifact one has to add “margins” of size at least twice filter radius and pad the source image to 0 there.
Meanwhile, as said above, in human vision simulation that radius is about the image
size, thus we apply FFT convolution to the thrice larger image, 120000ı75000 pixels in
the above example. And this is not the end. The usual FFT is radix-2, i.e. the resolution
must be an integer power of two. Although there do exist radix-3, 5 and other variants
[
Since the number of FFT operations is nearly linear in the total number of pixels, it still remains more or less admissible, and modern PCs can do that in half an hour or hour for 3 color components. The really critical thing is though memory because 131072ı131072 pixels, times 3 color channels is already 192 Gb. One must keep both the source image, the filtered image and the kernel at a time, all of that resolution. Plus several auxiliary arrays etc. So for our example this totals to more that 512 Gb which is rarely available nowadays.
Meanwhile, the source image is just about 12 Gb so all this huge memory is merely the FFT overhead. Happily, for symmetric kernels it is possible to design an alternative method that works in much less memory while still has decent speed. This paper describes two such algorithms (that can be also combined). 2
This is closely related to the well-known Level Of Detail or mipmap [
The kernel for human vision is rather sophisticated [
K(r) = where r px2 + y2 is the distance between the source and target pixels. The spread a can be both small (less than a pixel) for some components, and large for other components. As a result, for small r the kernel changes very steeply, while at large r it inevitably changes more and more slowly.
To utilize that let us decompose our filter into the sum of the “central” and “peripheral” parts, so that the former be distinct from zero for only, say, r 50 while the “periphery”, which can be not 0 everywhere, be smooth (slowly changing) enough.
The convolution with the original filter is then the sum of convolutions with the “central” and “peripheral” sub-kernels. The former, because of the small filter radius, can be calculated even directly, without any FFT and thus in its original resolution. It requires about 1002 operations per pixel which is a lot but still admissible for modern PCs.
As to the smooth peripheral part, it is insensitive to the fine-scale detail of the source image. And the filtered image is smooth. So we can calculate it in decreased resolution (say tenfold) and then interpolate the result to the original resolution. And in this highly reduced resolution we can already apply FFT because the overhead is admissible.
Alternatively, we can repeat the procedure, now decomposing the peripheral component into again the “central part of periphery” and “periphery of periphery”. Again the central part will be limited to say 50 pixels and the peripheral will be so smooth that admits already a 10-fold reduction of resolution (as compared to level 1, so already 102 = 100 as compared to the original resolution!). In our example this is merely 250ı400 so the filter, even if it covers the whole image area, can be applied directly (without FFT at all). So in this case we proceed without FFT at all levels.
The first central component is termed the first level of pyramid. The “central part of periphery” is the second level of pyramid. And “periphery of periphery” constitutes the third level (which in principle can be subdivided that sample manner leading to the fourth level and further). How this pyramid looks in reality is shown in Figure 1.
While convolving with the 1st level, we work in the original resolution and keep the source and the processed images which totals to about 24 Gb. The memory to keep the filter of 50 pixels radius is negligible.
At the second level, all is the same only for the tenfold reduced resolution, so it needs merely about 24=102 Gb = 240 Mb that is negligible. The third level takes even less.
So the total memory and speed are determined by the first level and are admissible. 2.2
In the above reasoning we spoke about tenfold reduction of resolution which is surely
imaginary. In reality reduction ratio is determined by the desired accuracy (which is
in turn determined by the kernel and requirements to the “smoothness” of periphery),
and is usually just 3. That is, resolution of the image processed at the next level is 3
times lower than that at the current level, and each pixel of the image at the next level
“comprises” a group of 3ı3=9 current level pixels. In the simplest case, the value stored
in that pixel is just the sum over that group (i.e. 9 pixels of the current level), though
later we shall also use another relations). The first 3 levels of pyramid for one of the
filters used in simulation of human vision acuity [
For the lower reduction ratio we shall need more levels to reach finally so low resolution that can be processed directly (without FFT). It can be already 7.
The radius of the central part which in our idea was said to be 50 pixels is, surely, not that round. It is calculated from the criterion: outside that area the function changes so gradually that the inaccuracy of the linear interpolation from the 3fold reduced resolution is below allowed threshold. For human vision kernel, it is usually somewhere between 7 thru 25 pixels depending on other parameters.
All that and other details are explained below. How the second level (of the kernel and image) is constructed? Let the reduction ratio be 2M + 1 where M is a positive integer. Than each pixel of the next level substitutes somehow the group (2M + 1) (2M + 1) pixels of current level “around” it.
It is convenient to take M = 1 so that reduction of resolution is 3fold, and in our calculations we used that value. Though M = 2 and larger is also possible.
Obviously, the convolution of the 1st level (i.e. the original) image with the 1st level kernel can be written as f (x; y) = X K((x0
Pyramid of Filters — Fast Image Filtering without FFT 5 where x; y; x0; y0 are the pixel coordinates i.e. the integer indices, K is the kernel as a function of the physical (continuous!) coordinates (meters, degrees, ...), while is the pixel step i.e. the change of that physical coordinate from one pixel to another.
Grouping the pixels of the 1st level into “clusters” of (2M + 1) (2M + 1) pixels, so that on the original resolution grid their centers are
((2M + 1)X0 + 1; (2M + 1)Y 0 + 1) where (X0; Y 0) are the pixels at the reduced resolution, we can write result of convolution at the point x = (2M + 1)X + 1; y = (2M + 1)Y + 1 (present in both resolutions!) as f (x; y) = X
M X
M
X
X0;Y 0 = M = M 0 (2M + 1)
K((X0 Near the origin the kernel usually changes steeply while far from it, for r r0 (below we shall calculate that r0), the kernel becomes so smooth that for it the original high resolution becomes superfluous. So according to our idea the kernel (which we shall assume isotropic just for convenience) is decomposed into the sum of the “central part” (vanishing for r > r0) and “periphery”:
K(r) = Kc(r) + Kp(r)
The peripheral part Kp(r) is known for r r0 where by construction it equals K(r), while for r r0 it is arbitrary. It is natural to require that there Kp(r) be maximally simple and smooth, e.g. a parabola seamlessly going into K(r) at r = r0: where
Kp(r) = (
+ r2; r K(r); r r0 r0 Kc(r) = = K(r0) =
K0(r0)
2r0 (K(r) 0; r0K0(r0) 2 r2; r r r0 r0 (4) (5) (6) (7) (8) (9) (10)
+ X So far we did not introduced any approximation but just regrouped the original expression and named the results of convolution with the central respectively peripheral part of the kernel as fc(x; y) and fp(x; y).
Now let us assume that r0 is so large that inside one pixel of 2nd level (i.e. for all (2M + 1) (2M + 1) pixels of the 1st level it comprises) Kp can be approximated by the linear function (effectively the first terms of the Taylor series)
Kp((X0 where
D(r)
0 Then introducing r 1 @Kp = 2(2M + 1) 2 @r M X
M
X = M = M Fm;n(X; Y ) m nf (X + ; Y + ) we have in the 2nd level pixel (X; Y ), fp(X; Y )
X) 0; (Y 0
Y ) 0)F0;0(X0; Y 0) X Kp((X0 X0;Y 0 + +
X (X0 X0;Y 0 X (Y 0 X0;Y 0
X)D((X0
X) 0; (Y 0
Y ) 0)F1;0(X0; Y 0) Y )D((X0
X) 0; (Y 0
Y ) 0)F0;1(X0; Y 0) (11) which is nothing but the convolution on the coarse pixel grid of 2nd level of resolution.
Deviation of the kernel from its linear approximation can be estimated through the second derivatives. For r r0 the relative error is
2 while within the central zone it is The use of the derivative in r2, not in r, is convenient when the kernel is given as an analytic function of r2, as it is in human vision. If the kernel is defined as a tabulated function of r,
Pyramid of Filters — Fast Image Filtering without FFT 7
r=r0 Usually @(r2) +2r2K@(r2) @(r2) decreases for r ! 1, so the relative error in any pixel @K @ @K is bounded by
K r=r0 Then r0 can be found as the radius for which the above expression equals the error of approximation. It can be done with e.g. bisection because of monotone decrease.
Construction of the third level i.e. decomposition of now Kp into its central and peripheral parts is done similarly, just taking Kp instead of K. The central part of Kp becomes level 2 in the pyramid of filters, while peripheral part of Kp is again decomposed, its central part becoming level 3 of the pyramid of filters and so on.
Notice the filter of the n-th (for n > 1) level is not scalar, but its value is a “tuple” n
Kp; r 1 @Kp o. The source image at this level also keeps a “triad” fF0;0; F0;1; F1;0g @r of sums over sub-pixels (i.e. of (2M + 1) (2M + 1) pixels of the previous level) in each its color channel. 3
These words mean a function of several arguments is a product of univariate function. Namely, let us look at the famous Gaussian filter kernel
K = C2e (r=a)2 :
K = Ce (x=a)2
Ce (y=a)2 k(x)k(y) so the convolution with it is just two successive one-dimensional convolutions (12) (13) (14) (15) (16) (17) For a filter of radius 100 pixels a direct calculation of a 2D convolution requires about 1002 operations per pixel, while two 1D convolutions need only 200 operations, i.e. 50 times less. For larger radius the gain can be giant. So it would be great if this method can be somehow extended to work with more general kernels besides Gaussian.
It happens that really a 2D convolution can be calculated as a series of 1D ones for
an isotropic kernel K(r) and even just a “Hermitian” kernel K(x; y) = K(y; x), this is
a generalization of “separable filtering” from [
Indeed, a Hermitian matrix K(x; y) = K(y; x) of dimension N N can be expanded as a sum of “eigenprojectors”
N 1 K(x; y) = X m=0 m m(x) m(y) (20) where m — is a normalized eigenvector, m being the corresponding eigenvalue. Notice that a matrix with nonnegative elements can have some negative eigenvalues.
The most contribution in (20) comes from eivenvectors with larger j j so that we order in the absolute value of the eigenvalues and terminate the series (20) when the sum of remaining 2m becomes less than the desired accuracy. Notice for a smooth kernel, the first eigenvectors are also smooth, while the last oscillate rapidly, changing the sign from pixel to pixel. This is a usual property.
For kernels used in simulation of the human vision, a good accuracy about 1% is achieved already for just 2 – 4 (first) eigenvectors. The convolution now is calculated as:
f (x; y) =
X k(x0 It is not necessary to keep all the images gm for all pixels. For a filter of size N is enough to keep just N rows of each of them. 4
We can combine the two approaches. Either we first construct the pyramid of filters and then each its level is factored.
In principle the opposite order is also possible, i.e. we first decompose the original kernel with (20) and then construct the pyramid for each of its m. But for this order (18) (19) (21) (22) N it
Pyramid of Filters — Fast Image Filtering without FFT 9 we should need to calculate eigenvectors (at least the main) of the matrix of size about 40000ı25000 which is too expensive.
For all levels beyond the 1st the filter is a tuple nKp; r 1 @@Krp o so one must factor both its components and so that the number of eigenvectors retained be the same. We find it from the obvious condition that Kp itself is approximated well enough, because the gradient component, r 1 @@Krp is used for improvement of accuracy to the 2nd order and is thus less essential than the main term Kp. 5
The method was applied to HDRI image of size 4000ı2500 shown in Figure 2 (left). Filter was K = 2 1a2 (1+(r=1a)2)3=2 with a = 12:5 pixels and was defined in area 2500ı2500 pixels.
We compared our methods with the exact filtering based on FFT, shown in Figure 2 (right). The results of our methods are shown in Figure 3. All methods and their combinations give very similar, acceptable results. The difference between images obtained by the 5 filtering methods is below 0.67%.
The memory and speed for computer with i7–4800MQ @ 2.7Ghz and 16Gb DDR3 RAM are shown in Table 1 pyramid 1st order pyramid 2nd order FFT with factorization 12.1 sec, 276 Mb 22.3 sec, 276 Mb n/a no factorization 27.8 sec, 196 Mb 39.3 sec, 218 Mb 18.8 sec, 2598 Mb
pyramid (no factorization, no gradient)
pyramid with gradient, no factorization pyramid with factorization, no gradient
pyramid with factorization and gradient We elaborated and tested the method that allows to filter images of very high resolutions without FFT. The new method needs much less memory while the speed is comparable with that achieved with FFT, and in case of pyramid with factorization and no gradient even higher. The deviation from the standard exact method is meanwhile very low, below 1% which is admissible for visual applications. It can be further improved by adding support of gradients (2nd order) and manipulating the settings of the method. This however increases memory and decreases speed, so the results shown in Section 5 are just a reasonable compromise.