Commit 39b6cc58 authored by Gerardo Semprum's avatar Gerardo Semprum
Browse files

Markdown de la entrega

parent 96d42a4f
#### Optimización y procesamiento de imágenes (numpy-scipy)
Gerardo Semprúm.
Universidad Central de Venezuela.
Módulo de datos.
El proyecto descrito es un primer paso a lo que seria la caracterización y procesamiento de datos e imagenes. El problema en cuestrión consta de una imagen que refleja muchisimas estrellas cada una con una luminosidad caracteristica. Este proyecto espera poder optimizar dichas imagenes y crear diferentes gaussianas de donde sea mas sencillo el procesamiento de datos.
Importamos las librerias que necesitaremos para trabajar en nuestro código. Dichas librerias fueron: **numpy**, **scipy** y **matplotlib**. Además de traer una herramienta de optimización contenida en scipy, llamada leastsq.
Por otro lado, aplicamos la libreria **sys** con la finalidad de darle a python un camino que seguir e importamos "cmap" de matplotlib que será util para procesar las imágenes.
```python
import numpy as np
import scipy
import matplotlib.pyplot as plt
from scipy.optimize import leastsq
```
```python
import sys
sys.path.append("./data/zapatocaImage.jpeg")
from matplotlib import cm
```
En la primera linea del codigo, traemos la imagen desde una carpeta destino, la cual automaticamente es procesada como un array. Ese array estará compuesto de 3 matrices, una tras otra, correspondiente a los canales R,G,B.
Podemos comprobar que es un array debido a que le pedimos a python que nos imprima la forma (.shape()) de nuestra imagen y obtenemos un array de 3 matrices cada una con 789 filas y 1184 columnas
```python
stars = plt.imread("./data/zapatocaImage.jpeg")
print(stars.shape)
plt.imshow(stars)
```
(789, 1184, 3)
<matplotlib.image.AxesImage at 0x7fc837567a58>
![png](output_6_2.png)
Separamos dicho array, asignado cada matriz a una variable "color" asignado asi los valores de rojo, verde y azul respectivamente. Por ejemplo, si pedimos la forma (.shape()) de cualquiera de nuestro nuevos arreglos, obtenemos una matriz 2D.
```python
#Asignando un cada color a un arreglo respectivo
red = stars[:,:,0]
green = stars[:,:,1]
blue = stars[:,:,2]
```
Para juntar el trio de matrices, podria establecerse una función "suma", en donde se multiplique cada color por una constante y despues se sumen entre sí.
**Uso de la constante**: Al transformar una imagen colorida a una escala de grises debe multiplicarse cada matriz correspondiente a los colores (R,G,B) por una constante propia de cada color para evitar la perdida de luminosidad de la imagen. Puede leerse mas al respecto en la referencia anexada:
[Convirtiendo color a escala de grises](https://en.wikipedia.org/wiki/Grayscale#Converting_color_to_grayscale)
Para sumar y multiplicar cada matriz "color" con su constante asignada, usamos la función de numpy.dot() la cual devuelve el producto punto entre nuestra matriz principal "stars" y una lista de con las constantes de cada color(R,G,B).
También es necesario aplicarle un mapeo de colores en escala de grises (cmap = "gray") para que la imagen adopte su escala de grises (este ultimo paso se repite en muchas imagenes a futuro por lo que no será explicado nuevamente).
```python
gray = np.dot(stars[...,:3], [0.299, 0.587, 0.114])
plt.imshow(gray, cmap="gray")
print(gray.shape)
plt.show()
```
(789, 1184)
![png](output_10_1.png)
Siguiendo el procedimiento, recortamos una sección de la imagen para aislar una estrella de nuestra preferencia para realizar una gaussiana de su luminosidad y la optimización de dicha gaussiana.
Para cortar la imagen solo es necesario crear una función (en este caso: **estrella1**) que sea un subconjunto de **gray** (nuestra imagen en escala de grises).
Para recortar una parte de la imagen, definimos la función de la siguiente manera:
$$estrellaX = gray[val_{y0}:val_{y1} , val_{x0}:Val_{x1}]$$
Siendo el eje "Y" el eje vertical y el eje "X" el eje horizontal. Finalmente solo pedimos imprimir la estrella deseada.
Para manejar con mas facilidad el array que contiene a nuestra estrella, usamos la variable "st1" (diminutivo de "star1"). Las estrellas siguientes usarán la misma etiqueta con diferente numeración
```python
#Recortando nuestra primera estrella
estrella1 = gray[540:570, 630:700]
plt.imshow(estrella1, cmap="gray")
st1 = np.asarray(estrella1)
print(st1.shape)
```
(30, 70)
![png](output_12_1.png)
Definimos una función gaussiana para modelar nuestro sistema de la forma:
$$Gauss = Ae^{-\frac{(x-x_0)^2 + (y-y_0)^2}{2C^2}}+ B $$
**Consejo para aplicar la función:** la función construida depende de varios parámetros además de Xe Y. El parametro 0 (A) es la amplitud deseada, el paramtro 1 (B) es una constante aditiva, el parametro 2 (C) es la desviación estandar, el parametro 3 (X0) es la posición inicial en X y el parametro 4 (Y0) es la posición inicial en Y.
```python
def gauss2D(params,x, y):
exponente = -((x-params[3])**2 + (y-params[4])**2) / (2*params[2]**2)
z = params[0]* np.exp(exponente) + params[1]
return z
#params[0]: constante de amplitud
#params[1]: constante aditiva
#params[2]: desviacion estandar
#params[3]: X0
#params[4]: Y0
```
Construmos un primer arreglo para modelar la primera estrella, se hace un arreglo de 70 números en X y 30 números en Y con la intención de recrear la imágen anterior.
Se realiza un .meshgrid() debido a que hace un arreglo de X-Y para evaluaciones vectorizadas de una escala X-Y sobre X-Y grids dando un arreglo de coordenada unidimensional.
```python
#Definiendo los valores modelo de nuestra primera estrella (Gaussiana 2D)
x = np.arange(0,70,1)
y = np.arange(0,30,1)
xx , yy = np.meshgrid(x,y)
```
Por aproximación damos unos parametros para tratar de simular nuestra estrella, asignamos la primera estrella a la variable zz haciendo uso de la función antes descrita e imprimimos
```python
params = [150,1,5,33,14]
zz = gauss2D(params,xx,yy)
plt.imshow(zz,cmap="gray")
```
<matplotlib.image.AxesImage at 0x7fc82e01d2b0>
![png](output_18_1.png)
Definiendo una función error, la cual es una función comparativa que trabajará enfrentando los valores reales y los valores modelados. Dicha función depende de tpl(que es una tupla con los parametros que se usaron para calcular el modelo), x, y e z.
Definimos dentro de la función a zmodel, que es nuestra función gauss2D antes descrita que depende de la tupla, x e y.
Finalmete los errores ("errors2D") serán la resta entre los valores modelo y los valores reales:
$$errors2D = Z_{model} - Z$$
Al realizar dicha resta, obtenemos una nueva matriz. De donde asignaremos a F(filas) y C(columnas) a los elementos de esa matriz respectivamente.
Luego de eso aplicaremos un reshape a la matriz del tipo (F*C) con el objetivo de que sea una función general para todas las estrellas a trabajar.
```python
#Definimos la función error:
def ErrorGauss2D(tpl,x,y,z):
zmodel = gauss2D(tpl,x,y)
errors2D = zmodel - z
F, C = errors2D.shape
return errors2D.reshape(F*C)
```
Damos un "chute inicial" (p1) y finalmente "best1" solo será una tupla con los mejores valores calculados entre el modelo y los valores reales. De donde, leastsq() es una función que minimiza la suma de cuadrados de un set de ecuaciones.
Notando que el "chute inicial" (p1) está definimo como los params usados en el modelo de la primera estrella. Estos chutes iran cambiando según la estrella.
Para obtener "best1", llamamos a nuestra función de Error, evaluada en p1 y con argumentos xx,yy(matriz modelo) y st1(primera estrella).
```python
p1 = params
best1,suss = leastsq(ErrorGauss2D, p1, args=(xx,yy,st1))
print(best1)
```
[149.49087839 116.36309968 4.88921373 32.41971607 12.85314581]
Solo queda graficar nuestra gaussiana optimizada ("st1_model_opt"), para eso, usamos la función modelo evaluada en los valores de "best1" y en el meshgrid de xx,yy.
```python
st1_model_opt = gauss2D(best1,xx,yy)
plt.imshow(st1_model_opt,cmap="gray")
```
<matplotlib.image.AxesImage at 0x7fc82dff9c88>
![png](output_24_1.png)
**Trabajando con la segunda estrella:**
Recortamos la estrella que deseamos trabajar como fué antes explicado y asignamos los pixeles de dicha estrella como un array a st2 (np.asarray())
```python
#Recortando la segunda estrella
estrella2 = gray[360:400,430:460]
plt.imshow(estrella2, cmap="gray")
st2 = np.asarray(estrella2)
print(st2.shape)
```
(40, 30)
![png](output_26_1.png)
Definiendo los valores modelo de nuestra segunda estrella (Gaussiana 2D). Creamos un arreglo de (30,40) que simule nuestra segunda estrella. Usando x2, y2, xx2, yy2 para no confudir valores con la estrella anterior.
```python
x2 = np.arange(0,30,1)
y2 = np.arange(0,40,1)
xx2 , yy2 = np.meshgrid(x2,y2)
```
También generamos una segunda lista de parametros para la estrella en cuestión (params2). Y evaluamos nuesta función gaussiana bajo el nombre de z2, en donde, aplicamos a la función los valores de interes para esta estrella.
```python
params2 = [5,5,3.5,20,17]
zz2 = gauss2D(params2,xx2,yy2)
plt.imshow(zz2,cmap="gray")
```
<matplotlib.image.AxesImage at 0x7fc82ea16630>
![png](output_30_1.png)
Definimos un "segundo chute" (p2) que estára definido por "params2". Calculamos los mejores valores posibles con "best2" e imprimimos dichos valores.
```python
p2 = params2
best2,suss = leastsq(ErrorGauss2D, p2, args=(xx2,yy2,st2))
print(best2)
```
[130.78142434 108.14028778 2.20724331 20.33084866 16.81819746]
Ahora graficamos nuestra segunda estrella ("st2_model_opt") aplicando los valores de "best2" a la función gaussiana junto con los parametros xx2,yy2.
```python
st2_model_opt = gauss2D(best2,xx2,yy2)
plt.imshow(st2_model_opt,cmap="gray")
```
<matplotlib.image.AxesImage at 0x7fc82ef02438>
![png](output_34_1.png)
**Trabajando con la tercera estrella:**
Los pasos explicados anteriormente aplican para la tercera estrella. De donde se hara un array de (30,20), definiremos un modelo con zz3 que dependerá de "params3" además de xx3,yy3.
```python
#Tercera estrella
estrella3 = gray[230:260,210:230]
plt.imshow(estrella3, cmap="gray")
st3 = np.asarray(estrella3)
print(st3.shape)
```
(30, 20)
![png](output_36_1.png)
```python
#Definiendo los valores modelo de nuestra tercera estrella (Gaussiana 2D)
x3 = np.arange(0,20,1)
y3 = np.arange(0,30,1)
xx3 , yy3 = np.meshgrid(x3,y3)
```
```python
params3 = [3,4,2,8,13]
zz3 = gauss2D(params3,xx3,yy3)
plt.imshow(zz3,cmap="gray")
```
<matplotlib.image.AxesImage at 0x7fc82ebb0f28>
![png](output_38_1.png)
Damos como siempre el "chute inicial" para la tercera estrella (p3), calculamos los mejores valores posibles ("best3") y finalmente graficamos nuestra gaussiana optimizada.
```python
p3 = params3
best3,suss = leastsq(ErrorGauss2D, p3, args=(xx3,yy3,st3))
print(best3)
```
[131.64909249 112.00306687 1.71925599 7.88240793 13.57465711]
```python
st3_model_opt = gauss2D(best3,xx3,yy3)
plt.imshow(st3_model_opt,cmap="gray")
```
<matplotlib.image.AxesImage at 0x7fc82eb1a9e8>
![png](output_41_1.png)
Como ultimo procedimiento de este proyecto, calcularemos la anchura a media altura(FWHM) para cada una de nuestras estrellas tratadas.
Para eso necesitamos la desviación estandar contenida en nuestras graficas optimizadas, para eso usaremos los valores contenido en las funciones "best".
Definimos la función FWHM:
```python
def FWHM(f):
FWHM = 2*np.sqrt(2*np.log(2))*f
return FWHM
```
Presentamos asi los FWHM de cada estrella 1, 2 y 3 respectivamente.
```python
FWHM(best1[2])
```
11.513218489642778
```python
FWHM(best2[2])
```
5.197660781112112
```python
FWHM(best3[2])
```
4.048538475508802
Como antesala a el proyecto antes mencionado, trabajamos con bandas de estrellas y generando gaussianas de una dimesión(1D). Se adjunta lo realiado pero no se explicará a detalle debido a que en escencia es el mismo codigo que se vio anteriormente.
**Para una banda de la primera estrella:**
```python
#Gaussiana de 1D. Intensidad Luminosa
#Cortando Las filas centrales de la estrella:
frag_star = estrella1[12:15,0:70]
print(frag_star.shape)
plt.imshow(frag_star, cmap="gray")
```
(3, 70)
<matplotlib.image.AxesImage at 0x7fc82e47a0f0>
![png](output_50_2.png)
Como ahora **frag_star** representa una matriz de (3,70) no es suficiente solo hacer un arreglo lineal de las variables x e y (como en el caso anteior). Por lo tanto, Y será un arreglo de 210 valores (producto de filas por columnas) y X fue diseñado como un arreglo de 3 matrices de 70 elementos con el fin de obtener 210 elementos para comparar con Y.
Graficamos xf vs Y y obtenemos todos los valores posibles contenidos en la banda de imagen.
```python
y = frag_star.reshape(210)
x = np.arange(70)
xf = np.array([x,x,x]).reshape(210)
plt.plot(xf,y,".")
```
[<matplotlib.lines.Line2D at 0x7fc82e30f400>]
![png](output_52_1.png)
```python
#Definimos la función modelo:
def gauss1D(params, x):
exponente = -((x-params[1])**2)/(2*(params[2])**2)
ymodel = params[0]*np.exp(exponente) + params[3]
return ymodel
#Params[0]: Constante de amplitud
#Params[1]: Xo
#Params[2]: desv. std
#Params[3]: constante aditiva
```
```python
#Definimos la función error:
def Error(tpl,x,y):
ymodel = gauss1D(tpl,x)
errors = (ymodel-y)
return errors
params1D =[140,32,4,120]
p0 = params1D
best,suss = leastsq(Error, p0, args=(xf,y))
print(best)
```
[145.52246647 32.43940401 4.98444533 117.26091299]
```python
ymodel_opt = gauss1D(best,xf)
plt.title("Gaussiana optimizada de la banda de la estrella 1")
plt.plot(xf,ymodel_opt,'r.')
plt.plot(xf,y,"b.")
```
[<matplotlib.lines.Line2D at 0x7fc82dfd9208>]
![png](output_55_1.png)
**Para una banda de la segunda estrella:**
```python
#Recortando las filas medias de la estrella 2
frag_star2 = estrella2[16:19,0:30]
print(frag_star2.shape)
plt.imshow(frag_star2, cmap="gray")
```
(3, 30)
<matplotlib.image.AxesImage at 0x7fc82de17a90>
![png](output_57_2.png)
```python
y2 = frag_star2.reshape(90)
x2 = np.arange(30)
xf2 = np.array([x2,x2,x2]).reshape(90)
plt.plot(xf2,y2,".")
```
[<matplotlib.lines.Line2D at 0x7fc82ddf19b0>]
![png](output_58_1.png)
```python
params1D_2 = [145,20,2.7,105]
p1 = params1D_2
best2,suss = leastsq(Error, p1, args=(xf2,y2))
print(best2)
```
[128.16071095 20.38057866 2.04737263 109.61442948]