operation/overlayng/OverlayNG.java

/*
 * Copyright (c) 2019 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 org.locationtech.jts.operation.overlayng;

import java.util.Collection;
import java.util.List;

import org.locationtech.jts.geom.Envelope;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.LineString;
import org.locationtech.jts.geom.Location;
import org.locationtech.jts.geom.Point;
import org.locationtech.jts.geom.Polygon;
import org.locationtech.jts.geom.PrecisionModel;
import org.locationtech.jts.geom.TopologyException;
import org.locationtech.jts.geomgraph.Label;
import org.locationtech.jts.noding.MCIndexNoder;
import org.locationtech.jts.noding.Noder;
import org.locationtech.jts.noding.snap.SnappingNoder;
import org.locationtech.jts.noding.snapround.SnapRoundingNoder;
import org.locationtech.jts.operation.overlay.OverlayOp;

/**
 * Computes the geometric overlay of two {@link Geometry}s,
 * using an explicit precision model to allow robust computation.
 * The overlay can be used to determine any of the
 * following set-theoretic operations (boolean combinations) of the geometries:
 * <ul>
 * <li>{@link INTERSECTION} - all points which lie in both geometries
 * <li>{@link UNION} - all points which lie in at least one geometry
 * <li>{@link DIFFERENCE} - all points which lie in the first geometry but not the second
 * <li>{@link SYMDIFFERENCE} - all points which lie in one geometry but not both
 * </ul>
 * Input geometries may have different dimension.
 * Input collections must be homogeneous (all elements must have the same dimension).
 * <p>
 * The precision model used for the computation can be supplied
 * independent of the precision model of the input geometry.
 * The main use for this is to allow using a fixed precision
 * for geometry with a floating precision model.
 * This does two things: ensures robust computation;
 * and forces the output to be validly rounded to the precision model.
 * <p>
 * For fixed precision models noding is performed using a {@link SnapRoundingNoder}.
 * This provides robust computation (as long as precision is limited to
 * around 13 decimal digits).
 * <p>
 * For floating precision an {@link MCIndexNoder} is used.
 * This is not fully robust, so can sometimes result in
 * {@link TopologyException}s being thrown.
 * For robust full-precision overlay see {@link OverlayNGRobust}.
 * <p>
 * A custom {@link Noder} can be supplied.
 * This allows using a more performant noding strategy in specific cases,
 * for instance in {@link CoverageUnion}.
 * <p>
 * <b>Note:</b If a {@link SnappingNoder} is used
 * it is best to specify a fairly small snap tolerance,
 * since the intersection clipping optimization can
 * interact with the snapping to alter the result.
 * <p>
 * Optionally the overlay computation can process using strict mode
 * (via {@link #setStrictMode(boolean)}.
 * In strict mode result semantics are:
 * <ul>
 * <li>Lines and Points resulting from topology collapses are not included in the result</li>
 * <li>Result geometry is homogeneous
 *     for the {@link #INTERSECTION} and {@link #DIFFERENCE} operations.
 * <li>Result geometry is homogeneous
 *     for the {@link #UNION} and {@link #SYMDIFFERENCE} operations
 *     if the inputs have the same dimension</li>
 * </ul>
 * Strict mode has the following benefits:
 * <ul>
 * <li>Results are simpler</li>
 * <li>Overlay operations are chainable
 *     without needing to remove lower-dimension elements</li>
 * </ul>
 * The original JTS overlay semantics corresponds to non-strict mode.
 * <p>
 * If a robustness error occurs, a {@link TopologyException} is thrown.
 * These are usually caused by numerical rounding causing the noding output
 * to not be fully noded.
 * For robust computation with full-precision {@link OverlayNGRobust} can be used.
 *
 * @author mdavis
 *
 * @see OverlayNGRobust
 *
 */
public class OverlayNG
{
  /**
   * The code for the Intersection overlay operation.
   */
  public static final int INTERSECTION  = OverlayOp.INTERSECTION;

  /**
   * The code for the Union overlay operation.
   */
  public static final int UNION         = OverlayOp.UNION;

  /**
   *  The code for the Difference overlay operation.
   */
  public static final int DIFFERENCE    = OverlayOp.DIFFERENCE;

  /**
   *  The code for the Symmetric Difference overlay operation.
   */
  public static final int SYMDIFFERENCE = OverlayOp.SYMDIFFERENCE;

  /**
   * The default setting for Strict Mode.
   *
   * The original JTS overlay semantics used non-strict result
   * semantics, including;
   * - An Intersection result can be mixed-dimension,
   *   due to inclusion of intersection components of all dimensions
   * - Results can include lines caused by Area topology collapse
   */
  static final boolean STRICT_MODE_DEFAULT = false;

  /**
   * Tests whether a point with a given topological {@link Label}
   * relative to two geometries is contained in
   * the result of overlaying the geometries using
   * a given overlay operation.
   * <p>
   * The method handles arguments of {@link Location#NONE} correctly
   *
   * @param label the topological label of the point
   * @param opCode the code for the overlay operation to test
   * @return true if the label locations correspond to the overlayOpCode
   */
  static boolean isResultOfOpPoint(OverlayLabel label, int opCode)
  {
    int loc0 = label.getLocation(0);
    int loc1 = label.getLocation(1);
    return isResultOfOp(opCode, loc0, loc1);
  }

  /**
   * Tests whether a point with given {@link Location}s
   * relative to two geometries would be contained in
   * the result of overlaying the geometries using
   * a given overlay operation.
   * This is used to determine whether components
   * computed during the overlay process should be
   * included in the result geometry.
   * <p>
   * The method handles arguments of {@link Location#NONE} correctly.
   *
   * @param overlayOpCode the code for the overlay operation to test
   * @param loc0 the code for the location in the first geometry
   * @param loc1 the code for the location in the second geometry
   *
   * @return true if a point with given locations is in the result of the overlay operation
   */
  static boolean isResultOfOp(int overlayOpCode, int loc0, int loc1)
  {
    if (loc0 == Location.BOUNDARY) loc0 = Location.INTERIOR;
    if (loc1 == Location.BOUNDARY) loc1 = Location.INTERIOR;
    switch (overlayOpCode) {
    case INTERSECTION:
      return loc0 == Location.INTERIOR
          && loc1 == Location.INTERIOR;
    case UNION:
      return loc0 == Location.INTERIOR
          || loc1 == Location.INTERIOR;
    case DIFFERENCE:
      return loc0 == Location.INTERIOR
          && loc1 != Location.INTERIOR;
    case SYMDIFFERENCE:
      return   (     loc0 == Location.INTERIOR &&  loc1 != Location.INTERIOR)
            || (     loc0 != Location.INTERIOR &&  loc1 == Location.INTERIOR);
    }
    return false;
  }

  /**
   * Computes an overlay operation for
   * the given geometry operands, with the
   * noding strategy determined by the precision model.
   *
   * @param geom0 the first geometry argument
   * @param geom1 the second geometry argument
   * @param opCode the code for the desired overlay operation
   * @param pm the precision model to use
   * @return the result of the overlay operation
   */
  public static Geometry overlay(Geometry geom0, Geometry geom1,
      int opCode, PrecisionModel pm)
  {
    OverlayNG ov = new OverlayNG(geom0, geom1, pm, opCode);
    Geometry geomOv = ov.getResult();
    return geomOv;
  }

  /**
   * Computes an overlay operation on the given geometry operands,
   * using a supplied {@link Noder}.
   *
   * @param geom0 the first geometry argument
   * @param geom1 the second geometry argument
   * @param opCode the code for the desired overlay operation
   * @param pm the precision model to use (which may be null if the noder does not use one)
   * @param noder the noder to use
   * @return the result of the overlay operation
   */
  public static Geometry overlay(Geometry geom0, Geometry geom1,
      int opCode, PrecisionModel pm, Noder noder)
  {
    OverlayNG ov = new OverlayNG(geom0, geom1, pm, opCode);
    ov.setNoder(noder);
    Geometry geomOv = ov.getResult();
    return geomOv;
  }

  /**
   * Computes an overlay operation on the given geometry operands,
   * using a supplied {@link Noder}.
   *
   * @param geom0 the first geometry argument
   * @param geom1 the second geometry argument
   * @param opCode the code for the desired overlay operation
   * @param noder the noder to use
   * @return the result of the overlay operation
   */
  public static Geometry overlay(Geometry geom0, Geometry geom1,
      int opCode, Noder noder)
  {
    OverlayNG ov = new OverlayNG(geom0, geom1, null, opCode);
    ov.setNoder(noder);
    Geometry geomOv = ov.getResult();
    return geomOv;
  }

  /**
   * Computes an overlay operation on
   * the given geometry operands,
   * using the precision model of the geometry.
   * and an appropriate noder.
   * <p>
   * The noder is chosen according to the precision model specified.
   * <ul>
   * <li>For {@link PrecisionModel#FIXED}
   * a snap-rounding noder is used, and the computation is robust.
   * <li>For {@link PrecisionModel#FLOATING}
   * a non-snapping noder is used,
   * and this computation may not be robust.
   * If errors occur a {@link TopologyException} is thrown.
   * </ul>
   *
   *
   * @param geom0 the first argument geometry
   * @param geom1 the second argument geometry
   * @param opCode the code for the desired overlay operation
   * @return the result of the overlay operation
   */
  public static Geometry overlay(Geometry geom0, Geometry geom1, int opCode)
  {
    OverlayNG ov = new OverlayNG(geom0, geom1, opCode);
    return ov.getResult();
  }

  /**
   * Computes a union operation on
   * the given geometry, with the supplied precision model.
   * <p>
   * The input must be a valid geometry.
   * Collections must be homogeneous.
   * <p>
   * To union an overlapping set of polygons in a more performant way use {@link UnaryUnionNG}.
   * To union a polyonal coverage or linear network in a more performant way,
   * use {@link CoverageUnion}.
   *
   * @param geom0 the geometry
   * @param pm the precision model to use
   * @return the result of the union operation
   *
   * @see OverlayMixedPoints
   */
  static Geometry union(Geometry geom, PrecisionModel pm)
  {
    OverlayNG ov = new OverlayNG(geom, pm);
    Geometry geomOv = ov.getResult();
    return geomOv;
  }

  /**
   * Computes a union of a single geometry using a custom noder.
   * <p>
   * The primary use of this is to support coverage union.
   * Because of this the overlay is performed using strict mode.
   *
   * @param geom the geometry to union
   * @param pm the precision model to use (maybe be null)
   * @param noder the noder to use
   * @return the result geometry
   *
   * @see CoverageUnion
   */
  static Geometry union(Geometry geom, PrecisionModel pm, Noder noder)
  {
    OverlayNG ov = new OverlayNG(geom, pm);
    ov.setNoder(noder);
    ov.setStrictMode(true);
    Geometry geomOv = ov.getResult();
    return geomOv;
  }

  private int opCode;
  private InputGeometry inputGeom;
  private GeometryFactory geomFact;
  private PrecisionModel pm;
  private Noder noder;
  private boolean isStrictMode = STRICT_MODE_DEFAULT;
  private boolean isOptimized = true;
  private boolean isAreaResultOnly = false;
  private boolean isOutputEdges = false;
  private boolean isOutputResultEdges = false;
  private boolean isOutputNodedEdges = false;

  /**
   * Creates an overlay operation on the given geometries,
   * with a defined precision model.
   * The noding strategy is determined by the precision model.
   *
   * @param geom0 the A operand geometry
   * @param geom1 the B operand geometry (may be null)
   * @param pm the precision model to use
   * @param opCode the overlay opcode
   */
  public OverlayNG(Geometry geom0, Geometry geom1, PrecisionModel pm, int opCode) {
    this.pm = pm;
    this.opCode = opCode;
    geomFact = geom0.getFactory();
    inputGeom = new InputGeometry( geom0, geom1 );
  }

  /**
   * Creates an overlay operation on the given geometries
   * using the precision model of the geometries.
   * <p>
   * The noder is chosen according to the precision model specified.
   * <ul>
   * <li>For {@link PrecisionModel#FIXED}
   * a snap-rounding noder is used, and the computation is robust.
   * <li>For {@link PrecisionModel#FLOATING}
   * a non-snapping noder is used,
   * and this computation may not be robust.
   * If errors occur a {@link TopologyException} is thrown.
   * </ul>
   *
   * @param geom0 the A operand geometry
   * @param geom1 the B operand geometry (may be null)
   * @param opCode the overlay opcode
   */
  public OverlayNG(Geometry geom0, Geometry geom1, int opCode) {
    this(geom0, geom1, geom0.getFactory().getPrecisionModel(), opCode);
  }

  /**
   * Creates a union of a single geometry with a given precision model.
   *
   * @param geom the geometry
   * @param pm the precision model to use
   */
  OverlayNG(Geometry geom, PrecisionModel pm) {
    this(geom, null, pm, UNION);
  }

  /**
   * Sets whether the overlay results are computed according to strict mode
   * semantics.
   * <ul>
   * <li>Lines resulting from topology collapse are not included
   * <li>Result geometry is homogeneous
   *     for the {@link #INTERSECTION} and {@link #DIFFERENCE} operations.
   * <li>Result geometry is homogeneous
   *     for the {@link #UNION} and {@link #SYMDIFFERENCE} operations
   *     if the inputs have the same dimension
   * </ul>
   *
   * @param isStrictMode true if strict mode is to be used
   */
  public void setStrictMode(boolean isStrictMode) {
    this.isStrictMode = isStrictMode;
  }

  /**
   * Sets whether overlay processing optimizations are enabled.
   * It may be useful to disable optimizations
   * for testing purposes.
   * Default is TRUE (optimization enabled).
   *
   * @param isOptimized whether to optimize processing
   */
  public void setOptimized(boolean isOptimized) {
    this.isOptimized = isOptimized;
  }

  /**
   * Sets whether the result can contain only {@link Polygon} components.
   * This is used if it is known that the result must be an (possibly empty) area.
   *
   * @param isAreaResultOnly true if the result should contain only area components
   */
  void setAreaResultOnly(boolean isAreaResultOnly) {
    this.isAreaResultOnly = isAreaResultOnly;
  }

  //------ Testing options -------

  /**
   *
   * @param isOutputEdges
   */
  public void setOutputEdges(boolean isOutputEdges ) {
    this.isOutputEdges = isOutputEdges;
  }

  public void setOutputNodedEdges(boolean isOutputNodedEdges ) {
    this.isOutputEdges = true;
    this.isOutputNodedEdges = isOutputNodedEdges;
  }

  public void setOutputResultEdges(boolean isOutputResultEdges ) {
    this.isOutputResultEdges = isOutputResultEdges;
  }
  //---------------------------------

  void setNoder(Noder noder) {
    this.noder = noder;
  }

  /**
   * Gets the result of the overlay operation.
   *
   * @return the result of the overlay operation.
   *
   * @throws IllegalArgumentException if the input is not supported (e.g. a mixed-dimension geometry)
   * @throws TopologyException if a robustness error occurs
   */
  public Geometry getResult() {
    // handle empty inputs which determine result
    if (OverlayUtil.isEmptyResult(opCode,
        inputGeom.getGeometry(0),
        inputGeom.getGeometry(1),
        pm)) {
      return createEmptyResult();
    }

    /**
     * The elevation model is only computed if the input geometries have Z values.
     */
    ElevationModel elevModel = ElevationModel.create(inputGeom.getGeometry(0), inputGeom.getGeometry(1));
    Geometry result;
    if (inputGeom.isAllPoints()) {
      // handle Point-Point inputs
      result = OverlayPoints.overlay(opCode, inputGeom.getGeometry(0), inputGeom.getGeometry(1), pm);
    }
    else if (! inputGeom.isSingle() &&  inputGeom.hasPoints()) {
      // handle Point-nonPoint inputs
      result = OverlayMixedPoints.overlay(opCode, inputGeom.getGeometry(0), inputGeom.getGeometry(1), pm);
    }
    else {
      // handle case where both inputs are formed of edges (Lines and Polygons)
      result = computeEdgeOverlay();
    }
    /**
     * This is a no-op if the elevation model was not computed due to Z not present
     */
    elevModel.populateZ(result);
    return result;
  }

  private Geometry computeEdgeOverlay() {

    List<Edge> edges = nodeEdges();

    OverlayGraph graph = buildGraph(edges);

    if (isOutputNodedEdges) {
      return OverlayUtil.toLines(graph, isOutputEdges, geomFact);
    }

    labelGraph(graph);
    //for (OverlayEdge e : graph.getEdges()) {  Debug.println(e);  }

    if (isOutputEdges || isOutputResultEdges) {
      return  OverlayUtil.toLines(graph, isOutputEdges, geomFact);
    }

    return extractResult(opCode, graph);
  }

  private List<Edge> nodeEdges() {
    /**
     * Node the edges, using whatever noder is being used
     */
    EdgeNodingBuilder nodingBuilder = new EdgeNodingBuilder(pm, noder);

    /**
     * Optimize Intersection and Difference by clipping to the
     * result extent, if enabled.
     */
    if ( isOptimized ) {
      Envelope clipEnv = OverlayUtil.clippingEnvelope(opCode, inputGeom, pm);
      if (clipEnv != null)
        nodingBuilder.setClipEnvelope( clipEnv );
    }

    List<Edge> mergedEdges = nodingBuilder.build(
        inputGeom.getGeometry(0),
        inputGeom.getGeometry(1));

    /**
     * Record if an input geometry has collapsed.
     * This is used to avoid trying to locate disconnected edges
     * against a geometry which has collapsed completely.
     */
    inputGeom.setCollapsed(0, ! nodingBuilder.hasEdgesFor(0) );
    inputGeom.setCollapsed(1, ! nodingBuilder.hasEdgesFor(1) );

    return mergedEdges;
  }

  private OverlayGraph buildGraph(Collection<Edge> edges) {
    OverlayGraph graph = new OverlayGraph();
    for (Edge e : edges) {
      graph.addEdge(e.getCoordinates(), e.createLabel());
    }
    return graph;
  }

  private void labelGraph(OverlayGraph graph) {
    OverlayLabeller labeller = new OverlayLabeller(graph, inputGeom);
    labeller.computeLabelling();
    labeller.markResultAreaEdges(opCode);
    labeller.unmarkDuplicateEdgesFromResultArea();
  }

  /**
   * Extracts the result geometry components from the fully labelled topology graph.
   * <p>
   * This method implements the semantic that the result of an
   * intersection operation is homogeneous with highest dimension.
   * In other words,
   * if an intersection has components of a given dimension
   * no lower-dimension components are output.
   * For example, if two polygons intersect in an area,
   * no linestrings or points are included in the result,
   * even if portions of the input do meet in lines or points.
   * This semantic choice makes more sense for typical usage,
   * in which only the highest dimension components are of interest.
   *
   * @param opCode the overlay operation
   * @param graph the topology graph
   * @return the result geometry
   */
  private Geometry extractResult(int opCode, OverlayGraph graph) {
    boolean isAllowMixedIntResult = ! isStrictMode;

    //--- Build polygons
    List<OverlayEdge> resultAreaEdges = graph.getResultAreaEdges();
    PolygonBuilder polyBuilder = new PolygonBuilder(resultAreaEdges, geomFact);
    List<Polygon> resultPolyList = polyBuilder.getPolygons();
    boolean hasResultAreaComponents = resultPolyList.size() > 0;

    List<LineString> resultLineList = null;
    List<Point> resultPointList = null;

    if (! isAreaResultOnly) {
      //--- Build lines
      boolean allowResultLines = ! hasResultAreaComponents
          || isAllowMixedIntResult
          || opCode == SYMDIFFERENCE
          || opCode == UNION;
      if ( allowResultLines ) {
        LineBuilder lineBuilder = new LineBuilder(inputGeom, graph, hasResultAreaComponents, opCode, geomFact);
        lineBuilder.setStrictMode(isStrictMode);
        resultLineList = lineBuilder.getLines();
      }
      /**
       * Operations with point inputs are handled elsewhere.
       * Only an Intersection op can produce point results
       * from non-point inputs.
       */
      boolean hasResultComponents = hasResultAreaComponents || resultLineList.size() > 0;
      boolean allowResultPoints = ! hasResultComponents || isAllowMixedIntResult;
      if ( opCode == INTERSECTION && allowResultPoints ) {
        IntersectionPointBuilder pointBuilder = new IntersectionPointBuilder(graph, geomFact);
        pointBuilder.setStrictMode(isStrictMode);
        resultPointList = pointBuilder.getPoints();
      }
    }

    if (isEmpty(resultPolyList)
        && isEmpty(resultLineList)
        && isEmpty(resultPointList))
      return createEmptyResult();

    Geometry resultGeom = OverlayUtil.createResultGeometry(resultPolyList, resultLineList, resultPointList, geomFact);
    return resultGeom;
  }

  private static boolean isEmpty(List list) {
    return list == null || list.size() == 0;
  }

  private Geometry createEmptyResult() {
    return OverlayUtil.createEmptyResult(
        OverlayUtil.resultDimension(opCode,
            inputGeom.getDimension(0),
            inputGeom.getDimension(1)),
          geomFact);
  }

}