Tatiana Acero Cuellar

Universidad Nacional de Colombia

Ejercicios Semana 5

Contexto: 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 (Ver imágenes en este link)

Para modelar este efecto, típicamente consideramos la acción de la atmósferacomo 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 el 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.

Ejercicio: Pasos

* Leer la imagen almacenada en la carpeta data como un array de numpy. Ese array estará compuesto de 3 matrices, una tras otra, correspondiente a los canales R,G,B * Combinar los 3 array para generar una versión blanco y negro de la imagen, en la cual ella consiste de una sola matriz 2D. Puede usar su intuición y prueba y error para combinar las 3 imágenes, explicando el procedimiento elegido. Esto será más interesante que usar un comando desconocido de una biblioteca sofisticada que haga las cosas como una caja negra (queremos practicar numpy) * Recorte un sector de la imagen conteniendo una estrella individual y ajuste una gaussiana 2D simétrica a la imagen de la estrella por mínimos cuadrados, incluyendo una constante aditiva (el cielo "vacio" brilla) * Repita este procedimiento para varias estrellas y presente alguna estadística sobre las medidas de la FWHM de las distintas gaussianas: histograma, media, mediana, desviación estándar * Repita el mismo ejercicio sobre cada una de las bandas R,G,B separadamente y comente si observa diferencias en los resultados

import numpy as np
import imageio as iio
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as patches

from scipy.optimize import leastsq 
# %matplotlib notebook

Importar la imagen almacenada como array de numpy.

La imagen se lee con imageio y se convierte en array con np.asarray

img = np.asarray(iio.imread('data/zapatocaImage.jpeg'))

# -- visualizar la imagen
fig, ax = plt.subplots(figsize=(20, 7))
ax.imshow(img)
plt.show()

# -- shape de la imagen
# -- el ultimo numero indica el numero de canales, al ser una imagen hace referencia a RedGreenBlue
img.shape

(789, 1184, 3)

Combinar los canales

Para combinar los canales, vamos a promediar los tres canales para así "crear" la imagen en escala de grises.

# -- promedio de la imagen en axis = 2 (la dimension de los canales), el tercer numero de la tupla
# -- ahora, no hay tercera dimension
gr_img = img.mean(axis=2)
gr_img.shape

(789, 1184)

# -- visualizar la imagen en escala de grises
fig, ax = plt.subplots(figsize=(20, 7))
# -- usar un mapa de color en grises
ax.imshow(gr_img, cmap ='gist_gray')
plt.show()

Las estrellas que voy a estar analizando las enmarco en la imagen a blanco y negro. En total seleccione 10 estrellas

fig, ax = plt.subplots(figsize=(20, 7))
ax.imshow(gr_img, cmap ='gist_gray')
rect1 = patches.Rectangle((889-5, 283-5), 15, 10, linewidth=1, edgecolor='r', facecolor='none')
ax.text(889-8,283-10,'1',fontsize=10, c='orange')
rect2 = patches.Rectangle((946-5, 250-5), 20, 20, linewidth=1, edgecolor='r', facecolor='none')
ax.text(946-8,250-10,'2',fontsize=10, c='orange')
rect3 = patches.Rectangle((617, 308-5), 20, 20, linewidth=1, edgecolor='r', facecolor='none')
ax.text(617,308-10,'3',fontsize=10, c='orange')
rect4 = patches.Rectangle((605, 16-5), 15, 15, linewidth=1, edgecolor='r', facecolor='none')
ax.text(605+18,16,'4',fontsize=10, c='orange')
rect5 = patches.Rectangle((542, 347-5), 15, 15, linewidth=1, edgecolor='r', facecolor='none')
ax.text(542,347-9,'5',fontsize=10, c='orange')
rect6 = patches.Rectangle((305, 295-5), 15, 15, linewidth=1, edgecolor='r', facecolor='none')
ax.text(305,295-9,'6',fontsize=10, c='orange')
rect7 = patches.Rectangle((386, 141-5), 15, 15, linewidth=1, edgecolor='r', facecolor='none')
ax.text(386,141-9,'7',fontsize=10, c='orange')
rect8 = patches.Rectangle((320, 458), 15, 15, linewidth=1, edgecolor='r', facecolor='none')
ax.text(320,458-9,'8',fontsize=10, c='orange')
rect9 = patches.Rectangle((1071-2, 277), 15, 15, linewidth=1, edgecolor='r', facecolor='none')
ax.text(1071-2,277-9,'9',fontsize=10, c='orange')
rect10 = patches.Rectangle((882, 42-9), 15, 15, linewidth=1, edgecolor='r', facecolor='none')
ax.text(882,42-9,'10',fontsize=10, c='orange')

ax.add_patch(rect1)
ax.add_patch(rect2)
ax.add_patch(rect3)
ax.add_patch(rect4)
ax.add_patch(rect5)
ax.add_patch(rect6)
ax.add_patch(rect7)
ax.add_patch(rect8)
ax.add_patch(rect9)
ax.add_patch(rect10)
plt.show()

Ajuste de gaussiana a imagen blanco y negro

Función gaussiana generalizada en dos dimensiones

$f(x,y) = A \exp\left(- \left(a(x - x_o)^2 + 2b(x-x_o)(y-y_o) + c(y-y_o)^2 \right)\right)$

${\displaystyle {\begin{aligned}a&={\frac {\cos ^{2}\theta }{2\sigma _{X}^{2}}}+{\frac {\sin ^{2}\theta }{2\sigma _{Y}^{2}}}\\[4pt]b&=-{\frac {\sin 2\theta }{4\sigma _{X}^{2}}}+{\frac {\sin 2\theta }{4\sigma _{Y}^{2}}}\\[4pt]c&={\frac {\sin ^{2}\theta }{2\sigma _{X}^{2}}}+{\frac {\cos ^{2}\theta }{2\sigma _{Y}^{2}}}\end{aligned}}}$

def function(params, x, y):
    '''
        Define función gaussiana generalizada en dos dimensiones.
    '''
    A = params[0]
    sigma_x = params[1]
    sigma_y = params[2]
    x0 = params[3]
    y0 = params[4]
    theta = params[5]
    K = params[6]
    
    a = np.cos(theta)**2/(2*sigma_x**2) + np.sin(theta)**2/(2*sigma_y**2)
    b = -np.sin(2*theta)/(4*sigma_x**2) + np.sin(2*theta)/(4*sigma_y**2)
    c = np.sin(theta)**2/(2*sigma_x**2) + np.cos(theta)**2/(2*sigma_y**2)

    z = A*np.exp( - (a*(x-x0)**2 + 2*b*(x-x0)*(y-y0) + c*(y-y0)**2)) + K
    return z

def Error(tpl,x,y,z,incer):
    '''
        Calcula el error entre el valor z dado por la imagen y el determinado por el método de mínimos cuadrados
        o cualquier otro. Al ser el resultado dos dimensional, fue necesario convertirlo en unidimensional.
    '''
    # x son las posiciones reportadas
    # y son los valores medidos en las posiciones x
    # tpl es una tupla con los parámetros para calcular el modelo

    zmodel = function(tpl,x,y)
    errors = (zmodel.flatten() - z.flatten()) / np.sqrt(abs(incer))
    return errors

def calculator(star,p0,Error,incer=1):
    '''
        Extrae las coordenadas x & y de la imagen y calcula la mejor gaussiana mediante el método de mínimos
        cuadrados
    '''
    x = np.arange(0,star.shape[0],1)
    y = np.arange(0,star.shape[1],1)
    x, y = np.meshgrid(y,x)
    best,suss = leastsq(Error, p0, args=(x,y,star,incer))
    return best,x,y

# -- listas para guardar los valore del Full Width at Half Maximum
FWHM_x = []
FWHM_y = []

Formúla para determinar el Full width at half maximum

${\mathrm {FWHM}}=2{\sqrt {2\ln 2}}\;\sigma$

def determinar(star,p0, FWHMx = FWHM_x, FWHMy = FWHM_y,incer=1):
    '''
        Determina la mejor gaussiana dado un punto inicial, usando la función calculator.
        Determina el valor del FWHM en dirección x & y, guarda el valor en una lista.
        Este es determinado usando la formula de https://en.wikipedia.org/wiki/Full_width_at_half_maximum
        Genera una gráfica donde se muestra el contorno de la gaussiana determinada encima de la imagen.
    '''
    best_star,x,y = calculator(star, p0, Error,incer)
    print(best_star)
    FWHM_star_x = 2*np.sqrt(2*np.log(2))*best_star[1]
    FWHM_star_y = 2*np.sqrt(2*np.log(2))*best_star[2]
    FWHMx.append(FWHM_star_x)
    FWHMy.append(FWHM_star_y)
    print(FWHM_star_x,FWHM_star_y)
    
    # -- plot el contorno de la gaussiana encontrada encima de la imagen de la estrella
    fig, ax = plt.subplots(figsize=(20, 7))
    ax.imshow(star,cmap='gist_gray')
    ax.contour(x, y, function(best_star,x,y), 3, colors='r')
    plt.show()

Para cada uno de los casos mostrados, los parámetros iniciales p0 son elegidos convenientemente para que la gaussiana determinada por el método sea la mejor posible.

Estrella 1

str_1 = gr_img.copy()
str_1 = str_1[283:288,890:895] 
plt.imshow(str_1, cmap ='gist_gray')
plt.show()

p0 = [70, 1.2, 11, 1, 1, 0, 0]
determinar(str_1,p0)

[186.07969768   0.92039049   1.01048277   2.00654537   2.0493886
  -0.300277    69.41303048]
2.1673539689071366 2.3795050880398305

Mirar como es para una sola dimesión

str_1

array([[ 76.33333333,  72.33333333,  92.        ,  89.        ,
         76.        ],
       [ 75.33333333, 123.33333333, 182.        , 136.        ,
         70.        ],
       [ 84.33333333, 176.33333333, 247.        , 176.        ,
         78.        ],
       [ 82.33333333, 143.33333333, 191.        , 141.        ,
         78.        ],
       [ 84.33333333,  80.33333333,  84.        ,  86.        ,
         83.        ]])

# -- extraer valores de cada columna para formar un vector con cada columna
# -- el tercer vector es el que parte la estrella en dos
wa = list(zip(*str_1))
list(wa[2])

[92.0, 182.0, 247.0, 191.0, 84.0]

# -- grafiquemos esos puntos
# -- no se ve tan linda como la del profe
plt.scatter(np.arange(0,len(wa[2]),1),wa[2],c='green')
plt.ylabel('Intensidad del pixel')
plt.xlabel('No del pixel fila')
plt.show()

def function1(params, x):
    '''
        Define función gaussiana generalizada en una dimension.
    '''
    A = params[0]
    sigma_x = params[1]
    #sigma_y = params[2]
    x0 = params[2]
    #y0 = params[4]
    theta = params[3]
    K = params[4]
    
    a = np.cos(theta)**2/(2*sigma_x**2) #+ np.sin(theta)**2/(2*sigma_y**2)
    #b = -np.sin(2*theta)/(4*sigma_x**2) + np.sin(2*theta)/(4*sigma_y**2)
    #c = np.sin(theta)**2/(2*sigma_x**2) + np.cos(theta)**2/(2*sigma_y**2)

    z = A*np.exp( - (a*(x-x0)**2 )) + K
    return z

def Error1(tpl,x,y):
    '''
        Calcula el error entre el valor y dado por la imagen y el determinado por el método de mínimos cuadrados
        o cualquier otro.
    '''
    # x son las posiciones reportadas
    # y son los valores medidos en las posiciones x
    # tpl es una tupla con los parámetros para calcular el modelo
    ymodel = function1(tpl,x)
    errors = (ymodel - y)
    #print(errors)
    return errors


p0 = [10, 10, 50, 4,0]
x = np.arange(0,len(wa[2]),1)
y = np.array(wa[2])
best,suss = leastsq(Error1, p0, args=(x,y))
plt.scatter(x,function1(best,x),c='orange',label='Fit')
plt.scatter(x,wa[2], c = 'green',label='Original')
plt.legend()

<matplotlib.legend.Legend at 0x7f63d48ae080>

Continuación en dos dimensiones para las demas estrellas

Estrella 2

str_2 = gr_img.copy()
str_2 = str_2[252:257,949:954]
plt.imshow(str_2, cmap='gist_gray')
plt.show()

p0 = [50, 1.3, 2.5, 4, 1, 0.5, 0]

#p0 = [60, 30.4, 11, 0.01, 1, 0, 0]
determinar(str_2,p0)

[199.17204186   1.08795697   1.01652482   2.3507139    2.25572616
  -1.63067875  71.03541102]
2.5619428907503394 2.393733033767761

Estrella 3

str_3 = gr_img.copy()
str_3 = str_3[309:324,619:636]
plt.imshow(str_3, cmap='gist_gray')
plt.show()

p0 = [50, 15, 10, 4, 1, 0, 0]
determinar(str_3,p0)

[167.1691       2.54202832   2.43955445   8.26104046   6.7936942
  -1.57329251  97.72363193]
5.986019246498141 5.744711731315182

Estrella 4

str_4 = gr_img.copy()
str_4 = str_4[18:23,607:612]
plt.imshow(str_4, cmap='gist_gray')
plt.show()

p0 = [10, 3, 6, 0.9, 1.4, 0, 0]
determinar(str_4,p0)

[171.67801287   0.72367521   0.95397152   1.9468729    2.26661796
  -3.22904682  85.47495094]
1.704124894388063 2.2464312597874616

Estrella 5

str_5 = gr_img.copy()
str_5 = str_5[347:356,542:552]
plt.imshow(str_5, cmap='gist_gray')
plt.show()

p0 = [10, 3.4, 6, 0.8, 1.4, 0, 0]
determinar(str_5,p0)

[164.57072408   1.12967336   1.30573147   4.56858863   4.23103384
   6.39463879 115.86139744]
2.660177475704662 3.0747626286144585

Estrella 6

str_6 = gr_img.copy()
str_6 = str_6[295:300,307:312]
plt.imshow(str_6, cmap='gist_gray')
plt.show()

p0 = [10, 3.3, 5.9, 0.78, 1.7, 0, 0]
determinar(str_6,p0)

[154.74939288   1.06584842   0.86557986   2.09066242   1.67091097
  -1.79900438 106.39500057]
2.509881229793308 2.038284809401981

Estrella 7

str_7 = gr_img.copy()
str_7 = str_7[142:154,387:399]
plt.imshow(str_7, cmap='gist_gray')
plt.show()

p0 = [10, 2, 5, 3, 5, 0, 0]
determinar(str_7,p0)

[140.94787706   1.79851031   2.03013514   6.17006805   5.63716518
   3.68873831 108.41886214]
4.235168131706809 4.780602919231583

Estrella 8

str_8 = gr_img.copy()
str_8 = str_8[459:465,322:327]
plt.imshow(str_8, cmap='gist_gray')
plt.show()

p0 = [10, 3.43, 6.5, 0.01, 0.5, 0, 0]
determinar(str_8,p0)

[133.76503279   1.29368943   0.97282715   2.12474538   2.37096544
  10.63544949 129.69305459]
3.0464058103768084 2.2908328681370387

Estrella 9

str_9 = gr_img.copy()
str_9 = str_9[277:284,1071:1076]
plt.imshow(str_9, cmap='gist_gray')
plt.show()

p0 = [10, 3, 6, 0.8, 1.4, 0, 0]
determinar(str_9,p0)

[188.16928656   0.91808841   1.25701937   2.10259749   3.14037841
   3.29910297  61.1276509 ]
2.161933002277811 2.96005440206742

Estrella 10

str_10 = gr_img.copy()
str_10 = str_10[43:49,883:889]
plt.imshow(str_10, cmap='gist_gray')
plt.show()

p0 = [700, 3.4, 6, 0.8, 1.4, 0, 0]
determinar(str_10,p0)

[2.12780258e+02 9.98780282e-01 1.13708459e+00 2.28269387e+00
 2.89718094e+00 1.60170734e-03 6.12198320e+01]
2.3519478278094055 2.677629596603565

Los valores para el Full Width at Half Maximum (FHWM)

Considere el valor en dirección x & en dirección y. Saque la desviación estándar, el promedio, la mediana y el histograma

FWHM_x = abs(np.array(FWHM_x))
print(FWHM_x)
stdx = FWHM_x.std()
print(stdx)
meanx = FWHM_x.mean()
print(meanx)
medianx = np.median(FWHM_x)
print(medianx)
plt.hist(FWHM_x.flatten())
plt.title('Valores FHWM en la dirección x para escala de grises')
plt.xlabel('Valor FHWM')
plt.ylabel('No de imagénes con ese valor')
plt.show()

[2.16735397 2.56194289 5.98601925 1.70412489 2.66017748 2.50988123
 4.23516813 3.04640581 2.161933   2.35194783]
1.2020310334918178
2.9384954478212486
2.5359120602718237

FWHM_y = abs(np.array(FWHM_y))
print(FWHM_y)
stdy = FWHM_y.std()
print(stdy)
meany = FWHM_y.mean()
print(meany)
mediany = np.median(FWHM_y)
print(mediany)
plt.hist(FWHM_y.flatten())
plt.title('Valores FHWM en la dirección y para escala de grises')
plt.xlabel('Valor FHWM')
plt.ylabel('No de imagénes con ese valor')
plt.show()

[2.37950509 2.39373303 5.74471173 2.24643126 3.07476263 2.03828481
 4.78060292 2.29083287 2.9600544  2.6776296 ]
1.1631400032813368
3.0586548336966284
2.5356813151856628

Ajuste de gaussiana para imagen: canal rojo

Para extraer el canal rojo, usamos todas las filas, todas las columnas, pero solo el primero de los tres canales
$img[:,:,0]$

FWHM_xr = []
FWHM_yr = []

Estrella 1

str_1r = img[:,:,0].copy()
str_1r = str_1r[283:288,890:895] 
plt.imshow(str_1r, cmap ='Reds')
plt.show()

p0 = [70, 1.2, 11, 1, 1, 0, 0]
determinar(str_1r,p0,FWHM_xr,FWHM_yr)

[178.76060595   0.9511361    1.04253055   2.00150507   2.05170422
  -0.30877778  88.35069654]
2.2397543569080045 2.454971833929507

Estrella 2

str_2r = img[:,:,0].copy()
str_2r = str_2r[252:257,949:955]
plt.imshow(str_2r, cmap='Reds')
plt.show()

p0 = [50, 1.3, 2.5, 4, 1, 0, 0]
determinar(str_2r,p0,FWHM_xr,FWHM_yr)

[190.6122258    1.09953422   0.9819177    2.34511946   2.27447217
  -1.62016404  85.47014377]
2.589205231289726 2.312239480960313

Estrella 3

str_3r = img[:,:,0].copy()
str_3r = str_3r[309:324,619:636]
plt.imshow(str_3r, cmap='Reds')
plt.show()

p0 = [50, 15, 10, 4, 1, 0, 0]
determinar(str_3r,p0,FWHM_xr,FWHM_yr)

[143.67174102   2.32328961   2.16090506   8.28633162   6.8475628
  -1.5779969  114.72814395]
5.470928949222188 5.0885425478166635

Estrella 4

str_4r = img[:,:,0].copy()
str_4r = str_4r[18:23,607:612]
plt.imshow(str_4r, cmap='Reds')
plt.show()

p0 = [10, 3, 6, 0.9, 1.4, 0, 0]
determinar(str_4r,p0,FWHM_xr,FWHM_yr)

[170.51074692   0.72404326   0.94899811   1.94940078   2.28406376
   6.190567   100.11931426]
1.704991580249038 2.2347197813928976

Estrella 5

str_5r = img[:,:,0].copy()
str_5r = str_5r[347:356,542:552]
plt.imshow(str_5r, cmap='Reds')
plt.show()

p0 = [10, 3.4, 6, 0.8, 1.4, 0, 0]
determinar(str_5r,p0,FWHM_xr,FWHM_yr)

[159.16045587   1.07357801   1.27234105   4.5756869    4.21360054
   6.28738116 128.40525837]
2.528083008621747 2.9961342032205103

Estrella 6

str_6r = img[:,:,0].copy()
str_6r = str_6r[295:300,307:312]
plt.imshow(str_6r, cmap='Reds')
plt.show()

p0 = [10, 3.3, 5.9, 0.78, 1.7, 0, 0]
determinar(str_6r,p0,FWHM_xr,FWHM_yr)

[149.35942882   1.0738311    0.87541102   2.08507806   1.68940448
  -1.81607981 124.90198335]
2.528679002803907 2.061435424259867

Estrella 7

str_7r = img[:,:,0].copy()
str_7r = str_7r[142:154,387:399]
plt.imshow(str_7r, cmap='Reds')
plt.show()

p0 = [10, 2, 5, 3, 5, 0, 0]
determinar(str_7r,p0,FWHM_xr,FWHM_yr)

[138.7387436    1.9386357   -2.22391119   6.19193906   5.64616817
   0.47421667 129.67231319]
4.565138195369105 -5.236910658145347

Estrella 8

str_8r = img[:,:,0].copy()
str_8r = str_8r[459:465,322:327]
plt.imshow(str_8r, cmap='Reds')
plt.show()

p0 = [10, 3.43, 6.4, 0.01, 0.5, 0, 0]
determinar(str_8r,p0,FWHM_xr,FWHM_yr)

[117.96132322   0.8918261    1.23916011   2.140237     2.33663264
   5.9503522  150.18634797]
2.1000899865145097 2.9179990671423646

Estrella 9

str_9r = img[:,:,0].copy()
str_9r = str_9r[277:284,1071:1076]
plt.imshow(str_9r, cmap='Reds')
plt.show()

p0 = [10, 3, 6, 0.8, 1.4, 0, 0]
determinar(str_9r,p0,FWHM_xr,FWHM_yr)

[182.9258356    0.94142815  -1.29220082   2.10748647   3.15148273
   3.294021    80.34963419]
2.2168938787115526 -3.042900394987331

Estrella 10

str_10r = img[:,:,0].copy()
str_10r = str_10r[43:49,883:889]
plt.imshow(str_10r, cmap='Reds')
plt.show()

p0 = [700, 3.4, 6, 0.8, 1.4, 0, 0]
determinar(str_10r,p0,FWHM_xr,FWHM_yr)

[2.09299250e+02 1.01842846e+00 1.15336547e+00 2.28048409e+00
 2.89212867e+00 1.67741306e-02 7.77637647e+01]
2.398215753935267 2.7159681216478253

Los valores para el Full Width at Half Maximum (FHWM)

Considere el valor en dirección x & en dirección y. Saque la desviación estándar, el promedio, la mediana y el histograma

FWHM_xr = abs(np.array(FWHM_xr))
print(FWHM_xr)
stdxr = FWHM_xr.std()
print(stdxr)
meanxr = FWHM_xr.mean()
print(meanxr)
medianxr = np.median(FWHM_xr)
print(medianxr)
plt.hist(FWHM_xr.flatten())
plt.title('Valores FHWM en la dirección x para la image roja')
plt.xlabel('Valor FHWM')
plt.ylabel('No de imagénes con ese valor')
plt.show()

[2.23975436 2.58920523 5.47092895 1.70499158 2.52808301 2.528679
 4.5651382  2.10008999 2.21689388 2.39821575]
1.1372820062135576
2.8341979943625044
2.463149381278507

FWHM_yr = abs(np.array(FWHM_yr))
print(FWHM_yr)
stdyr = FWHM_yr.std()
print(stdyr)
meanyr = FWHM_yr.mean()
print(meanyr)
medianyr = np.median(FWHM_yr)
print(medianyr)
plt.hist(FWHM_yr.flatten())
plt.title('Valores FHWM en la dirección y para la image roja')
plt.xlabel('Valor FHWM')
plt.ylabel('No de imagénes con ese valor')
plt.show()

[2.45497183 2.31223948 5.08854255 2.23471978 2.9961342  2.06143542
 5.23691066 2.91799907 3.04290039 2.71596812]
1.0760536200525121
3.1061821513502625
2.816983594395095

Ajuste de gaussiana para imagen: canal verde

Para extraer el canal verde, usamos todas las filas, todas las columnas, pero solo el segundo de los tres canales
$img[:,:,1]$

FWHM_xg = []
FWHM_yg = []

Estrella 1

str_1g = img[:,:,1].copy()
str_1g = str_1g[283:288,890:895] 
plt.imshow(str_1g, cmap ='Greens')
plt.show()

p0 = [70, 1.5, 11, 1, 1, 0, 0]
determinar(str_1g,p0,FWHM_xg,FWHM_yg)

[189.66980515   0.90810494   0.99594246   2.00139712   2.04842509
  -0.31654792  67.1450173 ]
2.138423711514054 2.3452652626833705

Estrella 2

str_2g = img[:,:,1].copy()
str_2g = str_2g[252:257,949:955]
plt.imshow(str_2g, cmap='Greens')
plt.show()

p0 = [50, 1.3, 2.5, 4, 1, 0, 0]
determinar(str_2g,p0,FWHM_xg,FWHM_yg)

[200.49323003   1.0742782    1.00605482   2.34595298   2.25284053
  -1.63474596  71.37063765]
2.5297318424914206 2.3690780599175865

Estrella 3

str_3g = img[:,:,1].copy()
str_3g = str_3g[309:324,619:636]
plt.imshow(str_3g, cmap='Greens')
plt.show()

p0 = [50, 15, 10, 4, 1, 0, 0]
determinar(str_3g,p0,FWHM_xg,FWHM_yg)

[170.3427937    2.50270952   2.42460793   8.20606974   6.75458075
  -1.59836573  96.7865555 ]
5.893430538924626 5.709515350529837

Estrella 4

str_4g = img[:,:,1].copy()
str_4g = str_4g[18:23,607:612]
plt.imshow(str_4g, cmap='Greens')
plt.show()

p0 = [10, 4, 6, 0.9, 1.4, 0, 0]
determinar(str_4g,p0,FWHM_xg,FWHM_yg)

[171.70162698   0.72207629   0.95112594   1.94675696   2.26657667
   9.33559274  84.76461655]
1.7003597200006297 2.2397304322878533

Estrella 5

str_5g = img[:,:,1].copy()
str_5g = str_5g[347:356,542:552]
plt.imshow(str_5g, cmap='Greens')
plt.show()

p0 = [10, 3.4, 6, 0.8, 1.4, 0, 0]
determinar(str_5g,p0,FWHM_xg,FWHM_yg)

[168.26372326   1.10033428   1.29450467   4.5902933    4.22695556
   9.45794583 113.51970684]
2.5910892262579335 3.048325549295308

Estrella 6

str_6g = img[:,:,1].copy()
str_6g = str_6g[295:300,307:312]
plt.imshow(str_6g, cmap='Greens')
plt.show()

p0 = [10, 3.5, 6, 0.78, 1.7, 0, 0]
determinar(str_6g,p0,FWHM_xg,FWHM_yg)

[156.84218534   1.06030104   0.85529269   2.10415204   1.66655508
  -1.79293087 104.3264753 ]
2.4968181499130977 2.0140603781124353

Estrella 7

str_7g = img[:,:,1].copy()
str_7g = str_7g[142:154,387:399]
plt.imshow(str_7g, cmap='Greens')
plt.show()

p0 = [15, 2, 5, 3, 5, 0, 0]
determinar(str_7g,p0,FWHM_xg,FWHM_yg)

[145.65432612   1.76651186   1.99543641   6.18715013   5.63414873
   0.57532504 107.4753585 ]
4.159817530728484 4.6988936558118475

Estrella 8

str_8g = img[:,:,1].copy()
str_8g = str_8g[459:465,322:327]
plt.imshow(str_8g, cmap='Greens')
plt.show()

p0 = [10, 3.43, 6.4, 1, 0.5, 0, 0]
determinar(str_8g,p0,FWHM_xg,FWHM_yg)

[132.62532774   0.95533301   1.27035965   2.15830574   2.3619008
   5.91555374 128.72993202]
2.249637328620205 2.9914683602476924

Estrella 9

str_9g = img[:,:,1].copy()
str_9g = str_9g[277:284,1071:1076]
plt.imshow(str_9g, cmap='Greens')
plt.show()

p0 = [10, 3, 6, 0.8, 1.4, 0, 0]
determinar(str_9g,p0,FWHM_xg,FWHM_yg)

[191.01057596   0.90731452   1.24170476   2.10161186   3.14651353
   3.29961415  61.2009776 ]
2.1365624274379273 2.923991268052027

Estrella 10

str_10g = img[:,:,1].copy()
str_10g = str_10g[43:49,883:889]
plt.imshow(str_10g, cmap='Greens')
plt.show()

p0 = [700, 3.4, 6, 0.8, 1.4, 0, 0]
determinar(str_10g,p0,FWHM_xg,FWHM_yg)

[ 2.16319665e+02  9.99190581e-01  1.12653091e+00  2.27946151e+00
  2.90122310e+00 -1.93991955e-02  6.06937104e+01]
2.3529140096920815 2.6527775769756388

Los valores para el Full Width at Half Maximum (FHWM)

Considere el valor en dirección x & en dirección y. Saque la desviación estándar, el promedio, la mediana y el histograma

FWHM_xg = abs(np.array(FWHM_xg))
print(FWHM_xg)
stdxg = FWHM_xg.std()
print(stdxr)
meanxg = FWHM_xg.mean()
print(meanxg)
medianxg = np.median(FWHM_xg)
print(medianxg)
plt.hist(FWHM_xg.flatten())
plt.title('Valores FHWM en la dirección x para la imagen verde')
plt.xlabel('Valor FHWM')
plt.ylabel('No de imagénes con ese valor')
plt.show()

[2.13842371 2.52973184 5.89343054 1.70035972 2.59108923 2.49681815
 4.15981753 2.24963733 2.13656243 2.35291401]
1.1372820062135576
2.824878448558046
2.4248660798025896

FWHM_yg = abs(np.array(FWHM_yg))
print(FWHM_yg)
stdyg = FWHM_yg.std()
print(stdyg)
meanyg = FWHM_yg.mean()
print(meanyg)
medianyg = np.median(FWHM_yg)
print(medianyg)
plt.hist(FWHM_yg.flatten())
plt.title('Valores FHWM en la dirección y para la imagen verde')
plt.xlabel('Valor FHWM')
plt.ylabel('No de imagénes con ese valor')
plt.show()

[2.34526526 2.36907806 5.70951535 2.23973043 3.04832555 2.01406038
 4.69889366 2.99146836 2.92399127 2.65277758]
1.124040162412368
3.0993105893913597
2.7883844225138326

Ajuste de gaussiana para imagen: canal azul

Para extraer el canal azul, usamos todas las filas, todas las columnas, pero solo el último de los tres canales
$img[:,:,2]$

FWHM_xb = []
FWHM_yb = []

Estrella 1

str_1b = img[:,:,2].copy()
str_1b = str_1b[283:288,890:895] 
plt.imshow(str_1b, cmap ='Blues')
plt.show()

p0 = [70, 1.7, 11, 1, 1, 0, 0]
determinar(str_1b,p0,FWHM_xb,FWHM_yb)

[190.2169832    0.90335067   0.99435099   2.01610059   2.04808311
  -0.27729555  52.66186116]
2.1272282649036094 2.3415176386132233

Estrella 2

str_2b = img[:,:,2].copy()
str_2b = str_2b[252:257,949:955]
plt.imshow(str_2b, cmap='Blues')
plt.show()

p0 = [50, 2.9, 2, 4, 1, 0, 0]
determinar(str_2b,p0,FWHM_xb,FWHM_yb)

[205.9355405    1.07443705   1.03024455   2.35673966   2.24362524
   1.41576639  58.80752854]
2.5301059117924476 2.42604050765156

Estrella 3

str_3b = img[:,:,2].copy()
str_3b = str_3b[309:324,619:636]
plt.imshow(str_3b, cmap='Blues')
plt.show()

p0 = [50, 15, 10, 4, 1, 0, 0]
determinar(str_3b,p0,FWHM_xb,FWHM_yb)

[191.02580433   2.73160061   2.64058181   8.29780833   6.77398083
  -1.59662047  81.16874448]
6.432427878391856 6.218094970738713

Estrella 4

str_4b = img[:,:,2].copy()
str_4b = str_4b[18:23,607:612]
plt.imshow(str_4b, cmap='Blues')
plt.show()

p0 = [9, 3.4, 6, 1.23, 1.46, 0, 0]
determinar(str_4b,p0,FWHM_xb,FWHM_yb)

[172.85830022   0.7251488    0.96189654   1.94448225   2.24924495
  -6.36452675  71.51882385]
1.707594921415306 2.265093259490133

Estrella 5

str_5b = img[:,:,2].copy()
str_5b = str_5b[347:356,542:552]
plt.imshow(str_5b, cmap='Blues')
plt.show()

p0 = [10, 12, 10, 2, 1, 0, 0]
determinar(str_5b,p0,FWHM_xb,FWHM_yb)

[165.9928307    1.21438496   1.37811936   4.53438415   4.24417785
   3.5362753  105.2004878 ]
2.8596580347727287 3.2452230957557635

Estrella 6

str_6b = img[:,:,2].copy()
str_6b = str_6b[295:300,307:312]
plt.imshow(str_6b, cmap='Blues')
plt.show()

p0 = [10, 3.3, 5.9, 0.78, 1.7, 0, 0]
determinar(str_6b,p0,FWHM_xb,FWHM_yb)

[158.14516944   1.06375236   0.86631248   2.08220315   1.65755899
  -1.78824612  89.94146728]
2.5049453718991064 2.0400099880751132

Estrella 7

str_7b = img[:,:,2].copy()
str_7b = str_7b[142:154,387:399]
plt.imshow(str_7b, cmap='Blues')
plt.show()

p0 = [11, 2, 5, 3, 5, 0, 0]
determinar(str_7b,p0,FWHM_xb,FWHM_yb)

[140.2375672    1.69057103   1.87880031   6.13700838   5.63296036
   0.59694216  87.61531461]
3.9809905380035873 4.424236622195258

Estrella 8

str_8b = img[:,:,2].copy()
str_8b = str_8b[459:465,322:327]
plt.imshow(str_8b, cmap='Blues')
plt.show()

p0 = [10, 3.43, 6.4, 0.01, 0.5, 0, 0]
determinar(str_8b,p0,FWHM_xb,FWHM_yb)

[152.46960483   1.37726272   1.08272843   2.08071456   2.39950791
  -1.95926817 107.90558619]
3.2432058673937094 2.5496306060688854

Estrella 9

str_9b = img[:,:,2].copy()
str_9b = str_9b[277:284,1071:1076]
plt.imshow(str_9b, cmap='Blues')
plt.show()

p0 = [13, 3, 6.4, 1, 1.1, 0, 0]
determinar(str_9b,p0,FWHM_xb,FWHM_yb)

[190.8831359    0.90622759   1.23752934   2.09884658   3.12378005
   3.30383272  41.79158734]
2.1340028974084926 2.9141588888311842

Estrella 10

str_10b = img[:,:,2].copy()
str_10b = str_10b[43:49,883:889]
plt.imshow(str_10b, cmap='Blues')
plt.show()

p0 = [50, 2, 4, 1, 1, 0, 0]
determinar(str_10b,p0,FWHM_xb,FWHM_yb)

[2.12869194e+02 9.79245215e-01 1.13164829e+00 2.28753325e+00
 2.89802157e+00 6.97495764e-03 4.51627862e+01]
2.305946260209862 2.6648280839211322

Los valores para el Full Width at Half Maximum (FHWM)

Considere el valor en dirección x & en dirección y. Saque la desviación estándar, el promedio, la mediana y el histograma

FWHM_xb = abs(np.array(FWHM_xb))
print(FWHM_xb)
stdxb = FWHM_xb.std()
print(stdxb)
meanxb = FWHM_xb.mean()
print(meanxb)
medianxb = np.median(FWHM_xb)
print(medianxb)
plt.hist(FWHM_xb.flatten())
plt.title('Valores FHWM en la dirección x para la image azul')
plt.xlabel('Valor FHWM')
plt.ylabel('No de imagénes con ese valor')
plt.show()

[2.12722826 2.53010591 6.43242788 1.70759492 2.85965803 2.50494537
 3.98099054 3.24320587 2.1340029  2.30594626]
1.3018950466931847
2.9826105946190706
2.517525641845777

FWHM_yb = abs(np.array(FWHM_yb))
print(FWHM_yb)
stdyb = FWHM_yb.std()
print(stdyb)
meanyb = FWHM_yb.mean()
print(meanyb)
medianyb = np.median(FWHM_yb)
print(medianyb)
plt.hist(FWHM_yb.flatten())
plt.title('Valores FHWM en la dirección y para la image azul')
plt.xlabel('Valor FHWM')
plt.ylabel('No de imagénes con ese valor')
plt.show()

[2.34151764 2.42604051 6.21809497 2.26509326 3.2452231  2.04000999
 4.42423662 2.54963061 2.91415889 2.66482808]
1.2195655968098804
3.1088833661340964
2.607229344995009

Comparemos los cuatro valores

En dirección x

print(stdx,stdxr,stdxg,stdxb)
print(meanx,meanxr,meanxg,meanxb)
print(medianx,medianxr,medianxg,medianxb)

1.2020310334918178 1.1372820062135576 1.1920847689204435 1.3018950466931847
2.9384954478212486 2.8341979943625044 2.824878448558046 2.9826105946190706
2.5359120602718237 2.463149381278507 2.4248660798025896 2.517525641845777

En dirección y

print(stdy,stdyr,stdyg,stdyb)
print(meany,meanyr,meanyg,meanyb)
print(mediany,medianyr,medianyg,medianyb)

1.1631400032813368 1.0760536200525121 1.124040162412368 1.2195655968098804
3.0586548336966284 3.1061821513502625 3.0993105893913597 3.1088833661340964
2.5356813151856628 2.816983594395095 2.7883844225138326 2.607229344995009

Algunas conclusiones sencillitas

El valor el FWHM es mas o menos el mismo para todas las estrellas, entre 2 y 3, en algunos casos es mayor y en otros menor.

Para la escala de grises y para los canales, debí cambiar el punto inicial p0 para que la gaussiana fuera mas apropiada, e.g. para la Estrella n, debí usar cuatro valores de p0, uno para cada canal y para la gris. También elegí los valores de p0 de tal manera que los valores de sigma_x y sigma_y fueran siempre positivos.

El valor para el FWHM si cambia dependiendo del canal del que se parta para hacer el ajuste. Mirando los valores de la desviación estándar, el promedio y la mediana se ve que hay diferencias, sin embargo, los valores podrian considerarse como iguales para los cuatro casos. Igual teniendo en cuenta que el valor inicial p0 para cada estrella debio ser modificado para cada canal.

Ahora, incluyendo las incertidumbres

Escala de grises

# -- listas para guardar los valore del Full Width at Half Maximum
FWHM_xi = []
FWHM_yi = []

Estrella 1

str_1 = gr_img.copy()
str_1 = str_1[283:288,890:895] 
plt.imshow(str_1, cmap ='gist_gray')
plt.show()

p0 = [71, 1.2, 11, 1, 1, 0, 0]
determinar(str_1,p0,FWHM_xi,FWHM_yi,incer=str_1.flatten())

[191.03306608   0.88416723   0.97198929   2.00370142   2.04447666
  -0.30033141  70.9737103 ]
2.0820547263026454 2.28885986684923

Estrella 2

str_2 = gr_img.copy()
str_2 = str_2[252:257,949:954]
plt.imshow(str_2, cmap='gist_gray')
plt.show()

p0 = [57, 1.3, 2.5, 4, 1, 0.5, 0]

#p0 = [60, 30.4, 11, 0.01, 1, 0, 0]
determinar(str_2,p0,FWHM_xi,FWHM_yi,incer=str_2.flatten())

[200.04931224   1.06771002   1.01254502   2.37072665   2.22891783
  -1.63770305  71.38762112]
2.514264960806871 2.3843613096910805

Estrella 3

str_3 = gr_img.copy()
str_3 = str_3[309:324,619:636]
plt.imshow(str_3, cmap='gist_gray')
plt.show()

p0 = [70, 7, 10, 4, 1, 0, 0]
determinar(str_3,p0,FWHM_xi,FWHM_yi,incer=str_3.flatten())

[165.09149595   2.56740779   2.48379149   8.25953555   6.79491264
   1.52624294  97.1682973 ]
6.045783318082368 5.848881983108898

Estrella 4

str_4 = gr_img.copy()
str_4 = str_4[18:23,607:612]
plt.imshow(str_4, cmap='gist_gray')
plt.show()

p0 = [11, 3, 6, 0.9, 1.4, 0, 0]
determinar(str_4,p0,FWHM_xi,FWHM_yi,incer=str_4.flatten())

[181.39720161   0.73469636   0.98245806   1.94925035   2.25033059
   3.05248272  78.45962331]
1.7300777056892003 2.3135119445798464

Estrella 5

str_5 = gr_img.copy()
str_5 = str_5[347:356,542:552]
plt.imshow(str_5, cmap='gist_gray')
plt.show()

p0 = [12, 3, 6, 0.8, 1.4, 0, 0]
determinar(str_5,p0,FWHM_xi,FWHM_yi,incer=str_5.flatten())

[166.83763898   1.12850905   1.28552333   4.56532092   4.21399213
   3.32990317 115.39600089]
2.657435737167189 3.0271761137992788

Estrella 6

str_6 = gr_img.copy()
str_6 = str_6[295:300,307:312]
plt.imshow(str_6, cmap='gist_gray')
plt.show()

p0 = [70, 2.3, 5.9, 0.78, 1.7, 0, 0]
determinar(str_6,p0,FWHM_xi,FWHM_yi,incer=str_6.flatten())

[155.79025228   0.85378903   1.05568553   2.08202815   1.70316379
  -0.22008333 106.48341521]
2.0105195108028022 2.4859494493339502

Estrella 7

str_7 = gr_img.copy()
str_7 = str_7[142:154,387:399]
plt.imshow(str_7, cmap='gist_gray')
plt.show()

p0 = [10, 2, 5, 3, 5, 0, 0]
determinar(str_7,p0,FWHM_xi,FWHM_yi,incer=str_7.flatten())

[138.57467342   1.84344675   2.08148352   6.17117586   5.62605279
   0.63536992 107.33799354]
4.340985359062775 4.901519126606775

Estrella 8

str_8 = gr_img.copy()
str_8 = str_8[459:465,322:327]
plt.imshow(str_8, cmap='gist_gray')
plt.show()

p0 = [10, 3.43, 6.5, 0.01, 0.5, 0, 0]
determinar(str_8,p0,FWHM_xi,FWHM_yi,incer=str_8.flatten())

[133.55992457   0.98677464   1.29549702   2.12033093   2.36053394
   2.75037941 129.10974029]
2.3236767077234046 3.050662357379447

Estrella 9

str_9 = gr_img.copy()
str_9 = str_9[277:284,1071:1076]
plt.imshow(str_9, cmap='gist_gray')
plt.show()

p0 = [10, 3, 6, 0.8, 1.4, 0, 0]
determinar(str_9,p0,FWHM_xi,FWHM_yi,incer=str_9.flatten())

[190.85299257   1.24748929   0.89798364   2.10782787   3.12535911
   1.72494314  61.26764694]
2.937612779019202 2.1145898720577216

Estrella 10

str_10 = gr_img.copy()
str_10 = str_10[43:49,883:889]
plt.imshow(str_10, cmap='gist_gray')
plt.show()

p0 = [43, 3.4, 6, 0.8, 1.4, 0, 0]
determinar(str_10,p0,FWHM_xi,FWHM_yi,incer=str_10.flatten())

[219.64233685   1.0717617    0.97269196   2.28430393   2.86515034
  -1.59071136  62.34114315]
2.523805945849836 2.2905145186106575

Los valores para el Full Width at Half Maximum (FHWM)

Considere el valor en dirección x & en dirección y. Saque la desviación estándar, el promedio, la mediana y el histograma

FWHM_xi = abs(np.array(FWHM_xi))
print(FWHM_xi)
stdxi = FWHM_xi.std()
print(stdxi)
meanxi = FWHM_xi.mean()
print(meanxi)
medianxi = np.median(FWHM_xi)
print(medianxi)
plt.hist(FWHM_xi.flatten())
plt.title('Valores FHWM en la dirección x para escala de grises considerando incertidumbre')
plt.xlabel('Valor FHWM')
plt.ylabel('No de imagénes con ese valor')
plt.show()

[2.08205473 2.51426496 6.04578332 1.73007771 2.65743574 2.01051951
 4.34098536 2.32367671 2.93761278 2.52380595]
1.2442375593713781
2.916621675050629
2.519035453328353

FWHM_yi = abs(np.array(FWHM_yi))
print(FWHM_yi)
stdyi = FWHM_yi.std()
print(stdyi)
meanyi = FWHM_yi.mean()
print(meanyi)
medianyi = np.median(FWHM_yi)
print(medianyi)
plt.hist(FWHM_yi.flatten())
plt.title('Valores FHWM en la dirección y para escala de grises considerando incertidumbre')
plt.xlabel('Valor FHWM')
plt.ylabel('No de imagénes con ese valor')
plt.show()

[2.28885987 2.38436131 5.84888198 2.31351194 3.02717611 2.48594945
 4.90151913 3.05066236 2.11458987 2.29051452]
1.2080118897663938
3.0706026542016884
2.4351553795125156

Canal rojo

FWHM_xir = []
FWHM_yir = []

Estrella 1

str_1r = img[:,:,0].copy()
str_1r = str_1r[283:288,890:895] 
plt.imshow(str_1r, cmap ='Reds')
plt.show()

p0 = [70, 1.2, 11, 1, 1, 0, 0]
determinar(str_1r,p0,FWHM_xir,FWHM_yir, incer = str_1r.flatten())

[182.91431936   0.91235058   1.00060875   1.99849318   2.04670046
  -0.31143976  90.21305084]
2.1484214221733677 2.3562535363489063

Estrella 2

str_2r = img[:,:,0].copy()
str_2r = str_2r[252:257,949:955]
plt.imshow(str_2r, cmap='Reds')
plt.show()

p0 = [24, 1.3, 2.5, 4, 1, 0, 0]
determinar(str_2r,p0,FWHM_xir,FWHM_yir, incer = str_2r.flatten())

[192.45744813   0.97887011   1.08961655   2.35993751   2.25608062
  -6.34642439  84.92551405]
2.3050629462193077 2.5658509037264716

Estrella 3

str_3r = img[:,:,0].copy()
str_3r = str_3r[309:324,619:636]
plt.imshow(str_3r, cmap='Reds')
plt.show()

p0 = [11, 15, 9, 4, 1, 0, 0]
determinar(str_3r,p0,FWHM_xir,FWHM_yir, incer = str_3r.flatten())

[141.9874119    2.34167953   2.19328524   8.28237019   6.86111281
  20.41751836 114.33748466]
5.514233897956943 5.16479203779562

Estrella 4

str_4r = img[:,:,0].copy()
str_4r = str_4r[18:23,607:612]
plt.imshow(str_4r, cmap='Reds')
plt.show()

p0 = [12, 3, 6, 0.9, 1.4, 0, 0]
determinar(str_4r,p0,FWHM_xir,FWHM_yir, incer = str_4r.flatten())

[179.47671572   0.73445207   0.9714805    1.95111908   2.27151171
   3.04648514  93.82977231]
1.729502460722425 2.287661764280814

Estrella 5

str_5r = img[:,:,0].copy()
str_5r = str_5r[347:356,542:552]
plt.imshow(str_5r, cmap='Reds')
plt.show()

p0 = [15, 3.4, 6, 0.8, 1.4, 0, 0]
determinar(str_5r,p0,FWHM_xir,FWHM_yir, incer = str_5r.flatten())

[162.50745202   1.0652327    1.2441567    4.57620325   4.20175448
   3.18754763 128.05035619]
2.5084313108756726 2.929765136330512

Estrella 6

str_6r = img[:,:,0].copy()
str_6r = str_6r[295:300,307:312]
plt.imshow(str_6r, cmap='Reds')
plt.show()

p0 = [10, 3.3, 5.9, 0.78, 1.7, 0, 0]
determinar(str_6r,p0,FWHM_xir,FWHM_yir, incer = str_6r.flatten())

[150.54027437   1.06288385   0.86245768   2.07589231   1.71659829
  -1.80244502 125.04225717]
2.5029001928382875 2.0309326365113836

Estrella 7

str_7r = img[:,:,0].copy()
str_7r = str_7r[142:154,387:399]
plt.imshow(str_7r, cmap='Reds')
plt.show()

p0 = [14, 10, 15, 1, 5, 0, 0]
determinar(str_7r,p0,FWHM_xir,FWHM_yir, incer = str_7r.flatten())

[137.08273591   1.97954232   2.26715815   6.20129932   5.63488281
   0.52724304 128.6534564 ]
4.66146594576438 5.338749467405802

Estrella 8

str_8r = img[:,:,0].copy()
str_8r = str_8r[459:465,322:327]
plt.imshow(str_8r, cmap='Reds')
plt.show()

p0 = [11, 3.43, 6.4, 0.01, 0.5, 0, 0]
determinar(str_8r,p0,FWHM_xir,FWHM_yir, incer = str_8r.flatten())

[118.37738938   1.2301783    0.88804317   2.13659917   2.33042575
   4.36580063 150.26652019]
2.896848514525203 2.091181859056181

Estrella 9

str_9r = img[:,:,0].copy()
str_9r = str_9r[277:284,1071:1076]
plt.imshow(str_9r, cmap='Reds')
plt.show()

p0 = [13, 3, 6, 0.8, 1.4, 0, 0]
determinar(str_9r,p0,FWHM_xir,FWHM_yir, incer = str_9r.flatten())

[185.59903905   1.27508233   0.91857923   2.11100885   3.13854179
   1.72086826  80.86325465]
3.0025894335728154 2.163088772629732

Estrella 10

str_10r = img[:,:,0].copy()
str_10r = str_10r[43:49,883:889]
plt.imshow(str_10r, cmap='Reds')
plt.show()

p0 = [700, 3.4, 6, 0.8, 1.4, 0, 0]
determinar(str_10r,p0,FWHM_xir,FWHM_yir, incer = str_10r.flatten())

[214.94414166   1.08939095   0.9932327    2.27855954   2.86276987
   1.57713906  79.07029388]
2.565319646114155 2.338884281994375

Los valores para el Full Width at Half Maximum (FHWM)

Considere el valor en dirección x & en dirección y. Saque la desviación estándar, el promedio, la mediana y el histograma

FWHM_xir = abs(np.array(FWHM_xir))
print(FWHM_xir)
stdxir = FWHM_xir.std()
print(stdxir)
meanxir = FWHM_xir.mean()
print(meanxir)
medianxir = np.median(FWHM_xir)
print(medianxir)
plt.hist(FWHM_xir.flatten())
plt.title('Valores FHWM en la dirección x para la image roja')
plt.xlabel('Valor FHWM')
plt.ylabel('No de imagénes con ese valor')
plt.show()

[2.14842142 2.30506295 5.5142339  1.72950246 2.50843131 2.50290019
 4.66146595 2.89684851 3.00258943 2.56531965]
1.1220401417419346
2.9834775770762554
2.536875478494914

FWHM_yir = abs(np.array(FWHM_yir))
print(FWHM_yir)
stdyir = FWHM_yir.std()
print(stdyir)
meanyir = FWHM_yir.mean()
print(meanyir)
medianyir = np.median(FWHM_yir)
print(medianyir)
plt.hist(FWHM_yir.flatten())
plt.title('Valores FHWM en la dirección y para la image roja')
plt.xlabel('Valor FHWM')
plt.ylabel('No de imagénes con ese valor')
plt.show()

[2.35625354 2.5658509  5.16479204 2.28766176 2.92976514 2.03093264
 5.33874947 2.09118186 2.16308877 2.33888428]
1.1882827213448535
2.9267160396079794
2.347568909171641

Canal verde

FWHM_xig = []
FWHM_yig = []

Estrella 1

str_1g = img[:,:,1].copy()
str_1g = str_1g[283:288,890:895] 
plt.imshow(str_1g, cmap ='Greens')
plt.show()

p0 = [71, 1.5, 11, 1, 1, 0, 0]
determinar(str_1g,p0,FWHM_xig,FWHM_yig, incer = str_1g.flatten())

[194.58735935   0.87499245   0.96168531   1.99756044   2.04355755
  -0.315852    68.48717297]
2.060449762326214 2.2645958531909964

Estrella 2

str_2g = img[:,:,1].copy()
str_2g = str_2g[252:257,949:955]
plt.imshow(str_2g, cmap='Greens')
plt.show()

p0 = [22, 1.3, 2.5, 4, 1, 0, 0]
determinar(str_2g,p0,FWHM_xig,FWHM_yig, incer = str_2g.flatten())

[202.16491005   1.00204564   1.05684175   2.36233933   2.22523495
  -3.24379486  71.25427427]
2.35963715103674 2.4886721282798

Estrella 3

str_3g = img[:,:,1].copy()
str_3g = str_3g[309:324,619:636]
plt.imshow(str_3g, cmap='Greens')
plt.show()

p0 = [50, 15, 10, 4, 1, 0, 0]
determinar(str_3g,p0,FWHM_xig,FWHM_yig, incer = str_3g.flatten())

[168.35593595   2.52518256   2.46583528   8.20589042   6.75978717
  -1.67205716  96.27653287]
5.946350509252893 5.806598337482345

Estrella 4

str_4g = img[:,:,1].copy()
str_4g = str_4g[18:23,607:612]
plt.imshow(str_4g, cmap='Greens')
plt.show()

p0 = [14, 3.4, 6, 0.9, 1.4, 0, 0]
determinar(str_4g,p0,FWHM_xig,FWHM_yig, incer = str_4g.flatten())

[181.55452689   0.73372479   0.98023727   1.94887108   2.25037477
   3.05065648  77.61482048]
1.7277898473442026 2.3082823783862074

Estrella 5

str_5g = img[:,:,1].copy()
str_5g = str_5g[347:356,542:552]
plt.imshow(str_5g, cmap='Greens')
plt.show()

p0 = [17, 3.4, 6, 0.8, 1.4, 0, 0]
determinar(str_5g,p0,FWHM_xig,FWHM_yig, incer = str_5g.flatten())

[171.08439474   1.09751253   1.27162261   4.58951748   4.21422199
   3.22642819 113.06046966]
2.5844445024818627 2.994442416716548

Estrella 6

str_6g = img[:,:,1].copy()
str_6g = str_6g[295:300,307:312]
plt.imshow(str_6g, cmap='Greens')
plt.show()

p0 = [32, 3.5, 6, 0.78, 1.7, 0, 0]
determinar(str_6g,p0,FWHM_xig,FWHM_yig, incer = str_6g.flatten())

[157.82945535   0.84466944   1.05013524   2.09746462   1.69957051
  -0.21827805 104.40660141]
1.9890445292987982 2.4728795039574236

Estrella 7

str_7g = img[:,:,1].copy()
str_7g = str_7g[142:154,387:399]
plt.imshow(str_7g, cmap='Greens')
plt.show()

p0 = [11, 2, 5, 3, 5, 0, 0]
determinar(str_7g,p0,FWHM_xig,FWHM_yig, incer = str_7g.flatten())

[142.44945667   1.81945909   2.06026126   6.18751642   5.62400824
   0.67907496 106.22773724]
4.284498726614627 4.851544509614211

Estrella 8

str_8g = img[:,:,1].copy()
str_8g = str_8g[459:465,322:327]
plt.imshow(str_8g, cmap='Greens')
plt.show()

p0 = [16, 3.43, 6.4, 1, 0.5, 0, 0]
determinar(str_8g,p0,FWHM_xig,FWHM_yig, incer = str_8g.flatten())

[132.27002448   0.97039729   1.27318881   2.15941699   2.35249593
  -0.40139803 128.14877035]
2.2851110001414745 2.9981305351356577

Estrella 9

str_9g = img[:,:,1].copy()
str_9g = str_9g[277:284,1071:1076]
plt.imshow(str_9g, cmap='Greens')
plt.show()

p0 = [10, 3, 6, 0.8, 1.4, 0, 0]
determinar(str_9g,p0,FWHM_xig,FWHM_yig, incer = str_9g.flatten())

[193.35932279   1.23687204   0.89029534   2.10755642   3.13488928
   1.72376695  61.15734226]
2.9126110808108847 2.0964853040996143

Estrella 10

str_10g = img[:,:,1].copy()
str_10g = str_10g[43:49,883:889]
plt.imshow(str_10g, cmap='Greens')
plt.show()

p0 = [10, 3.3, 3, 0.9, 1.4, 0, 0]
determinar(str_10g,p0,FWHM_xig,FWHM_yig, incer = str_10g.flatten())

[223.03027334   0.97687939   1.06282319   2.28196065   2.86747905
 -15.75893458  61.71980515]
2.300375167201353 2.5027573548739817

Los valores para el Full Width at Half Maximum (FHWM)

Considere el valor en dirección x & en dirección y. Saque la desviación estándar, el promedio, la mediana y el histograma

FWHM_xig = abs(np.array(FWHM_xig))
print(FWHM_xig)
stdxig = FWHM_xig.std()
print(stdxr)
meanxig = FWHM_xig.mean()
print(meanxig)
medianxig = np.median(FWHM_xig)
print(medianxig)
plt.hist(FWHM_xig.flatten())
plt.title('Valores FHWM en la dirección x para la imagen verde considerando incertidumbre')
plt.xlabel('Valor FHWM')
plt.ylabel('No de imagénes con ese valor')
plt.show()

[2.06044976 2.35963715 5.94635051 1.72778985 2.5844445  1.98904453
 4.28449873 2.285111   2.91261108 2.30037517]
1.1372820062135576
2.845031227650905
2.330006159119047

FWHM_yig = abs(np.array(FWHM_yig))
print(FWHM_yig)
stdyig = FWHM_yig.std()
print(stdyig)
meanyig = FWHM_yig.mean()
print(meanyig)
medianyig = np.median(FWHM_yig)
print(medianyig)
plt.hist(FWHM_yig.flatten())
plt.title('Valores FHWM en la dirección y para la imagen verde considerando incertidumbre')
plt.xlabel('Valor FHWM')
plt.ylabel('No de imagénes con ese valor')
plt.show()

[2.26459585 2.48867213 5.80659834 2.30828238 2.99444242 2.4728795
 4.85154451 2.99813054 2.0964853  2.50275735]
1.1775361127169395
3.0784388321736786
2.495714741576891

Canal azul

FWHM_xib = []
FWHM_yib = []

Estrella 1

str_1b = img[:,:,2].copy()
str_1b = str_1b[283:288,890:895] 
plt.imshow(str_1b, cmap ='Blues')
plt.show()

p0 = [71, 1.7, 11, 1, 1, 0, 0]
determinar(str_1b,p0,FWHM_xib,FWHM_yib, incer = str_1b.flatten())

[196.08431561   0.86622721   0.95511986   2.01460987   2.04323191
  -0.27430574  54.08557709]
2.039809207469812 2.249135380109503

Estrella 2

str_2b = img[:,:,2].copy()
str_2b = str_2b[252:257,949:955]
plt.imshow(str_2b, cmap='Blues')
plt.show()

p0 = [50, 2.9, 2, 4, 1, 0, 0]
determinar(str_2b,p0,FWHM_xib,FWHM_yib, incer = str_2b.flatten())

[207.39523454   1.05340756   1.02290869   2.3785664    2.2034983
   1.28470635  58.97679989]
2.480585247068367 2.4087658908107965

Estrella 3

str_3b = img[:,:,2].copy()
str_3b = str_3b[309:324,619:636]
plt.imshow(str_3b, cmap='Blues')
plt.show()

p0 = [14, 15, 10, 4, 1, 0, 0]
determinar(str_3b,p0,FWHM_xib,FWHM_yib, incer = str_3b.flatten())

[189.3566925    2.75095638   2.67312878   8.30245528   6.75026981
  -4.8848669   80.55741834]
6.47800722275355 6.294737245460569

Estrella 4

str_4b = img[:,:,2].copy()
str_4b = str_4b[18:23,607:612]
plt.imshow(str_4b, cmap='Blues')
plt.show()

p0 = [9, 3, 6, 1.23, 1.46, 0, 0]
determinar(str_4b,p0,FWHM_xib,FWHM_yib, incer = str_4b.flatten())

[183.38176381   0.73671432   0.99916499   1.94780626   2.22736734
   3.06024167  63.6038059 ]
1.7348296501661935 2.3528537457692438

Estrella 5

str_5b = img[:,:,2].copy()
str_5b = str_5b[347:356,542:552]
plt.imshow(str_5b, cmap='Blues')
plt.show()

p0 = [8, 1.7, 1, 3.5, 1, 0, 0]
determinar(str_5b,p0,FWHM_xib,FWHM_yib, incer = str_5b.flatten())

[166.04481043  -1.21990199   1.38650382   4.52448629   4.2147847
  -2.61241384 104.34156413]
-2.8726496561210997 3.2649669805187034

Estrella 6

str_6b = img[:,:,2].copy()
str_6b = str_6b[295:300,307:312]
plt.imshow(str_6b, cmap='Blues')
plt.show()

p0 = [123, 3, 5.9, 0.78, 1.7, 0, 0]
determinar(str_6b,p0,FWHM_xib,FWHM_yib, incer = str_6b.flatten())

[159.01793668   0.8552177    1.05546949   2.07193554   1.69511554
  -0.21091069  89.9142516 ]
2.0138837865360038 2.4854407042295907

Estrella 7

str_7b = img[:,:,2].copy()
str_7b = str_7b[142:154,387:399]
plt.imshow(str_7b, cmap='Blues')
plt.show()

p0 = [100, 2, 5, 3, 5, 0, 0]
determinar(str_7b,p0,FWHM_xib,FWHM_yib, incer = str_7b.flatten())

[137.6552333    1.92930042   1.73162599   6.12640336   5.62271455
  -0.84345339  86.564405  ]
4.543155312707578 4.077667599505995

Estrella 8

str_8b = img[:,:,2].copy()
str_8b = str_8b[459:465,322:327]
plt.imshow(str_8b, cmap='Blues')
plt.show()

p0 = [30, 3.43, 6.4, 0.01, 0.5, 0, 0]
determinar(str_8b,p0,FWHM_xib,FWHM_yib, incer = str_8b.flatten())

[152.36430635   1.12628322   1.39881391   2.06937349   2.37438637
  -0.44230575 105.37899408]
2.6521942943381274 3.293955036969304

Estrella 9

str_9b = img[:,:,2].copy()
str_9b = str_9b[277:284,1071:1076]
plt.imshow(str_9b, cmap='Blues')
plt.show()

p0 = [10, 3, 6.4, 7.3, 1.1, 0, 0]
determinar(str_9b,p0,FWHM_xib,FWHM_yib, incer = str_9b.flatten())

[193.86710559   0.88686474   1.23340715   2.10631332   3.10028245
 -15.54702541  41.51349323]
2.088406856621961 2.9044518727402315

Estrella 10

str_10b = img[:,:,2].copy()
str_10b = str_10b[43:49,883:889]
plt.imshow(str_10b, cmap='Blues')
plt.show()

p0 = [41, 12, 4, 1, 1, 0, 0]
determinar(str_10b,p0,FWHM_xib,FWHM_yib, incer = str_10b.flatten())

[ 2.21501512e+02  9.47560514e-01  1.06364334e+00  2.29364326e+00
  2.86413861e+00 -1.70292345e-02  4.60375723e+01]
2.231334493416455 2.5046886503061434

Los valores para el Full Width at Half Maximum (FHWM)

Considere el valor en dirección x & en dirección y. Saque la desviación estándar, el promedio, la mediana y el histograma

FWHM_xib = abs(np.array(FWHM_xib))
print(FWHM_xib)
stdxib = FWHM_xib.std()
print(stdxib)
meanxib = FWHM_xib.mean()
print(meanxib)
medianxib = np.median(FWHM_xib)
print(medianxib)
plt.hist(FWHM_xib.flatten())
plt.title('Valores FHWM en la dirección x para la image azul')
plt.xlabel('Valor FHWM')
plt.ylabel('No de imagénes con ese valor')
plt.show()

[2.03980921 2.48058525 6.47800722 1.73482965 2.87264966 2.01388379
 4.54315531 2.65219429 2.08840686 2.23133449]
1.4045244740696985
2.9134855727199147
2.355959870242411

FWHM_yib = abs(np.array(FWHM_yib))
print(FWHM_yib)
stdyib = FWHM_yib.std()
print(stdyib)
meanyib = FWHM_yib.mean()
print(meanyib)
medianyib = np.median(FWHM_yib)
print(medianyib)
plt.hist(FWHM_yib.flatten())
plt.title('Valores FHWM en la dirección y para la image azul')
plt.xlabel('Valor FHWM')
plt.ylabel('No de imagénes con ese valor')
plt.show()

[2.24913538 2.40876589 6.29473725 2.35285375 3.26496698 2.4854407
 4.0776676  3.29395504 2.90445187 2.50468865]
1.1699287444022204
3.183666310642008
2.7045702615231875

Comparemos los cuatro valores

En dirección x

print(stdxi,stdxir,stdxig,stdxib)
print(meanxi,meanxir,meanxig,meanxib)
print(medianxi,medianxir,medianxig,medianxib)

1.2442375593713781 1.1220401417419346 1.233085493598224 1.4045244740696985
2.916621675050629 2.9834775770762554 2.845031227650905 2.9134855727199147
2.519035453328353 2.536875478494914 2.330006159119047 2.355959870242411

En dirección y

print(stdyi,stdyir,stdyig,stdyib)
print(meanyi,meanyir,meanyig,meanyib)
print(medianyi,medianyir,medianyig,medianyib)

1.2080118897663938 1.1882827213448535 1.1775361127169395 1.1699287444022204
3.0706026542016884 2.9267160396079794 3.0784388321736786 3.183666310642008
2.4351553795125156 2.347568909171641 2.495714741576891 2.7045702615231875

Comparando los valores con y sin incertidumbre

Para x

# -- crear un array con todas las cantidades calculadas anteriormente para el FWHM
valuesx = np.array([[stdx,stdxr,stdxg,stdxb, meanx,meanxr,meanxg,meanxb, medianx,medianxr,medianxg,medianxb],
         [stdxi,stdxir,stdxig,stdxib, meanxi,meanxir,meanxig,meanxib, medianxi,medianxir,medianxig,medianxib]])
valuesx.shape

(2, 12)

# -- crear un data frame con esa información
FWHMxx = pd.DataFrame(data=valuesx,
             columns=['stdxL','stdxR','stdxG','stdxB','meanxL','meanxR','meanxG','meanxB','medianxL','medianxR','medianxG','medianxB'], 
             index=['Sin-incertidumbre','Con-incertidumbre'])
# -- crear un data frame con la diferencia entre valores con y sin incertidumbre
resta = pd.DataFrame(data=abs(FWHMxx.loc['Sin-incertidumbre'] - FWHMxx.loc['Con-incertidumbre']).values.reshape(1,12), 
                   columns=['stdxL','stdxR','stdxG','stdxB','meanxL','meanxR','meanxG','meanxB','medianxL','medianxR','medianxG','medianxB'],
                   index=['diferencia'])
# -- combinar ambos data frame
FWHMxx = FWHMxx.append(resta)

Para y

# -- lo mismo pero en dirección y
valuesy = np.array([[stdy,stdyr,stdyg,stdyb, meany,meanyr,meanyg,meanyb, mediany,medianyr,medianyg,medianyb],
         [stdyi,stdyir,stdyig,stdyib, meanyi,meanyir,meanyig,meanyib, medianyi,medianyir,medianyig,medianyib]])

FWHMyy = pd.DataFrame(data=valuesy,
             columns=['stdyL','stdyR','stdyG','stdyB','meanyL','meanyR','meanyG','meanyB','medianyL','medianyR','medianyG','medianyB'], 
             index=['Sin-incertidumbre','Con-incertidumbre'])

resta = pd.DataFrame(data=abs(FWHMyy.loc['Sin-incertidumbre'] - FWHMyy.loc['Con-incertidumbre']).values.reshape(1,12), 
                   columns=['stdyL','stdyR','stdyG','stdyB','meanyL','meanyR','meanyG','meanyB','medianyL','medianyR','medianyG','medianyB'],
                   index=['diferencia'])
FWHMyy = FWHMyy.append(resta)

FWHMxx

	stdxL	stdxR	stdxG	stdxB	meanxL	meanxR	meanxG	meanxB	medianxL	medianxR	medianxG	medianxB
Sin-incertidumbre	1.202031	1.137282	1.192085	1.301895	2.938495	2.834198	2.824878	2.982611	2.535912	2.463149	2.424866	2.517526
Con-incertidumbre	1.244238	1.122040	1.233085	1.404524	2.916622	2.983478	2.845031	2.913486	2.519035	2.536875	2.330006	2.355960
diferencia	0.042207	0.015242	0.041001	0.102629	0.021874	0.149280	0.020153	0.069125	0.016877	0.073726	0.094860	0.161566

FWHMyy

	stdyL	stdyR	stdyG	stdyB	meanyL	meanyR	meanyG	meanyB	medianyL	medianyR	medianyG	medianyB
Sin-incertidumbre	1.163140	1.076054	1.124040	1.219566	3.058655	3.106182	3.099311	3.108883	2.535681	2.816984	2.788384	2.607229
Con-incertidumbre	1.208012	1.188283	1.177536	1.169929	3.070603	2.926716	3.078439	3.183666	2.435155	2.347569	2.495715	2.704570
diferencia	0.044872	0.112229	0.053496	0.049637	0.011948	0.179466	0.020872	0.074783	0.100526	0.469415	0.292670	0.097341

Al incluir la incertidumbre los valores de FWHM si cambian. No se puede afirmar que el valor sea mas o menos preciso, sin embargo al darle un peso a cada pixel que conforma la imagen de la estrella, es mas facil manejar imagenes con mayor ruido como la Estrella 3 y 7.

El módelo es extremadamente sensible a cualquier cambio como el cambio de los parametros iniciales p0, cambio en los pixeles alrededor de la estrella.