17 const Eigen::MatrixXi &Q, Eigen::MatrixXi &edge_index,
18 std::vector<std::vector<int>> &
adj,
19 std::vector<std::pair<int, int>> *pairs_of_edges,
20 std::vector<std::pair<int, int>> *pairs_of_quads,
21 Eigen::MatrixXi *quad_index)
23 assert(Q.cols() == 4);
24 typedef std::pair<int, int>
Edge;
26 std::vector<std::pair<int, int>> edges;
28 edge_index.resizeLike(Q);
33 pairs_of_quads->clear();
37 quad_index->setConstant(Q.rows(), Q.cols(), -1);
41 std::vector<std::pair<Edge, std::pair<int, int>>> accu;
42 accu.reserve((
size_t)Q.size());
43 for (
int f = 0; f < Q.rows(); ++f)
45 for (
int lv = 0; lv < 4; ++lv)
48 int y = Q(f, (lv + 1) % 4);
49 accu.emplace_back(
Edge(std::min(
x,
y), std::max(
x,
y)), std::make_pair(f, lv));
52 std::sort(accu.begin(), accu.end());
53 edges.reserve(accu.size() / 2);
54 edge_index.setConstant(-1);
58 Edge prev_side(-1, -1);
59 for (
const auto &kv : accu)
63 int f = kv.second.first;
64 int lv = kv.second.second;
69 edges.emplace_back(e);
77 (*quad_index)(f, lv) = prev_quad;
78 (*quad_index)(prev_quad, prev_lv) = f;
80 pairs_of_quads->emplace_back(prev_quad, f);
83 edge_index(f, lv) = id;
88 adj.resize(edges.size());
89 for (
int f = 0; f < Q.rows(); ++f)
91 for (
int lv = 0; lv < 4; ++lv)
93 int e1 = edge_index(f, lv);
94 int e2 = edge_index(f, (lv + 2) % 4);
98 adj[e1].push_back(e2);
99 adj[e2].push_back(e1);
106 std::swap(*pairs_of_edges, edges);
115 enum class PatternType
135 PatternType compute_pattern_type(
136 const Eigen::MatrixXd &
V,
const Eigen::MatrixXi &
F,
137 std::array<Eigen::VectorXi, 4> &border_vertices)
139 Eigen::Vector2d lower =
V.colwise().minCoeff().head<2>();
140 Eigen::Vector2d upper =
V.colwise().maxCoeff().head<2>();
142 std::array<Eigen::MatrixXd, 2> border;
155 for (
int d = 0; d < 2; ++d)
157 std::vector<std::pair<double, int>> low, upp;
158 for (
int i = 0; i <
V.rows(); ++i)
160 if (
V(i, 1 - d) == lower[1 - d])
162 low.emplace_back(
V(i, d), i);
164 if (
V(i, 1 - d) == upper[1 - d])
166 upp.emplace_back(
V(i, d), i);
169 std::sort(low.begin(), low.end());
170 std::sort(upp.begin(), upp.end());
171 border_vertices[d ? 3 : 0].resize(low.size());
172 for (
size_t i = 0; i < low.size(); ++i)
174 border_vertices[d ? 3 : 0][i] = low[i].second;
176 border_vertices[d ? 1 : 2].resize(upp.size());
177 for (
size_t i = 0; i < upp.size(); ++i)
179 border_vertices[d ? 1 : 2][i] = upp[i].second;
184 border_vertices[i].reverseInPlace();
188 for (
int d = 0; d < 2; ++d)
190 std::vector<double> low, upp;
191 for (
int i = 0; i <
V.rows(); ++i)
193 if (
V(i, d) == lower[d])
195 low.push_back(
V(i, 1 - d));
197 if (
V(i, d) == upper[d])
199 upp.push_back(
V(i, 1 - d));
202 std::sort(low.begin(), low.end());
203 std::sort(upp.begin(), upp.end());
204 if (low.size() != upp.size())
206 return PatternType::NOT_PERIODIC;
208 size_t n = low.size();
209 border[d].resize(n, 2);
210 for (
size_t i = 0; i < n; ++i)
212 border[d](i, 0) = low[i];
213 border[d](i, 1) = upp[i];
215 if (!border[d].col(0).isApprox(border[d].col(1)))
217 return PatternType::NOT_PERIODIC;
222 for (
int d = 0; d < 2; ++d)
224 if (!(border[d].col(0).array() - lower[d]).isApprox(upper[d] - border[d].col(0).array().reverse()))
226 return PatternType::NOT_SYMMETRIC;
231 if (border[0].size() == border[1].size())
233 if (!border[0].col(0).isApprox(border[1].col(0)))
235 logger().warn(
"Pattern boundaries have the same number of vertices, but their position differ slighly.");
236 Eigen::MatrixXd
X(border[0].size(), 2);
237 X.col(0) = border[0].col(0);
238 X.col(1) = border[1].col(0);
240 return PatternType::DOUBLE_PERIODIC;
244 return PatternType::SIMPLE_PERIODIC;
247 return PatternType::SIMPLE_PERIODIC;
252 Eigen::VectorXi vertex_degree(
const Eigen::MatrixXd &
V,
const Eigen::MatrixXi &
F)
254 Eigen::VectorXi count = Eigen::VectorXi::Zero(
V.rows());
256 for (
unsigned i = 0; i <
F.rows(); ++i)
258 for (
unsigned j = 0; j <
F.cols(); ++j)
261 if (
F(i, j) <
F(i, (j + 1) %
F.cols()))
264 count(
F(i, (j + 1) %
F.cols())) += 1;
277 const Eigen::MatrixXd &IV,
const Eigen::MatrixXi &IQ,
278 const Eigen::MatrixXd &PV,
const Eigen::MatrixXi &PF,
279 Eigen::MatrixXd &OV, Eigen::MatrixXi &OF, Eigen::VectorXi *SF,
283 std::array<Eigen::VectorXi, 4> border_vertices;
285 auto pattern = compute_pattern_type(PV, PF, border_vertices);
286 auto valence = vertex_degree(IV, IQ);
292 case PatternType::NOT_PERIODIC:
293 logger().error(
"Pattern not periodic");
295 case PatternType::NOT_SYMMETRIC:
296 logger().error(
"Pattern not symmetric");
298 case PatternType::SIMPLE_PERIODIC:
299 logger().error(
"Pattern simple periodic");
301 case PatternType::DOUBLE_PERIODIC:
305 logger().error(
"Unknown patter type");
310 auto lower = PV.colwise().minCoeff();
311 auto upper = PV.colwise().maxCoeff();
312 Eigen::MatrixXd PVN = (PV.rowwise() - lower).array().rowwise() / (upper - lower).cwiseMax(1e-5).array();
317 evalFunc = [&](
const Eigen::MatrixXd &uv, Eigen::MatrixXd &mapped,
int q) {
318 const auto &u = uv.col(0).array();
319 const auto &v = uv.col(1).array();
320 Eigen::RowVector3d a = IV.row(IQ(q, 0));
321 Eigen::RowVector3d b = IV.row(IQ(q, 1));
322 Eigen::RowVector3d c = IV.row(IQ(q, 2));
323 Eigen::RowVector3d d = IV.row(IQ(q, 3));
324 mapped = ((1 - u) * (1 - v)).matrix() * a + (u * (1 - v)).matrix() * b + (u * v).matrix() * c + ((1 - u) * v).matrix() * d;
329 Eigen::MatrixXd mapped;
332 Eigen::MatrixXd P0 = PVN.row(0);
333 evalFunc(P0, mapped, 0);
334 Eigen::MatrixXd
V(IQ.rows() * PV.rows(), mapped.cols());
335 Eigen::MatrixXi
F(IQ.rows() * PF.rows(), PF.cols());
338 SF->resize(
V.rows());
342 Eigen::VectorXi remap(
V.rows());
343 remap.setConstant(-1);
345 for (
int q = 0; q < IQ.rows(); ++q)
347 evalFunc(PVN, mapped, q);
348 int v0 = (int)PVN.rows() * q;
349 V.middleRows(v0, mapped.rows()) = mapped;
350 int f0 = (int)PF.rows() * q;
351 F.middleRows(f0, PF.rows()) = PF.array() + v0;
354 SF->segment(v0, mapped.rows()).setConstant(q);
358 auto min_max = [](
int x,
int y) {
359 return std::make_pair(std::min(
x,
y), std::max(
x,
y));
366 for (
int q1 = 0; q1 < IQ.rows(); ++q1)
368 const int v1 = (int)PVN.rows() * q1;
369 for (
int lv1 = 0; lv1 < 4; ++lv1)
371 const auto res = getAdjLocalEdge(q1, lv1);
372 const int q2 = std::get<0>(res);
373 const int v2 = (int)PVN.rows() * q2;
374 const int lv2 = std::get<1>(res);
375 const bool rev = std::get<2>(res);
378 Eigen::VectorXi side1 = border_vertices[lv1].array() + v1;
379 Eigen::VectorXi side2 = border_vertices[lv2].array() + v2;
382 side2.reverseInPlace();
384 for (
int ii = 0; ii < (int)side1.size(); ++ii)
386 const int x1 = side1[ii];
387 const int x2 = side2[ii];
392 remap(x2) = remap(x1);
397 for (
int v = 0; v <
V.rows(); ++v)
408 Eigen::MatrixXi edge_index;
409 std::vector<std::vector<int>>
adj;
410 std::vector<std::pair<int, int>> pairs_of_quads;
413 for (
const auto &qq : pairs_of_quads)
421 const int v1 = (int)PVN.rows() * q1;
422 const int v2 = (int)PVN.rows() * q2;
425 for (; lv1 < 4; ++lv1)
427 int x1 = IQ(q1, lv1);
428 int y1 = IQ(q1, (lv1 + 1) % 4);
429 for (lv2 = 0; lv2 < 4; ++lv2)
431 int x2 = IQ(q2, lv2);
432 int y2 = IQ(q2, (lv2 + 1) % 4);
433 if (min_max(x1, y1) == min_max(x2, y2))
435 int e = edge_index(q1, lv1);
438 assert(edge_index(q2, lv2) == e);
439 Eigen::VectorXi side1 = border_vertices[lv1].array() + v1;
440 Eigen::VectorXi side2 = border_vertices[lv2].array() + v2;
443 side1.reverseInPlace();
447 side2.reverseInPlace();
449 for (
int ii = 0; ii < (int)side1.size(); ++ii)
451 remap(side2[ii]) = side1[ii];
466 int num_remapped = remap.maxCoeff() + 1;
467 OV.resize(num_remapped,
V.cols());
468 for (
int v = 0; v <
V.rows(); ++v)
470 OV.row(remap(v)) =
V.row(v);
472 for (
int f = 0; f <
F.rows(); ++f)
474 for (
int lv = 0; lv <
F.cols(); ++lv)
488 Eigen::VectorXi tmp = *SF;
489 SF->resize(OV.rows());
491 for (
int v = 0; v <
V.rows(); ++v)
493 (*SF)(remap(v)) = tmp(v);
503 const Eigen::MatrixXd &IV,
const Eigen::MatrixXi &IF,
504 Eigen::MatrixXd &OV, Eigen::MatrixXi &OF)
506 Eigen::MatrixXd PV(9, 3);
507 Eigen::MatrixXi PF(4, 4);
530 assert(IV.cols() == 2 || IV.cols() == 3);
531 Eigen::RowVector3d bary;
537 OV.resize(IV.rows() + 1, IV.cols());
538 OV.topRows(IV.rows()) = IV;
539 OV.bottomRows(1) = bary.head(IV.cols());
540 int n = (int)IV.rows();
541 for (
int i = 0; i < IV.rows(); ++i)
543 OF.push_back({i, (i + 1) % n, n});
548 throw std::invalid_argument(
"Value t should be >= 0.0.");
552 throw std::invalid_argument(
"Value t >= 1.0 would create degenerate elements.");
557 OV.resize(2 * IV.rows(), IV.cols());
558 OV.topRows(IV.rows()) = IV;
559 OV.bottomRows(IV.rows()) = (t * IV).rowwise() + (1.0 - t) * bary.head(IV.cols());
560 int n = (int)IV.rows();
562 for (
int i = 0; i < IV.rows(); ++i)
564 OF.front().push_back(i + n);
565 OF.push_back({i, (i + 1) % n, (i + 1) % n + n, i + n});
571 throw std::invalid_argument(
"Non star-shaped input polygon.");
579 assert(IV.cols() == 2 || IV.cols() == 3);
580 assert(IV.rows() % 2 == 0);
581 Eigen::RowVector3d bary;
585 OV.resize(IV.rows() + 1, IV.cols());
586 OV.topRows(IV.rows()) = IV;
587 OV.bottomRows(1) = bary.head(IV.cols());
588 int n = (int)IV.rows();
589 for (
int i = 1; i < IV.rows(); i += 2)
591 OF.push_back({i, (i + 1) % n, (i + 2) % n, n});
596 throw std::invalid_argument(
"Non star-shaped input polygon.");
607 OF.front().resize(IV.rows());
608 std::iota(OF.front().begin(), OF.front().end(), 0);
619 assert(&M_in != &M_out);
622 M_out.facets.clear();
623 GEO::Attribute<GEO::index_t> c2e(M_in.facet_corners.attributes(),
"edge_id");
626 std::vector<int> edge_to_midpoint(M_in.edges.nb(), -1);
627 for (index_t f = 0; f < M_in.facets.nb(); ++f)
629 index_t nv = M_in.facets.nb_vertices(f);
637 assert(idx.vertex == (
int)M_in.facets.vertex(f, lv));
638 if (edge_to_midpoint[idx.edge] == -1)
642 assert(edge_to_midpoint[idx.edge] + 1 == (
int)M_out.vertices.nb());
650 assert(vf + 1 == M_out.vertices.nb());
653 for (index_t lv = 0; lv < M_in.facets.nb_vertices(f); ++lv)
658 int v12 = edge_to_midpoint[idx.edge];
660 assert(v12 != -1 && v01 != -1);
661 M_out.facets.create_quad(v1, v12, vf, v01);
667 for (index_t f = 0; f < M_in.facets.nb(); ++f)
669 if (M_in.facets.nb_vertices(f) <= 4)
673 GEO::vector<index_t> hole;
676 for (index_t lv = 0; lv < M_in.facets.nb_vertices(f); ++lv)
678 hole.push_back(idx.vertex);
679 hole.push_back(edge_to_midpoint[idx.edge]);
684 if (M_in.vertices.dimension() != 2)
690 index_t n = hole.size();
691 Eigen::MatrixXd P(n, 2),
V;
692 std::vector<std::vector<int>>
F;
693 for (index_t k = 0; k < n; ++k)
696 P.row(k) << pk[0], pk[1];
699 assert(
V.rows() >= n);
700 std::vector<int> remap(
V.rows() - n);
701 for (index_t k = n, lk = 0; k <
V.rows(); ++k, ++lk)
703 GEO::vec3 qk = GEO::vec3(
V(k, 0),
V(k, 1), 0);
706 for (
const auto &poly :
F)
708 GEO::vector<index_t> vertices;
713 vertices.push_back(hole[vk]);
717 vertices.push_back(remap[vk - n]);
720 M_out.facets.create_polygon(vertices);
733 assert(&M_in != &M_out);
736 M_out.facets.clear();
737 GEO::Attribute<GEO::index_t> c2e(M_in.facet_corners.attributes(),
"edge_id");
740 std::vector<int> edge_to_midpoint(M_in.edges.nb(), -1);
741 for (index_t f = 0; f < M_in.facets.nb(); ++f)
743 index_t nv = M_in.facets.nb_vertices(f);
751 assert(idx.vertex == (
int)M_in.facets.vertex(f, lv));
752 if (edge_to_midpoint[idx.edge] == -1)
756 assert(edge_to_midpoint[idx.edge] + 1 == (
int)M_out.vertices.nb());
761 std::array<index_t, 3> e2v;
762 for (index_t lv = 0; lv < M_in.facets.nb_vertices(f); ++lv)
767 int v12 = edge_to_midpoint[idx.edge];
770 assert(v12 != -1 && v01 != -1);
771 M_out.facets.create_triangle(v1, v12, v01);
773 M_out.facets.create_triangle(e2v[0], e2v[1], e2v[2]);
vector< list< int > > adj
Index get_index_from_face(const GEO::Mesh &M, const GEO::Attribute< GEO::index_t > &c2e, int f, int lv)
Index switch_edge(const GEO::Mesh &M, const GEO::Attribute< GEO::index_t > &c2e, Index idx)
Index next_around_face(const GEO::Mesh &M, const GEO::Attribute< GEO::index_t > &c2e, Index idx)
Index switch_vertex(const GEO::Mesh &M, Index idx)
void polar_split(const Eigen::MatrixXd &IV, Eigen::MatrixXd &OV, std::vector< std::vector< int > > &OF, double t=0.5)
Split a polygon using polar refinement.
void catmul_clark_split(const Eigen::MatrixXd &IV, Eigen::MatrixXd &OV, std::vector< std::vector< int > > &OF)
Split a polygon using polar refinement.
std::function< void(Eigen::MatrixXd &IV, Eigen::MatrixXd &OV, std::vector< std::vector< int > > &OF)> SplitFunction
void no_split(const Eigen::MatrixXd &IV, Eigen::MatrixXd &OV, std::vector< std::vector< int > > &OF)
Don't split polygons.
void refine_triangle_mesh(const GEO::Mesh &M_in, GEO::Mesh &M_out)
Refine a triangle mesh.
void refine_quad_mesh(const Eigen::MatrixXd &IV, const Eigen::MatrixXi &IF, Eigen::MatrixXd &OV, Eigen::MatrixXi &OF)
GEO::vec3 mesh_vertex(const GEO::Mesh &M, GEO::index_t v)
Retrieve a 3D vector with the position of a given vertex.
GEO::index_t mesh_create_vertex(GEO::Mesh &M, const GEO::vec3 &p)
void edge_adjacency_graph(const Eigen::MatrixXi &Q, Eigen::MatrixXi &edge_index, std::vector< std::vector< int > > &adj, std::vector< std::pair< int, int > > *pairs_of_edges=nullptr, std::vector< std::pair< int, int > > *pairs_of_quads=nullptr, Eigen::MatrixXi *quad_index=nullptr)
GEO::vec3 facet_barycenter(const GEO::Mesh &M, GEO::index_t f)
bool instantiate_pattern(const Eigen::MatrixXd &IV, const Eigen::MatrixXi &IF, const Eigen::MatrixXd &PV, const Eigen::MatrixXi &PF, Eigen::MatrixXd &OV, Eigen::MatrixXi &OF, Eigen::VectorXi *SF=nullptr, EvalParametersFunc evalFunc=nullptr, GetAdjacentLocalEdge getAdjLocalEdge=nullptr)
std::function< void(const Eigen::MatrixXd &, Eigen::MatrixXd &, int)> EvalParametersFunc
void refine_polygonal_mesh(const GEO::Mesh &M_in, GEO::Mesh &M_out, Polygons::SplitFunction split_func)
Refine a polygonal mesh.
std::function< std::tuple< int, int, bool >(int, int)> GetAdjacentLocalEdge
bool is_star_shaped(const Eigen::MatrixXd &IV, Eigen::RowVector3d &bary)
Determine whether a polygon is star-shaped or not.
spdlog::logger & logger()
Retrieves the current logger.