Remove vector lines which are too close to each other?

Question

I'm currently using boost::polygon::detail::resize() function to enlarge or shrink the outline of a polygon by a specific size. This is working well and gives proper results.

Now in some cases, depending on the input shape, the resulting polygon (red) contains vector lines which are very close to each other (black line):

What I want to do is to have a minimum distance over all neighbouring lines of such a polygon. For this example this means, the resulting polygon needs to be modified and a vector line (blue) has to be added, replacing the red lines below of it. This of course would change the shape of the polygon but would keep it conform to the minimum-distance-rule. In this example picture given above, the new, blue line would have the length of the required minimum distance.

My questions: are there any boost-functions available which could do that job? If yes: which ones and how have they to be used?

Thanks!

Answer 1

Converted your image into code using:

That's:

Point
    A(-25.74,    2.5),
    B( 14.96,  12.96),
    C( 31.07,  -2.16),
    D(    37, -27.74),
    E( 29.37, -51.63),
    F( 46.62, -25.48),
    G( 48.17,   3.63),
    H( 41.53,  20.31),
    I( 77.57,  32.04);

String ls{A, B, C, D, E, F, G, H, I};

Using

using Point  = bgm::d2::point_xy<double>;
using String = bgm::linestring<Point>;

---

Let's reproduce the angle measurements:

Live On Coliru

#include <boost/geometry.hpp>
#include <boost/geometry/views/segment_view.hpp>
#include <boost/geometry/geometries/point_xyz.hpp>
#include <iostream>
#include <fstream>

namespace bg = boost::geometry;
namespace bgm = bg::model;

using Point  = bgm::d2::point_xy<double>;
using String = bgm::linestring<Point>;
using bg::wkt;

static inline auto angle_rad(Point a, Point b) {
    bg::subtract_point(a, b);
    return atan2(a.y(), a.x());
}

int main() 
{
    // clang-format off
    Point
        A(-25.74,    2.5),
        B( 14.96,  12.96),
        C( 31.07,  -2.16),
        D(    37, -27.74),
        E( 29.37, -51.63),
        F( 46.62, -25.48),
        G( 48.17,   3.63),
        H( 41.53,  20.31),
        I( 77.57,  32.04);
    // clang-format on

    String ls{A, B, C, D, E, F, G, H, I};
    if (std::string reason; !bg::is_valid(ls, reason)) {
        std::cout << reason << "\n";
        bg::correct(ls);
    }

    std::cout << wkt(ls) << "\n";

    auto desc_p = [&ls](Point const& p) {
        static auto names = "ABCDEFGHIJ";
        assert(not std::less<>{}(&ls.back(), &p));
        assert(not std::less<>{}(&p, &ls.front()));

        return std::string(names + (&p - &ls.front()), 1);
    };
    auto desc = [=](auto it) {
        return desc_p(*it->first) + desc_p(*it->second);
    };

    assert(ls.size() > 1);
    auto a = bg::segments_begin(ls), b = std::next(a), e = bg::segments_end(ls);

    for (; b != e; ++a, ++b) {
        // segment pair *a and *b
        auto pointing_angle = [](auto it) {
            return angle_rad(*it->first, *it->second);
        };
        auto aa = pointing_angle(a);
        auto ab = pointing_angle(b);
        auto deg = [](auto rad) { return rad / M_PI * 180; };
        std::cout << "angle " << deg(fmod(ab - aa + 3 * M_PI, 2 * M_PI))
                  << " between "                         //
                  << desc(a) << " " << deg(aa) << " vs " //
                  << desc(b) << " " << deg(ab) << "\n";
    }

    {
        // Declare a stream and an SVG mapper
        std::ofstream svg("output.svg");
        bg::svg_mapper<Point> mapper(svg, 400, 400);

        // Add geometries such that all these geometries fit on the map
        mapper.add(ls);

        mapper.map(ls, "fill-opacity:0.3;fill:rgb(51,0,0);stroke:rgb(51,0,0);stroke-width:1");
        for (auto& p : ls)
            mapper.text(p, desc_p(p), "");
    }
}

Prints

LINESTRING(-25.74 2.5,14.96 12.96,31.07 -2.16,37 -27.74,29.37 -51.63,46.62 -25.48,48.17 3.63,41.53 20.31,77.57 32.04)
angle 122.402 between AB -165.587 vs BC 136.816
angle 146.236 between BC 136.816 vs CD 103.052
angle 149.236 between CD 103.052 vs DE 72.2875
angle 344.301 between DE 72.2875 vs EF -123.411
angle 210.363 between EF -123.411 vs FG -93.0479
angle 204.754 between FG -93.0479 vs GH -68.2934
angle 86.322 between GH -68.2934 vs HI -161.971

With the side-effect of writing output.svg containing

Performing a cut-off

Detecting a sharp angle, we can check the cutoff length that would be required for the distance between the legs to grow beyond min_distance:

auto seg_angle = [](auto it) {
    return angle_rad(*it->first, *it->second);
};
auto rel      = fmod(seg_angle(b) - seg_angle(a) + 3 * M_PI, 2 * M_PI);
auto inner    = fabs(2 * M_PI - rel);
bool is_sharp = inner < M_PI / 2;

if (is_sharp) {
    std::cout << " ---- sharp angle, min_distance: " << min_distance
              << " (" << color << ")\n";
    auto deg = [](auto rad) { return rad / M_PI * 180; };
    std::cout << "angle " << deg(inner) << " between " << desc(a) << " and " << desc(b) << "\n";

    auto length_a = bg::length(*a);
    auto length_b = bg::length(*b);
    std::cout << "len(" << desc(a) << "): " << length_a << "\n";
    std::cout << "len(" << desc(b) << "): " << length_b << "\n";
    std::cout << "distance(" << desc(P1) << " - " << desc(P3) << "): " << bg::distance(P1, P3) << "\n";

    auto cutoff = min_distance / tan(inner);
    std::cout << "cutoff: " << cutoff << "\n";

Let's create a helper to manually interpolate the cutoff points:

static inline String do_cutoff(String s, long double amount)
{
    assert(s.size() >= 2);
    auto delta = s[1];
    bg::subtract_point(delta, s[0]);

    auto l = bg::length(s);
    if (l > 0) {
        bg::multiply_value(delta, std::min(amount, l) / l);
        bg::add_point(s[0], delta);
    }

    return s;
}

Now we can use it to get new mid-points instead of P2:

    String fs = do_cutoff({P2, P1}, cutoff);
    String ss = do_cutoff({P2, P3}, cutoff);
    std::reverse(fs.begin(), fs.end());

    // assert that end points didn't change
    assert(bg::equals(fs.front(), P1));
    assert(bg::equals(ss.back(), P3));

Now there's the degenerate case where both legs are too short for the cutoff, meaning the mid-point will disappear:

    if (bg::length(*a) < cutoff && bg::length(*b) < cutoff) {
        // both legs vanish, connect P1 to P3 directly
        auto it = ls.begin() + index_of_p2;
        std::cout << "Dropping " << desc(*it) << " (" << desc(P2) << " "
                  << wkt(P2) << ")\n";
        ls.erase(it);

        // note that distance (P1, P3) might still exceed min_distance
        // behaviour here is not specified in question

        break; // invalidated the loop iterators!
    }

Next up, there might be one leg that disappears:

    if (bg::length(*a) < cutoff) { // P1 P2 P3 becomes P1 P3' P3
        update(P2, ss.front());
    } else if (bg::length(*b) < cutoff) { // P1 P2 P3 becomes P1 P1' P3
        update(P2, fs.back());
    } else {

So far so good. If none of the legs disappear, we need to cut-off the sharp point by introducing an extra point (and updating the original mid-point):

        update(P2, ss.front());

        // adding a point (invalidating the segment iterators)
        std::cout << "Adding " << wkt(fs.back()) << " as X "
                  << " before  " << desc(P2) << "\n";

        names.insert(names.begin() + index_of_p2, "X");
        auto added = ls.insert(ls.begin() + index_of_p2, fs.back());
        assert(bg::equals(*added, fs.back()));

        break; // invalidated the loop iterators!
    }

Live Demo

Live On Coliru

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xyz.hpp>
#include <boost/geometry/views/segment_view.hpp>
#include <fstream>
#include <iostream>

namespace bg = boost::geometry;
namespace bgm = bg::model;

using Point   = bgm::d2::point_xy<double>;
using String  = bgm::linestring<Point>;
using bg::wkt;

static inline auto angle_rad(Point a, Point b) {
    bg::subtract_point(a, b);
    return atan2(a.y(), a.x());
}

static inline String do_cutoff(String s, long double amount)
{
    assert(s.size() >= 2);
    auto delta = s[1];
    bg::subtract_point(delta, s[0]);

    auto l = bg::length(s);
    if (l > 0) {
        bg::multiply_value(delta, std::min(amount, l) / l);
        bg::add_point(s[0], delta);
    }

    return s;
}

struct Test {
    String                   ls;
    std::vector<std::string> names;
    double const             min_distance;
    std::string              color, dash_pattern, other_style;

    std::string desc(Point const& p) const
    {
        return names.at(&p - &ls.front());
    };
    std::string desc(bg::segment_iterator<String const> it) const
    {
        return desc(*it->first) + desc(*it->second);
    };

    void perform_cutoff();
};

void Test::perform_cutoff()
{
    if (min_distance == 0)
        return; // keep original

    assert(ls.size() > 1);
    auto a = bg::segments_begin(ls), b = std::next(a), e = bg::segments_end(ls);

    for (; b != e; ++a, ++b) {
        // segments a,b form a joint, let's call it P1 P2 P3
        auto& P1 = *a->first;
        auto& P2 = *a->second;
        auto& P3 = *b->second;
        // keep in mind a->second and b->first alias the same point here
        assert(a->second == b->first);

        auto seg_angle = [](auto it) {
            return angle_rad(*it->first, *it->second);
        };
        auto rel      = fmod(seg_angle(b) - seg_angle(a) + 3 * M_PI, 2 * M_PI);
        auto inner    = fabs(2 * M_PI - rel);
        bool is_sharp = inner < M_PI / 2;

        if (is_sharp) {
            std::cout << " ---- sharp angle, min_distance: " << min_distance
                      << " (" << color << ")\n";
            auto deg = [](auto rad) { return rad / M_PI * 180; };
            std::cout << "angle " << deg(inner) << " between " << desc(a) << " and " << desc(b) << "\n";

            auto length_a = bg::length(*a);
            auto length_b = bg::length(*b);
            std::cout << "len(" << desc(a) << "): " << length_a << "\n";
            std::cout << "len(" << desc(b) << "): " << length_b << "\n";
            std::cout << "distance(" << desc(P1) << " - " << desc(P3) << "): " << bg::distance(P1, P3) << "\n";

            auto cutoff = min_distance / tan(inner);
            std::cout << "cutoff: " << cutoff << "\n";

            String fs = do_cutoff({P2, P1}, cutoff);
            String ss = do_cutoff({P2, P3}, cutoff);
            std::reverse(fs.begin(), fs.end());

            // assert that end points didn't change
            assert(bg::equals(fs.front(), P1));
            assert(bg::equals(ss.back(), P3));

            std::cout << wkt(ls) << "\n";
            auto update = [this](Point const& p, Point const& value) {
                std::cout << "Updating " << desc(p) << " from " << wkt(p)
                          << " to " << wkt(value) << "\n";
                bg::assign(const_cast<Point&>(p), value);
                names[&p - &ls.front()] += "'";
            };

            // For modifying we need the index to the middle point:
            auto const index_of_p2 = &P2 - &ls.front();

            if (bg::length(*a) < cutoff && bg::length(*b) < cutoff) {
                // both legs vanish, connect P1 to P3 directly
                auto it = ls.begin() + index_of_p2;
                std::cout << "Dropping " << desc(*it) << " (" << desc(P2) << " "
                          << wkt(P2) << ")\n";
                ls.erase(it);

                // note that distance (P1, P3) might still exceed min_distance
                // behaviour here is not specified in question

                break; // invalidated the loop iterators!
            }

            if (bg::length(*a) < cutoff) { // P1 P2 P3 becomes P1 P3' P3
                update(P2, ss.front());
            } else if (bg::length(*b) < cutoff) { // P1 P2 P3 becomes P1 P1' P3
                update(P2, fs.back());
            } else {
                update(P2, ss.front());

                // adding a point (invalidating the segment iterators)
                std::cout << "Adding " << wkt(fs.back()) << " as X "
                          << " before  " << desc(P2) << "\n";

                names.insert(names.begin() + index_of_p2, "X");
                auto added = ls.insert(ls.begin() + index_of_p2, fs.back());
                assert(bg::equals(*added, fs.back()));

                break; // invalidated the loop iterators!
            }
        }
    }
}
int main() 
{
    // clang-format off
    Point
        A(-25.74,    2.5),
        B( 14.96,  12.96),
        C( 31.07,  -2.16),
        D(    37, -27.74),
        E( 29.37, -51.63),
        F( 46.62, -25.48),
        G( 48.17,   3.63),
        H( 41.53,  20.31),
        I( 77.57,  32.04);
    // clang-format on

    {
        String const original_ls{A, B, C, D, E, F, G, H, I};

        if (std::string reason; !bg::is_valid(original_ls, reason)) {
            std::cout << reason << "\n";
            return 1;
        }
    }

    Test testcases[]{
        {
            {A, B, C, D, E, F, G, H, I},
                {"A", "B", "C", "D", "E", "F", "G", "H", "I", "J"},
                0, // keep original
                "grey", "1 1",
        },
            {
                {A, B, C, D, E, F, G, H, I},
                {"A", "B", "C", "D", "E", "F", "G", "H", "I", "J"},
                5,
                "green", "2 1",
            },
            {
                {A, B, C, D, E, F, G, H, I},
                {"A", "B", "C", "D", "E", "F", "G", "H", "I", "J"},
                8,
                "red", "1 2",
            },
            {
                {A, B, C, D, E, F, G, H, I},
                {"A", "B", "C", "D", "E", "F", "G", "H", "I", "J"},
                12,
                "blue", "2 3",
            },
    };

    std::ofstream svg("output.svg");
    bg::svg_mapper<Point> mapper(svg, 400, 400);

    for (auto& c : testcases) {
        c.perform_cutoff();
        std::cout << wkt(c.ls) << "\n";

        mapper.add(c.ls);
        mapper.map(
                c.ls,
                "stroke:" + c.color + ";stroke-width:1;stroke-dasharray:" + c.dash_pattern);
        for (auto& p : c.ls)
            mapper.text(p, c.desc(p), "fill:" + c.color + ";font-size:x-small");
    }
}

Prints

LINESTRING(-25.74 2.5,14.96 12.96,31.07 -2.16,37 -27.74,29.37 -51.63,46.62 -25.48,48.17 3.63,41.53 20.31,77.57 32.04)
 ---- sharp angle, min_distance: 5 (green)
angle 15.6986 between DE and EF
len(DE): 25.0789
len(EF): 31.3271
distance(D - F): 9.8819
cutoff: 17.7897
LINESTRING(-25.74 2.5,14.96 12.96,31.07 -2.16,37 -27.74,29.37 -51.63,46.62 -25.48,48.17 3.63,41.53 20.31,77.57 32.04)
Updating E from POINT(29.37 -51.63) to POINT(39.1658 -36.7802)
Adding POINT(34.7824 -34.6836) as X  before  E'
LINESTRING(-25.74 2.5,14.96 12.96,31.07 -2.16,37 -27.74,34.7824 -34.6836,39.1658 -36.7802,46.62 -25.48,48.17 3.63,41.53 20.31,77.57 32.04)
 ---- sharp angle, min_distance: 8 (red)
angle 15.6986 between DE and EF
len(DE): 25.0789
len(EF): 31.3271
distance(D - F): 9.8819
cutoff: 28.4636
LINESTRING(-25.74 2.5,14.96 12.96,31.07 -2.16,37 -27.74,29.37 -51.63,46.62 -25.48,48.17 3.63,41.53 20.31,77.57 32.04)
Updating E from POINT(29.37 -51.63) to POINT(45.0432 -27.8703)
LINESTRING(-25.74 2.5,14.96 12.96,31.07 -2.16,37 -27.74,45.0432 -27.8703,46.62 -25.48,48.17 3.63,41.53 20.31,77.57 32.04)
 ---- sharp angle, min_distance: 12 (blue)
angle 15.6986 between DE and EF
len(DE): 25.0789
len(EF): 31.3271
distance(D - F): 9.8819
cutoff: 42.6953
LINESTRING(-25.74 2.5,14.96 12.96,31.07 -2.16,37 -27.74,29.37 -51.63,46.62 -25.48,48.17 3.63,41.53 20.31,77.57 32.04)
Dropping E (E POINT(29.37 -51.63))
LINESTRING(-25.74 2.5,14.96 12.96,31.07 -2.16,37 -27.74,46.62 -25.48,48.17 3.63,41.53 20.31,77.57 32.04)

And the resulting output.svg:

Conclusion

Hopefully that helps.

Note that the code has limitations as written, because for simplicity I used segment_iterator.

• Due to the segment-iterator loop, currently only the first sharp corner might be treated, unless they can be treated without invalidating the iterators. You'd either have to repeat the loop or rewrite the loop in terms of indices instead of using the segment_iterators.

• Also due to the segment iterators, I resorted to an ugly const_cast to update existing points, instead of going throught the geometry interface. This would not pass my own code review :)

• Finally, nothing is decided when even the distance P1-P3 directly exceeded the min_distance (you didn't specify anything about this potentiality in the question).

---