com.google.common.geometry.S2Cap Maven / Gradle / Ivy
The newest version!
/*
* 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 + "]";
}
}
© 2015 - 2025 Weber Informatics LLC | Privacy Policy