Choix d’une bibliothèque de calculs matriciels pour un SIG

Les SIG font un usage intensif de matrices afin d’afficher leurs cartes ou transformer des coordonnées. On pourrait croire que le marché est suffisamment bien pourvu en excellentes bibliothèques de calculs matriciels, libres ou commerciales. Pourtant, les SIG ont des besoins spécifiques qui divergent un peu des objectifs de plusieurs bibliothèques existantes. Cet article présente les qualités que nous recherchons, ainsi que la façon dont MapFaces y répond.

Parmi les sortes d’opérations qu’un SIG doit effectuer sur les coordonnées géographiques, il en est une à la fois simple et très fréquente. Ce sont les opérations linéaires, constituées uniquement d’une combinaison d’additions et de certaines multiplications. Ces opérations n’effectuent pas de projections cartographiques, plus complexes, mais couvrent de nombreux autres cas:

• Changer l’ordre des axes, par exemple de (latitude, longitude) vers (longitude, latitude).
• Changer la direction des axes (par exemple l’axe des y des images pointe souvent vers le bas).
• Changer de méridien d’origine (par exemple de Paris vers Greenwich).
• Changer le nombre de dimensions (par exemple passer des coordonnées 3D vers 2D).
• Convertir des unités de mesures (par exemple convertir des pieds en mètres).
• Convertir les coordonnées pixels d’une image en coordonnées géographiques (par exemple la conversion exprimée dans les fichiers « .tfw » qui accompagnent certaines images TIFF).
• Prendre en charge une petite partie des projections cartographiques (par exemple les paramètres « False Easting », « False Northing » et « Scale factor »).
• Appliquer des rotations, translations, échelles ou cisaillements (des transformations dites affines).

Les opérations linéaires ont la propriété de toujours se combiner: peu importe le nombre d’opérations linéaires que l’on enchaîne, le résultat sera toujours exprimable par une seule opération linéaire. Cette propriété est plus facilement visible lorsque les opérations linéaires sont exprimées sous forme de matrices: pour les combiner, il suffit de multiplier les matrices.

Exemple: supposons que nous disposons d’une image selon une projection cartographique quelconque. Supposons que les coordonnées des pixels sont représentées par (i,j) et que nous souhaitons les convertir en coordonnées projetées (x,y). Cette transformation peut être représentée par la matrice suivante:

Concentrons notre attention sur la matrice du milieu. Une conversion des coordonnées des pixels vers les coordonnées de la projection cartographique pourrait s’exprimer par les valeurs ci-dessous, si nous n’interchangeons ni n’inversons la direction d’aucun axe:

Si l’on veut en outre appliquer les opérations suivantes:

• inverser la direction de l’axe des y (pour se conformer à la convention la plus courante appliquée aux images);
• interchanger l’ordre des axes (ce qui se produit souvent pour les systèmes de coordonnées géographiques exprimant la latitude avant la longitude)

alors on peut exprimer cette conversion par la matrice ci-dessous, qu’il suffira de multiplier avec la matrice ci-dessus:

Note: la matrice de gauche est un cas particulier de la matrice de rotation autour de l’origine de l’espace pour un angle de θ=90°. Notre exemple équivaut à une rotation, mais la conversion appliquée dans la pratique ne sera pas toujours une rotation.

L’élément clé est qu’il n’y a pas besoin d’écrire un code dédié à l’inversion des axes. Cette opération, et bien d’autres, est prise en compte naturellement par l’algèbre matricielle. On y gagne en généricité du code et en performance. Dans notre exemple, le résultat de la multiplication est:

Des manipulations de matrices comme l’exemple ci-haut interviennent dans quasiment toutes les opérations appliquées sur des coordonnées. Mais l’analyse de ces opérations révèle quelques patterns:

• Ces matrices sont presque toujours de petites tailles, dépassant rarement 5 lignes par 5 colonnes.
• Les opérations matricielles « lourdes » (multiplications ou inversions de matrices) ne surviennent pas dans des endroits où la performance est importante. Dans la quasi-totalité des cas, elles ne sont effectuées qu’une fois pour toute, à la lecture d’un fichier, ou lors des étapes de préparation avant de convertir des coordonnées. Elles ne surviennent quasiment jamais dans la boucle convertissant chacune des coordonnées.
• Dans une succession de multiplications et d’inversions de matrices, les erreurs d’arrondissement s’accumulent et grandissent rapidement au point de se confondre avec certaines opérations légitimes, notamment les changements de référentiel.

Exemple: un changement de référentiel s’exprime souvent par un changement de la taille, position et orientation de l’ellipsoïde choisie comme approximation de la forme de la Terre. Les changements de la taille s’expriment en parties par million et les rotations en arc-secondes. Retranscrites dans une matrice, ces valeurs sont donc assez petites.

• Il arrive fréquemment que des matrices s’annulent en tout ou en partie, c’est-à-dire que leurs multiplications ramènent des facteurs d’échelles à 1 et des translations à 0. Toutefois les erreurs d’arrondissements font que les valeurs obtenues sont rarement exactes, mais plutôt des valeurs s’en rapprochant telles que 0,9999…97 à la place de 1. Malheureusement, les erreurs d’arrondissement sont parfois telles qu’il est difficile de savoir si certains coefficients de la matrices sont des artefacts ou proviennent d’un réel changement de référentiel.

Ces points font que, pour un SIG, la précision d’une bibliothèque de calculs matriciels est plus importante que la performance. Paradoxalement, un bon moyen de gagner en performance est justement d’investir davantage de temps de CPU pour effectuer des opérations matricielles plus précises, car on augmente ainsi les chances de détecter correctement quelles opérations s’annulent. L’effort investit dans cette détection permet de sauver du temps là où ça compte: quand viendra le moment de boucler sur des millions de coordonnées à transformer. Ces concepts expérimentés pour le projet MapFaces ont permis de favoriser le traitement de données à la volée dans le cadre d’élaboration de clients SIG.

Mais les bibliothèques dédiées aux calculs matriciels sont souvent conçues pour opérer de manière très performante sur des matrices de grandes tailles, ayant par exemple des milliers de lignes et colonnes. Elles sont ainsi conçues pour être capable de résoudre efficacement des systèmes d’équations linéaires comportant des centaines d’inconnues. Les problèmes qu’elles résolvent sont certes difficiles, mais assez différents de ceux qui intéressent nos cas d’usage. Pour cette raison, et aussi à cause d’un autre besoin spécifique détaillé dans la section suivante, nos travaux utilisent leurs propres fonctions de calculs matriciels. Ces fonctions tentent de résoudre le problème de précision en utilisant l’arithmétique « double-double » (une technique permettant de simuler une précision d’environ 120 bits) au prix de la performance dans une partie du code où elle n’est pas jugée critique.

Que faire des matrices qui ne sont pas carrées (et pourquoi)

Nous avons très souvent besoin d’inverser des matrices, afin d’obtenir une conversion de coordonnées qui fasse le contraire de la conversion originale. Mais on n’inverse habituellement que des matrices carrées. Or, notre code a besoin d’effectuer des inversions de matrices non-carrées. Selon que l’on ait plus de lignes ou plus de colonnes:

• Pour nos cas d’usage, une matrice non-carrée est une conversion qui ajoute ou supprime une dimension aux coordonnées.
• Pour les bibliothèques d’algèbre linéaire, une matrice non-carrée est un système d’équations sous-déterminé ou surdéterminé.

Pour mieux comprendre les difficultés que causerait une transposition trop directe des bibliothèques d’algèbre linéaire aux SIG, imaginons une conversion qui projetterait les points d’un espace 3D vers une surface 2D:

Pour des bibliothèques d’algèbre linéaire, la matrice représentant cette conversion serait un système d’équations sous-déterminé, et donc insoluble. C’est-à-dire qu’on ne peut pas inverser cette conversion pour obtenir (λ₂, φ₂) → (λ₁, φ₁, z) puisqu’on ne sait pas quelle valeur donner à z, ce qui implique qu’on ne peut pas trouver (λ₁, φ₁) non-plus car ces valeurs dépendent peut-être de z. Toutefois, dans le cas des SIG, l’axe des z est très souvent perpendiculaire à la surface sur laquelle sont exprimées les coordonnées (λ, φ). Cette perpendicularité rend λ₁ et φ₁ indépendants de z. Dans ce cas particulier, et ce cas seulement, on peut encore sauver les meubles.

Notre programme procède en vérifiant si les coordonnées z sont indépendantes des coordonnées λ et φ. Nous reconnaissons ce cas en vérifiant quels coefficients de la matrice ont la valeur zéro. Si SIS arrive à identifier des dimensions indépendantes, il peut les exclure temporairement de manière à inverser sans ambiguïté la conversion dans les dimensions restantes. S’il ne trouve pas de dimension indépendante, alors une exception est levée.

Si une inversion a été possible, alors il reste à décider du sort des dimensions que nous avions temporairement exclues. Dans notre exemple, on assignera la valeur NaN (Not-a-Number) aux valeurs de z dans la conversion (λ₂, φ₂) → (λ₁, φ₁, z). Là encore, le choix du coefficient à mettre à NaN dans la matrice est basé sur la présomption qu’elle représente une conversion de coordonnées.

Le traitement particulier fait par ce code permet donc d’inverser des matrices que l’on rencontre couramment dans les SIG, même si en principe le système est sous-déterminé. Dans notre exemple la coordonnée z reste inconnue – nous ne faisons pas surgir de l’information du néant – mais au moins les coordonnées (λ, φ) ont pu être récupérées.

Le problème inverse, celui des systèmes surdéterminés, est plus subtil. Une approche classique des bibliothèques d’algèbre linéaire est de résoudre les systèmes surdéterminés par la méthode des moindres carrées. Transposée à notre exemple, cette approche proposerait une conversion (λ₂, φ₂, z) → (λ₁, φ₁) qui semble le meilleur compromis pour diverses valeurs de λ₂, φ₂ et z, tout en n’étant (sauf cas particuliers) une solution exacte pour personne. De plus, les éventuelles combinaisons linéaires entre ces trois variables sont délicates compte tenu de l’hétérogénéité des unités de mesures, où les z sont en mètres et (λ, φ) en degrés. Apache SIS procède plutôt comme pour les systèmes sous-déterminés: en exigeant que certaines dimensions soient indépendantes des autres, faute de quoi la matrice sera considérée non-inversible. Dans le cas des systèmes surdéterminés SIS refusera donc d’effectuer certaines opérations que les bibliothèques d’algèbre linéaire auraient faite, mais garantira que les conversions obtenues sont exactes (aux erreurs d’arrondissement prêts).

En résumé, les besoins qui ont amené MapFaces à fournir ses propres fonctions de calculs matriciels sont:

• Structure légère pour les petites matrices, particulièrement celles de taille 3×3.
• Précision accrue avec l’arithmétique « double-double », quitte à sacrifier un peu de performance dans des endroits où elle n’est pas critique.
• Traitement particulier de l’inversion des matrices non-carrées pour des conversions de coordonnées.