光流的hession矩阵 一些记录
上面的直观解释:
Hessian矩阵的特征值:就是形容其在该点附近特征向量方向的凹凸性,特征值越大,凸性越强。
对于二维图像的某点的hessian矩阵,其最大特征值和其对应的特征向量对应其邻域二维曲线最大曲率的强度和方向,即山坡陡的那面;最小特征值对应的特征向量对应与其垂直的方向,即平缓的方向。
简单来讲,图像某点的hessian矩阵特征值大小和符号决定了该点邻域内的几何结构
三维图像同理。由此原理设计的一些滤波器(最经典的是Frangi filter)在医学领域的MRA和CTA血管造影增强有很重要的应用。
基于Hessian矩阵,就可以判断多元函数的极值情况了,结论如下
如果是正定矩阵,则临界点处是一个局部极小值
如果是负定矩阵,则临界点处是一个局部极大值
如果是不定矩阵,则临界点处不是极值
对于实二次型矩阵正定负定的判定方法:实二次型矩阵为正定二次型的充要条件是的矩阵的特征值全大于零。为负定二次型的充要条件是的矩阵的特征值全小于零,否则是不定的。
1.矩阵的特征值之和等于矩阵的行列式
2.矩阵的特征值之积等于矩阵的迹
计算最小特征值的方法就用到上面的东西了。公式min_lambda + max_lambda = dttermiant; min_lambda * max_lambda = trace;用二次方程的求解公式很容易求得min_lambda就是下面的公式:min_lambda= trace / 2 - sqrt((trace * trace) / 4 - determinant);
float trace = hession(0, 0) + hession(1, 1);
float determinant = hession(0, 0) * hession(1, 1) - hession(0, 1) * hession(1, 0);
float min_lambda= trace / 2 - sqrt((trace * trace) / 4 - determinant);
result(0, 0) = determinant;
result(1, 0) = min_lambda;
光流中用到的具体代码位置:
hession计算特征值部分:
inline void CacluHessen(const MatrixXf &dx, const MatrixXf &dy, int x, int y,
MatrixXf &result, MatrixXf &hession) {
MatrixXf dx_square = MatrixXf::Zero(WinSize, WinSize);
MatrixXf dxy_square = MatrixXf::Zero(WinSize, WinSize);
MatrixXf dy_square = MatrixXf::Zero(WinSize, WinSize);
dx_square = dx.block(y - halfSize, x - halfSize, WinSize, WinSize).array() *
dx.block(y - halfSize, x - halfSize, WinSize, WinSize).array();
dy_square = dy.block(y - halfSize, x - halfSize, WinSize, WinSize).array() *
dy.block(y - halfSize, x - halfSize, WinSize, WinSize).array();
dxy_square = dx.block(y - halfSize, x - halfSize, WinSize, WinSize).array() *
dy.block(y - halfSize, x - halfSize, WinSize, WinSize).array();
hession(0, 0) = dx_square.sum();
hession(1, 1) = dy_square.sum();
hession(0, 1) = dxy_square.sum();
hession(1, 0) = dxy_square.sum();
float trace = hession(0, 0) + hession(1, 1);
float determinant = hession(0, 0) * hession(1, 1) - hession(0, 1) * hession(1, 0);
float min_lambda= trace / 2 - sqrt((trace * trace) / 4 - determinant);
result(0, 0) = determinant;
result(1, 0) = min_lambda;
}
应用部分:
if(result(1, 0) < LK_minEigThreshold || result(0, 0) < FLT_EPSILON)
{
if(0 == level)
{
status(i, 0) = 0;
}
continue;
}