// Copyright 2005 Google Inc. All Rights Reserved.

//

// 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.

//

// Author: ericv@google.com (Eric Veach)

#include "s2/s2cap.h"

#include <algorithm>

#include <cfloat>

#include <cmath>

#include <ostream>

#include <vector>

#include "absl/flags/flag.h"

#include "absl/log/absl_check.h"

#include "s2/util/coding/coder.h"

#include "s2/r1interval.h"

#include "s2/s1angle.h"

#include "s2/s1chord_angle.h"

#include "s2/s1interval.h"

#include "s2/s2cell.h"

#include "s2/s2cell_id.h"

#include "s2/s2debug.h"

#include "s2/s2edge_distances.h"

#include "s2/s2latlng.h"

#include "s2/s2latlng_rect.h"

#include "s2/s2metrics.h"

#include "s2/s2point.h"

#include "s2/s2pointutil.h"

#include "s2/s2region.h"

using std::fabs;

using std::max;

using std::vector;

double S2Cap::GetArea() const {

  return 2 * M_PI * max(0.0, height());

}

S2Point S2Cap::GetCentroid() const {

  // From symmetry, the centroid of the cap must be somewhere on the line

  // from the origin to the center of the cap on the surface of the sphere.

  // When a sphere is divided into slices of constant thickness by a set of

  // parallel planes, all slices have the same surface area. This implies

  // that the radial component of the centroid is simply the midpoint of the

  // range of radial distances spanned by the cap. That is easily computed

  // from the cap height.

  if (is_empty()) return S2Point();

  double r = 1.0 - 0.5 * height();

  return r * GetArea() * center_;

}

S2Cap S2Cap::Complement() const {

  // The complement of a full cap is an empty cap, not a singleton.

  // Also make sure that the complement of an empty cap is full.

  if (is_full()) return Empty();

  if (is_empty()) return Full();

  return S2Cap(-center_, S1ChordAngle::FromLength2(4 - radius_.length2()));

}

bool S2Cap::Contains(const S2Cap& other) const {

  if (is_full() || other.is_empty()) return true;

  return radius_ >= S1ChordAngle(center_, other.center_) + other.radius_;

}

bool S2Cap::Intersects(const S2Cap& other) const {

  if (is_empty() || other.is_empty()) return false;

  return radius_ + other.radius_ >= S1ChordAngle(center_, other.center_);

}

bool S2Cap::InteriorIntersects(const S2Cap& other) const {

  // Make sure this cap has an interior and the other cap is non-empty.

  if (radius_.length2() <= 0 || other.is_empty()) return false;

  return radius_ + other.radius_ > S1ChordAngle(center_, other.center_);

}

void S2Cap::AddPoint(const S2Point& p) {

  // Compute the squared chord length, then convert it into a height.

  ABSL_DCHECK(S2::IsUnitLength(p));

  if (is_empty()) {

    center_ = p;

    radius_ = S1ChordAngle::Zero();

  } else {

    // After calling cap.AddPoint(p), cap.Contains(p) must be true.  However

    // we don't need to do anything special to achieve this because Contains()

    // does exactly the same distance calculation that we do here.

    radius_ = max(radius_, S1ChordAngle(center_, p));

  }

}

void S2Cap::AddCap(const S2Cap& other) {

  if (is_empty()) {

    *this = other;

  } else if (!other.is_empty()) {

    // We round up the distance to ensure that the cap is actually contained.

    S1ChordAngle dist = S1ChordAngle(center_, other.center_) + other.radius_;

    dist = dist.PlusError(  //

        (2 * DBL_EPSILON + S1ChordAngle::kRelativeSumError) * dist.length2());

    radius_ = max(radius_, dist);

  }

}

S2Cap S2Cap::Expanded(S1Angle distance) const {

  ABSL_DCHECK_GE(distance.radians(), 0);

  if (is_empty()) return Empty();

  return S2Cap(center_, radius_ + S1ChordAngle(distance));

}

S2Cap S2Cap::Union(const S2Cap& other) const {

  if (radius_ < other.radius_) {

    return other.Union(*this);

  }

  if (is_full() || other.is_empty()) {

    return *this;

  }

  // This calculation would be more efficient using S1ChordAngles.

  S1Angle this_radius = GetRadius();

  S1Angle other_radius = other.GetRadius();

  S1Angle distance(center(), other.center());

  if (this_radius >= distance + other_radius) {

    return *this;

  } else {

    S1Angle result_radius = 0.5 * (distance + this_radius + other_radius);

    S2Point result_center =

        S2::GetPointOnLine(center(), other.center(),

                           0.5 * (distance - this_radius + other_radius));

    return S2Cap(result_center, result_radius);

  }

}

S2Region* S2Cap::Clone() const {

  return new S2Cap(*this);

}

S2Cap S2Cap::GetCapBound() const {

  return *this;

}

S2LatLngRect S2Cap::GetRectBound() const {

  if (is_empty()) return S2LatLngRect::Empty();

  // Convert the center to a (lat,lng) pair, and compute the cap angle.

  S2LatLng center_ll(center_);

  double cap_angle = GetRadius().radians();

  bool all_longitudes = false;

  double lat[2], lng[2];

  lng[0] = -M_PI;

  lng[1] = M_PI;

  // Check whether cap includes the south pole.

  lat[0] = center_ll.lat().radians() - cap_angle;

  if (lat[0] <= -M_PI_2) {

    lat[0] = -M_PI_2;

    all_longitudes = true;

  }

  // Check whether cap includes the north pole.

  lat[1] = center_ll.lat().radians() + cap_angle;

  if (lat[1] >= M_PI_2) {

    lat[1] = M_PI_2;

    all_longitudes = true;

  }

  if (!all_longitudes) {

    // 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 sin_a = sin(radius_);

    double sin_c = cos(center_ll.lat().radians());

    if (sin_a <= sin_c) {

      double angle_A = asin(sin_a / sin_c);

      lng[0] = remainder(center_ll.lng().radians() - angle_A, 2 * M_PI);

      lng[1] = remainder(center_ll.lng().radians() + angle_A, 2 * M_PI);

    }

  }

  return S2LatLngRect(R1Interval(lat[0], lat[1]),

                      S1Interval(lng[0], lng[1]));

}

// Computes a covering of the S2Cap.  In general the covering consists of at

// most 4 cells except for very large caps, which may need up to 6 cells.

// The output is not sorted.

void S2Cap::GetCellUnionBound(vector<S2CellId>* cell_ids) const {

  // TODO(ericv): The covering could be made quite a bit tighter by mapping

  // the cap to a rectangle in (i,j)-space and finding a covering for that.

  cell_ids->clear();

  // Find the maximum level such that the cap contains at most one cell vertex

  // and such that S2CellId::AppendVertexNeighbors() can be called.

  int level = S2::kMinWidth.GetLevelForMinValue(GetRadius().radians()) - 1;

  // If level < 0, then more than three face cells are required.

  if (level < 0) {

    cell_ids->reserve(6);

    for (int face = 0; face < 6; ++face) {

      cell_ids->push_back(S2CellId::FromFace(face));

    }

  } else {

    // The covering consists of the 4 cells at the given level that share the

    // cell vertex that is closest to the cap center.

    cell_ids->reserve(4);

    S2CellId(center_).AppendVertexNeighbors(level, cell_ids);

  }

}

bool S2Cap::Intersects(const S2Cell& cell, const S2Point* vertices) const {

  // 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 (radius_ >= S1ChordAngle::Right()) return false;

  // We need to check for empty caps due to the center check just below.

  if (is_empty()) return false;

  // Optimization: return true if the cell contains the cap center.  (This

  // allows half of the edge checks below to be skipped.)

  if (cell.Contains(center_)) return true;

  // At this point we know that the cell does not contain the cap center,

  // 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 sin2_angle = sin2(radius_);

  for (int k = 0; k < 4; ++k) {

    S2Point edge = cell.GetEdgeRaw(k);

    double dot = center_.DotProd(edge);

    if (dot > 0) {

      // The center 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 center is not contained with the cell).

      continue;

    }

    // The Norm2() factor is necessary because "edge" is not normalized.

    if (dot * dot > sin2_angle * 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.

    Vector3_d dir = edge.CrossProd(center_);

    if (dir.DotProd(vertices[k]) < 0 && dir.DotProd(vertices[(k+1)&3]) > 0)

      return true;

  }

  return false;

}

bool S2Cap::Contains(const S2Cell& cell) const {

  // 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[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);

}

bool S2Cap::MayIntersect(const S2Cell& cell) const {

  // If the cap contains any cell vertex, return true.

  S2Point vertices[4];

  for (int k = 0; k < 4; ++k) {

    vertices[k] = cell.GetVertex(k);

    if (Contains(vertices[k])) return true;

  }

  return Intersects(cell, vertices);

}

bool S2Cap::Contains(const S2Point& p) const {

  ABSL_DCHECK(S2::IsUnitLength(p));

  return S1ChordAngle(center_, p) <= radius_;

}

bool S2Cap::InteriorContains(const S2Point& p) const {

  ABSL_DCHECK(S2::IsUnitLength(p));

  return is_full() || S1ChordAngle(center_, p) < radius_;

}

bool S2Cap::operator==(const S2Cap& other) const {

  return (center_ == other.center_ && radius_ == other.radius_) ||

         (is_empty() && other.is_empty()) ||

         (is_full() && other.is_full());

}

bool S2Cap::ApproxEquals(const S2Cap& other, S1Angle max_error_angle) const {

  const double max_error = max_error_angle.radians();

  const double r2 = radius_.length2();

  const double other_r2 = other.radius_.length2();

  return (S2::ApproxEquals(center_, other.center_, max_error_angle) &&

          fabs(r2 - other_r2) <= max_error) ||

         (is_empty() && other_r2 <= max_error) ||

         (other.is_empty() && r2 <= max_error) ||

         (is_full() && other_r2 >= 2 - max_error) ||

         (other.is_full() && r2 >= 2 - max_error);

}

std::ostream& operator<<(std::ostream& os, const S2Cap& cap) {

  return os << "[Center=" << cap.center()

            << ", Radius=" << cap.GetRadius() << "]";

}

void S2Cap::Encode(Encoder* encoder) const {

  encoder->Ensure(4 * sizeof(double));

  encoder->putdouble(center_.x());

  encoder->putdouble(center_.y());

  encoder->putdouble(center_.z());

  encoder->putdouble(radius_.length2());

  ABSL_DCHECK_GE(encoder->avail(), 0);

}

bool S2Cap::Decode(Decoder* decoder) {

  if (decoder->avail() < 4 * sizeof(double)) return false;

  double x = decoder->getdouble();

  double y = decoder->getdouble();

  double z = decoder->getdouble();

  center_ = S2Point(x, y, z);

  radius_ = S1ChordAngle::FromLength2(decoder->getdouble());

  if (absl::GetFlag(FLAGS_s2debug)) {

    ABSL_CHECK(is_valid()) << "Invalid S2Cap: " << *this;

  }

  return true;

}