Commit f54778ba authored by Jennifer Lorena Ortega Aguilar's avatar Jennifer Lorena Ortega Aguilar
Browse files

BPF

parent d132a31b
%% Cell type:code id: tags:
``` python
import cv2
import numpy as np
from matplotlib import pyplot as plt
from PIL import Image
from PIL.ExifTags import TAGS
from numpy import asarray
from numpy.fft import fft2, fftshift, fftfreq, ifft2 #Python DFT
```
%% Cell type:code id: tags:
``` python
#imagename = "Data/crateres.jpg"
#imagename = "Data/Perseverance_ZR0_0074.png" -codigo adicional para recorte de bordes negros.
#imagename = "Data/mastcamz_2_prints.jpg"
imagename = "Data/Data_Mars_Mars8.png"
imagename = "Data_Mars/Mars8.png"
# read the image data using PIL
image = Image.open(imagename)
```
%% Cell type:code id: tags:
``` python
display (image)
```
%% Cell type:markdown id: tags:
crateres.jpg : Image downloaded from : https://mars.nasa.gov/mars2020-raw-images/pub/ods/surface/sol/00000/ids/edr/browse/edl/ELM_0000_0666952859_600ECM_N0000030LVS_04000_0000LUJ00_800.jpg
Perseverance_ZR=_0074.png : Image downloaded from :
%% Cell type:code id: tags:
``` python
img = cv2.imread(imagename, 0) # load an image
#3 = np.asarray(imagename, dtype = np.float32) # Image class instance, I1, to float32 Numpy array, a
img_float32 = np.float32(img)
```
%% Cell type:markdown id: tags:
Para procesar la imagen se lee como cv2.imread y luego se le aplica la DFT
DFT: Discrete Fourier Transform.
%% Cell type:code id: tags:
``` python
hist =cv2.calcHist([img],[0],None,[256],[0,256])
plt.plot(hist,color='blue')
plt.xlabel('Intensidad de Iluminacion')
plt.ylabel('Pixels Bins')
plt.show()
plt.hist(np.concatenate((img),axis=0),bins=250)
plt.xlabel('Intensidad de Iluminacion')
plt.ylabel('Pixels Bins')
plt.show()
```
%% Cell type:code id: tags:
``` python
#transform,ar el array de la imagen en punto flotante para ser procesado.
dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)
dft_shift =np.fft.fftshift(dft)
```
%% Cell type:code id: tags:
``` python
# convert image to numpy array
data = asarray(image)
print(type(data))
# caracteristicas de la imagen
print(data.shape)
```
%% Cell type:code id: tags:
``` python
#Muestra cada valor de cada pixel de la imagen como array en Numpy
#print(data)
```
%% Cell type:code id: tags:
``` python
#imagen en el espacio de las frecuencias :
f = np.fft.fft2(img) #fft2 para 2 dimensiones
fshift = np.fft.fftshift(f) # para centrar las frecuencias bajas.
# fft2 tiene valors imaginarios por ello se saca el espectro real o de magnitud.
magnitude_sprectrum = 20*np.log(np.abs(fshift))
```
%% Cell type:code id: tags:
``` python
f.shape
```
%% Cell type:code id: tags:
``` python
fft2(f)
```
%% Cell type:code id: tags:
``` python
#Matplotlib para crear 2 subplots en la misma linea para comparar la imagen original y
#la imagen transformada en el espectro de frecuencias.
plt.subplot(121),plt.imshow(img , cmap = 'gray')
plt.title('Original Image'), plt.xticks([]),plt.yticks([])
plt.subplot(122),plt.imshow(magnitude_sprectrum , cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]),plt.yticks([])
plt.show
```
%% Cell type:code id: tags:
``` python
#Filtro pasa alto HPF en el dominio de las frencuencias para demarcar contornos:
rows, cols = img.shape
crow, ccol = int(rows/2), int(cols/2)
#remover las frecuencias bajas usando una plantilla circular
mask=np.ones((rows,cols,2),np.uint8)
r = 40
center = [crow, ccol]
x, y = np.ogrid[:rows, :cols]
mask_area = (x - center[0]) ** 2 + (y - center[1]) ** 2 <= r*r
mask[mask_area] = 1
fshift=dft_shift*mask
fshift_mask_mag=2000*np.log(cv2.magnitude(fshift[:,:,0],fshift[:,:,1]))
f_ishift =np.fft.ifftshift(fshift)
img_back =cv2.idft(f_ishift)
img_back =cv2.magnitude(img_back[:,:,0],img_back[:,:,1])
#fshift[crow-30:crow+30,ccol-30:ccol+30]=0
f_ishift = np.fft.ifftshift(fshift)
img_back = np.fft.ifft2(f_ishift)
img_back = np.abs(img_back)
```
%% Cell type:code id: tags:
``` python
#Filtro pasa Banda BPF aplicando dos circulos concentricos
rows, cols = img.shape
crow, ccol = int(rows/2), int(cols/2) # circulo centro
#remover las frecuencias bajas usando una plantilla circular
rows, cols = img.shape
crow, ccol = int(rows/2), int(cols/2)
mask=np.zeros((rows,cols,2),np.uint8)
r_out = 40
r_in = 4
center = [crow, ccol]
x, y = np.ogrid[:rows, :cols]
mask_area = np.logical_and(((x-center[0])**2 +(y-center[1])**2>= r_in**2),
((x-center[0])**2 +(y-center[1])**2<= r_out**2))
mask[mask_area] = 1
fshift=dft_shift*mask
fshift_mask_mag= 2000 * np.log(cv2.magnitude(fshift[:, :, 0], fshift[:, :, 1]))
f_ishift =np.fft.ifftshift(fshift)
img_back =cv2.idft(f_ishift)
img_back =cv2.magnitude(img_back[:,:,0],img_back[:,:,1])
#fshift[crow-30:crow+30,ccol-30:ccol+30]=0
fig = plt.figure(figsize=(20,20))
fig.add_subplot(2, 2, 1), plt.imshow(img, cmap='gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
fig.add_subplot(2, 2, 2), plt.imshow(magnitude_sprectrum,cmap='gray')
plt.title('After FFT'), plt.xticks([]), plt.yticks([])
fig.add_subplot(2, 2, 3), plt.imshow(fshift_mask_mag, cmap='gray')
plt.title('FFT + Mask'), plt.xticks([]), plt.yticks([])
fig.add_subplot(2, 2, 4), plt.imshow(img_back, cmap='gray')
plt.title('After FFT Inverse'), plt.xticks([]), plt.yticks([])
plt.show()
```
%% Cell type:code id: tags:
``` python
#Filtro pasa Banda BPF aplicando dos circulos concentricos
rows, cols = img.shape
crow, ccol = int(rows/2), int(cols/2) # circulo centro
#remover las frecuencias bajas usando una plantilla circular
rows, cols = img.shape
crow, ccol = int(rows/2), int(cols/2)
mask=np.zeros((rows,cols,2),np.uint8)
#con diferentes radios tal que sean mayor que la zona central, tapando la mayor cantidad de bajas frencuencias.
r_out = 80
r_in = 4
center = [crow, ccol]
x, y = np.ogrid[:rows, :cols]
mask_area = np.logical_and(((x-center[0])**2 +(y-center[1])**2>= r_in**2),
((x-center[0])**2 +(y-center[1])**2<= r_out**2))
mask[mask_area] = 1
fshift=dft_shift*mask
fshift_mask_mag= 2000 * np.log(cv2.magnitude(fshift[:, :, 0], fshift[:, :, 1]))
f_ishift =np.fft.ifftshift(fshift)
img_back =cv2.idft(f_ishift)
img_back =cv2.magnitude(img_back[:,:,0],img_back[:,:,1])
fig = plt.figure(figsize=(20,20))
fig.add_subplot(121),plt.imshow(img , cmap = 'gray')
plt.title('Original Image'), plt.xticks([]),plt.yticks([])
fig.add_subplot(122),plt.imshow(img_back , cmap = 'gray')
plt.title('AFTER FFT Inverse image'), plt.xticks([]),plt.yticks([])
plt.show()
```
%% Cell type:code id: tags:
``` python
plt.subplot(131),plt.imshow(img, cmap='gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(132),plt.imshow(img_back , cmap = 'gray')
plt.title('Imagen con HPF'), plt.xticks([]),plt.yticks([])
plt.subplot(133),plt.imshow(img_back)
plt.title('Imagen Final'), plt.xticks([]),plt.yticks([])
plt.show
```
%% Cell type:code id: tags:
``` python
fig = plt.figure(figsize=(20,20))
fig.add_subplot(121),plt.imshow(img , cmap = 'gray')
plt.title('Original Image'), plt.xticks([]),plt.yticks([])
fig.add_subplot(122),plt.imshow(img_back , cmap = 'gray')
plt.title('HPF image'), plt.xticks([]),plt.yticks([])
plt.show
```
%% Cell type:code id: tags:
``` python
#Comparativa para la foto de las rocas
#recortado de foto
#Recortar una imagen
crop= img[200:400,200:500]
cv2.imshow('Seccion Imagen',crop)
```
%% Cell type:code id: tags:
``` python
#filtro Pasa bajos: LPF
dft_shift= np.fft.fftshift(dft)
rows, cols = img.shape
crow , ccol = int(rows/2), int(cols/2)
mask= np.zeros((rows,cols,2),np.uint8)
mask[crow-30:crow+30,ccol-30:ccol+30]=1
fshift = dft_shift*mask
f_ishift=np.fft.ifftshift(fshift)
img_back=cv2.idft(f_ishift)
img_back=cv2.magnitude(img_back[:,:,0],img_back[:,:,1])
```
%% Cell type:code id: tags:
``` python
fig = plt.figure(figsize=(20,20))
fig.add_subplot(121),plt.imshow(img , cmap = 'gray')
plt.title('Original Image'), plt.xticks([]),plt.yticks([])
fig.add_subplot(122),plt.imshow(img_back , cmap = 'gray')
plt.title('LPF image'), plt.xticks([]),plt.yticks([])
plt.show
```
%% Cell type:code id: tags:
``` python
```
......
%% Cell type:markdown id: tags:
# Data Analysis with images from Mars
## Nombres: Jennifer Ortega (EPN), Gerardo Semprum (UCV), Carla Gomez (USB)
## Segmentación de Imágenes de Marte
El presente proyecto realiza una segmentación de imágenes de Marte, adquiridadas por el rover Mars Perseverance de la NASA, mediante el algoritmo de watershed.
### Estructura del notebook
- Carga y visualización de la imagen
- Tranformación de la imagen de RGB a Escala de Grises
- Establecer umbral para la segmentación
- Cálculo de marcadores
- Aplicación del algoritmo Watershed
- Análisis
- Conclusiones
%% Cell type:markdown id: tags:
---
<font color='gray'>
### Imagen
Las imágenes digitales en color están hechas de píxeles. Un pixel es la menor unidad homogénea con color de la imagen, el pixel esta formado por combinaciones de colores primarios como son el **rojo**, **verde** y **azul** conocidos como canales **RGB**.
La imagen puede ser representada por una matriz de dimensiones NxM, donde cada elemento de la matriz es una matriz de 3 componentes, cada componente del elemento contiene los valores para los canales R, G y B de la imagen ([R G B]).
---
</font>
%% Cell type:markdown id: tags:
### Carga y visualización de la imagen
- Para cargar la imagen utilizamos la biblioteca **PIL** de Python diseñada para el manejo y procesamiento de imágenes.
- La imagen es guardada en formato **numpy.ndarray**, que representa una matriz multidimensional. La matriz es de dimensión 1200 x 1648, donde cada elemento de la matriz es un array de 3 componentes. Cada componente del elemento contiene los valores para los canales R, G y B de la imagen.
- La imagen se puede encontrar <a href="https://mars.nasa.gov/mars2020/multimedia/raw-images/ZL0_0072_0673337731_644EBY_N0032208ZCAM08032_1100LUJ">aquí.</a> El rover Mars Perseverance de la NASA adquirió esta imagen usando su cámara Left Mastcam-Z., Mastcam-Z que son un par de cámaras ubicadas en lo alto del mástil del rover. Esta imagen fue adquirida en mayo. 04, 2021 (Sol 72) a la hora solar media local de 14:00:34.
%% Cell type:code id: tags:
``` python
# Bibliotecas
import numpy as np
from PIL import Image
from matplotlib import pyplot as plt
import matplotlib.cm as cm
import cv2
from numpy import asarray
from matplotlib.pyplot import figure
```
%% Cell type:code id: tags:
``` python
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
#-------------------------------------------------------------------------------------------------------------
image = Image.open('Data_Mars/Mars8.png') # cargamos la imagen
matriz = asarray(image) # abrimos la imagen en formato numpy.ndarray
img = cv2.imread('Data_Mars/Mars8.png')
img = img[:,50:-50] # recortamos la imagen ya que existen bordes negros en la imagen
#-------------------------------------------------------------------------------------------------------------
print('La matriz es de la clase: ', type(matriz))
print('El número de elementos de la matriz es: ', matriz.size)
print('La dimensión de la matriz es de: ',matriz.shape)
```
%% Cell type:markdown id: tags:
- Visualización de la imagen original descargada de la página web de la NASA, y una representación de cada uno de los canales R,G,B de la imagen.
%% Cell type:code id: tags:
``` python
R, G, B = matriz[:, :, 0], matriz[:, :, 1], matriz[:, :, 2] # Separamos las matrices según R,G,B
zeros = np.zeros((B.shape), dtype=int)
# Procedemos a separar la imagen según los canales
# Canal Rojo --------------------------------------------------------------------------
matriz_canal_r = np.zeros((matriz.shape), dtype=int)
matriz_canal_r[:, :, 0], matriz_canal_r[:, :, 1], matriz_canal_r[:, :, 2] = R , zeros, zeros
# Canal Verde--------------------------------------------------------------------------
matriz_canal_g = np.zeros((matriz.shape), dtype=int)
matriz_canal_g[:, :, 0], matriz_canal_g[:, :, 1], matriz_canal_g[:, :, 2] = zeros, G, zeros
# Canal Azul---------------------------------------------------------------------------
matriz_canal_b = np.zeros((matriz.shape), dtype=int)
matriz_canal_b[:, :, 0], matriz_canal_b[:, :, 1], matriz_canal_b[:, :, 2] = zeros, zeros, B
# Gráficas-----------------------------------------------------------------------------
fig, ([ax1, ax2],[ ax3,ax4]) = plt.subplots(2,2, figsize=(12,10))
ax1.imshow(matriz)
ax1.set_title('Imagen Original')
ax2.imshow(matriz_canal_r)
ax2.set_title('Canal Rojo')
ax3.imshow(matriz_canal_g)
ax3.set_title('Canal Verde')
ax4.imshow(matriz_canal_b)
ax4.set_title('Canal Azul')
plt.show()
```
%% Cell type:markdown id: tags:
### Tranformación de la imagen de RGB a Escala de Grises
- Procedemos ha realizar la transformación de GRB a escala de Grises mediante la siguiente ecuación:
$$ Grayscale= (0.299\times R+0.587\times G+0.114\times B)$$
%% Cell type:code id: tags:
``` python
Grayscale =0.299*R+0.587*G+0.114*B
Grayscale = Grayscale[:,50:-50]
plt.imshow(Grayscale, 'gray') # imprimimos la imagen en escala de Grises.
plt.show()
```
%% Cell type:markdown id: tags:
---
<font color='gray'>
### Segmentación
La segmentación de imágenes divide la imagen en sus partes o regiones, esta técnica permite extraer información de los objetos. La división en partes se basa a menudo en las características de los píxeles de la imagen.
### Algoritmo Watershed
Es una técnica de segmentación de imágenes basada en marcadores, el algoritmo de **Watershed** se fundamenta en el concepto de visualizar una imagen como una superficie topográfica donde los valores de alta intensidad denotan picos y colinas mientras que los de baja intensidad denotan valles (mínimos locales). Esto se puede obtener trazando las coordenadas de la imagen (x, y) en función de sus valores de intensidad.
El algoritmo comienza "inundando de agua" alrededor de la superficie topográfica de la imagen. Cuando el agua fusiona los picos, se construyen barreras para evitar esta fusión, hasta que todos los picos estén bajo el agua. Las barreras resultantes dan la segmentación de picos y valles en la imagen. Este enfoque produce una sobresegmentación debido al ruido o cualquier otra irregularidad en la imagen, por lo que algunos algoritmos han implementado marcadores que especifican cuáles son todos los puntos de valle que se fusionarán y cuáles no.
<img src="imagenes/watershed.png">
</font>
%% Cell type:markdown id: tags:
<font color='gray'>
### Función Watershed de Open CV
Para realizar la segmentación se utilizará la función **cv2.watershed** de la biblioteca libre de visión artificial OpenCV (**cv2**), esta función implementa marcadores para evitar una segmentación excesiva.
Para establecer los marcadores tenemos que etiquetar la imagen en 3 regiones, la primera es la región que estamos seguros de ser fondo, la región que estamos seguros que es un objeto y finalmente la región que desconocemos que tipo de superficie sea, esta última región se la etiqueta con 0. El algoritmo retornará una matriz cuya región que sea etiquetada con -1 es la que corresponde a los límites del objeto.
---
</font>
%% Cell type:markdown id: tags:
### Establecer el umbral para la segmentación
Se establece un umbral en la imagen, esto se realiza para posteriormente determinar los marcadores, el umbral es el valor de píxel con el cual se realizará una división de la imagen en dos zonas, la primera corresponde a los píxeles menores a este umbral y la segunda los píxeles con valores superiores al umbral. Para determinar el umbral es necesario:
- Recortar una sección de la imagen donde se encuentre el objeto que quiere que el algoritmo segmente (divida).
- Realizar una cuantificación de los píxeles presentes en la imagen recortada, para poder elegir un buen umbral, esto se lo puede hacer mediante un histograma.
%% Cell type:code id: tags:
``` python
fig, (ax2, ax1) = plt.subplots(1,2, figsize=(12,4))
ax1.set_title('Histograma de la intensidad de los píxeles')
ax1.set_ylabel('Número de Píxeles')
ax1.set_xlabel('Intesidad del píxel')
ax1.hist(np.concatenate((Grayscale), axis=0), bins=250)
ax2.set_title('Imagen en escala de grises')
#ax2.xlabel('Intesidad del píxel')
ax2.imshow(Grayscale, 'gray')
plt.show()
```
%% Cell type:markdown id: tags:
- El histograma de las intensidades de luminosidad de los píxeles de la imagen nos da una idea de cual es la cantidad de objetos presentes en la imagen y su distribución en la misma.
- Los valores que se encuentran entre los valores [100,250] podrían ser considerados como fondo ya que aquí se encuentra concentrada la mayor cantidad de valores de luminosidad, mientras que los valores que son menores o mayores a este valor podrían ser los objetos.
Ahora concentrémonos en un objeto en específico que sería una roca, para eso recortamos la imagen y la analizamos, para determinar cuál es el umbral.
%% Cell type:code id: tags:
``` python
plt.imshow(Grayscale[200:400,200:600], 'gray')
plt.show()
```
%% Cell type:code id: tags:
``` python
gray_total = np.concatenate((Grayscale[200:400,200:600]), axis=0)
plt.title('Histograma de la intensidad de los píxeles en imagen recortada')
plt.ylabel('Número de Píxeles')
plt.xlabel('Intesidad del píxel')
plt.hist(gray_total, bins=250)
plt.show()
```
%% Cell type:markdown id: tags:
- El histograma nos muestra para la imagen recortada que los píxeles se dividen en 2 zonas, los píxeles con valores mayores a 125 y menores a este mismo número, por lo que tomar este valor como umbral, sería un buen comienzo para esta imagen, entonces el **umbral = 100**, este umbral nos da una perspectiva de como se encuentran distribuidas las luminosidades en la imagen y cuales píxeles se podrían considerar como fondo, y cuales no, además de lo que ya consideramos anteriormente en qué valores se encontraba el fondo, lo ideal sería elegir a los píxeles correspondientes a rocas como objeto y la arena como fondo, pero en este caso no se puede tomar de esta forma ya que las rocas tienen casi el mismo valor de píxel y algunas zonas corespondientes a fondo podría el algoritmo considerar como objeto.
%% Cell type:markdown id: tags:
### Cálculo de Marcadores
- Utilizaremos la función **cv2.threshold**, dependiendo del umbral ingresado la función genera una imagen de ceros y unos. Si el píxel de la imagen supera el umbral lo reescribe como 1 si el píxel no supera el umbral lo transforma en 0.
%% Cell type:code id: tags:
``` python
umbral = 110
_, thresh = cv2.threshold(Grayscale, umbral, 255,cv2.THRESH_BINARY_INV )
plt.imshow(thresh, 'gray')
plt.show()
```
%% Cell type:markdown id: tags:
- Ahora utilizaremos la función **Morphological Transformations**, ya que en la imagen anterior existen algunos agujeros o ruidos. Esta función es útil para cerrar pequeños agujeros dentro de los objetos en primer plano o pequeños puntos negros en el objeto.
%% Cell type:code id: tags:
``` python
kernel = np.ones((4,4))
opening = cv2.morphologyEx(thresh,cv2.MORPH_OPEN, kernel, iterations = 1)
closing = cv2.morphologyEx(opening, cv2.MORPH_CLOSE , kernel, iterations = 1)
plt.imshow(closing, 'gray')
plt.show
```
%% Cell type:markdown id: tags:
- Aquí utilizamos la función **cv2.dilate**, aumenta el límite del objeto al fondo, el objeto se hace más grande, de esta manera, podemos asegurarnos de que cualquier región de fondo en el resultado sea realmente un fondo, ya que se elimina la región de límite, que se encontro manualmente solamente determinando los valores de los píxeles. Este sera el primer marcador.
%% Cell type:code id: tags:
``` python
kernel = np.ones((1,1),np.uint8)
Mark_int = cv2.dilate(closing, kernel,iterations=1)
Mark_int = Mark_int.astype(np.uint8)
plt.imshow(Mark_int, 'gray')
plt.show()
```
%% Cell type:markdown id: tags:
- Ahora encontraremos el segundo marcador, para eso utilizaremos la función **distanceTransform**, en esta operación, las intensidades de nivel de gris de los puntos dentro de las regiones del objeto se cambian por las distancias entre estos y del valor 0 más cercano (límite), a continuación se muestra de una manera gráfica como opera esta función, por lo que la imagen se ve de una forma muy clara en el centro y va degradándose el color a medida que se acerca al límite.
<img src="imagenes/distance.png">
%% Cell type:code id: tags:
``` python
dist_transform = cv2.distanceTransform(Mark_int,cv2.DIST_L2, 3)
plt.imshow(dist_transform, 'gray')
```
%% Cell type:markdown id: tags:
- Para encontrar el segundo marcador, utilizamos la función threshold, por lo que todos los valores de píxeles que sobrepasen el umbral se identifican como objeto y todos los que no lo sobrepasen se consideran como zonas desconocidas (no se sabe si es objeto o fondo).
%% Cell type:code id: tags:
``` python
dist_transform = cv2.distanceTransform(Mark_int,cv2.DIST_L2, 3)
ret, Mark_fondo = cv2.threshold(dist_transform,0.1*dist_transform.max(),255,0)
plt.imshow(Mark_fondo, 'gray')