// 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/s1interval.h"

#include <algorithm>

#include <cfloat>

#include <cmath>

#include "absl/log/absl_check.h"

using std::fabs;

using std::max;

S1Interval S1Interval::FromPoint(double p) {

  if (p == -M_PI) p = M_PI;

  return S1Interval(p, p, ARGS_CHECKED);

}

double S1Interval::GetCenter() const {

  double center = 0.5 * (lo() + hi());

  if (!is_inverted()) return center;

  // Return the center in the range (-Pi, Pi].

  return (center <= 0) ? (center + M_PI) : (center - M_PI);

}

double S1Interval::GetLength() const {

  double length = hi() - lo();

  if (length >= 0) return length;

  length += 2 * M_PI;

  // Empty intervals have a negative length.

  return (length > 0) ? length : -1;

}

S1Interval S1Interval::Complement() const {

  if (lo() == hi()) return Full();   // Singleton.

  return S1Interval(hi(), lo(), ARGS_CHECKED);  // Handles empty and full.

}

double S1Interval::GetComplementCenter() const {

  if (lo() != hi()) {

    return Complement().GetCenter();

  } else {  // Singleton.

    return (hi() <= 0) ? (hi() + M_PI) : (hi() - M_PI);

  }

}

bool S1Interval::FastContains(double p) const {

  if (is_inverted()) {

    return (p >= lo() || p <= hi()) && !is_empty();

  } else {

    return p >= lo() && p <= hi();

  }

}

bool S1Interval::Contains(double p) const {

  // Works for empty, full, and singleton intervals.

  ABSL_DCHECK_LE(fabs(p), M_PI);

  if (p == -M_PI) p = M_PI;

  return FastContains(p);

}

bool S1Interval::InteriorContains(double p) const {

  // Works for empty, full, and singleton intervals.

  ABSL_DCHECK_LE(fabs(p), M_PI);

  if (p == -M_PI) p = M_PI;

  if (is_inverted()) {

    return p > lo() || p < hi();

  } else {

    return (p > lo() && p < hi()) || is_full();

  }

}

bool S1Interval::Contains(const S1Interval& y) const {

  // It might be helpful to compare the structure of these tests to

  // the simpler Contains(double) method above.

  if (is_inverted()) {

    if (y.is_inverted()) return y.lo() >= lo() && y.hi() <= hi();

    return (y.lo() >= lo() || y.hi() <= hi()) && !is_empty();

  } else {

    if (y.is_inverted()) return is_full() || y.is_empty();

    return y.lo() >= lo() && y.hi() <= hi();

  }

}

bool S1Interval::InteriorContains(const S1Interval& y) const {

  if (is_inverted()) {

    if (!y.is_inverted()) return y.lo() > lo() || y.hi() < hi();

    return (y.lo() > lo() && y.hi() < hi()) || y.is_empty();

  } else {

    if (y.is_inverted()) return is_full() || y.is_empty();

    return (y.lo() > lo() && y.hi() < hi()) || is_full();

  }

}

bool S1Interval::Intersects(const S1Interval& y) const {

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

  if (is_inverted()) {

    // Every non-empty inverted interval contains Pi.

    return y.is_inverted() || y.lo() <= hi() || y.hi() >= lo();

  } else {

    if (y.is_inverted()) return y.lo() <= hi() || y.hi() >= lo();

    return y.lo() <= hi() && y.hi() >= lo();

  }

}

bool S1Interval::InteriorIntersects(const S1Interval& y) const {

  if (is_empty() || y.is_empty() || lo() == hi()) return false;

  if (is_inverted()) {

    return y.is_inverted() || y.lo() < hi() || y.hi() > lo();

  } else {

    if (y.is_inverted()) return y.lo() < hi() || y.hi() > lo();

    return (y.lo() < hi() && y.hi() > lo()) || is_full();

  }

}

inline static double PositiveDistance(double a, double b) {

  // Compute the distance from "a" to "b" in the range [0, 2*Pi).

  // This is equivalent to (remainder(b - a - M_PI, 2 * M_PI) + M_PI),

  // except that it is more numerically stable (it does not lose

  // precision for very small positive distances).

  double d = b - a;

  if (d >= 0) return d;

  // We want to ensure that if b == Pi and a == (-Pi + eps),

  // the return result is approximately 2*Pi and not zero.

  return (b + M_PI) - (a - M_PI);

}

double S1Interval::GetDirectedHausdorffDistance(const S1Interval& y) const {

  if (y.Contains(*this)) return 0.0;  // this includes the case *this is empty

  if (y.is_empty()) return M_PI;  // maximum possible distance on S1

  double y_complement_center = y.GetComplementCenter();

  if (Contains(y_complement_center)) {

    return PositiveDistance(y.hi(), y_complement_center);

  } else {

    // The Hausdorff distance is realized by either two hi() endpoints or two

    // lo() endpoints, whichever is farther apart.

    double hi_hi = S1Interval(y.hi(), y_complement_center).Contains(hi()) ?

        PositiveDistance(y.hi(), hi()) : 0;

    double lo_lo = S1Interval(y_complement_center, y.lo()).Contains(lo()) ?

        PositiveDistance(lo(), y.lo()) : 0;

    ABSL_DCHECK(hi_hi > 0 || lo_lo > 0);

    return max(hi_hi, lo_lo);

  }

}

void S1Interval::AddPoint(double p) {

  ABSL_DCHECK_LE(fabs(p), M_PI);

  if (p == -M_PI) p = M_PI;

  if (FastContains(p)) return;

  if (is_empty()) {

    set_hi(p);

    set_lo(p);

  } else {

    // Compute distance from p to each endpoint.

    double dlo = PositiveDistance(p, lo());

    double dhi = PositiveDistance(hi(), p);

    if (dlo < dhi) {

      set_lo(p);

    } else {

      set_hi(p);

    }

    // Adding a point can never turn a non-full interval into a full one.

  }

}

double S1Interval::Project(double p) const {

  ABSL_DCHECK(!is_empty());

  ABSL_DCHECK_LE(fabs(p), M_PI);

  if (p == -M_PI) p = M_PI;

  if (FastContains(p)) return p;

  // Compute distance from p to each endpoint.

  double dlo = PositiveDistance(p, lo());

  double dhi = PositiveDistance(hi(), p);

  return (dlo < dhi) ? lo() : hi();

}

S1Interval S1Interval::FromPointPair(double p1, double p2) {

  ABSL_DCHECK_LE(fabs(p1), M_PI);

  ABSL_DCHECK_LE(fabs(p2), M_PI);

  if (p1 == -M_PI) p1 = M_PI;

  if (p2 == -M_PI) p2 = M_PI;

  if (PositiveDistance(p1, p2) <= M_PI) {

    return S1Interval(p1, p2, ARGS_CHECKED);

  } else {

    return S1Interval(p2, p1, ARGS_CHECKED);

  }

}

S1Interval S1Interval::Expanded(double margin) const {

  if (margin >= 0) {

    if (is_empty()) return *this;

    // Check whether this interval will be full after expansion, allowing

    // for a 1-bit rounding error when computing each endpoint.

    if (GetLength() + 2 * margin + 2 * DBL_EPSILON >= 2 * M_PI) return Full();

  } else {

    if (is_full()) return *this;

    // Check whether this interval will be empty after expansion, allowing

    // for a 1-bit rounding error when computing each endpoint.

    if (GetLength() + 2 * margin - 2 * DBL_EPSILON <= 0) return Empty();

  }

  S1Interval result(remainder(lo() - margin, 2*M_PI),

                    remainder(hi() + margin, 2*M_PI));

  if (result.lo() <= -M_PI) result.set_lo(M_PI);

  return result;

}

S1Interval S1Interval::Union(const S1Interval& y) const {

  // The y.is_full() case is handled correctly in all cases by the code

  // below, but can follow three separate code paths depending on whether

  // this interval is inverted, is non-inverted but contains Pi, or neither.

  if (y.is_empty()) return *this;

  if (FastContains(y.lo())) {

    if (FastContains(y.hi())) {

      // Either this interval contains y, or the union of the two

      // intervals is the Full() interval.

      if (Contains(y)) return *this;  // is_full() code path

      return Full();

    }

    return S1Interval(lo(), y.hi(), ARGS_CHECKED);

  }

  if (FastContains(y.hi())) return S1Interval(y.lo(), hi(), ARGS_CHECKED);

  // This interval contains neither endpoint of y.  This means that either y

  // contains all of this interval, or the two intervals are disjoint.

  if (is_empty() || y.FastContains(lo())) return y;

  // Check which pair of endpoints are closer together.

  double dlo = PositiveDistance(y.hi(), lo());

  double dhi = PositiveDistance(hi(), y.lo());

  if (dlo < dhi) {

    return S1Interval(y.lo(), hi(), ARGS_CHECKED);

  } else {

    return S1Interval(lo(), y.hi(), ARGS_CHECKED);

  }

}

S1Interval S1Interval::Intersection(const S1Interval& y) const {

  // The y.is_full() case is handled correctly in all cases by the code

  // below, but can follow three separate code paths depending on whether

  // this interval is inverted, is non-inverted but contains Pi, or neither.

  if (y.is_empty()) return Empty();

  if (FastContains(y.lo())) {

    if (FastContains(y.hi())) {

      // Either this interval contains y, or the region of intersection

      // consists of two disjoint subintervals.  In either case, we want

      // to return the shorter of the two original intervals.

      if (y.GetLength() < GetLength()) return y;  // is_full() code path

      return *this;

    }

    return S1Interval(y.lo(), hi(), ARGS_CHECKED);

  }

  if (FastContains(y.hi())) return S1Interval(lo(), y.hi(), ARGS_CHECKED);

  // This interval contains neither endpoint of y.  This means that either y

  // contains all of this interval, or the two intervals are disjoint.

  if (y.FastContains(lo())) return *this;  // is_empty() okay here

  ABSL_DCHECK(!Intersects(y));

  return Empty();

}

bool S1Interval::ApproxEquals(const S1Interval& y, double max_error) const {

  // Full and empty intervals require special cases because the "endpoints"

  // are considered to be positioned arbitrarily.

  if (is_empty()) return y.GetLength() <= 2 * max_error;

  if (y.is_empty()) return GetLength() <= 2 * max_error;

  if (is_full()) return y.GetLength() >= 2 * (M_PI - max_error);

  if (y.is_full()) return GetLength() >= 2 * (M_PI - max_error);

  // The purpose of the last test below is to verify that moving the endpoints

  // does not invert the interval, e.g. [-1e20, 1e20] vs. [1e20, -1e20].

  return (fabs(remainder(y.lo() - lo(), 2 * M_PI)) <= max_error &&

          fabs(remainder(y.hi() - hi(), 2 * M_PI)) <= max_error &&

          fabs(GetLength() - y.GetLength()) <= 2 * max_error);

}