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.

picture [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$$

image

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
In [4]:
# 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
In [5]:
# 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()
uint8
(789, 1184, 3)

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.

In [6]:
#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()
float64
(789, 1184)
<Figure size 1080x720 with 0 Axes>

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}$$
In [7]:
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()
float64
(789, 1184)

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.

image

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

In [8]:
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

In [9]:
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()
/usr/local/lib/python3.7/dist-packages/seaborn/distributions.py:2557: FutureWarning: `distplot` is a deprecated function and will be removed in a future version. Please adapt your code to use either `displot` (a figure-level function with similar flexibility) or `histplot` (an axes-level function for histograms).
  warnings.warn(msg, FutureWarning)
In [10]:
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)})
{'min': 46.0, 'max': 255.0, 'mean': 84.9697677777778, 'std': 13.678023265705608}

Crearemos una funcion que nos permitirá hayar las estrellas mediante detección de outliers

In [79]:
# 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]
Out[79]:
[[20, 359], [34, 83], [42, 474], [46, 635], [56, 205]]

Seleccionaremos 8 estrellas de las halladas por nuestro algoritmo con un corte de $ 3 \times 3$ pixeles para evaluar el resultado.

In [80]:
#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

In [81]:
#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

In [82]:
#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.

In [83]:
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()
In [84]:
# 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

In [85]:
# 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

In [86]:
# 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

In [87]:
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

In [108]:
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

In [110]:
parameters_df.head()
Out[110]:
height mean_x mean_y std_x std_y mean FWHM_x FWHM_y
0 181.334251 2.748455 2.640425 0.893804 0.693985 71.939352 2.100439 1.630865
1 172.334018 3.263542 2.938102 0.971330 0.752327 82.773592 2.282625 1.767968
2 171.822397 2.263954 2.941747 0.957741 0.737776 84.416809 2.250691 1.733774
3 136.628336 3.268683 3.139522 1.376193 0.862820 87.926760 3.234053 2.027627
4 202.479220 3.236553 3.675587 1.039083 0.854364 73.956121 2.441846 2.007756

Histograma de FWHM

In [111]:
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()
/usr/local/lib/python3.7/dist-packages/seaborn/distributions.py:2557: FutureWarning: `distplot` is a deprecated function and will be removed in a future version. Please adapt your code to use either `displot` (a figure-level function with similar flexibility) or `histplot` (an axes-level function for histograms).
  warnings.warn(msg, FutureWarning)
/usr/local/lib/python3.7/dist-packages/seaborn/distributions.py:2557: FutureWarning: `distplot` is a deprecated function and will be removed in a future version. Please adapt your code to use either `displot` (a figure-level function with similar flexibility) or `histplot` (an axes-level function for histograms).
  warnings.warn(msg, FutureWarning)

Minimo y máximo de FWHM

In [112]:
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"])})
{'MIN_FWHM_x': -0.14052292418427964, 'MIN_FWHM_y': -1.901729605128938, 'MAX_FWHM_x': 8.408491641993736, 'MAX_FWHM_y': 7.458131424313182}

Media de FWHM

In [113]:
print({"MEAN_FWHM_x":np.mean(parameters_df["FWHM_x"]),"MEAN_FWHM_y":np.mean(parameters_df["FWHM_y"])})
{'MEAN_FWHM_x': 3.095229626171772, 'MEAN_FWHM_y': 2.6116064423565195}

Mediana de FHWM

In [114]:
print({"MEDIAN_FWHM_x":np.median(parameters_df["FWHM_x"]),"MEDIAN_FWHM_y":np.median(parameters_df["FWHM_y"])})
{'MEDIAN_FWHM_x': 2.6158406674952857, 'MEDIAN_FWHM_y': 2.175254638094052}

Incertidumbre

In [115]:
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"]))})
{'INCERTIDUMBRE_FWHM_x': 0.09808421430866827, 'INCERTIDUMBRE_FWHM_y': 0.10376113432660834}

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

In [151]:
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')
In [153]:
weird_stars=[[147,143],[315,376],[376,201]]
In [154]:
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. image

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

In [98]:
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()
In [99]:
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])
In [100]:
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"]))})
FWHM FOR CHANNEL RED
{'MEDIAN_FWHM_x': 2.6991353415041006, 'MEDIAN_FWHM_y': 2.3302362332216235}
{'INCERTIDUMBRE_FWHM_x': 0.1385692010241145, 'INCERTIDUMBRE_FWHM_y': 0.13066209213728613}
FWHM FOR CHANNEL GREEN
{'MEDIAN_FWHM_x': 2.5587965710667557, 'MEDIAN_FWHM_y': 2.1251087172806953}
{'INCERTIDUMBRE_FWHM_x': 0.07221395934376766, 'INCERTIDUMBRE_FWHM_y': 0.0836238877251871}
FWHM FOR CHANNEL BLUE
{'MEDIAN_FWHM_x': 2.5142551846111694, 'MEDIAN_FWHM_y': 2.1050342234253527}
{'INCERTIDUMBRE_FWHM_x': 0.06972255049213588, 'INCERTIDUMBRE_FWHM_y': 0.06777412878176484}

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.

Referencias