Java tutorial
/* * 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 + "]"; } }