Calcul de la Transformée de RADON
Sommaire -Introduction générale …………………………………………………… …..1 -Chapitre I : définition du problème……… …………………………………. 2 I .1-Introduction………………………………………………………… 3 I .2-Principe de la transformée de RADON………………………… …. 3 I .3-Problème posé……………………………………………………… 4 I .3.1-Pourquoi choisir de la transformée de RADON I .4.2-Problème -Chapitre II : Approche mathématique…………………………………… …5 II. 1-Introduction au théorème de RADON ………………………… …5 II. 2-Définition de la TR d’une fonction à deux variables……………. ..6 II. 3-Propriétés de la TR……………………………………………..… .8 II .3.1-Linéarité II .3.2-Invariace à la translation II .3 .3-Invariance à la rotation II .3. 4-Changement d’échelle II .3. 5-Insensibilité aux bruits -Chapitre III :Formulation numérique…..…………………………………..…9 III. 1-Introduction et choix de la méthode…………………………….... ..9 III. 2-Discritisation et dimensionnement………………………………. . .9 III. 3-Résolution numérique…………………………………………….. 11 III .3.1-Optimisation du calcul……………………………… … 11 III .3.2-Optimisation sur la boucle ……………………………... .12 III .3.3-Considérations numériques..……………………………...13 III. 4-Algorithme de la transformée de RADON………………………. ..15 III. 5-Utilisation de l’approximation linéaire……………………………..17 -Chapitre IV : Résolution et Discussion………………………………………18 IV .1-Résultats obtenus pour différentes images.………………………...18 IV. 2-Influence du nombre de projections………………………………..19 IV. 3-Influence de la rotation ………………………………………….....21 IV .4-Comparaison avec les résultats obtenus sur Matlab……………… .22 IV .5-Conclusion ……………………..……………………………… .. .23
-Conclusion générale ……………………………………………………….. 24
2
-Bibliographie………………………………………………………………. .25 -Annexe………………………………………………………………………..26
Introduction générale : L’identification des cibles aériennes a toujours été l’une des préoccupations majeures de la surveillance aérienne. Les premiers systèmes développés pour identifier les cibles sont les ‘I FF’ (identification friend-foe) qui consiste en une série de questions /réponses entre la cible et la station de surveillance à travers des codes. La station prend alors des multi décisions sur la cible c’est-à-dire : cible coopérative ou hostile. Mais cette méthode s’est avérée périlleuse lorsque une cible hostile se fait passer pour un avion civil par exemple. Pour contourner ce problème, les recherches se sont orientées vers la reconnaissance des cibles sans leurs participations actives. Sachant que les ondes électromagnétiques se comportent suivant la géométrie de la cible et son mouvement de rotation, la signature radar peut nous fournir beaucoup d’informations sur une cible qu’elle soit coopérative ou non. La signature radar de la cible dépend directement de la surface vue par le radar qui est appelée SER (surface équivalente radar) qui représente le pouvoir réflecteur de la cible. L’une des difficultés de la SER est sa grande dépendance de l’angle de visée du radar. Les images ISAR (radar à ouverture synthétique inverse) représentent une distribution bidimensionnelle de la SER de la cible. Nous en parlons brièvement en disant qu’elle donne les informations concernant l’image en 2D et la puissance des centres de réflexion des ondes sur la cible .Mais la fiabilité de cette méthode nécessite des algorithmes sophistiqués et un grand espace mémoire pour stocker les images ISAR pour plusieurs cibles et divers angles d’orientation. C’est ici que commence vraiment notre problème qui consiste à exploiter les images ISAR de manière optimale. L’identification et la classification d’une cible non coopérative consistent à « comparer » l’image de la cible détectée à des images de référence stockées dans une banque de données. Mais l’orientation de la cible à classifier n’est pas connue à priori. Il convient donc de rendre l’image de la cible indépendante de l’angle d’observation. Ceci dit, un prétraitement des données, qui permet de réduire leur taille tout en conservant le maximum d’information devient nécessaire. Si on arrive à rendre l’image indépendante de l’angle d’observation, on pourra diminuer considérablement la taille des données. Afin d’atteindre ce but, nous allons utiliser la transformée de RADON qui fera l’objet de notre étude tout au long de ce travail.
3
Notre travail est constitué de quatre chapitres. Le premier chapitre « Définition du problème » introduit la transformée de RADON et la situe dans le problème. Le deuxième chapitre « approche mathématique » est consacré au formalisme mathématique concernant la TR, définition et propriétés. Le troisième chapitre « formulation numérique » consiste a construire un algorithme capable de calculer la transformée de RADON discrète. Enfin, dans le dernier chapitre, nous allons analyser les résultats obtenus, afin d’aboutir a des conclusions.
Chapitre I : Définition du problème : I.1. Introduction : La transformée de RADON a trouvé beaucoup d’application dans le domaine de la science et de l’ingénierie. Cette transformée a l’importante propriété de convertir des images 2-D en une série de projections 1-D. Cela permet de faciliter et d’accélérer énormément le traitement de l’information. De plus, la transformée de RADON d’un objet est invariante par rapport a la rotation et a la translation de ce dernier (voir chapitre 2), ce qui lui confère une place privilégiée dans le traitement d’images. Pour souligner l’importance de cette transformée nous allons citer quelques domaines d’applications : La tomographie a rayons X, la tomographie a émission de positrons, la tomographie ultrasonore, l’astronomie, la microscopie électronique, l’imagerie a résonance magnétique nucléaire et en géophysique.
I.2. Principe de la transformée de RADON : La transformée de RADON bidimensionnelle consiste, très brièvement, a faire des projections successives d’une fonction f(x,y) suivant différentes directions tel que schématisé dans la figure .1.
Figure.1 : représentation du domaine spatial et du domaine de RADON
Donc, pour chaque droite de direction donnée correspond une projection estimée par l’intégrale de la fonction sur la droite (voir définition, chapitre 2).Chaque rayon est indexé 4
par sa distance par rapport a l’origine et son angle. Une définition mathématique et plus rigoureuse sera donnée dans l’approche mathématique. Comme le voyons sur la figure 1, dans le domaine de RADON, une ligne est représentée par un point. Des caractéristiques géométriques de l’image originale produisent des pics ou des groupements de points dans le domaine de RADON. C’est ce qui constitue la force de cette transformée et sa large utilisation dans le traitement d’images.
I.3. Problème posé : I.3.1. Pourquoi choisir la transformée de RADON : La TR offre plusieurs avantages qui sont l’invariance à la rotation et à la translation, en plus de sa grande insensibilité par rapport aux bruits. Comme notre cible est sans cesse en mouvement, avec des changements climatiques permanents, nous avons pensé que la transformée de RADON était la mieux adaptée pour notre étude.
I.3.2. Problème : La base de données qui est en notre possession est constituée d’images de type ISAR. Qu’est-ce qu’une image ISAR ? C’est une image obtenue a partir de la SER d’un radar, qui fait apparaître les points les plus réfléchissants de la cible sous forme de points brillants. Voici représentées sur la figure 2 un exemple d’images ISAR d’avions de deux types différents. Chaque image est dimensionnée en pixels (dans notre cas : 96*96) et chaque pixel a un niveau de gris bien déterminé. Ceci permet de la représenter par une fonction a deux variables f(x, y) dont les valeurs seront stockées dans une matrice.
Figure 2 : Images ISAR de deux types d’avions Notre problème consiste a calculer la transformée de RADON de la fonction ainsi obtenue.
5
En entrée : Matrice g () représentant l’image a projeter. En sortie : Matrice Rg () représentant les différentes projections.
.Chapitre II
: Approche mathématique :
II.1. Introduction et Théorème de RADON : Depuis la mécanique céleste qui a conduit au développement du calcul différentiel, il est apparu que bon nombre de phénomènes naturels étaient régis par des équations différentielles, ceci justifie l’importance et l’étendue de cette théorie mathématique qui s’est alors développée pour elle-même. Voici un exemple de résultat appartenant a cette théorie. Rappelons d’abord que si f (t) est une fonction définie sur un intervalle [a, b] son intégrale : b
∫ a
f (t ) dt
représente l’aire comprise sous son graphe, et que le processus d’intégration est en quelque sorte inverse de la dérivation. _Considérons maintenant une fonction de deux variables f(x, y) définie sur un domaine du plan, par exemple un disque. A toute droite traversant ce disque, on associe l’intégrale de f sur l’intervalle constitué par l’intersection du disque et de la droite. Il s’agit donc de l’aire indiquée sur la figure 3.
6
Figure 3 : aire représentant l’intégrale de f sur un segment Supposons que la droite est d’équation y=A x+B, elle est donc déterminée par les nombres A et B, et on peut noter S (A, B) l’aire attachée a la droite. A la fonction f(x, y) on a ainsi associé une nouvelle fonction S (A, B).On peut alors se demander si connaissant la fonction S, on peut retrouver la fonction f. Il ne s’agit pas seulement d’une question de calcul pratique, mais il faut s’assurer que deux fonction f et g différentes ne peuvent pas donner la même fonction S, car sinon de S on ne saura pas s’il faut remonter a f ou a g. La réponse est oui, S détermine f .C’est le théorème de la transformée de RADON. Il a été démontré par RADON en 1917, simplement parce qu’il était intéressé par l’étude abstraite du calcul différentiel et intégral, sans souci des applications. Ce résultat qui était a priori purement théorique, a débouché sur une application inattendue : Le scanner médical. En effet, lorsqu’on fait une radiographie d’un corps, on envoie des rayons qui sont affaiblis lorsqu’ils rencontrent la matière, et la quantité soustraite mesure l’intégrale de la densité de matière sur le chemin parcouru.
II.2. Définition de la TR d’une fonction a deux variables : Soient (x, y) les coordonnées cartésiennes d’un point dans un espace a 2-D, g(x, y) une fonction continue a support compact, L un ensemble de droites dans ce plan défini par : ρ = x cos θ + y sin θ
(II.2.1)
La transformée de RADON est tout simplement donnée par l’intégrale sur la droite L : Rg ( x, y ) = ∫ g ( x, y ) dl
(II.2.2)
L
Avec le changement de variable suivant : x = ρ cos θ − q sin θ y = ρ sin θ + q cos θ
(II.2.3)
on obtient la représentation schématisée sur la figure 4, et l’équation mathématique de cette transformation devient : ∞ g ( ρ cos θ − q sin θ , ρ sin θ + q cos θ ) dq −∞
Rg ( ρ,θ ) = ∫
7
(II.2.4)
Figure 4 : ligne intégrale de RADON En introduisant un changement de variable approprié dans la formule précédente nous obtenons : Rg ( ρ,θ) =
∞ 1 ρ g ( x, − x cot θ) dx ∫ sin θ −∞ sin θ
(II.2.5)
En utilisant l’impulsion de Dirac, δ(.) (voir Annexe A), on obtient : Rg ( ρ,θ) =
∞ ∞
∫ ∫ g ( x, y )δ ( ρ − x cos θ − y sin θ)dxdy
−∞ −∞
(II.2.6)
L’intégrale étant évaluée sur la droite définie par l’équation II.2.1. Si maintenant, on considère l’équation de la droite L sous la forme, y = p * x +τ
(II.2.7)
On obtient la formule donnant la TR : Rg ( p,τ ) =
∞ ∞
∫ ∫ g ( x, y )δ ( y − px −τ)dxdy
−∞ −∞
(II.2.8)
Il existe d’autres définitions de la transformée de RADON (voir Annexe A), qui sont identiques et équivalentes aux précédentes.
II.3. Propriétés de la TR :
8
Dans cette partie nous allons exposer les propriétés de la TR pour mieux comprendre l’importance de cette transformée en particulier dans le traitement et la reconnaissance d’images.
II.3.1. Linéarité : R(∑α i g i ( x, y )) = ∑ α i g i ( ρ, θ ) i
(II.3.1)
i
II.3.2. Invariance a la translation : D’après la définition II.2.5 on a immédiatement :
g ( x − x0 , y − y0 ) → Rg ( ρ − x0 cosθ − y0 sinθ ,θ )
(II.3.2)
Donc, une translation de l’image de ( x0 , y0 ) se traduit par une translation ρ0 dans la direction de ρ . On dit alors que la TR est invariante.
II.3.3. Invariance a la rotation : Avec le même raisonnement on a : g ( x cos φ − y sin φ , x sin φ + y cos φ ) → Rg ( ρ ,θ + φ )
(II.3.3)
Une rotation de l’image d’un angle φ se traduit par une translation de φ dans la direction deθ. D’où l’invariance a la rotation de la TR.
II.3.4. Changement d’échelle : g (αx,αy ) =
1
α
Rg (αρ,θ)
(II.3.4)
Contrairement au changement à la translation et a la rotation, on remarque que la TR n’est pas invariante a l’échelle. Cependant l’information géométrique est toujours préservée.
II.3.5. Insensibilité aux bruits : Un autre grand avantage de la transformée de RADON est sa robustesse vis-à-vis des bruits blancs.
Chapitre III : Formulation numérique : III.1. Introduction et choix de la méthode : 9
Dans notre problème, nous allons appliquer la TR a des images qui sont représentées par des pixels avec différents niveau de gris. Donc la discrétisation dans ce cas est directement liée au nombre de pixels ainsi que la précision. D’où l’obligation de faire des approximations dues a la discontinuité des pixels et les résultats obtenus seront considérés du point de vue qualitatif que quantitatif. Ceci explique la méthode choisie qui a pour principe d’approcher par l’entier le plus proche les valeurs non entières L’image sera représentée par une fonction bidimensionnelle dont les valeurs seront stockées dans une matrice. Etant donné, d’après la définition II.2.7, que la transformée de RADON est calculée par l’intégrale suivant des droites, nous avons exposé une autre méthode qui utilise les même résultats que la première plus une interpolation linéaire .
III.2. Discrétisation et dimensionnement : L’image M * N pixels sera donnée par une matrice de dimension M * N représentant la fonction g ( x, y ) . Pour couvrir toute l’image, nous poserons :
0 <θ <π
et
− ρ min ≤ ρ ≤ ρ m ax
Discrétisons les paramètres :
x = xm = xmin + m∆ x y = y n = y min + n∆ y θ = θ k = θ min + k∆ θ ρ = ρ h = ρ m in + h∆ ρ
, , , ,
m = 0,1....., M −1 n = 0,1....., N −1
(III.2.2)
k = 0,1....., K −1 h = 0,1....., H −1
Le traitement se fera sur des images carrées (et c’est le cas pour les images ISAR), donc nous aurons : ∆x = ∆y
(III.2.3) (III.2.4)
M =N
De plus les images sont centrées d’où :
10
( M − 1) ∆x 2 ( M − 1) y min = xmin = − y max = − ∆x 2 ( H − 1) ρ min = − ρ max = − ∆ρ 2 xmin = − xmax = −
ainsi, discrétisons la fonction g ( x, y )
(III.2.5)
:
g ( m, n) = g ( xm , yn )
(III.2.6)
De même pour la transformée de RADON discrète on notera :
Rg ( h, k ) = Rg ( ρ h ,θ k ) Il est clair que
(III.2.7)
:
180 K = ∆θ
et θ min = 0
(III.2.8)
ou ∆θ est évalué en degré. le symbole [.] permet d’arrondir a l’entier le plus proche. Conformément a l’équation II.2.5 nous avons : ρ = x2 + y2
(III.2.9)
et puisque nos images comme on l’a dit sont centrées et de dimension M * M on a immédiatement : ρ max =
M 2
(III.2.10)
2ρ H = max + 1 ∆ρ
(III.2.11)
Il apparaît maintenant clairement que les paramètres dont nous aurons besoins se réduisent a : ∆θ, ∆ρ, ∆x, M .
III.3. Résolution numérique : D’après l’équation II.2.5 on a :
11
Rg ( ρh ,θ k ) =
∞ 1 g ( x, y = pk x +τ h )dx ∫ − ∞ sin θ k
(III.3.1)
avec
cosθ k sin θ k ρ τh = h sin θ k
pk = −
(III.3.2)
Approchons maintenant l’intégrale III.3.1 par la somme : Rg ( ρh ,θk ) =
avec :
∆x M −1 ∑g ( xm , yn ) sin θk m =0
(III.3.3)
y n = [ pk xm + τ h ]
(III.3.4)
Le symbole [.] signifie approcher a l’entier le plus proche. Maintenant utilisons les formules III.2.2 nous obtenons : p m∆ x + pk xmin + τ h − y min n= k ∆y
(III.3.5)
III.3.1-Optimisation de calcul : Pour minimiser les calculs dans notre algorithme nous poseront d’après III.3.5 :
n= [αm + β ] avec :
(III.3.6)
p ∆x α= k
∆y pk xmin + τ h − y min β= ∆y
(III.3.7)
D’où enfin notre somme III.3.3 prendra la forme Rg ( ρh ,θk ) =
:
∆x M −1 ∑g (m, [αm + β]) sin θ m =0
(III.3.8)
Des formules III.3.2 et III.3.6 nous déduisons les valeurs de α et β : ∆x cosθ k α =− ∆y sin θ k
12
(III.3.9)
ρ − x cosθ k − y min sinθ k β = h min ∆ y sinθ k
Maintenant sachant que : ∆x = ∆y et xmin = y min les formules III.3.9 se simplifient et prennent la forme :
α =−
cosθ k sin θ k
(III.3.10)
ρ − x min (cosθ k + sinθ k ) β= h ∆ x sinθ k
III.3.2-Optimisation sur la boucle : Ceci nous permettra de gagner dans le temps de la boucle, puisque nous savons que les valeurs de la fonction n’existe que pour : 0 ≤ n ≤ M − 1 ⇒ 0 ≤ [αm + β ] ≤ M − 1
D’où : β M −1 − β − ≤m ≤ , pour α > 0 α α M −1 − β β ≤m≤− , pour α < 0 α α pourα > 0 : β M −1− β mmin = max0, − et mmax = min M − 1, α α
Pour α < 0 : M −1− β β mmin = max0, Et mmax = min M − 1, − α α
(III.3.11)
Ainsi, nous avons pu diminuer de manière significative le nombre d’itérations dans une boucle.
III.3.3-Considérations numériques : Observons attentivement les formules III.3.9 donnant α et β , nous constatons que lorsque : θ k → 0 alors sin θ k → 0 , et ceci peut nous causer des problème du point de vue précision. Pour palier a ce problème nous allons procéder comme suit :
13
Pour θ k >
2 : nous considéreront les projections suivant l’axe des x , soit les droites : 2 y = pk x + τ h
Donc les formules à utiliser sont celles données par III.3.8 et III.3.10. Pour θ k >
2 : nous considéreront les projections suivant l’axe y, soit les droites : 2
x = pk′ y + τ h′
Et les formules utilisées sont obtenues de manière analogue en permutant introduisant quelque changements, ainsi on obtiendra les formules : Rg ( ρ h , θ k ) =
∆x cos θ
M −1
∑ g ([α ′n + β ′], n)
x
et y , et en
(III.3.12)
n =0
Avec les nouvelles valeurs de α ′ et β ′ données par : α′ =
1
α
(III.3.13)
ρ − x (cosθ k + sinθ k ) β ′ = h min ∆ x cosθ k Les valeurs de nmax et nmin sont données par : Pourα > 0 :
β ′ M − 1− β ′ nm in = m ax 0, − Et nm ax = m in M − 1, α ′ α ′ (III.3.14)
Pour α < 0 :
M − 1− β ′ β ′ nm in = m ax 0, n = m in M − 1 , − Et m ax ′ α α ′ Ainsi on aura contourné le problème et notre algorithme sera plus précis.
Dans notre problème, la fonction g a transformer n’est pas continue, on manipule des pixels, donc comme on l’a déjà dis c’est le coté information qui nous intéresse et comme ∆x .ne contient aucune information il sera négligés, c’est-à-dire nous poserons : ∆x = 1 . Dans le cas où on veut calculer la TR d’une fonction continue, autrement dit approximer quantitativement la TR, on doit prendre en compte ∆x car il représente la précision. Après ces simplifications, résumons les formules obtenues : 2 Pour θ k > : 2
14
m
max ∆x Rg ( ρh ,θk ) = g ( m, [αm + β ]) ∑ sin θ m min
cosθ k sin θ k ρ − x (cosθ k + sinθ k ) β = h min ∆ x sinθ k
Avec α = −
Pour θ k >
2 : 2
Rg ( ρ h , θ k ) =
Avec
∆x cosθ
n max
α′ =
∑ g ([α ′n + β ′], n)
n min
1
α
ρ − x (cosθ k + sinθ k ) β ′ = h min ∆ x cosθ k Les entiers mmin, mmax, nmin, nmax sont donnés par les formules III.3.11 et III.3.14 . Rappelons que dans notre cas : ∆x = 1 pour des raisons expliquée précédemment. (Voir considération numérique)
III.4. Algorithme de la TR : En entrée nous avons : g () , M , ∆ρ, ∆θ . En sortie on aura la matrice : g_radon. Subrotines utilisées dans l’algorithme : Radon_Donnée : subrotine qui permet de calculer H , K , xmin à partir Des formules III.2.5, III.2.8 , III.2.11. Radon_Remplire : permet de remplir les vecteurs θ() et ρ() en utilisant III.2.2. Radon_Limite : permet de calculer mmin, mm ax, nmin, nmaxsuivant Les formules III.3.11, III.3.14. Round : fonction qui permet d’approcher par l’entier le plus proche. Les codes donnant ces différentes subrotines sont donnés dans l’Annexe B.
15
L’algorithme de la TR se présentera comme suit : Début Calcul de K et H \\ Appel de Radon_Donnée Remplissage des vecteur θ() et ρ() \\ Appel de Radon_Remplire Pour k = 0 jusqu'à K −1 avec un pas 1 Costheta= cos θ k Sintheta= sin θk rho= xmin (cosθ k + sinθ k ) Si sin θ k > 2 / 2 alors Alpha= − cosθ k / sin θ k Bêta= ( ρh − rho)/sintheta Calcul de mmin et mmax \\ appel de Radon_Limite Somme = 0 Pour m = mmin jusqu'à mmax avec un pas 1 Somme=somme+ g (m, round (alpha* m +bêta )) Fin g_radon (h, k)=summe/abs (sintheta) Fin Sinon Alpha= − sin θ k / cosθ k Pour h = 0 jusqu'à H −1 avec un pas 1 Bêta= ( ρh − rho)/costheta
16
Calcul de nmin et nmax \\ appel de Radon_Limite Somme= 0 Pour n = nmin jusqu'à nmax avec un pas 1 Somme=somme+ g ( round (alpha* n +bêta ), n) Fin g_radon (h, k)=somme/abs (costheta) Fin Fin Si Fin
III.5. Utilisation de l’approximation linéaire: Par souci de précision, on se propose dans cette partie d’introduire une interpolation linéaire dans le calcul précédent. Le changement n’affectera que les formules III.3.8 et III.3.12. Définissons d’abord le paramètre ω : ∗ ω = n∗ − n , avec n =αm + β
ω
et n = [αm + β ]
représente l’écart entre la valeur de n∗et l’entier le plus proche
n.
D’où l’interpolation appliquée a III.3.8 : Rg ( ρh ,θk ) =
∆x M −1 ∑((1 −ω) g (m, n) +ωg (m, n +1)) sin θ m =0
(III.5.1)
De même posons pour III.3.12 : ∗ ω = m∗ − m , avec m =α′n + β′ et
m = [α ′n + β ′]
D’où l’interpolation :
Rg ( ρh ,θk ) =
∆x M −1 ∑((1 −ω) g (m, n) +ω(m +1, n)) cos θ n =0
(III.5.2)
L’algorithme est le même sauf qu’on utilisera les formules III.5.1 et III.5.2.
17
Chapitre IV : Résolution et discussion : IV.1. Résultat Obtenu pour différentes images : Voici représentées, figure 5, la transformée de Radon pour deux images différentes avec un nombre de projections égal à 175 .
18
Figure 5 : TR de deus images distinctes avec 175 projections On remarque que la signature de chaque image diffère de l’autre, ceci résulte de l’unicité de la transformée de RADON. Dans les deux cas on voit qu’il est possible de négliger la partie supérieure et inférieur de la TR obtenue car elle ne contient aucune information sur l’image originale.
IV.2. Influence du nombre de projections : Considérons l’image donnée par la figure 6 :
Figure 6 : Image ISAR d’un appareil
19
Appliquons la transformée de RADON avec différents nombre de projections, 31,60,121,et 175 projections, ce qui revient a faire varier le pas ∆θ , nous obtenons pour l’image précédente les signatures ci-dessous. Nous remarquons que plus le nombre de projections augmente et plus la signature devient plus claire et les discontinuités disparaissent.
60 projections
121 projections
20
175 projections
31 projections
IV.3.Influence de la rotation : Dans cette partie nous allons visualiser l’effet de la rotation de l’image originale sur la TR, pour cela, nous allons choisir un nombre donné de projections, soit 121, et effectuer une rotation donnée soit 45 degrés. Nous obtenons le résultat de la figure 7 :
21
22
Figure 7 : TR après une rotation de 45 degrés de l’image originale On voit clairement qu’une rotation dans le domaine spatial se traduit par une translation dans le domaine de RADON suivant l’axe θ. Donc la TR est bel et bien invariante pour la rotation.
IV.4.Comparaison avec les résultats obtenus sur Matlab : Une comparaison qualitative montre que les deux signature données par le VB et MATLAB sont presque identique.(voir figure 8 et 9) Néanmoins une comparaison quantitative et plus poussée s’impose. Pour cela nous allons procéder a la corrélation pour pouvoir confondre ou non les deux signatures.
TR donnée par VB TR donnée par MATLAB
Figure 8 : TR obtenues pour 30 projections
23
TR donnée par VB TR donnée par MATLAB Figure 9 : TR obtenues pour 120 projections
Résultats de la corrélation : Pour différent nombre de projections, nous avons calculé le coefficient de corrélation pour les deux méthodes utilisées, les résultats sont donnés dans le tableau figure 10. D’après le tableau, on remarque que les deux méthodes ont approximativement la meme précision avec une légère amélioration pour la deuxième méthode grâce a l’interpolation linéaire. 31 121 175 projections projections projections Sans interpolation linéaire
0.7028
0.7106
0.7144
Avec interpolation linéaire
0.7089
0.7110
0.7161
Figure 10 : Coefficient de corrélation des deux méthodes relativement a MATLAB La différence entre les résultats obtenus avec le VB et ceux de MATLAB est du au fait que MATLAB utilise des méthodes d’interpolation et d’approximation différentes, et du fait qu’il est mieux adapté pour le calcul matriciel.
IV.5-Conclusion : D’après les résultats obtenus précédemment, on a pu vérifier l’unicité de la signature ainsi que l’invariance par rapport a la rotation. Le nombre de projection effectué détermine la précision de la signature obtenue.
24
L’utilisation des méthodes d’approximation tel l’interpolation linéaire permettent de modifier relativement la signature finale, mais du point de vu qualitatif elle n’a pas de grande importance.
Conclusion générale : Il existe plusieurs méthodes de calcul de la Transformée de RADON, utilisant différentes approches. Dans ce travail, nous avons exposé une méthode permettant d’arriver a un résultat satisfaisant surtout du point de vue propriétés. L’importance de cette transformée réside dans le fait qu’elle permet de réduire sensiblement la taille de notre base de donnée, et d’avoir un traitement plus accéléré de l’information. Notons que cette transformée peut être inversée, et ainsi remonter d’une signature obtenue dans le domaine de RADON a l’image correspondante dans le domaine spatiale . Ceci constitue meme le principe du scanner médical. Notons qu’il existe une autre transformée déduite directement de celle de RADON appelée transformée de HOUGH, cette dernière est mieux adaptée dans le cas d’images « sparse » c’est-à-dire dont la matrice est pleine de zéros.
25
BIBLIOGRAPHIE
[1]
N. MEZHOUD, ‘‘Identification et classification de cibles radar non coopératives au moyen d’images de type ISAR’’, EMP, 30/ 12/ 2003.
[2]
Ali .Mohamed-Djafar, ‘‘Transformée de Radon et reconstruction d’image’’, Ecole supérieur d’électricité, 05/ 02/2002.
[3]
Luc Lemaire, ‘‘La recherche mathématique aujourd’hui’’, Université Libre de Bruxelles, 2000.
[4] ‘‘Formation traitement d’image avec interactive Data language (IDL)’’, copyright RSI, France, 2005 [6]
A. Averbuch, R.R. Coifman, D. L. Donoho, M. Israeli, J.Waldén, ‘‘a notion of Radon Transform for Data in a Cartesian Grid
[7] A.J.Higginson, “Analyzing Ship Wakes in Synthetic Aperture Radar Images (developpement of the RedBox Application)”, Defence research and developpement Canada-Ottawa, December 2002. [8]
Peter Toft. “Curve Parameter Detection using the Radon Transform. Part I”, 1996 26
[9]
Peter Toft. “Appendices, Part III ”,1996
[10] Stéphane Grognet, ‘‘ Calcul Symbolique et rigidité en courbure négative ’’,06/12 /2002 [11] Crlos A. Berenstein, Enrico casadio Tarabusi, and Arpad Kursa, ‘‘Radon Transform on spaces’’, Procedings of the American Mathematical Society, 1997.
ANNEXE A Impulsion de Dirac :
L’impulsion de Dirac est définie par : +∞
si x=0
δ(x) = 0
ailleurs
Autre définition de la transformée de Radon :
27
Transformée de Hough : La transformée de HOUGH d’une fonction g est définie par : g ( x, y ) =
∞ ∞
∫ ∫g(x
∗, y ∗)δ ( x − x∗)δ ( y − y ∗) dx ∗dy ∗
−∞−∞
ANNEXE B Programmes en VB : Function round(ByVal X As Double) As Integer Dim N As Integer '-------------------------------------------------------------------------------------------------'cette fonction permet de donner l'entier le plus proche du reel positif x '-------------------------------------------------------------------------------------------------N = Int(X) If X - N > 0.5 Then round = N + 1 Else round = N End Function ___________________________________________________________________________ Sub roundmat(ByRef X() As Double, _ ByVal N As Integer, _ ByVal m As Integer) Dim i As Integer, J As Integer '---------------------------------------------------------------------------------------------------'Cette subrotine permet d'arrondir chaque élément ' d 'une matrice a l 'entier le plus proche '---------------------------------------------------------------------------------------------------For i = 0 To N - 1 For J = 0 To m - 1 X(i, J) = round(X(i, J))
28
Next J Next i End Sub ___________________________________________________________________________ Sub Radon_Remplire(ByVal N As Integer, _ ByVal dx As Double, _ ByVal xmin As Double, _ ByRef X() As Double) Dim i As Integer '-----------------------------------------------------------------------------------' Cette subrotine permet de remplir un vecteur avec un pas constant dx '-----------------------------------------------------------------------------------ReDim X(0 To N - 1) For i = 0 To N - 1 X(i) = xmin + i * dx Next i End Sub ___________________________________________________________________________ ___________________________________________________________________________ Sub Radon_Donnée(ByVal m As Integer, _ ByVal dphi As Double, _ ByVal drho As Double, _ ByRef k As Integer, _ ByRef h As Integer, _ ByRef xmin As Double, _ ByRef rhomin As Double) Const pi As Double = 3.14159 ‘------------------------------------------------------------------------------------------------------‘cette subrotine permet de calculer des parametres utilisés dans la fonction principale ‘------------------------------------------------------------------------------------------------------k = round(pi / dphi) h = round((Sqr(2) * m) / drho + 1) xmin = -(m - 1) / 2 rhomin = (-(h - 1) / 2) * drho End Sub ___________________________________________________________________________ Sub Radon_limite(ByVal m As Integer, _ ByVal alpha As Double, _ ByVal beta As Double, _ ByRef mmin As Integer, _ ByRef mmax As Integer) ‘-------------------------------------------------------------------------------------------------‘cette subrotine calcule les limites de la boucle utilisée dans la fonction principale
29
‘-------------------------------------------------------------------------------------------------If alpha > 0 Then mmin = max(0, round(-beta / alpha)) mmax = min(m - 1, round((m - 1 - beta) / alpha)) ElseIf alpha < 0 Then mmin = max(0, round((m - 1 - beta) / alpha)) mmax = min(m - 1, round(-beta / alpha)) Else mmin = 0 mmax = m - 1 End If End Sub ___________________________________________________________________________
Fonctions donnant la TR: ___________________________________________________________________________ Function Radon_transform(ByRef g() As Double, _ ByVal m As Integer, _ ByVal dphi As Double, _ ByVal drho As Double, _ ByRef Rg() As Double, _ Optional ByRef k As Integer, _ Optional ByRef h As Integer) As Integer '-----------------------------------------------------------------------------------------------'Declaration de variables '-----------------------------------------------------------------------------------------------Dim rho() As Double, phi() As Double, rhof As Double, N As Double Dim xmin As Double, rhomin As Double, somme As Double Dim alpha As Double, beta As Double, mmin As Integer, mmax As Integer Dim i As Integer, J As Integer, t As Integer, R As Integer '----------------------------------------------------------------------------------------------' Radon_transform permet de calculer la transformée de Radon d'une fonction ' discrete donnée par la matrice g() et nous rend en sortie une matrice Rg ' Parametres d'entrée: ' M :dimension de la matrice g() ' dphi:pas de l'angle polaire ' drho:pas du rayon polaire '-------------------------------------------------------------------------------------------------
30
Call Radon_Donnée(m, dphi, drho, k, h, xmin, rhomin) 'Calcul de K,H,xmin et rhomin Call Radon_Remplire(k, dphi, 0, phi()) 'remplissage du vecteur phi() Call Radon_Remplire(h, drho, rhomin, rho()) 'remplissage du vecteur rho() ReDim Rg(h, k) For i = 0 To k - 1 rhof = xmin * (Cos(phi(i)) + Sin(phi(i))) If Sin(phi(i)) > 1 / Sqr(2) Then alpha = -Cos(phi(i)) / Sin(phi(i)) For J = 0 To h - 1 beta = (rho(J) - rhof) / Sin(phi(i)) Call Radon_limite(m, alpha, beta, mmin, mmax) somme = 0 For t = mmin To mmax N = alpha * t + beta If N >= 0 And N < m Then somme = somme + g(t, round(N)) End If Next t Rg(J, i) = somme / Abs(Sin(phi(i))) Next J Else alpha = -Sin(phi(i)) / Cos(phi(i)) For J = 0 To h - 1 beta = (rho(J) - rhof) / Cos(phi(i)) Call Radon_limite(m, alpha, beta, mmin, mmax) somme = 0 For t = mmin To mmax N = alpha * t + beta If N >= 0 And N < m Then somme = somme + g(round(N), t) End If Next t Rg(J, i) = somme / Abs(Cos(phi(i))) Next J End If Next i End Function ___________________________________________________________________________ Function Radon_Lineartrans(ByRef g() As Double, _ ByVal m As Integer, _ ByVal dphi As Double, _ ByVal drho As Double, _ ByRef Rg() As Double, _ Optional ByRef k As Integer, _ Optional ByRef h As Integer) As Integer '-----------------------------------------------------------------------------------------------------'Declaration de variables
31
'-----------------------------------------------------------------------------------------------------Dim rho() As Double, phi() As Double, rhof As Double, N As Double Dim xmin As Double, rhomin As Double, somme As Double Dim alpha As Double, beta As Double, mmin As Integer, mmax As Integer Dim i As Integer, J As Integer, t As Integer, R As Integer, w As Double '-----------------------------------------------------------------------------------------------------' Radon_Lineartrans utilise le meme principe que Radon_transform ' mais avec une interpolation linéaire. '------------------------------------------------------------------------------------------------------Call Radon_Donnée(m, dphi, drho, k, h, xmin, rhomin) 'Calcul de K,H,xmin et rhomin Call Radon_Remplire(k, dphi, 0, phi()) 'remplissage du vecteur phi() Call Radon_Remplire(h, drho, rhomin, rho()) 'remplissage du vecteur rho() ReDim Rg(h, k) For i = 0 To k - 1 rhof = xmin * (Cos(phi(i)) + Sin(phi(i))) If Sin(phi(i)) > 1 / Sqr(2) Then alpha = -Cos(phi(i)) / Sin(phi(i)) For J = 0 To h - 1 beta = (rho(J) - rhof) / Sin(phi(i)) Call Radon_limite(m, alpha, beta, mmin, mmax) somme = 0 For t = mmin To mmax N = alpha * t + beta w = N - round(N) If N >= 0 And N < m Then somme = somme + (1 - w) * g(t, round(N)) + w * g(t, round(N) + 1) End If Next t Rg(J, i) = somme / Abs(Sin(phi(i))) Next J Else alpha = -Sin(phi(i)) / Cos(phi(i)) For J = 0 To h - 1 beta = (rho(J) - rhof) / Cos(phi(i)) Call Radon_limite(m, alpha, beta, mmin, mmax) somme = 0 For t = mmin To mmax N = alpha * t + beta w = N - round(N) If N >= 0 And N < m Then somme = somme + (1 - w) * g(round(N), t) + w * g(round(N) + 1, t) End If Next t Rg(J, i) = somme / Abs(Cos(phi(i))) Next J
32
End If Next i End Function
Remarques de l’enseignant
33
