圖像分割中的自動求閾值法

撰文/熊氏阿回

圖像分析里，常常需要確定某一個閾值，把圖像中的不同類別分割開來，比如水體與陸地、土壤與植被等等。

比如上圖中，如果找到一個合適的閾值，就可以把水體比較完整地分割出來。打開直方圖，憑直覺我們知道閾值就存在於紅色區域中，只是我們不知道具體在哪裡：

2016，荷蘭科學家完成了全球地表水製圖，要解決的也就是這個問題。

每輻圖像的閥值會有所不同，如果要完成全球製圖，需要找到一個自動化程度高，用輸入的影像數據驅動的演算法，而不是在演算法里搞很多的判斷。這就需要用到圖像分割中的 otsu 方法，又稱「大津法」，是 1979 年日本學者大津展之提出的。為了找到這個閾值，他把問題定義成：

假定一張圖片共有n個像素，其中灰度值小於閾值的像素為 n1 個，大於等於閾值的像素為 n2 個（ n1 + n2 = n ）。w1 和 w2 表示這兩種像素各自的比重。而所有灰度值小於閾值的像素的平均值和方差分別為 μ1 和 σ1，所有灰度值大於等於閾值的像素的平均值和方差分別為 μ2 和 σ2。於是，可以得到

類內差異 = w1(σ1的平方) + w2(σ2的平方)

類間差異 = w1w2(μ1-μ2)^2

要找到合適的閾值，讓前者最小或者後者最大都可以，從計算角度出發，後者要好計算一些。又被稱為 BSS：

具體的做法，是可以先對圖像求一下直方圖，用一系列從小到大的閥值去取一下，分別代入上面的算式。使得」類內差異最小」或」類間差異最大」的那個值，就是最終的閾值。

用來驗測水體的數據是紅外波段，陸地上因為有植被，呈紅色，水體則差異較大。

在 Google Earth Engine 上編碼，分四個部分：

1.定義所用的遙感影像，計算這幅圖像的直方圖（ histogram）

// Compute the histogram of the NDVI band.nvar histogram = NDVI.reduceRegion({nreducer: ee.Reducer.histogram(50)n .combine(mean, null, true)n .combine(variance, null, true), n// geometry: polygon, nscale: 300,nbestEffort: truen});n

2.為了應用雲計算，把 otsu 的計算步驟用 map 和 reduce 包裹起來，便於 GEE 調度資源：

var otsu = function(histogram) {nvar counts = ee.Array(ee.Dictionary(histogram).get(histogram));nvar means = ee.Array(ee.Dictionary(histogram).get(bucketMeans));nvar size = means.length().get([0]);nvar total = counts.reduce(ee.Reducer.sum(), [0]).get([0]);nvar sum = means.multiply(counts).reduce(ee.Reducer.sum(), [0]).get([0]);nvar mean = sum.divide(total);nvar indices = ee.List.sequence(1, size);nn// Compute between sum of squares, where each mean partitions the data.nvar bss = indices.map(function(i) {n var aCounts = counts.slice(0, 0, i);n var aCount = aCounts.reduce(ee.Reducer.sum(), [0]).get([0]);n var aMeans = means.slice(0, 0, i);n var aMean = aMeans.multiply(aCounts)n .reduce(ee.Reducer.sum(), [0]).get([0])n .divide(aCount);n var bCount = total.subtract(aCount);n var bMean = sum.subtract(aCount.multiply(aMean)).divide(bCount);n return aCount.multiply(aMean.subtract(mean).pow(2)).add(n bCount.multiply(bMean.subtract(mean).pow(2)));n});n// Return the mean value corresponding to the maximum BSS.nreturn means.sort(bss).get([-1]);n};n

3.從伺服器取出閥值，應用到原圖的 NDVI 波段上

var threshold = otsu(histogram.get(nd_histogram)); nvar water = NDVI.lt(threshold);n

自動檢測的結果，既快又好。

參考資料

1. https://en.wikipedia.org/wiki/Otsu%27s_method
2. https://www.zhihu.com/question/19829234
3. https://medium.com/google-earth/otsus-method-for-image-segmentation-f5c48f405e
4. Donchyts, G.; Baart, F.; Winsemius, H.; Gorelick, N.; Kwadijk, J.; van de Giesen, N. Earths surface water change over the past 30 years. Nature Clim. Change 6, 810–813.

推薦閱讀：