数字图像处理期末复习

准备数字图像处理期末考试，根据老师重点整理。

一、深入理解并熟悉直方图均衡化过程，完成完整计算

image.png

1. 求出图像f的总体像素个数
   Nf = m*n (m,n分别为图像的长和宽）
2. 计算每个灰度级的分布概率，即每个像素在整个图像中所占的比例。
   hs(i)=h(i)/Nf (i=0,1,…,255)
   hs=[0.12,0,08,0.16,0.16,0.04,0.04,0.16,0.04,0.08,0.12]
3. 计算原图灰度的累计分布
   hp=[0.12, 0.20, 0.36, 0.52, 0.56, 0.60, 0.76, 0.80, 0.88, 1.00]
4. 计算原、新图灰度值的影射关系
   新图象的灰度值为g(i,j)=255*hp(k) (f(x,y)!==0) / 0(f(x,y)==0)

二、 熟悉模板卷积过程并能利用给出的模板实现滤波运算，熟悉各种模板实现的滤波功能。

模板运算

(1) 将模板在输入图象中漫游，将模板中心与图象中某个像素位置重合
(2) 将模板上各个系数与模板下各对象的灰度值相乘
(3) 将所有乘积相加/模板系数之和 赋予图象中对于模板中心位置的像素

线性平滑滤波

1. 邻域平均例如3x3: [[1,1,1],[1,1,1],[1,1,1]]
2. 加权平均中心系数大，周围系数小
3. 高斯平均根据高斯分布确定各模板系数

线性锐化滤波

1. 拉普拉斯算子[0,-1,0],[-1,-4,-1],[0,-1,0]或[-1,-1,-1],[-1,-8,-1],[-1,-1,-1]

噪声平滑&边缘检测

Priwitt 算子
[[-1,-1,-1],[0,0,0],[1,1,1]]或[[-1,0,1],[-1,0,1],[-1,0,1]]
Sobel 边缘检测算子
[[-1,-2,-1],[0,0,0],[1,2,1]]或[[-1,0,1],[-2,0,2],[-1,0,1]]

三、Canny算法

image.png
image.png
image.png

Canny算法关键函数

gaussian_smooth(image, rows, cols, sigma, &smoothedim);  //高斯平滑
derrivative_x_y(smoothedim, rows, cols, &delta_x, &delta_y);  //微分
magnitude_x_y(delta_x, delta_y, rows, cols, &magnitude);  //计算二维矢量的幅值
non_max_supp(magnitude, delta_x, delta_y, rows, cols, nms);  //非极大值抑制

熟悉Hough变换流程，理解Hough的应用领域(冈萨p475)

在xy平面上的一点(x0,y0)要转换为极坐标下的(θ,ρ)：
x0 = ρcosθ，y0 = ρsinθ，那么x0cosθ = ρcosθcosθ，y0sinθ = ρsinθsinθ
可得 x0cosθ + y0sinθ = ρ 那么xy坐标平面到极坐标θρ平面的转换公式为：xcosθ + ysinθ = ρ
可以看出xy平面上的一点在极坐标θρ平面上的投影为一条正弦曲线，那么n个点对应n条正弦曲线，这n条曲线如果有m条曲线线交于一点(θ0,ρ0)，那么这m条曲线对应在xy平面下的m个点就在同一直线上，而θ0就是这条直线的法线与x轴的夹角，ρ0为原点到直线的距离。
OpenCV自带的标准霍夫变换函数：

icvHoughLinesStandard( const CvMat* img, float rho, float theta, int threshold, CvSeq *lines, int linesMax )

步骤为：
1、将θ离散值的sin、cos的计算结果存在数组中，便于下一步使用。
2、将原图中灰度值不为0的点变换到θρ空间的正弦曲线，而θρ空间相当于一个累加器。
3、将累加器的值进行排序，使用的是最大邻域值法，可以看出邻域很小，邻域小计算速度就快。
4、输出累加器中值最大的几个值。

应用:

Works on Disconnected Edges处理不连续边缘
Relatively insensitive to occlusion对遮掩区域不敏感
Effective for simple shapes (lines, circles, etc)对简单形状很有效
Trade-off between work in Image Space and Parameter Space 图像空间与参数空间中工作的权衡
Handling inaccurate edge locations处理不准确的边缘位置:
Increment Patch in Accumulator rather than a single point累加器中的增量补片而不是单个点

四、熟悉链码跟踪原理及编码过程

1、freeman链码简介
首先给出 Freeman编码的定义描述: 任选一个像素点(通常对已细化的图像进行)作为参考点,与其相邻的像素分别在 8 个不 同 的 位 置上, 给它们赋予方向值 0～7(如图 1), 称为0～7 位链码方向值, 一个线条可以用Freeman 链码 的 码值串来表示称为该线条图形的链码。如图 2, 给出一个9 ×9 的点阵图, 其中一条线段, S 为起始点, E 为终点,此线段可做如表示为 L =43322100000066。

image.png

2. Freeman 链码基础上的矢量跟踪方法
   针对线条图像的链码识别问题,Freeman 提出了直线链码应满足的三个基本条件:1) 只包含最多两个方向的码值。2) 该两码值之一总是单独出现。3) 该单独出现的码值要尽可能均匀的分布在整个链码中。也就是说, 在绘制或编历一条直线的时候, 只需要进行两个方向的移动, 其中一个方向一定是独立出现的。在理想状态下, 当在一个单独出现的方向上识别到某像素后, 下一个像素点的方向就是另一个链码的方向。

Freeman 准则的直线检测算法通过跟踪线段子元，并根据两相邻线段子元间的偏转角度判断是否需要进行子线段合并。

Freeman链码算法设计：
参考论文： 尚振宏等《运用Freeman准则的直线检测算法》；史册《对一种快速边缘跟踪算法的讨论》
输入为二值图像,输出为图像中物体边界的各条线段·
整个算法流程如下所述,其中,Line Sub Cell、Line Cell和Line Segment分别表示存储线段子元、线段元和线段的结构·
Step1·从扫描二值图像,定位边界跟踪起点Pi,i=0·将Pi作为第一个像素保存到线段子元Line Sub Cell中·
Step2·i=i+1,利用文献[6]中提出的边界跟踪算法顺时针(或逆时针)跟踪得到下一边界点Pi,若跟踪回到起始点,则算法结束;否则,根据定义1确定Pi是否属于当前线段子元Line Sub Cell·若属于,则将Pi保存到Line SubCell中,转Step2;否则,当前线段子元Line Sub Cell即为一完整的线段子元,执行Step3·
Step3·若线段元Line Cell中无线段子元,则将LineSub Cell作为Line Cell的第一线段子元加入Line Cell;否则,根据定义2判断线段子元Line Sub Cell是否属于线段元Line Cell·若属于,将Line Sub Cell加入Line Cell,然后将Line Sub Cell置空并将Pi加入Line Sub Cell中,转Step2;否则,Line Cell为一完整的线段元,执行Step4·
Step4·若当前线段Line Segment中无线段元,则将Line Cell作为Line Segment的第一线段元加入LineSegment;否则根据下面的准则判断Line Cell是否属于Line Segment·若不属于,则Line Segment为一完整直线段,保存该直线段并置空Line Segment·然后,无论LineCell是否属于Line Segment均将Line Cell加入LineSegment,并分别用当前的线段子元Line Sub Cell和点Pi初始化线段元Line Cell和线段子元Line Sub Cell并转回Step2·
判断线段元Line Cell是否属于线段LineSegment的准则:见论文

五、熟悉形态学操作计算过程，掌握腐蚀、膨胀、开和闭运算（冈萨p420）

1. 腐蚀是一种消除边界点，使边界向内部收缩的过程。可以用来消除小且无意义的物体。
   腐蚀的算法：
   用3x3的结构元素，扫描图像的每一个像素
   用结构元素与其覆盖的二值图像做“与”操作
   如果都为1，结果图像的该像素为1。否则为0。
   结果：使二值图像减小一圈
2. 膨胀是将与物体接触的所有背景点合并到该物体中，使边界向外部扩张的过程。可以用来填补物体中的空洞。
   膨胀的算法：
   用3x3的结构元素，扫描图像的每一个像素
   用结构元素与其覆盖的二值图像做“与”操作
   如果都为0，结果图像的该像素为0。否则为1
   结果：使二值图像扩大一圈
3. 先腐蚀后膨胀的过程称为开运算。用来消除小物体、在纤细点处分离物体、平滑较大物体的边界的同时并不明显改变其面积。
4. 先膨胀后腐蚀的过程称为闭运算。用来填充物体内细小空洞、连接邻近物体、平滑其边界的同时并不明显改变其面积。

六、熟悉各种图像分割算法，讲得清算法原理和步骤，局限性。(冈萨p460)

p-参数法

对固定分辨率下的目标物，根据目标物在画面中所占的比例来选择阈值，进行二值化处理。
步骤：1. 先试探性地给出一个阈值（黄色） ，统计目标物的像素点数在整幅图中所占的比例是否满足要求，是则阈值合适；2. 否则，阈值则偏大（右）或者偏小（左），再进行调整，直到满足要求（白色）。

均匀性度量法

步骤： 1）给定一个初始阈值 Th=Th0 例如：可以默认为 1，或者是 128 等），则将原图分为 C1 和 C2 两类；2）分别计算两类的类内方差 3）分别计算两类像素在图像中的分布概率：4）选择最佳阈值 Th=Th*，使得下式成立

image.png

聚类方法

• 基于区域的分割

  • 区域增长算法：流程，会计算
    • 定义：区域生长是把图像分割成特征相似的若干小区域，比较相邻小区域的特征，若相似则合并为同一区域，如此进行直到不能合并为止，最后生成特征不同的各区域。也称为区域扩张法。
    • 选择合适的种子点
      确定相似性准则(生长准则)
      确定生长停止条件
      采用的判断准则是：如果所考虑的象素与种子象素灰度值差的绝对值小于某个门限 T，则将该象素包括进种子象素所在的区域
    • 生长准则：基于区域灰度差 基于区域灰度分布统计性质 基于区域形状

• 区域分裂与合并：会计算怎么分裂怎么合并

  • 另外一种分割的想法: 先从整幅图像开始通过不断分裂，得到各个区域（在实际中，先将图像分成任意大小且不重叠的区域，然后再合并或分裂这些区域，以满足分割的要求），在这类方法中，常根据图像的统计特性设定图像区域属性的一致性测度
  • 算法步骤：
    （1）对任一区域 Ri，如果 P(Ri)=FALSE, 就将其分裂成不重叠的四等分
    （2）对相邻的两个区域 Ri 和 Rj（它们可以大小不同，即不在同一层），如果条件 P(Ri∪Rj)=TRUE，就将它们合并
    （3）如果进一步的分裂或合并都不可能，则结束

• 分水岭分割算法：
  • 基于距离的分水岭分割算法：会计算二值图像距离图像
    • 二值图像的距离变换：每一个像素到最近非零像素的距离。
      值为 1 的像素距离变换为 0。

  - 局域图像梯度的分水岭分割算法：掌握梯度图像平滑的方法：开运算和闭运算。
        - 梯度幅度图像在沿对象边缘处有较高的值，而在其它地方值较小。

理想情况下，在沿对象边缘处产生分水岭。

• 基于图像标记的分水岭分割算法
  上述方法通常会由于噪声和其它诸如梯度的局部不规则性的影响造成过度分割。而过度分割足以令结果变得毫无用处。
  解决方法：通过一系列的预处理步骤，加入附加知识，限制允许存在的区域数目。
  标记符属于图像的连通分量。内部标记符标记感兴趣对象的内部；外部标记符标记背景区域。

基于不连续点的分割

*   点检测：如果一个孤立点与它周围的点不同，则可以使用上述模板进行检测。
*   点检测的另外一种方法是在m*n大小的邻域中，
    找到其最大像素值点和最小像素值点，其差值大于阈值的那些点则可认为是图像中的孤立点。
    g= ordfilter2(f，m*n，ones(m,n))-ordfilter2(f,1,ones(m*n));
    g= g>=T
*   线检测，4个模板 6个【-1】和3个【2】的排列组合；检测水平数值45度。

• Canny 算子的检测的流程及实现
  Step1：利用高斯滤波对图像平滑 （去噪）
  Step2：利用梯度算子的计算图像的梯度和梯度方向（求导）
  Step3: 进行非极大值抑制（初步得到边缘点）
  Step4: 双阈值连接边缘

七、熟悉自己做过的所有课内实验及所用的函数。

实验1 灰度直方图

import cv2
import numpy as np
import matplotlib.pyplot as plt

if __name__ == "__main__":
     
    img_path="3.jpg"
    img=cv2.imread(img_path)
    cv2.imshow("img",img)
    img_gray=cv2.cvtColor(img,cv2.COLOR_RGB2GRAY)//转成灰度图
    cv2.imshow("img_g",img_gray)
    equ=cv2.equalizeHist(img_gray)
    cv2.imshow("img_equalizeHist",equ)//直方图均衡化
    print(img_gray.ravel())
    print(equ.ravel())
    plt.hist(img_gray.ravel(), 256 )
    plt.show()
    plt.hist(equ.ravel(), 256)
    plt.show()
    cv2.imwrite('img_g.jpg',img_gray)
    cv2.imwrite('img_equ.jpg',equ)
    cv2.waitKey()

实验二 卷积，滤波

import cv2
import numpy as np
import matplotlib.pyplot as plt

if __name__ == "__main__":
     
    img_path="3.png"
    img=cv2.imread(img_path)
//设置卷积核
    kernel_sharpen_1 = np.array([
        [0,1,0],
        [1,-4,1],
        [0,1,0]])
    kernel_sharpen_2 = np.array([
            [-1,-2,-1],
            [0,0,0],
            [1,2,1]])
    kernel_sharpen_3 = np.array([
            [-1,0,1],
            [-1,0,1],
            [-1,0,1]])
    img_gauss = cv2.GaussianBlur(img,(3,3),1) //高斯滤波
    img_blur = cv2.blur(img,(3,3))//均值滤波
    img_medianblur = cv2.medianBlur(img,5)//中值滤波
 /* Priwitt 算子 [[-1,-1,-1],[0,0,0],[1,1,1]]或[[-1,0,1],[-1,0,1],[-1,0,1]] 
Sobel 边缘检测算子 [[-1,-2,-1],[0,0,0],[1,2,1]]或[[-1,0,1],[-2,0,2],[-1,0,1]]
拉普拉斯算子[0,-1,0],[-1,-4,-1],[0,-1,0]或[-1,-1,-1],[-1,-8,-1],[-1,-1,-1]*/
    output_1 = cv2.add(cv2.filter2D(img,-1,kernel_sharpen_1),img)
    output_2 = cv2.add(cv2.filter2D(img,-1,kernel_sharpen_2),img)
    output_3 = cv2.add(cv2.filter2D(img,-1,kernel_sharpen_3),img)

    cv2.imshow("img",img)
    cv2.imshow("img_gauss",img_gauss)
    cv2.imshow("img_blur",img_blur)
    cv2.imshow("img_medianblur",img_medianblur)
    cv2.imshow("img_laplace",output_1)
    cv2.imshow("img_Sobelx",output_2)
    cv2.imshow("img_Prewitty",output_3)

    cv2.waitKey()

实验三 分割、Hough变换

import cv2 
import numpy


def watershed_demo():
    # remove noise if any 消除噪声
    print(src.shape)
    # gray, binary image
    blurred = cv.pyrMeanShiftFiltering(src, 10, 100)
    gray = cv.cvtColor(blurred, cv.COLOR_BGR2GRAY)
    ret, binary = cv.threshold(gray, 0, 255, cv.THRESH_BINARY | cv.THRESH_OTSU)
#     cv.imshow("binary image", binary)

    # morphology operation 开操作去除噪点
    kernel = cv.getStructuringElement(cv.MORPH_RECT, (3, 3))
    # morphology binary,2次开操作
    mb = cv.morphologyEx(binary, cv.MORPH_OPEN, kernel, iterations=2)
    sure_bg = cv.dilate(mb, kernel, iterations=3)  # 3次膨胀
#     cv.imshow("mor-opt", sure_bg)

    # distance transform
    # distance transform
    # DIST_L1:曼哈顿距离，DIST_L2：欧氏距离, masksize:跟卷积一样
    # 这是我们获取的字段距离数值，对应每个像素都有，所以数组结构和图像数组一致
    dist = cv.distanceTransform(mb, cv.DIST_L2, 3)
    dist_output = cv.normalize(dist, 0, 1.0, cv.NORM_MINMAX)  # 归一化的距离图像数组 0-1之间标准化
    # cv.imshow("distance-t", dist_output*50)    # 这行代码运行不了

    ret, surface = cv.threshold(dist, dist.max()*0.5, 255, cv.THRESH_BINARY)
    # cv.imshow("surface_bin", surface)   # 个人运行不了
    # 计算marker
    surface_fg = numpy.uint8(surface)  # 计算前景
    unknown = cv.subtract(sure_bg, surface_fg)  # 计算未知区域
    ret, markers = cv.connectedComponents(surface_fg)
    print(ret)

    # watershed transform 分水岭变换
    markers = markers + 1   # 用label进行控制
    markers[unknown == 255] = 0
    markers = cv.watershed(src, markers=markers)  # 分水岭的地方就编程-1
    src[markers == -1] = [0, 255, 0]
    cv.imshow("img_watershed", src)


src = cv2.imread("4.jpg", cv.IMREAD_COLOR)
watershed_demo()

img = cv2.imread("4.jpg")
#初始化slic项，超像素平均尺寸20（默认为10），平滑因子20
slic = cv2.ximgproc.createSuperpixelSLIC(img,region_size=20,ruler = 20.0) 
slic.iterate(10)     #迭代次数，越大效果越好
mask_slic = slic.getLabelContourMask() #获取Mask，超像素边缘Mask==1
label_slic = slic.getLabels()        #获取超像素标签
number_slic = slic.getNumberOfSuperpixels()  #获取超像素数目
mask_inv_slic = cv2.bitwise_not(mask_slic)  
img_slic = cv2.bitwise_and(img,img,mask =  mask_inv_slic) #在原图上绘制超像素边界
cv2.imshow("img_slic",img_slic)

img = cv2.imread('4.jpg')
cv2.imshow("img", img)
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
edges = cv2.Canny(gray, 50, 120)
minLineLength = 100 # 最小直线长度，太小舍去
maxLineGap = 5 # 最大线段间隙， 太大舍去
lines = cv2.HoughLinesP( # P意味着概率，所谓“概率版本的Hough变换”
    edges, 1, np.pi/180, 100, minLineLength, maxLineGap
)
for each in lines[0: 30]: # 绘制直线，十条
    for x1, y1, x2, y2 in each:
        cv2.line(img, (x1, y1), (x2, y2), (255, 0, 0), 2)

        
cv2.imwrite('line-cycle-detect.jpg', img)
# cv2.imshow('edges', edges)
# cv2.imshow('lines', img)

planets = img
gray_img = cv2.cvtColor(planets, cv2.COLOR_BGR2GRAY)
img = cv2.medianBlur(gray_img, 5)
cimg = cv2.cvtColor(img,cv2.COLOR_GRAY2BGR)

circles = cv2.HoughCircles(img,cv2.HOUGH_GRADIENT,1,120,
                            param1=100,param2=30,minRadius=0,maxRadius=0)

circles = np.uint16(np.around(circles))

for i in circles[0,:]:
    # draw the outer circle
    cv2.circle(planets,(i[0],i[1]),i[2],(0,255,0),2)
    # draw the center of the circle
#     cv2.circle(planets,(i[0],i[1]),2,(0,0,255),3)

cv2.imwrite("planets_circles.jpg", planets)
cv2.imshow("img_hough", planets)
cv2.waitKey()
cv2.destroyAllWindows()

实验四

特征点提取算法

import cv2
import numpy as np
img = cv2.imread('5.jpg')
fast = cv2.FastFeatureDetector_create()
keypoints = fast.detect(img, None)
img_fast = cv2.drawKeypoints(img, keypoints, None, (0, 0, 255))
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
alg = cv2.xfeatures2d.SURF_create(10000)
keypoints, descriptor = alg.detectAndCompute(gray, None)
cv2.imshow("img", img)
img_surf = cv2.drawKeypoints(img, keypoints, img, (0, 0, 255), cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)

cv2.imshow('img_surf', img_surf)
cv2.imshow("img_fast", img_fast)
cv2.waitKey()
cv2.destroyAllWindows()

形态学验证实验

import cv2
import numpy as np
img = cv2.imread('5.jpg')
kernel = np.ones((5, 5), np.uint8)
erosion = cv2.erode(img, kernel)
dilation = cv2.dilate(img, kernel)
opening = cv2.morphologyEx(img, cv2.MORPH_OPEN, kernel)
closing = cv2.morphologyEx(img, cv2.MORPH_CLOSE, kernel)

cv2.imshow('img', img)
cv2.imshow("erosion", erosion)
cv2.imshow("dilation", dilation)
cv2.imshow("opening", opening)
cv2.imshow("closing", closing)
cv2.waitKey()
cv2.destroyAllWindows()

低通和高通滤波两种频域增强

import cv2
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import math
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False    # 用来正常显示负号

 
def low_pass_filtering(image, radius):
    """
    低通滤波函数
    :param image: 输入图像
    :param radius: 半径
    :return: 滤波结果
    """
    # 对图像进行傅里叶变换，fft是一个三维数组，fft[:, :, 0]为实数部分，fft[:, :, 1]为虚数部分
    fft = cv2.dft(np.float32(image), flags=cv2.DFT_COMPLEX_OUTPUT)
    # 对fft进行中心化，生成的dshift仍然是一个三维数组
    dshift = np.fft.fftshift(fft)
 
    # 得到中心像素
    rows, cols = image.shape[:2]
    mid_row, mid_col = int(rows / 2), int(cols / 2)
 
    # 构建掩模，256位，两个通道
    mask = np.zeros((rows, cols, 2), np.float32)
    mask[mid_row - radius:mid_row + radius, mid_col - radius:mid_col + radius] = 1
 
    # 给傅里叶变换结果乘掩模
    fft_filtering = dshift * mask
    # 傅里叶逆变换
    ishift = np.fft.ifftshift(fft_filtering)
    image_filtering = cv2.idft(ishift)
    image_filtering = cv2.magnitude(image_filtering[:, :, 0], image_filtering[:, :, 1])
    # 对逆变换结果进行归一化（一般对图像处理的最后一步都要进行归一化，特殊情况除外）
    cv2.normalize(image_filtering, image_filtering, 0, 1, cv2.NORM_MINMAX)
    return image_filtering
 
 
def high_pass_filtering(image, radius, n):
    """
    高通滤波函数
    :param image: 输入图像
    :param radius: 半径
    :param n: ButterWorth滤波器阶数
    :return: 滤波结果
    """
    # 对图像进行傅里叶变换，fft是一个三维数组，fft[:, :, 0]为实数部分，fft[:, :, 1]为虚数部分
    fft = cv2.dft(np.float32(image), flags=cv2.DFT_COMPLEX_OUTPUT)
    # 对fft进行中心化，生成的dshift仍然是一个三维数组
    dshift = np.fft.fftshift(fft)
 
    # 得到中心像素
    rows, cols = image.shape[:2]
    mid_row, mid_col = int(rows / 2), int(cols / 2)
 
    # 构建ButterWorth高通滤波掩模
 
    mask = np.zeros((rows, cols, 2), np.float32)
    for i in range(0, rows):
        for j in range(0, cols):
            # 计算(i, j)到中心点的距离
            d = math.sqrt(pow(i - mid_row, 2) + pow(j - mid_col, 2))
            try:
                mask[i, j, 0] = mask[i, j, 1] = 1 / (1 + pow(radius / d, 2*n))
            except ZeroDivisionError:
                mask[i, j, 0] = mask[i, j, 1] = 0
    # 给傅里叶变换结果乘掩模
    fft_filtering = dshift * mask
    # 傅里叶逆变换
    ishift = np.fft.ifftshift(fft_filtering)
    image_filtering = cv2.idft(ishift)
    image_filtering = cv2.magnitude(image_filtering[:, :, 0], image_filtering[:, :, 1])
    # 对逆变换结果进行归一化（一般对图像处理的最后一步都要进行归一化，特殊情况除外）
    cv2.normalize(image_filtering, image_filtering, 0, 1, cv2.NORM_MINMAX)
    return image_filtering
 
 
def bandpass_filter(image, radius, w, n=1):
    """
    带通滤波函数
    :param image: 输入图像
    :param radius: 带中心到频率平面原点的距离
    :param w: 带宽
    :param n: 阶数
    :return: 滤波结果
    """
    # 对图像进行傅里叶变换，fft是一个三维数组，fft[:, :, 0]为实数部分，fft[:, :, 1]为虚数部分
    fft = cv2.dft(np.float32(image), flags=cv2.DFT_COMPLEX_OUTPUT)
    # 对fft进行中心化，生成的dshift仍然是一个三维数组
    dshift = np.fft.fftshift(fft)
 
    # 得到中心像素
    rows, cols = image.shape[:2]
    mid_row, mid_col = int(rows / 2), int(cols / 2)
 
    # 构建掩模，256位，两个通道
    mask = np.zeros((rows, cols, 2), np.float32)
    for i in range(0, rows):
        for j in range(0, cols):
            # 计算(i, j)到中心点的距离
            d = math.sqrt(pow(i - mid_row, 2) + pow(j - mid_col, 2))
            if radius - w / 2 < d < radius + w / 2:
                mask[i, j, 0] = mask[i, j, 1] = 1
            else:
                mask[i, j, 0] = mask[i, j, 1] = 0
 
    # 给傅里叶变换结果乘掩模
    fft_filtering = dshift * np.float32(mask)
    # 傅里叶逆变换
    ishift = np.fft.ifftshift(fft_filtering)
    image_filtering = cv2.idft(ishift)
    image_filtering = cv2.magnitude(image_filtering[:, :, 0], image_filtering[:, :, 1])
    # 对逆变换结果进行归一化（一般对图像处理的最后一步都要进行归一化，特殊情况除外）
    cv2.normalize(image_filtering, image_filtering, 0, 1, cv2.NORM_MINMAX)
    return image_filtering
 
 
def bandstop_filter(image, radius, w, n=1):
    """
    带通滤波函数
    :param image: 输入图像
    :param radius: 带中心到频率平面原点的距离
    :param w: 带宽
    :param n: 阶数
    :return: 滤波结果
    """
    # 对图像进行傅里叶变换，fft是一个三维数组，fft[:, :, 0]为实数部分，fft[:, :, 1]为虚数部分
    fft = cv2.dft(np.float32(image), flags=cv2.DFT_COMPLEX_OUTPUT)
    # 对fft进行中心化，生成的dshift仍然是一个三维数组
    dshift = np.fft.fftshift(fft)
 
    # 得到中心像素
    rows, cols = image.shape[:2]
    mid_row, mid_col = int(rows / 2), int(cols / 2)
 
    # 构建掩模，256位，两个通道
    mask = np.zeros((rows, cols, 2), np.float32)
    for i in range(0, rows):
        for j in range(0, cols):
            # 计算(i, j)到中心点的距离
            d = math.sqrt(pow(i - mid_row, 2) + pow(j - mid_col, 2))
            if radius - w / 2 < d < radius + w / 2:
                mask[i, j, 0] = mask[i, j, 1] = 0
            else:
                mask[i, j, 0] = mask[i, j, 1] = 1
 
    # 给傅里叶变换结果乘掩模
    fft_filtering = dshift * np.float32(mask)
    # 傅里叶逆变换
    ishift = np.fft.ifftshift(fft_filtering)
    image_filtering = cv2.idft(ishift)
    image_filtering = cv2.magnitude(image_filtering[:, :, 0], image_filtering[:, :, 1])
    # 对逆变换结果进行归一化（一般对图像处理的最后一步都要进行归一化，特殊情况除外）
    cv2.normalize(image_filtering, image_filtering, 0, 1, cv2.NORM_MINMAX)
    return image_filtering
 
 
if __name__ == "__main__":
    image = cv2.imread("4.jpg", 0)
    
    image_low_pass_filtering1 = low_pass_filtering(image, 10)

    image_high_pass_filtering1 = high_pass_filtering(image, 10, 1)

    plt.subplot(221), plt.imshow(image, 'gray'), plt.title("原图"), plt.xticks([]), plt.yticks([])
    plt.subplot(222), plt.imshow(image_low_pass_filtering1, 'gray'), plt.title("半径为10像素的低通滤波"), plt.xticks([]), plt.yticks([])
    plt.subplot(223), plt.imshow(image_high_pass_filtering1, 'gray'), plt.title("半径为10像素的高通滤波"), plt.xticks([]), plt.yticks([])

    plt.show()
    cv2.imshow('img', image)
    cv2.imshow("img_low", image_low_pass_filtering1)
    cv2.imshow("img_high", image_high_pass_filtering1)

    cv2.waitKey()
    cv2.destroyAllWindows()