Python: encuentre el píxel más brillante en el archivo de ajustes e ignore los artefactos más brillantes

Tengo 4000 imágenes de archivos de ajuste de SPITZER, cada una de las cuales contiene la misma estrella durante un período de tiempo. Cada imagen es un cubo de fotometría de 32 píxeles * 32 píxeles * 64 marcos de tiempo consecuentes. Necesito mirar a través de todos los cuadros individuales y ordenar el píxel por brillo. No es difícil de hacer con Python, pero a veces tengo en los marcos eventos/artefactos de corta duración que son más brillantes que el píxel más brillante de la estrella (ver más abajo). Estos artefactos duran solo uno o dos marcos de tiempo consecutivos.

¿Cómo puedo corregir fácilmente los artefactos? ¿Existe un paquete astropy/pyfits que haga esto? ¿O algún tipo de marcación de píxeles defectuosos en SPITZER (sé que las imágenes de Kepler tienen marcas de píxeles defectuosos)?

Opcional: Sería conveniente tener el código en forma vectorizada, evitando loops, así ahorro tiempo de cómputo


Fig: En el cuadro de la izquierda todo es normal, en el cuadro de la derecha, en la esquina superior derecha se ve uno de los artefactos, que tiene un valor por debajo de un umbral de saturación

En el panel de la izquierda todo es normal, en el panel de la derecha, en la esquina superior derecha, se ve uno de los artefactos, que tiene un valor por debajo de un umbral de saturación

La función de dispersión de puntos de la estrella debe ser consistente en todo el conjunto de datos, y es probable que los artefactos tengan PSF diferentes a los de la estrella. En lugar de simplemente seleccionar el píxel más brillante en cada cuadro (que es lo que supongo de la descripción), identifique la estrella a partir de su PSF y luego seleccione los valores máximos dentro de ese dominio.
Parece que está pidiendo un kernel 3-D, es decir, x, y, tiempo, para aplicar en su búsqueda.
@Mick Sí, hasta ahora estoy buscando el píxel más brillante en los marcos. Y tienes razón, será mejor buscar el PSF, pero no veo una forma fácil de codificar esto. ¿Hay alguna función en astropy o pyfits que reconozca el PSF de la estrella y me dé las coordenadas del centro PSF?
@Nestak Bueno... buscar en Google "python astropy psf" me ha demostrado que hay una función de coincidencia de PSF en el paquete photutils de astropy.
@Mick Gracias, buscar en Google resolvió el problema, de hecho :) Basado en su sugerencia, escribí mi propia respuesta a la pregunta
¡Fresco! Me alegro de que hayas encontrado una solución, +1 por compartirla.

Respuestas (2)

Siguiendo el comentario de Mick, encontré una función astropía, DAOStarFinder, que hace lo que estaba buscando. Escanea una imagen en busca de fuentes de luz y se puede establecer un umbral para una detección, por ejemplo, la "redondez", la "nitidez" o el brillo de un objeto deben estar entre [-0.05,0.05], [0.6,0.8] o >5 *desviación_estándar, respectivamente. Debo decir que tengo la impresión de que las medidas de redondez y nitidez no funcionan de manera consistente y, a veces, los valores para el mismo objeto cambian enormemente de un cuadro a otro. Pero DAOStarFinder también le brinda una estimación aproximada de la densidad de flujo del objeto, que podría usar para clasificar las fuentes y filtrar la relevante. Aquí hay un enlace a un script de ejemplo ya la documentación de DAOStarFinder.

>>> from astropy.stats import sigma_clipped_stats
>>> from photutils import datasets
>>> from photutils import DAOStarFinder

>>> hdu = datasets.load_star_image()    
>>> data = hdu.data[0:400, 0:400]    
>>> mean, median, std = sigma_clipped_stats(data, sigma=3.0, iters=5)
>>> daofind = DAOStarFinder(fwhm=3.0, threshold=5.*std)    
>>> sources = daofind(data - median)    

>>> print(sources)

Además de la solución sugerida en los comentarios (ajuste PSF), que probablemente sea la forma adecuada de realizar análisis adicionales, se me ocurren varias soluciones:

  1. Verifique si los artefactos son de hecho píxeles saturados (a veces se resta el sesgo, por lo que es posible que deba verificar el valor de saturación menos el sesgo). Ignora cualquier píxel con ese valor.

Si eso no funciona:

  1. Obtenga el valor medio de cada píxel a lo largo de la serie temporal, rechace todos los píxeles que sean más de X (%) lejos de él.

  2. Para cada píxel de la serie temporal, compruebe si su brillo cambia en más de X (%). Tenga en cuenta el valor anterior, rechace siempre que los valores posteriores también estén por encima del umbral.

Todos estos métodos tendrán problemas si el píxel brillante está en la posición de la estrella.

Estaba pensando en hacer lo que sugieres, implementar algunas condiciones y cláusulas if, que mi código python tiene que verificar. El problema sería que esto ralentizará mi código, procesando 260 000 puntos de datos, que ya se ejecutan unas horas en un servidor computacional. Así que encontré una función astronómica, DAOStarFinder, que hace el trabajo.