Medusa  1.1
Coordinate Free Mehless Method implementation
PolygonShape.hpp
Go to the documentation of this file.
1 #ifndef MEDUSA_BITS_DOMAINS_POLYGONSHAPE_HPP_
2 #define MEDUSA_BITS_DOMAINS_POLYGONSHAPE_HPP_
3 
4 #include "PolygonShape_fwd.hpp"
8 
14 namespace mm {
15 
16 template <typename vec_t>
17 PolygonShape<vec_t>::PolygonShape(const std::vector<vec_t>& points) : points_(points) {
18  scalar_t area = 0;
19  int n = points_.size();
20  assert_msg(n >= 3, "At least three points are needed to form a polygon, got %d.", n);
21  for (int i = 0; i < n; ++i) {
22  int j = (i == n-1) ? 0 : (i+1); // next point
23  area += (points_[i][0] - points_[j][0]) * (points_[i][1] + points_[j][1]);
24  }
25  assert_msg(std::abs(area) > 1e-10, "Given polygon has no area.");
26  if (area < 0) {
27  std::reverse(points_.begin(), points_.end());
28  }
30 }
31 
32 template <typename vec_t>
34  base_t::setMargin(margin);
35  extendByMargin();
36 }
37 
38 template <typename vec_t>
40  int n = points_.size();
41  points_with_margin_ = points_;
42  for (int i = 0; i < n; ++i) {
43  int j = (i == n - 1) ? 0 : (i + 1); // next
44  int k = (i == 0) ? n - 1 : (i - 1); // prev
45  vec_t next_edge = (points_[j] - points_[i]).normalized();
46  vec_t prev_edge = (points_[i] - points_[k]).normalized();
47  vec_t normal = {next_edge[1] + prev_edge[1], -next_edge[0] - prev_edge[0]};
48  points_with_margin_[i] += margin_*normal.normalized();
49  }
50  points_with_margin_.push_back(points_with_margin_[0]);
51 }
52 
53 template <typename vec_t>
54 bool PolygonShape<vec_t>::contains(const vec_t& point) const {
55  // Winding number test loosely based on http://geomalgorithms.com/a03-_inclusion.html.
56  int n = points_.size(); // loop through all edges of the polygon
57  int wn = 0; // the winding number counter
58  auto& pts = points_with_margin_;
59  for (int i = 0; i < n; i++) { // edge from points[i] to points[i+1]
60  if (pts[i][1] <= point[1]) { // Start y <= point[1].
61  // An upward crossing and point left of edge.
62  if (pts[i+1][1] > point[1] && isLeft(pts[i], pts[i+1], point) > 0) {
63  ++wn; // Have a valid up intersect.
64  }
65  } else { // Start y > point[1] (no test needed).
66  // A downward crossing and point right of edge.
67  if (pts[i+1][1] <= point[1] && isLeft(pts[i], pts[i+1], point) < 0) {
68  --wn; // Have a valid down intersect.
69  }
70  }
71  }
72  return wn != 0;
73 }
74 
75 template <typename vec_t>
76 std::pair<vec_t, vec_t> PolygonShape<vec_t>::bbox() const {
77  scalar_t maxx = -INF, maxy = -INF, minx = INF, miny = INF;
78  int n = points_.size();
79  for (int i = 0; i < n; ++i) {
80  if (points_[i][0] > maxx) maxx = points_[i][0];
81  if (points_[i][0] < minx) minx = points_[i][0];
82  if (points_[i][1] > maxy) maxy = points_[i][1];
83  if (points_[i][1] < miny) miny = points_[i][1];
84  }
85  return {{minx, miny}, {maxx, maxy}};
86 }
87 
88 template <typename vec_t>
90  scalar_t step, int type) const {
92  int n = points_.size();
93  for (int i = 0; i < n; ++i) {
94  int j = (i == n-1) ? 0 : (i+1); // next
95  int k = (i == 0) ? n-1 : (i-1); // prev
96  vec_t edge = points_[j] - points_[i];
97  scalar_t len = edge.norm();
98 
99  // Add the vertex point first.
100  vec_t next_edge_n = edge / len;
101  vec_t prev_edge_n = (points_[i] - points_[k]).normalized();
102  vec_t p_normal = {next_edge_n[1] + prev_edge_n[1], -next_edge_n[0] - prev_edge_n[0]};
103  d.addBoundaryNode(points_[i], (type == 0) ? -i-1 : type, p_normal.normalized());
104 
105  // And then the interior points on the edge.
106  int count = iceil(len/step);
107  vec_t e_normal = {edge[1] / len, -edge[0] / len};
108  for (int c = 1; c < count; ++c) { // only internal points
109  d.addBoundaryNode(points_[i] + c / scalar_t(count) * edge,
110  (type == 0) ? -i-1 : type, e_normal);
111  }
112  }
113  return d;
114 }
115 
116 template <typename vec_t>
119  const std::function<scalar_t(vec_t)>& dr, int type) const {
121  int n = points_.size();
122  for (int i = 0; i < n; ++i) {
123  int j = (i == n-1) ? 0 : (i+1); // next
124  int k = (i == 0) ? n-1 : (i-1); // prev
125  vec_t edge = points_[j] - points_[i];
126  scalar_t len = edge.norm();
127 
128  // Add the vertex point first.
129  vec_t next_edge_n = edge / len;
130  vec_t prev_edge_n = (points_[i] - points_[k]).normalized();
131  vec_t p_normal = {next_edge_n[1] + prev_edge_n[1], -next_edge_n[0] - prev_edge_n[0]};
132  d.addBoundaryNode(points_[i], (type == 0) ? -i-1 : type, p_normal.normalized());
133 
134  auto points = discretization_helpers::discretizeLineWithDensity(points_[i], points_[j], dr);
135  vec_t e_normal = {points_[j][1] - points_[i][1], -points_[j][0] + points_[i][0]};
136  e_normal.normalize();
137  for (const auto& p : points) { // only internal points
138  d.addBoundaryNode(p, (type == 0) ? -i-1 : type, e_normal);
139  }
140  }
141  return d;
142 }
143 
144 } // namespace mm
145 
146 #endif // MEDUSA_BITS_DOMAINS_POLYGONSHAPE_HPP_
mm::PolygonShape::discretizeBoundaryWithDensity
DomainDiscretization< vec_t > discretizeBoundaryWithDensity(const std::function< scalar_t(vec_t)> &dr, int type) const override
Returns a discretization of the boundary of this shape with approximately uniform distance step betwe...
Definition: PolygonShape.hpp:118
mm
Root namespace for the whole library.
Definition: Gaussian.hpp:14
scalar_t
Scalar scalar_t
Type of the elements, alias of Scalar.
Definition: MatrixBaseAddons.hpp:16
mm::DomainDiscretization
Class representing domain discretization along with an associated shape.
Definition: DomainDiscretization_fwd.hpp:46
DomainDiscretization.hpp
mm::PolygonShape::discretizeBoundaryWithStep
DomainDiscretization< vec_t > discretizeBoundaryWithStep(scalar_t step, int type) const override
Returns a discretization of the boundary of this shape with approximately uniform distance step betwe...
Definition: PolygonShape.hpp:89
numutils.hpp
assert_msg
#define assert_msg(cond,...)
Assert with better error reporting.
Definition: assert.hpp:75
mm::PolygonShape::PolygonShape
PolygonShape(const std::vector< vec_t > &points)
Create polygon given a sequence of points.
Definition: PolygonShape.hpp:17
mm::iceil
int iceil(T x)
Ceils a floating point to an integer.
Definition: numutils.hpp:27
mm::INF
static const double INF
Infinite floating point value.
Definition: Config.hpp:45
mm::discretization_helpers::discretizeLineWithDensity
Range< vec_t > discretizeLineWithDensity(const vec_t &p, const vec_t &q, const func_t &delta_r)
Returns nodes lying on the line segment pq with approximate distances delta_r.
Definition: discretization_helpers.hpp:35
mm::PolygonShape::points_
std::vector< vec_t > points_
The points that define the polygon, stored in CCW order.
Definition: PolygonShape_fwd.hpp:31
discretization_helpers.hpp
mm::PolygonShape::extendByMargin
void extendByMargin()
Save a version of points extended by margin_ for later contains checks.
Definition: PolygonShape.hpp:39
mm::PolygonShape::contains
bool contains(const vec_t &point) const override
Winding number test for point in a polygon inclusion (paying respect to margin).
Definition: PolygonShape.hpp:54
mm::PolygonShape::setMargin
void setMargin(scalar_t margin) override
Computes the polygon extended by margin as well.
Definition: PolygonShape.hpp:33
mm::PolygonShape::bbox
std::pair< vec_t, vec_t > bbox() const override
Return the bounding box of the domain.
Definition: PolygonShape.hpp:76
PolygonShape_fwd.hpp
mm::DomainDiscretization::addBoundaryNode
int addBoundaryNode(const vec_t &point, int type, const vec_t &normal)
Adds a boundary node with given type and normal to the domain.
Definition: DomainDiscretization.hpp:56
mm::DomainShape::scalar_t
vec_t::Scalar scalar_t
Scalar data type used in computation.
Definition: DomainShape_fwd.hpp:55