GeoCircle.java
/*******************************************************************************
* Copyright (c) 2015 Voyager Search and MITRE
* All rights reserved. This program and the accompanying materials
* are made available under the terms of the Apache License, Version 2.0 which
* accompanies this distribution and is available at
* http://www.apache.org/licenses/LICENSE-2.0.txt
******************************************************************************/
package org.locationtech.spatial4j.shape.impl;
import org.locationtech.spatial4j.context.SpatialContext;
import org.locationtech.spatial4j.distance.DistanceUtils;
import org.locationtech.spatial4j.shape.Point;
import org.locationtech.spatial4j.shape.Rectangle;
import org.locationtech.spatial4j.shape.SpatialRelation;
import java.util.Formatter;
import java.util.Locale;
/**
* A circle as it exists on the surface of a sphere.
*/
public class GeoCircle extends CircleImpl {
private GeoCircle inverseCircle;//when distance reaches > 1/2 way around the world, cache the inverse.
private double horizAxisY;//see getYAxis
public GeoCircle(Point p, double radiusDEG, SpatialContext ctx) {
super(p, radiusDEG, ctx);
assert ctx.isGeo();
init();
}
@Override
public void reset(double x, double y, double radiusDEG) {
super.reset(x, y, radiusDEG);
init();
}
private void init() {
if (radiusDEG > 90) {
//--spans more than half the globe
assert enclosingBox.getWidth() == 360;
double backDistDEG = 180 - radiusDEG;
if (backDistDEG > 0) {
double backRadius = 180 - radiusDEG;
double backX = DistanceUtils.normLonDEG(getCenter().getX() + 180);
double backY = DistanceUtils.normLatDEG(getCenter().getY() + 180);
//Shrink inverseCircle as small as possible to avoid accidental overlap.
// Note that this is tricky business to come up with a value small enough
// but not too small or else numerical conditioning issues become a problem.
backRadius -= Math.max(Math.ulp(Math.abs(backY)+backRadius), Math.ulp(Math.abs(backX)+backRadius));
if (inverseCircle != null) {
inverseCircle.reset(backX, backY, backRadius);
} else {
inverseCircle = new GeoCircle(ctx.makePoint(backX, backY), backRadius, ctx);
}
} else {
inverseCircle = null;//whole globe
}
horizAxisY = getCenter().getY();//although probably not used
} else {
inverseCircle = null;
double _horizAxisY = ctx.getDistCalc().calcBoxByDistFromPt_yHorizAxisDEG(getCenter(), radiusDEG, ctx);
//some rare numeric conditioning cases can cause this to be barely beyond the box
if (_horizAxisY > enclosingBox.getMaxY()) {
horizAxisY = enclosingBox.getMaxY();
} else if (_horizAxisY < enclosingBox.getMinY()) {
horizAxisY = enclosingBox.getMinY();
} else {
horizAxisY = _horizAxisY;
}
//assert enclosingBox.relate_yRange(horizAxis,horizAxis,ctx).intersects();
}
}
@Override
protected double getYAxis() {
return horizAxisY;
}
/**
* Called after bounding box is intersected.
* @param bboxSect INTERSECTS or CONTAINS from enclosingBox's intersection
* @return DISJOINT, CONTAINS, or INTERSECTS (not WITHIN)
*/
@Override
protected SpatialRelation relateRectanglePhase2(Rectangle r, SpatialRelation bboxSect) {
if (inverseCircle != null) {
return inverseCircle.relate(r).inverse();
}
//if a pole is wrapped, we have a separate algorithm
if (enclosingBox.getWidth() == 360) {
return relateRectangleCircleWrapsPole(r, ctx);
}
//This is an optimization path for when there are no dateline or pole issues.
if (!enclosingBox.getCrossesDateLine() && !r.getCrossesDateLine()) {
return super.relateRectanglePhase2(r, bboxSect);
}
//Rectangle wraps around the world longitudinally creating a solid band; there are no corners to test intersection
if (r.getWidth() == 360) {
return SpatialRelation.INTERSECTS;
}
//do quick check to see if all corners are within this circle for CONTAINS
int cornersIntersect = numCornersIntersect(r);
if (cornersIntersect == 4) {
//ensure r's x axis is within c's. If it doesn't, r sneaks around the globe to touch the other side (intersect).
SpatialRelation xIntersect = r.relateXRange(enclosingBox.getMinX(), enclosingBox.getMaxX());
if (xIntersect == SpatialRelation.WITHIN)
return SpatialRelation.CONTAINS;
return SpatialRelation.INTERSECTS;
}
//INTERSECT or DISJOINT ?
if (cornersIntersect > 0)
return SpatialRelation.INTERSECTS;
//Now we check if one of the axis of the circle intersect with r. If so we have
// intersection.
/* x axis intersects */
if ( r.relateYRange(getYAxis(), getYAxis()).intersects() // at y vertical
&& r.relateXRange(enclosingBox.getMinX(), enclosingBox.getMaxX()).intersects() )
return SpatialRelation.INTERSECTS;
/* y axis intersects */
if (r.relateXRange(getXAxis(), getXAxis()).intersects()) { // at x horizontal
double yTop = getCenter().getY()+ radiusDEG;
assert yTop <= 90;
double yBot = getCenter().getY()- radiusDEG;
assert yBot >= -90;
if (r.relateYRange(yBot, yTop).intersects())//back bottom
return SpatialRelation.INTERSECTS;
}
return SpatialRelation.DISJOINT;
}
private SpatialRelation relateRectangleCircleWrapsPole(Rectangle r, SpatialContext ctx) {
//This method handles the case where the circle wraps ONE pole, but not both. For both,
// there is the inverseCircle case handled before now. The only exception is for the case where
// the circle covers the entire globe, and we'll check that first.
if (radiusDEG == 180)//whole globe
return SpatialRelation.CONTAINS;
//Check if r is within the pole wrap region:
double yTop = getCenter().getY() + radiusDEG;
if (yTop > 90) {
double yTopOverlap = yTop - 90;
assert yTopOverlap <= 90;
if (r.getMinY() >= 90 - yTopOverlap)
return SpatialRelation.CONTAINS;
} else {
double yBot = point.getY() - radiusDEG;
if (yBot < -90) {
double yBotOverlap = -90 - yBot;
assert yBotOverlap <= 90;
if (r.getMaxY() <= -90 + yBotOverlap)
return SpatialRelation.CONTAINS;
} else {
//This point is probably not reachable ??
assert yTop == 90 || yBot == -90;//we simply touch a pole
//continue
}
}
//If there are no corners to check intersection because r wraps completely...
if (r.getWidth() == 360)
return SpatialRelation.INTERSECTS;
//Check corners:
int cornersIntersect = numCornersIntersect(r);
// (It might be possible to reduce contains() calls within nCI() to exactly two, but this intersection
// code is complicated enough as it is.)
double frontX = getCenter().getX();
if (cornersIntersect == 4) {//all
double backX = frontX <= 0 ? frontX + 180 : frontX - 180;
if (r.relateXRange(backX, backX).intersects())
return SpatialRelation.INTERSECTS;
else
return SpatialRelation.CONTAINS;
} else if (cornersIntersect == 0) {//none
if (r.relateXRange(frontX, frontX).intersects())
return SpatialRelation.INTERSECTS;
else
return SpatialRelation.DISJOINT;
} else//partial
return SpatialRelation.INTERSECTS;
}
/** Returns either 0 for none, 1 for some, or 4 for all. */
private int numCornersIntersect(Rectangle r) {
//We play some logic games to avoid calling contains() which can be expensive.
boolean bool;//if true then all corners intersect, if false then no corners intersect
// for partial, we exit early with 1 and ignore bool.
bool = (contains(r.getMinX(),r.getMinY()));
if (contains(r.getMinX(),r.getMaxY())) {
if (!bool)
return 1;//partial
} else {
if (bool)
return 1;//partial
}
if (contains(r.getMaxX(),r.getMinY())) {
if (!bool)
return 1;//partial
} else {
if (bool)
return 1;//partial
}
if (contains(r.getMaxX(),r.getMaxY())) {
if (!bool)
return 1;//partial
} else {
if (bool)
return 1;//partial
}
return bool?4:0;
}
@Override
public String toString() {
//Add distance in km, which may be easier to recognize.
double distKm = DistanceUtils.degrees2Dist(radiusDEG, DistanceUtils.EARTH_MEAN_RADIUS_KM);
//instead of String.format() so that we get consistent output no matter the locale
String dStr = new Formatter(Locale.ROOT).format("%.1f\u00B0 %.2fkm", radiusDEG, distKm).toString();
return "Circle(" + point + ", d=" + dStr + ')';
}
}