More consistent definition of precision for double types

The double types (DCoord, DPoint, DVector ...) support a
more consistent equivalence criterion now. This supports
micron-unit objects without loss of precision.

Before this fix, the precision used for 0.01 which was
basically implying a resolution of 10nm for micrometer-unit
double objects.
This commit is contained in:
Matthias Koefferlein 2017-07-04 23:50:30 +02:00
parent 18f79dde67
commit 8910828726
5 changed files with 256 additions and 41 deletions

View File

@ -738,17 +738,19 @@ public:
bool res = true;
bool ends_on_edge = false;
area_type vxa = coord_traits::vprod (e.p2 ().x (), e.p2 ().y (), m_p1.x (), m_p1.y (), e.p1 ().x (), e.p1 ().y ());
if (vxa <= -coord_traits::prec_area ()) {
std::pair<area_type, int> vsa = coord_traits::vprod_with_sign (e.p2 ().x (), e.p2 ().y (), m_p1.x (), m_p1.y (), e.p1 ().x (), e.p1 ().y ());
area_type vxa = vsa.first;
if (vsa.second < 0) {
res = false;
} else if (vxa < coord_traits::prec_area ()) {
} else if (vsa.second == 0) {
ends_on_edge = true;
}
area_type vxb = -coord_traits::vprod (e.p2 ().x (), e.p2 ().y (), m_p2.x (), m_p2.y (), e.p1 ().x (), e.p1 ().y ());
if (vxb <= -coord_traits::prec_area ()) {
std::pair<area_type, int> vsb = coord_traits::vprod_with_sign (e.p2 ().x (), e.p2 ().y (), m_p2.x (), m_p2.y (), e.p1 ().x (), e.p1 ().y ());
area_type vxb = -vsb.first;
if (vsb.second > 0) {
res = !res;
} else if (vxb < coord_traits::prec_area ()) {
} else if (vsb.second == 0) {
ends_on_edge = true;
}
@ -833,14 +835,7 @@ public:
return 0;
} else {
// compute the side as the sign of the distance as in "distance"
area_type axb = coord_traits::vprod (m_p2.x (), m_p2.y (), p.x (), p.y (), m_p1.x (), m_p1.y ());
if (axb >= coord_traits::prec_area ()) {
return 1;
} else if (axb <= -coord_traits::prec_area ()) {
return -1;
} else {
return 0;
}
return coord_traits::vprod_sign (m_p2.x (), m_p2.y (), p.x (), p.y (), m_p1.x (), m_p1.y ());
}
}
@ -1007,17 +1002,17 @@ public:
bool res = true;
area_type vxa = coord_traits::vprod (m_p2.x (), m_p2.y (), e.p1 ().x (), e.p1 ().y (), m_p1.x (), m_p1.y ());
if (vxa <= -coord_traits::prec_area ()) {
int vsa = coord_traits::vprod_sign (m_p2.x (), m_p2.y (), e.p1 ().x (), e.p1 ().y (), m_p1.x (), m_p1.y ());
if (vsa < 0) {
res = false;
} else if (vxa < coord_traits::prec_area ()) {
} else if (vsa == 0) {
return true;
}
area_type vxb = -coord_traits::vprod (m_p2.x (), m_p2.y (), e.p2 ().x (), e.p2 ().y (), m_p1.x (), m_p1.y ());
if (vxb <= -coord_traits::prec_area ()) {
int vsb = coord_traits::vprod_sign (m_p2.x (), m_p2.y (), e.p2 ().x (), e.p2 ().y (), m_p1.x (), m_p1.y ());
if (vsb > 0) {
res = !res;
} else if (vxb < coord_traits::prec_area ()) {
} else if (vsb == 0) {
return true;
}
@ -1036,17 +1031,19 @@ public:
{
bool res = true;
area_type vxa = coord_traits::vprod (m_p2.x (), m_p2.y (), e.p1 ().x (), e.p1 ().y (), m_p1.x (), m_p1.y ());
if (vxa <= -coord_traits::prec_area ()) {
std::pair<area_type, int> vsa = coord_traits::vprod_with_sign (m_p2.x (), m_p2.y (), e.p1 ().x (), e.p1 ().y (), m_p1.x (), m_p1.y ());
area_type vxa = vsa.first;
if (vsa.second < 0) {
res = false;
} else if (vxa < coord_traits::prec_area ()) {
} else if (vsa.second == 0) {
return std::make_pair (true, e.p1 ());
}
area_type vxb = -coord_traits::vprod (m_p2.x (), m_p2.y (), e.p2 ().x (), e.p2 ().y (), m_p1.x (), m_p1.y ());
if (vxb <= -coord_traits::prec_area ()) {
std::pair<area_type, int> vsb = coord_traits::vprod_with_sign (m_p2.x (), m_p2.y (), e.p2 ().x (), e.p2 ().y (), m_p1.x (), m_p1.y ());
area_type vxb = -vsb.first;
if (vsb.second > 0) {
res = !res;
} else if (vxb < coord_traits::prec_area ()) {
} else if (vsb.second == 0) {
return std::make_pair (true, e.p2 ());
}
@ -1088,9 +1085,10 @@ public:
*/
std::pair <bool, db::point<C> > cut_point (const db::edge<C> &e2) const
{
double pr1 = double (coord_traits::vprod (e2.p1 ().x (), e2.p1 ().y (), this->p2 ().x (), this->p2 ().y (), this->p1 ().x (), this->p1 ().y ()));
double pr2 = double (coord_traits::vprod (e2.dx (), e2.dy (), this->dx (), this->dy (), 0, 0));
if (fabs (pr2) > double (coord_traits::prec_area ())) {
std::pair<typename coord_traits::area_type, int> vps = coord_traits::vprod_with_sign (e2.dx (), e2.dy (), this->dx (), this->dy (), 0, 0);
if (vps.second != 0) {
double pr1 = double (coord_traits::vprod (e2.p1 ().x (), e2.p1 ().y (), this->p2 ().x (), this->p2 ().y (), this->p1 ().x (), this->p1 ().y ()));
double pr2 = double (vps.first);
db::point<C> p = e2.p1 () - db::vector<C> ((e2.p2 () - e2.p1 ()) * (pr1 / pr2));
return std::make_pair (true, p);
} else {

View File

@ -128,6 +128,62 @@ compute_normals (const db::vector<C> &d, C dx, C dy, int nsign, db::DVector &ed,
}
}
/**
* @brief Provides a special DVector vprod sign for the purpose of representing integer-coordinate vectors
* The "zero" criterion is somewhat tigher than that of the normal integer value vectors.
* Hence, parallelity is somewhat more strict which makes the size function produce a
* better approximation to the desired target contour.
*/
static inline int
vprod_sign_for (const db::DVector &a, const db::DVector &b, const db::Vector &)
{
double vp = db::vprod (a, b);
if (vp <= -1e-2) {
return -1;
} else if (vp < 1e-2) {
return 0;
} else {
return 1;
}
}
/**
* @brief Fallback to the default vprod sign in the double-coordinate case
*/
static inline int
vprod_sign_for (const db::DVector &a, const db::DVector &b, const db::DVector &)
{
return db::vprod_sign (a, b);
}
/**
* @brief Provides a special DVector sprod sign for the purpose of representing integer-coordinate vectors
* The "zero" criterion is somewhat tigher than that of the normal integer value vectors.
* Hence, orthogonality is somewhat more strict which makes the size function produce a
* better approximation to the desired target contour.
*/
static inline int
sprod_sign_for (const db::DVector &a, const db::DVector &b, const db::Vector &)
{
double sp = db::sprod (a, b);
if (sp <= -1e-2) {
return -1;
} else if (sp < 1e-2) {
return 0;
} else {
return 1;
}
}
/**
* @brief Fallback to the default sprod sign in the double-coordinate case
*/
static inline int
sprod_sign_for (const db::DVector &a, const db::DVector &b, const db::DVector &)
{
return db::sprod_sign (a, b);
}
template <class C>
void polygon_contour<C>::size (C dx, C dy, unsigned int mode)
{
@ -192,7 +248,7 @@ void polygon_contour<C>::size (C dx, C dy, unsigned int mode)
db::DVector eed, nnd;
compute_normals (dd, dx, dy, nsign, eed, nnd);
int vpsign = db::vprod_sign (eed, ed) * nsign;
int vpsign = vprod_sign_for (eed, ed, dd) * nsign;
if (vpsign <= 0) {
@ -210,7 +266,7 @@ void polygon_contour<C>::size (C dx, C dy, unsigned int mode)
*pts++ = *pp + vector<C> (nd);
*pts++ = *pp;
} else if (vpsign == 0 && db::sprod_sign (nd, nnd) > 0) {
} else if (vpsign == 0 && sprod_sign_for (nd, nnd, dd) > 0) {
// colinear edges: simply shift the point
*pts++ = *pp + vector<C> (nd);

View File

@ -27,6 +27,7 @@
#include <stdint.h>
#include <math.h>
#include <stdio.h>
#include <algorithm>
namespace db {
@ -247,7 +248,25 @@ struct generic_coord_traits
}
}
/**
/**
* @brief A combination of sign and value of sprod.
*/
static std::pair<area_type, int> sprod_with_sign (coord_type ax, coord_type ay,
coord_type bx, coord_type by,
coord_type cx, coord_type cy)
{
area_type p1 = ((area_type) ax - (area_type) cx) * ((area_type) bx - (area_type) cx);
area_type p2 = -(((area_type) ay - (area_type) cy) * ((area_type) by - (area_type) cy));
if (p1 > p2) {
return std::make_pair (p1 - p2, 1);
} else if (p1 == p2) {
return std::make_pair (0, 0);
} else {
return std::make_pair (p1 - p2, -1);
}
}
/**
* @brief the vector product of two vectors.
*
* Computes the vector product of two vectors.
@ -303,6 +322,23 @@ struct generic_coord_traits
}
}
/**
* @brief A combination of sign and value of vprod.
*/
static std::pair<area_type, int> vprod_with_sign (coord_type ax, coord_type ay,
coord_type bx, coord_type by,
coord_type cx, coord_type cy)
{
area_type p1 = ((area_type) ax - (area_type) cx) * ((area_type) by - (area_type) cy);
area_type p2 = ((area_type) ay - (area_type) cy) * ((area_type) bx - (area_type) cx);
if (p1 > p2) {
return std::make_pair (p1 - p2, 1);
} else if (p1 == p2) {
return std::make_pair (0, 0);
} else {
return std::make_pair (p1 - p2, -1);
}
}
};
/**
@ -382,19 +418,39 @@ struct coord_traits<double>
return (ax - cx) * (bx - cx) + (ay - cy) * (by - cy);
}
static int sprod_sign (double ax, double ay, double bx, double by, double cx, double cy)
{
area_type p1 = (ax - cx) * (bx - cx);
area_type p2 = -(ay - cy) * (by - cy);
if (p1 <= p2 - prec_area ()) {
double dx1 = ax - cx, dy1 = ay - cy;
double dx2 = bx - cx, dy2 = by - cy;
double pa = (sqrt (dx1 * dx1 + dy1 * dy1) + sqrt (dx2 * dx2 + dy2 * dy2)) * prec ();
area_type p1 = dx1 * dx2;
area_type p2 = -dy1 * dy2;
if (p1 <= p2 - pa) {
return -1;
} else if (p1 < p2 + prec_area ()) {
} else if (p1 < p2 + pa) {
return 0;
} else {
return 1;
}
}
static std::pair<area_type, int> sprod_with_sign (double ax, double ay, double bx, double by, double cx, double cy)
{
double dx1 = ax - cx, dy1 = ay - cy;
double dx2 = bx - cx, dy2 = by - cy;
double pa = (sqrt (dx1 * dx1 + dy1 * dy1) + sqrt (dx2 * dx2 + dy2 * dy2)) * prec ();
area_type p1 = dx1 * dx2;
area_type p2 = -dy1 * dy2;
if (p1 <= p2 - pa) {
return std::make_pair (p1 - p2, -1);
} else if (p1 < p2 + pa) {
return std::make_pair (0, 0);
} else {
return std::make_pair (p1 - p2, 1);
}
}
static area_type vprod (coord_type ax, coord_type ay,
coord_type bx, coord_type by,
coord_type cx, coord_type cy)
@ -404,15 +460,34 @@ struct coord_traits<double>
static int vprod_sign (double ax, double ay, double bx, double by, double cx, double cy)
{
area_type p1 = (ax - cx) * (by - cy);
area_type p2 = (ay - cy) * (bx - cx);
if (p1 <= p2 - prec_area ()) {
double dx1 = ax - cx, dy1 = ay - cy;
double dx2 = bx - cx, dy2 = by - cy;
double pa = (sqrt (dx1 * dx1 + dy1 * dy1) + sqrt (dx2 * dx2 + dy2 * dy2)) * prec ();
area_type p1 = dx1 * dy2;
area_type p2 = dy1 * dx2;
if (p1 <= p2 - pa) {
return -1;
} else if (p1 < p2 + prec_area ()) {
} else if (p1 < p2 + pa) {
return 0;
} else {
return 1;
}
}
}
static std::pair<area_type, int> vprod_with_sign (double ax, double ay, double bx, double by, double cx, double cy)
{
double dx1 = ax - cx, dy1 = ay - cy;
double dx2 = bx - cx, dy2 = by - cy;
double pa = (sqrt (dx1 * dx1 + dy1 * dy1) + sqrt (dx2 * dx2 + dy2 * dy2)) * prec ();
area_type p1 = dx1 * dy2;
area_type p2 = dy1 * dx2;
if (p1 <= p2 - pa) {
return std::make_pair (p1 - p2, -1);
} else if (p1 < p2 + pa) {
return std::make_pair (0, 0);
} else {
return std::make_pair (p1 - p2, 1);
}
}
static bool equal (double c1, double c2)

View File

@ -582,6 +582,15 @@ int vprod_sign (const db::vector<C> &p, const db::vector<C> &q)
return db::coord_traits<C>::vprod_sign (p.x (), p.y (), q.x (), q.y (), 0, 0);
}
/**
* @brief Convenience wrappers for coord_traits functions: vector product with sign
*/
template <class C>
std::pair<typename db::coord_traits<C>::area_type, int> vprod_with_sign (const db::vector<C> &p, const db::vector<C> &q)
{
return db::coord_traits<C>::vprod_with_sign (p.x (), p.y (), q.x (), q.y (), 0, 0);
}
/**
* @brief Convenience wrappers for coord_traits functions: scalar product: 0->p x 0->q
*/
@ -600,6 +609,15 @@ int sprod_sign (const db::vector<C> &p, const db::vector<C> &q)
return db::coord_traits<C>::sprod_sign (p.x (), p.y (), q.x (), q.y (), 0, 0);
}
/**
* @brief Convenience wrappers for coord_traits functions: scalar product with sign
*/
template <class C>
std::pair<typename db::coord_traits<C>::area_type, int> sprod_with_sign (const db::vector<C> &p, const db::vector<C> &q)
{
return db::coord_traits<C>::sprod_with_sign (p.x (), p.y (), q.x (), q.y (), 0, 0);
}
/**
* @brief A generic conversion operator from double vector to any type
*/

View File

@ -103,4 +103,72 @@ TEST(5)
EXPECT_EQ (p1.to_string (), "150,-150")
}
TEST(6)
{
EXPECT_EQ (db::sprod (db::Vector (0, 1000), db::Vector (1000, 1)), 1000);
EXPECT_EQ (db::sprod_with_sign (db::Vector (0, 1000), db::Vector (1000, 1)).first, 1000);
EXPECT_EQ (db::sprod (db::DVector (0, 1000), db::DVector (1000, 0.5)), 500);
EXPECT_EQ (db::sprod_with_sign (db::DVector (0, 1000), db::DVector (1000, 0.5)).first, 500);
EXPECT_EQ (db::vprod (db::Vector (2, 1000), db::Vector (1000, 1)), -999998);
EXPECT_EQ (db::vprod_with_sign (db::Vector (2, 1000), db::Vector (1000, 1)).first, -999998);
EXPECT_EQ (db::vprod (db::DVector (0.5, 1000), db::DVector (1000, 2)), -999999);
EXPECT_EQ (db::vprod_with_sign (db::DVector (0.5, 1000), db::DVector (1000, 2)).first, -999999);
EXPECT_EQ (db::sprod_sign (db::Vector (0, 1000000000), db::Vector (1000000000, 0)), 0);
EXPECT_EQ (db::sprod_sign (db::Vector (0, 1000000000), db::Vector (1000000000, 1)), 1);
EXPECT_EQ (db::sprod_sign (db::Vector (0, 1000000000), db::Vector (1000000000, -1)), -1);
EXPECT_EQ (db::sprod_sign (db::Vector (1000000000, 0), db::Vector (0, 1000000000)), 0);
EXPECT_EQ (db::sprod_sign (db::Vector (1000000000, 0), db::Vector (1, 1000000000)), 1);
EXPECT_EQ (db::sprod_sign (db::Vector (1000000000, 0), db::Vector (-1, 1000000000)), -1);
EXPECT_EQ (db::vprod_sign (db::Vector (0, 1000000000), db::Vector (0, 1000000000)), 0);
EXPECT_EQ (db::vprod_sign (db::Vector (0, 1000000000), db::Vector (1, 1000000000)), -1);
EXPECT_EQ (db::vprod_sign (db::Vector (0, 1000000000), db::Vector (-1, 1000000000)), 1);
EXPECT_EQ (db::vprod_sign (db::Vector (1000000000, 0), db::Vector (1000000000, 0)), 0);
EXPECT_EQ (db::vprod_sign (db::Vector (1000000000, 0), db::Vector (1000000000, 1)), 1);
EXPECT_EQ (db::vprod_sign (db::Vector (1000000000, 0), db::Vector (1000000000, -1)), -1);
EXPECT_EQ (db::sprod_with_sign (db::Vector (0, 1000000000), db::Vector (1000000000, 0)).second, 0);
EXPECT_EQ (db::sprod_with_sign (db::Vector (0, 1000000000), db::Vector (1000000000, 1)).second, 1);
EXPECT_EQ (db::sprod_with_sign (db::Vector (0, 1000000000), db::Vector (1000000000, -1)).second, -1);
EXPECT_EQ (db::sprod_with_sign (db::Vector (1000000000, 0), db::Vector (0, 1000000000)).second, 0);
EXPECT_EQ (db::sprod_with_sign (db::Vector (1000000000, 0), db::Vector (1, 1000000000)).second, 1);
EXPECT_EQ (db::sprod_with_sign (db::Vector (1000000000, 0), db::Vector (-1, 1000000000)).second, -1);
EXPECT_EQ (db::vprod_with_sign (db::Vector (0, 1000000000), db::Vector (0, 1000000000)).second, 0);
EXPECT_EQ (db::vprod_with_sign (db::Vector (0, 1000000000), db::Vector (1, 1000000000)).second, -1);
EXPECT_EQ (db::vprod_with_sign (db::Vector (0, 1000000000), db::Vector (-1, 1000000000)).second, 1);
EXPECT_EQ (db::vprod_with_sign (db::Vector (1000000000, 0), db::Vector (1000000000, 0)).second, 0);
EXPECT_EQ (db::vprod_with_sign (db::Vector (1000000000, 0), db::Vector (1000000000, 1)).second, 1);
EXPECT_EQ (db::vprod_with_sign (db::Vector (1000000000, 0), db::Vector (1000000000, -1)).second, -1);
EXPECT_EQ (db::sprod_sign (db::DVector (0, 100000.0000), db::DVector (100000.0000, 0)), 0);
EXPECT_EQ (db::sprod_sign (db::DVector (0, 100000.0000), db::DVector (100000.0000, 1e-7)), 0);
EXPECT_EQ (db::sprod_sign (db::DVector (0, 100000.0000), db::DVector (100000.0000, 0.0001)), 1);
EXPECT_EQ (db::sprod_sign (db::DVector (0, 100000.0000), db::DVector (100000.0000, -1e-7)), 0);
EXPECT_EQ (db::sprod_sign (db::DVector (0, 100000.0000), db::DVector (100000.0000, -0.0001)), -1);
EXPECT_EQ (db::sprod_sign (db::DVector (100000.0000, 0), db::DVector (0, 100000.0000)), 0);
EXPECT_EQ (db::sprod_sign (db::DVector (100000.0000, 0), db::DVector (0.0001, 100000.0000)), 1);
EXPECT_EQ (db::sprod_sign (db::DVector (100000.0000, 0), db::DVector (-0.0001, 100000.0000)), -1);
EXPECT_EQ (db::vprod_sign (db::DVector (0, 100000.0000), db::DVector (0, 100000.0000)), 0);
EXPECT_EQ (db::vprod_sign (db::DVector (0, 100000.0000), db::DVector (0.0001, 100000.0000)), -1);
EXPECT_EQ (db::vprod_sign (db::DVector (0, 100000.0000), db::DVector (-0.0001, 100000.0000)), 1);
EXPECT_EQ (db::vprod_sign (db::DVector (100000.0000, 0), db::DVector (100000.0000, 0)), 0);
EXPECT_EQ (db::vprod_sign (db::DVector (100000.0000, 0), db::DVector (100000.0000, 0.0001)), 1);
EXPECT_EQ (db::vprod_sign (db::DVector (100000.0000, 0), db::DVector (100000.0000, -0.0001)), -1);
EXPECT_EQ (db::sprod_with_sign (db::DVector (0, 100000.0000), db::DVector (100000.0000, 0)).second, 0);
EXPECT_EQ (db::sprod_with_sign (db::DVector (0, 100000.0000), db::DVector (100000.0000, 1e-7)).second, 0);
EXPECT_EQ (db::sprod_with_sign (db::DVector (0, 100000.0000), db::DVector (100000.0000, 0.0001)).second, 1);
EXPECT_EQ (db::sprod_with_sign (db::DVector (0, 100000.0000), db::DVector (100000.0000, -1e-7)).second, 0);
EXPECT_EQ (db::sprod_with_sign (db::DVector (0, 100000.0000), db::DVector (100000.0000, -0.0001)).second, -1);
EXPECT_EQ (db::sprod_with_sign (db::DVector (100000.0000, 0), db::DVector (0, 100000.0000)).second, 0);
EXPECT_EQ (db::sprod_with_sign (db::DVector (100000.0000, 0), db::DVector (0.0001, 100000.0000)).second, 1);
EXPECT_EQ (db::sprod_with_sign (db::DVector (100000.0000, 0), db::DVector (-0.0001, 100000.0000)).second, -1);
EXPECT_EQ (db::vprod_with_sign (db::DVector (0, 100000.0000), db::DVector (0, 100000.0000)).second, 0);
EXPECT_EQ (db::vprod_with_sign (db::DVector (0, 100000.0000), db::DVector (0.0001, 100000.0000)).second, -1);
EXPECT_EQ (db::vprod_with_sign (db::DVector (0, 100000.0000), db::DVector (-0.0001, 100000.0000)).second, 1);
EXPECT_EQ (db::vprod_with_sign (db::DVector (100000.0000, 0), db::DVector (100000.0000, 0)).second, 0);
EXPECT_EQ (db::vprod_with_sign (db::DVector (100000.0000, 0), db::DVector (100000.0000, 0.0001)).second, 1);
EXPECT_EQ (db::vprod_with_sign (db::DVector (100000.0000, 0), db::DVector (100000.0000, -0.0001)).second, -1);
}