Entrega.md

Estudiante: Angie Nicole Hernández Durán - UIS

Numpy y optimizción con SciPy

Este ejercicio consiste en conseguir una medición de la resolución espacial a partir de una foto del cielo estrellado. Específicamente se calculará la

\text{anchura a media altura}

(FWHM), la cual es una medida derivada de uno de los datos que se puede obtener de un ajuste gaussiano, la desviación estándar. Para llegar a esta medición tomaremos los pasos que se presentan a continuación.

Primero se leyeron los datos correspondientes, en este caso, una imagen en formato jpeg.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import leastsq 

sky = plt.imread('data/zapatocaImage.jpeg')

type(sky)   #La imagen es efectivamente un array de numpy

numpy.ndarray

sky.shape   #Vemos las dimensiones del arreglo

(789, 1184, 3)

plt.imshow(sky)  #Y visualizamos la imagen

<matplotlib.image.AxesImage at 0x7f3cd5f6d208>

La imágen guardada en "sky" está a color, ahora se deben combinar los 3 arrays para generar una versión a blanco y negro. Una manera de hacer esto es añadir los 3 canales y luego dividirlos entre 3 para obtener un promedio de los valores de cada canal, que darán los valores de gris de los pixeles:

sky_bw = (sky[:,:,0] + sky[:,:,1] + sky[:,:,2]) / 3

plt.imshow(sky_bw, cmap = 'gray')    #Nuevamente visualizo

<matplotlib.image.AxesImage at 0x7f3cd5ed4e10>

Lo anterior no se ve como las imagenes a blanco y negro a las que estamos acostumbrados, esto es debido a que se suele usar un método que consiste en dar diferentes pesos a cada canal, ya que las diferentes longitudes de onda brindan diferentes sensaciones oculares al ojo humano. A continuación se utilizó dicho método.

sky_bw = sky[:,:,0] * 0.3 + sky[:,:,1] * 0.59 + sky[:,:,2] * 0.11

plt.imshow(sky_bw, cmap = 'gray')     #Visualizamos de nuevo

<matplotlib.image.AxesImage at 0x7f3cd5c6b208>

Dando diferentes contribuciones a los diferentes canales (30% para rojo, 59% para verde y 11% para azul) se obtuvo una imagen mucho más agradable a la vista y que recuerda más a la imagen original.

Ahora es necesario seleccionar una pequeña región en la imágen donde haya una estrella. las coordenadas de los pixeles que delimitan la estrella se hallaron utilizando una págia web externa que permite seleccionar manualmente una región rectangular de una imágen. Este método se utilizó igualmente con la selección de las regiones de todas las estrellas.

star1 = sky_bw[10:19,228:241]       #La estrella es un rectángulo obtenido de la imagen original

plt.imshow(star1, cmap = 'gray')    #Visusalizamos para ver que la región seleccionada corresponde a una estrella

<matplotlib.image.AxesImage at 0x7f3cd5e589e8>

El seguiente movimiento consiste en ajustar una gaussiana simétrica bidimensional a la imagen de la estrella, para lo cual se llevan a cabo los pasos a continuación.

#Defino una función que me de los valores de una gausiana bidimensional.

#tupla es un arreglo donde van los parámetros de la gausiana, la amplitud a, una constante aditiva b, 
#la desviación estándar c, x0 y y0 son las coordenadas donde está centrada la función.

def gauss2D(tupla, x, y):
    
    #tupla es un arreglo donde van los parámetros de la gausiana, la amplitud a, una constante aditiva b, 
    #la desviación estándar c, x0 y y0 son las coordenadas donde está centrada la función.
    
    a = tupla[0]
    b = tupla[1]
    c = tupla[2]
    x0 = tupla[3]
    y0 = tupla[4]
    
    exponente = -((x-x0)**2 + (y-y0)**2) / (2*c**2)
    z1 = a * np.exp(exponente) + b
    
    return z1

#Hago la malla donde pondré mi ajuste gausiano, esta debe ser del mismo tamaño de la imagen que quiero 
#ajustar, la de mi estrella

x = np.arange(0,star1.shape[1],1)
y = np.arange(0,star1.shape[0],1)

xx, yy = np.meshgrid(x, y)

z = star1  #Mis datos reales son los valores de gris de mi estrella

#Defino el error, el cual consiste en la diferencia entre los valores de mi ajuste gaussiano y los valores reales
#de mi estrella

def errormodel(tupla, x,y,z):
    
    #Entra la tupla de parámetros para la función gaussiana, los valores de x,y necesarios para hacer la malla
    #donde esta irá, y los valores reales (z) de los pixeles de la estrella

    m = np.ravel(gauss2D(tupla,x,y) - z)  #Estaba teniendo problemas con la salida de esta función cuando estraba
                                          #en leastsq, np.ravel aplana mi arreglo 2D y arregla el fomato de salida 
    return m

#Para realizar la optimización de parámetros voy a utilizar la función leastsq, que utiliza el método de 
#mínimos cuadrados.

p0 = [1.0, 0.0, 1.0, 0.0, 0.0]   #Esta es mi primera aleatoria elección de parámetros que entrará a la función
                                 #leastsq

best, suss = leastsq(errormodel, p0, args=(xx,yy,z))  #finalmente uso la función 

params = best                    #Extraigo los parámetros que mejor se ajustan a la estrella

zz = gauss2D(params, xx,yy)      #Utilizo los parpametros obtenidos

plt.imshow(zz,cmap="gray")       #Y visualizo la gausiana que mejor se ajusta a mi estrella
#plt.colorbar()

<matplotlib.image.AxesImage at 0x7f3cd5f77160>

Ahora debo repetir este proceso para varias estrellas, voy a implementar los anteriores pasos en una función para que el proceso no sea tedioso.

def get_param(xi,xf,yi,yf):
    
    #xi, xf, yi, yf son las coordenadas iniciales y finales que delimitan el área de mi estrella en la imágen
    #La función devuelve los parámetros que mejor se ajustan a la estrella, así como los valores necesarios que 
    #hacen la malla en la cual va el ajuste gaussiano, esto para fines de graficación.
    
    star = sky_bw[yi:yf,xi:xf]
    
    x = np.arange(0,star.shape[1],1)    #Todas mis estrellas son de diferentes tamaños
    y = np.arange(0,star.shape[0],1)
    
    xx, yy = np.meshgrid(x, y)
    
    p0 = [1.0, 0.0, 1.0, 0.0, 0.0]      #Dejo la misma elección inicial de parámetros para todas las estrellas
    
    best, suss = leastsq(errormodel, p0, args=(xx,yy,star))  
    
    return best, xx, yy                 

#El siguiente es un arreglo donde coloco las coordenadas de los pixeles que delimitan areas de varias estrellas.

#El orden es yi, yf, xi, xf

stars_coord = np.array([[305, 328, 620, 637],
                      [82, 96, 628, 637],
                      [115, 126, 726, 734],
                      [368, 386, 444, 458],
                      [264, 276, 749, 758],
                      [540, 564, 650, 676],
                      [452, 459, 205, 215],
                      [87, 100, 1080, 1092],
                      [177, 185, 1096, 1110],
                      [21, 30, 921, 928],
                      [307, 328, 617, 638],
                      [260, 268, 764, 773],
                      [15, 25, 673, 682],
                      [52, 63, 452, 460],
                      [10, 19, 228, 241],
                      [142, 156,387,400]])  

#Visualizo las areas seleccionadas para asegurarme que correspoden a estrellas

for i in range(0,len(stars_coord)):
    
    star = sky_bw[stars_coord[i,0]:stars_coord[i,1], stars_coord[i,2]:stars_coord[i,3]]
    plt.imshow(star, cmap = 'gray')
    plt.show()                       

#Paso el arreglo de estrellas a la función para obtener los parámetros que mejor se ajustan.

params = []  #la lista donde guardaré los parámetros óptimos de todas las estrellas.
x_ = []
y_ = []      #Arreglos de los tamaños de mis estrellas, para fines de graficación futuros.

for i in range(0,stars_coord.shape[0]):
    m , xm, ym= get_param(stars_coord[i,2],stars_coord[i,3],stars_coord[i,0],stars_coord[i,1])
    x_.append(xm)
    y_.append(ym)
    params.append(m)

#Visualizo las gaussianas que mejor se ajustan a mis estrellas, se ven bien

for i in range(0,len(FWHM)):
    zz = gauss2D(params[i], x_[i],y_[i]) 
    plt.imshow(zz, cmap = 'gray')
    plt.show() 

En la lista llamada 'params' están los parámetros, por columnas, a, b, c, x0, y0 que mejor se ajustan a cada una de mis estrellas. Con el parámetro c, la desviación estándar, puedo calcular la FWHM mediante la expresión:

FWHM = 2\sqrt(2 \ln 2)c\approx 2.355 c

Extraigo este parámetro para hacer el cálculo.

#Obtengo la FWHM usando la formula anterior

desvesta = [np.absolute(row[2]) for row in params]         #Extraigo el tercer elemento (c) de cada fila

FWHM = 2.0 * np.sqrt(2.0 * np.log(2.0)) * np.array(desvesta)

Tengo muy, muy pocos datos (16), sin embargo voy a derivar de estos algunas medidas estadísticas:

#Primero, un histograma 

n, bins, patches = plt.hist(FWHM, len(FWHM), facecolor = 'gray', alpha = 0.5)
plt.title('Histograma del ancho a media altura de mis estrellas seleccionadas')
plt.xlabel('FWHM [pixeles]')
plt.ylabel('frecuencia')
plt.show()

En el anterior histograma podemos ver que la mayor parte de las medidas se encuentran al rededor del valor de 2 pixeles. A continuación calcularmos algunas otras medidas.

#Ahora algunas otras medidas

average_FWHM_bw = np.average(FWHM)    #Calculo el promedio,
median_FWHM_bw = np.median(FWHM)      #mediana,
std_FWHM_bw = np.std(FWHM)            #desviacón estándar

print('El promedio de las medidas es: ' + str(average_FWHM_bw) + '\n'
     'La mediana es: ' + str(median_FWHM_bw) + '\n'
     'La desviación estándar es: ' + str(std_FWHM_bw))

El promedio de las medidas es: 3.432858986192058
La mediana es: 2.260080688320251
La desviación estándar es: 2.369829780288601

Para los diferentes canales

Ahora es necesario seleccionar las mismas regiones de los diferentes canales, para lo cual se modificó un poco la función get_params para incluir la posibilidad de seleccionar los diferentes canales de color.

def get_param_RGB(xi,xf,yi,yf,channel):
    
    if channel == 'R':       #Agrego estas condiciones para tomar los diferentes canales R,G o B
        n = 0
    elif channel == 'G':
        n = 1
    elif channel == 'B':
        n = 2
    
    star = sky[yi:yf,xi:xf,n]        #Se elecciona de la imágen original, el rectángulo del canal que corresponde
    
    x = np.arange(0,star.shape[1],1) #Todas mis estrellas son de diferentes tamaños
    y = np.arange(0,star.shape[0],1)
    
    xx, yy = np.meshgrid(x, y)
    
    p0 = [1.0, 0.0, 1.0, 0.0, 0.0]  #voy a dejar la misma elección inicial de parámetros para todas las estrellas
    
    best, suss = leastsq(errormodel, p0, args=(xx,yy,star))
    
    return best, xx, yy

#Obtengo los mejores ajustes repitiendo lo hecho anteriormente, no guardo el tamaño de los arreglos
# ya que no voy a visualizar más

params_R = []
params_B = [] 
params_G = []

for i in range(0,stars_coord.shape[0]):
    m , xm, ym= get_param_RGB(stars_coord[i,2],stars_coord[i,3],stars_coord[i,0],stars_coord[i,1], channel = 'R')
    params_R.append(m)
    
for i in range(0,stars_coord.shape[0]):
    m , xm, ym= get_param_RGB(stars_coord[i,2],stars_coord[i,3],stars_coord[i,0],stars_coord[i,1], channel = 'G')
    params_G.append(m)
    
for i in range(0,stars_coord.shape[0]):
    m , xm, ym= get_param_RGB(stars_coord[i,2],stars_coord[i,3],stars_coord[i,0],stars_coord[i,1], channel = 'B')
    params_B.append(m)

#Extraigo los valores de la desviación estándar c

desvesta_R = [np.absolute(row[2]) for row in params_R]
desvesta_G = [np.absolute(row[2]) for row in params_G]
desvesta_B = [np.absolute(row[2]) for row in params_B]

#Uso la formula para calcular el FWHM

FWHM_R = 2.0 * np.sqrt(2.0 * np.log(2.0)) * np.array(desvesta_R)
FWHM_G = 2.0 * np.sqrt(2.0 * np.log(2.0)) * np.array(desvesta_G)
FWHM_B = 2.0 * np.sqrt(2.0 * np.log(2.0)) * np.array(desvesta_B)

#Visualizo los histogramas
n_bins = len(FWHM)

fig, axes = plt.subplots(nrows=2, ncols=2, figsize = (15,9))  #Hago mi canvas
ax0, ax1, ax2, ax3 = axes.flatten()

ax0.hist(FWHM_R, n_bins, histtype='bar', color = 'r')
ax0.set_xlabel('FWHM [pixeles] (Rojo)')
ax0.set_ylabel('frecuencia')

ax1.hist(FWHM_G, n_bins, histtype='bar', color = 'g')
ax1.set_xlabel('FWHM [pixeles] (verde)')
ax1.set_ylabel('frecuencia')

ax2.hist(FWHM_B, n_bins ,histtype='bar', color = 'b')
ax2.set_xlabel('FWHM [pixeles] (azul)')
ax2.set_ylabel('frecuencia')

ax3.hist(FWHM, n_bins, histtype='bar', color = 'gray')
ax3.set_xlabel('FWHM [pixeles](blanco y negro)')
ax3.set_ylabel('frecuencia')

Text(0, 0.5, 'frecuencia')

En los histogramas anteriores podemos observar que la mayor cantidad de medidas del FWHM se encuentra cerca al valor de 2 pixeles, también se puede apreciar que hay valores bastante alejados de este, incluso mayores a 10 pixeles, derivadas de un par de estrellas que ocupan un area considerablemente mayor a las demás.

También obtengo las otras medidas estadísticas.

average_FWHM_R = np.average(FWHM_R)
median_FWHM_R = np.median(FWHM_R)    
std_FWHM_R = np.std(FWHM_R)

average_FWHM_G = np.average(FWHM_G)
median_FWHM_G = np.median(FWHM_G)    
std_FWHM_G = np.std(FWHM_G)

average_FWHM_B = np.average(FWHM_B)
median_FWHM_B = np.median(FWHM_B)    
std_FWHM_B = np.std(FWHM_B)

#Las acomodo en un dataframe para su fácil visualización 

medidas = pd.DataFrame([
  [average_FWHM_bw, average_FWHM_R, average_FWHM_G, average_FWHM_B],
  [median_FWHM_bw, median_FWHM_R, median_FWHM_G, median_FWHM_B],
  [std_FWHM_bw, std_FWHM_R, std_FWHM_G, std_FWHM_B],
],
  index = ['promedio', 'media', 'desv estándar'],
  columns = ['B&W', 'canal R', 'canal G', 'canal B']
)

medidas

.dataframe tbody tr th:only-of-type { vertical-align: middle; } <pre data-sourcepos="714:5-720:5"><code>.dataframe tbody tr th { vertical-align: top; } .dataframe thead th { text-align: right; } </code></pre>

	B&W	canal R	canal G	canal B
promedio	3.432859	3.355606	3.436947	3.518292
media	2.260081	2.234073	2.261304	2.295926
desv estándar	2.369830	2.180836	2.395190	2.592449

Voy a acomodarlos en una tabla para que se vean bien

Podemos observar que las medidas son diferentes entre todos los canales, lo cual tiene sentido, ya que cada canal aporta diferentes valores a la imágen final, en el caso de la imagen a color, si cada canal aportase el mismo valor, no se conseguirían los diferentes colores que se requieren. Así mismo, son diferentes a las medidas de la imagen en blanco y negro.

Teniendo en cuenta la incertidumbre

En los cálculos realizados anteriormente no se tuvo en cuenta la incertidumbre de las medidas, en esta sección se incluye de la siguiente manera: se obtiene la raíz cuadrada del valor de cada pixel, y se divide el error de cada pixel por este valor, en la función en la cual se define el error.

#Redefino el error 

def errormodel_uncert(tupla, x,y,z, sigma):
    
    #Entra la tupla de parámetros para la función gaussiana, los valores de x,y necesarios para hacer la malla
    #donde esta irá, y los valores reales (z) de los pixeles de la estrella. Adicionalmente la incertidumbre en
    #un arreglo llamado sigma

    m = np.ravel((gauss2D(tupla,x,y) - z) / sigma)
    
    return m

#Redefino mis funciones para obtener los parámetros óptimos, añadiendo la obtención de sigma y cambiando la
#función de error que uso por aquella que incluye la incertidumbre

def get_param_uncert(xi,xf,yi,yf):
    
    star = sky_bw[yi:yf,xi:xf]
    
    sigma = np.sqrt(star)               #Añado este cálculo
    
    x = np.arange(0,star.shape[1],1)    
    y = np.arange(0,star.shape[0],1)
    
    xx, yy = np.meshgrid(x, y)
    
    p0 = [1.0, 0.0, 1.0, 0.0, 0.0]      
    
    best, suss = leastsq(errormodel_uncert, p0, args=(xx,yy,star,sigma))  
    
    return best, xx, yy  

def get_param_RGB_uncert(xi,xf,yi,yf,channel):
    
    if channel == 'R':       
        n = 0
    elif channel == 'G':
        n = 1
    elif channel == 'B':
        n = 2
    
    star = sky[yi:yf,xi:xf,n]   
    
    sigma = np.sqrt(star)                  #Igualmente incluyo este cálculo
    
    x = np.arange(0,star.shape[1],1) 
    y = np.arange(0,star.shape[0],1)
    
    xx, yy = np.meshgrid(x, y)
    
    p0 = [1.0, 0.0, 1.0, 0.0, 0.0]  
    
    best, suss = leastsq(errormodel_uncert, p0, args=(xx,yy,star,sigma))
    
    return best, xx, yy

Continúo con los pasos previamente realizados, hallando el mejor ajuste y calculando el FWHM y los medidas estadísticas.

#Obtengo los parámetros que mejor se ajustan a cada estrella

params_bw = []

for i in range(0,stars_coord.shape[0]):
    m , xm, ym= get_param_uncert(stars_coord[i,2],stars_coord[i,3],stars_coord[i,0],stars_coord[i,1])
    x_.append(xm)
    y_.append(ym)
    params_bw.append(m)
    
params_R = []
params_B = [] 
params_G = []

for i in range(0,stars_coord.shape[0]):
    m , xm, ym= get_param_RGB_uncert(stars_coord[i,2],stars_coord[i,3],stars_coord[i,0],stars_coord[i,1], channel = 'R')
    params_R.append(m)
    
for i in range(0,stars_coord.shape[0]):
    m , xm, ym= get_param_RGB_uncert(stars_coord[i,2],stars_coord[i,3],stars_coord[i,0],stars_coord[i,1], channel = 'G')
    params_G.append(m)
    
for i in range(0,stars_coord.shape[0]):
    m , xm, ym= get_param_RGB_uncert(stars_coord[i,2],stars_coord[i,3],stars_coord[i,0],stars_coord[i,1], channel = 'B')
    params_B.append(m)

#Extraigo los valores de la desviación estándar c

desvesta_bw = [np.absolute(row[2]) for row in params_bw]
desvesta_R = [np.absolute(row[2]) for row in params_R]
desvesta_G = [np.absolute(row[2]) for row in params_G]
desvesta_B = [np.absolute(row[2]) for row in params_B]

#Uso la formula para calcular el FWHM

FWHM_bw = 2.0 * np.sqrt(2.0 * np.log(2.0)) * np.array(desvesta_bw)
FWHM_R = 2.0 * np.sqrt(2.0 * np.log(2.0)) * np.array(desvesta_R)
FWHM_G = 2.0 * np.sqrt(2.0 * np.log(2.0)) * np.array(desvesta_G)
FWHM_B = 2.0 * np.sqrt(2.0 * np.log(2.0)) * np.array(desvesta_B)

#Grafico los histogramas

n_bins= len(FWHM)

fig, axes = plt.subplots(nrows=2, ncols=2, figsize = (15,9))  #Hago mi canvas
ax0, ax1, ax2, ax3 = axes.flatten()

ax0.hist(FWHM_R, n_bins, histtype='bar', color = 'r')
ax1.set_xlabel('FWHM [pixeles] (rojo)')
ax1.set_ylabel('frecuencia')

ax1.hist(FWHM_G, n_bins, histtype='bar', color = 'g')
ax1.set_xlabel('FWHM [pixeles] (verde)')
ax1.set_ylabel('frecuencia')

ax2.hist(FWHM_B, n_bins ,histtype='bar', color = 'b')
ax1.set_xlabel('FWHM [pixeles] (azul)')
ax1.set_ylabel('frecuencia')

ax3.hist(FWHM, n_bins, histtype='bar', color = 'gray')
ax1.set_xlabel('FWHM [pixeles] (blanco y negro)')
ax1.set_ylabel('frecuencia')

Text(0, 0.5, 'frecuencia')

Estos histogramas no presentan diferencia respecto a los anteriores, las diferencias introducidas por la incertidumbre deben ser más finas y no evidentes en un histograma.

#Calculo las otras medidas

average_FWHM_bw = np.average(FWHM_bw)
median_FWHM_bw = np.median(FWHM_bw)    
std_FWHM_bw = np.std(FWHM_bw)

average_FWHM_R = np.average(FWHM_R)
median_FWHM_R = np.median(FWHM_R)    
std_FWHM_R = np.std(FWHM_R)

average_FWHM_G = np.average(FWHM_G)
median_FWHM_G = np.median(FWHM_G)    
std_FWHM_G = np.std(FWHM_G)

average_FWHM_B = np.average(FWHM_B)
median_FWHM_B = np.median(FWHM_B)    
std_FWHM_B = np.std(FWHM_B)

#Las acomodo en un dataframe para su fácil visualización 

medidas_uncert = pd.DataFrame([
  [average_FWHM_bw, average_FWHM_R, average_FWHM_G, average_FWHM_B],
  [median_FWHM_bw, median_FWHM_R, median_FWHM_G, median_FWHM_B],
  [std_FWHM_bw, std_FWHM_R, std_FWHM_G, std_FWHM_B],
],
  index = ['promedio', 'media', 'desv estándar'],
  columns = ['B&W', 'canal R', 'canal G', 'canal B']
)

medidas_uncert

.dataframe tbody tr th:only-of-type { vertical-align: middle; } <pre data-sourcepos="964:5-970:5"><code>.dataframe tbody tr th { vertical-align: top; } .dataframe thead th { text-align: right; } </code></pre>

	B&W	canal R	canal G	canal B
promedio	3.442413	3.362915	3.548846	3.515920
media	2.213994	2.197444	2.216574	2.192405
desv estándar	2.406436	2.223734	2.421031	2.609954

medidas    #Visualizo las medidas calculadas previamente para hacer una comparación visual de estas

.dataframe tbody tr th:only-of-type { vertical-align: middle; } <pre data-sourcepos="1024:5-1030:5"><code>.dataframe tbody tr th { vertical-align: top; } .dataframe thead th { text-align: right; } </code></pre>

	B&W	canal R	canal G	canal B
promedio	3.432859	3.355606	3.436947	3.518292
media	2.260081	2.234073	2.261304	2.295926
desv estándar	2.369830	2.180836	2.395190	2.592449

Podemos ver un pequeño cambio en las medidas estadísticas, todos los promedios aumentaron al tener en cuenta la incertidumbre, así como la desviación estándar, la media disminuyó en todos los casos.

Una vez realizado todo este procedimiento, podemos decir que la resolución espacial de la máquina utilizada para capturar la imágen es de 3.4 pixeles.