Ejercicios para practicar numpy y optimización con scipy

LA-CoNGA 2021

Ejercicio elaborado por Eduardo E. Escalante R.

TAREA 3 (ejercicios-clase-05-datos)

• 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

• 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)

1. La conversión de una imagen de color a escala de grises se puede hacer de varias formas.

El pixel $(1,1,1)$ representa el blanco y el $(0,0,0)$ el negro en RGB. Entonces, cuanto más grandes son los valores de los canales, más "blanco" es el pixel y viceversa.

Por eso, una forma simple para hacer la conversión consiste en sacar el promedio de los canales y usar eso como la intensidad. De hecho, esto es lo que implementaremos en nuestras celdas jupyter:

1. Las imágenes se codifican como matrices. En particular, las imágenes de intensidad o escala de grises se codifican como una matriz de dos dimensiones, donde cada número representa la intensidad de un pixel. Las imágenes se codifican como matrices. En particular, las imágenes de intensidad o escala de grises se codifican como una matriz de dos dimensiones, donde cada número representa la intensidad de un pixel. Lo que significa que cualquiera de estas matrices que generamos se puede visualizar como una matriz. Para visualizar imágenes, usamos el módulo pyplot de la librería matplotlib.

# importamos el modulo pyplot, y lo llamamos plt
import matplotlib.pyplot as plt

#el modulo matplotlib.colors nos permite transformar el espacio de color de una imagen
#import matplotlib.colors

#configuracion necesaria de pyplot para ver las imagenes en escala de grises
plt.rcParams['image.cmap'] = 'gray'

# comando jupyter para que las imagenes se muestren automaticamente
%matplotlib inline

# también importamos Numpy ya que la usaremos para crear y manipular matrices
import numpy as np

# Para leer imágenes utilizaremos el módulo io de la librería skimage. 
# En ese módulo tenemos la función imread para leer imágenes de un archivo
# (y también la función imsave para escribirlas a un archivo).
from skimage import io


from scipy import ndimage

from scipy.optimize import leastsq

import time

import cv2
# La función que nos permite visualizar matrices es imshow del modulo pyplot,
# que invocaremos como plt.imshow, y recibe la imagen como parámetro.

2. Busqueda y visualización de la imagen

2.1. Primero ubiquemos y visualicemos la imagen con la función imshow del modulo pyplot que nos permite visualizar matrices, que invocaremos como plt.imshow, y recibe la imagen como parámetro. Debido a que pyplot intenta ajustar los colores automaticamente de una imagen, vamos a pasarle como parámetros también vmin=0,vmax=1. Estos se conocen como keyword parameters o parámetros por palabra clave de python.

2.2. Para leer imágenes utilizaremos el módulo *io* de la librería `skimage`. En ese módulo disponemos de la función `imread` para leer imágenes de un archivo (y también la función `imsave` para escribirlas a un archivo) esta función nos permite leer imagenes tanto en blanco y negro como en color.

2.3. En este caso, tenemos tres dimensiones. La primera y segunda de tamaños 789, 1184 corresponden a la cantidad de pixels de la imagen. La última, de tamaño 3, contiene los tres canales *R (red), G (green) y B (blue)* de la imagen.

Entonces cada trio $(r,g,b)$ que representa un pixel se reemplaza por su promedio $(r+g+b)/3$

# imread lee las imagenes con los pixeles codificados como enteros 
# en el rango 0-255. Por lo que la convertiremos a flotante y en el rango 0-1

image = io.imread("./data/zapatocaImage.jpeg")/255.0
print("- Dimensiones de la imagen:")
print(image.shape)
plt.imshow(image)

- Dimensiones de la imagen:
(789, 1184, 3)

<matplotlib.image.AxesImage at 0x7fb72eb8ba90>

Podemos ver cada canal individualmente, como si el canal fuera una imagen en escala de grises, para ver la intensidad de cada canal.

#Anulo los canales rojo y azul para ver solo el verde
green=np.copy(image)
green[:,:,0]=0
green[:,:,2]=0
plt.figure()
plt.imshow(green)

#Anulo los canales azul y verde para ver solo el rojo
red=np.copy(image)
red[:,:,1:3]=0
plt.figure()
plt.imshow(red)

#Anulo los canales rojo y verde para ver solo el azul
blue=np.copy(image)
blue[:,:,0:2]=0
plt.figure()
plt.imshow(blue)

<matplotlib.image.AxesImage at 0x7fb72a757358>

#Sacar los promedios de cada uno de los pixels que componen esta imagen

h, w, c = image.shape

print(h,w,c)

789 1184 3

# Crearemos una matriz de las mismas dimensiones (negra al principio) para construir
# la figura en escala de grises

imageGris = np.zeros((h,w))

plt.imshow(imageGris)

<matplotlib.image.AxesImage at 0x7fb72a6f2240>

# Para convertir la imagen a escala de grises, recorreremos cada pixel de la matrix
# y tomaremos la intensidad de cada uno de los indices RGB en cada uno de los pixels
# para luego tomar el promedio y asignarlo a la matriz imageGris

for i in range(h):
    for j in range(w):
        imageGris[i,j] = (image[i,j,0]+image[i,j,1]+image[i,j,2])/3
plt.imshow(imageGris[:,:], cmap='gray')
print(imageGris.shape)

(789, 1184)

pyplot.figure crea una nueva figura con los atributos dados en los parámetros, donde figsize define el tamaño de la figura en pulgadas.

fig = plt.figure(figsize=(15,15))
ax1 = fig.add_subplot(2,2,1)
ax1.imshow(image)
ax2 = fig.add_subplot(2,2,2)
ax2.imshow(imageGris)

<matplotlib.image.AxesImage at 0x7fb72a691780>

#guardamos la imagen generada en escala de grises

io.imsave("./data/zapatocaImage_escala_de_grises.png",imageGris)

Lossy conversion from float64 to uint8. Range [0, 1]. Convert image to uint8 prior to saving to suppress this warning.

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

3. ¿Cómo recortar una imagen? Para recortar una imagen vamos a tratar a esta como una matriz, para ello emplearemos indexación de Numpy, de donde escogeremos las filas y columnas que deseemos recortar de la imagen.

Planteando algo como lo siguiente: imagen[fila inicial : fila final, columna inicial : columna final]

# Recortemos la imagen

imagen_recortada = imageGris[310:325,620:635]

print(imagen_recortada.shape)

#copia_Blanco_Negro = np.copy(imagen_recortada)

plt.imshow(imagen_recortada, cmap="gray")
copia_Blanco_Negro = np.copy(imagen_recortada)
print(imagen_recortada.shape)
plt.show()

(15, 15)
(15, 15)

io.imsave("./data/zapatocaRecortada.png",imagen_recortada)

Lossy conversion from float64 to uint8. Range [0, 1]. Convert image to uint8 prior to saving to suppress this warning.

#Se crea una gaussiana
def gauss2D(p,x,y):
    exponente = -((x-p[0])**2 + (y-p[1])**2) / (2*p[2]**2)
    z = p[3] * np.exp(exponente) + p[4]
    return z

#Se define el error que devuelve un array en 1 D

def Error(p,x,y,z):
    
    ymodel=gauss2D(p,x,y)
#    ymodel=gaussvec(p,x,y)
    error=ymodel-z
    return error.reshape(-1)

#Se toman los parametros iniciales

# x0 y0 c, a , b
p0=[12,12,10,1,0]

#Se crea un grilla de 15x15, correspondiente a las posiciones de cada pixel
x=np.arange(0,15,1)
y=np.arange(0,15,1)
yy,xx=np.meshgrid(x,y)

# Se realiza la optimizacion por minimos cuadrados
best,suss=leastsq(Error,p0,args=(xx,yy,copia_Blanco_Negro))

#Se verifica que se ajusto correctamente la gaussiana
fig, (ax1, ax2)=plt.subplots(1,2,figsize=(10,6))

ax1.imshow(gauss2D(best,xx,yy),cmap='gray')
ax1.set_title('Gaussiana')

ax2.imshow(copia_Blanco_Negro,cmap='gray')
ax2.set_title('Imagen')

Text(0.5, 1.0, 'Imagen')

4. Histograma de la imagen

Los histogramas son distribuciones que describen la frecuencia con las que se presentan los valores de intensidad (pixeles) de una imagen. Tomando como ejemplo una imagen típica a escala de grises de $8$ bits el histograma adoptará valores comprendidos entre $0$ y $255$ ya que $2^{8}=256$. Cada valor del histograma contiene el numero de pixeles correspondientes a esa intensidad.

Evidentemente el histograma no proporciona información acerca del origen de los pixeles que lo conforman. Esto significa que se pierde la informacion de la relación espacial que tenian los pixeles en la imagen, por lo que resulta imposible reconstruir una imagen a partir de su histograma.

El histograma muestra importantes caracteristicas de una imagen:

Por ejemplo el contraste y el rango dinamico, los cuales son problemas atribuidos a la captación de la imagen que es necesario resolver para que puedan ser analizados con mejor claridad propiedades de la imagen en bloques de procesamiento posterior; los errores de iluminación son reconocidos en el histograma, por lo que la region final o inicial de la escala de intensidad no es utilizada, mientras que las otras regiones del histograma son ocupadas por valores de intensidad de pixeles

Como contraste se entiende el campo de valores de intensidad de una imagen concreta que son utilizados, es decir, una imagen con contraste completo utiliza un intervalo completo de niveles de intensidad definidos para la imagen y bajo el termino de dinámica se entiende el número de pixeles diferentes que son utilizados en la imagen. El caso ideal resulta cuando es utilizado el rango completo de valores de intensidad disponibles, utilizando asi su maximo rango dinamico. Una dinamica alta representa una ventaja ya que se reduce el riesgo de perder calidad a traves de las siguientes etapas del procesamiento. Por esta razon las camaras digitales y los escaner profesionales tienen una resolucion mayor a los 8 bits, normalmente de 12 a 14 bits, aunque los elementos de despliegue tengan resoluciones normales con 256 valores de intensidad.

im = cv2.imread("./data/zapatocaRecortada.png")
# Calculamos el valor medio de los canales RGB
vals = im.mean(axis=2).flatten()
# Graficamos el histograma con 255 bins
b, bins, patches = plt.hist(vals, 255)
plt.xlim([0,255])
plt.show()

Tomemos una 2da estrella y hagamos el mismo tratamiento

# Recortemos la imagen

imagen_recortada = imageGris[540:565,655:670]

print(imagen_recortada.shape)

#copia_Blanco_Negro = np.copy(imagen_recortada)

plt.imshow(imagen_recortada, cmap="gray")
copia_Blanco_Negro = np.copy(imagen_recortada)
print(imagen_recortada.shape)
plt.show()

(25, 15)
(25, 15)

io.imsave("./data/zapatocaRecortada-2.png",imagen_recortada)

#Se crea una gaussiana
def gauss2D(p,x,y):
    exponente = -((x-p[0])**2 + (y-p[1])**2) / (2*p[2]**2)
    z = p[3] * np.exp(exponente) + p[4]
    return z

#Se define el error que devuelve un array en 1 D

def Error(p,x,y,z):
    
    ymodel=gauss2D(p,x,y)
#    ymodel=gaussvec(p,x,y)
    error=ymodel-z
    return error.reshape(-1)

#Se toman los parametros iniciales

# x0 y0 c, a , b
p0=[15,15,10,1,0]

#Se crea un grilla de 15x15, correspondiente a las posiciones de cada pixel
x=np.arange(0,15,1)
y=np.arange(0,25,1)
yy,xx=np.meshgrid(x,y)

# Se realiza la optimizacion por minimos cuadrados
best,suss=leastsq(Error,p0,args=(xx,yy,copia_Blanco_Negro))

#Se verifica que se ajusto correctamente la gaussiana
fig, (ax1, ax2)=plt.subplots(1,2,figsize=(10,6))

ax1.imshow(gauss2D(best,xx,yy),cmap='gray')
ax1.set_title('Gaussiana')

ax2.imshow(copia_Blanco_Negro,cmap='gray')
ax2.set_title('Imagen')

Text(0.5, 1.0, 'Imagen')

im = cv2.imread("./data/zapatocaRecortada-2.png")
# calculate mean value from RGB channels and flatten to 1D array
vals = im.mean(axis=2).flatten()
# plot histogram with 255 bins
b, bins, patches = plt.hist(vals, 255)
plt.xlim([0,255])
plt.show()