Détection des étoiles
Principe
D'abord, une image est une suite de valeurs dans un tableau. La valeur représente l'intensité de lumière dans une couleur.
Dans l'exemple suivant :
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | |
|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 8 | 10 | 8 | 7 | 8 | 10 | 11 | 7 | 9 | 10 |
| 2 | 9 | 7 | 8 | 11 | 9 | 10 | 8 | 11 | 9 | 9 |
| 3 | 7 | 11 | 9 | 8 | 12 | 6 | 5 | 12 | 8 | 8 |
| 4 | 8 | 12 | 10 | 7 | 20 | 70 | 80 | 8 | 10 | 7 |
| 5 | 9 | 8 | 9 | 12 | 50 | 208 | 60 | 1 | 9 | 12 |
| 6 | 10 | 7 | 8 | 7 | 30 | 80 | 20 | 10 | 11 | 10 |
| 7 | 9 | 10 | 9 | 8 | 10 | 30 | 7 | 8 | 9 | 10 |
| 8 | 11 | 6 | 10 | 7 | 7 | 8 | 9 | 7 | 8 | 12 |
| 9 | 9 | 8 | 12 | 9 | 11 | 10 | 9 | 8 | 9 | 6 |
| 10 | 8 | 8 | 10 | 8 | 9 | 8 | 7 | 6 | 8 | 8 |
La moyenne des valeurs de points est de 14.
Les étoiles sont des petites courbes de Gauss (courbe en cloche) qui forme des pics qui se distingue du “noir” ambiant. Le noir est représenté par la moyenne des valeurs de l'ensemble de l'image dans une couleur donnée. On peut choisir la couleur verte car il y a deux fois plus de pixel vert que de rouge et de bleu (Voir la matrice de Bayer).
Donc si on fixe un seuil qui est deux ou trois fois la valeur moyenne des intensités lumineuses, nous pouvons considérer que ce qui est au dessus sont les étoiles.
L'objectif de l'algorithme est d'isolé des rectangles de pixels qui contiennent des valeurs supérieures au seuil donnée.
Dans notre exemple, la grille devient en éliminant les points inférieurs à 28 (2 x 14) :
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | |
|---|---|---|---|---|---|---|---|---|---|---|
| 1 | ||||||||||
| 2 | ||||||||||
| 3 | ||||||||||
| 4 | 70 | 80 | ||||||||
| 5 | 50 | 208 | 60 | |||||||
| 6 | 30 | 80 | ||||||||
| 7 | 30 | |||||||||
| 8 | ||||||||||
| 9 | ||||||||||
| 10 |
Dans notre grille, une étoile apparait. Si on calcul le barycentre des points, on peut retrouver son centre.
Soit:
Le centre de l'étoile est donc aux coordonnées ( 6.1 , 5.03 ). On peut donc noter que même si les calculs sont fait à partir de valeurs entières, le résultat est en valeur réelle. On gardera cette précision pour l'alignement des images.
Algorithme
L'algorithme utilisé est le suivant :
Deux boucles imbriquées permettent de parcourir l'ensemble des lignes de l'images et pour chaque ligne, l'ensemble des points de l'image.
Les valeurs inférieurs au seuil défini sont ignorés.
Lecture d'une ligne de l'image
On défini une variable d'état qui bascule à chaque fois que l'on passe d'une valeur inférieur au seuil à une valeur supérieur. Cet état donne donc une information qui est le cœur du processus, à savoir si on est sur du noir ou si on est sur une étoile.
| Analyse de l'état | Conséquence |
|---|---|
| L'état passe du noir vers une étoile | On initialise une nouvelle étoile |
| tant que l'état est sur étoile | On collecte les informations de l'étoile |
| On passe de étoile vers noir | On termine l'étoile et on la range en mémoire |
L'algorithme range les étoiles dans une liste chaînée. Cette liste est rangée dans un tableau de longueur le nombre de lignes de l'image.
Donc à l'issu de la boucle, on aurait la liste des étoiles de chaque ligne. Si on reprend notre exemple, on obtient :
- Pour la ligne 4, on a une étoile qui commence en colonne 6 et qui se termine en 7.
- Pour la ligne 5, on a une étoile qui commence en colonne 5 et qui se termine en 7.
- Pour la ligne 6, on a une étoile qui commence en colonne 5 et qui se termine en 6.
- Pour la ligne 7, on a une étoile qui commence en colonne 6 et qui se termine en 6.
On peut remarquer que ces 4 étoiles sont en faite la même.
Structure étoile
Définissons une structure pour ranger les étoiles. Cette structure doit avoir les informations suivantes :
- Les coordonnées du pixel de départ sur la photo.
- Les coordonnées du pixel de fin sur la photo.
- Un pointeur vers l'étoile suivante puisque cette liste est une liste chainée
L'objectif est calculer à la fin du processus, le centre de l'étoile. Pour cela, il nous faut la liste des points composants l'étoile. Nous pourrons alors calculer le centre de l'étoile, son intensité (la somme de toutes les valeurs d'intensité des points composants cette étoile), l'intensité maximum atteint par un des points de l'étoile, une valeur FWHM et un indice. Ces deux dernières interviendront plus loin.
struct etoile { struct point debut; struct point fin; uint32 intensite; uint16 vMax; float centre_X; float centre_Y; float FWHM; int indice; struct point * points; struct etoile * suivant; };
La structure d'un point a pour objectif de ranger en mémoire, les coordonnées (x,y) de l'image, l'intesité du point (la valeur du pixel) et un pointeur vers le point suivant car cette structure est utilisée comme une liste chainée.
struct point { int x; int y; uint32 intensite; struct point * suivant; };
Terminaison d'une étoile
Lorsque l'état passe de “étoile” vers “noir”, l'algorithme “range l'étoile (donc la structure) dans le tableau des étoiles à l'indice “n° de ligne” mais avant, il regarde sur la ligne précédente si il n'y a pas une étoile qui touche celle en cours d'analyse. Si la réponse est oui, alors l'étoile de la ligne précédente est supprimée et ses points sont ajoutés à l'étoile courante.
Par exemple, lors du parcours de la ligne 5 de notre cas, le programme regarde la ligne 4. L'étoile de la ligne 4 touche celle de la ligne 5. Cette détection se fait en partant du principe que les segments définis par les rangs des colonnes doivent avoir au moins un point qui se recouvre. Dans notre cas, les segments sont [6,7] et [5,7], il y a bien un morceau en commun.
Lors du parcours de la ligne 4, le programme avait stocké une étoile qui débute en (6,4) et qui se termine en (7,4), avec la liste des points (6,4)→70 et (7,4)→80. Lors du parcours de la ligne 5, le programme a pour l'étoile de la ligne 5 un début en (5,5) et une fin en (7,5), avec la liste de points (5,5)→50, (6,5)→208 et (7,5)→60.
Le programme va:
- supprimer l'étoile de la ligne 4
- ajouter les points (6,4)→70 et (7,4)→80 à la liste (5,5)→50, (6,5)→208 et (7,5)→60
- modifier le point de début de l'étoile à (5,4) et
- modifier le point de fin de l'étoile à (7,5)
A la fin du processus, nous aurons qu'une étoile avec les coordonnées (5,4) (7,7) et la liste de points suivantes :
(6,4)→70 (7,4)→80 (5,5)→50, (6,5)→208 (7,5)→60 (5,6)→30 (6,6)→80 (6,7)→30
Réalisation en C
Le code suivant est l'application de ce qui est décrit ci-dessus. Il faut noter que la manipulation de liste chainée permet de faire un code rapide et concis mais la lecture n'est pas toujours simple.
La variable etat donne l'information de pilotage de l'algorithme et a pour valeur 0 pour “étoile” et 1 pour “noir”.
On note en début, l'initialisation de la structure étoile.
La variable etoiles est un tableau donnant le pointeur vers la liste chainée des étoiles détectées sur une ligne.
// On traite la ligne des verts int p = G[count]; if (p < vCoupure) p = 0 ; //vCoupure est la fréquence seuil de détection d'une étoile if ((p != 0) && (etat==1)) { // début d'une étoile etat = 0; struct point pp; e = malloc(sizeof(struct etoile)); pp.x=x; pp.y=y; e->debut=pp; e->fin.x=0; e->fin.y=0; e->intensite=0; e->points=NULL; e->suivant=NULL; e->indice = -1; e->vMax=0; //if (debug) printf("Debut x %i,y%i \n",e->debut.x,e->debut.y); } if ((p == 0) && (etat==0)) { // Fin d'une etoile struct point pp; pp.x=x; pp.y=y; e->fin=pp; //if (debug) printf(" fin x %i,y%i \n",e->fin.x,e->fin.y); // On cherche sur la ligne précédente une étoile if ((y > 0) && (etoiles[y-1]!=NULL)) { struct etoile *ee=etoiles[y-1]; struct etoile *pee=NULL; int ok=0; do { //struct etoile * ee =pe; if ((min(ee->fin.x,e->fin.x) - max(ee->debut.x,e->debut.x))>0) { // ok ok=1; e->debut.x= min(ee->debut.x,e->debut.x); e->fin.x= max(ee->fin.x,e->fin.x); e->debut.y=ee->debut.y; e->intensite+=ee->intensite; e->indice=ee->indice; e->vMax=max(e->vMax,ee->vMax); // Copie des points if (e->points==NULL) e->points=ee->points; else { struct point *pse=e->points; while (pse->suivant!=NULL) pse=pse->suivant; pse->suivant = ee->points; } // On supprime l'ancienne ligne if (pee==NULL) etoiles[y-1]=ee->suivant; else pee->suivant=ee->suivant; pee=ee; ee=ee->suivant; } while ((ee!=NULL)&&(ok==0)); } // Stockage des informations if (etoiles[y]!=NULL) { struct etoile *pe=etoiles[y]; while (pe->suivant!=NULL) pe=pe->suivant; pe->suivant=e; } else { etoiles[y]=e; } etat = 1; } if (etat==0) { // Etoile en cours e->intensite+=p; e->vMax=max(e->vMax,p); Point = malloc(sizeof(struct point)); Point->x=x; Point->y=y; Point->intensite=p; Point->suivant=NULL; if (e->points==NULL) e->points=Point; else { struct point *ppp=e->points; while (ppp->suivant!=NULL) ppp=ppp->suivant; ppp->suivant=Point; ppp=e->points; while (ppp!=NULL) ppp=ppp->suivant; } }
Calculs finaux
Une fois les étoiles détectées, il est possible de calculer l'intensité, le centre de l'étoile et une valeur FWHM.
L'intensité d'une étoile est donnée par la somme de tous les points supérieurs au seuil situé dans ce carré.
Le centre de l'étoile est donné par le barycentre (la moyenne pondérée) de ces points.
// Calcul barycentre float xM = 0; float yM = 0; float M = 0; struct point *ppp=pe->points; while (ppp!=NULL) { xM +=ppp->x*ppp->intensite; yM +=ppp->y*ppp->intensite; M +=ppp->intensite; //printf(" p (%i,y%i) %i \n",ppp->x,ppp->y,ppp->intensite); ; ppp=ppp->suivant; } pe->centre_X=xM/M; pe->centre_Y=yM/M;
Une valeur FWHM est calculée également. Le FWHM “Full width at half maximum” est la largeur de la courbe de gauss à mi-hauteur de son intensité maximum. Vous pouvez trouver plus de détail mathématique sur Wikipedia. Donc si on repère dans le carré de l'étoile le point maximum et que l'on divise par deux cette valeur maximale, on obtiendrait la mi-hauteur. La valeur FWHM serait donc la distance du point dont la valeur du pixel est au dessus de la mi-hauteur et qui est le plus éloigné du centre de l'étoile, multipliée par deux. Cette définition simpliste (qui est implémenté dans l'algorithme) ne marche pas pour les étoiles saturées mais donne quand même une idée de la qualité de la photo en faisant la moyenne de toutes les FWHM de toutes les étoiles trouvées.
Je précise bien que ma valeur FWHM est une approximation simpliste. Pour réaliser un calcul correct, il faudrait interpoler les points par une courbe de Gauss. Une fois ce travail effectué le calcul devient faisable.
A faire : comparer ma méthode avec Iris sur quelques étoiles.
// On prend le point superieur au FWMH le plus loin ... int iFWHM = pe->vMax/2; float d; float dMax=0; ppp=pe->points; while (ppp!=NULL) { if (ppp->intensite > iFWHM) { d = ((float)ppp->x - (float)pe->centre_X)* ((float)ppp->x - (float)pe->centre_X)+ ((float)ppp->y - (float)pe->centre_Y)* ((float)ppp->y - (float)pe->centre_Y); dMax=max(d,dMax); } ppp=ppp->suivant; } pe->FWHM=sqrt(d)*2;
