Mapeamento de pixel para RA / DEC em astrofotografia digitalizada

Tenho uma foto das estrelas em 1443x998 (tirada com uma câmera de 35 mm e depois digitalizada) com as seguintes estrelas nos seguintes locais de pixels:

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] 

Que função matemática converte a localização do pixel em RA / DEC e vice-versa? Notas:

• Estrelas brilhantes são bolhas na imagem; as coordenadas acima são aproximadamente o centro do blob, mas podem estar desativadas em + -2 pixels.

• Eu sei que posso girar a esfera celeste para que o centro da minha imagem tenha coordenadas polares 0,0. Portanto, a verdadeira questão é "como encontrar essa rotação" (mas veja o próximo ponto).

• Se a elevação / azimute fosse linear nas imagens, seria mais fácil, mas não é: Medir a distância angular com fotografias

• Posso fornecer locais de pixel de mais estrelas, se isso ajudar. Acredito que 3 deve ser suficiente, mas posso estar errado.

• Tentei escolher 3 estrelas que estavam "espalhadas" pela imagem (porque acho que isso reduz o erro, não tenho certeza), mas não tenho certeza se consegui.

• Estou fazendo isso para várias fotos e gostaria de um método geral.

• Fazer isso me ajudará a identificar estrelas mais fracas / objetos Messier / etc na imagem.

• Tenho certeza que muitos astrofotógrafos querem fazer isso, mas não encontraram nenhum software existente que faça isso.

EDIT: Obrigado, whuber! A projeção gnomônica é o que estava faltando. Eu já tinha feito isso assumindo uma transformação linear:

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

onde x, ye z eram, cada uma, listas de 4 elementos consistindo em RA de estrelas, declinação, coordenada x na imagem e coordenada y na imagem. xwid e ywid são a largura e a altura da imagem. Nesse 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} 
} 

Agora, referindo-se a "{-2250.51, -1182.52, 385.689}" como $ frow ", {-166.12, -543.746, -2376.73}" como $ srow e "{0.480698, -0.861509, 0.163497}" como $ null, esta sub-rotina PHP converte as coordenadas RA / DEC em 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)); 
    } 
} 

Infelizmente, eu não tenho mais idéia do por que isso funciona, mas usá-lo + adicionar posições em estrela conhecidas gera resultados toleráveis ​​(use "ver imagem" para ver em tamanho real):

No entanto, como você pode ver, os resultados não são perfeitos, me convencendo de que uma transformação linear não era a resposta certa. Eu acho que gnomônico pode ser o graal que eu estava procurando.

Vou descrever uma abordagem rigorosa e indicar qual software pode ajudar com ela. A maior parte disso será tangencial aos interesses do site de fotografia, mas como existem algumas informações úteis que se aplicam a qualquer circunstância em que os locais serão estimados a partir de medições em uma imagem, esse site parece um local razoável para essa análise.

Tirar uma imagem (com uma lente que foi corrigida para distorção) projeta a esfera celeste através do ponto focal da lente no plano do sensor. Este é um aspecto oblíquo de uma projeção gnomônica .

Matematicamente, a conversão de (RA, DEC) prossegue através de uma série de etapas:

1. Converter (RA, DEC) em coordenadas esféricas. O RA deve ser convertido de horas-minutos-segundos em graus (ou radianos) e o DEC deve ser convertido de graus-minutos segundos em graus (ou radianos), lembrando que é a elevação acima do plano, não o ângulo do pólo norte (que é a convenção usual de coordenadas esféricas). Ambas as conversões são aritméticas simples.

2. Cálculo (x, y, z) de coordenadas para as coordenadas esféricas das estrelas. Esta é uma conversão de coordenadas padrão (envolvendo trigonometria simples).

3. Gire a esfera celeste para alinhar seus pólos com o eixo da lente. Esta é uma transformação linear.

4. Gire a esfera celeste em torno de seus polos para se adaptar à orientação da câmera (outra transformação linear).

5. Colocando o plano de imagem a uma altura constante z acima do ponto focal, extraia raios de luz das estrelas em (x, y, z) através do ponto focal até que interceptem o plano. (Essa é a projeção gnomônica e, por sua natureza, é projetiva e não linear.)

[Na figura, que se destina a ser uma seção transversal plana através do eixo da lente,

• A é o ponto focal.
• O semicírculo BCD é a parte visível da esfera celeste.
• Pontos CA ao longo do eixo da lente.
• E, F e G são locais em estrela.
• EE, FF e GG são seus locais correspondentes na esfera celeste (invisível).
• E ', F' e G 'são suas imagens no sensor KL (de modo que EE', FF 'e GG' são caminhos de raios de luz das estrelas ao sensor).
• AD é o horizonte a partir do qual a declinação é medida.
• Alfa é a declinação da estrela E (ou, equivalentemente, uma coordenada angular de EE). As estrelas F e G têm declínios semelhantes (não mostrados).

Nossa tarefa é encontrar a relação matemática entre as coordenadas angulares de E, F e G - que são consideradas conhecidas com alta precisão - como alfa, e as coordenadas de suas imagens E ', F' e G ', medidas em pixels ao longo do sensor. Uma vez encontrada, essa relação pode ser invertida como descrito abaixo para estimar coordenadas angulares de objetos celestes a partir das posições de suas imagens no sensor. Não é mostrado, por simplicidade, a ampliação da lente. Com uma lente sem distorção, isso terá o efeito de redimensionar uniformemente as coordenadas de E ', F' e G 'em relação ao centro do sensor.]

Este procedimento descreve como a luz passa de uma estrela para o sensor para obter uma lente simples e perfeita. Envolve estes parâmetros (desconhecidos), que precisarão ser determinados:

• Três ângulos em (3) e (4) descrevendo a orientação da lente e da câmera.

• Um fator de escala em (5) que descreve os efeitos combinados do tamanho do sensor, distância do ponto focal e ampliação da lente.

Devido à projeção (5), essa é uma transformação complexa e não linear em geral, mas tem uma descrição matemática definida. Se deixarmos x = (RA, DEC) designar a posição de uma estrela, deixe que teta represente os quatro parâmetros para o processo de geração de imagens, e deixe que y = (coluna, linha) represente as coordenadas do pixel, então podemos abstrair, mas com mais simplicidade, escrever

y = f (x, teta).

Em seguida - e isso é muito importante - precisamos levar em conta erros. As estrelas fotografadas não estão em locais precisos. Portanto, temos que incluir um termo de erro em nossa fórmula e é convencional (desde cerca de 1800) modelar esse erro de forma probabilística. A nova fórmula é

y = f (x, teta) + e

Quando a lente está livre de distorção, o valor esperado de e é 0 e seu desvio padrão ( sigma ) mede o tamanho típico do erro. É razoável supor que os e são aproximadamente distribuídos normalmente, com desvios padrão aproximadamente iguais (o que não é verdade, mas para uma análise inicial é uma suposição razoável) e podemos esperar que esses erros sejam estatisticamente independentes um do outro (o que novamente não é verdade mas é uma boa suposição inicial). Isso justifica uma solução de mínimos quadrados usando a probabilidade máxima . Até uma constante universal, cujo valor não precisamos saber, a probabilidade logarítmica de qualquer observação em particular (x, y) é igual a

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

(As barras de valor absoluto denotam a distância euclidiana no plano da imagem, calculada como de costume com o Teorema de Pitágoras.)

Em virtude da suposta independência de erros, a probabilidade de log do conjunto de dados para uma imagem é a soma desses valores. Essa é a "probabilidade do log". As estimativas de máxima verossimilhança (ML) dos parâmetros theta e sigma (cinco números ao todo) são aqueles valores que maximizam a probabilidade logarítmica.

Podemos e devemos ir além. A teoria do BC também mostra como obter intervalos de confiança para as estimativas. Intuitivamente, os erros em nossas observações criam um pouco de incerteza nos valores conjuntos dos ângulos, no fator de escala e no desvio padrão. Precisamos desses valores para estimar RA e DEC para qualquer pixel em nossa imagem. Usando valores incertos, o que é inevitável, obteremos resultados incertos. Além disso, se identificarmos um pixel em nossa imagem observando uma mancha difusa de luz (espalhada por aproximadamente pi * sigma ^ 2 pixels no total), haverá incerteza adicional nas coordenadas do pixel. Coletivamente, essas duas formas de incerteza se combinam. Isso implicaa incerteza líquida na estimativa do RA e DEC de qualquer mancha de luz na imagem é maior do que você imagina.

Finalmente, ao tirar uma medida da imagem e usá-la para estimar as coordenadas reais de uma estrela ou objeto celeste, você está fazendo uma regressão inversa , que é uma forma de calibração do instrumento. A regressão inversa é um procedimento para explicar as incertezas que acabei de descrever. É claro que sua saída inclui as coordenadas estimadas em estrela para qualquer mancha de pixels na imagem. Ele também inclui um anel de coordenadas em torno dessa estimativa que também são consistentes com a localização desse blob. (Este é um "intervalo de previsão inversa" conjunto ou um conjunto de "limites fiduciais" para o RA e o DEC do blob.) Na prática, se você consultar um catálogo de objetos celestes, poderá usar esse anel para procurar todos os objetos conhecidos que são consistentes com as informações da sua imagem. Claramente, isso pode ser mais valioso do que um procedimento simplista que estima - às vezes incorretamente - apenas um único conjunto de coordenadas.

Em resumo, o que é necessário aqui é um software para

• Execute a otimização não linear exigida pelo ML.

• Estimar erros padrão nas estimativas.

• Execute regressão inversa.

A experiência com o software apropriado, como o comando ML da Stata ou o Mathematica , é essencial se você mesmo estiver codificando isso.

Independentemente de sua experiência, aqui estão algumas conclusões que você pode usar em suas estratégias de imagem:

• A precisão da imagem para obter uma correção em qualquer objeto nunca pode ser maior que a imprecisão inerente na imagem (medida por sigma , o tamanho típico de um ponto de luz na imagem).

• Você pode se aproximar desse nível de precisão identificando muitas estrelas conhecidas, não apenas três. Isso reduz a incerteza na transformação de céu em imagem quase para zero se você tiver estrelas conhecidas suficientes localizadas na imagem.

• É correto que você queira que as estrelas de referência sejam espalhadas pela imagem. Também é crucial que eles não estejam alinhados (o que, infelizmente, é o caso dos três locais indicados na pergunta). Se você puder localizar apenas três estrelas, coloque-as em um triângulo agradável. Quando as estrelas se alinham, a análise estatística indica que há uma enorme incerteza sobre os locais nas direções perpendiculares à linha. Neste exemplo em particular, o erro estimado ( sigma ) tem centenas de pixels de largura. Ter mais uma estrela para formar um bom triângulo deve reduzir esse erro a um ou dois pixels.

Alguns pensamentos de despedida:

• É possível detectar e até corrigir aberrações de lente, realizando uma análise estatística mais extensa. A idéia é traçar desvios entre a localização esperada e a real das estrelas na imagem. Isso é semelhante a "distorcer" ou " georreferenciar " dados do mapa. Como uma solução rápida e suja, você pode pressionar o GIS ou o software de processamento de imagem (como o ENVI ) em serviço para georreferenciar (ou astroreferenciar) qualquer imagem. Esse software geralmente não realiza estimativas de ML de transformações projetivas, mas pode fazer aproximações polinomiais de alta ordem. Uma transformação polinomial de ordem 2 ou ordem 3 pode fazer um trabalho suficientemente bom, dependendo do seu aplicativo.

• É possível melhorar a precisão combinando várias imagens dos mesmos objetos.

Fazer isso com o mesmo grau de precisão que os astrônomos profissionais fazem seria realmente difícil. Isso exigiria uma caracterização extremamente precisa das distorções produzidas pela lente e das imperfeições no sensor da câmera. No entanto, você provavelmente não precisa desse grau de precisão. Deveria ser suficiente para você supor que sua lente não apresenta grandes quantidades de distorção (o que é uma boa suposição para uma lente de qualidade) e que o sensor de sua câmera está muito próximo de uma grade perfeitamente regular (que é uma suposição muito boa para até uma câmera barata).

Tudo o que resta é descobrir a transformação de coordenadas que descreve a orientação da câmera, ou seja, a direção em que ela foi apontada e o grau em que foi girada.

O que você está procurando então é chamado de transformação afim ou mapa afim. Que é apenas um nome elegante para uma matriz pela qual você multiplicaria suas coordenadas de pixel para obter suas coordenadas astronômicas. No caso de um mapa afim, essa transformação pode incluir qualquer grau de rotação, escala, cisalhamento e translação.

O significado do componente de rotação é bastante óbvio. O fator de escala descreve simplesmente quanto do céu é coberto por cada pixel em termos de RA / Dec. O cisalhamento é uma transformação que faria a imagem de um retângulo se tornar um paralelogramo, mas não deveria haver esse efeito em uma imagem de objetos no infinito (como estrelas). Por fim, o componente de tradução simples adiciona um deslocamento para explicar o fato de que o pixel (x = 0, y = 0) na sua imagem provavelmente não corresponde a (RA = 0, dez = 0).

Como você tem três estrelas de referência em sua imagem, você tem informações suficientes para calcular o relacionamento entre as coordenadas de pixel e o RA / Dec que você está procurando. Isso seria feito pelo ajuste de mínimos quadrados lineares (não mínimos quadrados não lineares, como mencionado acima) para determinar os valores dos componentes da matriz que melhor correspondem às suas coordenadas de pixel com o RA / Dec conhecido das estrelas de referência. Uma vez que a matriz é estabelecida, você pode aplicá-la às coordenadas de pixel de outras estrelas para obter seu RA / Dec.

Embora eu pudesse fazer isso com relativa facilidade, infelizmente não tenho certeza de como ajudá-lo. Isso envolveria alguma habilidade matemática que está um pouco além do escopo do photo.SE. Sou engenheiro óptico, mas não sou muito fotógrafo; o software que eu usaria para isso foi projetado para os engenheiros fazerem computação numérica pesada e não é realmente uma ferramenta fotográfica. Pode haver maneiras de fazer isso usando pacotes de software voltados para fotógrafos, mas eu não os conheço.