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

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

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

  1. 2014年1月12日06:33 | #1

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

    • 2014年1月12日13:35 | #2

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

      • 2014年1月12日18:39 | #3

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

        • 2014年1月12日19:16 | #4

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

          • 2014年1月13日16:59 | #5

            工作好多年了

            • 2014年1月13日17:36 | #6

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

              • 2014年1月13日20:24 | #7

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

                • 2014年1月13日20:26 | #8

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

                  • 2014年1月13日22:09 | #9

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

                    • 2014年1月13日22:11 | #10

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

                    • 2014年1月13日22:19 | #11

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

                    • 2014年1月13日22:25 | #12

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

                    • 2014年1月13日22:35 | #13

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

                    • 2014年1月13日23:07 | #14

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

                    • 2014年1月14日09:18 | #15

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

                    • 2014年1月14日14:02 | #16

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

                    • 2014年1月14日14:48 | #17

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

                    • 2014年1月14日15:02 | #18

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

  2. 2014年1月12日08:44 | #19

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

    • 2014年1月12日13:43 | #20

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

      • 2014年1月12日13:56 | #21

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

        • 2014年1月12日14:04 | #22

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

  3. 2014年1月12日09:34 | #26

  4. 2014年1月12日16:34 | #28

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

    • 2014年1月12日16:39 | #29

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

    • 2014年1月12日16:46 | #31

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

  5. 2014年1月13日15:46 | #33

    高科技,看不懂

    • 2014年1月13日15:50 | #34

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

      • 2014年1月13日15:52 | #35

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

        • 2014年1月13日15:55 | #36

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

        • 2014年1月13日15:58 | #37

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

  6. 2014年1月13日16:12 | #44

    看不懂

  7. 2014年1月13日16:31 | #45

    不明觉得厉害

  8. 2014年1月13日17:36 | #46

    基本思想还是可以理解的

    • 2014年1月13日17:38 | #47

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

  9. 2014年1月14日22:14 | #48

  10. 2014年1月21日11:20 | #49

    好深奥的样子

    • 2014年1月21日12:01 | #50

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

  11. 2015年4月20日20:11 | #51

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

  12. 2015年4月29日22:43 | #53

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

  13. 2015年5月2日12:24 | #54

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

  14. Yode
    2016年10月31日15:58 | #55

    求原图,我想试一下

验证码:4 + 7 = ?

友情提示:留言可以使用大部分html标签和属性;

添加代码示例:[code lang="cpp"]your code...[/code]

添加公式请用Latex代码,前后分别添加两个$$