Procédure de calibration astrométrique

Alain KLOTZ
Mise à jour 20020301

Cette méthode est implantée dans la librairie TT. On expose ici les liens entre des coordonnées RA,DEC et les coordonnées x,y sur l'image dans le concept des mots clé FITS. Les sept premiers paragraphes concernent la calibration dite "linéaire". Les paragraphes suivants concernent la calibration d'un champ ayant des distorsions (l'exposé de cette dernière partie peut comporter des erreurs).

1. Paramètres instrumentaux :

(RA,DEC) les coordonnées célestes d'une étoile du champ
(x,y) ses coordonnées cartésiennes dans l'image.
(RAo,DECo) les coordonnées célestes approximatives du centre de l'image.
(F) la longueur ficale du télescope
(px, py) la dimension physique des pixels de la caméra CCD
(Nx, Ny) le nombre de pixels de la caméra CCD

2. La norme FITS définit :

CRPIX1 : est le pixel du centre de la projection sur NAXIS1 (axe x)
CRPIX2 : est le pixel du centre de la projection sur NAXIS2 (axe y)
CRVAL1 : est la coordonnée céleste exacte sur RA correspondant à CRPIX1. (doit être proche de RAo)
CRVAL2 : est la coordonnée céleste exacte sur DEC correspondant à CRPIX2. (doit être proche de DECo)
CDELT1 : est l'échantillonnage angulaire sur NAXIS1 (axe x)
CDELT2 : est l'échantillonnage angulaire sur NAXIS2 (axe y)
CROTA2 : est l'angle de position de l'axe des déclinaisons par rapport à NAXIS1 (axe x)

3. Le passage des paramètres du télescope aux paramètres FITS est :

(équations 3)
CDELT1=-atan (px/F)
CDELT2=atan (py/F)
CRVAL1=RAo
CRVAL2=DECo
CRPIX1= Nx/2
CRPIX2= Ny/2

A noter que dans le cas d'une image présentée nord en haut et ouest à droite, CDELT1 doit être négatif et CROTA2=0.

4. Matrice de passage CD :

(équations 4)
CD001001 = CDELT1*cos(CROTA2)
CD001002 = fabs(CDELT2)*CDELT1/fabs(CDELT1)*sin(CROTA2)
CD002001 = -fabs(CDELT1)*CDELT2/fabs(CDELT2)*sin(CROTA2)
CD002002 = CDELT2*cos(CROTA2)

5. Passage (RA,DEC) -> (x,y)

(équations 5)
H = sin(DEC)*sin(CRVAL2) + cos(DEC)*cos(CRVAL2)*cos(RA-CRVAL1)
dRA = cos(DEC)*sin(RA-CRVAL1) / H
dDEC = [ sin(DEC)*cos(CRVAL2) - cos(DEC)*sin(CRVAL2)*cos(RA-CRVAL1) ] / H
det=CD002002*CD001001-CD001002*CD002001
x = CRPIX1 - (CD001002*dDEC - CD002002*dRA) / det
y = CRPIX2 + (CD001001*dDEC - CD002001*dRA) / det

6. Passage (x,y) -> (RA,DEC)

(équations 6)
dRA = CD001001 * (x-CRPIX1) + CD001002 * (y-CRPIX2)
dDEC = CD002001 * (x-CRPIX1) + CD002002 * (y-CRPIX2)
delta= cos(CRVAL2) - dDEC*sin(CRVAL2)
gamma= sqrt( dRA*dRA + delta*delta )
RA = CRVAL1 + atan (dRA/delta)
DEC = atan ( [sin(CRVAL2)+dDEC*cos(CRVAL2)] / gamma )

7. Calibration d'un champ sans distorsion

On dispose de RAo,DECo (coordonnées des codeurs) , F, px, py, Nx, Ny et une liste d'étoiles (x,y,flux) mesurées sur l'image à calibrer.

On calcule (équations 3) et (équations 4)
On génère une liste (x',y',mag) d'étoiles du catalogue de référence (équations 5) correspondant au champ observé théorique.

On effectue l'appariement entre les listes (x,y) et (x',y') et on optimise le jeu de paramètres (a0, a1, a2, a3, a4, a5) :

x' = a0*x + a1*y + a2
y' = a3*x + a4*y + a5

On calcule les nouveaux paramètres de projection (suivi de ')

xc'= a0*CRPIX1 + a1*CRPIX2 + a2
yc'= a3*CRPIX1 + a4*CRPIX2 + a5

Equations (6) pour passer de (xc',yc') à (CRVAL1',CRVAL2')
CRPIX1'=CRPIX1
CRPIX2'=CRPIX2

CD001001' = CD001001*a0 + CD001002*a3
CD001002' = CD001001*a1 + CD001002*a4
CD002001' = CD002001*a0 + CD002002*a3
CD002002' = CD002001*a1 + CD002002*a4

Pour calculer CDROTA2', il faut tenir compte du fait qu'il peut y avoir des imprécisions sur les valeurs de la matrice CD. Il est donc conseillé de procéder comme suit. On commence à calculer l'angle crota2 de deux manières différentes (angles aa et bb) et de les moyenner.

aa=fmod(2*pi+atan2(CD002001',CD001001'),2*pi)
bb=fmod(2*pi+atan2(-CD001002',CD002002'),2*pi)

Avant de moyenner les angles aa et bb, il faut s'assurer qu'ils ne sont pas séparés de 180 degrés (cas qui arrive lorsque les angles aa et bb sont proches des axes x ou y). Il faut donc calculer l'angle de sépartion de aa et bb (angle ab) et ajouter 180 suivant le cas:

cosa=cos(aa)
sina=sin(aa)
cosb=cos(bb)
sinb=sin(bb)
cosab=cosa*cosb+sina*sinb
sinab=sina*cosb-cosa*sinb
ab=fabs(atan2(sinab,cosab))
if (ab>pisur2) {
    if (cosa>cosb) {
      bb=fmod(pi+bb,2*pi)
   } else {
      aa=fmod(pi+aa,2*pi)
   }
}

On peut maintenant effectuer la moyenne :
ab=bb-aa;
if (ab>pi) {
    aa=aa+2*pi
}
ab=aa-bb
if (ab>pi) {
    bb=bb+2*pi
}
CROTA2'=fmod(2*pi+(aa+bb)/2.,2*pi)

On peut alors calculer les échelles CDELT1' et CDELT2' :

cosr=fabs(cos(CROTA2'))
sinr=fabs(sin(CROTA2'))
if (cosr>sinr) {
    CDELT1'=CD001001'/cos(CROTA2')
    CDELT2'=CD002002'/cos(CROTA2')
} else {
    CDELT1'=fabs(-CD002001'/sin(CROTA2'))
    CDELT2'=fabs(CD001002'/sin(CROTA2'))
    signr=sinr/fabs(sinr)
    signc=CD001002'/fabs(CD001002')
    signd=signc/signr
    if (signd<0) { CDELT1'=-CDELT1' }
    signc=CD002001'/fabs(CD002001')
    signd=-signc/signr
    if (signd<0) { CDELT2'=-CDELT2' }
}

Dans l'entête FITS, on pensera aussi à ajouter les lignes suivantes pour assurer la compatibilité entre logiciels :

RADESYS='FK5     ' / Mean Place IAU 1984 system
EQUINOX=2000.0 / System of equatorial coordinates
CTYPE1='RA---TAN' /Gnomonic projection
CTYPE2='DEC--TAN' /Gnomonic projection
LONPOLE=180.0 / Long. of the celest.NP in native coor.syst.
CUNIT1='deg     ' / Angles are degrees always
CUNIT2='deg     ' / Angles are degrees always

7.1. Note pour les cas de transformations géométriques d'images déjà calibrées

Dans le cas où l'on souhaite mettre à jour les mots clé d'une image calibrée qui subit une tranformation géométrique de type "linéaire", il faut recalculer la valeur des mots clés. Si (x,y) est un point de l'image transformée et (x',y') le point correspondant sur l'image de départ, on utilise la relation suivante:

x' = a0*x + a1*y + a2
y' = a3*x + a4*y + a5

L'équation est similaire à celle concernant l'appariement pour la calibration astrométrique. Dans le cas d'une transformation géométrique, on doit considérer que le point de référence (CRVAL1,CRVAL2) est toujours le même et que ce sont les coordonnées (CRPIX1,CRPIX2) qu'il faut modifier. On utilisera donc le jeu de relations suivantes :

Calculer les coefficients (aa0,aa1,aa2,aa3,aa4,aa5) correspondant à la relation inverse :

x = aa0*x' + aa1*y' + aa2
y = aa3*x' + aa4*y' + aa5

delta=a1*a3-a0*e4
aa0=-a4/delta;
aa1= a1/delta;
aa2=-(a1*a5-a2*a4)/delta;
aa3= a3/delta;
aa4=-a0/delta;
aa5=-(a2*a3-a0*a5)/delta;

CRPIX1'=aa0*CRPIX1+aa1*CRPIX2+aa2
CRPIX2'=aa3*CRPIX1+aa4*CRPIX2+aa5

CRVAL1'=CRVAL1
CRVAL2'=CRVAL2

La suite des relations est rigoureusement la même que pour la calibration astrométrique:

CD001001' = CD001001*a0 + CD001002*a3
etc.

8. Calibration d'un champ avec distorsions

8.1. Exposé théorique

Soit (x,y) la liste des coordonnées des étoiles sur l'image
Soit (x',y',mag) la liste des étoiles du catalogue de référence correspondant au champ observé théorique, projetée selon une transformation gnomonique sans distorsion (équations 5).

On effectue l'appariement entre les listes (x,y) et (x',y') et on optimise le jeu de paramètres (b0 à b11) :

(équations 7)
x' = b0*x + b1*y + b2 + b3*x*y + b4*x*x + b5*y*y
y' = b6*x + b7*y + b8 + b9*x*y + b10*x*x + b11*y*y

Malheureusement, cette écriture ne permet pas de calculer des paramètres physiques tels que ceux décrits par les (équations 3) et les (équations 4). C'est pour cette raison que les coefficients (b0 à b6) ne sont pas intégré dans la norme FITS. Pour y arriver, il faut construire un système de coordonnées intermédiaires (x",y") correspondant aux positions (x,y) de façon à avoir :

(équations 8)
x' = a0*x" + a1*y" + a2
y' = a3*x" + a4*y" + a5

ainsi on conserve la notion de coefficients "linéaires" classiques (équations 3 et 4) auxquels ont ajoute les coefficients d'une transformation supplémentaire permettant de passer du système (x,y) à (x",y"). La transformation, à trouver, de (x,y) vers (x",y") sera du type :

(équations 9)
x" = c0*x + c1*y + c2 + c3*x*y + c4*x*x + c5*y*y
y" = c6*x + c7*y + c8 + c9*x*y + c10*x*x + c11*y*y

Les coefficients (c0 à c11) sont définits dans la norme FITS. Cherchons à déterminer leur valeur. En remplaçant les (équations 9) dans les (équations 8) on obtient :

x' = (a0*c0+ a1*c6)*x + (a0*c1+ a1*c7)*y + (a0*c2+ a1*c8+ a2) + (a0*c3+ a1*c9)*x*y + (a0*c4+ a1*c10)*x*x + (a0*c5+ a1*c11)*y*y
y' = (a3*c0+ a4*c6)*x + (a3*c1+ a4*c7)*y + (a3*c2+ a4*c8+ a5) + (a3*c3+ a4*c9)*x*y + (a3*c4+ a4*c10)*x*x + (a3*c5+ a4*c11)*y*y

Il reste à identifier les relations entre les coefficients b (connus) et les coefficients a et c (inconnus) :

(équations 10)
b0 = a0*c0+ a1*c6
b1 = a0*c1+ a1*c7
b2 = a0*c2+ a1*c8+ a2
b3 = a0*c3+ a1*c9
b4 = a0*c4+ a1*c10
b5 = a0*c5+ a1*c11
b6 = a3*c0+ a4*c6
b7 = a3*c1+ a4*c7
b8 = a3*c2+ a4*c8+ a5
b9 = a3*c3+ a4*c9
b10 = a3*c4+ a4*c10
b11 = a3*c5+ a4*c11

Nous sommes face à un système qui comporte 6 inconnues de trop par rapport au nombre d'équations. La valeur des coefficients a et c n'est donc pas univoque.

Une solution de bon sens consiste à calculer la valeur des coefficients de passage (a0 à a5) pour le centre du champ où les déformations ne sont pas trop importantes. Pour cela, on effecute une calibration sans distorsion sur la partie centrale du champ uniquement. Ainsi, les coefficients (a0 à a5) sont connus et les valeurs des coefficients (c0 à c11) prennent pour valeur :

(équations 11)
c0=(a1*b6-a4*b0)/(a1*a3-a0*a4)
c1=(a1*b7-a4*b1)/(a1*a3-a0*a4)
c2=(a1*b8-a4*(b2-a2))/(a1*a3-a0*a4)
c3=(a1*b9-a4*b3)/(a1*a3-a0*a4)
c4=(a1*b10-a4*b4)/(a1*a3-a0*a4)
c5=(a1*b11-a4*b5)/(a1*a3-a0*a4)
c6=(a3*b0-a0*b6)/(a1*a3-a0*a4)
c7=(a3*b1-a0*b7)/(a1*a3-a0*a4)
c8=(a3*b2-a0*(b8-a5))/(a1*a3-a0*a4)
c9=(a3*b3-a0*b9)/(a1*a3-a0*a4)
c10=(a3*b4-a0*b10)/(a1*a3-a0*a4)
c11=(a3*b5-a0*b11)/(a1*a3-a0*a4)

N.B. : La norme FITS permet de calculer les coefficients jusqu'au septième ordre. Nous avons volontairement limité notre étude au second ordre pour garder une écriture simple.

8.1. Exposé pratrique

Soit (x,y) la liste des coordonnées des étoiles sur l'image entière
Soit (x',y',mag) la liste des étoiles du catalogue de référence correspondant au champ observé théorique de l'image entière, projetée selon une transformation gnomonique sans distorsion (équations 5).

Extraire deux sous listes (X,Y) et (X',Y',MAG) correspondant à la partie centrale de l'image à calibrer.

Calculer les coefficients (a0 à a5) ainsi que les mots clés FITS en effectuant une calibration sans distorsion (cf. paragraphe 7).

On effectue l'appariement entre les listes (x,y) et (x',y') et on optimise le jeu de paramètres (b0 à b11):

(équations 7)
x' = b0*x + b1*y + b2 + b3*x*y + b4*x*x + b5*y*y
y' = b6*x + b7*y + b8 + b9*x*y + b10*x*x + b11*y*y

Calculer les valeurs numériques des coefficients (c0 à c11):

(équations 11)
c0=(a1*b6-a4*b0)/(a1*a3-a0*a4)
c1=(a1*b7-a4*b1)/(a1*a3-a0*a4)
c2=(a1*b8-a4*(b2-a2))/(a1*a3-a0*a4)
c3=(a1*b9-a4*b3)/(a1*a3-a0*a4)
c4=(a1*b10-a4*b4)/(a1*a3-a0*a4)
c5=(a1*b11-a4*b5)/(a1*a3-a0*a4)
c6=(a3*b0-a0*b6)/(a1*a3-a0*a4)
c7=(a3*b1-a0*b7)/(a1*a3-a0*a4)
c8=(a3*b2-a0*(b8-a5))/(a1*a3-a0*a4)
c9=(a3*b3-a0*b9)/(a1*a3-a0*a4)
c10=(a3*b4-a0*b10)/(a1*a3-a0*a4)
c11=(a3*b5-a0*b11)/(a1*a3-a0*a4)

Dans l'entête FITS, on ajoutera les mots clés suivants :

(équations 12)
PV1_0=c2
PV1_1=c0
PV1_2=c1
PV1_3=0.0
PV1_4=c4
PV1_5=c3
PV1_6=c5
PV2_0=c8
PV2_1=c7
PV2_2=c6
PV2_3=0.0
PV2_4=c11
PV2_5=c9
PV2_6=c10

Les coefficients PV*_3 correspondent à un coefficient en sqrt(x*x+y*y) que nous ne calculons pas ici.

9. Passage (x,y) -> (RA,DEC) avec distorsions

On effectue d'abord la transformation :

(équations 9)
x" = PV1_1*x + PV1_2*y + PV1_0 + PV1_5*x*y + PV1_4*x*x + PV1_6*y*y
y" = PV2_2*x + PV2_1*y + PV2_0 + PV2_5*x*y + PV2_6*x*x + PV2_4*y*y

puis on utilise (équations 6) en remplaçant (x,y) par (x",y") pour calculer (RA,DEC).

10. Passage (RA,DEC) -> (x,y) avec distorsions

Avant d'effectuer le passage, il faut calculer les coefficients (d0 à d11) de la transformation inverse de l'(équation 9).

(équations 13)
x = d0*x" + d1*y" + d2 + d3*x"*y" + d4*x"*x" + d5*y"*y"
y = d6*x" + d7*y" + d8 + d9*x"*y" + d10*x"*x" + d11*y"*y"

On pourra, par exemple, effectuer un ajustement polynomial sur un nuage de points calculés à partir de l'(équation 9). Ceci peut être réalisé lors de la lecture des mots clés FITS d'astrométrie de l'image.

A partir de (RA,DEC), on utilise d'abord les (équations 5) pour calculer (x",y") puis on utilise les (équations 13) pour déterminer (x,y).