Skip to content
Snippets Groups Projects
Select Git revision
  • master default protected
1 result
Blame History Permalink
Forked from an inaccessible project.
Entrega.md 14.30 KiB

Tareas Clase 5

Ejercicio - Resolución espacial

Nombre: Sasiri Juliana Vargas Urbano
Correo: sasiri.vargas@correounivalle.edu.co
Usuario Mattermost: vargass

Objetivo de la tarea: 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.

Se importan las librerias

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.optimize import leastsq    # Para optimizar las funciones
import statistics as stat
import math

NOTA:

Por facilidad, se creó una clase para el posterior análisis de todas las estrellas, las funciones dentro de la clase tienen como base el funcionamiento de las funciones presentadas a lo largo de este notebook

from Clase_estrella import *

Paso 1 - Leer la imagen almacenada

#Para abrir la información de la imagen

imagen = plt.imread('./data/zapatocaImage.jpeg')

#Son 3 matrices de dimensiones 789 x 1184
print(imagen.shape) 

(789, 1184, 3)

Imagen a tratar:

#Impresion de la imagen original

plt.imshow(imagen)
<matplotlib.image.AxesImage at 0x7f38c3d63940>

png

#Canales RGB

#Según la forma en la que se declaran los arreglos, se pueden formar diferentes matrices
#que representan los canales RGB

R = imagen[:,:,0] #Red
G = imagen[:,:,1] #Green
B = imagen[:,:,2] # Blue

Para entender la composición de colores RGB de la anterior imagen, se grafica cada una de las matrices (R,G y B), haciendo las demás cero, de esta forma se obtiene:

imagen_copy = np.copy(imagen)   #Se copia la imagen para no perder la original

#Para el canal Rojo ---------------------------------
imagen_copy[:,:,0]=imagen[:,:,0]
imagen_copy[:,:,1]=0
imagen_copy[:,:,2]=0

Rojo = imagen_copy

imagen_copy = np.copy(imagen)

#Para el canal Verde ---------------------------------
imagen_copy[:,:,0]=0
imagen_copy[:,:,1]=imagen[:,:,1]
imagen_copy[:,:,2]=0

Verde = imagen_copy

imagen_copy = np.copy(imagen)

#Para el canal Azul ---------------------------------
imagen_copy[:,:,0]=0
imagen_copy[:,:,1]=0
imagen_copy[:,:,2]=imagen[:,:,2]

Azul = imagen_copy

#Impresion de imágenes

fig, (ax1, ax2, ax3) = plt.subplots(1,3, figsize=(16,9))
ax1.imshow(Rojo)
ax1.set_title('Canal Rojo')
ax2.imshow(Verde)
ax2.set_title('Canal Verde')
ax3.imshow(Azul)
ax3.set_title('Canal Azul')
Text(0.5, 1.0, 'Canal Azul')

png

Paso 2 - Escala de grises

Existen distintos métodos por los cuales podemos encontrar la imagen a escala de grises

Método del promedio

Se realiza el promedio por cada pixel y se encuentra la luminosidad que luego se reescala con una escala de grises

#Método del promedio
imagen_gris=(R+G+B)/3

#NOTA: Por prueba y error, encontré que multiplicando por 0.1 a R, se encuentra una mejor visualización de la imagen
#(Esto está relacionado con el siguiente método)
imagen_griss=(R*0.1+G+B)/3

fig, (ax1, ax2) = plt.subplots(1,2, figsize=(16,9))
ax1.imshow(imagen_gris, cmap='gray')
ax1.set_title('Método del promedio')
ax2.imshow(imagen_griss, cmap='gray')
ax2.set_title('Modificación método del promedio')
Text(0.5, 1.0, 'Modificación método del promedio')

png

Método pesado o de luminosidad

Se pesa la cantidad de R,G y B, con parámetros tales que, cuando se tiene una superposición de los R, G, B pesados, el ojo humano ve en escala de grises.

#Método pesado
imagen_grisss = (0.299*R+0.587*G+0.114*B)

fig=plt.figure(figsize=(16,5))
plt.imshow(imagen_grisss, cmap='gray')
<matplotlib.image.AxesImage at 0x7f38c3c5b198>

png

Las imágenes cirrespondientes a R,G,B en escala de grises se verían como:

fig, (ax1, ax2, ax3) = plt.subplots(1,3, figsize=(16,9))
ax1.imshow(R, cmap='gray')
ax1.set_title('Canal Rojo')
ax2.imshow(G, cmap='gray')
ax2.set_title('Canal Verde')
ax3.imshow(B, cmap='gray')
ax3.set_title('Canal Azul')
Text(0.5, 1.0, 'Canal Azul')

png

Se pueden observar unas ciertas diferencias entre las intensidades de las estrellas y por lo tanto, para el análisis de las estrellas, pueden existir diferencias para los parámetros de ajuste.

Paso 3 - Recorte de estrellas y ajuste gausiano 2D

Se escoje el rango de los pixeles de una estrella con prueba y error En primera instancia se trata de escogerla centrada

estrella = imagen_grisss[535:570, 645:680]

plt.imshow(estrella, cmap='gray')
plt.title('Primera estrella')
Text(0.5, 1.0, 'Primera estrella')

png

Para efectos del análisis posterior, también se recorta un fragmento con varias estrellas.

estrella_2 = imagen_grisss[330:380, 520:570]

plt.imshow(estrella_2, cmap='gray')
plt.title('Segunda estrella')
Text(0.5, 1.0, 'Segunda estrella')

png

Para el primer análisis de optimización se toma la 'Primera estrella'.

Para utilizar la función leastsq de scipy para optimizar se necesitan de las siguientes funciones:

#Definición del error, para una función gaussiana
#Error = Valor teórico - Valor experimental

def residuos(p, yy, xx, zz):
    
    a, b, c, x0, y0 = p
    
    exponente = -((xx-x0)**2 + (yy-y0)**2) / (2*c**2)
    
    error  = zz - (a * np.exp(exponente) + b)
    
    return error.reshape(-1)    #Se retorna un arreglo 1D, ya que esto es lo que recibe leastsq (función de scipy)

#NOTA: Para mayor información de esta función, pedir help() dentro de la clase Estrella_a
#Definición de una función que calcula una Gaussiana 2D
def gauss2D(x, y, p):
    
    a, b, c, x0, y0 = p
    
    exponente = -((x-x0)**2 + (y-y0)**2) / (2*c**2)
    z = a * np.exp(exponente) + b
    
    return z

#NOTA: Para mayor información de esta función, pedir help() dentro de la clase Estrella_a

Para el ajuste, se hacen dos arreglos de prueba x y y de dimensiones igual a 'estrella' (la imagen a analizar)

#Arreglos de prueba
x = np.arange(0,255,7.3)
y = np.arange(0,255,7.3)

#Grilla
yy, xx = np.meshgrid(x,y)

#Parámetros de ajuste tentativos
pe=np.array([2,15,15,5,5])

Haciendo uso de todo lo anterior, se obtienen los valores ajustados por míminos cuadrados de los parámetros de ajuste, esto se hace con la función leastsq de scipy

ajuste, dim = leastsq(residuos, pe, args=(xx, yy, estrella))

print(ajuste)
[149.8213941  119.16039711 -34.31135287 127.16003106 130.26587992]

Visualización del ajuste y la estrella escogida

zz= gauss2D(xx,yy,ajuste)

fig, (ax1, ax2) = plt.subplots(1,2, figsize=(9,6))
ax1.imshow(estrella, cmap='gray')
ax1.set_title('Estrella a analizar')
ax2.imshow(zz,cmap="gray")
ax2.set_title('Ajuste')
Text(0.5, 1.0, 'Ajuste')

png

Estadística con más estrellas

Para el análisis de más estrellas se crea una función que busca estrellas analizando una imagen de entrada.

Para ello se crean dos funciones, miniframe_estrella y estrellas_recortadas

def miniframe_estrella(estrellita_, fila, z_):
    
    """
    
    Encuentra el arreglo de la imagen de la estrella recortada.
    
    ---------------------------------------------------
    
    FUNCIONAMIENTO:
    
    estrellita: Imagen a analizar
    
    fila: fila de estrellita donde se encuentra la estrella
    
    z_: Arreglo donde se guardan las posibles estrellas
    """
    
    ancho=[]

    nx, ny = np.shape(estrellita_)
    
    #El siguiente for recorre la imagen ingresada, en la fila seleccionada
    #Añade elementos a la lista ancho, mientras que la intensidad de la estrella sea mayor que 0.8
    #Cuando se encuentra el centro de la estrella, el for para
    #Conociendo la cantidad de elementos añadido a la lista ancho, se encuentra el radio de la estrella
    
    for j in range(ny):
        if (estrellita_[z_[0,fila],j] > 0.8) and (j != z_[1,fila]):
            ancho.append(1)
            
    radio = len(ancho)
            
    #Creacion de un miniframe que contiene solo la estrella a analizar

    star = estrellita_[z_[0,fila]-len(ancho)-5:z_[0,fila]+len(ancho)+5, z_[1,fila]-len(ancho)-5:z_[1,fila]+len(ancho)+5]
    
    #star = estrellita_[z_[0,fila]-10:z_[0,fila]+10, z_[1,fila]-10:z_[1,fila]+10]
    
    #Retorna el array correspondiente a la imagen del miniframe
    
    return radio, star

#Intensidad recomendada: 240

def estrellas_recortadas(img_es, intensidad):
    
    """
    
    Retorna un arreglo, tal que cada una de sus componentes es la imagen recortada de las estrellas identificadas en la imagen de ingreso
    
    ---------------------------------------------------
    
    FUNCIONAMIENTO:
    
    img_es: Imagen a analizar
    
    intensidad: Intensidad de estrellas deseada para analizar. Parámetro de 0 a 254. Intensidad recomendada: 240
    
    """
    
    #Primero se escala con la intensidad de estrellas deseada para buscar
    
    estrellita=img_es/intensidad
    
    #Se identifican las candidatas a estrellas
    
    z = np.where(estrellita > 1)

    z=np.array(z)

    n, nz = np.shape(z)
    
    k=0
    
    #Se crea una lista para las estrellas de salida

    lista_estrellas = []
    
    #El ciclo a continuación funciona de la siguiente manera:
    #Se hace la diferencia en valor absoluto de una candidata, con la candidata siguiente. Se obtiene un arreglo 2x2
    #Si la suma del anterior arreglo es menor o igual a dos, significa que las candidatas a estrellas son contiguas y por consiguiente se asume que es la misma estrella
    #Si la suma del anterior arreglo es mayor, entonces se continuan evaluando las demás estrellas con el mismo procedimiento
    #Los frames de estrellas se almacenan en la lista: 'lista_estrellas'

    while k in range(nz-2):
        
        # 2. Es decir que vuelve aquí*
        
        k += 1

        r = abs(z[:,k]-z[:,k+1])

        r_sum = r.sum()

        if r_sum <= 2:

            radio, miniestrella = miniframe_estrella(estrellita, k, z)
            
            lista_estrellas.append(miniestrella)

        else:
            
            continue # 1. Continua con el ciclo*
    
    #lista_estrellas = np.array(lista_estrellas)
    return lista_estrellas

Utilizando la imagen de 'Estrella 2' para el análisis, se encuentran las siguientes entrellas dentro.

lista_estrellas=estrellas_recortadas(estrella_2, 240)



fig, axs = plt.subplots(2, 2, figsize=(6,6))

axs[0, 0].imshow(lista_estrellas[0], cmap='gray')
axs[0, 0].set_title('Estrella 1')
axs[0, 1].imshow(lista_estrellas[1], cmap='gray')
axs[0, 1].set_title('Estrella 2')
axs[1, 0].imshow(lista_estrellas[2], cmap='gray')
axs[1, 0].set_title('Estrella 3')
axs[1, 1].imshow(lista_estrellas[3], cmap='gray')
axs[1, 1].set_title('Estrella 4')

# Hide x labels and tick labels for top plots and y ticks for right plots.
for ax in axs.flat:
    ax.label_outer()

png

Para el análisis de las anteriores estrellas se hace uso de la clase, se presenta además el ajuste gaussiano de dichas estrellas

#Conversion a la clase estrella

clase_estrella=[]

for item in lista_estrellas:
    clase_estrella.append(Estrella_a(item))

#Impresión de las imagenes

#Para recordar: p = [a, b, c, x0, y0]


p1=np.array([1,0,1,5,5])

graficas=[]

parametros=[]


for i in range(0,len(clase_estrella)):
    uno, dos = clase_estrella[i].ajusteGauss(p1)
    
    graficas.append(uno)
    
    parametros.append(dos)


#Impresion de imagenes

fig, axs = plt.subplots(2, 2, figsize=(6,6))

<<<<<<< HEAD
axs[0, 0].imshow(graficas[0], cmap='gray')
axs[0, 0].set_title('Estrella 1')
axs[0, 1].imshow(graficas[1], cmap='gray')
axs[0, 1].set_title('Estrella 2')
axs[1, 0].imshow(graficas[2], cmap='gray')
axs[1, 0].set_title('Estrella 3')
axs[1, 1].imshow(graficas[3], cmap='gray')
axs[1, 1].set_title('Estrella 4')

# Hide x labels and tick labels for top plots and y ticks for right plots.
for ax in axs.flat:
    ax.label_outer()

png

Análisis estadistico

Para encontrar la mediana, media, moda y desviacion estandar del ajuste hecho se hace lo siguiente:

evaluar=np.array(parametros)

evaluar=evaluar*2*math.sqrt(2*math.log(2))

zz=evaluar
#Media
media=np.mean(zz)

#Mediana
mediana=np.median(zz)

#Moda
moda=stat.mode(np.round(zz,3).reshape(-1))

#Desviacion estandar con scipy
desviacion=stat.stdev(zz.reshape(-1))

#Desviacion estandar con numpy
desviacion1=zz.std()

print('Mediana: ', mediana)
print('Media: ', media)
print('Moda: ', moda)
print('Desviacion estandar 1: ', desviacion)
print('Desviacion estandar 2: ', desviacion1)
Mediana:  2.9331335083105
Media:  6.958263965191911
Moda:  15.497
Desviacion estandar 1:  6.4437838006164565
Desviacion estandar 2:  6.280623550702009

Para hacer el histograma

#Histograma

plt.figure(figsize=(9,4))
plt.hist(zz, bins=5, histtype='bar', alpha=0.7, edgecolor = 'black', linewidth=0.2)
plt.show()

png

estrellas_para_estadistica = estrellas_recortadas(imagen_grisss, 230)

#Para guardar el notebook a .md
! jupyter nbconvert --to markdown Entrega.ipynb
[NbConvertApp] Converting notebook Entrega.ipynb to markdown
[NbConvertApp] Support files will be in Entrega_files/
[NbConvertApp] Making directory Entrega_files
[NbConvertApp] Making directory Entrega_files
[NbConvertApp] Making directory Entrega_files
[NbConvertApp] Making directory Entrega_files
[NbConvertApp] Making directory Entrega_files
[NbConvertApp] Making directory Entrega_files
[NbConvertApp] Making directory Entrega_files
[NbConvertApp] Making directory Entrega_files
[NbConvertApp] Making directory Entrega_files
[NbConvertApp] Making directory Entrega_files
[NbConvertApp] Making directory Entrega_files
[NbConvertApp] Writing 13781 bytes to Entrega.md

=======

0f2da7af