se añade el archivo .md

ENTREGA.md

+## Juan David Hernández ejercicio1-clase05
+
+#### La idea de este programa, es determinar la resolución espacial del cielo en Zapatoca (Santander). La resolución espacial, hace referencia a la distancia angular mínima que debe haber entre dos fuentes de luz puntuales tal que se puedan seguir reconociendo de forma individual.
+
+Primero importamos las librerías necesarias. Requerimos de numpy para interpretar las imágenes como matrices y operar con ellas. Matplotlib para graficar, scipy para aplicar método de mínimos cuadrados de forma optimizada y statistics para utilizar funciones de estadística.
+
+
+```python
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib
+import scipy
+from scipy.optimize import leastsq
+import statistics
+import math
+from tabulate import tabulate
+import pandas as pd
+```
+
+Primero importamos la imágen de la carpeta data y la convertimos en un array de numpy.
+
+
+```python
+estrella=plt.imread('data/zapatocaImage.jpeg')
+```
+
+
+```python
+estrella1=np.array(estrella)
+```
+
+Es importante "plotear" y utilizar la función .shape en la imagen para tener una idea de los colores y el tamaño del objeto que vamos a analizar.
+
+
+```python
+plt.imshow(estrella)
+```
+
+
+
+
+    <matplotlib.image.AxesImage at 0x7f79a0cfae80>
+
+
+
+
+    
+![png](output_8_1.png)
+    
+
+
+
+```python
+estrella1.shape
+```
+
+
+
+
+    (789, 1184, 3)
+
+
+
+Como la imagen está en formato RGB, es una superposición de 3 matrices las cuales podemos extraer de forma individual. Para el procesamiento de la imagen, queremos que esta está a escala de grises, es por esto que vamos a aplicar un mapa de colores 'gray' para entender cómo se ve cada matriz en esta tonalidad.
+
+
+```python
+R=estrella1[:,:,0]
+plt.imshow(R,cmap='gray')
+```
+
+
+
+
+    <matplotlib.image.AxesImage at 0x7f79a0c365f8>
+
+
+
+
+    
+![png](output_11_1.png)
+    
+
+
+
+```python
+G=estrella1[:,:,1]
+plt.imshow(G,cmap='gray')
+```
+
+
+
+
+    <matplotlib.image.AxesImage at 0x7f79a0bac1d0>
+
+
+
+
+    
+![png](output_12_1.png)
+    
+
+
+
+```python
+B=estrella1[:,:,2]
+plt.imshow(B,cmap='gray')
+```
+
+
+
+
+    <matplotlib.image.AxesImage at 0x7f79a0b907b8>
+
+
+
+
+    
+![png](output_13_1.png)
+    
+
+
+Para obtener la imagen inicial en escala de grises podemos pensar en promediar las matrices RGB y plotear el resultado con cmap='gray'.
+
+
+```python
+grisprom=(R+G+B)/3
+plt.imshow(grisprom,cmap='gray')
+```
+
+
+
+
+    <matplotlib.image.AxesImage at 0x7f79a0afbb00>
+
+
+
+
+    
+![png](output_15_1.png)
+    
+
+
+Claramente se nota una saturación en píxeles que deberían tener una tonalidad más oscura, pero, ¿qué ha ocurrido? Si nos fijamos en la fórmula que utilizamos para promediar las matrices, habrá puntos donde $$R+G+B>255$$ y como estamos ploteando el resultado, no pueden haber valores que sobrepasen 255, lo que genera que python convierta cualquier valor mayor a 255 en 255, generando así este efecto.
+
+Para darnos cuenta que la explicación anterior es verdadera corremos la siguiente celda.
+
+
+```python
+grissum=(R+G+B)
+grissum.max()
+```
+
+
+
+
+    255
+
+
+
+Efectivamente ningún valor de la matriz supera el valor 255. Es por esto que la imagen no se ve como deseamos. Para tomar el promedio debemos entonces hacerlo de la forma $$\frac{1}{3}R+\frac{1}{3}G+\frac{1}{3}B$$
+
+
+```python
+grisprom=(1/3)*R+(1/3)*G+(1/3)*B
+plt.imshow(grisprom,cmap='gray')
+```
+
+
+
+
+    <matplotlib.image.AxesImage at 0x7f79a0a684e0>
+
+
+
+
+    
+![png](output_20_1.png)
+    
+
+
+Aunque cumplimos con nuestro objetivo, vemos que hay partes donde las estrellas contrastan menos respecto al cielo en comparación con la imagen inicial. Pensamos entonces en hacer un promedio ponderado para aumentar el contraste.
+
+Nos damos cuenta que la matriz R posee más píxeles blancos que las otras en zonas donde no hay estrellas, ocultando así información sobre la posición de estas. Podemos pensar entonces que para generar la imagen en escala de grises, le vamos a dar menos peso a la matriz R para generar mayor contraste entre las estrellas y el cielo.
+
+
+```python
+gray1=(R*(1/10))+(G*(9/20))+(B*(9/20))
+gray1.mean()
+```
+
+
+
+
+    56.30558824033155
+
+
+
+Obtenemos entonces nuestra imagen en escala de grises a partir de un promedio ponderado
+
+
+```python
+plt.imshow(gray1,cmap='gray')
+```
+
+
+
+
+    <matplotlib.image.AxesImage at 0x7f79a0a4f898>
+
+
+
+
+    
+![png](output_25_1.png)
+    
+
+
+Queremos estudiar estrellas de forma individual, para esto seleccionamos secciones de la imagen para que sea más fácil dividir una parte de la imagen que contenga una única estrella.
+
+
+```python
+seccion1=gray1[300:400,450:650]
+seccion2=gray1[100:200,250:450]
+```
+
+Ahora mostramos las dos secciones que hemos escogido para el análisis. Las ploteamos en una malla de subplots para poder imprimir ambas imágenes en la misma celda.
+
+
+```python
+f, axarr = plt.subplots(2)
+plt.setp(axarr, xticks=np.arange(0,200,25),yticks=np.arange(0,100,20))
+axarr[0].imshow(seccion1,cmap='gray')
+axarr[1].imshow(seccion2,cmap='gray')
+```
+
+
+
+
+    <matplotlib.image.AxesImage at 0x7f799ceade10>
+
+
+
+
+    
+![png](output_29_1.png)
+    
+
+
+Ahora recortamos estrellas individuales de ambas secciones
+
+
+```python
+star1=seccion1[58:70,95:107]
+star2=seccion1[34:46,85:97]
+star3=seccion1[0:12,23:35]
+star4=seccion1[20:32,109:121]
+star5=seccion2[38:50,56:68]
+star6=seccion2[12:24,108:120]
+star7=seccion2[33:45,181:193]
+star8=seccion2[78:90,104:116]
+```
+
+
+```python
+estrellas=[star1,star2,star3,star4,star5,star6,star7,star8]
+```
+
+En total se recortaron 8 estrellas en imágenes de 12x12(144) píxeles.
+
+
+```python
+f, axarr = plt.subplots(2,4)
+axarr[0,0].imshow(star1,cmap='gray')
+axarr[0,1].imshow(star2,cmap='gray')
+axarr[0,2].imshow(star3,cmap='gray')
+axarr[0,3].imshow(star4,cmap='gray')
+axarr[1,0].imshow(star5,cmap='gray')
+axarr[1,1].imshow(star6,cmap='gray')
+axarr[1,2].imshow(star7,cmap='gray')
+axarr[1,3].imshow(star8,cmap='gray')
+```
+
+
+
+
+    <matplotlib.image.AxesImage at 0x7f799cccbda0>
+
+
+
+
+    
+![png](output_34_1.png)
+    
+
+
+Para determinar la resolución espacial, es necesario ajustar cada imagen a una función gaussiana 2D. Puesto que la resolución espacial se toma como el [FWHM](https://es.wikipedia.org/wiki/Anchura_a_media_altura#:~:text=La%20Anchura%20a%20media%20altura,mitad%20de%20su%20valor%20m%C3%A1ximo.) de la gaussiana que característica de la imagen. Donde FWHM viene dado por $$FWHM=2\sqrt{2\ln {2}}\sigma$$ siendo sigma la desviación estándar asociada a la gaussiana.
+
+Para este problema, ajustaremos por el método de mínimos cuadrados una función gaussiana 2D a cada una de las estrellas para determinar sigma y su FWHM correspondiente. Sabemos que una función gaussiana 2D es de la forma $$a\exp{\bigg(\frac{(x-x_{0})^{2}+(y-y_{0})^{2}}{2\sigma^2}\bigg)}+b.$$ Así que las constantes a determinar son $$a,b,x_{0},y_{0},\sigma$$. Donde b representa una constante aditiva que hace referencia al valor de los píxeles por fuera de la estrella.
+
+Procedemos a definir la función gaussiana y una función error la cual debe ser minimizada. En ambos casos añadiremos la función .ravel() para obtener una arreglo 1D como resultado de las funciones en lugar de uno 2D. Esto para facilitar el uso de la función optimize.leastsq.
+
+
+```python
+def gauss2D(params,x, y):
+    exponente = -((x-params[2])**2 + (y-params[3])**2) / (2*params[4]**2)
+    z = (params[0]*np.exp(exponente)) + params[1]
+    return z.ravel()
+```
+
+
+```python
+def error(tpl,x,y,z):
+    return (gauss2D(tpl,x,y)-z).ravel()
+```
+
+También es necesario aplicar reshape a las imágenes para volverlas arreglos 1D, para esto se utiliza la función reshape.
+
+
+```python
+estrellasshape=[np.reshape(star1,144),np.reshape(star2,144),np.reshape(star3,144),np.reshape(star4,144)
+               ,np.reshape(star5,144),np.reshape(star6,144),np.reshape(star7,144),np.reshape(star8,144)]
+```
+
+Creamos dos arreglos que tengan de longitud el eje x y y de las imágenes de las estrellas.
+
+
+```python
+x=np.arange(0,12,1)
+y=np.arange(0,12,1)
+```
+
+
+```python
+yy, xx = np.meshgrid(x,y)
+```
+
+Para el procedimiento de mínimos cuadrados, hay que proporcionar valores iniciales de las 5 constantes a determinar. La constante aditiva (b) tendrá un valor cercano al mínimo que toma la imagen (cielo apagado) ya que es el valor que toma la función cuando x y y son grandes (i.e por fuera de la estrella). La amplitud, será cercana a la diferencia entre el máximo (cerca al centro de la estrella) y el mínimo de la imagen. x_0 y y_0 serán cercanos al centro de la imagen y sigma a un cuarto de la longitud de la imagen. Entonces $$a=np.max(estrella)-np.min(estrella), b=np.min(estrella), x_0=\frac{len(estrella)}{2}, y_0=\frac{len(estrella)}{2}, \sigma^2=\frac{len(estrella)}{4}.$$
+
+
+```python
+p=[]
+a=0
+while a < len(estrellas):
+    p.append([np.max(estrellas[a])-np.min(estrellas[a]),np.min(estrellas[a]),
+              len(estrellas[a][0,:])/2,len(estrellas[a][:,0])/2,len(estrellas[a][0,:])/4])
+    a+=1
+```
+
+Ahora, aplicamos la función leastsq sobre las 8 imágenes para minimizar la función error.
+
+
+```python
+b=0
+best=[]
+while b < len(estrellas):
+    best.append(leastsq(error, p[b], args=(xx,yy,estrellasshape[b])))
+    b+=1
+```
+
+Podemos comprobar si el método que utilizamos fue eficiente graficando la superficie de la primera estrella y la gaussiana con los parámetros encontrados y haciendo una comparación. Para esto necesitamos la función gaussiana pero que esta vez retorne un array 2D.
+
+
+```python
+def gauss2Dsup(params,x, y):
+    exponente = -((x-params[2])**2 + (y-params[3])**2) / (2*params[4]**2)
+    z = (params[0]*np.exp(exponente)) + params[1]
+    return z
+```
+
+Ploteamos entonces la función gaussiana con los parámetros "best1"
+
+
+```python
+zz2=gauss2Dsup(best[0][0],xx,yy)
+fig2 = plt.figure()
+ax = fig2.gca(projection='3d')
+ax.plot_surface(xx,yy, zz2, rstride=3, cstride=3, linewidth=1, antialiased=True,cmap="rainbow")
+```
+
+
+
+
+    <mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7f79a0a15c18>
+
+
+
+
+    
+![png](output_52_1.png)
+    
+
+
+También graficamos la matriz que representa a la estrella 1 como una superficie.
+
+
+```python
+fig3 = plt.figure()
+ax = fig3.gca(projection='3d')
+ax.plot_surface(xx,yy, estrellas[0], rstride=3, cstride=3, linewidth=1, antialiased=True,cmap="rainbow")
+```
+
+
+
+
+    <mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7f799cbe12e8>
+
+
+
+
+    
+![png](output_54_1.png)
+    
+
+
+Ahora miramos de forma gráfica la diferencia entre ambas superficies, es decir el error.
+
+
+```python
+def error2(tpl,x,y,z):
+    return gauss2Dsup(tpl,x,y)-z
+```
+
+
+```python
+plt.imshow(error2(best[0][0],xx,yy,estrellas[0]))
+plt.colorbar()
+```
+
+
+
+
+    <matplotlib.colorbar.Colorbar at 0x7f799cac6748>
+
+
+
+
+    
+![png](output_57_1.png)
+    
+
+
+
+```python
+fig4 = plt.figure()
+ax = fig4.gca(projection='3d')
+ax.plot_surface(xx,yy, error2(best[0][0],xx,yy,estrellas[0]), rstride=3, cstride=3, linewidth=1, antialiased=True,cmap="rainbow")
+```
+
+
+
+
+    <mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7f799caa4048>
+
+
+
+
+    
+![png](output_58_1.png)
+    
+
+
+Y calculamos su promedio y desviación estándar.
+
+
+```python
+error2(best[0][0],xx,yy,estrellas[0]).mean()
+```
+
+
+
+
+    1.0861282243960584e-09
+
+
+
+
+```python
+statistics.stdev(error2(best[0][0],xx,yy,estrellas[0]).reshape(144))
+```
+
+
+
+
+    8.691719728230568
+
+
+
+Como se observa, los errores no sobrepasan el 10% del valor de las funciones, podemos entonces afirmar que se encontró una buena aproximación.
+
+### Cálculo del FWHM
+
+Creamos una lista con las raíces de $$\sigma ^{2}$$ para el cálculo del FWHM y la volvemos un array de numpy para poder multiplicarla por floats.
+
+
+```python
+sigmas=[]
+c=0
+while c < len(estrellas):
+    sigmas.append(math.sqrt(abs(best[c][0][4])))
+    c+=1
+```
+
+
+```python
+sigmasnp=np.array(sigmas)
+```
+
+Ahora aplicamos la fórmula de FWHM sobre el array y obtenemos la resolución espacial de cada estrella.
+
+
+```python
+FWHM=2*math.sqrt(2*math.log(2))*sigmasnp
+```
+
+Finalmente sacamos el promedio para obtener un valor aproximado de la resolución espacial en Zapatoca (Santander).
+
+
+```python
+FWHM.mean()
+```
+
+
+
+
+    2.0974170630756683
+
+
+
+### Ahora vamos a repetir el cálculo del FHWM pero utilizando las matrices R,G y B para ver si existe alguna diferencia.
+
+### Matriz R
+
+Seleccionamos las dos secciones en la matriz R.
+
+
+```python
+seccion1R=R[300:400,450:650]
+seccion2R=R[100:200,250:450]
+```
+
+Ahora recortamos las estrellas individuales
+
+
+```python
+star1R=seccion1R[58:70,95:107]
+star2R=seccion1R[34:46,85:97]
+star3R=seccion1R[0:12,23:35]
+star4R=seccion1R[20:32,109:121]
+star5R=seccion2R[38:50,56:68]
+star6R=seccion2R[12:24,108:120]
+star7R=seccion2R[33:45,181:193]
+star8R=seccion2R[78:90,104:116]
+estrellasR=[star1R,star2R,star3R,star4R,star5R,star6R,star7R,star8R]
+```
+
+Convertimos las imágenes en arreglos 1D.
+
+
+```python
+estrellasshapeR=[np.reshape(star1R,144),np.reshape(star2R,144),np.reshape(star3R,144),np.reshape(star4R,144)
+               ,np.reshape(star5R,144),np.reshape(star6R,144),np.reshape(star7R,144),np.reshape(star8R,144)]
+```
+
+Creamos la lista con los valores iniciales.
+
+
+```python
+pR=[]
+d=0
+while d < len(estrellas):
+    pR.append([np.max(estrellasR[d])-np.min(estrellasR[d]),
+               np.min(estrellasR[d]),len(estrellasR[d][0,:])/2,len(estrellasR[d][:,0])/2,len(estrellasR[d][0,:])/4])
+    d+=1
+```
+
+Aplicamos la función leastsq utilizando los valores iniciales y la función error.
+
+
+```python
+e=0
+bestR=[]
+while e < len(estrellasR):
+    bestR.append(leastsq(error, pR[e], args=(xx,yy,estrellasshapeR[e])))
+    e+=1
+```
+
+Como obtenemos el cuadrado de sigma, debemos tomar su raíz para calcular el FWHMR.
+
+
+```python
+sigmasR=[]
+f=0
+while f < len(estrellasR):
+    sigmasR.append(math.sqrt(abs(bestR[f][0][4])))
+    f+=1
+```
+
+Ahora aplicampos la ecuación de FWHMR y calculamos su promedio.
+
+
+```python
+sigmasnpR=np.array(sigmasR)
+FWHMR=2*math.sqrt(2*math.log(2))*sigmasnpR
+FWHMR.mean()
+```
+
+
+
+
+    2.3115564774230157
+
+
+
+Repetimos el proceso para la matriz G y B.
+
+## Matriz G
+
+
+```python
+seccion1G=G[300:400,450:650]
+seccion2G=G[100:200,250:450]
+```
+
+
+```python
+star1G=seccion1G[58:70,95:107]
+star2G=seccion1G[34:46,85:97]
+star3G=seccion1G[0:12,23:35]
+star4G=seccion1G[20:32,109:121]
+star5G=seccion2G[38:50,56:68]
+star6G=seccion2G[12:24,108:120]
+star7G=seccion2G[33:45,181:193]
+star8G=seccion2G[78:90,104:116]
+estrellasG=[star1G,star2G,star3G,star4G,star5G,star6G,star7G,star8G]
+```
+
+
+```python
+estrellasshapeG=[np.reshape(star1G,144),np.reshape(star2G,144),np.reshape(star3G,144),np.reshape(star4G,144)
+               ,np.reshape(star5G,144),np.reshape(star6G,144),np.reshape(star7G,144),np.reshape(star8G,144)]
+```
+
+
+```python
+pG=[]
+g=0
+while g < len(estrellasG):
+    pG.append([np.max(estrellasG[g])-np.min(estrellasG[g]),np.min(estrellasG[g]),
+               len(estrellasG[g][0,:])/2,len(estrellasG[g][:,0])/2,len(estrellasG[g][0,:])/4])
+    g+=1
+```
+
+
+```python
+bestG=[]
+h=0
+while h < len(estrellasG):
+    bestG.append(leastsq(error, pG[h], args=(xx,yy,estrellasshapeG[h])))
+    h+=1
+```
+
+
+```python
+sigmasG=[]
+i=0
+while i < len(estrellasG):
+    sigmasG.append(math.sqrt(abs(bestG[i][0][4])))
+    i+=1
+```
+
+
+```python
+sigmasnpG=np.array(sigmasG)
+FWHMG=2*math.sqrt(2*math.log(2))*sigmasnpG
+FWHMG.mean()
+```
+
+
+
+
+    2.320331194403334
+
+
+
+## Matriz B
+
+
+```python
+seccion1B=B[300:400,450:650]
+seccion2B=B[100:200,250:450]
+```
+
+
+```python
+star1B=seccion1B[58:70,95:107]
+star2B=seccion1B[34:46,85:97]
+star3B=seccion1B[0:12,23:35]
+star4B=seccion1B[20:32,109:121]
+star5B=seccion2B[38:50,56:68]
+star6B=seccion2B[12:24,108:120]
+star7B=seccion2B[33:45,181:193]
+star8B=seccion2B[78:90,104:116]
+estrellasB=[star1B,star2B,star3B,star4B,star5B,star6B,star7B,star8B]
+```
+
+
+```python
+estrellasshapeB=[np.reshape(star1B,144),np.reshape(star2B,144),np.reshape(star3B,144),np.reshape(star4B,144)
+               ,np.reshape(star5B,144),np.reshape(star6B,144),np.reshape(star7B,144),np.reshape(star8B,144)]
+```
+
+
+```python
+pB=[]
+j=0
+while j < len(estrellasB):
+    pB.append([np.max(estrellasB[j])-np.min(estrellasB[j]),
+               np.min(estrellasB[j]),len(estrellasB[j][0,:])/2,len(estrellasB[j][:,0])/2,len(estrellasB[j][0,:])/4])
+    j+=1
+```
+
+
+```python
+bestB=[]
+k=0
+while k < len(estrellasB):
+    bestB.append(leastsq(error, pB[k], args=(xx,yy,estrellasshapeB[k])))
+    k+=1
+```
+
+
+```python
+sigmasB=[]
+l=0
+while l < len(estrellasB):
+    sigmasB.append(math.sqrt(abs(bestB[l][0][4])))
+    l+=1
+```
+
+
+```python
+sigmasnpB=np.array(sigmasB)
+FWHMB=2*math.sqrt(2*math.log(2))*sigmasnpB
+FWHMB.mean()
+```
+
+
+
+
+    2.38171347904209
+
+
+
+### Resumiendo los valores de FWHM en las 4 matrices, obtenemos.
+
+
+```python
+df = pd.DataFrame({'Escala de Grises' : [FWHM.mean()],
+                   'R' : [FWHMR.mean()],
+                  'G' : [FWHMG.mean()],
+                 'B' : [FWHMB.mean()]})
+print(tabulate(df, headers='keys', tablefmt='psql',showindex="never"))
+```
+
+    +--------------------+---------+---------+---------+
+    |   Escala de Grises |       R |       G |       B |
+    |--------------------+---------+---------+---------|
+    |            2.09742 | 2.31156 | 2.32033 | 2.38171 |
+    +--------------------+---------+---------+---------+
+
+
+Se encuentra que el FWHM para el caso de la imagen en escala de grises es menor que al tomar cada una de las matrices de forma independiente.