search

google / s2geometry

Computational geometry and spatial indexing on the sphere
http://s2geometry.io/
Apache License 2.0
2.29k stars 302 forks source link

S2CellUnion.Normalize Test Failure Case #260

Closed mlptownsend closed 2 years ago

mlptownsend commented 2 years ago

I managed to come up with a test case for S2CellUnion.Normalize that fails. At least it fails in v.0.9.0, as I couldn't get the latest to build properly in either Windows or Ubuntu.

I found this by repeatedly running the unmodified Normalize over and over until it broke down. Then I extracted the input and expected values and placed them in the existing logic.

TEST(S2CellUnion, Normalize) {
  // Try a bunch of random test cases, and keep track of average
  // statistics for normalization (to see if they agree with the
  // analysis above).
  double in_sum = 0, out_sum = 0;
  static const int kIters = 1;
  for (int i = 0; i < kIters; ++i) {
    vector<S2CellId> input, expected;
      input.push_back(S2CellId(1152921504606846976)); input.push_back(S2CellId(3458764513820540928)); input.push_back(S2CellId(4683743612465315840)); input.push_back(S2CellId(4629700416936869888)); input.push_back(S2CellId(4665729213955833856)); input.push_back(S2CellId(6052837899185946624)); input.push_back(S2CellId(8070450532247928832)); input.push_back(S2CellId(9295429630892703744)); input.push_back(S2CellId(9439544818968559616)); input.push_back(S2CellId(9493588014497005568)); input.push_back(S2CellId(9529616811515969536)); input.push_back(S2CellId(9525113211888599040)); input.push_back(S2CellId(9528490911609126912)); input.push_back(S2CellId(9565645608534933504)); input.push_back(S2CellId(9601674405553897472)); input.push_back(S2CellId(9624192403690749952)); input.push_back(S2CellId(9633199602945490944)); input.push_back(S2CellId(9642206802200231936)); input.push_back(S2CellId(9651214001454972928)); input.push_back(S2CellId(9637703202572861440)); input.push_back(S2CellId(9673731999591825408)); input.push_back(S2CellId(9709760796610789376)); input.push_back(S2CellId(9701879497262891008)); input.push_back(S2CellId(9701035072332759040)); input.push_back(S2CellId(9704131297076576256)); input.push_back(S2CellId(9706383096890261504)); input.push_back(S2CellId(9708634896703946752)); input.push_back(S2CellId(9745789593629753344)); input.push_back(S2CellId(9781818390648717312)); input.push_back(S2CellId(10088063165309911040)); input.push_back(S2CellId(10664523917613334528)); input.push_back(S2CellId(11240984669916758016)); input.push_back(S2CellId(11024811887802974208)); input.push_back(S2CellId(11240984669916758016)); input.push_back(S2CellId(12682136550675316736));
      expected.push_back(S2CellId(1152921504606846976)); expected.push_back(S2CellId(3458764513820540928)); expected.push_back(S2CellId(4683743612465315840)); expected.push_back(S2CellId(6052837899185946624)); expected.push_back(S2CellId(8070450532247928832)); expected.push_back(S2CellId(10376293541461622784)); expected.push_back(S2CellId(12682136550675316736));
    //AddCells(S2CellId::None(), false, &input, &expected);
    in_sum += input.size();
    out_sum += expected.size();
    S2CellUnion cellunion(input);
    EXPECT_EQ(expected.size(), cellunion.size());
    for (int i = 0; i < expected.size(); ++i) {
      EXPECT_EQ(expected[i], cellunion[i]);
    }

    // Test GetCapBound().
    S2Cap cap = cellunion.GetCapBound();
    for (S2CellId id : cellunion) {
//Fails here
      EXPECT_TRUE(cap.Contains(S2Cell(id)));
    }

    // Test Contains(S2CellId) and Intersects(S2CellId).
    for (S2CellId input_id : input) {
      EXPECT_TRUE(cellunion.Contains(input_id));
      EXPECT_TRUE(cellunion.Contains(input_id.ToPoint()));
      EXPECT_TRUE(cellunion.Intersects(input_id));
      if (!input_id.is_face()) {
        EXPECT_TRUE(cellunion.Intersects(input_id.parent()));
        if (input_id.level() > 1) {
          EXPECT_TRUE(cellunion.Intersects(input_id.parent().parent()));
          EXPECT_TRUE(cellunion.Intersects(input_id.parent(0)));
        }
      }
      if (!input_id.is_leaf()) {
        EXPECT_TRUE(cellunion.Contains(input_id.child_begin()));
        EXPECT_TRUE(cellunion.Intersects(input_id.child_begin()));
        EXPECT_TRUE(cellunion.Contains(input_id.child_end().prev()));
        EXPECT_TRUE(cellunion.Intersects(input_id.child_end().prev()));
        EXPECT_TRUE(cellunion.Contains(
                        input_id.child_begin(S2CellId::kMaxLevel)));
        EXPECT_TRUE(cellunion.Intersects(
                        input_id.child_begin(S2CellId::kMaxLevel)));
      }
    }
    for (S2CellId expected_id : expected) {
      if (!expected_id.is_face()) {
        EXPECT_TRUE(!cellunion.Contains(expected_id.parent()));
        EXPECT_TRUE(!cellunion.Contains(expected_id.parent(0)));
      }
    }

    // Test Contains(S2CellUnion), Intersects(S2CellUnion), Union(),
    // Intersection(), and Difference().
    vector<S2CellId> x, y, x_or_y, x_and_y;
    for (S2CellId input_id : input) {
      bool in_x = rnd.OneIn(2);
      bool in_y = rnd.OneIn(2);
      if (in_x) x.push_back(input_id);
      if (in_y) y.push_back(input_id);
      if (in_x || in_y) x_or_y.push_back(input_id);
    }
    S2CellUnion xcells(std::move(x));
    S2CellUnion ycells(std::move(y));
    S2CellUnion x_or_y_expected(std::move(x_or_y));
    S2CellUnion x_or_y_cells = xcells.Union(ycells);
    EXPECT_TRUE(x_or_y_cells == x_or_y_expected);

    // Compute the intersection of "x" with each cell of "y",
    // check that this intersection is correct, and append the
    // results to x_and_y_expected.
    for (S2CellId yid : ycells) {
      S2CellUnion ucells = xcells.Intersection(yid);
      for (S2CellId xid : xcells) {
        if (xid.contains(yid)) {
          EXPECT_TRUE(ucells.size() == 1 && ucells[0] == yid);
        } else if (yid.contains(xid)) {
          EXPECT_TRUE(ucells.Contains(xid));
        }
      }
      for (S2CellId uid : ucells) {
        EXPECT_TRUE(xcells.Contains(uid));
        EXPECT_TRUE(yid.contains(uid));
      }
      x_and_y.insert(x_and_y.end(), ucells.begin(), ucells.end());
    }
    S2CellUnion x_and_y_expected(std::move(x_and_y));
    S2CellUnion x_and_y_cells = xcells.Intersection(ycells);
    EXPECT_TRUE(x_and_y_cells == x_and_y_expected);

    S2CellUnion x_minus_y_cells = xcells.Difference(ycells);
    S2CellUnion y_minus_x_cells = ycells.Difference(xcells);
    EXPECT_TRUE(xcells.Contains(x_minus_y_cells));
    EXPECT_TRUE(!x_minus_y_cells.Intersects(ycells));
    EXPECT_TRUE(ycells.Contains(y_minus_x_cells));
    EXPECT_TRUE(!y_minus_x_cells.Intersects(xcells));
    EXPECT_TRUE(!x_minus_y_cells.Intersects(y_minus_x_cells));

    S2CellUnion diff_intersection_union =
        x_minus_y_cells.Union(y_minus_x_cells).Union(x_and_y_cells);
    EXPECT_TRUE(diff_intersection_union == x_or_y_cells);

    vector<S2CellId> test, dummy;
    AddCells(S2CellId::None(), false, &test, &dummy);
    for (S2CellId test_id : test) {
      bool contains = false, intersects = false;
      for (S2CellId expected_id : expected) {
        if (expected_id.contains(test_id)) contains = true;
        if (expected_id.intersects(test_id)) intersects = true;
      }
      EXPECT_EQ(contains, cellunion.Contains(test_id));
      EXPECT_EQ(intersects, cellunion.Intersects(test_id));
    }
  }
  printf("avg in %.2f, avg out %.2f\n", in_sum / kIters, out_sum / kIters);
}
Failure
Value of: cap.Contains(S2Cell(id))
  Actual: false
Expected: true
avg in 35.00, avg out 7.00

image

jmr commented 2 years ago

I think this is expected.

https://github.com/google/s2geometry/blob/7a40135059545396237a0199c558d749fe3be0b1/src/s2/s2region.h#L81-L84

  // Returns true if the region completely contains the given cell, otherwise
  // either the region does not contain the cell or the containment relationship
  // could not be determined.
  virtual bool Contains(const S2Cell& cell) const = 0;

The cap will be tight to the corner of some cell, and containment won't be determinable due to numerical instability.

What's the S2ChordAngle of each vertex vs the cap center vs the chord angle of the cap?

mlptownsend commented 2 years ago

Yeah, you're right. It's throwing it out on the first vertex here. image

image 3.9874261770344996 image 3.9874261770344992

It's just a floating point comparison doing what floating point comparisons do then :)

// Test GetCapBound().
    S2Cap cap = cellunion.GetCapBound();
    S2Point center = cap.center();
    S2Point centerN = center.Normalize();
    S1ChordAngle centerAngle = cap.radius();
    char buf[200];

      sprintf(buf, "Cap: %.17f / %.17f", i, centerAngle.radians(), centerAngle.length2());
      puts(buf);
      memset(buf, 0, 200);

      for(int i = 0; i < input.size(); ++i) {
          S2CellId c = input[i];
          S2Point pc = c.ToPoint().Normalize();
          S1ChordAngle vsCenter = S1ChordAngle(centerN, pc);
          sprintf(buf, "V%i to Center: %.17f / %.17f", i, vsCenter.radians(), vsCenter.length2());
          puts(buf);
          memset(buf, 0, 200);
      }

      int j = 0;
    for (S2CellId id : cellunion) {
        S2Point pc = id.ToPoint().Normalize();
        S1ChordAngle vsCenter = S1ChordAngle(centerN, pc);
        sprintf(buf, "CU-V%i to Center: %.17f / %.17f", j++, vsCenter.radians(), vsCenter.length2());
        puts(buf);
        memset(buf, 0, 200);
      EXPECT_TRUE(cap.Contains(S2Cell(id)));
    }
Cap: 3.02940076358108490 / 3.98742617703449920
V0 to Center: 1.65004457710871555 / 2.15833065217771214
V1 to Center: 1.65004457710871555 / 2.15833065217771214
V2 to Center: 2.48244445722154428 / 3.58102839999664058
V3 to Center: 2.38355502666034313 / 3.45237307025370832
V4 to Center: 2.46511237321986565 / 3.55956214577003571
V5 to Center: 2.49691898001032797 / 3.59859182265260014
V6 to Center: 1.49154807648107801 / 1.84166934782228831
V7 to Center: 1.97582029277483429 / 2.78808157489446717
V8 to Center: 1.61375148565650628 / 2.08588390065200446
V9 to Center: 1.71782752487735846 / 2.29300402551243909
V10 to Center: 1.74908418543264177 / 2.35468966472460828
V11 to Center: 1.80730572343033846 / 2.46862125917851971
V12 to Center: 1.77805633212831760 / 2.41155864239179873
V13 to Center: 1.55677094755307266 / 1.97195016115732824
V14 to Center: 1.57259255779997487 / 2.00359246007834280
V15 to Center: 1.72370218603139924 / 2.30462145414183794
V16 to Center: 1.71526038248314472 / 2.28792417784496394
V17 to Center: 1.81935529331584345 / 2.49201493524619533
V18 to Center: 1.82846267105155436 / 2.50964926478571382
V19 to Center: 1.77162414590865458 / 2.39896116095157685
V20 to Center: 1.98229513488767850 / 2.79996693571606503
V21 to Center: 1.95436458746768915 / 2.74846363702268803
V22 to Center: 1.87317742760410821 / 2.59558823850561549
V23 to Center: 1.86249494934561266 / 2.57515901289319471
V24 to Center: 1.86395037577228018 / 2.57794629176132251
V25 to Center: 1.91445497573352763 / 2.67386812164724130
V26 to Center: 1.92430074130738826 / 2.69237529307332757
V27 to Center: 2.14823728819093684 / 3.09176320900185253
V28 to Center: 2.17870527130146785 / 3.14230455884868665
V29 to Center: 1.89329355565500546 / 2.63387207078301078
V30 to Center: 1.16308599361227838 / 1.20698323233063198
V31 to Center: 1.10028088631957144 / 1.09330844884037659
V32 to Center: 0.88547208090139673 / 0.73415183910062098
V33 to Center: 1.10028088631957144 / 1.09330844884037659
V34 to Center: 0.11219189000871095 / 0.01257382296550170
CU-V0 to Center: 1.65004457710871555 / 2.15833065217771214
CU-V1 to Center: 1.65004457710871555 / 2.15833065217771214
CU-V2 to Center: 2.48244445722154428 / 3.58102839999664058
CU-V3 to Center: 2.49691898001032797 / 3.59859182265260014
~/src/s2geometry/src/s2/s2cell_union_test.cc:234: Failure
Value of: cap.Contains(S2Cell(id))
  Actual: false
Expected: true
CU-V4 to Center: 1.49154807648107801 / 1.84166934782228831
CU-V5 to Center: 1.49154807648107801 / 1.84166934782228831
CU-V6 to Center: 0.11219189000871095 / 0.01257382296550170

Just chalk it up to that, thanks!