Skip to content
Snippets Groups Projects
Commit dce8b94a authored by Felipe Reyes Osorio's avatar Felipe Reyes Osorio
Browse files

More stars more files

parent a31aa159
No related branches found
No related tags found
No related merge requests found
This diff is collapsed.
Felipe Reyes Osorio (@reyesf)
# Ejercicio 1
Comenzamos importando tres librerías importantes: numpy, matplotlib, y PIL.Image. Esta última la usamos exclusivamente para leer la imagen como un arreglo de numpy de dimensión (789,1184,3). Las primeras dos denotan las coordenadas del pixel, mientras que la tercera, el canal o color RGB. Los elementos del arreglo codifican la intensidad en escala de 0 a 1.
Con `imshow` mostramos cada canal por separado
```python
import numpy as np
from PIL import Image
from matplotlib.pyplot import *
from mpl_toolkits.mplot3d import Axes3D
```
```python
sky = np.array(Image.open("./data/zapatocaImage.jpeg"))
sky.shape
```
(789, 1184, 3)
```python
fig, axes= subplots(nrows=1, ncols=3,figsize=(16,10))
channels = ["Reds","Greens", "Blues"]
for i in range(3):
axes[i].title.set_text(channels[i])
axes[i].imshow(sky[:,:,i], cmap=channels[i]+"_r")
```
![png](output_3_0.png)
Ahora, para convertir a escala de grises se podría realizar el promedio de los tres canales. Sin embargo, este método no preserva la intensidad lumínica percibida. Para lograr una conversión más fiel a la realidad debemos primero linealizar los datos de cada canal, pues al ser codificados se les aplica un filtro, tal que
$$
C_{\mathrm {lineal} }={\begin{cases}{\frac {C_{\mathrm {filtro} }}{12.92}},&{\text{si }}C_{\mathrm {filtro} }\leq 0.04045\\\left({\frac {C_{\mathrm {filtro} }+0.055}{1.055}}\right)^{2.4},&{\text{de lo contrario}}\end{cases}}
$$
donde $C$ representa cualquier canal RGB. Una vez linealizamos los datos, debemos realizar un promedio ponderado, pues la visión humana percibe con mayor intensidad el verde, y con menor intensidad el azul. De esta manera, la escala de grises $Y_{lineal}$ esta codificada por
$$
Y_{\mathrm {lineal} }=0.2126R_{\mathrm {lineal} }+0.7152G_{\mathrm {lineal} }+0.0722B_{\mathrm {lineal} }.
$$
Para guardar una imagen en escala de grises usando un formato de tres canales habría que aplicar la inversa del linealizador. Sin embargo, para nuestros propósitos es suficiente con graficar los datos a esta altura del proceso. Para más información, ver <a href="https://en.wikipedia.org/wiki/Grayscale#Converting_color_to_grayscale">este artículo</a>.
Teniendo todo esto en cuenta, definimos la función `linearizer` que quita el filtro de los canales individuales, y usamos `np.vectorize` para poder aplicársela a un arreglo de numpy. Luego definimos otra función, `grayScaler`, que linealiza cada canal, y realiza la suma ponderada apropiada. Adicionalmente, tiene un parámetro `naive` que permite realizar la transformación con el promedio simple (sin linealizar y sin ponderar). Graficamos ambas para mostrar la diferencia: se nota como el método más sofisticado da mejores resultados.
```python
def linearizer(x):
if x<=0.04045: return x/12.92
else: return ((x+0.055)/1.055)**2.4
linearizer = np.vectorize(linearizer)
```
```python
def grayScaler(mat, naive=False):
copyMat=np.zeros(mat[:,:,1].shape)
coefs=[0.2126,0.7152,0.0722]
if naive:
for chan in range(3):
copyMat += 0.333333*mat[:,:,chan]
return copyMat
else:
for chan in range(3):
copyMat += coefs[chan]*linearizer(mat[:,:,chan])
return copyMat
```
```python
fig1, ax1= subplots(nrows=1, ncols=2,figsize=(15,10))
ax1[0].title.set_text("Método sofisticado")
ax1[1].title.set_text("Método sencillo")
ax1[1].imshow(grayScaler(sky, naive=True), cmap="Greys_r")
ax1[0].imshow(grayScaler(sky), cmap="Greys_r")
```
<matplotlib.image.AxesImage at 0x7f68154c4be0>
![png](output_7_1.png)
Guardamos la información de la imagen en blanco y negro en la variable `graySky`, y seleccionamos una porción del cielo que contenga varias estrellas brillantes. El objetivo es estudiar ocho estrellas aquí guardadas pues todas estarán sujetas a condiciones aproximadamente iguales. Las estrellas seleccionadas se guardan en la lista `stars`, de manera que el punto más brillante quede cerca del centro de la imagen.
La porción de cielo escogida se muestra a continuación
```python
graySky = grayScaler(sky)
portion = graySky[280:420,400:650]
stars=[portion[85:105,40:60],portion[20:30,73:85],
portion[28:45,220:235],portion[114:123,218:226],
portion[55:65,135:145],portion[67:76,142:152],
portion[80:89,148:156],portion[42:50,161:170],
portion[9:17,190:200],portion[37:44,209:216],
portion[78:85,156:163],portion[20:28,232:240]]
fig2, ax2= subplots(nrows=1, ncols=2,figsize=(14,4))
ax2[0].title.set_text("Porción de cielo")
ax2[1].title.set_text("Estrella ejemplo stars[2]")
ax2[0].imshow(portion, "Greys_r")
ax2[1].imshow(stars[2], "Greys_r")
```
<matplotlib.image.AxesImage at 0x7f67f00c97f0>
![png](output_9_1.png)
Ahora definimos la función Gaussiana que se busca ajustar, y la función que calcula los errores. La primera, `gauss`, toma de argumentos coordenadas $x$ y $y$, además de una tupla de `params` organizada de esta manera:
* ($x_0$, $y_0$, $s$, $c$, $A$)
y retorna
$$
A\cdot gauss = e^{-(x^2+y^2)/2s^2}+c.
$$
La función `error` es sencilla, solo vale la pena mencionar que se debe reducir la dimensión del arreglo de errores con la función `flatten()`.
```python
def gauss(params,x, y):
x0,y0=params[1],params[0]
dis=(x-x0)**2+(y-y0)**2
return params[4]*np.exp(-0.5*dis**2/params[2]**2)+params[3]
def error(tpl, x,y, mat):
er = gauss(tpl, x,y)-mat
return er.flatten()
```
Ahora importamos la función `least_squares` de SciPy; esta es más recomendable que `leastsq` pues en funciona en casos más generales y complejos, además de ser una versión más reciente.
Hacemos un bucle que itera sobre cada estrella en la lista `stars`. Los detalles del loop se dan en los comentarios; al final se obtiene una lista `best` con los resultados de la optimización de cada estrella.
```python
from scipy.optimize import least_squares
```
```python
results = []
for item in stars: # Iteramos sobre cada estrella.
Y=np.arange(len(item[:,0]))
X=np.arange(len(item[0,:])) # Definimos un meshgrid con las
# mismas dimensiones de la estrella
X,Y=np.meshgrid(X,Y) # en cuestión.
inix,iniy = np.where(item == np.amax(item)) # Calculamos el centro de la imagen
# para usar las coordenadas como centro
p0 = (inix[0], iniy[0], 5.,0.1, np.amax(item)) # de la Gaussiana, y llenamos el vector
# de condiciones iniciales
best = least_squares(error, p0, args = (X,Y, item)) # Le pasamos al optimizador los errores
# comparando con los datos normalizados.
results.append(best)
```
Mostramos los resultados de la estrella `star[2]` como ejemplo: los contornos representan la función de mejor ajuste.
```python
Y=np.arange(len(stars[2][:,0]))
X=np.arange(len(stars[2][0,:]))
X,Y=np.meshgrid(X,Y)
fig3, ax3= subplots(figsize=(8,6))
im1 = ax3.imshow(stars[2], "Greys_r")
im2 = ax3.contour(X,Y, gauss(results[2].x, X,Y), cmap="summer_r")
cbar_im1 = colorbar(im1)
cbar_im2 = colorbar(im2)
```
![png](output_16_0.png)
Finalmente, convertimos los datos del ajuste para obtener el FWHM (**F**ull **W**idth at **H**alf **M**aximum), sabiendo que para la Gaussiana
$$
FWHM = 2\sqrt{2s\ln 2}.
$$
Con estos datos realizamos un histograma y calculamos la media, mediana, y desviación estándar.
```python
fwhm = np.array([item.x[2] for item in results])
fwhm = 2*np.sqrt(2*fwhm*np.log(2))
```
```python
fig4, ax4=subplots(figsize=(9,6))
ylabel("Frecuencia", fontsize=20)
xlabel("FWHM", fontsize=20)
n1, bins1, patches1=hist(fwhm, 10, alpha=0.6)
```
![png](output_19_0.png)
```python
avg = np.mean(fwhm)
med = np.median(fwhm)
dev = np.std(fwhm)
print("El promedio es "+ '{:.3f}'.format(avg)+"\nLa mediana es "+ "{:.3f}".format(med)
+ "\nLa desviación estándar es "+ "{:.3f}".format(dev))
```
El promedio es 3.310
La mediana es 3.051
La desviación estándar es 1.089
### Parte final
Realizamos todo de nuevo pero usando el canal azul, ya que es el que menos peso tiene en la conversión a escala de grises. En este caso obviamos la mayoría de figuras, y mostramos directamente los resultados estadísticos. Observamos en general Gaussianas más anchas, y una acentuación de los datos extremales.
```python
blueSky = sky[:,:,2]
portion = blueSky[280:420,400:650]
stars=[portion[85:105,40:60],portion[20:30,73:85],
portion[28:45,220:235],portion[114:123,218:226],
portion[55:65,135:145],portion[67:76,142:152],
portion[80:89,148:156],portion[42:50,161:170],
portion[9:17,190:200],portion[37:44,209:216],
portion[78:85,156:163],portion[20:28,232:240]]
results = []
for item in stars:
Y=np.arange(len(item[:,0]))
X=np.arange(len(item[0,:]))
X,Y=np.meshgrid(X,Y)
inix,iniy = np.where(item == np.amax(item))
p0 = (inix[0], iniy[0], 5.,0.1, np.amax(item))
best = least_squares(error, p0, args = (X,Y, item))
results.append(best)
fwhm = np.array([item.x[2] for item in results])
fwhm = 2*np.sqrt(2*fwhm*np.log(2))
fig5, ax5=subplots(figsize=(9,6))
ylabel("Frecuencia", fontsize=20)
xlabel("FWHM (Solo azul)", fontsize=20)
n1, bins1, patches1=hist(fwhm, 10, alpha=0.6)
avg = np.mean(fwhm)
med = np.median(fwhm)
dev = np.std(fwhm)
print("El promedio es "+ '{:.3f}'.format(avg)+"\nLa mediana es "+ "{:.3f}".format(med)
+ "\nLa desviación estándar es "+ "{:.3f}".format(dev))
```
El promedio es 3.871
La mediana es 3.626
La desviación estándar es 1.461
![png](output_22_1.png)
```python
```
```python
```
```python
```
output_16_0.png

32.3 KiB

output_19_0.png

7.27 KiB

output_22_1.png

7.56 KiB

output_3_0.png

185 KiB

output_7_1.png

157 KiB

output_9_1.png

60.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