Delaunay refinement, tests etc.

This commit is contained in:
Matthias Koefferlein 2023-08-16 00:05:46 +02:00
parent 26f1219cc5
commit e07d802500
5 changed files with 128 additions and 64 deletions

View File

@ -494,4 +494,45 @@ Triangle::contains (const db::DPoint &point) const
return res;
}
double
Triangle::min_edge_length () const
{
double lmin = edge (0)->d ().length ();
for (int i = 1; i < 3; ++i) {
lmin = std::min (lmin, edge (i)->d ().length ());
}
return lmin;
}
double
Triangle::b () const
{
double lmin = min_edge_length ();
auto cr = circumcircle ();
return lmin / cr.second;
}
bool
Triangle::has_segment () const
{
for (int i = 0; i < 3; ++i) {
if (edge (i)->is_segment ()) {
return true;
}
}
return false;
}
unsigned int
Triangle::num_segments () const
{
unsigned int n = 0;
for (int i = 0; i < 3; ++i) {
if (edge (i)->is_segment ()) {
++n;
}
}
return n;
}
}

View File

@ -497,6 +497,26 @@ public:
return mp_e1.get () == e || mp_e2.get () == e || mp_e3.get () == e;
}
/**
* @brief Returns the minimum edge length
*/
double min_edge_length () const;
/**
* @brief Returns the min edge length to circumcircle radius ratio
*/
double b () const;
/**
* @brief Returns a value indicating whether the triangle borders to a segment
*/
bool has_segment () const;
/**
* @brief Returns the number of segments the triangle borders to
*/
unsigned int num_segments () const;
private:
bool m_is_outside;
tl::weak_ptr<TriangleEdge> mp_e1, mp_e2, mp_e3;

View File

@ -26,6 +26,7 @@
#include "dbWriter.h"
#include "tlStream.h"
#include "tlLog.h"
#include "tlTimer.h"
#include <set>
@ -1389,16 +1390,12 @@ Triangles::create_constrained_delaunay (const db::Region &region, double dbu)
static bool is_skinny (db::Triangle *tri, const Triangles::TriangulateParameters &param)
{
if (param.b < db::epsilon) {
if (param.min_b < db::epsilon) {
return false;
} else {
auto cr = tri->circumcircle ();
double lmin = tri->edge (0)->d ().length ();
for (int i = 1; i < 3; ++i) {
lmin = std::min (lmin, tri->edge (i)->d ().length ());
}
double delta = (fabs (lmin / cr.second) + fabs (param.b)) * db::epsilon;
return lmin / cr.second < param.b - delta;
double b = tri->b ();
double delta = (b + param.min_b) * db::epsilon;
return b < param.min_b - delta;
}
}
@ -1410,43 +1407,28 @@ static bool is_invalid (db::Triangle *tri, const Triangles::TriangulateParameter
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) {
if (tri->has_segment ()) {
amax = param.max_area_border;
}
}
if (amax > db::epsilon) {
double a = tri->area ();
double delta = (fabs (a) + fabs (amax)) * db::epsilon;
double delta = (a + 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)
{
tl::SelfTimer timer (tl::verbosity () > parameters.base_verbosity, "Triangles::triangulate");
create_constrained_delaunay (region, dbu);
if (parameters.b < db::epsilon && parameters.max_area < db::epsilon && parameters.max_area_border < db::epsilon) {
if (parameters.min_b < db::epsilon && parameters.max_area < db::epsilon && parameters.max_area_border < db::epsilon) {
// no refinement requested - we're done.
remove_outside_triangles ();
@ -1464,7 +1446,9 @@ Triangles::triangulate (const db::Region &region, const TriangulateParameters &p
while (nloop < parameters.max_iterations) { // @@@
++nloop;
tl::info << "Iteration " << nloop << " ..";
if (tl::verbosity () >= parameters.base_verbosity + 10) {
tl::info << "Iteration " << nloop << " ..";
}
std::list<tl::weak_ptr<db::Triangle> > to_consider;
for (auto t = new_triangles.begin (); t != new_triangles.end (); ++t) {
@ -1472,12 +1456,15 @@ Triangles::triangulate (const db::Region &region, const TriangulateParameters &p
to_consider.push_back (*t);
}
}
tl::info << "@@@ to_consider " << to_consider.size();
if (to_consider.empty ()) {
break;
}
if (tl::verbosity () >= parameters.base_verbosity + 10) {
tl::info << to_consider.size() << " triangles to consider";
}
new_triangles.clear ();
for (auto t = to_consider.begin (); t != to_consider.end (); ++t) {
@ -1493,8 +1480,10 @@ Triangles::triangulate (const db::Region &region, const TriangulateParameters &p
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)}")
if (t->get ()->num_segments () <= 1) {
if (tl::verbosity () >= parameters.base_verbosity + 20) {
tl::info << "Inserting in-triangle center " << center.to_string () << " of " << (*t)->to_string (true);
}
insert_point (center, &new_triangles);
}
@ -1514,7 +1503,9 @@ Triangles::triangulate (const db::Region &region, const TriangulateParameters &p
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)}")
if (tl::verbosity () >= parameters.base_verbosity + 20) {
tl::info << "Inserting out-of-triangle center " << center << " of " << (*t)->to_string (true);
}
insert_point (center, &new_triangles);
} else {
@ -1522,7 +1513,9 @@ Triangles::triangulate (const db::Region &region, const TriangulateParameters &p
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)}")
if (tl::verbosity () >= parameters.base_verbosity + 20) {
tl::info << "Split edge " << edge->to_string (true) << " at " << pnew.to_string ();
}
db::Vertex *vnew = insert_point (pnew, &new_triangles);
auto vertexes_in_diametral_circle = find_points_around (vnew, sr);
@ -1537,7 +1530,9 @@ Triangles::triangulate (const db::Region &region, const TriangulateParameters &p
}
}
// @@@ print(f" -> found {len(to_delete)} vertexes to remove")
if (tl::verbosity () >= parameters.base_verbosity + 20) {
tl::info << " -> found " << to_delete.size () << " vertexes to remove";
}
for (auto v = to_delete.begin (); v != to_delete.end (); ++v) {
remove (*v, &new_triangles);
}
@ -1548,10 +1543,11 @@ Triangles::triangulate (const db::Region &region, const TriangulateParameters &p
}
// @@@ dump ("debug2.gds"); // @@@
}
if (tl::verbosity () >= parameters.base_verbosity + 20) {
tl::info << "Finishing ..";
}
remove_outside_triangles ();
}

View File

@ -43,16 +43,17 @@ public:
struct TriangulateParameters
{
TriangulateParameters ()
: b (1.0),
: min_b (1.0),
max_area (0.0),
max_area_border (0.0),
max_iterations (std::numeric_limits<size_t>::max ())
max_iterations (std::numeric_limits<size_t>::max ()),
base_verbosity (30)
{ }
/**
* @brief Min. readius-to-shortest edge ratio
*/
double b;
double min_b;
/**
* @brief Max area or zero for "no constraint"
@ -68,6 +69,11 @@ public:
* @brief Max number of iterations
*/
size_t max_iterations;
/**
* @brief The verbosity level above which triangulation reports details
*/
int base_verbosity;
};
typedef tl::shared_collection<db::Triangle> triangles_type;
@ -114,6 +120,9 @@ public:
/**
* @brief Creates a refined Delaunay triangulation for the given region
*
* The database unit should be chosen in a way that target area values are "in the order of 1".
* The algorithm becomes numerically unstable area constraints below 1e-4.
*/
void triangulate (const db::Region &region, const TriangulateParameters &parameters, double dbu = 1.0);

View File

@ -222,7 +222,7 @@ TEST(Triangles_illegal_edge2)
db::Triangle t2 (&ee1, &ee2, &e1);
EXPECT_EQ (db::Triangles::is_illegal_edge (&e1), true);
EXPECT_EQ (db::Triangles::is_illegal_edge (&e1), false);
}
}
@ -631,35 +631,33 @@ TEST(Triangles_test_triangulate)
r -= r2;
db::Triangles::TriangulateParameters param;
param.b = 1.21;
param.min_b = 1.2;
param.max_area = 1.0;
param.max_area_border = 0.0;
db::Triangles tri;
tri.triangulate (r, param, 0.001);
tri.dump ("debug.gds");
EXPECT_EQ (tri.check (), true);
for (auto t = tri.begin (); t != tri.end (); ++t) {
EXPECT_EQ (t->area () <= param.max_area, true);
EXPECT_EQ (t->b () >= param.min_b, true);
}
EXPECT_EQ (tri.num_triangles () > 100 && tri.num_triangles () < 150, true);
tl::info << tri.num_triangles ();
param.min_b = 1.0;
param.max_area = 0.1;
tri.triangulate (r, param, 0.001);
EXPECT_EQ (tri.check (), true);
for (auto t = tri.begin (); t != tri.end (); ++t) {
EXPECT_EQ (t->area () <= param.max_area, true);
EXPECT_EQ (t->b () >= param.min_b, true);
}
EXPECT_EQ (tri.num_triangles () > 900 && tri.num_triangles () < 1000, 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