# 卷积

# 一维卷积

连续域函数卷积

$g(x)=f(x)*k(x)=\int_{-\infty}^{+\infty}{f(t)·k(x-t)·dt}$

离散函数卷积

$g(i)=f(i)*k(i)=\sum\limits_{l}f(l)·k(i-l)$

# 二维卷积

$G(x,y)=I(x,y)*K(x,y)=\sum\limits_{i=1}^{H}\sum\limits_{j=1}^{W}I(i,j)·K(x-i,y-j)$

# 两种信号

# 冲激信号

a = [0, 1, 2, 3, 1, 0] # 输入信号
k = [0, 1, 0] # 卷积核

# 输入信号与冲激信号卷积后，输入信号在形状上不会发生变化，在位置上可能存在唯一

k_size = 15
k1 = np.zeros((k_size, k_size))
mid = k_size // 2
k1[mid][mid] = 1

# 输入信号（图像）与 k1 进行二维卷积后，图像不会发生变化

# 方波信号

a = [0, 1, 2, 3, 2, 1, 0] # 输入信号
k = [1, 1, 1, 1, 1] # 卷积核

# 输入信号与方波信号卷积后，输入信号的形状会变得平滑

k_size = 15
k1 = 1 / k_size / k_size * np.ones((k_size, k_size))
# 输入信号（图像）与 k1 进行二维卷积后，图像会变得平滑、模糊

# 图像噪声

# 椒盐噪声

像撒了盐和胡椒一样

def add_Salt(img, percent):
    # 生成噪声
    noise = np.random.uniform(0, 1, img[:, :, 0].shape)   # 获取图像的高和宽，不包括颜色通道

    # 添加盐噪声（白色 = [255, 255, 255]）
    salt_mask = noise < percent
    salt_mask = np.expand_dims(salt_mask, axis=2)
    salt_mask = np.repeat(salt_mask, 3, axis=2)
    img = img * (1 - salt_mask) + 255 * salt_mask

    # 添加胡椒噪声（黑色 = [0, 0, 0]）
    pepper_mask = noise > 1 - percent
    pepper_mask = np.expand_dims(pepper_mask, axis=2)
    pepper_mask = np.repeat(pepper_mask, 3, axis=2)
    img = img * (1 - pepper_mask) + 0 * pepper_mask

    # 确保输出是 uint8 类型
    return img.astype(np.uint8)

# 高斯噪声

感觉图片变模糊了

def add_Gaussian(img, sigma=20, mean=0):
    lab = cv2.cvtColor(img, cv2.COLOR_BGR2LAB)
    noise = np.random.normal(mean, sigma, lab[:, :, 0].shape)
    lab = lab.astype(float)
    lab[:, :, 0] = lab[:, :, 0] + noise
    lab[:, :, 0] = np.clip(lab[:, :, 0], 0, 255)
    lab = lab.astype(np.uint8)
    return cv2.cvtColor(lab, cv2.COLOR_LAB2BGR)

# 滤波技术

# 均值滤波

利用方波信号作为卷积核

$\boldsymbol{K}=\frac{1}{k^{2}}\left[\begin{array}{ccc} 1 &\cdots &1\\ \vdots &&\vdots\\ 1 &\cdots &1 \end{array}\right]_{k\times k}$

对上述添加了高斯噪声的图片进行均值滤波，得到

图片变得平滑，但也变得模糊，并且有 “块状效应”

# 高斯滤波

高斯滤波卷积核为

$G(i,j)=\frac{1}{2\pi \sigma^2}e^{-\frac{i^2+j^2}{2\sigma^2}}$

高斯滤波的效果更柔和，而且不会有块状效应，但图像会变的模糊

def gaussian_kernel(k_size, sigma):
    # k_size 核的尺寸，必须是奇数，这样才能找到 center 位置
    # sigma 是高斯分布的标准差，决定“模糊程度”，值越大，核越“平”
    # G(i, j) = 1 / (2 * pi * sigma^2) * e^{-(x^2 + y^2) / (2 * sigma^2)}

    k = np.zeros((k_size, k_size))
    center = k_size // 2    # 中心位置
    s = sigma ** 2  # sigma^2 方差
    sum_val = 0

    for i in range(k_size):
        for j in range(k_size):
            x, y = i - center, j - center
            k[i, j] = np.exp(-(x ** 2 + y ** 2) / (2 * s)) / (2 * np.pi * s)
            sum_val += k[i, j]

    return k / sum_val

# 双边滤波

图像 $I$ 在双边滤波后输出的图像 $\hat{I}$ ，他们之间的关系为

$\hat{I}(i,j)=\frac{1}{W_{sb}}\sum\limits_{m=i-k}^{i-1}\sum\limits_{n=j-k}^{j-1}I(m,n)·G(i-m,j-n)·G(I(m,n)-I(i,j))$

其中

G_s(i,j)=\frac{1}{2\pi\sigma_s^2}e^{-\frac{i^2+j^2}{2\sigma^2}}$$$$G_b(d)=\frac{1}{2\pi\sigma_b^2}e^{-\frac{d^2}{2\sigma_b^2}}

$W_{sb}=\sum\limits_{m=i-k}^{i-1}\sum\limits_{n=j-k}^{j-1}G_s(i-m,j-n)·G_b(|I(m,n)-I(i,j)|)$

经双边滤波后，得到图像

在去噪的同时，很好地保留了清晰度

def bilateral_filter_custom(img, k_size=5, sigma_s=5, sigma_r=25):
    img = img.astype(np.float32)    # 转换为浮点类型，防止数值溢出
    half = k_size // 2  # 核的"半径"
    h, w = img.shape    # 图像的高和宽

    # 对图像边界进行反射填充（防止边界黑边）
    img_pad = np.pad(img, ((half, half), (half, half)), mode='reflect')
    # 在图像上下左右都填充 half 个像素，确保卷积时，边缘像素也会被处理到

    # 输出图像
    output = np.zeros_like(img, dtype=np.float32)

    # 预计算空间高斯核 G_s(i,j)
    Gs = gaussian_kernel(k_size, sigma_s)

    for i in range(h):
        for j in range(w):
            region = img_pad[i : i + k_size, j : j + k_size]    # 取一块区域
            center_val = img_pad[i + half, j + half]    # 这块区域的中心像素

            # 颜色差权重 G_r
            Gr = np.exp(-((region - center_val) ** 2) / (2 * sigma_r ** 2))
            Gr /= (2 * np.pi * sigma_r ** 2)

            # 联合权重 W = Gs * Gr
            W = Gs * Gr

            # 加权平均
            output[i, j] = np.sum(W * region) / (np.sum(W) + 1e-8)

    return np.clip(output, 0, 255).astype(np.uint8)

# 中值滤波

祛痘是美颜时常用的操作，痘 $\to$ 椒盐噪声，高斯滤波不能有效地处理，可以用中值滤波
中值滤波器是一种统计排序滤波器，图像像素等于周围像素排序后的中值，可以有效抑制椒盐噪声

def median_filter(img, k_size):
    img = img.astype(np.float32)
    half = k_size // 2
    h, w = img.shape

    res = np.zeros_like(img, dtype=np.float32)
    img = np.pad(img, ((half, half), (half, half)), mode='reflect')

    for i in range(h):
        for j in range(w):
            region = img[i : i + k_size, j : j + k_size]
            res[i, j] = np.median(region)

    return np.clip(res, 0, 255).astype(np.uint8)

# 图像锐化

图像平滑是过滤掉高频信息，保留低频信息，图像锐化是反其道而行
如果用原来图减去滤波后的图像，会得到全黑图像（这里是完全去除低频信息，而图像锐化是增强高频信息）

k_size = 3
k0 = np.zeros((k_size, k_size))
center = k_size // 2
k0[center][center] = 1

kg = gaussian_kernel(k_size, 1)

k_sharpen = 1 * (k0 - kg)

img = cv2.imread('avatar_Gaussian.jpg')

conv = conv_2d(img, k_sharpen)
conv.plot_conv()

而我们不妨，在原图的基础上，加上上面的卷积核， k_sharpen = k0 + 1 * (k0 - kg)
这样我们就能使得图像锐化
由此可见，将原图减去高斯滤波后的图像，即图像的高频信息，反复叠加在原图上，可以不断增强锐化效果

# 图像匹配

# 匹配步骤

模板图像 $T~(w\times h)$ 和输入图像 $I~(W \times H)$

让模板图像在输入图像上滑动，对于每个滑动位置，计算相似度
滑动结束后，得到相似度矩阵 $R$ ，查找元素最大值位置 $(i,j)$ 即可

# 相似度度量

# 互相关

$R(i,j)=\sum\limits_{m=i}^{i+2-1}\sum\limits_{n=j}^{j+h-1}I(m,n)·T(m-i,n-j)$

# 计算互相关矩阵
def COR(img, tmp):
    w, h = tmp.shape
    W, H = img.shape

    img = np.array(img, dtype=np.float32)
    tmp = np.array(tmp, dtype=np.float32)
    res = np.zeros((W - w + 1, H - h + 1), dtype=np.float32)

    for i in range(W - w + 1):          # 行
        for j in range(H - h + 1):      # 列
            patch = img[i : i + w, j : j + h]
            res[i, j] = np.sum(tmp * patch)

    return res

# 模板匹配类
class Tmp_match():
    def __init__(self, img, tmp):
        self.img = img
        self.tmp = tmp

    def match(self):
        # 转为灰度图
        img_gray = cv2.cvtColor(self.img, cv2.COLOR_BGR2GRAY)
        tmp_gray = cv2.cvtColor(self.tmp, cv2.COLOR_BGR2GRAY)
        w, h = tmp_gray.shape

        # 计算互相关矩阵
        res = COR(img_gray, tmp_gray)

        # 找到最大匹配位置
        loc = np.unravel_index(np.argmax(res), res.shape)
        top = [int(loc[1]), int(loc[0])]                # [x, y]
        bottom = [top[0] + h, top[1] + w]               # 注意宽高是反过来的

        # 画出矩形框（在原彩色图上）
        cv2.rectangle(self.img, tuple(top), tuple(bottom), (0, 255, 0), 2)

        # 正确显示颜色：BGR → RGB
        plt.imshow(cv2.cvtColor(self.img, cv2.COLOR_BGR2RGB))
        plt.title("Template Matching Result")
        plt.axis('off')
        plt.show()

# 主程序部分
img = cv2.imread('avatar.jpg')
tmp = cv2.imread('img.png')

test = Tmp_match(img, tmp)
test.match()

我们发现，这个代码匹配很容易出问题，匹配到奇奇怪怪的地方
因为互相关是模板图像像素值和输入图像子图像素值直接相乘，那么模板一定的情况下，子图整体像素值越大（越亮），乘积越大，因此更容易找到亮度高的地方，而不是相似的地方

# 归一化互相关

$R(i,j)=\frac{\sum\limits_{m}\sum\limits_{n}I(m,n)·T(m-i,n-j)}{\sqrt{\sum\limits_{m}\sum\limits_{n}I^2(m,n)}·\sqrt{\sum\limits_{m}\sum\limits_{n}T^2(m-i,n-j)}}$

def COR(img, tmp):
    w, h = tmp.shape
    W, H = img.shape

    img = np.array(img, dtype=np.float32)
    tmp = np.array(tmp, dtype=np.float32)
    res = np.zeros((W - w + 1, H - h + 1), dtype=np.float32)
    t = np.sqrt(np.sum(tmp ** 2))

    for i in range(W - w + 1):          # 行
        for j in range(H - h + 1):      # 列
            patch = img[i : i + w, j : j + h]
            res[i, j] = np.sum(tmp * patch)

            res[i, j] = res[i, j] / t / np.sqrt(np.sum(patch ** 2))

    return res

# 多目标模板匹配

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

# 快速选择第 k 大的数，时间复杂度 O(n)
def findKth(s, k):
    def select(arr, k):
        idx = random.randint(0, len(arr) - 1)
        pivot = arr[idx]

        low = [x for x in arr if x > pivot]
        equal = [x for x in arr if x == pivot]
        high = [x for x in arr if x < pivot]

        if k < len(low):
            return select(low, k)
        elif k < len(low) + len(equal):
            return pivot
        else:
            return select(high, k - len(low) - len(equal))

    flat = s.flatten().tolist()
    if k < 1 or k > len(flat):
        raise ValueError("k 超出范围")
    return select(flat, k - 1)


# 多目标匹配类
class MultiMatch():
    def __init__(self, img, tmp, k):
        self.img = img
        self.tmp = tmp
        self.k = k

    def match(self):
        # 转为灰度图
        img_gray = cv2.cvtColor(self.img, cv2.COLOR_BGR2GRAY)
        tmp_gray = cv2.cvtColor(self.tmp, cv2.COLOR_BGR2GRAY)
        h, w = tmp_gray.shape

        # 模板匹配
        res = cv2.matchTemplate(img_gray, tmp_gray, cv2.TM_CCORR_NORMED)

        # 找到第 k 大的得分作为阈值
        threshold = findKth(res, self.k)
        loc = np.where(res >= threshold)

        # 画矩形框
        for pt in zip(*loc[::-1]):  # 注意坐标顺序 (x, y)
            cv2.rectangle(self.img, pt, (pt[0] + w, pt[1] + h), (0, 255, 0), 2)

        # 显示结果（BGR 转 RGB）
        plt.figure(figsize=(10, 8))
        plt.imshow(cv2.cvtColor(self.img, cv2.COLOR_BGR2RGB))
        plt.axis('off')
        plt.show()


img = cv2.imread('multi_match/all_mul_match.png')
tmp = cv2.imread('multi_match/smt.png')
matcher = MultiMatch(img, tmp, k=3)
matcher.match()

# 边缘检测

边缘是图像中局部亮度、纹理或颜色发生剧烈变化的位置（边缘蕴含图像中物体的形状信息和场景的结构信息）

# 数学模型

给定一幅图像 $I$ ，边缘检测的目标是计算出图像中每一点 $(i,j)$ 的

边缘强度 $E_s(i,j|I)$
边缘方向 $E_\theta(i,j|I)$

图像 $I$ 是一个矩阵，可以表示为一个二元函数 $f(x,y)$ 的输出，而二元函数的微分是

$\frac{\partial f(x,y)}{\partial x}=\lim_{x \to 0} \frac{f(x+\varepsilon,y)-f(x-\varepsilon,y)}{2\varepsilon}$

$\frac{\partial f(x,y)}{\partial y}=\lim_{\varepsilon \to 0} \frac{f(x,y+\varepsilon)-f(x,y-\varepsilon)}{2\varepsilon}$

由于数字图像是离散的，在数字图像中最小的位置距离单位为 $1$ ，即 $\varepsilon = 1$ ，那么在离散情况下有

$I_x(i,j)=\frac{I(i+1,j)-I(i-1,j)}{2}$

$I_y(i,j)=\frac{I(i,j+1)-I(i,j-1)}{2}$

其中， $I_x,I_y$ 分别是图像像素值关于像素 $(i,j)$ 沿 $x,y$ 轴方向的一阶导数，记为 $\nabla I(i,j)=[I_x(i,j),I_y(i,j)]$ ，表示图像在该点处的梯度

在数字图像中，对图像的一阶微分是通过一阶差分实现的

$I_x(i,j)=\frac{1}{2}\sum\limits_{m=-1}^{1}\sum\limits_{n=-1}^{1}I(i-m,j-n)·K_x(m,n)$

$I_y(i,j)=\frac{1}{2}\sum\limits_{m=-1}^{1}\sum\limits_{n=-1}^{1}I(i-m,j-n)·K_y(m,n)$

其中

$K_x=\begin{bmatrix} 0 & 0 & 0 \\ 1 & 0 & -1 \\ 0 & 0 & 0 \end{bmatrix},~K_y=\begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & -1 & 0 \end{bmatrix}$

边缘强度和边缘方向分别为

$E_s(i,j|I)=\sqrt{I_x^2(i,j)+I_y^2(i,j)}$

$E_\theta(i,j|I)=\tan^{-1}(\frac{I_y(i,j)}{I_x(i,j)})+\frac{\pi}{2}$

# 算法

边缘检测算法主要有三个步骤

图像平滑：边缘检测算法基于图像微分，但是图像微分对噪声很敏感
图像梯度计算：通过差分模板对平滑处理后的图像进行卷积，得到图像梯度
图像边缘定位：根据边缘强度和边缘方向进行定位

# Sobel 边缘检测算法

$K_x=\begin{bmatrix} 1 & 0 & -1 \\ 2 & 0 & -2 \\ 1 & 0 & -1 \end{bmatrix},~K_y=\begin{bmatrix} 1 & 2 & 1 \\ 0 & 0 & 0 \\ -1 & -2 & -1 \end{bmatrix}$

如何理解 $K_x$ 呢？

$K_x=\begin{bmatrix} 1 & 0 & -1 \\ 2 & 0 & -2 \\ 1 & 0 & -1 \end{bmatrix}=\begin{bmatrix}1 \\ 2 \\ 1\end{bmatrix}\begin{bmatrix}1 & 0 & -1\end{bmatrix}$

[1 2 -1]^T 可以理解成一个高斯滤波， [1 0 -1] 是一个差分算子

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

img = cv2.imread('scenery.jpg')

k_x = np.array([
    [1, 0, -1],
    [2, 0, -2],
    [1, 0, -1]
])

k_y = np.array([
    [1, 2, 1],
    [0, 0, 0],
    [-1, -2, -1]
])

conv_x = cv2.filter2D(img, -1, k_x)
conv_y = cv2.filter2D(img, -1, k_y)
plt.imshow(conv_x[:, :, ::-1])
plt.imshow(conv_y[:, :, ::-1])
plt.axis('off')
plt.show()

在计算边缘强度时，通常使用不开平方的近似值

$E_s(i,j|I)=||I_x(i,j)|| + ||I_y(i,j)||$

conv_x = cv2.filter2D(img, -1, k_x)
conv_y = cv2.filter2D(img, -1, k_y)

res = abs(conv_x) + abs(conv_y)

plt.imshow(res[:, :, ::-1])
plt.axis('off')
plt.show()

# Canny 边缘检测算法

# 图像平滑

$\sigma$ 越大，图像越平滑，但也越模糊，一般高斯核是 $\sigma$ 的 $3$ 倍时，滤波效果较好

k_size = 5
sigma = 1.5

img = cv2.imread('./edge_detect/a.jpg')

kernel = gaussian_kernel(k_size, sigma)

img = cv2.filter2D(img, -1, kernel)

plt.axis('off')
plt.imshow(img[:, :, ::-1])
plt.show()

# 梯度计算

使用和 $Sobel$ 边缘检测算法一样的方法，计算

梯度大小（边缘强度）

$E_s(i,j|I)=\sqrt{I_x^2(i,j)+I_y^2(i,j)}$

梯度方向

$E_\theta(i,j|I)=\tan^{-1}(\frac{I_y(i,j)}{I_x(i,j)})+\frac{\pi}{2}$

# 求沿 x 轴方向的梯度幅值
def partial_x(img):
    k = np.array([
        [0, 0, 0],
        [0.5, 0, -0.5],
        [0, 0, 0]
    ])

    return cv2.filter2D(img, -1, k)


# 求沿 y 轴方向的梯度幅值
def partial_y(img):
    k = np.array([
        [0, 0.5, 0],
        [0, 0, 0],
        [0, -0.5, 0]
    ])

    return cv2.filter2D(img, -1, k)

# 求梯度幅值及方向（处理彩色图像）
def gradient(img):
    # 分别计算每个通道的梯度
    channels = []
    thetas = []

    for i in range(3):  # 对BGR三个通道分别处理
        channel = img[:, :, i]
        dx = partial_x(channel)
        dy = partial_y(channel)

        G = np.sqrt(dx ** 2 + dy ** 2)
        theta = np.rad2deg(np.arctan2(dy, dx))
        theta %= 360

        channels.append(G)
        thetas.append(theta)

    # 合并三个通道的梯度幅值（可以取最大值或平均值）
    G = np.maximum.reduce(channels)
    # 或者使用平均值: G = np.mean(channels, axis=0)

    # 对于方向，我们通常只取一个通道的（比如绿色通道）
    theta = thetas[1]  # 取绿色通道的方向

    return G, theta

# 计算梯度
img = img.astype(np.float32)
G, theta = gradient(img)

# 显示梯度幅值
plt.subplot(1, 2, 2)
plt.axis('off')
plt.title('final img')
plt.imshow(G, cmap='gray')  # 用于显示灰度图
plt.show()

直接从梯度得到的边缘线很粗，显得很模糊，因此需要通过非极大值抑制将边缘细化
非极大值抑制沿每个梯度方向找到边缘强度局部极大的像素予以保留，去除该方向上边缘强度非局部极大的像素

# 边缘强度非极大值抑制

# 非极大值抑制
def non_max_suppress(G, theta):
    h, w = G.shape
    res = np.zeros((h, w))

    # 角度规范化至 [0, 180)，因为梯度方向是双向的
    theta = theta % 180

    for i in range(1, h - 1):
        for j in range(1, w - 1):
            a = theta[i, j]
            now = G[i, j]

            if 0 <= a <= 22.5 or 157.5 <= a <= 180:
                n1 = G[i, j + 1]
                n2 = G[i, j - 1]
            elif 22.5 <= a < 67.5:
                n1 = G[i - 1, j + 1]
                n2 = G[i + 1, j - 1]
            elif 67.5 <= a < 112.5:
                n1 = G[i - 1, j]
                n2 = G[i + 1, j]
            else:
                n1 = G[i - 1, j - 1]
                n2 = G[i + 1, j + 1]

            if now >= n1 and now >= n2:
                res[i, j] = now
            else:
                res[i, j] = 0

    return res

# 基于滞后的阈值化的边缘定位

经过非极大值抑制后，将边缘细化成边缘线，接下来我们要确定边缘的位置
我们对边缘强度 $E_s$ 进行二值化，将其转为 $E_b$ ，以 $E_b(i,j | I)=1$ 是否成立来确定 $(i,j)$ 是否是边缘点，可以利用单一阈值法进行二值化 0

$E_b(i,j|I)=\begin{cases}1, & E_s(i,j|I) \ge T \\ 0, & E_s(i,j | I) < T\end{cases}$

观察有，单一阈值法的图像容易造成断裂的边缘，我们此处提出双阈值法（基于滞后的阈值法）

$E_b(i,j|I)=\begin{cases}1, & E_s(i,j | I) \ge T_h \\ 0, & E_s(i,j | I) < T_l \\ unknow,& T_l \le E_s(i,j | I) < T_h \end{cases}$

强边缘点是已经确定的边缘点，而对于弱边缘点，当且仅当它们在已经确定的边缘点的领域内，才能被确定为边缘点

# 强弱边缘点的合并

def merge_edge(strong_res, weak_res):
    res = np.zeros(strong_res.shape)
    res[strong_res] = 255 # 强边缘保留

    keep = weak_res.copy()
    flag = True

    while flag:
        flag = False

        for i in range(1, res.shape[0] - 1):
            for j in range(1, res.shape[1] - 1):
                # 是弱边缘点且最终未被保留
                if keep[i, j] and not res[i, j]:
                    # 如果邻域有强边缘或已保留边缘
                    if np.any(res[i - 1 : i + 2, j - 1 : j + 2] == 255):
                        res[i, j] = 255
                        flag = True

    return res

可能效果还不如

Sobel

，但是边缘（明显的边缘）更精准

# 角点检测

前面提到，模板匹配是最基础、最简单的图像匹配（模板和待匹配子图必须有极高相似度）。那么一般情况下的图像匹配如何实现呢？这需要检测图像中的特征进行匹配，局部特征通常被称为关键点特征或者兴趣点特征，他们有以下四个性质

重复性
显著性
紧性
鲁棒性

# Harris 角点检测算法

对于图像每个位置，计算窗口移动前后其内部的像素值变化量
计算像素值变化量对应的角点响应函数，并对该函数进行阈值处理，提取角点

# 计算像素值变化量

在图像 $I$ 中某一位置放置一个窗口 $W$ ，将其移动一个微小位移 $(u,v)$ ，定义窗口移动前后的像素值变化量 $E(u,v)$ 如下

$E(u,v)=\sum\limits_{(i,j) \in W}(I(i+u, j+v)-I(i,j))^2$

如果窗口位于平坦区域，无论位置 $(u,v)$ 指向什么方向， $E(u,v)$ 的值都较小
位于边缘时，沿着边缘移动， $E(u,v)$ 也较小
当位于角点的时候，无论往哪移动， $E(u,v)$ 都较大

利用泰勒展开获取近似形式 $I(i+u,j+v) \approx I(i,j)+uI_x+vI_y$
其中， $I_{x,y}$ 是图像沿着 $x,y$ 轴的梯度幅值，原式可化为

$E(u,v)=\sum\limits_{(i,j) \in w}(u^2I_x^2+v^2I_y^2+2uvI_xI_y)$

提出 $(u,v)$

$E(u,v) = \begin{bmatrix}u & v\end{bmatrix} M \begin{bmatrix}u \\ v\end{bmatrix}$

矩阵 $M$ 为

$\boldsymbol{M} = \sum_{(i,j) \in \boldsymbol{W}} \begin{bmatrix} I_{x}^{2} & I_{x}I_{y} \\ I_{x}I_{y} & I_{y}^{2} \end{bmatrix}$

$M=\sum\limits_{(i,j) \in w}H$ ，其中 $H$ 是图像 $I$ 的黑塞矩阵

$M=R\begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix}R^{T}$

其中， $\lambda_1,\lambda_2$ 是矩阵 $M$ 的特征值， $R$ 是二维旋转矩阵，代入 $E(u,v)$ 有

$E(u,v) = \begin{bmatrix} u & v \end{bmatrix}\boldsymbol{R}\begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix}\boldsymbol{R}^\top\begin{bmatrix} u \\ v \end{bmatrix}$

令 $[m~n]=[u~v]R$ 可得

$E(u,v) = \left[\begin{array}{ll} m & n \end{array}\right] \left[\begin{array}{cc} \lambda_{1} & 0 \\ 0 & \lambda_{2} \end{array}\right] \left[\begin{array}{ll} m & n \end{array}\right]^{\top} = m^{2}\lambda_{1} + n^{2}\lambda_{2}$

由此，像素值变化量 $E(u,v)$ 可以直接与矩阵 $M$ 的两个特征值直接相关

# 计算角点响应函数

角点附近的 $E(u,v)$ 对于任意 $(u,v)$ 都很大，又有 $[m~n]=[u~v]R$ ，则有

$m^{2}+n^{2}= \left[\begin{array}{ll} m & n \end{array}\right] \left[\begin{array}{l} m \\ n \end{array}\right]= \left[\begin{array}{ll} u & v \end{array}\right] \boldsymbol{R}\boldsymbol{R}^{\top} \left[\begin{array}{l} u \\ v \end{array}\right]= u^{2}+v^{2}$

假设微小位移 $(u,v)$ 近似于单位向量，即 $u^2+v^2 \approx 1 \to u^2 + v^2 \approx 1$ ，可以得到

$\max_{u,v}E(u,v)=max(\lambda_1,\lambda_2)$

$\min_{u,v}E(u,v)=min(\lambda_1,\lambda_2)$

如上图

两个特征值都小，且近似相等：对应图像平坦区域
一个大、一个小，且相差较大：对应图像边缘区域
两个都大，且近似相等：对应图像角点

虽然可以根据大小关系判断角点，但是我们还是希望通过一个单一的度量指标，根据经验，提出下方角点响应函数

$\rho = \lambda_1 \lambda_2 - \alpha (\lambda_1 + \lambda_2)^2$

其中， $\alpha$ 是一个常数，通常取 $0.04 \to 0.06$

根据线性代数知识有

$\rho = \lambda_{1}\lambda_{2} - \alpha(\lambda_{1}+\lambda_{2})^{2} = \det(\boldsymbol{M}) - \alpha \operatorname{trace}(\boldsymbol{M})^{2}$

计算出 $\rho$ 后，与设置的阈值 $\gamma$ 进行比较，若 $\rho > \gamma$ ，则视为角点

# 代码实现

C++

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

from tools import *

# 角点响应函数
def response(M):
    # 超参数
    k = 0.04

    det = np.linalg.det(M)
    trace = np.trace(M)

    # 角点响应函数
    R = det - k * trace ** 2

    return R


# Harris 角点检测算法
def harris(img, NMS=True):
    h, w = img.shape[:2]
    img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

    point = np.zeros((h, w))

    grad = np.zeros((h, w, 2), dtype=np.float32)

    grad[:, :, 0] = cv2.Sobel(img, cv2.CV_16S, 1, 0)
    grad[:, :, 1] = cv2.Sobel(img, cv2.CV_16S, 0, 1)

    Ix = grad[:, :, 0] ** 2
    Iy = grad[:, :, 1] ** 2
    Ixy = grad[:, :, 0] * grad[:, :, 1]

    Ix = cv2.GaussianBlur(Ix, (3, 3), sigmaX=2)
    Iy = cv2.GaussianBlur(Iy, (3, 3), sigmaX=2)
    Ixy = cv2.GaussianBlur(Ixy, (3, 3), sigmaX=2)

    for i in range(img.shape[0]):
        for j in range(img.shape[1]):
            # 构建 M 矩阵
            sm = [[Ix[i, j], Ixy[i, j]], [Ixy[i, j], Iy[i, j]]]

            # 计算角点响应函数
            R = response(sm)

            point[i, j] = R

    # 非极大值抑制
    corner = np.zeros_like(img, dtype=np.float32)
    threshold = 0.01

    maxVal = np.max(point)

    for i in range(point.shape[0]):
        for j in range(point.shape[1]):
            if NMS:
                if point[i][j] > threshold * maxVal \
                        and point[i][j] == np.max(point[max(0, i - 1) : min(i + 1, h - 1), max(0, j - 1) : min(j + 1, w - 1)]):
                    corner[i, j] = 255
            else:
                if point[i][j] > threshold * maxVal:
                    corner[i, j] = 255

    return corner


if __name__ == "__main__":
    img = cv2.imread("./img/qp.png")
    
    # 可以先进行锐化（锐化可以提高角点检测的效果）
    k_size = 3
    k0 = np.zeros((k_size, k_size))
    center = k_size // 2
    k0[center][center] = 1
    kg = gaussian_kernel(k_size, 1)
    k_sharpen = k0 + 1 * (k0 - kg)
    conv = conv_2d(img, k_sharpen)
    for i in range(2):
        img = conv.conv(img, k_sharpen)
        conv = conv_2d(img, k_sharpen)

    plt.imshow(img[:, :, ::-1])
    plt.axis('off')
    plt.title("Original Image")

    response = harris(img)
    img[response == 255] = [0, 0, 255]
    plt.imshow(img[:, :, ::-1])
    plt.axis('off')
    plt.title("Harris Corner")
    plt.show()

# 基于尺度空间的块状区域检测的基本原理

由上一小章可以知道，角点对于尺度是不变的

随着高斯拉普拉斯算子标准差增大，当高斯拉普拉斯算子的标准差大小和块状区域大小相同时，他们的卷积结果会出现唯一极值点，此时记为 $\sigma^{*}$
但是，此极值点是在某一特定标准差 $\sigma$ （尺度）下 $x$ 取值范围内的唯一极值点，如果考虑多个尺度，即整个 $(x, \sigma)$ 空间，该极值并不是一个全局最小值

出现这种现象，是因为高斯拉普拉斯算子的响应会随着标准差的变大而衰减，因此需要对响应做归一化，即对其乘以 $\sigma^2$

在归一化后，此时特征尺度下的极值点是整个 $(x,\sigma)$ 空间上的全局最小值，只要检测到 $(x,\sigma)$ 空间归一化高斯拉普拉斯算子响应的全局最小值，就可以检测到块状区域，并确定其特征尺度

# 迁移到图像中

图像的尺度空间 $L(x,y,\sigma)$ 可通过一个不同尺度的高斯函数 $G(x,y,\sigma)$ 与输入图像 $I(x,y)$ 卷积得到，在该尺度空间 $(x,y,\sigma)$ 处的值 $L(x,y,\sigma)$ 为

$L(x,y,\sigma)=G(x,y,\sigma) * I(x,y)$

由滤波的知识可知

$\sigma$ 越小（即尺度越小），图像越清晰
$\sigma$ 越大（即尺度越大），图像越模糊

如果找到了两幅图像中对应的特征点，对特征点对应的特征尺度求比值，就可以得到两幅图像尺度缩放的比例

# SIFT 算法

计算机视觉