发明名称 一种基于GPU加速的海中场景建模与实时交互绘制方法
摘要 本发明公开了一种基于GPU加速的海中场景建模与实时交互绘制方法。其步骤为:(1)对海中场景声纳数据进行预处理:水体地层分离、目标检测、去噪增强;(2)利用GPU加速实现了可交互的基于光线投射方法的海中场景三维体数据可视化建模;(3)基于GPU加速实现了海中目标声纳数据表面提取和实时绘制。本发明解决了以往可视化算法不适用海中场景声纳数据快速建模和交互可视化分析的问题。减少了主观因素影响,使得数据的初始可视化结果更加准确。本发明使用曲面细分方法对marching cube算法的目标表面提取结果进行优化处理,得到平滑的渲染结果。最终,实现了海中场景建模与实时交互绘制。
申请公布号 CN103714574B 申请公布日期 2016.05.04
申请号 CN201310703038.0 申请日期 2013.12.19
申请人 浙江大学 发明人 王章野;王丰金;廖惟博
分类号 G06T17/00(2006.01)I;G06T5/00(2006.01)I 主分类号 G06T17/00(2006.01)I
代理机构 杭州求是专利事务所有限公司 33200 代理人 张法高
主权项 一种基于GPU加速的海中场景建模与实时交互绘制方法,其特征在于包括如下步骤:1)海中场景声纳数据针对声纳数据特点采用形态学方法和滤波对原始数据进行预处理,分为基于二值化的水体地层分离、基于梯度的目标检测和基于梯度统计与滤波的去噪增强;2)利用GPU加速光线投射算法,基于数据的统计数据生成传递函数,并加入实时交互功能,实现可交互的海中场景三维体数据可视化建模;3)使用基于范围检测的marchingCube方法和GPU加速方法,实现海中目标声纳数据表面提取和实时绘制;步骤1)所述的水体地层分离为:2.1)通过LOG变换改善数据范围,拓宽有效声纳数据范围;2.2)采用高斯滤波快速去噪减弱水体中噪音的影响,并保留了水体地层分界线的原有性质;2.3)通过灰度域上的形态学操作中的开运算去除目标与噪声物体,仅保留水体与地层数据;2.4)基于水体与地层对声纳信号反应的不同,采用自动化阈值分割水体地层,阈值计算公式为:<maths num="0001" id="cmaths0001"><math><![CDATA[<mrow><msubsup><mi>&sigma;</mi><mrow><mi>w</mi><mi>i</mi><mi>t</mi><mi>h</mi><mi>i</mi><mi>n</mi></mrow><mn>2</mn></msubsup><mrow><mo>(</mo><mi>T</mi><mo>)</mo></mrow><mo>=</mo><mfrac><msub><mi>N</mi><mrow><mi>F</mi><mi>g</mi><mi>r</mi><mi>n</mi><mi>d</mi><mrow><mo>(</mo><mi>T</mi><mo>)</mo></mrow></mrow></msub><mi>N</mi></mfrac><msubsup><mi>&sigma;</mi><mrow><mi>F</mi><mi>g</mi><mi>r</mi><mi>n</mi><mi>d</mi></mrow><mn>2</mn></msubsup><mrow><mo>(</mo><mi>T</mi><mo>)</mo></mrow><mo>+</mo><mfrac><mrow><msub><mi>N</mi><mrow><mi>B</mi><mi>g</mi><mi>r</mi><mi>n</mi><mi>d</mi></mrow></msub><mrow><mo>(</mo><mi>T</mi><mo>)</mo></mrow></mrow><mi>N</mi></mfrac><msubsup><mi>&sigma;</mi><mrow><mi>g</mi><mi>r</mi><mi>n</mi><mi>d</mi></mrow><mn>2</mn></msubsup><mrow><mo>(</mo><mi>T</mi><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000900358490000011.GIF" wi="1381" he="175" /></maths><maths num="0002" id="cmaths0002"><math><![CDATA[<mfenced open = "" close = ""><mtable><mtr><mtd><mrow><msubsup><mi>&sigma;</mi><mrow><mi>b</mi><mi>e</mi><mi>t</mi><mi>w</mi><mi>e</mi><mi>e</mi><mi>n</mi></mrow><mn>2</mn></msubsup><mrow><mo>(</mo><mi>T</mi><mo>)</mo></mrow><mo>=</mo><msup><mi>&sigma;</mi><mn>2</mn></msup><mo>-</mo><msubsup><mi>&sigma;</mi><mrow><mi>w</mi><mi>i</mi><mi>t</mi><mi>h</mi><mi>i</mi><mi>n</mi></mrow><mn>2</mn></msubsup><mrow><mo>(</mo><mi>T</mi><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mo>=</mo><mrow><mo>(</mo><mrow><munder><mo>&Sigma;</mo><mrow><mi>x</mi><mo>,</mo><mi>y</mi></mrow></munder><mrow><msup><mi>f</mi><mn>2</mn></msup><mrow><mo>&lsqb;</mo><mrow><mi>x</mi><mo>,</mo><mi>y</mi></mrow><mo>&rsqb;</mo></mrow></mrow></mrow><mo>)</mo></mrow><mo>-</mo><mfrac><msub><mi>N</mi><mrow><mi>F</mi><mi>g</mi><mi>r</mi><mi>n</mi><mi>d</mi></mrow></msub><mi>N</mi></mfrac><mrow><mo>(</mo><mrow><munder><mo>&Sigma;</mo><mrow><mi>x</mi><mo>,</mo><mi>y</mi><mo>&Element;</mo><mi>F</mi><mi>g</mi><mi>r</mi><mi>n</mi><mi>d</mi></mrow></munder><mrow><msup><mi>f</mi><mn>2</mn></msup><mrow><mo>&lsqb;</mo><mrow><mi>x</mi><mo>,</mo><mi>y</mi></mrow><mo>&rsqb;</mo></mrow></mrow><mo>-</mo><msubsup><mi>&mu;</mi><mrow><mi>F</mi><mi>g</mi><mi>r</mi><mi>n</mi><mi>d</mi></mrow><mn>2</mn></msubsup></mrow><mo>)</mo></mrow><mo>-</mo><mfrac><msub><mi>N</mi><mrow><mi>B</mi><mi>g</mi><mi>r</mi><mi>n</mi><mi>d</mi></mrow></msub><mi>N</mi></mfrac><mrow><mo>(</mo><mrow><munder><mi>&Sigma;</mi><mrow><mi>x</mi><mo>,</mo><mi>y</mi><mo>&Element;</mo><msub><mi>B</mi><mrow><mi>g</mi><mi>r</mi><mi>n</mi><mi>d</mi></mrow></msub></mrow></munder><mo>-</mo><msup><mi>&mu;</mi><mn>2</mn></msup></mrow><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mo>=</mo><mo>-</mo><msup><mi>&mu;</mi><mn>2</mn></msup><mo>+</mo><mfrac><msub><mi>N</mi><mrow><mi>F</mi><mi>g</mi><mi>r</mi><mi>n</mi><mi>d</mi></mrow></msub><mi>N</mi></mfrac><msubsup><mi>&mu;</mi><mrow><mi>F</mi><mi>g</mi><mi>r</mi><mi>n</mi><mi>d</mi></mrow><mn>2</mn></msubsup><mo>+</mo><mfrac><msub><mi>N</mi><mrow><mi>B</mi><mi>g</mi><mi>r</mi><mi>n</mi><mi>d</mi></mrow></msub><mi>N</mi></mfrac><msubsup><mi>&mu;</mi><mrow><mi>B</mi><mi>g</mi><mi>r</mi><mi>n</mi><mi>d</mi></mrow><mn>2</mn></msubsup></mrow></mtd></mtr><mtr><mtd><mrow><mo>=</mo><mfrac><msub><mi>N</mi><mrow><mi>F</mi><mi>g</mi><mi>r</mi><mi>n</mi><mi>d</mi></mrow></msub><mi>N</mi></mfrac><msup><mrow><mo>(</mo><mrow><msubsup><mi>&mu;</mi><mrow><mi>F</mi><mi>g</mi><mi>r</mi><mi>n</mi><mi>d</mi></mrow><mn>2</mn></msubsup><mo>-</mo><msup><mi>&mu;</mi><mn>2</mn></msup></mrow><mo>)</mo></mrow><mn>2</mn></msup><mo>+</mo><mfrac><msub><mi>N</mi><mrow><mi>B</mi><mi>g</mi><mi>r</mi><mi>n</mi><mi>d</mi></mrow></msub><mi>N</mi></mfrac><msup><mrow><mo>(</mo><mrow><msubsup><mi>&mu;</mi><mrow><mi>B</mi><mi>g</mi><mi>r</mi><mi>n</mi><mi>d</mi></mrow><mn>2</mn></msubsup><mo>-</mo><msup><mi>&mu;</mi><mn>2</mn></msup></mrow><mo>)</mo></mrow><mn>2</mn></msup></mrow></mtd></mtr><mtr><mtd><mrow><mo>=</mo><mfrac><mrow><msub><mi>N</mi><mrow><mi>F</mi><mi>g</mi><mi>r</mi><mi>n</mi><mi>d</mi></mrow></msub><mrow><mo>(</mo><mi>T</mi><mo>)</mo></mrow><msub><mi>N</mi><mrow><mi>B</mi><mi>r</mi><mi>n</mi><mi>d</mi></mrow></msub><mrow><mo>(</mo><mi>T</mi><mo>)</mo></mrow></mrow><mi>N</mi></mfrac><msup><mrow><mo>(</mo><mrow><msub><mi>&mu;</mi><mrow><mi>F</mi><mi>g</mi><mi>r</mi><mi>n</mi><mi>d</mi></mrow></msub><mrow><mo>(</mo><mi>T</mi><mo>)</mo></mrow><mo>-</mo><msub><mi>&mu;</mi><mrow><mi>B</mi><mi>g</mi><mi>r</mi><mi>n</mi><mi>d</mi></mrow></msub><mrow><mo>(</mo><mi>T</mi><mo>)</mo></mrow></mrow><mo>)</mo></mrow><mn>2</mn></msup></mrow></mtd></mtr></mtable></mfenced>]]></math><img file="FDA0000900358490000012.GIF" wi="1581" he="661" /></maths>式中σ代表方差,T代带表数据场,σ<sub>within</sub>(T)代表数据场分割后前景方差和背景方差加权求和,σ<sub>between</sub>(T)代表前景与背景的方差,即前景与背景的差别度量,N代表数据点个数,如N<sub>Fgrnd</sub>(T)代表数据场T中前景数据点个数,N<sub>Fgrnd</sub>(T)代表数据场T中背景数据点的个数,N代表数据场整体采样点个数;f[x,y]代表数据场中位于点[x,y]处的强度值,三维可拓展为f[x,y,z],μ代表全体数据平均值,μ<sub>Fgrnd</sub>和μ<sub>Bgrnd</sub>分别带便前景与背景数据点的平均值;步骤1)所述的目标检测和去噪增强为:3.1)计算三维数据场的梯度数据场;3.2)统计梯度直方图数据场,发现数据场具有双峰特性或多峰特性;3.3)基于梯度直方图进行门限分割出强度最大的波峰,即得到目标物体;3.4)对其他波峰数据进行高斯滤波模糊处理,即得到去噪效果;所述的步骤2)为:使用GPU优化以下步骤:4.1)基于梯度数据场和去噪后的声纳数据场计算加权求和得到新数据场C<sub>trans</sub>=αC<sub>volumedata</sub>+βC<sub>gradientdata</sub>其中C<sub>volumedata</sub>代表源数据场,C<sub>gradientdata</sub>代表梯度数据场,α与β为加权系数,C<sub>trans</sub>代表合成结果;4.2)基于新数据场的梯度域统计计算得到可交互更改的传递函数;4.3)以视点为起点,向RayCasting载体立方体表面投射射线,进行等距采样并进行颜色叠加混合;c<sub>o</sub>=α<sub>s</sub>c<sub>s</sub>+(1‑α<sub>s</sub>)c<sub>d</sub>其中α<sub>s</sub>代表体素透明度,由传递函数提供,C<sub>s</sub>为体素强度,C<sub>d</sub>代表已叠加强度值,C<sub>o</sub>代表新得到的强度值,直到采样强度值的不透明度超过1,或采样深度超过数据场停止计算;步骤3)所述的基于范围检测的marchingCube方法和GPU加速方法为:5.1)将立方体的八个点进行编码2<sup>0</sup>、2<sup>1</sup>、2<sup>2</sup>、2<sup>3</sup>、2<sup>4</sup>、2<sup>5</sup>、2<sup>6</sup>、2<sup>7</sup>,将立方体所有的模式进行编码索引,进而将所有的模式放在模式纹理texture_Pattern中进行GPU端的处理与查询;5.2)利用GPU的并行处理能力以及立方体之间相互独立的特性,对数据场的所有立方体进行分类,分类检索公式为:<img file="FDA0000900358490000021.GIF" wi="1622" he="175" />其中index为当前立方体类别的索引,vertex[i]为立方体顶点i的数据场的值,[isovalue1,isovalue2]是目标的可能阈值动态范围,利用index在texture_Pattern中查找立方体所能产生的顶点个数numofvertex.然后更新相应的状态变量:cubeVertexNumber[indexofcube]=numofvertex<maths num="0003" id="cmaths0003"><math><![CDATA[<mrow><mi>c</mi><mi>u</mi><mi>b</mi><mi>e</mi><mi>O</mi><mi>c</mi><mi>c</mi><mi>u</mi><mi>p</mi><mi>i</mi><mi>e</mi><mi>d</mi><mo>&lsqb;</mo><mi>i</mi><mi>n</mi><mi>d</mi><mi>e</mi><mi>x</mi><mi>o</mi><mi>f</mi><mi>c</mi><mi>u</mi><mi>b</mi><mi>e</mi><mo>&rsqb;</mo><mo>=</mo><mfenced open = "{" close = ""><mtable><mtr><mtd><mn>0</mn></mtd><mtd><mrow><mi>n</mi><mi>u</mi><mi>m</mi><mi>b</mi><mi>o</mi><mi>f</mi><mi>v</mi><mi>e</mi><mi>r</mi><mi>t</mi><mi>e</mi><mi>x</mi><mo>&lt;</mo><mo>=</mo><mn>0</mn></mrow></mtd></mtr><mtr><mtd><mn>1</mn></mtd><mtd><mrow><mi>n</mi><mi>u</mi><mi>m</mi><mi>b</mi><mi>e</mi><mi>r</mi><mi>o</mi><mi>f</mi><mi>v</mi><mi>e</mi><mi>r</mi><mi>t</mi><mi>e</mi><mi>x</mi><mo>&gt;</mo><mn>0</mn></mrow></mtd></mtr></mtable></mfenced></mrow>]]></math><img file="FDA0000900358490000031.GIF" wi="1302" he="174" /></maths>其中indexofcube为当前立方体在GPU上的编码;5.3)利用Prefix‑Sum算法对cubeOccupid数组进行处理,在cubeOccupidScan数组中保存当前包含顶点的立方体在所有包含顶点的立方体的排序;5.4)利用Prefix‑Sum处理的结果对立方体进行压缩处理,将不含有顶点的立方体剔除出去,压缩方法为:compactCube[cubeOccupidScan[i]]=i if(cubeOccupid[i]&gt;0)其中i为当前立方体的索引;处理结束后,compactCube数组中保存的是所有包含三角面片的立方体的索引;5.5)利用cube的256中模式分别对每个含有三角面片的立方体进行处理,产生最终的顶点、法向、以及三角面片,计算公式为:<maths num="0004" id="cmaths0004"><math><![CDATA[<mrow><mi>v</mi><mi>e</mi><mi>r</mi><mi>t</mi><mi>e</mi><mi>x</mi><mo>=</mo><msub><mi>v</mi><mn>1</mn></msub><mo>+</mo><mrow><mo>(</mo><mi>i</mi><mi>s</mi><mi>o</mi><mi>v</mi><mi>a</mi><mi>l</mi><mi>u</mi><mi>e</mi><mo>-</mo><mi>v</mi><mi>a</mi><mi>l</mi><mi>u</mi><mi>e</mi><mn>1</mn><mo>)</mo></mrow><mfrac><mrow><msub><mi>v</mi><mn>2</mn></msub><mo>-</mo><msub><mi>v</mi><mn>1</mn></msub></mrow><mrow><msub><mi>value</mi><mn>2</mn></msub><mo>-</mo><msub><mi>value</mi><mn>1</mn></msub></mrow></mfrac></mrow>]]></math><img file="FDA0000900358490000032.GIF" wi="1204" he="159" /></maths>其中vertex为产生的顶点坐标,v<sub>1</sub>、v<sub>2</sub>为两个端点的坐标,value<sub>1</sub>、value<sub>2</sub>为两个端点的体数据值,Isovalue为体数据的阈值范围的平均值:isovalue=(isovalue1+isovalue2)/2normal为三角面片的法向:normal=(v2‑v1)×(v3‑v2)5.6)利用openGL显示最终的结果,加入了Tessellation Shader,设置曲面细分的等级,得到平滑的最终显示结果;5.7)将marching cube获得的结果以obj格式进行保存。
地址 310027 浙江省杭州市西湖区浙大路38号