読者です 読者をやめる 読者になる 読者になる

島の海岸線の長さを推定してみる

この記事はOpenCV Advent Calendar 2015の23日目の記事です。

qiita.com

あらかじめ言います。しょぼいです。すみません!

目次

ネタ

島の海岸線の長さを計りたい

Google Mapsを見ているだけで余裕で1日過ごせる私にとっては、いくらでも調べたいことがあります。これはそのうちの1つです。皆さんも、気になりますよね?

手法

余裕がなくて全然自動化ができませんでした・・・。以降ローテクの連発です。

  • 入力は、ブラウザのスクリーンショット
  • 中央に映っている島の海岸線を抽出。
  • 長さを計算し、縮尺からkmに変換する。

ではまずは、この母なる大地北海道を入れてみましょう。

f:id:Schima:20151221225328p:plain

海領域の抽出

陸は様々な色が使われているので、青一色の海を抽出することで逆算して陸を求めます。

簡単に、RGBのBを見て判断することにします。

ちなみに、私はこれをやってみて初めて気づきました。Google Mapsの海の色は一様ではありません。左上が薄くて、だんだん濃くなっています。気付いた後は、もうなぜ今まで一様と思っていたのか意味が分からないほど、違いを感じます。

ただしいずれにせよ、RとGは見なくても"大体は"問題なく抽出できるようでした。

void extractBlueRegion(
    const cv::Mat &src, cv::Mat &dst, uchar bLow, uchar bHigh)
{
    CV_Assert(!src.empty() && src.type() == CV_8UC3);
    CV_Assert(bLow <= bHigh);

    std::vector<cv::Mat> planes;
    cv::split(src, planes);
    cv::Mat &b = planes[0];

    cv::threshold(b, dst, bLow, 255, cv::THRESH_TOZERO);
    cv::threshold(dst, dst, bHigh, 255, cv::THRESH_TOZERO_INV);
    cv::threshold(dst, dst, 0, 255, cv::THRESH_BINARY);
}
cv::Mat src = cv::imread("hokkaido.png");
    
// 海領域
cv::Mat blueRegion;
extractBlueRegion(src, blueRegion, 250, 255); // 適当!

まずまずですね。 f:id:Schima:20151221230334p:plain

輪郭抽出

お決まりの輪郭抽出です。findCountoursです。

輪郭は白い部分に対して求められますから、反転して陸地を白にしてから呼び出しましょう。

// 輪郭の型の別名(C++11)
using Contour = std::vector<cv::Point>;

// 陸に置き換える
cv::Mat land;
cv::bitwise_not(blueRegion, land);

// 最外郭のみの輪郭抽出
std::vector<Contour> contours;
cv::findContours(land, contours, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_NONE);

輪郭の描画には、drawContoursが便利です。

cv::Mat view = src.clone();
cv::drawContours(view, contours, -1, cv::Scalar(0, 0, 255), 3);

f:id:Schima:20151221231107p:plain

真ん中の輪郭のみを残す

今回は、「計測したい島は真ん中に映っている」という前提をおいていますので、一番真ん中にある輪郭を求めます。地味にコードが多い。

(これ、変な形の島だとたまにしくじります。)

double distance(cv::Point2d a, cv::Point2d b)
{
    return std::sqrt(std::pow(a.x - b.x, 2) + std::pow(a.y - b.y, 2));
}

cv::Point2d getCentroid(const Contour &contour)
{
    cv::Moments m = cv::moments(contour, false); 
    return cv::Point2d(m.m10 / m.m00, m.m01 / m.m00);
}

void orderContoursByDistanceToCenter(std::vector<Contour> &contours, const cv::Mat &image)
{
    cv::Point2d center(image.cols / 2.0, image.rows / 2.0);

    std::sort(contours.begin(), contours.end(), [center](const Contour &a, const Contour &b) -> int
    {
        double aDist = distance(center, getCentroid(a));
        double bDist = distance(center, getCentroid(b));
        return (aDist < bDist);
    });
}
// 最外郭のみの輪郭抽出
std::vector<Contour> contours;
cv::findContours(land, contours, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_NONE);

// 真ん中からの距離順に並べ替え、先頭(一番中心から近い)を取得
orderContoursByDistanceToCenter(contours, src);
const Contour &islandContour = contours[0];

f:id:Schima:20151221233046p:plain

輪郭の長さを求める

arcLengthで一発。 多分閉じた輪郭になっているので(そうなるようスクリーンショットを取りましょう)、第二引数はtrueにしています。

double perimeter = cv::arcLength(islandContour, true);

輪郭の長さをkmに変換する

上で求めたperimeterは画像におけるピクセル単位の長さですから、これを現実の距離に変換します。

ここがいちばんしょぼく、なんとかしたかったです。

Google Mapsは、右下隅に縮尺が出ています。これを手でピクセル長さを計り、ルックアップテーブルに入れておきます。

f:id:Schima:20151221232110p:plain

また、URLにzIndexが書かれており、これで今の表示の縮尺を定量的に得られます。 末尾の○○zという数値です。

f:id:Schima:20151221232418p:plain

Google Chormeでは、zは0.25単位で表示が変わるようだったので、下のコードはそれを基準で書きました。ただ、他のブラウザではそうとも限らないらしく、上の画像(Microsoft Edge)はさっそく違っていて、ダメダメ感がしています。OCRモジュールとかで読み取ってやればかっこよかったかもしれませんが・・・

double pixelToKilometer(double pixelLength, double z)
{
    // google mapsのz値(100倍)の値と、1ピクセルあたりのkm縮尺の対
    std::unordered_map<int, double> lut;
    lut[775] = 50 / 97.0; // 北海道画像での縮尺。 50kmで97px の縮尺だった。
    // ...
    // この調子でいっぱい登録しておくイメージ

    // zを100倍にし、25刻みにする (小数のキーを一応避ける)
    int zInteger = int(z);
    double zDecimal = z - zInteger;
    int key = (zInteger * 100) + (int((zDecimal + 0.125) * 4) * 100 / 4);

    if (lut.find(key) != lut.end())
        return pixelLength * lut[key];

    return -1;
}

z値は引数で与えます。これもすごいハードコーディングですが、動けばいいんです。

double kilometer = pixelToKilometer(perimeter, 7.75);
std::printf("周囲長: %.2fkm\n", kilometer);

以上、ごく初歩的な実装を書いてみました。

北海道の海岸線長

さて、結果は。

f:id:Schima:20151221233447p:plain

答え合わせ。

http://www1.kaiho.mlit.go.jp/KAN1/soudan/kyori.html

海岸線距離(km) 3,062

大敗北。500kmも違う。

敗因

しかしこれは、予想通りです。ぶっちゃけますと、ハナからだめだろうなと思いながら実装してました。

道東を見てみるだけでも、厚岸湖や風連湖のカウントがされてないなーとか、野付半島がちぎれてるなーとか、いろいろ見つけられます。逆に海抽出のしょぼさから根室や斜里の文字が輪郭に含まれてしまってもいますが、長さ挽回はできなかったのでしょう。

もっと拡大して計算すれば、より近い値が出そうです。そうすると画面に収まらないので、今回のコードでは苦しいですが。

f:id:Schima:20151221233921p:plain

佐渡島の海岸線長

だめとわかりつつ他にもやってみます。わりとなだらかそうにぱっと見は思える佐渡を入れてみましょう。

zは10.75、縮尺は5km/80pxです。 f:id:Schima:20151221234325p:plain

f:id:Schima:20151221235003p:plain

答え合わせ。

佐渡島 - Wikipedia

佐渡島(さどしま、さどがしま)[1]は、新潟県西部に位置する周囲262.7kmの島で、

またしても大敗北。まあ、だいたいは過少に出ますね。

このほか5島くらい試したのですが、惨めになるのでやめておきます。

まとめ

  • 縮尺が小さいと、島の輪郭の詳細は潰れてしまうため、正確に輪郭長は出せない。
  • 極限までズームしたのを貼り合わせるのが良さそう。
  • 意外と、海岸線長データが公開されている島があまりなく、答え合わせできないこともしばしば。
    • これの実装にかかった時間よりも、海岸線長の公式データを探しにネットをさまよっていた時間の方が長かったかもしれません。
    • 公式データがなかなか見つからないということは、これをちゃんと作れば価値があるのではないでしょうか。

結果は満足いくものではありませんが、いつも見るだけのGoogle Mapsで、GoogleAPIではなく自分の力で遊べた感が強く、楽しかったです。

面積への展開

contourAreaを使えば、輪郭の面積を求められます(輪郭内のピクセル数を直接数えても良いですね)。ということは、島の面積が求まります。

面積データは比較的入手しやすいので、答え合わせも容易でしょう。それゆえ動機に欠けますが。

コード

書き捨てですが上げておきます。 gist.github.com