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

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

声明：本文章中所涉及的代码并非我个人独立完成的成果。

        在撰写的过程中，我广泛地吸取了前一辈人，尤其是学长学姐们的宝贵学习经验。通过深入研究他们的学习轨迹，以及查阅和分析了众多相关的理论文献与资料，并在王老师的悉心指导下，经过反复的实验、调试与优化，最终得以总结完成本文所展现的代码。

        建议各位学弟学妹们，先根据王老师的授课内容，独立思考一下应该如何完成。如果实在是有理解上的困难，不知道从何下手，再来学习和参考本文。

        学长我是用Python写的，如果你使用MATLAB，同样可以参考此代码进行翻译。我在代码中加入了许多注释和测试环节，以便于理解。由于学长能力有限，代码中或许存在一些疏漏或错误，请批判性地参考。

实验五：相关性匹配

作业要求：

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

例子：

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

实验图片：

匹配图片：match.jpg

模板图片：match_window.jpg

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

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

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

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

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

Windows图片文件目录：

ImageSet文件目录下储存其他待处理图片

Python代码：

#导入库
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()

运行结果：