✋ 🔉 🔻 Calculation of the center of mass for O (1) using integrated images 👨🏿‍🎤 ✌🏾 🍫

An integrated image is an algorithm that allows you to efficiently calculate the sum of the values enclosed in a rectangular subset of a multidimensional array. His very idea dates back to the study of multidimensional probability distribution functions, and until now he has found successful application in those areas that directly use probability theory as the main toolkit. For example, in pattern recognition.

Today we will consider a curious case of how to apply integrated images in a radically different field - computational physics. Namely, let's see what happens if we calculate with them the center of mass of the pulse field, and what benefit can be derived from this symbiosis.

In this article I will tell:

What kind of task is this in question;
Read more about integrated images;
N (-);
;
, , .

N

A rigorous solution to the gravitational interaction of a system of N bodies refers to algorithms that have quadratic complexity: all other bodies in the system have a gravitational effect on each body ^[1] . Consequently, the force of the gravitational interaction need to calculate for each of the N ² /2 pair.

To illustrate how time-consuming a strict solution is, you can try to correlate the complexity of calculating a real system with the performance of modern computers.

To calculate the gravitational force between two bodies in three-dimensional space, you must perform 20 operations with a floating point (FLOP). In the solar system, there are about 100 bodies larger than 200 kilometers (including the Sun, 9 planets, dwarf planets and satellites) ^[2]. The approximate number of asteroids in near-Earth orbit is about 20 thousand ^[3] . In the asteroid belt, there are from 1.1 to 1.9 million bodies larger than 1 km ^[4] and about one million similar asteroids in the “Trojan” ^[5] and “Greek” groups of Jupiter. In the Kuiper belt, about 100 million objects larger than 20 km ^[6] and about 2 trillion more are in the Oort cloud ^[7] .

The computational power required to strictly solve the last gravitational problem is only an order of magnitude smaller than the computational power required to simulate the work of the human brain at the neural level (2.5 × 10 ²⁶ ) ^[8]. That is why its rigorous solution is usually replaced by an approximate one.

One of the most widely used algorithms for the approximate solution of the gravitational problem - Tree Code - has a quasilinear time complexity O (N * logN) ^[9] . The Tree Code builds a spatial tree and for each node in this tree calculates the total mass and center of mass of all the bodies that fall into it. Further, when calculating the gravitational forces acting on each body, the Tree Code can take into account not the direct influence of other bodies, but the influence of tree nodes, and depending on the distance, it selects ever larger nodes that contain parameters of an increasingly numerous subset of bodies.

Illustration from Wikipedia. On the body in the center are the masses of nodes indicated by green rectangles. Only the masses of the nearest bodies are taken into account directly.

Gravity problem for the momentum field

The momentum field is a discrete vector field mv (p) , which associates a mass-velocity pair with each point of the space p under consideration . The pulse field can be set for a space of any dimension. A mass-velocity pair characterizes the momentum and represents a vector of dimension R + 1, where R is the dimension of the space for which the pulse field is given. For example, for R = 2 this could be the vector { v _x , v _y , m }.

It is worth noting that the mathematically traditional version of the description of this system is a combination of a vector velocity field and a scalar mass field. In this case, we allow ourselves the freedom to combine speed and mass in one vector, since we do not pretend to be mathematical.

A small portion of the vector field of pulses: the masses are indicated by color, the directions by arrows. The

vector field of pulses of 729 × 729 elements in size. On the left are the masses. On the right are impulses. The pulses in this illustration are color coded in HSV format (Hue is the direction of the pulse, and Value linearly displays the order of its absolute value so that the least distinguishable (Value = 0.05) pulses are of the order of 10 ^-7 , and the brightest pulses (Value = 1 , 0) are of the order of 10 ³ and higher).

Similar to this - grid - methods for describing physical systems are widely used in astrophysics to simulate the evolution of gas clouds, the formation of stars and galaxies. These include the SAMR method ^[10] or the grid model of hydrodynamics ^[11] .

As in the case of the gravitational problem of N bodies, the evolution of a discrete field must take into account the gravitational influence of all the masses that make up this field on each of its discrete elements. It is necessary to take into account that the complexity of the problem indirectly depends on the dimension of the space R : for example, for a field on a plane consisting of 1000 × 1000 elements, the total number N will be 10 ⁶elements, and a field of the same size in three-dimensional space will already include 10 ⁹ elements.

Nevertheless, there are a number of tricks that make it possible to approximately solve the gravitational problem for the momentum field. The Multigrid method uses hierarchical discretization, supporting several grids of various sizes ^[12] . Multipole Expansion treats groups of sources located close to each other as one object ^[13] . Multigrid has quasilinear computational complexity and logarithmic memory complexity. One of the Multipole Expansion options - FMM - demonstrates linear computational complexity in exchange for reduced computational accuracy ^[14] .

Below we consider another method that allows us to solve the gravitational problem for a discrete field of pulses in quasilinear time and has linear complexity in memory.

Integrated image

The integral image allows calculating the constant time the sum of the elements of the original image in an arbitrary rectangular region ^[15] . In computer graphics, integrated images were originally proposed as an alternative to mipmapping and anisotropic filtering ^[16] . In addition, they are successfully used for digital image processing and in pattern recognition techniques.

An integrated image is one of the most obvious examples of a trade-off between computational complexity and memory complexity. A “naive” algorithm for calculating the sum of image elements has a time complexity of O (M * N) and a memory complexity of O (1). The integrated image allows you to perform the same calculations in O (1) time and has O (M * N) memory complexity.

Using an integrated image to calculate the sum of elements in a given region

The algorithm for calculating an integral image is very simple, characterized by linear computational complexity and is easily parallelized under GPGPU ^[17] . It is implemented like a two-pass Gaussian filter ^[18] : in the first pass, the partial sums for the rows are considered, in the second, these partial sums are accumulated in columns.

Calculation of the integrated image in two passes. On the left is the original image. In the center are the partial amounts for the lines. On the right is the final integrated image.

On the left is an image containing masses of elements. On the right is its integral image. The images on the left and right use different logarithmic scales.

An approximate solution to the gravitational problem using an integrated image

Having at its disposal an integral image of the masses, it is not difficult to adapt the technique characteristic of the Tree Code and Multipole Expansion to it. So, for each element of the discrete field:

We directly take into account the influence of the masses of only the neighboring eight elements;
We take into account the influence of eight neighboring regions, consisting of nine elements (3 × 3), by calculating the sum of their masses using an integrated image;
At each subsequent step, we increase the size of the region by 3 times (9 × 9, 27 × 27, 81 × 81, etc.) until this one exceeds the size of the entire vector field.

An approximate calculation of the forces acting on an element of the vector field of pulses using an integral image of masses

The complexity of the approximate solution of the gravitational problem for a field of pulses using an integral image grows quasilinearly, like O (N * 8 * log3 (sqrt (N))) for R = 2 and like O (N * 26 * log3 (cbrt (N))) for R = 3.

However, in this form, this solution has the same drawback as the Multigrid technique, namely, the perceptible discrete "rigidity" of the low-frequency components of the gravity function. The easiest way to demonstrate this problem is clearly.

The forces in this illustration are color coded in HSV format (Hue is the direction of the force, Value linearly displays the order of its absolute value in such a way that the least distinguishable (Value = 0.05) forces are of the order of 10 ^-7 , and the brightest forces (Value = 1, 0) are of the order of 10 ³ and higher).

In the illustration above, the effect of "rigidity" of the low-frequency components of the gravity function is noticeable to the naked eye. But if in “Multigrid” this “rigidity” occurs due to a decrease in the sampling frequency, then in the integrated image this is due to the lack of information about the nature of the spatial distribution of masses. The fact is that an algorithm that approximates forces using an integrated image considers that the center of mass coincides with the center of the region.

The center of the illustration shows the mass extremum with a relatively uniform mass distribution in the rest of the field. The

mass extremum falls into each of the four rectangular regions. Obviously, the position of the center of mass of each region should approximately coincide with the element indicated by white, the mass of which is three orders of magnitude greater than the mass of most elements of the region indicated by yellow. However, in the calculations of the low-frequency components of the gravity function for these regions, it is believed that their centers of mass coincide with the centers of the regions indicated by circles.

This problem can be solved by introducing additional information into the integrated image, which will make it possible to calculate not only the sum of the masses in a given region, but also the coordinates of the center of mass.

First, let's try to imagine an image, each element of which contains its own coordinates. The corresponding integral image will contain the sum of the coordinates. Therefore, the sum of an arbitrary region of this integrated image, divided by the total number of elements in this region, will be equal to the coordinates of the center of this region.

Integral image of coordinates

Let’s take for example the region in the lower left corner with coordinates {3; 1} and in the upper right with coordinates {7; 5}. The sum of the coordinates of this region {168; 120} - {18; 45} - {28; 0} + {3; 0} is {125; 75}, the total number of elements is 25, therefore, the coordinates of its center will be {5; 3}.

All that remains to be done is to generalize this particular case in which we consider that all elements of the image have the same weight coefficient equal to unity. To take into account different weights, we will multiply the coordinates of the elements by them. And we will get the weighted coordinates of the center of the region if we divide the sum of the weighted coordinates by the sum of the weighting factors.

Integrated image of weighted coordinates. Larger font indicates weights

Consider the region in the lower left corner with coordinates {2; 1} and in the upper right with coordinates {5; 4}. The sum of the weighting coefficients of this region 4.8 - 1.0 - 0.6 + 0.2 is 3.4, and the sum of its weighted coordinates {11.1; 13.2} - {0.5; 2.0} - {1.5; 0,0} + {0,1; 0,0} is equal to {9,2; 11.2}. Thus, the weighted coordinates of the center of this region are {2.7; 3.3}.

To make sure that this circuit really works, you can visually - by converting the integrated image with weighted coordinates to a distance field using distance transform ^[19] .

Converting an integrated image with weighted coordinates to a distance field. On the left is an image of the masses. The following image is the distance to the center of mass for a region of 3 × 3 elements. Next is the distance to the centers of mass for regions with sizes of 9 × 9 and 27 × 27 elements. The distance values in this illustration are normalized to the size of the region used for sampling.

Converts an integrated image with weighted coordinates to a directed distance field. The direction and distance are color coded in HSV format, where Hue is the direction, Value is the normalized distance. The distance values in this illustration are normalized to the size of the region used for sampling.

Summarizing all of the above:

, , ( , R = 3) ― , , .
.
O(1).
O(M*N).

. ― . ― .

Time to pay attention to one of the key features of integrated images, which still limits the scope of their application, namely, resource consumption and numerical stability. So, depending on the range of values that the original image contains, a longer data type may be required for its integral image. For example, when calculating the standard deviations, in addition to the original image (value range 0..255), an image containing squares of the corresponding values (value range 0..65535) is used. In this case, accuracy is not enough for calculating large-scale integrated images of 32 bits ^[20] .

A similar situation is observed in the case of integrated images with weighted coordinates. By itself, the value of the sum of coordinates that must be stored in the integrated image grows depending on the image size N: quadratically for the one-dimensional case 0 + 1 + 2 + ... + N - 1 = (N - 1) * N / 2 and cubically for two-dimensional N * (N - 1) * N / 2. For example, to store the sum of one integer coordinate for an image of 2048 × 2048 size (equal to 4292870144), a 32-bit unsigned integer (the maximum value of which is 4294967295) is barely enough. When calculating larger integrated images, overflow occurs.

The use of 32-bit floating-point numbers in integrated images delays the problem of overflow by an astronomical distance (10 ³⁸ - 10 ¹⁰), but at the same time, the accuracy of calculations with weighted coordinates is significantly reduced due to the truncation error accumulated during summation ^[21] .

The absolute value of the error in calculating the weighted center of mass of a region with a size of 4 × 4 elements. The integrated image contains single precision numbers (32 bits).

The absolute value of the error in calculating the weighted center of mass of a region with a size of 32 × 32 elements. The integrated image contains single precision numbers (32 bits).

The largest absolute value is achieved by the calculation of the weighted center of mass for small regions. An increase in the size of the region by two orders of magnitude leads to a decrease in the magnitude of the absolute error of calculations by four orders of magnitude. In this case, no dependence of the calculation error on the range of values of the weighting coefficients (which are used to calculate the weighted coordinates of the elements) was found.

An increase in the range of weights does not affect the absolute value of the error in calculating the weighted center of mass. The graph shows the errors for the “worst” region (in the upper right corner of the integrated image). The size of the integrated image is 256 × 256 elements.

The absolute value of the error in calculating the weighted center of mass decreases with increasing size of the region. The graph shows the errors for the “worst” region (in the upper right corner of the integrated image).

Based on the analysis described above, we can conclude that the use of single-precision numbers for calculating weighted centers of mass using integrated images makes practical sense only for images of about 512 × 512 elements in size. With a further increase in size, the magnitude of the error becomes comparable with the size of the image element. And although this error decreases with increasing size of the region, it is regions of a small size that have the greatest impact on the final result - the magnitude of the force of gravity - therefore, only the corresponding errors should be taken into account.

As for larger images, you can use double-precision weighted coordinates or add one more level of discretization: divide the original image into several parts and read the integral images separately for each of these parts. From the point of view of computational complexity, if more than one integral image is used to calculate the sum of the elements of a region, but this “more than one” is a constant, the complexity of the sum calculation algorithm does not change.

Conclusion

The above example of the use of an integrated image can probably be adapted to develop an optimal algorithm O (1) for sampling by importance (importance sampling). Existing - and extremely resource-intensive from the point of view of GPU - algorithms have linear complexity O (N) and find application in modern methods of global lighting ^{[22, 23]} .

In general, in my opinion, integrated images are one of the most underrated algorithms. They can be an excellent alternative to mipmapping and anisotropic filtering at the same time, and taking into account how the latter is implemented on the GPU ^{[24, 25]} , it is quite possible that they are already a more effective measure.

References

View list

Calculation of the center of mass for O (1) using integrated images

N

Gravity problem for the momentum field

Integrated image

An approximate solution to the gravitational problem using an integrated image

Conclusion

References

More articles: