[发明专利]区域灰度和方向成分算子的遥感影像中树冠直径提取方法

摘要

本发明公开一种区域灰度和方向成分算子的遥感影像中树冠直径提取方法，利用不断扩大的区域灰度和方向成分算子发现遥感影像一定区域范围内像元的构成、以及明暗方向的规律，并基于这一规律进行区域扩大来对遥感影像中树冠直径进行提取。本发明能够应对同一树木的像元颜色的不同，同时在能够容忍从树木中心到边沿树木形态的变化，使得提取的效率更高，根据同一树木的像元颜色的不同，同时在能够容忍从树木中心到边沿树木形态的变化，使得提取的效率更高。

主权项

1.一种区域灰度和方向成分算子的遥感影像中树冠直径提取方法，包括以下步骤：S1，输入由红、绿、蓝三个波段构成的遥感影像Image，将Image转换为256级灰度的灰度影像GrayImage，建立树木直径探测结构RStruct;输入区域变化容忍度最大阈值MaxRate、区域变化容忍度最小阈值MinRate，阈值阶梯Gap; 指定当前容忍度变量CurrentRate=MaxRate；S101，输入遥感影像Image；S102，将Image转换为256级灰度的灰度影像GrayImage；S103，对GrayImage的所有像元进行统计，获得GrayImage的所有像元的均值IMean；S104，建立树木直径探测结构RStruct，RStruct对应GrayImage上的一个位置，该结构包含以下属性；RX，该结构在GrayImage上对应的横坐标；RY，该结构在GrayImage上对应的纵坐标；RCircleVector圆区域灰度和方向成分矢量，为一个8个元素的矢量，初始值为[0,0,0,0,0,0,0,0]；RStop标志是否需要进行继续计算，值为TRUE或FALSE，初始值为FALSE；RScale圆形区域尺度，初始值为2；RRemoved删除标记，值为TRUE或FALSE，初始值为FALSE；S105,输入区域变化容忍度最大阈值MaxRate、区域变化容忍度最小阈值MinRate，阈值阶梯Gap；MaxRate的默认值为0.5,MinRate的默认值为0.1,Gap的默认值为0.01；S106，CurrentRate=MaxRate；S2，建立圆区域灰度和方向成分算子CircleOperator，该算子的输入为横坐标CX，纵坐标CY，尺度CScale，输出为圆区域灰度和方向成分矢量CircleVector；S201，建立圆区域成分算子CircleOperator，该算子的输入为横坐标CX，纵坐标CY，尺度CScale，S202, pixels=在GrayImage上取出以CX和CY为中心以CScale为半径的圆的范围内的所有像元；S203，pixelnum= pixels内像元的个数；S204，基于像元灰度值对pixels的所有像元进行K‑Means算法聚类，聚类中心点个数为4，K‑Means初始化中心点的灰度值分别为0、64、128、182，聚类的结果将pixels中的所有像元分配到4个像元分组内，分别为group1、group2、group3、group4；S205，计算如下灰度成分变量，a1=group1中像元个数/ pixelnum，a2=group2中像元个数/ pixelnum，a3=group3中像元个数/ pixelnum，a4=group4中像元个数/ pixelnum；S206，初始化方向成分数组b，该数组包含4个元素，初始值为0；S207，像元分组加和值sumvalue=0；S208，像元计数器pixelcounter=1；S209,tpixel=取出pixels中的第pixelcounter个元素，获取tpixel在GrayImage上的横坐标TCX和纵坐标TCY；S210,查找tpixel所处的像元分组，如果在group1内则tvalue=1，如果在group2内则tvalue=2,如果在group3内则tvalue=3，如果在group4内则tvalue=4；S211，建立向量v1, 该向量的起点为(CX，CY) 终点为（CX+1，CY）；建立向量v2，该向量的起点为(CX，CY)终点为（TCX，TCY）；计算以v1的方向为起始，以（CX，CY）为中心逆时针旋转到v2的方向对应的夹角angle；S212，计算位置变量position=round(angle/90)，其中round为进行四舍五入；S213, 设置方向成分数组b的第position个元素的值，b[position]= b[position]+tvalueS214,sumvalue=sumvalue+tvalue；S215,pixelcounter=pixelcounter+1；S216，如果pixelcounter>pixelnum则转到S217，否则转到S209；S218，计算方向成分数组b，对应公式为b[1]=b[1]/sumvalue, b[2]=b[2]/sumvalue, b[3]=b[3]/sumvalue, b[4]=b[4]/sumvalue；S219, 建立圆区域灰度和方向成分矢量CircleVector=[a1,a2,a3,a4,b[1],b[2],b[3],b[4]]S220, 将CircleVector作为结果返回；S3，对于GrayImage中的每一个像元，建立树木直径探测结构RStruct，放入RStruct列表RStructList中，RStructList的元素个数为RNum；S301，建立RStruct列表RStructList，RStructList初始状态下包含0个元素；S302，CTPixel=在GrayImage中取出一个像元；S303，CTRStruct=建立一个RStruct结构；S304，CTRStruct.RX=CTPixel在GrayImage中的横坐标，CTRStruct.RY=CTPixel在GrayImage中的纵坐标；S305，CTRStruct.RCircleVector=CircleOperator(CTRStruct.RX,CTRStruct.RY,CTRStruct.RScale)，利用圆区域灰度和方向成分算子计算CTRStruct对应位置的圆区域灰度和方向成分矢量；S306，将CTRStruct加入到RStructList中；S307，如果GrayImage还有项目没有处理那么放S302，否则转到S308；S308，RNum=RStructList的元素个数；S4，建立非树木与重合区域删除算子DeleteOperator，该算子的输入为一个在列表中的位置TX，判断RStructList的第TX个元素是否需要删除，如果需要则将第TX个元素的RRemoved属性标记为TRUE；S401，建立非树木与重合区域删除算子DeleteOperator，该算子的输入为一个在列表中的位置TX；S402，DORStruct=取出RStructList中第TX个元素；S403，DOPixels=在GrayImage上取出以DORStruct.RX和DORStruct.RY为中心，以DORStruct.RScale为半径的所有像元；S404，删除比对计数器DOCounter=1，初始像元数量DOPNum=统计DOPixels中的像元个数；S405，如果DORStruct与TX相等那么转到S410，否则转到S406；S406，DORStruct2=取出RStructList中第DOCounter个元素；S407，如果DORStruct2.RRemoved等于TRUE则转到S410，否则转到S408；S408，DOPixels2=在GrayImage上取出以DORStruct2.RX和DORStruct2.RY为中心，以DORStruct2.RScale为半径的所有像元；S409，基于像元位置计算DOPixels和DOPixels2的交集union，将union的所有像元从DOPixels中删除；S410，DOCounter= DOCounter+1；S411, 如果DOCounter<=RNum则转到S405，否则转到S412；S412，DOResult=统计DOPixels的元素个数/DOPNum；S413, 如果DOResult<0.75则转到S414，否则转到S415；S414，DORStruct.RRmoved=TRUE，DORStruct.RStoped=TRUE；S415, 该过程处理结束；S5，建立区域扩大算子EnlargeOperator，该算子的输入为一个在列表中的位置EX，判断RStructList的第EX个元素是否扩大其对应区域，如果不可以则将第EX个元素的RStoped属性标记为TRUE；S501, 建立区域扩大算子EnlargeOperator，该算子的输入为一个在列表中的位置EX；S502，EORstruct=取出RStructList中的第EX个元素；S503，vector1= EORstruct.RCircleVector；S504，vector2=CircleOperator(EORstruct.RX,EORstruct.RY, EORstruct.RScale+1)，计算EORstruct的RScale扩大1之后的圆区域灰度和方向成分矢量；S505，EOResult=l2norm(vector1‑vector2)/l2norm(vector1)，其中l2norm为l2范数；S506，如果EOResultRNum则转到S607，否则转到S602；S607, CurrentRate=CurrentRate‑Gap；S608, 如果CurrentRateRNum则转到S616，否则转到S611；S616，该处理过程结束。