Cálculo del centro de masa para O (1) utilizando imágenes integradas



Una imagen integrada es un algoritmo que le permite calcular eficientemente la suma de los valores encerrados en un subconjunto rectangular de una matriz multidimensional. Su idea se remonta al estudio de las funciones de distribución de probabilidad multidimensional, y hasta ahora ha encontrado una aplicación exitosa en aquellas áreas que utilizan directamente la teoría de probabilidad como el conjunto de herramientas principal. Por ejemplo, en reconocimiento de patrones.

Hoy consideraremos un caso curioso de cómo aplicar imágenes integradas en un campo radicalmente diferente: la física computacional. Es decir, veamos qué sucede si calculamos con ellos el centro de masa del campo de pulso, y qué beneficio se puede obtener de esta simbiosis.

En este artículo te diré:

  • ¿Qué tipo de tarea es esta en cuestión?
  • Lea más sobre imágenes integradas;
  • N (-);
  • ;
  • , , .

N


Una solución estricta de la interacción gravitacional de un sistema de N cuerpos se refiere a algoritmos que tienen complejidad cuadrática: todos los demás cuerpos en el sistema tienen un efecto gravitacional en cada cuerpo [1] . En consecuencia, la fuerza de la necesidad de interacción gravitatoria a calcular para cada uno de los N 2 /2 pares.

Para ilustrar cuán lenta es una solución estricta, puede intentar correlacionar la complejidad de calcular un sistema real con el rendimiento de las computadoras modernas.

Para calcular la fuerza gravitacional entre dos cuerpos en un espacio tridimensional, debe realizar 20 operaciones con un punto flotante (FLOP). En el sistema solar, hay alrededor de 100 cuerpos de más de 200 kilómetros (incluido el Sol, 9 planetas, planetas enanos y satélites) [2]. El número aproximado de asteroides en la órbita cercana a la Tierra es de aproximadamente 20 mil [3] . En el cinturón de asteroides, hay de 1,1 a 1,9 millones de cuerpos de más de 1 km [4] y alrededor de un millón de asteroides similares en los grupos "troyano" [5] y "griego" de Júpiter. En el cinturón de Kuiper, alrededor de 100 millones de objetos de más de 20 km [6] y alrededor de 2 billones más están en la nube de Oort [7] .



El poder computacional requerido para resolver estrictamente el último problema gravitacional es solo un orden de magnitud menor que el poder computacional requerido para simular el trabajo del cerebro humano a nivel neural (2.5 × 10 26 ) [8]. Es por eso que su solución rigurosa generalmente se reemplaza por una aproximada.

Uno de los algoritmos más utilizados para la solución aproximada del problema gravitacional, el Código de árbol, tiene una complejidad de tiempo cuasilineal O (N * logN) [9] . El Código de árbol construye un árbol espacial y para cada nodo en este árbol calcula la masa total y el centro de masa de todos los cuerpos que caen en él. Además, al calcular las fuerzas gravitacionales que actúan sobre cada cuerpo, el Código de árbol puede tener en cuenta no la influencia directa de otros cuerpos, sino la influencia de los nodos de los árboles y, según la distancia, selecciona nodos cada vez más grandes que contienen parámetros de un subconjunto de cuerpos cada vez más numeroso.


Ilustración de Wikipedia. En el cuerpo en el centro están las masas de nodos indicadas por rectángulos verdes. Solo las masas de los cuerpos más cercanos se tienen en cuenta directamente.

Problema de gravedad para el campo de impulso


El campo de momento es un campo de vector discreto mv (p) , que asocia un par de masa-velocidad con cada punto del espacio p en consideración . El campo de pulso se puede configurar para un espacio de cualquier dimensión. Un par masa-velocidad caracteriza el momento y representa un vector de dimensión R + 1, donde R es la dimensión del espacio para el cual se da el campo de pulso. Por ejemplo, para R = 2, este podría ser el vector { v x , v y , m }.

Vale la pena señalar que la versión matemáticamente tradicional de la descripción de este sistema es una combinación de un campo de velocidad vectorial y un campo de masa escalar. En este caso, nos permitimos la libertad de combinar velocidad y masa en un vector, ya que no pretendemos ser matemáticos.


Una pequeña porción del campo vectorial de pulsos: las masas se indican por color, las direcciones por flechas. El


campo vectorial de pulsos con un tamaño de 729 × 729 elementos. A la izquierda están las masas. A la derecha están los impulsos. Los pulsos en esta ilustración están codificados por colores en formato HSV (Hue es la dirección del pulso, y Value muestra linealmente el orden de su valor absoluto para que los pulsos menos distinguibles (Valor = 0.05) sean del orden de 10 -7 , y los pulsos más brillantes (Valor = 1 , 0) son del orden de 10 3 y superiores).

Similar a esto, la rejilla, los métodos para describir sistemas físicos son ampliamente utilizados en astrofísica para simular la evolución de las nubes de gas, la formación de estrellas y galaxias. Estos incluyen el método SAMR [10] o el modelo de cuadrícula de hidrodinámica [11] .

Como en el caso del problema gravitacional de N cuerpos, la evolución de un campo discreto debe tener en cuenta la influencia gravitacional de todas las masas que componen este campo en cada uno de sus elementos discretos. Es necesario tener en cuenta que la complejidad del problema depende indirectamente de la dimensión del espacio R : por ejemplo, para un campo en un plano que consta de 1000 × 1000 elementos, el número total N será 10 6elementos, y un campo del mismo tamaño en un espacio tridimensional ya incluirá 10 9 elementos.

Sin embargo, hay una serie de trucos que permiten resolver aproximadamente el problema gravitacional para el campo de impulso. El método Multigrid utiliza discretización jerárquica, que admite varias cuadrículas de varios tamaños [12] . La expansión multipolar trata los grupos de fuentes ubicados cerca uno del otro como un solo objeto [13] . Multigrid tiene complejidad computacional cuasilineal y complejidad de memoria logarítmica. Una de las opciones de expansión multipolar, FMM, demuestra la complejidad computacional lineal a cambio de una precisión computacional reducida [14] .

A continuación consideramos otro método que nos permite resolver el problema gravitacional para un campo discreto de pulsos en tiempo cuasilineal y tiene una complejidad lineal en la memoria.

Imagen integrada


La imagen integral permite calcular el tiempo constante de la suma de los elementos de la imagen original en una región rectangular arbitraria [15] . En gráficos de computadora, las imágenes integradas se propusieron originalmente como una alternativa al mipmapping y al filtrado anisotrópico [16] . Además, se utilizan con éxito para el procesamiento de imágenes digitales y en técnicas de reconocimiento de patrones.

Una imagen integrada es uno de los ejemplos más obvios de una compensación entre la complejidad computacional y la complejidad de la memoria. Un algoritmo "ingenuo" para calcular la suma de elementos de imagen tiene una complejidad temporal de O (M * N) y una complejidad de memoria de O (1). La imagen integrada le permite realizar los mismos cálculos en tiempo O (1) y tiene una complejidad de memoria O (M * N).


Uso de una imagen integrada para calcular la suma de elementos en una región determinada

El algoritmo para calcular una imagen integral es muy simple, se caracteriza por una complejidad computacional lineal y se paraleliza fácilmente bajo GPGPU [17] . Se implementa como un filtro gaussiano de dos pasos [18] : en el primer paso, las sumas parciales para las filas se consideran, en el segundo, estas sumas parciales se acumulan en columnas.


Cálculo de la imagen integrada en dos pasadas. A la izquierda está la imagen original. En el centro están las cantidades parciales para las líneas. A la derecha está la imagen integrada final.


A la izquierda hay una imagen que contiene masas de elementos. A la derecha está su imagen integral. Las imágenes de la izquierda y la derecha usan diferentes escalas logarítmicas.

Una solución aproximada al problema gravitacional utilizando una imagen integrada.


Teniendo a su disposición una imagen integral de las masas, no es difícil adaptar la técnica característica del Código de Árbol y la Expansión Multipolar. Entonces, para cada elemento del campo discreto:

  1. Tenemos en cuenta directamente la influencia de las masas de solo los ocho elementos vecinos;
  2. Tenemos en cuenta la influencia de ocho regiones vecinas, que constan de nueve elementos (3 × 3), calculando la suma de sus masas utilizando una imagen integrada;
  3. En cada paso posterior, aumentamos el tamaño de la región en 3 veces (9 × 9, 27 × 27, 81 × 81, etc.) hasta que este exceda el tamaño del campo vectorial completo.


Un cálculo aproximado de las fuerzas que actúan sobre un elemento del campo vectorial de pulsos usando una imagen integral de masas

La complejidad de la solución aproximada del problema gravitacional para un campo de pulsos usando una imagen integral crece de forma cuasilineal, como O (N * 8 * log3 (sqrt (N))) para R = 2 y como O (N * 26 * log3 (cbrt (N))) para R = 3.



Sin embargo, en esta forma, esta solución tiene el mismo inconveniente que la técnica Multigrid, es decir, la "rigidez" discreta perceptible de los componentes de baja frecuencia de la función de gravedad. La forma más fácil de demostrar este problema es claramente.


Las fuerzas en esta ilustración están codificadas por colores en formato HSV (Hue es la dirección de la fuerza, Value muestra linealmente el orden de su valor absoluto para que las fuerzas menos distinguibles (Value = 0.05) sean del orden de 10 -7 , y las fuerzas más brillantes (Value = 1, 0) son del orden de 10 3 y superiores).

En la ilustración anterior, el efecto de la "rigidez" de los componentes de baja frecuencia de la función de gravedad se nota a simple vista. Pero si en "Multigrid" esta "rigidez" se produce debido a una disminución en la frecuencia de muestreo, entonces en la imagen integrada esto se debe a la falta de información sobre la naturaleza de la distribución espacial de las masas. El hecho es que un algoritmo que aproxima fuerzas usando una imagen integrada considera que el centro de masa coincide con el centro de la región.


El centro de la ilustración muestra el extremo de la masa con una distribución de masa relativamente uniforme en el resto del campo. El extremo de la

masa cae en cada una de las cuatro regiones rectangulares. Obviamente, la posición del centro de masa de cada región debe coincidir aproximadamente con el elemento indicado por el blanco, cuya masa es tres órdenes de magnitud mayor que la masa de la mayoría de los elementos de la región indicada por el amarillo. Sin embargo, en los cálculos de los componentes de baja frecuencia de la función de gravedad para estas regiones, se cree que sus centros de masa coinciden con los centros de las regiones indicadas por los círculos.

Este problema se puede resolver mediante la introducción de información adicional en la imagen integrada, que permitirá calcular no solo la suma de las masas en una región determinada, sino también las coordenadas del centro de masa.

Primero, tratemos de imaginar una imagen, cada elemento del cual contiene sus propias coordenadas. La imagen integral correspondiente contendrá la suma de las coordenadas. Por lo tanto, la suma de una región arbitraria de esta imagen integrada, dividida por el número total de elementos en esta región, será igual a las coordenadas del centro de esta región.


Imagen integral de coordenadas

Tomemos, por ejemplo, la región en la esquina inferior izquierda con coordenadas {3; 1} y en la esquina superior derecha con las coordenadas {7; 5}. La suma de las coordenadas de esta región {168; 120} - {18; 45} - {28; 0} + {3; 0} es {125; 75}, el número total de elementos es 25, por lo tanto, las coordenadas de su centro serán {5; 3}.

Todo lo que queda por hacer es generalizar este caso particular en el que consideramos que todos los elementos de la imagen tienen el mismo coeficiente de peso igual a la unidad. Para tener en cuenta diferentes pesos, multiplicaremos las coordenadas de los elementos por ellos. Y obtendremos las coordenadas ponderadas del centro de la región si dividimos la suma de las coordenadas ponderadas por la suma de los factores de ponderación.


Imagen integrada de coordenadas ponderadas. La fuente más grande indica pesos

Considere la región en la esquina inferior izquierda con coordenadas {2; 1} y en la esquina superior derecha con las coordenadas {5; 4}. La suma de los coeficientes de ponderación de esta región 4.8 - 1.0 - 0.6 + 0.2 es 3.4, y la suma de sus coordenadas ponderadas {11.1; 13.2} - {0.5; 2.0} - {1.5; 0,0} + {0,1; 0,0} es igual a {9,2; 11.2}. Por lo tanto, las coordenadas ponderadas del centro de esta región son {2.7; 3.3}.

Para asegurarse de que este circuito realmente funcione, puede hacerlo visualmente al convertir la imagen integrada con coordenadas ponderadas en un campo de distancia utilizando la transformación de distancia [19] .


Convertir una imagen integrada con coordenadas ponderadas en un campo de distancia. A la izquierda hay una imagen de las masas. La siguiente imagen es la distancia al centro de masa para una región de 3 × 3 elementos. Luego está la distancia a los centros de masa para regiones con tamaños de elementos de 9 × 9 y 27 × 27. Los valores de distancia en esta ilustración están normalizados al tamaño de la región utilizada para el muestreo.


Convierte una imagen integrada con coordenadas ponderadas en un campo de distancia dirigido. La dirección y la distancia están codificadas por colores en formato HSV, donde Hue es la dirección, el valor es la distancia normalizada. Los valores de distancia en esta ilustración están normalizados al tamaño de la región utilizada para el muestreo.

Resumiendo todo lo anterior:

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


. ― . ― .


Es hora de prestar atención a una de las características clave de las imágenes integradas, que todavía limita el alcance de su aplicación, a saber, el consumo de recursos y la estabilidad numérica. Por lo tanto, dependiendo del rango de valores que contiene la imagen original, se puede requerir un tipo de datos más largo para su imagen integral. Por ejemplo, al calcular las desviaciones estándar, además de la imagen original (rango de valores 0..255), se usa una imagen que contiene cuadrados de los valores correspondientes (rango de valores 0..65535). En este caso, la precisión no es suficiente para calcular imágenes integradas a gran escala de 32 bits [20] .

Una situación similar se observa en el caso de imágenes integradas con coordenadas ponderadas. Por sí mismo, el valor de la suma de coordenadas que debe almacenarse en la imagen integrada aumenta dependiendo del tamaño de la imagen N: cuadráticamente para el caso unidimensional 0 + 1 + 2 + ... + N - 1 = (N - 1) * N / 2 y cúbicamente para bidimensional N * (N - 1) * N / 2. Por ejemplo, para almacenar la suma de una coordenada entera para una imagen de tamaño 2048 × 2048 (igual a 4292870144), un número entero sin signo de 32 bits (cuyo valor máximo es 4294967295) es apenas suficiente. Al calcular imágenes integradas más grandes, se produce un desbordamiento.

El uso de números de coma flotante de 32 bits en imágenes integradas retrasa el problema del desbordamiento por una distancia astronómica (10 38 - 10 10), pero al mismo tiempo, la precisión de los cálculos con coordenadas ponderadas se reduce significativamente debido al error de truncamiento acumulado durante la suma [21] .


El valor absoluto del error al calcular el centro de masa ponderado de una región con un tamaño de 4 × 4 elementos. La imagen integrada contiene números de precisión únicos (32 bits).


El valor absoluto del error al calcular el centro de masa ponderado de una región con un tamaño de 32 × 32 elementos. La imagen integrada contiene números de precisión únicos (32 bits).

El mayor valor absoluto se logra mediante el cálculo del centro de masa ponderado para regiones pequeñas. Un aumento en el tamaño de la región en dos órdenes de magnitud conduce a una disminución en la magnitud del error absoluto de los cálculos en cuatro órdenes de magnitud. En este caso, no se encontró dependencia del error de cálculo en el rango de valores de los coeficientes de ponderación (que se utilizan para calcular las coordenadas ponderadas de los elementos).


Un aumento en el rango de pesos no afecta el valor absoluto del error al calcular el centro de masa ponderado. El gráfico muestra los errores para la región "peor" (en la esquina superior derecha de la imagen integrada). El tamaño de la imagen integrada es de 256 × 256 elementos.


El valor absoluto del error al calcular el centro de masa ponderado disminuye al aumentar el tamaño de la región. El gráfico muestra los errores para la región "peor" (en la esquina superior derecha de la imagen integrada).

Con base en el análisis descrito anteriormente, podemos concluir que el uso de números de precisión simple para calcular centros de masa ponderados utilizando imágenes integradas tiene sentido práctico solo para imágenes de aproximadamente 512 × 512 elementos de tamaño. Con un aumento adicional en el tamaño, la magnitud del error se vuelve comparable con el tamaño del elemento de imagen. Y aunque este error disminuye al aumentar el tamaño de la región, son las regiones de pequeño tamaño las que tienen el mayor impacto en el resultado final, la magnitud de la fuerza de gravedad, por lo tanto, solo se deben tener en cuenta los errores correspondientes.

En cuanto a las imágenes más grandes, puede usar coordenadas ponderadas de doble precisión o agregar un nivel más de discretización: divida la imagen original en varias partes y lea las imágenes integrales por separado para cada una de estas partes. Desde el punto de vista de la complejidad computacional, si se usa más de una imagen integral para calcular la suma de los elementos de una región, pero este "más de uno" es una constante, la complejidad del algoritmo de cálculo de la suma no cambia.

Conclusión


El ejemplo anterior del uso de una imagen integrada probablemente se puede adaptar para desarrollar un algoritmo óptimo O (1) para el muestreo por importancia (muestreo de importancia). Los algoritmos existentes, y extremadamente intensivos en recursos desde el punto de vista de la GPU, tienen una complejidad lineal O (N) y encuentran aplicación en los métodos modernos de iluminación global [22, 23] .

En general, en mi opinión, las imágenes integradas son uno de los algoritmos más subestimados. Pueden ser una excelente alternativa al mipmapping y al filtrado anisotrópico al mismo tiempo, y teniendo en cuenta cómo se implementa este último en la GPU [24, 25] , es muy posible que ya sean una medida más efectiva.

Referencias



All Articles