Compression des images hyperspectrales



These en pdf

Chapitre 3
Compression des images hyperspectrales

REVENONS maintenant sur la compression avec ce chapitre qui présente l’adaptation des méthodes de transformée en ondelettes pour les propriétés spécifiques de l’hyperspectral. Une première partie développe les possibilités du standard de compression JPEG 2000. La seconde partie s’attache à trouver la transformation en ondelettes optimale pour les images hyperspectrales. Enfin, une adaptation des structures d’arbres de zéros est faite pour appliquer les algorithmes EZW et SPIHT aux images hyperspectrales. Cette partie de l’étude a fait l’objet de publications [Chr06a,Chr06b,Chr06c,Chr07a,Chr07b].

3.1 JPEG 2000 pour l’hyperspectral

3.1.1 La norme JPEG 2000

JPEG 2000 est la norme la plus récente pour la compression d’image fixe. Une documentation complète est disponible sur le sujet, on peut voir par exemple [Tau02,Sko01]. Cette partie propose uniquement une présentation des points clés de la norme.

Cette norme a été définie pour fournir un cadre à une grande variété d’applications compressant les images avec différentes caractéristiques (images naturelles, images scientifiques multicomposantes, texte codé sur deux niveaux…) pour différentes utilisations (transmission en temps réel, archive, ressources limitées…).

Les principales fonctionnalités de JPEG 2000 sont, entre autres, le tuilage d’images de grande taille, le codage de régions d’intérêt (region of interest ou ROI) tout en proposant le codage d’images de dynamiques variées avec différents nombres de composantes. La partie 1 du standard [ISO02] comprend la plupart des fonctionnalités de JPEG 2000 pouvant être utilisées sans payer de royalties ou de frais de licence. La partie 2 contient les extensions [ISO04]. La plupart des implémentations actuelles de JPEG 2000 correspondent uniquement la partie 1.

Le codeur JPEG 2000 est basé sur le modèle de codage par transformée. Pour coder une image à une composante, une transformation multirésolution en ondelettes est d’abord appliquée. L’annexe A concerne les principaux aspects de l’utilisation de la transformée en ondelettes. Dans le cas d’une compression sans pertes, c’est l’ondelette Le Gall 5/3 qui est utilisée, sinon c’est l’ondelette de Daubechies 9/7. Pour la compression sans pertes, c’est une transformation d’entiers en entiers basée sur une implémentation en lifting scheme de la transformée en ondelettes [Dau98]. Ces transformations ont été choisies pour leur capacité à effectuer une décorrélation efficace d’une image.

La partie codage entropique de JPEG 2000 est assez complexe pour respecter la flexibilité demandée par la norme. Pour permettre le codage progressif ainsi que la résistance aux erreurs de transmission, le codage est basé sur un codage par blocs avec points de troncatures optimaux (Embedded Block Coding with Optimal Truncation ou EBCOT) [Tau00]. EBCOT est la raison principale des excellentes performances en terme de débit-distorsion de cette norme. Chaque sous-bande issue de la transformée en ondelettes est partagée en blocs de taille relativement faible (64 × 64 pixels en général). Chacun de ces blocs est ensuite codé indépendamment pour produire un train binaire progressif. Une quantification scalaire est utilisée pour la partie 1 de la norme tandis que la quantification vectorielle est disponible dans les extensions. Le codage entropique est réalisé par un codeur arithmétique contextuel. Les symboles sont codés à partir de modèles de probabilités pour 18 contextes différents. La norme ne spécifie pas d’allocation de débit et donc sur ce point, l’utilisateur est libre de faire ce qu’il veut. Néanmoins, une méthode populaire et largement utilisée est l’algorithme PCRD-opt (Post Compression Rate Distortion optimization) qui détermine le point optimal de troncature pour chacun des blocs [Tau00].

Pour les images multicomposantes à trois bandes (une image couleur classique par exemple), une transformation standard entre les couleurs est appliquée avant le codage par JPEG 2000. Cette transformation convertit l’espace naturel RGB (Red Green Blue) dans un espace moins corrélé YCrCb (une luminance et deux chrominances). Ensuite, le codeur JPEG 2000 classique est appliqué à chaque composante. L’algorithme d’allocation de débit (PCRD-opt) est ensuite utilisé de manière globale pour répartir le débit entre toutes les composantes.

3.1.2 Adaptation aux images hyperspectrales

Les parties 1 et 2 de la norme JPEG 2000 sont destinées au codage des images fixes en niveau de gris ou des images avec 3 bandes de couleurs et éventuellement un quatrième canal alpha (pour spécifier la transparence). Dans ces parties, aucune transformation intercomposante n’est définie à l’exception de la transformation couleur. Cependant, la partie 2 prévoit l’utilisation de transformations intercomposantes pouvant être précisées par l’utilisateur, avec, entre autres, la transformation en ondelettes. La partie 10, également connue sous le terme de JP3D, vise plutôt les images tridimensionnelles aussi isotropiques que possible [wg104]. Ces spécifications ne conviennent pas aux images hyperspectrales qui ne sont pas isotropiques et ont, par exemple, une corrélation bien plus forte dans la direction spectrale que dans les directions spatiales. Par conséquent, la partie JP3D du standard ne convient pas à la compression des images hyperspectrales. Nous proposons donc plutôt d’utiliser les extensions de la norme (partie 2) en introduisant l’utilisation de transformées intercomposantes avant d’appliquer le codage JPEG 2000.

L’implémentation de référence de JPEG 2000 est le Verification Model (VM). Le VM est utilisé par le comité JPEG 2000 pour les expérimentations et évoluera au final en implémentation de la norme JPEG 2000 [wg101]. Pour cette étude, la version VM 9.1, la plus récente, est utilisée pour évaluer les performances de compression de la norme JPEG 2000 adaptée aux images hyperspectrales. Une correction est nécessaire pour pouvoir décomprimer les trains binaires obtenus en utilisant la transformée intercomposante, cette correction est détaillée dans l’annexe C. Une autre implémentation récente et populaire de JPEG 2000, Kakadu, version 5.0 est également utilisée pour confirmer les résultats. Le détail des lignes de commande et des options utilisées est donné dans l’annexe C.

Comme il a été précisé, les parties 1 et 2 de la norme ne donnent pas de spécifications pour les cas où l’image possède plus de 3 bandes (transformation couleur). Cependant, le VM permet à l’utilisateur de fournir une matrice de transformation permettant une transformée en cosinus discrète (DCT) ou une transformée de Karhunen-Loeve (KLT). Cela est fait en utilisant l’option -Mlin du VM. Une autre option (-Mtdt) permet d’appliquer une transformée en ondelettes (DWT 9/7) entre les différentes composantes. Dans ces deux cas, la transformation 1D est appliquée sur la dimension spectrale avant d’utiliser l’encodage classique JPEG2000 sur les images résultantes. Une optimisation débit-distorsion par lagrangien est une option utile car les composantes ont des propriétés statistiques différentes après la compaction d’énergie réalisée par la transformation (option -Flra). Sans cette option, le VM a du mal à atteindre le débit visé sur certaines images.

La figure 3.1 compare les performances de plusieurs transformations intercomposantes : DCT, KLT et DWT en terme de PSNR (Eq. 2.11) en fonction du nombre de bits par pixel par bande (bpppb).


PIC
Fig. 3.1: Performances de la compression JPEG 2000 en utilisant de la décorrélation interbande par DCT, KLT, DWT ou sans décorrélation.


Introduire une décorrélation intercomposante des données avant le codage JPEG 2000 permet d’augmenter de manière importante les performances. La KLT donne les meilleurs résultats car elle est adaptée aux données et est optimale en terme de décorrélation et de concentration d’énergie. La DWT augmente significativement les résultats tandis que la DCT présente des performances plus faibles.

La KLT est spécifique à chaque image et les matrices de transformation doivent donc être recalculées en fonction des données. Ce calcul est coûteux en terme de complexité : il est nécessaire de calculer la matrice de covariance ainsi que les vecteurs propres. Cette complexité est dissuasive pour la plupart des applications particulièrement lorsque le nombre de bandes augmente. Dans le cas des images multispectrales une approche utilisant une KLT générique, calculée sur plusieurs images est envisagée [Thi06]. Dans ce cas, il n’est plus nécessaire d’effectuer la diagonalisation de la matrice de covariance, il suffit d’effectuer la transformation. Avec 224 bandes, l’opération reste encore impossible en implémentation matérielle comme le montre le tableau 3.1. L’approche par transformation en ondelettes paraît plus adaptée, surtout dans le cas hyperspectral où le nombre de composantes est important.




# bandesSurface (mm2)


16 6
64 99
128 400
256 1500



Tab. 3.1: Surface de silicium nécessaire pour l’implémentation matérielle (ASIC) de la transformée KLT (sans le calcul de la matrice de transformation) en fonction du nombre de bandes spectrales. La limite actuelle est de l’ordre de 110mm2 [Alc06].

L’algorithme JPEG 2000 nécessite des ressources importantes (principalement pour le codeur arithmétique et l’allocation de débit) ce qui n’est pas possible particulièrement dans le contexte d’une implémentation spatiale où les contraintes sont très fortes. A notre connaissance, il n’existe qu’une seule implémentation de JPEG 2000 adaptée aux contraintes spatiales [VB05]. L’auteur de cette implémentation indique que les performances sont significativement moins bonnes que Kakadu dans le cas de la compression avec pertes à cause de l’impossibilité d’implémenter l’optimisation débit-distorsion. Le Consultative Committee for Space Data Systems (CCSDS), un groupe de travail regroupant les principales agences spatiales mondiales (NASA, JAXA, ESA, CNES) a émis des recommandations pour les systèmes de compression bord [Yeh05]. La recommandation consiste à regrouper les coefficients de la transformée en ondelettes dans une structure similaire à celle des arbres de zéros plutôt que d’utiliser la norme JPEG 2000. Le but de cette implémentation est d’alléger la complexité vis-à-vis de la norme JPEG 2000. Les performances de JPEG 2000 dans notre étude sont à voir plutôt comme une référence à atteindre, le but étant ici d’obtenir des performances proches de cette référence tout en gardant une complexité faible. La version de JPEG 2000 au fil de l’eau (scan-based mode), plus réaliste pour une implémentation spatiale donne elle-même des performances réduites par rapport à la référence.

La réduction de complexité par rapport à JPEG 2000 n’est pas obtenue directement par l’utilisation de la transformation en ondelettes selon les 3 directions, mais il apparaît que les ondelettes permettent l’utilisation d’outils performants et simples. C’est donc la transformée en ondelettes qui est choisie pour la décomposition.

3.2 Choix de la décomposition optimale

Avant d’adapter les arbres de zéros aux images hyperspectrales, il est nécessaire de définir une extension de la transformation en ondelettes performante pour les images hyperspectrales. La plupart des extensions actuelles sont basées sur des décompositions isotropiques [Tan03,Kim03], or, comme il a été montré dans les chapitres précédents, les données hyperspectrales ne sont clairement pas isotropiques. Dans le domaine de la compression vidéo, des structures anisotropiques sont utilisées avec succès [He03,Cho03]. Cependant, aucune justification théorique n’a été donnée concernant le choix de cette structure particulière et des décompositions plus efficaces pourraient exister. De toute façon, le meilleur choix pour les données vidéo n’est pas nécessairement le meilleur pour les images hyperspectrales à cause des propriétés statistiques différentes.

Le problème de la recherche de la décomposition en ondelettes optimale pour signaux à une dimension a été exploré dans plusieurs publications (par exemple [Coi90]). Pour les images naturelles 2D, les possibilités de décomposition ont souvent été restreintes à des quadtrees (conduisant à des sous-bandes carrées) mais ont évolué avec les décompositions anisotropiques [Xu03]. Plusieurs critères ont été utilisés pour choisir la décomposition optimale : sélection basée sur l’entropie [Coi92] ou sur un compromis débit-distorsion [Ram93] par exemple. L’avantage principal du dernier choix est qu’il propose simultanément l’allocation de débit entre les différentes sous-bandes [Sho88].

3.2.1 Décomposition anisotropique 3D en ondelettes

Classiquement, pour les images 2D, la transformée en ondelettes est isotropique i.e. pour une sous-bande donnée, le niveau de décomposition dans la direction horizontale est le même que dans la direction verticale. Cette alternance entre décomposition des lignes et des colonnes conduit à des sous-bandes carrées 1, l’équivalent en 3D étant des cubes. C’est le cas de la décomposition multirésolution définie par Mallat ou de la décomposition en paquets d’ondelettes pour les images [Mal89]. Un exemple de la décomposition multirésolution classique est illustré sur la figure 3.2. Les coefficients gris représentent des valeurs proches de zéros, tandis que les noirs et blancs correspondent respectivement à des valeurs fortement négatives ou positives.


PIC


Fig. 3.2: Décomposition multirésolution classique : on décompose successivement suivant les lignes et les colonnes. LF est la partie basse fréquence (LF = low frequency) et HF la partie haute fréquence (HF = high frequency).


Le terme anisotropique est plus général que l’utilisation courante du terme paquets d’ondelettes. Dans la plupart des cas, le terme paquets d’ondelettes désigne une transformation en quadtree conduisant à des sous-bandes carrées. Cette utilisation est justifiée par le fait que les images 2D classiques ont des propriétés statistiques similaires dans toutes les directions. Un exemple de décomposition anisotropique est donné sur la figure 3.3. On s’aperçoit ainsi qu’une telle décomposition non isotropique permet d’isoler une sous bande contenant beaucoup d’énergie (Fig. 3.3, sous-bande 1), tandis que les autres sous-bandes contiennent une majorité de coefficients nuls (la sous-bande 2 par exemple). Ceci permettra un codage efficace de l’image.


PIC


Fig. 3.3: Exemple de décomposition anisotropique.


On note Wi,j,kp,q,r une sous-bande de la décomposition 3D en ondelettes (Fig. 3.4) :

Cette notation est illustrée sur la décomposition multirésolution pour une image 2D dans l’annexe A.

Une relation peut être définie au niveau des espaces vectoriels des sous-bandes. Pour une décomposition selon les lignes, l’espace des ondelettes anisotropiques vérifie
W pi,,qj,,kr= W 2ip+,1q,,jr,k ⊕ W 2i+p1+,1j,q,k,r
(3.1)

est la somme directe de deux espaces vectoriels.

Pour une décomposition selon les colonnes, on a
W  p,q,r = W p,2q,r  ⊕ W p,2q+1,r.
  i,j,k     i,j+1,k     i,j+1,k
(3.2)

et selon les spectres
W p,q,r = W p,q,2r ⊕ W p,q,2r+1.
  i,j,k     i,j,k+1     i,j,k+1
(3.3)

À chaque étape de la décomposition, pour toutes les sous-bandes, il est possible de choisir la direction de la décomposition ce qui accroît la flexibilité de l’espace des décompositions. La décomposition multirésolution et les décompositions en paquets d’ondelettes sont toutes deux des cas particuliers de cette représentation.


PIC


Fig. 3.4: Décomposition anisotropique et notations.


Une représentation simple et bien adaptée pour ce type de décomposition peut être faite sous la forme d’arbre (Fig. 3.5 et 3.6). Chaque nœud porte un numéro, indiquant le sens de la décomposition (x, y ou λ), et deux branches, correspondant aux deux sous-bandes.


PIC PIC
Fig. 3.5: Exemple de décomposition en ondelettes d’une image 2D classique (décomposition de Mallat). La décomposition est itérée sur la bande basse fréquence (en haut à gauche). L’arbre associé correspondant à cette décomposition est représenté à droite.



PIC PIC


Fig. 3.6: Exemple de décomposition non isotropique d’une image hyperspectrale et son arbre de décomposition associé.


Il faut noter que la correspondance arbre-décomposition n’est pas bijective. Deux arbres différents peuvent donner la même décomposition. Par exemple, décomposer sur x puis sur chacune des deux sous-bandes décomposer sur y donne la même décomposition que décomposer d’abord sur y puis sur x pour chacune des deux sous-bandes. Les arbres correspondant à ces deux transformations seront différents. Par contre, un arbre donne bien une seule décomposition possible.

La représentation sous forme d’arbre se prête particulièrement bien à une programmation récursive.

3.2.2 Optimisation débit-distorsion

3.2.2.1 Le problème d’allocation

Le problème de l’allocation de débit, i.e. la distribution du budget de bits entre chaque sous-bande est un problème classique en compression de données. Shoham et Gersho ont traité le problème dans le cadre de la théorie débit-distorsion [Sho88]. Leur solution consiste à minimiser la distorsion sous une contrainte de débit en utilisant la méthode du lagrangien.

Dans le contexte de la décomposition en ondelettes, différentes types de quantificateurs peuvent être utilisés pour les différentes sous-bandes. Chaque quantificateur donnera une valeur de débit et une valeur de distortion pour chaque sous-bande (Fig. 3.7). On appelle S l’ensemble fini des combinaisons de quantificateurs pour les sous-bandes et B un élément de S. Un choix de B indiquera le quantificateur utilisé pour chacune des sous-bandes de la décomposition. Le problème est de minimiser la distorsion D(B) sous la contrainte de débit total R(B) dans le budget total Rc :
min{D (B )} sous R(B ) ≤ Rc.
B∈S
(3.4)

En utilisant la méthode du lagrangien, on transforme cette minimisation sous contrainte en minimisation d’une fonction de coût J sans contrainte mais avec un paramètre λJ
J(λJ) = D + λJR.
(3.5)

Sous une hypothèse de codage indépendant des différentes sous-bandes et d’additivité pour les mesures de distorsion et de débit, on montre que l’optimal global est atteint lorsque toutes les sous-bandes sont à un même point de fonctionnement λJ pour leur courbe débit-distorsion. Le problème devient alors :
min{D   + λ R  } pour chaque sous-bande k.
      k    J  k
(3.6)

La preuve de l’équivalence entre le problème sous contrainte et le problème sans contrainte est simple et peut être trouvée dans [Sho88].


PIC


Fig. 3.7: Calcul de la courbe débit-distorsion (R-D) pour une sous-bande. Pour différents pas de quantification, la distorsion et le débit sont mesurés conduisant à différents points de fonctionnement pour former la courbe débit-distorsion.


3.2.2.2 Algorithme

Le débit Rkq pour chaque sous-bande est évalué en utilisant le codeur arithmétique défini dans [Mof98]. Le choix du codeur n’est pas critique ici. Ce sont les positions relatives des différents choix possibles qui sont importants, pas leur performances absolues. Des simulations avec d’autres mesures de débits comme l’entropie des coefficients des sous-bandes ou par une combinaison de run length coding et de codage de Rice conduisent à des résultats similaires. Le fait que les sous-bandes sont codées de manière indépendante est une supposition implicite ici.

La courbe débit-distorsion est calculée pour les 3 décompositions suivantes possibles (correspondant aux 3 directions). Une représentation illustrée sur la figure 3.8 est obtenue. Pour chaque valeur de λJ, la fonction de coût J est calculée pour chacun des points de la courbe débit-distorsion. La décision de poursuivre ou non la décomposition est prise selon le coût le plus faible.


PIC


Fig. 3.8: Illustration de la décision de poursuivre ou non la décomposition. Pour plus de clarté, uniquement deux courbes débit-distorsion (sur les 4) sont représentées. Avec le λJ présenté ici, la décision de poursuivre la décomposition selon la direction x sera prise.


La recherche de type bottom-up est basée sur une fonction récursive et sur la propriété d’additivité de la fonction de coût J. On note Oi,j,kp,q,r la base optimale de Wi,j,kp,q,r. Soit Bi,j,kp,q,r la base de Wi,j,kp,q,r sans aucune transformation. On a alors une fonction de coût pour une base de représentation, J(λJ,B) = min{D + λJR} : R et D étant les points de fonctionnement de la sous-bande représentée sur la base B.

J est une fonction de coût additive : pour deux bases orthogonales B1 et B2, on a
J (λJ ,B1 ∪ B2 ) = J (λJ ,B1)+ J (λJ,B2)
(3.7)

en utilisant les propriétés d’additivité des débits et des distorsions.

Proposition 1 (Cas 1D : Coifman, Wickerhauser). Si J est une fonction additive de coût alors
      {   2p     2p+1            2p           2p+1          p
Op  =   O i+p1 ∪ O i+1   si J(λJ,O i+21p) + J(λJ,O i+21p+1) < J (λJ ,Bip)
  i     B i           si J(λJ,O i+1) + J(λJ,O i+1  ) ≥ J (λJ ,Bi)
(3.8)

   Preuve: La meilleure base Oip est soit égale à Bip ou à l’union B0 B1 des bases de Wi+12p et Wi+12p+1 respectivement. Dans ce dernier cas, la propriété d’additivité (Eq. 3.7) implique que le coût sur Oip est minimum si B0 et B1 minimisent le coût sur Wi+12p et Wi+12p+1. On a donc B0 = Oi+12p et B1 = Oi+12p+1. Cela montre que Oip est soit Bip, soit Oi+12p Oi+12p+1. La meilleure base est choisie en comparant le coût des deux possibilités.  __

Ce théorème, formulé à une dimension, peut être étendu au cas à 3 dimensions. A chaque étape on a alors 4 possibilités pour la meilleure base :

Proposition 2 (Extension à 3 dimensions). Si J est une fonction additive de coût alors
        (| O2p,q,r  ∪ O2p+1,q,r  si J =  min{J ,J ,J ,J }
        |||{   ip+,21q,j,r,k    ip+,12q,j+,k1,r     1        0  1  2  3
Op,q,r =   O i,j+1,k ∪ Oi,j+1,k   si J2 = min{J0,J1,J2,J3}
  i,j,k   ||| Op,iq,j,,k2r+1 ∪ Opi,,jq,,k2+r+11   si J3 = min{J0,J1,J2,J3}
        |( Bp,q,r               sinon
            i,j,k
(3.9)

En notant

               2p,q,r            2p+1,q,r
J1  =   J(λJ,O i+1,j,k)+  J(λJ,O i+1,j,k )
J   =   J(λ ,Op,2q,r ) + J(λ ,Op,2q+1,r)
  2        J   i,j+1,k       J   i,j+1,k
J3  =   J(λJ,Op,iq,j,,k2r+1)+  J(λJ,Op,i,qj,,2kr++11)
J   =   J(λ ,Bp,q,r)
  0        J   i,j,k

En pratique, l’algorithme est développé de manière récursive. On note J0 le coût de la sous-bande courante sans réaliser de décomposition supplémentaire. Les coûts J1, J2 et J3 sont les coûts correspondant à une poursuite de la décomposition de la sous-bande courante respectivement selon les directions x, y ou λ. Les coûts J1, J2 et J3 sont calculés par un appel récursif à la fonction de calcul.

Algorithme 1. Recherche de la décomposition optimale

Fonction récursive : cost(Wi,j,kp,q,rJ)

Fonction globale

Cet algorithme conduit à une décomposition optimale différente pour chaque image et pour chaque débit visé. Il est à noter que grâce à l’utilisation de la récursivité dans l’algorithme, les points de fonctionnement seront d’abord calculés pour les sous-bandes les plus petites et ensuite l’algorithme rassemblera ces valeurs pour prendre la décision de partager ou non la sous-bande. Cette recherche est exhaustive et ne conduit en aucun cas à un minimum local.

Cette recherche est similaire dans l’idée à ce qui est fait par Ramchandran dans [Ram93] avec une extension à un espace 3D anisotropique.

3.2.3 Résultats sur les images 2D

Pour illustration, la recherche de la décomposition optimale est d’abord appliquée à des images 2D classiques. Le débit est ici estimé avec un codeur arithmétique. Pour la plupart des images (comme Lena, Fig. 3.9), le gain apporté par la décomposition optimale ne justifie pas l’augmentation de complexité. Ce résultat est en accord avec le fait que la décomposition multirésolution est largement utilisée dans les standards de compression tel que JPEG 2000. Cependant, lorsque l’image a un contenu fréquentiel plus marqué (comme Barbara, Fig. 3.10), la décomposition optimale arrive à grouper les coefficients de forte amplitude dans une sous-bande (Fig. 3.10 sous-bande 1) tandis qu’une sous-bande de taille importante (Fig. 3.10 sous-bande 2) arrive à regrouper un grand nombre de coefficients de faible amplitude. Dans ce cas, le gain peut atteindre 1.5 dB par rapport à la décomposition multirésolution (pour des débits compris entre 0.2 et 2.0 bpp).


PIC PIC PIC


Fig. 3.9: Image Lena et la décomposition optimale obtenue pour un débit de 0.5 bit par pixel (bpp) (à gauche). Courbe débit-distorsion pour la décomposition classique et pour la décomposition optimale (droite). Dans le cas d’une image classique, la décomposition optimale n’apporte pas de gain significatif.

PIC PIC PIC


Fig. 3.10: Image Barbara et la décomposition optimale obtenue pour un débit de 0.9 bit par pixel (bpp) (à gauche). Courbe débit-distorsion pour la décomposition classique et pour la décomposition optimale (droite). Dans le cas d’une image avec des fréquences spatiales marquées, la décomposition optimale apporte un gain significatif.


3.2.4 Décomposition en ondelettes optimale pour l’hyperspectral

Les résultats présentés sur la figure 3.11 montrent que la décomposition anisotropique optimale amène une amélioration claire en terme de débit-distorsion. Ces résultats sont confirmés sur diverses images hyperspectrales. L’amélioration est d’environ 8 dB en terme de PSNR par rapport à une décomposition isotropique classique, pour des débits compris entre 0.1 et 4.0 bits par pixel par bande (bpppb). Si la comparaison est faite en terme de contrainte de qualité, par exemple pour un PSNR supérieur à 70 dB, le débit nécessaire passe de 1 bpppb à 0.5 bpppb, soit une réduction d’un facteur 2. L’importance de cette amélioration par rapport aux images classiques est due à la nature anisotropique des images hyperspectrales.


PIC
Fig. 3.11: Décomposition optimale sur des données hyperspectrales : la décomposition anisotropique optimale améliore clairement les performances (environ 8 dB), mais au prix d’une complexité importante.



PIC (a) Moffett4
PIC (b) Moffett3
PIC (c) Harvard1
PIC (d) Hawaii1
PIC (e) Hawaii2
PIC (f) Hyperion1

Fig. 3.12: Différentes images hyperspectrales utilisées pour les simulations. (a) et (b) sont différentes parties de la scène f970620t01p02_r03 du capteur AVIRIS sur le site de Moffett Field. (a) contient des zones uniformes avec des caractéristiques spectrales marquées. (b) est un mélange de ville (hautes fréquences spatiales). (c) provient de la scène f010903t01p01_r03 d’AVIRIS au dessus de Harvard Forest, cette image contient principalement de la végétation. (d) et (e) viennent de la scène f000414t01p03_r08 à Hawaï. (d) contient deux zones contrastées (minérale et océan) séparées par une ligne de côte. (e) a été sélectionnée pour illustrer une scène contenant des nuages, une forte dynamique et des contrastes importants. (f) provient du capteur spatial Hyperion (EO1H0440342002212110PY) et couvre une zone similaire à celle des images (a) et (b).


3.2.5 Décomposition fixe

Il y a deux inconvénients principaux à cette décomposition optimale :

Le coût calculatoire de la méthode est important. Par exemple, la recherche de la décomposition optimale sur un cube hyperspectral de dimension 256 × 256 × 224 avec une taille minimale de sous-bande de 8 × 8 × 7 (un maximum de 5 niveaux de décomposition) nécessite de calculer la courbe débit-distorsion totale pour 250047 sous-bandes. Cette valeur est obtenue en calculant les combinaisons de toutes les tailles de sous-bandes possibles dans les trois directions : pour une direction, on a 20 + 21 + + 25 = 63 tailles possibles, en combinant dans les 3 dimensions, on a 633 = 250047 sous-bandes possibles.

La dépendance de la décomposition par rapport à l’image et au débit pose des problèmes en terme d’implémentation de la transformation. En général, on préfère des transformations indépendantes des données.

Le but est donc de trouver une décomposition suffisamment proche de la décomposition optimale pour donner des performances quasi-optimales, tout en conservant une forme assez générale qui lui permette de rester valable pour une grande variété d’images. La structure de la décomposition optimale pour différentes images à différents débits, ainsi que la localisation de la corrélation résiduelle pour la décomposition isotropique conduit à essayer une décomposition particulière. Cette décomposition régulière consiste à appliquer une décomposition multirésolution standard sur les spectres suivie d’une décomposition multirésolution 2D sur les images résultantes (l’ordre des opérations étant réversible). Cette décomposition a montré de bonnes performances dans plusieurs travaux [Kim00,Wan04a], mais uniquement avec une justification empirique. Une justification théorique en terme d’entropie apparait dans [Pen06].

La décomposition obtenue est présentée sur la figure 3.13. Cette décomposition est comparée avec la décomposition isotropique classique, les coefficients en gris représentent les valeurs proches de 0, les coefficients en blanc sont les coefficients positifs tandis qu’en noir sont les coefficients négatifs. Comme on peut le voir sur la figure 3.14, cette décomposition fixe est presque aussi performante que la décomposition optimale. La régularité de cette structure ainsi que les faibles perspectives de gain possibles à chercher une autre décomposition conduisent à choisir cette décomposition.


PIC (a)
PIC (b)

Fig. 3.13: Décomposition isotropique classique (a) et décomposition anisotropique (b) avec 3 niveaux de décomposition (les simulations sont faites avec 5 niveaux). Pour la décomposition isotropique, on peut remarquer qu’il reste une corrélation importante entre les coefficients dans les basses fréquences spectrales et spatiales (lignes de coefficients consécutifs de valeur similaire). Cette corrélation résiduelle est plus faible pour la décomposition anisotropique.


PIC (a) moffett4
PIC (b) moffett3
PIC (c) Harvard1
PIC (d) Hawaii1
PIC (e) Hawaii2
PIC (f) hyperion1

Fig. 3.14: Résultats sur des données hyperspectrales : la décomposition anisotropique optimale améliore clairement les performances. La décomposition anisotropique fixe permet d’obtenir des performances proches de la décomposition optimale tout en gardant une complexité raisonnable.


3.3 Structures d’arbres

3.3.1 Idées principales

Une des faiblesses possibles de JPEG 2000 est qu’il n’utilise pas les relations existantes entre les localisations des coefficients significatifs entre les différentes sous-bandes. D’après Taubman [Tau02], le bénéfice venant du choix du point de troncature compense le fait que les relations entre les coefficients ne sont pas utilisées. Cette conclusion peut être différente dans le cas des images hyperspectrales, la corrélation entre les sous-bandes étant anormalement élevée.

Les arbres de zéros sur les coefficients d’ondelettes ont été développés pour utiliser la relation qui existe entre la position des singularités dans les différentes sous-bandes. Après la transformée en ondelettes, on observe que la localisation des coefficients significatifs est similaire entre les différentes sous-bandes même si leurs amplitudes sont décorrélées (Fig. 3.15). La propriété à l’origine des arbres de zéros est que si un coefficient n’est pas significatif dans une sous-bande alors le coefficient à la même position dans une sous-bande de plus haute fréquence sera lui aussi probablement non significatif.


PIC


Fig. 3.15: Illustration du principe des arbres de zéros : même si les sous-bandes sont décorrélées (les coefficients de corrélation entre les bandes HL, LH et HH sont très faibles), les endroits où les singularités apparaissent sont les mêmes dans les différentes sous-bandes.


Cette idée a été exploitée avec succès par Shapiro avec EZW (Embedded Zerotree of Wavelet coefficients) [Sha93]. Une amélioration remarquable a été apportée quelques années plus tard par Said et Pearlman avec SPIHT (Set Partitioning In Hierarchical Trees) [Sai96]. Toutefois, les performances des arbres de zéros sont supérieures à ce qu’on pouvait en attendre [Mar01] et ce succès n’est pas complètement expliqué. L’idée de SPIHT a ensuite été reprise pour donner des solutions progressives en résolution [Dan03], résistantes aux erreurs de transmission [Cho05a] ou généralisées pour des structures 3D, principalement pour la compression vidéo [Kim00,He03] ou pour les images médicales [Bil00,Xio03]. Sur le modèle de SPIHT, d’autres codeurs ont ensuite vu le jour comme SPECK [Isl99,Pea04] qui tire parti des similarités entre les coefficients adjacents dans une même sous-bande (blocs de zéros plutôt que des arbres de zéros). SPECK a ensuite été adapté en SBHP [Chr00] pour concurrencer EBCOT [Tau00] dans la norme JPEG 2000 ainsi qu’en EZBC [Hsi01] utilisant un codage arithmétique contextuel évolué pour augmenter les performances. LTW [Oli03], amélioré dans [Guo06] propose une approche différente pour réduire la complexité. Plus récemment, SPECK a été adapté pour fournir un codage progressif en résolution [Xie05].

3.3.2 Principes généraux de EZW et SPIHT

Au moment de sa publication, l’algorithme EZW de Shapiro, utilisant les arbres de zéros, produisait les meilleures performances en terme de débit-distorsion tout en ayant une complexité faible. Cette structure d’arbres de zéros a été généralisée quelques années après par Said et Pearlman pour donner l’algorithme SPIHT.

Ces deux algorithmes possèdent des propriétés qui les rendent particulièrement attractifs dans le domaine de la compression spatiale embarquée. Ils produisent tous deux un train binaire emboîté : le préfixe d’un train binaire produit par EZW (ou SPIHT) est lui-même un train binaire EZW (ou SPIHT) valide conduisant à une image décompressée avec une qualité plus faible. Ces deux algorithmes atteignent ce résultat avec un niveau modeste de complexité.

La question de la complexité algorithmique est un point critique pour la définition des algorithmes embarqués, la puissance de calcul étant très limitée. Cela est particulièrement vrai dans le cas des satellites. Le fait d’avoir un train binaire emboîté est intéressant pour éviter les dépassements des mémoires et garantir la meilleure qualité d’image étant données les ressources disponibles à bord (mémoire et calcul).

Pour assurer cette propriété de train binaire emboîté, la compression procède par plans de bits. On note cx,y,λ le coefficient de la transformée en ondelettes de la colonne x, de la ligne y et du plan spectral λ. On définit la suite de seuils TK-1,,T0, tels que Tk = Tk+12. Le seuil initial est choisi pour que |cx,y,λ| < 2TK-1 pour tous les coefficients d’ondelettes. Pour des raisons pratiques de représentation binaire, on choisit TK-1 comme une puissance de 2. Le coefficient cx,y,λ est significatif au plan de bits k si |cx,y,λ|≥ Tk. Les plans de bits sont codés les uns après les autres permettant de réduire la distorsion à chaque étape.

Un exemple de codage par plan de bits est présenté dans le tableau 3.2. On considère des coefficients d’ondelettes de valeurs 7, 30 et 180. 7 n’est pas significatif avant le plan de bits 2, 30 avant le 4 et 180 le devient dès le plan de bits 7. Après qu’un coefficient ait été marqué comme significatif, les bits suivants doivent être codés. Le signe des coefficients est codé séparément.

Pour chaque plan de bits, le codage suit une structure d’arbres de zéros, les arbres étant définis le long des différentes sous-bandes de la transformée en ondelettes.







Plan de bitsTk730180





7 1280 0 1
6 64 0 0 0
5 32 0 0 1
4 16 0 1 1
3 8 0 1 0
2 4 1 1 1
1 2 1 1 0
0 1 1 0 0






Tab. 3.2: Exemple de codage en plan de bits.

La définition d’un arbre de zéros varie en fonction des algorithmes. Dans [Cho05b] est développée l’idée des arbres de zéros de degré k. Un arbre de degré 0 est un arbre dont tous les coefficients sont égaux à zéros (Fig. 3.16), un arbre de degré 1 est un arbre dont tous les coefficients sauf la racine sont égaux à zéros (Fig. 3.17) et un arbre de degré 2 est un arbre dont tous les coefficients sauf la racine et ses enfants directs sont égaux à zéros (Fig. 3.18). EZW utilise des arbres de degré 0 tandis que SPIHT utilise des arbres de degré 1 et 2.


PIC
Fig. 3.16: Arbre de zéros de degré 0

PIC
Fig. 3.17: Arbre de zéros de degré 1

PIC
Fig. 3.18: Arbre de zéros de degré 2


Un exemple de la première passe des algorithmes EZW et SPIHT est détaillé dans l’annexe B.

3.3.3 Étude statistique pour le choix de la structure d’arbre

Étant donné la décomposition optimale précédente, plusieurs structures d’arbre peuvent être définies. On peut utiliser le lien entre les sous-bandes spatiales, entre les sous-bandes spectrales ou les deux. L’avantage principal des arbres de zéros est leur capacité de coder un grand nombre de coefficients à zéro (pour le plan de bits courant) en n’utilisant que très peu de symboles. La structure d’arbre optimale est donc celle qui maximise la longueur des arbres de zéros tout en ne laissant qu’un nombre très faible de zéros isolés. Quelle que soit la structure choisie, les coefficients significatifs seront les mêmes. La différence dans l’efficacité des structures d’arbres est uniquement due à leur capacité à rassembler les coefficients nuls.

Il y a trois structures régulières d’arbre possibles :

Le choix optimal n’est pas évident : mettre plus de coefficients dans un arbre permet d’encoder potentiellement plus de zéros avec un seul symbole, mais augmente également le risque d’avoir un coefficient significatif détruisant l’arbre. Des choix ont été fait précédemment mais sans justifications comme dans [He03].

Des statistiques sont calculées pour la transformée d’une image hyperspectrale avec 5 niveaux de décomposition et sont présentées dans les tableaux 3.3, 3.4 et 3.5. Ces tableaux donnent le nombre de coefficients significatifs pour chaque plan de bits (qui doivent être codés de toute façon et sont indépendants de la structure choisie), le nombre de zéros isolés (IZ) ainsi que le nombre d’arbres de zéros (ZTR). On considère ici les arbres de degré 0. Pour un plan de bits donné k, tous les coefficients en dessous du seuil Tk sont à zéro. Le but du codage par arbre de zéros est d’inclure le maximum de coefficients dans un arbre de 0 en utilisant un seul symbole (ZTR dans la terminologie EZW comme on le verra plus tard). Les coefficients en dessous du seuil ne pouvant pas être inclus dans un arbre de zéros sont codés avec un symbole (IZ). Les deux symboles, ZTR et IZ, sont utilisés pour coder les coefficients à zéros. Le but de la structure d’arbre est de minimiser le nombre de symboles utilisés ou de maximiser le nombre de 0 codés avec un seul symbole. Pour comparer l’efficacité des différentes structures, le nombre moyen de 0 codés par un symbole (ZTR ou IZ) est donné ci-dessous dans les trois tableaux 3.3, 3.4, 3.5.

Dans ces tableaux le terme moyenne désigne le nombre de 0 codés en moyenne par un des symboles IZ ou ZTR. Cette valeur est obtenue par :
moyenne =  nbr coeff---nbr-coeff-significatifs
                 nbr IZ + nbr ZTR
(3.10)

Par exemple, pour une structure d’arbre 3D (tab. 3.3), pour le plan de bits 16, 610 symboles sont significatifs, les autres ((256*256*224)-610=14679454) sont codés par 94+1960 = 2054 symboles, ce qui donne une moyenne de 7146.76 zéros codés par symbole. Cette valeur est à comparer avec l’arbre spatial pour lequel 15745 symboles sont nécessaires (932.32 zéros codés par symbole).







Plan de bitsSignificatif IZ ZTRMoyenne





19 99 0 74519704.65
18 183 0 99714724.05
17 249 2 119512263.84
16 610 94 1960 7146.76
15 1143 264 3595 3803.81
14 2975 975 9619 1385.41
13 9052 2678 28191 475.27
12 25525 6797 66564 199.76
11 64646 14587 122718 106.44
10 141869 29824 189039 66.43
9 276589 56712 268876 44.24
8 499196 104615 372197 29.74
7 860482 189221 517928 19.54
6 1461127 376246 718656 12.07
5 2439789 691967 926441 7.56
4 394055512993991015978 4.64
3 605323721303171098716 2.67
2 89501662860847 800219 1.57
1 117959232079228 341052 1.19
0 13673387 777902 156241 1.08






Tab. 3.3: Statistiques pour une structure d’arbre 3D.







Plan de bitsSignificatif IZ ZTRMoyenne





19 99 0 14534 1010.04
18 183 0 14702 998.50
17 249 2 14842 988.94
16 610 41 15704 932.32
15 1143 294 17256 836.41
14 2975 1877 22564 600.51
13 9052 10995 35115 318.17
12 25525 39488 53933