透過空拍圖進行稻田區域影像分割與秧苗偵測之研究

(1)

國立臺灣大學工學院工程科學及海洋工程學系碩士論文

Department of Engineering Science and Ocean Engineering College of Engineering

National Taiwan University Master Thesis

透過空拍圖進行稻田區域影像分割與秧苗偵測之研究 A Study on Paddy Field Segmentation and Rice Seedling

Detection from UAV Images

楊欣怡 HSIN-YI Yang

指導教授：黃乾綱博士

Advisor: CHIEN-KANG Huang, Ph.D.

中華民國 110 年 1 月

Jan 2021

(2)

口試委員會審定書

(3)

誌謝

研究就像走上一條充滿障礙物的道路，在過程中必須不斷的披荊斬棘，開拓新的前進方向。研究也像是身處一座迷宮森林，在尋找出口的過程中，走的路線不僅曲曲折折彎彎繞繞，常常走著走著碰到死路，就只好原路折返，再選擇下一條可能通往出口的道路，繼續摸索著前進。

感謝黃乾綱老師與劉力瑜老師的指導，使我在與老師們的無數次討論中獲益良多，學習新知的同時也受到很多啟發，並在研究過程中慢慢摸索解決問題的方法。

本研究為葉怡萱學姊論文之後續研究。感謝學姊先奠定了一定的基礎，讓我有可以參考的研究思路，也讓我對於相關文獻要找哪些領域的論文比較有概念。除此之外，

也要感謝這一路上好多好多願意幫助我的人，尤其是溫暖的實驗室同學們。在樸實無華且枯燥的研究生活中，能有一群戰友一起吃飯聊天、互相討論研究還有加油打氣，在心靈層面上給予我的支持真的很大，讓我得以堅持到現在。

(4)

摘要

近年來，由於人口逐漸高齡化，台灣農業面臨勞動力嚴重短缺的問題。為了減輕農事生產作業量，進而降低勞動力需求，運用物聯網、感測元件、大數據分析等技術的智慧農業在台灣開始被大力推廣。其中，空拍機結合影像資料分析提供農民更便捷的作物成長監控方式，在農業上應用的需求日益增加。

本研究針對早期水稻秧苗之空拍圖，提出之方法其目的指在分割影像中的稻田區域與秧苗植株偵測，除此之外，提出方法也根據偵測之秧苗，標記其所屬稻田劃分區域，並完成相鄰影像與多張影像之間的秧苗對齊。研究主要分為四大部分進

行實驗：（一）稻田區域影像分割、（二）秧苗偵測、（三）秧苗標記以及（四）秧

苗對齊。第一部分使用三種影像分割演算法將輸入影像分割為不同的子區域，再利用分類資料集將每塊子區域分類為稻田或者非稻田兩類，最後合併稻田子區域得到稻田區域遮罩影像。第二部分利用不同植被指數產生前處理影像，於第一部分得到的稻田區域進行斑點偵測，所得之偵測結果即為秧苗。其中，此部分設計一個系統化的流程，使得每種前處理方法之偵測參數能被自動決定。第三部分藉由將輸入影像轉換至鑲嵌圖座標平面，依偵測秧苗轉換後的座標決定其所屬稻田區域並進行標記。最後第四部份提出一種基於秧苗的影像拼接方法來完成相鄰影像之間的秧苗對齊，再延伸到多張影像之間的秧苗對齊。

研究成果方面，第一部分使用的三種影像分割演算法，於最後的稻田區域分割結果，Ｆ度量值皆可達到95% 以上，其中以使用 Felzenszwalb 演算法之 F 度量值 97.04% 為最佳。第二部分的秧苗偵測結果證明使用植被指數對於偵測結果有正面影響，而使用TGI 具有最佳平均偵測結果，可達成平均 F 度量值 93.84%，其中平均精確度為96.15%、平均召回率為 91.96%。至於第三部分的標記秧苗所在稻田區

域，其平均標記準確度可達99.2%。而第四部份在相鄰影像之間的秧苗對齊方面，

使用提出方法所得之秧苗對數量，與原始方法相比平均可提升 26.17%。在多張影

(5)

像之間的秧苗對齊方面，ABC 三張影像中，AC 透過 B 間接達成影像拼接的方式，

會比AC 直接進行影像拼接所得結果更加精確。

關鍵字：精準農業、空拍圖、影像分割、特徵偵測、秧苗偵測、影像拼接

(6)

ABSTRACT

Due to the aging of the population, Taiwan’s agricultural manpower has been in a large shortage in recent years. In order to reduce the labor demand of agricultural production, smart farming utilizing big data analysis, Internet of Things (IoTs), and sensors is of great interest. Particularly, combining with image data analysis, the inquiries of Unmanned Aerial Vehicles (UAVs) applied in agriculture is rising. It can provide farmers with better management on the process of growing crops.

The research proposed an approach to automatically label rice seedling from UAV images. There are 4 main parts in the proposed method: paddy field segmentation, rice seedling detection, rice seedling labeling, and rice seedling alignment. First, with the classification tool and three different segmentation algorithms: Watershed algorithm, Felzenszwalb’s algorithm and SLIC algorithm, paddy fields were extracted from the UAV images and the paddy mask image was computed. Second, by using blob detection with the paddy mask, rice seedlings were detected from different pre-processing images, which were generated by calculating vegetation indices (VIs) and image enhancement. Third, each detected rice seedling was labeled with an ID of paddy area by utilizing the homography matrix between the input UAV image and the mosaic. Last, the proposed method implemented an image stitching algorithm based on rice seedlings, and achieve

(7)

rice seedling alignment.

For paddy field segmentation, all the three image segmentation algorithm can achieve more than 95% f-measure, where Felzenszwalb’s algorithm got the best performance with 96.26% precision, 97.88% recall and 97.04% f-measure. For rice seedling detection, most of the pre-processing images using VIs obtained better result than the pre-processing image which was just converted from RGB to grayscale. Among the pre-processing images, the one using TGI gained the best detection result with 96.15%

precision, 91.96% recall and 93.84% f-measure. For rice seedling labeling, the result can achieve 99.20% accuracy. As for rice seedling alignment, the result shows that using seedling-based image stitching algorithm can increase 26.17% the numbers of rice seedling pairs between adjacent images. Besides, the result also shows that within multiple image stitching, the perspective transformation between the first image and the third image was more accurate by indirectly transforming with the second image.

Keywords: precision agriculture; UAV image; image segmentation; feature detection;

rice seedling detection; image stitching

(8)

圖目錄

圖 1 基於 NDVI 和 TGI 之辨識結果比較 [12] ... 6

圖 2 不同指數在植物像素分割方法上的優缺點比較 [17] ... 7

圖 3 Ahmed 等人提出方法之流程圖 [4] ... 9

圖 4 葉怡萱提出之秧苗標記方法流程圖 [10] ... 10

圖 5 以洪水氾濫為原則的分水嶺演算法 [40] ... 13

圖 6 k-NN 演算法 [31] ... 14

圖 7 RANSAC：內群值與離群值 [49] ... 15

圖 8 系統架構圖 ... 18

圖 9 稻田區域影像分割架構圖 ... 19

圖 10 稻田區域之特徵點偵測架構圖 ... 20

圖 11 相鄰影像之間的秧苗對齊架構圖 ... 21

圖 12 分水嶺演算法 ... 22

圖 13 分水嶺演算法使用不同閾值之分割結果 ... 23

圖 14 Felzenszwalb 演算法使用不同參數之分割結果 ... 24

圖 15 SLIC 演算法分割結果 ... 25

圖 16 利用 RAG 與不同閾值的子區域合併結果 ... 25

圖 17 使用不同演算法之影像分割結果 ... 26

圖 18 分類資料集之標記影像 ... 27

圖 19 分類資料集內的 HS 直方圖 ... 28

圖 20 待分類子區域與其 HS 直方圖 ... 29

圖 21 稻田區域影像分割之結果圖 ... 29

圖 22 使用不同方法之稻田區域影像分割結果圖 ... 30

圖 23 植被指數影像 ... 31

(12)

圖 24 調整後之植被指數影像 ... 31

圖 25 植被指數相關矩陣 ... 33

圖 26 選用之植被指數相關矩陣 ... 33

圖 27 使用 CLAHE 前後比較圖 ... 34

圖 28 前處理影像 ... 34

圖 29 影像二值化前後比較圖 ... 35

圖 30 稻田塊示意圖 ... 36

圖 31 使用不同 thresholdStep 之偵測結果示意圖 ... 36

圖 32 秧苗偵測結果示意圖 ... 37

圖 33 空間點透過相機中心的投影示意圖 [53] ... 39

圖 34 平面上一點映射至另一平面之示意圖 [53]... 40

圖 35 秧苗標記示意圖 ... 41

圖 36 秧苗之間的距離範例圖 ... 42

圖 37 對齊之秧苗對 𝑞, 𝑟 示意圖 ... 42

圖 38 SIFT 影像拼接示意圖 ... 44

圖 39 基於不同區域之影像拼接內群對 ... 46

圖 40 基於不同區域之影像拼接結果 ... 46

圖 41 對照組與實驗組之秧苗對齊結果 ... 47

圖 42 影像重新取樣示意圖 ... 48

圖 43 多張影像之間的拼接示意圖 ... 49

圖 44 多張影像之拼接結果對照圖 ... 52

圖 45 三張影像之間的重疊示意圖 ... 53

圖 46 單株秧苗資訊示意圖 ... 53

圖 47 空拍機飛行路線與高度 ... 54

(13)

圖 49 經緯航太提供之鑲嵌圖 ... 55

圖 50 稻田範圍與區域劃分圖 ... 56

圖 51 實驗影像示意圖 ... 57

圖 52 分類資料集影像示意圖 ... 57

圖 53 分類資料集之標記影像 ... 58

圖 54 分割結果與標記稻田區域之比較圖 ... 61

圖 55 稻田區域影像分割結果之盒鬚圖 ... 62

圖 56 Ws 稻田區域影像分割 ... 63

圖 57 使用 Fz 與 SL_R 之子區域數量盒鬚圖 ... 64

圖 58 各階段之執行時間 ... 64

圖 59 分割子區域中各階段之執行時間 ... 64

圖 60 為評估實驗結果之取樣區域 ... 65

圖 61 取樣區域之秧苗偵測結果 ... 66

圖 62 秧苗偵測結果之盒鬚圖 ... 66

圖 63 使用 CLAHE 之秧苗偵測結果 ... 67

圖 64 環境變化明顯區域之偵測結果 ... 69

圖 65 標記結果影像 ... 70

圖 66 稻田劃分區域中央與交界之標記結果 ... 72

圖 67 秧苗對齊結果影像一 ... 73

圖 68 秧苗對齊結果影像二 ... 74

圖 69 對照組 (A → C) 之秧苗對齊結果裁切影像... 76

圖 70 實驗組 (A → B → C) 之秧苗對齊結果裁切影像... 77

圖 71 對照組之秧苗對齊結果影像 ... 80

圖 72 實驗組之秧苗對齊結果影像 ... 81

(14)

表目錄

表 1 RGB 植被指數公式對照表 ... 8

表 2 使用範例影像產生之子區域數量 ... 26

表 3 分類資料集中的標記影像類別與張數 ... 27

表 4 各前處理方法之 thresholdStep 設定值 ... 36

表 5 秧苗資訊之儲存方式 ... 53

表 6 混淆矩陣表 ... 58

表 7 不同演算法之平均稻田區域影像分割結果 ... 62

表 8 分割所得平均子區域數量 ... 63

表 9 秧苗偵測之平均結果 ... 67

表 10 標記數字與顏色對照表 ... 70

表 11 秧苗標記實驗結果 ... 71

表 12 相鄰影像之間的秧苗對齊結果 ... 75

表 13 AB 秧苗對齊結果 ... 78

表 14 AC 秧苗對齊結果 ... 78

表 15 ABC 秧苗對齊結果 ... 79

表 16 Ws 之稻田區域影像分割結果 ... 89

表 17 Fz 之稻田區域影像分割結果 ... 90

表 18 SL_R 之稻田區域影像分割結果 ... 91

表 19 Gray 之秧苗偵測結果 ... 91

表 20 CLAHE 之秧苗偵測結果 ... 92

表 21 ExG 之秧苗偵測結果 ... 92

表 22 ExR 之秧苗偵測結果 ... 93

表 23 MExG 之秧苗偵測結果 ... 93

(15)

表 24 TGI 之秧苗偵測結果 ... 94 表 25 COM2 之秧苗偵測結果 ... 94

(16)

第一章緒論

1.1. 研究背景

臺灣為一個位於亞熱帶之海島國家，氣候適合農作物生長，但由於山多平原少等地理條件限制，使得耕地面積狹小，農戶平均耕地規模不到一公頃，屬於小農經營型態，生產成本較高。根據行政院所公布之農業經營現況，臺灣農作物種植面積為74.2 萬公頃，其中糧食作物以稻米為主，佔總作物種植面積中的 36%。而在全球皆面臨人口高齡化、勞動力短缺危機的現代，運用科技降低農業生產成本、提高產量、改良作物品質等成為農業領域研究之重點。

精準農業 (precision agriculture) 是一種使用 GPS (global positioning system) 獲得作物生長區域的詳細資訊，並依現場數據主導決策，從而管理農作物生長狀況以提高產量和作物品質的手段 [1, 2]。在精準農業中，由於無人飛行載具 (unmanned aerial vehicle, UAV) 的方便性、低使用成本與低飛行高度，常與相機組合成空拍機，

幫助人們獲取農作區域與作物的現場影像資訊。至今為止已有無數為了達成作物監測的目標，而依據空拍機所拍攝之空拍圖進行影像分析的研究，其研究內容涵蓋

各個領域，如UAV 影像拼接、作物偵測、作物生長階段分類、植物疾病辨識、作

物產量預測等等 [3-7]。

1.2. 研究動機

有鑒於現有的基於空拍圖之作物偵測研究，其偵測品種皆為果樹等體型較大的作物，又或者是以偵測作物生長區域為主 [4, 8, 9]，甚少看到針對水稻進行偵測的研究。若要監測水稻的生長狀況，則需掌握該區域的水稻數量以及植株位置。人工標記水稻會耗費大量人力與時間成本，而如果能利用程式將標記過程自動化，則

(17)

可大幅降低人工標記的所需成本。葉怡萱 [10] 的提出研究雖然是針對水稻秧苗植株的作物偵測研究，但偵測結果並不理想，因此本研究參考葉怡萱的提出方法，進行更多的延伸與發想，希望最後能藉由低空拍攝之空拍圖，從影像中有效偵測出水稻植株與其位置。

1.3. 研究目的

本研究最主要的目標為實現水稻植株標記。為了達成此目標，本研究將其分為四個研究目的：

(一) 稻田區域影像分割：於稻田空拍圖中分割出稻田區域。

(二) 秧苗偵測：偵測出影像中的秧苗植株所在位置。

(18)

(三) 秧苗標記：標記秧苗植株所在稻田區域。

(四) 秧苗對齊：找出涵蓋單一秧苗植株的所有影像編號與其座標。

1.4. 研究貢獻

本研究之實驗產出結果在稻田區域影像分割部分，最高可達到平均精確度 96.26%、召回率 97.88%以及 F 度量值 97.04%；在秧苗偵測部分之偵測結果，最高可達到平均精確度 96.15%、召回率 91.96%以及 F 度量值 93.84%；而秧苗標記部分之標記準確度平均有99.2%；至於秧苗對齊部分，使用提出方法相較於使用 SIFT- based 方法的秧苗對齊數量，平均可增加 26.17%。

而在實現秧苗偵測的過程中，本研究為測試不同植被指數對於偵測結果的影

響，比較了11 種不同植被指數之間的相關程度，這對於日後需挑選不同植被指數

進行比較的研究具有一定的參考價值。同時，本研究設計的參數決定流程，能自動

a

b

c

A B C

(19)

決定斑點偵測所用的參數數值。此流程可應用於當影像中欲偵測之目標，其大小、

形狀皆非常相近時的狀況。此外，秧苗標記的實現，可省去人工標記秧苗所需耗費的大量時間與人力成本。至於秧苗對齊部分，本研究為了找出涵蓋單一秧苗植株的所有影像編號與其座標，提出了一種基於秧苗的影像拼接方法，可使影像拼接時的特徵配對較為平均地分散於影像中。

1.5. 論文架構

論文第一章為緒論，概述本研究之研究背景、研究動機、研究目的與研究貢獻。

第二章相關文獻探討分為六個部分：第一部分介紹可用於輔助作物偵測的不同植被指數；第二部分介紹近年來與作物偵測相關之不同研究；第三部分講述本研究使

用的斑點偵測演算法；第四部分則闡述本研究使用的影像分割方法；第五部分為k-

NN 演算法進行簡單介紹；而第六部分則是針對影像拼接相關研究進行介紹。第三章研究方法首先對本研究所要解決之研究問題提出問題定義，並說明研究整體的系統架構，再解釋各階段的演算法細節。第四章為實驗結果與討論，說明各項實驗之評估方式與本研究所用的實驗資料，並以圖表與結果影像的方式呈現實驗結果，

再根據實驗結果做出比較與分析。至於第五章則是對本研究的研究成果與貢獻做出結論，除此之外也提出後續研究可延伸方向。

(20)

第二章相關文獻探討

2.1. 植被指數

植被指數 (vegetation index) 是依據植被的光譜特性，由至少兩種波段組合而成的數值指標，可以應用在評估一個地區的植被覆蓋率、生長狀況或是植被活力等。

本研究希望能利用植被指數增加被偵測到的秧苗數量，因此首先根據現有的植被指數進行調查。Xue 和 Su [11] 在研究中綜述了超過一百種植被指數，並依據不同環境情況討論了它們的適用性與代表性。每種植被指數都有特定的計算方式，適用在特定用途的同時也有其限制因素。他們表示簡單組合可見光 (visible) 和近紅外光 (NIR) 波段的植被指數就能顯著提高綠色植被偵測的靈敏度 (sensitivity)。然而，

他們的研究中所提到的植被指數大部分都是基於多光譜影像去計算的，而本研究

所使用的實驗影像皆為RGB 彩色影像，因此需要選用僅考慮紅色光、綠色光、藍

色光波段的RGB 植被指數。

Fuentes-Peailillo 等人 [12] 以通過多光譜相機獲得的 NDVI [13]，辨識空拍圖中的土壤和植被做為地表實況 (ground truth)，與 VARI [14]、TGI [15]、ExG [16]、

NGRDI [15] 四種由 RGB 影像得到的植被指數之影像辨識結果互相比較。他們的研究成果顯示，使用RGB 植被指數可辨識出與使用 NDVI 類似的圖樣 (pattern)，

其中以使用TGI 的辨識結果最佳 (如圖 1)。雖然兩者在辨識土壤和植被的像素數

量上的結果近似，但以肉眼觀察時能發現使用RGB 植被指數會在辨識植被時出現

錯誤，這表示若要精確界定土壤與植被的邊界，使用更複雜的分群 (clustering) 技術是必要的。

(21)

圖 1 基於 NDVI 和 TGI 之辨識結果比較 [12]

Hamuda、Glavin 以及 Jones [17] 在影像分割出植物像素與非植物像素的方法中，討論並分析使用不同RGB 植物指數的差異，如圖 2 所示，十種指數對於過曝或過暗的區域皆無法很好地處理。利用NDI [18] 進行分割雖然計算簡單，但分割結果會有大量偽陽性 (false positive)。VEG [19] 對色溫具有不變性，且對於亮度變化有穩健性，但VEG 的計算量較高。ExR [20] 計算簡單，但精確度不如 ExG [16]

高。ExG [16]、CIVE [21]、ExGR [22]、COM1 [23]、MExG [24]、COM2 [25] 對戶外環境皆有不錯的適應性，但CIVE 對於陰影的適應性很差，使用 ExGR 會讓陰影也被分割為植物像素，造成過度分割 (over-segmentation)，使用 COM1 的計算量較高，且由於 CIVE 為組成元素之一，同樣會造成過度分割的問題，而 COM2 的計算量也偏高，但相較於COM1 少。

本研究將上述所提到的RGB 植被指數與其公式，統一整理於表 1。

(22)

圖 2 不同指數在植物像素分割方法上的優缺點比較 [17]

(23)

表 1 RGB 植被指數公式對照表

植被指數公式參考

NDI NDI =G − R

G + R (式 1) [18]

NGRDI NGRDI =G − R

G + R (式 2) [15]

ExG

ExG = 2g − r − b r = ^R^∗

𝑅^∗+𝐺^∗+𝐵^∗, g = ^G^∗

𝑅^∗+𝐺^∗+𝐵^∗, b = ^𝐵^∗

𝑅^∗+𝐺^∗+𝐵^∗

R^∗= ^𝑅

𝑅𝑚𝑎𝑥, G^∗ = ^𝐺

𝐺𝑚𝑎𝑥, B^∗= ^𝐵

𝐵𝑚𝑎𝑥

(式 3) [16]

ExR ExR = 1.3R − G (式 4) [20]

CIVE CIVE = 0.441R − 0.811G + 0.385B + 18.78745 (式 5) [21]

MExG MExG = 1.262G − 0.884R − 0.311B (式 6) [24]

VEG VEG = G

R^aB^1−a a = 0.667

(式 7) [19]

TGI TGI = G − 0.39R − 0.61B (式 8) [15]

VARI VARI = G − R

𝐺 + 𝑅 − 𝐵 (式 9) [14]

ExGR ExGR = ExG − ExR (式 10) [22]

COM1 COM1 = ExG + CIVE + ExGR + VEG (式 11) [23]

COM2 COM2 = 0.36ExG + 0.47CIVE + 0.17VEG (式 12) [25]

(24)

2.2. 作物偵測

在作物偵測方面，Arango Quiroz、Pereira Guidotti 和 Bedoya [8] 根據芒果園的

空拍圖，提出一種自動辨識影像中果樹行列的演算法。提出方法首先將影像由RGB

色彩空間轉換至 YCrCb 色彩空間，接著對 Cr 通道影像進行高斯自適應二值化 (adaptive Gaussian thresholding) 以找出影像中的果樹區域，然後利用霍夫轉換 (Hough transform) 尋找影像中的果樹行列，最後以各行列與果樹區域的交集，標記出所有果樹的中心點位置。

Ahmed 等人 [4] 提出的演算法則是透過多光譜相機所拍攝之空拍圖，自動偵測並分割出金麥豌的育種樣區。他們表示提出方法對於各生長階段的金麥豌皆有不錯的分割成果。圖 3 為提出方法之流程圖：首先計算植被指數影像，經過對比度受限的自適應直方圖均衡化 (contrast limited adaptive histogram equalization, CLAHE) [26] 強化對比度後，以盒狀濾波器 (box filter) 去除影像雜訊。接著，使用 LoG 斑點偵測 [27] 偵測作物所在位置，然後利用隨機漫步 (random walker) 演算法 [28] 於所偵測之斑點周圍分割出作物區域，最後產生作物分割區域遮罩影像。

圖 3 Ahmed 等人提出方法之流程圖 [4]

Chebrolu、Läbe 以及 Stachniss [29] 根據跨時期的作物土地空拍圖，提出一種穩健的跨時期影像對準 (image registration) 演算法。他們發現作物缺少點相較於作

(25)

取代一般的特徵偵測演算法。提出方法中利用植被指數ExG 與 Otsu 演算法 [30]，

產生植物遮罩影像，以此尋找作物點與作物缺少點。至於描述符則是利用k-NN 演

算法 [31]，以作物缺少點與其最近的 4 個缺少點之相對方向與距離建構而成。

不同於上述研究針對的作物植株在影像中皆有一定大小，葉怡萱 [10] 研究中所針對的作物為早期水稻秧苗，其植株密集分布於稻田區域，且單一植株在影像中所佔面積甚至不到十萬分之一。葉怡萱根據空拍稻田圖進行秧苗偵測。圖 4 為其提出方法中秧苗標記之流程圖。研究中使用 SLIC 超像素分割演算法 [32] 將輸入影像分割為多塊大小相近的子區域後，藉由建立的分類資料集將每個子區域分類為稻田與非稻田，然後合併被分為稻田類別的所有子區域作為影像中的稻田區域。

接著，將輸入影像轉換為灰階影像後，於稻田區域進行斑點偵測，最後紀錄所得斑點座標完成秧苗標記。

圖 4 葉怡萱提出之秧苗標記方法流程圖 [10]

(26)

2.3. 斑點偵測

斑點偵測 (blob detection) 中的斑點，指的是影像中擁有與周遭不同性質的區

域，該區域通常與周遭有顏色或亮暗上的差異。常見的斑點偵測演算法會利用LoG

(Laplacian of Gaussian) [33]、DoG (Difference of Gaussian) [34]、DoH (Determinant of Hessian) [27] 等方法來進行，整體來說皆為透過尋找局部極值的方式找出斑點位置。

本研究使用 OpenCV [35] 提供函式進行斑點偵測。以下根據函式的運算原理進行介紹：演算法主要由四個步驟構成。首先，輸入影像依照給定的閾值間隔 𝑡ℎ𝑟𝑒𝑠ℎ𝑜𝑙𝑑𝑆𝑡𝑒𝑝，依低至高閾值進行迭代式影像二值化，經過 k 次迭代後產生 k 張二值化影像。若令輸入影像為 𝐼(𝑥, 𝑦)，第 k 次迭代產生之二值化影像為 𝐼_𝑘(𝑥, 𝑦)，

𝑡_𝑘 為該次二值化使用之閾值，則第一個步驟可以 (式 13) 表示：

𝐼_𝑘(𝑥, 𝑦) = {1, 𝐼(𝑥, 𝑦) ≥ 𝑡_𝑘

0, 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒 (式 13) 其中， 𝑡_𝑘 = 𝑚𝑖𝑛𝑇ℎ𝑟𝑒𝑠ℎ𝑜𝑙𝑑 + 𝑘 ∗ 𝑡ℎ𝑟𝑒𝑠ℎ𝑜𝑙𝑑𝑆𝑡𝑒𝑝

𝑘 = [0,⌊𝑚𝑎𝑥𝑇ℎ𝑟𝑒𝑠ℎ𝑜𝑙𝑑 − 𝑚𝑖𝑛𝑇ℎ𝑟𝑒𝑠ℎ𝑜𝑙𝑑 𝑡ℎ𝑟𝑒𝑠ℎ𝑜𝑙𝑑𝑆𝑡𝑒𝑝 ⌋]

完成影像二值化後，可於各二值化影像找出連通體 (connected component) 與其中心位置，接著依據不同二值化影像之間的連通體距離來決定他們是否屬於同一個斑點，最後重新計算斑點中心位置與大小作為輸出結果。

(27)

2.4. 影像分割

影像分割 (image segmentation) 指的是將一張數位影像分割成多個子區域的過程，每個子區域內的像素之間具有相似的性質，如顏色、明暗度、紋理 (texture) 等 [36]。Masmoudi、Zennouhi 和 El Ansari [37] 提出一種基於二維直方圖的彩色影像分割演算法。他們比較了利用RGB 色彩空間二維直方圖、HSV 色彩空間一維直方圖以及 HSV 二維直方圖在影像分割的結果差異，其研究成果顯示利用 HSV 二維直方圖的彩色影像分割演算法能得到更好的分割結果，且該方法可使用在某些農業應用上以分離植被與土壤。本研究透過影像分割演算法與直方圖比較 (histogram comparison) [38] 將空拍圖分割成稻田區與非稻田區。以下依據本研究所使用的三種影像分割演算法進行介紹。

2.4.1. 分水嶺演算法

分水嶺演算法 (watershed algorithm) [39, 40] 的核心概念源於地理學的分水嶺，

係指一地形中將河流流域分隔開來的界線。演算法將灰階影像視為一種地形圖，像素值代表地勢高度，如此一來，影像中的極小值 (local minimum) 會轉化為盆地。

現在假設該地形因暴雨而逐漸氾濫，最後形成數個湖泊。水淹沒的高度會決定湖泊數量，而區隔相鄰湖泊的山脊即為分水嶺。其中，湖泊即為分割之子區域，而分水嶺為子區域邊界，如圖 5 所示。

(28)

圖 5 以洪水氾濫為原則的分水嶺演算法 [40]

2.4.2. Felzenszwalb-Huttenlocher 演算法

Felzenszwalb 和 Huttenlocher 提出了一種基於圖的影像分割演算法 [41]，以下簡稱為Felzenszwalb 演算法。影像在演算法中被視為無向圖 (graph)，像素為圖中的節點 (node)，鄰近像素以邊 (edge) 連接，而邊上的權重 (weight) 則由像素間的相異程度決定。每個像素在初始化階段都視為一個子區域，然後在尋找最小生成樹 (minimal spanning tree) 的過程中逐漸合併，最後所得之分割區域都是由各個最小生成樹中的節點集合之像素所構成。

2.4.3. SLIC 超像素分割演算法

SLIC (Simple Linear Iterative Clustering) 超像素分割演算法 [32] 是一種改良版本的 k-means 分群演算法 [42]，透過不斷更新群集中心位置來決定分割的子區域邊界。SLIC 將影像由 RGB 色彩空間轉換至 CIELAB 色彩空間進行考慮，並透過參數 𝑘 決定希望得到的分割子區域數量。SLIC 在初始化階段，𝑘 個群集中心點 C_𝑘 會被平均分布於影像中，然後將每個像素分配到與其距離最近的群集。當像

(29)

素皆被分配完成時，則更新群集中心，並重新分配像素所屬群集，以此迭代執行。

迭代完成後的像素分配結果即為分割之子區域。

2.5. k-NN 演算法

k-NN 演算法 (k-nearest neighbors algorithm) [31] 在資料科學中是一種常被使用的分類演算法。圖 6 為演算法之運作概念：以待分類資料點為中心，尋找與其

距離最近的 k 個資料點，再依據這些點的所屬類別，以投票形式決定待分類資料

點的類別。

圖 6 k-NN 演算法 [31]

2.6. 影像拼接

由於本研究在秧苗對齊部分以影像拼接 (image stitching) 的角度進行思考，為此也參考了影像拼接的相關研究。一般來說，影像拼接主要由特徵點偵測 (feature detection) 與描述 (description)、特徵配對 (feature matching)、離群值 (outlier) 拋棄與產生單應性矩陣，以及影像混和 (image blending) 等幾個步驟組成 [43]。

目前已有非常多種特徵點偵測與描述的演算法，其中 Tareen 和 Saleem [43]

針對SIFT [34]、SURF [44]、KAZE [45]、AKAZE [46]、ORB [47] 以及 BRISK [48]

六種特徵偵測演算法進行比較。他們在研究中提到，SIFT、SURT 以及 BRISK 是

(30)

尺度不變性最佳的偵測器，能適應大幅的尺度變化，而ORB 的尺度不變性為六者之中最低；SIFT、KAZE、AKAZE 以及 BRISK 具有較高的影像旋轉準確度；ORB 以及BRISK 較其他四者對仿射變換 (affine transformation) 的穩定性高。雖然六種演算法各有優缺，但對於任意類型的幾何變換，以 SIFT 和 BRISK 整體上的準確度最佳，其中又以SIFT 最為準確。

特徵配對完成後的離群值拋棄與產生單應性矩陣部分，經常使用的演算法是 RANSAC (random sample consensus) [49]。其核心概念為，隨機選擇一些特徵配對產生單應性矩陣，依此計算特徵配對之間的距離，並分出內群值 (inlier) 與離群值，

接著不斷重複這個過程，直到產生的單應性矩陣使內群值數量最多，如圖 7 所示。

圖 7 RANSAC：內群值與離群值 [49]

(31)

第三章研究方法

3.1. 問題定義

本研究依據前述所提之研究目的提出以下三個研究問題：

（１）從影像中抓取稻田區可被定義為將影像簡易區隔為稻田與非稻田區的影像分割問題。此類問題的深入應用為將影像中的每個像素都劃分到一個類別的語義分割 (semantic segmentation) 問題。

（２）從局部來看，秧苗與周圍泥土在肉眼中的差異十分明顯，在影像中應該也會被認為是一個特徵點。而一張資訊豐富的影像會擁有大量且不一的特徵點，

因此本研究先假設秧苗為特徵點，並將秧苗偵測問題定義為尋找具有特定性質的特徵點偵測問題。

（３）為了得知單株秧苗出現的影像與位置，勢必要先對齊相鄰影像之間的秧苗，

再延伸到處理多張影像的對齊問題。由於本研究將秧苗視為一種特徵點，所以秧苗對齊可視為對齊指定類型特徵點的問題。而影像之間的對齊屬於影像拼接的範圍，因此定義秧苗對齊問題為：基於指定類型特徵點之影像拼接問題。這可以歸類為影像拼接的特殊應用之一。

3.2. 系統架構

本研究之系統架構圖如圖 8 所示，系統架構分為秧苗偵測、秧苗標記與秧苗對齊三大部分。在秧苗偵測部分，本研究參考葉怡萱 [10] 提出之研究方法思路，

由稻田區域影像分割與稻田區域之特徵點偵測這兩大步驟所構成，如圖 9、圖 10 所示。首先利用不同影像分割演算法將輸入影像分割成多個子區域，再將子區域分類為稻田或者非稻田，接著合併所有被分類為稻田的子區域，以此獲得稻田區域。

(32)

進行特徵點偵測之前，輸入影像會經過多種不同的前處理方法，並以一個系統化的評估方式調整每種前處理方法適合的偵測參數，再於前述所提之稻田區域進行特徵點偵測，以此得到秧苗偵測結果。

秧苗標記部分則是使用影像拼接方法，尋找輸入影像與鑲嵌圖 (mosaic) 之間的轉換矩陣，再以所得之秧苗偵測結果將秧苗座標轉換至鑲嵌圖的座標平面上。如此，即可依照事先標記的稻田區域劃分區域與秧苗轉換完的座標位置，標記秧苗所屬的稻田區域。

至於秧苗對齊的系統架構，包含相鄰影像之間的秧苗對齊與多張影像之間的秧苗對齊兩大步驟。在相鄰影像之間的秧苗對齊部分如圖 11 所示，首先將輸入影像等比例分割為 2 × 2 的四個子區域，並分別基於各子區域找出相對應的轉換矩陣，再將輸入影像轉換到相鄰影像的座標平面上，然後進行秧苗對齊得到秧苗對。

由於有四個子區域，因此總共會得到四批秧苗對。接著合併四批秧苗對，並將它們視為影像拼接中的特徵配對，重新計算轉換矩陣後再次尋找新的秧苗對做為相鄰影像之間對齊的秧苗。而在多張影像之間的秧苗對齊部分，主要利用前述相鄰影像拼接所得的轉換矩陣，對於任三張彼此相鄰的影像皆統一轉換到第二張影像的座標平面上，再進行秧苗對齊尋找秧苗對，最後給定秧苗編號並記錄每個編號所在的影像編號與座標。

(33)

圖 8 系統架構圖秧苗偵測稻田區域影像分割

稻田區域之特徵點偵測

秧苗對齊基於秧苗之相鄰影像拼接

多張影像之間的秧苗對齊秧苗標記

座標平面轉換

標記秧苗之稻田區域

(34)

圖 9 稻田區域影像分割架構圖影像分割演算法

輸入影像

子區域與分類資料集進行直方圖比較與 3-NN 分類

合併被分類為稻田之子區域

分割之稻田區域分割之稻田區域

分割之稻田區域 Felzenszwalb

演算法分水嶺演算法

SLIC 演算法

利用區域相鄰圖的子區域合併

形態學運算

(35)

圖 10 稻田區域之特徵點偵測架構圖輸入影像

灰階

灰階 COM2ExR MExG TGI

ExG

前處理影像

前處理影像前處理影像

前處理影像前處理影像前處理影像

前處理影像

特徵點偵測

特徵點偵測特徵點偵測

特徵點偵測特徵點偵測特徵點偵測

特徵點偵測

補數處理 (Complement)

對比度受限的自適應直方圖均衡化 (CLAHE) 提升影像對比度

(36)

圖 11 相鄰影像之間的秧苗對齊架構圖輸入影像

相鄰影像

子區域 1 子區域 2 子區域 3 子區域 4

基於子區域尋找轉換矩陣

輸入影像轉換到相鄰影像之座標平面

尋找秧苗對

尋找秧苗對尋找秧苗對尋找秧苗對

基於秧苗對尋找轉換矩陣

輸入影像轉換到相鄰影像之座標平面合併秧苗對

尋找秧苗對

(37)

3.3. 稻田區域影像分割

由於秧苗皆被種植在稻田中，因此本研究首先使用影像分割演算法將影像分割成數個子區域，再藉由分類子區域將影像分割為稻田像素與非稻田像素。

3.3.1. 影像分割演算法使用與比較

本研究使用 scikit-image [50] 所提供之模組，以 2.4. 介紹的分水嶺演算法、

Felzenszwalb 演算法與 SLIC 演算法，比較影像分割出的子區域特性，並在第四章討論使用不同演算法所分割出的稻田區域有何差異。分水嶺演算法如圖 12 所示，

藉由產生梯度影像將一張影像視為地形圖，再根據給定的閾值 (threshold) 假想該地形的水平面高度，以此標記盆地所在位置並尋找分水嶺作為影像分割之子區域邊界。

圖 12 分水嶺演算法

(38)

如圖 13 所示，給定不同的閾值會有不一樣的分割結果。本研究經過測試不同閾值後將使用於實驗中的閾值設定為30。

圖 13 分水嶺演算法使用不同閾值之分割結果

Felzenszwalb 演算法有 3 個可調控參數：scale 的值會間接影響分割的子區域大小與數量，值越大代表分割出的子區域面積越大且總數量越少；sigma 的值會影響高斯核 (Gaussian kernel) 大小，用於分割前先對影像進行平滑處理 (smoothing)，

值越高代表分割前的影像越平滑；min_size 則是用來限制分割出的子區域最小面積。本研究將sigma 維持預設值，藉由調整另外兩個參數得到不同分割結果，如圖 14 所示，最後設定實驗所用之參數為：scale = 500，min_size = 3000。

(39)

圖 14 Felzenszwalb 演算法使用不同參數之分割結果

至於 SLIC 演算法中的參數，由 n_segments 和 compactness 主導分割結果：

n_segments 可決定所分割子區域大致的總數量；compactness 則會影響分割時顏色接近度 (color proximity) 與空間接近度 (space proximity) 之間的權重，值越高代表空間接近度的權重越大，分割出的子區域形狀會越接近正方形。本研究將 compactness 維持預設值，並設定 n_segments 為 1500，使分割出的子區域大小適當，能大致分出田埂與稻田的邊界，如圖 15 所示。

(40)

(a) 原始影像 (b) 分割影像，n_segments = 1500 圖 15 SLIC 演算法分割結果

產生多個子區域之後，計算每個子區域的顏色平均值並建立區域相鄰圖 (region adjacency graph, RAG)，即可藉由閾值，將顏色相近的子區域合併為較大的子區域 [51]。圖 16 為利用 RAG 與不同閾值合併子區域後的分割結果，由肉眼觀察最後選擇將此處閾值設為15。

(a) threshold = 10 (b) threshold = 15 (c) threshold = 20 圖 16 利用 RAG 與不同閾值的子區域合併結果

(41)

為了方便稱呼，論文中將（１）利用分水嶺演算法（２）利用Felzenszwalb 演算法與（３）利用SLIC 演算法加上 RAG 子區域合併，三種稻田區域影像分割方法分別簡稱為Ws、Fz、SL_R。圖 17 為三種方法之影像分割結果，表 2 為產生之子區域數量。觀察圖 17 與表 2 可看出，使用 Ws 所分割出的子區域數量是最龐大的，這是由於田埂、秧苗等處像素值的梯度變化相較於泥土區明顯地多，因此在這

些區域很容易被分割為大量且面積非常小的子區域。使用Fz 所分割出的子區域數

量次中，但缺點是無法直接從參數得知最後分割出的子區域數量，需經多次測試才能決定參數。而使用SL_R 所得分割結果的子區域數量是最少的，其優點是 SLIC 能直接從參數得知大略的子區域數量，從而能計算合併完後的子區域的最小面積為何。

(a) Ws (b) Fz (c) SL_R

圖 17 使用不同演算法之影像分割結果表 2 使用範例影像產生之子區域數量

使用方法 Ws Fz SL_R

子區域數量 25103 275 142

(42)

3.3.2. 直方圖比較與 k-NN 演算法

本階段之目的在於將輸入影像分類為稻田像素或者非稻田像素。在前一階段

獲得分割結果後，將每個子區域與建立好的分類資料集進行直方圖比較，依k-NN

演算法將子區域分類為稻田子區域或者非稻田子區域。分類資料集之中包含資料集影像與標記影像，其中表 3 列出標記影像之類別與張數：稻田類別 7 張，非稻田類別19 張，一共 26 張標記影像。考量到稻田類別只有稻田一種元素，而非稻田類別中有建築物、田埂、道路等多種元素，因此非稻田類別需要較多張標記影像作為資料集。標記影像皆為二元影像 (binary image)，如圖 18 所示，白色像素代表標記為該分類之像素。

表 3 分類資料集中的標記影像類別與張數

類別稻田非稻田

標記影像張數 7 19

圖 18 分類資料集之標記影像

本研究利用HSV 色彩空間之色相 (Hue) 與飽和度 (Saturation) 進行二維直方圖比較。因此資料集影像會先被轉成HSV 影像，再將標記影像作為遮罩 (mask)，

計算資料集影像於遮罩區域內的HS 直方圖，如圖 19 所示。資料集建立完畢後，

於待分類子區域計算 HS 直方圖 (圖 20) 並與資料集進行直方圖比較，可得到子

(43)

區域與資料集之相關係數排序。接著即可依據子區域與資料集中各標記的相似程度，使用k-NN 演算法將子區域分類為稻田或者非稻田。本研究將 k 設定為 3，也就是從與子區域相關係數前三高之對應標記類別中，依佔多數的類別為該子區域之分類結果。若子區域被分類為稻田，則將該區域作為遮罩加入結果，最後會得到輸入影像中所有被分類為稻田像素的遮罩影像作為初步結果。

圖 19 分類資料集內的 HS 直方圖

(44)

圖 20 待分類子區域與其 HS 直方圖

3.3.3. 影像分割後處理

獲得初步的稻田區域遮罩後，本研究對該區域進行形態學開運算清除可能分割不理想的邊緣後，尋找遮罩的凸多邊形讓稻田區域的形狀更為完整。圖 21 (a) 為初步之稻田區域影像分割結果，圖 21 (b) 為後處理完畢之最終結果。

(a) 稻田區域影像分割之初步結果 (b) 經過後處理之最終結果圖 21 稻田區域影像分割之結果圖

圖 22 為 Ws、Fz 以及 SL_R 的最終稻田區域影像分割結果圖。詳細的結果比較與分析將於第四章進行討論。

(45)

(a) Ws (b) Fz (c) SL_R 圖 22 使用不同方法之稻田區域影像分割結果圖

3.4. 稻田區域之特徵點偵測

為了確認若於特徵點偵測前將影像進行不同前處理，對偵測結果是否產生正面影響，本研究以僅將輸入影像轉換成灰階影像之前處理作為對照組，以其他有經過對比度強化或者使用植被指數影像之前處理作為實驗組。

3.4.1. 分析植被指數間的相關程度

如第二章所提，現有植被指數的數量十分龐大，僅考慮 R、G、B 三個通道的 RGB 植被指數也在 2.1. 提到了十二種。其中，NDI 與 NGRDI 在計算公式上完全一樣，因此可將兩者視為相同，並在兩者以 NDI 作為代表。圖 23 為將原始影像

邊長縮小為0.1 倍，並經過灰階處理以及經過不同植被指數計算後所得影像。為了

方便與對照組比較，NDI、ExG、MExG、VEG、TGI、VARI、ExGR 以及 COM1 之植被指數影像又再進行補數處理 (式 14)，其中 I 為輸入影像，I_C 為經過補數處理之影像，(x, y) 為影像中的像素位置。

I_C(x, y) = 255 − I(x, y) (式 14)

(46)

圖 23 植被指數影像

圖 24 調整後之植被指數影像

(47)

調整後的結果如圖 24 所示。接著，根據植被指數影像每個位置的像素值，可以計算出植被指數兩兩之間的皮爾森相關係數 (Pearson’s correlation coefficient)，

然後建立相關矩陣 (correlation matrix) 來衡量植被指數之間的線性相依程度。以下為變數X 與 Y 之間的皮爾森相關係數定義 [52]：

ρ_X,Y =cov(X, Y)

𝜎_𝑋𝜎_𝑌 (式 15)

其中，cov(X, Y)為 X 和 Y 之共變異數 (covariance)，σ為標準差 (standard deviation)。

公式如下，給定n 對觀測值 (observations) (𝑋₁, 𝑌₁), ⋯ , (𝑋_𝑛, 𝑌_𝑛)：

cov(X, Y) = 1

n − 1∑(𝑋_𝑖− 𝜇_𝑋)

𝑛

𝑖=1

(𝑌_𝑖 − 𝜇_𝑌) (式 16)

𝜎_𝑋 = √ 1

n − 1∑(𝑋_𝑖 − 𝜇_𝑋)²

𝑛

𝑖=1

, 𝜎_𝑌 = √ 1

n − 1∑(𝑌_𝑖 − 𝜇_𝑌)²

𝑛

𝑖=1

(式 17)

𝜇_𝑋 = 1 𝑛∑ 𝑋_𝑖

𝑛

𝑖=1

, 𝜇_𝑌 = 1 𝑛∑ 𝑌_𝑖

𝑛

𝑖=1

(式 18)

𝜇為變數之平均數。由以上兩式可將 (式 15) 改寫為 (式 19)：

ρ_X,Y= ∑^𝑛_𝑖=1(𝑋_𝑖 − 𝜇_𝑋)(𝑌_𝑖− 𝜇_𝑌)

√∑^𝑛_𝑖=1(𝑋_𝑖 − 𝜇_𝑋)²√∑^𝑛_𝑖=1(𝑌_𝑖 − 𝜇_𝑌)²

(式 19) 套用在影像中，則n 代表總像素數量，𝑋_𝑖、𝑌_𝑖 分別代表影像 X 與影像 Y 的第 i 個像素值。相關係數的值會落在 −1 與 1 之間，越靠近 1 代表兩變數正相關程度越高，反之越靠近 −1 代表兩變數負相關程度越高，而相關係數的值若等於 0 代表兩變數無線性相依。

圖 25 為計算出植被指數彼此間相關係數後所建立的相關矩陣。由圖中框起處可發現，ExG 與 VEG 完全相關，與 COM1 高度相關；ExR 與 ExGR 幾乎完全相關；MExG 與 NDI、VARI 高度相關；TGI 與 CIVE 幾乎完全相關；而 COM2 與其

(48)

他植被指數的相關程度都不高。經過以上分析，本研究決定挑選彼此相關程度相對較低的ExG、ExR、MExG、TGI 與 COM2 作為實驗所使用的植被指數，如圖 26 所示，相關係數的值最高不超過0.95，最低幾乎等於 0。

圖 25 植被指數相關矩陣

圖 26 選用之植被指數相關矩陣

(49)

3.4.2. 產生前處理影像

選定使用於實驗的 5 種植被指數後，本階段對灰階影像與植被指數影像皆使

用對比度受限的自適應直方圖均衡化 (CLAHE) [26] 來增強影像對比度，如圖 27

所示。加上未做其它處理的灰階影像，本階段產生的前處理影像總共有7 種，如圖

28，以下將 7 種前處理方法分別簡稱為 Gray、CLAHE、ExG、ExR、MExG、TGI 以及COM2。其中 Gray 為對照組，其餘六者為實驗組。

(a) 強化對比度前 (b) 強化對比度後

圖 27 使用 CLAHE 前後比較圖

圖 28 前處理影像

(50)

3.4.3. 斑點偵測與參數設定

本研究使用 OpenCV [35] 所提供之 SimpleBlobDetector 針對前處理影像的稻田區域執行斑點偵測，以此偵測秧苗。函式可以依參數來決定是否對偵測斑點的顏色亮暗、面積大小、形狀圓度、形狀凸包度，以及形狀長短比進行篩選。為了決定這些參數，本階段首先使用Otsu 演算法 [30] 將稻田灰階影像二值化，如圖 29 所

示。觀察結果可發現秧苗為黑色像素，且所占面積幾乎都在15 個像素以上，此外

由於秧苗葉片方向不一，二值化完後會呈不規則狀。

(a) 稻田灰階影像 (b) 二值化影像

圖 29 影像二值化前後比較圖

由此設定實驗所用之斑點偵測參數為：

parameters. filterByColor = True parameters. blobColor = 0 parameters. filterByArea = True parameters. minArea = 15 parameters. maxArea = 300

parameters. filterByCircularity = False parameters. filterByConvexity = False parameters. filterByInertia = False

(51)

除了以上依肉眼觀察設定完成的參數，仍需針對每種前處理影像挑選適當的 parameters.thresholdStep 數值。此項參數控制輸入影像在斑點偵測過程中，數次二值化期間的閾值間隔。本階段設計一個系統化的流程用以決定參數。首先於原始影像切出一塊 500 × 500 像素的稻田塊 (圖 30)，人工計算稻田塊中的秧苗數量後，

再依thresholdStep 預設值的偵測結果，校正秧苗數量作為比較基準。接著針對稻田塊的每種前處理影像，設定thresholdStep 由 50 至 1，測試使用不同參數值所得斑點數量 (圖 31)，最後選擇使偵測所得斑點數量與實際秧苗數量最為接近的 thresholdStep 值，作為該前處理影像的參數值。

秧苗數量 = 252

thresholdStep = 30 秧苗數量 = 90

thresholdStep = 11 秧苗數量 = 248

圖 30 稻田塊示意圖圖 31 使用不同 thresholdStep 之偵測結果示意圖經過上述流程即可為各前處理方法選定斑點偵測中的 thresholdStep 設定值，

如表 4 所列。

表 4 各前處理方法之 thresholdStep 設定值

前處理影像 Gray CLAHE ExG ExR MExG TGI COM2 thresholdStep 11 35 1 21 17 15 15

(52)

斑點偵測之各項參數皆設定完成後，本階段利用 3.3. 所得之稻田區域影像分割，於各前處理影像中的稻田區域進行斑點偵測，並將偵測結果作為秧苗偵測結果。

同時從結果影像擷取 500 × 500 像素大小之區域作為結果示意圖，如圖 32。

(a) 原始影像 (b) Gray (c) CLAHE

(d) ExG (e) ExR (f) MExG

(g) TGI (h) COM2

圖 32 秧苗偵測結果示意圖

(53)

3.5. 秧苗標記

此階段目的在於標記所偵測之秧苗位於哪一個稻田區域，為此會涉及影像平面轉換的問題 [53, 54]。圖 33 為空間點透過相機中心的投影示意圖。影像平面上的點 x_i 為空間中的點 X_i 指向相機中心 c 之射線與影像平面的交點，如圖 33 之 a。若空間點共面，則可以稱空間中的平面與影像平面之間存在投影轉換 (projective transformation) 的關係，亦稱為單應性 (homography)。圖 33 之 b 可以想像成實際場景透過相機投影到影像平面，而所有相同相機中心的影像都藉由投影轉換進行連結，如圖 33 之 c 為影像平面之間的映射 (mapping)。空間中的點若透過不同相機中心分別投影到影像平面上，則影像平面之間通常無法由投影轉換進行連結，如圖 33 之 d，但若空間中的點皆共平面，如圖 33 之 e，則影像平面之間會存在單應性。

(54)

圖 33 空間點透過相機中心的投影示意圖 [53]

(55)

根據定義，平面投影轉換 (圖 34) 是由一 3 × 3 非奇異 (non-singular) 矩陣所構成之齊次座標系3 維向量的線性變換，可以表示為：

[ 𝑥₁^′ 𝑥₂^′ 𝑥₃^′

] = [

ℎ₁₁ ℎ₁₂ ℎ₁₃ ℎ₂₁ ℎ₂₂ ℎ₂₃ ℎ₃₁ ℎ₃₂ ℎ₃₃

] [ 𝑥₁ 𝑥₂

𝑥₃] (式 20)

或是以向量形式表示：

𝒙^′= 𝐻𝒙 (式 21)

圖 34 平面上一點映射至另一平面之示意圖 [53]

本階段藉由尋找輸入影像與鑲嵌圖之間的單應性矩陣 (homography matrix)，

將輸入影像轉換到鑲嵌圖之影像平面上，以判斷秧苗落在哪個稻田區域。為了尋找單應性矩陣，首先使用SIFT 演算法對輸入影像與鑲嵌圖進行影像拼接 (圖 35 (a))，

以此獲得單應性矩陣。接著將 3.4. 偵測出的秧苗座標轉換到鑲嵌圖的座標平面上，

如此即可依據在鑲嵌圖上已劃分好的稻田區域範圍，為秧苗進行標記，如圖 35 (b) 所示。現在令 𝐼 為輸入影像之像素座標集合，𝐻 為輸入影像與鑲嵌圖之間的單應性矩陣，則 𝑛 株秧苗 𝑟 之座標轉換可以表示成 (式 22)：

𝒓_𝒊^′= 𝐻𝒓_𝒊，其中 𝑟_𝑖 ∈ 𝐼，𝑖 = 1, … , 𝑛 (式 22)

(56)

(a) 輸入影像轉換到鑲嵌圖上 (b) 依稻田區域劃分標記秧苗點圖 35 秧苗標記示意圖

3.6. 秧苗對齊

為了得知單株秧苗出現的影像與位置，此小節先討論相鄰影像之間的秧苗對齊問題，再延伸到多張影像之間的秧苗對齊問題。其中，秧苗對齊以影像拼接的角度進行思考。以下針對如何判定秧苗是否為一組配對進行說明：

基於計算方便，本研究使用無窮範數 (infinity norm) 來定義兩張影像的秧苗在拼接影像上的距離。現在令 I_input 為輸入影像，I_adj 為相鄰影像，𝑞、𝑟 分別為輸入影像與相鄰影像上的秧苗，I_input 轉換至 I_adj 之單應性矩陣為 H。依 (式 22) 轉換 𝑞 之座標為 𝑞^′，使兩秧苗位於同一座標平面後，即可計算兩秧苗的距離。

定義秧苗 𝑞^′、𝑟 之間的距離 𝑑(𝑞^′, 𝑟) 為：

∀ 𝑞 ∈ I_input，𝑟 ∈ I_adj，

𝑞^′= 𝐻𝑞

𝑑(𝑞^′, 𝑟) = ‖𝒓 − 𝒒^′‖_∞∶= max(|𝑟_𝑥− 𝑞_𝑥^′|, |𝑟_𝑦− 𝑞_𝑦^′|) (式 23)

(57)

，其中 𝒒 = (𝑞_𝑥, 𝑞_𝑦)，𝒓 = (𝑟_𝑥, 𝑟_𝑦)。圖 36 為秧苗之間的距離範例圖。當 𝒓 落在圖中的正方形邊上時，𝑑(𝑞^′, 𝑟) = 1，代表輸入影像與相鄰影像上的秧苗彼此之間的 𝑥 座標差和 𝑦 座標差之最大值為 1。

圖 36 秧苗之間的距離範例圖

而一組秧苗對是否對齊，由它們彼此之間是否足夠靠近決定。給定一閾值 𝑡，

若 𝑞^′ 與 𝑟 彼此互為距離最近的秧苗，且距離小於 𝑡 時，則稱 (𝑞, 𝑟) 為一對齊之秧苗對，如圖 37 所示。

if 𝑑(𝑞^′, 𝑟) < 𝑡:

𝑟 = 𝑓𝑖𝑛𝑑_𝑐𝑙𝑜𝑠𝑒𝑠𝑡(𝑞) = arg min

𝑟𝑖∈𝐼𝑎𝑑𝑗

(𝑑(𝑞^′, 𝑟_𝑖))

圖 37 對齊之秧苗對 (𝑞, 𝑟) 示意圖 𝑑(𝑞, 𝑟) = 1 𝑞^′

𝑟 𝑥 𝑦

1 1

𝑞 𝑟

𝑥 𝑦

𝑡

(58)

3.6.1. 相鄰影像之間的影像拼接與秧苗對齊

一般常見的影像拼接演算法透過特徵比對尋找單應性矩陣來完成。圖 38 為使用SIFT 演算法的影像拼接示意圖。由圖 38 (d) 可看出經 RANSAC 篩選後留下的內群對集中在田埂上，這會造成影像拼接結果在肉眼觀察下，靠近田埂處較為理想，

遠離田埂處則會產生重影，如圖 38 (e) 所示。

(59)

(a) 輸入影像 (b) 相鄰影像

(c) 特徵比對結果，藍線為對應的特徵配對

(d) 經 RANSAC 篩選後之內群對

(e) 影像拼接結果圖 38 SIFT 影像拼接示意圖

(60)

若要使影像各區域的拼接理想程度不會差異太大，則需要令特徵配對盡可能

地均勻分布於影像上，如此一來才不會讓RANSAC 篩選後的內群對分布過於密集，

導致影像拼接以內群對的密集區為主。由以上想法，且考慮到此階段目的在於使對

齊的秧苗越多越好，本研究提出一種基於秧苗的影像拼接方法，核心概念是以SIFT

演算法進行影像拼接，先得到單應性矩陣 H 並以此找出秧苗對後，將秧苗對作為新的特徵配對，再進行一次影像拼接演算法，校正先前得出的單應性矩陣，最後再依校正過的單應性矩陣尋找秧苗對，作為秧苗對齊的結果。

為了使找到的秧苗對分佈不要過於集中，本階段先將輸入影像等比例劃分為左上、右上、左下、右下四個區域，基於不同區域進行影像拼接，如圖 39、圖 40 所示，再分別根據所獲得的單應性矩陣尋找秧苗對，此處將 𝑡 設為 5 像素，以減少秧苗配對錯誤的可能性。接著，合併所有秧苗對並刪除重複配對，然後以合併後的秧苗對作為新的特徵配對，得出新的單應性矩陣 H^′。最後重新尋找秧苗對做為秧苗對齊結果，此處將 𝑡 設為 10 像素。

本階段將秧苗對齊方法中，使用SIFT 演算法的影像拼接方法作為對照組，簡稱為 SIFT-based 方法，而以上述所提使用秧苗對做為特徵配對的影像拼接方法作為實驗組，簡稱為seedling-based 方法。

(61)

(a) 基於左上區域 (b) 基於右上區域

(c) 基於左下區域 (d) 基於右下區域

圖 39 基於不同區域之影像拼接內群對

(a) 基於左上區域 (b) 基於右上區域

(c) 基於左下區域 (d) 基於右下區域

圖 40 基於不同區域之影像拼接結果

(62)

(a) SIFT-based 之秧苗對齊結果 (對照組)

(b) seedling-based 之秧苗對齊結果 (實驗組) 圖 41 對照組與實驗組之秧苗對齊結果

圖 41 為對照組與實驗組之秧苗對齊結果範例影像，其中紅點與藍點分別屬於不同影像的秧苗點，而結果範例影像中只顯示對齊之秧苗對。

(63)

3.6.2. 多張影像之間的影像拼接與秧苗對齊

當影像拼接涉及到多張影像之間的拼接，就勢必要考慮累積誤差的問題了。影像經過平移 (translation)、縮放 (scaling)、旋轉 (rotation)、推移 (shearing) 等平面線性變換時，由於影像是由像素所構成，因此需要對影像進行重新取樣 (resampling)。如圖 42 所示，當影像在重新取樣時，落在同一格像素內的紅點座標皆會以黑點之像素座標表示。這代表影像中的點在轉換過程中，與實際狀況相比，

最多可在 𝑥 方向及 𝑦 方向各產生 1 像素的誤差。

圖 42 影像重新取樣示意圖

現在針對多張影像之間的拼接，考慮如圖 43 之範例。A、B、C 為三張大小相等的影像，寬度皆為300 像素。其中，AB、BC 之間彼此有三分之二的面積重疊，

AC 之間有三分之一的面積重疊。若要對 AC 進行影像拼接，可經由兩種方式達成：

一種是 AC 直接進行影像拼接，另一種則是透過 AB、BC 之間的影像拼接達成。

若將影像座標平面轉換至另一影像座標平面的過程視為一個路徑 (path)，而影像座 1px

1px

(64)

標平面視為節點，則上述兩種方式可以表示成 (A → C) 與 (A → B → C)。(式 24)、

(式 25) 表示藉由兩種路徑將 A 轉換至 C 座標平面的過程。其中，H_AB 代表 A 轉換至 B 座標平面的單應性矩陣，以此類推。

路徑一 (A → C) A^′ = H_ACA (式 24)

路徑二 (A → B → C) A^′= H_BCH_ABA (式 25)

圖 43 多張影像之間的拼接示意圖

理想情況下的影像拼接，會使特徵配對完美重疊。假設每次進行影像拼接時，

皆會於重疊區兩端的 𝑥、𝑦 方向產生 1 像素的誤差，則影像拼接於一張影像兩端所產生的誤差，會與重疊區所佔比例有關，而多次影像拼接所產生的總誤差為每次拼接誤差的乘積。以 𝑥 方向為例，(式 26) 表示單次影像拼接在各方向所產生的誤差：

1 𝑝𝑥

重疊區寬度= 𝑥 𝑝𝑥

影像寬度⇒ 𝑥 = 影像寬度

重疊區寬度 (式 26)

A C B

100px 200px

200px

300px

(65)

以下分別為路徑一 (A → C) 與路徑二 (A → B → C) 在 A 之兩端所產生的 𝑥 方向誤差：

路徑一 1

100= 𝑥

300⇒ 𝑥 =300 100= 3

路徑二 (A → B)：

1

200= 𝑥_𝑎𝑏 300 (B → C)：

1

200= 𝑥_𝑏𝑐 300 (A → B → C)：

𝑥 = 𝑥_𝑎𝑏× 𝑥_𝑏𝑐 = 300

200×300

200= 2.25

而 𝑦 方向亦同，因此選擇路徑二對 AC 進行影像拼接所得誤差較小。

由上述例子可發現影像拼接時，誤差與所經路徑之重疊率乘積倒數成正比，即誤差與所經路徑之重疊率乘積成反比，因此當 AC 之間的影像拼接有 A → C 與 A → B → C 兩種路徑可達成，且 AB 之間重疊率為 r_AB、BC 之間的重疊率為 r_BC、 AC 之間的重疊率為 r_AC 時，應選擇所經路徑之重疊率乘積較大的路徑作為影像拼接之轉換路徑。

為了測試影像拼接時的路徑選擇重疊率乘積較大的路徑，於實際情況中是否真的會使得拼接結果較為理想，本研究將相鄰的三張影像 ABC 中，以路徑一 (A → C) 進行影像拼接的方式作為對照組，以路徑二 (A → B → C) 進行影像拼接的方式作為實驗組，比較兩種路徑在 A ∩ B ∩ C 區域中的秧苗對齊數量。

而有了估算誤差與選擇轉換路徑的方法後，即可以圖的概念考慮影像拼接問