Convolución 2D¶

Métodos potenciales de prospección, FCAG, 2023.

import numpy as np
from numpy.fft import fft,ifft
import matplotlib.pyplot as plt
plt.xkcd()
plt.rcParams.update({'font.size': 14}) # tamaño de fuente

Convolución lineal 2D $\ast$¶

Para calcular la convolución lineal 2D entre $\mathbf{A}$ y $\mathbf{B}$ procedemos de la siguiente manera. Hacemos zero-padding para embebeber los operandos en matrices de dimensiones iguales a la dimensión de la convolución lineal 2D. Operamos en frecuencia, volvemos al dominio original y finalmente removemos efectos de borde:

from scipy.signal import convolve2d # para utilizar la subrutina de la convolución 2D

np.random.seed(2020) # para poder reproducir los mismos resultados

A = np.random.random((5,5))
B = np.random.random((10,10))

Nx,Ny = B.shape
Mx,My = A.shape

Lx    = Nx+Mx-1
Ly    = Ny+My-1

# Relleno con ceros para obtener la dimensión deseada:

A0            = np.zeros((Lx,Ly))
A0[:Mx,:My]   = A

B0            = np.zeros((Lx,Ly))
B0[:Nx,:Ny]   = B

conv_fft = np.real ( np.fft.ifft2( np.fft.fft2(A0)*np.fft.fft2(B0) ) )

# Otra forma: haciendo el zero-padding dentro de la FFT
#
#conv_fft = np.real ( np.fft.ifft2( np.fft.fft2(A,s=(Lx,Ly))*np.fft.fft2(B,s=(Lx,Ly)) ) )

# Fig.

cmap="inferno"

av = A0 #.reshape(1,N+N-1) # para graficar
bv = B0 #.reshape(1,N+M-1) # para graficar
cv = conv_fft

vmin,vmax   = np.min(av), np.max(av)
vmin2,vmax2 = np.min(cv), np.max(cv)

extent = [0,cv.shape[0],0,cv.shape[1]]


plt.figure(figsize=(15,5))
plt.subplot(1,3,1)
plt.title("$\mathbf{A_{pad}}$")
plt.imshow(av,interpolation="None",cmap=cmap,vmin=vmin,vmax=vmax,extent=extent);
plt.grid(color="black",linestyle='-',linewidth=1.5)
plt.xticks([]);plt.yticks([])
plt.subplot(1,3,2)
plt.title("$\mathbf{B_{pad}}$")
plt.imshow(bv,interpolation="None",cmap=cmap,vmin=vmin,vmax=vmax,extent=extent);
plt.xticks([]);plt.yticks([])
plt.subplot(1,3,3)
plt.title("$\mathbf{A*B}$")
plt.imshow(cv,interpolation="None",cmap=cmap,vmin=vmin2,vmax=vmax2,extent=extent);
plt.xticks([]);plt.yticks([])
plt.show()

Verificamos el resultado de la convolución lineal por medio de scipy.signal.convolve2d:

from scipy.signal import convolve2d
conv_full = convolve2d(image, kernel, mode="full",boundary="fill", fillvalue=0)
...

#conv_full     = convolve2d(A, B, mode='full',boundary='fill', fillvalue=0, )
conv_full     = convolve2d(A, B, mode="full")

print(np.allclose(conv_full,conv_fft))

True

Removemos efectos de borde del resultado considerando que la imagen a procesar es $\mathbf{B}$. Es decir, la salida tiene la dimensión de $\mathbf{B}$ y suponemos que la muestra central de $\mathbf{A}$ es la muestra en el orígen:

# Removiendo efectos de borde:

offset_x   = Mx//2  # válido cuando Mx es impar
offset_y   = My//2  # válido cuando My es impar


conv_fft_same = conv_fft[offset_x:-offset_x,offset_y:-offset_y]

# Fig.
plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
plt.title("$\mathbf{A*B}$")
plt.imshow(cv,interpolation="None",cmap=cmap,vmin=vmin2,vmax=vmax2,extent=extent);
plt.xticks([]);plt.yticks([])

plt.subplot(1,2,2)
plt.title("$\mathbf{A*B}$ sin bordes")
plt.imshow(conv_fft_same,interpolation="None",cmap=cmap,vmin=vmin2,vmax=vmax2,extent=extent);
plt.xticks([]);plt.yticks([])

# plt.subplot(1,3,3)
# plt.title("$\mathbf{a*b}$ same ")
# plt.imshow(conv_same,interpolation="None",cmap=cmap,vmin=vmin2,vmax=vmax2,extent=extent);
# plt.xticks([]);plt.yticks([])
plt.show()

Verificamos el resultado de la convolución lineal sin efectos de borde por medio de scipy.signal.convolve2d:

from scipy.signal import convolve2d
conv_same = convolve2d(image, kernel, mode="same",boundary="fill", fillvalue=0)
...

conv_same     = convolve2d(B, A, mode="same") # Notar que de usar "same", entonces el tamaño de salida es de la primer entrada
print(np.allclose(conv_same,conv_fft_same))

True

Convolución circular 2D $\circledast$¶

np.random.seed(2020) # para poder reproducir los mismos resultados

A = np.random.random((3,3))
B = np.random.random((5,5))

Nx,Ny = B.shape
Mx,My = A.shape

AC = np.zeros_like(B)
# Centramos el operador para el ejemplo considerado
AC[0,0]   = A[1,1]
AC[0,1]   = A[1,2]
AC[0,-1]  = A[1,0]
AC[1,0]   = A[2,1]
AC[1,1]   = A[2,2]
AC[1,-1]  = A[2,0]
AC[-1,0]  = A[0,1]
AC[-1,1]  = A[0,2]
AC[-1,-1] = A[0,0]

convc_fft = np.real ( np.fft.ifft2( np.fft.fft2(AC)*np.fft.fft2(B) ) )

# Fig.

cmap="inferno"

av = AC 
bv = B 
cv = convc_fft

vmin,vmax   = np.min(av), np.max(av)
vmin2,vmax2 = np.min(cv), np.max(cv)

extent = [0,cv.shape[0],0,cv.shape[1]]


plt.figure(figsize=(15,5))
plt.subplot(1,3,1)
plt.title("$\mathbf{A_{C}}$")
plt.imshow(av,interpolation="None",cmap=cmap,vmin=vmin,vmax=vmax,extent=extent);
plt.grid(color="black",linestyle='-',linewidth=1.5)
plt.xticks([]);plt.yticks([])
plt.subplot(1,3,2)
plt.title("$\mathbf{B}$")
plt.imshow(bv,interpolation="None",cmap=cmap,vmin=vmin,vmax=vmax,extent=extent);
plt.xticks([]);plt.yticks([])
plt.subplot(1,3,3)
plt.title("$\mathbf{A_C {\circ} B}$")
plt.imshow(cv,interpolation="None",cmap=cmap,vmin=vmin2,vmax=vmax2,extent=extent);
plt.xticks([]);plt.yticks([])
plt.show()

Verificamos el resultado de la convolución circular por medio de scipy.signal.convolve2d:

from scipy.signal import convolve2d
convc = convolve2d(image, kernel, mode="same",boundary="wrap")
...

convc     = convolve2d(B, A, mode="same", boundary="wrap") # "same" espera que la imagen a procesar sea el primer elemento
print(np.allclose(convc,convc_fft))

True

Nota¶

De no hacer el ordenamiento de $\mathbf{A}$, obtendremos (por el teorema del corrimiento de fase) el resultado desplazado considerando el dominio cíclico. El resultado aparece corrido $i_0$ muestras en una dimensión y $j_0$ en la otra dimensión, donde los valores $i_0$ y $j_0$ son las coordenadas de la muestra central. Calculamos para observar esto, tanto de forma numérica como en una figura:

A0          = np.zeros_like(B)
A0[:Mx,:My] = A
conv_fft    = np.real ( np.fft.ifft2( np.fft.fft2(A0)*np.fft.fft2(B) ) )

print(np.round(convc_fft,2))
print("")
print(np.round(conv_fft,2))

[[1.72 1.97 2.13 2.24 1.64]
 [2.44 1.88 2.12 2.8  2.51]
 [2.16 1.88 2.03 1.47 2.03]
 [2.38 2.26 2.54 1.94 2.19]
 [1.24 2.03 2.61 2.03 1.76]]

[[1.76 1.24 2.03 2.61 2.03]
 [1.64 1.72 1.97 2.13 2.24]
 [2.51 2.44 1.88 2.12 2.8 ]
 [2.03 2.16 1.88 2.03 1.47]
 [2.19 2.38 2.26 2.54 1.94]]

# Fig.

cmap="inferno"

av = convc_fft #.reshape(1,N+N-1) # para graficar
bv = conv_fft #.reshape(1,N+M-1) # para graficar

vmin=np.min(conv_fft)
vmax =np.max(conv_fft)
extent = [0,B.shape[0],0,B.shape[1]]


plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
plt.title("$\mathbf{A_C \circ B}$")
plt.imshow(av,interpolation="None",cmap=cmap,vmin=vmin,vmax=vmax,extent=extent);
plt.xticks([]);plt.yticks([])
plt.subplot(1,2,2)
plt.title("$\mathbf{A \circ B}$")
plt.imshow(bv,interpolation="None",cmap=cmap,vmin=vmin,vmax=vmax,extent=extent);
plt.xticks([]);plt.yticks([])
plt.show()

Eso es todo por hoy.