Gaussian smoothing and new geometry densifier

modules/core/src/main/java/org/locationtech/jts/densify/LittleThumblingDensifier.java

@@ -0,0 +1,154 @@
+/*
+ * Copyright (c) 2016 Vivid Solutions.
+ *
+ * All rights reserved. This program and the accompanying materials
+ * are made available under the terms of the Eclipse Public License v1.0
+ * and Eclipse Distribution License v. 1.0 which accompanies this distribution.
+ * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v10.html
+ * and the Eclipse Distribution License is available at
+ *
+ * http://www.eclipse.org/org/documents/edl-v10.php.
+ */
+package org.locationtech.jts.densify;
+
+import org.locationtech.jts.geom.Coordinate;
+import org.locationtech.jts.geom.CoordinateSequence;
+import org.locationtech.jts.geom.Geometry;
+import org.locationtech.jts.geom.GeometryFactory;
+import org.locationtech.jts.geom.LineString;
+import org.locationtech.jts.geom.MultiPolygon;
+import org.locationtech.jts.geom.Polygon;
+import org.locationtech.jts.geom.PrecisionModel;
+import org.locationtech.jts.geom.util.GeometryTransformer;
+
+/**
+ * Densify a geometry using the "Little Thumbling" strategy: Add a vertex for each step along the line.
+ * 
+ * @author julien Gaffuri
+ *
+ */
+public class LittleThumblingDensifier {
+
+	/**
+	 * Densify a geometry: Walk along the line step by step and record the positions at each step
+	 * (like the "Little Thumbling" would do).
+	 * Return the line composed of the recorded positions.
+	 * The vertices of the input geometry are not kept, which can result in an output line with shorter length.
+	 * The result is different from org.locationtech.jts.densify.Densifier
+	 * 
+	 * @param geom
+	 * @param stepLength the step length
+	 * @return the densified geometry
+	 */
+	public static Geometry densify(Geometry geom, double stepLength) {
+		LittleThumblingDensifier densifier = new LittleThumblingDensifier(geom);
+		densifier.setStepLength(stepLength);
+		return densifier.getResultGeometry();
+	}
+
+	/**
+	 * Densifies a coordinate sequence.
+	 * 
+	 * @param pts
+	 * @param stepLength the step length
+	 * @return the densified coordinate sequence
+	 */
+	private static Coordinate[] densifyPoints(Coordinate[] pts, double stepLength, PrecisionModel precModel) {
+		//out coords
+		LineString line = new GeometryFactory(precModel).createLineString(pts);
+		int nb = (int) (line.getLength()/stepLength);
+		Coordinate[] out = new Coordinate[nb+1];
+
+		double d=0.0, a=0.0, dTot;
+		int densIndex=0;
+		for(int i=0; i<pts.length-1; i++) {
+			Coordinate c0=pts[i], c1=pts[i+1];
+			dTot = c0.distance(c1);
+			if (d<=dTot) a = Math.atan2( c1.y-c0.y, c1.x-c0.x );
+			while(d <= dTot){
+				//use LineSegment.pointAlong instead ?
+				Coordinate c = new Coordinate(c0.x+d*Math.cos(a), c0.y+d*Math.sin(a));
+				precModel.makePrecise(c);
+				out[densIndex] = c;
+				densIndex++;
+				d+=stepLength;
+			}
+			d-=dTot;
+		}
+		out[nb] = pts[pts.length-1];
+		return out;
+	}
+
+
+	private Geometry inputGeom;
+	private double stepLength;
+
+	public LittleThumblingDensifier(Geometry inputGeom) {
+		this.inputGeom = inputGeom;
+	}
+
+	public void setStepLength(double stepLength) {
+		if (stepLength <= 0.0)
+			throw new IllegalArgumentException("Step length must be positive");
+		this.stepLength = stepLength;
+	}
+
+	/**
+	 * Gets the densified geometry.
+	 * 
+	 * @return the densified geometry
+	 */
+	public Geometry getResultGeometry() {
+		return (new LittleThumblingDensifyTransformer(stepLength)).transform(inputGeom);
+	}
+
+
+	static class LittleThumblingDensifyTransformer extends GeometryTransformer {
+		double stepLength;
+
+		LittleThumblingDensifyTransformer(double stepLength) {
+			this.stepLength = stepLength;
+		}
+
+		protected CoordinateSequence transformCoordinates(CoordinateSequence coords, Geometry parent) {
+			Coordinate[] inputPts = coords.toCoordinateArray();
+			Coordinate[] newPts = LittleThumblingDensifier.densifyPoints(inputPts, stepLength, parent.getPrecisionModel());
+			// prevent creation of invalid linestrings
+			if (parent instanceof LineString && newPts.length == 1) {
+				newPts = new Coordinate[0];
+			}
+			return factory.getCoordinateSequenceFactory().create(newPts);
+		}
+
+		protected Geometry transformPolygon(Polygon geom, Geometry parent) {
+			Geometry roughGeom = super.transformPolygon(geom, parent);
+			// don't try and correct if the parent is going to do this
+			if (parent instanceof MultiPolygon) {
+				return roughGeom;
+			}
+			return createValidArea(roughGeom);
+		}
+
+		protected Geometry transformMultiPolygon(MultiPolygon geom, Geometry parent) {
+			Geometry roughGeom = super.transformMultiPolygon(geom, parent);
+			return createValidArea(roughGeom);
+		}
+
+		/**
+		 * Creates a valid area geometry from one that possibly has bad topology
+		 * (i.e. self-intersections). Since buffer can handle invalid topology, but
+		 * always returns valid geometry, constructing a 0-width buffer "corrects"
+		 * the topology. Note this only works for area geometries, since buffer
+		 * always returns areas. This also may return empty geometries, if the input
+		 * has no actual area.
+		 * 
+		 * @param roughAreaGeom
+		 *          an area geometry possibly containing self-intersections
+		 * @return a valid area geometry
+		 */
+		private Geometry createValidArea(Geometry roughAreaGeom) {
+			return roughAreaGeom.buffer(0.0);
+		}
+	}
+
+}

modules/core/src/main/java/org/locationtech/jts/simplify/GaussianLineSmoothing.java

@@ -0,0 +1,154 @@
+/*
+ * Copyright (c) 2016 Vivid Solutions.
+ *
+ * All rights reserved. This program and the accompanying materials
+ * are made available under the terms of the Eclipse Public License v1.0
+ * and Eclipse Distribution License v. 1.0 which accompanies this distribution.
+ * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v10.html
+ * and the Eclipse Distribution License is available at
+ *
+ * http://www.eclipse.org/org/documents/edl-v10.php.
+ */
+package org.locationtech.jts.simplify;
+
+import org.locationtech.jts.densify.LittleThumblingDensifier;
+import org.locationtech.jts.geom.Coordinate;
+import org.locationtech.jts.geom.LineString;
+
+/**
+ * Line gaussian smoothing.
+ * 
+ * @author julien Gaffuri
+ *
+ */
+public class GaussianLineSmoothing {
+
+	/**
+	 * @param line
+	 * @param sigmaM
+	 * @return
+	 */
+	public static LineString get(LineString line, double sigmaM){ return get(line, sigmaM, -1); }
+
+	/**
+	 * Line gaussian smoothing.
+	 * The position of each point is the average position of its neighbors, weighted by a gaussian kernel.
+	 * For non-closed lines, the initial and final points are preserved.
+	 * 
+	 * @param line The input line
+	 * @param sigmaM The standard deviation of the gaussian kernel. The larger, the more smoothed.
+	 * @param resolution The target resolution of the geometry. This parameter is used to filter the final geometry.
+	 * @return
+	 */
+	public static LineString get(LineString line, double sigmaM, double resolution){
+		if(line.getCoordinates().length <= 2) return (LineString) line.copy();
+
+		boolean isClosed = line.isClosed();
+		double length = line.getLength();
+		double densifiedResolution = sigmaM/3;
+
+		//handle extreme cases: too large sigma resulting in too large densified resolution.
+		if(densifiedResolution > 0.25*length ) {
+			if(isClosed){
+				//return tiny triangle nearby center point
+				Coordinate c = line.getCentroid().getCoordinate();
+				length *= 0.01;
+				return line.getFactory().createLineString(new Coordinate[]{ new Coordinate(c.x-length,c.y-length), new Coordinate(c.x,c.y+length), new Coordinate(c.x+length,c.y-length), new Coordinate(c.x-length,c.y-length) });
+			} else {
+				//return segment
+				return line.getFactory().createLineString(new Coordinate[]{ line.getCoordinateN(0), line.getCoordinateN(line.getNumPoints()-1) });
+			}
+		}
+
+		//compute densified line
+		Coordinate[] densifiedCoords = LittleThumblingDensifier.densify(line, densifiedResolution).getCoordinates();
+
+		//build ouput line structure
+		int nb = (int) (length/densifiedResolution);
+		Coordinate[] out = new Coordinate[nb+1];
+
+		//prepare gaussian coefficients
+		int n = 7*3; //it should be: E(7*sigma/densifiedResolution), which is 7*3;
+		double gcs[] = new double[n+1];
+		{
+			double a = sigmaM*Math.sqrt(2*Math.PI);
+			double b = sigmaM*sigmaM*2;
+			double d = densifiedResolution*densifiedResolution;
+			for(int i=0; i<n+1; i++) gcs[i] = Math.exp(-i*i*d/b) /a;
+		}
+
+		Coordinate c0 = densifiedCoords[0];
+		Coordinate cN = densifiedCoords[nb];
+		for(int i=0; i<nb; i++) {
+			if(!isClosed && i==0) continue;
+
+			//compute coordinates of point i of the smoothed line (gauss mean)
+			double x=0.0, y=0.0;
+			for(int j=-n; j<=n; j++) {
+				//index of the point to consider on the original densified line
+				int q = i+j;
+				//find coordinates (xq,yq) of point q
+				double xq, yq;
+				if(q<0) {
+					if(isClosed) {
+						//make loop to get the right point
+						q = q%nb; if(q<0) q+=nb;
+						Coordinate c = densifiedCoords[q];
+						xq = c.x;
+						yq = c.y;
+					} else {
+						//get symetric point
+						q = (-q)%nb; if(q==0) q=nb;
+						Coordinate c = densifiedCoords[q];
+						xq = 2*c0.x-c.x;
+						yq = 2*c0.y-c.y;
+					}
+				} else if (q>nb) {
+					if(isClosed) {
+						//make loop to get the right point
+						q = q%nb; if(q==0) q=nb;
+						Coordinate c = densifiedCoords[q];
+						xq = c.x;
+						yq = c.y;
+					} else {
+						//get symetric point
+						q = nb-q%nb; if(q==nb) q=0;
+						Coordinate c = densifiedCoords[q];
+						xq = 2*cN.x-c.x;
+						yq = 2*cN.y-c.y;
+					}
+				} else {
+					//general case (most frequent)
+					Coordinate c = densifiedCoords[q];
+					xq = c.x;
+					yq = c.y;
+				}
+				//get gaussian coefficient
+				double gc = gcs[j>=0?j:-j];
+				//add contribution of point q to new position of point i
+				x += xq*gc;
+				y += yq*gc;
+			}
+			//assign smoothed position of point i
+			out[i] = new Coordinate(x*densifiedResolution, y*densifiedResolution);
+		}
+
+		//handle start and end points
+		if(isClosed) {
+			//ensure start and end locations are the same
+			out[nb] = out[0];
+		} else {
+			//ensure start and end points are at the same position as the initial geometry
+			out[0] = densifiedCoords[0];
+			out[nb] = densifiedCoords[densifiedCoords.length-1];
+		}
+
+		//prepare final line, applying some filtering
+		LineString lsOut = line.getFactory().createLineString(out);
+		if(resolution<0) resolution = densifiedResolution /3;
+		lsOut = (LineString) DouglasPeuckerSimplifier.simplify( lsOut , resolution);
+
+		return lsOut;
+	}
+
+}

modules/core/src/test/java/org/locationtech/jts/densify/LittleThumblingDensifierTest.java

@@ -0,0 +1,61 @@
+/*
+ * Copyright (c) 2016 Vivid Solutions.
+ *
+ * All rights reserved. This program and the accompanying materials
+ * are made available under the terms of the Eclipse Public License v1.0
+ * and Eclipse Distribution License v. 1.0 which accompanies this distribution.
+ * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v10.html
+ * and the Eclipse Distribution License is available at
+ *
+ * http://www.eclipse.org/org/documents/edl-v10.php.
+ */
+package org.locationtech.jts.densify;
+
+import java.util.Collection;
+
+import org.locationtech.jts.geom.Geometry;
+import org.locationtech.jts.io.WKTFileReader;
+import org.locationtech.jts.io.WKTReader;
+
+import junit.framework.TestCase;
+
+/**
+ * @author Julien Gaffuri
+ *
+ */
+public class LittleThumblingDensifierTest extends TestCase {
+	private final WKTReader wr = new WKTReader();
+
+	public LittleThumblingDensifierTest(String name) { super(name); }
+
+
+	public static void main(String[] args) {
+		junit.textui.TestRunner.run(LittleThumblingDensifierTest.class);
+	}
+
+	public void test0() throws Exception{
+		Geometry g = wr.read("LINESTRING(0 0, 1 0)");
+		for(double step : new double[] {0.1, 0.05, 0.5, 0.7, 0.005}) {
+			Geometry g_ = LittleThumblingDensifier.densify(g, step);
+			assertEquals(g.getLength(), g_.getLength());
+			assertEquals((int)(1.0/step)+1, g_.getNumPoints());
+		}
+	}
+
+	public void test1() throws Exception{
+		Geometry g = LittleThumblingDensifier.densify(wr.read("LINESTRING(0 0, 1 1)"), 0.1);
+		assertEquals(Math.sqrt(2), g.getLength());
+		assertEquals(15, g.getNumPoints());
+	}
+
+	public void test2() throws Exception{
+		WKTFileReader wfr = new WKTFileReader("src/test/resources/testdata/plane.wkt", new WKTReader());
+		Collection<?> gs = wfr.read();
+		Geometry g = (Geometry) gs.iterator().next();
+		Geometry g_ = LittleThumblingDensifier.densify(g, 0.1);
+
+		assertTrue(g.getLength()>g_.getLength());
+		assertEquals(2612, g_.getNumPoints());
+	}
+
+}