/* This file contains liberal use of asserts to assist code development and debugging.
* Standard Matplotlib builds disable asserts so they cause no performance reduction. To
* enable the asserts, you need to undefine the NDEBUG macro, which is achieved by
* passing ``b_ndebug=false`` to the Meson configuration.
*/
#include "../mplutils.h"
#include "_tri.h"
#include
#include
#include
TriEdge::TriEdge()
: tri(-1), edge(-1)
{}
TriEdge::TriEdge(int tri_, int edge_)
: tri(tri_), edge(edge_)
{}
bool TriEdge::operatorx;
*segs_ptr++ = point->y;
*codes_ptr++ = (point == line.cbegin() ? MOVETO : LINETO);
}
if (line.size() > 1) {
*(codes_ptr-1) = CLOSEPOLY;
}
}
py::list vertices_list(1);
vertices_list[0] = segs;
py::list codes_list(1);
codes_list[0] = codes;
return py::make_tuple(vertices_list, codes_list);
}
py::tuple TriContourGenerator::create_contour(const double& level)
{
clear_visited_flags(false);
Contour contour;
find_boundary_lines(contour, level);
find_interior_lines(contour, level, false);
return contour_line_to_segs_and_kinds(contour);
}
py::tuple TriContourGenerator::create_filled_contour(const double& lower_level,
const double& upper_level)
{
if (lower_level >= upper_level)
throw std::invalid_argument("filled contour levels must be increasing");
clear_visited_flags(true);
Contour contour;
find_boundary_lines_filled(contour, lower_level, upper_level);
find_interior_lines(contour, lower_level, false);
find_interior_lines(contour, upper_level, true);
return contour_to_segs_and_kinds(contour);
}
XY TriContourGenerator::edge_interp(int tri, int edge, const double& level)
{
return interp(_triangulation.get_triangle_point(tri, edge),
_triangulation.get_triangle_point(tri, (edge+1)%3),
level);
}
void TriContourGenerator::find_boundary_lines(Contour& contour,
const double& level)
{
// Traverse boundaries to find starting points for all contour lines that
// intersect the boundaries. For each starting point found, follow the
// line to its end before continuing.
const Triangulation& triang = _triangulation;
const Boundaries& boundaries = get_boundaries();
for (const auto & boundary : boundaries) {
bool startAbove, endAbove = false;
for (auto itb = boundary.cbegin(); itb != boundary.cend(); ++itb) {
if (itb == boundary.cbegin())
startAbove = get_z(triang.get_triangle_point(*itb)) >= level;
else
startAbove = endAbove;
endAbove = get_z(triang.get_triangle_point(itb->tri,
(itb->edge+1)%3)) >= level;
if (startAbove && !endAbove) {
// This boundary edge is the start point for a contour line,
// so follow the line.
ContourLine& contour_line = contour.emplace_back();
TriEdge tri_edge = *itb;
follow_interior(contour_line, tri_edge, true, level, false);
}
}
}
}
void TriContourGenerator::find_boundary_lines_filled(Contour& contour,
const double& lower_level,
const double& upper_level)
{
// Traverse boundaries to find starting points for all contour lines that
// intersect the boundaries. For each starting point found, follow the
// line to its end before continuing.
const Triangulation& triang = _triangulation;
const Boundaries& boundaries = get_boundaries();
for (Boundaries::size_type i = 0; i < boundaries.size(); ++i) {
const Boundary& boundary = boundaries[i];
for (Boundary::size_type j = 0; j < boundary.size(); ++j) {
if (!_boundaries_visited[i][j]) {
// z values of start and end of this boundary edge.
double z_start = get_z(triang.get_triangle_point(boundary[j]));
double z_end = get_z(triang.get_triangle_point(
boundary[j].tri, (boundary[j].edge+1)%3));
// Does this boundary edge's z increase through upper level
// and/or decrease through lower level?
bool incr_upper = (z_start < upper_level && z_end >= upper_level);
bool decr_lower = (z_start >= lower_level && z_end < lower_level);
if (decr_lower || incr_upper) {
// Start point for contour line, so follow it.
ContourLine& contour_line = contour.emplace_back();
TriEdge start_tri_edge = boundary[j];
TriEdge tri_edge = start_tri_edge;
// Traverse interior and boundaries until return to start.
bool on_upper = incr_upper;
do {
follow_interior(contour_line, tri_edge, true,
on_upper ? upper_level : lower_level, on_upper);
on_upper = follow_boundary(contour_line, tri_edge,
lower_level, upper_level, on_upper);
} while (tri_edge != start_tri_edge);
// Close polygon.
contour_line.push_back(contour_line.front());
}
}
}
}
// Add full boundaries that lie between the lower and upper levels. These
// are boundaries that have not been touched by an internal contour line
// which are stored in _boundaries_used.
for (Boundaries::size_type i = 0; i < boundaries.size(); ++i) {
if (!_boundaries_used[i]) {
const Boundary& boundary = boundaries[i];
double z = get_z(triang.get_triangle_point(boundary[0]));
if (z >= lower_level && z < upper_level) {
ContourLine& contour_line = contour.emplace_back();
for (auto edge : boundary) {
contour_line.push_back(triang.get_point_coords(
triang.get_triangle_point(edge)));
}
// Close polygon.
contour_line.push_back(contour_line.front());
}
}
}
}
void TriContourGenerator::find_interior_lines(Contour& contour,
const double& level,
bool on_upper)
{
const Triangulation& triang = _triangulation;
int ntri = triang.get_ntri();
for (int tri = 0; tri < ntri; ++tri) {
int visited_index = (on_upper ? tri+ntri : tri);
if (_interior_visited[visited_index] || triang.is_masked(tri))
continue; // Triangle has already been visited or is masked.
_interior_visited[visited_index] = true;
// Determine edge via which to leave this triangle.
int edge = get_exit_edge(tri, level, on_upper);
assert(edge >= -1 && edge < 3 && "Invalid exit edge");
if (edge == -1)
continue; // Contour does not pass through this triangle.
// Found start of new contour line loop.
ContourLine& contour_line = contour.emplace_back();
TriEdge tri_edge = triang.get_neighbor_edge(tri, edge);
follow_interior(contour_line, tri_edge, false, level, on_upper);
// Close line loop
contour_line.push_back(contour_line.front());
}
}
bool TriContourGenerator::follow_boundary(ContourLine& contour_line,
TriEdge& tri_edge,
const double& lower_level,
const double& upper_level,
bool on_upper)
{
const Triangulation& triang = _triangulation;
const Boundaries& boundaries = get_boundaries();
// Have TriEdge to start at, need equivalent boundary edge.
int boundary, edge;
triang.get_boundary_edge(tri_edge, boundary, edge);
_boundaries_used[boundary] = true;
bool stop = false;
bool first_edge = true;
double z_start, z_end = 0;
while (!stop)
{
assert(!_boundaries_visited[boundary][edge] && "Boundary already visited");
_boundaries_visited[boundary][edge] = true;
// z values of start and end points of boundary edge.
if (first_edge)
z_start = get_z(triang.get_triangle_point(tri_edge));
else
z_start = z_end;
z_end = get_z(triang.get_triangle_point(tri_edge.tri,
(tri_edge.edge+1)%3));
if (z_end > z_start) { // z increasing.
if (!(!on_upper && first_edge) &&
z_end >= lower_level && z_start < lower_level) {
stop = true;
on_upper = false;
} else if (z_end >= upper_level && z_start < upper_level) {
stop = true;
on_upper = true;
}
} else { // z decreasing.
if (!(on_upper && first_edge) &&
z_start >= upper_level && z_end < upper_level) {
stop = true;
on_upper = true;
} else if (z_start >= lower_level && z_end < lower_level) {
stop = true;
on_upper = false;
}
}
first_edge = false;
if (!stop) {
// Move to next boundary edge, adding point to contour line.
edge = (edge+1) % (int)boundaries[boundary].size();
tri_edge = boundaries[boundary][edge];
contour_line.push_back(triang.get_point_coords(
triang.get_triangle_point(tri_edge)));
}
}
return on_upper;
}
void TriContourGenerator::follow_interior(ContourLine& contour_line,
TriEdge& tri_edge,
bool end_on_boundary,
const double& level,
bool on_upper)
{
int& tri = tri_edge.tri;
int& edge = tri_edge.edge;
// Initial point.
contour_line.push_back(edge_interp(tri, edge, level));
while (true) {
int visited_index = tri;
if (on_upper)
visited_index += _triangulation.get_ntri();
// Check for end not on boundary.
if (!end_on_boundary && _interior_visited[visited_index])
break; // Reached start point, so return.
// Determine edge by which to leave this triangle.
edge = get_exit_edge(tri, level, on_upper);
assert(edge >= 0 && edge < 3 && "Invalid exit edge");
_interior_visited[visited_index] = true;
// Append new point to point set.
assert(edge >= 0 && edge < 3 && "Invalid triangle edge");
contour_line.push_back(edge_interp(tri, edge, level));
// Move to next triangle.
TriEdge next_tri_edge = _triangulation.get_neighbor_edge(tri,edge);
// Check if ending on a boundary.
if (end_on_boundary && next_tri_edge.tri == -1)
break;
tri_edge = next_tri_edge;
assert(tri_edge.tri != -1 && "Invalid triangle for internal loop");
}
}
const TriContourGenerator::Boundaries& TriContourGenerator::get_boundaries() const
{
return _triangulation.get_boundaries();
}
int TriContourGenerator::get_exit_edge(int tri,
const double& level,
bool on_upper) const
{
assert(tri >= 0 && tri < _triangulation.get_ntri() &&
"Triangle index out of bounds.");
unsigned int config =
(get_z(_triangulation.get_triangle_point(tri, 0)) >= level) |
(get_z(_triangulation.get_triangle_point(tri, 1)) >= level) = level) = 0 && point < _triangulation.get_npoints() &&
"Point index out of bounds.");
return _z.data()[point];
}
XY TriContourGenerator::interp(int point1,
int point2,
const double& level) const
{
assert(point1 >= 0 && point1 < _triangulation.get_npoints() &&
"Point index 1 out of bounds.");
assert(point2 >= 0 && point2 < _triangulation.get_npoints() &&
"Point index 2 out of bounds.");
assert(point1 != point2 && "Identical points");
double fraction = (get_z(point2) - level) / (get_z(point2) - get_z(point1));
return _triangulation.get_point_coords(point1)*fraction +
_triangulation.get_point_coords(point2)*(1.0 - fraction);
}
TrapezoidMapTriFinder::TrapezoidMapTriFinder(Triangulation& triangulation)
: _triangulation(triangulation),
_points(nullptr),
_tree(nullptr)
{}
TrapezoidMapTriFinder::~TrapezoidMapTriFinder()
{
clear();
}
bool
TrapezoidMapTriFinder::add_edge_to_tree(const Edge& edge)
{
std::vector trapezoids;
if (!find_trapezoids_intersecting_edge(edge, trapezoids))
return false;
assert(!trapezoids.empty() && "No trapezoids intersect edge");
const Point* p = edge.left;
const Point* q = edge.right;
Trapezoid* left_old = nullptr; // old trapezoid to the left.
Trapezoid* left_below = nullptr; // below trapezoid to the left.
Trapezoid* left_above = nullptr; // above trapezoid to the left.
// Iterate through trapezoids intersecting edge from left to right.
// Replace each old trapezoid with 2+ new trapezoids, and replace its
// corresponding nodes in the search tree with new nodes.
size_t ntraps = trapezoids.size();
for (size_t i = 0; i < ntraps; ++i) {
Trapezoid* old = trapezoids[i]; // old trapezoid to replace.
bool start_trap = (i == 0);
bool end_trap = (i == ntraps-1);
bool have_left = (start_trap && edge.left != old->left);
bool have_right = (end_trap && edge.right != old->right);
// Old trapezoid is replaced by up to 4 new trapezoids: left is to the
// left of the start point p, below/above are below/above the edge
// inserted, and right is to the right of the end point q.
Trapezoid* left = nullptr;
Trapezoid* below = nullptr;
Trapezoid* above = nullptr;
Trapezoid* right = nullptr;
// There are 4 different cases here depending on whether the old
// trapezoid in question is the start and/or end trapezoid of those
// that intersect the edge inserted. There is some code duplication
// here but it is much easier to understand this way rather than
// interleave the 4 different cases with many more if-statements.
if (start_trap && end_trap) {
// Edge intersects a single trapezoid.
if (have_left)
left = new Trapezoid(old->left, p, old->below, old->above);
below = new Trapezoid(p, q, old->below, edge);
above = new Trapezoid(p, q, edge, old->above);
if (have_right)
right = new Trapezoid(q, old->right, old->below, old->above);
// Set pairs of trapezoid neighbours.
if (have_left) {
left->set_lower_left(old->lower_left);
left->set_upper_left(old->upper_left);
left->set_lower_right(below);
left->set_upper_right(above);
}
else {
below->set_lower_left(old->lower_left);
above->set_upper_left(old->upper_left);
}
if (have_right) {
right->set_lower_right(old->lower_right);
right->set_upper_right(old->upper_right);
below->set_lower_right(right);
above->set_upper_right(right);
}
else {
below->set_lower_right(old->lower_right);
above->set_upper_right(old->upper_right);
}
}
else if (start_trap) {
// Old trapezoid is the first of 2+ trapezoids that the edge
// intersects.
if (have_left)
left = new Trapezoid(old->left, p, old->below, old->above);
below = new Trapezoid(p, old->right, old->below, edge);
above = new Trapezoid(p, old->right, edge, old->above);
// Set pairs of trapezoid neighbours.
if (have_left) {
left->set_lower_left(old->lower_left);
left->set_upper_left(old->upper_left);
left->set_lower_right(below);
left->set_upper_right(above);
}
else {
below->set_lower_left(old->lower_left);
above->set_upper_left(old->upper_left);
}
below->set_lower_right(old->lower_right);
above->set_upper_right(old->upper_right);
}
else if (end_trap) {
// Old trapezoid is the last of 2+ trapezoids that the edge
// intersects.
if (left_below->below == old->below) {
below = left_below;
below->right = q;
}
else
below = new Trapezoid(old->left, q, old->below, edge);
if (left_above->above == old->above) {
above = left_above;
above->right = q;
}
else
above = new Trapezoid(old->left, q, edge, old->above);
if (have_right)
right = new Trapezoid(q, old->right, old->below, old->above);
// Set pairs of trapezoid neighbours.
if (have_right) {
right->set_lower_right(old->lower_right);
right->set_upper_right(old->upper_right);
below->set_lower_right(right);
above->set_upper_right(right);
}
else {
below->set_lower_right(old->lower_right);
above->set_upper_right(old->upper_right);
}
// Connect to new trapezoids replacing prevOld.
if (below != left_below) {
below->set_upper_left(left_below);
if (old->lower_left == left_old)
below->set_lower_left(left_below);
else
below->set_lower_left(old->lower_left);
}
if (above != left_above) {
above->set_lower_left(left_above);
if (old->upper_left == left_old)
above->set_upper_left(left_above);
else
above->set_upper_left(old->upper_left);
}
}
else { // Middle trapezoid.
// Old trapezoid is neither the first nor last of the 3+ trapezoids
// that the edge intersects.
if (left_below->below == old->below) {
below = left_below;
below->right = old->right;
}
else
below = new Trapezoid(old->left, old->right, old->below, edge);
if (left_above->above == old->above) {
above = left_above;
above->right = old->right;
}
else
above = new Trapezoid(old->left, old->right, edge, old->above);
// Connect to new trapezoids replacing prevOld.
if (below != left_below) { // below is new.
below->set_upper_left(left_below);
if (old->lower_left == left_old)
below->set_lower_left(left_below);
else
below->set_lower_left(old->lower_left);
}
if (above != left_above) { // above is new.
above->set_lower_left(left_above);
if (old->upper_left == left_old)
above->set_upper_left(left_above);
else
above->set_upper_left(old->upper_left);
}
below->set_lower_right(old->lower_right);
above->set_upper_right(old->upper_right);
}
// Create new nodes to add to search tree. Below and above trapezoids
// may already have owning trapezoid nodes, in which case reuse them.
Node* new_top_node = new Node(
&edge,
below == left_below ? below->trapezoid_node : new Node(below),
above == left_above ? above->trapezoid_node : new Node(above));
if (have_right)
new_top_node = new Node(q, new_top_node, new Node(right));
if (have_left)
new_top_node = new Node(p, new Node(left), new_top_node);
// Insert new_top_node in correct position or positions in search tree.
Node* old_node = old->trapezoid_node;
if (old_node == _tree)
_tree = new_top_node;
else
old_node->replace_with(new_top_node);
// old_node has been removed from all of its parents and is no longer
// needed.
assert(old_node->has_no_parents() && "Node should have no parents");
delete old_node;
// Clearing up.
if (!end_trap) {
// Prepare for next loop.
left_old = old;
left_above = above;
left_below = below;
}
}
return true;
}
void
TrapezoidMapTriFinder::clear()
{
delete [] _points;
_points = nullptr;
_edges.clear();
delete _tree;
_tree = nullptr;
}
TrapezoidMapTriFinder::TriIndexArray
TrapezoidMapTriFinder::find_many(const CoordinateArray& x,
const CoordinateArray& y)
{
if (x.ndim() != 1 || y.ndim() != 1 || x.shape(0) != y.shape(0))
throw std::invalid_argument(
"x and y must be array-like with same shape");
// Create integer array to return.
auto n = x.shape(0);
TriIndexArray tri_indices_array(n);
auto tri_indices = tri_indices_array.mutable_unchecked();
auto x_data = x.data();
auto y_data = y.data();
// Fill returned array.
for (py::ssize_t i = 0; i < n; ++i)
tri_indices(i) = find_one(XY(x_data[i], y_data[i]));
return tri_indices_array;
}
int
TrapezoidMapTriFinder::find_one(const XY& xy)
{
const Node* node = _tree->search(xy);
assert(node != nullptr && "Search tree for point returned null node");
return node->get_tri();
}
bool
TrapezoidMapTriFinder::find_trapezoids_intersecting_edge(
const Edge& edge,
std::vector& trapezoids)
{
// This is the FollowSegment algorithm of de Berg et al, with some extra
// checks to deal with simple colinear (i.e. invalid) triangles.
trapezoids.clear();
Trapezoid* trapezoid = _tree->search(edge);
if (trapezoid == nullptr) {
assert(trapezoid != nullptr && "search(edge) returns null trapezoid");
return false;
}
trapezoids.push_back(trapezoid);
while (edge.right->is_right_of(*trapezoid->right)) {
int orient = edge.get_point_orientation(*trapezoid->right);
if (orient == 0) {
if (edge.point_below == trapezoid->right)
orient = +1;
else if (edge.point_above == trapezoid->right)
orient = -1;
else {
assert(0 && "Unable to deal with point on edge");
return false;
}
}
if (orient == -1)
trapezoid = trapezoid->lower_right;
else if (orient == +1)
trapezoid = trapezoid->upper_right;
if (trapezoid == nullptr) {
assert(0 && "Expected trapezoid neighbor");
return false;
}
trapezoids.push_back(trapezoid);
}
return true;
}
py::list
TrapezoidMapTriFinder::get_tree_stats()
{
NodeStats stats;
_tree->get_stats(0, stats);
py::list ret(7);
ret[0] = stats.node_count;
ret[1] = stats.unique_nodes.size(),
ret[2] = stats.trapezoid_count,
ret[3] = stats.unique_trapezoid_nodes.size(),
ret[4] = stats.max_parent_count,
ret[5] = stats.max_depth,
ret[6] = stats.sum_trapezoid_depth / stats.trapezoid_count;
return ret;
}
void
TrapezoidMapTriFinder::initialize()
{
clear();
const Triangulation& triang = _triangulation;
// Set up points array, which contains all of the points in the
// triangulation plus the 4 corners of the enclosing rectangle.
int npoints = triang.get_npoints();
_points = new Point[npoints + 4];
BoundingBox bbox;
for (int i = 0; i < npoints; ++i) {
XY xy = triang.get_point_coords(i);
// Avoid problems with -0.0 values different from 0.0
if (xy.x == -0.0)
xy.x = 0.0;
if (xy.y == -0.0)
xy.y = 0.0;
_points[i] = Point(xy);
bbox.add(xy);
}
// Last 4 points are corner points of enclosing rectangle. Enclosing
// rectangle made slightly larger in case corner points are already in the
// triangulation.
if (bbox.empty) {
bbox.add(XY(0.0, 0.0));
bbox.add(XY(1.0, 1.0));
}
else {
const double small = 0.1; // Any value > 0.0
bbox.expand( (bbox.upper - bbox.lower)*small );
}
_points[npoints ] = Point(bbox.lower); // SW point.
_points[npoints+1] = Point(bbox.upper.x, bbox.lower.y); // SE point.
_points[npoints+2] = Point(bbox.lower.x, bbox.upper.y); // NW point.
_points[npoints+3] = Point(bbox.upper); // NE point.
// Set up edges array.
// First the bottom and top edges of the enclosing rectangle.
_edges.emplace_back(&_points[npoints], &_points[npoints+1], -1, -1,
nullptr, nullptr);
_edges.emplace_back(&_points[npoints+2], &_points[npoints+3], -1, -1,
nullptr, nullptr);
// Add all edges in the triangulation that point to the right. Do not
// explicitly include edges that point to the left as the neighboring
// triangle will supply that, unless there is no such neighbor.
int ntri = triang.get_ntri();
for (int tri = 0; tri < ntri; ++tri) {
if (!triang.is_masked(tri)) {
for (int edge = 0; edge < 3; ++edge) {
Point* start = _points + triang.get_triangle_point(tri,edge);
Point* end = _points +
triang.get_triangle_point(tri,(edge+1)%3);
Point* other = _points +
triang.get_triangle_point(tri,(edge+2)%3);
TriEdge neighbor = triang.get_neighbor_edge(tri,edge);
if (end->is_right_of(*start)) {
const Point* neighbor_point_below = (neighbor.tri == -1) ?
nullptr : _points + triang.get_triangle_point(
neighbor.tri, (neighbor.edge+2)%3);
_edges.emplace_back(start, end, neighbor.tri, tri,
neighbor_point_below, other);
}
else if (neighbor.tri == -1) {
_edges.emplace_back(end, start, tri, -1, other, nullptr);
}
// Set triangle associated with start point if not already set.
if (start->tri == -1)
start->tri = tri;
}
}
}
// Initial trapezoid is enclosing rectangle.
_tree = new Node(new Trapezoid(&_points[npoints], &_points[npoints+1],
_edges[0], _edges[1]));
_tree->assert_valid(false);
// Randomly shuffle all edges other than first 2.
std::mt19937 rng(1234); // NOLINT(bugprone-random-generator-seed)
std::shuffle(_edges.begin()+2, _edges.end(), rng);
// Add edges, one at a time, to tree.
size_t nedges = _edges.size();
for (size_t index = 2; index < nedges; ++index) {
if (!add_edge_to_tree(_edges[index]))
throw std::runtime_error("Triangulation is invalid");
_tree->assert_valid(index == nedges-1);
}
}
void
TrapezoidMapTriFinder::print_tree()
{
assert(_tree != nullptr && "Null Node tree");
_tree->print();
}
TrapezoidMapTriFinder::Edge::Edge(const Point* left_,
const Point* right_,
int triangle_below_,
int triangle_above_,
const Point* point_below_,
const Point* point_above_)
: left(left_),
right(right_),
triangle_below(triangle_below_),
triangle_above(triangle_above_),
point_below(point_below_),
point_above(point_above_)
{
assert(left != nullptr && "Null left point");
assert(right != nullptr && "Null right point");
assert(right->is_right_of(*left) && "Incorrect point order");
assert(triangle_below >= -1 && "Invalid triangle below index");
assert(triangle_above >= -1 && "Invalid triangle above index");
}
int
TrapezoidMapTriFinder::Edge::get_point_orientation(const XY& xy) const
{
double cross_z = (xy - *left).cross_z(*right - *left);
return (cross_z > 0.0) ? +1 : ((cross_z < 0.0) ? -1 : 0);
}
double
TrapezoidMapTriFinder::Edge::get_slope() const
{
// Divide by zero is acceptable here.
XY diff = *right - *left;
return diff.y / diff.x;
}
double
TrapezoidMapTriFinder::Edge::get_y_at_x(const double& x) const
{
if (left->x == right->x) {
// If edge is vertical, return lowest y from left point.
assert(x == left->x && "x outside of edge");
return left->y;
}
else {
// Equation of line: left + lambda*(right - left) = xy.
// i.e. left.x + lambda(right.x - left.x) = x and similar for y.
double lambda = (x - left->x) / (right->x - left->x);
assert(lambda >= 0 && lambda y + lambda*(right->y - left->y);
}
}
bool
TrapezoidMapTriFinder::Edge::has_point(const Point* point) const
{
assert(point != nullptr && "Null point");
return (left == point || right == point);
}
bool
TrapezoidMapTriFinder::Edge::operator==(const Edge& other) const
{
return this == &other;
}
void
TrapezoidMapTriFinder::Edge::print_debug() const
{
std::cout remove_parent(this))
delete _union.xnode.right;
break;
case Type_YNode:
if (_union.ynode.below->remove_parent(this))
delete _union.ynode.below;
if (_union.ynode.above->remove_parent(this))
delete _union.ynode.above;
break;
case Type_TrapezoidNode:
delete _union.trapezoid;
break;
}
}
void
TrapezoidMapTriFinder::Node::add_parent(Node* parent)
{
assert(parent != nullptr && "Null parent");
assert(parent != this && "Cannot be parent of self");
assert(!has_parent(parent) && "Parent already in collection");
_parents.push_back(parent);
}
void
TrapezoidMapTriFinder::Node::assert_valid(bool tree_complete) const
{
#ifndef NDEBUG
// Check parents.
for (const auto & parent : _parents) {
assert(parent != this && "Cannot be parent of self");
assert(parent->has_child(this) && "Parent missing child");
}
// Check children, and recurse.
switch (_type) {
case Type_XNode:
assert(_union.xnode.left != nullptr && "Null left child");
assert(_union.xnode.left->has_parent(this) && "Incorrect parent");
assert(_union.xnode.right != nullptr && "Null right child");
assert(_union.xnode.right->has_parent(this) && "Incorrect parent");
_union.xnode.left->assert_valid(tree_complete);
_union.xnode.right->assert_valid(tree_complete);
break;
case Type_YNode:
assert(_union.ynode.below != nullptr && "Null below child");
assert(_union.ynode.below->has_parent(this) && "Incorrect parent");
assert(_union.ynode.above != nullptr && "Null above child");
assert(_union.ynode.above->has_parent(this) && "Incorrect parent");
_union.ynode.below->assert_valid(tree_complete);
_union.ynode.above->assert_valid(tree_complete);
break;
case Type_TrapezoidNode:
assert(_union.trapezoid != nullptr && "Null trapezoid");
assert(_union.trapezoid->trapezoid_node == this &&
"Incorrect trapezoid node");
_union.trapezoid->assert_valid(tree_complete);
break;
}
#endif
}
void
TrapezoidMapTriFinder::Node::get_stats(int depth,
NodeStats& stats) const
{
stats.node_count++;
if (depth > stats.max_depth)
stats.max_depth = depth;
bool new_node = stats.unique_nodes.insert(this).second;
if (new_node)
stats.max_parent_count = std::max(stats.max_parent_count,
static_cast(_parents.size()));
switch (_type) {
case Type_XNode:
_union.xnode.left->get_stats(depth+1, stats);
_union.xnode.right->get_stats(depth+1, stats);
break;
case Type_YNode:
_union.ynode.below->get_stats(depth+1, stats);
_union.ynode.above->get_stats(depth+1, stats);
break;
default: // Type_TrapezoidNode:
stats.unique_trapezoid_nodes.insert(this);
stats.trapezoid_count++;
stats.sum_trapezoid_depth += depth;
break;
}
}
int
TrapezoidMapTriFinder::Node::get_tri() const
{
switch (_type) {
case Type_XNode:
return _union.xnode.point->tri;
case Type_YNode:
if (_union.ynode.edge->triangle_above != -1)
return _union.ynode.edge->triangle_above;
else
return _union.ynode.edge->triangle_below;
default: // Type_TrapezoidNode:
assert(_union.trapezoid->below.triangle_above ==
_union.trapezoid->above.triangle_below &&
"Inconsistent triangle indices from trapezoid edges");
return _union.trapezoid->below.triangle_above;
}
}
bool
TrapezoidMapTriFinder::Node::has_child(const Node* child) const
{
assert(child != nullptr && "Null child node");
switch (_type) {
case Type_XNode:
return (_union.xnode.left == child || _union.xnode.right == child);
case Type_YNode:
return (_union.ynode.below == child ||
_union.ynode.above == child);
default: // Type_TrapezoidNode:
return false;
}
}
bool
TrapezoidMapTriFinder::Node::has_no_parents() const
{
return _parents.empty();
}
bool
TrapezoidMapTriFinder::Node::has_parent(const Node* parent) const
{
return (std::find(_parents.begin(), _parents.end(), parent) !=
_parents.end());
}
void
TrapezoidMapTriFinder::Node::print(int depth /* = 0 */) const
{
for (int i = 0; i < depth; ++i) std::cout search(xy);
}
default: // Type_TrapezoidNode:
return this;
}
}
TrapezoidMapTriFinder::Trapezoid*
TrapezoidMapTriFinder::Node::search(const Edge& edge)
{
switch (_type) {
case Type_XNode:
if (edge.left == _union.xnode.point)
return _union.xnode.right->search(edge);
else {
if (edge.left->is_right_of(*_union.xnode.point))
return _union.xnode.right->search(edge);
else
return _union.xnode.left->search(edge);
}
case Type_YNode:
if (edge.left == _union.ynode.edge->left) {
// Coinciding left edge points.
if (edge.get_slope() == _union.ynode.edge->get_slope()) {
if (_union.ynode.edge->triangle_above ==
edge.triangle_below)
return _union.ynode.above->search(edge);
else if (_union.ynode.edge->triangle_below ==
edge.triangle_above)
return _union.ynode.below->search(edge);
else {
assert(0 &&
"Invalid triangulation, common left points");
return nullptr;
}
}
if (edge.get_slope() > _union.ynode.edge->get_slope())
return _union.ynode.above->search(edge);
else
return _union.ynode.below->search(edge);
}
else if (edge.right == _union.ynode.edge->right) {
// Coinciding right edge points.
if (edge.get_slope() == _union.ynode.edge->get_slope()) {
if (_union.ynode.edge->triangle_above ==
edge.triangle_below)
return _union.ynode.above->search(edge);
else if (_union.ynode.edge->triangle_below ==
edge.triangle_above)
return _union.ynode.below->search(edge);
else {
assert(0 &&
"Invalid triangulation, common right points");
return nullptr;
}
}
if (edge.get_slope() > _union.ynode.edge->get_slope())
return _union.ynode.below->search(edge);
else
return _union.ynode.above->search(edge);
}
else {
int orient =
_union.ynode.edge->get_point_orientation(*edge.left);
if (orient == 0) {
// edge.left lies on _union.ynode.edge
if (_union.ynode.edge->point_above != nullptr &&
edge.has_point(_union.ynode.edge->point_above))
orient = -1;
else if (_union.ynode.edge->point_below != nullptr &&
edge.has_point(_union.ynode.edge->point_below))
orient = +1;
else {
assert(0 && "Invalid triangulation, point on edge");
return nullptr;
}
}
if (orient < 0)
return _union.ynode.above->search(edge);
else
return _union.ynode.below->search(edge);
}
default: // Type_TrapezoidNode:
return _union.trapezoid;
}
}
TrapezoidMapTriFinder::Trapezoid::Trapezoid(const Point* left_,
const Point* right_,
const Edge& below_,
const Edge& above_)
: left(left_), right(right_), below(below_), above(above_),
lower_left(nullptr), lower_right(nullptr),
upper_left(nullptr), upper_right(nullptr),
trapezoid_node(nullptr)
{
assert(left != nullptr && "Null left point");
assert(right != nullptr && "Null right point");
assert(right->is_right_of(*left) && "Incorrect point order");
}
void
TrapezoidMapTriFinder::Trapezoid::assert_valid(bool tree_complete) const
{
#ifndef NDEBUG
assert(left != nullptr && "Null left point");
assert(right != nullptr && "Null right point");
if (lower_left != nullptr) {
assert(lower_left->below == below &&
lower_left->lower_right == this &&
"Incorrect lower_left trapezoid");
assert(get_lower_left_point() == lower_left->get_lower_right_point() &&
"Incorrect lower left point");
}
if (lower_right != nullptr) {
assert(lower_right->below == below &&
lower_right->lower_left == this &&
"Incorrect lower_right trapezoid");
assert(get_lower_right_point() == lower_right->get_lower_left_point() &&
"Incorrect lower right point");
}
if (upper_left != nullptr) {
assert(upper_left->above == above &&
upper_left->upper_right == this &&
"Incorrect upper_left trapezoid");
assert(get_upper_left_point() == upper_left->get_upper_right_point() &&
"Incorrect upper left point");
}
if (upper_right != nullptr) {
assert(upper_right->above == above &&
upper_right->upper_left == this &&
"Incorrect upper_right trapezoid");
assert(get_upper_right_point() == upper_right->get_upper_left_point() &&
"Incorrect upper right point");
}
assert(trapezoid_node != nullptr && "Null trapezoid_node");
if (tree_complete) {
assert(below.triangle_above == above.triangle_below &&
"Inconsistent triangle indices from trapezoid edges");
}
#endif
}
XY
TrapezoidMapTriFinder::Trapezoid::get_lower_left_point() const
{
double x = left->x;
return XY(x, below.get_y_at_x(x));
}
XY
TrapezoidMapTriFinder::Trapezoid::get_lower_right_point() const
{
double x = right->x;
return XY(x, below.get_y_at_x(x));
}
XY
TrapezoidMapTriFinder::Trapezoid::get_upper_left_point() const
{
double x = left->x;
return XY(x, above.get_y_at_x(x));
}
XY
TrapezoidMapTriFinder::Trapezoid::get_upper_right_point() const
{
double x = right->x;
return XY(x, above.get_y_at_x(x));
}
void
TrapezoidMapTriFinder::Trapezoid::print_debug() const
{
std::cout