二値画像からのGeoJSON生成を実装してみたがうまくいかなかったメモ
Qiitaへの投稿「Tellus画像処理の練習メモ:緑地を切り出してGeoJSONを生成してみる - Qiita」では、potraceとsvg2geojsonで二値画像からからGeoJSONを生成しましたが、この処理を自力でできないか週末にJavaで試してみました。
うまく生成できたらQiitaに投稿するところなのですが、現段階ではきれいなGeoJSONは生成できていません。GeoJSONをうまく生成するには、境界線のトレースだけではなく、Pathの統合や曲線化等が必要な感じです。
後学のため、メモを記録しておきます。
二値画像JSON化の処理手順
今回は、次の手順でJSON生成を試してみました。 なお、極力、他のライブラリは使用しない方針としました。
- 二値画像から境界線を抽出する
- 境界線をトレースし、Pathに変換する
- Pathの座標情報からGeoJSONを生成する
二値画像から境界線を抽出し、境界線をトレース
二値画像のピクセルを検査して境界線を抽出した後、以下のアルゴリズムで境界線をトレースしました。
実装したクラスは以下のとおりです。
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を結合・統合する必要がある
- ジクザクな境界線を一度曲線化して、その後、座標を再配置する必要がある