前言

好吧，我是强迫症发作来刷出在感的，因为我发现当年写博文的时候，乱开Categories，导致现在博客很多Categories下只有一篇博文，强迫症患者表示，必须至少两篇，所以我最近就来慢慢填这些坑好了。。所以今天填的是第一个坑——计算机视觉。。。【好吧，其实我不搞DIP很久了。。。

作为一个从大四毕设开始就基本没研究过图像处理的人来说，要写一篇算法的科普文，虽然不是不可以，但是你会发现和网上很多人写的差不多【好吧，一般都是别人写的比我好。。】，所以我就决定写一篇关于去年(好吧，前年)研究生数模比赛的某道题的博文好了【一直觉得那道题很适合拿来水一篇博文。。】；

当年赛题发下来时，要做的第一件事自然就是选题，我把几道题都瞄了一下，然后基本马上下定了“嗯，就做这道题吧！！”的决心；因为那道题其中最核心的一问我已经有十足的信心可以“秒杀”了！！

那道题目的大概意思我也就不说了，我说说那核心的一问吧，假如你有以下三幅卫星云图(为了让大家看的比较清楚区别，我将其转为gif了)，然后你就需要估算出各个地方的风速和风向【也叫风矢场或者云导风】；

你觉得这个不好做是吧，但是如果你研究过光流法的话，这个东西会在你看到题目的瞬间马上闪现在你脑海里！！

光流法，简而言之一句话光流法是干嘛的呢？就是： 根据相邻两帧图片，然后算出各个部分的运动方向和速度；说起来好像很厉害的样子。。是吧。。【所以我感觉这道题就是为了光流法而设计的，不知道其他没碰过这个方法的人最后使用什么办法解决这道题的呢？】

光流法

下面“简单”介绍一下光流法原理吧，不用数学就讲清楚。。。应该是不大可能的。。另外，光流法版本诸多，这里只讲最最最最最最简单的那种，懂了这个的话，其他也挺好懂的。。。如果懒得看数学，那么请无视我。。。

光流法的思想就是把整个图片分成一个一个的块，然后计算出这一刻这一小块图出现在下一帧的什么位置，从而获得这一小块的速度和方向，使用的前提条件是以下三个：

• 亮度恒定：每个分块的颜色前后帧之间要恒定；
• 连续性：两帧之间图像相差不大，或者说分块前后帧之间位置偏差不大；
• 一致性：假设每一小块内各个像素点运动方向一致

假设我们用\(H(x,y,t)\)来表示第t帧图片中\((x,y)\)这个像素点的值，那么因为亮度恒定，可以有下面式子成立:

\(H(x,y,t)=H(x+\delta x,y+\delta y,t+\delta t)\)

也就是说这一帧中\((x,y)\)这个点在下一帧中移动到了\((x+\delta x,y+\delta y)\)这个位置；其实这个时候速度就是:

\(v=(v_x,v_y)=(\dfrac{\delta x}{\delta t},\dfrac{\delta y}{\delta t})\)

接下来将\(H(x+\delta x,y+\delta y,t+\delta t)\)做一阶泰勒展开【\(\varepsilon\)为无穷小项】：

\(H(x+\delta x,y+\delta y,t+\delta t)=H(x,y,t)+\dfrac{\partial H}{\partial x}\delta x+\dfrac{\partial H}{\partial y}\delta y+\dfrac{\partial H}{\partial t}\delta t+\varepsilon\)

于是乎就得到：

\(\dfrac{\partial H}{\partial x}\delta x+\dfrac{\partial H}{\partial y}\delta y+\dfrac{\partial H}{\partial t}\delta t=-\varepsilon\)

为了描述简单，我们定义\(\dfrac{\partial H}{\partial x}=H_x\)，\(\dfrac{\partial H}{\partial y}=H_y\)，然后上面式子两边同除\(\delta t\)，再把\(\dfrac{\delta x}{\delta t}\)用速度\(v_x\)代替，\(\dfrac{\delta y}{\delta t}\)用\(v_y\)代替，就得到：

\(H_xv_x+H_yv_y+H_t=-\varepsilon\)

然而我们这里一直思考的都是一个点，我们之前说过，我们是一个小分块一个小分块来考虑的，假设我们考虑的小分块大小为\(a\times a=m\)，那么这m个点每个点的速度都是一样的，于是就可以得到下面m个方程：

\(\left\{\begin{array}{llcl}H_{x1}v_x+H_{y1}v_y+H_{t1}=-\varepsilon_1\\H_{x2}v_x+H_{y2}v_y+H_{t2}=-\varepsilon_2\\\ldots \ldots\\H_{xm}v_x+H_{ym}v_y+H_{tm}=-\varepsilon_m\end{array}\right.\)

这里提醒一点，为什么上面泰勒展开的时候要保留\(\varepsilon\)这一项，原因很简单，如果它是0的话，那么上面的方程岂不是就变成m个方程两个未知数的形式了？！

我们所要计算的就是要让误差最小，也就是令\(\varepsilon_k\)的平方和最小，令：

\(t=\sum\limits_{n=1}^{m}\varepsilon_n^2=\sum\limits_{n=1}^{m}(H_{xn}v_x+H_{yn}v_y+H_{tn})^2\)

为了使\(t\)最小，很简单，也就是对变量偏导为0：

\(\left\{\begin{array}{ll}\dfrac{\partial t}{\partial v_x}=0\\\dfrac{\partial t}{\partial v_y}=0\end{array}\right.\)

本来这里可以直接给出结论的，我还是简单推一下吧。。。
\(\begin{eqnarray*}\dfrac{\partial t}{\partial v_x} & = &\sum\limits_{n=1}^{m}2(H_{xn}v_x+H_{yn}v_y+H_{tn})H_{xn} \\& = & 2\sum\limits_{n=1}^{m}(H_{xn}^2v_x+H_{yn}H_{xn}v_y+H_{tn}H_{xn})\\& = & 2(v_x\sum\limits_{n=1}^{m}H_{xn}^2+v_y\sum\limits_{n=1}^{m}H_{yn}H_{xn}+\sum\limits_{n=1}^{m}H_{tn}H_{xn})=0\end{eqnarray*}\)
\(\begin{eqnarray*}\dfrac{\partial t}{\partial v_y} & = &2(v_x\sum\limits_{n=1}^{m}H_{xn}H_{yn}+v_y\sum\limits_{n=1}^{m}H_{yn}^2+\sum\limits_{n=1}^{m}H_{tn}H_{yn})=0\end{eqnarray*}\)

也就是说要解的方程就是：
\(\left\{\begin{array}{ll}v_x\sum\limits_{n=1}^{m}H_{xn}^2+v_y\sum\limits_{n=1}^{m}H_{yn}H_{xn}+\sum\limits_{n=1}^{m}H_{tn}H_{xn}=0\\v_x\sum\limits_{n=1}^{m}H_{xn}H_{yn}+v_y\sum\limits_{n=1}^{m}H_{yn}^2+\sum\limits_{n=1}^{m}H_{tn}H_{yn}=0\end{array}\right.\)

于是乎，用矩阵装逼就可以写成：

\(\underbrace{\left(\begin{array}{cc}\sum H_x^2 & \sum H_xH_y \\\sum H_xH_y & \sum H_y^2\end{array}\right)}_{A}\left(\begin{array}{cc}v_x \\v_y\end{array}\right)=-\underbrace{\left(\begin{array}{cc} \sum H_xH_t \\\sum H_yH_t\end{array}\right)}_B\)

所以只要A可逆，那么方程有解~

其实我一开始的计划是，用最少的数学来完成这一篇博文的，嘛~算了，虽然我觉得这篇里面的数学都很简单，至少我已经详细到了每一步推导都写出来了。。。

实用结论

我觉得有人估计想跳上面的数学了，那么我这里就直接给出实用方法吧。。上面的那个矩阵中，\(H_x\)就是每个小分块内每个点在x方向上的颜色变化，\(\sum H_x\)自然就是m个点的x方向颜色变化之和；

\(H_y\)自然就是y方向上的变化；\(H_t\)是这个点在前后两帧之间的颜色变化；

使用的时候，我们有两幅图片是吧，然后将其划分成很多个\(a \times a\)的小格，对于第一帧每个格子内的\(a \times a\)个点，计算出上述矩阵方程中的A和B，得到一个二元二次方程：
\(\left\{\begin{array}{ll}A_{11}v_x+A_{12}v_y=-B_1\\A_{21}v_x+A_{22}v_y=-B_2\end{array}\right.\)
这样解出来的\(v_x,v_y\)就是这个小块的偏移速度了~

我们需要做的就是对每个小块都这么算就是了~

结果

下图是我们去年算出来的结果，最后是二等奖。。。

不过还是要说明一下，上面那个算法的出来的结果和下图还是有点差别的，因为我们还做了一些别的滤波处理什么的，但是由于跟本文没关系，所以我就不讲了。。。另外代码一般贪快都是用的OPENCV的API，而API里面不是简单地上面的算法，而是金字塔多尺度的LK光流法。

事实上这一问是这道题最难的一问，但是由于以前接触过光流法和写过相关代码，所以当时这道题直接把以前视频流部分换成读取两张图片，就完成了。。。

代码

#ifdef _DEBUG
#pragma comment ( lib, "cxcore200d.lib" )
#pragma comment ( lib, "cv200d.lib" )
#pragma comment ( lib, "highgui200d.lib" )
#else
#pragma comment ( lib, "cxcore200.lib" )
#pragma comment ( lib, "cv200.lib" )
#pragma comment ( lib, "highgui200.lib" )
#endif

#include "cv.h"
#include "highgui.h"
#include <math.h>
#include <time.h>
#include <iostream>
#include <iomanip>
#include <fstream>
using namespace std;
static const double pi = 3.14159265358979323846;

double square(double a)
{
	return a * a;
}

void allocateOnDemand(IplImage **img,
					  CvSize size,int depth,int channels)
{
	*img = cvCreateImage( size, depth, channels );
}

int main(int argc,char* argv[])
{
	const int corner_num = 4000;

	//8点半云图
	IplImage* frame = cvLoadImage("2030.bmp");
	CvSize frame_size = cvSize(frame->width,frame->height);

	cvNamedWindow("Optical Flow", 0);

	IplImage *frame1_1C = NULL,*frame2_1C = NULL,*eig_image = NULL, 
			 *temp_image = NULL,*pyramid1 = NULL,*pyramid2 = NULL;
	allocateOnDemand( &frame1_1C, frame_size, IPL_DEPTH_8U, 1);
	allocateOnDemand( &frame2_1C, frame_size, IPL_DEPTH_8U, 1);
	allocateOnDemand( &eig_image, frame_size, IPL_DEPTH_32F, 1);
	allocateOnDemand( &temp_image, frame_size, IPL_DEPTH_32F, 1);
	allocateOnDemand( &pyramid1, frame_size, IPL_DEPTH_8U, 1);
	allocateOnDemand( &pyramid2, frame_size, IPL_DEPTH_8U, 1);

	cvConvertImage(frame, frame1_1C, CV_CVTIMG_FLIP);
	//9点云图
	frame = cvLoadImage("2100.bmp");
	cvConvertImage(frame, frame2_1C, CV_CVTIMG_FLIP);
	

	CvPoint2D32f frame1_features[corner_num];
	int number_of_features;
	number_of_features = corner_num;

	//获取特征点
	cvGoodFeaturesToTrack(frame1_1C, eig_image, temp_image,
		                  frame1_features, &number_of_features, 
						  .01, 15, NULL);

	//去除无云区特征点
	for(int i = 0;i < corner_num;i++)
	{
		int w = frame1_features[i].x;
		int h = frame_size.height - frame1_features[i].y;

		if(CV_IMAGE_ELEM(frame,uchar,h,3*w)< 120)
			frame1_features[i].x = frame1_features[i].y = 1;
	}

	//计算光流
	CvPoint2D32f frame2_features[corner_num];
	char optical_flow_found_feature[corner_num];
	float optical_flow_feature_error[corner_num];

	CvSize optical_flow_window = cvSize(64,64);
	CvTermCriteria optical_flow_termination_criteria 
	    = cvTermCriteria(CV_TERMCRIT_ITER|CV_TERMCRIT_EPS,20, .3);
	cvCalcOpticalFlowPyrLK(frame1_1C, frame2_1C, pyramid1,
		                    pyramid2, frame1_features,
							frame2_features,
							number_of_features,optical_flow_window, 
							5,optical_flow_found_feature, 
							optical_flow_feature_error, 
							optical_flow_termination_criteria,0);

	//下面部分为绘图画箭头部分
	int none_zero_vector_num = 0;
	int unstable_vector_num = 0;
	for(int i = 0; i < number_of_features; i++)
	{
		if ( optical_flow_found_feature[i] == 0 )	
			continue;
		int line_thickness;				
		CvScalar line_color;			
		line_color = CV_RGB(255,0,0);

		CvPoint2D32f p,q;
		p.x =  frame1_features[i].x;
		p.y =  frame1_features[i].y;
		q.x =  frame2_features[i].x;
		q.y =  frame2_features[i].y;


		double angle;		
		angle = atan2((double)p.y - q.y, (double)p.x - q.x);
		double hypotenuse;	
		hypotenuse = sqrt(square(p.y - q.y) + square(p.x - q.x));

		if(hypotenuse > 25 || hypotenuse < 0.8) 
		{
			unstable_vector_num++;
			continue;
		}
		none_zero_vector_num++;
		if(hypotenuse*6.1 < 15)//4.61
			line_color = CV_RGB(255,0,0);
		else if(hypotenuse*6.1 < 25)
			line_color = CV_RGB(0,0,255);
		else
			line_color = CV_RGB(0,255,0);

		hypotenuse = 14;
		q.x = (p.x - 4 * hypotenuse * cos(angle));
		q.y = (p.y - 4 * hypotenuse * sin(angle));

		//draw arrow	
		cvLine(frame, cvPoint(p.x,p.y), cvPoint(q.x,q.y),
			   line_color, 2, CV_AA, 0 );
		p.x = (int) (q.x + 9 * cos(angle + pi / 4));
		p.y = (int) (q.y + 9 * sin(angle + pi / 4));
		cvLine(frame, cvPoint(p.x,p.y), cvPoint(q.x,q.y), 
			   line_color, 2, CV_AA, 0 );
		p.x = (int) (q.x + 9 * cos(angle - pi / 4));
		p.y = (int) (q.y + 9 * sin(angle - pi / 4));
		cvLine(frame, cvPoint(p.x,p.y), cvPoint(q.x,q.y), 
			   line_color, 2, CV_AA, 0 );
	}
	cvConvertImage(frame, frame, CV_CVTIMG_FLIP);
	cvShowImage("Optical Flow", frame);
	cvSaveImage("temp1.bmp",frame);
	cvWaitKey(0);

	return 0;
}

【完】

本文内容遵从CC版权协议，转载请注明出自http://www.kylen314.com