Hi,
I have question regarding proper use of Polygonizer class (I'm using NTS, but it shouldn't matter I hope) I'm working on utility to create convex hull of a point set. Result of first step of that utility is whole bunch of line strings representing convex hull or point set. Then I'm using Polygonizer to transform list of line strings to polygon. Problem is that sometime result of polygonizer is not one 'Big' outlining polygon (having potentially some holes in it), but instead lot of small polygons. I suspect that problem is when some line segments lie on the boundary and as well form smaller hole in polygon that somehow create probelm. It's a bit difficult to describe, but I'll try to simplify it as much as possible. Imagine chess board (8x8 square) having line segments of lenght 1 all arround. If it happens that there are line segments arroung some of little squares arround the edges (e.g. field C1) then instead or returning big polygon (8x8) polygonizer return only small one (1x1 C1 square). My guess is that line segment on the edge get assigned to that small polygon, and then rest of the line segments do not form closed ring. That is of course totaly different of what I need, as I want to get back largest possible polygon, potentially discarding all those smaller but it seem that they take precedence and mess up desired output. Is there any way to control how polygonizer work so that i get my big polygon back. Thanks Dragan Blagojevic
Senior Developer
GeoDigital International Inc. 775 Topaz Avenue, Unit 104 Victoria, BC V8T 4Z7 Tel: 250.388.0500 |
Hi Dragan,
Nice question :-) > Problem is that sometime result of polygonizer is not one 'Big' outlining polygon (having potentially some holes in it), but instead lot of small polygons. > When you get back a collection fo small polys instead of one large one, will unioning the polys give you the desired result ? Also I don't quite understand about the potential for holes, are you trying to generate a variation on the general convex hull ? Michael _______________________________________________ jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
Hi Michael,
What I'm trying to get is Concave hull of set of points. Result of that particular concave hull algorithm (Alpha Shapes) is bunch of line segments. That part is fine, and I get bunch of line segments which outline my points just fine, however they are line segments, not polygon which is what I need. When you visually inspect those line segments, you can see big outlining polygon, lot of holes in it, lot of dangling lines, cut lines and what not. Generally polygonizer munch that lines and create polygon with holes just fine. However there is quite often situation that along the border of that (not yet formed) big polygon group of line segments form a small polygon which would in ideal case be interpreted as a hole in big polygon. Problem generally arises if one or more of line segments is shared by big and small polygon. I don't know how polygonizer works, but it would appear that those small holes get somehow created first, and line segment that belong both to hole and big polygon get 'stolen' by that small polygon, which leave big outlining polygon without one of line segments. That in turn make it not closed and basically it get thrown away. Now, what I would like to get is exact opposite, I wouldn't mind to have some small inside discarded as long as I get big outline. I hope that I explained problem better this time. If you wish I can send you .shp files with line segments and resulting polygons where is much easier to see what the problem is. Another question, is there any good user guide for JTS / NTS? I didn't find anything except class reference and some very basic user guide. Regards Dragan
Senior Developer
GeoDigital International Inc. 775 Topaz Avenue, Unit 104 Victoria, BC V8T 4Z7 Tel: 250.388.0500 |
Just realized that in my first post I said by mistake that I want to get convex hull of a point set.
I'm actually after concave hull Regards Dragan
Senior Developer
GeoDigital International Inc. 775 Topaz Avenue, Unit 104 Victoria, BC V8T 4Z7 Tel: 250.388.0500 |
2009/4/2 Dragan Blagojevic wrote:
> > Just realized that in my first post I said by mistake that I want to get convex hull of a point set. > > I'm actually after concave hull OK - it all makes a lot more sense now :-) I'm actually very interested in finding a concave hull implementation myself. I contacted the authors of this paper recently... http://repositorium.sdum.uminho.pt/handle/1822/6429?locale=en ...and asked for a copy. Unfortunately, they told me that they are patenting their algorithm (not just their code) and any use of it will require commercial licensing. Given that, I decided that ignorance was legally safer and it was best not to take up their offer to send me the paper ! But back to your problem... yes, if you could send me an example shapefile that would be great. I'm not one of the real gurus here so I can't promise to find the solution. But perhaps with help from others we will work it out. cheers Michael _______________________________________________ jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
Yes, there is quite a lot of papers mentioning concave hull but it apear that it's kind of secret :-) It took me a while to find description of Alpha Shape algorithm, I believe that is basically the same thing as one in the link you mentioned above, main difference being that this one was published open. I'll send you a paper with algorithm description I found tomorrow when I get to the office, together with shape files. If you wait a while, I'll send you code of my implementation as well - at the moment code is in total mess as I was experimenting with various options and optimizations, so it won't do you any good in it's current state. It's in C#, but it will be trivial to change it to Java. Regards Dragan
Senior Developer
GeoDigital International Inc. 775 Topaz Avenue, Unit 104 Victoria, BC V8T 4Z7 Tel: 250.388.0500 |
In reply to this post by Dragan Blagojevic
Two possibilities:
1. bug in Polygonizer 2. your input is not fully noded. I would suspect #2. Can you post sample data? (to the list or to me directly). I'd like to see the concave hull paper you mention as well - I'm quite interested in an implementation of this. Although as pointed out there is no single interpretation of how concave hull should be defined - so there could be several ways to compute something reasonable) As for a user guide to JTS, I realize this would be very nice to have. I'd very much like to provide something for this. Unfortunately it's quite time-consuming to create, and it's been hard for me to find the time to put to it. For the time being I try and keep enhancing the Javadoc whenever it looks like it could be improved. Martin Dragan Blagojevic wrote: > Hi, > > I have question regarding proper use of Polygonizer class (I'm using NTS, but it shouldn't matter I hope) > > I'm working on utility to create convex hull of a point set. > > Result of first step of that utility is whole bunch of line strings representing convex hull or point set. > > Then I'm using Polygonizer to transform list of line strings to polygon. > > Problem is that sometime result of polygonizer is not one 'Big' outlining polygon (having potentially some holes in it), but instead lot of small polygons. > > I suspect that problem is when some line segments lie on the boundary and as well form smaller hole in polygon that somehow create probelm. > > It's a bit difficult to describe, but I'll try to simplify it as much as possible. > > Imagine chess board (8x8 square) having line segments of lenght 1 all arround. > If it happens that there are line segments arroung some of little squares arround the edges (e.g. field C1) then instead or returning big polygon (8x8) polygonizer return only small one (1x1 C1 square). > > My guess is that line segment on the edge get assigned to that small polygon, and then rest of the line segments do not form closed ring. > > That is of course totaly different of what I need, as I want to get back largest possible polygon, potentially discarding all those smaller but it seem that they take precedence and mess up desired output. > > Is there any way to control how polygonizer work so that i get my big polygon back. > > > Thanks > Dragan Blagojevic > > > -- Martin Davis Senior Technical Architect Refractions Research, Inc. (250) 383-3022 _______________________________________________ jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
This post has NOT been accepted by the mailing list yet.
Hi Martin, Paper describing algorithm can be found here: http://www.isprs.org/congresses/beijing2008/proceedings/3b_pdf/37.pdf Input set is guaranteed to be fully noded, it's a set of simple line segments , 2 poins each and they do not intersect. Those line segments are result of runing alpha shape algorithm on set of LIDAR points. Input shape file and resulting polygons are attached. You can see that few small polygons arround the edge get formed first and then 'big' one don't have those line segments to close contour - or at least that's my guess work on how this happens as I don't have a clue how polygonizer does his work. I tried runing LineMerger first before polygonizer, to no effect. Regards Dragan shapes.zip
Senior Developer
GeoDigital International Inc. 775 Topaz Avenue, Unit 104 Victoria, BC V8T 4Z7 Tel: 250.388.0500 |
In reply to this post by Martin Davis
Well done Dragan - you've got Martin interested :-) I suspect he will
find a solution slightly faster than I will ! Michael 2009/4/4 Martin Davis <[hidden email]>: > Two possibilities: > > 1. bug in Polygonizer > 2. your input is not fully noded. > > I would suspect #2. Can you post sample data? (to the list or to me > directly). > I'd like to see the concave hull paper you mention as well - I'm quite > interested in an implementation of this. Although as pointed out there is > no single interpretation of how concave hull should be defined - so there > could be several ways to compute something reasonable) > jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
Hi Dragan and everyone,
I've written some code (see below) that demonstrates the concave hull algorithm described in the Shen Wei paper and I think I now have a better idea of what can go wrong with it. The paper is a bit naughty really :) None of the point clouds depicted in the figures are very testing for the algorithm. In my demo I generate a sparser cloud within an irregular polygon to give the algorithm a bit more challenge. If you run the program a few times you will see that while it does a generally good job of defining a hull, only a small number of replicate runs is required to find at least one solution with a degenerate polygon. I suspect this is what is happening with your data Dragan. Still, the simplicity of the algorithm is really appealing and there are various tweaks that could be tried. One obvious one is to make the alpha parameter locally adaptive. At present, the algorithm will fail if there is any point further than 2*alpha from its nearest neighbour, but if the alpha parameter was treated as a minimum, and allowed to increase when necessary, this problem would be avoided and the algorithm could be used to generate a more detailed hull for those regions of the point cloud that supported it. Michael package mxb.concavehull; import com.vividsolutions.jts.geom.Coordinate; import com.vividsolutions.jts.geom.Envelope; import com.vividsolutions.jts.geom.Geometry; import com.vividsolutions.jts.geom.GeometryFactory; import com.vividsolutions.jts.geom.LineString; import com.vividsolutions.jts.geom.Point; import com.vividsolutions.jts.geom.Polygon; import com.vividsolutions.jts.index.ItemVisitor; import com.vividsolutions.jts.index.strtree.STRtree; import java.awt.Color; import java.awt.Graphics; import java.awt.Graphics2D; import java.awt.Insets; import java.util.ArrayList; import java.util.Arrays; import java.util.List; import java.util.Random; import javax.swing.JFrame; import javax.swing.JPanel; public class ConcaveHullBuilder { private GeometryFactory gf = new GeometryFactory(); private Random rand = new Random(); private double alpha = 50d; public static void main(String[] args) { ConcaveHullBuilder hull = new ConcaveHullBuilder(); hull.demo(); } private void demo() { /** * Make an irregular polygon and generate random points within it * to test the concave hull algorithm */ final int numPoints = 200; Coordinate[] vertices = { new Coordinate(0, 0), new Coordinate(100, 100), new Coordinate(0, 250), new Coordinate(200, 300), new Coordinate(200, 450), new Coordinate(350, 450), new Coordinate(350, 200), new Coordinate(250, 200), new Coordinate(350, 100), new Coordinate(0, 0) }; Polygon poly = gf.createPolygon(gf.createLinearRing(vertices), null); Envelope env = poly.getEnvelopeInternal(); int n = 0; Point[] points = new Point[numPoints]; while (n < numPoints) { Coordinate c = new Coordinate(); c.x = rand.nextDouble() * env.getWidth() + env.getMinX(); c.y = rand.nextDouble() * env.getHeight() + env.getMinY(); Point p = gf.createPoint(c); if (poly.contains(p)) { points[n++] = p; } } List<LineString> edges = getConcaveHull(points, alpha); display((int)(env.getWidth()*1.1), (int)(env.getHeight()*1.1), poly, Color.GRAY, points, Color.RED, edges, Color.RED); } /** * Identify a concave hull using the simple alpha-shape algorithm * described in: * <blockquote> * Shen Wei (2003) Building boundary extraction based on LIDAR point clouds data. * The International Archives of the Photogrammetry, Remote Sensing and Spatial * Information Sciences. Vol. XXXVII. Part B3b. Beijing 2008 * </blockquote> * @param points the point cloud * @param alpha the single parameter for the algorithm * @return a List of LineString boundary segments for the concave hull */ private List<LineString> getConcaveHull(Point[] points, double alpha) { double alpha2 = 2 * alpha; /* * Index the points for faster processing */ STRtree index = new STRtree(); for (Point p : points) { index.insert(p.getEnvelopeInternal(), p); } index.build(); List<LineString> edges = new ArrayList<LineString>(); List<Point> candidates = new ArrayList<Point>(); candidates.addAll(Arrays.asList(points)); while (!candidates.isEmpty()) { Point p1 = candidates.remove(0); Envelope qEnv = new Envelope(p1.getX() - alpha2, p1.getX() + alpha2, p1.getY() - alpha2, p1.getY() + alpha2); PointVisitor visitor = new PointVisitor(p1, alpha2); index.query(qEnv, visitor); if (visitor.plist.size() < 2) { System.out.println("isolated point"); break; } visitor.plist.remove(p1); boolean[] used = new boolean[visitor.plist.size()]; Arrays.fill(used, false); int numPts = visitor.plist.size(); while (numPts > 0) { Point p2; while (true) { int pindex = rand.nextInt(visitor.plist.size()); if (!used[pindex]) { p2 = visitor.plist.get(pindex); used[pindex] = true; numPts--; break; } } Point pcentre = createCircle(p1, p2, alpha); boolean onBoundary = true; for (Point vp : visitor.plist) { if (vp != p2) { if (pcentre.distance(vp) <= alpha) { onBoundary = false; break; } } } if (onBoundary) { edges.add(gf.createLineString(new Coordinate[]{p1.getCoordinate(), p2.getCoordinate()})); } } } return edges; } /** * Calculate the centre coordinates of a circle of radius alpha that has * point coordinates c1 and c2 on its circumference. * * @param c1 first circumference point * @param c2 second circumference point * @param alpha radius * @return a Coordinate representing the circle centre */ private Point createCircle(Point p1, Point p2, double alpha) { Coordinate centre = new Coordinate(); double dx = (p2.getX() - p1.getX()); double dy = (p2.getY() - p1.getY()); double s2 = dx * dx + dy * dy; double h = Math.sqrt(alpha * alpha / s2 - 0.25d); centre.x = p1.getX() + dx / 2 + h * dy; centre.y = p1.getY() + dy / 2 + h * (p1.getX() - p2.getX()); return gf.createPoint(centre); } /** * Display demo results * @param poly the polygon used to generate the point cloud * @param polyCol display colour for the polygon * @param points the point cloud * @param pointCol display colour for the points * @param edges concave hull boundary segments as a List of LineStrings * @param edgeCol display colour for the segments */ private void display(int w, int h, final Polygon poly, final Color polyCol, final Point[] points, final Color pointCol, final List<LineString> edges, final Color edgeCol) { JFrame frame = new JFrame("Concave hull demo"); JPanel panel = new JPanel() { final int ow = 2; final int ow2 = 4; @Override protected void paintComponent(Graphics g) { super.paintComponent(g); Graphics2D g2 = (Graphics2D) g; g2.setColor(polyCol); LineString ring = poly.getExteriorRing(); Coordinate clast = ring.getCoordinateN(0); for (int i = 1; i < ring.getNumPoints(); i++) { Coordinate c = ring.getCoordinateN(i); g2.drawLine((int) clast.x, (int) clast.y, (int) c.x, (int) c.y); clast = c; } g2.setColor(pointCol); for (int i = 0; i < points.length; i++) { int x = (int) Math.round(points[i].getX()); int y = (int) Math.round(points[i].getY()); g2.fillOval(x - ow, y - ow, ow2, ow2); } g2.setColor(edgeCol); Geometry polyEnv = poly.getEnvelope(); for (LineString l : edges) { Point p0 = l.getStartPoint(); Point p1 = l.getEndPoint(); int x0 = (int) l.getStartPoint().getX(); int y0 = (int) l.getStartPoint().getY(); int x1 = (int) l.getEndPoint().getX(); int y1 = (int) l.getEndPoint().getY(); g2.drawLine(x0, y0, x1, y1); } } }; panel.setBackground(Color.WHITE); frame.getContentPane().add(panel); Insets insets = frame.getInsets(); frame.setSize(w + insets.left + insets.right, h + insets.bottom + insets.top); frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE); frame.setVisible(true); } } class PointVisitor implements ItemVisitor { public List<Point> plist = new ArrayList<Point>(); private double maxDist; private Point refP; PointVisitor(Point refP, double maxDist) { this.refP = refP; this.maxDist = maxDist; } public void visitItem(Object o) { if (refP.isWithinDistance((Point) o, maxDist)) { plist.add((Point) o); } } } _______________________________________________ jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
This post has NOT been accepted by the mailing list yet.
Hi Michael, Limitations of the algorithm are obvious, but I don't have problem with that. If you however come up with some idea how to treat alpha as minimum and to have it increase if needed depending on local conditions I will be very interested. For my particular purpose loss of few points on the edges is quite irrelevant. We are doing aerial LIDAR surveys for utility companies, and task at hand is to find outline of LIDAR points cloud and compare it with predefined area we are supposed to cover. Due to enormous amount of data I'm already discarding quite lot of them before I run algorithm - we get about 1 Billion (with B) points per day of flying. My problem however is not when alpha shape do not return good enough list of edges to represent polygon, but when it does as in sample shapes I uploaded. When there is a 'big' polygon that can be made, but with some edges on the outline of it forming smaller polygons, polygonizer sometime fails to make big one that is most important for my particular purpose. If polygonizer can somehow be made to do the exact oposite of example I uploaded I would quite happy... Regards Dragan
Senior Developer
GeoDigital International Inc. 775 Topaz Avenue, Unit 104 Victoria, BC V8T 4Z7 Tel: 250.388.0500 |
Free forum by Nabble | Edit this page |