Resolución espacial en Observaciones astronómicas: Ajuste Gaussiano y detección automática de objetos celestes¶
Siria Sadeddin¶
Introducción: Resolución espacial¶
En observaciones astronómicas e imágenes en general, llamamos resolución espacial a la distancia angular minima a la que pueden estar dos fuentes puntuales de luz y aun poder ser reconocidas como objetos individuales. En el caso de la astronomía, este efecto tiene que ver con la dispersión de la luz al atravezar la atmósfera, la cual hace que una estrella, que debería en principio aparecer como una fuente puntual (pues las estrellas están muy lejos), aparezca en cambio como una mancha. Así, si dos estrellas están demasiado cerca sus manchas se superpondrán hasta el punto en que sea imposible distinguirlas como fuentes individuales.
[1]
Para modelar este efecto, típicamente consideramos la acción de la atmósfera como la convolución de la imagen "perfecta" (como se vería desde el espacio) con un kernel gaussiano. El ancho de esa función gaussiana 2D caracteriza las condiciones de observación, varía con las condiciones climáticas y para cada sitio de la Tierra.
La resolución espacial normalmente se toma como la "Anchura a media altura" (FWHM) de la gaussiana caracteristica registrada durante una observación. Es decir, si dos estrellas están a una distancia aparente en el cielo menor que el FWHM del efecto atmosférico, la luz de ambas fuentes se mezclará después de la convolución hasta el punto de impedir reconocerlas de modo individual. Además, la atmósfera puede interactuar de maneras distintas con la luz de distintas longitudes de onda, de manera que el ancho de la gaussiana puede ser distinto para observaciones con diferentes filtros. El objetivo de esta tarea es medir de forma aproximada la resolución espacial en una noche de observación en Zapatoca, Santander (Colombia), a partir de una foto del cielo estrellado.
Anchura a media altura (FWHM)¶
Para una Gaussiana en 1D¶
La Anchura a media altura, abreviada FWHM (del inglés Full Width at Half Maximum) es una medida de la extensión de una función, que viene dada por la diferencia entre los dos valores extremos de la variable independiente en los que la variable dependiente es igual a la mitad de su valor máximo. [2]
La FWHM se aplica a fenómenos tales como la duración de un pulso y la anchura espectral de fuentes utilizadas para comunicaciones ópticas y la resolución de espectrómetros.
Cuando la función considerada es una distribución normal de la forma:
$$ \frac{1}{\sigma\sqrt{2\pi}}\exp{-\frac{(x-x_0)^2}{2\sigma^2}}$$donde $\sigma$ es la desviación típica y $x_{0}$ puede ser cualquier valor (la anchura de la función no cambia con una traslación), la relación entre FWHM y la desviación típica es:
$$ FWHM = 2\sqrt{2\ln(2)}\sigma \approx 2.35482\sigma$$
Para una Gaussiana en 2D¶
En dos dimensiones, el exponente de la potencia de e dentro de la función de Gauss es cualquier valor negativo y definido en forma cuadrática. Como consecuencia, los niveles de la función siempre serán elipses. [3]
Un ejemplo de la función de dos dimensiones es:
$$f(x,y)=A\exp \left(-\left({\frac {(x-x_{0})^{2}}{2\sigma _{X}^{2}}}+{\frac {(y-y_{0})^{2}}{2\sigma _{Y}^{2}}}\right)\right)$$Para este experimento tomaremos el maximo entre $\sigma_{x}$ y $\sigma_{y}$ como resultado final para el cáculo del FWHM.
Exploración de los datos¶
Lectura de la imagen¶
Se leerá la imagen almacenada en la carpeta "data" llamada "zapatocaImage.jpeg" como un array de numpy. Ese array estará compuesto de 3 matrices, una tras otra, correspondiente a los canales R,G,B. Lo que significa:
- R: Rojo
- G: Verde
- B: Azul
# load libraries
import numpy as np
import math
from matplotlib import image
from matplotlib import pyplot
import matplotlib.patches as patches
import pandas as pd
from mpl_toolkits.mplot3d import Axes3D
from scipy import optimize
import seaborn as sns
from scipy.stats import norm
import statistics as st
# load image
img = np.asarray(image.imread('./zapatocaImage.jpeg'))
# summarize shape of the pixel array
print(img.dtype)
print(img.shape)
# display the array of pixels as an image
figure, plot = pyplot.subplots(figsize=(15,10))
plot.imshow(img)
plot.set_axis_off()
pyplot.show()
Pretratamiento de la imagen¶
Conversion de la imagen a blanco y negro¶
Para hacer un estudio primario de la imagen, la convertiremos a una imagen de un solo canal (escala de grises) que contienne para cada nuevo pixel el promedio de valores de los pixeles en la imagen RGB. En principio, podemos hacer cualquier tipo de combinación lineal de los tres canales RGB para obtener la imagen en escala de grises.
#compute the mean of the image pixels in the channel axis
image_gray_m=np.mean(img,axis=2)
pyplot.figure(figsize=(15,10))
# summarize shape of the pixel array
print(image_gray_m.dtype)
print(image_gray_m.shape)
# display the array of pixels as an image
figure, plot = pyplot.subplots(figsize=(15,10))
plot.imshow(image_gray_m,cmap='gray')
plot.set_axis_off()
pyplot.show()
Como el color es una percepción psicológica, se ha encontrado que una mejor aproximación al problema es hacer la combinación lineal:
$$y_{gray}=0.2126R_{lineal}+0.7252G_{lineal}+0.0722B_{lineal}$$image_gray_l=img[:,:,0]*0.2126+img[:,:,1]*0.7252+img[:,:,2]*0.0722
# summarize shape of the pixel array RGB BGR
print(image_gray_l.dtype)
print(image_gray_l.shape)
# display the array of pixels as an image
fig, (ax1, ax2) = pyplot.subplots(1, 2,figsize=(20,40))
ax1.imshow(image_gray_m,cmap='gray')
ax1.set_axis_off()
ax2.imshow(image_gray_l,cmap='gray')
ax2.set_axis_off()
pyplot.show()
Podriamos usar cualquera de los dos resultados, pero me parece que la promedio presenta una imagen mas clara (image a la izquierda)
Extracción de las estrellas en la imagen mediante el método de detección de outliers¶
Valores atípicos o outliers¶
En estadística, tales como muestras estratificadas, un valor atípico (en inglés outlier) es una observación que es numéricamente distante del resto de los datos. En una distribución normal, los valores que se encuentran mas por encima o por debajo de 3 desviaciones estandar representan menos del 0.1% de la población, por esta razón podemos decir que estos valores son extraños para esa distribución.
Para hallar la posición de las estrellas en el cielo utilizaremos un método de deteccion de outliers en imagenes, el cielo es obscuro, por lo tanto las estrellas son outliers en la imagen. lo que se hará es calcular la media y desviación estandar de valores en los pixeles de la imagen, para luego encontrar los valores que se alejen de la media mas de 3 desviaciones estandar, estos valores son las estrellas. Además de esto, la selección de estrellas se hará para la sección de la imagen $y<400$ y $250<x<1000$, esto elimina los objetos que no pertenecen al cielo en obsevacion. Se nota la precencia de un satelite (presuntamente de Elon Musk :D) pasando por el cielo, que no corresponde a las estrellas en estudio (linea punteada brillante en la parte inferior de la imagen).
Como hemos mencionado, se tomará como estrella todos los valores que superen 3 desviaciones estandar por encima de la media (los valores altos en los pixeles son brillantes y los cercanos a cero son obscuros o negros)
$$estrella=mean_{pixels}+3*std_{pixels}$$Seleccionaremos la zona de la foto que solo incluya el cielo nocturno¶
image_gray_sky = image_gray_m[:400,250:1000]
image_sky = img[:400,250:1000]
# display the array of pixels as an image
fig, (ax1, ax2) = pyplot.subplots(1, 2,figsize=(20,40))
ax1.imshow(image_gray_sky,cmap="gray")
ax1.set_axis_off()
ax1.set_title("Cielo en escala de grises",fontsize=16)
ax2.imshow(image_sky,cmap="gray")
ax2.set_axis_off()
ax2.set_title("Cielo con canales RGB",fontsize=16)
pyplot.show()
Demostraremos que los datos del cielo siguen la distribución normal, por ende es natural hacer la detección de outliers como hemos descrito antes¶
h,w=image_gray_sky.shape
fig, (ax1) = pyplot.subplots(1, 1,figsize=(15,5))
sns.set_color_codes()
sns.distplot(image_gray_sky.reshape(h*w),ax=ax1,fit=norm,color="y")
ax1.legend(labels=['Distribución Normal','Distribución real'])
ax1.set_title("Distribucion de valores en los pixeles")
ax1.set_ylabel("frecuencia",fontweight ='bold',fontsize=14)
ax1.set_xlabel("Valor del pixel",fontweight ='bold',fontsize=14)
pyplot.show()
print({"min":np.min(image_gray_sky),"max":np.max(image_gray_sky),"mean":np.mean(image_gray_sky),"std":np.std(image_gray_sky)})
Crearemos una funcion que nos permitirá hayar las estrellas mediante detección de outliers¶
# searching outliers
def star_search(image):
'''
Esta función identifica las estrellas en el
cielo nocturno identificando valores atípicos
en los pixeles de la imagen
Input:
image: array en 2D
output:
list: posiciones de los valores atípicos encontrados
'''
mean=np.mean(image)
std=np.std(image)
h,w=image.shape
stars=[]
for i in range(h):
for j in range(w):
if (image[i,j]>=mean+10*std) and (i+15<h) and (j+15<w) and (i-15>0) and (j-15>0): ## dont take stars near the borders
stars.append([i,j])
return stars
stars=star_search(image_gray_sky)
stars[0:5]
Seleccionaremos 8 estrellas de las halladas por nuestro algoritmo con un corte de $ 3 \times 3$ pixeles para evaluar el resultado.¶
#corte de la seccion donde se encuentra la estrella
fig, axs = pyplot.subplots(2, 4, figsize=(15,10))
for i in range(2):
for j in range(4):
image_unit=image_gray_sky[stars[i+j][0]-6:stars[i+j][0]+6,stars[i+j][1]-6:stars[i+j][1]+6]
axs[i, j].imshow(image_unit,cmap='gray')
axs[i, j].set_axis_off()
pyplot.show()
Vamos a marcar el centroide de cada estrella solo por curiosidad¶
#corte de la seccion donde se encuentra la estrella
def star_centroid(image,stars):
fig, axs = pyplot.subplots(2, 4, figsize=(15,10))
for i in range(2):
for j in range(4):
image_unit=image[stars[i+j][0]-6:stars[i+j][0]+6,stars[i+j][1]-6:stars[i+j][1]+6]
center_x=math.ceil(np.sum(np.sum(image_unit,axis=1)*range(12))/np.sum(np.sum(image_unit,axis=1)))
center_y=math.ceil(np.sum(np.sum(image_unit,axis=0)*range(12))/np.sum(np.sum(image_unit,axis=0)))
axs[i, j].plot(center_x,center_y,c='r',marker="+", markersize=40)
axs[i, j].imshow(image_unit,cmap='gray')
axs[i, j].set_axis_off()
pyplot.show()
star_centroid(image_gray_sky,stars)
Ahora veamos el resultado marcando todas las estrellas halladas rodeandolas con una caja de color azul en la imagen original¶
#dentro de la caja roja esta la estrella mas brillante detectada
def star_boxing(image,stars,color):
# Create figure and axes
fig, ax = pyplot.subplots(figsize=(14,8))
# Display the image
ax.imshow(image)
for i in range(len(stars)):
# Create a Rectangle patch
rect = patches.Rectangle((stars[i][1]-5,stars[i][0]-5), 10, 10, linewidth=1, edgecolor=color, facecolor='none')
# Add the patch
ax.add_patch(rect)
ax.set_axis_off()
star_boxing(image_sky,stars,'b')
Ajuste Gaussiano de la intensidad luminosa de las estrellas¶
Veremos que la luminosidad de la estrella se distribuye en forma gaussiana, por lo tanto es natural hacer un ajuste gaussiano a los datos de esta.
Vamos a graficar la estrella en un espacio tridimensional.¶
def surface_plot (matrix, **kwargs):
'''
Hace un grafico de superficie a partir de una matriz de valores para el eje z
'''
# acquire the cartesian coordinate matrices from the matrix
# x is cols, y is rows
(x, y) = np.meshgrid(np.arange(matrix.shape[0]), np.arange(matrix.shape[1]))
fig = pyplot.figure(figsize=(10,10))
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(x, y, matrix, **kwargs)
return (fig, ax, surf)
# m.shape must be (6,6)
def star3D(image,title):
'''
Toma la imagen recortada de una estrella y
hace un grafico de superficie para la intensidad de luz
'''
m = np.fromfunction(lambda x, y: image, (6, 6))
(fig, ax, surf) = surface_plot(m, cmap = pyplot.cm.coolwarm)
pyplot.title(title)
pyplot.show()
# estrella seleccionada
i=0
#corte de la seccion donde se encuentra la estrella
image_unit=image_gray_sky[stars[i][0]-8:stars[i][0]+8,stars[i][1]-8:stars[i][1]+8]
#centroide de la estrella
center_x=math.ceil(np.sum(np.sum(image_unit,axis=1)*range(16))/np.sum(np.sum(image_unit,axis=1)))
center_y=math.ceil(np.sum(np.sum(image_unit,axis=0)*range(16))/np.sum(np.sum(image_unit,axis=0)))
#recorte alrededor del centroide
image_center = image_unit[center_y-3:center_y+3,center_x-3:center_x+3]
# grafico de la superficie recortada
star3D(image_center,"Representación en 3D de la intesidad luminosa de una estrella")
Procedemos entonces a hacer un ajuste gaussiano de los datos para una estrella¶
# calcula los parametros de la funcion gaussiana de acuerdo a los datos entregados
def moments(data):
"""Returns (height, x, y, width_x, width_y)
the gaussian parameters of a 2D distribution by calculating its
moments """
total = data.sum()
X, Y = np.indices(data.shape)
##gaussian mean position
x = (X*data).sum()/total
y = (Y*data).sum()/total
col = data[:, int(y)]
##std x
width_x = np.sqrt(np.abs((np.arange(col.size)-x)**2*col).sum()/col.sum())
row = data[int(x), :]
##std y
width_y = np.sqrt(np.abs((np.arange(row.size)-y)**2*row).sum()/row.sum())
height = data.max()
mean=np.min(image_gray_sky)
return height, x, y, width_x, width_y, mean
# con los parametros calculados por la funcion "moments"
# se predicen los datos de una funcion gaussiana que tenga esos parametros
def gaussian(height, center_x, center_y, width_x, width_y,mean):
"""Returns a gaussian function with the given parameters"""
width_x = float(width_x)
width_y = float(width_y)
return lambda x,y: mean + height*np.exp(-(((center_x-x)/width_x)**2+((center_y-y)/width_y)**2)/2)
# se optimizan los parametros por minimos cuadrados
def fitgaussian(data):
"""Returns (height, x, y, width_x, width_y)
los parámetros optimizados por minimos cuadrados"""
params = moments(data)
#la funcion errorfunction se aplicará sobre los parametros p
#se calcula la diferencia de los valores predichos y los reales
errorfunction = lambda p: np.ravel(gaussian(*p)(*np.indices(data.shape)) - data)
#se optimiza la funcion de error, hallando los parametros que
# proveen el minimo error cuadrático
p, success = optimize.leastsq(errorfunction, params)
return p
Graficaremos el resultado de la optimización¶
# Create the gaussian data
Xin, Yin = image_center.shape
data = image_center
params = fitgaussian(data)
fit = gaussian(*params)
pyplot.matshow(data, cmap="Reds")
pyplot.contour(fit(*np.indices(data.shape)), cmap=pyplot.cm.copper)
ax = pyplot.gca()
(height, x, y, width_x, width_y,mean) = params
pyplot.text(0.95, 0.05, """
height : %.1f
x : %.1f
y : %.1f
width_x : %.1f
width_y : %.1f"""
%(height,x, y, width_x, width_y),
fontsize=12, horizontalalignment='right',
verticalalignment='bottom', transform=ax.transAxes)
pyplot.show()
Observemos graficamente el proceso de ajuste gaussiano paso a paso¶
star3D(fit(*np.indices(data.shape)),"Gaussian Fit")
star3D(image_center,"Cropped data")
star3D(image_unit,"Original Data")
Ajuste gaussiano para todas las estrellas¶
Vamos a hacer el ajuste gaussiano para todas las estrellas detectadas, con los parametros calculados podemos generar estadísticas interesantes
def star_parameters(image,stars):
#initialize parameters list
parameters=[]
## iteration over stars
for star in stars:
image_unit = image[star[0]-5:star[0]+5,star[1]-5:star[1]+5]
#centroide de la estrella
center_x=math.ceil(np.sum(np.sum(image_unit,axis=1)*range(10))/np.sum(np.sum(image_unit,axis=1)))
center_y=math.ceil(np.sum(np.sum(image_unit,axis=0)*range(10))/np.sum(np.sum(image_unit,axis=0)))
#recorte alrededor del centroide
image_center = image_unit[center_y-3:center_y+3,center_x-3:center_x+3]
## gaussian fit parameters
param = fitgaussian(image_center)
parameters.append(param)
parameters_df = pd.DataFrame(parameters,columns=["height", "mean_x", "mean_y", "std_x", "std_y","mean"])
parameters_df["FWHM_x"] = 2.35*parameters_df["std_x"]
parameters_df["FWHM_y"] = 2.35*parameters_df["std_y"]
return parameters_df
parameters_df=star_parameters(image_gray_sky,stars)
Obtenemos un data set de parametros para el estudio¶
parameters_df.head()
Histograma de FWHM¶
fig, (ax1,ax2) = pyplot.subplots(1, 2,figsize=(15,5))
sns.set_color_codes()
sns.distplot(parameters_df["FWHM_x"],ax=ax1,color="y")
ax1.legend(labels=['Distribución FWHM x','Distribución real'])
ax1.set_title("Distribucion de valores FWHM")
ax1.set_ylabel("frecuencia",fontweight ='bold',fontsize=14)
ax1.set_xlabel("FWHM",fontweight ='bold',fontsize=14)
sns.distplot(parameters_df["FWHM_y"],ax=ax2,color="y")
ax2.legend(labels=['Distribución FWHM y','Distribución real'])
ax2.set_title("Distribucion de valores FWHM")
ax2.set_ylabel("frecuencia",fontweight ='bold',fontsize=14)
ax2.set_xlabel("FWHM",fontweight ='bold',fontsize=14)
pyplot.show()
Minimo y máximo de FWHM¶
print({"MIN_FWHM_x":np.min(parameters_df["FWHM_x"]),"MIN_FWHM_y":np.min(parameters_df["FWHM_y"]),"MAX_FWHM_x":np.max(parameters_df["FWHM_x"]),"MAX_FWHM_y":np.max(parameters_df["FWHM_y"])})
Media de FWHM¶
print({"MEAN_FWHM_x":np.mean(parameters_df["FWHM_x"]),"MEAN_FWHM_y":np.mean(parameters_df["FWHM_y"])})
Mediana de FHWM¶
print({"MEDIAN_FWHM_x":np.median(parameters_df["FWHM_x"]),"MEDIAN_FWHM_y":np.median(parameters_df["FWHM_y"])})
Incertidumbre¶
print({"INCERTIDUMBRE_FWHM_x":np.std(parameters_df["FWHM_x"])/math.sqrt(len(parameters_df["FWHM_x"])),"INCERTIDUMBRE_FWHM_y":np.std(parameters_df["FWHM_y"])/math.sqrt(len(parameters_df["FWHM_y"]))})
Tomaremos el promedi de los resultados en MEAN_FWHM y diremos que: $$FWHM\approx 2.4 \pm 0.1 pixeles$$
Es interesante observar un grupo de estrellas que se alejan de la media en el histograma, exploremos la razon de esto
weird_stars=[stars[i] for i in parameters_df[(parameters_df.FWHM_x>4) | (parameters_df.FWHM_y>4)].index.values]
star_boxing(image_sky,weird_stars,'b')
weird_stars=[[147,143],[315,376],[376,201]]
fig, axs = pyplot.subplots(2, 2, figsize=(15,10))
for i in range(2):
for j in range(2):
image_unit=image_gray_sky[weird_stars[i+j][0]-6:weird_stars[i+j][0]+6,weird_stars[i+j][1]-6:weird_stars[i+j][1]+6]
axs[i, j].imshow(image_unit,cmap='gray')
axs[i, j].set_axis_off()
pyplot.show()
Notamos 2 de estas 3 estrellas se encuentran en la constelación de Orion, y las identificamos como Betelgeuse y Rigel, dos de las estrellas mas brillantes del cielo.
Por ende esta rareza en el histograma se debe a lo brillante de estas estrellas, lo que haria mas dificil distinguir una estrella cercana a ellas.
Ajuste gaussiano para cada canal RGB¶
Dado que ya tenemos funciones que permiten hacer los calculos de el ajuste gaussiano y detección de estrellas rápidamente, procederemos a realizar el mismo ejercicio, pero esta vez para cada canal RGB de la imagen
figure, plots = pyplot.subplots(ncols=3, nrows=1,figsize=(15,10))
for i, subplot in zip(range(3), plots):
temp = np.zeros(image_sky.shape, dtype='uint8')
temp[:,:,i] = image_sky[:,:,i]
subplot.imshow(temp)
subplot.set_axis_off()
pyplot.show()
color=['b','r','g']
for i in range(3):
temp = np.zeros(image_sky.shape, dtype='uint8')
temp[:,:,i] = image_sky[:,:,i]
stars = star_search(temp[:,:,i])
pyplot.axis=None
star_boxing(temp,stars,color[i])
color=['RED','GREEN','BLUE']
for i in range(3):
class bcolors:
BLUE = '\033[94m'
GREEN = '\033[92m'
RED = '\033[91m'
ENDC = '\033[0m'
temp = np.zeros(image_sky.shape, dtype='uint8')
temp[:,:,i] = image_sky[:,:,i]
#busca las estrellas
stars = star_search(temp[:,:,i])
# calcula parametros
parameters_df = star_parameters(temp[:,:,i],stars)
parameters_df=parameters_df[(parameters_df["height"]<255)]
parameters_df=parameters_df[(parameters_df["FWHM_x"]<np.mean(parameters_df["FWHM_x"])+2.5*np.std(parameters_df["FWHM_x"])) & (parameters_df["FWHM_y"]<np.median(parameters_df["FWHM_y"])+2.5*np.std(parameters_df["FWHM_y"]))]
#muestra resultados
print(bcolors.RED +"FWHM FOR CHANNEL "+color[i] + bcolors.ENDC)
print({"MEDIAN_FWHM_x":np.median(parameters_df["FWHM_x"]),"MEDIAN_FWHM_y":np.median(parameters_df["FWHM_y"])})
print({"INCERTIDUMBRE_FWHM_x":np.std(parameters_df["FWHM_x"])/math.sqrt(len(parameters_df["FWHM_x"])),"INCERTIDUMBRE_FWHM_y":np.std(parameters_df["FWHM_y"])/math.sqrt(len(parameters_df["FWHM_y"]))})
Encontramos que el mejor canal es el azul con $$FWHM = 2.3 \pm 0.1 $$
Resultados¶
Se encuentra que el peor canal para hacer observaciones de estrellas es el rojo por su mayor FWHM. El mejor canal sería el azul, con menor FWHM. En el caso de la imagen en escala de grises la distancia mínima de los centroides para distinguir dos estrellas es de 2 pixeles.