diff --git a/Ejercicios-clase-05-datosF.md b/Ejercicios-clase-05-datosF.md new file mode 100644 index 0000000000000000000000000000000000000000..b0d1a80db746e445fa22119568b6353bb0ce3bea --- /dev/null +++ b/Ejercicios-clase-05-datosF.md @@ -0,0 +1,729 @@ +## Alumno: Christian Solis Calero. +### Universidad Nacional Mayor de San Marcos, Lima. Perú + +### Ejercicio 05-01: 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ósfera como la convolución de la imagen "perfecta" (como se verÃa desde el espacio) +con un kernel gaussiano. El ancho de esa función gaussiana 2D caracteriza las condiciones de observación, varÃa con las condiciones climáticas y para cada sitio de la Tierra. +La resolución espacial normalmente se toma como 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 +1. 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 + +2. 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) + +3. 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) + +4. 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 + +5. Repita el mismo ejercicio sobre cada una de las bandas R,G,B separadamente y comente si observa diferencias en los resultados + +### 1. Manipulando la imagen + + +```python +import numpy as np +import matplotlib.pyplot as plt +from matplotlib import image as img +%matplotlib inline + +# load the image +imagen = img.imread('data/zapatocaImage.jpeg') +# convert image to numpy array +imagen_array = np.asarray(imagen) +plt.imshow(imagen_array) +plt.grid(None) +plt.show() +``` + + + + + + + + +```python +#In RGB (Red-Green-Blue) channels, the intensity of the color is presented on a 0–255 scale, at first we rescaled to the [0,1] range. +imagen_array = imagen_array/255 + +#Creating a function for generate gray image, combinig channels, multiplying the red matrix by 0.2989, +#the green matrix by 0.587 and the blue matrix by 0.114 +#and summing: 0.299 R + 0.587 G + 0.114 B to get each pixel of the grayscale image + +def gray_generator(rgb): + ''' + Define function for generating gray image fromo RGB data. + ''' + r, g, b = rgb[:,:,0], rgb[:,:,1], rgb[:,:,2] + gray = 0.2989 * r + 0.5870 * g + 0.1140 * b + return gray + +imagen_gris=gray_generator(imagen_array) + +plt.imshow(imagen_gris, cmap = 'gray') +plt.grid(None) +plt.show() +``` + + + + + + + +### 2. Procesando la información de la primera estrella seleccionada +Generación de arrays conteniendo las variables de posisción (x e y) y de la intensidad de pixeles (z) + + +```python +#Working with a FIRST STAR +#based in coordinates of gray image crop the stars +star01 = imagen_gris[540:565,650:675] +plt.imshow(star01, cmap ='gray') +plt.grid(None) +plt.show() + +``` + + + + + + + + +```python +#convert its negative image: +y=np.shape(star01) +negative_star01=np.zeros(y) +#convert the image into its negative value. +negative_star01=1-star01 +plt.xlabel("Value") +plt.ylabel("pixels Frequency") +plt.title("Negative image ") +plt.imshow(negative_star01,cmap="gray") +plt.grid(None) +plt.show() +``` + + + + + + + +Usando la información de los arrays x, y, z visualizamos la data a través de plots en 3D + + +```python +# Generate a 3D Plot +# get coordinates (y,x) --- alternately see below for (x,y) +yx_coords = np.column_stack(np.where(star01 >= 0)) + +#Generating vectors of x and y coordinates from x,y array +y=[item[0] for item in yx_coords] +x=[item[1] for item in yx_coords] + +#Generating a vector from matrix of pixeles +zcoords = star01.flatten() +``` + + +```python +#Creating plot with the data +# Import libraries +from mpl_toolkits import mplot3d +import numpy as np +import matplotlib.pyplot as plt + +# Creating dataset +z =zcoords + +# Creating figure +fig = plt.figure(figsize = (10, 7)) +ax = plt.axes(projection ="3d") + +# Creating plot +ax.scatter3D(x, y, z, color = "green") +plt.title("simple 3D scatter plot") + + +# show plot +plt.show() + +``` + + + + + + + +### Distribución Gaussiana +Trabajamos primero considerando solo los pixeles de la lÃnea que pasa por la mitad de la estrella, obteniendo un vector de valores de intensidad luminosa solo para esa linea. Lo graficamos y lo ajustamos a una una función gaussiana 1D + + + +```python +#tomar solo los pixeles de la lÃnea que pasa por la mitad de la estrella +#Valores de intensidad de pixeles correpondiente al valor 12 del array X +import numpy as np +import matplotlib.pyplot as plt +plt.style.use('seaborn-whitegrid') +%matplotlib inline + +Zlinea=[] +Ylinea=[] +j=0 +for i in range(z.size-1): + if x[i]==12: + Zlinea =np.insert(Zlinea, j, z[i]) + j=j+1 +j=0 +for i in range(z.size-1): + if x[i]==12: + Ylinea =np.insert(Ylinea, j, y[i]) + j=j+1 + +plt.plot(Ylinea, Zlinea, 'o', color='blue'); +plt.xlabel('Position in image') +plt.ylabel('Intensity of pixels') + + +``` + + + + + Text(0, 0.5, 'Intensity of pixels') + + + + + + + + + + +```python +import numpy as np +import matplotlib.pyplot as plt +from numpy import asarray as ar, exp, sqrt +from scipy.optimize import curve_fit + +YLa = ar(Ylinea) +ZLa = ar(Zlinea) + +n = len(ZLa) ## <--- +mean = sum(ZLa*YLa)/n +sigma = sqrt(sum(ZLa*(YLa-mean)**2)/n) + +def gaus(x,a,mu,sigma): + return a*exp(-(x-mu)**2/(2*sigma**2)) + +popt,pcov = curve_fit(gaus,YLa,ZLa) # first estimation of the parameters +print('Fitted parameters:') +print(popt) + +xx = np.linspace( 0, 25, 100 ) ## <--- calculate against a continuous variable + +fig = plt.figure() +plt.plot(YLa, ZLa, "ob", label = "Measured") +plt.plot(xx,gaus(xx,*popt),'r',label='Fit') ## <--- plot against the contious variable +plt.xlim(0, 25) +plt.ylim(0.4, 1.1) +plt.xticks(YLa) +plt.title("Fitting gaussian on pixel data of star image") +plt.xlabel("Position in image") +plt.ylabel("Intensity of pixels") +plt.grid() +plt.legend() +plt.savefig('Fitting2D.png') +plt.show() +``` + + Fitted parameters: + [ 0.97409685 12.68315045 8.98659285] + + + + + + + + + +```python +#Generating an histogram of the pixel values +plt.hist(Zlinea, bins = 'auto') +plt.xlabel('Intervals of Intensity of pixels') +plt.ylabel('Number of pixels') +plt.show() +``` + + + + + + + +### Distribución Gaussiana multivariante +Trabajamos ahora considerando todos los pixeles de la imagen de la estrella, usando los valores de los vectores de posición (x, y) y de intensidad lumninosa: pixeles (z). Lo graficamos y lo ajustamos a una una función gaussiana 2D. + + + +```python +import numpy as np +from scipy.optimize import curve_fit +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +from numpy import asarray as ar, exp, sqrt + +# The two-dimensional domain of the fit. +xmin, xmax, nx = np.min(x), np.max(x), np.size(x) +ymin, ymax, ny = np.min(y), np.max(y), np.size(y) +#xTa, yTa = np.linspace(xmin, xmax, nx), np.linspace(ymin, ymax, ny) +xTa = ar(x) +yTa = ar(y) + +XTa, YTa = np.meshgrid(xTa, yTa) + +#sigma_YTa = sqrt(sum(z*(yTa-12)**2)/(np.size) +sigma_YTa = sqrt(sum(z*(yTa-12)**2)/625) + + +# Our function to fit is a two-dimensional Gaussians +def gaussian(x, y, x0, y0, xalpha, yalpha, A): + return A * np.exp( -0.5*((x-x0)/xalpha)**2 -0.5*((y-y0)/yalpha)**2) + +# A list of the Gaussian parameters: x0, y0, xalpha, yalpha, A + +gprms=[12,12,sigma_YTa,sigma_YTa,1] + +# Standard deviation of normally-distributed noise to add in generating +# our test function to fit. +noise_sigma = 0.1 + +# The function to be fit is Z. +ZTa = gaussian(XTa, YTa, 12,12,sigma_YTa,sigma_YTa,1) + +# Plot the 3D figure of the gaussian function. +fig = plt.figure() +ax = fig.gca(projection='3d') +ax.plot_surface(XTa, YTa, ZTa, cmap=cm.viridis,rstride=3, cstride=3, linewidth=1, antialiased=True) +ax.set_zlim(0,np.max(ZTa)) + +cset = ax.contourf(XTa, YTa, ZTa, zdir='z', offset=-0.95, cmap=cm.viridis) + +# Adjust the limits, ticks and view angle +ax.set_zlim(-0.95,1.0) +ax.set_zticks(np.linspace(0,1,5)) +ax.view_init(27, -21) + +plt.show() +``` + + + + + + + +Usamos los datos de intensidad luminosa generados por la función gaussiana, asi como los datos de posición (x e y) para ajustar la gaussiana a los datos iniciales de intensidad luminosa (z) +Usamos el ajuste no lineal leastsq de scipy, y obtenemos imagnese reconstruidas de la estrella sin y con ajuste + + +```python +from scipy.optimize import leastsq +#print(ZTa) +#print(np.size(ZTa)) +ZTr = [z] +for i in range(z.size-1): + ZTr = np.append(ZTr, [z], 0) + +#print(ZTr) +#print(np.size(ZTr)) + +parameters=[12,12,sigma_YTa,sigma_YTa,1] + +def gaussian3D(params,x, y): + #z = params[4] * np.exp( -((x-params[0])/params[2])**2 -((y-params[1])/params[3])**2) + z = params[4] * np.exp( -0.5*((x-params[0])/(params[2]))**2 -0.5*((y-params[1])/(params[3]))**2) + return z +zOriginal=gaussian3D(parameters,xTa,yTa) + + +def error3D(params,x,y,z, cols, rows): + errors = (gaussian3D(params,x,y)-z) + return errors.reshape(cols*rows) + +obtained_error=error3D(parameters,xTa,yTa, z, 25, 25) + +print('The average error between gaussian model and data is:') +print(np.average(obtained_error)) + +best3D,sus3D = leastsq(error3D, parameters, args = (xTa,yTa, z, 25, 25)) + +print('Fitted parameters:') +print(best3D) +Fitted_parameters=[12.31,12.31,16.59,16.59,0.84] + +#Using the fitted parameters for obtaining pixel values +zFitted=gaussian3D(best3D,xTa,yTa) + +#Showing the imagen of star with original values, values according to gaussian function and values according to fitted parameters +plt.subplot(1, 3, 1) +star01_01 = z.reshape(25, 25) +plt.title('original') +plt.imshow(star01_01,cmap="gray") +plt.grid(None) +plt.subplot(1, 3, 2) +star01_02 = zOriginal.reshape(25, 25) +plt.title('Gaussian') +plt.imshow(star01_02,cmap="gray") +plt.grid(None) +plt.subplot(1, 3, 3) +star01_03 = zFitted.reshape(25, 25) +plt.title('Fitted') +plt.imshow(star01_03,cmap="gray") +plt.grid(None) +plt.show() + +#Calculando algunos parámetros estadÃsticos para los 3 casos +print("The original data has an average of: ",np.average(z)) +print("The original data has an standard deviation of: ",np.std(z)) + +print("The data generated directly from gaussian model has an average of: ",np.average(zOriginal)) +print("The data generated directly from gaussian model has an standard deviation of: ",np.std(zOriginal)) + +print("The fitted data has an average of: ",np.average(zFitted)) +print("The fitted data has an standard deviation of: ",np.std(zFitted)) +``` + + The average error between gaussian model and data is: + -0.33308824487283 + Fitted parameters: + [12.31003544 12.64869953 11.72961613 11.65211515 0.84817359] + + + + + + + + + The original data has an average of: 0.5972459312941176 + The original data has an standard deviation of: 0.1470436454930245 + The data generated directly from gaussian model has an average of: 0.2641576864212876 + The data generated directly from gaussian model has an standard deviation of: 0.25792679949611685 + The fitted data has an average of: 0.5953058814848675 + The fitted data has an standard deviation of: 0.13721730657930556 + + +### 2. Procesando la información de la segunda estrella seleccionada +Trabajamos con los arrays conteniendo los valores para las variables de posición (x2 e y2) y de la intensidad de pixeles (z2) provenientes de la imagen de tercera estrella seleccionada. + + +```python +#Working with a SECOND STAR +#based in coordinates of gray image crop the stars +star02 = imagen_gris[307:326,618:637] +plt.imshow(star02, cmap ='gray') +plt.grid(None) +plt.show() +``` + + + + + + + + +```python +# Generate a 3D Plot +# Import libraries +from mpl_toolkits import mplot3d +import numpy as np +import matplotlib.pyplot as plt + +# get coordinates (y,x) --- alternately see below for (x,y) +yx_coords2 = np.column_stack(np.where(star02 >= 0)) + +#Generating vectors of x and y coordinates from x,y array +y2=[item[0] for item in yx_coords2] +x2=[item[1] for item in yx_coords2] + +#Generating a vector from matrix of pixeles +zcoords2 = star02.flatten() + +# Creating dataset +z2 =zcoords2 + +# Creating figure +fig = plt.figure(figsize = (10, 7)) +ax = plt.axes(projection ="3d") + +# Creating plot +ax.scatter3D(x2, y2, z2, color = "green") +plt.title("Simple 3D scatter plot for Star 2") + +# show plot +plt.show() + +``` + + + + + + + + +```python +from scipy.optimize import leastsq + +xTa2 = ar(x2) +yTa2 = ar(y2) + +#Calculating parameters +xTa2_average=np.average(xTa2) +yTa2_average=np.average(yTa2) + +sigma_yTa2 = sqrt(sum(z2*(yTa2-9)**2)/361) + +parameters2=[9,9,sigma_yTa2,sigma_yTa2,1] + +def gaussian3D(params,x, y): + z = params[4] * np.exp( -0.5*((x-params[0])/(params[2]))**2 -0.5*((y-params[1])/(params[3]))**2) + return z +zOriginal2=gaussian3D(parameters2,xTa2,yTa2) + + +def error3D(params,x,y,z, cols, rows): + errors = (gaussian3D(params,x,y)-z) + return errors.reshape(cols*rows) + +obtained_error2=error3D(parameters2,xTa2,yTa2, z2, 19, 19) + +print('The average error between gaussian model and data is:') +print(np.average(obtained_error2)) + + +best3D2,sus3D2 = leastsq(error3D, parameters2, args = (xTa2,yTa2, z2, 19, 19)) + +print('Fitted parameters:') +print(best3D2) + + +#Using the fitted parameters for obtaining pixel values +zFitted2=gaussian3D(best3D2,xTa2,yTa2) + +#Showing the imagen of star with original values, values according to gaussian function and values according to fitted parameters +plt.subplot(1, 3, 1) +star02_01 = z2.reshape(19, 19) +plt.title('original') +plt.imshow(star02_01,cmap="gray") +plt.grid(None) +plt.subplot(1, 3, 2) +star02_02 = zOriginal2.reshape(19, 19) +plt.title('Gaussian') +plt.imshow(star02_02,cmap="gray") +plt.grid(None) +plt.subplot(1, 3, 3) +star02_03 = zFitted2.reshape(19, 19) +plt.title('Fitted') +plt.imshow(star02_03,cmap="gray") +plt.grid(None) +plt.show() + +#Calculando algunos parámetros estadÃsticos para los 3 casos +print("The original data has an average of: ",np.average(z2)) +print("The original data has an standard deviation of: ",np.std(z2)) + +print("The data generated directly from gaussian model has an average of: ",np.average(zOriginal2)) +print("The data generated directly from gaussian model has an standard deviation of: ",np.std(zOriginal2)) + +print("The fitted data has an average of: ",np.average(zFitted2)) +print("The fitted data has an standard deviation of: ",np.std(zFitted2)) +``` + + The average error between gaussian model and data is: + -0.2504618299460217 + Fitted parameters: + [9.15491427 8.91137258 8.54230266 8.65676436 0.66111467] + + + + + + + + + The original data has an average of: 0.45683074900874476 + The original data has an standard deviation of: 0.1332971573460381 + The data generated directly from gaussian model has an average of: 0.20636891906272306 + The data generated directly from gaussian model has an standard deviation of: 0.24861526646143148 + The fitted data has an average of: 0.454608092485463 + The fitted data has an standard deviation of: 0.11042951535837742 + + +### 2. Procesando la información de la Tercera estrella seleccionada +Trabajamos con los arrays conteniendo los valores para las variables de posición (x3 e y3) y de la intensidad de pixeles (z3) provenientes de la imagen de tercera estrella seleccionada. + + +```python +#Working with a THIRD STAR +#based in coordinates of gray image crop the stars +star03 = imagen_gris[369:388,441:460] +plt.imshow(star03, cmap ='gray') +plt.grid(None) +plt.show() +``` + + + + + + + + +```python +# Generate a 3D Plot +# Import libraries +from mpl_toolkits import mplot3d +import numpy as np +import matplotlib.pyplot as plt + +# get coordinates (y,x) --- alternately see below for (x,y) +yx_coords3 = np.column_stack(np.where(star03 >= 0)) + +#Generating vectors of x and y coordinates from x,y array +y3=[item[0] for item in yx_coords3] +x3=[item[1] for item in yx_coords3] + +#Generating a vector from matrix of pixeles +zcoords3 = star03.flatten() + +# Creating dataset +z3 =zcoords3 + +# Creating figure +fig = plt.figure(figsize = (10, 7)) +ax = plt.axes(projection ="3d") + +# Creating plot +ax.scatter3D(x3, y3, z3, color = "green") +plt.title("Simple 3D scatter plot for Star 3") + +# show plot +plt.show() + +``` + + + + + + + + +```python +from scipy.optimize import leastsq + +xTa3 = ar(x3) +yTa3 = ar(y3) + +#Calculating parameters +xTa3_average=np.average(xTa3) +yTa3_average=np.average(yTa3) + +sigma_yTa3 = sqrt(sum(z3*(yTa3-9)**2)/361) + +parameters3=[9,9,sigma_yTa3,sigma_yTa3,1] + +def gaussian3D(params,x, y): + z = params[4] * np.exp( -0.5*((x-params[0])/(params[2]))**2 -0.5*((y-params[1])/(params[3]))**2) + return z +zOriginal3=gaussian3D(parameters3,xTa3,yTa3) + + +def error3D(params,x,y,z, cols, rows): + errors = (gaussian3D(params,x,y)-z) + return errors.reshape(cols*rows) + +obtained_error3=error3D(parameters3,xTa3,yTa3, z3, 19, 19) + +print('The average error between gaussian model and star 3 data is:') +print(np.average(obtained_error2)) + + +best3D3,sus3D3 = leastsq(error3D, parameters3, args = (xTa3,yTa3, z3, 19, 19)) + +print('Fitted parameters:') +print(best3D3) + + +#Using the fitted parameters for obtaining pixel values +zFitted3=gaussian3D(best3D3,xTa3,yTa3) + +#Showing the imagen of star with original values, values according to gaussian function and values according to fitted parameters +plt.subplot(1, 3, 1) +star03_01 = z3.reshape(19, 19) +plt.title('original') +plt.imshow(star03_01,cmap="gray") +plt.grid(None) +plt.subplot(1, 3, 2) +star03_02 = zOriginal3.reshape(19, 19) +plt.title('Gaussian') +plt.imshow(star03_02,cmap="gray") +plt.grid(None) +plt.subplot(1, 3, 3) +star03_03 = zFitted3.reshape(19, 19) +plt.title('Fitted') +plt.imshow(star03_03,cmap="gray") +plt.grid(None) +plt.show() + +#Calculando algunos parámetros estadÃsticos para los 3 casos +print("The original data has an average of: ",np.average(z3)) +print("The original data has an standard deviation of: ",np.std(z3)) + +print("The data generated directly from gaussian model has an average of: ",np.average(zOriginal3)) +print("The data generated directly from gaussian model has an standard deviation of: ",np.std(zOriginal3)) + +print("The fitted data has an average of: ",np.average(zFitted3)) +print("The fitted data has an standard deviation of: ",np.std(zFitted3)) + + + +``` + + The average error between gaussian model and star 3 data is: + -0.2504618299460217 + Fitted parameters: + [ 9.18996877 8.30415795 11.50308609 11.38848193 0.58975201] + + + + + + + + + The original data has an average of: 0.4736430003802075 + The original data has an standard deviation of: 0.09551903614251601 + The data generated directly from gaussian model has an average of: 0.22400209120618317 + The data generated directly from gaussian model has an standard deviation of: 0.2523793578898535 + The fitted data has an average of: 0.4730893645054868 + The fitted data has an standard deviation of: 0.06740127908186569 + + + +```python + +```