二变量函数的傅里叶变换

二维冲激及其取样性质

两个连续变量的冲激函数定义为：
$$\delta (t,z)=\{\begin{array}{ll}1,& t=z=0\\ 0,& \text{others}\end{array}\text{\hspace{1em}(4.52)}$$
$${\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\delta (t,z)\text{d}t\text{d}z\text{\hspace{1em}(4.53)}$$

二维冲激在积分下展现了取样性质
$${\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(t,z)\delta (t,z)\text{d}t\text{d}z=f(0,0)\text{\hspace{1em}(4.54)}$$
一般的情况
$${\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(t,z)\delta (t-t_0,z-z_0)\text{d}t\text{d}z=f(t_0,z_0)\text{\hspace{1em}(4.55)}$$

二维离散单位冲激定义为
$$\delta (x,y)=\{\begin{array}{ll}1,& x=y=0\\ 0,& \text{others}\end{array}\text{\hspace{1em}(4.56)}$$
取样性质为
$$\sum _{-\infty }^{\infty }\sum _{-\infty }^{\infty }f(x,y)\delta (x-x_0,y-y_0)\text{d}x\text{d}y=f(x_0,y_0)\text{\hspace{1em}(4.58)}$$

二维连续傅里叶变换对

$$F(\mu ,v)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(t,z)e^{-j2\pi (\mu t+vz)}\text{d}t\text{d}z\text{\hspace{1em}(4.59)}$$

$$f(t,z)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }F(\mu ,v)e^{j2\pi (\mu t+vz)}\text{d}\mu \text{d}v\text{\hspace{1em}(4.60)}$$

$\mu$和$v$是频率变量，涉及图像时，$t$和$z$解释为连续空间变量。变量$\mu$和$v$的域定义了连续频率域

# 二维盒式函数的傅里叶变换
height, width = 128, 128
m = int((height - 1) / 2)
n = int((width - 1) / 2)
f = np.zeros([height, width])
# T 控制方格的大小
T = 5
f[m-T:m+T, n-T:n+T] = 1fft = np.fft.fft2(f)
shift_fft = np.fft.fftshift(fft)
amp = np.log(1 + np.abs(shift_fft))plt.figure(figsize=(16, 16))
plt.subplot(2, 2, 1), plt.imshow(f, 'gray'), plt.title('Box filter'), plt.xticks([]), plt.yticks([])
plt.subplot(2, 2, 2), plt.imshow(amp, 'gray'), plt.title('FFT Spectrum'), plt.xticks([]), plt.yticks([])# 不同的盒式函数对应的傅里叶变换
height, width = 128, 128
m = int((height - 1) / 2)
n = int((width - 1) / 2)
f = np.zeros([height, width])
T = 20
f[m-T:m+T, n-T:n+T] = 1fft = np.fft.fft2(f)
shift_fft = np.fft.fftshift(fft)
amp = np.abs(shift_fft)
plt.subplot(2, 2, 3), plt.imshow(f, 'gray'), plt.title('Box filter'), plt.xticks([]), plt.yticks([])
plt.subplot(2, 2, 4), plt.imshow(amp, 'gray'), plt.title('FFT Spectrum'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

二维取样和二维取样定理

$$s_{\Delta T\Delta Z}(t,z)=\sum _{m=-\infty }^{\infty }\sum _{n=-\infty }^{\infty }\delta (t-m\Delta T,z-n\Delta Z)\text{\hspace{1em}(4.61)}$$

在区间$[-{\mu }_{max},{\mu }_{max}]$和$[-v_{max},v_{max}]$建立的频率域矩形之外，函数$f(t,z)$的傅里叶变换为零，即时，
$$F(\mu ,v)=0,|\mu |\ge {\mu }_{max}且|v|\ge v_{max}\text{\hspace{1em}(4.62)}$$
称该函数为带限函数。

二维取样定理：
$$\Delta T<\frac{1}{2{\mu }_{max}}\text{\hspace{1em}(4.63)}$$
和
$$\Delta Z<\frac{1}{2v_{max}}\text{\hspace{1em}(4.64)}$$

或者是：
$$\frac{1}{\Delta T}>2{\mu }_{max}\text{\hspace{1em}(4.65)}$$
和
$$\frac{1}{\Delta Z}>2v_{max}\text{\hspace{1em}(4.66)}$$

则连续带限函数$f(t,z)$可由基一组样本无误地复原。

图像中的混叠

def get_check(height, width, check_size=(5, 5), lower=130, upper=255):"""create check pattern imageheight: input, height of the image you wantwidth: input, width of the image you wantcheck_size: the check size you want, default is 5x5lower: dark color of the check, default is 130, which is dark gray, 0 is black, 255 is whiteupper: light color of the check, default is 255, which is white, 0 is blackreturn uint8[0, 255] grayscale check pattern image"""m, n = check_sizeblack = np.zeros((m, n), np.uint8)white = np.zeros((m, n), np.uint8)black[:] = lower  # darkwhite[:] = upper  # whiteblack_white = np.concatenate([black, white], axis=1)white_black = np.concatenate([white, black], axis=1)black_white_black_white = np.vstack((black_white, white_black))tile_times_h = int(np.ceil(height / m / 2))tile_times_w = int(np.ceil(width / n / 2))img_temp = np.tile(black_white_black_white, (tile_times_h, tile_times_w))img_dst = np.zeros([height, width])img_dst = img_temp[:height, :width]return img_dst

# 混叠
img_16 = get_check(512, 800, check_size=(16, 16), lower=10, upper=255)
img_16_show = img_16[:48, :96]
img_6 = get_check(512, 800, check_size=(6, 6), lower=10, upper=255)
img_6_show = img_6[:48, :96]# 16 * 0.95 = 15.2
img_095 = img_16[::15, ::15]
img_095 = np.concatenate((img_095, img_095[:, 3:]), axis=1)
img_095 = np.concatenate((img_095, img_095[:, 3:]), axis=1)
img_095 = np.concatenate((img_095, img_095[:, 3:]), axis=1)
img_095 = np.concatenate((img_095, img_095[3:, :]), axis=0)
img_095 = np.concatenate((img_095, img_095[3:, :]), axis=0)
img_095 = img_095[2:50, 2:98]# 16 * 0.48 = 7.68
img_05 = img_16[::2, ::2]  # 为了显示这里的步长选了2
img_05 = img_05[:48, :96]fig = plt.figure(figsize=(15, 8))
plt.subplot(2, 2, 1), plt.imshow(img_16_show, 'gray', vmin=0, vmax=255), plt.xticks([]), plt.yticks([])
plt.subplot(2, 2, 2), plt.imshow(img_6_show, 'gray', vmin=0, vmax=255), plt.xticks([]), plt.yticks([])
plt.subplot(2, 2, 3), plt.imshow(img_095, 'gray', vmin=0, vmax=255), plt.xticks([]), plt.yticks([])
plt.subplot(2, 2, 4), plt.imshow(img_05, 'gray', vmin=0, vmax=255), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

图像重取样和刚插

在图像被取样之前，必须在前端进行搞混叠滤波，但对于违反取样定理导致的混叠效应，事实上不存在事后能降低它的抗混叠滤波的软件。多数抗混叠的功能主要是模糊数字图像，进行降低由重取样导致的其他混叠伪影，而不能降低原取样图像中的混叠。

对数字图像降低混叠，在重取样之前 需要使用低通滤波器来平滑，以衰减数字图像的高频分量。但实际上只是平滑图像，而减少了那些令人讨厌的混叠现象。

def nearest_neighbor_interpolation(img, new_h, new_w):"""get nearest_neighbor_interpolation for image, can up or down scale image into any ratioparam: img: input image, grady image, 1 channel, shape like [512, 512]param: new_h: new image height param: new_w: new image widthreturn a nearest_neighbor_interpolation up or down scale image"""new_img = np.zeros([new_h, new_w])src_height, src_width = img.shape[:2]r = new_h / src_heightl = new_w / src_widthfor i in range(new_h):for j in range(new_w):x0 = int(i / r)y0 = int(j / l)new_img[i, j] = img[x0, y0]return new_img

def box_filter(image, kernel):""":param image: input image:param kernel: input kernel:return: image after convolution"""img_h = image.shape[0]img_w = image.shape[1]m = kernel.shape[0]n = kernel.shape[1]# paddingpadding_h = int((m -1)/2)padding_w = int((n -1)/2)image_pad = np.zeros((image.shape[0]+padding_h*2, image.shape[1]+padding_w*2), np.uint8)image_pad[padding_h:padding_h+img_h, padding_w:padding_w+img_w] = imageimage_convol = image.copy()for i in range(padding_h, img_h + padding_h):for j in range(padding_w, img_w + padding_w):temp = np.sum(image_pad[i-padding_h:i+padding_h+1, j-padding_w:j+padding_w+1] * kernel)image_convol[i - padding_h][j - padding_w] = temp # 1/(m * n) * tempreturn image_convol

# 抗混叠
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0417(a)(barbara).tif', 0)# 先缩小图像，再放大图像
img_down = nearest_neighbor_interpolation(img_ori, int(img_ori.shape[0]*0.33), int(img_ori.shape[1]*0.33))
img_up = nearest_neighbor_interpolation(img_down, img_ori.shape[0], img_ori.shape[1])# 先对原图像进行5x5的平均滤波，再缩小图像，再放大图像
kernel = np.ones([5, 5])
kernel = kernel / kernel.size
img_box_filter = box_filter(img_ori, kernel=kernel)
img_down_1 = nearest_neighbor_interpolation(img_box_filter, int(img_ori.shape[0]*0.33), int(img_ori.shape[1]*0.33))
img_up_1 = nearest_neighbor_interpolation(img_down_1, img_ori.shape[0], img_ori.shape[1])fig = plt.figure(figsize=(15, 10))
plt.subplot(1, 3, 1), plt.imshow(img_ori, 'gray'), plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 2), plt.imshow(img_up, 'gray'), plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 3), plt.imshow(img_up_1, 'gray'), plt.xticks([]), plt.yticks([])plt.tight_layout()
plt.show()

混叠和莫尔模式

莫尔模式是由近似等间隔的两个光栅叠加所产生的一种视觉现象。

def rotate_image(img, angle=45):height, width = img.shape[:2]if img.ndim == 3:channel = 3else:channel = Noneif int(angle / 90) % 2 == 0:reshape_angle = angle % 90else:reshape_angle = 90 - (angle % 90)reshape_radian = np.radians(reshape_angle)  # 角度转弧度# 三角函数计算出来的结果会有小数，所以做了向上取整的操作。new_height = int(np.ceil(height * np.cos(reshape_radian) + width * np.sin(reshape_radian)))new_width = int(np.ceil(width * np.cos(reshape_radian) + height * np.sin(reshape_radian)))if channel:new_img = np.zeros((new_height, new_width, channel), dtype=np.uint8)else:new_img = np.zeros((new_height, new_width), dtype=np.uint8)radian = np.radians(angle)cos_radian = np.cos(radian)sin_radian = np.sin(radian)# dx = 0.5 * new_width + 0.5 * height * sin_radian - 0.5 * width * cos_radian# dy = 0.5 * new_height - 0.5 * width * sin_radian - 0.5 * height * cos_radian# ---------------前向映射--------------------# for y0 in range(height):#     for x0 in range(width):#         x = x0 * cos_radian - y0 * sin_radian + dx#         y = x0 * sin_radian + y0 * cos_radian + dy#         new_img[int(y) - 1, int(x) - 1] = img[int(y0), int(x0)]  # 因为整体映射的结果会比偏移一个单位，所以这里x,y做减一操作。# ---------------后向映射--------------------dx_back = 0.5 * width - 0.5 * new_width * cos_radian - 0.5 * new_height * sin_radiandy_back = 0.5 * height + 0.5 * new_width * sin_radian - 0.5 * new_height * cos_radianfor y in range(new_height):for x in range(new_width):x0 = x * cos_radian + y * sin_radian + dx_backy0 = y * cos_radian - x * sin_radian + dy_backif 0 < int(x0) <= width and 0 < int(y0) <= height:  # 计算结果是这一范围内的x0，y0才是原始图像的坐标。new_img[int(y), int(x)] = img[int(y0) - 1, int(x0) - 1]  # 因为计算的结果会有偏移，所以这里做减一操作。#     # ---------------双线性插值--------------------
#     if channel:
#         fill_height = np.zeros((height, 2, channel), dtype=np.uint8)
#         fill_width = np.zeros((2, width + 2, channel), dtype=np.uint8)
#     else:
#         fill_height = np.zeros((height, 2), dtype=np.uint8)
#         fill_width = np.zeros((2, width + 2), dtype=np.uint8)
#     img_copy = img.copy()
#     # 因为双线性插值需要得到x+1，y+1位置的像素，映射的结果如果在最边缘的话会发生溢出，所以给图像的右边和下面再填充像素。
#     img_copy = np.concatenate((img_copy, fill_height), axis=1)
#     img_copy = np.concatenate((img_copy, fill_width), axis=0)
#     for y in range(new_height):
#         for x in range(new_width):
#             x0 = x * cos_radian + y * sin_radian + dx_back
#             y0 = y * cos_radian - x * sin_radian + dy_back
#             x_low, y_low = int(x0), int(y0)
#             x_up, y_up = x_low + 1, y_low + 1
#             u, v = np.modf(x0)[0], np.modf(y0)[0]  # 求x0和y0的小数部分
#             x1, y1 = x_low, y_low
#             x2, y2 = x_up, y_low
#             x3, y3 = x_low, y_up
#             x4, y4 = x_up, y_up
#             if 0 < int(x0) <= width and 0 < int(y0) <= height:
#                 pixel = (1 - u) * (1 - v) * img_copy[y1, x1] + (1 - u) * v * img_copy[y2, x2] + u * (1 - v) * img_copy[y3, x3] + u * v * img_copy[y4, x4]  # 双线性插值法，求像素值。
#                 new_img[int(y), int(x)] = pixelreturn new_img

# 混叠和莫尔模式
# 竖线原图
img_lines = np.ones([129, 129]) * 255
img_lines[:, ::3] = 0# 旋转时使用刚内插，产生了混叠
rotate_matrix = cv2.getRotationMatrix2D((int(img_lines.shape[0]*0.5), int(img_lines.shape[1]*0.5)), -5, 1)
img_lines_r = cv2.warpAffine(img_lines, rotate_matrix, dsize=(img_lines.shape[0], img_lines.shape[1]), flags=cv2.INTER_CUBIC, borderValue=255)# 相加后，产生的混叠更明显
img_add = img_lines + img_lines_r
img_add = np.uint8(normalize(img_add) * 255)fig = plt.figure(figsize=(15, 10))
plt.subplot(2, 3, 1), plt.imshow(img_lines, 'gray', vmin=0, vmax=255), plt.xticks([]), plt.yticks([])
plt.subplot(2, 3, 2), plt.imshow(img_lines_r, 'gray', vmin=0, vmax=255), plt.xticks([]), plt.yticks([])
plt.subplot(2, 3, 3), plt.imshow(img_add, 'gray', vmin=0, vmax=255), plt.xticks([]), plt.yticks([])# 格子原图
img_check = np.ones([129, 129]) * 255
img_check[::2, ::2] = 0# 旋转时使用刚内插，产生了混叠
rotate_matrix = cv2.getRotationMatrix2D((int(img_check.shape[0]*0.5), int(img_check.shape[1]*0.5)), -5, 1)
img_check_r = cv2.warpAffine(img_check, rotate_matrix, dsize=(img_check.shape[0], img_check.shape[1]), flags=cv2.INTER_CUBIC, borderValue=255)# 相加后，产生的混叠更明显
img_add = img_check + img_check_r
img_add = np.uint8(normalize(img_add) * 255)plt.subplot(2, 3, 4), plt.imshow(img_check, 'gray', vmin=0, vmax=255), plt.xticks([]), plt.yticks([])
plt.subplot(2, 3, 5), plt.imshow(img_check_r, 'gray', vmin=0, vmax=255), plt.xticks([]), plt.yticks([])
plt.subplot(2, 3, 6), plt.imshow(img_add, 'gray', vmin=0, vmax=255), plt.xticks([]), plt.yticks([])plt.tight_layout()
plt.show()

# 印刷采用欠取样时产生莫尔模式效应
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0421(car_newsprint_sampled_at_75DPI).tif', 0)fig = plt.figure(figsize=(15, 10))
plt.subplot(1, 3, 1), plt.imshow(img_ori, 'gray'), plt.xticks([]), plt.yticks([])plt.tight_layout()
plt.show()

二维离散傅里叶变换及其反变换

二维DFT
$$F_{u,v}=\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}f(x,y)e^{-j2\pi (ux/M+vy/N)}\text{\hspace{1em}(4.67)}$$

二维IDFT
$$f(x,y)=\frac{1}{MN}\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}F(u,v)e^{j2\pi (ux/M+vy/N)}\text{\hspace{1em}(4.68)}$$

上式，$u=0,1,2,\cdots ,M-1$和$v=0,1,2,\cdots ,N-1$，$x=0,1,2,\cdots ,M-1$和$y=0,1,2,\cdots ,N-1$