基本操作
from PIL import Image
import matplotlib.pyplot as plt
import cv2
import numpy as np
import platform
import math
from mpl_toolkits.mplot3d import axes3d
if "Ubuntu" in platform.version():
plt.rcParams['font.sans-serif'] = ['AR PL UKai CN']
elif "Windows" in platform.platform():
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
%config InlineBackend.figure_formats = ['svg','jpg']
读取图像并显示图像
img = Image.open("lena.bmp")
img
图像加噪声
arr = np.array(img)
nosie = np.random.random(arr.shape)
shape = (arr.shape[0], arr.shape[1])
w = 0.3
img_nosie = arr/255+nosie*w
img_nosie = img_nosie/(1+w)
img_nosie = np.array(img_nosie*255, np.uint8)
Image.fromarray(img_nosie)
图像插值
最近邻内插、双线性内插、双三次内插
a = 1.2 # 放大倍数
new_shape = (int(shape[0]*a), int(shape[1]*a))
names = ["最近邻内插", "双线性内插", "双三次内插"]
methods = [cv2.INTER_NEAREST, cv2.INTER_LINEAR, cv2.INTER_CUBIC]
fig = plt.figure(figsize=(20, 30))
axes = fig.subplots(3, 1)
for ax, name, method in zip(axes, names, methods):
out = cv2.resize(arr, new_shape, interpolation=method)
ax.imshow(out)
ax.set_title(name, fontsize=20)
ax.axis('off')
fig.tight_layout()
颜色空间
分别在RGB、HSI、CIE Lab颜色空间，显示各个通道；
RGB
arr_trans = arr.transpose(2, 0, 1)
color = "RGB"
for i, j in enumerate(arr_trans, 1):
plt.subplot(1, 4, i)
plt.axis('off')
plt.title(color[i-1])
plt.imshow(j, cmap="gray")
plt.subplot(1, 4, 4)
plt.axis('off')
plt.title(color)
plt.imshow(arr)
plt.tight_layout()
plt.show()
HSL
HSL = cv2.cvtColor(arr, cv2.COLOR_RGB2HLS)
HSL = HSL.transpose(2, 0, 1)
HSL[[1, 2]] = HSL[[2, 1]]
color = "HSL"
for i, j in enumerate(HSL, 1):
plt.subplot(1, 3, i)
plt.axis('off')
plt.title(color[i-1])
plt.imshow(j, cmap="gray")
plt.tight_layout()
plt.show()
ICE lab
Lab_cv = cv2.cvtColor(arr, cv2.COLOR_RGB2Lab)
Lab_cv = Lab_cv.transpose(2, 0, 1)
Lab_cv[[1, 2]] = Lab_cv[[2, 1]]
color = "Lab"
for i, j in enumerate(Lab_cv, 1):
plt.subplot(1, 3, i)
plt.axis('off')
plt.title(color[i-1])
plt.imshow(j, cmap="gray")
plt.tight_layout()
plt.show()
灰度图像
加权平均法是根据人的亮度感知系统调节出来的参数，是个广泛使用的标准化参数。
$$
I = 0.3 \times R + 0.59 \times G + 0.11 \times B
$$
grey_img = (arr_trans[0]*0.3 + arr_trans[1]*0.59 + arr_trans[2]*0.11+0.5).astype(np.uint8)
plt.axis('off')
plt.title('灰度图像')
plt.imshow(grey_img, cmap="gray")
plt.show()
图像变换
缩放、旋转、左右翻转、上下翻转
使用Python自带函数
缩放
img.resize((100, 100))
旋转
img.rotate(30, expand=True)
左右翻转
img.transpose(Image.FLIP_LEFT_RIGHT)
上下翻转
img.transpose(Image.FLIP_TOP_BOTTOM)
构造变换矩阵
缩放
$$
\begin{bmatrix}
c_x &0 &0 \\
0 &c_y &0 \\
0 &0 &1
\end{bmatrix}
$$
k = 0.5
m = np.array([[k, 0, 0], [0, k, 0]])
Image.fromarray(cv2.warpAffine(arr, m, (int(shape[0]*k), int(shape[1]*k))))
旋转
$$
\begin{bmatrix}
\cos\theta &-\sin\theta &0 \\
\sin\theta &\cos\theta &0 \\
0 &0 &1
\end{bmatrix}
$$
a = 30/180*math.pi
m = np.array([[np.cos(a), -np.sin(a), shape[0]*0.5], [np.sin(a), np.cos(a), 0]])
Image.fromarray(cv2.warpAffine(arr, m, (int(shape[0]*1.4), int(shape[1]*1.4))))
左右翻转
放射变换是不足以完成左右翻转和上下翻转，所以需要对图片做透视变换
$$
\begin{bmatrix}x',y',w'\end{bmatrix}=
\begin{bmatrix}u,v,w\end{bmatrix}
\begin{bmatrix}
a_{11} &a_{12} &a_{13} \\
a_{21} &a_{22} &a_{23} \\
a_{31} &a_{32} &1
\end{bmatrix}
$$
$$
x=\frac{x'}{w'}=\frac{a_{11}u+a_{21}v+a_{31}}{a_{13}u+a_{23}v+1}\\
y=\frac{y'}{w'}=\frac{a_{12}u+a_{22}v+a_{32}}{a_{13}u+a_{23}v+1}
$$
8个未知数，两条方程，所以我们需要4个点才能确定透视变换矩阵中所有的参数，所以我们选择图片的四个角的点作为参考点
r, c = shape
point1 = [[0, 0], [c, 0], [0, r], [c, r]]
point2 = [[c, 0], [0, 0], [c, r], [0, r]]
m = cv2.getPerspectiveTransform(np.float32(point1), np.float32(point2))
Image.fromarray(cv2.warpPerspective(arr, m, (r, c), borderMode=cv2.BORDER_REFLECT))
上下翻转
r, c = shape
point1 = [[0, 0], [c, 0], [0, r], [c, r]]
point2 = [[0, r], [c, r], [0, 0], [c, 0]]
m = cv2.getPerspectiveTransform(np.float32(point1), np.float32(point2))
Image.fromarray(cv2.warpPerspective(arr, m, (c, r), borderMode=cv2.BORDER_REFLECT))
直方图均衡化
def curve(x):
if type(x) == np.ndarray:
x = Image.fromarray(x)
hist = x.histogram()
curve = [hist[0]]
for i in range(1, 256):
curve.append(curve[i-1]+hist[i])
curve = np.array(curve)
return curve/curve[255]
color = "rgb"
fig = plt.figure(figsize=(10, 4))
ax = fig.add_subplot(1, 2, 1, projection='3d')
for i in range(3):
ax.bar(np.arange(256), Image.fromarray(arr_trans[i]).histogram(), zs=i, zdir='y', color=color[i], alpha=0.8)
ax.set_title("直方图")
ax.view_init(elev=10, azim=-75)
ax.set_xlim(0, 255)
ax.set_yticks([2, 1, 0])
ax.set_yticklabels(['B', 'G', 'R'])
ax = fig.add_subplot(1, 2, 2)
arr_trans_curve = [curve(i) for i in arr_trans]
for i in range(3):
ax.plot(arr_trans_curve[i], color=color[i])
ax.set_title("曲线")
ax.set_xlim(0, 255)
ax.set_ylim(0, 1.0+1e-2)
plt.show()
自己编程实现
基于直方图均衡化原理，自己编程实现
out = []
for i in range(3):
channel = np.array([arr_trans_curve[i][j]*255+0.5 for j in arr_trans[i].flatten()],np.uint8).reshape(arr_trans[i].shape)
out.append(channel)
out = np.array(out, np.uint8).transpose(1, 2, 0)
out_trans = out.transpose(2, 0, 1)
fig = plt.figure(figsize=(15, 4))
fig.add_artist
ax = fig.add_subplot(1, 3, 1)
ax.axis('off')
ax.imshow(out)
ax = fig.add_subplot(1, 3, 2, projection='3d')
for i in range(3):
ax.bar(np.arange(256), Image.fromarray(out_trans[i]).histogram(), zs=i, zdir='y', color=color[i], alpha=0.8)
ax.set_title("直方图")
ax.view_init(elev=10, azim=-60)
ax.set_xlim(0, 255)
ax.set_yticks([2, 1, 0])
ax.set_yticklabels(['B', 'G', 'R'])
out_curve = [curve(i) for i in out_trans]
ax = fig.add_subplot(1, 3, 3)
arr_trans_curve = [curve(i) for i in arr_trans]
for i in range(3):
ax.plot(out_curve[i], color=color[i])
ax.set_title("曲线")
ax.set_xlim(0, 255)
ax.set_ylim(0, 1.0+1e-2)
plt.show()
调用Python自带函数
调用python/Matlab自带直方图均衡化
out = np.array([cv2.equalizeHist(i) for i in arr_trans]).transpose(1, 2, 0)
out_trans = out.transpose(2, 0, 1)
fig = plt.figure(figsize=(15, 4))
ax = fig.add_subplot(1, 3, 1)
ax.axis('off')
ax.imshow(out)
ax = fig.add_subplot(1, 3, 2, projection='3d')
for i in range(3):
ax.bar(np.arange(256), Image.fromarray(out_trans[i]).histogram(), zs=i, zdir='y', color=color[i], alpha=0.8)
ax.set_title("直方图")
ax.view_init(elev=10, azim=-60)
ax.set_xlim(0, 255)
ax.set_yticks([2, 1, 0])
ax.set_yticklabels(['B', 'G', 'R'])
out_curve = [curve(i) for i in out_trans]
ax = fig.add_subplot(1, 3, 3)
arr_trans_curve = [curve(i) for i in arr_trans]
for i in range(3):
ax.plot(out_curve[i], color=color[i])
ax.set_title("曲线")
ax.set_xlim(0, 255)
ax.set_ylim(0, 1.0+1e-2)
plt.show()
对比
将两种方式的结果进行对比（给出均衡化前后的直方图）
fig = plt.figure(figsize=(10, 4))
ax = fig.add_subplot(1, 2, 1, projection='3d')
for i in range(3):
ax.bar(np.arange(256), Image.fromarray(arr_trans[i]).histogram(), zs=i, zdir='y', color=color[i], alpha=0.8)
ax.set_title("均衡化前直方图")
ax.view_init(elev=10, azim=-60)
ax.set_xlim(0, 255)
ax.set_yticks([2, 1, 0])
ax.set_yticklabels(['B', 'G', 'R'])
ax = fig.add_subplot(1, 2, 2, projection='3d')
for i in range(3):
ax.bar(np.arange(256), Image.fromarray(out_trans[i]).histogram(), zs=i, zdir='y', color=color[i], alpha=0.8)
ax.set_title("均衡化后直方图")
ax.view_init(elev=10, azim=-60)
ax.set_xlim(0, 255)
ax.set_yticks([2, 1, 0])
ax.set_yticklabels(['B', 'G', 'R'])
plt.show()
图像傅里叶变换
对图像进行傅里叶变换，并显示频谱和相位谱
def fft2_show(x):
plt.figure(figsize=(8, 8))
plt.subplot(2, 2, 1)
plt.imshow(x, cmap="gray")
plt.axis('off')
plt.title("原始图像")
plt.subplot(2, 2, 2)
f = np.fft.fft2(x)
f_img = np.log(np.abs(f)+1e-9)
f_img = (f_img-f_img.min())/(f_img.max()-f_img.min())
plt.title("相位谱")
plt.axis('off')
plt.imshow(f_img, cmap="gray")
plt.subplot(2, 2, 3)
fshift = np.fft.fftshift(f)
fshift_img = np.log(np.abs(fshift)+1e-9)
fshift_img = (fshift_img-fshift_img.min())/(fshift_img.max()-fshift_img.min())
plt.title("频谱")
plt.axis('off')
plt.imshow(fshift_img, cmap="gray")
plt.tight_layout()
plt.show()
fft2_show(grey_img)
进行傅里叶逆变换，计算恢复后图像与原始图像之间的L2距离（即MSE）
inv_f = np.abs(np.fft.ifft2(np.fft.fft2(grey_img)))
mse = np.power(inv_f-grey_img, 2).mean()
print(f"MSE: {mse:.4f}")
MSE: 0.0000
傅里叶变换定理
- 分别对图像做平移、旋转、尺度变换，计算变换后图像的傅里叶变换
- 显示其频谱和相位谱，并与原始图像的结果进行比较。
平移
m = np.array([[1, 0, 0.3*shape[0]], [0, 1, 0]])
out = cv2.warpAffine(grey_img, m, shape)
fft2_show(out)
旋转
a = 30/180*math.pi
m = np.array([[np.cos(a), -np.sin(a), shape[0]*0.5], [np.sin(a), np.cos(a), 0]])
out = cv2.warpAffine(grey_img, m, (int(shape[0]*1.4), int(shape[1]*1.4)))
fft2_show(out)
尺度变换
k = 1/2
m = np.array([[k, 0, shape[0]*(1-k)/2], [0, k, shape[1]*(1-k)/2]])
out = cv2.warpAffine(grey_img, m, shape)
fft2_show(out)
低通滤波
- 在傅里叶变换的基础上进行低通滤波；
- 采用理想低通、巴特沃斯低通；
- 对滤波后结果进行逆变换，得到滤波后的图像；
- 计算滤波后图像与原始图像之间的信噪比，并分析实验结果；
理想低通
d0 = 40
mat = np.zeros(shape, np.uint8)
cv2.circle(mat, (r // 2, c // 2), d0, (1, 1, 1), thickness=-1)
fig = plt.figure()
ax = fig.add_subplot(1, 2, 1, projection='3d')
X = np.arange(0, shape[1], 1)
Y = np.arange(0, shape[0], 1)
X, Y = np.meshgrid(X, Y)
ax.set_title("理想低通滤波器透视图")
ax.plot_wireframe(X, Y, mat)
ax.view_init(elev=10, azim=45)
ax.grid(False)
ax.set_zticks([])
ax.set_xticks([])
ax.set_yticks([])
ax = fig.add_subplot(1, 2, 2)
ax.imshow(mat, cmap="gray")
ax.axis('off')
ax.set_title("理想低通滤波器图像")
plt.show()
f = np.fft.fft2(grey_img)
fshift = np.fft.fftshift(f)
fshift_new = fshift*mat
fshift_new = np.fft.ifftshift(fshift_new) # 对新的进行逆变换
img_new = np.fft.ifft2(fshift_new)
img_new = np.abs(img_new)
img_new = (img_new-img_new.min())/(img_new.max()-img_new.min())
img_new = (img_new*255+0.5).astype(np.uint8)
Image.fromarray(img_new)
计算信噪比
$$
PSNR = 10\cdot \log_{10}{(\frac{MAX^2_I}{MSE})}=20 \cdot \log_{10}{(\frac{MAX_I}{\sqrt{MSE}})}
$$
mse = np.power((img_new-grey_img)/255, 2).mean()
psnr = 10*np.log10(1/mse)
print(f"信噪比为：{psnr:.2f}")
信噪比为：1.13
巴特沃斯低通
def fft_distances(m, n):
u = np.array([i if i <= m / 2 else m - i for i in range(m)], dtype=np.float32)
v = np.array([i if i <= n / 2 else n - i for i in range(n)], dtype=np.float32)
u.shape = m, 1
ret = np.sqrt(u * u + v * v)
return np.fft.fftshift(ret)
d0 = 50
n = 3
duv = fft_distances(*fshift.shape)
filter_mat = 1 / (1 + np.power(duv / d0, 2 * n))
fig = plt.figure()
ax = fig.add_subplot(1, 2, 1, projection='3d')
X = np.arange(0, shape[1], 1)
Y = np.arange(0, shape[0], 1)
X, Y = np.meshgrid(X, Y)
ax.set_title("巴特沃斯低通滤波器透视图")
ax.plot_wireframe(X, Y, filter_mat)
ax.view_init(elev=10, azim=45)
ax.grid(False)
ax.set_zticks([])
ax.set_xticks([])
ax.set_yticks([])
ax = fig.add_subplot(1, 2, 2)
ax.imshow(filter_mat, cmap="gray")
ax.axis('off')
ax.set_title("巴特沃斯低通滤波器图像")
plt.show()
fshift_new = filter_mat*fshift
fshift_new = np.fft.ifftshift(fshift_new) # 对新的进行逆变换
img_new = np.fft.ifft2(fshift_new)
img_new = np.abs(img_new)
img_new = (img_new-img_new.min())/(img_new.max()-img_new.min())
img_new = (img_new*255+0.5).astype(np.uint8)
Image.fromarray(img_new)
计算信噪比
mse = np.power((img_new-grey_img)/255, 2).mean()
psnr = 10*np.log10(1/mse)
print(f"信噪比为：{psnr:.2f}")
信噪比为：2.42
高通滤波
- 在傅里叶变换的基础上进行高通滤波；
- 分别采用理想高通、巴特沃斯高通；
- 对滤波后结果进行逆变换，得到滤波后的图像；
- 计算滤波后图像与原始图像之间的信噪比，并分析实验结果；
理想高通
d0 = 10
mat = np.ones(shape, np.uint8)
cv2.circle(mat, (r // 2, c // 2), d0, (0, 0, 0), thickness=-1)
fig = plt.figure()
ax = fig.add_subplot(1, 2, 1, projection='3d')
X = np.arange(0, shape[1], 1)
Y = np.arange(0, shape[0], 1)
X, Y = np.meshgrid(X, Y)
ax.set_title("理想高通滤波器透视图")
ax.plot_wireframe(X, Y, mat)
ax.view_init(elev=10, azim=45)
ax.grid(False)
ax.set_zticks([])
ax.set_xticks([])
ax.set_yticks([])
ax = fig.add_subplot(1, 2, 2)
ax.imshow(mat, cmap="gray")
ax.set_xticks([])
ax.set_yticks([])
ax.set_title("理想高通滤波器图像")
plt.show()
f = np.fft.fft2(grey_img)
fshift = np.fft.fftshift(f)
fshift_new = fshift*mat
fshift_new = np.fft.ifftshift(fshift_new) # 对新的进行逆变换
img_new = np.fft.ifft2(fshift_new)
img_new = np.abs(img_new)
img_new = (img_new-img_new.min())/(img_new.max()-img_new.min())
img_new = (img_new*255+0.5).astype(np.uint8)
Image.fromarray(img_new)
巴特沃斯高通
d0 = 4
n = 2
duv = fft_distances(*fshift.shape)
duv[r // 2, c // 2] = 1e-8
filter_mat = 1 / (1 + np.power(d0 / duv, 2 * n))
fig = plt.figure()
ax = fig.add_subplot(1, 2, 1, projection='3d')
X = np.arange(0, shape[1], 1)
Y = np.arange(0, shape[0], 1)
X, Y = np.meshgrid(X, Y)
ax.set_title("巴特沃斯高通滤波器透视图")
ax.plot_wireframe(X, Y, filter_mat)
ax.view_init(elev=10, azim=45)
ax.grid(False)
ax.set_zticks([])
ax.set_xticks([])
ax.set_yticks([])
ax = fig.add_subplot(1, 2, 2)
ax.imshow(filter_mat, cmap="gray")
ax.set_xticks([])
ax.set_yticks([])
ax.set_title("巴特沃斯高通滤波器图像")
plt.show()
fshift_new = filter_mat*fshift
fshift_new = np.fft.ifftshift(fshift_new) # 对新的进行逆变换
img_new = np.fft.ifft2(fshift_new)
img_new = np.abs(img_new)
img_new = (img_new-img_new.min())/(img_new.max()-img_new.min())
img_new = (img_new*255+0.5).astype(np.uint8)
Image.fromarray(img_new)
