GEEでOtsuの二値化 - 海老ノート

できるはずだけど面倒くさそうと思っていた

画像処理でおなじみ大津の二値化。考え方はシンプルだし、GoogleEarthEngineでできるんだろうけどコード書くの面倒だなと。以前八郎潟の湛水の有無を見るのにSARデータ使った時にも、二値化したいと思っていたけど適当にごまかした。が、先日検索してみると随分前にコードが紹介されてました、しかも公式Nick先生。というわけで、こちらをお手本に試します。

medium.com

お手本コードのメモ

このコードのポイントは ee.Reducer.histogram()を使うところ。お手本では2~256クラスに分けるよう引数で指定して、ヒストグラムの各階級ごとに平均とデータ数をまとめてBSSと呼んでいる分割の良さを表す指標を計算している。

これまで面倒と思ってたのは長いリストを作って処理しないといけないところ。0~255のDNくらいだとそこまで大変ではないだろうけど、NDWIを0.001刻みでとなると処理に時間かかりそうだなと。そこのところ、このお手本ではhistogramの関数をうまく使うことでざっくりと最大256段階に分けている。

お手本へのリンクhttps://code.earthengine.google.com/6f3abb73b6642198485e36000024827a

もっと詳細にみたい

なるほどと思いつつ、手を加えたくなっていじってみる。Nick先生もざっくりだよ、とおっしゃっているし：

We’re not going to search exhaustively. We’re just going to search over the thresholds that are represented by the bins in a histogram. The advantage of that approach is that it only requires a single pass over the data.

それなら、ざっくりから細かくしてみようかと。蛇足だったけど。

力ずくのコード

単純にお手本コードのee.Reducer.histogram()の一つ目の引数を大きくすれば閾値のレベルは細かくなる。ただ、それじゃあやった気がしないというわけで、お手本で求められた閾値の前後を細かく区切ってBSSがどうなるかを調べることにします。リストにした閾値候補ごとにmaskを作って、地道にBSSを計算してみる。大体のつくりはお手本コードをお借りしてます。nClassを小さくして、計算を数回繰り返せば早くて細かい計算ができるはず。

var nClass = 30 // 閾値候補の数。30くらいなら我慢のいらない速さで回る。

var bucketWidth = ee.Dictionary(histogram.get('B5_histogram')).get('bucketWidth');
var thresholdList = ee.List.sequence({start: threshold.subtract(bucketWidth), 
                                      end: threshold.add(bucketWidth), 
                                      count: nClass}); //閾値候補のリスト
print(thresholdList);

//今回メインの関数
var detailedOtsu = function(image){
  var globalMean = image.reduceRegion({
                    reducer: ee.Reducer.mean(),
                    geometry: polygon,
                    scale: 30,
                    bestEffort: true }).values().getNumber(0) // 全体の平均

 var funcBss = thresholdList.map(function(i){
  var greaterThanMask = image.gt(ee.Image.constant(i))
  var greaterThan = image.mask(greaterThanMask)  //閾値候補超のピクセル抽出
  var lessThanEq = image.mask(greaterThanMask.eq(0)) //閾値候補以下のピクセル抽出

  var funcStats = function(image){
    var intraMean = image.reduceRegion({
                     reducer:ee.Reducer.mean(),
                     geometry:polygon,
                     scale:30,
                     bestEffort: true }).values().getNumber(0) //マスクしたピクセルの平均
    var intraCount = image.reduceRegion({
                     reducer: ee.Reducer.count(),
                     geometry: polygon,
                     scale: 30,
                     bestEffort: true }).values().getNumber(0) //マスクしたピクセル数
  return intraCount.multiply(intraMean.subtract(globalMean).pow(2))}
    
 return funcStats(greaterThan).add(funcStats(lessThanEq))}
 )
 print(ui.Chart.array.values(ee.Array(funcBss), 0, thresholdList)) //BSSのグラフ作成
 
return ee.Array(thresholdList).sort(funcBss).get([-1]) //BSS最大の閾値候補の選択
}

どっちでもいい感じ

お手本コードでの閾値は8127.50。一方細かく見た方では8184.88。グラフだとこういう感じ。

閾値が60違うとどう変わるのか？

あまり違わない。画像で見てみるとほぼ同じ。

左：元画像（LS8 False Colour）、中：お手本の閾値、右：詳細な閾値

結局水と陸が混ざっているピクセルの扱いの問題みたいから、もっと解像度が高い画像か現地データがないと精度の評価は難しそう。

k-meansでもできるのでは？

GEEでも使えるお手軽分類、k-meansを2クラスでやってもできるんじゃないかと思ってやってみたら、これもほぼ同じ結果になった。

左：元画像（LS8 False Colour）、中：お手本の閾値、右：2クラスのk-means

k-meansがいいのはGEEで書くコードが断然短くて楽。悪いのはどちらが水でどちらが陸のクラスなのかはランダムに決まるから、時系列で多くの画像から水域を抽出する処理をするような場合はひと手間が必要になるはずというところか。

まとめ：どこまで頑張ってどこで妥協するか

そういうことです。Landsat8の画像でざっと見る限りではどっちでも、というところ。値を取得する範囲となるポリゴンの位置によっても結構変わるだろうし。より詳細に突き詰めれば細かくなるんだろうけど、それがどの程度意味あるかっていうのはわからない、正解データがあれば数値で評価できるんだろうけど。バランスは大事です。

今回のコードはこちら： https://code.earthengine.google.com/b95c53a32bca944425f0f7d3b08e2f5d