◽️ 👩🏿‍💻 🔇 Calcul du centre de masse pour O (1) à l'aide d'images intégrées 🥐 👩🏻‍🎨 ☦️

Une image intégrée est un algorithme qui vous permet de calculer efficacement la somme des valeurs contenues dans un sous-ensemble rectangulaire d'un tableau multidimensionnel. Son idée même remonte à l'étude des fonctions de distribution de probabilité multidimensionnelle, et jusqu'à présent, il a trouvé une application réussie dans les domaines qui utilisent directement la théorie des probabilités comme boîte à outils principale. Par exemple, dans la reconnaissance des formes.

Aujourd'hui, nous allons examiner un cas curieux de la façon d'appliquer des images intégrées dans un domaine radicalement différent - la physique computationnelle. À savoir, voyons ce qui se passe si nous calculons avec eux le centre de masse du champ d'impulsions et quels avantages peuvent être tirés de cette symbiose.

Dans cet article, je dirai:

De quel type de tâche s'agit-il?
En savoir plus sur les images intégrées;
N (-);
;
, , .

N

Une solution rigoureuse à l'interaction gravitationnelle d'un système de N corps fait référence à des algorithmes qui ont une complexité quadratique: tous les autres corps du système ont un effet gravitationnel sur chaque corps ^[1] . Par conséquent, la force de la nécessité de l' interaction gravitationnelle de calculer pour chacun des N ² paires / 2.

Pour illustrer à quel point une solution stricte prend du temps, vous pouvez essayer de corréler la complexité du calcul d'un système réel avec les performances des ordinateurs modernes.

Pour calculer la force gravitationnelle entre deux corps dans un espace tridimensionnel, vous devez effectuer 20 opérations avec une virgule flottante (FLOP). Dans le système solaire, il y a environ 100 corps de plus de 200 kilomètres (y compris le Soleil, 9 planètes, planètes naines et satellites) ^[2]. Le nombre approximatif d'astéroïdes en orbite proche de la Terre est d'environ 20 000 ^[3] . Dans la ceinture d'astéroïdes, il y a de 1,1 à 1,9 million de corps de plus de 1 km ^[4] et environ un million d'astéroïdes similaires dans les groupes «troyen» ^[5] et «grec» de Jupiter. Dans la ceinture de Kuiper, environ 100 millions d'objets de plus de 20 km ^[6] et environ 2 trillions de plus se trouvent dans le nuage d'Oort ^[7] .

La puissance de calcul requise pour résoudre strictement le dernier problème gravitationnel n'est que d'un ordre de grandeur inférieur à la puissance de calcul requise pour simuler le travail du cerveau humain au niveau neuronal (2,5 × 10 ²⁶ ) ^[8]. C'est pourquoi sa solution rigoureuse est généralement remplacée par une solution approximative.

L'un des algorithmes les plus largement utilisés pour la solution approximative du problème gravitationnel - Tree Code - a une complexité temporelle quasi-linéaire O (N * logN) ^[9] . Le code d'arbre crée un arbre spatial et pour chaque nœud de cet arbre calcule la masse totale et le centre de masse de tous les corps qui y tombent. De plus, lors du calcul des forces gravitationnelles agissant sur chaque corps, le code d'arbre peut prendre en compte non pas l'influence directe d'autres corps, mais l'influence des nœuds d'arbre, et en fonction de la distance, il sélectionne des nœuds de plus en plus grands qui contiennent les paramètres d'un sous-ensemble de plus en plus nombreux de corps.

Illustration de Wikipédia. Sur le corps au centre se trouvent les masses de nœuds indiquées par des rectangles verts. Seules les masses des corps les plus proches sont directement prises en compte.

Problème de gravité pour le champ de quantité de mouvement

Le champ de quantité de mouvement est un champ vectoriel discret mv (p) , qui associe une paire masse-vitesse à chaque point de l'espace p considéré . Le champ d'impulsions peut être réglé pour un espace de n'importe quelle dimension. Une paire masse-vitesse caractérise l'impulsion et représente un vecteur de dimension R + 1, où R est la dimension de l'espace pour lequel le champ d'impulsions est donné. Par exemple, pour R = 2, cela pourrait être le vecteur { v _x , v _y , m }.

Il convient de noter que la version mathématique traditionnelle de la description de ce système est une combinaison d'un champ de vitesse vectorielle et d'un champ de masse scalaire. Dans ce cas, nous nous accordons la liberté de combiner vitesse et masse dans un même vecteur, car nous ne prétendons pas être mathématiques.

Une petite partie du champ vectoriel d'impulsions: les masses sont indiquées par la couleur, les directions par des flèches Le

champ vectoriel d'impulsions d'une taille de 729 × 729 éléments. À gauche, les masses. À droite, des impulsions. Les impulsions de cette illustration sont codées par couleur au format HSV (la teinte est la direction de l'impulsion et la valeur affiche linéairement l'ordre de sa valeur absolue de sorte que les impulsions les moins distinctes (valeur = 0,05) soient de l'ordre de 10 ^-7 et les impulsions les plus brillantes (valeur = 1 , 0) sont de l'ordre de 10 ³ et plus).

Semblable à ceci - grille - les méthodes pour décrire les systèmes physiques sont largement utilisées en astrophysique pour simuler l'évolution des nuages de gaz, la formation des étoiles et des galaxies. Il s'agit notamment de la méthode SAMR ^[10] ou du modèle de grille d'hydrodynamique ^[11] .

Comme dans le cas du problème gravitationnel des N corps, l'évolution d'un champ discret doit prendre en compte l'influence gravitationnelle de toutes les masses qui composent ce champ sur chacun de ses éléments discrets. Il faut tenir compte du fait que la complexité du problème dépend indirectement de la dimension de l'espace R : par exemple, pour un champ sur un plan constitué de 1000 × 1000 éléments, le nombre total N sera de 10 ⁶éléments, et un champ de même taille dans l'espace tridimensionnel comprendra déjà 10 ⁹ éléments.

Néanmoins, il existe un certain nombre de trucs qui permettent de résoudre approximativement le problème gravitationnel pour le champ de quantité de mouvement. La méthode Multigrid utilise une discrétisation hiérarchique, supportant plusieurs grilles de différentes tailles ^[12] . L'extension multipolaire traite les groupes de sources proches les uns des autres comme un seul objet ^[13] . Multigrid a une complexité de calcul quasi-linéaire et une complexité de mémoire logarithmique. L'une des options d'extension multipolaire - FMM - démontre une complexité de calcul linéaire en échange d'une précision de calcul réduite ^[14] .

Ci-dessous, nous considérons une autre méthode qui nous permet de résoudre le problème gravitationnel pour un champ discret d'impulsions en temps quasi-linéaire et a une complexité linéaire en mémoire.

Image intégrée

L'image intégrale permet de calculer à temps constant la somme des éléments de l'image d'origine dans une région rectangulaire arbitraire ^[15] . En infographie, les images intégrées ont été proposées à l'origine comme une alternative au mipmapping et au filtrage anisotrope ^[16] . De plus, ils sont utilisés avec succès pour le traitement d'images numériques et dans les techniques de reconnaissance de formes.

Une image intégrée est l'un des exemples les plus évidents d'un compromis entre la complexité de calcul et la complexité de la mémoire. Un algorithme «naïf» pour calculer la somme des éléments d'image a une complexité temporelle de O (M * N) et une complexité de mémoire de O (1). L'image intégrée vous permet d'effectuer les mêmes calculs en O (1) temps et a une complexité de mémoire O (M * N).

Utiliser l'image intégrée pour calculer la somme des éléments dans une région donnée

L'algorithme de calcul de l'image intégrale est très simple, caractérisé par une complexité de calcul linéaire et est facilement parallélisé sous GPGPU ^[17] . Il est implémenté comme un filtre gaussien à deux passes ^[18] : dans la première passe, les sommes partielles pour les lignes sont considérées, dans la seconde, ces sommes partielles sont accumulées en colonnes.

Calcul de l'image intégrée en deux passes. À gauche, l'image d'origine. Au centre, les montants partiels des lignes. À droite, l'image finale intégrée.

À gauche, une image contenant des masses d'éléments. À droite, son image intégrale. Les images de gauche et de droite utilisent différentes échelles logarithmiques.

Une solution approximative au problème gravitationnel en utilisant une image intégrée

Ayant à sa disposition une image intégrale des masses, il n'est pas difficile d'y adapter la technique caractéristique de l'Arbre Code et de l'Expansion Multipolaire. Donc, pour chaque élément du champ discret:

Nous prenons directement en compte l'influence des masses des huit éléments voisins seulement;
Nous prenons en compte l'influence de huit régions voisines, constituées de neuf éléments (3 × 3), en calculant la somme de leurs masses à l'aide d'une image intégrée;
À chaque étape suivante, nous augmentons la taille de la région de 3 fois (9 × 9, 27 × 27, 81 × 81, etc.) jusqu'à ce que celle-ci dépasse la taille de l'ensemble du champ vectoriel.

Un calcul approximatif des forces agissant sur un élément du champ vectoriel d'impulsions à l'aide d'une image intégrale des masses

La complexité de la solution approximative du problème gravitationnel pour un champ d'impulsions utilisant une image intégrale croît de façon quasi linéaire, comme O (N * 8 * log3 (sqrt (N))) pour R = 2 et comme O (N * 26 * log3 (cbrt (N))) pour R = 3.

Cependant, sous cette forme, cette solution présente le même inconvénient que la technique Multigrid, à savoir la "rigidité" discrète perceptible des composantes basse fréquence de la fonction de gravité. La façon la plus simple de démontrer ce problème est clairement.

Les forces dans cette illustration sont codées par couleur au format HSV (la teinte est la direction de la force, la valeur affiche linéairement l'ordre de sa valeur absolue de telle sorte que les forces les moins distinctes (valeur = 0,05) sont de l'ordre de 10 ^-7 et les forces les plus brillantes (valeur = 1, 0) sont de l'ordre de 10 ³ et plus).

Dans l'illustration ci-dessus, l'effet de la "rigidité" des composantes basse fréquence de la fonction de gravité est perceptible à l'œil nu. Mais si dans «Multigrid» cette «rigidité» se produit en raison d'une diminution de la fréquence d'échantillonnage, alors dans l'image intégrée cela est dû au manque d'informations sur la nature de la distribution spatiale des masses. Le fait est qu'un algorithme qui rapproche les forces à l'aide d'une image intégrée considère que le centre de masse coïncide avec le centre de la région.

Le centre de l'illustration montre l'extremum de masse avec une distribution de masse relativement uniforme dans le reste du champ,

l' extremum de masse se situant dans chacune des quatre régions rectangulaires. Évidemment, la position du centre de masse de chaque région doit coïncider approximativement avec l'élément indiqué par le blanc, dont la masse est de trois ordres de grandeur supérieure à la masse de la plupart des éléments de la région indiquée par le jaune. Cependant, dans les calculs des composantes à basse fréquence de la fonction de gravité pour ces régions, on pense que leurs centres de masse coïncident avec les centres des régions indiquées par des cercles.

Ce problème peut être résolu en introduisant des informations supplémentaires dans l'image intégrée, ce qui permettra de calculer non seulement la somme des masses dans une région donnée, mais aussi les coordonnées du centre de masse.

Tout d'abord, essayons d'imaginer une image, dont chaque élément contient ses propres coordonnées. L'image intégrale correspondante contiendra la somme des coordonnées. Par conséquent, la somme d'une région arbitraire de cette image intégrée, divisée par le nombre total d'éléments dans cette région, sera égale aux coordonnées du centre de cette région.

Image intégrale des coordonnées

Prenons par exemple la région dans le coin inférieur gauche avec les coordonnées {3; 1} et en haut à droite avec les coordonnées {7; 5}. La somme des coordonnées de cette région {168; 120} - {18; 45} - {28; 0} + {3; 0} est {125; 75}, le nombre total d'éléments est de 25, donc les coordonnées de son centre seront {5; 3}.

Il ne reste plus qu'à généraliser ce cas particulier où l'on considère que tous les éléments de l'image ont le même coefficient de pondération égal à l'unité. Pour prendre en compte différents poids, nous multiplions par eux les coordonnées des éléments. Et nous obtenons les coordonnées pondérées du centre de la région si nous divisons la somme des coordonnées pondérées par la somme des facteurs de pondération.

Image intégrée des coordonnées pondérées. Une police plus grande indique les poids

Considérez la région dans le coin inférieur gauche avec les coordonnées {2; 1} et en haut à droite avec les coordonnées {5; 4}. La somme des coefficients de pondération de cette région 4,8 - 1,0 - 0,6 + 0,2 est 3,4, et la somme de ses coordonnées pondérées {11,1; 13,2} - {0,5; 2,0} - {1,5; 0,0} + {0,1; 0,0} est égal à {9,2; 11.2}. Ainsi, les coordonnées pondérées du centre de cette région sont {2,7; 3.3}.

Pour vous assurer que ce circuit fonctionne vraiment, vous pouvez le faire visuellement - en convertissant l'image intégrée avec les coordonnées pondérées en un champ de distance en utilisant la transformation de distance ^[19] .

Conversion d'une image intégrée avec des coordonnées pondérées en champ de distance. À gauche, une image des masses. L'image suivante est la distance au centre de masse pour une région de 3 × 3 éléments. Vient ensuite la distance aux centres de masse pour les régions de tailles 9 × 9 et 27 × 27 éléments. Les valeurs de distance dans cette illustration sont normalisées en fonction de la taille de la région utilisée pour l'échantillonnage.

Convertit une image intégrée avec des coordonnées pondérées en un champ de distance dirigé. La direction et la distance sont codées par couleur au format HSV, où Hue est la direction, Value est la distance normalisée. Les valeurs de distance dans cette illustration sont normalisées en fonction de la taille de la région utilisée pour l'échantillonnage.

Résumant tout ce qui précède:

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

. ― . ― .

Il est temps de prêter attention à l'une des principales caractéristiques des images intégrées, qui limite encore le champ d'application de celles-ci, à savoir la consommation des ressources et la stabilité numérique. Ainsi, selon la plage de valeurs que contient l'image d'origine, un type de données plus long peut être requis pour son image intégrale. Par exemple, lors du calcul des écarts-types, en plus de l'image d'origine (plage de valeurs 0..255), une image contenant des carrés des valeurs correspondantes (plage de valeurs 0..65535) est utilisée. Dans ce cas, la précision n'est pas suffisante pour calculer des images intégrées à grande échelle de 32 bits ^[20] .

Une situation similaire est observée dans le cas d'images intégrées avec des coordonnées pondérées. En soi, la valeur de la somme des coordonnées qui doivent être stockées dans l'image intégrée augmente en fonction de la taille de l'image N: quadratique pour le cas unidimensionnel 0 + 1 + 2 + ... + N - 1 = (N - 1) * N / 2 et cubiquement pour bidimensionnel N * (N - 1) * N / 2. Par exemple, pour stocker la somme d'une coordonnée entière pour une image de taille 2048 × 2048 (égale à 4292870144), un entier non signé 32 bits (dont la valeur maximale est 4294967295) suffit à peine. Lors du calcul d'images intégrées plus grandes, un débordement se produit.

L'utilisation de nombres à virgule flottante 32 bits dans les images intégrées retarde le problème du débordement d'une distance astronomique (10 ³⁸ - 10 ¹⁰), mais en même temps, la précision des calculs avec les coordonnées pondérées est considérablement réduite en raison de l'erreur de troncature accumulée lors de la sommation ^[21] .

La valeur absolue de l'erreur dans le calcul du centre de masse pondéré d'une région d'une taille de 4 × 4 éléments. L'image intégrée contient des nombres à simple précision (32 bits).

Valeur absolue de l'erreur dans le calcul du centre de masse pondéré d'une région d'une taille de 32 × 32 éléments. L'image intégrée contient des nombres à simple précision (32 bits).

La plus grande valeur absolue est obtenue par le calcul du centre de masse pondéré pour les petites régions. Une augmentation de la taille de la région de deux ordres de grandeur entraîne une diminution de l'ampleur de l'erreur absolue des calculs de quatre ordres de grandeur. Dans ce cas, aucune dépendance de l'erreur de calcul sur la plage de valeurs des coefficients de pondération (qui sont utilisés pour calculer les coordonnées pondérées des éléments) n'a été trouvée.

Une augmentation de la plage de poids n'affecte pas la valeur absolue de l'erreur dans le calcul du centre de masse pondéré. Le graphique montre les erreurs pour la «pire» région (dans le coin supérieur droit de l'image intégrée). La taille de l'image intégrée est de 256 × 256 éléments.

La valeur absolue de l'erreur dans le calcul du centre de masse pondéré diminue avec l'augmentation de la taille de la région. Le graphique montre les erreurs pour la «pire» région (dans le coin supérieur droit de l'image intégrée).

Sur la base de l'analyse décrite ci-dessus, nous pouvons conclure que l'utilisation de nombres à simple précision pour calculer les centres de masse pondérés à l'aide d'images intégrées n'a de sens pratique que pour des images d'environ 512 × 512 éléments. Avec une nouvelle augmentation de la taille, l'ampleur de l'erreur devient comparable à la taille de l'élément d'image. Et bien que cette erreur diminue avec l'augmentation de la taille de la région, ce sont les régions de petite taille qui ont le plus grand impact sur le résultat final - l'ampleur de la force de gravité - par conséquent, seules les erreurs correspondantes doivent être prises en compte.

Comme pour les images plus grandes, vous pouvez utiliser des coordonnées pondérées double précision ou ajouter un niveau de discrétisation supplémentaire: divisez l'image d'origine en plusieurs parties et lisez les images intégrales séparément pour chacune de ces parties. Du point de vue de la complexité de calcul, si plus d'une image intégrale est utilisée pour calculer la somme des éléments d'une région, mais que «plus d'un» est une constante, la complexité de l'algorithme de calcul de somme ne change pas.

Conclusion

L'exemple ci-dessus de l'utilisation d'une image intégrée peut probablement être adapté pour développer un algorithme optimal O (1) pour l'échantillonnage par importance (échantillonnage d'importance). Les algorithmes existants - et extrêmement gourmands en ressources du point de vue du GPU - ont une complexité linéaire O (N) et trouvent une application dans les méthodes modernes d'éclairage global ^{[22, 23]} .

En général, à mon avis, les images intégrées sont l'un des algorithmes les plus sous-estimés. Ils peuvent être une excellente alternative au mipmapping et au filtrage anisotrope en même temps, et compte tenu de la façon dont ce dernier est implémenté sur le GPU ^{[24, 25]} , il est fort possible qu'ils soient déjà une mesure plus efficace.

Références

Afficher la liste

Calcul du centre de masse pour O (1) à l'aide d'images intégrées

N

Problème de gravité pour le champ de quantité de mouvement

Image intégrée

Une solution approximative au problème gravitationnel en utilisant une image intégrée

Conclusion

Références

More articles: