### 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
importnumpyasnp
importmatplotlib.pyplotasplt
frommatplotlibimportimageasimg
%matplotlibinline
# 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
defgray_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
returngray
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]foriteminyx_coords]
x=[item[1]foriteminyx_coords]
#Generating a vector from matrix of pixeles
zcoords=star01.flatten()
```
```python
#Creating plot with the data
# Import libraries
frommpl_toolkitsimportmplot3d
importnumpyasnp
importmatplotlib.pyplotasplt
# 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
importnumpyasnp
importmatplotlib.pyplotasplt
plt.style.use('seaborn-whitegrid')
%matplotlibinline
Zlinea=[]
Ylinea=[]
j=0
foriinrange(z.size-1):
ifx[i]==12:
Zlinea=np.insert(Zlinea,j,z[i])
j=j+1
j=0
foriinrange(z.size-1):
ifx[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
importnumpyasnp
importmatplotlib.pyplotasplt
fromnumpyimportasarrayasar,exp,sqrt
fromscipy.optimizeimportcurve_fit
YLa=ar(Ylinea)
ZLa=ar(Zlinea)
n=len(ZLa)## <---
mean=sum(ZLa*YLa)/n
sigma=sqrt(sum(ZLa*(YLa-mean)**2)/n)
defgaus(x,a,mu,sigma):
returna*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
importnumpyasnp
fromscipy.optimizeimportcurve_fit
importmatplotlib.pyplotasplt
frommpl_toolkits.mplot3dimportAxes3D
fromnumpyimportasarrayasar,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
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
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
frommpl_toolkitsimportmplot3d
importnumpyasnp
importmatplotlib.pyplotasplt
# 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
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
frommpl_toolkitsimportmplot3d
importnumpyasnp
importmatplotlib.pyplotasplt
# 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
