From a8213fe2fca2322f671f421e721cf879c2726226 Mon Sep 17 00:00:00 2001 From: Jes Wulfsberg Nielsen Date: Thu, 10 Feb 2022 02:12:16 +0100 Subject: [PATCH] CatmullRomSpline interpolation --- .../jts/shape/CatmullRomSpline.java | 384 ++++++++++++++++++ .../jts/shape/CatmullRomSplineTest.java | 140 +++++++ 2 files changed, 524 insertions(+) create mode 100644 modules/core/src/main/java/org/locationtech/jts/shape/CatmullRomSpline.java create mode 100644 modules/core/src/test/java/org/locationtech/jts/shape/CatmullRomSplineTest.java diff --git a/modules/core/src/main/java/org/locationtech/jts/shape/CatmullRomSpline.java b/modules/core/src/main/java/org/locationtech/jts/shape/CatmullRomSpline.java new file mode 100644 index 0000000000..d6bb27edc9 --- /dev/null +++ b/modules/core/src/main/java/org/locationtech/jts/shape/CatmullRomSpline.java @@ -0,0 +1,384 @@ +/* + * Copyright (c) 2022 Jes Wuilfsberg Nielsen. + * + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License 2.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-v20.html + * and the Eclipse Distribution License is available at + * + * http://www.eclipse.org/org/documents/edl-v10.php. + */ +package org.locationtech.jts.shape; + +import org.locationtech.jts.geom.*; +import org.locationtech.jts.geom.util.GeometryMapper; + +/** + * Creates a curved geometry by replacing segments of the input with Catmull-Rom splines. + * The Coordinates of the input geometry are used as control points for the spline, rolling a window of four coordinates + * along the geometry to calculate the spline for each line segment (either looping around at the ends of a ring, or + * extrapolating endpoints linearly from linestrings to get the first and last control point). + *

+ * The implementation uses the "2D coordinates with additional values" approach which JTS generally uses, in the sense + * that it treats the geometries as 2D, but allows for them to carry Z and/or M values. The Z and M values will be + * interpolated if they exist, but they are not used to define the shape of the spline. + *

+ * The implementation is based on the example given at https://stackoverflow.com/a/23980479 + * / https://ideone.com/NoEbVM + * with a nod to https://github.com/mrdoob/three.js for further inspiration. + */ +public class CatmullRomSpline { + private final boolean includeZ; + private final boolean includeM; + + /** + * As {@link #catmullRomSpline(Geometry, int, double, double)} with alpha = 0.5 (centripetal) and + * endpointExtrapolationWeight = 1.0. + */ + public static Geometry catmullRomSpline(Geometry g, int numberOfSegments) { + return catmullRomSpline(g, numberOfSegments, 0.5, 1); + } + + /** + * Creates a geometry of linearized Catmull-Rom splines defined by the coordinates of each line segment in the input + * geometry, with the immediately preceding and following coordinates used as control points. In other words, it + * uses the input geometry as the control curve for a Catmull-Rom spline. + *

+ * In the case of LinearRings (including interior and exterior rings of polygons), the "preceding" and "following" + * coordinates at the end of the ring are calculated by looping around to the end/start of the ring. + *

+ * In the case of LineStrings, additional points are linearly extrapolated at the start and end of the geometry to be + * used as control points. + * + * @param geom The geometry defining the control curve. (Multigeometries and polygons are + * valid input. Each constituent geometry will be calculated separately). + * @param numberOfSegments The number of spline segments on each segment of the input geometry + * @param alpha Describes the "looseness" of the result, from 0.0 to 1.0. + *

+ * @param endpointExtrapolationWeight When extra control points are linearly extrapolated at the start and end of + * linestrings, this parameter adjusts the length of that extrapolation, affecting + * the "lead-in" to the actual curve, and hence the shape of the first and last + * segment (though the impact is fairly minimal). + * @return the linearized curved geometry + */ + public static Geometry catmullRomSpline(Geometry geom, int numberOfSegments, double alpha, double endpointExtrapolationWeight) { + // Pick the best Coordinate type amongst known JTS types: + Coordinate c = geom.getCoordinate(); + boolean includeZ = !(c instanceof CoordinateXY || c instanceof CoordinateXYM); + boolean includeM = c instanceof CoordinateXYM || c instanceof CoordinateXYZM; + return new CatmullRomSpline(includeZ, includeM).catmullRomGeometry(geom, numberOfSegments, alpha, endpointExtrapolationWeight); + } + + /** + * As {@link #CatmullRomSpline(boolean, boolean)}, with includeZ = true, includeM = false. + */ + public CatmullRomSpline() { + this(true, false); + } + + /** + * Instantiate a version of the {@link CatmullRomSpline} to use specific functionality from it if the static + * {@link #catmullRomSpline(Geometry, int, double, double)} does not suffice. + * + * @param includeZ Whether the generated coordinates include the Z ordinate. ({@link Coordinate} or {@link CoordinateXYZM}) + * @param includeM Whether the generated coordinates include the M ordinate. ({@link CoordinateXYM} or {@link CoordinateXYZM}) + */ + public CatmullRomSpline(boolean includeZ, boolean includeM) { + this.includeZ = includeZ; + this.includeM = includeM; + } + + /** + * Creates a geometry of linearized Catmull-Rom splines defined by the coordinates of each line segment in the input + * geometry. + *

+ * Parameters as per {@link #catmullRomSpline(Geometry, int, double, double)} + */ + public Geometry catmullRomGeometry(Geometry g, int numberOfSegments, double alpha, double endpointExtrapolationWeight) { + return GeometryMapper.flatMap(g, 1, geom -> { + if (geom instanceof LinearRing) { + return catmullRomLinearRing((LinearRing) geom, numberOfSegments, alpha); + } + if (geom instanceof LineString) { + return catmullRomLineString((LineString) geom, numberOfSegments, alpha, endpointExtrapolationWeight); + } + if (geom instanceof Polygon) { + return catmullRomPolygon((Polygon) geom, numberOfSegments, alpha); + } + return geom.copy(); + }); + } + + /** + * Creates a LineString of linearized Catmull-Rom splines defined by the coordinates of each line segment in the input + * geometry. This is a specialized form with stronger type signature, but parameters are otherwise as + * {@link #catmullRomSpline(Geometry, int, double, double)} + */ + public LineString catmullRomLineString(LineString l, int numberOfSegments, double alpha, double endpointExtrapolationWeight) { + return l.getFactory().createLineString(this.calculateSplineCoordinates(l.getCoordinates(), numberOfSegments, alpha, false, endpointExtrapolationWeight)); + } + + /** + * Creates a LinearRing of linearized Catmull-Rom splines defined by the coordinates of each line segment in the input + * geometry. The created geometry will loop around, with the last point of the input connecting smoothly to the first. + *

+ * The input does not need to be a valid LinearRing. If the last point of the input equals the first (per + * {@link Coordinate#equals2D(Coordinate, double)}), the last point will be skipped to avoid duplicate coordinates. + *

+ * Parameters as per {@link #catmullRomSpline(Geometry, int, double, double)}, except endpointExtrapolationWeight is + * not needed for loops. + */ + public LinearRing catmullRomLinearRing(LineString l, int numberOfSegments, double alpha) { + return l.getFactory().createLinearRing(this.calculateSplineCoordinates(l.getCoordinates(), numberOfSegments, alpha, true, Double.NaN)); + } + + /** + * Creates a Polygon of linearized Catmull-Rom splines defined by the coordinates of each line segment in the input + * geometry. + *

+ * This method will calculate a {@link #catmullRomLinearRing(LineString, int, double)} for the exterior and + * interior rings of the polygon. Note that there is no guarantee that the generated polygon will be valid, as e.g. + * the spline interpolation of holes could make them overlap or cross outside the exterior hull. + */ + public Polygon catmullRomPolygon(Polygon poly, int numberOfSegments, double alpha) { + LinearRing shell = catmullRomLinearRing(poly.getExteriorRing(), numberOfSegments, alpha); + LinearRing[] holes = null; + if (poly.getNumInteriorRing() > 0) { + holes = new LinearRing[poly.getNumInteriorRing()]; + for (int i = 0; i < poly.getNumInteriorRing(); i++) { + holes[i] = catmullRomLinearRing(poly.getInteriorRingN(i), numberOfSegments, alpha); + } + } + return poly.getFactory().createPolygon(shell, holes); + } + + /** + * The "raw" Catmull-Rom spline interpolation of coordinates. The various geometry-generating methods use this + * internally. + * + * @param coords An array of coordinates defining the control points of the spline. + * @param numberOfSegments The number of spline segments between each coordinate. + * @param alpha Describes the "looseness" of the result, from 0.0 to 1.0. + *

+ * @param loop Whether the spline forms a closed loop (ring) or a linestring. + * @param endpointExtrapolationWeight If loop is false, a "phantom" control point will be linearly extrapolated at + * each end of the input line, to function as control points for the lead-in and + * -out to the spline. The endpointExtrapolationWeight determines the length of + * this extrapolation, as a factor of the distance between the two points used for + * the extrapolation. 1 is a sensible default, keeping the same distance as the + * existing points. A larger number will make the spline run at a flatter angle + * to the control line for a longer while before it turns away towards the next + * control point. + * @return An array of interpolated coordinates. + */ + public Coordinate[] calculateSplineCoordinates(Coordinate[] coords, int numberOfSegments, double alpha, boolean loop, double endpointExtrapolationWeight) { + CoordinateList coordList = this.cleanInput(coords, loop, endpointExtrapolationWeight); + if (coordList == null) { + return coords; + } + + // Initialize the polynomial structure for each ordinate: + CubicPolynomial polyX, polyY, polyZ = null, polyM = null; + polyX = new CubicPolynomial(); + polyY = new CubicPolynomial(); + if (this.includeZ) { + polyZ = new CubicPolynomial(); + } + if (this.includeM) { + polyM = new CubicPolynomial(); + } + + // Prep the result array with the known length: + Coordinate[] resultArray = new Coordinate[(coordList.size() - 3) * numberOfSegments + 1]; + int resultIndex = 0; + + // Spin over each segment in the input, possibly looping back to the start point: + int len = coordList.size() - 2; + for (int i = 1; i < len; i++) { + Coordinate c0 = coordList.get(i - 1); + Coordinate c1 = coordList.get(i); + Coordinate c2 = coordList.get(i + 1); + Coordinate c3 = coordList.get(i + 2); + // A somewhat non-Java-idiomatic way, reusing allocated structures rather than creating new instances: + this.initializePolynomialsForSegment(c0, c1, c2, c3, alpha, polyX, polyY, polyZ, polyM); + + // Spin over the wanted number of steps on each input segment: + for (int j = 0; j < numberOfSegments; j++) { + double t = j / (double) numberOfSegments; + resultArray[resultIndex++] = newCoord( + polyX.eval(t), + polyY.eval(t), + polyZ == null ? Double.NaN : polyZ.eval(t), + polyM == null ? Double.NaN : polyM.eval(t) + ); + } + } + // And drop the final control point in: + Coordinate lastCoord = coordList.get(coordList.size() - 2); + resultArray[resultIndex] = newCoord(lastCoord.x, lastCoord.y, lastCoord.getZ(), lastCoord.getM()); + return resultArray; + } + + /** + * Cleans the input of duplicate coordinates, as well as inserts the "phantom" control points at the start and end of + * the coordinate list. + * + * @return A CoordinateList with the full list of control points to iterate over, or null if cleaning it resulted in + * too few coordinates to be able to interpolate. + */ + private CoordinateList cleanInput(Coordinate[] coords, boolean loop, double endpointExtrapolationWeight) { + // Can't interpolate less than a line: + if (coords == null || coords.length < 2) { + return null; + } + + CoordinateList coordList = new CoordinateList(); + + // Placeholder for the first coord. We defer actually picking it until we have weeded out duplicates, since we may + // need to extrapolate from input points. + coordList.add(new Coordinate(Double.NaN, Double.NaN, Double.NaN)); + + // The input coords, non-repeating, and possibly skipping the last one if it's equal to the first. + // If we have a ring (touching endpoints), ignore the last, duplicated point. + int numberOfPoints = loop && coords[0].equals2D(coords[coords.length - 1]) + ? coords.length - 1 + : coords.length; + for (int i = 0; i < numberOfPoints; i++) { + coordList.add(coords[i], false); + } + + // If we do not have at least two unique coordinates (plus the placeholder), we cannot interpolate: + if (coordList.size() < 3) { + return null; + } + + if (loop) { + // The first control point in a loop is the last point in the (cleaned) input + coordList.set(0, coordList.get(coordList.size() - 1)); + // And for the end, we loop back to the first (input) point, and pick the next as control point: + coordList.add(coordList.getCoordinate(1)); + coordList.add(coordList.get(2)); + } else { + // Extrapolate backwards along the line to find first control point + coordList.set(0, this.extrapolate(coordList.get(2), coordList.get(1), endpointExtrapolationWeight)); + // ...and forward beyond the end: + coordList.add(this.extrapolate(coordList.get(coordList.size() - 2), coordList.get(coordList.size() - 1), endpointExtrapolationWeight)); + } + return coordList; + } + + /** + * Conceptually, the algorithm is based on distance, but since it's only ever used in an exponential function, we just + * keep it as squared and divide the exponent by two later. + * If we want true 3D or 4D interpolation, we simply need to include those coordinates in the distance function here. + */ + double distSquared(Coordinate p, Coordinate q) { + double dx = q.x - p.x; + double dy = q.y - p.y; + return dx * dx + dy * dy; + } + + void initializePolynomialsForSegment(Coordinate p0, + Coordinate p1, + Coordinate p2, + Coordinate p3, + double alpha, + CubicPolynomial polyX, + CubicPolynomial polyY, + CubicPolynomial polyZ, + CubicPolynomial polyM) { + double dt0 = Math.pow(distSquared(p0, p1), alpha * 0.5); + double dt1 = Math.pow(distSquared(p1, p2), alpha * 0.5); + double dt2 = Math.pow(distSquared(p2, p3), alpha * 0.5); + + polyX.initNonuniformCatmullRom(p0.x, p1.x, p2.x, p3.x, dt0, dt1, dt2); + polyY.initNonuniformCatmullRom(p0.y, p1.y, p2.y, p3.y, dt0, dt1, dt2); + if (polyZ != null) { + polyZ.initNonuniformCatmullRom(p0.getZ(), p1.getZ(), p2.getZ(), p3.getZ(), dt0, dt1, dt2); + } + if (polyM != null) { + polyM.initNonuniformCatmullRom(p0.getM(), p1.getM(), p2.getM(), p3.getM(), dt0, dt1, dt2); + } + } + + Coordinate extrapolate(Coordinate c1, Coordinate c2, double weight) { + return newCoord(c2.x + (c2.x - c1.x) * weight, c2.y + (c2.y - c1.y) * weight, c2.getZ() + (c2.getZ() - c1.getZ()) * weight, c2.getM() + (c2.getM() - c1.getM()) * weight); + } + + private Coordinate newCoord(double x, double y, double z, double m) { + if (includeZ) { + if (includeM) { + return new CoordinateXYZM(x, y, z, m); + } else { + return new Coordinate(x, y, z); + } + } else { + if (includeM) { + return new CoordinateXYM(x, y, m); + } else { + return new CoordinateXY(x, y); + } + } + } + + /** + * Reuse the data structure and formula for each ordinate. This includes pre-calculations for each segment, making the + * interpolation over t simpler. + */ + private static class CubicPolynomial { + double c0, c1, c2, c3; + boolean finite; + + public void init(double c0, double c1, double tangent0, double tangent1) { + this.c0 = c0; + this.c1 = tangent0; + this.c2 = -3 * c0 + 3 * c1 - 2 * tangent0 - tangent1; + this.c3 = 2 * c0 - 2 * c1 + tangent0 + tangent1; + } + + public void initNonuniformCatmullRom(double c0, double c1, double c2, double c3, double dt0, double dt1, double dt2) { + // I don't know if this optimization is worth it, but it seems very likely that the most common use will be + // "standard" Coordinate objects without actual Z values, so it will probably be a win to early-out those. + this.finite = Double.isFinite(c0) && Double.isFinite(c1) && Double.isFinite(c2) && Double.isFinite(c3) && Double.isFinite(dt0) && Double.isFinite(dt1) && Double.isFinite(dt2); + if (this.finite) { + // Compute tangents + double t1 = (c1 - c0) / dt0 - (c2 - c0) / (dt0 + dt1) + (c2 - c1) / dt1; + double t2 = (c2 - c1) / dt1 - (c3 - c1) / (dt1 + dt2) + (c3 - c2) / dt2; + + // Normalize into 0-1 range + t1 *= dt1; + t2 *= dt1; + + this.init(c1, c2, t1, t2); + } + } + + double eval(double t) { + if (!this.finite) { + return Double.NaN; + } + double t2 = t * t; + double t3 = t2 * t; + return c0 + c1 * t + c2 * t2 + c3 * t3; + } + } +} diff --git a/modules/core/src/test/java/org/locationtech/jts/shape/CatmullRomSplineTest.java b/modules/core/src/test/java/org/locationtech/jts/shape/CatmullRomSplineTest.java new file mode 100644 index 0000000000..dbad234454 --- /dev/null +++ b/modules/core/src/test/java/org/locationtech/jts/shape/CatmullRomSplineTest.java @@ -0,0 +1,140 @@ +package org.locationtech.jts.shape; + +import junit.textui.TestRunner; +import org.locationtech.jts.geom.*; +import test.jts.GeometryTestCase; + +public class CatmullRomSplineTest extends GeometryTestCase { + + public static void main(String[] args) { + TestRunner.run(CatmullRomSplineTest.class); + } + + public CatmullRomSplineTest(String name) { + super(name); + } + + public void testSimpleLineUniform() { + checkCurve("LINESTRING(-1 0, -2 1, 2 1, 1 0)", 10, 0.0, 1, + "LINESTRING (-1 0, -1.12 0.1, -1.28 0.22, -1.46 0.33, -1.64 0.45, -1.81 0.56, -1.96 0.67, -2.07 0.77, -2.12 0.86, -2.1 0.94, -2 1, -1.78 1.05, -1.44 1.08, -1.01 1.11, -0.52 1.12, 0 1.13, 0.52 1.12, 1.01 1.11, 1.44 1.08, 1.78 1.05, 2 1, 2.1 0.94, 2.12 0.86, 2.07 0.77, 1.96 0.67, 1.81 0.56, 1.64 0.45, 1.46 0.33, 1.28 0.22, 1.12 0.1, 1 0)"); + } + + public void testSimpleCentripetal() { + checkCurve("LINESTRING(-1 0, -2 1, 2 1, 1 0)", 10, 0.5, 1, + "LINESTRING (-1 0, -1.11 0.1, -1.24 0.21, -1.38 0.32, -1.52 0.44, -1.66 0.55, -1.78 0.65, -1.89 0.75, -1.96 0.85, -2 0.93, -2 1, -1.86 1.09, -1.54 1.17, -1.1 1.22, -0.57 1.25, 0 1.26, 0.57 1.25, 1.1 1.22, 1.54 1.17, 1.86 1.09, 2 1, 2 0.93, 1.96 0.85, 1.89 0.75, 1.78 0.65, 1.66 0.55, 1.52 0.44, 1.38 0.32, 1.24 0.21, 1.11 0.1, 1 0)"); + } + + public void testSimpleLineChordal() { + checkCurve("LINESTRING(-1 0, -2 1, 2 1, 1 0)", 10, 1.0, 1, + "LINESTRING (-1 0, -1.11 0.1, -1.22 0.21, -1.34 0.32, -1.46 0.43, -1.58 0.53, -1.69 0.64, -1.79 0.74, -1.88 0.83, -1.95 0.92, -2 1, -1.96 1.19, -1.68 1.33, -1.22 1.44, -0.64 1.5, 0 1.52, 0.64 1.5, 1.22 1.44, 1.68 1.33, 1.96 1.19, 2 1, 1.95 0.92, 1.88 0.83, 1.79 0.74, 1.69 0.64, 1.58 0.53, 1.46 0.43, 1.34 0.32, 1.22 0.21, 1.11 0.1, 1 0)"); + } + + public void testPolygonWithHoles() { + checkCurve("POLYGON ((0 0, 0 10, 6 10, 6 0, 0 0), (3 4, 1 8, 3 9, 5 8, 3 4), (3 1, 2 2, 3 3, 4 2, 3 1))", 10, 0.5, 1, + "POLYGON ((-0.39 0.59, -0.7 1.46, -0.92 2.53, -1.05 3.73, -1.09 5, -1.05 6.27, -0.92 7.47, -0.7 8.54, -0.39 9.41, 0 10, 0.41 10.3, 0.95 10.54, 1.58 10.71, 2.27 10.81, 3 10.85, 3.73 10.81, 4.42 10.71, 5.05 10.54, 5.59 10.3, 6 10, 6.39 9.41, 6.7 8.54, 6.92 7.47, 7.05 6.27, 7.09 5, 7.05 3.73, 6.92 2.53, 6.7 1.46, 6.39 0.59, 6 0, 5.59 -0.3, 5.05 -0.54, 4.42 -0.71, 3.73 -0.81, 3 -0.85, 2.27 -0.81, 1.58 -0.71, 0.95 -0.54, 0.41 -0.3, 0 0, -0.39 0.59), (3.23 4.09, 3.49 4.34, 3.78 4.71, 4.07 5.17, 4.35 5.69, 4.61 6.23, 4.82 6.77, 4.96 7.27, 5.03 7.69, 5 8, 4.91 8.17, 4.78 8.33, 4.61 8.47, 4.4 8.61, 4.18 8.72, 3.94 8.82, 3.69 8.89, 3.45 8.95, 3.21 8.99, 3 9, 2.79 8.99, 2.55 8.95, 2.31 8.89, 2.06 8.82, 1.82 8.72, 1.6 8.61, 1.39 8.47, 1.22 8.33, 1.09 8.17, 1 8, 0.97 7.69, 1.04 7.27, 1.18 6.77, 1.39 6.23, 1.65 5.69, 1.93 5.17, 2.22 4.71, 2.51 4.34, 2.77 4.09, 3 4, 3.23 4.09), (3.11 1.02, 3.23 1.07, 3.36 1.15, 3.5 1.26, 3.63 1.38, 3.74 1.5, 3.85 1.64, 3.93 1.77, 3.98 1.89, 4 2, 3.98 2.11, 3.93 2.23, 3.85 2.36, 3.74 2.5, 3.63 2.63, 3.5 2.74, 3.36 2.85, 3.23 2.93, 3.11 2.98, 3 3, 2.89 2.98, 2.77 2.93, 2.64 2.85, 2.5 2.74, 2.38 2.63, 2.26 2.5, 2.15 2.36, 2.07 2.23, 2.02 2.11, 2 2, 2.02 1.89, 2.07 1.77, 2.15 1.64, 2.26 1.5, 2.38 1.38, 2.5 1.26, 2.64 1.15, 2.77 1.07, 2.89 1.02, 3 1, 3.11 1.02))"); + } + + public void testSimpleCentripetalWithLongExtrapolation() { + checkCurve("LINESTRING(-1 0, -2 1, 2 1, 1 0)", 10, 0.5, 10, + "LINESTRING (-1 0, -1.15 0.15, -1.31 0.28, -1.46 0.4, -1.6 0.51, -1.72 0.61, -1.83 0.7, -1.92 0.79, -1.98 0.86, -2.01 0.93, -2 1, -1.86 1.09, -1.54 1.17, -1.1 1.22, -0.57 1.25, 0 1.26, 0.57 1.25, 1.1 1.22, 1.54 1.17, 1.86 1.09, 2 1, 2.01 0.93, 1.98 0.86, 1.92 0.79, 1.83 0.7, 1.72 0.61, 1.6 0.51, 1.46 0.4, 1.31 0.28, 1.15 0.15, 1 0)"); + } + + public void testIgnoresRepeatedCoordsAtStart() { + checkCurve("LINESTRING (-2 0, -2 0, 0 1, 2 0)", 10, 0.5, 1, + "LINESTRING (-2 0, -1.8 0.11, -1.6 0.23, -1.4 0.36, -1.2 0.5, -1 0.63, -0.8 0.74, -0.6 0.85, -0.4 0.93, -0.2 0.98, 0 1, 0.2 0.98, 0.4 0.93, 0.6 0.85, 0.8 0.74, 1 0.63, 1.2 0.5, 1.4 0.36, 1.6 0.23, 1.8 0.11, 2 0)"); + } + + public void testIgnoresRepeatedCoordsInMiddle() { + checkCurve("LINESTRING (-2 0, 0 1, 0 1, 2 0)", 10, 0.5, 1, + "LINESTRING (-2 0, -1.8 0.11, -1.6 0.23, -1.4 0.36, -1.2 0.5, -1 0.63, -0.8 0.74, -0.6 0.85, -0.4 0.93, -0.2 0.98, 0 1, 0.2 0.98, 0.4 0.93, 0.6 0.85, 0.8 0.74, 1 0.63, 1.2 0.5, 1.4 0.36, 1.6 0.23, 1.8 0.11, 2 0)"); + } + + public void testIgnoresRepeatedCoordsAtEnd() { + checkCurve("LINESTRING (-2 0, 0 1, 2 0, 2 0)", 10, 0.5, 1, + "LINESTRING (-2 0, -1.8 0.11, -1.6 0.23, -1.4 0.36, -1.2 0.5, -1 0.63, -0.8 0.74, -0.6 0.85, -0.4 0.93, -0.2 0.98, 0 1, 0.2 0.98, 0.4 0.93, 0.6 0.85, 0.8 0.74, 1 0.63, 1.2 0.5, 1.4 0.36, 1.6 0.23, 1.8 0.11, 2 0)"); + } + + public void testZ() { + Geometry geom = read("LINESTRINGZ(0 0 0, 1 1 2, 2 0 3)"); + Geometry actual = CatmullRomSpline.catmullRomSpline(geom, 10); + Coordinate[] actualCoords = actual.getCoordinates(); + // Sanity-check that we have filled out the Z-value + assertEquals(1.0625, actualCoords[5].getZ(), 0.0001); + // Check that we handle the "last point" special case correct as well + assertEquals(3, actualCoords[actualCoords.length - 1].getZ(), 0.0001); + } + + public void testM() { + CatmullRomSpline catmullRom = new CatmullRomSpline(false, true); + Coordinate[] coords = new Coordinate[3]; + coords[0] = new CoordinateXYM(0, 0, 0); + coords[1] = new CoordinateXYM(1, 1, 2); + coords[2] = new CoordinateXYM(2, 0, 3); + Geometry geom = new GeometryFactory().createLineString(coords); + Geometry actual = CatmullRomSpline.catmullRomSpline(geom, 10); + Coordinate[] actualCoords = actual.getCoordinates(); + // Sanity-check that we have filled out the M-value + assertEquals(1.0625, actualCoords[5].getM(), 0.0001); + // Check that we handle the "last point" special case correct as well + assertEquals(3, actualCoords[actualCoords.length - 1].getM(), 0.0001); + } + + public void testZMRing() { + Coordinate[] coords = new Coordinate[4]; + coords[0] = new CoordinateXYZM(0, 0, 0, 0); + coords[1] = new CoordinateXYZM(1, 1, 10, 10); + coords[2] = new CoordinateXYZM(2, 0, 20, 0); + coords[3] = new CoordinateXYZM(0, 0, 0, 0); + + Geometry geom = new GeometryFactory().createLinearRing(coords); + Geometry actual = CatmullRomSpline.catmullRomSpline(geom, 10); + Coordinate[] actualCoords = actual.getCoordinates(); + + // Sanity check that the Z and M are filled out (with values which match 10 times the corresponding X and Y) + assertEquals(3.4687, actualCoords[5].getZ(), 0.001); + assertEquals(5.6790, actualCoords[5].getM(), 0.001); + } + + public void testXY() { + Geometry geom = read("LINESTRING(0 0, 1 1, 2 0)"); + // Large lead-in weight to exaggerate the effect into something which is visible within the test tolerance: + Geometry actual = new CatmullRomSpline(false, false) + .catmullRomGeometry(geom, 10, 0.5, 10000); + + // This is the curve as it would look with weight 1, and it does indeed differ from the actual. + // Geometry expected = read( "LINESTRING (0 0, 0.1 0.109, 0.2 0.232, 0.3 0.363, 0.4 0.4960000000000001, 0.5 0.625, 0.6 0.744, 0.7 0.847, 0.8 0.928, 0.9 0.9809999999999999, 1 1, 1.1 0.981, 1.2 0.9279999999999999, 1.3 0.8470000000000001, 1.4 0.744, 1.5 0.625, 1.6 0.496, 1.7 0.363, 1.8 0.2319999999999999, 1.9 0.109, 2 0)"); + // This is the expected with weight 10000: + Geometry expected = read("LINESTRING (0 0, 0.1793960396039605 0.1883960396039605, 0.3254653465346536 0.3574653465346537, 0.4440891089108913 0.5070891089108913, 0.5411485148514853 0.6371485148514854, 0.6225247524752477 0.7475247524752477, 0.6940990099009903 0.8380990099009903, 0.7617524752475247 0.9087524752475247, 0.8313663366336633 0.9593663366336634, 0.9088217821782177 0.9898217821782177, 1 1, 1.0911782178217821 0.9898217821782178, 1.1686336633663368 0.9593663366336632, 1.2382475247524756 0.9087524752475244, 1.3059009900990102 0.8380990099009896, 1.377475247524753 0.747524752475247, 1.4588514851485157 0.6371485148514845, 1.5559108910891095 0.5070891089108904, 1.6745346534653471 0.3574653465346527, 1.8206039603960398 0.18839603960396, 2 0)"); + checkEqual(expected, actual, 0.01); + assertEquals(CoordinateXY.class, actual.getCoordinate().getClass()); + } + + public void testCopiesPoints() { + Geometry geom = read("POINT(10 20)"); + Geometry actual = new CatmullRomSpline().catmullRomGeometry(geom, 10, 0.5, 1); + checkEqual(geom, actual); + assertNotSame(geom, actual); + } + + public void testEarlyOutWhenTooFewCoords() { + Coordinate[] coords = new Coordinate[1]; + coords[0] = new Coordinate(1, 2, 3); + Coordinate[] actual = new CatmullRomSpline().calculateSplineCoordinates(coords, 10, 0.5, false, 1); + assertSame(coords, actual); + } + + public void testEarlyOutWhenCollapsingIntoTooFewCoords() { + Coordinate[] coords = new Coordinate[2]; + coords[0] = new Coordinate(1, 2, 3); + coords[1] = new Coordinate(1, 2, 3); + Coordinate[] actual = new CatmullRomSpline().calculateSplineCoordinates(coords, 10, 0.5, false, 1); + assertSame(coords, actual); + } + + private void checkCurve(String wkt, int numberOfSegments, double alpha, double extrapolationWeight, String wktExpected) { + Geometry geom = read(wkt); + Geometry actual = CatmullRomSpline.catmullRomSpline(geom, numberOfSegments, alpha, extrapolationWeight); + Geometry expected = read(wktExpected); + checkEqual(expected, actual, 0.01); + } +}