二値画像からのGeoJSON生成を実装してみたがうまくいかなかったメモ

Qiitaへの投稿「Tellus画像処理の練習メモ:緑地を切り出してGeoJSONを生成してみる - Qiita」では、potracesvg2geojsonで二値画像からからGeoJSONを生成しましたが、この処理を自力でできないか週末にJavaで試してみました。
うまく生成できたらQiitaに投稿するところなのですが、現段階ではきれいなGeoJSONは生成できていません。GeoJSONをうまく生成するには、境界線のトレースだけではなく、Pathの統合や曲線化等が必要な感じです。
後学のため、メモを記録しておきます。

二値画像JSON化の処理手順

今回は、次の手順でJSON生成を試してみました。 なお、極力、他のライブラリは使用しない方針としました。

  1. 二値画像から境界線を抽出する
  2. 境界線をトレースし、Pathに変換する
  3. Pathの座標情報からGeoJSONを生成する

二値画像から境界線を抽出し、境界線をトレース

二値画像のピクセルを検査して境界線を抽出した後、以下のアルゴリズムで境界線をトレースしました。

f:id:termat:20190407221933p:plain 実装したクラスは以下のとおりです。

import java.awt.geom.GeneralPath;
import java.awt.image.BufferedImage;
import java.io.IOException;
import java.util.ArrayList;
import java.util.List;
import java.util.Stack;

public class Bitmap {
    public int[] data;
    public int[] mark;
    private int w;
    private int h;

    /**
     * BufferedImageからBitmapインスタンスを生成する
     *
     * @param is
     * @return
     * @throws IOException
     */
    public static Bitmap createBitmap(BufferedImage im){
        return new Bitmap(im.getWidth(),im.getHeight(),getImageData(im));
    }

    /*
     * 二値のピクセルデータを生成する
     */
    private static int[] getImageData(BufferedImage im){
        int w=im.getWidth();
        int h=im.getHeight();
        int[] bimg=new int[w*h];
        int it=0;
        for(int i=0;i<h;i++){
            for(int j=0;j<w;j++){
                int cc=im.getRGB(j, i);
                int r=(cc & 0x000000FF);
                int g=((cc & 0x0000FF00) >> 8);
                int b=((cc & 0x00FF0000) >> 16);
                int col=(r+g+b)/3;
                bimg[it++]=(col < 128 ? 1 : 0);
            }
        }
        return bimg;
    }

    /*
     * コンストラクタ
     */
    private Bitmap(int w,int h,int[] d){
        this.w=w;
        this.h=h;
        data=inspection(d);
        mark=new int[data.length];
    }

    /*
     * ピクセル(x,y)の値を取得
     */
    private int get(int x,int y){
        return data[w*y+x];
    }

    /*
     * ピクセル(x,y)に値を入力
     */
    private void set(int x,int y,int v){
        data[w*y+x]=v;
    }

    /*
     * ピクセル(x,y)の検査回数を取得
     */
    private int num_mark(int x,int y){
        return mark[w*y+x];
    }

    /*
     * ピクセル(x,y)の検査回数を1回増やす
     */
    private void mark(int x, int y){
        mark[w*y+x]++;
    }

    /*
     * 配列番号iをピクセル位置(x,y)に変換する
     */
    private Vec index(int i){
        Vec p = new Vec();
        p.y = (int)Math.floor(i/w);
        p.x = i - p.y * w;
        return p;
    }

    /*
     * 画像を検査し、境界線を抽出する
     */
    private int[] inspection(int[] d){
        int[] ret=new int[d.length];
        for(int x=1;x<w-1;x++){
            for(int y=1;y<h-1;y++){
                if(d[w*y+x]==0)continue;
                boolean flg=false;
                for(int i=x-1;i<=x+1;i++){
                    for(int j=y-1;j<=y+1;j++){
                        if(d[w*j+i]==0){
                            flg=true;
                            break;
                        }
                    }
                    if(flg)break;
                }
                if(flg&&d[w*y+x]==1)ret[w*y+x]=1;
            }
        }
        return ret;
    }

    /*
     * 境界線をトレースした座標を返す
     */
    private Vec trace(Vec p){
        int x=p.x;
        int y=p.y;
        for(int i=x-1;i<=x+1;i++){
            for(int j=y-1;j<=y+1;j++){
                if(i==x&&j==y)continue;
                if(get(i,j)==1&&num_mark(i,j)<2){
                    mark(i,j);
                    return new Vec(i,j);
                }
            }
        }
        return null;
    }

    /*
     * 境界トレースの始点を検索する
     */
    private Vec findNext(Vec p){
        int i=w*p.y+p.x;
        while(i<data.length&&data[i]==0)i++;
        return i< w*h ? index(i) : null;
    }

    /*
     * 境界をトレースしたGeneralPathを生成する
     */
    private GeneralPath getPath(Vec p){
        List<Vec> list=new ArrayList<Vec>();
        Vec cursor=new Vec(p.x,p.y);
        mark(p.x,p.y);
        while((cursor=trace(cursor))!=null){
            list.add(cursor);
        }
        list=reduction(list);
        int n=list.size();
        if(n>2){
            GeneralPath gp=new GeneralPath();
            gp.moveTo(p.x, p.y);
            set(p.x,p.y,0);
            for(Vec c : list){
                gp.lineTo(c.x, c.y);
                if(c.x==0&&c.y==0)System.out.println("113");
                set(c.x,c.y,0);
            }
            if(p.distSq(list.get(list.size()-1))<=2)gp.closePath();
            return gp;
        }else{
            set(p.x,p.y,0);
            return null;
        }
    }

    /*
     * トレースした境界の座標を縮減する
     */
    private List<Vec> reduction(List<Vec> list){
        Stack<Vec> stack=new Stack<Vec>();
        Stack<Vec> tmp=new Stack<Vec>();
        stack.addAll(list);
        while(stack.size()>=3){
            Vec v1=stack.pop();
            Vec v2=stack.pop();
            Vec v3=stack.pop();
            if(isLinear(v1,v2,v3)){
                stack.push(v3);
                stack.push(v1);
            }else{
                tmp.push(v1);
                stack.push(v3);
                stack.push(v2);
            }
        }
        while(!stack.isEmpty())tmp.push(stack.pop());
        List<Vec> ret=new ArrayList<Vec>();
        while(!tmp.isEmpty())ret.add(tmp.pop());
        return ret;
    }

    /*
     * 3つの座標が直線上にあるか否かを判定する
     */
    private boolean isLinear(Vec v1,Vec v2,Vec v3){
        double x1=v2.x-v1.x;
        double y1=v2.y-v1.y;
        double l1=Math.sqrt(x1*x1+y1*y1);
        x1=x1/l1;
        y1=y1/l1;
        double x2=v3.x-v1.x;
        double y2=v3.y-v1.y;
        double l2=Math.sqrt(x2*x2+y2*y2);
        x2=x2/l2;
        y2=y2/l2;
        return ((Math.abs(x1-x2)<1e-8)&&(Math.abs(y1-y2)<1e-8));
    }

    /**
     * 画像の境界線をトレースし、GeneralPathのリストを取得する
     *
     * @return
     */
    public List<GeneralPath> getPathList(){
        List<GeneralPath> ret=new ArrayList<GeneralPath>();
        Vec p=new Vec();
        while((p=findNext(p))!=null){
            GeneralPath gp=getPath(p);
            if(gp!=null)ret.add(gp);
        }
        return ret;
    }
}
public class Vec {
    public int x;
    public int y;

    public Vec(){
        x=0;
        y=0;
    }

    public Vec(int x,int y){
        this.x=x;
        this.y=y;
    }

    public double distSq(Vec v){
        return (x-v.x)*(x-v.x)+(y-v.y)*(y-v.y);
    }
}

トレースした境界線からGeoJSON生成

 GeneralPathからGeometryに変換し、GeoJSONを生成するため、以下のユーティリティクラスを作成しました。

import java.awt.geom.AffineTransform;
import java.awt.geom.GeneralPath;
import java.awt.geom.PathIterator;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

import com.google.gson.Gson;
import com.google.gson.GsonBuilder;

public class Geotrace {

    /**
     * GeneralPathのリストからGeoJSONを生成
     *
     * @param list:List<GeneralPath>
     * @param af: AffineTransform
     * @return JSON
     */
    public static String createJSON(List<GeneralPath> list,
        AffineTransform af){
        List<Map<String,Object>> fes=new ArrayList<Map<String,Object>>();
        for(GeneralPath gp : list){
            boolean poly=false;
            List<double[]> tmp=new ArrayList<double[]>();
            PathIterator it=gp.getPathIterator(af);
            while(!it.isDone()){
                double[] co=new double[2];
                int ck=it.currentSegment(co);
                switch(ck){
                case PathIterator.SEG_MOVETO:
                case PathIterator.SEG_LINETO:
                    tmp.add(co);
                    break;
                case PathIterator.SEG_CLOSE:
                    poly=true;
                    break;
                }
                it.next();
            }
            Map<String,Object> fea=createFeature(
                createGeometry(tmp,poly),new HashMap<String,Object>());
            fes.add(fea);
        }
        Map<String,Object> col=createFeatureCollection(fes);
        Gson gson=new GsonBuilder().setPrettyPrinting().create();
        return gson.toJson(col);
    }

    /* FeatureCollection */
    private static Map<String,Object> createFeatureCollection(
        List<Map<String,Object>> list){
        Map<String,Object> ret=new HashMap<String,Object>();
        ret.put("type", "FeatureCollection");
        ret.put("features", list);
        return ret;
    }

    /* Feature */
    private static Map<String,Object> createFeature(
        Map<String,Object> geo,Map<String,Object> prop){
        Map<String,Object> ret=new HashMap<String,Object>();
        ret.put("type", "Feature");
        ret.put("geometry", geo);
        ret.put("properties", prop);
        return ret;
    }

    /* Geometry */
    private static Map<String,Object> createGeometry(
        List<double[]> coord,boolean isPoly){
        Map<String,Object> ret=new HashMap<String,Object>();
        if(isPoly){
            ret.put("type", "Polygon");
        }else{
            ret.put("type", "LineString");
        }
        ret.put("coordinates", coord);
        return ret;
    }

    /**
     * ワールドファイルのパラメータからAffineTransformを生成
     *
     * @param param:Map<String,Double>
     * @return
     */
    public static AffineTransform createTransformByWorldfile(
        Map<String,Double> param){
        double m00=param.get("m00");
        double m01=param.get("m01");
        double m02=param.get("m02");
        double m10=param.get("m10");
        double m11=param.get("m11");
        double m12=param.get("m12");
        return new AffineTransform(m00,m10,m01,m11,m02,m12);
    }

    /**
     * 矩形領域の四隅のピクセル座標、緯度経度を指定してAffineTransformを生成
     *
     * @param param:Map<String,Double>
     * @return
     */
    public static AffineTransform createTransformByRect(
        Map<String,Double> param){
        double min_x=param.get("min_x");
        double min_y=param.get("min_y");
        double max_x=param.get("max_x");
        double max_y=param.get("max_y");
        double min_lon=param.get("min_lon");
        double min_lat=param.get("min_lat");
        double max_lon=param.get("max_lon");
        double max_lat=param.get("max_lat");
        double m00=(max_lon-min_lon)/(max_x-min_x);
        double m01=0;
        double m02=min_lon;
        double m10=0;
        double m11=(min_lat-max_lat)/(max_y-min_y);
        double m12=max_lat;
        return new AffineTransform(m00,m10,m01,m11,m02,m12);
    }
}

GeoJSON生成結果

GeoJSONの生成結果を以下に示します。 境界線を抽出した段階では行けそうな感じですが、GeoJSONにして拡大したところ、かなりダメな感じです。

■元座標
■境界線抽出(GeneralPath)
■GeoJSON
■GeoJSON(拡大図)

修正が必要な部分

使えるものにするには、以下の処理が必要そうです。 このまま開発を続けるか、potraceとsvg2geojsonを移植・統合するか、さて、どうしようか考え中です。

  • 抽出したGeneralPathを結合・統合する必要がある
  • ジクザクな境界線を一度曲線化して、その後、座標を再配置する必要がある