Mapeo de píxel a RA / DEC en astrofotografía digitalizada

Tengo una imagen de 1443x998 de las estrellas (tomada con una cámara de 35 mm y luego escaneada) con las siguientes estrellas en las siguientes ubicaciones de píxeles:

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] 

¿Qué función matemática convierte la ubicación de píxeles en RA / DEC y viceversa? Notas:

• Las estrellas brillantes son manchas en la imagen; las coordenadas anteriores son aproximadamente el centro de la burbuja, pero pueden estar apagadas en + -2 píxeles.

• Sé que puedo rotar la esfera celeste para que el centro de mi imagen tenga coordenadas polares 0,0. Entonces, la verdadera pregunta es "cómo encontrar esta rotación" (pero vea el siguiente punto).

• Si la elevación / acimut fuera lineal en las imágenes, esto sería fácil (er), pero no lo son: medir la distancia angular con fotografías

• Puedo proporcionar ubicaciones de píxeles de más estrellas si eso ayuda. Creo que 3 debería ser suficiente, pero podría estar equivocado.

• Traté de elegir 3 estrellas que se "extendieron" en la imagen (porque creo que eso reduce el error, no estoy seguro), pero no estoy seguro de haber tenido éxito.

• Estoy haciendo esto para varias imágenes y me gustaría un método general.

• Hacer esto me ayudará a identificar estrellas más débiles / objetos Messier / etc en la imagen.

• Estoy seguro de que muchos astrofotógrafos quieren hacer esto, pero no han encontrado ningún software existente que lo haga.

EDITAR: ¡Gracias, whuber! La proyección gnomónica es lo que me faltaba. Ya había hecho esto suponiendo una transformación lineal:

(* 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}]; 
]; 

donde x, y y z eran cada una de las listas de 4 elementos que consisten en una estrella RA, declinación, coordenada x en la imagen y coordenada y en la imagen. xwid e ywid son el ancho y la altura de la imagen. En este caso:

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} 
} 

Ahora, refiriéndose a "{-2250.51, -1182.52, 385.689}" como $ frow, "{-166.12, -543.746, -2376.73}" como $ srow y "{0.480698, -0.861509, 0.163497}" como $ null, Esta subrutina PHP traduce RA / DEC a coordenadas 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)); 
    } 
} 

Lamentablemente, ya no tengo idea de por qué esto funciona, pero usarlo + agregar posiciones conocidas de estrellas produce resultados tolerables (use "ver imagen" para verlo a tamaño completo):

Sin embargo, como puede ver, los resultados no son perfectos, convenciéndome de que una transformación lineal no era la respuesta correcta. Creo que gnomonic podría ser el grial que estaba buscando.

Esbozaré un enfoque riguroso e indicaré qué software puede ayudarlo. La mayor parte de esto será tangencial a los intereses del sitio de fotografía, pero debido a que hay algunas ideas útiles que se aplican a cualquier circunstancia en la que se estimarán ubicaciones a partir de mediciones en una imagen, este sitio parece un lugar razonable para tal análisis.

Tomar una imagen (con una lente que ha sido corregida por distorsión) proyecta la esfera celeste a través del punto focal de la lente en el plano del sensor. Este es un aspecto oblicuo de una proyección gnomónica .

Matemáticamente, la conversión de (RA, DEC) continúa a través de una serie de pasos:

1. Convierta (RA, DEC) a coordenadas esféricas. RA se debe convertir de horas-minutos-segundos a grados (o radianes) y DEC se debe convertir de grados-minutos segundos a grados (o radianes), recordando que es una elevación sobre el plano, no un ángulo desde el polo norte (que es la convención de coordenadas esférica habitual). Ambas conversiones son aritméticas simples.

2. Calcule las coordenadas (x, y, z) para las coordenadas esféricas de las estrellas. Esta es una conversión de coordenadas estándar (que involucra trigonometría simple).

3. Gire la esfera celeste para alinear sus polos con el eje de la lente. Esta es una transformación lineal.

4. Gire la esfera celeste alrededor de sus polos para ajustarse a la orientación de la cámara (otra transformación lineal).

5. Colocando el plano de imagen a una altura constante z sobre el punto focal, dibuje rayos de luz de las estrellas en (x, y, z) a través del punto focal hasta que intercepten el plano. (Esta es la proyección gnomónica y, por su naturaleza, es proyectiva y no lineal).

[En la figura, que pretende ser una sección transversal plana a través del eje de la lente,

• A es el punto focal.
• El semicírculo BCD es la parte visible de la esfera celeste.
• Puntos de CA a lo largo del eje de la lente.
• E, F y G son ubicaciones de estrellas.
• EE, FF y GG son sus ubicaciones correspondientes en la esfera celeste (invisible).
• E ', F' y G 'son sus imágenes en el sensor KL (de modo que EE', FF 'y GG' son trayectorias de rayos de luz desde las estrellas hasta el sensor).
• AD es el horizonte desde el cual se mide la declinación.
• Alfa es la declinación de la estrella E (o, equivalentemente, una coordenada angular de EE). Las estrellas F y G tienen declinaciones similares (no se muestran).

Nuestra tarea es encontrar la relación matemática entre las coordenadas angulares para E, F y G, que se supone que se conocen con gran precisión, como alfa, y las coordenadas de sus imágenes E ', F' y G ', medidas en píxeles a lo largo del sensor. Una vez encontrada, esta relación puede invertirse como se describe a continuación para estimar las coordenadas angulares de los objetos celestes a partir de las posiciones de sus imágenes en el sensor. No se muestra, por simplicidad, el aumento de la lente. Con una lente sin distorsión, esto tendrá el efecto de reescalar uniformemente las coordenadas de E ', F' y G 'en relación con el centro del sensor.]

Este procedimiento describe cómo la luz pasa de una estrella al sensor para obtener una lente simple perfecta. Implica estos parámetros (desconocidos), que deberán determinarse:

• Tres ángulos en (3) y (4) que describen la lente y la orientación de la cámara.

• Un factor de escala en (5) que describe los efectos combinados del tamaño del sensor, la distancia desde el punto focal y el aumento de la lente.

Debido a la proyección (5), esta es una transformación compleja, no lineal en general, pero tiene una descripción matemática definida. Si dejamos que x = (RA, DEC) designe la posición de una estrella, dejemos que theta represente los cuatro parámetros para el proceso de obtención de imágenes y dejemos que y = (columna, fila) represente las coordenadas de los píxeles, entonces podemos escribir de manera abstracta pero más simple

y = f (x, theta).

A continuación, y esto es muy importante, debemos tener en cuenta los errores. Las estrellas fotografiadas no están en ubicaciones precisas. Por lo tanto, tenemos que incluir un término de error en nuestra fórmula y es convencional (desde aproximadamente 1800) modelar este error probabilísticamente. La nueva fórmula es

y = f (x, theta) + e

Cuando la lente está libre de distorsión, el valor esperado de e es 0 y su desviación estándar ( sigma ) mide el tamaño típico de error. Es razonable suponer que las e están distribuidas aproximadamente de manera normal, con desviaciones estándar aproximadamente iguales (lo cual no es cierto, pero para un análisis inicial es una suposición razonable) y podemos esperar que estos errores sean estadísticamente independientes entre sí (lo que de nuevo no es cierto pero es una buena suposición inicial). Esto justifica una solución de mínimos cuadrados usando la máxima probabilidad . Hasta una constante universal, cuyo valor no necesitamos saber, la probabilidad logarítmica de cualquier observación particular (x, y) es igual a

- | f (x, theta) - y | ^ 2 / (2 sigma ^ 2) - 2 log (sigma).

(Las barras de valor absoluto denotan la distancia euclidiana en el plano de imagen, calculada como de costumbre con el teorema de Pitágoras).

En virtud de la supuesta independencia de los errores, la probabilidad de registro del conjunto de datos para una imagen es la suma de estos valores. Esta es la "probabilidad de registro". Las estimaciones de máxima verosimilitud (ML) de los parámetros theta y sigma (cinco números en total) son aquellos valores que maximizan la probabilidad logarítmica.

Podemos y debemos ir más allá. La teoría de ML también muestra cómo obtener intervalos de confianza para las estimaciones. Intuitivamente, los errores en nuestras observaciones crean una pequeña incertidumbre en los valores conjuntos de los ángulos, el factor de escala y la desviación estándar. Necesitamos estos valores para estimar RA y DEC para cualquier píxel en nuestra imagen. Al usar valores inciertos, lo cual es inevitable, obtendremos resultados inciertos. Además, si identificamos un píxel en nuestra imagen al observar una gota de luz difusa (dispersa en aproximadamente pi * sigma ^ 2 píxeles en total), habrá incertidumbre adicional en las coordenadas de los píxeles. Colectivamente, estas dos formas de incertidumbre se combinan. Esto implicaLa incertidumbre neta en la estimación de la RA y la DEC de cualquier gota de luz en la imagen es mayor de lo que se podría suponer.

Finalmente, cuando tomas una medida de la imagen y la usas para estimar las coordenadas verdaderas de una estrella u objeto celeste, estás haciendo una regresión inversa , que es una forma de calibración del instrumento. La regresión inversa es un procedimiento para explicar las incertidumbres que acabo de describir. Su salida, por supuesto, incluye las coordenadas de estrella estimadas para cualquier burbuja de píxeles en la imagen. También incluye un anillo de coordenadas alrededor de esa estimación que también son consistentes con la ubicación de esa burbuja. (Este es un "intervalo de predicción inversa" conjunto o un conjunto de "límites fiduciales" para la RA y la DEC del blob.) En la práctica, si consulta un catálogo de objetos celestes, puede usar este anillo para buscar todos los objetos conocidos que son consistentes con la información en su imagen. Claramente, esto puede ser más valioso que un procedimiento simplista que estima, a veces incorrectamente, solo un conjunto único de coordenadas.

En resumen, lo que se necesita aquí es software para

• Realice la optimización no lineal requerida por ML.

• Estimar errores estándar en las estimaciones.

• Realizar regresión inversa.

La experiencia con el software apropiado, como el comando ML de Stata o Mathematica , es esencial si está codificando esto usted mismo.

Independientemente de su experiencia, aquí hay algunas conclusiones que puede usar en sus estrategias de imagen:

• La precisión de la imagen para obtener una solución en cualquier objeto nunca puede ser mayor que la imprecisión inherente en la imagen (medida por sigma , el tamaño típico de un punto de luz en la imagen).

• Puede acercarse a este nivel de precisión identificando muchas estrellas conocidas, no solo tres. Esto reduce la incertidumbre en la transformación del cielo a la imagen casi a cero si tiene suficientes estrellas conocidas ubicadas en la imagen.

• Es correcto que desee que las estrellas de referencia se extiendan por la imagen. También es crucial que no estén alineados (lo cual, desafortunadamente, es el caso con las tres ubicaciones dadas en la pregunta). Si puede permitirse localizar solo tres estrellas, consígalas en un bonito triángulo. Cuando las estrellas se alinean, el análisis estadístico indica que existe una gran incertidumbre sobre las ubicaciones en direcciones perpendiculares a la línea. En este ejemplo particular, el error estimado ( sigma ) es de cientos de píxeles de ancho. Tener una estrella más para hacer un buen triángulo debería reducir este error a uno o dos píxeles.

Algunos pensamientos de despedida:

• Es posible detectar e incluso corregir las aberraciones de la lente llevando a cabo un análisis estadístico más extenso. La idea es trazar las desviaciones entre las ubicaciones esperadas y reales de las estrellas en la imagen. Esto es similar a los datos del mapa "deformación" o " georreferenciación ". Como solución rápida y sucia, puede presionar GIS o software de procesamiento de imágenes (como ENVI ) en servicio para georreferenciar (o astroreferenciar) cualquier imagen. Tal software generalmente no realiza estimaciones ML de transformaciones proyectivas, pero puede hacer aproximaciones polinómicas de alto orden. Una transformación polinómica de orden 2 u orden 3 podría hacer un buen trabajo, dependiendo de su aplicación.

• Es posible mejorar la precisión combinando múltiples imágenes de los mismos objetos.

Hacer esto con el mismo grado de precisión que hacen los astrónomos profesionales sería realmente difícil. Requeriría que tenga una caracterización extremadamente precisa de las distorsiones producidas por su lente y las imperfecciones en el sensor de su cámara. Sin embargo, probablemente no necesite ese grado de precisión. Debería ser suficiente para suponer que su lente no introduce grandes cantidades de distorsión (lo cual es una buena suposición para una lente de calidad) y que el sensor de su cámara está bastante cerca de una cuadrícula perfectamente regular (lo cual es una muy buena suposición para incluso una cámara barata).

Todo lo que queda es resolver la transformación de coordenadas que describe la orientación de la cámara, es decir, la dirección en la que apuntaba y el grado en que giraba.

Entonces, lo que está buscando se llama transformación afín o mapa afín. El cual es solo un nombre elegante para una matriz por la cual multiplicarías tus coordenadas de píxeles para obtener tus coordenadas astronómicas. En el caso de un mapa afín, esta transformación puede incluir cualquier grado de rotación, escala, corte y traslación.

El significado del componente de rotación es bastante obvio. El factor de escala simplemente describe la cantidad de cielo que cubre cada píxel en términos de RA / Dec. La cizalladura es una transformación que haría que la imagen de un rectángulo se convirtiera en un paralelogramo, pero no debería haber ninguno de estos efectos en una imagen de objetos en el infinito (como las estrellas). Por último, el componente de traducción simple agrega un desplazamiento para tener en cuenta el hecho de que el píxel (x = 0, y = 0) en su imagen probablemente no corresponde a (RA = 0, Dec = 0).

Debido a que tiene 3 estrellas de referencia en su imagen, tiene suficiente información para calcular la relación entre sus coordenadas de píxeles y el RA / Dec que está buscando. Esto se haría mediante el ajuste de mínimos cuadrados lineales (no mínimos cuadrados no lineales como se mencionó anteriormente) para determinar los valores de los componentes de la matriz que mejor coincidan con sus coordenadas de píxeles con la RA / Dec conocida de las estrellas de referencia. Una vez establecida la matriz, puede aplicarla a las coordenadas de píxeles de otras estrellas para obtener su RA / Dec.

Si bien podría hacerlo con relativa facilidad, lamentablemente no estoy seguro de cómo ayudarlo a hacerlo. Implicaría alguna habilidad matemática que está un poco más allá del alcance de la foto. Soy ingeniero óptico, pero no soy muy fotógrafo; El software que usaría para esto está diseñado para que los ingenieros realicen cálculos numéricos de alta resistencia, y no es realmente una herramienta fotográfica. Puede haber maneras de hacer esto usando paquetes de software dirigidos a fotógrafos, pero no los conozco.