google / s2geometry

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

S2Builder adding extra points to overlapping polygons #351

Open MattOslin opened 6 months ago

MattOslin commented 6 months ago

Hello, I'm having some issues using the S2Builder to correct input geometry.

The input is as follows (in e6 lat/lng, the input has repeated terminal point and can be oriented either way): PolygonA = (37.65026, -122.400647), (37.649969, -122.400399), (37.650271, -122.399859), (37.650581, -122.400093), (37.65026, -122.400647) PolygonB = (37.65058, -122.400095), (37.650271, -122.39986), (37.650429, -122.399543), (37.650714, -122.399845), (37.65058, -122.400095)

I would like to

To do this, I first create a builder with the e7 snap function and EdgeType::UNDIRECTED. However, when I add each polygon as its own layer and run, I get

PolygonA = (37.65026, -122.400647), (37.649969, -122.400399), (37.650271, -122.399859), (37.650581, -122.400093), (37.650578, -122.400095) PolygonB = (37.65058, -122.400095), (37.650271, -122.39986), (37.650271, -122.399859), (37.650429, -122.399543), (37.650714, -122.399845), (37.650581, -122.400093)

Why is the builder adding extra points to each polygon? I have my own edge welding logic I can do after running the builder to identify edges that are close to be joined together, but this introduction of extra points is pesky. Why are they added, when they only complicate the individual polygons? I can't see any reason in the API that this would occur.

I have also tried increasing the snap radius of the snap function, but these extra points still get added. Setting simplify_edge_chains to true removes the introduced points, but it doesn't seem like I should have to go to such extremes.

Also, is there a convenient way to weld these points together in the builder, without having to identify and weld them later? Thanks!

smcallis commented 6 months ago

Can you post the code you're using to call S2Builder?

MattOslin commented 6 months ago

Can you post the code you're using to call S2Builder?

You could imagine that I'm doing something like the following:

using Polygons = absl::flat_hash_map<std::string_view, std::unique_ptr<S2Polygon>>;

struct InputPolygon {
    std::string_view id = {};
    std::vector<std::array<double, 2>> points = {};
};

std::vector<S2Point> MakeLoop(const std::vector<std::array<double, 2>>& polygon) {
  return polygon | std::views::transform([](const auto& pt) {
           return S2LatLng::FromDegrees(pt[0], pt[1]).ToPoint();
         }) | ToVector{};
}

Polygons ExtractPolygons(const std::vector<InputPolygon>& input_polgons) {
  auto builder = S2Builder{S2Builder::Options{s2builderutil::IntLatLngSnapFunction(7)}};

  auto ret = Polygons{};
  ret.reserve(input_polgons.size());
  for (const auto& input_polygon : input_polgons) {
    const auto [iter, _] = ret.emplace(input_polygon.id, std::make_unique<S2Polygon>());
    builder.StartLayer(std::make_unique<s2builderutil::S2PolygonLayer>(
        iter->second.get(), s2builderutil::S2PolygonLayer::Options{S2Builder::EdgeType::UNDIRECTED}));
    builder.AddLoop(MakeLoop(input_polygon.points));
  }

  auto error = S2Error{};
  if (!builder.Build(&error)) std::cout << error.text() << std::endl;

  absl::erase_if(ret, [](const auto& key_value) { return key_value.second->is_empty(); });

  return ret;
}
smcallis commented 6 months ago

Hmmm doesn't seem to reproduce for me. If I run this:

...

using s2builderutil::S2PolygonLayer;

int main() {
  std::vector<S2Point> poly0 = s2textformat::ParsePointsOrDie(
      "37.65026:-122.400647, 37.649969:-122.400399, 37.650271:-122.399859, "
      "37.650581:-122.400093, 37.65026:-122.400647");
  std::vector<S2Point> poly1 = s2textformat::ParsePointsOrDie(
      "37.65058:-122.400095, 37.650271:-122.39986, 37.650429:-122.399543, "
      "37.650714:-122.399845, 37.65058:-122.400095");

  absl::PrintF("input0: %s\n", s2textformat::ToString(poly0));
  absl::PrintF("input1: %s\n", s2textformat::ToString(poly1));

  S2Builder builder(
      S2Builder::Options(s2builderutil::IntLatLngSnapFunction(7)));

  S2Polygon out0, out1;

  builder.StartLayer(std::make_unique<S2PolygonLayer>(
      &out0, S2PolygonLayer::Options(S2Builder::EdgeType::UNDIRECTED)));
  builder.AddLoop(poly0);

  builder.StartLayer(std::make_unique<S2PolygonLayer>(
      &out1, S2PolygonLayer::Options(S2Builder::EdgeType::UNDIRECTED)));
  builder.AddLoop(poly1);

  S2Error error;
  if (!builder.Build(&error)) {
    LOG(FATAL) << error.text();
  }

  absl::PrintF("output0: %s\n", s2textformat::ToString(out0));
  absl::PrintF("output1: %s\n", s2textformat::ToString(out1));
}

The output is:

input0: 37.65026:-122.400647, 37.649969:-122.400399, 37.650271:-122.399859, 37.650581:-122.400093, 37.65026:-122.400647
input1: 37.65058:-122.400095, 37.650271:-122.39986, 37.650429:-122.399543, 37.650714:-122.399845, 37.65058:-122.400095
output0: 37.65026:-122.400647, 37.649969:-122.400399, 37.650271:-122.399859, 37.650581:-122.400093
output1: 37.65058:-122.400095, 37.650271:-122.39986, 37.650429:-122.399543, 37.650714:-122.399845, 37.650581:-122.400093  

These are the two input polygons:

image

And the corresponding outputs:

image

The only change I see is that the topmost blue edge was free of the red polygon before:

image

And was snapped after building:

image

If you can put together an MRE I'm happy to take another look.

MattOslin commented 6 months ago

Hmmm doesn't seem to reproduce for me. If I run this:

...

using s2builderutil::S2PolygonLayer;

int main() {
  std::vector<S2Point> poly0 = s2textformat::ParsePointsOrDie(
      "37.65026:-122.400647, 37.649969:-122.400399, 37.650271:-122.399859, "
      "37.650581:-122.400093, 37.65026:-122.400647");
  std::vector<S2Point> poly1 = s2textformat::ParsePointsOrDie(
      "37.65058:-122.400095, 37.650271:-122.39986, 37.650429:-122.399543, "
      "37.650714:-122.399845, 37.65058:-122.400095");

  absl::PrintF("input0: %s\n", s2textformat::ToString(poly0));
  absl::PrintF("input1: %s\n", s2textformat::ToString(poly1));

  S2Builder builder(
      S2Builder::Options(s2builderutil::IntLatLngSnapFunction(7)));

  S2Polygon out0, out1;

  builder.StartLayer(std::make_unique<S2PolygonLayer>(
      &out0, S2PolygonLayer::Options(S2Builder::EdgeType::UNDIRECTED)));
  builder.AddLoop(poly0);

  builder.StartLayer(std::make_unique<S2PolygonLayer>(
      &out1, S2PolygonLayer::Options(S2Builder::EdgeType::UNDIRECTED)));
  builder.AddLoop(poly1);

  S2Error error;
  if (!builder.Build(&error)) {
    LOG(FATAL) << error.text();
  }

  absl::PrintF("output0: %s\n", s2textformat::ToString(out0));
  absl::PrintF("output1: %s\n", s2textformat::ToString(out1));
}

The output is:

input0: 37.65026:-122.400647, 37.649969:-122.400399, 37.650271:-122.399859, 37.650581:-122.400093, 37.65026:-122.400647
input1: 37.65058:-122.400095, 37.650271:-122.39986, 37.650429:-122.399543, 37.650714:-122.399845, 37.65058:-122.400095
output0: 37.65026:-122.400647, 37.649969:-122.400399, 37.650271:-122.399859, 37.650581:-122.400093
output1: 37.65058:-122.400095, 37.650271:-122.39986, 37.650429:-122.399543, 37.650714:-122.399845, 37.650581:-122.400093  

These are the two input polygons:

image

And the corresponding outputs:

image

The only change I see is that the topmost blue edge was free of the red polygon before:

image

And was snapped after building:

image

If you can put together an MRE I'm happy to take another look.

Thanks for taking a look! That makes sense that the edge snapping is what is causing extra points to be added. I will try setting min_edge_vertex_separation to zero and see if that helps.