使用等值面绘制和体绘制2种可视化方法对电磁环境进行融合绘制是本文研究的主要内容。本文在分析2种电磁环境可视化方法特点的基础上，对等值面绘制进行了适当修改，不仅便于其融合到体绘制结果中，还能够增强显示等值面的深度信息，有利于用户直接观察电磁环境，从而理解其中的信息。
等值面绘制和体绘制方法的计算量均较大，通常需要基于硬件进行加速实现[4-6]。因此，为保证融合后的算法执行效率，同样需要进行基于硬件的加速实现。近年来，统一计算设备架构(Compute Unified Device Architecture，CUDA)由于其对通用计算的支持，在并行算法的实现方面应用广泛[7-8]。因此，研究了融合绘制算法在CUDA架构下的并行实现，以获得高执行效率。
1 电磁环境体可视化 现有的电磁环境可视化研究中，多使用模拟数据，如文献[3-4]中的数据生成方法，以三维数据场表示电磁环境。因此，可视化电磁环境即是对电磁环境数据场进行可视化，将此三维数据场绘制为二维图像以展示给用户。绘制三维数据场效果最好的是光线投射体绘制算法[9]，本文基于此算法对电磁环境数据场进行绘制。
但是，光线投射体绘制算法将三维数据场映射为二维图像后，在图像中丢失了相关的深度信息[10]，虽然能够绘制出电磁环境数据场的整体效果，但是缺少对其边界等形状信息的展示。因此，在光线投射体绘制算法的绘制结果中添加更多显示信息能够提高电磁环境的绘制效果。
2 融合等值面绘制 等值面绘制是电磁环境的另一种三维可视化方法，其主要显示电磁环境数据场中某些等值面的分布情况。等值面绘制与体绘制方法中，对电磁环境的表现侧重点不同，本文将二者融合，以显示更多信息，从而提高电磁环境的绘制效果。
2.1 网格化等值面 在等值面绘制中，通常以不同颜色表示不同的等值面[11]，而体绘制中同样通过不同颜色来表现不同的数据信息[12]。因此，融合后的算法中应该注意二者在颜色上的区分。本文对传统的等值面绘制进行了适当修改，不仅利于等值面绘制的融合，还能够增强其信息表现能力。
将等值面进行网格化，即使用覆盖等值面的网格代替等值面，称此网格为等值面网格。如图 1所示，均匀网格平行投影到空间曲面后，网格将附着在曲面表面而变形，从侧面观察变形后的网格，能够得出曲面起伏变化的相关信息。因此，投影后的网格能够反映出空间曲面的信息，可以使用网格代替等值面，即将等值面替换为如图 1所示的投影到此等值面后的网格。
 |
| 图 1 网格投影到空间曲面 Fig. 1 Projecting grid to a curved space surface |
| 图选项 |
2.2 等值面网格提取 提取等值面[13]是绘制等值面的前提条件，因此，网格化等值面后，首先需要提取等值面网格。如图 1所示，均匀网格投影到空间表面后，网格将沿此表面展开，即投影后的网格与空间表面重合。因此，使用均匀网格投影的方法提取等值面网格。
将均匀网格平行投影到电磁环境数据场中，如图 2所示，其中，网格的所有格孔处于同一平面且均为相等的正方形。投影前需要将连续网格离散为均匀点阵，然后将均匀点阵投影到指定等值面。
 |
| 图 2 均匀网格投影到数据场 Fig. 2 Projecting inerratic grid to data volume |
| 图选项 |
网格均匀离散过程中，需要保证网格中任意2条网格线的交叉点处均有离散后的点；均匀离散即保证任意2个相邻(同一条网格线上)离散点间距离相等。对于一条网格线，其将被离散成多少个点由式(1) 确定：
 | (1) |
式中：I为生成图像的分辨率；D为体数据的大小，下标x、y、z表示对应的维度；L为网格线长度(网格为正方形长宽相等)；D和L的度量单位相同。
离散后的所有网格点沿同一方向向数据场进行投影，网格点最终投影到其投影方向上的第1个值为ε的数据点，ε为指定的等值面值。如此，所有网格点的投影位置均是值为ε的数据点，即网格覆盖在值为ε的等值面上，完成了对体数据中等值面网格的提取。
2.3 与光线投射的融合 网格化等值面后，需要将网格投影融合到光线投射操作中，已完成融合绘制。离散网格投影与光线投射均是对数据场的操作，融合后二者共同对数据场进行操作的关系如图 3所示。
 |
| 图 3 离散网格投影与光线投射的关系 Fig. 3 Relationship between discrete grid projecting and ray-casting |
| 图选项 |
融合绘制的最终结果必然是生成一幅图像，而图 3中的光线投射已经定义了一幅图像，可以将投影后的网格绘制到光线投射生成的图像中，完成融合绘制。做一条直线连接网格投影点与光线投射的视点，则此直线穿过图像所在平面。当直线与平面的交点在图像内部时，选择距离此交点最近的像素点作为网格投影点对应的像素点，使用标准颜色混合算法将网格投影点的颜色C1(R1, G1, B1, A1)与已生成的像素点颜色值C0(R0, G0, B0, A0)合成为颜色C(R, G, B, A)作为此像素点的新颜色值。
颜色合成时，网格的颜色被融合到光线投射生成的图像中。当需要同时绘制多层等值面网格时，使用深度优先原则，从离视点最远的等值面网格开始，依次将每层网格融合到图像中。融合多层等值面网格为半透明融合，采用深度优先原则能够避免不同深度网格之间的相互覆盖。
3 算法描述 如图 3所示，网格投影和光线投射均对数据场进行处理，且影响生成图像的像素点值。两者能够独立生成相应的像素点颜色值后再进行颜色融合，因此，两者的先后执行顺序无需严格限定。本文由于是将多层投影后的网格融合到体绘制结果中，指定执行顺序为：首先进行光线投射生成一幅图像，然后按照等值面网格的深度依次执行对应的网格投影，每层网格的投影结果均绘制在已生成的图像当中。
3.1 并行线程 由于光线投射独立执行，为提高其执行效率, 在CUDA下并行实现。使用文献[14]中的方法实现电磁环境数据场的并行光线投射绘制算法，本文重点关注网格投影的并行执行。由图 2可知，离散后的网格投影中，每个网格点的操作是相同的，易于并行实现。CUDA下，每个离散网格点投影的主要步骤如下：
步骤1??计算网格点V沿投影方向穿越数据场的起点S和穿出点E。
步骤2??网格点V沿投影方向从点S开始以步长α对数据场进行采样。
步骤3??判断当前采样点P的值是否为要提取的等值面值ε，是则转到步骤5，否则继续步骤4。
步骤4??沿投影方向前进一个步长α对数据场进行采样，若采样点穿出数据场，转到步骤6，否则，转到步骤3。
步骤5??作连接采样点P与光线投射中视点O的直线l，若l与图像相交，获取离交点最近的像素点M，将网格点V对应的颜色值融合到像素点M中。
步骤6??结束。
上述步骤中，采样步长α与光线投射时的采样步长相同，设置为数据场中同一维度上相邻两点间的距离。步骤3中，判断采样点值是否为等值面值使用如下方法：当采样点值第1次不小于等值面值时，即认为二者相等。如此进行判断是因为电磁环境数据场中，数据值由连续数据场采样得到，当网格从数据场上方投影时，网格点遇到的数据点值将递增。采用上述方法判断能够避免由于投影时采样步长的影响使网格点错过等值面的情况。后续的实验结果中能够说明此判断方法是可行的。
在实际地理环境中，观察视点通常不会位于电磁环境数据场下方，而电磁辐射源处于数据场内部下方。因此，离观察视点较近(数据场外侧)的数据一般数据值较小，投影后网格相对视点的距离由网格对应的等值面值确定，值越大则距离越大。因此，根据等值面值的大小确定对应网格的投影顺序。
整个算法在CUDA下执行时的关键代码为：主机端对光线投射和网格投影的控制流程；每层网格投影时在设备端的执行流程。2个不同流程对应的伪代码如下所示。
1) 主机端：控制流程
1. ??g_CastRay < < < raySize, blockSize>>>(Image);
2. ??for i=1 to the number of isosurface grids
3.????g_ProjectGrid < < < gridSize, blockSize>>>(i, Image);
4. ??end for
2) 设备端：网格投影——g_ProjectGrid()函数
5. ??X=blockIdx.x*blockDim.x+threadIdx.x; //线程的X坐标
6. ??Y=blockIdx.y*blockDim.y+threadIdx.y; //线程的Y坐标
7. ??gridX=getGridX(X, Y); //获取网格点X坐标
8. ??gridY=getGridY(X, Y); //获取网格点Y坐标
9. ??if (gridX, gridY) is in the Grid then
10. ???isoValue=getGridValue(i); //获取等值面值
11. ????gridColor=getGridColor(i); //获取网格颜色
12. ????(S, E, d)=getProjectParameters(gridX, gridY);
13. ??for P=S to E
14.????v=getVolumeDataValue(P); //获取数据值
15.????if v is not less than isoValue then
16. ????????gridPosition=P;
17. ????????mark=true;
18. ????????break;
19. ????end if
20.????P=P+d*α; //沿投影方向向前采样
21. ??end for
22. ??if mark is true then
23.????addGridColor2Image(gridPosition, gridColor);
24. ??end if
25.end if
上述流程中，第(1) 行中函数g_CastRay()完成光线投射操作并生成图像Image。第1行和第3行中的参数raySize、blockSize、gridSize为开启线程个数的参数，blockSize固定为16×16，raySize和gridSize分别根据光线和网格点的数量进行设置。网格投影流程中，第12行中的函数getProjectParameters()根据当前投影网格点的位置(gridX, gridY)获取其投影穿越数据场时的参数：起点S、终点E和前进方向d。第23行中的函数addGridColor2Image()根据网格点在数据场中的投影位置gridPosition计算其对应的像素点，并进行颜色融合。
3.2 交互控制 并行执行算法的效率高，能够及时响应用户的交互更新生成图像，从而使用户能够及时观察到调整的结果，利于算法的使用。因此，算法支持用户的交互控制操作，使绘制结果能够被用户控制，满足用户的观察需求。算法提供给用户交互控制的参数主要有：光线投射中的颜色设置；网格投影中各层网格的大小、格孔大小、颜色设置、投影位置等。同时，用户对数据场的观察视角可以随时进行调整，获得对应的融合绘制结果，便于全面观察数据场。
4 实验验证与分析 本文模拟计算电磁环境数据场的方法与文献[4]中相同，由于光线投射体绘制算法独立于网格投影，使用直接体绘制(DVR)算法进行光线投射，以验证融合绘制的效果。实验中，重点关注使用网格投影将等值面绘制融合到光线投射体绘制结果中的效果，分别从网格的层数、交互控制绘制参数等方面进行了展示。统计说明了算法在CUDA下的并行实现效率。
4.1 绘制结果与分析 实验中，投影两层等值面网格，并生成分辨率为800×600的图像。借鉴文献[15]中对使用双层网格时的结论，设置投影后离视点较近的网格格孔大于较远的网格格孔。同时，投影的两层网格之间网格线互相交叉：平移一个网格到另一网格平面时，网格线不平行，夹角取45°。为更好地表现电磁环境数据场，对投影的两层网格进行如下设定：在数据场最外侧固定投影一层等值面网格，以表现整个数据场在空间的分别范围；另一层投影到数据场内部，且其投影深度可随用户交互进行变化，以观察数据场内部不同的等值面分布情况。
图 4中对比显示了融合绘制等值面网格前后的结果。图 4(a)为只使用体绘制生成的图像；图 4(b)为在图 4(a)基础上融合绘制两层等值面网格的结果；图 4(c)显示了只融合绘制外层网格的结果；图 4(d)则显示了只融合绘制内层网格的结果。
 |
| 图 4 融合等值面绘制的结果对比 Fig. 4 Comparing results of fusion isosurface rendering |
| 图选项 |
通过图 4中不同绘制结果的对比可以看出，融合绘制等值面后的结果立体感增强了：投影后的绿色网格能够反映出数据在空间中的覆盖范围，红色网格能够反映出数据场内部指定等值面的分布情况。以上分析说明，将等值面进行网格化并融合到体绘制算法后，生成的绘制结果增强了对原数据场中的信息表现能力，便于用户感知和理解更多电磁环境数据场中的信息。
算法在融合绘制的过程中支持用户对绘制参数的交互调整，图 5显示了在用户交互调整情况下生成的不同绘制结果。图 5 (a)为调整内层投影的等值面网格值后对应的不同绘制结果，可以观察到电磁环境数据场中各个等值面的分布情况；图 5 (b)为不同观察视角下的融合绘制结果，用户可以根据需要调整对电磁环境数据场的观察视角，便于使用。
 |
| 图 5 不同参数下的融合绘制结果 Fig. 5 Results of fusion rendering with different parameters |
| 图选项 |
图 5中的红色网格格孔比图 4中的要小，这是由于进行了交互调整的结果。除了图 5中使用的可交互调整的参数外，算法还允许对各层网格的颜色及其不透明度等进行调整。上述参数的调整均能够控制融合绘制的网格效果，使绘制结果尽量满足用户的期望。同时，网格投影与光线投射相互独立，改变光线投射的算法时不影响网格投影的操作，仍可对网格参数进行调整以获得较好的融合绘制结果。
将上述融合绘制电磁环境的结果与文献[3-4]进行对比可知，融合绘制结果在体绘制结果的基础上能够表现出更多的数据场信息，更利于用户对电磁环境的观察和感知。
4.2 并行执行验证 算法在CUDA架构下编程实现，验证算法能够并行执行从而获得较高的执行效率。实验中主要使用的软件为NVIDIA GPU Computing SDK 4.2和Microsoft Visual Studio 2005，硬件支持环境为Intel Core i5-2400 CPU、4GB内存、NVIDIA GeForce GTX470显卡。
图 6为对融合绘制算法并行执行时间的统计结果，主要统计了视角变化情况下的算法执行时间。图 6(a)为实验中不同视角的位置示意图；图 6(b)为统计的200次算法执行时间。实验中，网格位于数据场上方并向下投影，所有观察视点均分布于数据场外部的一个半球面上，使用鼠标移动视点并统计了200次算法的执行时间。图 6(b)中有3条曲线，分别对应不同情况下的时间统计结果：下面曲线是只进行体绘制的情况，对应图 4(a)；中间曲线是只投影内层网格的情况，对应图 4(d)；上面曲线是投影两层网格的情况，对应图 4(b)。
 |
| 图 6 不同观察视角下算法的执行时间统计 Fig. 6 Algorithm execution time statistics under different viewpoints |
| 图选项 |
从图 6(b)中上面那条曲线能够得出结论为：融合绘制算法的并行执行时间在6ms以内，完全满足实时绘制电磁环境数据场的需要。同时，由于光线投射与网格投影两者相互独立，用图 6(b)中上面曲线减去下面曲线可得出投影两层网格的执行时间小于4ms。同理可得，投影一层网格的执行时间小于2.5ms(图 6(b)的中间曲线减去下面曲线)。可见，投影网格并行执行的效率很高，在投影两层网格(图 4(b)和图 5中所示)且保持25帧/s的绘制速度情况下，能够为光线投射体绘制算法提供超过36ms的执行时间。因此，使用本文算法进行融合绘制不仅效果较好，还具有很高的执行效率，为支持更多的光线投射体绘制算法奠定了基础。
将实验中使用的光线投射体绘制算法在CPU中串行实现后，绘制一幅图像的时间在1s以上，加入网格投影后图像生成时间将更久，不利于算法使用。同时与图 6中的执行时间进行对比，充分说明了融合绘制算法在CUDA下并行实现的效率优势。
5 结论 1) 通过网格化等值面，将电磁环境数据场的等值面绘制与体绘制进行了融合，使融合绘制结果兼具2种绘制方法的效果。
2) 设计了算法的并行实现并进行了验证，算法具有较高的执行效率，可支持实时交互控制。
本文提出的算法只是对2种绘制方法进行融合的初步研究，讨论情况较简单，在将来的研究中将讨论更加复杂和精细的融合绘制情况。
参考文献
| [1] | 吴迎年, 张霖, 张利芳, 等. 电磁环境仿真与可视化研究综述[J].系统仿真学报, 2009, 21(20): 6332–6338. WU Y N, ZHANG L, ZHANG L F, et al. Survey on electromagnetic environment simulation and visualization[J].Journal of System Simulation, 2009, 21(20): 6332–6338.(in Chinese) |
| [2] | CHEN P, WU L D. 3D representation of radar coverage in complicated environment[J].Simulation Modeling Practice and Theory, 2008, 16(9): 1190–1199.DOI:10.1016/j.simpat.2008.06.005 |
| [3] | 杨超, 徐江斌, 赵健, 等. 基于多层等值面的电磁环境三维可视化研究[J].系统工程与电子技术, 2009, 31(11): 2767–2772. YANG C, XU J B, ZHAO J, et al. Research on 3D visualization of electromagnetic environment based on multi-isosurface[J].Systems Engineering and Electronics, 2009, 31(11): 2767–2772.(in Chinese) |
| [4] | 杨超, 徐江斌, 吴玲达. 硬件加速的虚拟电磁环境体可视化[J].北京邮电大学学报, 2011, 34(1): 55–59. YANG C, XU J B, WU L D. Hardware accelerated volume visualization in virtual electromagnetic environment[J].Journal of Beijing University of Posts and Telecommunications, 2011, 34(1): 55–59.(in Chinese) |
| [5] | 于荣欢, 邓宝松, 吴玲达, 等. 三维标量场并行等值面提取与绘制技术[J].计算机辅助设计与图形学学报, 2012, 24(2): 244–251. YU R H, DENG B S, WU L D, et al. Parallel isosurface extracting and rendering of 3D scalar fields[J].Journal of Computer-Aided Design & Computer Graphics, 2012, 24(2): 244–251.(in Chinese) |
| [6] | BEYER J, HADWIGER M, PFISTER H.A survey of GPU-based large-scale volume visualization[C]//Proceedings of the Eurographics Conference on Visualization 2014, Swansea, 2014:105-123. |
| [7] | ROSEN P. A visual approach to investigating shared and global memory behavior of CUDA kernels[J].Computer Graphics Forum, 2013, 32(3): 161–170. |
| [8] | ZHANG Y, DONG Z, MA K L. Real-time volume rendering in dynamic lighting environments using precomputed photon mapping[J].IEEE Transactions on Visualization and Computer Graphics, 2013, 19(8): 1317–1330.DOI:10.1109/TVCG.2013.17 |
| [9] | 马千里, 李思昆, 白晓征, 等. CFD非结构化网格格心格式数据高质量体绘制方法[J].计算机学报, 2011, 34(3): 508–516. MA Q L, LI S K, BAI X Z, et al. High-quality volume rendering of unstructured-grid cell-centered data in CFD[J].Chinese Journal of Computers, 2011, 34(3): 508–516.(in Chinese) |
| [10] | BRAMON R, RUIZ M, BARDERA A, et al. An information-theoretic observation channel for volume visualization[J].Computer Graphics Forum, 2013, 32(3): 411–420. |
| [11] | MARTIN T, COHEN E, KIRBY R M. Direct isosurface visualization of hex-based high-order geometry and attribute representations[J].IEEE Transactions on Visualization and Computer Graphics, 2012, 18(5): 753–766.DOI:10.1109/TVCG.2011.103 |
| [12] | WANG L, GIESEN J, MCDONNELL K T, et al. Color design for illustrative visualization[J].IEEE Transactions on Visualization and Computer Graphics, 2008, 14(6): 1739–1746.DOI:10.1109/TVCG.2008.118 |
| [13] | SCHOLZ M, BENDER J, DACHSBACHER C. Real-time isosurface extraction with view-dependent level of detail and applications[J].Computer Graphics Forum, 2015, 34(1): 103–115.DOI:10.1111/cgf.2015.34.issue-1 |
| [14] | 冯晓萌, 吴玲达, 董士伟. CUDA加速的动态电磁环境数据场实时绘制[J].系统仿真学报, 2014, 26(9): 2044–2049. FENG X M, WU L D, DONG S W. CUDA accelerated real-time rendering for dynamic electromagnetic environment volume data[J].Journal of System Simulation, 2014, 26(9): 2044–2049.(in Chinese) |
| [15] | BAIR A, HOUSE D. A grid with a view:Optimal texturing for perception of layered surface shape[J].IEEE Transactions on Visualization and Computer Graphics, 2007, 13(6): 1656–1663.DOI:10.1109/TVCG.2007.70559 |
