diff --git a/ENTREGA.md b/ENTREGA.md new file mode 100644 index 0000000000000000000000000000000000000000..80ecc8f9e0e60fd41d0c4dd9b6c033352eef9e02 --- /dev/null +++ b/ENTREGA.md @@ -0,0 +1,773 @@ +## 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> + + + + + + + + + + +```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> + + + + + + + + + + +```python +G=estrella1[:,:,1] +plt.imshow(G,cmap='gray') +``` + + + + + <matplotlib.image.AxesImage at 0x7f79a0bac1d0> + + + + + + + + + + +```python +B=estrella1[:,:,2] +plt.imshow(B,cmap='gray') +``` + + + + + <matplotlib.image.AxesImage at 0x7f79a0b907b8> + + + + + + + + + +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> + + + + + + + + + +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> + + + + + + + + + +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> + + + + + + + + + +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> + + + + + + + + + +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> + + + + + + + + + +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> + + + + + + + + + +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> + + + + + + + + + +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> + + + + + + + + + + +```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> + + + + + + + + + +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. diff --git a/output_11_1.png b/output_11_1.png new file mode 100644 index 0000000000000000000000000000000000000000..c9ff4d5d72ad598de4eb4fe9fd91d28b62e8fa88 Binary files /dev/null and b/output_11_1.png differ diff --git a/output_12_1.png b/output_12_1.png new file mode 100644 index 0000000000000000000000000000000000000000..21e22cd013d7431a6c0b2bc0792bf43b1fb0e30b Binary files /dev/null and b/output_12_1.png differ diff --git a/output_13_1.png b/output_13_1.png new file mode 100644 index 0000000000000000000000000000000000000000..de7c22121bd81a8af4d51e0bc5adde39d2d3e169 Binary files /dev/null and b/output_13_1.png differ diff --git a/output_15_1.png b/output_15_1.png new file mode 100644 index 0000000000000000000000000000000000000000..836bd689112d5727c85196f8235a27f99c896f04 Binary files /dev/null and b/output_15_1.png differ diff --git a/output_20_1.png b/output_20_1.png new file mode 100644 index 0000000000000000000000000000000000000000..56c6d2764f057c07607bce84d324c4f0e92e4516 Binary files /dev/null and b/output_20_1.png differ diff --git a/output_25_1.png b/output_25_1.png new file mode 100644 index 0000000000000000000000000000000000000000..820b0241bbf421aa9b0c74f515a001ddac6ae396 Binary files /dev/null and b/output_25_1.png differ diff --git a/output_29_1.png b/output_29_1.png new file mode 100644 index 0000000000000000000000000000000000000000..a7a88608eda67cf6c68ccf80aa7c32c2d8153cb5 Binary files /dev/null and b/output_29_1.png differ diff --git a/output_34_1.png b/output_34_1.png new file mode 100644 index 0000000000000000000000000000000000000000..7bdfafbac87402921c15add506ad1d20358e76fb Binary files /dev/null and b/output_34_1.png differ diff --git a/output_52_1.png b/output_52_1.png new file mode 100644 index 0000000000000000000000000000000000000000..c81b7209ddb5794a9c7f380107f84246ed5d7ffd Binary files /dev/null and b/output_52_1.png differ diff --git a/output_54_1.png b/output_54_1.png new file mode 100644 index 0000000000000000000000000000000000000000..440d82bfac20df71ed9521f32b8f120c31b01acb Binary files /dev/null and b/output_54_1.png differ diff --git a/output_57_1.png b/output_57_1.png new file mode 100644 index 0000000000000000000000000000000000000000..183e1235ec775222955b38c452d165a9b06d13f2 Binary files /dev/null and b/output_57_1.png differ diff --git a/output_58_1.png b/output_58_1.png new file mode 100644 index 0000000000000000000000000000000000000000..77a4f7dc00ba2570a4483168accd9705704a5993 Binary files /dev/null and b/output_58_1.png differ diff --git a/output_8_1.png b/output_8_1.png new file mode 100644 index 0000000000000000000000000000000000000000..c15b76aa0b3a254e8f60902b7296bf8b72d8e749 Binary files /dev/null and b/output_8_1.png differ