common/geometry/S2Cap.java

/*
 * Copyright 2005 Google Inc.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package com.google.common.geometry;


/**
 * This class represents a spherical cap, i.e. a portion of a sphere cut off by
 * a plane. The cap is defined by its axis and height. This representation has
 * good numerical accuracy for very small caps (unlike the (axis,
 * min-distance-from-origin) representation), and is also efficient for
 * containment tests (unlike the (axis, angle) representation).
 *
 * Here are some useful relationships between the cap height (h), the cap
 * opening angle (theta), the maximum chord length from the cap's center (d),
 * and the radius of cap's base (a). All formulas assume a unit radius.
 *
 * h = 1 - cos(theta) = 2 sin^2(theta/2) d^2 = 2 h = a^2 + h^2
 *
 */
public final strictfp class S2Cap implements S2Region {

  /**
   * Multiply a positive number by this constant to ensure that the result of a
   * floating point operation is at least as large as the true
   * infinite-precision result.
   */
  private static final double ROUND_UP = 1.0 + 1.0 / (1L << 52);

  private final S2Point axis;
  private final double height;

  // Caps may be constructed from either an axis and a height, or an axis and
  // an angle. To avoid ambiguity, there are no public constructors
  private S2Cap() {
    axis = new S2Point();
    height = 0;
  }

  private S2Cap(S2Point axis, double height) {
    this.axis = axis;
    this.height = height;
    // assert (isValid());
  }

  /**
   * Create a cap given its axis and the cap height, i.e. the maximum projected
   * distance along the cap axis from the cap center. 'axis' should be a
   * unit-length vector.
   */
  public static S2Cap fromAxisHeight(S2Point axis, double height) {
    // assert (S2.isUnitLength(axis));
    return new S2Cap(axis, height);
  }

  /**
   * Create a cap given its axis and the cap opening angle, i.e. maximum angle
   * between the axis and a point on the cap. 'axis' should be a unit-length
   * vector, and 'angle' should be between 0 and 180 degrees.
   */
  public static S2Cap fromAxisAngle(S2Point axis, S1Angle angle) {
    // The height of the cap can be computed as 1-cos(angle), but this isn't
    // very accurate for angles close to zero (where cos(angle) is almost 1).
    // Computing it as 2*(sin(angle/2)**2) gives much better precision.

    // assert (S2.isUnitLength(axis));
    double d = Math.sin(0.5 * angle.radians());
    return new S2Cap(axis, 2 * d * d);

  }

  /**
   * Create a cap given its axis and its area in steradians. 'axis' should be a
   * unit-length vector, and 'area' should be between 0 and 4 * M_PI.
   */
  public static S2Cap fromAxisArea(S2Point axis, double area) {
    // assert (S2.isUnitLength(axis));
    return new S2Cap(axis, area / (2 * S2.M_PI));
  }

  /** Return an empty cap, i.e. a cap that contains no points. */
  public static S2Cap empty() {
    return new S2Cap(new S2Point(1, 0, 0), -1);
  }

  /** Return a full cap, i.e. a cap that contains all points. */
  public static S2Cap full() {
    return new S2Cap(new S2Point(1, 0, 0), 2);
  }


  // Accessor methods.
  public S2Point axis() {
    return axis;
  }

  public double height() {
    return height;
  }

  public double area() {
    return 2 * S2.M_PI * Math.max(0.0, height);
  }

  /**
   * Return the cap opening angle in radians, or a negative number for empty
   * caps.
   */
  public S1Angle angle() {
    // This could also be computed as acos(1 - height_), but the following
    // formula is much more accurate when the cap height is small. It
    // follows from the relationship h = 1 - cos(theta) = 2 sin^2(theta/2).
    if (isEmpty()) {
      return S1Angle.radians(-1);
    }
    return S1Angle.radians(2 * Math.asin(Math.sqrt(0.5 * height)));
  }

  /**
   * We allow negative heights (to represent empty caps) but not heights greater
   * than 2.
   */
  public boolean isValid() {
    return S2.isUnitLength(axis) && height <= 2;
  }

  /** Return true if the cap is empty, i.e. it contains no points. */
  public boolean isEmpty() {
    return height < 0;
  }

  /** Return true if the cap is full, i.e. it contains all points. */
  public boolean isFull() {
    return height >= 2;
  }

  /**
   * Return the complement of the interior of the cap. A cap and its complement
   * have the same boundary but do not share any interior points. The complement
   * operator is not a bijection, since the complement of a singleton cap
   * (containing a single point) is the same as the complement of an empty cap.
   */
  public S2Cap complement() {
    // The complement of a full cap is an empty cap, not a singleton.
    // Also make sure that the complement of an empty cap has height 2.
    double cHeight = isFull() ? -1 : 2 - Math.max(height, 0.0);
    return S2Cap.fromAxisHeight(S2Point.neg(axis), cHeight);
  }

  /**
   * Return true if and only if this cap contains the given other cap (in a set
   * containment sense, e.g. every cap contains the empty cap).
   */
  public boolean contains(S2Cap other) {
    if (isFull() || other.isEmpty()) {
      return true;
    }
    return angle().radians() >= axis.angle(other.axis)
      + other.angle().radians();
  }

  /**
   * Return true if and only if the interior of this cap intersects the given
   * other cap. (This relationship is not symmetric, since only the interior of
   * this cap is used.)
   */
  public boolean interiorIntersects(S2Cap other) {
    // Interior(X) intersects Y if and only if Complement(Interior(X))
    // does not contain Y.
    return !complement().contains(other);
  }

  /**
   * Return true if and only if the given point is contained in the interior of
   * the region (i.e. the region excluding its boundary). 'p' should be a
   * unit-length vector.
   */
  public boolean interiorContains(S2Point p) {
    // assert (S2.isUnitLength(p));
    return isFull() || S2Point.sub(axis, p).norm2() < 2 * height;
  }

  /**
   * Increase the cap height if necessary to include the given point. If the cap
   * is empty the axis is set to the given point, but otherwise it is left
   * unchanged. 'p' should be a unit-length vector.
   */
  public S2Cap addPoint(S2Point p) {
    // Compute the squared chord length, then convert it into a height.
    // assert (S2.isUnitLength(p));
    if (isEmpty()) {
      return new S2Cap(p, 0);
    } else {
      // To make sure that the resulting cap actually includes this point,
      // we need to round up the distance calculation. That is, after
      // calling cap.AddPoint(p), cap.Contains(p) should be true.
      double dist2 = S2Point.sub(axis, p).norm2();
      double newHeight = Math.max(height, ROUND_UP * 0.5 * dist2);
      return new S2Cap(axis, newHeight);
    }
  }

  // Increase the cap height if necessary to include "other". If the current
  // cap is empty it is set to the given other cap.
  public S2Cap addCap(S2Cap other) {
    if (isEmpty()) {
      return new S2Cap(other.axis, other.height);
    } else {
      // See comments for FromAxisAngle() and AddPoint(). This could be
      // optimized by doing the calculation in terms of cap heights rather
      // than cap opening angles.
      double angle = axis.angle(other.axis) + other.angle().radians();
      if (angle >= S2.M_PI) {
        return new S2Cap(axis, 2); //Full cap
      } else {
        double d = Math.sin(0.5 * angle);
        double newHeight = Math.max(height, ROUND_UP * 2 * d * d);
        return new S2Cap(axis, newHeight);
      }
    }
  }

  // //////////////////////////////////////////////////////////////////////
  // S2Region interface (see {@code S2Region} for details):
  @Override
  public S2Cap getCapBound() {
    return this;
  }

  @Override
  public S2LatLngRect getRectBound() {
    if (isEmpty()) {
      return S2LatLngRect.empty();
    }

    // Convert the axis to a (lat,lng) pair, and compute the cap angle.
    S2LatLng axisLatLng = new S2LatLng(axis);
    double capAngle = angle().radians();

    boolean allLongitudes = false;
    double[] lat = new double[2], lng = new double[2];
    lng[0] = -S2.M_PI;
    lng[1] = S2.M_PI;

    // Check whether cap includes the south pole.
    lat[0] = axisLatLng.lat().radians() - capAngle;
    if (lat[0] <= -S2.M_PI_2) {
      lat[0] = -S2.M_PI_2;
      allLongitudes = true;
    }
    // Check whether cap includes the north pole.
    lat[1] = axisLatLng.lat().radians() + capAngle;
    if (lat[1] >= S2.M_PI_2) {
      lat[1] = S2.M_PI_2;
      allLongitudes = true;
    }
    if (!allLongitudes) {
      // Compute the range of longitudes covered by the cap. We use the law
      // of sines for spherical triangles. Consider the triangle ABC where
      // A is the north pole, B is the center of the cap, and C is the point
      // of tangency between the cap boundary and a line of longitude. Then
      // C is a right angle, and letting a,b,c denote the sides opposite A,B,C,
      // we have sin(a)/sin(A) = sin(c)/sin(C), or sin(A) = sin(a)/sin(c).
      // Here "a" is the cap angle, and "c" is the colatitude (90 degrees
      // minus the latitude). This formula also works for negative latitudes.
      //
      // The formula for sin(a) follows from the relationship h = 1 - cos(a).

      double sinA = Math.sqrt(height * (2 - height));
      double sinC = Math.cos(axisLatLng.lat().radians());
      if (sinA <= sinC) {
        double angleA = Math.asin(sinA / sinC);
        lng[0] = Math.IEEEremainder(axisLatLng.lng().radians() - angleA,
          2 * S2.M_PI);
        lng[1] = Math.IEEEremainder(axisLatLng.lng().radians() + angleA,
          2 * S2.M_PI);
      }
    }
    return new S2LatLngRect(new R1Interval(lat[0], lat[1]), new S1Interval(
      lng[0], lng[1]));
  }

  @Override
  public boolean contains(S2Cell cell) {
    // If the cap does not contain all cell vertices, return false.
    // We check the vertices before taking the Complement() because we can't
    // accurately represent the complement of a very small cap (a height
    // of 2-epsilon is rounded off to 2).
    S2Point[] vertices = new S2Point[4];
    for (int k = 0; k < 4; ++k) {
      vertices[k] = cell.getVertex(k);
      if (!contains(vertices[k])) {
        return false;
      }
    }
    // Otherwise, return true if the complement of the cap does not intersect
    // the cell. (This test is slightly conservative, because technically we
    // want Complement().InteriorIntersects() here.)
    return !complement().intersects(cell, vertices);
  }

  @Override
  public boolean mayIntersect(S2Cell cell) {
    // If the cap contains any cell vertex, return true.
    S2Point[] vertices = new S2Point[4];
    for (int k = 0; k < 4; ++k) {
      vertices[k] = cell.getVertex(k);
      if (contains(vertices[k])) {
        return true;
      }
    }
    return intersects(cell, vertices);
  }

  /**
   * Return true if the cap intersects 'cell', given that the cap vertices have
   * alrady been checked.
   */
  public boolean intersects(S2Cell cell, S2Point[] vertices) {
    // Return true if this cap intersects any point of 'cell' excluding its
    // vertices (which are assumed to already have been checked).

    // If the cap is a hemisphere or larger, the cell and the complement of the
    // cap are both convex. Therefore since no vertex of the cell is contained,
    // no other interior point of the cell is contained either.
    if (height >= 1) {
      return false;
    }

    // We need to check for empty caps due to the axis check just below.
    if (isEmpty()) {
      return false;
    }

    // Optimization: return true if the cell contains the cap axis. (This
    // allows half of the edge checks below to be skipped.)
    if (cell.contains(axis)) {
      return true;
    }

    // At this point we know that the cell does not contain the cap axis,
    // and the cap does not contain any cell vertex. The only way that they
    // can intersect is if the cap intersects the interior of some edge.

    double sin2Angle = height * (2 - height); // sin^2(capAngle)
    for (int k = 0; k < 4; ++k) {
      S2Point edge = cell.getEdgeRaw(k);
      double dot = axis.dotProd(edge);
      if (dot > 0) {
        // The axis is in the interior half-space defined by the edge. We don't
        // need to consider these edges, since if the cap intersects this edge
        // then it also intersects the edge on the opposite side of the cell
        // (because we know the axis is not contained with the cell).
        continue;
      }
      // The Norm2() factor is necessary because "edge" is not normalized.
      if (dot * dot > sin2Angle * edge.norm2()) {
        return false; // Entire cap is on the exterior side of this edge.
      }
      // Otherwise, the great circle containing this edge intersects
      // the interior of the cap. We just need to check whether the point
      // of closest approach occurs between the two edge endpoints.
      S2Point dir = S2Point.crossProd(edge, axis);
      if (dir.dotProd(vertices[k]) < 0
          && dir.dotProd(vertices[(k + 1) & 3]) > 0) {
        return true;
      }
    }
    return false;
  }

  public boolean contains(S2Point p) {
    // The point 'p' should be a unit-length vector.
    // assert (S2.isUnitLength(p));
    return S2Point.sub(axis, p).norm2() <= 2 * height;

  }


  /** Return true if two caps are identical. */
  @Override
  public boolean equals(Object that) {

    if (!(that instanceof S2Cap)) {
      return false;
    }

    S2Cap other = (S2Cap) that;
    return (axis.equals(other.axis) && height == other.height)
        || (isEmpty() && other.isEmpty()) || (isFull() && other.isFull());

  }

  @Override
  public int hashCode() {
    if (isFull()) {
      return 17;
    } else if (isEmpty()) {
      return 37;
    }
    int result = 17;
    result = 37 * result + axis.hashCode();
    long heightBits = Double.doubleToLongBits(height);
    result = 37 * result + (int) ((heightBits >>> 32) ^ heightBits);
    return result;
  }

  // /////////////////////////////////////////////////////////////////////
  // The following static methods are convenience functions for assertions
  // and testing purposes only.

  /**
   * Return true if the cap axis and height differ by at most "max_error" from
   * the given cap "other".
   */
  boolean approxEquals(S2Cap other, double maxError) {
    return (axis.aequal(other.axis, maxError) && Math.abs(height - other.height) <= maxError)
      || (isEmpty() && other.height <= maxError)
      || (other.isEmpty() && height <= maxError)
      || (isFull() && other.height >= 2 - maxError)
      || (other.isFull() && height >= 2 - maxError);
  }

  boolean approxEquals(S2Cap other) {
    return approxEquals(other, 1e-14);
  }

  @Override
  public String toString() {
    return "[Point = " + axis.toString() + " Height = " + height + "]";
  }
}