(十)OpenCV-Python学习—频率域滤波

　　之前的图像处理，都是再原图上进行；而频率域滤波，是在图像的傅里叶谱上进行处理，最后再通过傅里叶逆变换得到处理后的图像，则是因为图片的傅里叶谱包含图片的频率信息，方便对其频率进行处理。对于图像，低频信息表示图像中灰度值缓慢变化的区域，如背景信息等；而高频信息则表示灰度值迅速变化的区域，如边缘处等细节信息。

　　在经过中心化后的傅立叶谱(幅度谱)，其中心位置的幅度值最大，频率最低，随着离中心位置的距离增加频率会越来越大，所以，中心化后的傅里叶谱，中心位置为低频区域，四个角落处为高频区域。频率域滤波通常的处理步骤如下：

 　　常用的滤波器有四种：低通滤波器，高通滤波器，带通滤波器，带阻滤波器

低通滤波器

　　低通滤波器，即保留傅里叶变换的低频信息，过滤掉高频信息，会使图片变得更模糊。常用的低通滤波器包括理想低通滤波器，巴特沃斯低通滤波器，高斯低通滤波器。假设图像傅里叶变换的高，宽为H、W，傅里叶谱的最大值在中心点位置(maxR, maxC)， D(r, c)代表点(r, c)到中心点的距离:

　　那么三种滤波器可以表示为：

　　理想低通滤波器：

 　　巴特沃斯低通滤波器：

 　　高斯低通滤波器：

 　　低通滤波器的使用代码及结果如下：

#coding:utf-8

import cv2
import numpy as np


def createLPFilter(shape, center, radius, lpType=2, n=2):
    rows, cols = shape[:2]
    r, c = np.mgrid[0:rows:1, 0:cols:1]
    c -= center[0]
    r -= center[1]
    d = np.power(c, 2.0) + np.power(r, 2.0)
    lpFilter_matrix = np.zeros(shape, np.float32)
    if lpType == 0:  # 理想低通滤波器
        lpFilter = np.copy(d)
        lpFilter[lpFilter < pow(radius, 2.0)] = 1
        lpFilter[lpFilter >= pow(radius, 2.0)] = 0
    elif lpType == 1: #巴特沃斯低通滤波器
        lpFilter = 1.0 / (1 + np.power(np.sqrt(d)/radius, 2*n))
    elif lpType == 2: # 高斯低通滤波器
        lpFilter = np.exp(-d/(2*pow(radius, 2.0)))
    lpFilter_matrix[:, :, 0] = lpFilter
    lpFilter_matrix[:, :, 1] = lpFilter
    return lpFilter_matrix

def stdFftImage(img_gray, rows, cols):
    fimg = np.copy(img_gray)
    fimg = fimg.astype(np.float32)   #注意这里的类型转换
    # 1.图像矩阵乘以（-1）^(r+c), 中心化
    for r in range(rows):
        for c in range(cols):
            if (r+c) % 2:
                fimg[r][c] = -1 * img_gray[r][c]
    img_fft = fftImage(fimg, rows, cols)
    return img_fft


def fftImage(img_gray, rows, cols):
    rPadded = cv2.getOptimalDFTSize(rows)
    cPadded = cv2.getOptimalDFTSize(cols)
    imgPadded = np.zeros((rPadded, cPadded), dtype=np.float32)
    imgPadded[:rows, :cols] = img_gray
    img_fft = cv2.dft(imgPadded, flags=cv2.DFT_COMPLEX_OUTPUT)
    return img_fft


def graySpectrum(fft_img):
    real = np.power(fft_img[:, :, 0], 2.0)
    imaginary = np.power(fft_img[:, :, 1], 2.0)
    amplitude = np.sqrt(real+imaginary)
    spectrum = np.log(amplitude+1.0)
    spectrum = cv2.normalize(spectrum, 0, 1, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_32F)
    spectrum *= 255
    return amplitude, spectrum


def nothing(args):
    pass

if __name__ == "__main__":
    img_file = r"C:Userssilence_choDesktopMessi.jpg"
    # img_file = r"D:data
eceipt_rotate.jpg"
    img_gray = cv2.imread(img_file, 0)
    # 1.快速傅里叶变换
    rows, cols = img_gray.shape[:2]
    img_fft = stdFftImage(img_gray, rows, cols)
    amplitude, _ = graySpectrum(img_fft)
    minValue, maxValue, minLoc, maxLoc = cv2.minMaxLoc(amplitude)  #中心化后频谱的最大值在图片中心位置处

    cv2.namedWindow("tracks")
    max_radius = np.sqrt(pow(rows, 2) + pow(cols, 2))/2
    cv2.createTrackbar("Radius", "tracks", 0, int(max_radius), nothing)
    cv2.createTrackbar("Filter type", "tracks", 0, 2, nothing)

    while True:
        # 2.构建低通滤波器
        radius = cv2.getTrackbarPos("Radius", "tracks")
        lpType = cv2.getTrackbarPos("Filter type", "tracks")
        nrows, ncols = img_fft.shape[:2]

        # x, y = int(ncols/2), int(nrows/2)  # 注意这里是坐标
        # ilpFilter = createLPFilter(img_fft.shape, (x, y), radius, lpType)
        ilpFilter = createLPFilter(img_fft.shape, maxLoc, radius, lpType)

        # 3.低通滤波器滤波
        img_filter = ilpFilter*img_fft

        _, gray_spectrum = graySpectrum(img_filter)  #观察滤波器的变化

        # 4. 傅里叶反变换，并取实部进行裁剪, 并去中心化
        img_ift = cv2.dft(img_filter, flags=cv2.DFT_INVERSE+cv2.DFT_REAL_OUTPUT+cv2.DFT_SCALE)
        ori_img = np.copy(img_ift[:rows, :cols])
        for r in range(rows):
            for c in range(cols):
                if(r+c)%2:
                    ori_img[r][c] = -1*ori_img[r][c]
                # 截断高低值
                if ori_img[r][c] < 0:
                    ori_img[r][c] = 0
                if ori_img[r][c] > 255:
                    ori_img[r][c] = 255
        # ori_img[ori_img < 0] = 0
        # ori_img[ori_img > 255] = 255
        ori_img = ori_img.astype(np.uint8)

        cv2.imshow("img_gray", img_gray)
        cv2.imshow("ori_img", ori_img)
        cv2.imshow("gray_spectrum", gray_spectrum)
        key = cv2.waitKey(1)
        if key == 27:
            break
    cv2.destroyAllWindows()

 　　因为更多的高频信息被过滤掉了，从上图也可以发现，低通滤波器对图片起到了模糊作用。

高通滤波器

　　高通滤波器保留图片的高频信息，过滤掉低频信息，会使图片的细节部分更加明显。常用的高通滤波器包括理想高通滤波器，巴特沃斯高通滤波器，高斯高通滤波器。三种高通滤波器可以表示为：

　　理性高通滤波器：

 　　巴特沃斯高通滤波器：

 　　高斯高通滤波器：

 　　高通滤波器的使用代码及结果如下：

#coding:utf-8

import cv2
import numpy as np


def createHPFilter(shape, center, radius, lpType=2, n=2):
    rows, cols = shape[:2]
    r, c = np.mgrid[0:rows:1, 0:cols:1]
    c -= center[0]
    r -= center[1]
    d = np.power(c, 2.0) + np.power(r, 2.0)
    lpFilter_matrix = np.zeros(shape, np.float32)
    if lpType == 0:  # 理想高通滤波器
        lpFilter = np.copy(d)
        lpFilter[lpFilter < pow(radius, 2.0)] = 0
        lpFilter[lpFilter >= pow(radius, 2.0)] = 1
    elif lpType == 1: #巴特沃斯高通滤波器
        lpFilter = 1.0 - 1.0 / (1 + np.power(np.sqrt(d)/radius, 2*n))
    elif lpType == 2: # 高斯高通滤波器
        lpFilter = 1.0 - np.exp(-d/(2*pow(radius, 2.0)))
    lpFilter_matrix[:, :, 0] = lpFilter
    lpFilter_matrix[:, :, 1] = lpFilter
    return lpFilter_matrix

def stdFftImage(img_gray, rows, cols):
    fimg = np.copy(img_gray)
    fimg = fimg.astype(np.float32)   #注意这里的类型转换
    # 1.图像矩阵乘以（-1）^(r+c), 中心化
    for r in range(rows):
        for c in range(cols):
            if (r+c) % 2:
                fimg[r][c] = -1 * img_gray[r][c]
    img_fft = fftImage(fimg, rows, cols)
    return img_fft

def fftImage(img_gray, rows, cols):
    rPadded = cv2.getOptimalDFTSize(rows)
    cPadded = cv2.getOptimalDFTSize(cols)
    imgPadded = np.zeros((rPadded, cPadded), dtype=np.float32)
    imgPadded[:rows, :cols] = img_gray
    img_fft = cv2.dft(imgPadded, flags=cv2.DFT_COMPLEX_OUTPUT)
    return img_fft

def graySpectrum(fft_img):
    real = np.power(fft_img[:, :, 0], 2.0)
    imaginary = np.power(fft_img[:, :, 1], 2.0)
    amplitude = np.sqrt(real+imaginary)
    spectrum = np.log(amplitude+1.0)
    spectrum = cv2.normalize(spectrum, 0, 1, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_32F)
    spectrum *= 255
    return amplitude, spectrum

def nothing(args):
    pass

if __name__ == "__main__":
    img_file = r"C:Userssilence_choDesktopMessi.jpg"
    # img_file = r"D:data
eceipt_rotate.jpg"
    img_gray = cv2.imread(img_file, 0)
    # 1.快速傅里叶变换
    rows, cols = img_gray.shape[:2]
    img_fft = stdFftImage(img_gray, rows, cols)
    amplitude, _ = graySpectrum(img_fft)
    minValue, maxValue, minLoc, maxLoc = cv2.minMaxLoc(amplitude)  # 中心化后频谱的最大值在图片中心位置处

    cv2.namedWindow("tracks")
    max_radius = np.sqrt(pow(rows, 2) + pow(cols, 2))
    cv2.createTrackbar("Radius", "tracks", 0, int(max_radius), nothing)
    cv2.createTrackbar("Filter type", "tracks", 0, 2, nothing)

    while True:
        # 2.构建高通滤波器
        radius = cv2.getTrackbarPos("Radius", "tracks")
        lpType = cv2.getTrackbarPos("Filter type", "tracks")
        nrows, ncols = img_fft.shape[:2]
        # x, y = int(ncols / 2), int(nrows / 2)  # 注意这里是坐标
        # ilpFilter = createHPFilter(img_fft.shape, (x, y), radius, lpType)
        ilpFilter = createHPFilter(img_fft.shape, maxLoc, radius, lpType)

        # 3.高通滤波器滤波
        img_filter = ilpFilter * img_fft

        _, gray_spectrum = graySpectrum(img_filter)  # 观察滤波器的变化

        # 4. 傅里叶反变换，并取实部进行裁剪, 并去中心化
        img_ift = cv2.dft(img_filter, flags=cv2.DFT_INVERSE + cv2.DFT_REAL_OUTPUT + cv2.DFT_SCALE)
        ori_img = np.copy(img_ift[:rows, :cols])
        for r in range(rows):
            for c in range(cols):
                if (r + c) % 2:
                    ori_img[r][c] = -1 * ori_img[r][c]
                # 截断高低值
                if ori_img[r][c] < 0:
                    ori_img[r][c] = 0
                if ori_img[r][c] > 255:
                    ori_img[r][c] = 255
        # ori_img[ori_img < 0] = 0
        # ori_img[ori_img > 255] = 255
        ori_img = ori_img.astype(np.uint8)

        cv2.imshow("img_gray", img_gray)
        cv2.imshow("ori_img", ori_img)
        cv2.imshow("gray_spectrum", gray_spectrum)
        key = cv2.waitKey(1)
        if key == 27:
            break
    cv2.destroyAllWindows()

　　因为高通滤波器过滤掉了低频信息，从上图发现，高通滤波器对图片起到了锐化的作用，仅保留了图片中物体边缘信息。

带通滤波器

　　带通滤波器是只保留某一范围区域的频率带，频率信息过滤掉，能选择性的图片的部分信息。常用的带通滤波器包括理想带通滤波器，巴特沃斯带通滤波器，高斯带通滤波器。假设BW代表带宽，D0代表带宽的径向中心，则三种带通滤波器可以表示为：

　　理想带通滤波器：

 　　巴特沃斯带通滤波器：

 　　高斯带通滤波器：

　　 带通滤波器的使用代码及效果如下：

#coding:utf-8

import cv2
import numpy as np


def createBPFilter(shape, center, bandCenter, bandWidth, lpType=2, n=2):
    rows, cols = shape[:2]
    r, c = np.mgrid[0:rows:1, 0:cols:1]
    c -= center[0]
    r -= center[1]
    d = np.sqrt(np.power(c, 2.0) + np.power(r, 2.0))
    lpFilter_matrix = np.zeros(shape, np.float32)
    if lpType == 0:  # 理想带通滤波器
        lpFilter = np.copy(d)
        lpFilter[:, :] = 1
        lpFilter[d > (bandCenter+bandWidth/2)] = 0
        lpFilter[d < (bandCenter-bandWidth/2)] = 0
    elif lpType == 1: #巴特沃斯带通滤波器
        lpFilter = 1.0 - 1.0 / (1 + np.power(d*bandWidth/(d - pow(bandCenter,2)), 2*n))
    elif lpType == 2: # 高斯带通滤波器
        lpFilter = np.exp(-pow((d-pow(bandCenter,2))/(d*bandWidth), 2))
    lpFilter_matrix[:, :, 0] = lpFilter
    lpFilter_matrix[:, :, 1] = lpFilter
    return lpFilter_matrix

def stdFftImage(img_gray, rows, cols):
    fimg = np.copy(img_gray)
    fimg = fimg.astype(np.float32)   #注意这里的类型转换
    # 1.图像矩阵乘以（-1）^(r+c), 中心化
    for r in range(rows):
        for c in range(cols):
            if (r+c) % 2:
                fimg[r][c] = -1 * img_gray[r][c]
    img_fft = fftImage(fimg, rows, cols)
    return img_fft

def fftImage(img_gray, rows, cols):
    rPadded = cv2.getOptimalDFTSize(rows)
    cPadded = cv2.getOptimalDFTSize(cols)
    imgPadded = np.zeros((rPadded, cPadded), dtype=np.float32)
    imgPadded[:rows, :cols] = img_gray
    img_fft = cv2.dft(imgPadded, flags=cv2.DFT_COMPLEX_OUTPUT)
    return img_fft

def graySpectrum(fft_img):
    real = np.power(fft_img[:, :, 0], 2.0)
    imaginary = np.power(fft_img[:, :, 1], 2.0)
    amplitude = np.sqrt(real+imaginary)
    spectrum = np.log(amplitude+1.0)
    spectrum = cv2.normalize(spectrum, 0, 1, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_32F)
    spectrum *= 255
    return amplitude, spectrum

def nothing(args):
    pass

if __name__ == "__main__":
    img_file = r"C:Userssilence_choDesktopMessi.jpg"
    # img_file = r"D:data
eceipt_rotate.jpg"
    img_gray = cv2.imread(img_file, 0)
    # 1.快速傅里叶变换
    rows, cols = img_gray.shape[:2]
    img_fft = stdFftImage(img_gray, rows, cols)
    amplitude, _ = graySpectrum(img_fft)
    minValue, maxValue, minLoc, maxLoc = cv2.minMaxLoc(amplitude)  # 中心化后频谱的最大值在图片中心位置处

    cv2.namedWindow("tracks")
    max_radius = np.sqrt(pow(rows, 2) + pow(cols, 2))
    cv2.createTrackbar("BandCenter", "tracks", 0, int(max_radius), nothing)
    cv2.createTrackbar("BandWidth", "tracks", 0, int(max_radius), nothing)
    cv2.createTrackbar("Filter type", "tracks", 0, 2, nothing)

    while True:
        # 2.构建带通滤波器
        bandCenter = cv2.getTrackbarPos("BandCenter", "tracks")
        bandWidth = cv2.getTrackbarPos("BandWidth", "tracks")
        lpType = cv2.getTrackbarPos("Filter type", "tracks")
        nrows, ncols = img_fft.shape[:2]
        # x, y = int(ncols / 2), int(nrows / 2)  # 注意这里是坐标
        # ilpFilter = createBPFilter(img_fft.shape, (x, y), bandCenter, bandWidth, lpType)
        ilpFilter = createBPFilter(img_fft.shape, maxLoc, bandCenter, bandWidth, lpType)

        # 3.带通滤波器滤波
        img_filter = ilpFilter * img_fft

        _, gray_spectrum = graySpectrum(img_filter)  # 观察滤波器的变化

        # 4. 傅里叶反变换，并取实部进行裁剪, 并去中心化
        img_ift = cv2.dft(img_filter, flags=cv2.DFT_INVERSE + cv2.DFT_REAL_OUTPUT + cv2.DFT_SCALE)
        ori_img = np.copy(img_ift[:rows, :cols])
        for r in range(rows):
            for c in range(cols):
                if (r + c) % 2:
                    ori_img[r][c] = -1 * ori_img[r][c]
                # 截断高低值
                if ori_img[r][c] < 0:
                    ori_img[r][c] = 0
                if ori_img[r][c] > 255:
                    ori_img[r][c] = 255
        # ori_img[ori_img < 0] = 0
        # ori_img[ori_img > 255] = 255
        ori_img = ori_img.astype(np.uint8)

        cv2.imshow("img_gray", img_gray)
        cv2.imshow("ori_img", ori_img)
        cv2.imshow("gray_spectrum", gray_spectrum)
        key = cv2.waitKey(1)
        if key == 27:
            break
    cv2.destroyAllWindows()

带阻滤波器

　　与带通滤波器相反，带阻滤波器指过滤掉或者削弱指定范围区域的频率带，常用的带阻滤波器包括理想带阻滤波器，巴特沃斯带阻滤波器，高斯带阻滤波器。三种带阻滤波器表示如下：

　　理想带阻滤波器：

 　　巴特沃斯带阻滤波器：

　　高斯带阻滤波器：

　　带阻滤波器使用代码及效果如下：

#coding:utf-8

import cv2
import numpy as np


def createBRFilter(shape, center, bandCenter, bandWidth, lpType=2, n=2):
    rows, cols = shape[:2]
    r, c = np.mgrid[0:rows:1, 0:cols:1]
    c -= center[0]
    r -= center[1]
    d = np.sqrt(np.power(c, 2.0) + np.power(r, 2.0))
    lpFilter_matrix = np.zeros(shape, np.float32)
    if lpType == 0:  # 理想带阻滤波器
        lpFilter = np.copy(d)
        lpFilter[:, :] = 0
        lpFilter[d > (bandCenter+bandWidth/2)] = 1
        lpFilter[d < (bandCenter-bandWidth/2)] = 1
    elif lpType == 1: #巴特沃斯带阻滤波器
        lpFilter = 1.0 / (1 + np.power(d*bandWidth/(d - pow(bandCenter,2)), 2*n))
    elif lpType == 2: # 高斯带阻滤波器
        lpFilter = 1 - np.exp(-pow((d-pow(bandCenter,2))/(d*bandWidth), 2))
    lpFilter_matrix[:, :, 0] = lpFilter
    lpFilter_matrix[:, :, 1] = lpFilter
    return lpFilter_matrix

def stdFftImage(img_gray, rows, cols):
    fimg = np.copy(img_gray)
    fimg = fimg.astype(np.float32)   #注意这里的类型转换
    # 1.图像矩阵乘以（-1）^(r+c), 中心化
    for r in range(rows):
        for c in range(cols):
            if (r+c) % 2:
                fimg[r][c] = -1 * img_gray[r][c]
    img_fft = fftImage(fimg, rows, cols)
    return img_fft

def fftImage(img_gray, rows, cols):
    rPadded = cv2.getOptimalDFTSize(rows)
    cPadded = cv2.getOptimalDFTSize(cols)
    imgPadded = np.zeros((rPadded, cPadded), dtype=np.float32)
    imgPadded[:rows, :cols] = img_gray
    img_fft = cv2.dft(imgPadded, flags=cv2.DFT_COMPLEX_OUTPUT)
    return img_fft

def graySpectrum(fft_img):
    real = np.power(fft_img[:, :, 0], 2.0)
    imaginary = np.power(fft_img[:, :, 1], 2.0)
    amplitude = np.sqrt(real+imaginary)
    spectrum = np.log(amplitude+1.0)
    spectrum = cv2.normalize(spectrum, 0, 1, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_32F)
    spectrum *= 255
    return amplitude, spectrum

def nothing(args):
    pass

if __name__ == "__main__":
    img_file = r"C:Userssilence_choDesktopMessi.jpg"
    # img_file = r"D:data
eceipt_rotate.jpg"
    img_gray = cv2.imread(img_file, 0)
    # 1.快速傅里叶变换
    rows, cols = img_gray.shape[:2]
    img_fft = stdFftImage(img_gray, rows, cols)
    amplitude, _ = graySpectrum(img_fft)
    minValue, maxValue, minLoc, maxLoc = cv2.minMaxLoc(amplitude)  # 中心化后频谱的最大值在图片中心位置处

    cv2.namedWindow("tracks")
    max_radius = np.sqrt(pow(rows, 2) + pow(cols, 2))
    cv2.createTrackbar("BandCenter", "tracks", 0, int(max_radius), nothing)
    cv2.createTrackbar("BandWidth", "tracks", 0, int(max_radius), nothing)
    cv2.createTrackbar("Filter type", "tracks", 0, 2, nothing)

    while True:
        # 2.构建带阻滤波器
        bandCenter = cv2.getTrackbarPos("BandCenter", "tracks")
        bandWidth = cv2.getTrackbarPos("BandWidth", "tracks")
        lpType = cv2.getTrackbarPos("Filter type", "tracks")
        nrows, ncols = img_fft.shape[:2]
        # x, y = int(ncols / 2), int(nrows / 2)  # 注意这里是坐标
        # ilpFilter = createBRFilter(img_fft.shape, (x, y), bandCenter, bandWidth, lpType)
        ilpFilter = createBRFilter(img_fft.shape, maxLoc, bandCenter, bandWidth, lpType)

        # 3.带阻滤波器滤波
        img_filter = ilpFilter * img_fft

        _, gray_spectrum = graySpectrum(img_filter)  # 观察滤波器的变化

        # 4. 傅里叶反变换，并取实部进行裁剪, 并去中心化
        img_ift = cv2.dft(img_filter, flags=cv2.DFT_INVERSE + cv2.DFT_REAL_OUTPUT + cv2.DFT_SCALE)
        ori_img = np.copy(img_ift[:rows, :cols])
        for r in range(rows):
            for c in range(cols):
                if (r + c) % 2:
                    ori_img[r][c] = -1 * ori_img[r][c]
                    # 截断高低值
                if ori_img[r][c] < 0:
                    ori_img[r][c] = 0
                if ori_img[r][c] > 255:
                    ori_img[r][c] = 255
        # ori_img[ori_img < 0] = 0
        # ori_img[ori_img > 255] = 255
        ori_img = ori_img.astype(np.uint8)

        cv2.imshow("img_gray", img_gray)
        cv2.imshow("ori_img", ori_img)
        cv2.imshow("gray_spectrum", gray_spectrum)
        key = cv2.waitKey(1)
        if key == 27:
            break
    cv2.destroyAllWindows()