函数处理流程
提取圆形网格有两种方法，第一种是基于聚类的方法，第二种是基于距离向量的方法，此外，圆形网格还有对称和非对称，本篇文章将介绍第二种方法提取对称圆形网格。提取流程如下：
源码分析
本文实验时是以下面这张图片作为输入，Debug分析函数的处理流程（关于opencv源码调试的配置网上由许多教程）。后面展示的代码会将源码中部分无关的代码给删除，便于算法的理解。
findCirclesGrid源码
该函数内部直接调用findCirclesGrid2函数
bool findCirclesGrid(InputArray _image, Size patternSize,
OutputArray _centers, int flags, const Ptr<FeatureDetector> &blobDetector)
{
return cv::findCirclesGrid2(_image, patternSize, _centers, flags, blobDetector, CirclesGridFinderParameters2());
}
findCirclesGrid2 函数源码
提取网格前，会先调用FeatureDetector类的detect方法检测斑点，对于该函数有许多文章介绍，本文不分析，从输入图片中共检测到如下57个斑点(不同的参数检测到的斑点可能会不同)。
{x=450.351746 y=541.933289 }
{x=702.884644 y=541.117859 }
{x=660.913635 y=541.265869 }
{x=618.978760 y=541.394836 }
{x=576.906372 y=541.419922 }
{x=534.894287 y=541.474426 }
{x=492.848663 y=541.467468 }
{x=450.299591 y=500.272919 }
{x=703.152405 y=499.412323 }
{x=661.153748 y=499.559418 }
{x=619.142029 y=499.691437 }
{x=577.047241 y=499.749542 }
{x=534.983032 y=499.820282 }
{x=492.865265 y=499.785675 }
{x=619.321289 y=457.898041 }
{x=577.148376 y=457.953430 }
{x=535.045349 y=458.008728 }
{x=492.900116 y=457.997284 }
{x=450.287537 y=458.460907 }
{x=703.406006 y=457.627686 }
{x=661.354248 y=457.777527 }
{x=619.425903 y=415.991364 }
{x=577.232239 y=416.060608 }
{x=535.069641 y=416.080322 }
{x=492.906128 y=416.133026 }
{x=450.271027 y=416.573059 }
{x=703.605286 y=415.771484 }
{x=661.530212 y=415.896912 }
{x=450.200745 y=374.489288 }
{x=703.747742 y=373.774719 }
{x=661.654114 y=373.868896 }
{x=619.513733 y=373.977051 }
{x=577.312500 y=374.025116 }
{x=535.106689 y=374.059265 }
{x=492.884338 y=374.070251 }
{x=450.173920 y=332.365631 }
{x=703.903564 y=331.690369 }
{x=661.787842 y=331.757385 }
{x=619.618530 y=331.834778 }
{x=577.357910 y=331.877655 }
{x=535.137329 y=331.917389 }
{x=492.884796 y=331.933594 }
{x=704.000427 y=289.565308 }
{x=661.844299 y=289.628479 }
{x=619.650513 y=289.658173 }
{x=577.386047 y=289.713257 }
{x=535.139709 y=289.759460 }
{x=492.869232 y=289.746613 }
{x=450.145630 y=290.128418 }
{x=445.127472 y=102.963287 }
{x=475.568237 y=90.3130951 }
{x=508.711945 y=92.3338394 }
{x=509.194763 y=80.1243744 }
{x=504.870880 y=607.987488 }
{x=367.777588 y=106.329857 }
{x=802.892395 y=30.8430290 }
{x=485.435120 y=609.837585 }
//InputArray _image:输入图片; Size patternSize:网格大小; _centers:提取的网格点中心;
//flags:是否为对称网格; blobDetector:斑点检测器; parameters:参数
bool findCirclesGrid2(InputArray _image, Size patternSize,
OutputArray _centers, int flags, const Ptr<FeatureDetector> &blobDetector,
CirclesGridFinderParameters2 parameters)
{
//判断是否为对称网格
bool isAsymmetricGrid = (flags & CALIB_CB_ASYMMETRIC_GRID) ? true : false;
bool isSymmetricGrid = (flags & CALIB_CB_SYMMETRIC_GRID ) ? true : false;
//圆形网格点中心
std::vector<Point2f> centers;
//斑点检测得到的斑点中心
std::vector<Point2f> points;
if (blobDetector)
{
std::vector<KeyPoint> keypoints;
//斑点检测
blobDetector->detect(_image, keypoints);
for (size_t i = 0; i < keypoints.size(); i++)
{
points.push_back(keypoints[i].pt);
}
}
bool isValid = false;
const int attempts = 2;
const size_t minHomographyPoints = 4;
Mat H;
for (int i = 0; i < attempts; i++)
{
centers.clear();
//对称圆形网格点检测器
CirclesGridFinder boxFinder(patternSize, points, parameters);
//提取圆形网格点
bool isFound = boxFinder.findHoles();
//提取成功
if (isFound)
{
//获取圆形网格点
boxFinder.getHoles(centers);
isValid = true;
break; // done, return result
}
//提取失败时,获取提到的部分圆形网格点
boxFinder.getHoles(centers);
//若不是在最后一次提取
if (i != attempts - 1)
{
//minHomographyPoints：表示计算单应性矩阵需要最少特征点数量
if (centers.size() < minHomographyPoints)
break;
/*计算已提到的部分圆形网格点到一个标准网格的单应性矩阵，同时将斑点检测到的所有斑点通过这个单应性矩阵投影到
标x准网格。这样做的原因是该方法基于距离向量来提取，若是拍摄图片倾斜严重，会使得网格中相邻点间的距
离变化较大，导致提取失败，投影到标准网格后，相邻点间的距离变化就比较小，更容易提取成功。
*/
H = CirclesGridFinder::rectifyGrid(boxFinder.getDetectedGridSize(), centers, points, points);
}
}
//若是做了单应性投影提取，则要将提取到的圆形网格点投影回来
if (!centers.empty() && !H.empty()) // undone rectification
{
Mat orgPointsMat;
transform(centers, orgPointsMat, H.inv());
convertPointsFromHomogeneous(orgPointsMat, centers);
}
Mat(centers).copyTo(_centers);
return isValid;
}
CirclesGridFinder::findHoles 函数源码
该函数会从前面检测到的斑点中提取圆形网格。
bool CirclesGridFinder::findHoles()
{
switch (parameters.gridType)
{
case CirclesGridFinderParameters::SYMMETRIC_GRID:
{
std::vector<Point2f> vectors, filteredVectors, basis;
Graph rng(0);
//计算所有邻居点间的距离向量，并在邻居间添加一条边，保存到rng图中
computeRNG(rng, vectors);
//过滤向量
filterOutliersByDensity(vectors, filteredVectors);
std::vector<Graph> basisGraphs;
//将向量聚类为4类，并基于其中的两类生成两个图
findBasis(filteredVectors, basis, basisGraphs);
//从生成的两个图中得到圆形网格点
findMCS(basis, basisGraphs);
break;
}
}
return (isDetectionCorrect());
}
CirclesGridFinder::computeRNG 函数源码
该函数会计算所有邻居间的距离向量。判断两个点A、B是否为邻居，首先计算A、B间的距离 r = dAB ,然后分别以A、B为中心 r 为半径画圆，若两圆的公共区域（也即下图中的红色区域）没有其它点，则这两个点互为邻居，同时将距离向量PA-PB 与 PB-PA保存。
void CirclesGridFinder::computeRNG(Graph &rng, std::vector<cv::Point2f> &vectors, Mat *drawImage) const
{
rng = Graph(keypoints.size());
vectors.clear();
//TODO: use more fast algorithm instead of naive N^3
//判断 i 于 j 是否为邻居
for (size_t i = 0; i < keypoints.size(); i++)
{
for (size_t j = 0; j < keypoints.size(); j++)
{
if (i == j)
continue;
Point2f vec = keypoints[i] - keypoints[j];
double dist = norm(vec);
bool isNeighbors = true;
//检查除 i 与 j 外的所有点
for (size_t k = 0; k < keypoints.size(); k++)
{
if (k == i || k == j)
continue;
double dist1 = norm(keypoints[i] - keypoints[k]);
double dist2 = norm(keypoints[j] - keypoints[k]);
//若以i j 为圆心，i j距离为半径的两圆相交的公共区域内存在一个点， 那么i 与 j 不是邻居
if (dist1 < dist && dist2 < dist)
{
isNeighbors = false;
break;
}
}
// 若i 与 j 是邻居，在i j间增加一条边，同时计算它们间的距离向量
if (isNeighbors)
{
rng.addEdge(i, j);
vectors.push_back(keypoints[i] - keypoints[j]);
}
}
}
}
该函数还可以被优化 ，事实上圆形网格上的邻居点在图像上的距离 d 会有一个取值范围，在图像中只有圆形网格时，d取得最大值 dmax= max(img.width，img.height) / min(circleGrid.width，circleGrid.height)，img.width与img.height表示图像宽高，circleGrid.width与circleGrid.height表示圆形网格横向与竖向的点数目。求点P的邻居时可以先将斑点按照 x坐标排序，然后在 这个范围搜索。
CirclesGridFinder::filterOutliersByDensity 函数源码
对于前面计算的每一个向量，在以向量的终点为中心的窗口内统计向量个数，若数量小于阈值，该向量会被剔除。如下图中的P向量，由于在红色窗口内只有一个向量，则可能会被舍弃；q向量由于窗口内有四个向量，则可能会被保留下来。
void CirclesGridFinder::filterOutliersByDensity(const std::vector<Point2f> &samples,
std::vector<Point2f> &filteredSamples)
{
filteredSamples.clear();
for (size_t i = 0; i < samples.size(); i++)
{
// 生成以samples[i]为中心,parameters.densityNeighborhoodSize 大小的矩形
//parameters.densityNeighborhoodSiz为Size类型，默认值为(16,16)
Rect_<float> rect(samples[i] - Point2f(parameters.densityNeighborhoodSize) * 0.5,
parameters.densityNeighborhoodSize);
int neighborsCount = 0;
//统计矩形内向量的个数
for (size_t j = 0; j < samples.size(); j++)
{
if (rect.contains(samples[j]))
neighborsCount++;
}
//当个数>=parameters.minDensity时保留下来, parameters.minDensity默认值为10
if (neighborsCount >= parameters.minDensity)
filteredSamples.push_back(samples[i]);
}
}
CirclesGridFinder::findBasis 源码
该函数执行后basis.size()与basisGraphs.szie()都为2，basis[0]值为{x=42.2135696 y=-0.125327706 }，basis[1]的值为{x=-0.0775655136 y=41.9493866 }。basisGraphs[0]与basisGraphs[1]建立的图可以用下面的图来近似模拟。
void CirclesGridFinder::findBasis(const std::vector<Point2f> &samples,
std::vector<Point2f> &basis, std::vector<Graph> &basisGraphs)
{
basis.clear();
Mat bestLabels;
TermCriteria termCriteria;
Mat centers;
const int clustersCount = 4;
//将向量聚类为4类
kmeans(Mat(samples).reshape(1, 0), clustersCount, bestLabels, termCriteria, parameters.kmeansAttempts,
KMEANS_RANDOM_CENTERS, centers);
std::vector<int> basisIndices;
//TODO: only remove duplicate
//对于每一个聚类中心(x,y)，若x y中绝对值最大的那个变量大于0，则将该聚类中心加入到basis
for (int i = 0; i < clustersCount; i++)
{
int maxIdx = (fabs(centers.at<float> (i, 0)) < fabs(centers.at<float> (i, 1)));
if (centers.at<float> (i, maxIdx) > 0)
{
Point2f vec(centers.at<float> (i, 0), centers.at<float> (i, 1));
basis.push_back(vec);
basisIndices.push_back(i);
}
}
//若保留下来的聚类中心不只两个(正确情况下是保留2个聚类中心，且他们接近垂直)
if (basis.size() != 2)
CV_Error(0, "Basis size is not 2");
//将x取值较大的聚类中心放在vector的第一个位置
if (basis[1].x > basis[0].x)
{
std::swap(basis[0], basis[1]);
std::swap(basisIndices[0], basisIndices[1]);
}
const float minBasisDif = 2;
//两个聚类中心距离要>=2
if (norm(basis[0] - basis[1]) < minBasisDif)
CV_Error(0, "degenerate basis" );
//将属于basis[0]类的向量放在cluster[0]中，属于basis[1]类的向量放在cluster[1]中
std::vector<std::vector<Point2f> > clusters(2), hulls(2);
for (int k = 0; k < (int)samples.size(); k++)
{
int label = bestLabels.at<int> (k, 0);
int idx = -1;
if (label == basisIndices[0])
idx = 0;
if (label == basisIndices[1])
idx = 1;
if (idx >= 0)
{
clusters[idx].push_back(basis[idx] + parameters.convexHullFactor * (samples[k] - basis[idx]));
}
}
//计算cluster[0]与cluster[1]的凸包，保存在hull[0]与hull[1]中
for (size_t i = 0; i < basis.size(); i++)
{
convexHull(clusters[i], hulls[i]);
}
//初始化两个图basisGraphs[0]与basisGraphs[1]，图中结点为斑点检测出的所有斑点，无边
basisGraphs.resize(basis.size(), Graph(keypoints.size()));
//向两个图中添加边
for (size_t i = 0; i < keypoints.size(); i++)
{
for (size_t j = 0; j < keypoints.size(); j++)
{
if (i == j)
continue;
Point2f vec = keypoints[i] - keypoints[j];
//如果vec在hull[0]中，则在图basisGraphs[0]中添加一条i到j的边
//如果vec在hull[1]中，则在图basisGraphs[1]中添加一条i到j的边
for (size_t k = 0; k < hulls.size(); k++)
{
if (pointPolygonTest(hulls[k], vec, false) >= 0)
{
basisGraphs[k].addEdge(i, j);
}
}
}
}
}
在函数的后面向basisGraphs[0]与basisGraphs[1]添加边时是遍历任意两个点之间的边，事实上只需要检查所有邻居间的边就足够，由于computeRNG 函数中已经构建了一个邻居图，所以可以通过检查已构建的邻居图中的边来进行优化。
CirclesGridFinder::findMCS 源码
该函数先是在前一步建立的两个图中寻找最长路径并将该路径作为初始行，然后以初始行为基准向两边扩充剩余的行。其扩充过程可以用下图模拟（这里假设最长路径在第一个图中，若在第二个图中，最长路径是竖直）。
void CirclesGridFinder::findMCS(const std::vector<Point2f> &basis, std::vector<Graph> &basisGraphs)
{
holes.clear();
Path longestPath;
/*
在上一步建立的两个图中寻找最长路径，保存在longestPath中，并返回该路径所在图的索引。
Path是一个struct，有四个成员分别是int firstVertex, int lastVertex, int length;
vector<size_t> vertices;前面三个成员默认值为-1。
*/
size_t bestGraphIdx = findLongestPath(basisGraphs, longestPath);
//将寻找出来的最长路径作为初始行。
std::vector<size_t> holesRow = longestPath.vertices;
//如果最长路径长度超过网格最大边的长度，剔除路径中的首尾点。注意，如果patternSize=（7，7），
//holesRow.size()=8，处理后holesRow.size()=6，这样不会对结果造成影响。
while (holesRow.size() > std::max(patternSize.width, patternSize.height))
{
holesRow.pop_back();
holesRow.erase(holesRow.begin());
}
//如果最长路径在第一个图中，则在第二个图中向两边扩充行，为何如此做，可以看后续对于addHolesByGraph函数源码的分析。
if (bestGraphIdx == 0)
{
//holes是一个向量，存放最终提取到的网格点。
holes.push_back(holesRow);
size_t w = holes[0].size();
size_t h = holes.size();
//扩充行的最小置信度，如果新扩充的行的置信度大于该值，将该行加入到最终网格中。parameters.existingVertexGain=10000
parameters.minGraphConfidence = holes[0].size() * parameters.existingVertexGain;
//扩充patternSize.height-1行。在basisGraphs[1]中寻找。
for (size_t i = h; i < patternSize.height; i++)
{
addHolesByGraph(basisGraphs, true, basis[1]);
}
//如果w小于patternSize.width,则需要扩充列。在basisGraphs[0]中寻找列。
parameters.minGraphConfidence = holes.size() * parameters.existingVertexGain;
for (size_t i = w; i < patternSize.width; i++)
{
addHolesByGraph(basisGraphs, false, basis[0]);
}
}
//如果最长路径在第二个图中，则在第一个图中向两边扩充行
else
{
holes.resize(holesRow.size());
for (size_t i = 0; i < holesRow.size(); i++)
holes[i].push_back(holesRow[i]);
size_t w = holes[0].size();
size_t h = holes.size();
parameters.minGraphConfidence = holes.size() * parameters.existingVertexGain;
for (size_t i = w; i < patternSize.width; i++)
{
addHolesByGraph(basisGraphs, false, basis[0]);
}
parameters.minGraphConfidence = holes[0].size() * parameters.existingVertexGain;
for (size_t i = h; i < patternSize.height; i++)
{
addHolesByGraph(basisGraphs, true, basis[1]);
}
}
}
findLongestPath 源码分析
函数搜寻一条最长路径，并返回该路径所在图的索引，若有多条路径长度相等，则输出最先搜到的路径。若最长路径在第一个图中，函数保证路径从左到右；若在第二个图中，函数保证路径从上到下。
size_t CirclesGridFinder::findLongestPath(std::vector<Graph> &basisGraphs, Path &bestPath)
{
std::vector<Path> longestPaths(1);
std::vector<int> confidences;
size_t bestGraphIdx = 0;
const int infinity = -1;
//遍历所有的图(2个)
for (size_t graphIdx = 0; graphIdx < basisGraphs.size(); graphIdx++)
{
const Graph &g = basisGraphs[graphIdx];
Mat distanceMatrix;
//使用Floyd-Warshall算法求任意两点间的最短距离
g.floydWarshall(distanceMatrix, infinity);
Mat predecessorMatrix;
//计算i到j的最长路径中j前面的一个点。
computePredecessorMatrix(distanceMatrix, (int)g.getVerticesCount(), predecessorMatrix);
double maxVal;
Point maxLoc;
//计算最大路径的长度maxVal，maxLoc.x记录路径的起始点，maxLoc.y记录路径的终点。
minMaxLoc(distanceMatrix, 0, &maxVal, 0, &maxLoc);
//新的最长路径最长，则将longestPaths清空。
if (maxVal > longestPaths[0].length)
{
longestPaths.clear();
confidences.clear();
bestGraphIdx = graphIdx;
}
//将新的最长路径加入到longestPaths中
if (longestPaths.empty() || (maxVal == longestPaths[0].length && graphIdx == bestGraphIdx))
{
Path path = Path(maxLoc.x, maxLoc.y, cvRound(maxVal));
size_t id1 = static_cast<size_t> (maxLoc.x);
size_t id2 = static_cast<size_t> (maxLoc.y);
//根据predecessorMatrix得到最长路径上的所有点
computeShortestPath(predecessorMatrix, id1, id2, path.vertices);
longestPaths.push_back(path);
int conf = 0;
for (int v2 = 0; v2 < (int)path.vertices.size(); v2++)
{
conf += (int)basisGraphs[1 - (int)graphIdx].getDegree(v2);
}
confidences.push_back(conf);
}
}
int maxConf = -1;
int bestPathIdx = -1;
//寻找最好的那条最长路径，事实上longestPaths中只有一条路径。
for (int i = 0; i < (int)confidences.size(); i++)
{
if (confidences[i] > maxConf)
{
maxConf = confidences[i];
bestPathIdx = i;
}
}
bestPath = longestPaths.at(bestPathIdx);
//后续的步骤是保证最长路径是从左到右或者从上到下。
bool needReverse = (bestGraphIdx == 0 && keypoints[bestPath.lastVertex].x < keypoints[bestPath.firstVertex].x)
|| (bestGraphIdx == 1 && keypoints[bestPath.lastVertex].y < keypoints[bestPath.firstVertex].y);
if (needReverse)
{
std::swap(bestPath.lastVertex, bestPath.firstVertex);
std::reverse(bestPath.vertices.begin(), bestPath.vertices.end());
}
return bestGraphIdx;
}
该函数求最长路径时采用的Floyd-Warshall算法，该算法时间复杂度为O(N^3)，太过暴力。事实上可以用图的深度优先遍历来寻找最长路径，遍历时从度为1的点开始遍历即可。
addHolesByGraph源码
该函数向holes（holes中按序保存已找到的网格行）中新增加一行或一列，当addRow为true时，向holes中增加一行；当addRow为false时，向holes中增加一列。
void CirclesGridFinder::addHolesByGraph(const std::vector<Graph> &basisGraphs, bool addRow, Point2f basisVec)
{
std::vector<size_t> above, below, aboveSeeds, belowSeeds;
/*
分别向上和向下扩展。。
若addRow为true，aboveSeeds保存holes()中第一行，above保存向最上扩展的行；否则，aboveSeeds保存holes的第一列，
above保存向最右扩展的列，且有above[index]=aboveSeeds-basisVec。
若addRow为false，belowSeeds保存holes()中最后一行，below保存向最下扩展的行；否则，belowSeeds保存holes的最后一列，
below保存向最左扩展的列，且有above[index]=aboveSeeds+basisVec。
*/
findCandidateHoles(above, below, addRow, basisVec, aboveSeeds, belowSeeds);
//计算扩展above的置信值，该置信值是根据新扩展的点在basisGraphs图中是否有边与其相连来计算。
float aboveConfidence = computeGraphConfidence(basisGraphs, addRow, above, aboveSeeds);
//计算扩展below的置信值。
float belowConfidence = computeGraphConfidence(basisGraphs, addRow, below, belowSeeds);
//根据返回的置信值判断新扩展的点是否可以加入到最终网格中。在两个新扩展都可加入时，只加入above扩展。
//above扩展加入到holes的第一行或第一列前面，below扩展加入到holes的最后一行或最后一列后面。
insertWinner(aboveConfidence, belowConfidence, parameters.minGraphConfidence, addRow, above, below, holes);
}
rectifyGrid 函数源码
Mat CirclesGridFinder::rectifyGrid(Size detectedGridSize, const std::vector<Point2f>& centers,
const std::vector<Point2f> &keypoints, std::vector<Point2f> &warpedKeypoints)
{
//标准网格的圆点中心距为30
const float edgeLength = 30;
const Point2f offset(150, 150);
std::vector<Point2f> dstPoints;
//计算 (centers[centers.size()-centers[0])与(centers[detectedGridSize.width - 1]-centers[0])的向量积
bool isClockwiseBefore =
getDirection(centers[0], centers[detectedGridSize.width - 1], centers[centers.size() - 1]) < 0;
int iStart = isClockwiseBefore ? 0 : detectedGridSize.height - 1;
int iEnd = isClockwiseBefore ? detectedGridSize.height : -1;
int iStep = isClockwiseBefore ? 1 : -1;
//生成标准网格上的点
for (int i = iStart; i != iEnd; i += iStep)
{
for (int j = 0; j < detectedGridSize.width; j++)
{
dstPoints.push_back(offset + Point2f(edgeLength * j, edgeLength * i));
}
}
//计算单应性矩阵
Mat H = findHomography(centers, dstPoints, RANSAC)；
if (H.empty())
{
H = Mat::zeros(3, 3, CV_64FC1);
warpedKeypoints.clear();
return H;
}
std::vector<Point2f> srcKeypoints;
//keypoints保存的是斑点检测的点
for (size_t i = 0; i < keypoints.size(); i++)
{
srcKeypoints.push_back(keypoints[i]);
}
Mat dstKeypointsMat;
//将斑点检测到的所有点投影到标准网格上
transform(srcKeypoints, dstKeypointsMat, H);
std::vector<Point2f> dstKeypoints;
//将Mat点集转为vector<Point2f> 点集
convertPointsFromHomogeneous(dstKeypointsMat, dstKeypoints);
warpedKeypoints.clear();
for (size_t i = 0; i < dstKeypoints.size(); i++)
{
Point2f pt = dstKeypoints[i];
warpedKeypoints.push_back(pt);
}
return H;
}
结尾：从函数源码我们可以知道该函数并没有考虑偏心误差。
文章出处登录后可见！
已经登录？立即刷新
