Skip to content
Snippets Groups Projects
Commit b1d4b64a authored by sirias's avatar sirias
Browse files

version final tarea

parent 1b238a2f
No related branches found
No related tags found
No related merge requests found
Showing
with 3094 additions and 881 deletions
source diff could not be displayed: it is too large. Options to address this: view the blob.
0.11 0 → 100644
Requirement already up-to-date: seaborn in c:\users\ecf0124a\anaconda3\lib\site-packages (0.11.1)
Requirement already satisfied, skipping upgrade: scipy>=1.0 in c:\users\ecf0124a\anaconda3\lib\site-packages (from seaborn) (1.5.0)
Requirement already satisfied, skipping upgrade: matplotlib>=2.2 in c:\users\ecf0124a\appdata\roaming\python\python38\site-packages (from seaborn) (3.3.2)
Requirement already satisfied, skipping upgrade: numpy>=1.15 in c:\users\ecf0124a\anaconda3\lib\site-packages (from seaborn) (1.18.5)
Requirement already satisfied, skipping upgrade: pandas>=0.23 in c:\users\ecf0124a\anaconda3\lib\site-packages (from seaborn) (1.0.5)
Requirement already satisfied, skipping upgrade: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.3 in c:\users\ecf0124a\anaconda3\lib\site-packages (from matplotlib>=2.2->seaborn) (2.4.7)
Requirement already satisfied, skipping upgrade: kiwisolver>=1.0.1 in c:\users\ecf0124a\anaconda3\lib\site-packages (from matplotlib>=2.2->seaborn) (1.2.0)
Requirement already satisfied, skipping upgrade: cycler>=0.10 in c:\users\ecf0124a\anaconda3\lib\site-packages (from matplotlib>=2.2->seaborn) (0.10.0)
Requirement already satisfied, skipping upgrade: pillow>=6.2.0 in c:\users\ecf0124a\anaconda3\lib\site-packages (from matplotlib>=2.2->seaborn) (7.2.0)
Requirement already satisfied, skipping upgrade: certifi>=2020.06.20 in c:\users\ecf0124a\anaconda3\lib\site-packages (from matplotlib>=2.2->seaborn) (2020.6.20)
Requirement already satisfied, skipping upgrade: python-dateutil>=2.1 in c:\users\ecf0124a\anaconda3\lib\site-packages (from matplotlib>=2.2->seaborn) (2.8.1)
Requirement already satisfied, skipping upgrade: pytz>=2017.2 in c:\users\ecf0124a\anaconda3\lib\site-packages (from pandas>=0.23->seaborn) (2020.1)
Requirement already satisfied, skipping upgrade: six in c:\users\ecf0124a\anaconda3\lib\site-packages (from cycler>=0.10->matplotlib>=2.2->seaborn) (1.15.0)
source diff could not be displayed: it is too large. Options to address this: view the blob.
# Ejercicios para practicar numpy y optimización con scipy
# Resolución espacial en Observaciones astronómicas: Ajuste Gaussiano y detección automática de objetos celestes
## ***Siria Sadeddin***
......@@ -15,7 +15,7 @@ demasiado cerca sus manchas se superpondrán hasta el punto en que sea imposible
distinguirlas como fuentes individuales.
![picture](http://www.carlostapia.es/fisica/resolucion_criterios_practica_files/resolucion_rayleigh.png)
![picture](http://www.carlostapia.es/fisica/resolucion_criterios_practica_files/resolucion_rayleigh.png) [1]
Para modelar este efecto, típicamente consideramos la acción de la atmósfera
......@@ -25,7 +25,7 @@ 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
La resolución espacial normalmente se toma como la "Anchura a media altura" (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
......@@ -37,16 +37,47 @@ 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.
## Ejercicios:
## Anchura a media altura (FWHM)
### Para una Gaussiana en 1D
La Anchura a media altura, abreviada FWHM (del inglés Full Width at Half Maximum) es una medida de la extensión de una función, que viene dada por la diferencia entre los dos valores extremos de la variable independiente en los que la variable dependiente es igual a la mitad de su valor máximo. [2]
La FWHM se aplica a fenómenos tales como la duración de un pulso y la anchura espectral de fuentes utilizadas para comunicaciones ópticas y la resolución de espectrómetros.
Cuando la función considerada es una distribución normal de la forma:
$$ \frac{1}{\sigma\sqrt{2\pi}}\exp{-\frac{(x-x_0)^2}{2\sigma^2}}$$
donde $\sigma$ es la desviación típica y $x_{0}$ puede ser cualquier valor (la anchura de la función no cambia con una traslación), la relación entre FWHM y la desviación típica es:
$$ FWHM = 2\sqrt{2\ln(2)}\sigma \approx 2.35482\sigma$$
![image](https://upload.wikimedia.org/wikipedia/commons/thumb/c/cb/FWHM.svg/250px-FWHM.svg.png)
### Para una Gaussiana en 2D
En dos dimensiones, el exponente de la potencia de e dentro de la función de Gauss es cualquier valor negativo y definido en forma cuadrática. Como consecuencia, los niveles de la función siempre serán elipses. [3]
Un ejemplo de la función de dos dimensiones es:
$$f(x,y)=A\exp \left(-\left({\frac {(x-x_{0})^{2}}{2\sigma _{X}^{2}}}+{\frac {(y-y_{0})^{2}}{2\sigma _{Y}^{2}}}\right)\right)$$
Para este experimento tomaremos el maximo entre $\sigma_{x}$ y $\sigma_{y}$ como resultado final para el cáculo del FWHM.
## Exploración de los datos
### Lectura de la imagen
* 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
Se leerá la imagen almacenada en la carpeta **"data"** llamada **"zapatocaImage.jpeg"** como un array de numpy. Ese array estará compuesto de 3 matrices, una tras otra, correspondiente a los canales **R,G,B**. Lo que significa:
* <span style="color:red">R</span>: Rojo
* <span style="color:green">G</span>: Verde
* <span style="color:blue">B</span>: Azul
```python
```
# load libraries
import numpy as np
import math
......@@ -56,20 +87,24 @@ import matplotlib.patches as patches
import pandas as pd
from mpl_toolkits.mplot3d import Axes3D
from scipy import optimize
import os
import seaborn as sns
from scipy.stats import norm
import statistics as st
```
```python
```
# load image
image = image.imread('./data/zapatocaImage.jpeg')
img = np.asarray(image.imread('./data/zapatocaImage.jpeg'))
# summarize shape of the pixel array
print(image.dtype)
print(image.shape)
print(img.dtype)
print(img.shape)
# display the array of pixels as an image
pyplot.imshow(image)
figure, plot = pyplot.subplots(figsize=(15,10))
plot.imshow(img)
plot.set_axis_off()
pyplot.show()
```
......@@ -78,27 +113,30 @@ pyplot.show()
![png](output_8_1.png)
![png](output_10_1.png)
## Pretratamiento de la imagen
### Conversion de la imagen a blanco y negro
* 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)
Para hacer un estudio primario de la imagen, la convertiremos a una imagen de un solo canal (escala de grises) que contienne para cada nuevo pixel el promedio de valores de los pixeles en la imagen RGB. En principio, podemos hacer cualquier tipo de combinacion lineal de los tres canales RGB para obtener la imagen en escala de grises.
```python
```
#compute the mean of the image pixels in the channel axis
image_gray_m=np.mean(image,axis=2)
image_gray_m=np.mean(img,axis=2)
pyplot.figure(figsize=(15,10))
# summarize shape of the pixel array
print(image_gray_m.dtype)
print(image_gray_m.shape)
# display the array of pixels as an image
pyplot.imshow(image_gray_m,cmap='gray')
figure, plot = pyplot.subplots(figsize=(15,10))
plot.imshow(image_gray_m,cmap='gray')
plot.set_axis_off()
pyplot.show()
```
......@@ -107,16 +145,20 @@ pyplot.show()
![png](output_11_1.png)
<Figure size 1080x720 with 0 Axes>
![png](output_14_2.png)
Como el color es una percepcion psicologica, se ha encontrado que una mejor aproximacion al problema es hacer la combinación lineal:
$y_{gray}=0.2126R_{lineal}+0.7252G_{lineal}+0.0722B_{lineal}$
$$y_{gray}=0.2126R_{lineal}+0.7252G_{lineal}+0.0722B_{lineal}$$
```python
image_gray_l=image[:,:,0]*0.2126+image[:,:,1]*0.7252+image[:,:,2]*0.0722
```
image_gray_l=img[:,:,0]*0.2126+img[:,:,1]*0.7252+img[:,:,2]*0.0722
# summarize shape of the pixel array RGB BGR
print(image_gray_l.dtype)
......@@ -125,7 +167,9 @@ print(image_gray_l.shape)
# display the array of pixels as an image
fig, (ax1, ax2) = pyplot.subplots(1, 2,figsize=(20,40))
ax1.imshow(image_gray_m,cmap='gray')
ax1.set_axis_off()
ax2.imshow(image_gray_l,cmap='gray')
ax2.set_axis_off()
pyplot.show()
```
......@@ -134,70 +178,149 @@ pyplot.show()
![png](output_13_1.png)
![png](output_16_1.png)
Podriamos usar cualquera de los dos resultados, pero me parece que la combinación linear presenta una imagen mas clara (image a la derecha)
Podriamos usar cualquera de los dos resultados, pero me parece que la promedio presenta una imagen mas clara (image a la izquierda)
### Extracción de las estrellas en la imagen mediante el método de detección de outliers
* 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)
### Valores atípicos o outliers
En estadística, tales como muestras estratificadas, un valor atípico (en inglés outlier) es una observación que es numéricamente distante del resto de los datos. En una distribucion normal, los valores que se encuentran mas por encima o por debajo de 3 desviaciones estandar representan menos del 0.1% de la población, por esta razon podemos decir que estos valores son extraños para esa distribución.
![image](https://upload.wikimedia.org/wikipedia/commons/thumb/8/8c/Standard_deviation_diagram.svg/350px-Standard_deviation_diagram.svg.png)
Para hallar la posición de las estrellas en el cielo utilizaremos un método de deteccion de outliers en imagenes, el cielo es obscuro, por lo tanto las estrellas son outliers en la imagen. lo que se hará es calcular la media y desviación estandar de valores en los pixeles de la imagen, para luego encontrar los valores que se alejen de la media mas de 3 desviaciones estandar, estos valores son las estrellas. Además de esto, la selección de estrellas se hará para la sección de la imagen $y<400$ y $250<x<1000$, esto elimina los objetos que no pertenecen al cielo en obsevacion. Se nota la precencia de un satelite (presuntamente de Elon Musk :D) pasando por el cielo, que no corresponde a las estrellas en estudio (linea punteada brillante en la parte inferior de la imagen). Además, la estrella que se nota mas abajo parece ser un planeta o algún otro objeto celeste (opinion personal).
Como hemos mencionado, se tomará como estrella todos los valores que superen 3 desviaciones estandar por encima de la media (los valores altos en los pixeles son brillantes y los cercanos a cero son obscuros o negros)
$$estrella=mean_{pixels}+3*std_{pixels}$$
Para hallar la posición de las estrellas en el cielo utilizaremos un método de deteccion de outliers en imagenes, el cielo es obscuro, por lo tanto las estrellas son outliers en la imagen, lo que se hará es calcular la media y desviación estandar de valores en los pixeles de la imagen, para luego encontrar los valores que se alejen de la media mas de 4 desviaciones estandar, estos valores son las estrellas. Además de esto, la selección de estrellas se hará para la sección de la imagen $y<500$ ya que se nota la precencia de un satelite pasando por el cielo que no corresponde a las estrellas en estudio, además la estrella que se nota mas abajo parece ser un planeta (opinion personal).
```python
#### Seleccionaremos la zona de la foto que solo incluya el cielo nocturno
```
image_gray_sky = image_gray_m[:400,250:1000]
image_sky = img[:400,250:1000]
# display the array of pixels as an image
fig, (ax1, ax2) = pyplot.subplots(1, 2,figsize=(20,40))
ax1.imshow(image_gray_sky,cmap="gray")
ax1.set_axis_off()
ax1.set_title("Cielo en escala de grises",fontsize=16)
ax2.imshow(image_sky,cmap="gray")
ax2.set_axis_off()
ax2.set_title("Cielo con canales RGB",fontsize=16)
pyplot.show()
```
![png](output_21_0.png)
#### Demostraremos que los datos del cielo siguen la distribución normal, por ende es natural hacer la detección de outliers como hemos descrito antes
```
h,w=image_gray_sky.shape
fig, (ax1) = pyplot.subplots(1, 1,figsize=(15,5))
sns.set_color_codes()
sns.distplot(image_gray_sky.reshape(h*w),ax=ax1,fit=norm,color="y")
ax1.legend(labels=['Distribución Normal','Distribución real'])
ax1.set_title("Distribucion de valores en los pixeles")
ax1.set_ylabel("frecuencia",fontweight ='bold',fontsize=14)
ax1.set_xlabel("Valor del pixel",fontweight ='bold',fontsize=14)
pyplot.show()
```
C:\Users\ECF0124A\Anaconda3\lib\site-packages\seaborn\distributions.py:2557: FutureWarning: `distplot` is a deprecated function and will be removed in a future version. Please adapt your code to use either `displot` (a figure-level function with similar flexibility) or `histplot` (an axes-level function for histograms).
warnings.warn(msg, FutureWarning)
![png](output_23_1.png)
```
print({"min":np.min(image_gray_sky),"max":np.max(image_gray_sky),"mean":np.mean(image_gray_sky),"std":np.std(image_gray_sky)})
```
{'min': 46.0, 'max': 255.0, 'mean': 84.9697677777778, 'std': 13.678023265705608}
#### Crearemos una funcion que nos permitirá hayar las estrellas mediante deteccion de outliers
```
# searching outliers
def star_search(image):
'''
Esta función identifica las estrellas en el
cielo nocturno identificando valores atípicos
en los pixeles de la imagen
Input:
image: array en 2D
output:
list: posiciones de los valores atípicos encontrados
'''
image_sec=image[:500,:]
mean=np.mean(image_sec)
std=np.std(image_sec)
mean=np.mean(image)
std=np.std(image)
h,w=image_sec.shape
h,w=image.shape
stars=[]
for i in range(h):
for j in range(w):
if image_sec[i,j]>=mean+4.5*std:
if (image[i,j]>=mean+6*std) and (i+15<h) and (j+15<w) and (i-15>0) and (j-15>0): ## dont take stars near the borders
stars.append([i,j])
return stars
stars=star_search(image_gray_l)
stars=star_search(image_gray_sky)
stars[0:5]
```
[[19, 677], [20, 609], [21, 609], [34, 333], [42, 723]]
[[18, 427], [19, 357], [19, 426], [19, 427], [20, 227]]
#### Seleccionaremos 8 estrellas de las halladas por nuestro algoritmo con un corte de $ 3 \times 3$ pixeles para evaluar el ressultado.
#### Seleccionaremos 8 estrellas de las halladas por nuestro algoritmo con un corte de $ 3 \times 3$ pixeles para evaluar el resultado.
```python
```
#corte de la seccion donde se encuentra la estrella
fig, axs = pyplot.subplots(2, 4, figsize=(15,10))
for i in range(2):
for j in range(4):
image_unit=image_gray_l[stars[i+j][0]-6:stars[i+j][0]+6,stars[i+j][1]-6:stars[i+j][1]+6]
image_unit=image_gray_sky[stars[i+j][0]-6:stars[i+j][0]+6,stars[i+j][1]-6:stars[i+j][1]+6]
axs[i, j].imshow(image_unit,cmap='gray')
axs[i, j].set_axis_off()
pyplot.show()
```
![png](output_20_0.png)
![png](output_28_0.png)
### Vamos a marcar el centroide de cada estrella solo por curiosidad
#### Vamos a marcar el centroide de cada estrella solo por curiosidad
```python
```
#corte de la seccion donde se encuentra la estrella
def star_centroid(image,stars):
......@@ -206,94 +329,111 @@ def star_centroid(image,stars):
for i in range(2):
for j in range(4):
image_unit=image[stars[i+j][0]-8:stars[i+j][0]+8,stars[i+j][1]-8:stars[i+j][1]+8]
image_unit=image[stars[i+j][0]-6:stars[i+j][0]+6,stars[i+j][1]-6:stars[i+j][1]+6]
image_min_max=(image_unit-np.min(image_unit))/(np.max(image_unit)-np.min(image_unit))
center_x=math.ceil(np.sum(np.sum(image_min_max,axis=1)*range(16))/np.sum(np.sum(image_min_max,axis=1)))
center_y=math.ceil(np.sum(np.sum(image_min_max,axis=0)*range(16))/np.sum(np.sum(image_min_max,axis=0)))
center_x=math.ceil(np.sum(np.sum(image_unit,axis=1)*range(12))/np.sum(np.sum(image_unit,axis=1)))
center_y=math.ceil(np.sum(np.sum(image_unit,axis=0)*range(12))/np.sum(np.sum(image_unit,axis=0)))
axs[i, j].plot(center_x,center_y,c='r',marker="+", markersize=40)
axs[i, j].imshow(image_min_max,cmap='gray')
axs[i, j].imshow(image_unit,cmap='gray')
axs[i, j].set_axis_off()
pyplot.show()
star_centroid(image_gray_l,stars)
star_centroid(image_gray_sky,stars)
```
![png](output_22_0.png)
![png](output_30_0.png)
#### Ahora veamos el resultado marcando todas las estrellas halladas rodeandolas con una caja de color rojo en la imagen original
```python
```
#dentro de la caja roja esta la estrella mas brillante detectada
def star_boxing(image,stars,color):
# Create figure and axes
fig, ax = pyplot.subplots()
fig, ax = pyplot.subplots(figsize=(14,8))
# Display the image
ax.imshow(image)
for i in range(len(stars)):
# Create a Rectangle patch
rect = patches.Rectangle((stars[i][1]-3,stars[i][0]-3), 10, 10, linewidth=1, edgecolor=color, facecolor='none')
rect = patches.Rectangle((stars[i][1]-5,stars[i][0]-5), 10, 10, linewidth=1, edgecolor=color, facecolor='none')
# Add the patch
ax.add_patch(rect)
ax.set_axis_off()
star_boxing(image,stars,'r')
star_boxing(image_sky,stars,'b')
```
![png](output_24_0.png)
![png](output_32_0.png)
## Ajuste Gaussiano de la intensidad luminosa de las estrellas
Veremos que la luminosidad de la estrella se distribuye en forma gaussiana, por lo tanto es natural hacer un ajuste gaussiano a los datos de esta.
#### Vamos a graficar la estrella en un espacio tridimensional. Veremos que la luminosidad de la estrella se distribuye en forma gaussiana, por lo tanto es natural hacer un ajuste gaussiano a los datos de esta.
### Vamos a graficar la estrella en un espacio tridimensional.
```python
```
def surface_plot (matrix, **kwargs):
'''
Hace un grafico de superficie a partir de una matriz de valores para el eje z
'''
# acquire the cartesian coordinate matrices from the matrix
# x is cols, y is rows
(x, y) = np.meshgrid(np.arange(matrix.shape[0]), np.arange(matrix.shape[1]))
fig = pyplot.figure(figsize=(15,5))
fig = pyplot.figure(figsize=(10,10))
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(x, y, matrix, **kwargs)
return (fig, ax, surf)
#corte de la seccion donde se encuentra la estrella
i=0
image_unit=image_gray_l[stars[i][0]-8:stars[i][0]+8,stars[i][1]-8:stars[i][1]+8]
image_min_max=(image_unit-np.min(image_unit))/(np.max(image_unit)-np.min(image_unit))
center_x=math.ceil(np.sum(np.sum(image_min_max,axis=1)*range(16))/np.sum(np.sum(image_min_max,axis=1)))
center_y=math.ceil(np.sum(np.sum(image_min_max,axis=0)*range(16))/np.sum(np.sum(image_min_max,axis=0)))
image_center = image_min_max[center_y-3:center_y+3,center_x-3:center_x+3]
# m.shape must be (6,6)
def star3D(image,title):
'''
Toma la imagen recortada de una estrella y
hace un grafico de superficie para la intensidad de luz
'''
m = np.fromfunction(lambda x, y: image, (6, 6))
(fig, ax, surf) = surface_plot(m, cmap = pyplot.cm.coolwarm)
pyplot.title(title)
pyplot.show()
```
```
# estrella seleccionada
i=0
#corte de la seccion donde se encuentra la estrella
image_unit=image_gray_sky[stars[i][0]-8:stars[i][0]+8,stars[i][1]-8:stars[i][1]+8]
#centroide de la estrella
center_x=math.ceil(np.sum(np.sum(image_unit,axis=1)*range(16))/np.sum(np.sum(image_unit,axis=1)))
center_y=math.ceil(np.sum(np.sum(image_unit,axis=0)*range(16))/np.sum(np.sum(image_unit,axis=0)))
star3D(image_center,"3D representation of a star intensity")
#recorte alrededor del centroide
image_center = image_unit[center_y-3:center_y+3,center_x-3:center_x+3]
# grafico de la superficie recortada
star3D(image_center,"Representación en 3D de la intesidad luminosa de una estrella")
```
![png](output_26_0.png)
![png](output_36_0.png)
#### Procedemos entonces a hacer un ajuste gaussiano de los datos para una estrella
### Procedemos entonces a hacer un ajuste gaussiano de los datos para una estrella
```python
```
# calcula los parametros de la funcion gaussiana de acuerdo a los datos entregados
def moments(data):
"""Returns (height, x, y, width_x, width_y)
the gaussian parameters of a 2D distribution by calculating its
......@@ -317,8 +457,8 @@ def moments(data):
return height, x, y, width_x, width_y
#once we have the parameters of the Gaussian curve from the provided data we make a function that calculates
#the predicted values
# con los parametros calculados por la funcion "moments"
# se predicen los datos de una funcion gaussiana que tenga esos parametros
def gaussian(height, center_x, center_y, width_x, width_y):
"""Returns a gaussian function with the given parameters"""
......@@ -327,85 +467,88 @@ def gaussian(height, center_x, center_y, width_x, width_y):
return lambda x,y: height*np.exp(
-(((center_x-x)/width_x)**2+((center_y-y)/width_y)**2)/2)
#optimize de parameters given by moments function
# se optimizan los parametros por minimos cuadrados
def fitgaussian(data):
"""Returns (height, x, y, width_x, width_y)
the gaussian parameters of a 2D distribution found by a fit"""
los parámetros optimizados por minimos cuadrados"""
params = moments(data)
#function errorfunction will be applied over paramereters p
#the function finds the difference between the actual data points and the predicted by the
## gaussian function
#la funcion errorfunction se aplicará sobre los parametros p
#se calcula la diferencia de los valores predichos y los reales
errorfunction = lambda p: np.ravel(gaussian(*p)(*np.indices(data.shape)) - data)
#coputes mean square error of the parameters and find the minimun error
#se optimiza la funcion de error, hallando los parametros que
# proveen el minimo error cuadrático
p, success = optimize.leastsq(errorfunction, params)
return p
```
### Graficaremos el resultado de la optimizacion
```python
```
# Create the gaussian data
Xin, Yin = image_center.shape
data = image_center
pyplot.matshow(data, cmap="Reds")
params = fitgaussian(data)
fit = gaussian(*params)
pyplot.matshow(data, cmap="Reds")
pyplot.contour(fit(*np.indices(data.shape)), cmap=pyplot.cm.copper)
ax = pyplot.gca()
(height, x, y, width_x, width_y) = params
pyplot.text(0.95, 0.05, """
height : %.1f
x : %.1f
y : %.1f
width_x : %.1f
width_y : %.1f""" %(x, y, width_x, width_y),
fontsize=12, horizontalalignment='right',
verticalalignment='bottom', transform=ax.transAxes)
width_y : %.1f"""
%(height,x, y, width_x, width_y),
fontsize=12, horizontalalignment='right',
verticalalignment='bottom', transform=ax.transAxes)
pyplot.show()
```
![png](output_29_0.png)
![png](output_40_0.png)
#### Observemos graficamente el proceso de ajuste gaussiano paso a paso
### Observemos graficamente el proceso de ajuste gaussiano paso a paso
```python
```
star3D(fit(*np.indices(data.shape)),"Gaussian Fit")
star3D(image_center,"Cropped data")
star3D(image_min_max,"Original Data")
star3D(image_unit,"Original Data")
```
![png](output_31_0.png)
![png](output_42_0.png)
![png](output_31_1.png)
![png](output_42_1.png)
![png](output_31_2.png)
![png](output_42_2.png)
### Vamos a hacer el ajuste gaussiano para todas las estrellas detectadas, con los parametros calculados podemos generar estadísticas interesantes
## Ajuste gaussiano para todas las estrellas
* 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
Vamos a hacer el ajuste gaussiano para todas las estrellas detectadas, con los parametros calculados podemos generar estadísticas interesantes
```python
```
def star_parameters(image,stars):
#initialize parameters list
......@@ -414,18 +557,17 @@ def star_parameters(image,stars):
## iteration over stars
for star in stars:
image_unit = image[star[0]-5:star[0]+5,star[1]-5:star[1]+5]
#centroide de la estrella
center_x=math.ceil(np.sum(np.sum(image_unit,axis=1)*range(10))/np.sum(np.sum(image_unit,axis=1)))
center_y=math.ceil(np.sum(np.sum(image_unit,axis=0)*range(10))/np.sum(np.sum(image_unit,axis=0)))
#recorte alrededor del centroide
image_center = image_unit[center_y-3:center_y+3,center_x-3:center_x+3]
image_unit = image[star[0]-6:star[0]+6,star[1]-6:star[1]+6]
image_min_max = (image_unit-np.min(image_unit))/(np.max(image_unit)-np.min(image_unit))
## centroid of the star
center_x = math.ceil(np.sum(np.sum(image_min_max,axis=1)*range(12))/np.sum(np.sum(image_min_max,axis=1)))
center_y = math.ceil(np.sum(np.sum(image_min_max,axis=0)*range(12))/np.sum(np.sum(image_min_max,axis=0)))
# cropped star
image_center = image_min_max[center_y-3:center_y+3,center_x-3:center_x+3]
## gaussian fit parameters
param = fitgaussian(image_center)
parameters.append(param)
......@@ -433,20 +575,24 @@ def star_parameters(image,stars):
parameters_df = pd.DataFrame(parameters,columns=["height", "mean_x", "mean_y", "std_x", "std_y"])
parameters_df["FWHM_x"]=2.35*parameters_df["std_x"]
parameters_df["FWHM_y"]=2.35*parameters_df["std_y"]
parameters_df
parameters_df["FWHM_x"] = 2.35*parameters_df["std_x"]
parameters_df["FWHM_y"] = 2.35*parameters_df["std_y"]
return parameters_df
parameters_df=star_parameters(image_gray_l,stars)
parameters_df=star_parameters(image_gray_sky,stars)
```
C:\Users\ECF0124A\Anaconda3\lib\site-packages\scipy\optimize\minpack.py:475: RuntimeWarning: Number of calls to function has reached maxfev = 1200.
warnings.warn(errors[info][0], RuntimeWarning)
#### Obtenemos un data set de parametros para el estudio
```python
parameters_df
```
parameters_df.sort_values('FWHM_x',ascending=False)
```
......@@ -481,54 +627,54 @@ parameters_df
</thead>
<tbody>
<tr>
<th>0</th>
<td>1.159396</td>
<td>2.745116</td>
<td>2.662124</td>
<td>1.003062</td>
<td>0.805998</td>
<td>2.357196</td>
<td>1.894096</td>
</tr>
<tr>
<th>1</th>
<td>0.827333</td>
<td>3.112776</td>
<td>2.698797</td>
<td>1.268333</td>
<td>1.302106</td>
<td>2.980583</td>
<td>3.059948</td>
</tr>
<tr>
<th>2</th>
<td>0.822568</td>
<td>2.149674</td>
<td>2.706100</td>
<td>1.294210</td>
<td>1.297260</td>
<td>3.041393</td>
<td>3.048560</td>
</tr>
<tr>
<th>3</th>
<td>0.888378</td>
<td>3.303022</td>
<td>3.154786</td>
<td>1.664099</td>
<td>1.073939</td>
<td>3.910632</td>
<td>2.523758</td>
</tr>
<tr>
<th>4</th>
<td>1.103022</td>
<td>2.232578</td>
<td>3.679309</td>
<td>1.136575</td>
<td>0.950258</td>
<td>2.670952</td>
<td>2.233107</td>
<th>90</th>
<td>14409.841451</td>
<td>-222.262629</td>
<td>2.343065</td>
<td>74.058552</td>
<td>2.488918</td>
<td>174.037597</td>
<td>5.848958</td>
</tr>
<tr>
<th>516</th>
<td>267271.868327</td>
<td>115.166949</td>
<td>2.769692</td>
<td>29.110955</td>
<td>2.450235</td>
<td>68.410745</td>
<td>5.758053</td>
</tr>
<tr>
<th>477</th>
<td>107493.971955</td>
<td>96.620752</td>
<td>3.101467</td>
<td>26.096038</td>
<td>3.044724</td>
<td>61.325689</td>
<td>7.155100</td>
</tr>
<tr>
<th>288</th>
<td>169.559297</td>
<td>-16.954423</td>
<td>2.975433</td>
<td>25.863356</td>
<td>3.902506</td>
<td>60.778886</td>
<td>9.170890</td>
</tr>
<tr>
<th>441</th>
<td>21979.573625</td>
<td>66.797420</td>
<td>2.603415</td>
<td>20.616066</td>
<td>2.744205</td>
<td>48.447755</td>
<td>6.448882</td>
</tr>
<tr>
<th>...</th>
......@@ -541,129 +687,547 @@ parameters_df
<td>...</td>
</tr>
<tr>
<th>205</th>
<td>0.919854</td>
<td>3.328575</td>
<td>2.181515</td>
<td>1.610132</td>
<td>1.405021</td>
<td>3.783809</td>
<td>3.301799</td>
</tr>
<tr>
<th>206</th>
<td>0.958172</td>
<td>2.357734</td>
<td>3.100645</td>
<td>1.499172</td>
<td>1.318109</td>
<td>3.523054</td>
<td>3.097557</td>
</tr>
<tr>
<th>207</th>
<td>1.069374</td>
<td>3.486341</td>
<td>3.118657</td>
<td>1.169026</td>
<td>0.924830</td>
<td>2.747210</td>
<td>2.173351</td>
</tr>
<tr>
<th>208</th>
<td>1.098067</td>
<td>2.498427</td>
<td>4.089789</td>
<td>1.169563</td>
<td>0.865817</td>
<td>2.748472</td>
<td>2.034670</td>
</tr>
<tr>
<th>209</th>
<td>1.007167</td>
<td>4.316177</td>
<td>1.770797</td>
<td>1.709760</td>
<td>1.290679</td>
<td>4.017935</td>
<td>3.033095</td>
<th>21</th>
<td>202.186462</td>
<td>2.289328</td>
<td>2.644426</td>
<td>2.005263</td>
<td>1.836719</td>
<td>4.712369</td>
<td>4.316290</td>
</tr>
<tr>
<th>29</th>
<td>210.588381</td>
<td>1.978485</td>
<td>3.167382</td>
<td>1.943865</td>
<td>1.831859</td>
<td>4.568083</td>
<td>4.304869</td>
</tr>
<tr>
<th>30</th>
<td>217.389734</td>
<td>1.980785</td>
<td>2.331454</td>
<td>1.931921</td>
<td>1.708335</td>
<td>4.540015</td>
<td>4.014588</td>
</tr>
<tr>
<th>26</th>
<td>216.603135</td>
<td>2.821754</td>
<td>3.173432</td>
<td>1.804289</td>
<td>1.820440</td>
<td>4.240079</td>
<td>4.278035</td>
</tr>
<tr>
<th>27</th>
<td>222.917653</td>
<td>2.823663</td>
<td>2.336885</td>
<td>1.803631</td>
<td>1.702241</td>
<td>4.238533</td>
<td>4.000267</td>
</tr>
</tbody>
</table>
<p>210 rows × 7 columns</p>
<p>524 rows × 7 columns</p>
</div>
#### Promedios de FWHM
#### El ajuste diverge para todos los height>255
```python
```
parameters_df=parameters_df[(parameters_df["height"]<255)]
parameters_df.sort_values('FWHM_x',ascending=True)
```
<div>
<style scoped>
.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}
.dataframe tbody tr th {
vertical-align: top;
}
.dataframe thead th {
text-align: right;
}
</style>
<table border="1" class="dataframe">
<thead>
<tr style="text-align: right;">
<th></th>
<th>height</th>
<th>mean_x</th>
<th>mean_y</th>
<th>std_x</th>
<th>std_y</th>
<th>FWHM_x</th>
<th>FWHM_y</th>
</tr>
</thead>
<tbody>
<tr>
<th>27</th>
<td>222.917653</td>
<td>2.823663</td>
<td>2.336885</td>
<td>1.803631</td>
<td>1.702241</td>
<td>4.238533</td>
<td>4.000267</td>
</tr>
<tr>
<th>26</th>
<td>216.603135</td>
<td>2.821754</td>
<td>3.173432</td>
<td>1.804289</td>
<td>1.820440</td>
<td>4.240079</td>
<td>4.278035</td>
</tr>
<tr>
<th>30</th>
<td>217.389734</td>
<td>1.980785</td>
<td>2.331454</td>
<td>1.931921</td>
<td>1.708335</td>
<td>4.540015</td>
<td>4.014588</td>
</tr>
<tr>
<th>29</th>
<td>210.588381</td>
<td>1.978485</td>
<td>3.167382</td>
<td>1.943865</td>
<td>1.831859</td>
<td>4.568083</td>
<td>4.304869</td>
</tr>
<tr>
<th>21</th>
<td>202.186462</td>
<td>2.289328</td>
<td>2.644426</td>
<td>2.005263</td>
<td>1.836719</td>
<td>4.712369</td>
<td>4.316290</td>
</tr>
<tr>
<th>...</th>
<td>...</td>
<td>...</td>
<td>...</td>
<td>...</td>
<td>...</td>
<td>...</td>
<td>...</td>
</tr>
<tr>
<th>309</th>
<td>151.959876</td>
<td>6.048875</td>
<td>3.048254</td>
<td>6.335821</td>
<td>3.106507</td>
<td>14.889180</td>
<td>7.300292</td>
</tr>
<tr>
<th>351</th>
<td>195.582213</td>
<td>7.582757</td>
<td>2.703236</td>
<td>6.548991</td>
<td>3.505496</td>
<td>15.390128</td>
<td>8.237915</td>
</tr>
<tr>
<th>425</th>
<td>251.594399</td>
<td>7.994551</td>
<td>1.593031</td>
<td>6.646625</td>
<td>3.216397</td>
<td>15.619568</td>
<td>7.558533</td>
</tr>
<tr>
<th>47</th>
<td>115.774855</td>
<td>4.630637</td>
<td>2.210311</td>
<td>7.270685</td>
<td>5.813427</td>
<td>17.086110</td>
<td>13.661553</td>
</tr>
<tr>
<th>288</th>
<td>169.559297</td>
<td>-16.954423</td>
<td>2.975433</td>
<td>25.863356</td>
<td>3.902506</td>
<td>60.778886</td>
<td>9.170890</td>
</tr>
</tbody>
</table>
<p>462 rows × 7 columns</p>
</div>
```
fig, (ax1,ax2) = pyplot.subplots(1, 2,figsize=(15,5))
sns.set_color_codes()
sns.distplot(parameters_df["FWHM_x"],ax=ax1,color="y")
ax1.legend(labels=['Distribución FWHM x','Distribución real'])
ax1.set_title("Distribucion de valores FWHM")
ax1.set_ylabel("frecuencia",fontweight ='bold',fontsize=14)
ax1.set_xlabel("FWHM",fontweight ='bold',fontsize=14)
sns.distplot(parameters_df["FWHM_y"],ax=ax2,color="y")
ax2.legend(labels=['Distribución FWHM y','Distribución real'])
ax2.set_title("Distribucion de valores FWHM")
ax2.set_ylabel("frecuencia",fontweight ='bold',fontsize=14)
ax2.set_xlabel("FWHM",fontweight ='bold',fontsize=14)
pyplot.show()
```
C:\Users\ECF0124A\Anaconda3\lib\site-packages\seaborn\distributions.py:2557: FutureWarning: `distplot` is a deprecated function and will be removed in a future version. Please adapt your code to use either `displot` (a figure-level function with similar flexibility) or `histplot` (an axes-level function for histograms).
warnings.warn(msg, FutureWarning)
C:\Users\ECF0124A\Anaconda3\lib\site-packages\seaborn\distributions.py:2557: FutureWarning: `distplot` is a deprecated function and will be removed in a future version. Please adapt your code to use either `displot` (a figure-level function with similar flexibility) or `histplot` (an axes-level function for histograms).
warnings.warn(msg, FutureWarning)
![png](output_49_1.png)
#### Minimo y máximo de FWHM
```
print({"MIN_FWHM_x":np.min(parameters_df["FWHM_x"]),"MIN_FWHM_y":np.min(parameters_df["FWHM_y"]),"MAX_FWHM_x":np.max(parameters_df["FWHM_x"]),"MAX_FWHM_y":np.max(parameters_df["FWHM_y"])})
```
{'MIN_FWHM_x': 4.238532797637339, 'MIN_FWHM_y': 4.000266867579268, 'MAX_FWHM_x': 60.778886236256064, 'MAX_FWHM_y': 31.31075303962154}
#### Moda de FWHM
```
print({"MODE_FWHM_x":st.mode(parameters_df["FWHM_x"]),"MODE_FWHM_y":st.mode(parameters_df["FWHM_y"])})
```
{'MODE_FWHM_x': 7.738015242592363, 'MODE_FWHM_y': 5.074958049968591}
#### Media de FWHM
```
print({"MEAN_FWHM_x":np.mean(parameters_df["FWHM_x"]),"MEAN_FWHM_y":np.mean(parameters_df["FWHM_y"])})
```
{'MEAN_FWHM_x': 3.4649369689063634, 'MEAN_FWHM_y': 3.112886739637155}
{'MEAN_FWHM_x': 7.940323515966517, 'MEAN_FWHM_y': 7.337040419590766}
#### Mediana de FHWM
```
print({"MEDIAN_FWHM_x":np.median(parameters_df["FWHM_x"]),"MEDIAN_FWHM_y":np.median(parameters_df["FWHM_y"])})
```
{'MEDIAN_FWHM_x': 7.396254460793973, 'MEDIAN_FWHM_y': 6.900833212383647}
#### Incertidumbre
```python
```
print({"INCERTIDUMBRE_FWHM_x":np.std(parameters_df["FWHM_x"])/math.sqrt(len(parameters_df["FWHM_x"])),"INCERTIDUMBRE_FWHM_y":np.std(parameters_df["FWHM_y"])/math.sqrt(len(parameters_df["FWHM_y"]))})
```
{'INCERTIDUMBRE_FWHM_x': 0.07533479193502118, 'INCERTIDUMBRE_FWHM_y': 0.10722958941431958}
{'INCERTIDUMBRE_FWHM_x': 0.15200550254891215, 'INCERTIDUMBRE_FWHM_y': 0.12480846830393083}
### Dado que ya tenemos funciones que permiten hacer los calculos de el ajuste gaussiano y detección de estrellas rápidamente, procederemos a realizar el mismo ejercicio, pero esta vez para cada canal RGB de la imagen
Vemos que hay una gran diferencia entre la media y la mediana de la distribución, esto se debe a que nuestros resultados para el FWHM tiene unos valores extremos inesperados y muy altos. Por ejemplo, uno de los FWHM es 60, lo cual significaría que no podemos distinguir dos estrellas en una imagen de 60x60, lo cual sabemos que es falso. Una manera de evitar este problema es tomar la mediana como referente en lugar de la media, ya que esta ordena los valores y toma el valor que está en la posición central, ignorando los valores extremos.
```
parameters_df=parameters_df[(parameters_df["FWHM_x"]<np.median(parameters_df["FWHM_x"])+3*np.std(parameters_df["FWHM_x"])/math.sqrt(len(parameters_df["FWHM_x"]))) & (parameters_df["FWHM_y"]<np.median(parameters_df["FWHM_y"])+3*np.std(parameters_df["FWHM_y"])/math.sqrt(len(parameters_df["FWHM_y"])))]
parameters_df.sort_values('FWHM_x',ascending=True)
```
Tenemos 2 resultados, uno para el eje x y otro para y, se tomará el maximo de los dos resultados y diremos que nuestra resolucion para diferenciar dos estrellas en una imagen seria de 4 pixeles. Si los centros de dos estrellas no se encuentran separados al menos por 4 pixeles no será posible distinguirlas.
* Repita el mismo ejercicio sobre cada una de las bandas R,G,B separadamente
y comente si observa diferencias en los resultados
```python
<div>
<style scoped>
.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}
.dataframe tbody tr th {
vertical-align: top;
}
.dataframe thead th {
text-align: right;
}
</style>
<table border="1" class="dataframe">
<thead>
<tr style="text-align: right;">
<th></th>
<th>height</th>
<th>mean_x</th>
<th>mean_y</th>
<th>std_x</th>
<th>std_y</th>
<th>FWHM_x</th>
<th>FWHM_y</th>
</tr>
</thead>
<tbody>
<tr>
<th>27</th>
<td>222.917653</td>
<td>2.823663</td>
<td>2.336885</td>
<td>1.803631</td>
<td>1.702241</td>
<td>4.238533</td>
<td>4.000267</td>
</tr>
<tr>
<th>26</th>
<td>216.603135</td>
<td>2.821754</td>
<td>3.173432</td>
<td>1.804289</td>
<td>1.820440</td>
<td>4.240079</td>
<td>4.278035</td>
</tr>
<tr>
<th>30</th>
<td>217.389734</td>
<td>1.980785</td>
<td>2.331454</td>
<td>1.931921</td>
<td>1.708335</td>
<td>4.540015</td>
<td>4.014588</td>
</tr>
<tr>
<th>29</th>
<td>210.588381</td>
<td>1.978485</td>
<td>3.167382</td>
<td>1.943865</td>
<td>1.831859</td>
<td>4.568083</td>
<td>4.304869</td>
</tr>
<tr>
<th>21</th>
<td>202.186462</td>
<td>2.289328</td>
<td>2.644426</td>
<td>2.005263</td>
<td>1.836719</td>
<td>4.712369</td>
<td>4.316290</td>
</tr>
<tr>
<th>...</th>
<td>...</td>
<td>...</td>
<td>...</td>
<td>...</td>
<td>...</td>
<td>...</td>
<td>...</td>
</tr>
<tr>
<th>400</th>
<td>148.694932</td>
<td>3.251515</td>
<td>3.323809</td>
<td>3.311014</td>
<td>3.089162</td>
<td>7.780883</td>
<td>7.259531</td>
</tr>
<tr>
<th>32</th>
<td>120.009106</td>
<td>2.582832</td>
<td>2.669800</td>
<td>3.323717</td>
<td>2.781075</td>
<td>7.810735</td>
<td>6.535527</td>
</tr>
<tr>
<th>150</th>
<td>140.401733</td>
<td>2.814937</td>
<td>2.461857</td>
<td>3.330074</td>
<td>3.020962</td>
<td>7.825674</td>
<td>7.099261</td>
</tr>
<tr>
<th>281</th>
<td>192.304218</td>
<td>3.922074</td>
<td>2.380541</td>
<td>3.330268</td>
<td>2.254765</td>
<td>7.826129</td>
<td>5.298697</td>
</tr>
<tr>
<th>228</th>
<td>123.573138</td>
<td>2.781581</td>
<td>2.546152</td>
<td>3.334843</td>
<td>2.962802</td>
<td>7.836880</td>
<td>6.962584</td>
</tr>
</tbody>
</table>
<p>206 rows × 7 columns</p>
</div>
```
fig, (ax1,ax2) = pyplot.subplots(1, 2,figsize=(15,5))
sns.set_color_codes()
sns.distplot(parameters_df["FWHM_x"],ax=ax1,color="y")
ax1.legend(labels=['Distribución FWHM x','Distribución real'])
ax1.set_title("Distribucion de valores FWHM")
ax1.set_ylabel("frecuencia",fontweight ='bold',fontsize=14)
ax1.set_xlabel("FWHM",fontweight ='bold',fontsize=14)
sns.distplot(parameters_df["FWHM_y"],ax=ax2,color="y")
ax2.legend(labels=['Distribución FWHM y','Distribución real'])
ax2.set_title("Distribucion de valores FWHM")
ax2.set_ylabel("frecuencia",fontweight ='bold',fontsize=14)
ax2.set_xlabel("FWHM",fontweight ='bold',fontsize=14)
pyplot.show()
```
C:\Users\ECF0124A\Anaconda3\lib\site-packages\seaborn\distributions.py:2557: FutureWarning: `distplot` is a deprecated function and will be removed in a future version. Please adapt your code to use either `displot` (a figure-level function with similar flexibility) or `histplot` (an axes-level function for histograms).
warnings.warn(msg, FutureWarning)
C:\Users\ECF0124A\Anaconda3\lib\site-packages\seaborn\distributions.py:2557: FutureWarning: `distplot` is a deprecated function and will be removed in a future version. Please adapt your code to use either `displot` (a figure-level function with similar flexibility) or `histplot` (an axes-level function for histograms).
warnings.warn(msg, FutureWarning)
![png](output_62_1.png)
```
print({"MEDIAN_FWHM_x":np.median(parameters_df["FWHM_x"]),"MEDIAN_FWHM_y":np.median(parameters_df["FWHM_y"])})
```
{'MEDIAN_FWHM_x': 6.361155779828058, 'MEDIAN_FWHM_y': 5.865753036777117}
```
print({"INCERTIDUMBRE_FWHM_x":np.std(parameters_df["FWHM_x"])/math.sqrt(len(parameters_df["FWHM_x"])),"INCERTIDUMBRE_FWHM_y":np.std(parameters_df["FWHM_y"])/math.sqrt(len(parameters_df["FWHM_y"]))})
```
{'INCERTIDUMBRE_FWHM_x': 0.05897470200906251, 'INCERTIDUMBRE_FWHM_y': 0.0583083389179374}
Tomaremos el promedi de los resultados en MEAN_FWHM y diremos que:
$$FWHM\approx 6 pixeles$$
### Ajuste gaussiano para cada canal RGB
Dado que ya tenemos funciones que permiten hacer los calculos de el ajuste gaussiano y detección de estrellas rápidamente, procederemos a realizar el mismo ejercicio, pero esta vez para cada canal RGB de la imagen
```
figure, plots = pyplot.subplots(ncols=3, nrows=1,figsize=(15,10))
for i, subplot in zip(range(3), plots):
temp = np.zeros(image.shape, dtype='uint8')
temp[:,:,i] = image[:,:,i]
temp = np.zeros(image_sky.shape, dtype='uint8')
temp[:,:,i] = image_sky[:,:,i]
subplot.imshow(temp)
subplot.set_axis_off()
pyplot.show()
```
![png](output_44_0.png)
![png](output_67_0.png)
```python
```
color=['b','r','g']
for i in range(3):
temp = np.zeros(image.shape, dtype='uint8')
temp[:,:,i] = image[:,:,i]
temp = np.zeros(image_sky.shape, dtype='uint8')
temp[:,:,i] = image_sky[:,:,i]
stars = star_search(temp[:,:,i])
pyplot.axis=None
star_boxing(temp,stars,color[i])
```
![png](output_45_0.png)
![png](output_68_0.png)
![png](output_45_1.png)
![png](output_68_1.png)
![png](output_45_2.png)
![png](output_68_2.png)
```python
```
color=['RED','GREEN','BLUE']
for i in range(3):
......@@ -674,31 +1238,45 @@ for i in range(3):
RED = '\033[91m'
ENDC = '\033[0m'
temp = np.zeros(image.shape, dtype='uint8')
temp[:,:,i] = image[:,:,i]
temp = np.zeros(image_sky.shape, dtype='uint8')
temp[:,:,i] = image_sky[:,:,i]
#busca las estrellas
stars = star_search(temp[:,:,i])
# calcula parametros
parameters_df = star_parameters(temp[:,:,i],stars)
parameters_df=parameters_df[(parameters_df["height"]<255)]
parameters_df=parameters_df[(parameters_df["FWHM_x"]<np.median(parameters_df["FWHM_x"])+3*np.std(parameters_df["FWHM_x"])/math.sqrt(len(parameters_df["FWHM_x"]))) & (parameters_df["FWHM_y"]<np.median(parameters_df["FWHM_y"])+3*np.std(parameters_df["FWHM_y"])/math.sqrt(len(parameters_df["FWHM_y"])))]
#muestra resultados
print(bcolors.RED +"FWHM FOR CHANNEL "+color[i] + bcolors.ENDC)
print({"MEAN_FWHM_x":np.mean(parameters_df["FWHM_x"]),"MEAN_FWHM_y":np.mean(parameters_df["FWHM_y"])})
print({"INCERTIDUMBRE_FWHM_x":np.std(parameters_df["FWHM_x"])/math.sqrt(len(parameters_df["FWHM_x"])),"INCERTIDUMBRE_FWHM_y":np.std(parameters_df["FWHM_y"])/math.sqrt(len(parameters_df["FWHM_y"]))})
```
/usr/local/lib/python3.6/dist-packages/scipy/optimize/minpack.py:454: RuntimeWarning: Number of calls to function has reached maxfev = 1200.
C:\Users\ECF0124A\Anaconda3\lib\site-packages\scipy\optimize\minpack.py:475: RuntimeWarning: Number of calls to function has reached maxfev = 1200.
warnings.warn(errors[info][0], RuntimeWarning)
FWHM FOR CHANNEL RED
{'MEAN_FWHM_x': 3.4680630487569712, 'MEAN_FWHM_y': 7.0390829738856135}
{'INCERTIDUMBRE_FWHM_x': 0.09250146545235435, 'INCERTIDUMBRE_FWHM_y': 3.8420845224178906}
{'MEAN_FWHM_x': 6.85040177317909, 'MEAN_FWHM_y': 6.286265858550066}
{'INCERTIDUMBRE_FWHM_x': 0.07397839012804115, 'INCERTIDUMBRE_FWHM_y': 0.06519770408323534}
FWHM FOR CHANNEL GREEN
{'MEAN_FWHM_x': 3.4312513275072645, 'MEAN_FWHM_y': 3.081207982746502}
{'INCERTIDUMBRE_FWHM_x': 0.07271127428241228, 'INCERTIDUMBRE_FWHM_y': 0.09999192643863575}
{'MEAN_FWHM_x': 6.245164160492783, 'MEAN_FWHM_y': 5.852337602343364}
{'INCERTIDUMBRE_FWHM_x': 0.05644305840123783, 'INCERTIDUMBRE_FWHM_y': 0.05873235620560445}
FWHM FOR CHANNEL BLUE
{'MEAN_FWHM_x': 3.9344262901840414, 'MEAN_FWHM_y': 5.221441969716197}
{'INCERTIDUMBRE_FWHM_x': 0.17003334171988463, 'INCERTIDUMBRE_FWHM_y': 1.1152480701320464}
{'MEAN_FWHM_x': 5.656095446703643, 'MEAN_FWHM_y': 5.248410485649535}
{'INCERTIDUMBRE_FWHM_x': 0.0511572917957581, 'INCERTIDUMBRE_FWHM_y': 0.04996212650827945}
Encontramos que el mejor canal es el azul con $$FWHM \approx 6 $$
### Resultados
Se encuentra que el peor canal para hacer observaciones de estrellas es el rojo por su mayor FWHM y mayor incertidumbre. El mejor canal sería el verde, con menor FWHM y menor incertidumbre en la medida. En el caso de la imagen en escala de grises la distancia minima de los centroides para distinguir dos estrellas es de 4 pixeles.
Se encuentra que el peor canal para hacer observaciones de estrellas es el rojo por su mayor FWHM. El mejor canal sería el azul, con menor FWHM. En el caso de la imagen en escala de grises la distancia minima de los centroides para distinguir dos estrellas es de 6 pixeles.
[1] http://www.carlostapia.es/fisica/resolucion_criterios_practica.html
[2] https://es.wikipedia.org/wiki/Anchura_a_media_altura#:~:text=La%20Anchura%20a%20media%20altura,mitad%20de%20su%20valor%20m%C3%A1ximo.
[3] https://es.wikipedia.org/wiki/Funci%C3%B3n_gaussiana
Ejercicios_clase_5/output_10_1.png

556 KiB

Ejercicios_clase_5/output_14_2.png

357 KiB

Ejercicios_clase_5/output_16_1.png

293 KiB

Ejercicios_clase_5/output_21_0.png

298 KiB

Ejercicios_clase_5/output_23_1.png

28.5 KiB

Ejercicios_clase_5/output_28_0.png

6.99 KiB

Ejercicios_clase_5/output_30_0.png

7.99 KiB

Ejercicios_clase_5/output_32_0.png

321 KiB

Ejercicios_clase_5/output_36_0.png

73.7 KiB

Ejercicios_clase_5/output_40_0.png

26.2 KiB

Ejercicios_clase_5/output_42_0.png

71.2 KiB

Ejercicios_clase_5/output_42_1.png

71.7 KiB

Ejercicios_clase_5/output_42_2.png

96.4 KiB

Ejercicios_clase_5/output_49_1.png

26.2 KiB

Ejercicios_clase_5/output_62_1.png

31.5 KiB

Ejercicios_clase_5/output_67_0.png

61.4 KiB

0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment