用 hough 变换检测抛物线
告诉你某张图中很可能有某个东西长得很像开口向上的抛物线，想把它检测出来，嘿嘿，看上去很难吧？
其实呢，这种事情很适合用 hough 变换 来做：
假设抛物线方程为 y = a ( x - b )2 + c
先把图像边缘检测出来，对于图像边缘上的每一个点 (x,y), 可以有很多条抛物线穿过它，即有很多个 a, b, c 可以使上面的方程成立，于是， 你很暴力的把这些 a, b, c 都记下来！ ，对于每个 x,y 点都记下一次，然后看看哪些 a, b, c 的组合被命中的次数比较多就好了。
当然，上面的那个过程有不少地方是可以优化的，比如说对于图片，b 和 c 肯定是整数值，而且 b 很可能在比较偏中心的位置，而 a 的取值范围可以预先设定，比如说如果你想检测下巴，保守的说 a 肯定就是小于 1 的（不然那下巴也太尖了吧！），如此这般等等——而 OpenCV 中用 hough 变换检查圆的代码中显然还有些其他的不那么容易理解的优化，反正我也用不着把它做得很高效，就没有研究了。
下面上图：
（图片均来源于 pixdaus.com ，如果碰巧它侵犯了您的版权，请与我联系）
再上代码——
（你确定要看吗？？！！）
（警告！魔数出没！）
（警告！没有注释！）
（警告！诡异的变量名出没！）
（警告！你的大脑可能因为看到这种乱七八糟的代码受到永久性的损伤！！）
（你确定依然要看吗？？！！！！）
// y = a/150*(x-b)^2+c;
struct Parabola {
int a;
int b;
int c;
Parabola():a(0),b(0),c(0){}
Parabola(int a, int b, int c):a(a),b(b),c(c){}
Parabola(Parabola const& p):a(p.a),b(p.b),c(p.c){}
};
struct Record {
Parabola para;
int count;
Record():para(),count(0){}
Record(int a, int b, int c, int cnt):para(a,b,c),count(cnt){}
};
struct RecordCmp {
inline bool operator()(Record const& a, Record const& b) {
return a.count>b.count;
}
};
std::vector<Parabola> detectParabolas( cv::Mat const& img,
int minA=-1,
int maxA=-1,
int minB=-1,
int maxB=-1,
int minC=-1,
int maxC=-1 )
{
cv::Mat edges;
// cv::imshow("debug", img);
// cv::waitKey(3000);
cv::Canny( img, edges, 50, 150 );
cv::Mat_<uint8_t> img_(edges);
// cv::imshow("debug", edges);
// cv::waitKey(3000);
int ***acc = 0;
if(minA < 0)
minA = 1;
if(minB < 0)
minB = img.cols/4;
if(minC < 0)
minC = img.rows/4;
if(maxA < 0)
maxA = 60;
if(maxB < 0)
maxB = img.cols/2+img.cols/4;
if(maxC < 0)
maxC = img.rows;
acc = new int**[maxA+1];
for(int i=0; i<maxA+1; ++i) {
acc[i] = new int*[maxB+1];
for(int j=0; j<maxB+1; ++j) {
acc[i][j] = new int[maxC+1];
for(int k=minC; k<maxC+1; ++k) {
acc[i][j][k]=0;
}
}
}
printf("calculating ...\n");
int percent=0;
for(int y=0; y<img.rows; ++y) {
for(int x=0; x<img.cols; ++x) {
if(img_(y,x)) {
for(int a=minA; a<=maxA; ++a){
for(int b=minB; b<maxB; ++b) {
int c=y+a*(x-b)*(x-b)/150;
if(c>=minC && c<=maxC) {
++acc[a][b][c];
}
}
}
}
}
int p=int(round(double(y)/img.rows*100.));
if(p!=percent) {
percent = p;
printf("> %02d%%\n", percent);
}
}
printf("done!\n");
sizedpq<Record, RecordCmp> pq(3);
for(int i=minA; i<=maxA; ++i) {
for(int j=minB; j<=maxB; ++j) {
for(int k=minC; k<=maxC; ++k) {
pq.insert(Record(-i,j,k, acc[i][j][k]));
}
}
}
for(int i=0; i<=maxA; ++i) {
for(int j=0; j<=maxB; ++j) {
delete[] acc[i][j];
}
delete[] acc[i];
}
delete[] acc;
std::vector<Parabola> res;
for(int i=0; i<3; ++i) {
res.push_back(pq[i].para);
}
return res;
}
上面那段代码最诡异而需要说明的地方是 a ，我希望检测的是开口向上的抛物线，于是在计算机图像的坐标系统下 a 的取值应该是负数，但作为数组下标时它当然得要是正数，所以这里多绕了几个弯，下次我自己看到肯定也会忘记的，不过还好这算法相当简单，下次需要用到我肯定更愿意重写。
sizedpq 类是俺以前做马赛克拼接图像时折腾出来的一个固定大小的优先队列类，其实设计得并不好，不过接口挺好用，所以.. 就拿来用了。
PS. 速度与图像边缘数量关系很密切，不过通常情况下处理一张上面展示的那种大小的图像也就需要 3 秒钟左右吧，其实还可以接受。
PPS. 人果然还是逼出来的啊，今早老宋很客气的问我节日快乐，弄得我很不好意思的在两个小时之内把这个做好了……
PPPS. 上周末和室友去西湖边逛了一圈，在柳浪闻莺景区里，我刚拍完一张照，发现旁边飘过了一个长发白衣飘逸出尘的美女背影，我说哇美女，他说其实不是美女我已经观察很久了...... 唉，男人。
PPPPS. 谁能教我玩 DotA 啊！！！
