AreaPrecisionPerfTest.java

/*
 * Copyright (c) 2016 Martin Davis.
 *
 * 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 test.jts.perf.algorithm;

import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.LinearRing;
import org.locationtech.jts.geom.Polygon;

public class AreaPrecisionPerfTest
{
  public static void main(String[] args) throws Exception
  {

    double originX = 1000000;
    double originY = 5000000;
    long start = System.currentTimeMillis();

    for (int nrVertices = 4; nrVertices <= 1000000; nrVertices *= 2) {
      Coordinate[] coordinates = new Coordinate[nrVertices + 1];

      Coordinate vertex;
      for (int i = 0; i <= nrVertices; i++) {
        vertex = new Coordinate(originX
            + (1 + Math.sin((float) i / (float) nrVertices * 2 * Math.PI)),
            originY
                + (1 + Math.cos((float) i / (float) nrVertices * 2 * Math.PI)));
        coordinates[i] = vertex;
      }
      // close ring
      coordinates[nrVertices] = coordinates[0];
      
      Geometry g1 = new GeometryFactory().createLinearRing(coordinates);
      LinearRing[] holes = new LinearRing[] {};
      Polygon polygon = (Polygon) new GeometryFactory().createPolygon(
          (LinearRing) g1, holes);
      System.out.println(polygon);
      
      double area = originalSignedArea(coordinates);
      double area2 = accurateSignedArea(coordinates);
      double exactArea = 0.5 * nrVertices * Math.sin(2 * Math.PI / nrVertices);
      
      double eps = exactArea - area;
      double eps2 = exactArea - area2;
      
      System.out.println(nrVertices + "   orig err: " + eps 
          + "    acc err: " + eps2);
    }
    System.out.println("Time: " + (System.currentTimeMillis() - start) / 1000.0);
  }

  public static double originalSignedArea(Coordinate[] ring)
  {
    if (ring.length < 3)
      return 0.0;
    double sum = 0.0;
    for (int i = 0; i < ring.length - 1; i++) {
      double bx = ring[i].x;
      double by = ring[i].y;
      double cx = ring[i + 1].x;
      double cy = ring[i + 1].y;
      sum += (bx + cx) * (cy - by);
    }
    return -sum / 2.0;
  }

  public static double accurateSignedArea(Coordinate[] ring)
  {
    if (ring.length < 3)
      return 0.0;
    double sum = 0.0;
    // http://en.wikipedia.org/wiki/Shoelace_formula
    double x0 = ring[0].x;
    for (int i = 1; i < ring.length - 1; i++) {
      double x = ring[i].x - x0;
      double y1 = ring[i + 1].y;
      double y2 = ring[i == 0 ? ring.length - 1 : i - 1].y;
      sum += x * (y2 - y1);
    }
    return sum / 2.0;
  }

}