Cartographie pixel à RA / DEC en astrophotographie numérisée

J'ai une photo 1443x998 des étoiles (prise avec un appareil photo 35 mm puis scannée) avec les étoiles suivantes aux emplacements de pixels suivants:

Altair x=782, y=532 [19h50m46.9990s RA, +08 52'05.959'' DEC] 
Sualocin, x=311, y=146 [20h 39m 38.287s +15 54'43.49'' DEC] 
Denebokab, x=1023, y=815 [19h25m29.9005s +03 06' 53.191'' DEC] 

Quelle fonction mathématique convertit l'emplacement des pixels en RA / DEC et vice versa? Remarques:

• Les étoiles brillantes sont des taches sur l'image; les coordonnées ci-dessus sont à peu près le centre de la goutte, mais peuvent être décalées de + -2 pixels.

• Je sais que je peux faire pivoter la sphère céleste de sorte que le centre de mon image ait des coordonnées polaires 0,0. La vraie question est donc "comment trouver cette rotation" (mais voir le point suivant).

• Si l'élévation / azimut était linéaire sur les images, ce serait plus facile (euh), mais ce n'est pas le cas: mesurer la distance angulaire avec des photographies

• Je peux fournir des emplacements de pixels de plus d'étoiles si cela aide. Je pense que 3 devraient suffire, mais je peux me tromper.

• J'ai essayé de choisir 3 étoiles qui étaient "réparties" sur l'image (car je pense que cela réduit les erreurs, je ne sais pas), mais je ne suis pas sûr d'avoir réussi.

• Je fais cela pour plusieurs photos et je voudrais une méthode générale.

• Faire cela m'aidera à identifier les étoiles plus faibles / les objets Messier / etc. dans l'image.

• Je suis sûr que beaucoup d'astrophotographes veulent le faire, mais n'ont pas trouvé de logiciel existant qui le fasse.

EDIT: Merci, whuber! La projection gnomonique est ce qui me manquait. J'avais déjà fait cela en supposant une transformation linéaire:

(* convert RA/DEC to xyz coords on celestial psuedo-sphere of radius 1 *) 
radecxyz[ra_,dec_] = 
{Cos[ra/12*Pi]*Cos[dec/180*Pi],Sin[ra/12*Pi]*Cos[dec/180*Pi],Sin[dec/180*Pi]}; 

(* I no longer have any idea how this works *) 
astrosolve[x_,y_,z_,xwid_,ywid_] := Module[{a,m,ans,nullans}, 
m=Array[a,{2,3}]; 
temp=Solve[{ 
m.radecxyz[x[[1]],x[[2]]]=={x[[3]]-xwid/2,x[[4]]-ywid/2}, 
m.radecxyz[y[[1]],y[[2]]]=={y[[3]]-xwid/2,y[[4]]-ywid/2}, 
m.radecxyz[z[[1]],z[[2]]]=={z[[3]]-xwid/2,z[[4]]-ywid/2} 
}]; 
ans = m /. Flatten[temp]; 
nullans=Flatten[NullSpace[ans]]; 
If[nullans.radecxyz[x[[1]],x[[2]]]<0,nullans=-nullans]; 
Return[{ans,nullans}]; 
]; 

où x, y et z étaient chacun des listes à 4 éléments comprenant une étoile RA, une déclinaison, une coordonnée x sur l'image et une coordonnée y sur l'image. xwid et ywid sont la largeur et la hauteur de l'image. Dans ce cas:

astrosolve[ 
 {19.8463886110, 8.8683219443, 782, 532}, 
 {20.6606352777, 15.9120805555, 311, 146}, 
 {19.4249723610, 3.1147752777, 1023, 815}, 
 1443, 998] 

{ 
 {{-2250.51, -1182.52, 385.689},  {-166.12, -543.746, -2376.73}},  
 {0.480698, -0.861509, 0.163497} 
} 

Maintenant, en faisant référence à "{-2250.51, -1182.52, 385.689}" comme $ frow, "{-166.12, -543.746, -2376.73}" comme $ srow et "{0.480698, -0.861509, 0.163497}" comme $ null, ce sous-programme PHP traduit RA / DEC en coordonnées xy:

# radecxy(ra,dec): converts ra/dec to x,y using a quasi-linear transformation 

function radecxy($ra,$dec) { 
    global $null,$frow,$srow,$xwid,$ywid; 
    list($x,$y,$z)=array(cos($dec)*cos($ra),cos($dec)*sin($ra),sin($dec)); 

    $dotprod=$null[0]*$x+$null[1]*$y+$null[2]*$z; 
    if ($dotprod<0) {return(array(-1,-1));}

 list($fx,$fy)  = array($frow[0]*$x+$frow[1]*$y+$frow[2]*$z,$srow[0]*$x+$srow[1]*$y+$srow[2]*$z); 
    $fx+=$xwid/2; 
    $fy+=$ywid/2; 
    if ($fx<0 || $fy<0 || $fx>$xwid || $fy>$ywid) { 
        return(array(-1,-1)); 
    } else { 
        return(array($fx,$fy)); 
    } 
} 

Malheureusement, je ne sais plus pourquoi cela fonctionne, mais l'utiliser + ajouter des positions d'étoiles connues donne des résultats tolérables (utilisez "voir l'image" pour la voir en taille réelle):

Cependant, comme vous pouvez le voir, les résultats ne sont pas parfaits, ce qui me convainc qu'une transformation linéaire n'était pas la bonne réponse. Je pense que gnomonic pourrait être le Graal que je cherchais.

Je décrirai une approche rigoureuse et indiquerai quel logiciel peut l'aider. La plupart de ces informations seront tangentielles aux intérêts du site de photographie, mais comme il existe des informations utiles qui s'appliquent à toutes les circonstances dans lesquelles les emplacements seront estimés à partir des mesures sur une image, ce site semble un endroit raisonnable pour une telle analyse.

Prendre une image (avec un objectif corrigé pour la distorsion) projette la sphère céleste à travers le point focal de l'objectif sur le plan du capteur. Il s'agit d' un aspect oblique d'une projection gnomonique .

Mathématiquement, la conversion de (RA, DEC) passe par une série d'étapes:

1. Convertissez (RA, DEC) en coordonnées sphériques. RA doit être converti des heures-minutes-secondes en degrés (ou radians) et DEC doit être converti des degrés-minutes secondes en degrés (ou radians), en se rappelant que c'est l'élévation au-dessus du plan, pas l'angle du pôle nord (qui est la convention de coordonnées sphériques habituelle). Les deux conversions sont arithmétiques simples.

2. Calculez les coordonnées (x, y, z) pour les coordonnées sphériques des étoiles. Il s'agit d'une conversion de coordonnées standard (impliquant une trigonométrie simple).

3. Faites pivoter la sphère céleste pour aligner ses pôles avec l'axe de l'objectif. Il s'agit d'une transformation linéaire.

4. Faites pivoter la sphère céleste autour de ses pôles pour vous conformer à l'orientation de la caméra (une autre transformation linéaire).

5. En plaçant le plan d'imagerie à une hauteur constante z au-dessus du point focal, attirer les rayons lumineux des étoiles à (x, y, z) à travers le point focal jusqu'à ce qu'ils interceptent le plan. (Il s'agit de la projection gnomonique et, de par sa nature, elle est projective et non linéaire.)

[Sur la figure, qui est destinée à être une coupe transversale plane à travers l'axe de la lentille,

• A est le point focal.
• Le demi-cercle BCD est la partie visible de la sphère céleste.
• Points CA le long de l'axe de la lentille.
• E, F et G sont des emplacements étoiles.
• EE, FF et GG sont leurs emplacements correspondants sur la sphère céleste (invisible).
• E ', F' et G 'sont leurs images sur le capteur KL (de sorte que EE', FF 'et GG' sont des trajets de rayons lumineux des étoiles vers le capteur).
• AD est l'horizon à partir duquel la déclinaison est mesurée.
• Alpha est la déclinaison de l'étoile E (ou, de manière équivalente, une coordonnée angulaire de EE). Les étoiles F et G ont des déclinaisons similaires (non représentées).

Notre tâche consiste à trouver la relation mathématique entre les coordonnées angulaires pour E, F et G - qui sont supposées connues avec une grande précision - telles que l'alpha, et les coordonnées de leurs images E ', F' et G ', mesurées en pixels le long du capteur. Une fois trouvée, cette relation peut être inversée comme décrit ci-dessous pour estimer les coordonnées angulaires des objets célestes à partir des positions de leurs images sur le capteur. Pour plus de simplicité, le grossissement de l'objectif n'est pas illustré. Avec un objectif sans distorsion, cela aura pour effet de redimensionner uniformément les coordonnées de E ', F' et G 'par rapport au centre du capteur.]

Cette procédure décrit comment la lumière passe d'une étoile sur le capteur pour une lentille simple parfaite. Il implique ces paramètres (inconnus), qui devront être déterminés:

• Trois angles en (3) et (4) décrivant l'orientation de l'objectif et de la caméra.

• Un facteur d'échelle dans (5) décrivant les effets combinés de la taille du capteur, de la distance du point focal et du grossissement de l'objectif.

En raison de la projection (5), il s'agit d'une transformation complexe et non linéaire en général, mais elle a une description mathématique définie. Si nous laissons x = (RA, DEC) désigner la position d'une étoile, que thêta représente les quatre paramètres du processus d'imagerie, et que y = (colonne, ligne) représente les coordonnées des pixels, alors nous pouvons écrire de manière abstraite mais plus simplement

y = f (x, thêta).

Ensuite - et c'est très important - nous devons tenir compte des erreurs. Les étoiles imagées ne sont pas dans des emplacements précis. Nous devons donc inclure un terme d'erreur dans notre formule et il est classique (depuis environ 1800) de modéliser cette erreur de manière probabiliste. La nouvelle formule est

y = f (x, thêta) + e

Lorsque l'objectif est exempt de distorsion, la valeur attendue de e est 0 et son écart-type ( sigma ) mesure la taille d'erreur typique. Il est raisonnable de supposer que les e sont à peu près normalement distribués, avec des écarts-types approximativement égaux (ce qui n'est pas vrai, mais pour une analyse initiale, c'est une hypothèse raisonnable) et nous pouvons espérer que ces erreurs sont statistiquement indépendantes les unes des autres (ce qui n'est pas encore vrai mais c'est une bonne hypothèse de départ). Cela justifie une solution des moindres carrés utilisant le maximum de vraisemblance . Jusqu'à une constante universelle, dont nous n'avons pas besoin de connaître la valeur, la probabilité log d'une observation particulière (x, y) est égale à

- | f (x, thêta) - y | ^ 2 / (2 sigma ^ 2) - 2 log (sigma).

(Les barres de valeur absolue indiquent la distance euclidienne dans le plan d'imagerie, calculée comme d'habitude avec le théorème de Pythagore.)

En raison de l'indépendance supposée des erreurs, la probabilité logarithmique de l'ensemble de données pour une image est la somme de ces valeurs. Il s'agit de la «vraisemblance logarithmique». Les estimations du maximum de vraisemblance (ML) des paramètres thêta et sigma (cinq nombres en tout) sont les valeurs qui maximisent la vraisemblance logarithmique.

Nous pouvons et devons aller plus loin. La théorie du ML montre également comment obtenir des intervalles de confiance pour les estimations. Intuitivement, les erreurs dans nos observations créent une petite incertitude dans les valeurs conjointes des angles, du facteur d'échelle et de l'écart-type. Nous avons besoin de ces valeurs pour estimer RA et DEC pour tous les pixels de notre image. En utilisant des valeurs incertaines, ce qui est inévitable, nous obtiendrons des résultats incertains. De plus, si nous identifions un pixel dans notre image en regardant une goutte de lumière diffuse (dispersée sur environ pi * sigma ^ 2 pixels au total), il y aura une incertitude supplémentaire dans les coordonnées des pixels. Ensemble, ces deux formes d'incertitude se combinent. Cela impliquel'incertitude nette dans l'estimation de la RA et de la DEC de toute goutte de lumière sur l'image est plus grande que vous ne le pensez.

Enfin, lorsque vous prenez une mesure de l'image et que vous l'utilisez pour estimer les véritables coordonnées d'une étoile ou d'un objet céleste, vous effectuez une régression inverse , qui est une forme d'étalonnage de l'instrument. La régression inverse est une procédure pour tenir compte des incertitudes que je viens de décrire. Sa sortie comprend bien sûr les coordonnées estimées des étoiles pour n'importe quel blob de pixels sur l'image. Il comprend également un anneau de coordonnées autour de cette estimation qui sont également cohérentes avec l'emplacement de cette goutte. (Il s'agit d'un "intervalle de prédiction inverse" ou d'un ensemble de "limites fiduciales" pour les RA et DEC du blob.) En pratique, si vous consultez un catalogue d'objets célestes, vous pouvez utiliser cet anneau pour rechercher tous les objets connus qui sont cohérentes avec les informations de votre image. De toute évidence, cela peut être plus utile qu'une procédure simpliste qui estime - parfois incorrectement - un seul ensemble de coordonnées.

En résumé, ce qu'il faut ici, c'est un logiciel pour

• Effectuez l'optimisation non linéaire requise par ML.

• Estimer les erreurs types dans les estimations.

• Effectuez une régression inverse.

Une expertise avec un logiciel approprié, comme la commande ML de Stata ou Mathematica , est essentielle si vous codez vous-même.

Quelle que soit votre expertise, voici quelques conclusions que vous pouvez utiliser dans vos stratégies d'imagerie:

• La précision de l'image pour obtenir une correction sur n'importe quel objet ne peut jamais être supérieure à l'imprécision inhérente à l'imagerie (telle que mesurée par sigma , la taille typique d'un point de lumière sur l'image).

• Vous pouvez vous rapprocher de ce niveau de précision en identifiant de nombreuses étoiles connues, pas seulement trois. Cela réduit l'incertitude dans la transformation du ciel en image presque à zéro si vous avez suffisamment d'étoiles connues situées dans l'image.

• Il est vrai que vous souhaitez que les étoiles de référence soient réparties sur l'image. Il est également crucial de ne pas les aligner (ce qui est malheureusement le cas pour les trois lieux cités dans la question). Si vous pouvez vous permettre de localiser seulement trois étoiles, placez-les dans un joli triangle. Lorsque les étoiles s'alignent, l'analyse statistique indique qu'il existe une énorme incertitude sur les emplacements dans les directions perpendiculaires à la ligne. Dans cet exemple particulier, l'erreur estimée ( sigma ) fait des centaines de pixels de large. Avoir une étoile de plus pour faire un bon triangle devrait réduire cette erreur à un ou deux pixels.

Quelques pensées d'adieu:

• Il est possible de détecter et même de corriger les aberrations des lentilles en effectuant une analyse statistique plus approfondie. L'idée est de tracer les écarts entre les emplacements prévus et réels des étoiles sur l'image. Cela s'apparente à des données cartographiques de "déformation" ou de " géoréférencement ". En tant que solution rapide et sale, vous pouvez mettre en service un SIG ou un logiciel de traitement d'image (tel que ENVI ) pour géoréférencer (ou astroréférencer) n'importe quelle image. Un tel logiciel n'effectue généralement pas d'estimations ML des transformations projectives, mais il peut effectuer des approximations polynomiales d'ordre élevé. Une transformation polynomiale d'ordre 2 ou d'ordre 3 peut faire assez bien le travail, selon votre application.

• Il est possible d'améliorer la précision en combinant plusieurs images des mêmes objets.

Pour ce faire, avec le même degré de précision que les astronomes professionnels, il serait en effet difficile. Cela nécessiterait une caractérisation extrêmement précise des distorsions produites par votre objectif et des imperfections du capteur de votre appareil photo. Cependant, vous n'avez probablement pas besoin de ce degré de précision. Il devrait suffire que vous supposiez que votre objectif n'introduit pas de grandes quantités de distorsion (ce qui est une bonne hypothèse pour un objectif de qualité) et que le capteur de votre appareil photo est assez proche d'une grille parfaitement régulière (ce qui est une très bonne hypothèse pour même un appareil photo bon marché).

Il ne reste plus qu'à déterminer la transformation des coordonnées qui décrit l'orientation de la caméra, c'est-à-dire la direction dans laquelle elle a été pointée et le degré de rotation.

Ce que vous recherchez alors s'appelle une transformation affine ou une carte affine. Ce qui est juste un nom de fantaisie pour une matrice par laquelle vous multiplieriez vos coordonnées de pixels pour obtenir vos coordonnées astronomiques. Dans le cas d'une carte affine, cette transformation peut inclure n'importe quel degré de rotation, d'échelle, de cisaillement et de translation.

La signification du composant de rotation est assez évidente. Le facteur d'échelle décrit simplement la quantité de ciel couverte par chaque pixel en termes de RA / Dec. Le cisaillement est une transformation qui ferait que l'image d'un rectangle devienne un parallélogramme, mais il ne devrait y avoir aucun de ces effets dans une image d'objets à l'infini (comme les étoiles). Enfin, le composant de traduction simple ajoute un décalage pour tenir compte du fait que le pixel (x = 0, y = 0) dans votre image ne correspond probablement pas à (RA = 0, Dec = 0).

Parce que vous avez 3 étoiles de référence dans votre image, vous avez suffisamment d'informations pour calculer la relation entre vos coordonnées de pixels et le RA / Dec que vous recherchez. Cela se ferait par ajustement des moindres carrés linéaires (pas les moindres carrés non linéaires comme mentionné ci-dessus) pour déterminer les valeurs des composants de la matrice qui correspondent le mieux à vos coordonnées de pixels à la RA / Dec connue des étoiles de référence. Une fois la matrice établie, vous pouvez ensuite l'appliquer aux coordonnées en pixels des autres étoiles pour obtenir leur RA / Dec.

Bien que je puisse le faire relativement facilement, je ne suis malheureusement pas sûr de savoir comment vous aider à le faire. Cela impliquerait des compétences mathématiques qui dépassent un peu le cadre de photo.SE. Je suis ingénieur optique, mais je ne suis pas très photographe; le logiciel que j'utiliserais pour cela est conçu pour que les ingénieurs effectuent des calculs numériques intensifs, et n'est pas vraiment un outil photographique du tout. Il peut y avoir des moyens de le faire en utilisant des logiciels destinés aux photographes, mais je ne les connais pas.