设为首页 加入收藏

1. 引言

边缘检测是图像处理和计算机视觉中的基本方法，构成了图像特征提取的重要研究领域。边缘检测的目的是标识数字图像中亮度变化明显的点。图像属性中的显著变化通常反映了属性的重要事件和变化。例如：深度上的不连续、表面方向上的不连续、物质属性的变化和场景照明的变化等。因此，物质属性的变化必然在探地雷达GPR (Ground Penetrating Radar)成像中有所反映^[1]。本文探讨应用近年来新发展起来的GPU-CUDA^[2]技术实现Prewitt边缘检测算法，提出将实测探地雷达数据转换成灰度图像数据后，再实施基于GPU的Prewitt算法处理的方法思路，对实测探地雷达数据进行处理^[3]和基于GPU的快速成像，以求取得好的应用效果。试验结果表明该算法不仅运行效率高，较好的解决了探地雷达一体化处理解释中的图形冻结问题，而且在突出有效异常，提高目标体的识别能力方面取得实效。

2. Prewitt算法原理

多年以来，人们提出了许多以微分法^[4]为基础的边缘检测算法，并以数据块与模板卷积的手段提取图像边缘。常用边缘检测算法包括Roberts Cross算法、Prewitt算法、Sobel算法、Canny算法和Laplace算法。Prewitt算法作为一阶微分算子的边缘检测算法^[5]，它利用像素点左右、上下相邻点的灰度差在边缘处达到极值的重要特征来检测边缘和去掉部分伪边缘，对噪声具有平滑作用。其原理是在图像空间中，将左右、上下相邻点构成的3 × 3数据块(图1)与水平、垂直两个方向模板进行卷积运算来完成。其中：水平模板(图2)用于检测水平边缘，垂直模板(图3)用于检测垂直边缘。

对于像素点p(i,j)，Prewitt算法定义如下：

水平方向梯度：

图1. P(i,j)像素点及其3 × 3像素邻域

图2. 水平模板

图3. 垂直模板

垂直方向梯度：

梯度：

其中：x为像素点p(i,j)邻域的3 × 3数据块与水平模板的卷积；y为像素点p(i,j)邻域的3 × 3数据块与垂直模板的卷积；计算出梯度g(i,j)后，把g(i,j)值大于255的点标记为边界点，其像素值设定为255(即白点)，否则设定像素值为g(i,j)的值。Prewitt算法的串行实现简单，但是，由于需要做大量的卷积运算，当图像分辨率较高时，计算量很大，需要较长的运算时间，难以满足实际工作要求。

3. Prewitt算法的GPU-CUDA实现

通过Prewitt算法分析可知：计算过程中不需要利用已经完成的像素处理结果参与，符合GPU-CUDA (Compute Unified Device Architecture)数据并行处理要求。因此，可以将Prewitt算法中的所有单线程循环改为GPU-CUDA多线程并行处理，从而提高算法的效率。根据CUDA程序设计思想(即CPU进行数据准备、初始化GPU设备以及执行串行代码，GPU进行并行计算并向CPU返回计算结果)，按如下算法过程完成Prewitt算法在GPU-CUDA的实现：

1) 由主机端程序完成数据的初始化；

2) 申请两块与图像数据相等容量的设备内存空间，其中一块用于存放由主机传送过来的图像数据，另一块用于存放设备中内核并行程序进行Prewitt算法处理后的结果；

3) 将主机端图像数据拷贝回至设备端；

4) 调用设备内核并行程序；

5) 将处理后的图像数据拷贝回主机端。

Prewitt并行算法的关键是并行处理网格的构造以及内核并行程序实现，构造并行处理网格源代码如下^[6]：

int lineByte = screenWidth*sizeof(float);//计算图像数据大小

cudaMalloc((void **)&d_DataIn,lineByte*screen High);//申请设备内存空间

cudaMalloc((void **)&d_DataOut,lineByte*screen High);//申请设备内存空间

cudaMemcpy(d_DataIn,h_DataIn, lineByte*screen High,cudaMemcpyHostToDevice);//图像数据拷贝至设备端

dim3 dimBlock1(256,1);//构造并行处理网格

dim3 dimGrid1(((screenWidth-2)*(screenHigh-2) + 255)/256,1);//构造并行处理网格

prewittProcess<<<dimGrid1,dimBlock1>>>(d_Data In,d_DataOut,screenWidth,screenHigh);//调用内核并行处理函数

cudaMemcpy(h_DataOut,d_DataOut,lineByte*screen High,cudaMemcpyDeviceToHost);//处理结果拷贝回主机端

内核并行程序实现源代码如下：

__global__ void prewittProcess(float *d_DataIn, float*d_DataOut,int screenWidth,int screenHigh)

{unsigned int id = 256*gridDim.x + threadIdx.x;//计算并行处理网格的索引号id if(id<(screenWidth-2)*(screenHigh-2))//控制图像处理范围

{int row=id/(screenWidth-2) + 1;//图像矩阵行编号row int col=id%(screenWidth-2) + 1; //图像矩阵列编号col int x,y,result;

x=*(d_DataIn+(row-1)*screenWidth+(col-1))+*(d_DataIn+(row-1)*screenWidth+(col))+*(d_DataIn+(row-1)*screenWidth + (col+1))-*(d_DataIn+(row+1)*screen Width+(col-1))-*(d_DataIn+(row+1)*screenWidth+(col))-*(d_DataIn+(row+1)*screenWidth+(col+1));//计算水平方向梯度

y=*(d_DataIn+(row-1)*screenWidth+(col-1))+*(d_DataIn+row*screenWidth+(col-1))+*(d_DataIn + (row + 1)*screenWidth + (col-1))-*(d_DataIn + (row-1)*screen Width+(col+1))-*(d_DataIn+row*screenWidth+(col+1))-*(d_DataIn + (row + 1)*screenWidth+(col+1));//计算垂直方向梯度

result=x*x+y*y;//计算梯度值g(i,j)的平方值

*(d_DataOut+row*screenWidth+col)=result;//保存Prewitt并行处理结果}}

实验在Acer ASPIRE 4736G(含GEFORCE G105M CUDA 512M图形卡)手提电脑上完成，试验计算结果如表1所示，获得了较好的加速比，满足了对探地雷达数据进行实时交互处理解释的实际需要，圆满地解决了操作界面冻结问题。

4. Prewitt算法在探地雷达中的成功应用

4.1. 实测雷达数据转化为灰度图像数据

应用Prewitt边缘检测技术处理探地雷达数据的前提是需要把实测雷达数据转化为灰度图像数据。实测雷达数据是按道存储，必须在转化为灰度图像数据的同时进行行列变换。具体实现方法为^[7]：遍历所有实测雷达数据求取最大、最小值，将实测雷达按数据最大值对应255，最小值对应0的准则进行线性变换求取灰度图像数据。与此同时，完成行列变换后成为真正的灰度图像数据。主要代码如下：

for(int i=0;i<screenHigh;i++)//遍历所有实测雷达数据求取最大、最小值

for(int j=0;j<screenWidth;j++)

{if(min>*(h_DataIn+i*screenWidth+j)) min=* (h_DataIn+i*screenWidth+j);

if(max<*(h_DataIn+i*screenWidth+j)) max=* (h_DataIn+i*screenWidth+j);}

for(int i=0;i<screenHigh;i++)//每点数据进行行列变换，同时转换为灰度图像数据。

for(int j=0;j<screenWidth;j++)

{*(h_DataIn+i*screenWidth+j)=((255**(h_DataIn+ i*screenWidth+j)/(max-min)-255*min/(max-min)));}

表1. 不同分辨率图像CUP与GPU处理时间对比

以上预处理工作完成后，调用基于GPU-CUDA实现的Prewitt算法对灰度图像数据进行边缘检测。边缘检测结果再按如下实现的GPU快速彩色成像方法成像。

4.2. Prewitt处理结果的GPU快速彩色成像

探地雷达GPU快速彩色成像算法的关键是在色标数据基础上，如何实现色标数据的转换及彩色成像的GPU多线程处理^[8]？文中通过GPU多线程并行并读写OpenGL成图缓冲区方法直接成像，极大地提高了探地雷达成像速度，解决了当前探地雷达彩色成像速度慢造成程序操作界面冻结等问题。探地雷达GPU快速彩色成像算法流程图如图4所示^[9]。

根据图4所示流程图，设计快速彩色成像算法内核并行程序源代码如下：

__global__void prepareMap(float4 *positions,float *seismic, float *tmp_color, int screenWidth, int screenHigh, float min, float max, int colorTableLength, float TranspSliderValue)//每个线程负责生成一个像素点数据

{unsigned int i=blockIdx.x*blockDim.x + thread

图4. GPU快速彩色成像算法流程图

Idx.x;//计算并行处理网格的索引号i unsigned int j=blockIdx.y*blockDim.y+threadIdx.y; //计算并行处理网格的索引号j if(i<screenWidth && j<screenHigh)//成像数据范围：screenWidth*screenHigh

{int index=(int)(((colorTableLength-1)**(seismic +i*screenHigh+j)/(max-min)-(colorTableLength-1)*min/(max-min)));//根据定义的色标数据，计算成像数据

if(index<0)index=0; if(index>colorTableLength-1) index=colorTableLength-1; *(positions+i*screenHigh+j) =make_float4(*(tmp_color+index*3 + 0),*(tmp_color + index*3 + 1), *(tmp_color + index*3 + 2),TranspSlider Value);//生成彩色成像数据}}

4.3. 实验结果与分析

图5是实测探地雷达数据采用GPU快速彩色成像算法的直接成像，并包括探地雷达软件系统的总控界面；图6是实测探地雷达数据经过基于GPU-CUDA的Prewitt边缘检测处理后再采用GPU快速彩色成像算法的成像；由于Prewitt算法和彩色成像均采用GPU多线程并行处理，两者组合处理能达到近10倍的加速比，克服了常规CPU单进程处理速度过慢引起的窗口界面冻结问题。从图6看出，目标体周围的电磁波绕射干扰^[10]受到明显压制，目标体的边界更加清晰。位于A、B、C三处的地下管线埋设情况及管线边界更加准确，其中，A、B两处地下管是距地面埋藏深度只有25 CM的塑料管，C处是人工砖砌排水管，埋藏深度1.1 M(注：为了便于观察目标异常，已对上述成像在垂直方向进行了1 M的向下平移)。由此可见，应用Prewitt算法处理实测探地雷达数据对突出有效异常，提高目标体的识别能力效果明显，是一种提高探地雷达分辨率的有效方法。

5. 结语

在分析Prewitt边缘检测算法的基础上，首先设计并实现了基于GPU-CUDA的Prewitt快速算法。为了能够应用Prewitt快速算法于探地雷达的数据处理，文中提出了先将实测探地雷达数据转化为灰度图像数据，然后再调用Prewitt快速算法处理，最后应用GPU快速彩色成像算法成像的工作流程。实验结果表明：

图5. 野外实测探地雷达成像

图6. Prewitt算法处理后的探地雷达成像

该工作流程及实现方法既有效解决了常规CPU算法的界面冻结问题，又在突出有效异常，提高目标体的识别能力方面取得了良好的应用效果。

参考文献 (References)

[1]       粟毅. 探地雷达理论与应用[M]. 北京: 科学出版社, 2006.

[2]       NVIDIA. CUDA programming guide, 2012. http://www.nvidia.com/object/cuda_home.html

[3]       杨峰, 彭苏萍. 地质雷达探测原理与方法研究[M]. 北京: 科学出版社, 2010.

[4]       李荣华, 刘播. 微分方程数值解法(第四版)[M]. 北京: 高等教育出版社, 2010.

[5]       曾胜田, 刘羽等. 基于CUDA的Prewitt算子并行实现[J]. 微计算机应用, 2011, 11: 71-75.

[6]       蒋长锦. 科学计算与C程序集[M]. 北京: 中国水利水电出版社, 2010.

[7]       J. E. Lucius, M. H. Powers. GPR data processing computer software for the PC. Reston: United States Geological Survey (USGS), 2002: 151-158.

[8]       T. Y. Peng, J. X. Cao. The implementation and fast visualization application of Linux-based QT-GPU parallel architecture. American Journal of Engineering and Technology Research, 2011, 11 (12): 2063-2068.

[9]    彭土有, 叶亚平等. 意大利RIS-2K探地雷达数据解密及其快速可视化方法[J]. 计算机时代, 2012, 7: 17-20.

[10]    孟陆波, 李天斌等. 地质雷达超前预测不良地质体图像的智能识别[J]. 煤田地质与勘探, 2009, 2: 86-89.

NOTES

^*基金项目：基于Linux地质雷达一体化处理解释技术及应用基金项目：2011年省部产学研结合项目(2011B090400120)。