天津中医药大学-20级医学信息工程 教师：王翌 学生：邓集亲
声明：本文章中所涉及的代码并非我个人独立完成的成果。
在撰写的过程中，我广泛地吸取了前一辈人，尤其是学长学姐们的宝贵学习经验。通过深入研究他们的学习轨迹，以及查阅和分析了众多相关的理论文献与资料，并在王老师的悉心指导下，经过反复的实验、调试与优化，最终得以总结完成本文所展现的代码。
建议各位学弟学妹们，先根据王老师的授课内容，独立思考一下应该如何完成。如果实在是有理解上的困难，不知道从何下手，再来学习和参考本文。
学长我是用Python写的，如果你使用MATLAB，同样可以参考此代码进行翻译。我在代码中加入了许多注释和测试环节，以便于理解。由于学长能力有限，代码中或许存在一些疏漏或错误，请批判性地参考。
实验六：低通滤波器
作业要求：
- 参考频率域处理和空间域处理课程的内容，分别采用频率域和空间域的若干种低通滤波器对图片进行滤波操作，取半径为5、15、30、80和230，分别输出结果图片。
例子：
不同半径的高斯低通滤波器的滤波效果：
实验图片：
filter.jpg，再自选其它至少20张图片（包括灰度图片和彩色图片）。
低通滤波是一种常见的图像处理方法，其原理是通过去除高频信号和噪声，保留低频信号来平滑图像。低通滤波器可以应用于多种图像处理任务，如降噪、模糊、平滑等。
低通滤波器的基本原理是在空间域或频率域中将高频分量抑制或消除，从而实现平滑效果。在空间域中，低通滤波器可以通过卷积运算实现。常见的低通滤波器有均值滤波器和高斯滤波器。
在频率域中，低通滤波器可以通过将图像进行傅里叶变换，然后将高频分量截断或抑制来实现。这种方法通常被称为频域滤波。在频率域中，低通滤波器的效果可以通过掩模函数来可视化。
Windows图片文件目录：
ImageSet文件目录下储存其他待处理图片
Python代码：
import cv2 as cv
import numpy as np
from matplotlib import pyplot
from numpy.core import integer, asarray, roll,zeros, swapaxes, take
from numpy.core.multiarray import normalize_axis_index
from numpy.fft import _pocketfft_internal as pfi
def imagefft2(a, s=None, axes=(-2, -1)): #图片二维傅里叶变换
return raw_fftnd(a, s, axes, fft)
def imageifft2(a, s=None, axes=(-2, -1)):
return raw_fftnd(a, s, axes, ifft)
def fft_shift(x, axes=None): #将图像低频部分转移到图像的中心
x = asarray(x)
if axes is None:
axes = tuple(range(x.ndim))
shift = [dim // 2 for dim in x.shape]
elif isinstance(axes, integer):
shift = x.shape[axes] // 2
else:
shift = [x.shape[ax] // 2 for ax in axes]
return roll(x, shift, axes)
def ifft_shift(x, axes=None): #将图像低频部分移到原来位置
x = asarray(x)
if axes is None:
axes = tuple(range(x.ndim))
shift = [-(dim // 2) for dim in x.shape]
elif isinstance(axes, integer):
shift = -(x.shape[axes] // 2)
else:
shift = [-(x.shape[ax] // 2) for ax in axes]
return roll(x, shift, axes)
def fft(a, n=None, axis=-1): #定义一维傅里叶变换函数
a = asarray(a)
if n is None:
n = a.shape[axis]
inv_norm = 1
output = raw_fft(a, n, axis, False, True, inv_norm)
return output
def ifft(a, n=None, axis=-1):
a = asarray(a)
if n is None:
n = a.shape[axis]
inv_norm = n
output = raw_fft(a, n, axis, False, False, inv_norm)
return output
def raw_fft(a, n, axis, is_real, is_forward, inv_norm): #为了变换结果，这里需要定义浮点数inv_norm
axis = normalize_axis_index(axis, a.ndim)
if n is None:
n = a.shape[axis]
fct = 1/inv_norm
if a.shape[axis] != n:
s = list(a.shape)
index = [slice(None)]*len(s)
if s[axis] > n:
index[axis] = slice(0, n)
a = a[tuple(index)]
else:
index[axis] = slice(0, s[axis])
s[axis] = n
z = zeros(s, a.dtype.char)
z[tuple(index)] = a
a = z
if axis == a.ndim-1:
r = pfi.execute(a, is_real, is_forward, fct)
else:
a = swapaxes(a, axis, -1)
r = pfi.execute(a, is_real, is_forward, fct)
r = swapaxes(r, axis, -1)
return r
def cook_nd_args(a, s=None, axes=None, invreal=0):
if s is None:
shapeless = 1
if axes is None:
s = list(a.shape)
else:
s = take(a.shape, axes)
else:
shapeless = 0
s = list(s)
if axes is None:
axes = list(range(-len(s), 0))
if invreal and shapeless:
s[-1] = (a.shape[axes[-1]] - 1) * 2
return s, axes
def raw_fftnd(a, s=None, axes=None, function=fft):
a = asarray(a)
s, axes = cook_nd_args(a, s, axes)
itl = list(range(len(axes)))
itl.reverse()
for ii in itl:
a = function(a, n=s[ii], axis=axes[ii])
return a
#读取图像
original = cv.imread(r'D:\deng\smalljob\ImageSet\filter.jpg')
original_gray = cv.cvtColor(original, cv.COLOR_BGR2GRAY)
#傅里叶变换
f = imagefft2(original_gray)
fshift = fft_shift(f)
#创建存储图像的数组
all_images=[[]]*5
r=[5,15,30,80,230]
#低通蒙板
rows,cols = fshift.shape
mid_x,mid_y = int(rows/2),int(cols/2)
mask = np.zeros((rows,cols),dtype=np.uint8)
m=0
for n in r:
mask[mid_x-n:mid_x+n,mid_y-n:mid_y+n] = 1
fshift2 = mask*fshift
#傅里叶逆变换
ishift = ifft_shift(fshift2)
all_images[m] = imageifft2(ishift)
all_images[m] = np.abs(all_images[m])
m=m+1
#展示结果
pyplot.subplot(231),pyplot.imshow(original_gray,cmap='gray')
pyplot.title('original')
pyplot.axis('off')
pyplot.subplot(232),pyplot.imshow(all_images[0],cmap='gray')
pyplot.title('r=5')
pyplot.axis('off')
pyplot.subplot(233),pyplot.imshow(all_images[1],cmap='gray')
pyplot.title('r=15')
pyplot.axis('off')
pyplot.subplot(234),pyplot.imshow(all_images[2],cmap='gray')
pyplot.title('r=30')
pyplot.axis('off')
pyplot.subplot(235),pyplot.imshow(all_images[3],cmap='gray')
pyplot.title('r=80')
pyplot.axis('off')
pyplot.subplot(236),pyplot.imshow(all_images[4],cmap='gray')
pyplot.title('r=230')
pyplot.axis('off')
pyplot.show()
运行结果：（图像变灰是因为高频被过滤掉了）
以上是处理灰度图，以下是处理彩色图像
import cv2 as cv
import numpy as np
from matplotlib import pyplot
from numpy.core import integer, asarray, roll,zeros, swapaxes, take
from numpy.core.multiarray import normalize_axis_index
from numpy.fft import _pocketfft_internal as pfi
def fft_shift(x, axes=None): #将图像低频部分转移到图像的中心
x = asarray(x)
if axes is None:
axes = tuple(range(x.ndim))
shift = [dim // 2 for dim in x.shape]
elif isinstance(axes, integer):
shift = x.shape[axes] // 2
else:
shift = [x.shape[ax] // 2 for ax in axes]
return roll(x, shift, axes)
def ifft_shift(x, axes=None): #将图像低频部分移到原来位置
x = asarray(x)
if axes is None:
axes = tuple(range(x.ndim))
shift = [-(dim // 2) for dim in x.shape]
elif isinstance(axes, integer):
shift = -(x.shape[axes] // 2)
else:
shift = [-(x.shape[ax] // 2) for ax in axes]
return roll(x, shift, axes)
def raw_fft(a, n, axis, is_real, is_forward, inv_norm): #为了变换结果，这里需要定义浮点数inv_norm
axis = normalize_axis_index(axis, a.ndim)
if n is None:
n = a.shape[axis]
fct = 1/inv_norm
if a.shape[axis] != n:
s = list(a.shape)
index = [slice(None)]*len(s)
if s[axis] > n:
index[axis] = slice(0, n)
a = a[tuple(index)]
else:
index[axis] = slice(0, s[axis])
s[axis] = n
z = zeros(s, a.dtype.char)
z[tuple(index)] = a
a = z
if axis == a.ndim-1:
r = pfi.execute(a, is_real, is_forward, fct)
else:
a = swapaxes(a, axis, -1)
r = pfi.execute(a, is_real, is_forward, fct)
r = swapaxes(r, axis, -1)
return r
def fft(a, n=None, axis=-1): #定义一维傅里叶变换函数
a = asarray(a)
if n is None:
n = a.shape[axis]
inv_norm = 1
output = raw_fft(a, n, axis, False, True, inv_norm)
return output
def ifft(a, n=None, axis=-1):
a = asarray(a)
if n is None:
n = a.shape[axis]
inv_norm = n
output = raw_fft(a, n, axis, False, False, inv_norm)
return output
def cook_nd_args(a, s=None, axes=None, invreal=0):
if s is None:
shapeless = 1
if axes is None:
s = list(a.shape)
else:
s = take(a.shape, axes)
else:
shapeless = 0
s = list(s)
if axes is None:
axes = list(range(-len(s), 0))
if invreal and shapeless:
s[-1] = (a.shape[axes[-1]] - 1) * 2
return s, axes
def raw_fftnd(a, s=None, axes=None, function=fft):
a = asarray(a)
s, axes = cook_nd_args(a, s, axes)
itl = list(range(len(axes)))
itl.reverse()
for ii in itl:
a = function(a, n=s[ii], axis=axes[ii])
return a
def imagefft2(a, s=None, axes=(-2, -1)): #图片二维傅里叶变换
return raw_fftnd(a, s, axes, fft)
def imageifft2(a, s=None, axes=(-2, -1)):
return raw_fftnd(a, s, axes, ifft)
#读取图像
original = cv.imread(r'D:\deng\smalljob\ImageSet\ying.jpg')
#original_gray = cv.cvtColor(original, cv.COLOR_BGR2GRAY) #将图像从RGB颜色空间转换到灰度颜色空间
original = cv.cvtColor(original, cv.COLOR_BGR2RGB)
b,g,r=cv.split(original)
b_1=imagefft2(b)
g_1=imagefft2(g)
r_1=imagefft2(r)
b_fshift=fft_shift(b_1)
g_fshift=fft_shift(g_1)
r_fshift=fft_shift(r_1)
#傅里叶变换
#f = imagefft2(original_gray)
#fshift = fft_shift(f) #将变换后图像的低频部分转移到图像的中心
#创建存储图像的数组
all_images=[[]]*5
allimg=[[]]*5
r=[5,15,30,80,230]
#低通蒙板
rows,cols = b.shape
mid_x,mid_y = int(rows/2),int(cols/2)
mask = np.zeros((rows,cols),dtype=np.uint8)
m=0
for n in r:
mask[mid_x-n:mid_x+n,mid_y-n:mid_y+n] = 1
#fshift2 = mask*fshift
b_fshift2=mask*b_fshift
g_fshift2=mask*g_fshift
r_fshift2=mask*r_fshift
b_ishift = ifft_shift(b_fshift2)
g_ishift = ifft_shift(g_fshift2)
r_ishift = ifft_shift(r_fshift2)
# 傅里叶逆变换后，对结果进行归一化处理
b_2 = np.clip(imageifft2(b_ishift).real, 0, 255).astype(np.uint8)
g_2 = np.clip(imageifft2(g_ishift).real, 0, 255).astype(np.uint8)
r_2 = np.clip(imageifft2(r_ishift).real, 0, 255).astype(np.uint8)
allimg[m] = np.dstack((b_2, g_2, r_2))
allimg[m]=np.abs(allimg[m]).astype(np.float32)/255.0
#傅里叶逆变换
#ishift = ifft_shift(fshift2)
#all_images[m] = imageifft2(ishift)
#all_images[m] = np.abs(all_images[m])
m=m+1
#展示结果
pyplot.subplot(231),pyplot.imshow(original)
pyplot.title('original')
pyplot.axis('off')
pyplot.subplot(232),pyplot.imshow(allimg[0])
pyplot.title('r=5')
pyplot.axis('off')
pyplot.subplot(233),pyplot.imshow(allimg[1])
pyplot.title('r=15')
pyplot.axis('off')
pyplot.subplot(234),pyplot.imshow(allimg[2])
pyplot.title('r=30')
pyplot.axis('off')
pyplot.subplot(235),pyplot.imshow(allimg[3])
pyplot.title('r=80')
pyplot.axis('off')
pyplot.subplot(236),pyplot.imshow(allimg[4])
pyplot.title('r=230')
pyplot.axis('off')
pyplot.show()
运行结果：
747
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
