mongo/src/third_party/s2/s2cell_test.cc

// Copyright 2005 Google Inc. All Rights Reserved.

#include "s2cell.h"

#include <cstdio>
#include <map>
using std::map;
using std::multimap;

#include <vector>
using std::vector;


#include "base/commandlineflags.h"
#include "base/logging.h"
#include "testing/base/public/gunit.h"
#include "s2.h"
#include "s2cap.h"
#include "s2latlngrect.h"
#include "s2testing.h"

TEST(S2Cell, TestFaces) {
  map<S2Point, int> edge_counts;
  map<S2Point, int> vertex_counts;
  for (int face = 0; face < 6; ++face) {
    S2CellId id = S2CellId::FromFacePosLevel(face, 0, 0);
    S2Cell cell(id);
    EXPECT_EQ(id, cell.id());
    EXPECT_EQ(face, cell.face());
    EXPECT_EQ(0, cell.level());
    // Top-level faces have alternating orientations to get RHS coordinates.
    EXPECT_EQ(face & S2::kSwapMask, cell.orientation());
    EXPECT_FALSE(cell.is_leaf());
    for (int k = 0; k < 4; ++k) {
      edge_counts[cell.GetEdgeRaw(k)] += 1;
      vertex_counts[cell.GetVertexRaw(k)] += 1;
      EXPECT_DOUBLE_EQ(0.0, cell.GetVertexRaw(k).DotProd(cell.GetEdgeRaw(k)));
      EXPECT_DOUBLE_EQ(0.0,
                       cell.GetVertexRaw((k+1)&3).DotProd(cell.GetEdgeRaw(k)));
      EXPECT_DOUBLE_EQ(1.0,
                       cell.GetVertexRaw(k)
                       .CrossProd(cell.GetVertexRaw((k+1)&3))
                       .Normalize().DotProd(cell.GetEdge(k)));
    }
  }
  // Check that edges have multiplicity 2 and vertices have multiplicity 3.
  for (map<S2Point, int>::iterator i = edge_counts.begin();
       i != edge_counts.end(); ++i) {
    EXPECT_EQ(2, i->second);
  }
  for (map<S2Point, int>::iterator i = vertex_counts.begin();
       i != vertex_counts.end(); ++i) {
    EXPECT_EQ(3, i->second);
  }
}

struct LevelStats {
  double count;
  double min_area, max_area, avg_area;
  double min_width, max_width, avg_width;
  double min_edge, max_edge, avg_edge, max_edge_aspect;
  double min_diag, max_diag, avg_diag, max_diag_aspect;
  double min_angle_span, max_angle_span, avg_angle_span;
  double min_approx_ratio, max_approx_ratio;
  LevelStats()
    : count(0), min_area(100), max_area(0), avg_area(0),
      min_width(100), max_width(0), avg_width(0),
      min_edge(100), max_edge(0), avg_edge(0), max_edge_aspect(0),
      min_diag(100), max_diag(0), avg_diag(0), max_diag_aspect(0),
      min_angle_span(100), max_angle_span(0), avg_angle_span(0),
      min_approx_ratio(100), max_approx_ratio(0) {}
};
static vector<LevelStats> level_stats(S2CellId::kMaxLevel+1);

static void GatherStats(S2Cell const& cell) {
  LevelStats* s = &level_stats[cell.level()];
  double exact_area = cell.ExactArea();
  double approx_area = cell.ApproxArea();
  double min_edge = 100, max_edge = 0, avg_edge = 0;
  double min_diag = 100, max_diag = 0;
  double min_width = 100, max_width = 0;
  double min_angle_span = 100, max_angle_span = 0;
  for (int i = 0; i < 4; ++i) {
    double edge = cell.GetVertexRaw(i).Angle(cell.GetVertexRaw((i+1)&3));
    min_edge = min(edge, min_edge);
    max_edge = max(edge, max_edge);
    avg_edge += 0.25 * edge;
    S2Point mid = cell.GetVertexRaw(i) + cell.GetVertexRaw((i+1)&3);
    double width = M_PI_2 - mid.Angle(cell.GetEdgeRaw(i^2));
    min_width = min(width, min_width);
    max_width = max(width, max_width);
    if (i < 2) {
      double diag = cell.GetVertexRaw(i).Angle(cell.GetVertexRaw(i^2));
      min_diag = min(diag, min_diag);
      max_diag = max(diag, max_diag);
      double angle_span = cell.GetEdgeRaw(i).Angle(-cell.GetEdgeRaw(i^2));
      min_angle_span = min(angle_span, min_angle_span);
      max_angle_span = max(angle_span, max_angle_span);
    }
  }
  s->count += 1;
  s->min_area = min(exact_area, s->min_area);
  s->max_area = max(exact_area, s->max_area);
  s->avg_area += exact_area;
  s->min_width = min(min_width, s->min_width);
  s->max_width = max(max_width, s->max_width);
  s->avg_width += 0.5 * (min_width + max_width);
  s->min_edge = min(min_edge, s->min_edge);
  s->max_edge = max(max_edge, s->max_edge);
  s->avg_edge += avg_edge;
  s->max_edge_aspect = max(max_edge / min_edge, s->max_edge_aspect);
  s->min_diag = min(min_diag, s->min_diag);
  s->max_diag = max(max_diag, s->max_diag);
  s->avg_diag += 0.5 * (min_diag + max_diag);
  s->max_diag_aspect = max(max_diag / min_diag, s->max_diag_aspect);
  s->min_angle_span = min(min_angle_span, s->min_angle_span);
  s->max_angle_span = max(max_angle_span, s->max_angle_span);
  s->avg_angle_span += 0.5 * (min_angle_span + max_angle_span);
  double approx_ratio = approx_area / exact_area;
  s->min_approx_ratio = min(approx_ratio, s->min_approx_ratio);
  s->max_approx_ratio = max(approx_ratio, s->max_approx_ratio);
}

static void TestSubdivide(S2Cell const& cell) {
  GatherStats(cell);
  if (cell.is_leaf()) return;

  S2Cell children[4];
  CHECK(cell.Subdivide(children));
  S2CellId child_id = cell.id().child_begin();
  double exact_area = 0;
  double approx_area = 0;
  double average_area = 0;
  for (int i = 0; i < 4; ++i, child_id = child_id.next()) {
    exact_area += children[i].ExactArea();
    approx_area += children[i].ApproxArea();
    average_area += children[i].AverageArea();

    // Check that the child geometry is consistent with its cell ID.
    EXPECT_EQ(child_id, children[i].id());
    EXPECT_TRUE(S2::ApproxEquals(children[i].GetCenter(), child_id.ToPoint()));
    S2Cell direct(child_id);
    EXPECT_EQ(direct.face(), children[i].face());
    EXPECT_EQ(direct.level(), children[i].level());
    EXPECT_EQ(direct.orientation(), children[i].orientation());
    EXPECT_EQ(direct.GetCenterRaw(), children[i].GetCenterRaw());
    for (int k = 0; k < 4; ++k) {
      EXPECT_EQ(direct.GetVertexRaw(k), children[i].GetVertexRaw(k));
      EXPECT_EQ(direct.GetEdgeRaw(k), children[i].GetEdgeRaw(k));
    }

    // Test Contains() and MayIntersect().
    EXPECT_TRUE(cell.Contains(children[i]));
    EXPECT_TRUE(cell.MayIntersect(children[i]));
    EXPECT_FALSE(children[i].Contains(cell));
    EXPECT_TRUE(cell.Contains(children[i].GetCenterRaw()));
    EXPECT_TRUE(cell.VirtualContainsPoint(children[i].GetCenterRaw()));
    for (int j = 0; j < 4; ++j) {
      EXPECT_TRUE(cell.Contains(children[i].GetVertexRaw(j)));
      if (j != i) {
        EXPECT_FALSE(children[i].Contains(children[j].GetCenterRaw()));
        EXPECT_FALSE(
            children[i].VirtualContainsPoint(children[j].GetCenterRaw()));
        EXPECT_FALSE(children[i].MayIntersect(children[j]));
      }
    }

    // Test GetCapBound and GetRectBound.
    S2Cap parent_cap = cell.GetCapBound();
    S2LatLngRect parent_rect = cell.GetRectBound();
    if (cell.Contains(S2Point(0, 0, 1)) || cell.Contains(S2Point(0, 0, -1))) {
      EXPECT_TRUE(parent_rect.lng().is_full());
    }
    S2Cap child_cap = children[i].GetCapBound();
    S2LatLngRect child_rect = children[i].GetRectBound();
    EXPECT_TRUE(child_cap.Contains(children[i].GetCenter()));
    EXPECT_TRUE(child_rect.Contains(children[i].GetCenterRaw()));
    EXPECT_TRUE(parent_cap.Contains(children[i].GetCenter()));
    EXPECT_TRUE(parent_rect.Contains(children[i].GetCenterRaw()));
    for (int j = 0; j < 4; ++j) {
      EXPECT_TRUE(child_cap.Contains(children[i].GetVertex(j)));
      EXPECT_TRUE(child_rect.Contains(children[i].GetVertex(j)));
      EXPECT_TRUE(child_rect.Contains(children[i].GetVertexRaw(j)));
      EXPECT_TRUE(parent_cap.Contains(children[i].GetVertex(j)));
      EXPECT_TRUE(parent_rect.Contains(children[i].GetVertex(j)));
      EXPECT_TRUE(parent_rect.Contains(children[i].GetVertexRaw(j)));
      if (j != i) {
        // The bounding caps and rectangles should be tight enough so that
        // they exclude at least two vertices of each adjacent cell.
        int cap_count = 0;
        int rect_count = 0;
        for (int k = 0; k < 4; ++k) {
          if (child_cap.Contains(children[j].GetVertex(k)))
            ++cap_count;
          if (child_rect.Contains(children[j].GetVertexRaw(k)))
            ++rect_count;
        }
        EXPECT_LE(cap_count, 2);
        if (child_rect.lat_lo().radians() > -M_PI_2 &&
            child_rect.lat_hi().radians() < M_PI_2) {
          // Bounding rectangles may be too large at the poles because the
          // pole itself has an arbitrary fixed longitude.
          EXPECT_LE(rect_count, 2);
        }
      }
    }

    // Check all children for the first few levels, and then sample randomly.
    // Also subdivide one corner cell, one edge cell, and one center cell
    // so that we have a better chance of sample the minimum metric values.
    bool force_subdivide = false;
    S2Point center = S2::GetNorm(children[i].face());
    S2Point edge = center + S2::GetUAxis(children[i].face());
    S2Point corner = edge + S2::GetVAxis(children[i].face());
    for (int j = 0; j < 4; ++j) {
      S2Point p = children[i].GetVertexRaw(j);
      if (p == center || p == edge || p == corner)
        force_subdivide = true;
    }
    if (force_subdivide || cell.level() < (DEBUG_MODE ? 5 : 6) ||
        S2Testing::rnd.OneIn(DEBUG_MODE ? 5 : 4)) {
      TestSubdivide(children[i]);
    }
  }

  // Check sum of child areas equals parent area.
  //
  // For ExactArea(), the best relative error we can expect is about 1e-6
  // because the precision of the unit vector coordinates is only about 1e-15
  // and the edge length of a leaf cell is about 1e-9.
  //
  // For ApproxArea(), the areas are accurate to within a few percent.
  //
  // For AverageArea(), the areas themselves are not very accurate, but
  // the average area of a parent is exactly 4 times the area of a child.

  EXPECT_LE(fabs(log(exact_area / cell.ExactArea())), fabs(log(1 + 1e-6)));
  EXPECT_LE(fabs(log(approx_area / cell.ApproxArea())), fabs(log(1.03)));
  EXPECT_LE(fabs(log(average_area / cell.AverageArea())), fabs(log(1 + 1e-15)));
}

template <int dim>
static void CheckMinMaxAvg(
    char const* label, int level, double count, double abs_error,
    double min_value, double max_value, double avg_value,
    S2::Metric<dim> const& min_metric,
    S2::Metric<dim> const& max_metric,
    S2::Metric<dim> const& avg_metric) {

  // All metrics are minimums, maximums, or averages of differential
  // quantities, and therefore will not be exact for cells at any finite
  // level.  The differential minimum is always a lower bound, and the maximum
  // is always an upper bound, but these minimums and maximums may not be
  // achieved for two different reasons.  First, the cells at each level are
  // sampled and we may miss the most extreme examples.  Second, the actual
  // metric for a cell is obtained by integrating the differential quantity,
  // which is not constant across the cell.  Therefore cells at low levels
  // (bigger cells) have smaller variations.
  //
  // The "tolerance" below is an attempt to model both of these effects.
  // At low levels, error is dominated by the variation of differential
  // quantities across the cells, while at high levels error is dominated by
  // the effects of random sampling.
  double tolerance = (max_metric.GetValue(level) - min_metric.GetValue(level)) /
                     sqrt(min(count, 0.5 * double(1 << level)));
  if (tolerance == 0) tolerance = abs_error;

  double min_error = min_value - min_metric.GetValue(level);
  double max_error = max_metric.GetValue(level) - max_value;
  double avg_error = fabs(avg_metric.GetValue(level) - avg_value);
  printf("%-10s (%6.0f samples, tolerance %8.3g) - min (%9.3g : %9.3g) "
         "max (%9.3g : %9.3g), avg (%9.3g : %9.3g)\n",
         label, count, tolerance,
         min_error / min_value, min_error / tolerance,
         max_error / max_value, max_error / tolerance,
         avg_error / avg_value, avg_error / tolerance);

  EXPECT_LE(min_metric.GetValue(level), min_value + abs_error);
  EXPECT_GE(min_metric.GetValue(level), min_value - tolerance);
  EXPECT_LE(max_metric.GetValue(level), max_value + tolerance);
  EXPECT_GE(max_metric.GetValue(level), max_value - abs_error);
  EXPECT_NEAR(avg_metric.GetValue(level), avg_value, 10 * tolerance);
}

TEST(S2Cell, TestSubdivide) {
  for (int face = 0; face < 6; ++face) {
    TestSubdivide(S2Cell::FromFacePosLevel(face, 0, 0));
  }

  // The maximum edge *ratio* is the ratio of the longest edge of any cell to
  // the shortest edge of any cell at the same level (and similarly for the
  // maximum diagonal ratio).
  //
  // The maximum edge *aspect* is the maximum ratio of the longest edge of a
  // cell to the shortest edge of that same cell (and similarly for the
  // maximum diagonal aspect).

  printf("Level    Area      Edge          Diag          Approx       Average\n");
  printf("        Ratio  Ratio Aspect  Ratio Aspect    Min    Max    Min    Max\n");
  for (int i = 0; i <= S2CellId::kMaxLevel; ++i) {
    LevelStats* s = &level_stats[i];
    if (s->count > 0) {
      s->avg_area /= s->count;
      s->avg_width /= s->count;
      s->avg_edge /= s->count;
      s->avg_diag /= s->count;
      s->avg_angle_span /= s->count;
    }
    printf("%5d  %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\n",
           i, s->max_area / s->min_area,
           s->max_edge / s->min_edge, s->max_edge_aspect,
           s->max_diag / s->min_diag, s->max_diag_aspect,
           s->min_approx_ratio, s->max_approx_ratio,
           S2Cell::AverageArea(i) / s->max_area,
           S2Cell::AverageArea(i) / s->min_area);
  }

  // Now check the validity of the S2 length and area metrics.
  for (int i = 0; i <= S2CellId::kMaxLevel; ++i) {
    LevelStats const* s = &level_stats[i];
    if (s->count == 0) continue;

    printf("Level %2d - metric (error/actual : error/tolerance)\n", i);

    // The various length calculations are only accurate to 1e-15 or so,
    // so we need to allow for this amount of discrepancy with the theoretical
    // minimums and maximums.  The area calculation is accurate to about 1e-15
    // times the cell width.
    CheckMinMaxAvg("area", i, s->count, 1e-15 * s->min_width,
                   s->min_area, s->max_area, s->avg_area,
                   S2::kMinArea, S2::kMaxArea, S2::kAvgArea);
    CheckMinMaxAvg("width", i, s->count, 1e-15,
                   s->min_width, s->max_width, s->avg_width,
                   S2::kMinWidth, S2::kMaxWidth, S2::kAvgWidth);
    CheckMinMaxAvg("edge", i, s->count, 1e-15,
                   s->min_edge, s->max_edge, s->avg_edge,
                   S2::kMinEdge, S2::kMaxEdge, S2::kAvgEdge);
    CheckMinMaxAvg("diagonal", i, s->count, 1e-15,
                   s->min_diag, s->max_diag, s->avg_diag,
                   S2::kMinDiag, S2::kMaxDiag, S2::kAvgDiag);
    CheckMinMaxAvg("angle span", i, s->count, 1e-15,
                   s->min_angle_span, s->max_angle_span, s->avg_angle_span,
                   S2::kMinAngleSpan, S2::kMaxAngleSpan, S2::kAvgAngleSpan);

    // The aspect ratio calculations are ratios of lengths and are therefore
    // less accurate at higher subdivision levels.
    EXPECT_LE(s->max_edge_aspect, S2::kMaxEdgeAspect + 1e-15 * (1 << i));
    EXPECT_LE(s->max_diag_aspect, S2::kMaxDiagAspect + 1e-15 * (1 << i));
  }
}

static int const kMaxLevel = DEBUG_MODE ? 6 : 11;

static void ExpandChildren1(S2Cell const& cell) {
  S2Cell children[4];
  CHECK(cell.Subdivide(children));
  if (children[0].level() < kMaxLevel) {
    for (int pos = 0; pos < 4; ++pos) {
      ExpandChildren1(children[pos]);
    }
  }
}

static void ExpandChildren2(S2Cell const& cell) {
  S2CellId id = cell.id().child_begin();
  for (int pos = 0; pos < 4; ++pos, id = id.next()) {
    S2Cell child(id);
    if (child.level() < kMaxLevel) ExpandChildren2(child);
  }
}

TEST(S2Cell, TestPerformance) {
  double subdivide_start = S2Testing::GetCpuTime();
  ExpandChildren1(S2Cell::FromFacePosLevel(0, 0, 0));
  double subdivide_time = S2Testing::GetCpuTime() - subdivide_start;
  fprintf(stderr, "Subdivide: %.3f seconds\n", subdivide_time);

  double constructor_start = S2Testing::GetCpuTime();
  ExpandChildren2(S2Cell::FromFacePosLevel(0, 0, 0));
  double constructor_time = S2Testing::GetCpuTime() - constructor_start;
  fprintf(stderr, "Constructor: %.3f seconds\n", constructor_time);
}