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

权利要求书

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]；

RStoped标志是否需要进行继续计算，值为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]+tvalue

S214,sumvalue=sumvalue+tvalue；

S215,pixelcounter=pixelcounter+1；

S216，如果pixelcounterpixelnum则转到S217，否则转到S209；

S217，计算方向成分数组b，对应公式为b[1]=b[1]/sumvalue, b[2]=b[2]/sumvalue,b[3]=b[3]/sumvalue, b[4]=b[4]/sumvalue；

S218, 建立圆区域灰度和方向成分矢量CircleVector=[a1,a2,a3,a4,b[1],b[2],b[3],b[4]]

S219, 将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，如果DOCounter与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, 如果DOResult0.75则转到S414，否则转到S415；

S414，DORStruct.RRmoved=TRUE，DORStruct.RStoped=TRUE；

S415, 第S4步骤处理结束；

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，如果EOResultCurrentRate则转到S507，否则转到S508；

S507，EORstruct.RSCale= EORstruct.RSCale+1, 转到S509；

S508，EORstruct.RStoped=TRUE；

S509, 第S5步骤处理结束；

S6，利用非树木与重合区域删除算子DeleteOperator和区域扩大算子EnlargeOperator处理RStructList，输出树冠直径检测结果；

S601，处理计数器DTCounter=1；

S602,DTRstruct=取出RStructList中的第DTCounter个元素；

S603, 如果DTRstruct.RRemoved=FALSE则执行DeleteOperator(DTCounter)；

S604, 如果DTRstruct.RStoped=FALSE则执行EnlargeOperator(DTCounter)；

S605, DTCounter=DTCounter+1；

S606, 如果DTCounterRNum则转到S607，否则转到S602；

S607, CurrentRate=CurrentRate-Gap；

S608, 如果CurrentRateMinRate则执行CurrentRate=MinRate；

S609，如果RStructList包含的元素中存在RStoped属性为FALSE的元素，则转到S601；否则转到S610；

S610，DTCounter=1；

S611,DTRstruct2=取出RStructList中的第DTCounter个元素；

S612, 如果DTRstruct.RRemoved=FALSE则转到S613，否则转到S614；

S613，检测到树冠，以 (DTRstruct2.RX, DTRstruct2.RY)为中心，以DTRstruct2.RScalse为半径绘制圆形，并输出该圆形的直径为2×DTRstruct2.RScale；

S614，DTCounter= DTCounter+1；

S615，如果DTCounterRNum则转到S616，否则转到S611；

S616，第S6步骤处理结束。