Fӧrstner特征提取算法

### Förstner 特征提取算法概述

Förstner特征提取是一种基于局部灰度变化分析的方法，在计算机视觉领域广泛应用。此方法不仅考虑了图像平面内的噪声特性，还通过点的观测将相机坐标系中的协方差传递至下一个坐标系中[^1]。

#### 原理

核心在于构建并求解所谓的“误差椭圆”，即对于每一个候选角点位置，估计其作为实际角点的可能性大小。具体来说：

- **二阶矩矩阵**：类似于Harris算子，Förstner同样依赖于计算图像梯度及其外积形成的结构张量（Structure Tensor）。这一过程涉及到两个尺度参数的选择——一个是用于抑制高频噪音和平滑数据的空间尺度；另一个则是用来聚合局部信息的时间/空间尺度[^2]。

- **不确定性评估**：不同于其他仅返回响应强度值的传统角点探测器，Förstner进一步量化了每个潜在兴趣点处测量不确定性的程度。这使得即使是在存在显著成像畸变的情况下也能获得较为可靠的定位精度。

#### 实现细节

以下是简化版Python代码片段展示如何实现基本功能：

import numpy as np
from scipy.ndimage import gaussian_filter, sobel

def compute_structure_tensor(image, sigma_space=0.7, sigma_time=1.5):
    """Compute the structure tensor at each pixel."""
    
    Ix = sobel(image, axis=1, mode='constant')
    Iy = sobel(image, axis=0, mode='constant')

    # Smooth derivatives with Gaussian kernel
    Ixx = gaussian_filter(Ix * Ix, sigma=sigma_space)
    Iyy = gaussian_filter(Iy * Iy, sigma=sigma_space)
    Ixy = gaussian_filter(Ix * Iy, sigma=sigma_space)

    return (Ixx, Iyy, Ixy)


def forstner_detector(image, mask_size=3, min_grad=0.1, min_apparentness=0.01):
    """Detect corners using Forstner's method."""

    rows, cols = image.shape[:2]

    # Compute elements of the structure tensor field
    A_xx, A_yy, A_xy = compute_structure_tensor(image)

    # Initialize output arrays
    certainty_map = np.zeros_like(image, dtype=float)
    row_coords = []
    col_coords = []

    half_mask = int(mask_size / 2)

    for r in range(half_mask, rows-half_mask):
        for c in range(half_mask, cols-half_mask):

            # Extract local neighborhood values from structure tensors
            M_ii = A_xx[r,c]
            M_ij = A_xy[r,c]
            M_jj = A_yy[r,c]

            detM = M_ii*M_jj - M_ij**2
            traceM = M_ii + M_jj
            
            if detM > 0 and traceM != 0:

                w = detM/(traceM*traceM)  # Certainty measure
                
                if w >= min_apparentness:
                    certainty_map[r][c] = w
                    
                    if abs(M_ii-M_jj)*np.sqrt(detM)/abs(traceM)**(3./2.) >= min_grad:
                        row_coords.append(r)
                        col_coords.append(c)

    return certainty_map, row_coords, col_coords