傅里叶图像相关性匹配-《医学图像处理》小作业五-Python代码/Matlab代码

天津中医药大学-20级医学信息工程   教师：王翌   学生：邓集亲

学长我是用的python写的，matlab同样可以参考

实验五：相关性匹配

作业要求：

1. 参考“傅里叶变换”课的内容，采用快速傅里叶变换(FFT)进行相关性匹配，如下图示例输出结果图片。

例子：

通过相关性匹配找到感兴趣的物体区域：

实验图片：

匹配图片：match.jpg

模板图片：match_window.jpg

再自选其它至少20张图片（包括灰度图片和彩色图片）。

  二维傅里叶变换用于分析图像的频率特性，对采样的离散信号进行傅里叶变换，得到周期性的N点DFT。可以将图像视为在两个方向上采样的信号。因此，在X和Y方向都进行傅立叶变换，可以得到图像的频率表示。

  更直观地说，对于正弦信号，如果幅度在短时间内变化很快，则可以说它是高频信号。如果变化缓慢，则为低频信号。

  将相同的想法扩展到图像。图像中的振幅在哪里急剧变化？在边缘点或噪声。因此，可以说边缘和噪声是图像中的高频内容。如果幅度没有太大变化，则它是低频分量。

  模板匹配的原理就是利用目标的边缘信息用于搜索目标图像的模板所在位置，然后将模板图像在输入图像上滑动匹配，找到最大匹配度，模板匹配就完成了。

#导入库
import numpy as np
import cv2 as cv
from numpy.core import asarray, zeros, swapaxes, take
from numpy.core.multiarray import normalize_axis_index
from numpy.fft import _pocketfft_internal as pfi

#定义函数进行fft变换计算
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, shape=None, axes=None, invreal=0):
    if shape is None:
        shapeless = 1
        if axes is None:
            shape = list(a.shape)
        else:
            shape = take(a.shape, axes)
    else:
        shapeless = 0
    shape = list(shape)
    if axes is None:
        axes = list(range(-len(shape), 0))
    if invreal and shapeless:
        shape[-1] = (a.shape[axes[-1]] - 1) * 2
    return shape, axes


def raw_fftnd(a, shape=None, axes=None, function=fft):
    a = asarray(a)
    shape, axes = cook_nd_args(a, shape, axes)
    itl = list(range(len(axes)))
    itl.reverse()
    for ii in itl:
        a = function(a, n=shape[ii], axis=axes[ii])
    return a
def imagefft2(a, shape=None, axes=(-2, -1)):   #图片二维傅里叶变换
    return raw_fftnd(a, shape, axes, fft)


def imageifft2(a, shape=None, axes=(-2, -1)):
    return raw_fftnd(a, shape, axes, ifft)

#读取图像

original = cv.imread(r"D:\deng\smalljob\ImageSet\match.jpg")   #待匹配的原图像
template = cv.imread(r"D:\deng\smalljob\ImageSet\match_window.jpg")   #需要匹配的部分

original_gray = cv.cvtColor(original, cv.COLOR_BGR2GRAY)   #将图像从RGB颜色空间转换到灰度颜色空间
template_gray = cv.cvtColor(template, cv.COLOR_BGR2GRAY)

t_height, t_width = template_gray.shape
o_height, o_width = original_gray.shape
#傅里叶变换

t_fft = imagefft2(template_gray, shape=(o_height, o_width))
o_fft = imagefft2(original_gray)

c = imageifft2(np.multiply(o_fft, t_fft.conj()) / np.abs(np.multiply(o_fft, t_fft.conj())))
c = c.real

result = np.where(c == np.amax(c))

#压缩 2 个数组以获得确切的坐标
max_coordinate = list(zip(result[0], result[1]))[0]
print(max_coordinate)

start_point = (max_coordinate[1], max_coordinate[0] )
end_point = (max_coordinate[1] + t_height, max_coordinate[0] + t_width)

#用绿色矩形方框囊括匹配结果

color = (0, 255, 0)
thickness = 2
image = cv.rectangle(original, start_point, end_point, color, thickness)
cv.imshow("fft matching", image)
cv.waitKey(0)
cv.destroyAllWindows()

运行结果