PIV解析时用的算法

好吧这个部分比较复杂所以本来没有打算多说的（事实上期中发表的时候这个地方也不是重点，所以时间有限决定跳过），8过既然conanshang问道，这里就简单说一说把w

其实对比相似度的算法有很多，比如PIV通常用的标准“直接相互相关法”的算法的话是这样：Rfg = A / (B * C)，其中：
A = sum[]sum[](f(x + i, y + j) – fbar) * (g(x + ζ + i, y + η + j) – gbar)
B = sqrt(sum[]sum[](f(x + i, y + j) – fbar)2)
C = sqrt(sum[]sum[](g(x + ζ + i, y + η + j) – gbar)2)
fbar = sum[]sum[](f(x + i, y + j)) / (n * m)
gbar = sum[]sum[](g(x + ζ + i, y + η + j)) / (n * m)
（sum[]sum[]是sum[j = 0 to n]sum[i = 0 to m]，相当于for (j = 0; j < n; j++){ for (i = 0; i < m; i++) }}） Rfg是上一篇文章中的红色正方形区域和白色正方形区域的相关度，x, y分别是小区域的起点坐标，对应上一篇文章的红色正方形的左上角坐标；n, m分别是长宽，i, j分别是正方形内坐标；ζ, η分别是白色正方形相对红色正方形的相对坐标偏差，一点点移动ζ和η，得出每一个位置的R值，然后R值最大的位置的ζ和η就是上一篇文章红色正方形范围的流体的移动方向距离

8过这个算法的一个问题是需要进行两次for j for i的loop（第一次算出fbar和gbar，第二次才能算出ABC得出R），所以实际使用的算法是一个简略版的算法：Rfg2 = (A * A) / (B * C)，其中：
A = sum[]sum[](f(x + i, y + j) * g(x + ζ + i, y + η + j))
B = sum[]sum[](f(x + i, y + j)2)
C = sum[]sum[](g(x + ζ + i, y + η + j)2)

这个算法的话相对来说计算量就要小很多，而且最大的方便之处是在于可以只用一遍for循环就可以w

有兴趣的童鞋的话，CUDA的程序咱已经写好了，可以下下来试试w
http://cl.ly/3T1Q3S1w1Z3G
解压之后里面有4个文件：
2ms2000018.raw
2ms2000019.raw
occ.c
occ.cu
上面两个.raw文件是解析用的两张照片的raw文件，occ.c是C语言程序，occ.cu是CUDA程序w

使用方法：
1、确认自己的显卡是N卡并且支持CUDA环境
2、参照网上教程安装CUDA环境（如果已安装，请跳过w）
3、下载上面的文件，解压
4、打开命令工具（Mac或UNIX/Linux的话是Terminal，Windows的话是Command prompt（cmd））
5、cd到被解压路径
6、使用nvcc编译occ.cu（比如nvcc occ.cu -o occ）
7、运行被打包程序（比如./occ）
8、如果运行成功会看到一大堆的summary然后最下面是CUDA success，此时文件夹里面应该会多出一个result.dat文件
9、如果需要显示.dat文件的解析结果，需要安装GNUPlot（如果已安装，请跳过）
10、运行GNUPlot（gnuplot）
11、设置x轴和y轴（set xrange [0:1024]，set yrange [1024:0]）
12、设置图表为正方形（set size square）
13、描图（plot “result.dat” with vector）
这样的话应该就能和看到上一篇文章的那张图表差不多的结果了w

当然如果无法支持CUDA环境（比如A卡）的话，可以使用C语言的版本来测试效果，使用方法和上面差不多，只是不需要CUDA环境，外加使用cc或gcc打包occ.c就行；另外C语言版本只是用来做CUDA移植前的算法和写法测试用的所以和现在的CUDA版有一些微妙的不同（比如运行结束之后没有那么多的summary）

此外如果有兴趣的童鞋可以自己找两张.raw照片，然后修改一下代码里面导入文件的大小和文件名，再有能力的童鞋也可以修改一下其他的参数试试，现在的参数是照片大小1024×1024，Interrogation Window（红色/白色正方形）是16×16，Search Window算起点坐标大小是(8+8)×(8+8)，算范围大小（白色虚线正方形）是（8+16+8)×(8+16+8），Step Size（红色正方形的每一次计算完毕之后的移动距离）是(16/2)×(16/2)，CUDA的Block Size是8×8×1，恩大概就是这样？（只是需要是8bit黑白的.raw文件，因为研究室和实验室用的是8bit黑白raw文件，所以完全没有考虑24bit彩色raw文件等其他的文件格式；另外也同样因为只是CUDA的运行测试所以也没有装switch的option，只能手动改代码里面的文件名和文件大小，相关参数在代码的comment里面都写出来了w）

p.s. 这个代码还有很多问题，比如咱女仆本的9400m的话经常会告诉偶out of memory（哈！？），外加时不时还会死掉ww而且就算不死，运行速度也非常慢（快的时候大约1秒钟，慢的时候和CPU几乎一样要7秒多种w），当然这和9400m的架构太老估计也有关系，比如实验室的GTX 560 Ti的话就没有上述问题，运行起来都是0.5秒就能搞定w当然，上面也说了，里面很多参数（比如Interrogation Window的大小啥的），修改那些参数运行下来的时间也不一样，当然得出来的结果也有可能会不太一样w（因为附送的raw文件没有除噪点，外加偶也没有把求sub pixel移动量和除误差的代码加进去）

再p.s. 一不小心就发现原来距离发表只剩下一个星期了((((；ﾟДﾟ)))

再再p.s. 话说入手了蓝牙的小型joystick——之所以入手这种东西完全是因为咱找不到更好的红外线遥控器的替代品了w恩演讲的时候比起神马该死的激光棒果然还是普通的遥控器最好用了——可是该死的是水果那个红外线遥控器毕竟是红外线的操作的时候有各种局限性，所以一直向找一个蓝牙的——可是居然蓝牙的遥控器还真是没几个……而且一个二个都是要么做成鼠标状（真想不通这种糟糕企划是怎么通过会议的= =）要么做成激光棒= =相对来说只有这货的手感最好而且价格也比较公道了——喂水果乃敢不敢把乃们发布会的时候用的那个clicker也作为一个Accesory发售啊ノ ﾟДﾟ)ノ ＝＝＝＝ ┻━━┻

再再再p.s. 好吧说到水果发布会，今天半夜终于要发布新iPhone了——8过比起那货咱更期待的其实是新iMac——如果出新的绝对让教授给偶换一台XD

再再再再p.s. Wow! Retina!

Author: 星野恵瑠

Mac user, Niji-Ota, Chinese, Now working in Japan at MAGES. Inc., Future's aim is that one day my name can be listed in Wikipedia

15 thoughts on “PIV解析时用的算法”

  3. 前面全部忽略（拖走

    话说为啥不考虑iOS上的keynote remote程序？…那个貌似也是蓝牙wifi支持的而且还能看到幻灯内容.. (OoO)

    Reply

  4. 下载下来膜拜一下，虽然有bbb可以使用N显卡，但是没装CUDA编译环境的说。最后一张图偶要吐槽，画这个图的肯定没用过linux (……)

    Reply

  5. 我的垃圾网那个连接没反应…另,看了大概实现逻辑,感觉要是我做的话会加入个图片分割的逻辑过程,将图片分块然后交给多CPU进行处理,最后feedback. 得到结果汇总.另….这种真心需要个服务器来玩…估计本子的显卡CUDA也没好到哪里去吧…我记得我同学一个GT530M开CUDA视频转码,CPU花了12min,CUDA辅助下花了10min,也不是十分惊艳的结果吧…估计NV Tesla能得出非常好的效果,但是咱真心没见过实物….以上.

    Reply

  6. 又想了下, 看到上面的举例算法描述,浮点运算还是很多的. CUDA也许会有比较理想的效果?
    另,想到另外一个参考就是关于开方的问题,如果涉及到开方运算的话,不要用sqrt(),把常用的开方值写到缓存文件里,比如sqrt(10) = 3.16227766….把这组值按需要精度做个索引表,每次用到的时候去查询,查询不到再sqrt(),这样速度就会得到嗑药般的提升,另外C++啥的SSE优化也有显著效果,不过目前的游戏开发全部倾向于第一种.缓存sqrt()结果.
    又另,这种纯计算的算法的一般的参考都是用FORTRAN编写的….所以效率有天生优势…有没有FORTRAN2CU的编译器啥的…. (;;)

    Reply

  13. @Decmes
    梦梦别这样 (;;)

    @人好哇！
    喂偶不写水果发布会的文很久了啊好吧 (^^;)

    @小明猪
    肥猪啥时候有了偶逢水果必写的错觉 (!!!!)

    @Linfcstmr
    问题就是不存在“常用的开方数字”啊 (￣‥￣) 所以简略版的算法干脆直接求R2不开方了——反正如果R小于1的话本来就不太可能是他了

    @maplebeats
    CUDA的编译环境其实很容易装的啊 (OoO) 偶研究室里的GTX 560 Ti用的就是CentOS

    @bi119aTe5hXk
    iPhone太大了不好用而且对于偶来说那个遥控器的屏幕完全没用 (=v=o)

    @Phoenix
    其实现在写习惯了觉得C也就这么回事了，还是调教CUDA更麻烦 (;;) （主要是不像C可以随时printf调变数的值来查看函数是否正确啥的

    @conanshang
    可恶咱也好想要iPhone 5 (TAT)

    Reply

Leave a Reply

Your email address will not be published. Required fields are marked *

To create code blocks or other preformatted text, indent by four spaces:

    This will be displayed in a monospaced font. The first four 
    spaces will be stripped off, but all other whitespace
    will be preserved.
    
    Markdown is turned off in code blocks:
     [This is not a link](http://example.com)

To create not a block, but an inline code span, use backticks: 

Here is some inline `code`.

For more help see http://daringfireball.net/projects/markdown/syntax

(;;) (:D) (!!!!) (……) (^o^;) (==) (OoO) (=v=o) more »Note: Commenter is allowed to use '@User+blank' to automatically notify your reply to other commenter. e.g, if ABC is one of commenter of this post, then write '@ABC '(exclude ') will automatically send your comment to ABC. Using '@all ' to notify all previous commenters. Be sure that the value of User should exactly match with commenter's name (case sensitive).