CN105719306B - 一种高分辨率遥感影像中的建筑物快速提取方法 - Google Patents

一种高分辨率遥感影像中的建筑物快速提取方法 Download PDF

Info

Publication number
CN105719306B
CN105719306B CN201610050661.4A CN201610050661A CN105719306B CN 105719306 B CN105719306 B CN 105719306B CN 201610050661 A CN201610050661 A CN 201610050661A CN 105719306 B CN105719306 B CN 105719306B
Authority
CN
China
Prior art keywords
pixel
point
angle
straight line
array
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201610050661.4A
Other languages
English (en)
Other versions
CN105719306A (zh
Inventor
班瑞
郑延召
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Zhengzhou Hengzheng Electronic Technology Co Ltd
Original Assignee
Zhengzhou Hengzheng Electronic Technology Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Zhengzhou Hengzheng Electronic Technology Co Ltd filed Critical Zhengzhou Hengzheng Electronic Technology Co Ltd
Priority to CN201610050661.4A priority Critical patent/CN105719306B/zh
Publication of CN105719306A publication Critical patent/CN105719306A/zh
Application granted granted Critical
Publication of CN105719306B publication Critical patent/CN105719306B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G06T5/70—
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20024—Filtering details

Abstract

本发明提供了一种高分辨率遥感影像中的建筑物快速提取方法，先对图像进行灰度提取，然后通过高斯滤波进行去噪，接着使用基于Otsu方法的Canny自适应处理方法提取到最符合实际情况的边缘二值图像，再使用基于Freeman链码的方法快速提取二值图像中的直线，最后取到角点，连接形成建筑物轮廓，最终完成高分辨率遥感影像中建筑物的快速自动提取；本发明提供的高分辨率遥感影像中的建筑物快速提取方法，可实现快速提取遥感影像中的建筑物，处理过程中不需人为干预，提取速度快，精度高。

Description

一种高分辨率遥感影像中的建筑物快速提取方法
技术领域
本发明属于高分辨率遥感影像信息提取领域，涉及对高分辨遥感影像的中人造地物的识别和提取，特别是涉及对那些边缘由直线、弧线组成的建筑物提取的技术领域。
背景技术
随着遥感影像的分辨率不断提高，影像中的信息更加复杂，地物的纹理形状信息更加多样化，而建筑物在大小形状上各有区别，再加上周边其他地物的混淆、阴影的干扰等原因容易造成在影像上建筑物形状的变化，因此要对建筑物的形状做精确的描述是件非常困难的事情。近年来，许多研究人员在利用遥感影像进行建筑物提取方面作了大量研究工作，各种新方法、新理论层出不穷，其中基于建筑物轮廓大多是由矩形或是弧线构成的特点，通过特定处理方法，寻找图像中的边缘、线条、角点等特征来识别建筑物是主流研究方向之一。
发明内容
本发明提供了一种高分辨率遥感影像中的建筑物快速提取方法，可快速寻找图像中的边缘、线条、角点等特征来识别建筑物；
一种高分辨率遥感影像中的建筑物快速提取方法，其中，包括如下步骤：
步骤1)、读取遥感影像，提取灰度图像：
遥感图像中每个像素点的色彩具有R、G、B三个分量,P(r,g,b)代表一像素点,其中r代表该像素点的R分量的值，g代表该像素点的G分量的值，b代表该像素点的B分量的值；
灰度函数H(P)的值表示像素点P(r,g,b)的灰度值：
H(P)＝0.3r+0.59g+0.11b 式(1)
用式(1)按照从左到右、从上到下的顺序依次扫描遥感图像中的每个像素点，得到每个像素点的灰度值，按照式(1)扫描时同样的顺序组成第一灰度图像；
步骤2)、使用高斯滤波对第一灰度图像进行降噪处理：
2a)、计算大小为(2k+1)×(2k+1)的高斯模板U(x，y)：
其中，k为正整数且k≥1，x，y分别代表高斯模板中元素的横坐标和纵坐标；σ＞0，为第一灰度图像的平滑程度参数；
2b)、取k＝1，大小为3×3的高斯模板U(x，y)中9个元素的模板值分别为m1、m2、m3、m4、m5、m6、m7、m8、m9；
第一灰度图像中，除左、上、右、下四条宽度为一像素的边缘外任一像素点P点的8邻像素，即紧邻像素点P的左上、上、右上、左、右、左下、下、右下8个方向上的像素；像素点P与8邻像素的灰度值按照第一灰度图像中从左至右、从上到下的顺序依次为g1、g2、g3、g4、g5、g6、g7、g8、g9；
用该高斯模板U(x，y)按照从左到右、从上到下的顺序依次扫描第一灰度图像中除左、上、右、下四条宽度为一像素的边缘外的每一个像素P，各像素点P由高斯模板U(x，y)做高斯滤波的方法见式(3)；
所述的U0表示以高斯模板U(m，n)做高斯滤波后所到的像素点P的灰度值；S为中间变量，表示像素点P和其8邻像素的灰度值的总和；
按照式(3)依从左到右、从上到下的顺序对第一灰度图像中像素点P进行依次扫描时，对当前所扫描的像素点P由高斯模板U(x，y)所确定的邻域内所有像素的加权平均灰度值替代当前像素点P的灰度值，扫描结束后得到第二灰度图像；
步骤3)、使用基于Otsu方法的自适应Canny算法提取图像中的边缘；
31)：用一阶偏导的有限差分来计算梯度的幅值和方向；
采用Canny算法,其中所采用的卷积模板为式(4)；
其中M1、M2分别为第二灰度图像x、y轴方向上的一阶偏导数矩阵,依照所述一阶偏导数矩阵M1、M2按照从左到右、从上到下的顺序依次扫描第二灰度图像中的每一个像素点，得到每个像素点x、y方向上的梯度幅值Gx、Gy：
Gx＝[f(x+1,y)-f(x,y)+f(x+1,y+1)-f(x,y+1)]/2 式(5)
Gy＝[f(x,y)-f(x,y+1)+f(x+1,y)-f(x+1,y+1)]/2 式(6)
x、y表示所处理像素点的横坐标、纵坐标值，f(x,y)表示第二灰度图像中坐标为(x,y)的像素点的灰度值，该像素点的梯度幅值G(x,y)和方向角θ(x,y)由式(7)和式(8)处理得到；
由式(7)和式(8)对第二灰度图像中每一个像素点按照从左到右、从上到下的顺序依次扫描，得到第二灰度图像中每一个像素点梯度幅值和方向角；
32)：非极大值抑制；
经过上述步骤31)的处理，得到了每个像素的梯度幅值G(x,y)和方向角θ(x,y)，为了得到最真实的边缘，须保留每个像素点在其梯度方向上的极大值，而删掉当前像素点梯度方向上的其他值，即非极大值抑制；对第二灰度图像中每个像素点按照从左到右、从上到下的顺序进行非极大值抑制，即处理得到第三灰度图像，第三灰度图像显示了遥感影像中建筑物影像的边缘；
非极大值抑制的具体处理方法如下；
判断是否是极大值，需要根据当前扫描像素点梯度方向所在的区间，在x轴和y轴方向分别进行插值计算，然后再进行比较,以找到最大值，梯度方向所在的区间指的是梯度方向α所在的角度范围，插值是为了更合理地判断当前像素点的梯度方向：
①当以及时，令Gp＝G3w+G2(1-w),Gn＝G7w+G8(1-w)；所述的w为中间变量；
②当以及时，令Gp＝G7w+G4(1-w),Gn＝G3w+G6(1-w)；
③当以及时，令Gp＝G6w+G9(1-w),Gn＝G4w+G1(1-w)；
④当以及时，令Gp＝G2w+G1(1-w),Gn＝G8w+G9(1-w)；
其中，G5为当前处理像素的梯度幅值，Gk(1≤k≤4,6≤k≤9)表示8邻像素的梯度幅值，Gp表示梯度方向上前一个像素的梯度幅值，Gn表示梯度方向上后一个像素的梯度幅值；
当梯度方向α的值在上述相应的区间内时，当且仅当目前所处理像素的梯度幅值G5大于前述①～④中各插值方法计算出的Gp和Gn时，G5为极大值，也就表明当前处理像素是边缘点；
33)、计算最优高低阈值
读取第三灰度图像的宽度为W像素、高度为H像素，则第三灰度图像中的像素总数N为：
N＝W×H 式(9)
设第三灰度图像中每一像素的梯度幅值G(x,y)的范围为[0,G]，读取每一像素的梯度幅值G(x,y)，找梯度幅值等于确定值为i的像素，得到梯度幅值等于i的像素个数为Ni个，i∈[0,G]；设T1为最优低阈值，T2为最优高阈值；
非边缘点总个数N1(T1,T2)为：
其中，Ni是梯度幅值等于i的像素的个数；
非边缘点所占比例值ω1(T1,T2)为：
潜在边缘点总个数N2(T1,T2)为：
潜在边缘点所占比例值ω2(T1,T2)为：
真边缘点总个数N3(T1,T2)为：
真边缘点所占比例值ω3(T1,T2)为：
非边缘点的平均梯度μ1(T1,T2)、潜在边缘点的平均梯度μ2(T1,T2)和真边缘点的平均梯度μ3(T1,T2)如下：
第三灰度图像总的梯度均值μs为：
μs＝μ1(T1,T2)ω1(T1,T2)+μ2(T1,T2)ω2(T1,T2)+μ3(T1,T2)ω3(T1,T2) 式(17)
非边缘点、潜在边缘点和真边缘点这三类像素的类间方差σ2(T1,T2)为：
σ2(T1,T2)＝ω1(T1,T2)(μ1(T1,T2)-μS)2+ω2(T1,T2)(μ2(T1,T2)-μS)2+ω3(T1,T2)(μ3(T1,T2)-μS)2
式(18)
使得σ2(T1,T2)的输入参数T1、T2在总范围[1,G-1]依次取值，当类间方差σ2(T1,T2)为最大值时，此时的T1、T2的值即为最优低阈值T1和最优高阈值T2；
34)：用双阈值算法检测和边缘连接；
根据上一步33)中得到的最优低阈值T1和最优高阈值T2对第三灰度图像进行过滤，像素点的梯度幅值如果大于高阈值则认为当前像素点必然是边缘点，即真边缘点，如果当前处理像素点的梯度幅值小于低阈值，则认为当前像素点必然不是边缘点；处于低阈值和高阈值之间的像素点是潜在的边缘点，即假边缘；将第三灰度图像根据高阈值过滤得到的边缘点链接成轮廓，当到达轮廓的端点时，在端点的8邻点寻找潜在的边缘点，再根据所找到的潜在的边缘点,在其8邻点中循环处理收集新的边缘点，如此处理整个第三灰度图像，直至边缘闭合，得到第四灰度图像；
步骤4)、基于Freeman链码的直线提取方法：
单条直线L与Y轴的夹角θ范围为0～π,共有八种情况的直线L：
4)-①θ＝0或θ＝π；
4)-②
4)-③
4)-④
4)-⑤
4)-⑥
4)-⑦
4)-⑧
此外还有两条直线有交叉的情况；
设po表示左右均无目标像素的单像素；
ps表示连续像素个数≥2的像素组；
整型常量R表示目标之间可以重叠的像素数；
整型常量L表示一条直线内不包含重叠像素的最少像素数；
整型变量W(P0)表示单像素po在行或列内的位置，行内顺序为从左至右，列内顺序为从上到下；
整型变量S(Ps)表示ps在行或列内的起始位置，行内顺序为从左至右，列内顺序为从上到下；
整型变量E(Ps)表示ps在行或列内的终止位置；
数组Ax用于存储x轴方向上的所有po和ps信息，结构包含：目标类型、单像素位置，其中目标类型为po或ps；
数组Ay用于存储y轴方向上的所有po和ps信息，结构包含：目标类型、单像素位置、像素组在列内的起始点位置，其中目标类型为po或ps；
整型变量R(Po)、R(Ps)分别表示单像素po、像素组ps在数组Ax内的行号；
整型变量C(Po)、C(Ps)分别表示po、ps在Ay内的列号；
数组AL表示存储一条直线内的所有po和ps信息；
数组GL用于以(θ,ρ)的形式存储检索到的所有直线信息；
第四灰度图像中的直线具备以下特点：
4)-(一)：直线都是由po和ps组成；
4)-(二)：θ角越接近或其线条中的po和ps越多；θ角等于或时，直线仅由po组成；θ角越接近0、或π，直线中po和ps越少，但ps内的像素个数越多；θ角等于0、或π时，直线仅由ps组成；
4)-(三)：交叉的直线，共用其中一个po或ps；
4)-(四)：直线都是由po和ps同时沿x轴、y轴两个方向连接形成；x轴方向上直线的长度逐渐增加，y轴方向上逐渐产生一定的斜率，或者，y轴方向上直线的长度逐渐增加，x轴方向上逐渐产生一定的斜率；
4)-(五)：当θ角为0、或π时，直线仅包含一个ps；
令R＝1，以检索的直线，第四灰度图像中的直线提取通过下述步骤实现：
4)-⑴从第四灰度图像的左上角像素开始，按照逐行、逐列的顺序，即从左到右、从上到下的顺序，依次扫描第四灰度图像中的每一个像素点，来检索第四灰度图像中的目标点，目标点即边缘点；提取到的结果保存在Ax和Ay中；
4)-⑵取出Ax中的一个目标At，将At插入数组AL，其中，第一次抽取的目标At是随机抽取，以后抽取依照下步即第4)-⑶步所示，将抽取的目标At按照抽取的先后顺序插入数组AL；
4)-⑶沿x轴负方向，查找At的下一个目标At+1，At+1必须同时符合以下条件：
4)-⑶-①R(At+1)＝R(At)+1；
4)-⑶-②At为po时：
4)-⑶-②-Ⅰ：如果At+1为单像素，W(At+1)＝W(At)-R；
4)-⑶-②-Ⅱ：如果At+1为像素组，W(At)-R≤E(At+1)≤W(At)；
4)-⑶-③At为像素组时：
4)-⑶-③-Ⅰ：如果At+1为单像素，W(At+1)＝S(At)-R；
4)-⑶-③-Ⅱ：如果At+1为像素组，S(At)-R≤E(At+1)≤S(At)；
4)-⑶-④At+1未访问过；
如果找到At+1，将At+1标记为已访问且将其插入数组AL；令At＝At+1，递归执行该步骤4)-⑶处理过程；如果找不到At+1，结束查找；
4)-⑷令数组AL中的第一个元素为AL0，最后一个元素为ALn，如果数组AL中只有一个元素时：若该元素长度大于等于L，该元素长度即像素数，此条直线(此条直线为AL中该元素(θ,ρ)所表示的一条直线)起点为S(AL0)，终点为E(AL0)；若长度小于L,此次检索失败；如果数组AL中不少于2个元素，此条直线起点为E(AL0)，终点为S(ALn)；
4)-⑸沿x轴正方向，查找At的下一个目标At+1，At+1必须同时符合以下条件：
4)-⑸-①R(At+1)＝R(At)+1；
4)-⑸-②At为单像素时：
4)-⑸-②-Ⅰ：如果At+1为单像素，W(At+1)＝W(At)+R；
4)-⑸-②-Ⅱ：如果At+1为像素组，W(At)≤S(At+1)≤W(At)+R；
4)-⑸-③At为像素组时：
4)-⑸-③-Ⅰ：如果At+1为单像素，W(At+1)＝E(At)+R；
4)-⑸-③-Ⅱ：如果At+1为像素组，E(At)≤S(At+1)≤E(At)+R；
4)-⑸-④At+1未访问过；
如果找到At+1，将At+1标记为已访问且将其插入数组AL；令At＝At+1，递归执行该步骤4)-⑸的处理过程；如果找不到At+1，结束查找，执行下一步4)-⑹；
4)-⑹令数组AL中的第一个元素为AL0，最后一个元素为ALn；最终提取的直线分别以下述三种情况存储：
4)-⑹-①如果数组AL中只有一个元素且其长度大于等于L，该元素长度即为像素数，此条直线起点为S(AL0)，终点为E(AL0)；
4)-⑹-②如果数组AL中的元素数量≥2且E(AL0)-S(AL0)≥L，则此条直线起点为S(AL0)，终点为E(ALn)；
4)-⑹-③如果数组AL中的元素数量≥2且E(AL0)-S(AL0)＜L，放弃此条直线；
4)-⑺将检索到的直线插入到数组GL中,清空数组AL；
4)-⑻抽取数组Ax中一个未访问过的目标，令其等于At,递归执行上述对第四灰度图像中的直线提取步骤中的第4)-⑶步，直至数组Ax中所有目标都被标记为已访问；
4)-⑼扫描数组Ay中的目标，检索和内的直线；取出数组Ay中的一个目标At，将At插入数组AL；
4)-⑽沿y轴正方向，查找At的下一个目标At+1，At+1必须同时符合以下条件：
4)-⑽-①C(At+1)＝C(At)+1；
4)-⑽-②At为po时：
4)-⑽-②-Ⅰ：如果At+1为单像素，W(At+1)＝W(At)-R；
4)-⑽-②-Ⅱ：如果At+1为像素组，W(At)+R≤E(At+1)≤W(At)；
4)-⑽-③At为像素组时：
4)-⑽-③-Ⅰ：如果At+1为单像素，W(At+1)＝S(At)-R；
4)-⑽-③-Ⅱ：如果At+1为像素组，S(At)+R≤E(At+1)≤S(At)；
4)-⑽-④At+1未访问过；
如果找到At+1，将At+1标记为已访问且将其插入数组AL；令At＝At+1，递归执行该步4)-⑽；如果找不到At+1，结束查找，执行下一步4)-⑾；
4)-⑾令数组AL中的第一个元素为AL0，最后一个元素为ALn，如果数组AL中只有一个元素且该元素长度小于等于L，该元素长度即为像素数，则此条直线起点为S(AL0)，终点为E(AL0)；否则，如果数组AL中有至少2个元素，此条直线起点为E(AL0)，终点为S(ALn)；
4)-⑿沿y轴负方向，查找At的下一个目标At+1，At+1必须同时符合以下条件：
4)-⑿-①C(At+1)＝C(At)+1；
4)-⑿-②At为po时：
4)-⑿-②-Ⅰ：如果At+1为单像素，W(At+1)＝W(At)+R；
4)-⑿-②-Ⅱ：如果At+1为像素组，W(At)-R≤S(At+1)≤W(At)；
4)-⑿-③At为像素组时：
4)-⑿-③-Ⅰ：如果At+1为单像素，W(At+1)＝S(At)+R；
4)-⑿-③-Ⅱ：如果At+1为像素组，E(At)-R≤S(At+1)≤E(At)；
4)-⑿-④At+1未访问过；
如果找到At+1，将At+1标记为已访问且将其插入AL；令At＝At+1，递归执行该步4)-⑿；如果找不到At+1，结束查找，执行下一步4)-⒀；
4)-⒀令数组AL中的第一个元素为AL0，最后一个元素为ALn；最终提取的直线以下述三种情况存储：
4)-⒀-①如果数组AL中只有一个元素且其长度大于等于L，该元素长度即为像素数，此条直线起点为S(AL0)，终点为E(AL0)；
4)-⒀-②如果数组AL中的元素数量≥2且E(AL0)-S(AL0)≥L，则此条直线起点为S(AL0)，终点为E(ALn)；
4)-⒀-③如果AL中的元素数量≥2且E(AL0)-S(AL0)＜L，放弃此条直线；
4)-⒁将AL中的内容以Hough变换方法处理求出θ角和ρ值，插入到数组GL中,清空数组AL；
4)-⒂抽取Ay中一个未访问过的目标，令其等于At,递归执行上述第4)-⑽步，直至Ay中所有目标都被标记为已访问；
至此，直线的提取全部完成，数组GL中存储了第四灰度图像中0≤θ≤π的所有直线；
步骤5)、直线的归并
设阈值LR、θR、ρR分别表示两条直线的端点、θ角、ρ值的容差，令直线L1的起点、终点、θ角和ρ值分别为S1、E1、θ1和ρ1，直线L2的起点、终点、θ角、ρ值分别为S2、E2、θ2和ρ2；从第四灰度图像中提取到的直线会有部分具备以下特点之一：
5)-⑴θ1＝θ2，ρ1＝ρ2，MIN{|S1-E2|,|S2-E1|}≤LR；
5)-⑵|θ1-θ2|≤θR，|ρ1-ρ2|≤ρR；
5)-⑶ρ1＝ρ2，|θ1-θ2|≤θR；
5)-⑷θ1＝θ2，|ρ1-ρ2|≤ρR；
这四类直线需要进行归并，以减少重复的直线；
设合并后的直线为L3，其起点、终点、θ角和ρ值为S3、E3、θ3和ρ3，针对以上四种情况，直线归并如下处理：
5)-①当MIN{|S1-E2|,|S2-E1|}＝|S1-E2|时，取E3＝E1，S3＝S2；否则，当MIN{|S1-E2|,|S2-E1|}≠|S1-E2|时，E3＝E2，S3＝S1；
5)-②当ρ1＝ρ2，|θ1-θ2|≤θR时，取ρ3＝ρ1＝ρ2,
5)-③当L1‖L2时，取θ3＝θ1＝θ2,符号“‖”表示直线间的平行；
四类直线归并前，两条直线分别表示为L1(S1,E1)和L2(S2,E2)，归并后的直线表示为L3(S3,E3)；
步骤6)、图像中的角点提取
(61)：角度阈值过滤
定义角度阈值Tθ，判断两条直线的交叉角度在之间时，认为是直角，直线的交点即为建筑物的角点；
(62)：间断直线过滤
在背景复杂且噪声较多的遥感影像中，部分边缘直线不会被完整的获取到，出现间断的线段，设置阈值TL，确定在间断长度为阈值TL以内时，认为此角点可用，否则，间断的直线不做为直角边使用；
(63)：角分裂；
根据直线交叉的情况，分裂成多个直角；
(64)：计算角点；
已知两条直线L1:A1x+B1y+C1＝0；L2:A2x+B2y+C2＝0，两直线交点是坐标为(x,y)的像素点；
通过式(19)求出两直线的交点即角点的坐标x值、y值；
为了便于后续的角点归并和排序，在记录角点时，还需要考虑到角的开口方向；在得到角点的坐标值后，令：角的两条边为e1和e2，以角点为起点、远离角点的端点为终点的两个向量分别为和角方向如下表示：
的方向即为角的方向，所述角以θr表示，角的范围为：0≤θr≤2π；
(65)：角点过滤和归并；
由于遥感影像中的复杂背景，提取目标周围会出现一些伪边缘，需要进行角点过滤和归并：
如果两角点的开口方向相同：
令角点P1的两边长和夹角θr分别为E1、E2、θr1，角点P2的两边长和夹角分别为E3、E4、θr2，设置距离差阈值TD、角度容差Tθ；当角点P1和角点P2的距离D1≤TD且|θr1-θr2|≤Tθ时，对两角点进行归并；如果E1E2≥E3E4，则删除角点P2；如果E1E2＜E3E4，删除角点P1；
如果两角点的开口方向相对，且两角点之间出现了第三个角点P3：
令角点P3的两边长、夹角分别为E5、E6、θr3，当角点P1和角点P3的距离D2＞TD，且|θr1-θr3|≤Tθ时，以角点P2的θr2方向为检索方向，在θr2方向上寻找θr2-Tθ≤θr≤θr2+Tθ的角点，找到角点P1和角点P3后，如果E1E2≥E5E6，则删除角点P3；否则，当E1E2小于E5E6时，删除角点P1；
步骤7)、归属同一建筑物的角点匹配
建筑物角点的检索可以使用如下步骤完成：
(71)从角点集A中取一个角点As,如果As已访问过，则重新取；角点集A为经过角点过滤和归并后得到的角点集合；
(72)：沿角点As的一条边的方向寻找下一个与角点As匹配的角点As+1，匹配的条件是：
角点As+1的一条边与的方向相反，且误差在设定的角度容差θE以内；
(73)：将当前角点As标记为已访问，并将其保存在数组B内；
(74)：令As＝As+1，e1＝e4，递归执行步骤(72),直到匹配到最初的As角点为止；
(75)：遍历角点集A中所有角点，直到所有角点均为已访问；匹配结束；
对于未完成循环匹配的一组角点，由如下步骤处理推断出其余角点，构成由至少4个角点组成的矩形建筑：
⑴单独角点：由角点As在角点集A中匹配不到As+1时，以角点As的两条边以及角方向的三个端点，构建一个四边形的建筑物轮廓；
⑵两个角点：角点As在角点集A中匹配到As+1，而As+1匹配不到其它角点时，以As和As+1的两条边的端点为另两个角点，构建一个四边形的建筑物轮廓；
⑶多于两个角点：匹配到两个以上角点，但没有形成环路时，以最初的角点As的边与最后一个角点的边的交点做为最后一个角点；
步骤8)、扫描结束，输出处理后的图像，结束程序。
本发明提供的高分辨率遥感影像中的建筑物快速提取方法，步骤1)是提取遥感影像的灰度图像；步骤2)是通过高斯滤波，去除图像中的一些噪点，使图像变的更加平滑，为后面的边缘提取做好准备；步骤3)是使用基于Otsu方法的canny自适应处理方法进行边缘提取；步骤4)是采用基于Freeman链码的自主创新的快速直线提取处理方法，获取边缘图像中的所有直线，利用Hough变换中的点-正弦曲线对偶原理，将直线以二元参数形式存储和处理；步骤5)是根据直线的二元参数值，对直线进行归并和过滤，找到最优直线集；步骤6)是根据提取到的直线线段，检索图像中的角点；步骤7)是根据角方向、直角边、阈值θE，对归属同一建筑物的角点进行匹配，最终绘制出完整的建筑物轮廓；本发明提供的高分辨率遥感影像中的建筑物快速提取方法，可实现快速提取遥感影像中的建筑物，处理过程中不需人为干预，提取速度快，精度高。
附图说明
图1为本发明一种高分辨率遥感影像中的建筑物快速提取方法的处理步骤流程图；
图2为高、低阈值获取过程处理步骤流程图：
图3为直线检索处理步骤流程图；
图4为像素点P的8邻像素以及梯度方向角的示意图；
图5为两条直线L1、L2分别与Y轴的夹角示意图；
图6为θ＝0或θ＝π的直线；
图7为的直线；
图8为的直线；
图9为的直线；
图10为的直线；
图11为的直线；
图12为的直线；
图13为的直线；
图14为两条直线交叉情况下的图示；
图15为直线形成时像素连接方式示意图；
图16为At为单像素时与At+1的位置关系示意图；
图17为At为像素组时与At+1的位置关系示意图；
图18为数组AL中直线的起点和终点示意图；
图19为间断直线的交点三种情况示意图；
图20为三种直线交叉情况的角分裂示意图；
图21为角方向示意图；
图22为最典型的有3类角点归并情况示意图；
图23为凸型建筑物的角点分布及角方向示意图；
图24为样本图像经过本发明方法步骤2)高斯滤波处理后示意图；
图25为样本图像经过本发明方法步骤3)边缘提取处理后的示意图；
图26为样本图像经过本发明方法步骤6)角点提取处理后的示意图。
具体实施方式
本发明提供了一种遥感影像中的建筑物快速提取方法，具体流程如图1所示，下面提供一具体实施例：
步骤1)、读取遥感影像，提取灰度图像；
遥感图像中每个像素点的色彩由R、G、B三个分量共同决定，P(r,g,b)代表遥感图像中的一像素点,其中r代表该像素点的R分量的值，g代表该像素点的G分量的值，b代表该像素点的B分量的值，每个分量在内存所占的位数共同决定了图像深度，即每个像素点所占的字节数；以常见的24位深度彩色RGB图来说，其三个分量各占1个字节，这样每个分量取值范围为[0,255],令H(P)为灰度函数,其值表示像素点P(r,g,b)的灰度值，灰度函数H(P)的计算方法如下：
H(P)＝0.3r+0.59g+0.11b 式(1)
用式(1)按照从左到右、从上到下的顺序依次扫描遥感图像中的每个像素点，得到每个像素点的灰度值，按照式(1)扫描时同样的顺序组成第一灰度图像；
步骤2)、使用高斯滤波对第一灰度图像进行降噪处理；
噪声是由于拍摄成像环境、设备、天气等随机因素造成；可以通过高斯滤波减弱其对原图的影响；
高斯滤波是一种线性平滑滤波，适用于消除高斯噪声，广泛应用于图像处理的降噪过程；
大小为(2k+1)×(2k+1)的高斯模板U(x，y)为：
其中，k为正整数且k≥1，代表高斯模板的大小；x、y分别代表高斯模板中元素的横坐标和纵坐标；σ＞0，为第一灰度图像的平滑程度参数；
高斯滤波的具体操作是：用一个由式(2)计算出的高斯模板(或称模板、卷积、掩模)按照从左至右、从上到下的顺序依次扫描第一灰度图像中的每一个像素，用高斯模板确定的邻域内像素的加权平均灰度值去替代高斯模板中心像素点的灰度值；
第一灰度图像中，除左、上、右、下四条宽度为一像素的边缘外任一像素点P点的8邻像素，即紧邻像素点P的左上、上、右上、左、右、左下、下、右下8个方向上的像素；
令k＝1，根据式(2)输入阈值σ计算得出大小为3×3的高斯模板U(x，y)，该高斯模板U(x，y)中9个元素的模板值分别为m1、m2、m3、m4、m5、m6、m7、m8、m9；如矩阵1所示的一大小为3×3的高斯模板：
矩阵1中9个元素的模板值分别为1、2、1、2、4、2、1、2、1；
在对第一灰度图像做滤波时，高斯模板的中心点(即值为4的位置)对应的是第一灰度图像中的像素点P；
第一灰度图像中，除左、上、右、下四条宽度为一像素的边缘外任一像素点P点的8邻像素，即紧邻像素点P的左上、上、右上、左、右、左下、下、右下8个方向上的像素，见图4，像素点P与8邻像素的灰度值按照第一灰度图像中从左至右、从上到下的顺序依次为g1、g2、g3、g4、g5、g6、g7、g8、g9；像素点P位于g5位置，g1、g2、g3、g4、g6、g7、g8、g9位置分别为P点的8邻像素；
用该高斯模板U(x，y)按照从左到右、从上到下的顺序依次扫描第一灰度图像中除左、上、右、下四条宽度为一像素的边缘外的每一个像素P，各像素点P由高斯模板U(x，y)做高斯滤波的方法见式(3)；
所述的U0表示以高斯模板U(m，n)做高斯滤波后所到的像素点P的灰度值；S为中间变量，表示像素点P和其8邻像素的灰度值的总和；
按照式(3)依从左到右、从上到下的顺序对第一灰度图像中像素点P进行依次扫描时，对当前所扫描的像素点P由高斯模板U(x，y)所确定的邻域内所有像素(9个像素)的加权平均灰度值替代当前像素点P的灰度值，扫描结束后得到第二灰度图像；
高斯滤波相当于对图像做平滑操作，使单独的一个像素噪声在经过高斯滤波平滑处理后的图像上变得几乎没有影响；
步骤3)、使用基于Otsu的自适应Canny算法提取第二灰度图像中的边缘，如图2所示流程图，主要可分为四个步骤：
31)：用一阶偏导的有限差分来计算梯度的幅值和方向；
采用Canny算法,其中所采用的卷积模板表示为式(4)
其中M1、M2分别为第二灰度图像x、y轴方向上的一阶偏导数矩阵,依照所述一阶偏导数矩阵M1、M2按照从左到右、从上到下的顺序依次扫描第二灰度图像中的每一个像素点，得到每个像素点x、y方向上的梯度幅值Gx、Gy：
Gx＝[f(x+1,y)-f(x,y)+f(x+1,y+1)-f(x,y+1)]/2 式(5)
Gy＝[f(x,y)-f(x,y+1)+f(x+1,y)-f(x+1,y+1)]/2 式(6)
x、y表示所处理像素点的横坐标、纵坐标值，f(x,y)表示第二灰度图像中坐标为(x,y)的像素点的灰度值，该像素点的梯度幅值G(x,y)和方向角θ(x,y)由式(7)和式(8)处理得到；
由式(7)和式(8)对第二灰度图像中每一个像素点按照从左到右、从上到下的顺序依次扫描，得到第二灰度图像中每一个像素点梯度幅值和方向角；
32)：非极大值抑制；
经过上述步骤31)的处理，得到了每个像素的G(x,y)和方向角θ(x,y)，由方向角θ(x,y)即可得到梯度方向α,梯度方向α如图4所示，方向角θ(x,y)表示一个角度、度数，而梯度方向α是一个向量，起点为g5终点为Gp或Gn,Gp表示梯度方向上前一个像素的梯度幅值，Gn表示梯度方向上后一个像素的梯度幅值；
为了得到最真实的边缘，须保留每个像素点在其梯度方向上的极大值，而删掉当前像素点梯度方向上的其他值，对第二灰度图像中每个像素点按照从左到右、从上到下的顺序进行非极大值抑制，即处理得到第三灰度图像，第三灰度图像显示了遥感影像中建筑物影像的边缘，非极大值抑制的具体处理方法如下；
判断是否是极大值，需要根据当前扫描像素点梯度方向所在的区间，在x和y轴方向分别进行插值计算，然后再进行比较,以找到最大值，梯度方向所在的区间指的是梯度方向α所在的角度范围，插值是为了更合理地判断当前像素点的梯度方向：
①当以及时，令Gp＝G3w+G2(1-w),Gn＝G7w+G8(1-w)；所述的w为中间变量；
②当以及时，令Gp＝G7w+G4(1-w),Gn＝G3w+G6(1-w)；
③当以及时，令Gp＝G6w+G9(1-w),Gn＝G4w+G1(1-w)；
④当以及时，令Gp＝G2w+G1(1-w),Gn＝G8w+G9(1-w)；
其中，G5为当前处理像素的梯度幅值，Gk(1≤k≤4,6≤k≤9)表示8邻像素的梯度幅值，Gp表示梯度方向上前一个像素的梯度幅值，Gn表示梯度方向上后一个像素的梯度幅值；
当梯度方向α的值在上述相应的区间内时，当且仅当目前所处理像素的梯度幅值G5大于前述①～④中各插值方法计算出的Gp和Gn时，G5为极大值，也就表明当前处理像素是边缘点；
33)、计算最优高低阈值
本步骤采用了Canny算法中一种叫双阈值的技术，即设定一个高阈值和低阈值,图像中的像素点的梯度值如果大于高阈值则认为必然是边缘(真边缘)，如果小于低阈值，则认为必然不是边缘点；处于低阈值和高阈值之间的像素点是潜在的边缘点(假边缘)，需要在后面进行轮廓跟踪时进行取舍；
计算最优高低阈值:
读取第三灰度图像的宽度为W像素、高度为H像素，则第三灰度图像中的像素总数N为：
N＝W×H 式(9)
设第三灰度图像中每一像素的梯度幅值G(x,y)的范围为[0,G]，读取每一像素的梯度幅值G(x,y)，找梯度幅值等于确定值为i的像素，得到梯度幅值等于i的像素个数为Ni个，i∈[0,G]；设T1为最优低阈值，T2为最优高阈值；
非边缘点总个数N1(T1,T2)为：
其中，Ni是梯度幅值等于i的像素的个数；
非边缘点所占比例值ω1(T1,T2)为：
潜在边缘点总个数N2(T1,T2)为：
潜在边缘点所占比例值ω2(T1,T2)为：
真边缘点总个数N3(T1,T2)为：
真边缘点所占比例值ω3(T1,T2)为：
非边缘点的平均梯度μ1(T1,T2)、潜在边缘点的平均梯度μ2(T1,T2)和真边缘点的平均梯度μ3(T1,T2)如下：
第三灰度图像总的梯度均值μs为：
μs＝μ1(T1,T2)ω1(T1,T2)+μ2(T1,T2)ω2(T1,T2)+μ3(T1,T2)ω3(T1,T2) 式(17)
非边缘点、潜在边缘点和真边缘点这三类像素的类间方差σ2(T1,T2)为：
σ2(T1,T2)＝ω1(T1,T2)(μ1(T1,T2)-μS)2+ω2(T1,T2)(μ2(T1,T2)-μS)2+ω3(T1,T2)(μ3(T1,T2)-μS)2
式(18)
使得σ2(T1,T2)的输入参数T1、T2在总范围[1,G-1]依次取值,当类间方差σ2(T1,T2)为最大值时，此时的T1、T2的值即为最优低阈值T1和最优高阈值T2；
34)：用双阈值算法检测和边缘连接；
根据上一步33)中得到的最优低阈值T1和最优高阈值T2对第三灰度图像图像进行过滤，像素点的梯度幅值如果大于高阈值则认为当前像素点必然是边缘点，即真边缘点，如果当前处理像素点的梯度幅值小于低阈值，则认为当前像素点必然不是边缘点；处于低阈值和高阈值之间的像素点是潜在的边缘点，即假边缘；将第三灰度图像根据高阈值过滤得到的边缘点链接成轮廓，当到达轮廓的端点时，在端点的8邻点寻找潜在的边缘点，再根据所找到的潜在的边缘点,在其8邻点中循环处理收集新的边缘点，如此处理整个第三灰度图像，直至边缘闭合，得到第四灰度图像；
步骤4)、基于Freeman链码的直线提取方法，见图3所示流程图：
设单条直线L与Y轴的夹角θ范围为0～π,如图5，直线L1与Y轴的夹角为θ1,直线L2与Y轴的夹角为θ2；
按夹角θ的范围划分共有八种情况的直线L：
4)-①θ＝0或θ＝π；见图6；
4)-②见图7；
4)-③见图8；
4)-④见图9；
4)-⑤见图10；
4)-⑥见图11；
4)-⑦见图12；
4)-⑧见图13；
此外还有两条直线有交叉的情况，见图14；
设po表示左右均无目标像素的单像素；(如：p3(x,y)点像素，如果p1(x-1,y)和p2(x+1,y)均不是边缘点,p3点为po)
ps表示连续像素个数≥2的像素组；如：p3(x,y)像素是边缘点，如果p1(x-1,y)或p2(x+1,y)均是边缘点,p3点为ps；
整型常量R表示目标之间可以重叠的像素数；目标是指po或ps，例如：第一行的第2、3、4个像素是边缘点，第二行的4、5、6也是边缘点，那么从x方向上来说，重叠数就是1，就是R＝1，重叠的是第一行和第二行的第4个像素；
整型常量L表示一条直线内不包含重叠像素的最少像素数；
整型变量W(P0)表示单像素po在行或列内的位置，行内顺序为从左至右，列内顺序为从上到下；例如某行的第7个像素是po，那么W(P0)的值就是6，数组标号均从0开始，下同；
整型变量S(Ps)表示ps在行或列内的起始位置，行内顺序为从左至右，列内顺序为从上到下；例如某行的第7、8、9个像素表示一个ps，那么S(Ps)的值就是6；
整型变量E(Ps)表示ps在行或列内的终止位置；例如某行的第7、8、9个像素表示一个ps，那么E(Ps)的值就是8；
数组Ax用于存储x轴方向上的所有po和ps信息，结构包含：目标类型、单像素位置，目标是指po或ps；例如：Ax的一条记录为0,3,5:表示一个单像素，位置在第4行的第6列；如，Ax另一条记录为1,7,2,7:表示一个像素组，位置在第8行的第3列～第8列；
数组Ay用于存储y轴方向上的所有po和ps信息，结构包含：目标类型、单像素位置、像素组在列内的起始点位置，目标类型是指po或ps；例如：Ay的一条记录为0,3,5:表示一个单像素，位置在第4列的第6行；如：Ay另一条记录为1,7,2,7:表示一个像素组，位置在第8列的第3行～第8行；
整型变量R(Po)、R(Ps)分别表示单像素po、像素组ps在数组Ax内的行号；(参考Ax结构描述)
整型变量C(Po)、C(Ps)分别表示po、ps在Ay内的列号；(参考Ay结构描述)
数组AL表示存储一条直线内的所有po和ps信息；(数组AL结构与数组Ax、Ax相同)
数组GL用于以(θ,ρ)的形式存储检索到的所有直线信息；
第四灰度图像中的直线具备以下特点：
4)-(一)直线都是由po和ps组成；
4)-(二)θ角越接近或其线条中的po和ps越多；θ角等于或时，直线仅由po组成；θ角越接近0、或π，直线中po和ps越少，但ps内的像素个数越多；θ角等于0、或π时，直线仅由ps组成；
4)-(三)交叉的直线，共用其中一个po或ps；
4)-(四)直线都是由po和ps同时沿x轴、y轴两个方向连接形成；x轴方向上直线的长度逐渐增加，y轴方向上逐渐产生一定的斜率，或者，y轴方向上直线的长度逐渐增加，x轴方向上逐渐产生一定的斜率；见图15所示的直线形成时像素连接方式示意图；
4)-(五)当θ角为0、或π时，直线仅包含一个ps；
令R＝1，以检索的直线，第四灰度图像中的直线提取通过下述步骤实现：
4)-⑴从第四灰度图像的左上角像素开始，按照逐行、逐列的顺序，即从左到右、从上到下的顺序，依次扫描第四灰度图像中的每一个像素点，来检索第四灰度图像中的目标点，所述目标点即边缘点；提取到的结果保存在Ax和Ay中；
4)-⑵取出Ax中的一个目标At，将At插入数组AL，其中，第一次抽取的目标At是随机抽取，以后抽取依照下步即第(3)步所示，将抽取的目标At按照抽取的先后顺序插入数组AL；
4)-⑶沿x轴负方向，查找At的下一个目标At+1，At+1必须同时符合以下条件：
4)-⑶-①R(At+1)＝R(At)+1；
4)-⑶-②At为po时：
4)-⑶-②-Ⅰ：如果At+1为单像素，W(At+1)＝W(At)-R；
4)-⑶-②-Ⅱ：如果At+1为像素组，W(At)-R≤E(At+1)≤W(At)；
图16为At为单像素时与At+1的位置关系示意图；
4)-⑶-③At为像素组时：
4)-⑶-③-Ⅰ：如果At+1为单像素，W(At+1)＝S(At)-R；
4)-⑶-③-Ⅱ：如果At+1为像素组，S(At)-R≤E(At+1)≤S(At)；
图17为At为像素组时与At+1的位置关系示意图；
4)-⑶-④At+1未访问过；
如果找到At+1，将At+1标记为已访问且将其插入AL；令At＝At+1，递归执行该步骤4)-⑶的处理过程；如果找不到At+1，结束查找；
4)-⑷令数组AL中的第一个元素为AL0，最后一个元素为ALn：
如果数组AL中只有一个元素时：若该元素长度大于等于L，元素长度即像素数，此条直线(此条直线为AL中该元素(θ,ρ)所表示的一条直线)起点为S(AL0)，终点为E(AL0)；若长度小于L,此次检索失败；
如果数组AL中不少于2个元素，此条直线起点为E(AL0)，终点为S(ALn)；
见图18为数组AL中直线的起点和终点示意图；
4)-⑸沿x轴正方向，查找At的下一个目标At+1，At+1必须同时符合以下条件：
4)-⑸-①R(At+1)＝R(At)+1；
4)-⑸-②At为单像素时：
4)-⑸-②-Ⅰ：如果At+1为单像素，W(At+1)＝W(At)+R；
4)-⑸-②-Ⅱ：如果At+1为像素组，W(At)≤S(At+1)≤W(At)+R；
4)-⑸-③At为像素组时：
4)-⑸-③-Ⅰ：如果At+1为单像素，W(At+1)＝E(At)+R；
4)-⑸-③-Ⅱ：如果At+1为像素组，E(At)≤S(At+1)≤E(At)+R；
4)-⑸-④At+1未访问过；
如果找到At+1，将At+1标记为已访问且将其插入AL；令At＝At+1，递归执行该步骤4)-⑸的处理过程；如果找不到At+1，结束查找，执行下一步4)-⑹；
4)-⑹令数组AL中的第一个元素为AL0，最后一个元素为ALn；最终提取的直线以下述三种情况存储：
4)-⑹-①如果数组AL中只有一个元素且其长度(像素数)大于等于L，此条直线起点为S(AL0)，终点为E(AL0)；
4)-⑹-②如果数组AL中的元素数量≥2且E(AL0)-S(AL0)≥L，则此条直线起点为S(AL0)，终点为E(ALn)；
4)-⑹-③如果数组AL中的元素数量≥2且E(AL0)-S(AL0)＜L，放弃此条直线；
4)-⑺将检索到的直线插入到数组GL中,清空数组AL；
4)-⑻抽取数组Ax中一个未访问过的目标，令其等于At,递归执行上述对第四灰度图像中的直线提取步骤中的第⑶步，直至Ax中所有目标都被标记为已访问；
4)-⑼扫描数组Ay中的目标，检索和内的直线；取出数组Ay中的一个目标At，将At插入数组AL；
4)-⑽沿y轴正方向，查找At的下一个目标At+1，At+1必须同时符合以下条件：
4)-⑽-①C(At+1)＝C(At)+1；
4)-⑽-②At为po时：
4)-⑽-②-Ⅰ：如果At+1为单像素，W(At+1)＝W(At)-R；
4)-⑽-②-Ⅱ：如果At+1为像素组，W(At)+R≤E(At+1)≤W(At)；
4)-⑽-③At为像素组时：
4)-⑽-③-Ⅰ：如果At+1为单像素，W(At+1)＝S(At)-R；
4)-⑽-③-Ⅱ：如果At+1为像素组，S(At)+R≤E(At+1)≤S(At)；
4)-⑽-④At+1未访问过；
如果找到At+1，将At+1标记为已访问且将其插入数组AL；令At＝At+1，递归执行该步4)-⑽；如果找不到At+1，结束查找,执行下一步4)-⑾；
4)-⑾令数组AL中的第一个元素为AL0，最后一个元素为ALn，如果数组AL中只有一个元素且长度小于等于L，元素长度为像素数，此条直线起点为S(AL0)，终点为E(AL0)；否则，即如果数组AL中有至少2个元素，则此条直线起点为E(AL0)，终点为S(ALn)；
4)-⑿沿y轴负方向，查找At的下一个目标At+1，At+1必须同时符合以下条件：
4)-⑿-①C(At+1)＝C(At)+1；
4)-⑿-②At为po时：
4)-⑿-②-Ⅰ：如果At+1为单像素，W(At+1)＝W(At)+R；
4)-⑿-②-Ⅱ：如果At+1为像素组，W(At)-R≤S(At+1)≤W(At)；
4)-⑿-③At为像素组时：
4)-⑿-③-Ⅰ：如果At+1为单像素，W(At+1)＝S(At)+R；
4)-⑿-③-Ⅱ：如果At+1为像素组，E(At)-R≤S(At+1)≤E(At)；
4)-⑿-④At+1未访问过；
如果找到At+1，将At+1标记为已访问且将其插入数组AL；令At＝At+1，递归执行该步4)-⑿；如果找不到At+1，结束查找，执行下一步4)-⒀；
4)-⒀令数组AL中的第一个元素为AL0，最后一个元素为ALn；最终提取的直线以下述三种情况存储：
4)-⒀-①如果AL中只有一个元素且其长度大于等于L，元素长度为像素数，此条直线起点为S(AL0)，终点为E(AL0)；
4)-⒀-②如果AL中的元素数量≥2且E(AL0)-S(AL0)≥L，则此条直线起点为S(AL0)，终点为E(ALn)；
4)-⒀-③如果AL中的元素数量≥2且E(AL0)-S(AL0)＜L，放弃此条直线；
4)-⒁将数组AL中的内容以Hough变换方法处理求出θ角和ρ值，插入到数组GL中,清空数组AL；
4)-⒂抽取数组Ay中一个未访问过的目标，令其等于At,递归执行上述第4)-⑽步，直至数组Ay中所有目标都被标记为已访问；
至此，直线的提取全部完成，数组GL中存储了第四灰度图像中0≤θ≤π的所有直线；
步骤5)、直线的归并；
设阈值LR、θR、ρR分别表示两条直线的端点、θ角、ρ值的容差，令直线L1的起点、终点、θ角和ρ值分别为S1、E1、θ1和ρ1，直线L2的起点、终点、θ角、ρ值分别为S2、E2、θ2和ρ2；从第四灰度图像中提取到的直线会有部分具备以下特点之一：
5)-⑴θ1＝θ2，ρ1＝ρ2，MIN{|S1-E2|,|S2-E1|}≤LR；
5)-⑵|θ1-θ2|≤θR，|ρ1-ρ2|≤ρR；
5)-⑶ρ1＝ρ2，|θ1-θ2|≤θR；
5)-⑷θ1＝θ2，|ρ1-ρ2|≤ρR；
这四类直线需要进行归并，以减少重复的直线；
设合并后的直线为L3，其起点、终点、θ角和ρ值为S3、E3、θ3和ρ3，针对以上四种情况，直线归并如下处理：
5)-①当MIN{|S1-E2|,|S2-E1|}＝|S1-E2|时，取E3＝E1，S3＝S2；否则，当MIN{|S1-E2|,|S2-E1|}≠|S1-E2|时，E3＝E2，S3＝S1；
5)-②当ρ1＝ρ2，|θ1-θ2|≤θR时，取ρ3＝ρ1＝ρ2,
5)-③当L1‖L2时，取θ3＝θ1＝θ2,符号“‖”表示直线间的平行；
四类直线归并前，两条直线分别表示为L1(S1,E1)和L2(S2,E2)，归并后的直线表示为L3(S3,E3)；
步骤6)、图像中的角点提取
图像中的角点提取共分为5个步骤；
(61)：角度阈值过滤
定义角度阈值Tθ，直线的交叉角度在之间时，认为是直角，直线的交点即为建筑物的角点；
(62)：间断直线过滤
在背景复杂且噪声较多的遥感影像中，部分边缘直线不会被完整的获取到，出现间断的线段；需要设置阈值TL，确定在间断长度为阈值TL以内时，认为此角点可用，否则，间断的直线不做为直角边使用；
图19为间断直线的交点三种情况示意图，图19中虚线部分为间断部分，图19中仅标出了一条边间断的情况，此类间断可能在所有直角边上；
(63)：角分裂；
需要根据直线交叉的情况不同，分裂成多个直角，见图20的三种直线交叉情况的角分裂示意图；
(64)：计算角点；
已知两条直线L1:A1x+B1y+C1＝0；L2:A2x+B2y+C2＝0，两直线交点是坐标为(x,y)的像素点；
通过式(19)求出两直线的交点即角点的坐标x值、y值；
为了便于后续的角点归并和排序，在记录角点时，还需要考虑到角的开口方向；在得到角点的坐标值后，令：角的两条边为e1和e2，以角点为起点、远离角点的端点为终点的两个向量分别为和角方向如下表示：
的方向即为角的方向，所述角以θr表示，角的范围为：0≤θr≤2π；
见图21角方向示意图；
(65)：角点过滤和归并；
由于遥感影像中的复杂背景，提取目标周围会出现一些伪边缘，这些边缘主要来自影像中的植被、河流等自然地物，或者来自其它建筑物的阴影，在仅依靠灰度变化来获取边缘的前提下，这些边缘是无法与建筑物边缘区分开来的，这些边缘直线构成的角点也会被认为是建筑物的角点，如此，便需要进行角点过滤和归并，见图22为最典型的有3类角点归并情况，；
如果两角点的开口方向相同(见图22的Ⅰ、Ⅱ两种情况)：
令角点P1的两边长和夹角θr分别为E1、E2、θr1，角点P2的两边长和夹角分别为E3、E4、θr2，设置距离差阈值TD、角度容差Tθ；当角点P1和角点P2的距离D1≤TD且|θr1-θr2|≤Tθ时，对两角点进行归并；如果E1E2≥E3E4，则删除角点P2；如果E1E2＜E3E4，删除角点P1；
如果两角点的开口方向相对，且两角点之间出现了第三个角点P3：
令角点P3的两边长、夹角分别为E5、E6、θr3，当角点P1和角点P3的距离D2＞TD，且|θr1-θr3|≤Tθ时，以角点P2的θr2方向为检索方向，在θr2方向上寻找θr2-Tθ≤θr≤θr2+Tθ的角点，找到角点P1和角点P3后，如果E1E2≥E5E6，则删除角点P3；否则，当E1E2小于E5E6时，删除角点P1；
步骤7)、归属同一建筑物的角点匹配
以凸型建筑物为例，其角点分布及角方向见图23；
建筑物角点的检索可以使用如下步骤完成：
(71)：从角点集A中取一个角点As,如果As已访问过，则重新取；角点集A为经过角点过滤和归并后得到的角点集合；
(72)：沿角点As的一条边的方向寻找下一个与角点As匹配的角点As+1，匹配的条件是：
角点As+1的一条边与的方向相反，且误差在设定的角度容差θE以内；
(73)：将当前角点As标记为已访问，并将其保存在数组B内；
(74)：令As＝As+1，e1＝e4，递归执行步骤(72),直到匹配到最初的As角点为止；
(75)：遍历角点集A中所有角点，直到所有角点均为已访问；匹配结束；
对于未完成循环匹配的一组角点，由如下步骤处理推断出其余角点，构成由至少4个角点组成的矩形建筑：
⑴单独角点：由角点As在角点集A中匹配不到As+1时，以角点As的两条边以及角方向的三个端点，构建一个四边形的建筑物轮廓；
⑵两个角点：角点As在角点集A中匹配到As+1，而As+1匹配不到其它角点时，以As和As+1的两条边的端点为另两个角点，构建一个四边形的建筑物轮廓；
⑶多于两个角点：匹配到两个以上角点，但没有形成环路时，以最初的角点As的边与最后一个角点的边的交点做为最后一个角点；
步骤8)、扫描结束，输出处理后的图像，结束程序。
图24为样本图像经过本发明方法步骤2)高斯滤波处理后示意图；
图25为样本图像经过本发明方法步骤3)边缘提取处理后的示意图；
图26为样本图像经过本发明方法步骤6)角点提取处理后的示意图。
经过实测，一幅由无人机航拍的5784*4985大小、比例尺1:50000的高分辨率遥感影像,处理过程中不需要人工干预(不需要输入任何阈值)，整个处理过程在5-7秒内可完成，得到二值图像格式的处理结果。
综上可见，本发明提供的高分辨率遥感影像中的建筑物快速提取方法，可实现快速提取遥感影像中的建筑物，处理过程中不需人为干预，提取速度快，精度高；
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内；因此，本发明的保护范围应所述以权利要求的保护范围为准。

Claims (1)

1.一种高分辨率遥感影像中的建筑物快速提取方法，其特征在于，包括如下步骤：
步骤1)、读取遥感影像，提取灰度图像：
遥感图像中每个像素点的色彩具有R、G、B三个分量,P(r,g,b)代表一像素点,其中r代表该像素点的R分量的值，g代表该像素点的G分量的值，b代表该像素点的B分量的值；
灰度函数H(P)的值表示像素点P(r,g,b)的灰度值：
H(P)＝0.3r+0.59g+0.11b 式(1)
用式(1)按照从左到右、从上到下的顺序依次扫描遥感图像中的每个像素点，得到每个像素点的灰度值，按照式(1)扫描时同样的顺序组成第一灰度图像；
步骤2)、使用高斯滤波对第一灰度图像进行降噪处理：
2a)、计算大小为(2k+1)×(2k+1)的高斯模板U(x，y)：
其中，k为正整数且k≥1，x，y分别代表高斯模板中元素的横坐标和纵坐标；σ＞0，为第一灰度图像的平滑程度参数；
2b)、取k＝1，大小为3×3的高斯模板U(x，y)中9个元素的模板值分别为m1、m2、m3、m4、m5、m6、m7、m8、m9；
第一灰度图像中，除左、上、右、下四条宽度为一像素的边缘外任一像素点P点的8邻像素，即紧邻像素点P的左上、上、右上、左、右、左下、下、右下8个方向上的像素；像素点P与8邻像素的灰度值按照第一灰度图像中从左至右、从上到下的顺序依次为g1、g2、g3、g4、g5、g6、g7、g8、g9；
用该高斯模板U(x，y)按照从左到右、从上到下的顺序依次扫描第一灰度图像中除左、上、右、下四条宽度为一像素的边缘外的每一个像素P，各像素点P由高斯模板U(x，y)做高斯滤波的方法见式(3)；
所述的U0表示以高斯模板U(m，n)做高斯滤波后所到的像素点P的灰度值；S为中间变量，表示像素点P和其8邻像素的灰度值的总和；
按照式(3)依从左到右、从上到下的顺序对第一灰度图像中像素点P进行依次扫描时，对当前所扫描的像素点P由高斯模板U(x，y)所确定的邻域内所有像素的加权平均灰度值替代当前像素点P的灰度值，扫描结束后得到第二灰度图像；
步骤3)、使用基于Otsu方法的自适应Canny算法提取图像中的边缘；
31)：用一阶偏导的有限差分来计算梯度的幅值和方向；
采用Canny算法,其中所采用的卷积模板为式(4)；
其中M1、M2分别为第二灰度图像x、y轴方向上的一阶偏导数矩阵,依照所述一阶偏导数矩阵M1、M2按照从左到右、从上到下的顺序依次扫描第二灰度图像中的每一个像素点，得到每个像素点x、y方向上的梯度幅值Gx、Gy：
Gx＝[f(x+1,y)-f(x,y)+f(x+1,y+1)-f(x,y+1)]/2 式(5)
Gy＝[f(x,y)-f(x,y+1)+f(x+1,y)-f(x+1,y+1)]/2 式(6)
x、y表示所处理像素点的横坐标、纵坐标值，f(x,y)表示第二灰度图像中坐标为(x,y)的像素点的灰度值，该像素点的梯度幅值G(x,y)和方向角θ(x,y)由式(7)和式(8)处理得到；
由式(7)和式(8)对第二灰度图像中每一个像素点按照从左到右、从上到下的顺序依次扫描，得到第二灰度图像中每一个像素点梯度幅值和方向角；
32)：非极大值抑制；
经过上述步骤31)的处理，得到了每个像素的梯度幅值G(x,y)和方向角θ(x,y)，为了得到最真实的边缘，须保留每个像素点在其梯度方向上的极大值，而删掉当前像素点梯度方向上的其他值，即非极大值抑制；对第二灰度图像中每个像素点按照从左到右、从上到下的顺序进行非极大值抑制，即处理得到第三灰度图像，第三灰度图像显示了遥感影像中建筑物影像的边缘；
非极大值抑制的具体处理方法如下；
判断是否是极大值，需要根据当前扫描像素点梯度方向所在的区间，在x轴和y轴方向分别进行插值计算，然后再进行比较,以找到最大值，梯度方向所在的区间指的是梯度方向α所在的角度范围，插值是为了更合理地判断当前像素点的梯度方向：
①当以及时，令Gp＝G3w+G2(1-w),Gn＝G7w+G8(1-w)；所述的w为中间变量；
②当以及时，令Gp＝G7w+G4(1-w),Gn＝G3w+G6(1-w)；
③当以及时，令Gp＝G6w+G9(1-w),Gn＝G4w+G1(1-w)；
④当以及时，令Gp＝G2w+G1(1-w),Gn＝G8w+G9(1-w)；
其中，G5为当前处理像素的梯度幅值，Gk(1≤k≤4,6≤k≤9)表示8邻像素的梯度幅值，Gp表示梯度方向上前一个像素的梯度幅值，Gn表示梯度方向上后一个像素的梯度幅值；
当梯度方向α的值在上述相应的区间内时，当且仅当目前所处理像素的梯度幅值G5大于前述①～④中各插值方法计算出的Gp和Gn时，G5为极大值，也就表明当前处理像素是边缘点；
33)、计算最优高低阈值
读取第三灰度图像的宽度为W像素、高度为H像素，则第三灰度图像中的像素总数N为：
N＝W×H 式(9)
设第三灰度图像中每一像素的梯度幅值G(x,y)的范围为[0,G]，读取每一像素的梯度幅值G(x,y)，找梯度幅值等于确定值为i的像素，得到梯度幅值等于i的像素个数为Ni个，i∈[0,G]；设T1为最优低阈值，T2为最优高阈值；
非边缘点总个数N1(T1,T2)为：
其中，Ni是梯度幅值等于i的像素的个数；
非边缘点所占比例值ω1(T1,T2)为：
潜在边缘点总个数N2(T1,T2)为：
潜在边缘点所占比例值ω2(T1,T2)为：
真边缘点总个数N3(T1,T2)为：
真边缘点所占比例值ω3(T1,T2)为：
非边缘点的平均梯度μ1(T1,T2)、潜在边缘点的平均梯度μ2(T1,T2)和真边缘点的平均梯度μ3(T1,T2)如下：
第三灰度图像总的梯度均值μs为：
μs＝μ1(T1,T2)ω1(T1,T2)+μ2(T1,T2)ω2(T1,T2)+μ3(T1,T2)ω3(T1,T2) 式(17)
非边缘点、潜在边缘点和真边缘点这三类像素的类间方差σ2(T1,T2)为：
σ2(T1,T2)＝ω1(T1,T2)(μ1(T1,T2)-μS)2+ω2(T1,T2)(μ2(T1,T2)-μS)2+ω3(T1,T2)(μ3(T1,T2)-μS)2
式(18)
使得σ2(T1,T2)的输入参数T1、T2在总范围[1,G-1]依次取值，当类间方差σ2(T1,T2)为最大值时，此时的T1、T2的值即为最优低阈值T1和最优高阈值T2；
34)：用双阈值算法检测和边缘连接；
根据上一步33)中得到的最优低阈值T1和最优高阈值T2对第三灰度图像进行过滤，像素点的梯度幅值如果大于高阈值则认为当前像素点必然是边缘点，即真边缘点，如果当前处理像素点的梯度幅值小于低阈值，则认为当前像素点必然不是边缘点；处于低阈值和高阈值之间的像素点是潜在的边缘点，即假边缘；将第三灰度图像根据高阈值过滤得到的边缘点链接成轮廓，当到达轮廓的端点时，在端点的8邻点寻找潜在的边缘点，再根据所找到的潜在的边缘点,在其8邻点中循环处理收集新的边缘点，如此处理整个第三灰度图像，直至边缘闭合，得到第四灰度图像；
步骤4)、基于Freeman链码的直线提取方法：
单条直线L与Y轴的夹角θ范围为0～π,共有八种情况的直线L：
4)-①θ＝0或θ＝π；
4)-②
4)-③
4)-④
4)-⑤
4)-⑥
4)-⑦
4)-⑧
此外还有两条直线有交叉的情况；
设po表示左右均无目标像素的单像素；
ps表示连续像素个数≥2的像素组；
整型常量R表示目标之间可以重叠的像素数；
整型常量L表示一条直线内不包含重叠像素的最少像素数；
整型变量W(P0)表示单像素po在行或列内的位置，行内顺序为从左至右，列内顺序为从上到下；
整型变量S(Ps)表示ps在行或列内的起始位置，行内顺序为从左至右，列内顺序为从上到下；
整型变量E(Ps)表示ps在行或列内的终止位置；
数组Ax用于存储x轴方向上的所有po和ps信息，结构包含：目标类型、单像素位置，其中目标类型为po或ps；
数组Ay用于存储y轴方向上的所有po和ps信息，结构包含：目标类型、单像素位置、像素组在列内的起始点位置，其中目标类型为po或ps；
整型变量R(Po)、R(Ps)分别表示单像素po、像素组ps在数组Ax内的行号；
整型变量C(Po)、C(Ps)分别表示po、ps在Ay内的列号；
数组AL表示存储一条直线内的所有po和ps信息；
数组GL用于以(θ,ρ)的形式存储检索到的所有直线信息；
第四灰度图像中的直线具备以下特点：
4)-(一)：直线都是由po和ps组成；
4)-(二)：θ角越接近或其线条中的po和ps越多；θ角等于或时，直线仅由po组成；θ角越接近0、或π，直线中po和ps越少，但ps内的像素个数越多；θ角等于0、或π时，直线仅由ps组成；
4)-(三)：交叉的直线，共用其中一个po或ps；
4)-(四)：直线都是由po和ps同时沿x轴、y轴两个方向连接形成；x轴方向上直线的长度逐渐增加，y轴方向上逐渐产生一定的斜率，或者，y轴方向上直线的长度逐渐增加，x轴方向上逐渐产生一定的斜率；
4)-(五)：当θ角为0、或π时，直线仅包含一个ps；
令R＝1，以检索的直线，第四灰度图像中的直线提取通过下述步骤实现：
4)-⑴从第四灰度图像的左上角像素开始，按照逐行、逐列的顺序，即从左到右、从上到下的顺序，依次扫描第四灰度图像中的每一个像素点，来检索第四灰度图像中的目标点，目标点即边缘点；提取到的结果保存在Ax和Ay中；
4)-⑵取出Ax中的一个目标At，将At插入数组AL，其中，第一次抽取的目标At是随机抽取，以后抽取依照下步即第4)-⑶步所示，将抽取的目标At按照抽取的先后顺序插入数组AL；
4)-⑶沿x轴负方向，查找At的下一个目标At+1，At+1必须同时符合以下条件：
4)-⑶-①R(At+1)＝R(At)+1；
4)-⑶-②At为po时：
4)-⑶-②-Ⅰ：如果At+1为单像素，W(At+1)＝W(At)-R；
4)-⑶-②-Ⅱ：如果At+1为像素组，W(At)-R≤E(At+1)≤W(At)；
4)-⑶-③At为像素组时：
4)-⑶-③-Ⅰ：如果At+1为单像素，W(At+1)＝S(At)-R；
4)-⑶-③-Ⅱ：如果At+1为像素组，S(At)-R≤E(At+1)≤S(At)；
4)-⑶-④At+1未访问过；
如果找到At+1，将At+1标记为已访问且将其插入数组AL；令At＝At+1，递归执行该步骤4)-⑶处理过程；如果找不到At+1，结束查找；
4)-⑷令数组AL中的第一个元素为AL0，最后一个元素为ALn，如果数组AL中只有一个元素时：若该元素长度大于等于L，该元素长度即像素数，此条直线(此条直线为AL中该元素(θ,ρ)所表示的一条直线)起点为S(AL0)，终点为E(AL0)；若长度小于L,此次检索失败；如果数组AL中不少于2个元素，此条直线起点为E(AL0)，终点为S(ALn)；
4)-⑸沿x轴正方向，查找At的下一个目标At+1，At+1必须同时符合以下条件：
4)-⑸-①R(At+1)＝R(At)+1；
4)-⑸-②At为单像素时：
4)-⑸-②-Ⅰ：如果At+1为单像素，W(At+1)＝W(At)+R；
4)-⑸-②-Ⅱ：如果At+1为像素组，W(At)≤S(At+1)≤W(At)+R；
4)-⑸-③At为像素组时：
4)-⑸-③-Ⅰ：如果At+1为单像素，W(At+1)＝E(At)+R；
4)-⑸-③-Ⅱ：如果At+1为像素组，E(At)≤S(At+1)≤E(At)+R；
4)-⑸-④At+1未访问过；
如果找到At+1，将At+1标记为已访问且将其插入数组AL；令At＝At+1，递归执行该步骤4)-⑸的处理过程；如果找不到At+1，结束查找，执行下一步4)-⑹；
4)-⑹令数组AL中的第一个元素为AL0，最后一个元素为ALn；最终提取的直线分别以下述三种情况存储：
4)-⑹-①如果数组AL中只有一个元素且其长度大于等于L，该元素长度即为像素数，此条直线起点为S(AL0)，终点为E(AL0)；
4)-⑹-②如果数组AL中的元素数量≥2且E(AL0)-S(AL0)≥L，则此条直线起点为S(AL0)，终点为E(ALn)；
4)-⑹-③如果数组AL中的元素数量≥2且E(AL0)-S(AL0)＜L，放弃此条直线；
4)-⑺将检索到的直线插入到数组GL中,清空数组AL；
4)-⑻抽取数组Ax中一个未访问过的目标，令其等于At,递归执行上述对第四灰度图像中的直线提取步骤中的第4)-⑶步，直至数组Ax中所有目标都被标记为已访问；
4)-⑼扫描数组Ay中的目标，检索和内的直线；取出数组Ay中的一个目标At，将At插入数组AL；
4)-⑽沿y轴正方向，查找At的下一个目标At+1，At+1必须同时符合以下条件：
4)-⑽-①C(At+1)＝C(At)+1；
4)-⑽-②At为po时：
4)-⑽-②-Ⅰ：如果At+1为单像素，W(At+1)＝W(At)-R；
4)-⑽-②-Ⅱ：如果At+1为像素组，W(At)+R≤E(At+1)≤W(At)；
4)-⑽-③At为像素组时：
4)-⑽-③-Ⅰ：如果At+1为单像素，W(At+1)＝S(At)-R；
4)-⑽-③-Ⅱ：如果At+1为像素组，S(At)+R≤E(At+1)≤S(At)；
4)-⑽-④At+1未访问过；
如果找到At+1，将At+1标记为已访问且将其插入数组AL；令At＝At+1，递归执行该步4)-⑽；如果找不到At+1，结束查找，执行下一步4)-⑾；
4)-⑾令数组AL中的第一个元素为AL0，最后一个元素为ALn，如果数组AL中只有一个元素且该元素长度小于等于L，该元素长度即为像素数，则此条直线起点为S(AL0)，终点为E(AL0)；否则，如果数组AL中有至少2个元素，此条直线起点为E(AL0)，终点为S(ALn)；
4)-⑿沿y轴负方向，查找At的下一个目标At+1，At+1必须同时符合以下条件：
4)-⑿-①C(At+1)＝C(At)+1；
4)-⑿-②At为po时：
4)-⑿-②-Ⅰ：如果At+1为单像素，W(At+1)＝W(At)+R；
4)-⑿-②-Ⅱ：如果At+1为像素组，W(At)-R≤S(At+1)≤W(At)；
4)-⑿-③At为像素组时：
4)-⑿-③-Ⅰ：如果At+1为单像素，W(At+1)＝S(At)+R；
4)-⑿-③-Ⅱ：如果At+1为像素组，E(At)-R≤S(At+1)≤E(At)；
4)-⑿-④At+1未访问过；
如果找到At+1，将At+1标记为已访问且将其插入AL；令At＝At+1，递归执行该步4)-⑿；如果找不到At+1，结束查找，执行下一步4)-⒀；
4)-⒀令数组AL中的第一个元素为AL0，最后一个元素为ALn；最终提取的直线以下述三种情况存储：
4)-⒀-①如果数组AL中只有一个元素且其长度大于等于L，该元素长度即为像素数，此条直线起点为S(AL0)，终点为E(AL0)；
4)-⒀-②如果数组AL中的元素数量≥2且E(AL0)-S(AL0)≥L，则此条直线起点为S(AL0)，终点为E(ALn)；
4)-⒀-③如果AL中的元素数量≥2且E(AL0)-S(AL0)＜L，放弃此条直线；
4)-⒁将AL中的内容以Hough变换方法处理求出θ角和ρ值，插入到数组GL中,清空数组AL；
4)-⒂抽取Ay中一个未访问过的目标，令其等于At,递归执行上述第4)-⑽步，直至Ay中所有目标都被标记为已访问；
至此，直线的提取全部完成，数组GL中存储了第四灰度图像中0≤θ≤π的所有直线；
步骤5)、直线的归并
设阈值LR、θR、ρR分别表示两条直线的端点、θ角、ρ值的容差，令直线L1的起点、终点、θ角和ρ值分别为S1、E1、θ1和ρ1，直线L2的起点、终点、θ角、ρ值分别为S2、E2、θ2和ρ2；从第四灰度图像中提取到的直线会有部分具备以下特点之一：
5)-⑴θ1＝θ2，ρ1＝ρ2，MIN{|S1-E2|,|S2-E1|}≤LR；
5)-⑵|θ1-θ2|≤θR，|ρ1-ρ2|≤ρR；
5)-⑶ρ1＝ρ2，|θ1-θ2|≤θR；
5)-⑷θ1＝θ2，|ρ1-ρ2|≤ρR；
这四类直线需要进行归并，以减少重复的直线；
设合并后的直线为L3，其起点、终点、θ角和ρ值为S3、E3、θ3和ρ3，针对以上四种情况，直线归并如下处理：
5)-①当MIN{|S1-E2|,|S2-E1|}＝|S1-E2|时，取E3＝E1，S3＝S2；否则，当MIN{|S1-E2|,|S2-E1|}≠|S1-E2|时，E3＝E2，S3＝S1；
5)-②当ρ1＝ρ2，|θ1-θ2|≤θR时，取ρ3＝ρ1＝ρ2,
5)-③当L1‖L2时，取θ3＝θ1＝θ2,符号“‖”表示直线间的平行；
四类直线归并前，两条直线分别表示为L1(S1,E1)和L2(S2,E2)，归并后的直线表示为L3(S3,E3)；
步骤6)、图像中的角点提取
(61)：角度阈值过滤
定义角度阈值Tθ，判断两条直线的交叉角度在之间时，认为是直角，直线的交点即为建筑物的角点；
(62)：间断直线过滤
在背景复杂且噪声较多的遥感影像中，部分边缘直线不会被完整的获取到，出现间断的线段，设置阈值TL，确定在间断长度为阈值TL以内时，认为此角点可用，否则，间断的直线不做为直角边使用；
(63)：角分裂；
根据直线交叉的情况，分裂成多个直角；
(64)：计算角点；
已知两条直线L1:A1x+B1y+C1＝0；L2:A2x+B2y+C2＝0，两直线交点是坐标为(x,y)的像素点；
通过式(19)求出两直线的交点即角点的坐标x值、y值；
为了便于后续的角点归并和排序，在记录角点时，还需要考虑到角的开口方向；在得到角点的坐标值后，令：角的两条边为e1和e2，以角点为起点、远离角点的端点为终点的两个向量分别为和角方向如下表示：
的方向即为角的方向，所述角以θr表示，角的范围为：0≤θr≤2π；
(65)：角点过滤和归并；
由于遥感影像中的复杂背景，提取目标周围会出现一些伪边缘，需要进行角点过滤和归并：
如果两角点的开口方向相同：
令角点P1的两边长和夹角θr分别为E1、E2、θr1，角点P2的两边长和夹角分别为E3、E4、θr2，设置距离差阈值TD、角度容差Tθ；当角点P1和角点P2的距离D1≤TD且|θr1-θr2|≤Tθ时，对两角点进行归并；如果E1E2≥E3E4，则删除角点P2；如果E1E2＜E3E4，删除角点P1；
如果两角点的开口方向相对，且两角点之间出现了第三个角点P3：
令角点P3的两边长、夹角分别为E5、E6、θr3，当角点P1和角点P3的距离D2＞TD，且|θr1-θr3|≤Tθ时，以角点P2的θr2方向为检索方向，在θr2方向上寻找θr2-Tθ≤θr≤θr2+Tθ的角点，找到角点P1和角点P3后，如果E1E2≥E5E6，则删除角点P3；否则，当E1E2小于E5E6时，删除角点P1；
步骤7)、归属同一建筑物的角点匹配
建筑物角点的检索可以使用如下步骤完成：
(71)从角点集A中取一个角点As,如果As已访问过，则重新取；角点集A为经过角点过滤和归并后得到的角点集合；
(72)：沿角点As的一条边的方向寻找下一个与角点As匹配的角点As+1，匹配的条件是：
角点As+1的一条边与的方向相反，且误差在设定的角度容差θE以内；
(73)：将当前角点As标记为已访问，并将其保存在数组B内；
(74)：令As＝As+1，e1＝e4，递归执行步骤(72),直到匹配到最初的As角点为止；
(75)：遍历角点集A中所有角点，直到所有角点均为已访问；匹配结束；
对于未完成循环匹配的一组角点，由如下步骤处理推断出其余角点，构成由至少4个角点组成的矩形建筑：
⑴单独角点：由角点As在角点集A中匹配不到As+1时，以角点As的两条边以及角方向的三个端点，构建一个四边形的建筑物轮廓；
⑵两个角点：角点As在角点集A中匹配到As+1，而As+1匹配不到其它角点时，以As和As+1的两条边的端点为另两个角点，构建一个四边形的建筑物轮廓；
⑶多于两个角点：匹配到两个以上角点，但没有形成环路时，以最初的角点As的边与最后一个角点的边的交点做为最后一个角点；
步骤8)、扫描结束，输出处理后的图像，结束程序。
CN201610050661.4A 2016-01-26 2016-01-26 一种高分辨率遥感影像中的建筑物快速提取方法 Active CN105719306B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610050661.4A CN105719306B (zh) 2016-01-26 2016-01-26 一种高分辨率遥感影像中的建筑物快速提取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610050661.4A CN105719306B (zh) 2016-01-26 2016-01-26 一种高分辨率遥感影像中的建筑物快速提取方法

Publications (2)

Publication Number Publication Date
CN105719306A CN105719306A (zh) 2016-06-29
CN105719306B true CN105719306B (zh) 2018-09-11

Family

ID=56154130

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610050661.4A Active CN105719306B (zh) 2016-01-26 2016-01-26 一种高分辨率遥感影像中的建筑物快速提取方法

Country Status (1)

Country Link
CN (1) CN105719306B (zh)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108304840B (zh) * 2017-08-31 2022-11-11 腾讯科技（深圳）有限公司 一种图像数据处理方法以及装置
CN107730518B (zh) * 2017-09-30 2021-06-01 惠州华阳通用电子有限公司 一种基于fpga的canny算法阈值获取方法及装置
CN110388919B (zh) * 2019-07-30 2023-05-23 上海云扩信息科技有限公司 增强现实中基于特征图和惯性测量的三维模型定位方法
CN110569745B (zh) * 2019-08-20 2022-04-12 北方工业大学 一种遥感图像建筑区检测方法
CN111523391B (zh) * 2020-03-26 2021-11-05 上海刻羽信息科技有限公司 建筑物的识别方法、系统、电子设备及可读存储介质
CN111754536B (zh) * 2020-06-29 2024-04-16 上海商汤智能科技有限公司 图像标注方法、装置、电子设备及存储介质
CN112164142A (zh) * 2020-10-21 2021-01-01 江苏科技大学 一种基于智能手机的楼盘采光模拟方法
CN113743351B (zh) * 2021-09-14 2023-07-04 北京石油化工学院 一种基于边缘方向语义信息的遥感影像场景识别的方法
CN116310447B (zh) * 2023-05-23 2023-08-04 维璟(北京)科技有限公司 基于计算机视觉的遥感图像变化智能检测方法及系统

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101599120A (zh) * 2009-07-07 2009-12-09 华中科技大学 一种遥感影像建筑物识别方法
KR100963651B1 (ko) * 2009-03-12 2010-06-15 세종대학교산학협력단 항공라이다자료를 이용한 건물외곽선 자동 추출 방법
CN102565810A (zh) * 2011-12-30 2012-07-11 武汉大学 一种遥感影像上土地利用地物边界轮廓提取方法
CN103020342A (zh) * 2012-12-04 2013-04-03 南京大学 一种从地面LiDAR数据中提取建筑物轮廓和角点的方法
CN103218618A (zh) * 2013-01-09 2013-07-24 重庆交通大学 一种基于遥感数字图像的公路路线自动提取方法
CN104182754A (zh) * 2014-08-19 2014-12-03 山东临沂烟草有限公司 一种基于高分辨率遥感影像的农村居民点信息提取方法
CN104794723A (zh) * 2015-05-04 2015-07-22 福建师范大学 一种基于概率的遥感影像建筑物位置检测方法
CN104915672A (zh) * 2014-03-13 2015-09-16 北京大学 一种基于高分辨率遥感图像的矩形建筑物提取方法及系统

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR100963651B1 (ko) * 2009-03-12 2010-06-15 세종대학교산학협력단 항공라이다자료를 이용한 건물외곽선 자동 추출 방법
CN101599120A (zh) * 2009-07-07 2009-12-09 华中科技大学 一种遥感影像建筑物识别方法
CN102565810A (zh) * 2011-12-30 2012-07-11 武汉大学 一种遥感影像上土地利用地物边界轮廓提取方法
CN103020342A (zh) * 2012-12-04 2013-04-03 南京大学 一种从地面LiDAR数据中提取建筑物轮廓和角点的方法
CN103218618A (zh) * 2013-01-09 2013-07-24 重庆交通大学 一种基于遥感数字图像的公路路线自动提取方法
CN104915672A (zh) * 2014-03-13 2015-09-16 北京大学 一种基于高分辨率遥感图像的矩形建筑物提取方法及系统
CN104182754A (zh) * 2014-08-19 2014-12-03 山东临沂烟草有限公司 一种基于高分辨率遥感影像的农村居民点信息提取方法
CN104794723A (zh) * 2015-05-04 2015-07-22 福建师范大学 一种基于概率的遥感影像建筑物位置检测方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
基于Freeman编码的遥感影像直角型房屋半自动提取方法研究;秦永 等;《测绘科学》;20091130;第34卷(第6期);第203-205页 *
基于改进Canny算子的遥感图像边缘检测;石桂名 等;《大连交通大学学报》;20150630;第36卷(第3期);第87-90页 *
航空影像的直线特征检测算法;孙爱蓉 等;《武汉理工大学学报 交通科学与工程版》;20031231;第27卷(第6期);第808-810页 *

Also Published As

Publication number Publication date
CN105719306A (zh) 2016-06-29

Similar Documents

Publication Publication Date Title
CN105719306B (zh) 一种高分辨率遥感影像中的建筑物快速提取方法
CN111931538B (zh) 一种Micro QR二维码的定位方法
CN102930277B (zh) 一种基于识别反馈的字符图像验证码识别方法
CN107092871B (zh) 基于多尺度多特征融合的遥感影像建筑物检测方法
O'Gorman k× k thinning
CN110309687A (zh) 一种二维码图像的校正方法及校正装置
CN109858325B (zh) 一种表格检测方法和装置
CN104933434A (zh) 一种结合LBP特征提取和surf特征提取方法的图像匹配方法
CN104123554B (zh) 基于mmtd的sift图像特征提取方法
CN109559324A (zh) 一种线阵图像中的目标轮廓检测方法
CN103198319A (zh) 用于矿山井筒环境下的模糊图像角点提取方法
CN103632137A (zh) 一种人眼虹膜图像分割方法
CN109559273A (zh) 一种面向车底图像的快速拼接方法
CN110544300A (zh) 一种基于二维手绘图像特征自动生成三维模型的方法
CN109615598A (zh) 一种基于边缘绘图参数自由算法的输电线识别方法
CN111861866A (zh) 一种变电站设备巡检图像全景重建方法
CN109509151A (zh) 图像及视频拼接方法、计算机可读存储介质及计算机设备
Kumar et al. Comparative analysis for edge detection techniques
CN114092491A (zh) 一种建筑户型语义分割图矢量化方法及装置
Manandhar et al. Segmentation based building detection in high resolution satellite images
CN110569711B (zh) 面向人体动作识别方法
Shahab et al. A modified 2D chain code algorithm for object segmentation and contour tracing.
Unsalan Gradient-magnitude-based support regions in structural land use classification
CN106228553A (zh) 高分辨率遥感图像阴影检测装置与方法
CN108171771A (zh) 一种结合外部边缘信息和内部聚合笔道的线描画生成算法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant