[發明專利]區域灰度和方向成分算子的遙感影像中樹冠直徑提取方法

有效

• 申請號: 201810895227.5
• 申請日: 2018-08-08
• 公開（公告）號: CN109166148B
• 公開（公告）日: 2021-08-27
• 發明（設計）人: 許駿;付浩海;潘欣;張敏;張妳妳
• 申請（專利權）人: 長春工程學院
• 主分類號: G06T7/62
• 分類號: G06T7/62
• 代理公司: 吉林長春新紀元專利代理有限責任公司 22100
• 代理人: 陳宏偉

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步驟處理結束。