How to find two points that form closest distance between two rectangles?

Question

I'm trying to find algorithm that will find two points that represent closest distance between two rectangles. Like points C and J that form smallest distance on image below:

I'm trying to not reinvent the wheel here and use something that is already battletested like boost::geometry::distance, but it only return distance and not also the points.

Answer 1

Not making things overly generic (by assuming floating point coordinates and cartesian coordinate system), here's an implementation of point-to-line-segment distance that returns the projected point as well as the distance:

struct DistancePoint {
    double distance;
    P projected_point;
};

template <typename Strategy = bg::strategy::distance::pythagoras<> >
DistancePoint point_to_segment(P const& p, P const& p1, P const& p2) {
    P v = p2, w = p;
    bg::subtract_point(v, p1);
    bg::subtract_point(w, p1);

    auto const c1 = bg::dot_product(w, v);
    if (c1 <= 0)  return { Strategy::apply(p, p1), p1 };

    auto const c2 = bg::dot_product(v, v);
    if (c2 <= c1) return { Strategy::apply(p, p2), p2 };

    P prj = p1;
    bg::multiply_value(v, c1/c2);
    bg::add_point(prj, v);

    return { Strategy::apply(p, prj), prj };
}

Now you can use it with your geometries. I didn't want to satisfy all the criteria for the distance-strategy concept, so you can't use the above with boost::geometry::distance.

But since your input geometries have low point counts, you might get away with "brute-forcing" (without requiring the library internal closest_feature selection):

R a = gen_rect(),
  b = gen_rect();

// make sure a and b don't overlap (distance > 0)
while (!bg::disjoint(a,b)) { b = gen_rect(); }

std::cout
    << wkt(a) << "\n"
    << wkt(b) << "\n"
    << bg::distance(a, b) << " apart\n";

DistancePoint nearest;
P const* which = nullptr;

for (auto& [a,b] : { std::tie(a,b), std::tie(b,a) } ) {
    auto segments = boost::make_iterator_range(bg::segments_begin(a), bg::segments_end(a));
    auto points   = boost::make_iterator_range(bg::points_begin(b), bg::points_end(b));

    for (auto&& pq : segments) {
        for (auto&& r : points) {
            auto d = point_to_segment(r, *pq.first, *pq.second);

            if (!which || d.distance < nearest.distance) {
                which = &r;
                nearest = d;
            }
        }
    }
}

std::cout << wkt(which) << " at " << nearest.distance << " from " << wkt(nearest.projected_point) << "\n";

DEMO

Live On Coliru

#include <boost/geometry.hpp>
#include <boost/geometry/util/range.hpp>
#include <boost/geometry/io/io.hpp>
#include <boost/geometry/algorithms/assign.hpp>
#include <boost/geometry/algorithms/convex_hull.hpp>
#include <boost/geometry/algorithms/point_on_surface.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/geometry/geometries/box.hpp>
#include <boost/range/adaptors.hpp>
#include <iostream>
#include <fstream>

namespace bg = boost::geometry;
namespace bgm = boost::geometry::model;
using namespace boost::adaptors;
using bg::wkt;

using P = bgm::d2::point_xy<double>;
using B = bgm::box<P>;
using S = bgm::segment<P>;
using R = bgm::polygon<P>;

R gen_rect(); // generates a random rectangle

namespace {
    struct DistancePoint {
        double distance;
        P projected_point;
    };

    // after strategy::distance::projected_point<>
    template <typename Strategy = bg::strategy::distance::pythagoras<> >
    DistancePoint point_to_segment(P const& p, P const& p1, P const& p2) {
        P v = p2, w = p;
        bg::subtract_point(v, p1);
        bg::subtract_point(w, p1);

        auto const c1 = bg::dot_product(w, v);
        if (c1 <= 0)  return { Strategy::apply(p, p1), p1 };

        auto const c2 = bg::dot_product(v, v);
        if (c2 <= c1) return { Strategy::apply(p, p2), p2 };

        P prj = p1;
        bg::multiply_value(v, c1/c2);
        bg::add_point(prj, v);

        return { Strategy::apply(p, prj), prj };
    }
}

int main() {
    std::cout << std::setprecision(2);

    for (auto i = 0; i<10; ++i) {
        R a = gen_rect(),
          b = gen_rect();

        // make sure a and b don't overlap (distance > 0)
        while (!bg::disjoint(a,b)) { b = gen_rect(); }

        std::cout
            << wkt(a) << "\n"
            << wkt(b) << "\n"
            << bg::distance(a, b) << " apart\n";

        DistancePoint nearest;
        P const* which = nullptr;

        for (auto& [a,b] : { std::tie(a,b), std::tie(b,a) } ) {
            auto segments = boost::make_iterator_range(bg::segments_begin(a), bg::segments_end(a));
            auto points   = boost::make_iterator_range(bg::points_begin(b), bg::points_end(b));

            for (auto&& pq : segments) {
                for (auto&& r : points) {
                    auto d = point_to_segment(r, *pq.first, *pq.second);

                    if (!which || d.distance < nearest.distance) {
                        which = &r;
                        nearest = d;
                    }
                }
            }
        }

        std::cout << wkt(which) << " at " << nearest.distance << " from " << wkt(nearest.projected_point) << "\n";

        {
            std::ofstream svg("output" + std::to_string(i) + ".svg");
            boost::geometry::svg_mapper<P> mapper(svg, 400, 400, "style='fill-opacity:1;fill:rgb(255,255,255)'");
            mapper.add(a);
            mapper.add(b);
            S dline {*which, nearest.projected_point};
            mapper.add(dline);

            mapper.map(a, "fill-opacity:0.5;fill:rgb(153,204,0);stroke:rgb(153,204,0);stroke-width:2");
            mapper.map(b, "fill-opacity:0.5;fill:rgb(204,153,0);stroke:rgb(202,153,0);stroke-width:2");
            mapper.map(dline, "stroke-dasharray:1,1;stroke:rgb(255,0,0);stroke-width:1");
        }
    }
}

// details for generating the rectangles
#include <boost/geometry/strategies/transform/matrix_transformers.hpp>
#include <random>
std::mt19937 prng { std::random_device{}() };

static auto rand(double b, double e) { 
    return std::uniform_real_distribution<double>(b, e)(prng);
}

R gen_rect() {
    B initial {{0, 0}, { rand(0.1, 1), rand(0.1, 1) } };
    R raw, rect; // todo rotate and stuff
    bg::assign(raw, initial);

    using namespace bg::strategy::transform;
    auto rot   = rand(-M_PI, +M_PI);
    auto scale = rand(1, 3);
    auto x     = rand(-5, 5),
         y     = rand(-5, 5);

    matrix_transformer<double, 2, 2> xfrm(
         scale* cos(rot), scale*sin(rot), x,
         scale*-sin(rot), scale*cos(rot), y,
                       0,              0, 1);

    bg::transform(raw, rect, xfrm);
    bg::correct(rect);
    return rect;
}

This generates some randomized scenes like: