AnyMath 文档
Notebook

关键点检测与匹配

在现代计算机视觉中,关键任务是找到图像之间的几何对应关系。 无论是创建全景图,跟踪物体的运动,还是从照片中重建3D场景,整个过程都基于检测独特点并匹配它们。

所呈现的代码是经典特征检测和匹配管道的教学性但功能性的实现。 与使用现成库(例如OpenCV)的"黑匣子"不同,Julia中的这个例子演示了手动编程。 它不仅展示了如何应用算法,还揭示了每个阶段的数学本质算法细节

接下来,连接库:

*Images:用于处理图像的基本软件包,提供数据类型和基本输入/输出功能。
*TestImages:标准测试图像库。
*ImageFeatures:包含用于检测关键点的现成算法,尽管此示例实现了自己版本的Harris检测器。
*ImageDraw:专为渲染图形原语而设计,其代码中的功能被用于教育目的的手动实现所取代。
*ImageFiltering:提供Sobel算子计算梯度所需的滤波和卷积工具。
*CoordinateTransformations:描述和组成坐标空间几何变换的框架。
*Rotations:在计算机视觉任务中处理旋转和旋转矩阵的专用软件包。
*StaticArrays:为高效的矢量数学提供高性能的固定大小数组。
*LinearAlgebra:用于求解单应性估计中的方程组的标准线性代数库。
*Statistics:用于计算统计指标的标准模块,例如描述符的均值和标准差。

In [ ]:
# Pkg.add(["Images", "TestImages", "ImageFeatures", "ImageDraw", "ImageFiltering", "CoordinateTransformations", "Rotations", "StaticArrays"])
using Images, TestImages, CoordinateTransformations, Rotations, StaticArrays, LinearAlgebra, Statistics

接下来,我们将对辅助渲染函数进行详细分析,这两个函数都实现了经典的Bresenham算法-仅使用整数算术的几何图元的有效光栅化方法。

功能 draw_line! 线绘制的Bresenham算法。 该算法在点之间建立一条线 在每个步骤中选择最接近理想直线的像素。

直线的方程:

而不是计算 (这需要除法和浮动算术),Bresenham使用误差项函数,该函数显示当前像素与理想线的偏差。

此外,在每个步骤中,算法评估当前点偏离理想直线的程度。 差异函数根据规则进行更新:

乘以2(e2 = 2 * err)允许您避免使用小数,同时保持比较的准确性。

功能 draw_circle!:中点圆算法。 该算法建立一个半径圆 与中心 使用八向对称:如果点 如果它位于一个圆上,那么与它对称的点也属于它:

圆的方程:

**决策参数评估两个候选之间中点的位置:

公式****解决方案参数更新

参数 它会定期更新,而无需在每个步骤中计算正方形。:

这些公式是通过分解得到的 和重复计算的排除。

In [ ]:
function draw_line!(img, p1, p2, color)
    x1, y1 = round(Int, p1[1]), round(Int, p1[2])
    x2, y2 = round(Int, p2[1]), round(Int, p2[2])
    dx = abs(x2 - x1)
    dy = abs(y2 - y1)
    sx = x1 < x2 ? 1 : -1
    sy = y1 < y2 ? 1 : -1
    err = dx - dy
    x, y = x1, y1
    while true
        if 1 <= x <= size(img,2) && 1 <= y <= size(img,1)
            img[y, x] = color
        end
        if x == x2 && y == y2
            break
        end
        e2 = 2 * err
        if e2 > -dy
            err -= dy
            x += sx
        end
        if e2 < dx
            err += dx
            y += sy
        end
    end
end

function draw_circle!(img, center, radius, color)
    cx, cy = round(Int, center[1]), round(Int, center[2])
    r = radius
    x = 0
    y = r
    d = 3 - 2r
    while x <= y
        for (dx, dy) in [(x, y), (y, x), (-x, y), (-y, x), (x, -y), (y, -x), (-x, -y), (-y, -x)]
            px = cx + dx
            py = cy + dy
            if 1 <= px <= size(img,2) && 1 <= py <= size(img,1)
                img[py, px] = color
            end
        end
        if d < 0
            d += 4x + 6
        else
            d += 4(x - y) + 10
            y -= 1
        end
        x += 1
    end
end
Out[0]:
draw_circle! (generic function with 1 method)

哈里斯角探测器: harris_corners

该函数实现了经典的Harris-Stevens算法(1988)**来检测角度,这些角度是抵抗旋转,光照变化和轻微几何扭曲的特征图像点。

角度是图像强度在两个方向上同时显着变化的点。 该算法分析每个像素周围的局部窗口,并评估它的"角度"。

第二矩的矩阵(结构张量),为每个像素构建一个矩阵 :

哪里:

  • -沿X和Y轴的图像渐变
    -沿着像素周围的窗口进行求和(在代码中,一个3×3窗口)

特征值分析 矩阵 :

/区域类型||/视觉|
|:---|:---|:---|:---|
//小/小/均匀面积|
/边框/大/小/对象的边缘|
/角度/大/大/跨越边界|

而不是计算特征值(一个昂贵的操作),使用近似值。:

哪里:

  • -经验系数(通常为0.04-0.06)

回应的解释:

  • 角度(两个特征值都很大)
  • 边框(一大一小)
  • 平坦区域(均小)

,探测器操作的可视化

Исходное изображение → Градиенты → Матрица M → Отклик R → Порог → Локальные максимумы
     [512×512]        [Ix, Iy]     [Sxx,Sxy,Syy]  [R]    [>1% max]   [5×5 окно]

Пример откликов в разных областях:

Область Тип
Плоская 0.1 0.1 0.0 ~0.01
Граница 10.0 0.1 0.5 ~−0.4
Угол 10.0 10.0 2.0 ~+96

⚙️ Параметры и их влияние

Параметр Значение по умолчанию Эффект увеличения Эффект уменьшения
k 0.04 Меньше углов, только самые выраженные Больше углов, выше чувствительность к шуму
threshold 0.01 Только сильные углы (10% от максимума) Больше углов, включая слабые
Окно сглаживания 3×3 Более устойчиво к шуму Лучшая локализация, но шумнее
Окно NMS 5×5 Более разреженные точки Плотнее точки, возможны дубликаты
In [ ]:
function harris_corners(img; k=0.04, threshold=0.01)
    img_f = Float64.(img)
    img_f = replace(img_f, NaN => 0.0, Inf => 0.0, -Inf => 0.0)
    Ix = imfilter(img_f, Kernel.sobel()[1])
    Iy = imfilter(img_f, Kernel.sobel()[2])
    Ixx = Ix.^2
    Iyy = Iy.^2
    Ixy = Ix .* Iy
    window = ones(3,3) ./ 9
    Sxx = imfilter(Ixx, window)
    Syy = imfilter(Iyy, window)
    Sxy = imfilter(Ixy, window)
    det = Sxx .* Syy - Sxy.^2
    trace = Sxx + Syy
    response = det - k * trace.^2
    max_resp = maximum(response)
    corners = findall(response .> threshold * max_resp)
    local_max = []
    for c in corners
        i, j = c.I
        val = response[i,j]
        r_min = max(1, i-2)
        r_max = min(size(img,1), i+2)
        c_min = max(1, j-2)
        c_max = min(size(img,2), j+2)
        window_vals = response[r_min:r_max, c_min:c_max]
        if all(val .>= window_vals)
            push!(local_max, (j, i))  
        end
    end
    return local_max
end
Out[0]:
harris_corners (generic function with 1 method)

Пошаговый разбор кода

Шаг 1: Подготовка изображения

img_f = Float64.(img)
img_f = replace(img_f, NaN => 0.0, Inf => 0.0, -Inf => 0.0)

Зачем: Конвертация в Float64 необходима для точных вычислений градиентов. Замена NaN/Inf защищает от артефактов после геометрических трансформаций (warp).

Шаг 2: Вычисление градиентов (Оператор Собеля)

Ix = imfilter(img_f, Kernel.sobel()[1])
Iy = imfilter(img_f, Kernel.sobel()[2])

Ядра Собеля:

Формула свёртки:

Градиенты показывают скорость изменения интенсивности в каждом направлении.

步骤3:渐变的正方形和乘积

Ixx = Ix.^2
Iyy = Iy.^2
Ixy = Ix .* Iy

物理意义:

  • -x中梯度的"能量"
  • -y中梯度的"能量"
  • -梯度之间的相关性(显示变化的一致性)

步骤4:使用窗口过滤器平滑

window = ones(3,3) ./ 9
Sxx = imfilter(Ixx, window)
Syy = imfilter(Iyy, window)
Sxy = imfilter(Ixy, window)

**为什么:**在3×3窗口上取平均相当于矩阵公式中的求和 . 如果没有抗混叠,探测器将对噪声过于敏感。

权重窗口(在经典版本-高斯):

此代码使用统一窗口。 .

步骤5:计算Harris响应

det = Sxx .* Syy - Sxy.^2
trace = Sxx + Syy
response = det - k * trace.^2

每个像素的完整公式:

参数 k=0.04:

-较低值(0.04)→更多可检测角度(更高灵敏度)
-更高的值(0.06–0.08)→只有最明显的角度

步骤6:阈值过滤

max_resp = maximum(response)
corners = findall(response .> threshold * max_resp)

自适应阈值:

threshold=0.01 选择响应高于最大值1%的点。 这使得检测器能够抵抗图像对比度的变化。

步骤7:非最大抑制

for c in corners
    i, j = c.I
    val = response[i,j]
    r_min = max(1, i-2)
    r_max = min(size(img,1), i+2)
    c_min = max(1, j-2)
    c_max = min(size(img,2), j+2)
    window_vals = response[r_min:r_max, c_min:c_max]
    if all(val .>= window_vals)
        push!(local_max, (j, i))  # (x, y)
    end
    end

**原因:经过阈值滤波后,几个相邻像素可能会保留在一个角的旁边。 抑制非最大值只会在5x5窗口中留下局部最大值

算法:

  1. 对于每个候选者,在其周围采取5×5窗口。
  2. 如果中心的值不小于所有邻居→点被保存
  3. 否则,→被丢弃为"重复"

**结果:**一组稀疏的局部角,没有聚类。

关键点的过滤: filter_valid_keypoints -此功能对关键点进行空间验证,筛选出图像中不正确区域的点。 虽然该功能在算法上很简单,但它解决了确保后续流水线阶段数据质量的关键任务。

滤波的数学基础。 该函数将图像的有效区域定义为满足两个条件的像素子集:

1. 几何约束(边距):

哪里:

  • -图像的宽度和高度
  • -参数 margin (默认为10像素)

2. 强度限制:

哪里 -坐标中的像素强度 .

最终有效性条件

它被认为是有效的,如果:

参数及其影响

/参数/默认值/放大效果/缩小效果|
|:---|:---|:---|:---|
|margin/10/有效积分更少,可靠性更高/积分更多,越界风险/

选择保证金的建议:

/任务/建议保证金|理由/
|:---|:---|:---|
/提取7×7补丁|≥4|补丁大小的一半|
/提取补丁15x15/≥8/补丁大小的一半|
/强力转弯后/≥20/占大黑区|
/精确定位/≥3/最小边缘保护|

最小保证金公式:

In [ ]:
function filter_valid_keypoints(img, keypoints; margin=10)
    h, w = size(img)
    valid = []
    for (x, y) in keypoints
        if margin < x <= w-margin && margin < y <= h-margin
            if img[round(Int,y), round(Int,x)] != 0
                push!(valid, (x, y))
            end
        end
    end
    return valid
end
Out[0]:
filter_valid_keypoints (generic function with 1 method)

分步代码分析

步骤1:获取图像尺寸

h, w = size(img)

**重要:**在朱莉娅 size(img) 申报表 (высота, ширина) 那就是 , (rows, columns). 这对应于矩阵表示法,其中第一个索引是行(Y),第二个是列(X)。

步骤2:检查图像的边框

if margin < x <= w-margin && margin < y <= h-margin

为什么我需要保证金?:

-在提取补丁(描述符)时防止越界
-考虑几何变换的误差
-排除边缘过滤伪影

步骤3:检查像素强度

if img[round(Int,y), round(Int,x)] != 0
    push!(valid, (x, y))
end

**为什么重要:**几何变换(旋转,缩放)后,图像中出现黑色区域(用零填充):

**数学原因:**当图像旋转时,新坐标计算为:

在后向映射期间未落入原始图像的像素用值0(或NaN,其先前替换为0)填充。

步骤4:坐标转换

img[round(Int,y), round(Int,x)]  # Обратите внимание: [y, x]!

**关键时刻:**在计算机视觉中,坐标通常写成 ,在哪里:

  • -水平坐标(列)
  • -垂直坐标(行)

但在Julia中,数组被索引为 [строка, столбец] 这就是 , [y, x]. 函数正确执行此转换。

基于补丁提取描述符: extract_patch_descriptors

该函数实现了最简单的基于补丁的描述符,一种通过提取和归一化每个点周围的局部图像区域来生成关键点的矢量特征的方法。 虽然这种方法不如现代描述符(SIFT,ORB,深度学习),但它展示了不变特征形成的基本原则。

描述符的数学基础

什么是描述符?

****描述符是图像局部区域的数字表示,应该是:

  1. 独特-不同的点有不同的描述符
  2. 不变—不同条件下的相同点给出相似的描述符
  3. 紧凑-足够小,以便进行有效的比较

基于补丁的方法

而不是梯度的复杂直方图(如在SIFT中),使用raw patch-围绕关键点的像素矩形区域:

哪里 -补丁的一半大小。

描述符大小:

`

两个拒绝标准:

/标准/公式/阈值/原因|
|----------|---------|-------|---------|
/低方差| |1e-6/均匀区域(无纹理)|
/太多零| |30%/变换后黑色区域|

标准差公式:

为什么重要:
-❌补丁与 不包含信息(所有像素相同)
-具有>30%的零的补丁包含变换伪影

步骤6:描述符的归一化(Z-score)

desc = (desc .- mean(desc)) / (std(desc) + 1e-8)

归一化公式:

归一化描述符的属性:

为什么正常化?:

/问题|无归一化/有归一化|
|----------|-----------------|-----------------|
/光照变化/描述符不同|描述符相似|
/对比度/取决于亮度/不变|
|比较/欧几里德距离不正确|正确比较|

例子:

Исходный патч (разная яркость):
[100, 110, 120]     [200, 210, 220]
[105, 115, 125][205, 215, 225]  ← Разные значения!
[110, 120, 130]     [210, 220, 230]

После нормализации:
[-1.22, 0.00, 1.22]     [-1.22, 0.00, 1.22]
[-0.61, 0.61, 1.83]  =  [-0.61, 0.61, 1.83]  ← Одинаковые!
[0.00, 1.22, 2.44]      [0.00, 1.22, 2.44]

**增加ε=1e-8:**当 (虽然这样的补丁已经被过滤掉)。

步骤7:组装描述符矩阵

push!(descriptors, desc)
push!(valid_indices, idx)

return hcat(descriptors...), valid_indices

输出格式:

哪里:
-每个是一个关键点的描述符

  • valid_indices -已过滤的原始关键点的索引

例子:

Исходные keypoints: 297 точек
Валидные дескрипторы: 227 точек
Отсеяно: 70 точек (23.6%)

Матрица дескрипторов: 49 × 227

步骤8:处理空结果

if length(descriptors) == 0
    return zeros(0, 0), Int[]
end

**为什么:**如果所有补丁都被拒绝,则在后续阶段(映射,RANSAC)进行错误保护。

、参数及其影响

/参数/默认值/缩放效果/缩放效果|
|----------|----------------------|-------------------|-------------------|
| patch_size|7/更多的上下文,更高的不变性/更少的计算,更好的本地化|
| std threshold/1e-6/更少的剪报/更严格的过滤同质区域|
| zero ratio/0.3/更多描述符|更严格的工件过滤|

选择patch_size的建议:

/任务/建议大小|理由/
|--------|---------------------|-------------|
/快速匹配|5×5(25dim)/速度的最小尺寸|
|标准映射|7×7(49dim)/质量和速度的平衡|
/复杂转型|9×9(81dim)/更多可持续性背景|
/深度学习描述符/32×32(1024dim)|神经网络需要|

尺寸选择公式:

In [ ]:
function extract_patch_descriptors(img, keypoints; patch_size=7)
    img_f = Float64.(img)
    img_f = replace(img_f, NaN => 0.0, Inf => 0.0, -Inf => 0.0)
    descriptors = []
    offset = div(patch_size, 2)
    valid_indices = []
    
    for (idx, (x,y)) in enumerate(keypoints)
        x_int = clamp(round(Int, x), offset+1, size(img,2)-offset)
        y_int = clamp(round(Int, y), offset+1, size(img,1)-offset)
        patch = img_f[y_int-offset:y_int+offset, x_int-offset:x_int+offset]
        if std(patch) < 1e-6 || count(patch .== 0) / length(patch) > 0.3
            continue
        end
        desc = patch[:]
        desc = (desc .- mean(desc)) / (std(desc) + 1e-8)
        push!(descriptors, desc)
        push!(valid_indices, idx)
    end
    
    if length(descriptors) == 0
        return zeros(0, 0), Int[]
    end
    
    return hcat(descriptors...), valid_indices
end
Out[0]:
extract_patch_descriptors (generic function with 1 method)

分步代码分析

步骤1:图像准备

img_f = Float64.(img)
img_f = replace(img_f, NaN => 0.0, Inf => 0.0, -Inf => 0.0)

为什么:

  • Float64 精确计算统计数据(平均值,std)所必需的
    -更换 NaN/Inf 防止转换工件

步骤2:计算居中的偏移量

offset = div(patch_size, 7)  # Для patch_size=7 → offset=3

补丁边界公式:

步骤3:坐标限制(夹紧)

x_int = clamp(round(Int, x), offset+1, size(img,2)-offset)
y_int = clamp(round(Int, y), offset+1, size(img,1)-offset)

钳位功能:

**为什么:**确保补丁完全适合图像,即使关键点靠近边界。

可接受范围:

步骤4:提取补丁

patch = img_f[y_int-offset:y_int+offset, x_int-offset:x_int+offset]

补丁大小:

矢量化:

desc = patch[:]  # Преобразование 7×7 → 49×1

步骤5:过滤低质量补丁

if std(patch) < 1e-6 || count(patch .== 0) / length(patch) > 0.3
    continue
end
``

Сопоставление дескрипторов с Ratio Test: match_descriptors_euclidean

Эта функция реализует алгоритм сопоставления дескрипторов на основе евклидова расстояния с применением Ratio Test Лоу (Lowe's Ratio Test) — стандартного метода фильтрации ложных совпадений в компьютерном зрении. Хотя функция компактна, она содержит критически важную логику для обеспечения качества сопоставления.

Математическая основа сопоставления

Задача сопоставления (Matching Problem)

Даны два набора дескрипторов:

Цель: Найти пары , где и описывают одну и ту же физическую точку на разных изображениях.

Евклидово расстояние как метрика схожести

Для сравнения двух дескрипторов используется L2-норма (евклидово расстояние):

Где — размерность дескриптора (в нашем случае для патча 7×7).

Свойства:

  • → Идентичные дескрипторы
  • → Совершенно разные дескрипторы
  • Меньшее расстояние → Более вероятное совпадение

Ratio Test Лоу (1999)

Проблема: Простой выбор ближайшего соседа (nearest neighbor) даёт много ложных совпадений, особенно когда в сцене есть повторяющиеся текстуры.

Решение Lowe's Ratio Test: Сравнивать расстояние до первого и второго ближайших соседей:

Где:

  • — расстояние до ближайшего соседа (best match)
  • — расстояние до второго ближайшего соседа (second best)
  • — пороговый коэффициент (обычно 0.7–0.8)

直观的解释:

/情况| | /比率/产量|
|----------|-----------|-----------|-------|-------|
/绝对匹配 | 0.5 | 2.0 | 0.25 | ✅ 我们接受|
/不确定的匹配 | 0.5 | 0.6 | 0.83 | ❌ 拒绝|
/重复纹理 | 0.3 | 0.35 | 0.86 | ❌ 拒绝|

为什么有效:
-如果比例很小→第一个邻居明显优于第二个→肯定匹配
-如果比率接近1→两个邻居都同样好→模糊性(很可能是假匹配)

、参数及其影响

/参数/默认值/缩放效果/缩放效果|
|----------|----------------------|-------------------|-------------------|
| ratio_thresh|0.7/更多匹配,更低的精度/更少的匹配,更高的精度|

选择门槛的建议:

/任务/推荐问题|理由/
|--------|-----------------|-------------|
/高精度(RANSAC)| 0.6–0.7 | 减少RANSAC的排放量|
/最大保障范围| 0.8–0.85 | 更多的巧合,但更多的噪音|
/复杂场景(重复纹理)| 0.5–0.6 | 严格的歧义过滤|
/简单的场景(独特的功能)| 0.75–0.85 | 您可以松开过滤器|

劳的经验法则:

In [ ]:
function match_descriptors_euclidean(desc1, desc2, idx1, idx2; ratio_thresh=0.7)
    matches = []
    n1 = size(desc1, 2)
    n2 = size(desc2, 2)
    if n1 == 0 || n2 == 0
        return matches, Int[], Int[]
    end
    for i in 1:n1
        dists = [norm(desc1[:,i] - desc2[:,j]) for j in 1:n2]
        best_idx = argmin(dists)
        sorted_dists = sort(dists)
        if length(sorted_dists) >= 2 && sorted_dists[1] < ratio_thresh * sorted_dists[2]
            push!(matches, (idx1[i], idx2[best_idx]))
        end
    end
    return matches
end
Out[0]:
match_descriptors_euclidean (generic function with 1 method)

分步代码分析

步骤1:检查空描述符

n1=大小(desc1,2)
n2=大小(desc2,2)

如果n1==0|/n2==0
    返回匹配项,Int[],Int[]
结束

Зачем: Защита от ошибок, если предыдущие этапы не извлекли ни одного валидного дескриптора.

Формат дескрипторов:

  • size(desc, 2) возвращает количество дескрипторов (колонок)
  • size(desc, 1) возвращает размерность дескриптора (строк)

Шаг 2: Перебор всех дескрипторов первого изображения

对于我在1:n1
    #对于第一个集合中的每个描述符,我们在第二个集合中搜索匹配

Сложность: — наивный перебор всех пар

Пример:

1:227描述符
描述符2:227
距离检查:227×227=51.529

Шаг 3: Вычисление расстояний до всех кандидатов

dists=[norm(desc1[:,i]-desc2[:,j])for j in1:n2]

Математическая запись:

Векторизованная операция в Julia:

  • desc1[:,i] — i-й дескриптор из первого набора (вектор 49×1)
  • desc2[:,j] — j-й дескриптор из второго набора (вектор 49×1)
  • norm() — вычисляет L2-норму разности

Шаг 4: Поиск ближайшего соседа

best_idx=argmin(dists)

Функция argmin: Возвращает индекс минимального элемента:

Результат: Индекс дескриптора во втором наборе, наиболее похожего на текущий.

Шаг 5: Сортировка расстояний для Ratio Test

sorted_dists=排序(dists)

分类后:

**重要:**我们只需要比率测试的前两个值。

步骤6:应用比率测试

如果长度(sorted_dists)>=2&&sorted_dists[1]<ratio_thresh*sorted_dists[2]
    推!(匹配,(idx1[i],idx2[best_idx]))
结束

Условие принятия совпадения:

Дополнительная проверка: length(sorted_dists) >= 2 — требуется минимум 2 кандидата для вычисления ratio.

Сохранение соответствия:

(idx1[i],idx2[best_idx])
  • idx1[i] — индекс исходной ключевой точки в первом изображении
  • idx2[best_idx] — индекс соответствующей точки во втором изображении

Оценка гомографии через RANSAC: ransac_homography

Эта функция реализует алгоритм RANSAC (RANdom SAmple Consensus) для робастной оценки матрицы гомографии между двумя изображениями. Это финальный этап пайплайна, который очищает совпадения от выбросов и восстанавливает геометрическое преобразование, связывающее два кадра.

Математическая основа

Что такое гомография?

Гомография (Homography) — это проективное преобразование плоскости, которое описывает соответствие между точками на двух изображениях одной и той же плоской сцены:

Нормализованные координаты:

Степени свободы: 8 (матрица 3×3 с точностью до масштаба, обычно )

Минимальное число точек

Для однозначного определения гомографии требуется минимум 4 пары точек:

Пар точек Уравнений Неизвестных Система
1 пара 2 8 Недоопределена
2 пары 4 8 Недоопределена
3 пары 6 8 Недоопределена
4 пары 8 8 Определена
5+ пар 10+ 8 Переопределена (МНК)

Линейная система для гомографии

Для каждой пары точек получаем два линейных уравнения:

После домножения на знаменатель:

При :

Матричная форма для 4 точек (8 уравнений):

Где:


RANSAC算法

**问题:**在匹配中存在异常值-与真变换不对应的错误比较。

**RANSAC解决方案:**迭代地选择随机点子集,构建模型,并计算一致点(inliers)。

 RANSAC算法:
1. 重复两次:
   A.随机选择4对点(最小集合)
b.从这些点计算单应性H
C.统计inliers的个数(误差<threshold的点)
d.如果inliers大于最佳→保持H和inliers
2. 带回最好的模型

Вероятность успеха RANSAC:

Где:

  • — доля inliers среди всех совпадений
  • — минимальное число точек (4 для гомографии)
  • — количество итераций

Пример расчёта:

Доля inliers (w) Итераций для 99% успеха
50% 163
60% 106
70% 72
80% 50
90% 35

⚙️ Параметры и их влияние

Параметр Значение по умолчанию Эффект увеличения Эффект уменьшения
iter 500 Выше вероятность успеха Быстрее выполнение
thresh 3.0 Больше inliers, меньше точность Меньше inliers, выше точность

Рекомендации по выбору порога:

Задача Рекомендуемый thresh Обоснование
Точная геометрия 1.0–2.0 px Строгая проверка соответствия
Стандартное сопоставление 2.0–3.0 px Баланс точности и полноты
Шумные данные 4.0–5.0 px Учёт погрешности детекции
Панорамная склейка 3.0–4.0 px Допустимы небольшие несовпадения

Формула выбора порога:

Где — точность детектора ключевых точек (обычно 1-2 пикселя).

In [ ]:
function ransac_homography(kp1, kp2, matches; iter=500, thresh=3.0)
    best_inliers = []
    best_H = nothing
    n = length(matches)
    if n < 4
        return matches, nothing
    end
    for _ in 1:iter
        sample_idx = rand(1:n, 4)
        sample = matches[sample_idx]
        A = zeros(8, 8)
        b = zeros(8)
        for (row, (idx1, idx2)) in enumerate(sample)
            x1, y1 = kp1[idx1]
            x2, y2 = kp2[idx2]
            A[2*row-1, :] = [x1, y1, 1, 0, 0, 0, -x2*x1, -x2*y1]
            b[2*row-1] = x2
            A[2*row, :] = [0, 0, 0, x1, y1, 1, -y2*x1, -y2*y1]
            b[2*row] = y2
        end
        try
            h = A \ b
            H = [h[1] h[2] h[3]; h[4] h[5] h[6]; h[7] h[8] 1.0]
            inliers = []
            for (idx1, idx2) in matches
                x1, y1 = kp1[idx1]
                x2, y2 = kp2[idx2]
                p = [x1, y1, 1.0]
                pred = H * p
                pred = pred[1:2] ./ pred[3]
                if norm(pred - [x2, y2]) < thresh
                    push!(inliers, (idx1, idx2))
                end
            end
            if length(inliers) > length(best_inliers)
                best_inliers = inliers
                best_H = H
            end
        catch
            continue
        end
    end
    return best_inliers, best_H
end
Out[0]:
ransac_homography (generic function with 1 method)

Пошаговый разбор кода

Шаг 1: Проверка минимального числа совпадений

n=长度(匹配)
如果n<4
    返回匹配,什么都没有
结束

Зачем: Для гомографии требуется минимум 4 пары точек. Если меньше — система недоопределена.

Шаг 2: Случайная выборка 4 точек

对于_in1:iter
    sample_idx=rand(1:n,4)
    样本=匹配[样本_idx]

Важно: Выборка с возвращением (может выбрать одни и те же точки дважды). В продакшене лучше использовать sample(1:n, 4, replace=false).

Вероятность выбрать 4 inliers:

При (84% inliers в нашем примере): (50% chance per iteration)

Шаг 3: Построение матрицы системы

A=零(8,8)
b=零(8)

for(row,(idx1,idx2))in enumerate(sample)
    x1,y1=kp1[idx1]
    x2,y2=kp2[idx2]
    
    #X2的方程
    A[2*row-1,:]=[x1,y1, 1, 0, 0, 0, -x2*x1,-x2*y1]
    b[2*row-1]=x2
    
    #Y2的方程
    A[2*行, :] = [0, 0, 0, x1,y1,1,-y2*x1,-y2*y1]
    b[2*行]=y2
结束

Структура матрицы A (8×8):

Строка
1 (x₁) x₁ y₁ 1 0 0 0 -x₂x₁ -x₂y₁
2 (y₁) 0 0 0 x₁ y₁ 1 -y₂x₁ -y₂y₁
3 (x₂) x₁ y₁ 1 0 0 0 -x₂x₁ -x₂y₁
4 (y₂) 0 0 0 x₁ y₁ 1 -y₂x₁ -y₂y₁
... ... ... ... ... ... ... ... ...

Шаг 4: Решение системы уравнений

试试
    h=A\b
    H=[h[1]h[2]h[3];h[4]h[5]h[6];h[7]h[8]1.0]

Оператор \ в Julia: Решает систему методом наименьших квадратов (для переопределённых систем) или точным решением (для квадратных).

Восстановление матрицы 3×3:

Обработка ошибок: try-catch защищает от вырожденных матриц (например, когда 4 точки коллинеарны).

Шаг 5: Подсчёт inliers

inliers=[]
对于(idx1,idx2)在匹配
    x1,y1=kp1[idx1]
    x2,y2=kp2[idx2]
    p=[x1,y1,1.0]
    pred=H*p
    pred=pred[1:2]。/pred[3]
    
    如果norm(pred-[x2,y2])<thresh
        推!(inliers,(idx1,idx2))
    结束
结束

Проективное преобразование точки:

Ошибка репроекции (Reprojection Error):

Критерий inlier:

Шаг 6: Обновление лучшей модели

如果长度(inliers)>长度(best_inliers)
    best_inliers=inliers
    best_H=H
结束

Критерий качества: Максимальное количество inliers (не минимальная ошибка!).

Почему: RANSAC оптимизирует consensus (согласие), а не точность модели. Точность можно улучшить потом через МНК по всем inliers.

Шаг 7: Возврат результата

返回best_inliers,best_H
```</span>

**输出数据:**
- `best_inliers` -与模型一致的匹配索引数组
- `best_H` -估计3×3单应矩阵


所有阶段的启动和可视化

这部分代码是对整个计算机视觉管道的集成测试。 在这里,所有六个阶段都按顺序执行:从图像加载到通过RANSAC进行最终单应性评估。 每个步骤都伴随着统计输出和结果的可视化,这使您可以直观地监控每个阶段的算法质量和数据丢失。

,这段代码演示了什么

/阶段/输入数据/输出数据|关键度量/
|------|---------------|-----------------|------------------|
| 1. 上传/图片文件|512×512数组|—/
| 2. 变换/原始图像/旋转图像/角度:π/15(12°)|
| *3. 检测||2图像|421→297点|损失:29%/
| 4. 描述符|297分|227有效|损失:24%|
| *5. 匹配|/227×227描述符|73对/精度:~84%|
| 6. RANSAC/73匹配/61更早|准确性:95%+|

**共同目标:**展示如何从421个最初检测到的点获得61个可靠匹配,足以准确地重建图像之间的几何变换。

In [ ]:
println("\N[1]上传图片。..")
img_original = testimage("lena_gray_512")
img_original = Float64.(img_original)  
display(Gray.(img_original))
[1] Загрузка изображения...
Warning: Usage of "lena" is not recommended, and the image may be removed in a later release. See https://womenlovetech.com/losing-lena-why-we-need-to-remove-one-image-and-end-techs-original-sin/ for more information.
@ TestImages /usr/local/ijulia-demos/packages/TestImages/UtVkc/src/TestImages.jl:124
No description has been provided for this image
In [ ]:
println("\n[2]变换(旋转+缩放)。..")
angle = π/15
scale = 1.0
center = SVector(size(img_original, 2)/2, size(img_original, 1)/2)
tfm = Translation(center)  LinearMap(RotMatrix(angle))  LinearMap(SDiagonal(scale, scale))  Translation(-center)

img_transformed = warp(img_original, tfm, axes(img_original))
img_transformed = replace(img_transformed, NaN => 0.0, Inf => 0.0, -Inf => 0.0)
img_transformed_gray = Gray.(clamp01.(img_transformed))
display(img_transformed_gray)
[2] Трансформация (поворот + масштаб)...
No description has been provided for this image
In [ ]:
println("\n[3]关键点检测(Harris角点)。..")
keypoints_original = harris_corners(img_original)
keypoints_transformed = harris_corners(img_transformed)
keypoints_transformed = filter_valid_keypoints(img_transformed, keypoints_transformed)

println("   原文上的要点:$(length(keypoints_original))")
println("   转换(有效)上的关键点:$(长度(keypoints_transformed))")

img_kp_original = RGB.(img_original)
for (x,y) in keypoints_original
    draw_circle!(img_kp_original, (x, y), 3, RGB(1,0,0))
end
display(img_kp_original)

img_kp_transformed = RGB.(img_transformed)
for (x,y) in keypoints_transformed
    draw_circle!(img_kp_transformed, (x, y), 3, RGB(0,0,1))
end
display(img_kp_transformed)
[3] Детекция ключевых точек (Harris corners)...
   Ключевых точек на исходном: 421
   Ключевых точек на трансформированном (валидных): 297
No description has been provided for this image
No description has been provided for this image
In [ ]:
println("\n[4]提取描述符(基于补丁)。..")
desc_orig, valid_idx_orig = extract_patch_descriptors(img_original, keypoints_original)
desc_trans, valid_idx_trans = extract_patch_descriptors(img_transformed, keypoints_transformed)

println("   原文中的有效描述符:$(length(valid_idx_orig))")
println("   转换后的有效描述符:$(length(valid_idx_trans))")
[4] Извлечение дескрипторов (patch-based)...
   Валидных дескрипторов на исходном: 421
   Валидных дескрипторов на трансформированном: 227
In [ ]:
println("\n[5]用比率检验匹配描述符。..")
matches = match_descriptors_euclidean(desc_orig, desc_trans, valid_idx_orig, valid_idx_trans, ratio_thresh=0.7)
println("   比率测试后的匹配:$(长度(匹配))")

h, w = size(img_original)
combined = fill(Gray(0), h, 2w)
combined[:, 1:w] = Gray.(img_original)
combined[:, w+1:2w] = img_transformed_gray
combined_rgb = RGB.(combined)

for (i1, i2) in matches
    x1, y1 = keypoints_original[i1]
    x2, y2 = keypoints_transformed[i2]
    p1 = (x1, y1)
    p2 = (x2 + w, y2)  
    draw_line!(combined_rgb, p1, p2, RGB(0,1,0))
end
display(combined_rgb)
[5] Сопоставление дескрипторов с ratio test...
   Совпадений после ratio test: 73
No description has been provided for this image
In [ ]:
println("\n[6] RANSAC...")
inliers, H = ransac_homography(keypoints_original, keypoints_transformed, matches, iter=1000, thresh=3.0)
println("   RANSAC后的匹配项(inliers):$(length(inliers))")

combined_ransac = RGB.(combined)
for (i1, i2) in inliers
    x1, y1 = keypoints_original[i1]
    x2, y2 = keypoints_transformed[i2]
    p1 = (x1, y1)
    p2 = (x2 + w, y2)
    draw_line!(combined_ransac, p1, p2, RGB(0,1,0))
    draw_circle!(combined_ransac, p1, 5, RGB(1,0,0))
    draw_circle!(combined_ransac, p2, 5, RGB(1,0,0))
end
display(combined_ransac)
[6] RANSAC...
   Совпадений после RANSAC (inliers): 61
No description has been provided for this image

本文以Julia语言呈现一个完整的计算机视觉训练示例,它实现了从零开始检测和匹配关键点的经典管道。 与使用现成的库(OpenCV,VLFeat)不同,每个算法都是手动实现的,并详细解释了数学基础。

实现的组件:

/组件/算法/代码行/复杂性|
|-----------|----------|------------|-----------|
/渲染/画布(线,圆)/~60| |
/检测/哈里斯-史蒂文斯/~40| |
/描述符/基于补丁的+Z-分数/~35| |
/匹配/欧几里德+比率测试/~20| |
/模型的评估/RANSAC+单应性|~50| |

总计:~200行代码,演示从像素到几何模型的整个周期。

结论

这段代码表明,理解计算机视觉的基本算法比使用现成库的能力更重要。
使用此代码作为您自己实验的基础-更改参数,尝试不同的图像,添加尺度不变性,或与现成的描述符集成。 在这里获得的理解将成为使用任何计算机视觉系统的基础。