数字图像处理

数字图像处理

2021/4/24: 解决(?)公式渲染问题,卸了重装了渲染器,然后参考配置了服务器渲染

空间&频率域滤波

用于记录学习数字图像处理的历程(不完全同步学校课程,不定期摸鱼)

空间域滤波

直方图均衡化

直方图的形状往往与图像的外观有关。

  • 在暗图像的直方图中,大多数直方图容器集中在灰度级较低的一端(靠近 0)
  • 亮的图像的直方图中,多数直方图容器集中在灰度级的高端
  • 低对比度的图像的直方图中,直方图容器基本位于灰度级的中间
  • 高对比度图像的直方图,容器覆盖了较宽的范围,并且像素的分布是基本均匀的 直方图均衡化的数学原理我们暂且不提,假如我们把直方图每个部分出现频率和总统计量做除法,得到一个概率分布函数$p_r(r)$,均衡化变换就是将一个概率分布不规则,不均匀的$p_r(r)$变化为分布均匀的新的函数$p_s(s)$。

Python 代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
from PIL import Image
import matplotlib.pyplot as plt
import numpy as np

# 划分不同的空间来绘图
img = Image.open("./IMG_2546.JPG").convert("L") # 转为灰度图
img_array = np.array(img)
x_num, y_num = img.size
# 查看直方图
plt.subplot(121)
plt.xticks([])
plt.yticks([])
plt.hist(img_array.flatten(), 256) # 变换成 1 维再绘制直方图


# 使用直方图均衡再查看直方图
# 每个 bins 数组的区间值对应一个 im_hist 数组中的强度值
im_hist, bins = np.histogram(img_array.flatten(), 256, density=True) # 转化为直方图
cdf = im_hist.cumsum() # 计算累计分布函数 cumulative distribution function
# cdf(a)=P(x<=a)
cdf = 255 * (cdf / cdf[-1]) # 归一化到 0-1 之后转化到 0-255
# interp(x, xp, yp)输入原函数的一系列点(xp, yp)使用线性插值方法模拟函数并计算 f(x)
# s = cdf(r)
img_array_transformed = np.interp(img_array.flatten(), bins[:256], cdf)
# 变换回原来的数组结构后转为图片
plt.subplot(122)
plt.hist(img_array_transformed, 256)
plt.xticks([])
plt.yticks([])
plt.savefig("exp.png")
plt.show()

原直方图和均衡化后的图片的直方图:

直方图均衡化.png

原图和均衡化后的图的对比:

原图.png

均衡化后.png

如果我们对局部直方图均衡化,还可以显示出全局直方图均衡化无法显示的灰度细节。

低通空间滤波器

  • 高斯核

3x3高斯核

有现成的库可以调用:

1
2
from PIL import Image, ImageFilter
img = img.filter(ImageFilter.GaussianBlur(radius=2))

高通空间滤波器

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
from scipy import signal

def sharpening_f(img_array, kernel, c):
# 卷积锐化图像

new_array = signal.convolve2d(img_array, kernel,
boundary='symm',
mode='same')
new_array[new_array > 255] = 255
new_array[new_array < 0] = 0

new_array = img_array + c*new_array
new_array[new_array > 255] = 255
new_array[new_array < 0] = 0
return new_array
  • Sobel

横方向的 Sobel 算子,用于提取出水平方向的边界:

  • 拉普拉斯

拉普拉斯是导数算子,因此会突出图像中急剧的灰度过渡,并不强调缓慢变化的灰度区域,这会使得原图像产生灰色边缘线和其他不连续的特征,因此将原图像与拉普拉斯变换后的图像相加就能够恢复背景特征,并且保留拉普拉斯锐化的效果。

一种拉普拉斯核如下:

又或者是:

  • 钝化掩蔽

分为三个步骤:

  1. 模糊原图像
  2. 从原图中减去模糊后的图像(产生的差称为模板)
  3. 将模板与原图像相加

令 $\bar{f}(x, y)$ 表示模糊后的图像,钝化掩蔽的过程可以用公式表示为:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
from PIL import Image, ImageFilter
import numpy as np
import matplotlib.pyplot as plt

image = Image.open("./IMG_2546.JPG").convert("L")
origin_array = np.array(image)
blurred_image = image.filter(ImageFilter.BLUR)
blurred_array = np.array(blurred_image)
mask_array = origin_array - blurred_array

weight_k = 2
new_array = origin_array + weight_k * mask_array
fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.hist(origin_array.flatten(), 256)
ax1.set_title("Origin")
ax2.hist(new_array.flatten(), 256)
ax2.set_title("shielding Sharpened")
plt.savefig("exp6_sharping2.png")
plt.show()

原图像:
原图

增强后的图像:
钝化掩蔽后的图像

可以看到,我们仅仅是把 k 设置为 2,就出现了边界线。当 $k = 1$ 称为钝化掩蔽,$k > 1$ 时被称为高提升滤波,选择 $k < 1$ 可以减少钝化模板的贡献。

频率域滤波

在深度学习出来前,频率域滤波一直是数字图像处理的比较热门的研究点(我们老师说的),频率域滤波的功能还是挺强大的,同样关于原理不过多介绍。

FFT&频谱

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
import matplotlib.pyplot as plt
import numpy as np
import os
from PIL import Image


def fft_transform(in_path, name, out_path="."):
image = Image.open(in_path)
plt.subplot(221)
plt.imshow(image)
# 关闭坐标轴
plt.axis("off")
plt.title("Origin Image")
image = image.convert("L")
plt.subplot(222)
plt.imshow(image, "gray")
plt.axis("off")
plt.title("Gray Image")
# 进行傅立叶变换,并显示结果
fft2 = np.fft.fft2(image)
# 取 log
log_fft2 = np.log(1 + np.abs(fft2))
# 将图像变换的原点移动到频域矩形的中心,并显示效果
shift2center = np.fft.fftshift(fft2)
# 对中心化后的结果进行对数变换
log_shift2center = np.log(1 + np.abs(shift2center))
plt.subplot(223)
plt.imshow(np.absolute(shift2center), "gray")
plt.axis("off")
plt.title("shift2center")
plt.subplot(224)
plt.imshow(log_shift2center, "gray")
plt.axis("off")
plt.title("log_shift2center")
plt.savefig(os.path.join(out_path, ("fft_transformed_" + name)))
# clear all figure
plt.clf()

fft_transformed.png

需要注意的是要对变换后的图像取 log,不然值太大会无法在正常的灰度级下显示。

高低通滤波器

课本上介绍的 3 个高通滤波器:

  • 理想高通滤波器
  • 高斯高通滤波器
  • 巴德沃斯高通滤波器

其中 $D_0$ 表示截止频率到矩阵中心的距离,$D(u,v)$ 表示频率矩阵中心到矩阵中任意一点的距离。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
from PIL import Image
import numpy as np
import matplotlib.pyplot as plt


def cal_distance(pa, pb):
from math import sqrt
dis = sqrt((pa[0] - pb[0]) ** 2 + (pa[1] - pb[1]) ** 2)
return dis


class Filter:

@classmethod
def generate_filter(cls, d, shape, *args, **kwargs):
transfer_matrix = np.zeros(shape)
center_point = tuple(map(lambda x: (x - 1) // 2, shape))
for i in range(transfer_matrix.shape[0]):
for j in range(transfer_matrix.shape[1]):
dist = cal_distance(center_point, (i, j))
transfer_matrix[i, j] = cls.get_one(d, dist, *args, **kwargs)
return transfer_matrix

@classmethod
def get_one(cls, d, dist, *args, **kwargs) -> float:
return 1


class ILPFLowPassFilter(Filter):

@classmethod
def get_one(cls, d, dist, *args, **kwargs) -> float:
if dist <= d:
return 1
else:
return 0


class ILPFHighPassFilter(Filter):

@classmethod
def get_one(cls, d, dist, *args, **kwargs) -> float:
if dist <= d:
return 0
else:
return 1


class GaussianHighPassFilter(Filter):

@classmethod
def get_one(cls, d, dist, *args, **kwargs) -> float:
return 1 - np.exp(-(dist ** 2) / (2 * (d ** 2)))


class GaussianLowPassFilter(Filter):

@classmethod
def get_one(cls, d, dist, *args, **kwargs) -> float:
return np.exp(-(dist ** 2) / (2 * (d ** 2)))


class ButterworthFilter(Filter):

@classmethod
def get_one(cls, d, dist, *args, **kwargs) -> float:
n = kwargs["n"]
return 1 / ((1 + dist / d) ** (2 * n))

使用滤波器进行滤波:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

choice = "material/1.jpg"
img = Image.open(choice).convert("L")

sq = min(img.size[0], img.size[1])
img = img.resize((sq, sq))
f = np.fft.fft2(img)
f_shift = np.fft.fftshift(f)
# gauss_filter_matrix = GaussianLowPassFilter.generate_filter(50, img.size)
butter_filter_matrix = ButterworthFilter.generate_filter(30, img.size, n=2)
filter_matrix = np.fft.ifftshift(f_shift*butter_filter_matrix)
transform_img = np.abs(np.fft.ifft2(filter_matrix))
plt.axis("off")
plt.imshow(transform_img, cmap="gray")
plt.savefig("gauss_filter.png")
plt.show()

扩展:人脸磨皮算法

Bilateral

双边滤波(Bilateral filter)结合图像的空间邻近度和像素值相似度,同时考虑空域信息和灰度相似性,达到保边去噪的目的。

它的滤波器核由两个函数生成:空间域核和值域核。

空间域核是由像素位置欧式距离决定的模板权值,公式:

其中i,j 代表的是当前坐标点的位置 k,l 为中心坐标点,$\delta_d$ 代表高斯函数的标准差。很明显 $w_d$ 是计算临近点 ij 到中心点的临近程度,因此空间域核是用于衡量空间临近的程度。这代表空间域的高斯函数。

值域核是由灰度像素值的差值决定模板的权值的:

$f(i,j)$代表每个点的灰度像素值,$f(k,l)$代表中点的像素值,$\delta_r$也是值域核下高斯函数的标准差。将两者相乘就能得到双边滤波的模板权值:

化简:

1
2
3
4
5
6
7
8
9
10
11
# 公式看起来实现起来好麻烦,还好有库
import cv2
file_name = './material/4.jpg'
image = cv2.imread(file_name)

dst = cv2.bilateralFilter(src=image, d=0, sigmaColor=100, sigmaSpace=15)
cv2.imshow("Source", image)
cv2.imshow("Filter", dst)

cv2.waitKey(0)
cv2.destroyAllWindows()

bilateral.png

表面模糊

图像的表面模糊处理,其作用是在保留图像边缘的情况下,对图像的表面进行模糊处理。在对人物皮肤处理上,比高斯模糊更有效。(高斯模糊在使人物皮肤光洁的同时,也将一些边缘特征给模糊了)

在处理手法上,表面模糊也与前面提到的卷积处理手段不同,表面模糊是每一个像素点都有自己的卷积矩阵,而且还是 3 套,用以对应于像素的 R、G、B 分量。

表面模糊有 2 个参数,即模糊半径 r 和模糊阈值 T,模糊半径确定模糊的范围,而模糊范围确定的是卷积矩阵的大小,模糊矩阵是一个长宽相等的矩阵,长度 $l=2r+1$。

矩阵的中间元素是当前的像素点,其余的元素按照下面的方法计算:

$I_{ij}$ 是图像值,$I_0$是模板矩阵中心的图像值

一般来说,会有预处理: $w_{ij}=max(0, w_{ij})$

根据卷积运算,每个像素通过表面模糊之后的值为:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
# 表面模糊算法
import matplotlib.pyplot as plt
import numpy
from skimage import io


def sur_blur(origin, threshold, r):
transformed = origin * 1.0
row, col = origin.shape
w_size = r * 2 + 1
for i in range(r, row - 1 - r):
for j in range(r, col - 1 - r):
iij = origin[i-r: i+r+1, j-r: j+r+1]
i0 = numpy.ones([w_size, w_size]) * origin[i, j]
wij = 1 - abs(iij - i0) / (2.5 * threshold)
wij[wij < 0] = 0
tmp = iij * wij
transformed[i, j] = tmp.sum() / wij.sum()
return transformed


file_name = './material/4.jpg'
img = io.imread(file_name)
img_out = img * 1.0
boundary = 20
half_size = 10
img_out[:, :, 0] = sur_blur(img[:, :, 0], boundary, half_size)
img_out[:, :, 1] = sur_blur(img[:, :, 1], boundary, half_size)
img_out[:, :, 2] = sur_blur(img[:, :, 2], boundary, half_size)

img_out = img_out / 255

plt.subplot(121)
plt.imshow(img)
plt.axis('off')

plt.subplot(122)
plt.imshow(img_out)
plt.axis('off')
plt.savefig("surface_blur.png")
plt.show()

效果挺好,就是算的有点慢:
surface_blur.png

参考链接

Author

Ctwo

Posted on

2021-04-23

Updated on

2021-04-24

Licensed under

Comments