Triangles: Better numerical stability

This commit is contained in:
Matthias Koefferlein 2023-08-15 19:57:11 +02:00
parent 7351810e55
commit bfccd24016
5 changed files with 371 additions and 46 deletions

View File

@ -119,7 +119,7 @@ Vertex::in_circle (const DPoint &point, const DPoint &center, double radius)
double dy = point.y () - center.y ();
double d2 = dx * dx + dy * dy;
double r2 = radius * radius;
double delta = std::max (1.0, fabs (d2 + r2)) * db::epsilon;
double delta = fabs (d2 + r2) * db::epsilon;
if (d2 < r2 - delta) {
return 1;
} else if (d2 < r2 + delta) {
@ -427,15 +427,15 @@ Triangle::circumcircle () const
db::DVector n1 = db::DVector (v1.y (), -v1.x ());
db::DVector n2 = db::DVector (v2.y (), -v2.x ());
double p1s = vertex(0)->sq_distance (db::DPoint ());
double p2s = vertex(1)->sq_distance (db::DPoint ());
double p3s = vertex(2)->sq_distance (db::DPoint ());
double p1s = v1.sq_length ();
double p2s = v2.sq_length ();
double s = db::vprod (v1, v2);
tl_assert (fabs (s) > db::epsilon);
db::DPoint center = db::DPoint () + (n2 * (p1s - p2s) - n1 * (p1s - p3s)) * (0.5 / s);
double radius = (*vertex (0) - center).length ();
db::DVector r = (n1 * p2s - n2 * p1s) * (0.5 / s);
db::DPoint center = *vertex (0) + r;
double radius = r.length ();
return std::make_pair (center, radius);
}

View File

@ -395,8 +395,8 @@ private:
bool m_is_segment;
// no copying
TriangleEdge &operator= (const TriangleEdge &);
TriangleEdge (const TriangleEdge &);
// @@@ TriangleEdge &operator= (const TriangleEdge &);
// @@@ TriangleEdge (const TriangleEdge &);
void set_left (Triangle *t);
void set_right (Triangle *t);

View File

@ -314,19 +314,19 @@ Triangles::find_points_around (db::Vertex *vertex, double radius)
}
db::Vertex *
Triangles::insert_point (const db::DPoint &point, std::vector<db::Triangle *> *new_triangles)
Triangles::insert_point (const db::DPoint &point, std::list<tl::weak_ptr<db::Triangle> > *new_triangles)
{
return insert (create_vertex (point), new_triangles);
}
db::Vertex *
Triangles::insert_point (db::DCoord x, db::DCoord y, std::vector<db::Triangle *> *new_triangles)
Triangles::insert_point (db::DCoord x, db::DCoord y, std::list<tl::weak_ptr<db::Triangle> > *new_triangles)
{
return insert (create_vertex (x, y), new_triangles);
}
db::Vertex *
Triangles::insert (db::Vertex *vertex, std::vector<db::Triangle *> *new_triangles)
Triangles::insert (db::Vertex *vertex, std::list<tl::weak_ptr<db::Triangle> > *new_triangles)
{
std::vector<db::Triangle *> tris = find_triangle_for_point (*vertex);
@ -452,7 +452,7 @@ Triangles::find_closest_edge (const db::DPoint &p, db::Vertex *vstart, bool insi
}
void
Triangles::insert_new_vertex (db::Vertex *vertex, std::vector<db::Triangle *> *new_triangles_out)
Triangles::insert_new_vertex (db::Vertex *vertex, std::list<tl::weak_ptr<db::Triangle> > *new_triangles_out)
{
if (mp_triangles.empty ()) {
@ -547,7 +547,7 @@ Triangles::add_more_triangles (std::vector<db::Triangle *> &new_triangles,
}
void
Triangles::split_triangle (db::Triangle *t, db::Vertex *vertex, std::vector<db::Triangle *> *new_triangles_out)
Triangles::split_triangle (db::Triangle *t, db::Vertex *vertex, std::list<tl::weak_ptr<db::Triangle> > *new_triangles_out)
{
t->unlink ();
@ -577,7 +577,7 @@ Triangles::split_triangle (db::Triangle *t, db::Vertex *vertex, std::vector<db::
}
void
Triangles::split_triangles_on_edge (const std::vector<db::Triangle *> &tris, db::Vertex *vertex, db::TriangleEdge *split_edge, std::vector<db::Triangle *> *new_triangles_out)
Triangles::split_triangles_on_edge (const std::vector<db::Triangle *> &tris, db::Vertex *vertex, db::TriangleEdge *split_edge, std::list<tl::weak_ptr<db::Triangle> > *new_triangles_out)
{
TriangleEdge *s1 = create_edge (split_edge->v1 (), vertex);
TriangleEdge *s2 = create_edge (split_edge->v2 (), vertex);
@ -654,7 +654,7 @@ Triangles::find_inside_circle (const db::DPoint &center, double radius) const
}
void
Triangles::remove (db::Vertex *vertex, std::vector<db::Triangle *> *new_triangles)
Triangles::remove (db::Vertex *vertex, std::list<tl::weak_ptr<db::Triangle> > *new_triangles)
{
if (vertex->begin_edges () == vertex->end_edges ()) {
// removing an orphan vertex -> ignore
@ -666,7 +666,7 @@ Triangles::remove (db::Vertex *vertex, std::vector<db::Triangle *> *new_triangle
}
void
Triangles::remove_outside_vertex (db::Vertex *vertex, std::vector<db::Triangle *> *new_triangles_out)
Triangles::remove_outside_vertex (db::Vertex *vertex, std::list<tl::weak_ptr<db::Triangle> > *new_triangles_out)
{
auto to_remove = vertex->triangles ();
@ -689,7 +689,7 @@ Triangles::remove_outside_vertex (db::Vertex *vertex, std::vector<db::Triangle *
}
void
Triangles::remove_inside_vertex (db::Vertex *vertex, std::vector<db::Triangle *> *new_triangles_out)
Triangles::remove_inside_vertex (db::Vertex *vertex, std::list<tl::weak_ptr<db::Triangle> > *new_triangles_out)
{
std::set<db::Triangle *, TriangleLessFunc> triangles_to_fix;
@ -789,7 +789,7 @@ Triangles::remove_inside_vertex (db::Vertex *vertex, std::vector<db::Triangle *>
}
int
Triangles::fix_triangles (const std::vector<db::Triangle *> &tris, const std::vector<db::TriangleEdge *> &fixed_edges, std::vector<db::Triangle *> *new_triangles)
Triangles::fix_triangles (const std::vector<db::Triangle *> &tris, const std::vector<db::TriangleEdge *> &fixed_edges, std::list<tl::weak_ptr<db::Triangle> > *new_triangles)
{
int flips = 0;
@ -1387,4 +1387,172 @@ Triangles::create_constrained_delaunay (const db::Region &region, double dbu)
constrain (edge_contours);
}
static bool is_skinny (db::Triangle *tri, const Triangles::TriangulateParameters &param)
{
if (param.b < db::epsilon) {
return false;
} else {
auto cr = tri->circumcircle ();
double lmin = tri->edge (0)->d ().sq_length ();
for (int i = 1; i < 3; ++i) {
lmin = std::min (lmin, tri->edge (i)->d ().sq_length ());
}
double delta = (fabs (lmin / cr.second) + fabs (param.b)) * db::epsilon;
return lmin / cr.second < param.b - delta;
}
}
static bool is_invalid (db::Triangle *tri, const Triangles::TriangulateParameters &param)
{
if (is_skinny (tri, param)) {
return true;
}
double amax = param.max_area;
if (param.max_area_border > db::epsilon) {
bool at_border = false;
for (int i = 0; i < 3; ++i) {
if (tri->edge (i)->is_segment ()) {
at_border = true;
}
}
if (at_border) {
amax = param.max_area_border;
}
}
if (amax > db::epsilon) {
double a = tri->area ();
double delta = (fabs (a) + fabs (amax)) * db::epsilon;
return tri->area () > amax + delta;
} else {
return false;
}
}
static unsigned int num_segments (db::Triangle *tri)
{
unsigned int n = 0;
for (int i = 0; i < 3; ++i) {
if (tri->edge (i)->is_segment ()) {
++n;
}
}
return n;
}
void
Triangles::triangulate (const db::Region &region, const TriangulateParameters &parameters, double dbu)
{
create_constrained_delaunay (region, dbu);
if (parameters.b < db::epsilon && parameters.max_area < db::epsilon && parameters.max_area_border < db::epsilon) {
// no refinement requested - we're done.
remove_outside_triangles ();
return;
}
unsigned int nloop = 0;
std::list<tl::weak_ptr<db::Triangle> > new_triangles;
for (auto t = mp_triangles.begin (); t != mp_triangles.end (); ++t) {
new_triangles.push_back (t.operator-> ());
}
// @@@ TODO: break if iteration gets stuck
while (nloop < 20) { // @@@
++nloop;
tl::info << "Iteration " << nloop << " ..";
std::list<tl::weak_ptr<db::Triangle> > to_consider;
for (auto t = new_triangles.begin (); t != new_triangles.end (); ++t) {
if (t->get () && ! (*t)->is_outside () && is_invalid (t->get (), parameters)) {
to_consider.push_back (*t);
}
}
tl::info << "@@@ to_consider " << to_consider.size();
if (to_consider.empty ()) {
break;
}
new_triangles.clear ();
for (auto t = to_consider.begin (); t != to_consider.end (); ++t) {
if (! t->get ()) {
// triangle got removed during loop
continue;
}
auto cr = (*t)->circumcircle();
auto center = cr.first;
if ((*t)->contains (center) >= 0) {
// heuristics #1: never insert a point into a triangle with more than two segments
if (num_segments (t->get ()) <= 1) {
// @@@ print(f"Inserting in-triangle center {repr(center)} of {repr(tri)}")
insert_point (center, &new_triangles);
}
} else {
db::Vertex *vstart = 0;
for (int i = 0; i < 3; ++i) {
db::TriangleEdge *edge = (*t)->edge (i);
vstart = (*t)->opposite (edge);
if (edge->side_of (*vstart) * edge->side_of (center) < 0) {
break;
}
}
db::TriangleEdge *edge = find_closest_edge (center, vstart, true /*inside only*/);
tl_assert (edge != 0);
if (! edge->is_segment () || edge->side_of (*vstart) * edge->side_of (center) >= 0) {
// @@@ print(f"Inserting out-of-triangle center {repr(center)} of {repr(tri)}")
insert_point (center, &new_triangles);
} else {
double sr = edge->d ().length () * 0.5;
db::DPoint pnew = *edge->v1 () + edge->d () * 0.5;
// @@@ print(f"split edge {repr(edge)} at {repr(pnew)}")
db::Vertex *vnew = insert_point (pnew, &new_triangles);
auto vertexes_in_diametral_circle = find_points_around (vnew, sr);
std::vector<db::Vertex *> to_delete;
for (auto v = vertexes_in_diametral_circle.begin (); v != vertexes_in_diametral_circle.end (); ++v) {
bool has_segment = false;
for (auto e = (*v)->begin_edges (); e != (*v)->end_edges () && ! has_segment; ++e) {
has_segment = e->is_segment ();
}
if (! has_segment) {
to_delete.push_back (*v);
}
}
// @@@ print(f" -> found {len(to_delete)} vertexes to remove")
for (auto v = to_delete.begin (); v != to_delete.end (); ++v) {
remove (*v, &new_triangles);
}
}
}
}
// @@@ tris.dump_as_gdstxt("debug2.txt")
}
remove_outside_triangles ();
}
}

View File

@ -40,6 +40,30 @@ class Layout;
class DB_PUBLIC Triangles
{
public:
struct TriangulateParameters
{
TriangulateParameters ()
: b (1.0),
max_area (0.0),
max_area_border (0.0)
{ }
/**
* @brief Min. readius-to-shortest edge ratio
*/
double b;
/**
* @brief Max area or zero for "no constraint"
*/
double max_area;
/**
* @brief Max area for border triangles or zero for "use max_area"
*/
double max_area_border;
};
typedef tl::shared_collection<db::Triangle> triangles_type;
typedef triangles_type::const_iterator triangle_iterator;
@ -82,7 +106,13 @@ public:
*/
void clear ();
// exposed for testing purposes:
/**
* @brief Creates a refined Delaunay triangulation for the given region
*/
void triangulate (const db::Region &region, const TriangulateParameters &parameters, double dbu = 1.0);
// -- exposed for testing purposes --
/**
* @brief Checks the triangle graph for consistency
@ -113,7 +143,7 @@ public:
* If "new_triangles" is not null, it will receive the list of new triangles created during
* the remove step.
*/
db::Vertex *insert_point (const db::DPoint &point, std::vector<db::Triangle *> *new_triangles = 0);
db::Vertex *insert_point (const db::DPoint &point, std::list<tl::weak_ptr<db::Triangle> > *new_triangles = 0);
/**
* @brief Inserts a new vertex as the given point
@ -121,7 +151,7 @@ public:
* If "new_triangles" is not null, it will receive the list of new triangles created during
* the remove step.
*/
db::Vertex *insert_point (db::DCoord x, db::DCoord y, std::vector<db::Triangle *> *new_triangles = 0);
db::Vertex *insert_point (db::DCoord x, db::DCoord y, std::list<tl::weak_ptr<db::Triangle> > *new_triangles = 0);
/**
* @brief Removes the given vertex
@ -129,7 +159,7 @@ public:
* If "new_triangles" is not null, it will receive the list of new triangles created during
* the remove step.
*/
void remove (db::Vertex *vertex, std::vector<db::Triangle *> *new_triangles = 0);
void remove (db::Vertex *vertex, std::list<tl::weak_ptr<db::Triangle> > *new_triangles = 0);
/**
* @brief Flips the given edge
@ -180,6 +210,11 @@ public:
*/
void create_constrained_delaunay (const db::Region &region, double dbu = 1.0);
/**
* @brief Returns a value indicating whether the edge is "illegal" (violates the Delaunay criterion)
*/
static bool is_illegal_edge (db::TriangleEdge *edge);
// NOTE: these functions are SLOW and intended to test purposes only
std::vector<db::Vertex *> find_touching (const db::DBox &box) const;
std::vector<db::Vertex *> find_inside_circle (const db::DPoint &center, double radius) const;
@ -198,21 +233,20 @@ private:
db::Triangle *create_triangle (db::TriangleEdge *e1, db::TriangleEdge *e2, db::TriangleEdge *e3);
void remove (db::Triangle *tri);
void remove_outside_vertex (db::Vertex *vertex, std::vector<db::Triangle *> *new_triangles = 0);
void remove_inside_vertex (db::Vertex *vertex, std::vector<db::Triangle *> *new_triangles_out = 0);
void remove_outside_vertex (db::Vertex *vertex, std::list<tl::weak_ptr<db::Triangle> > *new_triangles = 0);
void remove_inside_vertex (db::Vertex *vertex, std::list<tl::weak_ptr<db::Triangle> > *new_triangles_out = 0);
std::vector<db::Triangle *> fill_concave_corners (const std::vector<TriangleEdge *> &edges);
int fix_triangles(const std::vector<db::Triangle *> &tris, const std::vector<db::TriangleEdge *> &fixed_edges, std::vector<db::Triangle *> *new_triangles);
static bool is_illegal_edge (db::TriangleEdge *edge);
int fix_triangles(const std::vector<db::Triangle *> &tris, const std::vector<db::TriangleEdge *> &fixed_edges, std::list<tl::weak_ptr<db::Triangle> > *new_triangles);
std::vector<db::Triangle *> find_triangle_for_point (const db::DPoint &point);
db::TriangleEdge *find_closest_edge (const db::DPoint &p, db::Vertex *vstart = 0, bool inside_only = false);
db::Vertex *insert (db::Vertex *vertex, std::vector<db::Triangle *> *new_triangles = 0);
void split_triangle (db::Triangle *t, db::Vertex *vertex, std::vector<db::Triangle *> *new_triangles_out);
void split_triangles_on_edge (const std::vector<db::Triangle *> &tris, db::Vertex *vertex, db::TriangleEdge *split_edge, std::vector<db::Triangle *> *new_triangles_out);
db::Vertex *insert (db::Vertex *vertex, std::list<tl::weak_ptr<db::Triangle> > *new_triangles = 0);
void split_triangle (db::Triangle *t, db::Vertex *vertex, std::list<tl::weak_ptr<db::Triangle> > *new_triangles_out);
void split_triangles_on_edge (const std::vector<db::Triangle *> &tris, db::Vertex *vertex, db::TriangleEdge *split_edge, std::list<tl::weak_ptr<db::Triangle> > *new_triangles_out);
void add_more_triangles (std::vector<Triangle *> &new_triangles,
db::TriangleEdge *incoming_edge,
db::Vertex *from_vertex, db::Vertex *to_vertex,
db::TriangleEdge *conn_edge);
void insert_new_vertex(db::Vertex *vertex, std::vector<db::Triangle *> *new_triangles_out);
void insert_new_vertex(db::Vertex *vertex, std::list<tl::weak_ptr<Triangle> > *new_triangles_out);
std::vector<db::TriangleEdge *> ensure_edge_inner (db::Vertex *from, db::Vertex *to);
void join_edges (std::vector<TriangleEdge *> &edges);
};

View File

@ -29,7 +29,7 @@
#include <cstdlib>
#include <cmath>
TEST(Triangle_basic)
TEST(Triangles_basic)
{
db::Triangles tris;
tris.init_box (db::DBox (1, 0, 5, 4));
@ -40,7 +40,7 @@ TEST(Triangle_basic)
EXPECT_EQ (tris.check (), true);
}
TEST(Triangle_flip)
TEST(Triangles_flip)
{
db::Triangles tris;
tris.init_box (db::DBox (0, 0, 1, 1));
@ -61,7 +61,7 @@ TEST(Triangle_flip)
EXPECT_EQ (tris.check (), true);
}
TEST(Triangle_insert)
TEST(Triangles_insert)
{
db::Triangles tris;
tris.init_box (db::DBox (0, 0, 1, 1));
@ -71,7 +71,7 @@ TEST(Triangle_insert)
EXPECT_EQ (tris.check (), true);
}
TEST(Triangle_split_segment)
TEST(Triangles_split_segment)
{
db::Triangles tris;
tris.init_box (db::DBox (0, 0, 1, 1));
@ -81,7 +81,7 @@ TEST(Triangle_split_segment)
EXPECT_EQ (tris.check(), true);
}
TEST(Triangle_insert_vertex_twice)
TEST(Triangles_insert_vertex_twice)
{
db::Triangles tris;
tris.init_box (db::DBox (0, 0, 1, 1));
@ -93,7 +93,7 @@ TEST(Triangle_insert_vertex_twice)
EXPECT_EQ (tris.check(), true);
}
TEST(Triangle_insert_vertex_convex)
TEST(Triangles_insert_vertex_convex)
{
db::Triangles tris;
tris.insert_point (0.2, 0.2);
@ -105,7 +105,7 @@ TEST(Triangle_insert_vertex_convex)
EXPECT_EQ (tris.check(), true);
}
TEST(Triangle_insert_vertex_convex2)
TEST(Triangles_insert_vertex_convex2)
{
db::Triangles tris;
tris.insert_point (0.25, 0.1);
@ -116,7 +116,7 @@ TEST(Triangle_insert_vertex_convex2)
EXPECT_EQ (tris.check(), true);
}
TEST(Triangle_insert_vertex_convex3)
TEST(Triangles_insert_vertex_convex3)
{
db::Triangles tris;
tris.insert_point (0.25, 0.5);
@ -127,7 +127,7 @@ TEST(Triangle_insert_vertex_convex3)
EXPECT_EQ (tris.check(), true);
}
TEST(Triangle_search_edges_crossing)
TEST(Triangles_search_edges_crossing)
{
db::Triangles tris;
db::Vertex *v1 = tris.insert_point (0.2, 0.2);
@ -147,6 +147,85 @@ TEST(Triangle_search_edges_crossing)
EXPECT_EQ (std::find (xedges.begin (), xedges.end (), s2) != xedges.end (), true);
}
TEST(Triangles_illegal_edge1)
{
db::Vertex v1 (0, 0);
db::Vertex v2 (1.6, 1.6);
db::Vertex v3 (1, 2);
db::Vertex v4 (2, 1);
{
db::TriangleEdge e1 (&v1, &v3);
db::TriangleEdge e2 (&v3, &v4);
db::TriangleEdge e3 (&v4, &v1);
db::Triangle t1 (&e1, &e2, &e3);
db::TriangleEdge ee1 (&v2, &v3);
db::TriangleEdge ee2 (&v4, &v2);
db::Triangle t2 (&ee1, &e2, &ee2);
EXPECT_EQ (db::Triangles::is_illegal_edge (&e2), true);
}
{
// flipped
db::TriangleEdge e1 (&v1, &v2);
db::TriangleEdge e2 (&v2, &v3);
db::TriangleEdge e3 (&v3, &v1);
db::Triangle t1 (&e1, &e2, &e3);
db::TriangleEdge ee1 (&v1, &v4);
db::TriangleEdge ee2 (&v4, &v2);
db::Triangle t2 (&ee1, &ee2, &e1);
EXPECT_EQ (db::Triangles::is_illegal_edge (&e2), false);
}
}
TEST(Triangles_illegal_edge2)
{
// numerical border case
db::Vertex v1 (773.94756216690905, 114.45875269431208);
db::Vertex v2 (773.29574734131643, 113.47402096138073);
db::Vertex v3 (773.10652961562653, 114.25497975904504);
db::Vertex v4 (774.08856345337881, 113.60495072750861);
{
db::TriangleEdge e1 (&v1, &v2);
db::TriangleEdge e2 (&v2, &v4);
db::TriangleEdge e3 (&v4, &v1);
db::Triangle t1 (&e1, &e2, &e3);
db::TriangleEdge ee1 (&v2, &v3);
db::TriangleEdge ee2 (&v3, &v4);
db::Triangle t2 (&ee1, &ee2, &e2);
EXPECT_EQ (db::Triangles::is_illegal_edge (&e2), false);
}
{
// flipped
db::TriangleEdge e1 (&v1, &v2);
db::TriangleEdge e2 (&v2, &v3);
db::TriangleEdge e3 (&v3, &v1);
db::Triangle t1 (&e1, &e2, &e3);
db::TriangleEdge ee1 (&v1, &v4);
db::TriangleEdge ee2 (&v4, &v2);
db::Triangle t2 (&ee1, &ee2, &e1);
EXPECT_EQ (db::Triangles::is_illegal_edge (&e1), true);
}
}
// Returns a random float number between 0.0 and 1.0
inline double flt_rand ()
{
@ -163,7 +242,7 @@ namespace {
};
}
TEST(Triangle_test_heavy_insert)
TEST(Triangles_test_heavy_insert)
{
tl::info << "Running test_heavy_insert " << tl::noendl;
@ -219,7 +298,7 @@ TEST(Triangle_test_heavy_insert)
tl::info << tl::endl << "done.";
}
TEST(Triangle_test_heavy_remove)
TEST(Triangles_test_heavy_remove)
{
tl::info << "Running test_heavy_remove " << tl::noendl;
@ -273,7 +352,7 @@ TEST(Triangle_test_heavy_remove)
tl::info << tl::endl << "done.";
}
TEST(Triangle_test_ensure_edge)
TEST(Triangles_test_ensure_edge)
{
srand (0);
@ -339,7 +418,7 @@ bool safe_inside (const db::DBox &b1, const db::DBox &b2)
(ct::less (b1.top (), b2.top ()) || ct::equal (b1.top (), b2.top ()));
}
TEST(Triangle_test_constrain)
TEST(Triangles_test_constrain)
{
srand (0);
@ -398,7 +477,7 @@ TEST(Triangle_test_constrain)
EXPECT_EQ (tl::to_string (area_in), "0.25");
}
TEST(Triangle_test_heavy_constrain)
TEST(Triangles_test_heavy_constrain)
{
tl::info << "Running test_heavy_constrain " << tl::noendl;
@ -468,7 +547,7 @@ TEST(Triangle_test_heavy_constrain)
tl::info << tl::endl << "done.";
}
TEST(Triangle_test_heavy_find_point_around)
TEST(Triangles_test_heavy_find_point_around)
{
tl::info << "Running Triangle_test_heavy_find_point_around " << tl::noendl;
@ -514,7 +593,7 @@ TEST(Triangle_test_heavy_find_point_around)
tl::info << tl::endl << "done.";
}
TEST(Triangle_test_create_constrained_delaunay)
TEST(Triangles_test_create_constrained_delaunay)
{
db::Region r;
r.insert (db::Box (0, 0, 1000, 1000));
@ -528,6 +607,8 @@ TEST(Triangle_test_create_constrained_delaunay)
tri.create_constrained_delaunay (r);
tri.remove_outside_triangles ();
EXPECT_EQ (tri.check (), true);
EXPECT_EQ (tri.to_string (),
"((1000, 0), (0, 0), (200, 200)), "
"((0, 1000), (200, 200), (0, 0)), "
@ -538,3 +619,45 @@ TEST(Triangle_test_create_constrained_delaunay)
"((0, 1000), (800, 800), (200, 800)), "
"((0, 1000), (200, 800), (200, 200))");
}
TEST(Triangles_test_triangulate)
{
db::Region r;
r.insert (db::Box (0, 0, 1000, 1000));
db::Region r2;
r2.insert (db::Box (200, 200, 800, 800));
r -= r2;
db::Triangles::TriangulateParameters param;
param.max_area = 0.1;
db::Triangles tri;
tri.triangulate (r, param);
tri.dump ("debug.gds");
EXPECT_EQ (tri.check (), true);
}
#if 0
assert tris.check(check_delaunay = False)
amax = 0.0
l2rmin = 2.0
for tri in tris.triangles:
if not tri.is_outside:
_, radius = tri.circumcircle()
lmin = min([math.sqrt(t.square(s.d())) for s in tri.edges()])
l2rmin = min(l2rmin, lmin / radius)
amax = max(amax, tri.area())
print(f"max. area = {'%.5g'%amax}")
print(f"l/R min = {'%.5g'%l2rmin}")
tris.dump_as_gdstxt("out.txt")
#endif