2. 福建师范大学 医学光电科学与技术教育部重点实验室, 福州 350007;
3. 福建师范大学 福建省光子技术重点实验室, 福州 350007
2. Key Laboratory of OptoElectronic Science and Technology for Medicine of Ministry of Education, Fujian Normal University, Fuzhou Fujian 350007, China;
3. Fujian Provincial Key Laboratory of Photonics Technology, Fujian Normal University, Fuzhou Fujian 350007, China
建筑物是城市下垫面主要的地物类型,是城市大比例尺基础地理底图中必须重点表现的专题要素。通过自动识别与提取高分辨率影像上的建筑物,以此提高城市基础地理数据库建设与更新的效率、降低人工作业的工作量与劳动强度具有重要意义。
目前,各类高分辨率卫星遥感影像以及航空遥感影像得到普及,成为建筑物提取的主要数据源,以下对相关的建筑物提取方法进行综述。
1) 基于分类的方法。分类的主要目的是先进行建筑物或建筑物区域的粗提取,再利用先验知识或建筑物特征进行筛选、提高检测性能[1-3]。通过分类可以排除一部分非建筑物区域的干扰,提高方法的效率。但分类器的选择和设计是关键,通常需要根据具体场景对分类器的参数进行调校,是一个较为繁琐且复杂的过程。
2) 基于分割的提取方法。为了获取更为准确的建筑物对象,通常先经过分割并计算分割区域内的梯度、纹理、方向和形状等特征信息,获得候选建筑物,再根据空间格局和结构语义信息,如建筑物和道路的空间关系、建筑物之间的方向一致性等提取建筑物对象,使得准确轮廓信息的获取成为可能[4-6]。但缺点是准确性过分依赖于分割算法的性能,过分割或欠分割都会造成建筑物的误检或漏检。
3) 基于阴影的方法。阴影包含了建筑物结构的大量信息,而无需另外描述其模型,因此被认为是用于建筑物提取的一种重要特征。一般的做法是根据阴影区域的饱和度和强度遵循某种特征,通过阈值方案检测阴影区域[7-8]。理论上,准确的阴影可以用于任意形状建筑物的自动提取,然而,在目前的大多数研究中,阴影没得到充分的提取,一般只是提取建筑物周围的暗区域,另外,空间关系如何确定以及阴影的自动准确提取仍有待进一步研究。
4) 基于边缘的提取方法。该类方法先利用边缘检测算子提取边缘,根据空间关系对边缘进行分组,搜索符合建筑物空间结构的边缘组合,得到建筑物[9]。此类方法由底层的边缘特征提取和搜索开始,逐层往上获得更高层的目标信息,在遮蔽少、建筑物独立和结构简单的情况下效果较好,但在边缘提取及其拟合等步骤中使用了较多阈值,这些阈值难以同时适应不同的建筑物类型和影像条件,实用性受到限制。
近年来,出现了诸如振荡器神经网络[10]、视觉图形主题模型[11]、时空推论[12]等新模型和方法,同时各种地理信息系统(Geographic Information System, GIS)数据库和先验知识库也不断更新和完善,虽然可以在一定程度上提高建筑物检测的准确性,但以上方法还存在如下问题:1) 需要大量的数据集作为训练样本;2) 需要复杂的模型建立过程;3) 需要庞大的数据库和专家系统的支持;4) 需要多种不同类型的遥感数据。总之,从实用性来看,受到较大的限制。
针对基于二维高分辨率遥感影像提取建筑物现有相关研究成果的不足以及数据限制,本文利用势直方图函数(Potential Histogram Function, PHF)和邻域总变分(Neighborhood Total Variation, NTV),提出一种无需多光谱及高程信息的方法,可以高效、鲁棒地自动识别和提取遥感影像中的建筑物。
1 算法原理本文的算法流程如图 1所示,主要包括基于NTV的区域生长分割、基于PHF的阴影提取和基于数学形态学计算邻接关系的建筑物提取三个部分。
以建筑物顶部具有均匀的光谱特性为前提条件,提出一种基于NTV的区域生长分割方法,实现了种子点的自动选取,并通过对分割对象进行后处理,提取候选建筑物。
设遥感影像为I(x, y),则总变分定义为:
$ {J_T}(I) = \int_{{D_I}} {\sqrt {\left| {\nabla I} \right|} \;{\rm{d}}x{\rm{d}}y} = \int_{{D_I}} {\sqrt {I_x^2 + I_y^2} \;{\rm{d}}x{\rm{d}}y} $ | (1) |
其中:
总变分实质上是函数I的梯度的模的积分。函数梯度反映了函数在某点的变化程度;梯度的模的积分反映了函数变化的激烈程度;总变分表示影像整体能量的变化,取值越小,代表着支持域内影像的变化越小,影像越平滑。由于建筑物顶部一般具有均匀光谱特征,使得遥感影像I中建筑物区域的总变分相对非建筑物区域要小得多,利用上述特征,定义如下的NTV:
$ \begin{array}{l} {J_{LT}}(I(x, y)) = \int_{D_{I(x, y)}^r} {\sqrt {{{\left| {\nabla I} \right|}^2} + \beta } \;{\rm{d}}u{\rm{d}}v} = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\int_{D_{I(x, y)}^r} {\sqrt {I_u^2 + I_v^2 + \beta } \;{\rm{d}}u{\rm{d}}v} \end{array} $ | (2) |
其中:DI(x, y)r是以I(x, y)为圆心、r为半径的支持域;β是大于0的可调参数,其取值与测量噪声的方差成正比,其作用是避免JLT(I(x, y))在Ix=Iy=0处导数不存在。NTV反映了影像中某一点附近的平滑程度,其值越小则越平滑。
考虑到邻域中不同位置的点可能具有不同的贡献,定义如下的加权NTV:
$ {J_{WLT}}(I(x, y)) = \int_{D_{I(x, y)}^r} {\lambda (u, v)\sqrt {I_u^2 + I_v^2 + \beta } \;{\rm{d}}u{\rm{d}}v} $ | (3) |
其中:λ(u, v)为对应点的权值,本文研究算法中,采用高斯核函数模板。
为了表征影像中的每一点的局部均质程度,根据加权NTV定义如下的似然函数:
$ LS(I({x_0}, {y_0})) = \frac{{\int_{D_{I({x_0}, {y_0})}^r} {\lambda (u, v)\;{\rm{d}}u{\rm{d}}v} }}{{\int_{D_{I({x_0}, {y_0})}^r} {\lambda (u, v)\sqrt {I_u^2 + I_v^2 + \beta } \;{\rm{d}}u{\rm{d}}v} }} $ | (4) |
LS(I(x0, y0))越大,表明I(x0, y0)属于建筑物区域的可能性越大。在计算各点的似然函数取值后,通过线性变换将其转换至区间[0, 255],如图 2(a)所示。考虑后续步骤中区域生长分割的准确性是以种子点的正确选取为前提的,因此对LS的二值化需要选取较大的阈值以保留最有可能属于建筑物区域的像素,本文通过对20幅遥感影像进行基于NTV的分割,比较阈值取值范围从100~250、间隔为10的效果,通过目视判断,发现阈值取值在220~240最为合适。这里选取了最佳范围的中间值230,二值化结果如图 2(b)所示。由于二值化后的结果中存在大量的孤立像素或小面积区域,将影响候选建筑物的提取精度,因此需要对其进一步处理,本文算法删除面积小于10的区域并对保留区域进行标注,结果如图 2(c)所示。图 2(d)为每个保留区域的质心,以每个质心为种子点进行区域生长分割,结果如图 2(e)所示。最后,利用矩形度RLW=L/W和长宽比RU=S1/S2作为约束条件对区域生长分割后的对象进行后处理,其中L和W分别代表分割对象的最小外接矩形的长和宽;S1/S2分别代表分割对象的面积和最小外接矩形的面积,最后保留同时满足矩形度大于阈值且长宽比小于阈值的分割对象(阈值的设置见表 2),得到候选建筑物(其中可能包含阴影),如图 2(f)所示。
阴影和非阴影的分割关键在于确定一个合适的阈值,而传统的分割方法大多采用人工设置。本文在文献[11]的基础上,提出一种基于PHF的自动阴影提取算法。
设I=[f(i, j)]M×N为M×N的灰度影像,G为灰度级,I的灰度统计直方图函数H(k)计算公式为:
$ H(k) = \frac{1}{{MN}}\sum\limits_{i = 0}^{M-1} {\sum\limits_{j = 0}^{N-1} {{\delta _{ij}}(k)} } $ | (5) |
其中:
${\delta _{ij}}(k) = \left\{ \begin{array}{l} 1, \;\;\;\;当\;f(i, j) = k\\ 0, \;\;\;\;其他 \end{array} \right.;\;\;k \in \{ 0, 1, \cdots, G - 1\} \;$ |
H(k)可以表征灰度值为k的像元在影像I中出现的概率。
定义PHF为:
$ {P_H}(k) = \frac{1}{{{P_{{\rm{max}}}}}}\sum\limits_{g = 0}^{G-1} {\frac{{H(g)}}{{1 + \alpha {{(g-k)}^2}}}} \; $ | (6) |
其中,
对图 2(f)分别计算灰度统计直方图函数H(k)和PHF,得到的结果如图 3(a)和图 3(b)所示,可以看出,PHF相当于对灰度统计直方图函数进行低通滤波。
根据阴影的灰度级较低的特征,通过构造PHF的一阶差分自动确定分割阈值,如式(7) 所示:
$ P(k) = {P_H}(k + 1)-{P_H}(k)\;;\;\;\;\;\;\;k \in \{ 0, 1, \cdots .G-1\} $ | (7) |
图 3(b)中的PHF曲线的一阶差分结果如图 4(a)所示,图 4(b)为图 4(a)的局部放大(椭圆包含部分)。将P(k)中第一个由负变正的阶跃点作为阴影提取的自动分割阈值T。图 4(b)中,P(63) < 0且P(64)>0,因此阈值T=63。
从I=[f(i, j)]M×N中提取的阴影为:
$SH(i, j) = \left\{ \begin{array}{l} 1, \;\;\;\;\;\;f(i, j) \le T\\ 0, \;\;\;\;\;其他 \end{array} \right.\;$ | (8) |
由图 5可见,本文的阴影提取算法对于图 5中的#1和#2测试影像均能得到较为可靠的阴影结果。
1.3 数学形态学计算基于形态学计算邻接关系进行建筑物提取的处理流程如图 6所示。
图 6中邻接关系的计算式为:
$ {B_c} = \left\{ \begin{array}{l} 1, \;\;\;\;\;{B_c} \cap S{H_m} \ne \emptyset \\ 0, \;\;\;\;其他 \end{array} \right. $ | (9) |
其中:Bc为任一候选建筑物;SHm为经过数学形态学处理后的所有阴影,取值为1表示保留Bc,取值为0表示删除Bc。当SHm为膨胀后的阴影对象时,保留的Bc为精细候选建筑物;当SHm为腐蚀后的阴影对象时,保留的Bc为精细阴影。
开运算、膨胀和腐蚀三种数学形态学方法均采用圆形结构元素,对应的半径分别记为R1、R2和R3。圆形结构元素ESTR可定义为:
$ {E_{{\rm{STR}}}} = {\rm{round}}(r) $ | (10) |
其中:round(*)为圆形生成函数; r为半径。
在候选建筑物和阴影的基础上,首先用数学形态学方法对阴影进行处理:
1) 开运算处理。目的是为了减小与建筑物阴影具有相似光谱特征的植被阴影及小面积人工地物的干扰,可以提高建筑物对象的查准率。
2) 膨胀处理。对开运算处理后的阴影进行膨胀处理,目的是为了弥补候选分割对象提取的不准确性,保证阴影和对应的建筑物有一定的重叠部分,可以提高建筑物检测的查全率。
3) 腐蚀处理。对开运算处理后的阴影进行腐蚀处理,目的是为了防止阴影本身被当作建筑物,可以提高建筑物提高的查准率。
然后,将与膨胀处理后的阴影有重叠部分的候选建筑物判定为精细候选建筑物;将与腐蚀处理后的阴影有重叠部分的候选建筑物判定为精细阴影。将精细候选建筑物减去精细阴影得到建筑物对象。
2 实验结果与分析 2.1 数据和评价方法本文的实验数据为深圳市2012年11月6日的PLEIADES影像,经过数据融合预处理,分辨率为0.5 m,该区域内建筑物类型较为典型,具有尺寸、方向和分布等特征的多样性。选择其中的红色波段进行处理,从中选取9幅具有代表性的测试影像,尺寸均为400×400像素,如图 8所示。#1测试影像中的建筑物呈多个角度分布且尺寸相差较大;#2测试影像中的建筑物的排列较为规则,但顶部光谱特性与周围的道路较为一致,容易受到道路的干扰;#3和#9测试影像中的建筑物阴影混叠在一起;#4和#6测试影像中的建筑物尺寸较大,且顶部包含小面积的规则同质区域;#7~#9测试影像中建筑物周围几乎被植被所覆盖(植被容易对阴影检测产生干扰);#5测试影像中的建筑物阴影不明显。
分别采用式(11)~(13) 的查准率Pre、查全率Rec和F1(将查准率和查全率联合量化为一个数值,取值越大,表明算法性能越好),对每幅测试影像进行评价。
$ Pre = \frac{{TP}}{{TP + FP}} \times 100\% \; $ | (11) |
$ Rec = \frac{{TP}}{{TP + FN}} \times 100\% \; $ | (12) |
$ {F_1} = \frac{{2 \times Pre \times Rec}}{{Pre + Rec}} \times 100\; $ | (13) |
其中:TP(True Positive)表示本文算法和人工标注法均提取到的建筑物数量(如果建筑物的部分区域被检测到,则认为建筑物是被正确提取的);FP(False Positive)表示本文算法检测到而人工标注法未检测到的建筑物数量;FN(False Negative)表示本文算法未检测到而人工标注法检测到的建筑物数量。
2.2 参数分析与取值本文算法共涉及8个参数,参数描述如表 1所示。
本文算法的综合平均查准率、查全率和F1分别随8个参数的变化曲线如图 7所示,图 7曲线是以图 8的#1~#5测试影像作为分析对象得到。
1) 参数α和RU均是先使F1曲线达到一个最大值,然后急剧下降至0,其原因在于α达到一定数值(图 7(a)中α=2.5) 后,将导致阴影完全消失,结果是提取不到任何建筑物;同样,当RU达到一定数值(图 7(d)中RU=0.75) 后,候选建筑物的数量为0,也提取不到任何建筑物。
2) 根据对式(8) 的分析,LS越大,属于建筑物区域的可能性越大,LS取值大于TBW的像素点被保留,因此,随着TBW由小到大的变化过程,建筑物被正确提取的可能性越来越大,但超过一定的取值(图 7(b)中TBW=230) 后,正确提取的建筑物数量越来越少,导致查准率和查全率下降。
3) 参数Tseg控制区域生长分割算法中一个像素的8邻域内与其相似的像素数量,Tseg越大,更多的像素被合并成同一对象。当Tseg取值小于2时,建筑物被过分割,每个分割对象的长宽比和矩形度无法满足约束条件,算法性能较低,随着Tseg取值增大,分割后的对象越来越接近于建筑物,建筑物顶部的光谱同质特性得到最大的利用,算法性能最好,当Tseg取值继续增加,出现欠分割现象,即相邻的多个建筑物或建筑物与相邻地物可能被检测为一个对象,导致查准率和查全率急剧下降。
4) 参数R1从小到大的变化,从一定程度上消除了与建筑物阴影具有相似光谱特征的植被阴影及小面积人工地物的干扰,到一定取值(图 7(f)中R1=5),效果达到最佳,但若再继续增大R1,阴影面积将不断被削弱,最终导致查全率下降得很快。
5) 参数RLW和R2使得F1曲线达到最大值后保持稳定,其原因在于RLW越大,约束条件越宽松,越不容易产生漏检和误检;此外,R2越大,保留的候选建筑物数量越多,漏检数量越少,因此当这两个参数达到一定数值后,查准率、查全率和F1均保持稳定。
6) 参数R3使得F1曲线上升到最大值后下降,其原因在于R3较小时,阴影被误检为建筑物的可能性大,导致查准率低,随着R3增大到一定的取值(图 7(h)中,R3=10),腐蚀后的阴影完全包含于候选建筑物中,对应的候选建筑物被删除,得到最大的查准率,当R3继续增大到足以将阴影完全腐蚀,与阴影所对应的候选建筑物得到保留,导致查准率下降。
综上分析可知,上述8个参数对算法性能的影响均有一定的规律性,即参数在变化过程中有某一取值或区间使算法总体性能(F1)达到最高或保持稳定,选择F1达到最大时对应的参数取值如表 2所示。
2.3 实验结果与分析结合表 2的参数取值,利用本文算法对9幅测试影像进行建筑物提取,结果如图 8所示,其中:矩形(即提取的建筑物的最小外接矩形)包含的区域代表TP建筑物的数量;实心圆代表FP建筑物的数量,即误检数;空心圆代表FN建筑物的数量,即漏检数。计算得到的性能结果如表 3所示。
如图 8所示,9幅测试影像中的大多数建筑物都能被正确提取,表明本文算法能适用于不同尺寸、分布及具有周围地物干扰的建筑物提取。表 3中列举的数据也证明了这一点。基于对象的性能统计数据中,9幅测试影像的平均查准率和查全率分别为97.71%和84.21%,F1为90.46。多数的测试影像均能达到100%的查准率,其中,#1、#3和#7测试影像各存在一个将道路片段误检为建筑物的对象,导致查准率有所下降,原因是部分道路片段具有均匀的光谱特性且满足一定的长宽比和矩形度。9幅测试影像都存在漏检建筑物,其原因在于:1) 建筑物被影像边界分割成不完整的片段,无法满足长宽比和矩形度的约束条件,如#2测试影像右上方的漏检建筑物;2) 建筑物顶部不具备光谱均一性,如#3、#4和#6测试影像中的漏检建筑物;3) 阴影没有被准确提取,导致在计算阴影和周围建筑物的重叠关系时造成漏检,如#3、#7和#9测试影像中的漏检建筑物。
需要重点说明的是:1) 阴影是本文算法利用的关键特征,当阴影漏检时,将直接导致建筑物的漏检,当阴影误检时,由于后续步骤利用数学形态学的腐蚀和膨胀分别对阴影进行处理,因此只要不把建筑物误检为阴影,将对建筑物提取结果无影响;2) 本文方法还无法完整地提取建筑物,存在如#1测试影像中建筑物漏检较为严重的现象,这主要取决于分割算法的性能,将在以后的研究中,综合颜色、边缘和像元之间的空间关系设计性能更优的分割算法,提高建筑物提取的准确性。
综合以上分析可知,本文方法对于多数的场景具有较好的提取效果,但对顶部光谱不均匀、周围阴影不明显、轮廓难以用矩形表示的建筑物提取还存在一定的局限性,同时,值得注意的是,表 3中9幅测试影像的性能数据是基于表 2的同一配置参数下获得的,如果分别对每幅测试影像配置各自最佳的参数,总体性能将会有一定程度的提高。
2.4 方法对比鉴于水平集是目前广泛应用且效果较好的分割方法,而基于不变颜色特征和阴影的建筑物提取方法与本文的方法结构基本一致,因此选择这两种方法进行对比具有一定的代表性和可比性。基于#1~#9测试影像,对比结果如表 4所示, 其中:方法1为利用水平集进行分割的建筑物提取方法[13],方法2为基于不变颜色特征和阴影的建筑物提取方法[14]。
由表 4可知:相比方法1和方法2,本文方法在总体性能F1方面均提高了10%以上。其中,方法1具有较高的查全率,但查准率相对较低,其原因在于方法1的水平集初始曲线由均匀分布在整幅影像上的圆形组成,当建筑物边界较为明显及分布较为稀疏时,能够较好地检测到建筑物,相反,活动轮廓线难以准确收敛至建筑物边界;方法2先利用不变颜色特征进行分割,提取建筑物(为了保证满足此条件,随机将图 8中9幅测试影像中的两处建筑物进行红色填充),方法2达到88.65%的查准率,但条件要求较为苛刻。
3 结语相对于现有多数建筑物提取方法存在的需要多光谱数据、高程数据以及人工设置大量参数等问题,本文方法只需利用遥感影像的灰度信息,降低了对遥感影像维度和波段的要求,通过基于阴影和候选建筑物之间的邻接关系筛选并提取建筑物,且该方法所涉及的参数容易确定,为建筑物提取提供了一种有效方案。实验结果表明,本文方法能够很好地适应于具有顶部光谱一致性的矩形或类矩形建筑物的提取,查准率高达97.71 %,通过目视解译,提取的建筑物位置和轮廓较为准确,不足之处是对于边缘模糊的、形状不规则和顶部光谱特征不均匀的建筑物提取的查全率相对较低。下一步研究工作将考虑建筑物位置和边缘的准确性评价问题并扩展至大幅面影像中建筑物提取,提高查全率,使本文算法能够应用于城市建筑物提取的工程中。
[1] | LEE D S, SHAN J, BETHEL J S. Class-guided building extraction from Ikonos imagery[J]. Photogrammetric Engineering & Remote Sensing, 2003, 69(2): 143-150. |
[2] | SUMER E, TURKER M. An adaptive fuzzy-genetic algorithm approach for building detection using high-resolution satellite images[J]. Computers, Environment and Urban Systems, 2013, 39: 48-62. doi: 10.1016/j.compenvurbsys.2013.01.004 |
[3] | VNSALAN C, BOYER K L. A system to detect houses and residential street networks in multispectral satellite images[J]. Computer Vision and Image Understanding, 2005, 98(3): 423-461. doi: 10.1016/j.cviu.2004.10.006 |
[4] | SHRIVASTAVA N, RAI P K. Automatic building extraction based on multiresolution segmentation using remote sensing data[J]. Geographia Polonica, 2015, 88(3): 407-421. doi: 10.7163/GPol |
[5] | GRINIAS I, PANAGIOTAKIS C, TZIRITAS G. MRF-based segmentation and unsupervised classification for building and road detection in peri-urban areas of high-resolution satellite images[J]. ISPRS Journal of Photogrammetry & Remote Sensing, 2016, 122: 145-166. |
[6] | 施文灶, 毛政元. 基于图分割的高分辨率遥感影像建筑物变化检测研究[J]. 地球信息科学学报, 2016, 18(3): 423-432. ( SHI W Z, MAO Z Y. The research on building change detection from high resolution remotely sensed imagery based on graph-cut segmentation[J]. Journal of Geo-Information Science, 2016, 18(3): 423-432. ) |
[7] | YU B, WANG L, NIU Z. A novel algorithm in buildings/shadow detection based on Harris detector[J]. Optik-International Journal for Light and Electron Optics, 2014, 125(2): 741-744. doi: 10.1016/j.ijleo.2013.07.043 |
[8] | NGO T T, COLLET C, MAZET V. Automatic rectangular building detection from VHR aerial imagery using shadow and image segmentation[C]//Proceedings of the 2015 IEEE International Conference on Image Processing. Piscataway, NJ:IEEE, 2015:1483-1487. |
[9] | DHULE J J, DHURPATE N B, GONGE S S, et al. Edge detection technique used for identification of cracks on vertical walls of the building[C]//Proceedings of the 2015 International Conference on Computing and Network Communications. Piscataway, NJ:IEEE, 2015:263-268. |
[10] | LIU C, SHI B Q, YANG X, et al. Automatic buildings extraction from LiDAR data in urban area by neural oscillator network of visual cortex[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2013, 6(4): 2008-2019. doi: 10.1109/JSTARS.2012.2234726 |
[11] | LI Y S, TAN Y H, TIAN J W. Urban building extraction via visual graphical topic model[C]//Proceedings of the 2014 IEEE International Geoscience and Remote Sensing Symposium. Piscataway, NJ:IEEE, 2014:4746-4749. |
[12] | QIN R J, TIAN J J, REINARTZ P. Spatiotemporal inferences for use in building detection using series of very-high-resolution space-borne stereo images[J/OL]. International Journal of Remote Sensing, 2015.[2016-10-09]. http://rongjunqin.weebly.com/uploads/2/3/2/3/23230064/spatiotemporal_inferences_for_use_in_building_detection_using_series_of_very-high-resolution_space-borne_stereo_imagesauthorcopy.pdf. |
[13] | AHMADI S, ZOEJ M J V, EBADI H, et al. Automatic urban building boundary extraction from high resolution aerial images using an innovative model of active contours[J]. International Journal of Applied Earth Observation & Geoinformation, 2010, 12(3): 150-157. |
[14] | SIRMACEK B, UNSALAN C. Building detection from aerial images using invariant color features and shadow information[C]//Proceedings of the 2008 23rd IEEE International Symposium on Computer and Information Sciences. Piscataway, NJ:IEEE, 2008:1-5. |