带显式仿射项FastVFC

Posted on 2026-05-20 In Tech , Vision Views: Disqus:

一、统一符号定义（维度 + 含义 + 代码对应）

数学符号	矩阵维度	物理含义	代码变量
\(N\)	标量	匹配点对总数	`N = Y.rows()`
\(M\)	标量	核矩阵低秩近似主成分数	`M = Q.cols()`
\(D=2\)	标量	二维平面坐标	固定值
\(X\)	\(\mathbb{R}^{N \times 2}\)	归一化参考点坐标 \((x_i,y_i)\)	入参 `const Eigen::MatrixXf& X`
\(Y\)	\(\mathbb{R}^{N \times 2}\)	观测位移向量场	入参 `const Eigen::MatrixXf& Y`
\(P\)	\(\mathbb{R}^{N}\)	内点后验概率向量（对角矩阵 \(\mathrm{diag}(P)\)）	入参 `const Eigen::VectorXf& P`
\(Q\)	\(\mathbb{R}^{N \times M}\)	高斯核低秩特征向量矩阵	入参 `const Eigen::MatrixXf& Q`
\(S\)	\(\mathbb{R}^{M}\)	高斯核低秩特征值向量	入参 `const Eigen::VectorXf& S`
\(C\)	\(\mathbb{R}^{N \times 2}\)	输出：局部非刚性核系数	出参 `Eigen::MatrixXf& C`
\(\Theta\)	\(\mathbb{R}^{3 \times 2}\)	输出：全局仿射参数 \([平移, A_x, A_y]^T\)	出参 `Eigen::MatrixXf& Affine`
\(\lambda\)	标量	正则化系数	入参 `float lambda`
\(\sigma^2\)	标量	噪声方差	入参 `float sigma2`
\(\lambda_s\)	标量	联合正则系数 \(\lambda \cdot \sigma^2\)	`lambda_sigma2`
\(P_{aff}\)	\(\mathbb{R}^{N \times 3}\)	仿射基矩阵 \([\mathbf{1}, X]\)	`P_affine`

二、核心数学模型

将 VFC 向量场解耦为全局仿射变换 + 局部非刚性形变，同时加入正交约束（强制非刚性分量不含全局仿射趋势，与 TPS 数学同源）： ### 1. 向量场定义 \[ f(x) = \underbrace{P_{aff} \Theta}_{全局仿射} + \underbrace{K C}_{局部非刚性} \] ### 2. 扩展线性系统（VFC M 步核心） \[ \begin{cases} \left(K + \lambda \sigma^2 P^{-1}\right) C + P_{aff} \Theta = Y \quad (1) \quad // 拟合方程 \\ P_{aff}^T C = 0 \quad\quad\quad\quad\quad\quad\quad\quad\quad\ (2) \quad // 正交约束：非刚性系数群仿射力矩为0 \end{cases} \]

三、分步数学推导（核心修正：正交约束代入）

步骤 1：方程等价变形（消去拟合逆矩阵）

对公式 (1) 左乘概率对角矩阵 \(P\)，化简得： \[ \left(P K + \lambda \sigma^2 I_N\right) C + P P_{aff} \Theta = P Y \tag{3} \]

步骤 2：代入核矩阵低秩近似

高斯核低秩分解：\(K \approx Q S Q^T\)，代入公式 (3)： \[ \left(P Q S Q^T + \lambda \sigma^2 I_N\right) C + P P_{aff} \Theta = P Y \tag{4} \]

步骤 3：定义辅助变量（FastVFC 降维核心）

定义低维辅助变量 \(T = S Q^T C \in \mathbb{R}^{M \times 2}\)，代入公式 (4) 并移项，直接得到核系数闭式解： \[ \boldsymbol{C = \frac{1}{\lambda \sigma^2} \Big( P Y - P Q T - P P_{aff} \Theta \Big)} \tag{5} \]

步骤 4：正定子系统的反向构建

将显式关系 (5) 代回物理性质方程，构建极小系统：

(A) 代回辅助变量定义式 \(T = S Q^T C\)： \[ T = S Q^T \left[ \frac{1}{\lambda \sigma^2} \Big( P Y - P Q T - P P_{aff} \Theta \Big) \right] \] 整理化简可得第一方程： \[ \left( Q^T P Q + \lambda \sigma^2 S^{-1} \right) T + \left( Q^T P P_{aff} \right) \Theta = Q^T P Y \tag{6} \]

(B) [严格修正] 代回正交约束 \(P_{aff}^T C = 0\)： \[ P_{aff}^T \left[ \frac{1}{\lambda \sigma^2} \Big( P Y - P Q T - P P_{aff} \Theta \Big) \right] = 0 \] 消去 \(\frac{1}{\lambda \sigma^2}\)，并将待求解变量移项。此时伴生在 \(\Theta\) 上的系数矩阵浮现： \[ \left( P_{aff}^T P Q \right) T + \boldsymbol{\left( P_{aff}^T P P_{aff} \right)} \Theta = P_{aff}^T P Y \tag{7} \]

步骤 5：构建低维对称正定系统 (SPD)

将 (6) 与 (7) 联立为 \((M+3) \times (M+3)\) 块矩阵：

\[ \underbrace{ \begin{bmatrix} Q^T P Q + \lambda \sigma^2 S^{-1} & Q^T P P_{aff} \\ P_{aff}^T P Q & \boldsymbol{P_{aff}^T P P_{aff}} \end{bmatrix}}_{K_{ext}} \underbrace{\begin{bmatrix} T \\ \Theta \end{bmatrix}}_{T_{ext}} = \underbrace{\begin{bmatrix} Q^T P Y \\ P_{aff}^T P Y \end{bmatrix}}_{B_{ext}} \tag{8} \]

左侧系数大矩阵 \(K_{ext}\) 是一个绝对的对称正定矩阵 (SPD)。由于它是 SPD 矩阵，我们可以安全地只计算下三角部分并使用最高效的 LDLT 分解求逆，彻底避免鞍点系统带来的数值崩溃。

四、C++ 实现代码

利用上一步论证得出的 SPD 对称性，剔除多余乘法、转置操作：

/**
 * @brief 带全局仿射项的VFC控制点更新 (数学严格版，利用 SPD 半区域加速)
 * @param X 归一化参考点坐标 N×2
 * @param Y 观测位移向量场 N×2
 * @param P 内点概率向量 N×1
 * @param Q 核低秩特征向量 N×M
 * @param S 核低秩特征值 M×1
 * @param C 输出：非刚性核系数 N×2
 * @param Affine 输出：仿射参数 3×2 [平移, x线性系数, y线性系数]
 * @param sigma2 噪声方差
 * @param lambda 正则化系数
 * @return 求解是否成功
 */
static bool UpdateControlPointsWithAffine(const Eigen::MatrixXf& X,
                                          const Eigen::MatrixXf& Y,
                                          const Eigen::VectorXf& P,
                                          const Eigen::MatrixXf& Q,
                                          const Eigen::VectorXf& S,
                                          Eigen::MatrixXf& C,
                                          Eigen::MatrixXf& Affine,
                                          float sigma2,
                                          float lambda) {
    const int N = Y.rows();
    const int M = Q.cols();
    const float eps = 1e-6f;
    const float lambda_sigma2 = lambda * sigma2;

    // ====================== 1. 构建仿射基矩阵 P_aff = [1, x, y] (N×3) ======================
    Eigen::MatrixXf P_affine(N, 3);
    P_affine.col(0).setOnes();
    P_affine.col(1) = X.col(0);
    P_affine.col(2) = X.col(1);

    // ====================== 2. 概率加权矩阵 (N 维度运算) ======================
    Eigen::MatrixXf PQ(N, M); PQ.noalias() = P.asDiagonal() * Q;
    Eigen::MatrixXf PY(N, 2); PY.noalias() = P.asDiagonal() * Y;
    Eigen::MatrixXf PP(N, 3); PP.noalias() = P.asDiagonal() * P_affine;

    // ====================== 3. 构建低维系统子矩阵 (极限加速区) ======================
    // 因为最终大矩阵是对称正定的，大块矩阵中针对自身的内积我们只算 Lower 下半角，减少近一半计算量
    Eigen::MatrixXf QtPQ(M, M);
    QtPQ.triangularView<Eigen::Lower>() = Q.transpose() * PQ;
    Eigen::MatrixXf PtPP(3, 3);
    PtPP.triangularView<Eigen::Lower>() = P_affine.transpose() * PP; 
    
    Eigen::MatrixXf QtPP(M, 3); QtPP.noalias() = Q.transpose() * PP;
    
    // 生成右侧目标 B_ext 块
    Eigen::MatrixXf QtPY(M, 2); QtPY.noalias() = Q.transpose() * PY;
    Eigen::MatrixXf PtPY(3, 2); PtPY.noalias() = P_affine.transpose() * PY;

    // ====================== 4. 添加正则化项 ======================
    Eigen::VectorXf inv_S = S.cwiseMax(eps).cwiseInverse();
    QtPQ.diagonal() += lambda_sigma2 * inv_S;

    // ====================== 5. 拼装分块矩阵 K_ext (M+3 × M+3) ======================
    Eigen::MatrixXf K_ext(M + 3, M + 3);
    // 因为 LDLT 使用 Lower 模式，我们只需填充下三角区域
    K_ext.topLeftCorner(M, M).triangularView<Eigen::Lower>() = QtPQ; 
    K_ext.bottomRightCorner(3, 3).triangularView<Eigen::Lower>() = PtPP;
    K_ext.bottomLeftCorner(3, M) = QtPP.transpose();   

    Eigen::MatrixXf B_ext(M + 3, 2);
    B_ext.topRows(M) = QtPY;
    B_ext.bottomRows(3) = PtPY;

    // ====================== 6. 下三角 LDLT 分解求解 ======================
    Eigen::LDLT<Eigen::MatrixXf> ldlt;
    ldlt.compute(K_ext.selfadjointView<Eigen::Lower>());
    if (ldlt.info() != Eigen::Success) return false;
    
    Eigen::MatrixXf T_ext(M + 3, 2);
    T_ext.noalias() = ldlt.solve(B_ext);

    // ====================== 7. 还原最终系数 ======================
    Eigen::MatrixXf T = T_ext.topRows(M);
    Affine = T_ext.bottomRows(3);
    C.noalias() = (PY - PQ * T - PP * Affine) / lambda_sigma2;

    return true;
}

配套的向量场合成函数（严禁使用隐式会导致失真的 .transpose()）：

static void UpdateVectorFieldWithAffine(const Eigen::MatrixXf& Q,
                                        const Eigen::VectorXf& S,
                                        const Eigen::MatrixXf& C,
                                        const Eigen::MatrixXf& Affine,
                                        const Eigen::MatrixXf& X,
                                        Eigen::MatrixXf& V) {
    // 1. 局部非刚性形变
    V.noalias() = Q * (S.asDiagonal() * (Q.transpose() * C));
    // 2. 全局仿射形变 (严格坐标方向映射相加)
    V.rowwise() += Affine.row(0);            // 广播相加平移量
    V.noalias() += X * Affine.bottomRows(2); // x,y 线性矩阵点乘，输出直接为 N×2 尺寸坐标系
}

五、实验结果

VFC NORMAL	VFC FAST	VFC FAST AFFINE

可以看到，FastVFC加上带显式仿射项后，并没有使得算法快速收敛。