http://www.cnblogs.com/ronny/p/3895883.html
http://www.cnblogs.com/ronny/p/4045979.html
http://www.cnblogs.com/ronny/p/4048213.html
1. 什么是斑点
斑点通常是指与周围有着颜色和灰度差别的区域。在实际地图中,往往存在着大量这样的斑点,如一颗树是一个斑点,一块草地是一个斑点,一栋房子也可以是一个斑点。由于斑点代表的是一个区域,相比单纯的角点,它的稳定性要好,抗噪声能力要强,所以它在图像配准上扮演了很重要的角色。
同时有时图像中的斑点也是我们关心的区域,比如在医学与生物领域,我们需要从一些X光照片或细胞显微照片中提取一些具有特殊意义的斑点的位置或数量。
比如下图中天空的飞机、向日葵的花盘、X线断层图像中的两个斑点。
在视觉领域,斑点检测的主要思路都是检测出图像中比它周围像素灰度值大或比周围灰度值小的区域。一般有两种方法来实现这一目标:
- 基于求导的微分方法,这类的方法称为微分检测器;
- 基于局部极值的分水岭算法。
这里我们重点介绍第一种方法,主要检测LOG斑点。而OpenCV中SimpleBlobDetector斑点检测算子就实现了第二种方法,我们这里也会介绍它的接口使用方法。
2. LOG斑点检测
2.1 基本原理
利用高斯拉普通拉斯(Laplace of Gaussian,LOG)算子检测图像斑点是一种十分常用的方法,对于二维高斯函数:
2. LOG斑点检测
2.1 基本原理
利用高斯拉普通拉斯(Laplace of Gaussian,LOG)算子检测图像斑点是一种十分常用的方法,对于二维高斯函数:
它的拉普拉斯变换为:
规范化的高斯拉普变换为:
规范化算法子在二维图像上显示是一个圆对称函数,如下图所示。我们可以用这个算子来检测图像中的斑点,并且可以通过改变σ的值,可以检测不同尺寸的二维斑点。
2.2 LOG原理解释
其实从更直观的角度去解释为什么LOG算子可以检测图像中的斑点是:
图像与某一个二维函数进行卷积运算实际就是求取图像与这一函数的相似性。同理,图像与高斯拉普拉斯函数的卷积实际就是求取图像与高斯拉普拉斯函数的相似性。当图像中的斑点尺寸与高斯拉普拉斯函数的形状趋近一致时,图像的拉普拉斯响应达到最大。
从概率的角度解释为:假设原图像是一个与位置有关的随机变量X的密度函数,而LOG为随机变量Y的密度函数,则随机变量X+Y的密度分布函数即为两个函数的卷积形式(这一部分的理论,可以参见本博客概率与统计相关文章)。如果想让X+Y能取到最大值,则X与Y能保持步调一致最好,即X上升时,Y也上升,X最大时,Y也最大。
那么LOG算子是怎么被构想出来的呢?
事实上我们知道Laplace可以用来检测图像中的局部极值点,但是对噪声敏感,所以在我们对图像进行Laplace卷积之前,我们用一个高斯低通滤波对图像进行卷积,目标是去除图像中的噪声点。这一过程 可以描述为:
先对图像f(x,y)
用方差为σ的高斯核进行高斯滤波,去除图像中的噪点。
然后对图像的拉普拉斯图像则为:
而实际上有下面等式:
所以,我们可以先求高斯核的拉普拉斯算子,再对图像进行卷积。也就是一开始描述的步骤。
2.3 LOG算子的实现
Mat Feat::getHOGKernel(Size& ksize, double sigma)
{
Mat kernel(ksize, CV_64F);
Point centPoint = Point((ksize.width -1)/2, ((ksize.height -1)/2));
// first calculate Gaussian
for (int i=0; i < kernel.rows; i++)
{
double* pData = kernel.ptr<double>(i);
for (int j = 0; j < kernel.cols; j++)
{
double param = -((i - centPoint.y) * (i - centPoint.y) + (j - centPoint.x) * (j - centPoint.x)) / (2*sigma*sigma);
pData[j] = exp(param);
}
}
double maxValue;
minMaxLoc(kernel, NULL, &maxValue);
for (int i=0; i < kernel.rows; i++)
{
double* pData = kernel.ptr<double>(i);
for (int j = 0; j < kernel.cols; j++)
{
if (pData[j] < EPS* maxValue)
{
pData[j] = 0;
}
}
}
double sumKernel = sum(kernel)[0];
if (sumKernel != 0)
{
kernel = kernel / sumKernel;
}
// now calculate Laplacian
for (int i=0; i < kernel.rows; i++)
{
double* pData = kernel.ptr<double>(i);
for (int j = 0; j < kernel.cols; j++)
{
double addition = ((i - centPoint.y) * (i - centPoint.y) + (j - centPoint.x) * (j - centPoint.x) - 2*sigma*sigma)/(sigma*sigma*sigma*sigma);
pData[j] *= addition;
}
}
// make the filter sum to zero
sumKernel = sum(kernel)[0];
kernel -= (sumKernel/(ksize.width * ksize.height));
return kernel;
}
2.4 多尺度检测
我们注意到当σ尺度一定时,只能检测对应半径的斑点,那么检测的是多大半径的斑点呢,我们可以通过对规范化的二维拉普拉斯高斯算子求导:
规范化的高斯拉普拉斯函数为:
求的极点值等价于求取下式:
得到:
对于图像中的斑点,在尺度σ=r/√2时,高斯拉普拉斯响应值达到最大。同理,如果图像中的圆形斑点黑白反向,那么,它的高斯拉普拉斯响应值在σ=r/√2时达到最小。将高斯拉普拉斯响应达到峰值时的尺度σ值,称为特征尺度。
那么在多尺度的情况下,同时在空间和尺度上达到最大值(或最小值)的点就是我们所期望的斑点。对于二维图像I(x,y)
,计算图像在不同尺度下的离散拉普拉斯响应值,然后检查位置空间中的每个点;如果该点的拉普拉斯响应值都大小于或小于其他26个立方空间领域(9+8+9)的值,那么该点就是被检测到的图像斑点。
3 OpenCV进行斑点检测
opencv中检测Blobs的类为SimpleBlobDetector,这个类在opencv中的定义如下:
class SimpleBlobDetector : public FeatureDetector
{
public:
struct Params
{
Params();
float thresholdStep;
float minThreshold;
float maxThreshold;
size_t minRepeatability;
float minDistBetweenBlobs;
bool filterByColor;
uchar blobColor;
bool filterByArea;
float minArea, maxArea;
bool filterByCircularity;
float minCircularity, maxCircularity;
bool filterByInertia;
float minInertiaRatio, maxInertiaRatio;
bool filterByConvexity;
float minConvexity, maxConvexity;
};
SimpleBlobDetector(const SimpleBlobDetector::Params ¶meters = SimpleBlobDetector::Params());
protected:
...
};
算法的大致步骤如下:
- 对[minThreshold,maxThreshold)区间,以thresholdStep为间隔,做多次二值化。
- 对每张二值图片,使用findContours()提取连通域并计算每一个连通域的中心。
- 根据2得到的中心,全部放在一起。一些很接近的点[由theminDistBetweenBlobs控制多少才算接近]被归为一个group,对应一个bolb特征..
- 从3得到的那些点,估计最后的blob特征和相应半径,并以key points返回。
同时该支持提取特征的方法,一共有5个选项,这里就不多加描述了,默认是提取黑色圆形的Blob特征。下面是一个示例
int main(int argc, char** argv)
{
Mat image = imread(argv[1]);
vector<KeyPoint> keyPoints;
SimpleBlobDetector::Params params;
SimpleBlobDetector blobDetect(params);
blobDetect.create("SimpleBlob");
blobDetect.detect(image, keyPoints);
cout << keyPoints.size() << endl;
drawKeypoints(image, keyPoints, image, Scalar(255,0,0));
namedWindow("blobs");
imshow("blobs", image);
waitKey();
return 0;
}
总体来说,OpenCV的斑点检测效果还算不错,但是在有些图像的效果上明显不如LOG算子检测的检测效果。
4. 扩展阅读
一个与LOG滤波核近似的是高斯差分DOG滤波核,它的定义为:
其中k
为两个相邻尺度间的比例因子。
DOG可以看作为LOG的一个近似,但是它比LOG的效率更高。
前面介绍的微分算子在近圆的斑点检测方面效果很好,但是这些检测算子被限定于只能检测圆形斑点,而且不能估计斑点的方向,因为LOG算子等都是中心对称的。如果我们定义一种二维高斯核的变形,记它在X方向与Y方向上具有不同的方差,则这种算子可以用来检测带有方向的斑点。
其中A是规一性因子。
5. 参考资料
1. 《现代数字图像 -- 处理技术提高与应用案例详解》
2. 《图像局部不变性特征与描述》
3. Lindeberg, T. Feature Detection with Automatic Scale Selection
4. Hui Kong. A Generalized Laplacian Of Gaussian Filter for Blob Detection and Its Applications.
5. OpenCV2马拉松第20圈——blob特征检测原理与实现
6. 积分图像
SURF算法中要用到积分图像的概念。借助积分图像,图像与高斯二阶微分模板的滤波转化为对积分图像的加减运算。
积分图像中任意一点(i,j)的值ii(i,j),为原图像左上角到点(i,j)相应的对角线区域灰度值的总和.
OpenCV中提供了用于计算积分图像的接口
/*
* src :输入图像,大小为M*N
* sum: 输出的积分图像,大小为(M+1)*(N+1)
* sdepth:用于指定sum的类型,-1表示与src类型一致
*/
void integral(InputArray src, OutputArray sum, int sdepth = -1);
值得注意的是OpenCV里的积分图大小比原图像少一行一列,那是因为OpenCV中积分图的计算公式为:
7. DoH近似
前述斑点检测已经提到过,我们可以利用Hessian矩阵行列式的极大值检测斑点。下面我们给出Hessian矩阵的定义。
给定图像I中的一个点x(i,j),在点x处,尺度为σ的Hessian矩阵H(x,σ)定义如下:
式中,Lxx(x,σ)是高斯二阶微分在点x处与图像I的卷积,Lx,y(x,σ)和Lyy(x,σ)具有类似的含义。
下面显示的是上面三种高斯微分算子的图形。
但是利用Hessian行列式进行图像斑点检测时,有一个缺点。由于二阶高斯微分被离散化和裁剪的原因,导致了图像在旋转奇数倍的π/4时,即转换到模板的对角线方向时,特征点检测的重复性降低(也就是说,原来特征点的地方,可能检测不到特征点了)。而在π/2时,特征点检测的重现率真最高。但这一小小的不足不影响我们使用Hessian矩阵进行特征点的检测。
为了将模板与图产像的卷积转换为盒子滤波运算,我们需要对高斯二阶微分模板进行简化,使得简化后的模板只是由几个矩形区域组成,矩形区域内填充同一值,如下图所示,在简化模板中白色区域的值为正数,黑色区域的值为负数,灰度区域的值为0。
对于σ=1.2的高斯二阶微分滤波器,我们设定模板的尺寸为9×9的大小,并用它作为最小尺度空间值对图像进行滤波和斑点检测。我们使用Dxx、Dxy和Dyy表示模板与图像进行卷积的结果。这样,便可以将Hessian矩阵的行列式作如下的简化。
滤波器响应的相关权重w是为了平衡Hessian行列式的表示式。这是为了保持高斯核与近似高斯核的一致性。
其中|…|F为Frobenius范数。理论上来说对于不同的σ的值和对应尺寸的模板尺寸,w值是不同的,但为了简化起见,可以认为它是同一个常数。
使用近似的Hessian矩阵行列式来表示图像中某一点x
处的斑点响应值,遍历图像中所有的像元点,便形成了在某一尺度下琉璃点检测的响应图像。使用不同的模板尺寸,便形成了多尺度斑点响应的金字塔图像,利用这一金字塔图像,就可以进行斑点响应极值点的搜索,其过程完全与SIFT算法类同。
8. 尺度空间表示
通常想要获取不同尺度的斑点,必须建立图像的尺度空间金字塔。一般的方法是通过不同σ
的高斯函数,对图像进行平滑滤波,然后重采样图像以获得更高一层的金字塔图像。SIFT特征检测算法中就是通过相邻两层图像金字塔相减得到DoG图像,然后再在DoG图像上进行斑点和边缘检测工作的。
由于采用了盒子滤波和积分图像,所以,我们并不需要像SIFT算法那样去直接建立图像金字塔,而是采用不断增大盒子滤波模板的尺寸的间接方法。通过不同尺寸盒子滤波模板与积分图像求取Hessian矩阵行列式的响应图像。然后在响应图像上采用3D非最大值抑制,求取各种不同尺度的斑点。
如前所述,我们使用9×9
的模板对图像进行滤波,其结果作为最初始的尺度空间层(此时,尺度值为s=1.2,近似σ=1.2
的高斯微分),后续的层将通过逐步放大滤波模板尺寸,以及放大后的模板不断与图像进行滤波得到。由于采用盒子滤波和积分图像,滤波过程并不随着滤波模板尺寸的增加而使运算工作量增加。
与SIFT算法类似,我们需要将尺度空间划分为若干组(Octaves)。一个组代表了逐步放大的滤波模板对同一输入图像进行滤波的一系列响应图。每个组又由若干固定的层组成。由于积分图像离散化的原因,两个层之间的最小尺度变化量是由高斯二阶微分滤波器在微分方向上对正负斑点响应长度l0
决定的,它是盒子滤波器模板尺寸的1/3。对于9×9的模板,它的l0=3。一下层的响应长度至少应该在l0的基础上增加2个像素,以保证一边一个像素,即l0=5。这样模板的尺寸就为15×15。以此类推,我们可以得到一个尺寸增大模板序列,它们的尺寸分别为:9×9,15×15,21×21,27×27,黑色、白色区域的长度增加偶数个像素,以保证一个中心像素的存在。
采用类似的方法来处理其他几组的模板序列。其方法是将滤波器尺寸增加量翻倍(6,12,24,38)。这样,可以得到第二组的滤波器尺寸,它们分别为15,27,39,51。第三组的滤波器尺寸为27,51,75,99。如果原始图像的尺寸仍然大于对应的滤波器尺寸,尺度空间的分析还可以进行第四组,其对应的模板尺寸分别为51,99,147和195。下图显示了第一组至第三组的滤波器尺寸变化。
在通常尺度分析情况下,随着尺度的增大,被检测到的斑点数量迅速衰减。所以一般进行3-4组就可以了,与此同时,为了减少运算量,提高计算的速度,可以考虑在滤波时,将采样间隔设为2。
对于尺寸为L的模板,当用它与积分图运算来近似二维高斯核的滤波时,对应的二维高斯核的参数σ=1.2×L//9
,这一点至关重要,尤其是在后面计算描述子时,用于计算邻域的半径时。
9. 兴趣点的定位
为了在图像及不同尺寸中定位兴趣点,我们用了3×3×3邻域非最大值抑制。具体的步骤基本与SIFT一致,而且Hessian矩阵行列式的最大值在尺度和图像空间被插值。
下面显示了我们用的快速Hessian检测子检测到的兴趣点。
10. SURF源码解析
这份源码来自OpenCV nonfree模块。
10.1 主干函数 fastHessianDetector
特征点定位的主干函数为fastHessianDetector,该函数接受一个积分图像,以及尺寸相关的参数,组数与每组的层数,检测到的特征点保存在vector<KeyPoint>类型的结构中。
static void fastHessianDetector(const Mat& sum, const Mat& msum, vector<KeyPoint>& keypoints,
int nOctaves, int nOctaveLayers, float hessianThreshold)
{
/*first Octave图像采样的步长,第二组的时候加倍,以此内推
增加这个值,将会加快特征点检测的速度,但是会让特征点的提取变得不稳定*/
const int SAMPLE_STEP0 = 1;
int nTotalLayers = (nOctaveLayers + 2)*nOctaves; // 尺度空间的总图像数
int nMiddleLayers = nOctaveLayers*nOctaves; // 用于检测特征点的层的 总数,也就是中间层的总数
vector<Mat> dets(nTotalLayers); // 每一层图像 对应的 Hessian行列式的值
vector<Mat> traces(nTotalLayers); // 每一层图像 对应的 Hessian矩阵的迹的值
vector<int> sizes(nTotalLayers); // 每一层用的 Harr模板的大小
vector<int> sampleSteps(nTotalLayers); // 每一层用的采样步长
vector<int> middleIndices(nMiddleLayers); // 中间层的索引值
keypoints.clear();
// 为上面的对象分配空间,并赋予合适的值
int index = 0, middleIndex = 0, step = SAMPLE_STEP0;
for (int octave = 0; octave < nOctaves; octave++)
{
for (int layer = 0; layer < nOctaveLayers + 2; layer++)
{
/*这里sum.rows - 1是因为 sum是积分图,它的大小是原图像大小加1*/
dets[index].create((sum.rows - 1) / step, (sum.cols - 1) / step, CV_32F); // 这里面有除以遍历图像用的步长
traces[index].create((sum.rows - 1) / step, (sum.cols - 1) / step, CV_32F);
sizes[index] = (SURF_HAAR_SIZE0 + SURF_HAAR_SIZE_INC*layer) << octave;
sampleSteps[index] = step;
if (0 < layer && layer <= nOctaveLayers)
middleIndices[middleIndex++] = index;
index++;
}
step *= 2;
}
// Calculate hessian determinant and trace samples in each layer
for (int i = 0; i < nTotalLayers; i++)
{
calcLayerDetAndTrace(sum, sizes[i], sampleSteps[i], dets[i], traces[i]);
}
// Find maxima in the determinant of the hessian
for (int i = 0; i < nMiddleLayers; i++)
{
int layer = middleIndices[i];
int octave = i / nOctaveLayers;
findMaximaInLayer(sum, msum, dets, traces, sizes, keypoints, octave, layer, hessianThreshold, sampleSteps[layer]);
}
std::sort(keypoints.begin(), keypoints.end(), KeypointGreater());
}
10.2 计算Hessian矩阵的行列式与迹calcLayerDetAndTrace
这个函数首先定义了尺寸为9的第一层图像的三个模板。模板分别为一个3×5
、3×5、4×5的二维数组表示,数组的每一行表示一个黑白块的位置参数。函数里只初始化了第一层图像的模板参数,后面其他组其他层的Harr模板参数都是用resizeHaarPattern这个函数来计算的。这个函数返回的是一个SurfHF的结构体,这个结构体由两个点及一个权重构成。
struct SurfHF
{
int p0, p1, p2, p3;
float w;
SurfHF() : p0(0), p1(0), p2(0), p3(0), w(0) {}
};
resizeHaarPattern这个函数非常的巧妙,它把模板中的点坐标。转换到在积分图中的相对(模板左上角点)坐标。
static void
resizeHaarPattern(const int src[][5], SurfHF* dst, int n, int oldSize, int newSize, int widthStep)
{
float ratio = (float)newSize / oldSize;
for (int k = 0; k < n; k++)
{
int dx1 = cvRound(ratio*src[k][0]);
int dy1 = cvRound(ratio*src[k][1]);
int dx2 = cvRound(ratio*src[k][2]);
int dy2 = cvRound(ratio*src[k][3]);
/*巧妙的坐标转换*/
dst[k].p0 = dy1*widthStep + dx1; // 转换为一个相对距离,距离模板左上角点的 在积分图中的距离 !!important!!
dst[k].p1 = dy2*widthStep + dx1;
dst[k].p2 = dy1*widthStep + dx2;
dst[k].p3 = dy2*widthStep + dx2;
dst[k].w = src[k][4] / ((float)(dx2 - dx1)*(dy2 - dy1));// 原来的+1,+2用 覆盖的所有像素点平均。
}
}
在用积分图计算近似卷积时,用的是calcHaarPattern函数。这个函数比较简单,只用知道左上与右下角坐标即可。
inline float calcHaarPattern(const int* origin, const SurfHF* f, int n)
{
/*orgin即为积分图,n为模板中 黑白 块的个数 */
double d = 0;
for (int k = 0; k < n; k++)
d += (origin[f[k].p0] + origin[f[k].p3] - origin[f[k].p1] - origin[f[k].p2])*f[k].w;
return (float)d;
}
最终我们可以看到了整个calcLayerDetAndTrack的代码
static void calcLayerDetAndTrace(const Mat& sum, int size, int sampleStep,
Mat& det, Mat& trace)
{
const int NX = 3, NY = 3, NXY = 4;
const int dx_s[NX][5] = { { 0, 2, 3, 7, 1 }, { 3, 2, 6, 7, -2 }, { 6, 2, 9, 7, 1 } };
const int dy_s[NY][5] = { { 2, 0, 7, 3, 1 }, { 2, 3, 7, 6, -2 }, { 2, 6, 7, 9, 1 } };
const int dxy_s[NXY][5] = { { 1, 1, 4, 4, 1 }, { 5, 1, 8, 4, -1 }, { 1, 5, 4, 8, -1 }, { 5, 5, 8, 8, 1 } };
SurfHF Dx[NX], Dy[NY], Dxy[NXY];
if (size > sum.rows - 1 || size > sum.cols - 1)
return;
resizeHaarPattern(dx_s, Dx, NX, 9, size, sum.cols);
resizeHaarPattern(dy_s, Dy, NY, 9, size, sum.cols);
resizeHaarPattern(dxy_s, Dxy, NXY, 9, size, sum.cols);
/* The integral image 'sum' is one pixel bigger than the source image */
int samples_i = 1 + (sum.rows - 1 - size) / sampleStep; // 最大能遍历到的 行坐标,因为要减掉一个模板的尺寸
int samples_j = 1 + (sum.cols - 1 - size) / sampleStep; // 最大能遍历到的 列坐标
/* Ignore pixels where some of the kernel is outside the image */
int margin = (size / 2) / sampleStep;
for (int i = 0; i < samples_i; i++)
{
/*坐标为(i,j)的点是模板左上角的点,所以实际现在模板分析是的i+margin,j+margin点处的响应*/
const int* sum_ptr = sum.ptr<int>(i*sampleStep);
float* det_ptr = &det.at<float>(i + margin, margin); // 左边空隙为 margin
float* trace_ptr = &trace.at<float>(i + margin, margin);
for (int j = 0; j < samples_j; j++)
{
float dx = calcHaarPattern(sum_ptr, Dx, 3);
float dy = calcHaarPattern(sum_ptr, Dy, 3);
float dxy = calcHaarPattern(sum_ptr, Dxy, 4);
sum_ptr += sampleStep;
det_ptr[j] = dx*dy - 0.81f*dxy*dxy;
trace_ptr[j] = dx + dy;
}
}
}
10.3 局部最大值搜索findMaximaInLayer
这里算法思路很简单,值得注意的是里面的一些坐标的转换很巧妙,里面比较重的函数就是interpolateKeypoint函数,通过插值计算最大值点。
/*
* Maxima location interpolation as described in "Invariant Features from
* Interest Point Groups" by Matthew Brown and David Lowe. This is performed by
* fitting a 3D quadratic to a set of neighbouring samples.
*
* The gradient vector and Hessian matrix at the initial keypoint location are
* approximated using central differences. The linear system Ax = b is then
* solved, where A is the Hessian, b is the negative gradient, and x is the
* offset of the interpolated maxima coordinates from the initial estimate.
* This is equivalent to an iteration of Netwon's optimisation algorithm.
*
* N9 contains the samples in the 3x3x3 neighbourhood of the maxima
* dx is the sampling step in x
* dy is the sampling step in y
* ds is the sampling step in size
* point contains the keypoint coordinates and scale to be modified
*
* Return value is 1 if interpolation was successful, 0 on failure.
*/
static int
interpolateKeypoint(float N9[3][9], int dx, int dy, int ds, KeyPoint& kpt)
{
Vec3f b(-(N9[1][5] - N9[1][3]) / 2, // Negative 1st deriv with respect to x
-(N9[1][7] - N9[1][1]) / 2, // Negative 1st deriv with respect to y
-(N9[2][4] - N9[0][4]) / 2); // Negative 1st deriv with respect to s
Matx33f A(
N9[1][3] - 2 * N9[1][4] + N9[1][5], // 2nd deriv x, x
(N9[1][8] - N9[1][6] - N9[1][2] + N9[1][0]) / 4, // 2nd deriv x, y
(N9[2][5] - N9[2][3] - N9[0][5] + N9[0][3]) / 4, // 2nd deriv x, s
(N9[1][8] - N9[1][6] - N9[1][2] + N9[1][0]) / 4, // 2nd deriv x, y
N9[1][1] - 2 * N9[1][4] + N9[1][7], // 2nd deriv y, y
(N9[2][7] - N9[2][1] - N9[0][7] + N9[0][1]) / 4, // 2nd deriv y, s
(N9[2][5] - N9[2][3] - N9[0][5] + N9[0][3]) / 4, // 2nd deriv x, s
(N9[2][7] - N9[2][1] - N9[0][7] + N9[0][1]) / 4, // 2nd deriv y, s
N9[0][4] - 2 * N9[1][4] + N9[2][4]); // 2nd deriv s, s
Vec3f x = A.solve(b, DECOMP_LU);
bool ok = (x[0] != 0 || x[1] != 0 || x[2] != 0) &&
std::abs(x[0]) <= 1 && std::abs(x[1]) <= 1 && std::abs(x[2]) <= 1;
if (ok)
{
kpt.pt.x += x[0] * dx;
kpt.pt.y += x[1] * dy;
kpt.size = (float)cvRound(kpt.size + x[2] * ds);
}
return ok;
}
static void findMaximaInLayer(const Mat& sum, const Mat& mask_sum,
const vector<Mat>& dets, const vector<Mat>& traces,
const vector<int>& sizes, vector<KeyPoint>& keypoints,
int octave, int layer, float hessianThreshold, int sampleStep)
{
// Wavelet Data
const int NM = 1;
const int dm[NM][5] = { { 0, 0, 9, 9, 1 } };
SurfHF Dm;
int size = sizes[layer];
// 当前层图像的大小
int layer_rows = (sum.rows - 1) / sampleStep;
int layer_cols = (sum.cols - 1) / sampleStep;
// 边界区域大小,考虑的下一层的模板大小
int margin = (sizes[layer + 1] / 2) / sampleStep + 1;
if (!mask_sum.empty())
resizeHaarPattern(dm, &Dm, NM, 9, size, mask_sum.cols);
int step = (int)(dets[layer].step / dets[layer].elemSize());
for (int i = margin; i < layer_rows - margin; i++)
{
const float* det_ptr = dets[layer].ptr<float>(i);
const float* trace_ptr = traces[layer].ptr<float>(i);
for (int j = margin; j < layer_cols - margin; j++)
{
float val0 = det_ptr[j]; // 中心点的值
if (val0 > hessianThreshold)
{
// 模板左上角的坐标
int sum_i = sampleStep*(i - (size / 2) / sampleStep);
int sum_j = sampleStep*(j - (size / 2) / sampleStep);
/* The 3x3x3 neighbouring samples around the maxima.
The maxima is included at N9[1][4] */
const float *det1 = &dets[layer - 1].at<float>(i, j);
const float *det2 = &dets[layer].at<float>(i, j);
const float *det3 = &dets[layer + 1].at<float>(i, j);
float N9[3][9] = { { det1[-step - 1], det1[-step], det1[-step + 1],
det1[-1], det1[0], det1[1],
det1[step - 1], det1[step], det1[step + 1] },
{ det2[-step - 1], det2[-step], det2[-step + 1],
det2[-1], det2[0], det2[1],
det2[step - 1], det2[step], det2[step + 1] },
{ det3[-step - 1], det3[-step], det3[-step + 1],
det3[-1], det3[0], det3[1],
det3[step - 1], det3[step], det3[step + 1] } };
/* Check the mask - why not just check the mask at the center of the wavelet? */
if (!mask_sum.empty())
{
const int* mask_ptr = &mask_sum.at<int>(sum_i, sum_j);
float mval = calcHaarPattern(mask_ptr, &Dm, 1);
if (mval < 0.5)
continue;
}
/* 检测val0,是否在N9里极大值,??为什么不检测极小值呢*/
if (val0 > N9[0][0] && val0 > N9[0][1] && val0 > N9[0][2] &&
val0 > N9[0][3] && val0 > N9[0][4] && val0 > N9[0][5] &&
val0 > N9[0][6] && val0 > N9[0][7] && val0 > N9[0][8] &&
val0 > N9[1][0] && val0 > N9[1][1] && val0 > N9[1][2] &&
val0 > N9[1][3] && val0 > N9[1][5] &&
val0 > N9[1][6] && val0 > N9[1][7] && val0 > N9[1][8] &&
val0 > N9[2][0] && val0 > N9[2][1] && val0 > N9[2][2] &&
val0 > N9[2][3] && val0 > N9[2][4] && val0 > N9[2][5] &&
val0 > N9[2][6] && val0 > N9[2][7] && val0 > N9[2][8])
{
/* Calculate the wavelet center coordinates for the maxima */
float center_i = sum_i + (size - 1)*0.5f;
float center_j = sum_j + (size - 1)*0.5f;
KeyPoint kpt(center_j, center_i, (float)sizes[layer],
-1, val0, octave, CV_SIGN(trace_ptr[j]));
/* 局部极大值插值,用Hessian,类似于SIFT里的插值,里面没有迭代5次,只进行了一次查找,why? */
int ds = size - sizes[layer - 1];
int interp_ok = interpolateKeypoint(N9, sampleStep, sampleStep, ds, kpt);
/* Sometimes the interpolation step gives a negative size etc. */
if (interp_ok)
{
/*printf( "KeyPoint %f %f %d\n", point.pt.x, point.pt.y, point.size );*/
keypoints.push_back(kpt);
}
}
}
}
}
}
11. SURF特征点方向分配
为了保证特征矢量具有旋转不变性,与SIFT特征一样,需要对每个特征点分配一个主方向。为些,我们需要以特征点为中心,以6s(s=1.2∗L/9为特征点的尺度)为半径的圆形区域,对图像进行Haar小波响应运算。这样做实际就是对图像进行梯度运算只不过是我们需要利用积分图像,提高计算图像梯度的效率。在SIFT特征描述子中我们在求取特征点主方向时,以是特征点为中心,在以4.5σ
为半径的邻域内计算梯度方向直方图。事实上,两种方法在求取特征点主方向时,考虑到Haar小波的模板带宽,实际计算梯度的图像区域是相同的。用于计算梯度的Harr小波的尺度为4s。
与SIFT类似,使用σ=2s
的高斯加权函数对Harr小波的响应值进行高斯加权。为了求取主方向值,需要设计一个以特征点为中心,张角为π/3的扇形滑动窗口。以步长为0.2弧度左右,旋转这个滑动窗口,并对滑动窗口内的图像Harr小波响应值dx、dy进行累加,得到一个矢量(mw,θw):
主方向为最大Harr响应累加值所对应的方向,也就是最长矢量所对应的方向,即
可以依照SIFT求方方向时策略,当存在另一个相当于主峰值80%能量的峰值时,则将这个方向认为是该特征点的辅方向。一个特征点可能会被指定具有多个方向(一个主方向,一个以上辅方向),这可以增强匹配的鲁棒性。和SIFT的描述子类似,如果在mw中出现另一个大于主峰能量max{mw}80时的次峰,可以将该特征点复制成两个特征点。一个主的方向为最大响应能量所对应的方向,另一个主方向为次大响应能量所对应的方向。
图 1 求取主方向时扇形滑动窗口围绕特征点转动,统计Haar小波响应值,并计算方向角
12. 特征点特征矢量生成
生成特征点描述子与确定特征点方向有些类似,它需要计算图像的Haar小波响应。不过,与主方向的确定不同的是,这次我们不是使用一个圆形区域,而是在一个矩形区域来计算Haar小波响应。以特征点为中心,沿上一节讨论得到的主方向,沿主方向将s20s×20s
的图像划分为4×4个子块,每个子块利用尺寸2s的Harr模板进行响应值进行响应值计算,然后对响应值进行统计∑dx、∑|dx|、∑dy、∑|dy|形成特征矢量。如下图2所示。图中,以特征点为中心,以20s为边长的矩形窗口为特征描述子计算使用的窗口,特征点到矩形边框的线段表示特征点的主方向。
图2 特征描述子表示
将20s的窗口划分成4×4子窗口,每个子窗口有5s×5s个像素。使用尺寸为2s的Harr小波对子窗口图像进行其响应值计算,共进行25次采样,分别得到沿主方向的dy和垂直于主方向的dx。然后,以特征点为中心,对dy和dx进行高斯加权计算,高斯核的参数为σ=3.3s(即20s/6)。最后,分别对每个子块的响应值进行统计,得到每个子块的矢量:
由于共有4×4个子块,因此,特征描述子共由4×4×4=64维特征矢量组成。SURF描述子不仅具有尺度和旋转不变性,而且对光照的变化也具有不变性。使小波响应本身就具有亮度不变性,而对比度的不变性则是通过将特征矢量进行归一化来实现。图3 给出了三种不同图像模式的子块得到的不同结果。对于实际图像的描述子,我们可以认为它们是由这三种不同模式图像的描述子组合而成的。
图3 不同的图像密度模式得到的不同的描述子结果
为了充分利用积分图像进行Haar小波的响应计算,我们并不直接旋转Haar小波模板求得其响应值,而是在积图像上先使用水平和垂直的Haar模板求得响应值dy和dx,然后根据主方向旋转dx和dy与主方向操持一致,如下图4所示。为了求得旋转后Haar小波响应值,首先要得到旋转前图像的位置。旋转前后图偈的位置关系,可以通过点的旋转公式得到:
在得到点(j,i)在旋转前对应积分图像的位置(x,y)后,利用积分图像与水平、垂直Harr小波,求得水平与垂直两个方向的响应值dx和dy。对dx和dy进行高斯加权处理,并根据主方向的角度,对dx和dy进行旋转变换,从而,得到旋转后的dx’和dy’。其计算公式如下:
图4 利用积分图像进行Haar小波响应计算示意图,左边为旋转后的图像,右边为旋转前的图像
13. 特征描述子的维数
一般而言,特征矢量的长度越长,特征矢量所承载的信息量就越大,特征描述子的独特性就越好,但匹配时所付出的时间代价就越大。对于SURF描述子,可以将它扩展到用128维矢量来表示。具体方法是在求∑dx
、∑|dx|时区分dy<0和dy≥0情况。同时,在求取∑dy、∑|dy|时区分dx<0和dx≥0情况。这样,每个子块就产生了8个梯度统计值,从而使描述子特征矢量的长度增加到8×4×4=128维。
为了实现快速匹配,SURF在特征矢量中增加了一个新的变量,即特征点的拉普拉斯响应正负号。在特征点检测时,将Hessian矩阵的迹的正负号记录下来,作为特征矢量中的一个变量。这样做并不增加运算量,因为特征点检测进已经对Hessian矩阵的迹进行了计算。在特征匹配时,这个变量可以有效地节省搜索的时间,因为只有两个具有相同正负号的特征点才有可能匹配,对于正负号不同的特征点就不进行相似性计算。
简单地说,我们可以根据特征点的响应值符号,将特征点分成两组,一组是具有拉普拉斯正响应的特征点,一组是具有拉普拉斯负响应的特征点,匹配时,只有符号相同组中的特征点才能进行相互匹配。显然,这样可以节省特征点匹配的时间。如下图5所示。
图5 黑背景下的亮斑和白背景下的黑斑 因为它们的拉普拉斯响应正负号不同,不会对它们进行匹配
14. 源码解析
特征点描述子的生成这一部分的代码主要是通过SURFInvoker这个类来实现。在主流程中,通过一个parallel_for_()函数来并发计算。
struct SURFInvoker
{
enum{ORI_RADIUS = 6, ORI_WIN = 60, PATCH_SZ = 20};
// Parameters
const Mat* img;
const Mat* sum;
vector<KeyPoint>* keypoints;
Mat* descriptors;
bool extended;
bool upright;
// Pre-calculated values
int nOriSamples;
vector<Point> apt; // 特征点周围用于描述方向的邻域的点
vector<float> aptw; // 描述 方向时的 高斯 权
vector<float> DW;
SURFInvoker(const Mat& _img, const Mat& _sum,
vector<KeyPoint>& _keypoints, Mat& _descriptors,
bool _extended, bool _upright)
{
keypoints = &_keypoints;
descriptors = &_descriptors;
img = &_img;
sum = &_sum;
extended = _extended;
upright = _upright;
// 用于描述特征点的 方向的 邻域大小: 12*sigma+1 (sigma =1.2) 因为高斯加权的核的参数为2sigma
// nOriSampleBound为 矩形框内点的个数
const int nOriSampleBound = (2 * ORI_RADIUS + 1)*(2 * ORI_RADIUS + 1); // 这里把s近似为1 ORI_DADIUS = 6s
// 分配大小
apt.resize(nOriSampleBound);
aptw.resize(nOriSampleBound);
DW.resize(PATCH_SZ*PATCH_SZ); // PATHC_SZ为特征描述子的 区域大小 20s(s 这里初始为1了)
/* 计算特征点方向用的 高斯分布 权值与坐标 */
Mat G_ori = getGaussianKernel(2 * ORI_RADIUS + 1, SURF_ORI_SIGMA, CV_32F); // SURF_ORI_SIGMA = 1.2 *2 =2.5
nOriSamples = 0;
for (int i = -ORI_RADIUS; i <= ORI_RADIUS; i++)
{
for (int j = -ORI_RADIUS; j <= ORI_RADIUS; j++)
{
if (i*i + j*j <= ORI_RADIUS*ORI_RADIUS) // 限制在圆形区域内
{
apt[nOriSamples] = cvPoint(i, j);
// 下面这里有个坐标转换,因为i,j都是从-ORI_RADIUS开始的。
aptw[nOriSamples++] = G_ori.at<float>(i + ORI_RADIUS, 0) * G_ori.at<float>(j + ORI_RADIUS, 0);
}
}
}
CV_Assert(nOriSamples <= nOriSampleBound); // nOriSamples为圆形区域内的点,nOriSampleBound是正方形区域的点
/* 用于特征点描述子的高斯 权值 */
Mat G_desc = getGaussianKernel(PATCH_SZ, SURF_DESC_SIGMA, CV_32F); // 用于生成特征描述子的 高斯加权 sigma = 3.3s (s初取1)
for (int i = 0; i < PATCH_SZ; i++)
{
for (int j = 0; j < PATCH_SZ; j++)
DW[i*PATCH_SZ + j] = G_desc.at<float>(i, 0) * G_desc.at<float>(j, 0);
}
/* x与y方向上的 Harr小波,参数为4s */
const int NX = 2, NY = 2;
const int dx_s[NX][5] = { { 0, 0, 2, 4, -1 }, { 2, 0, 4, 4, 1 } };
const int dy_s[NY][5] = { { 0, 0, 4, 2, 1 }, { 0, 2, 4, 4, -1 } };
float X[nOriSampleBound], Y[nOriSampleBound], angle[nOriSampleBound]; // 用于计算特生点主方向
uchar PATCH[PATCH_SZ + 1][PATCH_SZ + 1];
float DX[PATCH_SZ][PATCH_SZ], DY[PATCH_SZ][PATCH_SZ]; // 20s * 20s区域的 梯度值
CvMat matX = cvMat(1, nOriSampleBound, CV_32F, X);
CvMat matY = cvMat(1, nOriSampleBound, CV_32F, Y);
CvMat _angle = cvMat(1, nOriSampleBound, CV_32F, angle);
Mat _patch(PATCH_SZ + 1, PATCH_SZ + 1, CV_8U, PATCH);
int dsize = extended ? 128 : 64;
int k, k1 = 0, k2 = (int)(*keypoints).size();// k2为Harr小波的 模板尺寸
float maxSize = 0;
for (k = k1; k < k2; k++)
{
maxSize = std::max(maxSize, (*keypoints)[k].size);
}
// maxSize*1.2/9 表示最大的尺度 s
int imaxSize = std::max(cvCeil((PATCH_SZ + 1)*maxSize*1.2f / 9.0f), 1);
Ptr<CvMat> winbuf = cvCreateMat(1, imaxSize*imaxSize, CV_8U);
for (k = k1; k < k2; k++)
{
int i, j, kk, nangle;
float* vec;
SurfHF dx_t[NX], dy_t[NY];
KeyPoint& kp = (*keypoints)[k];
float size = kp.size;
Point2f center = kp.pt;
/* s是当前层的尺度参数 1.2是第一层的参数,9是第一层的模板大小*/
float s = size*1.2f / 9.0f;
/* grad_wav_size是 harr梯度模板的大小 边长为 4s */
int grad_wav_size = 2 * cvRound(2 * s);
if (sum->rows < grad_wav_size || sum->cols < grad_wav_size)
{
/* when grad_wav_size is too big,
* the sampling of gradient will be meaningless
* mark keypoint for deletion. */
kp.size = -1;
continue;
}
float descriptor_dir = 360.f - 90.f;
if (upright == 0)
{
// 这一步 是计算梯度值,先将harr模板放大,再根据积分图计算,与前面求D_x,D_y一致类似
resizeHaarPattern(dx_s, dx_t, NX, 4, grad_wav_size, sum->cols);
resizeHaarPattern(dy_s, dy_t, NY, 4, grad_wav_size, sum->cols);
for (kk = 0, nangle = 0; kk < nOriSamples; kk++)
{
int x = cvRound(center.x + apt[kk].x*s - (float)(grad_wav_size - 1) / 2);
int y = cvRound(center.y + apt[kk].y*s - (float)(grad_wav_size - 1) / 2);
if (y < 0 || y >= sum->rows - grad_wav_size ||
x < 0 || x >= sum->cols - grad_wav_size)
continue;
const int* ptr = &sum->at<int>(y, x);
float vx = calcHaarPattern(ptr, dx_t, 2);
float vy = calcHaarPattern(ptr, dy_t, 2);
X[nangle] = vx*aptw[kk];
Y[nangle] = vy*aptw[kk];
nangle++;
}
if (nangle == 0)
{
// No gradient could be sampled because the keypoint is too
// near too one or more of the sides of the image. As we
// therefore cannot find a dominant direction, we skip this
// keypoint and mark it for later deletion from the sequence.
kp.size = -1;
continue;
}
matX.cols = matY.cols = _angle.cols = nangle;
// 计算邻域内每个点的 梯度角度
cvCartToPolar(&matX, &matY, 0, &_angle, 1);
float bestx = 0, besty = 0, descriptor_mod = 0;
for (i = 0; i < 360; i += SURF_ORI_SEARCH_INC) // SURF_ORI_SEARCH_INC 为扇形区域扫描的步长
{
float sumx = 0, sumy = 0, temp_mod;
for (j = 0; j < nangle; j++)
{
// d是 分析到的那个点与 现在主方向的偏度
int d = std::abs(cvRound(angle[j]) - i);
if (d < ORI_WIN / 2 || d > 360 - ORI_WIN / 2)
{
sumx += X[j];
sumy += Y[j];
}
}
temp_mod = sumx*sumx + sumy*sumy;
// descriptor_mod 是最大峰值
if (temp_mod > descriptor_mod)
{
descriptor_mod = temp_mod;
bestx = sumx;
besty = sumy;
}
}
descriptor_dir = fastAtan2(-besty, bestx);
}
kp.angle = descriptor_dir;
if (!descriptors || !descriptors->data)
continue;
/* 用特征点周围20*s为边长的邻域 计算特征描述子 */
int win_size = (int)((PATCH_SZ + 1)*s);
CV_Assert(winbuf->cols >= win_size*win_size);
Mat win(win_size, win_size, CV_8U, winbuf->data.ptr);
if (!upright)
{
descriptor_dir *= (float)(CV_PI / 180); // 特征点的主方向 弧度值
float sin_dir = -std::sin(descriptor_dir); // - sin dir
float cos_dir = std::cos(descriptor_dir);
float win_offset = -(float)(win_size - 1) / 2;
float start_x = center.x + win_offset*cos_dir + win_offset*sin_dir;
float start_y = center.y - win_offset*sin_dir + win_offset*cos_dir;
uchar* WIN = win.data;
int ncols1 = img->cols - 1, nrows1 = img->rows - 1;
size_t imgstep = img->step;
for (i = 0; i < win_size; i++, start_x += sin_dir, start_y += cos_dir)
{
double pixel_x = start_x;
double pixel_y = start_y;
for (j = 0; j < win_size; j++, pixel_x += cos_dir, pixel_y -= sin_dir)
{
int ix = cvFloor(pixel_x), iy = cvFloor(pixel_y);
if ((unsigned)ix < (unsigned)ncols1 &&
(unsigned)iy < (unsigned)nrows1)
{
float a = (float)(pixel_x - ix), b = (float)(pixel_y - iy);
const uchar* imgptr = &img->at<uchar>(iy, ix);
WIN[i*win_size + j] = (uchar)
cvRound(imgptr[0] * (1.f - a)*(1.f - b) +
imgptr[1] * a*(1.f - b) +
imgptr[imgstep] * (1.f - a)*b +
imgptr[imgstep + 1] * a*b);
}
else
{
int x = std::min(std::max(cvRound(pixel_x), 0), ncols1);
int y = std::min(std::max(cvRound(pixel_y), 0), nrows1);
WIN[i*win_size + j] = img->at<uchar>(y, x);
}
}
}
}
else
{
float win_offset = -(float)(win_size - 1) / 2;
int start_x = cvRound(center.x + win_offset);
int start_y = cvRound(center.y - win_offset);
uchar* WIN = win.data;
for (i = 0; i < win_size; i++, start_x++)
{
int pixel_x = start_x;
int pixel_y = start_y;
for (j = 0; j < win_size; j++, pixel_y--)
{
int x = MAX(pixel_x, 0);
int y = MAX(pixel_y, 0);
x = MIN(x, img->cols - 1);
y = MIN(y, img->rows - 1);
WIN[i*win_size + j] = img->at<uchar>(y, x);
}
}
}
// Scale the window to size PATCH_SZ so each pixel's size is s. This
// makes calculating the gradients with wavelets of size 2s easy
resize(win, _patch, _patch.size(), 0, 0, INTER_AREA);
// Calculate gradients in x and y with wavelets of size 2s
for (i = 0; i < PATCH_SZ; i++)
for (j = 0; j < PATCH_SZ; j++)
{
float dw = DW[i*PATCH_SZ + j]; // 高斯加权系数
float vx = (PATCH[i][j + 1] - PATCH[i][j] + PATCH[i + 1][j + 1] - PATCH[i + 1][j])*dw;
float vy = (PATCH[i + 1][j] - PATCH[i][j] + PATCH[i + 1][j + 1] - PATCH[i][j + 1])*dw;
DX[i][j] = vx;
DY[i][j] = vy;
}
// Construct the descriptor
vec = descriptors->ptr<float>(k);
for (kk = 0; kk < dsize; kk++)
vec[kk] = 0;
double square_mag = 0;
if (extended)
{
// 128维描述子,考虑dx与dy的正负号
for (i = 0; i < 4; i++)
for (j = 0; j < 4; j++)
{
// 每个方块内是一个5s * 5s的区域,每个方法由8个特征描述
for (int y = i * 5; y < i * 5 + 5; y++)
{
for (int x = j * 5; x < j * 5 + 5; x++)
{
float tx = DX[y][x], ty = DY[y][x];
if (ty >= 0)
{
vec[0] += tx;
vec[1] += (float)fabs(tx);
}
else {
vec[2] += tx;
vec[3] += (float)fabs(tx);
}
if (tx >= 0)
{
vec[4] += ty;
vec[5] += (float)fabs(ty);
}
else {
vec[6] += ty;
vec[7] += (float)fabs(ty);
}
}
}
for (kk = 0; kk < 8; kk++)
square_mag += vec[kk] * vec[kk];
vec += 8;
}
}
else
{
// 64位描述子
for (i = 0; i < 4; i++)
for (j = 0; j < 4; j++)
{
for (int y = i * 5; y < i * 5 + 5; y++)
{
for (int x = j * 5; x < j * 5 + 5; x++)
{
float tx = DX[y][x], ty = DY[y][x];
vec[0] += tx; vec[1] += ty;
vec[2] += (float)fabs(tx); vec[3] += (float)fabs(ty);
}
}
for (kk = 0; kk < 4; kk++)
square_mag += vec[kk] * vec[kk];
vec += 4;
}
}
// 归一化 描述子 以满足 光照不变性
vec = descriptors->ptr<float>(k);
float scale = (float)(1. / (sqrt(square_mag) + DBL_EPSILON));
for (kk = 0; kk < dsize; kk++)
vec[kk] *= scale;
}
}
};
15. 总结
总体来说,如果理解了SIFT算法,再来看SURF算法会发现思路非常简单。尤其是局部最大值查找方面,基本一致。关键还是一个用积分图来简化卷积的思路,以及怎么用不同的模板来近似原来尺度空间中的高斯滤波器。
实际上有文献指出,SURF比SIFT工作更出色。他们认为主要是因为SURF在求取描述子特征矢量时,是对一个子块的梯度信息进行求和,而SIFT则是依靠单个像素梯度的方向。
16. surf特征+FLANN特征匹配+knn筛选匹配点+单应性矩阵映射
https://blog.csdn.net/panda1234lee/article/details/10896099
FlannBasedMatcher.match效果不如FlannBasedMatcher.knnMatch后再进行比较筛选的的方法。
为了排除因为图像遮挡和背景混乱而产生的无匹配关系的关键点,SIFT的作者Lowe提出了比较最近邻距离与次近邻距离的SIFT匹配方式:取一幅图像中的一个SIFT关键点,并找出其与另一幅图像中欧式距离最近的前两个关键点,在这两个关键点中,如果最近的距离除以次近的距离得到的比率ratio少于某个阈值T,则接受这一对匹配点。因为对于错误匹配,由于特征空间的高维性,相似的距离可能有大量其他的错误匹配,从而它的ratio值比较高。显然降低这个比例阈值T,SURF匹配点数目会减少,但更加稳定,反之亦然。
Lowe推荐ratio的阈值为0.8,但作者对大量任意存在尺度、旋转和亮度变化的两幅图片进行匹配,结果表明ratio取值在0. 4~0. 6 之间最佳,小于0. 4的很少有匹配点,大于0. 6的则存在大量错误匹配点,所以建议ratio的取值原则如下:
ratio=0. 4:对于准确度要求高的匹配;
ratio=0. 6:对于匹配点数目要求比较多的匹配;
ratio=0. 5:一般情况下。
17. findFundamentalMat
使用该函数效果不好?Nister的五点法实现效果更好。
Even if your algorithm is correct, 8 point F matrix computation is very error prone due to image noise. The lesser correspondences you use the better. The best you can do is doing 5 point Essential (E) matrix computation, but that would require you to pre-calibrate the camera and convert the detected pixel image points after SIFT/SURF to normalized pixels (metric pixel locations). Then apply Nister's 5-point algorithm either from the freely available Matlab implementation or from Bundler (c++ implementation by Noah Snavely). In my experience with SfM, 5-point E matrix is much much better/stable than 7 or 8 point F matrix computation. And ofcourse do RANSAC after 5 point to get more robust estimates. Hope this helps.
即使你的算法是正确的,8点F矩阵计算由于图像噪声。使用较少的对应性越好。最好的做法是做5点基本(E)矩阵计算,但这将需要您预先校准摄像机,并将检测到的像素图像点SIFT / SURF后转换为归一化像素(度量像素位置)。然后应用Nister的5点算法,从免费的Matlab实现或Bundler(由Noah Snavely的c ++实现)。在我的SfM经验中,5点E矩阵比7或8点F矩阵计算好得多/稳定。而且在5点后做RANSAC得到更稳健的估计。
18. findHomography
使用findHomography接口获取两张图之间的单应性矩阵H,获得的结果不具有可逆性。举例来说,两张图M1和M2,那么调用findHomography(M1, M2, CV_RANSAC, 4)得到的矩阵H12和findHomography(M2, M1, CV_RANSAC, 4)得到的矩阵H21不具有可逆关系,即:H12的逆不等于H21,H12和H21的逆存在一定的误差,这个误差的产生是因为这个单应性矩阵的求取本身就是采用了RANSAC算法,得到的估算矩阵,是个估算值,所以存在误差。