path: root/src/mongo/db/geo/shapes.cpp

/**
*    Copyright (C) 2008-2012 10gen Inc.
*
*    This program is free software: you can redistribute it and/or  modify
*    it under the terms of the GNU Affero General Public License, version 3,
*    as published by the Free Software Foundation.
*
*    This program is distributed in the hope that it will be useful,
*    but WITHOUT ANY WARRANTY; without even the implied warranty of
*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*    GNU Affero General Public License for more details.
*
*    You should have received a copy of the GNU Affero General Public License
*    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*
*    As a special exception, the copyright holders give permission to link the
*    code of portions of this program with the OpenSSL library under certain
*    conditions as described in each individual source file and distribute
*    linked combinations including the program with the OpenSSL library. You
*    must comply with the GNU Affero General Public License in all respects for
*    all of the code used other than as permitted herein. If you modify file(s)
*    with this exception, you may extend this exception to your version of the
*    file(s), but you are not obligated to do so. If you do not wish to do so,
*    delete this exception statement from your version. If you delete this
*    exception statement from all source files in the program, then also delete
*    it in the license file.
*/

#include "mongo/db/jsobj.h"
#include "mongo/db/geo/shapes.h"
#include "mongo/util/mongoutils/str.h"

using std::abs;

// So we can get at the str namespace.
using namespace mongoutils;

namespace mongo {

////////////// Point

Point::Point() : x(0), y(0) {}

Point::Point(double x, double y) : x(x), y(y) {}

Point::Point(const BSONElement& e) {
    BSONObjIterator i(e.Obj());
    x = i.next().number();
    y = i.next().number();
}

Point::Point(const BSONObj& o) {
    BSONObjIterator i(o);
    x = i.next().number();
    y = i.next().number();
}

string Point::toString() const {
    StringBuilder buf;
    buf << "(" << x << "," << y << ")";
    return buf.str();
}

////////////// Circle

Circle::Circle() {}
Circle::Circle(double radius, Point center) : radius(radius), center(center) {}


////////////// Box

Box::Box() {}

Box::Box(double x, double y, double size) : _min(x, y), _max(x + size, y + size) {}

Box::Box(const Point& ptA, const Point& ptB) {
    init(ptA, ptB);
}

void Box::init(const Point& ptA, const Point& ptB) {
    _min.x = min(ptA.x, ptB.x);
    _min.y = min(ptA.y, ptB.y);
    _max.x = max(ptA.x, ptB.x);
    _max.y = max(ptA.y, ptB.y);
}

void Box::init(const Box& other) {
    init(other._min, other._max);
}

BSONArray Box::toBSON() const {
    return BSON_ARRAY(BSON_ARRAY(_min.x << _min.y) << BSON_ARRAY(_max.x << _max.y));
}

string Box::toString() const {
    StringBuilder buf;
    buf << _min.toString() << " -->> " << _max.toString();
    return buf.str();
}

bool Box::between(double min, double max, double val, double fudge) const {
    return val + fudge >= min && val <= max + fudge;
}

bool Box::onBoundary(double bound, double val, double fudge) const {
    return (val >= bound - fudge && val <= bound + fudge);
}

bool Box::mid(double amin, double amax, double bmin, double bmax, bool min, double* res) const {
    verify(amin <= amax);
    verify(bmin <= bmax);

    if (amin < bmin) {
        if (amax < bmin)
            return false;
        *res = min ? bmin : amax;
        return true;
    }
    if (amin > bmax)
        return false;
    *res = min ? amin : bmax;
    return true;
}

bool Box::intersects(const Box& other) const {
    bool intersectX = between(_min.x, _max.x, other._min.x)  // contain part of other range
        || between(_min.x, _max.x, other._max.x)             // contain part of other range
        || between(other._min.x, other._max.x, _min.x);      // other range contains us

    bool intersectY = between(_min.y, _max.y, other._min.y) ||
        between(_min.y, _max.y, other._max.y) || between(other._min.y, other._max.y, _min.y);

    return intersectX && intersectY;
}

double Box::legacyIntersectFraction(const Box& other) const {
    Point boundMin(0, 0);
    Point boundMax(0, 0);

    if (!mid(_min.x, _max.x, other._min.x, other._max.x, true, &boundMin.x) ||
        !mid(_min.x, _max.x, other._min.x, other._max.x, false, &boundMax.x) ||
        !mid(_min.y, _max.y, other._min.y, other._max.y, true, &boundMin.y) ||
        !mid(_min.y, _max.y, other._min.y, other._max.y, false, &boundMax.y))
        return 0;

    Box intersection(boundMin, boundMax);
    return intersection.area() / area();
}

double Box::area() const {
    return (_max.x - _min.x) * (_max.y - _min.y);
}

double Box::maxDim() const {
    return max(_max.x - _min.x, _max.y - _min.y);
}

Point Box::center() const {
    return Point((_min.x + _max.x) / 2, (_min.y + _max.y) / 2);
}

void Box::truncate(double min, double max) {
    if (_min.x < min)
        _min.x = min;
    if (_min.y < min)
        _min.y = min;
    if (_max.x > max)
        _max.x = max;
    if (_max.y > max)
        _max.y = max;
}

void Box::fudge(double error) {
    _min.x -= error;
    _min.y -= error;
    _max.x += error;
    _max.y += error;
}

void Box::expandToInclude(const Point& pt) {
    _min.x = min(_min.x, pt.x);
    _min.y = min(_min.y, pt.y);
    _max.x = max(_max.x, pt.x);
    _max.y = max(_max.y, pt.y);
}

bool Box::onBoundary(Point p, double fudge) const {
    return onBoundary(_min.x, p.x, fudge) || onBoundary(_max.x, p.x, fudge) ||
        onBoundary(_min.y, p.y, fudge) || onBoundary(_max.y, p.y, fudge);
}

bool Box::inside(Point p, double fudge) const {
    bool res = inside(p.x, p.y, fudge);
    return res;
}

bool Box::inside(double x, double y, double fudge) const {
    return between(_min.x, _max.x, x, fudge) && between(_min.y, _max.y, y, fudge);
}

bool Box::contains(const Box& other, double fudge) const {
    return inside(other._min, fudge) && inside(other._max, fudge);
}

////////////// Polygon

Polygon::Polygon() {}

Polygon::Polygon(const vector<Point>& points) {
    init(points);
}

void Polygon::init(const vector<Point>& points) {
    _points.clear();
    _bounds.reset();
    _centroid.reset();

    _points.insert(_points.begin(), points.begin(), points.end());
}

void Polygon::init(const Polygon& other) {
    init(other._points);
}

int Polygon::size(void) const {
    return _points.size();
}

bool Polygon::contains(const Point& p) const {
    return contains(p, 0) > 0;
}

/*
 * Return values:
 * -1 if no intersection
 * 0 if maybe an intersection (using fudge)
 * 1 if there is an intersection
 *
 * A ray casting intersection method is used.
 */
int Polygon::contains(const Point& p, double fudge) const {
    Box fudgeBox(Point(p.x - fudge, p.y - fudge), Point(p.x + fudge, p.y + fudge));

    int counter = 0;
    Point p1 = _points[0];
    for (int i = 1; i <= size(); i++) {
        // XXX: why is there a mod here?
        Point p2 = _points[i % size()];

        // We need to check whether or not this segment intersects our error box
        if (fudge > 0 &&
            // Points not too far below box
            fudgeBox._min.y <= std::max(p1.y, p2.y) &&
            // Points not too far above box
            fudgeBox._max.y >= std::min(p1.y, p2.y) &&
            // Points not too far to left of box
            fudgeBox._min.x <= std::max(p1.x, p2.x) &&
            // Points not too far to right of box
            fudgeBox._max.x >= std::min(p1.x, p2.x)) {
            // If our box contains one or more of these points, we need to do an exact
            // check.
            if (fudgeBox.inside(p1)) {
                return 0;
            }
            if (fudgeBox.inside(p2)) {
                return 0;
            }

            // Do intersection check for vertical sides
            if (p1.y != p2.y) {
                double invSlope = (p2.x - p1.x) / (p2.y - p1.y);

                double xintersT = (fudgeBox._max.y - p1.y) * invSlope + p1.x;
                if (fudgeBox._min.x <= xintersT && fudgeBox._max.x >= xintersT) {
                    return 0;
                }

                double xintersB = (fudgeBox._min.y - p1.y) * invSlope + p1.x;
                if (fudgeBox._min.x <= xintersB && fudgeBox._max.x >= xintersB) {
                    return 0;
                }
            }

            // Do intersection check for horizontal sides
            if (p1.x != p2.x) {
                double slope = (p2.y - p1.y) / (p2.x - p1.x);

                double yintersR = (p1.x - fudgeBox._max.x) * slope + p1.y;
                if (fudgeBox._min.y <= yintersR && fudgeBox._max.y >= yintersR) {
                    return 0;
                }

                double yintersL = (p1.x - fudgeBox._min.x) * slope + p1.y;
                if (fudgeBox._min.y <= yintersL && fudgeBox._max.y >= yintersL) {
                    return 0;
                }
            }
        } else if (fudge == 0) {
            // If this is an exact vertex, we won't intersect, so check this
            if (p.y == p1.y && p.x == p1.x)
                return 1;
            else if (p.y == p2.y && p.x == p2.x)
                return 1;

            // If this is a horizontal line we won't intersect, so check this
            if (p1.y == p2.y && p.y == p1.y) {
                // Check that the x-coord lies in the line
                if (p.x >= std::min(p1.x, p2.x) && p.x <= std::max(p1.x, p2.x))
                    return 1;
            }
        }

        // Normal intersection test.
        // TODO: Invert these for clearer logic?
        if (p.y > std::min(p1.y, p2.y)) {
            if (p.y <= std::max(p1.y, p2.y)) {
                if (p.x <= std::max(p1.x, p2.x)) {
                    if (p1.y != p2.y) {
                        double xinters = (p.y - p1.y) * (p2.x - p1.x) / (p2.y - p1.y) + p1.x;
                        // Special case of point on vertical line
                        if (p1.x == p2.x && p.x == p1.x) {
                            // Need special case for the vertical edges, for example:
                            // 1) \e   pe/----->
                            // vs.
                            // 2) \ep---e/----->
                            //
                            // if we count exact as intersection, then 1 is in but 2 is out
                            // if we count exact as no-int then 1 is out but 2 is in.

                            return 1;
                        } else if (p1.x == p2.x || p.x <= xinters) {
                            counter++;
                        }
                    }
                }
            }
        }

        p1 = p2;
    }

    if (counter % 2 == 0) {
        return -1;
    } else {
        return 1;
    }
}

const Point& Polygon::centroid() const {
    if (_centroid) {
        return *_centroid;
    }

    _centroid.reset(new Point());

    double signedArea = 0.0;
    double area = 0.0;  // Partial signed area

    /// For all vertices except last
    int i = 0;
    for (i = 0; i < size() - 1; ++i) {
        area = _points[i].x * _points[i + 1].y - _points[i + 1].x * _points[i].y;
        signedArea += area;
        _centroid->x += (_points[i].x + _points[i + 1].x) * area;
        _centroid->y += (_points[i].y + _points[i + 1].y) * area;
    }

    // Do last vertex
    area = _points[i].x * _points[0].y - _points[0].x * _points[i].y;
    _centroid->x += (_points[i].x + _points[0].x) * area;
    _centroid->y += (_points[i].y + _points[0].y) * area;
    signedArea += area;
    signedArea *= 0.5;
    _centroid->x /= (6 * signedArea);
    _centroid->y /= (6 * signedArea);

    return *_centroid;
}

const Box& Polygon::bounds() const {
    if (_bounds) {
        return *_bounds;
    }

    _bounds.reset(new Box(_points[0], _points[0]));

    for (int i = 1; i < size(); i++) {
        _bounds->expandToInclude(_points[i]);
    }

    return *_bounds;
}

R2Annulus::R2Annulus() : _inner(0.0), _outer(0.0) {}

R2Annulus::R2Annulus(const Point& center, double inner, double outer)
    : _center(center), _inner(inner), _outer(outer) {}

const Point& R2Annulus::center() const {
    return _center;
}

double R2Annulus::getInner() const {
    return _inner;
}

double R2Annulus::getOuter() const {
    return _outer;
}

bool R2Annulus::contains(const Point& point) const {
    // See if we're inside the inner radius
    if (distanceCompare(point, _center, _inner) < 0) {
        return false;
    }

    // See if we're outside the outer radius
    if (distanceCompare(point, _center, _outer) > 0) {
        return false;
    }

    return true;
}

Box R2Annulus::getR2Bounds() const {
    return Box(
        _center.x - _outer, _center.y - _outer, 2 * _outer);  // Box(_min.x, _min.y, edgeLength)
}

bool R2Annulus::fastContains(const Box& other) const {
    return circleContainsBox(Circle(_outer, _center), other) &&
        !circleInteriorIntersectsWithBox(Circle(_inner, _center), other);
}

bool R2Annulus::fastDisjoint(const Box& other) const {
    return !circleIntersectsWithBox(Circle(_outer, _center), other) ||
        circleInteriorContainsBox(Circle(_inner, _center), other);
}

string R2Annulus::toString() const {
    return str::stream() << "center: " << _center.toString() << " inner: " << _inner
                         << " outer: " << _outer;
}

/////// Other methods

double S2Distance::distanceRad(const S2Point& pointA, const S2Point& pointB) {
    S1Angle angle(pointA, pointB);
    return angle.radians();
}

double S2Distance::minDistanceRad(const S2Point& point, const S2Polyline& line) {
    int tmp;
    S1Angle angle(point, line.Project(point, &tmp));
    return angle.radians();
}

double S2Distance::minDistanceRad(const S2Point& point, const S2Polygon& polygon) {
    S1Angle angle(point, polygon.Project(point));
    return angle.radians();
}

double S2Distance::minDistanceRad(const S2Point& point, const S2Cap& cap) {
    S1Angle angleToCenter(point, cap.axis());
    return (angleToCenter - cap.angle()).radians();
}

/**
 * Distance method that compares x or y coords when other direction is zero,
 * avoids numerical error when distances are very close to radius but axis-aligned.
 *
 * An example of the problem is:
 * (52.0 - 51.9999) - 0.0001 = 3.31965e-15 and 52.0 - 51.9999 > 0.0001 in double arithmetic
 * but:
 * 51.9999 + 0.0001 <= 52.0
 *
 * This avoids some (but not all!) suprising results in $center queries where points are
 * (radius + center.x, center.y) or vice-versa.
 */
bool distanceWithin(const Point& p1, const Point& p2, double radius) {
    return distanceCompare(p1, p2, radius) <= 0.0;
}

// Compare the distance between p1 and p2 with the radius.
// Float-number comparison might be inaccurate.
//
// > 0: distance is greater than radius
// = 0: distance equals radius
// < 0: distance is less than radius
double distanceCompare(const Point& p1, const Point& p2, double radius) {
    double a = p2.x - p1.x;
    double b = p2.y - p1.y;

    if (a == 0) {
        //
        // Note:  For some, unknown reason, when a 32-bit g++ optimizes this call, the sum is
        // calculated imprecisely.  We need to force the compiler to always evaluate it
        // correctly, hence the weirdness.
        //
        // On some 32-bit linux machines, removing the volatile keyword or calculating the sum
        // inline will make certain geo tests fail.  Of course this check will force volatile
        // for all 32-bit systems, not just affected systems.
        if (sizeof(void*) <= 4) {
            volatile double sum = p2.y > p1.y ? p1.y + radius : p2.y + radius;
            return p2.y > p1.y ? p2.y - sum : p1.y - sum;
        } else {
            // Original math, correct for most systems
            return p2.y > p1.y ? p2.y - (p1.y + radius) : p1.y - (p2.y + radius);
        }
    }

    if (b == 0) {
        if (sizeof(void*) <= 4) {
            volatile double sum = p2.x > p1.x ? p1.x + radius : p2.x + radius;
            return p2.x > p1.x ? p2.x - sum : p1.x - sum;
        } else {
            return p2.x > p1.x ? p2.x - (p1.x + radius) : p1.x - (p2.x + radius);
        }
    }

    return sqrt((a * a) + (b * b)) - radius;
}

// note: multiply by earth radius for distance
double spheredist_rad(const Point& p1, const Point& p2) {
    // this uses the n-vector formula: http://en.wikipedia.org/wiki/N-vector
    // If you try to match the code to the formula, note that I inline the cross-product.

    double sinx1(sin(p1.x)), cosx1(cos(p1.x));
    double siny1(sin(p1.y)), cosy1(cos(p1.y));
    double sinx2(sin(p2.x)), cosx2(cos(p2.x));
    double siny2(sin(p2.y)), cosy2(cos(p2.y));

    double cross_prod =
        (cosy1 * cosx1 * cosy2 * cosx2) + (cosy1 * sinx1 * cosy2 * sinx2) + (siny1 * siny2);

    if (cross_prod >= 1 || cross_prod <= -1) {
        // fun with floats
        verify(fabs(cross_prod) - 1 < 1e-6);
        return cross_prod > 0 ? 0 : M_PI;
    }

    return acos(cross_prod);
}

// @param p1 A point on the sphere where x and y are degrees.
// @param p2 A point on the sphere where x and y are degrees.
// @return The distance between the two points in RADIANS.  Multiply by radius to get arc
// length.
double spheredist_deg(const Point& p1, const Point& p2) {
    return spheredist_rad(Point(deg2rad(p1.x), deg2rad(p1.y)), Point(deg2rad(p2.x), deg2rad(p2.y)));
}

// Technically lat/long bounds, not really tied to earth radius.
bool isValidLngLat(double lng, double lat) {
    return abs(lng) <= 180 && abs(lat) <= 90;
}

double distance(const Point& p1, const Point& p2) {
    double a = p1.x - p2.x;
    double b = p1.y - p2.y;

    // Avoid numerical error if possible...
    if (a == 0)
        return abs(b);
    if (b == 0)
        return abs(a);

    return sqrt((a * a) + (b * b));
}

static inline Vector2_d toVector2(const Point& p) {
    return Vector2_d(p.x, p.y);
}

// Given a segment (A, B) and a segment (C, D), check whether they intersect.
bool linesIntersect(const Point& pA, const Point& pB, const Point& pC, const Point& pD) {
    Vector2_d a = toVector2(pA);
    Vector2_d b = toVector2(pB);
    Vector2_d c = toVector2(pC);
    Vector2_d d = toVector2(pD);

    // The normal of line AB
    Vector2_d normalAB = (b - a).Ortho();

    // Dot products of AC and the normal of AB
    // = 0 : C is on the line AB
    // > 0 : C is on one side
    // < 0 : C is on the other side
    double dotProdNormalAB_AC = normalAB.DotProd(c - a);
    double dotProdNormalAB_AD = normalAB.DotProd(d - a);

    // C and D can not on the same side of line AB
    if (dotProdNormalAB_AC * dotProdNormalAB_AD > 0)
        return false;

    // AB and CD are on the same line
    if (dotProdNormalAB_AC == 0 && dotProdNormalAB_AD == 0) {
        // Test if C or D is on segment AB.
        return (c - a).DotProd(c - b) <= 0 || (d - a).DotProd(d - b) <= 0;
    }

    // Check if A and B are on different sides of line CD.
    Vector2_d normalCD = (d - c).Ortho();
    double dotProdNormalCD_CA = normalCD.DotProd(a - c);
    double dotProdNormalCD_CB = normalCD.DotProd(b - c);
    return dotProdNormalCD_CA * dotProdNormalCD_CB <= 0;  // Perhaps A or B is on line CD
}

static bool circleContainsBoxInternal(const Circle& circle,
                                      const Box& box,
                                      bool includeCircleBoundary) {
    // NOTE: a circle of zero radius is a point, and there are NO points contained inside a
    // zero-radius circle, not even the point itself.

    const Point& a = box._min;
    const Point& b = box._max;
    double compareLL = distanceCompare(circle.center, a, circle.radius);  // Lower left
    double compareUR = distanceCompare(circle.center, b, circle.radius);  // Upper right
    // Upper Left
    double compareUL = distanceCompare(circle.center, Point(a.x, b.y), circle.radius);
    // Lower right
    double compareLR = distanceCompare(circle.center, Point(b.x, a.y), circle.radius);
    if (includeCircleBoundary) {
        return compareLL <= 0 && compareUR <= 0 && compareUL <= 0 && compareLR <= 0;
    } else {
        return compareLL < 0 && compareUR < 0 && compareUL < 0 && compareLR < 0;
    }
}

bool circleContainsBox(const Circle& circle, const Box& box) {
    return circleContainsBoxInternal(circle, box, true);
}

bool circleInteriorContainsBox(const Circle& circle, const Box& box) {
    return circleContainsBoxInternal(circle, box, false);
}

// Check the intersection by measuring the distance between circle center and box center.
static bool circleIntersectsWithBoxInternal(const Circle& circle,
                                            const Box& box,
                                            bool includeCircleBoundary) {
    // NOTE: a circle of zero radius is a point, and there are NO points to intersect inside a
    // zero-radius circle, not even the point itself.
    if (circle.radius == 0.0 && !includeCircleBoundary)
        return false;

    /* Collapses the four quadrants down into one.
     *   ________
     * r|___B___ \  <- a quarter round corner here. Let's name it "D".
     *  |       | |
     * h|       | |
     *  |   A   |C|
     *  |_______|_|
     *      w    r
     */

    Point boxCenter = box.center();
    double dx = abs(circle.center.x - boxCenter.x);
    double dy = abs(circle.center.y - boxCenter.y);
    double w = (box._max.x - box._min.x) / 2;
    double h = (box._max.y - box._min.y) / 2;
    const double& r = circle.radius;

    // Check if circle.center is in A, B or C.
    // The circle center could be above the box (B) or right to the box (C), but close enough.
    if (includeCircleBoundary) {
        if ((dx <= w + r && dy <= h) || (dx <= w && dy <= h + r))
            return true;
    } else {
        if ((dx < w + r && dy < h) || (dx < w && dy < h + r))
            return true;
    }

    // Now check if circle.center is in the round corner "D".
    double compareResult = distanceCompare(Point(dx, dy), Point(w, h), r);
    return compareResult < 0 || (compareResult == 0 && includeCircleBoundary);
}

bool circleIntersectsWithBox(const Circle& circle, const Box& box) {
    return circleIntersectsWithBoxInternal(circle, box, true);
}

bool circleInteriorIntersectsWithBox(const Circle& circle, const Box& box) {
    return circleIntersectsWithBoxInternal(circle, box, false);
}

bool lineIntersectsWithBox(const Point& a, const Point& b, const Box& box) {
    Point upperLeft(box._min.x, box._max.y);
    Point lowerRight(box._max.x, box._min.y);

    return linesIntersect(a, b, upperLeft, box._min) ||
        linesIntersect(a, b, box._min, lowerRight) || linesIntersect(a, b, lowerRight, box._max) ||
        linesIntersect(a, b, box._max, upperLeft);
}

// Doc: The last point specified is always implicitly connected to the first.
// [[ 0 , 0 ], [ 3 , 6 ], [ 6 , 0 ]]
bool edgesIntersectsWithBox(const vector<Point>& vertices, const Box& box) {
    for (size_t i = 0; i < vertices.size() - 1; i++) {
        if (lineIntersectsWithBox(vertices[i], vertices[i + 1], box))
            return true;
    }
    // The last point and first point.
    return lineIntersectsWithBox(vertices[vertices.size() - 1], vertices[0], box);
}

bool polygonContainsBox(const Polygon& polygon, const Box& box) {
    // All vertices of box have to be inside the polygon.
    if (!polygon.contains(box._min) || !polygon.contains(box._max) ||
        !polygon.contains(Point(box._min.x, box._max.y)) ||
        !polygon.contains(Point(box._max.x, box._min.y)))
        return false;

    // No intersection between the polygon edges and the box.
    return !edgesIntersectsWithBox(polygon.points(), box);
}

bool polygonIntersectsWithBox(const Polygon& polygon, const Box& box) {
    // 1. Polygon contains the box.
    // Check the relaxed condition that whether the polygon include any vertex of the box.
    if (polygon.contains(box._min) || polygon.contains(box._max) ||
        polygon.contains(Point(box._min.x, box._max.y)) ||
        polygon.contains(Point(box._max.x, box._min.y)))
        return true;

    // 2. Box contains polygon.
    // Check the relaxed condition that whether the box include any vertex of the polygon.
    for (vector<Point>::const_iterator it = polygon.points().begin(); it != polygon.points().end();
         it++) {
        if (box.inside(*it))
            return true;
    }

    // 3. Otherwise they intersect on a portion of both shapes.
    // Edges intersects
    return edgesIntersectsWithBox(polygon.points(), box);
}

bool ShapeProjection::supportsProject(const PointWithCRS& point, const CRS crs) {
    // Can always trivially project or project from SPHERE->FLAT
    if (point.crs == crs || point.crs == SPHERE)
        return true;

    invariant(point.crs == FLAT);
    // If crs is FLAT, we might be able to upgrade the point to SPHERE if it's a valid SPHERE
    // point (lng/lat in bounds).  In this case, we can use FLAT data with SPHERE predicates.
    return isValidLngLat(point.oldPoint.x, point.oldPoint.y);
}

bool ShapeProjection::supportsProject(const PolygonWithCRS& polygon, const CRS crs) {
    return polygon.crs == crs || (polygon.crs == STRICT_SPHERE && crs == SPHERE);
}

void ShapeProjection::projectInto(PointWithCRS* point, CRS crs) {
    dassert(supportsProject(*point, crs));

    if (point->crs == crs)
        return;

    if (FLAT == point->crs) {
        // Prohibit projection to STRICT_SPHERE CRS
        invariant(SPHERE == crs);

        // Note that it's (lat, lng) for S2 but (lng, lat) for MongoDB.
        S2LatLng latLng = S2LatLng::FromDegrees(point->oldPoint.y, point->oldPoint.x).Normalized();
        dassert(latLng.is_valid());
        point->point = latLng.ToPoint();
        point->cell = S2Cell(point->point);
        point->crs = SPHERE;
        return;
    }

    // Prohibit projection to STRICT_SPHERE CRS
    invariant(SPHERE == point->crs && FLAT == crs);
    // Just remove the additional spherical information
    point->point = S2Point();
    point->cell = S2Cell();
    point->crs = FLAT;
}

void ShapeProjection::projectInto(PolygonWithCRS* polygon, CRS crs) {
    if (polygon->crs == crs)
        return;

    // Only project from STRICT_SPHERE to SPHERE
    invariant(STRICT_SPHERE == polygon->crs && SPHERE == crs);
    polygon->crs = SPHERE;
}

}  // namespace mongo