pyViewFactor v1.1.0 — noyau SA-30 & accélération BVH  

Un an après la dernière mise à jour significative de pyViewFactor , la version 1.1.0 marque une refonte significative des algorithmes de calcul. Au programme : un noyau d’intégration semi-analytique qui élimine les biais numériques sur les maillages fermés, et une structure hiérarchique d’accélération des tests d’obstruction — le tout compilé à la volée avec Numba pour un gain mesuré de ~33× sur un cas urbain de référence.

Rappel : le problème

Le facteur de forme $F_{i \to j}$ quantifie la fraction du rayonnement émis par la surface $i$ qui atteint la surface $j$ :

$$ F_{i \to j} = \frac{1}{A_i} \iint_{A_i} \iint_{A_j} \frac{\cos\theta_i \cos\theta_j}{\pi r^2} , dA_i , dA_j $$

Pour des polygones plans arbitraires, Mazumder & Ravishankar (2012) montrent qu’une double application du théorème de Green réduit cette intégrale de surface en une intégrale de contour :

$$ A_i F_{i \to j} = \frac{1}{2\pi} \oint_{\partial A_i} \oint_{\partial A_j} \ln(r_{kl}) ; d\mathbf{k} \cdot d\mathbf{l} $$

En discrétisant les contours en arêtes dirigées, on obtient une somme sur toutes les paires d’arêtes $(k, l)$ :

$$ A_i F_{i \to j} = \frac{1}{4\pi} \sum_{k} \sum_{l} (\mathbf{e}_k \cdot \mathbf{e}l) \int_0^1 \int_0^1 \ln r^2{kl}(t, s) ; dt , ds $$

où $r^2_{kl}(t,s)$ est un polynôme de degré 2 en $t$ et $s$. C’est sur cette forme que travaillent tous les noyaux de la bibliothèque.

Le problème des noyaux v1.0

En v1.0.3, le calcul de la matrice commençait par un tri des paires : les faces partageant un sommet ou une arête (touching) étaient séparées des faces disjointes. Les deux groupes étaient ensuite traités différemment :

• Paires touchantes → dblquad + décalage epsilon : la face $j$ était légèrement déplacée dans la direction centroïde-à-centroïde avant l’intégration, pour éloigner la singularité de $\ln(r^2)$. L’epsilon était un paramètre à ajuster : trop petit, la singularité restait problématique ; trop grand, la géométrie était déformée et le résultat biaisé. L’intégrale elle-même était ensuite évaluée par scipy.integrate.dblquad (quadrature adaptative sur les deux niveaux).
• Paires disjointes → GL-30×30 : règle de Gauss-Legendre à 30 points sur chacune des deux intégrales, sans décalage.

Le résultat : sur un maillage fermé, la somme des facteurs de forme d’une face est systématiquement sous-estimée d’environ 0,16 % — biais structurel lié au décalage epsilon, impossible à éliminer sans changer d’approche.

Le noyau v1.1.0 résout ce problème à la racine.

Précision et vitesse : le diagramme de Pareto

Précision vs temps de calcul par paire de faces (cas disjoint). En bas à gauche = optimal. Le noyau SA-30 (triangles orange) domine nettement GL-30 (cercles bleus, noyau v1.0) : erreur relative ~10⁻¹¹ contre ~10⁻⁴, à temps comparable.

Le noyau SA-30 : intégrale interne en forme fermée

Le noyau SA-30 (v1.1.0) évalue l’intégrale intérieure exactement, pour chaque point de quadrature extérieur. L’intégrale à résoudre est de la forme :

$$ I(x) = \int_0^1 \ln(A y^2 + B y + C) , dy $$

avec $A$, $B$, $C$ fonctions du point de quadrature extérieur $x$ et de la paire d’arêtes. Le discriminant $\Delta = 4AC - B^2$ détermine la branche :

Le cas $\Delta \leq 0$ est précisément celui des faces adjacentes ou coplanaires — traité exactement, sans décalage epsilon, sans biais. Les deux directions $F(i \to j)$ et $F(j \to i)$ sont calculées indépendamment (pas de remplissage par réciprocité), ce qui apporte un bénéfice supplémentaire de précision sur les maillages fermés.

import pyviewfactor as pvf

# Un seul appel — noyau SA-30 automatique
F = pvf.compute_viewfactor(face_i, face_j)

# Matrice complète N×N avec obstruction BVH
F_matrix = pvf.compute_viewfactor_matrix(mesh, obstacles=obstacle_mesh)

Accélération BVH pour les tests d’obstruction

Pour chaque paire de faces $(i, j)$, il faut vérifier si le segment centroïde-à-centroïde est bloqué par un triangle de l’obstacle. En v1.0, un premier filtre AABB plat (boîte englobante commune des deux faces) éliminait déjà les triangles manifestement hors de portée — mais les candidats restants étaient ensuite testés séquentiellement, un par un.

L’idée : une hiérarchie de boîtes englobantes (AABB)

Une Bounding Volume Hierarchy (BVH) organise les triangles de l’obstacle en un arbre binaire de boîtes englobantes alignées sur les axes (AABB). Chaque nœud de l’arbre est une AABB qui contient tous les triangles de son sous-arbre.

Construction (une seule fois, à l’initialisation) :

1. Calculer l’AABB de l’ensemble courant de triangles.
2. Découper selon le plus long axe, au centroïde médian.
3. Récurser jusqu’à ce que chaque feuille contienne un seul triangle.

Traversée pour un rayon :

1. Tester le rayon contre l’AABB du nœud racine (slab test : 6 comparaisons).
2. Raté → élaguer tout le sous-arbre, zéro triangle testé.
3. Touché → descendre dans les deux fils.
4. À la feuille → test de Möller–Trumbore sur le triangle.

Chaque niveau de l’arbre élimine environ la moitié des candidats restants : coût O(log N_tri) par rayon au lieu de O(N_tri).

L’arbre est stocké sous forme de cinq tableaux numpy plats (node_mins, node_maxs, node_left, node_right, node_tri) — sans pointeurs, sans objets Python — pour une traversée compatible Numba.

Résultats sur le cas de référence

Maillage urbain de 382 faces, 72 771 paires candidates, auto-obstruction complète :

Les deux effets se cumulent :

• BVH : réduit le coût d’obstruction de O(N_tri) à O(log N_tri) par rayon
• Numba prange : parallélise obstruction et intégration sur l’ensemble des paires

Nouvelle documentation : migration vers MkDocs

La v1.1.0 inaugure également une documentation entièrement réécrite. La génération automatique via pdoc a été remplacée par MkDocs + mkdocstrings, avec le thème Material.

La documentation est structurée en quatre sections :

• Théorie : dérivation de l’intégrale double, noyau SA-30 (forme fermée, cas $\Delta > 0$ et $\Delta \leq 0$), moteur BVH,
• Implémentation : architecture du pipeline, sémantique visibilité / obstruction, guide de performance,
• Exemples : de la paire de faces à l’analyse d’une scène urbaine complète,
• Référence API : générée automatiquement à partir des docstrings au format NumPy.

Elle est disponible à l’adresse : arep-dev.gitlab.io/pyViewFactor

Installation

pip install pyviewfactor

Numba est optionnel mais fortement recommandé. Sans lui, tous les noyaux fonctionnent en Python pur séquentiel avec des résultats identiques.