首页 > 计算机视觉 > 卫星云图定风向——光流法

卫星云图定风向——光流法

2014年1月11日 发表评论 阅读评论

前言

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

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

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

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

of

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

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

光流法

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

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

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

假设我们用\(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光流法。
OF result

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

代码


【完】

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

  • 知其然,知其所以然,不知其如何然

    • 。。。。你想表达什么。。。【我都不知道要怎么回你了。。。】

      • 简单来说就是以前知道怎么实现,但是现在我高数已经忘到连符号都几乎看不懂的境界了(用境界这两个字很合适啊)

        • 你现实是工作了还是还在读书?

          • 工作好多年了

            • 其实在你那聊VPN余额的时候我已经知道了。。

              • 为什么前面你这条12号的留言没有提醒啊。我是进来找新文章时才看到的

                • 我也好奇为什么我在你那边聊了半天,然后你会突然来回这一条。。。是不是好几条一块提醒然后你手滑关掉了。。

                  • 可能吧。最近鼠标经常单击变双击

                    • 我经常前一天到处逛逛留言,然后第二天一起床一看n条回复。。。总感觉点的时候搞不好会漏掉一些。。还有多说同一个人在一篇博文里面同一个人回复你多条的话,多说只会提示一条,经常就漏看。。

                    • 恩,这是个麻烦。等我主机迁移后,准备使用原生评论。另外你还没有睡觉啊,我已经躺下了

                    • 还在实验室赶论文。。。正准备撤。。。晚安好梦~等我主机搬了可以发邮件之后我肯定也是要把多说撤了的节奏!!

                    • 我现在的主机太慢,必须缓存,因此只能用第三方评论服务。你不能通过smtp发送邮件么

                    • 为什么?缓存会导致评论不能实时显示么?我试过了,smtp插件总是返回连接不上smtp host,我后来就没搞了。。虽然google了各种方法都不行。。虽然官方说了这个主机支持mail函数的。。我也不是很明白。突然想到不用多说的话版聊好像就没那么方便了。。每次要看邮件才知道别人回复了。。

                    • 支持mail的话,可以直接用mail()发送,或者用wp内置的wp_mail()发送。开启了缓存自然不能实时显示了。缓存和实时是个冲突的概念

                    • 我当时弄了一阵子没弄出来,就放弃转多说了。。好吧。。。。缓存不能局部缓存是吧。。

                    • 缓存时把整个页面另存一份html静态文件,然后通过.htaccess规则调用

                    • 哦,好吧。。我也想弄缓存。。但是缓存插件不用用于破WIN主机。。所以我要换主机!!

  • 跟用 BMA 之类的做 Motion Evaluation 好像_________________

    • BMA是用于做视频压缩那个的么?后面下斜杠是什么意思。。。我看你多说历史评论好多都是加这么一下杠的

      • 块匹配,在运动向量评估中用到过,就是描述前后帧物体的运动,类似于你 PO 的这篇说的,视频压缩应该可以用这个_____________________(ps下划线是我的习惯,可以说癖好,以前主要是应对百度贴吧的15字,^_^)_______________

        • 我瞄了一下块匹配的文献,好像主要是用于视频压缩的吧;感觉上块匹配和光流相比,块匹配应该还具有预测功能,根据本帧图片和运动向量应该可以大致预测出下一帧长啥样,但是光流法就完全不行。。。好吧。。。就和我打省略号一个性质。。。只不过下划线很少见。。。

          • (*^__^*) 嘻嘻……

          • 一看就是动漫控,so am i____________

            • 看头像懂PO主系列~

  • 逆天啊啊啊啊!!!!这神一样的公式不用图片肿么做到的。。

    • 其实对于行内人士看来这就是一篇没啥水准的水文。。公式的话,可以用插件LaTeX for WordPress,或者Mathjax Latex,我现在用的是后者但是实际上也有不用插件更简单的方法,在你需要公式的地方插入这篇里面的代码就可以了。。http://www.cnblogs.com/ilogic/archive/2012/08/05/latex.html

    • 我才发现你的头像不是静态图片。。

  • 高科技,看不懂

    • 水文,别在意~换头像啦?

      • 你的水文都能写这么多字,我的水文才两三行字

        • 我不习惯写短文。。。。短文就感觉变微博了。。所以写出来都是有点长。。最后导致没人看。。。

        • 虽然我自己不习惯写短文,但是我很喜欢看别人短文。。看几秒就可以留言了~

          • 我也是看到很长的文就头晕

            • 长文是写给有这方面需求的人细细研究用的。。我看你微博说是室内设计师?那你可以写这方面的博文啊~

              • 。。。谁说室内设计的就能写博文??

                • 总有心得可以写的嘛。。。

                  • 这个方面的心得写了没意思

                    • 好吧~那就当一个折腾的地方好了~

  • 看不懂

  • 不明觉得厉害

  • 基本思想还是可以理解的

    • 其实仔细慢慢看的话,我感觉我应该说清楚了。。

  • 好深奥的样子

    • 个人觉得仔细看下来的话应该不难理解。。因为我推导得很详细了。。

  • 讲的好清晰,看着就是一种享受!

    • 谢谢哈,这种评价听着好开心!!(逃

  • 看了好文章江南留个脚印是美德 http://www.timi520.com 三亚婚纱摄影哪家好

  • 第一次来楼主博客,欢迎互动我的小站,呼叫中心www.fuziseo.com

  • Yode

    求原图,我想试一下