opencv3中camshift详解（三）opencv库中CamShift()函数源码探究
在opencv3中camshift详解（一）中，已经对opencv自带的camshiftdemo进行了分析，对camshift的实际使用有了感性认识。
在opencv3中camshift详解（二）中，已经对论文中的camshift原理进行了梳理，终于到了最后一部分，能够进到CamShift函数源码中一探究竟了。demo中使用的camshift函数的具体实现源码位置在…\opencv\sources\modules\video\src\camshift.cpp，3.2版本的实现代码如下：可以先粗略地看一下，重点并不是代码实现。
#include "precomp.hpp"
int cv::meanShift( InputArray _probImage, Rect& window, TermCriteria criteria )
{
CV_INSTRUMENT_REGION()
Size size;
int cn;
Mat mat;
UMat umat;
bool isUMat = _probImage.isUMat();
if (isUMat)
umat = _probImage.getUMat(), cn = umat.channels(), size = umat.size();
else
mat = _probImage.getMat(), cn = mat.channels(), size = mat.size();
Rect cur_rect = window;//保存当前搜索窗window
CV_Assert( cn == 1 );//检测图像通道数是不是1，若不是会返回错误信息
//搜索窗面积为负会报错
if( window.height <= 0 || window.width <= 0 )
CV_Error( Error::StsBadArg, "Input window has non-positive sizes" );
window = window & Rect(0, 0, size.width, size.height);//确保window在图像矩阵内
double eps = (criteria.type & TermCriteria::EPS) ? std::max(criteria.epsilon, 0.) : 1.;
eps = cvRound(eps*eps);
int i, niters = (criteria.type & TermCriteria::MAX_ITER) ? std::max(criteria.maxCount, 1) : 100;
for( i = 0; i < niters; i++ )
{
cur_rect = cur_rect & Rect(0, 0, size.width, size.height);
if( cur_rect == Rect() )
{ //如果当前窗口是空的，则把它默认放在图像中间
cur_rect.x = size.width/2;
cur_rect.y = size.height/2;
}
cur_rect.width = std::max(cur_rect.width, 1);
cur_rect.height = std::max(cur_rect.height, 1);
Moments m = isUMat ? moments(umat(cur_rect)) : moments(mat(cur_rect));//计算图像的矩
// 计算质心
if( fabs(m.m00) < DBL_EPSILON )
break;
//计算质心到搜索窗口中心的距离向量（漂移向量）
//值得注意的是，此处的质心坐标是以cur_rect（ROI）为标准的，也就是说，以cur_rect左上角为原点
//而不是以整个图像为标准，从moments(mat(cur_rect)也可以看出，质心计算只和cur_rect有关系
//cur_rect.x和cur_rect.y是以整个图像为标准，两套坐标轴不要搞混了
int dx = cvRound( m.m10/m.m00 - window.width*0.5 );
int dy = cvRound( m.m01/m.m00 - window.height*0.5 );
//若cur_rect.x + dx<0则cur_rect左上角的横坐标nx到图片外了，要限制nx=0；就是这个表达式max(cur_rect.x + dx, 0)
//若cur_rect.x + dx>0,则要保证cur_rext左上角的横坐标nx加上cur_rect.width后不能超过图片的宽度size.width，即nx要小于size.width - cur_rect.width
//两个边界条件合成一个表达式了
int nx = std::min(std::max(cur_rect.x + dx, 0), size.width - cur_rect.width);
int ny = std::min(std::max(cur_rect.y + dy, 0), size.height - cur_rect.height);
//因为可能遇上边界条件，再一次计算漂移向量
dx = nx - cur_rect.x;
dy = ny - cur_rect.y;
cur_rect.x = nx;
cur_rect.y = ny;
// 计算精度是否满足条件，满足了就break
if( dx*dx + dy*dy < eps )
break;
}
window = cur_rect;//windows用的是Rect&，是引用；更新搜索窗
return i;
}
cv::RotatedRect cv::CamShift( InputArray _probImage, Rect& window,
TermCriteria criteria )
{
//应该是OpenCV相关算法表现性能测试框架，测量函数执行时间，在函数内部追踪函数执行状况
CV_INSTRUMENT_REGION()
const int TOLERANCE = 10;
Size size;
Mat mat;
UMat umat;
bool isUMat = _probImage.isUMat();
if (isUMat)
umat = _probImage.getUMat(), size = umat.size();
else
mat = _probImage.getMat(), size = mat.size();
meanShift( _probImage, window, criteria );//使用meanshift
window.x -= TOLERANCE;
if( window.x < 0 )
window.x = 0;
window.y -= TOLERANCE;
if( window.y < 0 )
window.y = 0;
window.width += 2 * TOLERANCE;
if( window.x + window.width > size.width )
window.width = size.width - window.x;
window.height += 2 * TOLERANCE;
if( window.y + window.height > size.height )
window.height = size.height - window.y;
// 计算新位置上的搜索框图像的矩
Moments m = isUMat ? moments(umat(window)) : moments(mat(window));
double m00 = m.m00, m10 = m.m10, m01 = m.m01;//一阶原点矩
double mu11 = m.mu11, mu20 = m.mu20, mu02 = m.mu02;//二阶中心矩，注意还有几何矩，注意区别，此处的u代表拉丁字母 https://blog.youkuaiyun.com/keith_bb/article/details/70197104
//m11,m20这些不带u的是几何矩
//目标矩形质量太小就退出
if( fabs(m00) < DBL_EPSILON )
return RotatedRect();
double inv_m00 = 1. / m00;//m00取倒数，1.表示double(1.0),这里并不是matlab里的点除
int xc = cvRound( m10 * inv_m00 + window.x );//因为计算区域大了一些，要算目标质心的新位置x centre，此处质心以整幅图像为坐标
int yc = cvRound( m01 * inv_m00 + window.y );
double a = mu20 * inv_m00, b = mu11 * inv_m00, c = mu02 * inv_m00;
//
double square = std::sqrt( 4 * b * b + (a - c) * (a - c) );
// 计算目标主轴方向角度
double theta = atan2( 2 * b, a - c + square ); //theta是与x轴的夹角
//实际上和论文的公式是一样的
// 计算宽度和长度
double cs = cos( theta );
double sn = sin( theta );
double rotate_a = cs * cs * mu20 + 2 * cs * sn * mu11 + sn * sn * mu02;//转动惯量
double rotate_c = sn * sn * mu20 - 2 * cs * sn * mu11 + cs * cs * mu02;
//https://fr.wikipedia.org/wiki/Camshift
double length = std::sqrt( rotate_a * inv_m00 ) * 4;
double width = std::sqrt( rotate_c * inv_m00 ) * 4;
// 若长度小于宽度，交换两个值，对角度值也要处理
if( length < width )
{
std::swap( length, width );
std::swap( cs, sn );
theta = CV_PI*0.5 - theta;
}
//根据length和width的大小对window进行调整，因为meanshift只能接受横平竖直的矩形
// 储存结果
int _xc = cvRound( xc );//xc既是矩形质心坐标，也是矩形中心坐标，因为meanshift已经把矩形中心移动到质心位置了
int _yc = cvRound( yc );
//可以尝试以下函数
//返回包含旋转矩形的最小矩形
//Rect boundingRect() const;
int t0 = cvRound( fabs( length * cs ));//长轴在X轴的投影长度
int t1 = cvRound( fabs( width * sn ));//短轴在X轴的投影长度
//选择长短轴在x轴上的投影中最长的一个加上2作为矩形的宽度
t0 = MAX( t0, t1 ) + 2;
window.width = MIN( t0, (size.width - _xc) * 2 );
t0 = cvRound( fabs( length * sn ));//长轴在Y轴的投影长度
t1 = cvRound( fabs( width * cs ));//短轴在Y轴的投影长度
//选择长短轴在Y轴上的投影中最长的一个加上2作为矩形的高度
t0 = MAX( t0, t1 ) + 2;
window.height = MIN( t0, (size.height - _yc) * 2 );
window.x = MAX( 0, _xc - window.width / 2 );//保证左上角在图像内
window.y = MAX( 0, _yc - window.height / 2 );
window.width = MIN( size.width - window.x, window.width );//保证右下角在图像内，从质心延伸出的长度不能超过边界，最大是size.width-x，画图比较好理解
window.height = MIN( size.height - window.y, window.height );
RotatedRect box;
box.size.height = (float)length;
box.size.width = (float)width;
box.angle = (float)((CV_PI*0.5+theta)*180./CV_PI);//顺时针旋转确实要+90度
while(box.angle < 0)
box.angle += 360;
while(box.angle >= 360)
box.angle -= 360;
if(box.angle >= 180)
box.angle -= 180;
box.center = Point2f( window.x + window.width*0.5f, window.y + window.height*0.5f);
return box;
}
/* 结束 */
代码已经写了一些注释注释，理解起来应该会容易一些，先把代码实现放到一边，来研究下CamShift()函数中对下一帧搜索框的确定。
由代码可以发现和论文中的一样，CamShift()函数先进行了meanshift()，再利用和论文不太一样的公式计算了主轴角度，长轴长度、短轴长度，之后再根据这些数据确定了下一帧搜索框的大小和返回的RotateRect box。
（1）主轴角度的计算
代码：
原理：
来自于camshift详解（二）中的参考文献图像处理基础知识（二）—— 中心矩求主轴方向
可以看出来和论文中求主轴角度的公式有所不同，但本质其实是一样的，无非是论文使用的是原点矩表示，代码中使用的是中心矩，感兴趣的可以自己化简一下，我化简到一半已经一头雾水了哈哈。
（2）主轴长短轴的长度计算
代码：
原理：
不意外地发现和camshift论文中长短轴计算公式又不一样，毕竟a、b、c的定义都不一样。经过一段时间地查找资料后，我认为代码中求主轴长度length和width的公式实际上是利用了转动惯量作为桥梁。
在参考文献基于几何代数理论的医学图像配准研究的第54页中，写有不规则物体的绕任意轴旋转的转动惯量公式：
可以发现式子（3-17）已经很接近代码中所使用的求rotate_a的公式了，但是还是有一个正负号的差别，原因是参考文献中y轴是向上的，而在图像矩阵中，已经提到过很多次了，建系是y轴向下的，因此表示角度的绝对值相等，而正负性相反。将（-α）代入参考文献的公式，就得到了代码中求绕长轴转动的转动惯量的公式。
再将（1）主轴角度的计算求得的角度θ代入转动惯量公式，就得到了转动惯量rotate_a，而长轴和短轴是垂直关系，角度相差90°，绕短轴的转动惯量只要把公式中的角度加上90°就行，得到绕短轴旋转的转动惯量totate-c。
参考文献中的求转动惯量的公式是针对不规则物体绕任意轴旋转的通用公式，而实际上我们求的是规则的椭圆绕自身长轴短轴旋转的结果，有没有更加简洁的，不依赖于角度和中心矩的公式呢？答案是肯定的，
椭圆绕长轴旋转的转动惯量的公式为：其中a为长半轴长度，b为短半轴长度，k为密度。
I y = 1 4 ∗ a 3 ∗ b ∗ p i ∗ k I_y=\frac{1}{4}*a^3*b*pi*k Iy=41∗a3∗b∗pi∗k
公式来源为椭圆转动惯量计算
a* b* pi为椭圆的面积，乘以密度k之后就是椭圆的重量，椭圆的重量是什么呢？就是零阶矩M00啊！M00的公式不就是求全部像素点质量之和吗？（回头来看这句话不太严谨，这两者相等是巧合还是有什么物理数学原理在里面，我已经不清楚了，网络上关于二阶矩椭圆表示图像的内容也没搜索到多少）
转动惯量rotate_a已知，0阶矩M00已知，就可以求得长半轴长度a了。
a = 2 ∗ r o t a t e a M 00 a=2*\sqrt{\frac{rotatea}{M_{00}}} a=2∗M00rotatea
a再乘2就得到了椭圆长轴的长度length，同理短轴的长度也可以求出来。
（3）计算新搜索窗大小和旋转矩形
代码中已经注释过了，因为meanshift()只能接受横平竖直的矩形，所以将长短轴分别投影于X轴和Y轴上，取较长的值作为新的搜索窗的长和宽。
接下来是旋转矩形，矩形高度为长轴长度length，宽度为短轴长度width就不说了，为什么旋转角度弧度转角度后要加上90°呢？
旋转矩形的产生可以分解，第一步是画出未经旋转的矩形，从x轴出发顺时针旋转，先碰到的轴默认作为width，第二步是旋转。既然将椭圆短轴长width赋给旋转矩形宽度width，将椭圆长轴长length赋给旋转矩形高度height，为了能准确表示椭圆，确实是要多旋转90°的。下面画了一张图帮助理解。
最后旋转角度被限制在了0-180°之间，和弧度制的（-pi/2，pi/2）的范围也是相吻合的。opencv里没有绘制旋转矩形的函数，所以在demo主程序中又绘制了矩形的内接椭圆，按道理来讲这个绘制的椭圆就是就是上文计算出的图像椭圆。
遗留问题和总结
camshift到此应该告一段落了，代码实现，论文原理，代码原理都走了一遍，我觉得已经没什么决定性的问题了，还遗留了一些“知其然不知其所以然”的问题，比如：
（1）下一帧搜索窗的大小为什么要牵扯上二阶矩，用一系列投影来确定下一帧搜索框，论文的零阶矩有什么不合适吗？
（2）图像椭圆为什么就能覆盖追踪的目标呢？虽然能理解二阶中心矩（方差）能表示离平均值的距离远近，但是对于由二阶中心矩求出来的椭圆长短轴能覆盖目标这件事，我还是希望有更确信的解答。
（3）图像椭圆的面积和M00是相等关系吗？
其他没什么好总结的了，只是想到这个算法是二十多年前的了就有种无力感，就算搞懂了也已经落后于时代了，现在是深度学习的天下了。还有深入去探究具体实现真的很耗时间和精力，如果是为了做工程的话有点得不偿失。
扫一扫
1527
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
