4#include <geogram/numerics/predicates.h>
5#include <igl/barycenter.h>
11 inline GEO::Sign point_is_in_half_plane(
12 const Eigen::RowVector2d &p,
const Eigen::RowVector2d &q1,
const Eigen::RowVector2d &q2)
14 return GEO::PCK::orient_2d(q1.data(), q2.data(), p.data());
17 inline bool intersect_segments(
18 const Eigen::RowVector2d &p1,
const Eigen::RowVector2d &p2,
19 const Eigen::RowVector2d &q1,
const Eigen::RowVector2d &q2,
20 Eigen::RowVector2d &result)
23 Eigen::RowVector2d Vp = p2 - p1;
24 Eigen::RowVector2d Vq = q2 - q1;
25 Eigen::RowVector2d pq = q1 - p1;
32 double delta = a * d - b * c;
38 double tp = (d * pq(0) - b * pq(1)) / delta;
40 result << (1.0 - tp) * p1(0) + tp * p2(0),
41 (1.0 - tp) * p1(1) + tp * p2(1);
51 const Eigen::RowVector2d &q1,
const Eigen::RowVector2d &q2, Eigen::MatrixXd &P_out)
54 assert(P_in.cols() == 2);
55 std::vector<Eigen::RowVector2d> result;
65 if (point_is_in_half_plane(P_in.row(0), q1, q2))
77 Eigen::RowVector2d prev_p = P_in.row(P_in.rows() - 1);
78 Sign prev_status = point_is_in_half_plane(prev_p, q1, q2);
80 for (
unsigned int i = 0; i < P_in.rows(); ++i)
82 Eigen::RowVector2d p = P_in.row(i);
83 Sign status = point_is_in_half_plane(p, q1, q2);
84 if (status != prev_status && status != ZERO && prev_status != ZERO)
86 Eigen::RowVector2d intersect;
87 if (intersect_segments(prev_p, p, q1, q2, intersect))
89 result.push_back(intersect);
108 prev_status = status;
111 P_out.resize((
int)result.size(), 2);
112 for (
size_t i = 0; i < result.size(); ++i)
114 P_out.row((
int)i) = result[i];
122 assert(IV.cols() == 2 || IV.cols() == 3);
125 Eigen::MatrixXd src, dst;
126 const auto &minV = IV.colwise().minCoeff().array();
127 const auto &maxV = IV.colwise().maxCoeff().array();
129 src.row(0) << minV(0), minV(1);
130 src.row(1) << maxV(0), minV(1);
131 src.row(2) << maxV(0), maxV(1);
132 src.row(3) << minV(0), maxV(1);
137 for (
unsigned int i = 0; i < IV.rows(); ++i)
139 unsigned int j = ((i + 1) % (
int)IV.rows());
140 const Eigen::RowVector2d &q1 = IV.row(i).head<2>();
141 const Eigen::RowVector2d &q2 = IV.row(j).head<2>();
160 Eigen::MatrixXi F(1, OV.rows());
161 for (
int i = 0; i < OV.rows(); ++i)
166 igl::barycenter(OV, F, BC);
169 int n = std::min(3, (
int)BC.cols());
171 bary.head(n) = BC.row(0).head(n);
180#ifdef POLYFEM_WITH_CLIPPER
181 using namespace ClipperLib;
186 co.AddPath(PolygonClipping::toClipperPolygon(IV), jtSquare, etClosedPolygon);
190 co.Execute(solution, cInt(eps * DOUBLE_TO_INT_SCALE_FACTOR));
191 assert(solution.size() == 1);
194 OV = PolygonClipping::fromClipperPolygon(solution.front());
196 throw "Compile with clipper!";
205 typedef Eigen::RowVector2d Vec2d;
207 double inline det(Vec2d u, Vec2d v)
209 return u[0] * v[1] - u[1] * v[0];
213 bool intersect_segment(
const Vec2d &a,
const Vec2d &b, Vec2d c, Vec2d d, Vec2d &ans)
215 const double eps = 1e-10;
216 double x = det(c - a, d - c);
217 double y = det(b - a, a - c);
218 double z = det(b - a, d - c);
219 if (std::abs(
z) < eps ||
x * z < 0 || x * z >
z *
z ||
y * z < 0 || y * z >
z *
z)
221 ans = c + (d - c) *
y /
z;
225 bool is_point_inside(
const Eigen::MatrixXd &poly,
const Vec2d &outside,
const Vec2d &query)
228 bool tmp, ans =
false;
229 for (
long i = 0; i < poly.rows(); ++i)
232 tmp = intersect_segment(query, outside, poly.row(i), poly.row((i + 1) % n), m);
244 assert(IV.cols() == 2);
245 Eigen::RowVector2d minV = IV.colwise().minCoeff().array();
246 Eigen::RowVector2d maxV = IV.colwise().maxCoeff().array();
247 Eigen::RowVector2d center = 0.5 * (maxV + minV);
248 Eigen::RowVector2d outside(2.0 * minV(0) - center(0), center(1));
249 inside.resize(Q.rows());
251 for (
long i = 0; i < Q.rows(); ++i)
253 inside[i] = is_point_inside(IV, outside, Q.row(i));
266 assert(poly.rows() >= 2);
269 auto length = [&](
int i) {
270 return (poly.row(i) - poly.row((i + 1) % n)).norm();
275 double max_length = 0;
276 double total_length = 0;
277 for (
int i = 0; i < n; ++i)
279 double len = length(i);
281 if (len > max_length)
290 S.resize(num_samples, poly.cols());
291 double spacing = total_length / num_samples;
292 double offset = length(i0) / 2.0;
293 double distance_to_next = length(i0);
294 for (
int s = 0, i = i0; s < num_samples; ++s)
296 double distance_to_sample = s * spacing + offset;
297 while (distance_to_sample > distance_to_next)
300 distance_to_next += length(i);
302 double t = (distance_to_next - distance_to_sample) / length(i);
303 S.row(s) = t * poly.row(i) + (1.0 - t) * poly.row((i + 1) % n);
void clip_polygon_by_half_plane(const Eigen::MatrixXd &P, const Eigen::RowVector2d &q1, const Eigen::RowVector2d &q2, Eigen::MatrixXd &result)
Clip a polygon by a half-plane.
void sample_polygon(const Eigen::MatrixXd &IV, int num_samples, Eigen::MatrixXd &S)
Sample points on a polygon, evenly spaced from each other.
void offset_polygon(const Eigen::MatrixXd &IV, Eigen::MatrixXd &OV, double eps)
Compute offset polygon.
void compute_visibility_kernel(const Eigen::MatrixXd &IV, Eigen::MatrixXd &OV)
Determine the kernel of the given polygon.
bool is_star_shaped(const Eigen::MatrixXd &IV, Eigen::RowVector3d &bary)
Determine whether a polygon is star-shaped or not.
int is_inside(const Eigen::MatrixXd &IV, const Eigen::MatrixXd &Q, std::vector< bool > &inside)
Compute whether points are inside a polygon.