(十)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:\Users\silence_cho\Desktop\Messi.jpg"
    # img_file = r"D:\data\receipt_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:\Users\silence_cho\Desktop\Messi.jpg"
    # img_file = r"D:\data\receipt_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:\Users\silence_cho\Desktop\Messi.jpg"
    # img_file = r"D:\data\receipt_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:\Users\silence_cho\Desktop\Messi.jpg"
    # img_file = r"D:\data\receipt_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()