• 二维离散傅里叶变换DFT
• 二维图像卷积
• Canny边缘检测

# 二维离散傅里叶变换DFT

F ( u , v ) = ∑ x = 0 M − 1 ∑ y = 0 N − 1 f ( x , y ) e − j 2 π ( u x M + v y N ) (1.1) F(u, v) = \sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x, y)e^{-j2\pi(\frac{ux}{M}+\frac{vy}{N})} \tag{1.1}

for u in range(M):
for v in range(N):
for x in range(M):
for y in range(N):
comp = np.exp(-2j * np.pi * (u * x / M + v * y / N))


G 1 = [ e − j 2 π 0 × 0 M e − j 2 π 0 × 1 M . . . e − j 2 π 0 × ( M − 1 ) M e − j 2 π 1 × 0 M e − j 2 π 1 × 1 M . . . e − j 2 π 1 × ( M − 1 ) M . . . . . . . . . . . . e − j 2 π ( M − 1 ) × 0 M e − j 2 π ( M − 1 ) × 1 M . . . e − j 2 π ( M − 1 ) × ( M − 1 ) M ] G 2 = [ e − j 2 π 0 × 0 N e − j 2 π 0 × 1 N . . . e − j 2 π 0 × ( N − 1 ) N e − j 2 π 1 × 0 N e − j 2 π 1 × 1 N . . . e − j 2 π 1 × ( N − 1 ) N . . . . . . . . . . . . e − j 2 π ( N − 1 ) × 0 N e − j 2 π ( N − 1 ) × 1 N . . . e − j 2 π ( N − 1 ) × ( N − 1 ) N ] G_1 = \left[ \begin{matrix} e^{-j2\pi\frac{0\times0}{M}}&e^{-j2\pi\frac{0\times1}{M}}&...&e^{-j2\pi\frac{0\times(M-1)}{M}}\\e^{-j2\pi\frac{1\times0}{M}}&e^{-j2\pi\frac{1\times1}{M}}&...&e^{-j2\pi\frac{1\times(M-1)}{M}}\\ ... & ... & ... & ...\\e^{-j2\pi\frac{(M-1)\times0}{M}}&e^{-j2\pi\frac{(M-1)\times1}{M}}&...&e^{-j2\pi\frac{(M-1)\times(M-1)}{M}}\end{matrix} \right] \\ \\ G_2 = \left[ \begin{matrix} e^{-j2\pi\frac{0\times0}{N}}&e^{-j2\pi\frac{0\times1}{N}}&...&e^{-j2\pi\frac{0\times(N-1)}{N}}\\e^{-j2\pi\frac{1\times0}{N}}&e^{-j2\pi\frac{1\times1}{N}}&...&e^{-j2\pi\frac{1\times(N-1)}{N}}\\ ... & ... & ... & ...\\e^{-j2\pi\frac{(N-1)\times0}{N}}&e^{-j2\pi\frac{(N-1)\times1}{N}}&...&e^{-j2\pi\frac{(N-1)\times(N-1)}{N}}\end{matrix} \right]

F ( u , v ) = G 1 f ( x , y ) G 2 (1.2) F(u, v) = G_1f(x, y)G_2 \tag{1.2}

m, n = img.shape[:2]
img_fft = np.zeros((m, n), dtype=np.complex)    # 经过计算得到的图像的像素是复数
G1 = dftMatrix(m)
G2 = dftMatrix(n)
img_fft = G1.dot(img).dot(G2)
img_fft = np.abs(img_fft)   # 取绝对值, 将复数映射到实数
img_fft = np.log(img_fft)   # 将像素值映射到0~255


def dftMatrix(N):
i, j = np.meshgrid(np.arange(N), np.arange(N))
return np.power(np.exp(-2j * np.pi / N), i * j)


# 将图像四个象限对换
def transform(img):
(n, m) = img.shape
nmid, mmid = n // 2, m // 2
if (n & 1) and (m & 1):  # n, m都为奇数

piece1 = img[0:nmid, :][:, 0:mmid]
piece2 = img[0:nmid, :][:, mmid + 1:]
piece3 = img[nmid + 1:, :][:, mmid + 1:]
piece4 = img[nmid + 1:, :][:, 0:mmid]
block1 = img[0:nmid, mmid].reshape((mmid, 1))
block2 = img[nmid + 1:, mmid].reshape((mmid, 1))
block3 = np.array([img[nmid, :]])

part1 = np.hstack((piece3, block1))
part1 = np.hstack((part1, piece4))
part2 = np.hstack((piece2, block2))
part2 = np.hstack((part2, piece1))
part1 = np.concatenate((part1, block3), axis=0)
whole = np.concatenate((part1, part2), axis=0)
return whole

elif (n & 1) and not (m & 1):  # n为奇数, m为偶数
piece1 = img[0:nmid, :][:, 0:mmid]
piece2 = img[0:nmid, :][:, mmid:]
piece3 = img[nmid + 1:, :][:, mmid:]
piece4 = img[nmid + 1:, :][:, 0:mmid]
block1 = np.array([img[nmid, :]])

part1 = np.hstack((piece3, piece4))
part2 = np.hstack((piece2, piece1))
part1 = np.concatenate((part1, block1), axis=0)
whole = np.concatenate((part1, part2), axis=0)
return whole

elif not (n & 1) and (m & 1):  # n为偶数, m为奇数
piece1 = img[0:nmid, :][:, 0:mmid]
piece2 = img[0:nmid, :][:, mmid + 1:]
piece3 = img[nmid:, :][:, mmid + 1:]
piece4 = img[nmid:, :][:, 0:mmid]
block1 = img[:, mmid].reshape((n, 1))

part1 = np.concatenate((piece3, piece2), axis=0)
part2 = np.concatenate((piece4, piece1), axis=0)
part1 = np.hstack((part1, block1))
whole = np.hstack((part1, part2))
return whole

else:
piece1 = img[0:nmid, :][:, 0:mmid]
piece2 = img[0:nmid, :][:, mmid:]
piece3 = img[nmid:, :][:, mmid:]
piece4 = img[nmid:, :][:, 0:mmid]

part1 = np.concatenate((piece3, piece2), axis=0)
part2 = np.concatenate((piece4, piece1), axis=0)
whole = np.hstack((part1, part2))
return whole


def myFFT(img):
"""
对图像img做DFT变换, 返回傅里叶频谱
:param img:
:return:
"""
# 先将图片转换成灰度图
if 2 == len(img.shape):  # 长度为2说明ff_img已经是灰度图
pass
else:
img = cv.cvtColor(img, cv.COLOR_BGR2GRAY)
m, n = img.shape[:2]

img_fft = np.zeros((m, n), dtype=np.complex)  # 经过计算得到的图像的像素是复数
G1 = dftMatrix(m)
G2 = dftMatrix(n)
img_fft = G1.dot(img).dot(G2)
img_fft = np.abs(img_fft)  # 取绝对值, 将复数映射到实数
img_fft = np.log(img_fft)  # 将像素值映射到0~255

# 将四个象限对调
return transform(img_fft)

def testFFT(img):
if 2 == len(img.shape):  # 长度为2说明ff_img已经是灰度图
pass
else:
img = cv.cvtColor(img, cv.COLOR_BGR2GRAY)
img_orig = deepcopy(img)
img_fshift = np.fft.fftshift(np.fft.fft2(img))  # 转移像素做幅度普
img_fft_np = np.log(np.abs(img_fshift))  # 取绝对值,将复数变化成实数,取对数是为了将数据变化到0-255
img_fft_my = myFFT(deepcopy(img))
images = [img_orig, img_fft_np, img_fft_my]
titles = ['原图', 'NumPy傅里叶频谱', '自定义傅里叶频谱']
for i in range(3):
plt.subplot(1, 3, i + 1)
plt.imshow(images[i], cmap='gray')
plt.title(titles[i])
plt.xticks([]), plt.yticks([])  # 隐藏X,Y轴
plt.show()


1. 计算加速矩阵 G 1 , G 2 G_1, G_2
2. 计算原始傅里叶频谱
3. 将原始傅里叶频谱中心化

# 二维图像卷积

A [ i ] [ j ] = ∑ p r o d u c t A[i][j] = \sum product

1. 将卷积核 K e r n e l Kernel 上下翻转
3. K e r n e l Kernel 作为滑动窗口, 在被卷积矩阵 A A 上滑动, 同时计算卷积结果

def myConvol2D(convol_img, kernel, padding_values=0):
"""

:param convol_img:      被卷积的图像,必须为灰度图
:param kernel:          卷积核,必须为方阵,且长宽为奇数
:return convol_matrix:  卷积后的矩阵, shape同convol_img
"""
if 2 != len(convol_img.shape):  # 判断convol_img是否为灰度图
raise TypeError('wrong shape of parameter img:{}'.format(convol_img.shape))
elif (kernel.shape[0] != kernel.shape[1]) or (1 != kernel.shape[0] % 2):  # 判断kernel是否为奇数尺寸的方阵
raise TypeError('wrong shape of parameter kernel:{}'.format(kernel.shape))
else:
# 使用原图像的复制矩阵, 避免在处理时改变原图像
convol_matrix = np.zeros(convol_img.shape, dtype=np.uint8)
kernel = kernel[::-1]  # 上下翻转卷积核
kernel_rows_cols = kernel.shape[0]  # 卷积核的尺寸
pad_num = kernel_rows_cols // 2  # 填充深度
(convol_matrix_rows, convol_matrix_cols) = convol_matrix.shape  # 原矩阵形状

# 开始卷积

# np.multiply(a,b)返回矩阵a和矩阵b所有元素的乘积之和
temp_sum = np.sum(np.multiply(x_y_filtered_matrix, kernel))

# 防止溢出, 并将像素值数据类型转换为8位无符号整数0~255
if temp_sum < 0:
temp_sum = np.uint8(0)
elif temp_sum > 255:
temp_sum = np.uint8(255)
else:
temp_sum = np.uint8(temp_sum)
# 将卷积结果存入对应的convol_matrix[x][y]

return convol_matrix

def testConvol2D(img):
if 2 == len(img.shape):  # 长度为2说明ff_img已经是灰度图
pass
else:
img = cv.cvtColor(img, cv.COLOR_BGR2GRAY)
kernel = np.array([[-1, -1, -1, -1, -1],
[-1, 1, 2, 1, -1],
[-1, 2, 4, 2, -1],
[-1, 1, 2, 1, -1],
[-1, -1, -1, -1, -1]])
img_orig = deepcopy(img)
img_con_cv = cv.filter2D(img, -1, kernel)
img_con_my = myConvol2D(img, kernel)
images = [img_orig, img_con_my, img_con_cv]
titles = ['原图', '自定义卷积', 'OpenCV卷积']

for i in range(3):
plt.subplot(1, 3, i + 1)
plt.imshow(images[i], cmap='gray')
plt.title(titles[i])
plt.xticks([]), plt.yticks([])  # 隐藏X,Y轴
plt.show()


# Canny边缘检测

canny边缘检测的算法过于复杂, 这里我引用另一博主的文章

def myCanny(img):
# 先将图片转换成灰度图
if 2 == len(img.shape):  # 长度为2说明ff_img已经是灰度图
pass
else:
img = cv.cvtColor(img, cv.COLOR_BGR2GRAY)
sigma1 = sigma2 = 1
sum = 0
gaussian = np.zeros([5, 5])
for i in range(5):
for j in range(5):
gaussian[i, j] = math.exp(-1 / 2 * (np.square(i - 3) / np.square(sigma1)  # 生成二维高斯分布矩阵
+ (np.square(j - 3) / np.square(sigma2)))) / (
2 * math.pi * sigma1 * sigma2)
sum = sum + gaussian[i, j]

gaussian = gaussian / sum

# step1.高斯滤波
W, H = img.shape
new_gray = np.zeros([W - 5, H - 5])
for i in range(W - 5):
for j in range(H - 5):
new_gray[i, j] = np.sum(img[i:i + 5, j:j + 5] * gaussian)  # 与高斯矩阵卷积实现滤波

# step2.增强 通过求梯度幅值
W1, H1 = new_gray.shape
dx = np.zeros([W1 - 1, H1 - 1])
dy = np.zeros([W1 - 1, H1 - 1])
d = np.zeros([W1 - 1, H1 - 1])
for i in range(W1 - 1):
for j in range(H1 - 1):
dx[i, j] = new_gray[i, j + 1] - new_gray[i, j]
dy[i, j] = new_gray[i + 1, j] - new_gray[i, j]
d[i, j] = np.sqrt(np.square(dx[i, j]) + np.square(dy[i, j]))  # 图像梯度幅值作为图像强度值

# setp3.非极大值抑制 NMS
W2, H2 = d.shape
NMS = np.copy(d)
NMS[0, :] = NMS[W2 - 1, :] = NMS[:, 0] = NMS[:, H2 - 1] = 0
for i in range(1, W2 - 1):
for j in range(1, H2 - 1):
if d[i, j] == 0:
NMS[i, j] = 0
else:
# 如果Y方向幅度值较大
grad2 = d[i - 1, j]
grad4 = d[i + 1, j]
# 如果x,y方向梯度符号相同
grad1 = d[i - 1, j - 1]
grad3 = d[i + 1, j + 1]
# 如果x,y方向梯度符号相反
else:
grad1 = d[i - 1, j + 1]
grad3 = d[i + 1, j - 1]
# 如果X方向幅度值较大
else:
grad2 = d[i, j - 1]
grad4 = d[i, j + 1]
# 如果x,y方向梯度符号相同
grad1 = d[i + 1, j - 1]
grad3 = d[i - 1, j + 1]
# 如果x,y方向梯度符号相反
else:
grad1 = d[i - 1, j - 1]
grad3 = d[i + 1, j + 1]
else:
NMS[i, j] = 0

# step4. 双阈值算法检测、连接边缘
W3, H3 = NMS.shape
DT = np.zeros([W3, H3])
# 定义高低阈值
TL = 0.2 * np.max(NMS)
TH = 0.3 * np.max(NMS)
for i in range(1, W3 - 1):
for j in range(1, H3 - 1):
if (NMS[i, j] < TL):
DT[i, j] = 0
elif (NMS[i, j] > TH):
DT[i, j] = 1
elif ((NMS[i - 1, j - 1:j + 1] < TH).any() or (NMS[i + 1, j - 1:j + 1]).any()
or (NMS[i, [j - 1, j + 1]] < TH).any()):
DT[i, j] = 1

return DT

def testCanny(img):
if 2 == len(img.shape):  # 长度为2说明ff_img已经是灰度图
pass
else:
img = cv.cvtColor(img, cv.COLOR_BGR2GRAY)
img_orig = deepcopy(img)
img_canny_cv = cv.Canny(deepcopy(img), 200, 300)
img_canny_my = myCanny(deepcopy(img))
images = [img_orig, img_canny_cv, img_canny_my]
titles = ['原图', 'OpenCV Canny', '自定义Canny']
for i in range(3):
plt.subplot(1, 3, i + 1)
plt.imshow(images[i], cmap='gray')
plt.title(titles[i], loc='center')
plt.xticks([]), plt.yticks([])  # 隐藏X,Y轴
plt.show()


# 本篇完整代码

06-26 1947

01-26 3428
10-05 351
07-31 5179
12-28
01-23
©️2020 CSDN 皮肤主题: Age of Ai 设计师:meimeiellie