Boost Polygon: Issue with euclidean_distance

Question

I have the following code which is supposed to compute the Euclidean distance between two rectangles. I compiled using GCC 4.7.3 and Boost v1.58.0

#include <iostream>
#include <cmath>
#include <boost/polygon/polygon.hpp>
#include <boost/geometry.hpp>

namespace gtl = boost::polygon;
using namespace boost::polygon::operators;

typedef gtl::rectangle_data<int> LayoutRectangle;

int main(int argc, char** argv)
{
    LayoutRectangle t(16740130,29759232,16740350,29760652);
    LayoutRectangle n(16808130,29980632,16808350,29982052);

    std::cout << gtl::euclidean_distance(t, n) << std::endl;

    std::cout << gtl::euclidean_distance(t, n, gtl::HORIZONTAL) << " "
              << gtl::euclidean_distance(t, n, gtl::VERTICAL) << std::endl;

    std::cout << gtl::square_euclidean_distance(t, n) << std::endl;
    std::cout << std::sqrt(gtl::square_euclidean_distance(t, n)) << std::endl;
    std::cout << (int) std::sqrt(gtl::square_euclidean_distance(t, n)) << std::endl;

    return 0;
}

The code above produced the following output:

38022.6
67780 219980
52985328800
230185
230185

The correct answer is 230185. Now if I go look at the implementation of euclidean_distance() in the boost polygon library, I see this:

  template <typename rectangle_type, typename rectangle_type_2>
  typename enable_if< typename gtl_and_3<y_r_edist2, typename is_rectangle_concept<typename geometry_concept<rectangle_type>::type>::type,
                                                          typename is_rectangle_concept<typename geometry_concept<rectangle_type_2>::type>::type>::type,
                       typename rectangle_distance_type<rectangle_type>::type>::type
  euclidean_distance(const rectangle_type& lvalue, const rectangle_type_2& rvalue) {
    double val = (int)square_euclidean_distance(lvalue, rvalue);
    return std::sqrt(val);
  }

This looks identical to the std::sqrt(gtl::square_eclidean_distance(t,n)) line in my code which gives the correct answer (230185). So why am I getting 38022.6 with gtl::euclidean_distance()? What am I not seeing here?

Answer 1

Looks like the internal computation is overflowing. I don't think this is a library bug, the library is used incorrectly with the underlying (unchecked) int type. (However, there is a different bug in the library that I mention at the end.)

Try using a smaller "integer representation" of the problem:

For example:

LayoutRectangle t(167402,297592,167404,297607);
LayoutRectangle n(168082,299806,168084,299821);

Unfortunately there is no general solution of the problem in integer arithmetic, except 0) using higher precision can buy you something, 1) scaling the problem 2) using multiprecision, 3) using rational arithmetic and integer part

(For floating point the solution is simply normalizing the components, this is how std::abs for std::complex<double> works to avoid floating point overflow)

It is good to use large integers to represent a geometric problem BUT for this reason, as a workaround, use coordinates that span distance of at most (int)std::sqrt((double)std::numeric_limits<int>::max()/2) = 2^15 = 32768. Which is a surprisingly small number.

Complete code:

#include <iostream>
#include <cmath>
#include <boost/polygon/polygon.hpp>
#include <boost/geometry.hpp>

int main(){

namespace gtl = boost::polygon;
using namespace boost::polygon::operators;

typedef gtl::rectangle_data<int> LayoutRectangle;

    LayoutRectangle t(167401,297592,167403,297606);
    LayoutRectangle n(168081,299806,168083,299820);

    std::cout << gtl::euclidean_distance(t, n) << std::endl;

    std::cout << gtl::euclidean_distance(t, n, gtl::HORIZONTAL) << " "
              << gtl::euclidean_distance(t, n, gtl::VERTICAL) << std::endl;

    std::cout << gtl::square_euclidean_distance(t, n) << std::endl;
    std::cout << std::sqrt(gtl::square_euclidean_distance(t, n)) << std::endl;
    std::cout << (int) std::sqrt(gtl::square_euclidean_distance(t, n)) << std::endl;

}

Output:

2302.1
678 2200
5299684
2302.1
2302

Which is the expected result.

---

Looking at the code, it seems that there is a bug in the library, not because it gives overflow but because an internal computation is casted to int and not the the underlying generic integer data type. This means that probably even if you use multiprecision integers the results will overflow.