diff --git a/Data/data/polylines_2/Archimedean_spiral.svg b/Data/data/polylines_2/Archimedean_spiral.svg new file mode 100644 index 00000000000..c2481821d07 --- /dev/null +++ b/Data/data/polylines_2/Archimedean_spiral.svg @@ -0,0 +1,40 @@ + + + + + + diff --git a/Data/data/polylines_2/nano.svg b/Data/data/polylines_2/nano.svg new file mode 100644 index 00000000000..15635a16f74 --- /dev/null +++ b/Data/data/polylines_2/nano.svg @@ -0,0 +1,27 @@ + + + + + + diff --git a/Documentation/doc/Documentation/packages.txt b/Documentation/doc/Documentation/packages.txt index 9b1fa73107e..e8a1738d1ee 100644 --- a/Documentation/doc/Documentation/packages.txt +++ b/Documentation/doc/Documentation/packages.txt @@ -118,6 +118,7 @@ \package_listing{Surface_mesh_shortest_path} \package_listing{Surface_mesh_skeletonization} \package_listing{Surface_mesh_approximation} +\package_listing{Vector_graphics_on_surfaces} \package_listing{Ridges_3} \package_listing{Jet_fitting_3} \package_listing{Point_set_3} diff --git a/Installation/cmake/modules/CGAL_NanoSVG_support.cmake b/Installation/cmake/modules/CGAL_NanoSVG_support.cmake new file mode 100644 index 00000000000..b1db19f3152 --- /dev/null +++ b/Installation/cmake/modules/CGAL_NanoSVG_support.cmake @@ -0,0 +1,6 @@ +if(NanoSVG_FOUND AND NOT TARGET CGAL::NanoSVG_support) + add_library(CGAL::NanoSVG_support INTERFACE IMPORTED) + set_target_properties(CGAL::NanoSVG_support PROPERTIES + INTERFACE_COMPILE_DEFINITIONS "CGAL_NANOSVG_ENABLED") + target_link_libraries(CGAL::NanoSVG_support INTERFACE NanoSVG::nanosvg) +endif() diff --git a/Installation/include/CGAL/license/Vector_graphics_on_surfaces.h b/Installation/include/CGAL/license/Vector_graphics_on_surfaces.h new file mode 100644 index 00000000000..32ca4b440c5 --- /dev/null +++ b/Installation/include/CGAL/license/Vector_graphics_on_surfaces.h @@ -0,0 +1,54 @@ +// Copyright (c) 2016 GeometryFactory SARL (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org) +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Andreas Fabri +// +// Warning: this file is generated, see include/CGAL/license/README.md + +#ifndef CGAL_LICENSE_VECTOR_GRAPHICS_ON_SURFACES_H +#define CGAL_LICENSE_VECTOR_GRAPHICS_ON_SURFACES_H + +#include +#include + +#ifdef CGAL_VECTOR_GRAPHICS_ON_SURFACES_COMMERCIAL_LICENSE + +# if CGAL_VECTOR_GRAPHICS_ON_SURFACES_COMMERCIAL_LICENSE < CGAL_RELEASE_DATE + +# if defined(CGAL_LICENSE_WARNING) + + CGAL_pragma_warning("Your commercial license for CGAL does not cover " + "this release of the Vector_graphics_on_surfaces package.") +# endif + +# ifdef CGAL_LICENSE_ERROR +# error "Your commercial license for CGAL does not cover this release \ + of the Vector_graphics_on_surfaces package. \ + You get this error, as you defined CGAL_LICENSE_ERROR." +# endif // CGAL_LICENSE_ERROR + +# endif // CGAL_VECTOR_GRAPHICS_ON_SURFACES_COMMERCIAL_LICENSE < CGAL_RELEASE_DATE + +#else // no CGAL_VECTOR_GRAPHICS_ON_SURFACES_COMMERCIAL_LICENSE + +# if defined(CGAL_LICENSE_WARNING) + CGAL_pragma_warning("\nThe macro CGAL_VECTOR_GRAPHICS_ON_SURFACES_COMMERCIAL_LICENSE is not defined." + "\nYou use the CGAL Vector_graphics_on_surfaces package under " + "the terms of the GPLv3+.") +# endif // CGAL_LICENSE_WARNING + +# ifdef CGAL_LICENSE_ERROR +# error "The macro CGAL_VECTOR_GRAPHICS_ON_SURFACES_COMMERCIAL_LICENSE is not defined.\ + You use the CGAL Vector_graphics_on_surfaces package under the terms of \ + the GPLv3+. You get this error, as you defined CGAL_LICENSE_ERROR." +# endif // CGAL_LICENSE_ERROR + +#endif // no CGAL_VECTOR_GRAPHICS_ON_SURFACES_COMMERCIAL_LICENSE + +#endif // CGAL_LICENSE_VECTOR_GRAPHICS_ON_SURFACES_H diff --git a/Installation/include/CGAL/license/gpl_package_list.txt b/Installation/include/CGAL/license/gpl_package_list.txt index 75ff9d5c9ed..6d7bd6ea978 100644 --- a/Installation/include/CGAL/license/gpl_package_list.txt +++ b/Installation/include/CGAL/license/gpl_package_list.txt @@ -106,5 +106,6 @@ Triangulation_2 2D Triangulation Triangulation_3 3D Triangulations Triangulation dD Triangulations Triangulation_on_sphere_2 2D Triangulation on Sphere +Vector_graphics_on_surfaces Visibility_2 2D Visibility Computation Voronoi_diagram_2 2D Voronoi Diagram Adaptor diff --git a/Lab/demo/Lab/CMakeLists.txt b/Lab/demo/Lab/CMakeLists.txt index eb2a5dc56b6..d94c08e37bb 100644 --- a/Lab/demo/Lab/CMakeLists.txt +++ b/Lab/demo/Lab/CMakeLists.txt @@ -274,6 +274,11 @@ if(CGAL_Qt6_FOUND AND Qt6_FOUND) add_item(scene_edit_box_item Plugins/PCA/Scene_edit_box_item.cpp) + if (CGAL_Eigen3_support) + add_item(locally_shortest_path_item Plugins/Bsurf/Locally_shortest_path_item.cpp) + target_link_libraries(locally_shortest_path_item PUBLIC scene_surface_mesh_item scene_polylines_item CGAL::Eigen3_support) + endif() + add_item(scene_image_item Scene_image_item.cpp) add_item(scene_surface_mesh_item Scene_surface_mesh_item.cpp) diff --git a/Lab/demo/Lab/Plugins/Bsurf/CMakeLists.txt b/Lab/demo/Lab/Plugins/Bsurf/CMakeLists.txt new file mode 100644 index 00000000000..9c3c4f1dbab --- /dev/null +++ b/Lab/demo/Lab/Plugins/Bsurf/CMakeLists.txt @@ -0,0 +1,9 @@ +include(CGALlab_macros) + +if(TARGET CGAL::Eigen3_support) + cgal_lab_plugin(locally_shortest_path_plugin Locally_shortest_path_plugin) + target_link_libraries(locally_shortest_path_plugin PUBLIC locally_shortest_path_item + scene_surface_mesh_item + scene_polylines_item + CGAL::Eigen3_support) +endif() diff --git a/Lab/demo/Lab/Plugins/Bsurf/Locally_shortest_path_item.cpp b/Lab/demo/Lab/Plugins/Bsurf/Locally_shortest_path_item.cpp new file mode 100644 index 00000000000..224d4aaf838 --- /dev/null +++ b/Lab/demo/Lab/Plugins/Bsurf/Locally_shortest_path_item.cpp @@ -0,0 +1,869 @@ +#include "Locally_shortest_path_item.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include + +using namespace CGAL::Three; + +namespace VGoS = CGAL::Vector_graphics_on_surfaces; +namespace PMP = CGAL::Polygon_mesh_processing; + +typedef Viewer_interface Vi; +typedef Triangle_container Tc; +typedef Edge_container Ec; +typedef Locally_shortest_path_item_priv Priv; + + +typedef Scene_surface_mesh_item::Face_graph Mesh; +typedef CGAL::AABB_face_graph_triangle_primitive Primitive; +typedef CGAL::AABB_traits_3 Traits; +typedef CGAL::AABB_tree Tree; + +// TODO: update +struct Locally_shortest_path_item::vertex{ + int id; + double x; + double y; + double z; + + CGAL::Point_3 position()const + { + return CGAL::Point_3(x,y,z); + } + double operator[](int i) + { + switch(i) + { + case 0: + return x; + case 1: + return y; + case 2: + return z; + default: + return 0; + } + } + + template + void set(const P& p, int id_) + { + x=p.x(); + y=p.y(); + z=p.z(); + id=id_; + } +}; + +struct Locally_shortest_path_item_priv{ + typedef CGAL::Simple_cartesian Kernel; + enum Face_containers{ + Faces = 0, + S_Faces, + Spheres, + S_Spheres, + P_Spheres, + P_Faces, + Nbf + }; + + enum Line_containers{ + Edges = 0, + S_Edges, + P_Edges, + Nbe + }; + + enum HL_Primitive{ + VERTEX=0, + EDGE, + FACE, + NO_TYPE + }; + + Locally_shortest_path_item_priv(const CGAL::Three::Scene_interface* scene_interface, + const Scene_surface_mesh_item* sm_item, + Scene_polylines_item* polyline_item, + std::size_t nb_pts, + Locally_shortest_path_item* ebi) + { + const CGAL::qglviewer::Vec offset = static_cast(CGAL::QGLViewer::QGLViewerPool().first())->offset(); + ready_to_hl = true; + scene = scene_interface; + mesh_item = sm_item; + spath_item=polyline_item; + spath_item->polylines.resize(1); + const Mesh& mesh = *mesh_item->face_graph(); + aabb_tree = Tree(faces(mesh).first, + faces(mesh).second, + mesh); + item = ebi; + selection_on = false; + Scene_item::Bbox bb = scene->bbox(); + double x=(bb.xmin()+bb.xmax())/2; + double y=(bb.ymin()+bb.ymax())/2; + double z=(bb.zmin()+bb.zmax())/2; + center_ = CGAL::qglviewer::Vec(x,y,z); + relative_center_ = CGAL::qglviewer::Vec(0,0,0); + remodel_frame = new Scene_item::ManipulatedFrame(); + remodel_frame->setTranslationSensitivity(1.0); + frame = new Scene_item::ManipulatedFrame(); + frame->setPosition(center_+offset); + frame->setSpinningSensitivity(100.0); //forbid spinning + constraint.setRotationConstraintType(CGAL::qglviewer::AxisPlaneConstraint::AXIS); + constraint.setTranslationConstraintType(CGAL::qglviewer::AxisPlaneConstraint::FREE); + constraint.setRotationConstraintDirection(CGAL::qglviewer::Vec(.0,.0,.1)); + frame->setConstraint(&constraint); + //create the sphere model + vertex_spheres.resize(0); + normal_spheres.resize(0); + create_flat_sphere(1.0f, vertex_spheres, normal_spheres,10); + + + // TODO: change the default + if (nb_pts==2) + { + //shortest_path + vertices.resize(2); + vertices[0].set( mesh.point(*mesh.vertices().begin()), 0 ); + vertices[1].set( mesh.point(*std::next(mesh.vertices().begin())), 1 ); + locations.resize(2); + } + else if (nb_pts==4) + { + // bezier + vertices.resize(4); + vertices[0].set( mesh.point(*mesh.vertices().begin()), 0 ); + vertices[1].set( mesh.point(*std::next(mesh.vertices().begin())), 1 ); + vertices[2].set( mesh.point(*std::next(mesh.vertices().begin(),2)), 2 ); + vertices[3].set( mesh.point(*std::next(mesh.vertices().begin(),3)), 3 ); + locations.resize(4); + } + + + vertex_faces.resize(0); + normal_faces.resize(0); + + reset_selection(); + last_picked_id = -1; + last_picked_type = -1; + QPixmap pix(":/cgal/cursors/resources/rotate_around_cursor.png"); + rotate_cursor = QCursor(pix); + +#ifndef CGAL_BSURF_USE_DIJKSTRA_SP + VGoS::init_geodesic_dual_solver(geodesic_solver, mesh); +#endif + } + + ~Locally_shortest_path_item_priv(){ + delete frame; + delete remodel_frame; + } + + mutable std::vector vertex_edges; + mutable std::vector color_edges; + mutable std::vector vertex_spheres; + mutable std::vector normal_spheres; + mutable std::vector center_spheres; + mutable std::vector color_spheres; + mutable std::vector vertex_faces; + mutable std::vector normal_faces; + mutable std::vector color_faces; + mutable std::vector hl_vertex; + mutable std::vector hl_normal; + + bool ready_to_hl; + + CGAL::qglviewer::ManipulatedFrame* frame; + CGAL::qglviewer::ManipulatedFrame* remodel_frame; + CGAL::qglviewer::Vec rf_last_pos; + CGAL::qglviewer::LocalConstraint constraint; + CGAL::qglviewer::Vec center_; + CGAL::qglviewer::Vec relative_center_; + + mutable std::vector vertices; + mutable std::vector> locations; + std::vector selected_vertices; + + void reset_selection(); + bool selection_on; + void picking(int& type, int& id, Viewer_interface *viewer); + + void initializeBuffers(Viewer_interface *viewer)const; + + void computeElements() const; + void draw_picking(Viewer_interface*); + void update_points(const QVector3D &dir); + + const Scene_interface* scene; + const Scene_surface_mesh_item* mesh_item; + Scene_polylines_item* spath_item; + Locally_shortest_path_item* item; + Tree aabb_tree; + QPoint picked_pixel; + HL_Primitive hl_type; + int last_picked_id; + int last_picked_type; + QCursor rotate_cursor; + bool path_invalidated=true; +#ifndef CGAL_BSURF_USE_DIJKSTRA_SP + VGoS::Dual_geodesic_solver geodesic_solver; +#endif +}; + +Locally_shortest_path_item::Locally_shortest_path_item(const CGAL::Three::Scene_interface* scene_interface, + const Scene_surface_mesh_item *sm_item, + Scene_polylines_item* polyline_item, + std::size_t nb_pts) +{ + d = new Locally_shortest_path_item_priv(scene_interface, sm_item, polyline_item, nb_pts, this); + + are_buffers_filled = false; + + for(CGAL::QGLViewer* v : CGAL::QGLViewer::QGLViewerPool()) + { + v->setMouseTracking(true); + } + connect(Three::mainWindow(), SIGNAL(newViewerCreated(QObject*)), + this, SLOT(connectNewViewer(QObject*))); + + setTriangleContainer(Priv::P_Faces , new Tc(Vi::PROGRAM_NO_SELECTION, false)); + setTriangleContainer(Priv::Faces , new Tc(Vi::PROGRAM_WITH_LIGHT, false)); + setTriangleContainer(Priv::S_Faces , new Tc(Vi::PROGRAM_WITH_LIGHT, false)); + setTriangleContainer(Priv::Spheres , new Tc(Vi::PROGRAM_SPHERES, false)); + setTriangleContainer(Priv::S_Spheres, new Tc(Vi::PROGRAM_SPHERES, false)); + setTriangleContainer(Priv::P_Spheres, new Tc(Vi::PROGRAM_DARK_SPHERES, false)); + + + for(int i=Priv::Nbe-1; i>=0; --i) + { + setEdgeContainer(i, + new Ec(Three::mainViewer()->isOpenGL_4_3() + ? Vi::PROGRAM_SOLID_WIREFRAME + : Vi::PROGRAM_NO_SELECTION, + false)); + } + contextMenu(); +} +QString Locally_shortest_path_item::toolTip() const +{ + std::stringstream ss; + ss << "Face locations:
"; + ss << std::setprecision(17); + for (auto fl : d->locations) + { + ss << " - " << fl.first << " (" << fl.second[0] << "," << fl.second[1] << "," << fl.second[2] << ")
"; + } + return QString::fromStdString(ss.str()); +} + + +QMenu* Locally_shortest_path_item::contextMenu() +{ + // disable "Alpha slider" in menu + QMenu* resMenu = Scene_item::contextMenu(); + bool prop = property("menu_changed").toBool(); + if(!prop) + { + setProperty("menu_changed", true); + } + return resMenu; +} + +void Locally_shortest_path_item::drawSpheres(Viewer_interface *viewer, const QMatrix4x4 f_matrix ) const +{ + GLdouble d_mat[16]; + QMatrix4x4 mv_mat; + viewer->camera()->getModelViewMatrix(d_mat); + for (int i=0; i<16; ++i) + mv_mat.data()[i] = GLfloat(d_mat[i]); + mv_mat = mv_mat*f_matrix; + // TODO: must depend on the mesh (and zoom?) + double radius = 0.01 ; + + Tc* tc = getTriangleContainer(Priv::Spheres); + tc->setFrameMatrix(f_matrix); + tc->setMvMatrix(mv_mat); + tc->setClipping(false); + tc->getVao(viewer)->bind(); + tc->getVao(viewer)->program->setAttributeValue("radius",radius); + tc->getVao(viewer)->release(); + tc->setColor(QColor(Qt::red)); + tc->draw(viewer, true); +} + +void Locally_shortest_path_item::drawPath() const +{ + if (!d->path_invalidated) return; + typedef PMP::Edge_location Edge_location; + typedef PMP::Face_location Face_location; + + const Mesh& mesh = *d->mesh_item->face_graph(); + + if (d->vertices.size()==2) + { + std::vector edge_locations; + CGAL::Epick::Point_3 src_pt(d->vertices[0].x,d->vertices[0].y,d->vertices[0].z), + tgt_pt(d->vertices[1].x,d->vertices[1].y,d->vertices[1].z); + + //TODO store that in the vector vertices + d->locations[0] = PMP::locate_with_AABB_tree(src_pt, d->aabb_tree, mesh); + d->locations[1] = PMP::locate_with_AABB_tree(tgt_pt, d->aabb_tree, mesh); + + VGoS::locally_shortest_path(d->locations[0], d->locations[1], mesh, edge_locations, d->geodesic_solver); + d->spath_item->polylines.back().clear(); + d->spath_item->polylines.back().push_back(src_pt); + for (auto el : edge_locations) + d->spath_item->polylines.back().push_back(PMP::construct_point(el, mesh)); + d->spath_item->polylines.back().push_back(tgt_pt); + d->spath_item->setRenderingMode(Wireframe); + d->spath_item->invalidateOpenGLBuffers(); + } + else if (d->vertices.size()==4) + { + std::vector edge_locations; + CGAL::Epick::Point_3 c1_pt(d->vertices[0].x,d->vertices[0].y,d->vertices[0].z), + c2_pt(d->vertices[1].x,d->vertices[1].y,d->vertices[1].z), + c3_pt(d->vertices[2].x,d->vertices[2].y,d->vertices[2].z), + c4_pt(d->vertices[3].x,d->vertices[3].y,d->vertices[3].z); + + //TODO store that in the vector vertices + d->locations[0] = PMP::locate_with_AABB_tree(c1_pt, d->aabb_tree, mesh); + d->locations[1] = PMP::locate_with_AABB_tree(c2_pt, d->aabb_tree, mesh); + d->locations[2] = PMP::locate_with_AABB_tree(c3_pt, d->aabb_tree, mesh); + d->locations[3] = PMP::locate_with_AABB_tree(c4_pt, d->aabb_tree, mesh); + + VGoS::Bezier_segment control_points=CGAL::make_array(d->locations[0], + d->locations[1], + d->locations[2], + d->locations[3]); + + std::vector face_locations = + VGoS::recursive_de_Casteljau(mesh, control_points, 8, d->geodesic_solver); + + // TODO: we should connect points with geodesics and not segments + d->spath_item->polylines.back().clear(); + for (auto fl : face_locations) + d->spath_item->polylines.back().push_back(PMP::construct_point(fl, mesh)); + d->spath_item->setRenderingMode(Wireframe); + d->spath_item->invalidateOpenGLBuffers(); + } + d->path_invalidated=false; +} + +void Locally_shortest_path_item::draw(Viewer_interface *viewer) const +{ + if(!isInit(viewer)) + initGL(viewer); + if ( getBuffersFilled() && + ! getBuffersInit(viewer)) + { + initializeBuffers(viewer); + setBuffersInit(viewer, true); + } + if(!getBuffersFilled()) + { + computeElements(); + initializeBuffers(viewer); + } + QMatrix4x4 f_matrix; + for (int i=0; i<16; ++i){ + f_matrix.data()[i] = (float)d->frame->matrix()[i]; + } + + drawSpheres(viewer, f_matrix); + drawHl(viewer); + + drawPath(); +} + +void Locally_shortest_path_item::compute_bbox() const +{ + double xmin=d->vertices[0].x, xmax=xmin; + double ymin=d->vertices[0].y, ymax=ymin; + double zmin=d->vertices[0].z, zmax=zmin; + for(std::size_t i=0; ivertices.size(); ++i) + { + xmin=(std::min)(xmin, d->vertices[i].x); + ymin=(std::min)(ymin, d->vertices[i].y); + zmin=(std::min)(zmin, d->vertices[i].z); + xmax=(std::max)(xmax, d->vertices[i].x); + ymax=(std::max)(ymax, d->vertices[i].y); + zmax=(std::max)(zmax, d->vertices[i].z); + } + + setBbox(Scene_item::Bbox(xmin, ymin, zmin,xmax, ymax, zmax)); +} + +void Locally_shortest_path_item_priv::computeElements() const +{ + vertex_edges.clear(); + vertex_faces.clear(); + normal_faces.clear(); + center_spheres.clear(); + color_edges.clear(); + color_faces.clear(); + color_spheres.clear(); + + for (std::size_t i=0; isetMouseTracking(false); + + delete d; +} + +// Indicate if rendering mode is supported +bool Locally_shortest_path_item::supportsRenderingMode(RenderingMode m) const { + return m==FlatPlusEdges; +} + +Scene_item::ManipulatedFrame* +Locally_shortest_path_item::manipulatedFrame() +{ + return d->frame; +} + +void Locally_shortest_path_item::highlight(Viewer_interface *viewer) +{ + d->ready_to_hl = true; + viewer->makeCurrent(); + int type = -1, id = -1; + //pick + if(!d->selection_on) + { + d->picking(type, id, viewer); + d->last_picked_id = id; + d->last_picked_type = type; + } + //highlight + d->hl_normal.clear(); + d->hl_vertex.clear(); + if(type !=-1) + { + switch(d->last_picked_type) + { + case 0: + { + //compute + d->hl_vertex.push_back(d->vertices[d->last_picked_id].x); + d->hl_vertex.push_back(d->vertices[d->last_picked_id].y); + d->hl_vertex.push_back(d->vertices[d->last_picked_id].z); + //fill buffers + Tc* tc = getTriangleContainer(Priv::S_Spheres); + tc->reset_vbos(ALL); + tc->allocate( + Tc::Flat_vertices, + d->vertex_spheres.data(), + static_cast(d->vertex_spheres.size()*sizeof(float))); + + tc->allocate( + Tc::Flat_normals, + d->normal_spheres.data(), + static_cast(d->normal_spheres.size()*sizeof(float))); + + tc->allocate( + Tc::Facet_centers, + d->hl_vertex.data(), + static_cast(d->hl_vertex.size()*sizeof(float))); + + tc->setFlatDataSize(d->vertex_spheres.size()); + tc->setCenterSize(d->hl_vertex.size()); + //draw + d->hl_type = Locally_shortest_path_item_priv::VERTEX; + for(CGAL::QGLViewer* v : CGAL::QGLViewer::QGLViewerPool()) + { + CGAL::Three::Viewer_interface* viewer = + static_cast(v); + tc->initializeBuffers(viewer); + } + break; + } + default: + d->hl_type = Locally_shortest_path_item_priv::NO_TYPE; + break; + } + } + else + clearHL(); + redraw(); + + d->ready_to_hl = false; +} + +void Locally_shortest_path_item::clearHL() +{ + Viewer_interface* viewer = dynamic_cast(*CGAL::QGLViewer::QGLViewerPool().begin()); + viewer->makeCurrent(); + d->hl_normal.clear(); + d->hl_vertex.clear(); + + Tc* tc = getTriangleContainer(Priv::S_Spheres); + tc->reset_vbos(ALL); + tc->allocate(Tc::Flat_vertices, d->vertex_spheres.data(), + static_cast(d->vertex_spheres.size()*sizeof(float))); + tc->allocate(Tc::Flat_normals, + d->normal_spheres.data(), + static_cast(d->normal_spheres.size()*sizeof(float))); + + tc->allocate(Tc::Facet_centers, nullptr, 0); + + tc->initializeBuffers(viewer); + tc->setFlatDataSize(0); + tc->setCenterSize(0); + //draw + Ec* ec = getEdgeContainer(Priv::S_Edges); + ec->reset_vbos(ALL); + ec->allocate(Ec::Vertices, nullptr, 0); + ec->initializeBuffers(viewer); + ec->setFlatDataSize(0); + + tc = getTriangleContainer(Priv::S_Faces); + tc->reset_vbos(ALL); + tc->allocate(Tc::Flat_vertices, nullptr, 0); + tc->allocate(Tc::Flat_normals, nullptr, 0); + tc->initializeBuffers(viewer); + tc->setFlatDataSize(0); + d->hl_type = Locally_shortest_path_item_priv::NO_TYPE; + + itemChanged(); + +} +void Locally_shortest_path_item_priv::reset_selection() +{ + selection_on = false; + CGAL::QGLViewer* viewer = *CGAL::QGLViewer::QGLViewerPool().begin(); + viewer->setManipulatedFrame(frame); + viewer->setMouseBinding(Qt::ShiftModifier, Qt::LeftButton, CGAL::qglviewer::SELECT); + constraint.setTranslationConstraintType(CGAL::qglviewer::AxisPlaneConstraint::FREE); + selected_vertices.clear(); +} + +//intercept events for picking +bool Locally_shortest_path_item::eventFilter(QObject *obj, QEvent *event) +{ + if(!visible()) + return false; + Viewer_interface* viewer = qobject_cast(obj); + if(!viewer) + return false; + if(event->type() == QEvent::MouseButtonPress) + { + QMouseEvent* e = static_cast(event); + if(e->modifiers() == Qt::NoModifier) + { + //pick + int type, picked; + d->picked_pixel = e->pos(); + d->picking(type, picked, viewer); + viewer->makeCurrent(); + if(type == 0) + { + bool found = false; + QApplication::setOverrideCursor(Qt::DragMoveCursor); + CGAL::qglviewer::Vec pos = viewer->camera()->pointUnderPixel(d->picked_pixel, found); + if(found) + { + d->rf_last_pos = pos; + d->remodel_frame->setPosition(pos); + } + + d->selection_on = true; + d->selected_vertices.push_back(&d->vertices[picked]); + d->constraint.setTranslationConstraintType(CGAL::qglviewer::AxisPlaneConstraint::FREE); + d->remodel_frame->setConstraint(&d->constraint); + + viewer->setManipulatedFrame(d->remodel_frame); + viewer->setMouseBinding( + Qt::NoModifier, + Qt::LeftButton, + CGAL::qglviewer::FRAME, + CGAL::qglviewer::TRANSLATE); + } + else + { + d->reset_selection(); + } + } + return false; + } + else if(event->type() == QEvent::MouseMove) + { + QMouseEvent* e = static_cast(event); + if(e->modifiers() == Qt::NoModifier) + { + if(d->selection_on) + { + d->remodel_frame->setOrientation(d->frame->orientation()); + CGAL::qglviewer::Vec td(d->remodel_frame->position() - d->rf_last_pos); + QVector3D dir(td.x, td.y, td.z); + d->update_points(dir); + d->rf_last_pos=d->remodel_frame->position(); + d->path_invalidated=true; + } + d->ready_to_hl= true; + d->picked_pixel = e->pos(); + QTimer::singleShot(0, this, + [this, viewer](){ + highlight(viewer); + }); + } + else if(e->modifiers() == Qt::ControlModifier && + e->buttons() == Qt::LeftButton) + { + QApplication::setOverrideCursor(d->rotate_cursor); + } + else if(d->selection_on) + { + d->reset_selection(); + } + d->picked_pixel = e->pos(); + return false; + } + else if(event->type() == QEvent::MouseButtonRelease) + { + d->reset_selection(); + QApplication::setOverrideCursor(QCursor()); + viewer->setMouseBinding( + Qt::NoModifier, + Qt::LeftButton, + CGAL::qglviewer::CAMERA, + CGAL::qglviewer::ROTATE); + } + else if(event->type() == QEvent::KeyRelease) + { + QKeyEvent* e = static_cast(event); + if(e->key() == Qt::Key_Control) + { + QApplication::setOverrideCursor(QCursor()); + } + } + return false; +} + +void Locally_shortest_path_item_priv::draw_picking(Viewer_interface* viewer) +{ + + QMatrix4x4 f_matrix; + for (int i=0; i<16; ++i){ + f_matrix.data()[i] = (float)frame->matrix()[i]; + } + GLdouble d_mat[16]; + QMatrix4x4 mv_mat; + viewer->camera()->getModelViewMatrix(d_mat); + for (int i=0; i<16; ++i) + mv_mat.data()[i] = GLfloat(d_mat[i]); + mv_mat = mv_mat*f_matrix; + + // TODO: radius + double radius = 0.01 ; + Tc* tc = item->getTriangleContainer(P_Spheres); + tc->setFrameMatrix(f_matrix); + tc->setClipping(false); + tc->getVao(viewer)->bind(); + tc->getVao(viewer)->program->setAttributeValue("radius", (float)radius); + tc->getVao(viewer)->release(); + tc->draw(viewer, false); +} + +void Locally_shortest_path_item_priv::update_points(const QVector3D &dir) +{ + CGAL::qglviewer::AxisPlaneConstraint::Type prev_cons = constraint.translationConstraintType(); + constraint.setTranslationConstraintType(CGAL::qglviewer::AxisPlaneConstraint::FREE); + for(Locally_shortest_path_item::vertex* selected_vertex: selected_vertices ) + { + int id = selected_vertex->id; + CGAL_assume(id<8 && id >=0); + double x = selected_vertex->x + dir.x(); + double y = selected_vertex->y + dir.y(); + double z = selected_vertex->z + dir.z(); + + // project point onto the mesh + auto p = aabb_tree.closest_point(CGAL::Epick::Point_3(x,y,z)); + selected_vertex->x = p.x(); + selected_vertex->y = p.y(); + selected_vertex->z = p.z(); + } + item->invalidateOpenGLBuffers(); + constraint.setTranslationConstraintType(prev_cons); +} + +//type : 0 = vertex, 1 = edge, 2 = face +void Locally_shortest_path_item_priv::picking(int& type, int& id, Viewer_interface *viewer) +{ + viewer->makeCurrent(); + type = -1; + id = -1; + int deviceWidth = viewer->camera()->screenWidth(); + int deviceHeight = viewer->camera()->screenHeight(); + QOpenGLFramebufferObject* fbo = new QOpenGLFramebufferObject(deviceWidth, deviceHeight,QOpenGLFramebufferObject::Depth); + fbo->bind(); + viewer->glEnable(GL_DEPTH_TEST); + viewer->glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); + QColor bgColor(viewer->backgroundColor()); + //draws the image in the fbo + viewer->setBackgroundColor(::Qt::white); + draw_picking(viewer); + + const auto buffer = read_pixel_as_ubyte_rgba(picked_pixel, viewer, viewer->camera()); + //decode ID and pick (don't forget the case nothing is picked + if(!(buffer[0]==buffer[1] && buffer[1]==buffer[2])) + { + int r(std::ceil((buffer[0]-10)/20)), g(std::ceil((buffer[1]-10)/20)), b(std::ceil((buffer[2]-10)/20)); + id = (std::max)(r,g); + id = (std::max)(id,b); + if(buffer[0] > 0) + { + if(id <8) + type = 0 ; + } + else if(buffer[1] > 0) + { + if(id <12) + { + type = 1; + } + } + else if(buffer[2] > 0) + { + if(id <6) + { + type = 2; + } + } + } + viewer->setBackgroundColor(bgColor); + fbo->release(); + delete fbo; +} + +void Locally_shortest_path_item::drawHl(Viewer_interface* viewer)const +{ + QMatrix4x4 f_matrix; + for (int i=0; i<16; ++i){ + f_matrix.data()[i] = (float)d->frame->matrix()[i]; + } + GLdouble d_mat[16]; + QMatrix4x4 mv_mat; + viewer->camera()->getModelViewMatrix(d_mat); + for (int i=0; i<16; ++i) + mv_mat.data()[i] = GLfloat(d_mat[i]); + mv_mat = mv_mat*f_matrix; + + if(d->hl_type == Locally_shortest_path_item_priv::VERTEX) + { + Tc* tc = getTriangleContainer(Priv::S_Spheres); + + tc->setFrameMatrix(f_matrix); + tc->setMvMatrix(mv_mat); + tc->setColor(QColor(Qt::yellow)); + + // TODO radius + double radius = 0.01 ; + tc->setClipping(false); + tc->getVao(viewer)->bind(); + tc->getVao(viewer)->program->setUniformValue("radius", (float)radius); + tc->getVao(viewer)->release(); + tc->draw(viewer, true); + } +} + +void Locally_shortest_path_item::invalidateOpenGLBuffers() +{ + compute_bbox(); + setBuffersFilled(false); + getTriangleContainer(Priv::Spheres)->reset_vbos(ALL); + getTriangleContainer(Priv::P_Spheres)->reset_vbos(ALL); +} + +void Locally_shortest_path_item::computeElements() const +{ + d->computeElements(); + + Tc* tc = getTriangleContainer(Priv::Spheres); + tc->allocate( + Tc::Flat_vertices, + d->vertex_spheres.data(), + static_cast(d->vertex_spheres.size()*sizeof(float))); + + tc->allocate( + Tc::Flat_normals, + d->normal_spheres.data(), + static_cast(d->normal_spheres.size()*sizeof(float))); + tc->allocate( + Tc::Facet_centers, + d->center_spheres.data(), + static_cast(d->center_spheres.size()*sizeof(float))); + + tc = getTriangleContainer(Priv::P_Spheres); + tc->allocate( + Tc::Flat_vertices, + d->vertex_spheres.data(), + static_cast(d->vertex_spheres.size()*sizeof(float))); + + tc->allocate( + Tc::Facet_centers, + d->center_spheres.data(), + static_cast(d->center_spheres.size()*sizeof(float))); + + tc->allocate( + Tc::FColors, + d->color_spheres.data(), + static_cast(d->color_spheres.size()*sizeof(float))); + setBuffersFilled(true); +} + +void Locally_shortest_path_item::initializeBuffers(Viewer_interface *v) const +{ + getTriangleContainer(Priv::Spheres)->initializeBuffers(v); + getTriangleContainer(Priv::Spheres)->initializeBuffers(v); + getTriangleContainer(Priv::P_Spheres)->initializeBuffers(v); + getTriangleContainer(Priv::P_Spheres)->initializeBuffers(v); + + getTriangleContainer(Priv::Spheres)->setFlatDataSize(d->vertex_spheres.size()); + getTriangleContainer(Priv::Spheres)->setCenterSize(d->center_spheres.size()); + getTriangleContainer(Priv::P_Spheres)->setFlatDataSize(d->vertex_spheres.size()); + getTriangleContainer(Priv::P_Spheres)->setCenterSize(d->center_spheres.size()); +} + +void Locally_shortest_path_item::connectNewViewer(QObject *o) +{ + Vi* viewer = qobject_cast(o); + if(!viewer) + return; + viewer->setMouseTracking(true); +} diff --git a/Lab/demo/Lab/Plugins/Bsurf/Locally_shortest_path_item.h b/Lab/demo/Lab/Plugins/Bsurf/Locally_shortest_path_item.h new file mode 100644 index 00000000000..29a131b37d4 --- /dev/null +++ b/Lab/demo/Lab/Plugins/Bsurf/Locally_shortest_path_item.h @@ -0,0 +1,70 @@ +#ifndef LOCALLY_SHORTEST_PATH_ITEM_H +#define LOCALLY_SHORTEST_PATH_ITEM_H + +#include +#include +#include "Scene_surface_mesh_item.h" +#include "Scene_polylines_item.h" +#include "create_sphere.h" +#include "Locally_shortest_path_item_config.h" + +struct Locally_shortest_path_item_priv; +class LOCALLY_SHORTEST_PATH_ITEM_EXPORT Locally_shortest_path_item: + public CGAL::Three::Scene_item_rendering_helper +{ + Q_OBJECT + public: + typedef CGAL::Simple_cartesian Kernel; + struct vertex; + struct edge; + struct face; + + Locally_shortest_path_item(const CGAL::Three::Scene_interface* scene_interface, + const Scene_surface_mesh_item* sm_item, + Scene_polylines_item* polyline_item, + std::size_t nb_pts); + ~Locally_shortest_path_item(); + bool isFinite() const { return true; } + bool isEmpty() const { return false; } + void compute_bbox() const; + + bool manipulatable() const { return true; } + ManipulatedFrame* manipulatedFrame(); + Locally_shortest_path_item* clone() const { + return nullptr; + } + + QString toolTip() const; + QMenu* contextMenu(); + + bool eventFilter(QObject *, QEvent *); + // Indicate if rendering mode is supported + bool supportsRenderingMode(RenderingMode m) const; + void draw(CGAL::Three::Viewer_interface *) const; + void drawHl(CGAL::Three::Viewer_interface *) const; + //~ void drawEdges(CGAL::Three::Viewer_interface* viewer) const; + void drawSpheres(CGAL::Three::Viewer_interface* viewer, const QMatrix4x4 f_matrix) const; + void drawPath() const; + void invalidateOpenGLBuffers(); + + // 5-----6 + // . | . | + // 4------7 | + // | | | | + // | 1-|---2 + // | . |. + // 0------3 + + double point(short i, short j) const; + void initializeBuffers(CGAL::Three::Viewer_interface *) const; + void computeElements() const; +public Q_SLOTS: + void highlight(CGAL::Three::Viewer_interface* viewer); + void clearHL(); + void connectNewViewer(QObject* o); + +protected: + friend struct Locally_shortest_path_item_priv; + Locally_shortest_path_item_priv* d; +}; +#endif // SCENE_EDIT_BOX_ITEM_H diff --git a/Lab/demo/Lab/Plugins/Bsurf/Locally_shortest_path_item_config.h b/Lab/demo/Lab/Plugins/Bsurf/Locally_shortest_path_item_config.h new file mode 100644 index 00000000000..06796b33eba --- /dev/null +++ b/Lab/demo/Lab/Plugins/Bsurf/Locally_shortest_path_item_config.h @@ -0,0 +1,10 @@ +#ifndef LOCALLY_SHORTEST_PATH_ITEM_CONFIG_H +#define LOCALLY_SHORTEST_PATH_ITEM_CONFIG_H + +#ifdef locally_shortest_path_item_EXPORTS +# define LOCALLY_SHORTEST_PATH_ITEM_EXPORT Q_DECL_EXPORT +#else +# define LOCALLY_SHORTEST_PATH_ITEM_EXPORT Q_DECL_IMPORT +#endif + +#endif // LOCALLY_SHORTEST_PATH_ITEM_CONFIG_H diff --git a/Lab/demo/Lab/Plugins/Bsurf/Locally_shortest_path_plugin.cpp b/Lab/demo/Lab/Plugins/Bsurf/Locally_shortest_path_plugin.cpp new file mode 100644 index 00000000000..986975987b4 --- /dev/null +++ b/Lab/demo/Lab/Plugins/Bsurf/Locally_shortest_path_plugin.cpp @@ -0,0 +1,151 @@ +#include + +#include +#include +#include "Locally_shortest_path_item.h" +#include "Scene_surface_mesh_item.h" +#include "Scene_polylines_item.h" +#include +#include +#include +#include +#include +#include + + +#include +#include +using namespace CGAL::Three; +class Locally_shortest_path_plugin : + public QObject, + public CGAL_Lab_plugin_interface +{ + Q_OBJECT + Q_INTERFACES(CGAL::Three::CGAL_Lab_plugin_interface) + Q_PLUGIN_METADATA(IID "com.geometryfactory.PolyhedronDemo.PluginInterface/1.0") + +public: + void init(QMainWindow* mainWindow, CGAL::Three::Scene_interface* scene_interface, Messages_interface*); + QList actions() const { + return QList() << actionTracePath << actionTraceBezier; + } + + bool applicable(QAction*) const { + if(scene->numberOfEntries() > 0) + { + int item_id = scene->mainSelectionIndex(); + return qobject_cast(scene->item(item_id)); + } + return false; + } +public Q_SLOTS: + + void trace_path(); + void trace_bezier(); + void enableAction(); + void connectNewViewer(QObject* o) + { + for(int i=0; inumberOfEntries(); ++i) + { + Locally_shortest_path_item* item = qobject_cast( + scene->item(i)); + if(item) + o->installEventFilter(item); + } + } + +private: + CGAL::Three::Scene_interface* scene; + QMainWindow* mw; + QAction* actionTracePath; + QAction* actionTraceBezier; +}; // end Locally_shortest_path_plugin + +void Locally_shortest_path_plugin::init(QMainWindow* mainWindow, CGAL::Three::Scene_interface* scene_interface, Messages_interface*) +{ + scene = scene_interface; + mw = mainWindow; + actionTracePath = new QAction(tr("Create Locally Shortest Path"), mainWindow); + actionTraceBezier= new QAction(tr("Trace Bezier Curve"), mainWindow); + connect(actionTracePath, SIGNAL(triggered()), + this, SLOT(trace_path())); + connect(actionTraceBezier, SIGNAL(triggered()), + this, SLOT(trace_bezier())); + connect(mw, SIGNAL(newViewerCreated(QObject*)), + this, SLOT(connectNewViewer(QObject*))); +} + +void Locally_shortest_path_plugin::trace_path() +{ + for(int i = 0, end = scene->numberOfEntries(); + i < end; ++i) + { + if(qobject_cast(scene->item(i))) + return; + } + QApplication::setOverrideCursor(Qt::WaitCursor); + + Scene_surface_mesh_item* sm_item = qobject_cast(scene->item(scene->mainSelectionIndex())); + + Scene_polylines_item* polyline_item = new Scene_polylines_item(); + + polyline_item->setName(tr("Locally Shortest Path")); + polyline_item->setColor(Qt::red); + scene->addItem(polyline_item); + polyline_item->invalidateOpenGLBuffers(); + + Locally_shortest_path_item* item = new Locally_shortest_path_item(scene, sm_item, polyline_item,2); + connect(item, SIGNAL(destroyed()), + this, SLOT(enableAction())); + item->setName("Source and Target Points"); + item->setRenderingMode(FlatPlusEdges); + for(CGAL::QGLViewer* viewer : CGAL::QGLViewer::QGLViewerPool()) + viewer->installEventFilter(item); + + scene->addItem(item); + actionTracePath->setEnabled(false); + actionTraceBezier->setEnabled(false); + + QApplication::restoreOverrideCursor(); +} + +void Locally_shortest_path_plugin::trace_bezier() +{ + for(int i = 0, end = scene->numberOfEntries(); + i < end; ++i) + { + if(qobject_cast(scene->item(i))) + return; + } + QApplication::setOverrideCursor(Qt::WaitCursor); + + Scene_surface_mesh_item* sm_item = qobject_cast(scene->item(scene->mainSelectionIndex())); + + Scene_polylines_item* polyline_item = new Scene_polylines_item(); + + polyline_item->setName(tr("Locally Shortest Path")); + polyline_item->setColor(Qt::red); + scene->addItem(polyline_item); + polyline_item->invalidateOpenGLBuffers(); + + Locally_shortest_path_item* item = new Locally_shortest_path_item(scene, sm_item, polyline_item,4); + connect(item, SIGNAL(destroyed()), + this, SLOT(enableAction())); + item->setName("Control Points"); + item->setRenderingMode(FlatPlusEdges); + for(CGAL::QGLViewer* viewer : CGAL::QGLViewer::QGLViewerPool()) + viewer->installEventFilter(item); + + scene->addItem(item); + actionTracePath->setEnabled(false); + actionTraceBezier->setEnabled(false); + + QApplication::restoreOverrideCursor(); +} + +void Locally_shortest_path_plugin::enableAction() { + actionTracePath->setEnabled(true); + actionTraceBezier->setEnabled(true); +} + +#include "Locally_shortest_path_plugin.moc" diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/locate.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/locate.h index c8653c0196d..184c66b2782 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/locate.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/locate.h @@ -59,6 +59,12 @@ template using descriptor_variant = std::variant::vertex_descriptor, typename boost::graph_traits::halfedge_descriptor, typename boost::graph_traits::face_descriptor>; +/// \ingroup PMP_locate_grp +/// +/// A variant used in the function `get_descriptor_from_location()` overload for `Edge_location`. +template +using vertex_or_halfedge_variant = std::variant::vertex_descriptor, + typename boost::graph_traits::halfedge_descriptor>; /// \ingroup PMP_locate_grp /// @@ -81,6 +87,109 @@ template using Face_location = std::pair::face_descriptor, Barycentric_coordinates >; +/// \ingroup PMP_locate_grp +/// +/// If `pm` is the input surface mesh and given the pair (`e`, `bc`) +/// such that `bc` is a pair of barycentric coordinates `(w0, w1)`, the correspondence +/// between the coordinates in `bc` and the vertices of the edge `e` is the following: +/// - `w0` corresponds to `source(e, pm)` +/// - `w1` corresponds to `target(e, pm)` +template +using Edge_location = + std::pair::edge_descriptor, + std::array>; + +/// \ingroup PMP_locate_grp +/// returns `true` if `loc1` and `loc2` indicate the same point on `tm`, and `false` otherwise, +template +bool +are_locations_identical(const Face_location& loc1, + const Face_location& loc2, + const TriangleMesh& tm) +{ + if (loc1==loc2) return true; + + using vertex_descriptor = typename boost::graph_traits::vertex_descriptor; + using halfedge_descriptor = typename boost::graph_traits::halfedge_descriptor; + using face_descriptor = typename boost::graph_traits::face_descriptor; + using Barycentric_coordinates = CGAL::Polygon_mesh_processing::Barycentric_coordinates; + + const face_descriptor fd1 = loc1.first; + const Barycentric_coordinates& bar1 = loc1.second; + const face_descriptor fd2 = loc2.first; + const Barycentric_coordinates& bar2 = loc2.second; + + // the first barycentric coordinate corresponds to source(halfedge(fdX, tm), tm) + halfedge_descriptor hd1 = prev(halfedge(fd1, tm), tm); + halfedge_descriptor hd2 = prev(halfedge(fd2, tm), tm); + + // check if the point is a vertex + for(int i=0; i<3; ++i) + { + if(bar1[i] == FT(1)) + { + vertex_descriptor vd1 = target(hd1, tm); + for(int k=0; k<3; ++k) + { + if(bar2[k] == FT(1)) + return target(hd2, tm) == vd1; + hd2 = next(hd2, tm); + } + } + hd1 = next(hd1, tm); + } + CGAL_assertion(hd1 == prev(halfedge(fd1, tm), tm)); + + // check loc2 is not a vertex + for(int k=0; k<3; ++k) + { + if(bar2[k] == FT(1)) + return false; + hd2 = next(hd2, tm); + } + CGAL_assertion(hd2 == prev(halfedge(fd2, tm), tm)); + + + // check if the point is on an edge + for(int i=0; i<3; ++i) + { + if(bar1[i] == FT(0)) // coordinate at target(hd1, tm) + { + for(int k=0; k<3; ++k) + { + if(bar2[k] == FT(0)) // coordinate at target(hd2, tm) + return prev(hd2, tm) == opposite(prev(hd1, tm), tm); + hd2 = next(hd2, tm); + } + } + hd1 = next(hd1, tm); + } + + return false; +} + +template +Face_location +to_face_location(Edge_location loc, + const TriangleMesh& tm) +{ + auto h = halfedge(loc.first, tm); + if (is_border(h, tm)) + { + h=opposite(h, tm); + std::swap(loc.second[0], loc.second[1]); + } + + auto f = face(h, tm); + auto hf = halfedge(f, tm); + if (hf == h) + return Face_location(f, CGAL::make_array(loc.second[0], loc.second[1], FT(0))); + if (next(hf, tm)==h) + return Face_location(f, CGAL::make_array(FT(0),loc.second[0], loc.second[1])); + CGAL_assertion(prev(hf, tm)==h); + return Face_location(f, CGAL::make_array(loc.second[1],FT(0),loc.second[0])); +} + namespace internal { // The Ray must have the same ambient dimension as the property map's value type (aka, the point type) @@ -339,6 +448,18 @@ struct Barycentric_point_constructor // 3D version return P(x, y, z); } + + P operator()(const P& p, const FT wp, const P& q, const FT wq, + const K& /*k*/) const + { + FT sum = wp + wq; + CGAL_assertion(sum != FT(0)); + FT x = (wp * p.x() + wq * q.x()) / sum; + FT y = (wp * p.y() + wq * q.y()) / sum; + FT z = (wp * p.z() + wq * q.z()) / sum; + + return P(x, y, z); + } }; } // namespace internal @@ -503,11 +624,7 @@ random_location_on_mesh(const TriangleMesh& tm, /// template descriptor_variant -#ifdef DOXYGEN_RUNNING // just for convenience because template alias do not allow template deduction get_descriptor_from_location(const Face_location& loc, -#else -get_descriptor_from_location(const Face_location& loc, -#endif const TriangleMesh& tm) { typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; @@ -547,7 +664,123 @@ get_descriptor_from_location(const Face_location& loc, /// \ingroup PMP_locate_grp /// -/// \brief Given a location in a face, returns the geometric position described +/// \brief Given a location, if it designates a vertex, returns the corresponding vertex descriptor and `std::nullopt` otherwise. +/// +/// \tparam FT must be a model of `FieldNumberType` +/// \tparam TriangleMesh must be a model of `FaceGraph` +/// +/// \param loc a location with `loc.first` a face of `tm` +/// \param tm a triangulated surface mesh +/// +/// \pre `loc.first` is a face descriptor corresponding to a face of `tm`. +/// \pre `loc` describes the barycentric coordinates of a point that lives within the face (boundary included), +/// meaning the barycentric coordinates are all positive. +/// +template +std::optional::vertex_descriptor> +vertex_descriptor_from_location(const Face_location& loc, + const TriangleMesh& tm) +{ + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + + typedef Barycentric_coordinates Barycentric_coordinates; + + const face_descriptor fd = loc.first; + const Barycentric_coordinates& bar = loc.second; + + CGAL_precondition(is_valid_face_descriptor(fd, tm)); + CGAL_precondition(is_triangle(halfedge(fd, tm), tm)); + CGAL_precondition(is_in_face(loc, tm)); + + // the first barycentric coordinate corresponds to source(halfedge(fd, tm), tm) + halfedge_descriptor hd = prev(halfedge(fd, tm), tm); + + // check if the point is a vertex + for(int i=0; i<3; ++i) + { + if(bar[i] == FT(1)) // coordinate at target(hd, tm) + return target(hd, tm); + hd = next(hd, tm); + } + return std::nullopt; +} + +/// \ingroup PMP_locate_grp +/// +/// \brief Given a location, if it designates a vertex, returns the corresponding vertex descriptor and `std::nullopt` otherwise. +/// +/// \tparam FT must be a model of `FieldNumberType` +/// \tparam TriangleMesh must be a model of `FaceGraph` +/// +/// \param loc a location with `loc.first` a face of `tm` +/// \param tm a triangulated surface mesh +/// +/// \pre `loc.first` is a face descriptor corresponding to a face of `tm`. +/// \pre `loc` describes the barycentric coordinates of a point that lives within the face (boundary included), +/// meaning the barycentric coordinates are all positive. +/// +template +std::optional::vertex_descriptor> +vertex_descriptor_from_location(const Edge_location& loc, + const TriangleMesh& tm) +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + vertex_descriptor src = source(loc.first, tm), tgt = target(loc.first, tm); + const auto& bar = loc.second; + + if (bar[1]==0) return src; + if (bar[0]==0) return tgt; + + return std::nullopt; +} + + +/// \ingroup PMP_locate_grp +/// +/// \brief Given a location, returns a descriptor to the simplex of smallest dimension +/// on which the point corresponding to the location lies. +/// +/// \details In other words: +/// - if the point lies on a vertex, this function returns a `boost::graph_traits::%vertex_descriptor` `v`; +/// - if the point lies on a halfedge, this function returns a `boost::graph_traits::%halfedge_descriptor` `hd` +/// (note that in that case, `loc.first == face(hd, tm)` holds). +/// +/// \tparam FT must be a model of `FieldNumberType` +/// \tparam TriangleMesh must be a model of `FaceGraph` +/// +/// \param loc a location with `loc.first` a face of `tm` +/// \param tm a triangulated surface mesh +/// +/// \pre `loc.first` is an edge descriptor corresponding to an edge of `tm`. +/// \pre `loc` describes the barycentric coordinates of a point that lives within the edge (boundary included), +/// meaning the barycentric coordinates are all positive. +/// +template +vertex_or_halfedge_variant +get_descriptor_from_location(const Edge_location& loc, + const TriangleMesh& tm) +{ + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + + typedef std::array Barycentric_coordinates; + + const edge_descriptor ed = loc.first; + const Barycentric_coordinates& bar = loc.second; + + if (bar[0]==FT(0)) + return source(ed, tm); + + if (bar[1]==FT(0)) + return target(ed, tm); + + return halfedge(ed, tm); +} + +/// \ingroup PMP_locate_grp +/// +/// \brief returns, given a location in a face, the geometric position described /// by these coordinates, as a point. /// /// \tparam FT must be a model of `FieldNumberType` @@ -616,6 +849,83 @@ construct_point(const Face_location& loc, return bp_constructor(p0, loc.second[0], p1, loc.second[1], p2, loc.second[2], gt); } +/// \ingroup PMP_locate_grp +/// +/// \brief returns, given a location in an edge, the geometric position described +/// by these coordinates, as a point. +/// +/// \tparam FT must be a model of `FieldNumberType` +/// \tparam TriangleMesh must be a model of `FaceGraph` +/// \tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" +/// +/// \param loc the location from which a point is constructed +/// \param tm a triangulated surface mesh +/// \param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below +/// +/// \cgalNamedParamsBegin +/// \cgalParamNBegin{vertex_point_map} +/// \cgalParamDescription{a property map associating points to the vertices of `tm`} +/// \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits::%vertex_descriptor` +/// as key type and `%Point_3` as value type} +/// \cgalParamDefault{`boost::get(CGAL::vertex_point, tm)`} +/// \cgalParamNEnd +/// +/// \cgalParamNBegin{geom_traits} +/// \cgalParamDescription{an instance of a geometric traits class} +/// \cgalParamType{a class model of `Kernel`} +/// \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`} +/// \cgalParamExtra{The geometric traits class must be compatible with the vertex point type.} +/// \cgalParamExtra{If such traits class is provided, its type `FT` must be identical +/// to the template parameter `FT` of this function.} +/// \cgalParamNEnd +/// \cgalNamedParamsEnd +/// +/// \pre `loc.first` is a face descriptor corresponding to a face of `tm`. +/// +/// \returns a point whose type is the same as the value type of the vertex point property map +/// provided by the user or via named parameters, or the internal point map of the mesh `tm`. +/// +template +#ifdef DOXYGEN_RUNNING +Point construct_point( + const Edge_location &loc, +#else +typename internal::Location_traits::Point +construct_point(const Edge_location &loc, +#endif + const TriangleMesh &tm, + const NamedParameters &np = parameters::default_values()) { + typedef typename boost::graph_traits::edge_descriptor + edge_descriptor; + typedef + typename GetGeomTraits::type Geom_traits; + + typedef typename GetVertexPointMap::const_type + VertexPointMap; + typedef typename boost::property_traits::value_type Point; + typedef typename boost::property_traits::reference + Point_reference; + + using parameters::choose_parameter; + using parameters::get_parameter; + + CGAL_precondition(CGAL::is_triangle_mesh(tm)); + + VertexPointMap vpm = parameters::choose_parameter( + parameters::get_parameter(np, internal_np::vertex_point), + get_const_property_map(boost::vertex_point, tm)); + Geom_traits gt = choose_parameter( + get_parameter(np, internal_np::geom_traits)); + + edge_descriptor ed = loc.first; + const Point_reference p0 = get(vpm, source(ed, tm)); + const Point_reference p1 = get(vpm, target(ed, tm)); + + internal::Barycentric_point_constructor bp_constructor; + return bp_constructor(p0, loc.second[0], p1, loc.second[1], gt); +} + /// \name Location Predicates /// @{ diff --git a/Vector_graphics_on_surfaces/doc/Vector_graphics_on_surfaces/Doxyfile.in b/Vector_graphics_on_surfaces/doc/Vector_graphics_on_surfaces/Doxyfile.in new file mode 100644 index 00000000000..e1d1b533f78 --- /dev/null +++ b/Vector_graphics_on_surfaces/doc/Vector_graphics_on_surfaces/Doxyfile.in @@ -0,0 +1,3 @@ +@INCLUDE = ${CGAL_DOC_PACKAGE_DEFAULTS} + +PROJECT_NAME = "CGAL ${CGAL_DOC_VERSION} - Vector Graphics on Triangulated Surface Meshes" diff --git a/Vector_graphics_on_surfaces/doc/Vector_graphics_on_surfaces/PackageDescription.txt b/Vector_graphics_on_surfaces/doc/Vector_graphics_on_surfaces/PackageDescription.txt new file mode 100644 index 00000000000..c4ce5e4b288 --- /dev/null +++ b/Vector_graphics_on_surfaces/doc/Vector_graphics_on_surfaces/PackageDescription.txt @@ -0,0 +1,68 @@ +/// \defgroup PkgVGoS_Ref Vector Graphics on Triangulated Surface Meshes Reference +/// \defgroup VGSConcepts Concepts +/// \ingroup PkgVGoS_Ref + +/// \defgroup VGSAlgorithmClasses Algorithm Classes +/// \ingroup PkgVGoS_Ref + +/// \defgroup VGSFunctions Functions +/// \ingroup PkgVGoS_Ref + +/// \defgroup VGSTraitsClasses Traits Classes +/// \ingroup PkgVGoS_Ref + +/// \defgroup VGSMiscellaneous Miscellaneous +/// \ingroup PkgVGoS_Ref + +/*! +\addtogroup PkgVGoS_Ref +\cgalPkgDescriptionBegin{Vector Graphics on Triangulated Surface Meshes,PkgVGoS} +\cgalPkgPicture{vgos-small.png} +\cgalPkgSummaryBegin +\cgalPkgAuthors{Claudio Mancinelli and Sébastien Loriot} +\cgalPkgDesc{The package provides functions to draw and manipulate basic geometric primitives such as geodesics curves and Bézier splines + on triangulated surface meshes.} +\cgalPkgManuals{Chapter_VGoS,PkgVGoS_Ref} +\cgalPkgSummaryEnd +\cgalPkgShortInfoBegin +\cgalPkgSince{6.2} +\cgalPkgDependsOn{\ref PkgPolygonMeshProcessing} +\cgalPkgBib{cgal:ml-vgos} +\cgalPkgLicense{\ref licensesGPL "GPL"} +\cgalPkgDemo{CGAL Lab,CGALlab.zip} +\cgalPkgShortInfoEnd +\cgalPkgDescriptionEnd + + +\cgalClassifedRefPages + +\cgalCRPSection{Base Functions} +- `CGAL::Vector_graphics_on_surfaces::locally_shortest_path()` +- `CGAL::Vector_graphics_on_surfaces::recursive_de_Casteljau()` +- `CGAL::Vector_graphics_on_surfaces::straightest_geodesic()` + +\cgalCRPSection{Higher Level Tracing Functions} +- `CGAL::Vector_graphics_on_surfaces::trace_geodesic_polygon()` +- `CGAL::Vector_graphics_on_surfaces::trace_geodesic_polygons()` +- `CGAL::Vector_graphics_on_surfaces::trace_geodesic_label()` +- `CGAL::Vector_graphics_on_surfaces::trace_geodesic_label_along_curve` +- `CGAL::Vector_graphics_on_surfaces::trace_bezier_curves()` +- `CGAL::Vector_graphics_on_surfaces::trace_polyline_of_bezier_curves()` + +\cgalCRPSection{Utilities} +- `CGAL::Vector_graphics_on_surfaces::path_length()` +- `CGAL::Vector_graphics_on_surfaces::convert_path_to_polyline()` +- `CGAL::Vector_graphics_on_surfaces::convert_to_polar_coordinates` + +\cgalCRPSection{Miscellaneous} +- `CGAL::Vector_graphics_on_surfaces::Bezier_segment` +- `CGAL::Vector_graphics_on_surfaces::Dual_geodesic_solver` +- `CGAL::Vector_graphics_on_surfaces::init_geodesic_dual_solver()` +- `CGAL::Vector_graphics_on_surfaces::approximate_geodesic_distance_field()` +- `CGAL::Vector_graphics_on_surfaces::refine_mesh_along_paths()` + + +\todo tracing functions should probably be merged in only one or two functions, with parameters for the placement of centers + parallel transport, as well as an option for joining with locally shortest paths + + +*/ diff --git a/Vector_graphics_on_surfaces/doc/Vector_graphics_on_surfaces/Vector_graphics_on_surfaces.txt b/Vector_graphics_on_surfaces/doc/Vector_graphics_on_surfaces/Vector_graphics_on_surfaces.txt new file mode 100644 index 00000000000..0a50533331f --- /dev/null +++ b/Vector_graphics_on_surfaces/doc/Vector_graphics_on_surfaces/Vector_graphics_on_surfaces.txt @@ -0,0 +1,22 @@ +namespace CGAL { +/*! + +\mainpage User Manual +\anchor Chapter_VGoS +\cgalAutoToc +\author Claudio Mancinelli and Sébastien Loriot + +This chapter describes the ... + +\section vgos_definitions Definitions + +Section on definitions here ... + +\section vgos_examples Examples + +\subsection vgos_examples_first First Example + +The following example shows ... + +*/ +} /* namespace CGAL */ diff --git a/Vector_graphics_on_surfaces/doc/Vector_graphics_on_surfaces/dependencies b/Vector_graphics_on_surfaces/doc/Vector_graphics_on_surfaces/dependencies new file mode 100644 index 00000000000..26c05c82d78 --- /dev/null +++ b/Vector_graphics_on_surfaces/doc/Vector_graphics_on_surfaces/dependencies @@ -0,0 +1,8 @@ +Manual +Kernel_23 +STL_Extension +Algebraic_foundations +Circulator +Stream_support +Polygon_mesh_processing +BGL diff --git a/Vector_graphics_on_surfaces/doc/Vector_graphics_on_surfaces/examples.txt b/Vector_graphics_on_surfaces/doc/Vector_graphics_on_surfaces/examples.txt new file mode 100644 index 00000000000..e69de29bb2d diff --git a/Vector_graphics_on_surfaces/doc/Vector_graphics_on_surfaces/fig/vgos-small.png b/Vector_graphics_on_surfaces/doc/Vector_graphics_on_surfaces/fig/vgos-small.png new file mode 100644 index 00000000000..aab567220e5 Binary files /dev/null and b/Vector_graphics_on_surfaces/doc/Vector_graphics_on_surfaces/fig/vgos-small.png differ diff --git a/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/CMakeLists.txt b/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/CMakeLists.txt new file mode 100644 index 00000000000..a7bb19178b5 --- /dev/null +++ b/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/CMakeLists.txt @@ -0,0 +1,41 @@ +# Created by the script cgal_create_CMakeLists +# This is the CMake script for compiling a set of CGAL applications. + +cmake_minimum_required(VERSION 3.12...3.29) +project(Vector_graphics_on_surfaces_Examples) + +# CGAL and its components +find_package(CGAL REQUIRED) + +find_package(NanoSVG) +include(CGAL_NanoSVG_support) + +find_package(Eigen3 3.2.0 QUIET) #(requires 3.2.0 or greater) +include(CGAL_Eigen3_support) +if(TARGET CGAL::Eigen3_support) + create_single_source_cgal_program("trace_bezier_segment_sm_example.cpp") + target_link_libraries(trace_bezier_segment_sm_example PUBLIC CGAL::Eigen3_support) + create_single_source_cgal_program("geodesic_circles_sm_example.cpp") + target_link_libraries(geodesic_circles_sm_example PUBLIC CGAL::Eigen3_support) + create_single_source_cgal_program("straightest_geodesic_sm_example.cpp") + target_link_libraries(straightest_geodesic_sm_example PUBLIC CGAL::Eigen3_support) + create_single_source_cgal_program("trace_polygon_example.cpp") + target_link_libraries(trace_polygon_example PUBLIC CGAL::Eigen3_support) + create_single_source_cgal_program("trace_regular_polygons_example.cpp") + target_link_libraries(trace_regular_polygons_example PUBLIC CGAL::Eigen3_support) + create_single_source_cgal_program("trace_polygon_on_polygonal_curve_example.cpp") + target_link_libraries(trace_polygon_on_polygonal_curve_example PUBLIC CGAL::Eigen3_support) + create_single_source_cgal_program("locally_shortest_path_sm_example.cpp") + target_link_libraries(locally_shortest_path_sm_example PUBLIC CGAL::Eigen3_support) + if (NanoSVG_FOUND) + create_single_source_cgal_program("trace_svg_example.cpp") + target_link_libraries(trace_svg_example PUBLIC CGAL::Eigen3_support CGAL::NanoSVG_support) + create_single_source_cgal_program("trace_polygon_on_svg_curve_example.cpp") + target_link_libraries(trace_polygon_on_svg_curve_example PUBLIC CGAL::Eigen3_support CGAL::NanoSVG_support) + else() + message(STATUS "NOTICE: trace_svg_example requires NanoSVG and will not be compiled.") + endif() +else() + message(STATUS "NOTICE: Examples that use Eigen will not be compiled.") +endif() + diff --git a/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/geodesic_circles_sm_example.cpp b/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/geodesic_circles_sm_example.cpp new file mode 100644 index 00000000000..c24d173bc38 --- /dev/null +++ b/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/geodesic_circles_sm_example.cpp @@ -0,0 +1,74 @@ +#include +#include +#include + +#include + +namespace VGoS = CGAL::Vector_graphics_on_surfaces; +namespace PMP = CGAL::Polygon_mesh_processing; + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Surface_mesh Mesh; +typedef PMP::Face_location Face_location; +typedef PMP::Edge_location Edge_location; + + +int main(int argc, char** argv) +{ + std::string filename = (argc > 1) ? std::string(argv[1]) + : CGAL::data_file_path("meshes/elephant.off"); + + Mesh mesh; + if(!CGAL::IO::read_polygon_mesh(filename, mesh) || !CGAL::is_triangle_mesh(mesh)) + { + std::cerr << "Invalid input." << std::endl; + return 1; + } + + std::size_t nb_faces = faces(mesh).size(); + + // take two random faces and pick the centroid + CGAL::Random rnd = CGAL::get_default_random(); + // CGAL::Random rnd(1695720148); + // CGAL::Random rnd(1695724381); + // CGAL::Random rnd(1695813638); + + std::cout << "seed " << rnd.get_seed() << std::endl; + Mesh::Face_index f = *std::next(faces(mesh).begin(), rnd.get_int(0, nb_faces)); + // or pick specific faces + + Face_location src(f, CGAL::make_array(0.3,0.3,0.4)); + std::vector radii = {0.001, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 10}; + + std::cout << "src = " << PMP::construct_point(src, mesh) << "\n"; + + auto distance_map = mesh.add_property_map("v:dstmap").first; + VGoS::approximate_geodesic_distance_field(src, distance_map, mesh); + + std::ofstream out("circles.polylines.txt"); + + for (double r : radii) + { + for (Mesh::Face_index f : faces(mesh)) + { + std::vector pts; + Mesh::Halfedge_index h = halfedge(f, mesh); + for (int i=0; i<3; ++i) + { + Mesh::Vertex_index src = source(h, mesh), tgt = target(h, mesh); + double ds = get(distance_map, src); + double dt = get(distance_map, tgt); + if ((ds < r) != (dt < r)) + { + double alpha = (r - dt) / (ds - dt); + pts.push_back( CGAL::barycenter(mesh.point(src), alpha, mesh.point(tgt), 1-alpha) ); + } + h=next(h, mesh); + } + if (pts.size()==2) + out << "2 " << pts[0] << " " << pts[1] << "\n"; + } + } + + return 0; +} diff --git a/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/locally_shortest_path_sm_example.cpp b/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/locally_shortest_path_sm_example.cpp new file mode 100644 index 00000000000..0347c57b78a --- /dev/null +++ b/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/locally_shortest_path_sm_example.cpp @@ -0,0 +1,82 @@ +#include +#include +#include + +#include + +namespace VGoS = CGAL::Vector_graphics_on_surfaces; +namespace PMP = CGAL::Polygon_mesh_processing; + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Surface_mesh Mesh; +typedef PMP::Face_location Face_location; +typedef PMP::Edge_location Edge_location; + + +int main(int argc, char** argv) +{ + std::string filename = (argc > 1) ? std::string(argv[1]) + : CGAL::data_file_path("meshes/elephant.off"); + + Mesh mesh; + if(!CGAL::IO::read_polygon_mesh(filename, mesh) || !CGAL::is_triangle_mesh(mesh)) + { + std::cerr << "Invalid input." << std::endl; + return 1; + } + + std::size_t nb_faces = faces(mesh).size(); + + // take two random faces and pick the centroid + CGAL::Random rnd = CGAL::get_default_random(); + // CGAL::Random rnd(1695720148); + // CGAL::Random rnd(1695724381); + // CGAL::Random rnd(1695813638); + + std::cout << "seed " << rnd.get_seed() << std::endl; + Mesh::Face_index f1 = *std::next(faces(mesh).begin(), rnd.get_int(0, nb_faces)); + Mesh::Face_index f2 = *std::next(faces(mesh).begin(), rnd.get_int(0, nb_faces)); + // or pick specific faces + + // TODO: add in testsuite: special case + // Mesh::Face_index f1( 3268 ); + // Mesh::Face_index f2( 3014 ); + + // TODO: add in testsuite: special case + // Mesh::Face_index f2( 3265 ); + // Mesh::Face_index f1( 3014 ); + + // TODO: add in testsuite: special case + // Mesh::Face_index f2( 3268 ); + // Mesh::Face_index f1( 3014 ); + + // TODO: add in testsuite: special case + // Mesh::Face_index f1( 3265 ); + // Mesh::Face_index f2( 3014 ); + + // TODO: add in testsuite: special case + // Mesh::Face_index f1( 3543 ); + // Mesh::Face_index f2( 4356 ); + + Face_location src(f1, CGAL::make_array(0.3,0.3,0.4)); + Face_location tgt(f2, CGAL::make_array(0.3,0.3,0.4)); + + std::cout << "src = " << PMP::construct_point(src, mesh) << "\n"; + std::cout << "tgt = " << PMP::construct_point(tgt, mesh) << "\n"; + + std::vector edge_locations; + VGoS::locally_shortest_path(src, tgt, mesh, edge_locations); + + + std::vector poly; + poly.reserve(edge_locations.size()+2); + VGoS::convert_path_to_polyline(src, edge_locations, tgt, mesh, std::back_inserter(poly)); + + std::ofstream out("locally_shortest_path.polylines.txt"); + out << poly.size(); + for (auto p : poly) + out << " " << p; + out << "\n"; + + return 0; +} diff --git a/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/straightest_geodesic_sm_example.cpp b/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/straightest_geodesic_sm_example.cpp new file mode 100644 index 00000000000..0a7390c08bc --- /dev/null +++ b/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/straightest_geodesic_sm_example.cpp @@ -0,0 +1,65 @@ +#include +#include +#include + +#include + +namespace VGoS = CGAL::Vector_graphics_on_surfaces; +namespace PMP = CGAL::Polygon_mesh_processing; + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Surface_mesh Mesh; +typedef PMP::Face_location Face_location; +typedef PMP::Edge_location Edge_location; + + +int main(int argc, char** argv) +{ + std::string filename = (argc > 1) ? std::string(argv[1]) + : CGAL::data_file_path("meshes/elephant.off"); + + Mesh mesh; + if(!CGAL::IO::read_polygon_mesh(filename, mesh) || !CGAL::is_triangle_mesh(mesh)) + { + std::cerr << "Invalid input." << std::endl; + return 1; + } + + std::size_t nb_faces = faces(mesh).size(); + + // take two random faces and pick the centroid + CGAL::Random rnd = CGAL::get_default_random(); + //CGAL::Random rnd(1706709591); + + std::cout << "seed " << rnd.get_seed() << std::endl; + Mesh::Face_index f = *std::next(faces(mesh).begin(), rnd.get_int(0, nb_faces)); + Face_location src(f, CGAL::make_array(0.3,0.3,0.4)); + +//case opposite edge (for testsuite) +// Mesh::Face_index f = *std::next(faces(mesh).begin(), 1031); +// Face_location src(f, CGAL::make_array(0, 0.65258992669550031, 0.34741007330449963)); + + K::Point_3 src_pt = PMP::construct_point(src, mesh); + std::cout << "src = " << src_pt << "\n"; + +//case opposite edge (for testsuite) +// double target_distance = 0.0072766352463858128; +// K::Vector_2 dir(0.079665117730984697, 0.99682168366107904); + + double target_distance = 0.5; + K::Vector_2 dir(1,1); + std::vector path = VGoS::straightest_geodesic(src, dir, target_distance, mesh); + + std::vector poly; + poly.reserve(path.size()); + VGoS::convert_path_to_polyline(path, mesh, std::back_inserter(poly)); + + std::ofstream out("straightest_geodesic_path.polylines.txt"); + out << path.size() << " "; + + for (auto p : poly) + out << " " << p; + out << "\n"; + + return 0; +} diff --git a/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/trace_bezier_segment_sm_example.cpp b/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/trace_bezier_segment_sm_example.cpp new file mode 100644 index 00000000000..e585fb8b137 --- /dev/null +++ b/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/trace_bezier_segment_sm_example.cpp @@ -0,0 +1,64 @@ +#include +#include +#include + +#include + +namespace VGoS = CGAL::Vector_graphics_on_surfaces; +namespace PMP = CGAL::Polygon_mesh_processing; + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Surface_mesh Mesh; +typedef PMP::Face_location Face_location; +typedef PMP::Edge_location Edge_location; + + +int main(int argc, char** argv) +{ + std::string filename = (argc > 1) ? std::string(argv[1]) + : CGAL::data_file_path("meshes/elephant.off"); + + Mesh mesh; + if(!CGAL::IO::read_polygon_mesh(filename, mesh) || !CGAL::is_triangle_mesh(mesh)) + { + std::cerr << "Invalid input." << std::endl; + return 1; + } + + std::size_t nb_faces = faces(mesh).size(); + + // take two random faces and pick the centroid + CGAL::Random rnd = CGAL::get_default_random(); + // CGAL::Random rnd(1695795785); + + std::cout << "seed " << rnd.get_seed() << std::endl; + Mesh::Face_index f1 = *std::next(faces(mesh).begin(), rnd.get_int(0, nb_faces)); + Mesh::Face_index f2 = *std::next(faces(mesh).begin(), rnd.get_int(0, nb_faces)); + Mesh::Face_index f3 = *std::next(faces(mesh).begin(), rnd.get_int(0, nb_faces)); + Mesh::Face_index f4 = *std::next(faces(mesh).begin(), rnd.get_int(0, nb_faces)); + + + Face_location a(f1, CGAL::make_array(0.3,0.3,0.4)), + b(f2, CGAL::make_array(0.3,0.3,0.4)), + c(f3, CGAL::make_array(0.3,0.3,0.4)), + d(f4, CGAL::make_array(0.3,0.3,0.4)); + VGoS::Bezier_segment control_points=CGAL::make_array(a, b, c, d); + + std::cout << "Using the following control points:\n"; + std::cout << PMP::construct_point(a, mesh) << "\n"; + std::cout << PMP::construct_point(b, mesh) << "\n"; + std::cout << PMP::construct_point(c, mesh) << "\n"; + std::cout << PMP::construct_point(d, mesh) << "\n"; + + std::vector face_locations = + VGoS::recursive_de_Casteljau(mesh, control_points, 8); + + + std::ofstream out("bezier.polylines.txt"); + out << face_locations.size(); + for (auto fl : face_locations) + out << " " << PMP::construct_point(fl, mesh); // TODO: we should connect points with geodesics and not segments + out << "\n"; + + return 0; +} diff --git a/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/trace_polygon_example.cpp b/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/trace_polygon_example.cpp new file mode 100644 index 00000000000..f1a42ab43e1 --- /dev/null +++ b/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/trace_polygon_example.cpp @@ -0,0 +1,299 @@ +#include +#include +#include +#include +#include +#include + +#include + +#if 0 +#include +#else + +#endif + +namespace VGoS = CGAL::Vector_graphics_on_surfaces; +namespace PMP = CGAL::Polygon_mesh_processing; + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Surface_mesh Mesh; +typedef PMP::Face_location Face_location; +typedef PMP::Edge_location Edge_location; + + +int main(int argc, char** argv) +{ + std::string filename = (argc > 1) ? std::string(argv[1]) + : CGAL::data_file_path("meshes/elephant.off"); + + std::string filename_poly = (argc > 2) ? std::string(argv[2]) + : CGAL::data_file_path("XXXXX"); + + Mesh mesh; + if(!CGAL::IO::read_polygon_mesh(filename, mesh) || !CGAL::is_triangle_mesh(mesh)) + { + std::cerr << "Invalid input." << std::endl; + return 1; + } + + std::vector> polygons; + std::ifstream in(filename_poly); + if (!in) + { + std::cerr << "Error cannot open " << filename_poly << "\n"; + return 1; + } + + int nb_pt; + K::Point_3 pt; + CGAL::Bbox_2 bb2; + while (in >> nb_pt) + { + polygons.emplace_back(); + polygons.back().reserve(nb_pt-1); + for (int i=0; i> pt; + polygons.back().emplace_back(pt.x(), pt.y()); + bb2+=polygons.back().back().bbox(); + } + in >> pt; + if (!in) + { + std::cerr << "Error reading input polygons\n"; + return 1; + } + // check if last point is duplicated + if (polygons.back().back().x()!=pt.x() || polygons.back().back().y()!=pt.y()) + { + polygons.back().emplace_back(pt.x(), pt.y()); + bb2+=polygons.back().back().bbox(); + } + if (!in) break; + } + + std::cout << polygons.size() << " polygons read\n"; + + // tracing center + std::size_t nb_faces = faces(mesh).size(); + Mesh::Face_index f = *std::next(faces(mesh).begin(), (2154)%nb_faces); + Face_location center(f, CGAL::make_array(0.3,0.3,0.4)); + + // convert polygons to polar coordinates + typename K::Point_2 center_2((bb2.xmax()+bb2.xmin())/2., (bb2.ymax()+bb2.ymin())/2.); + double diag = std::sqrt( CGAL::square(bb2.xmin()-bb2.xmax()) + CGAL::square(bb2.xmin()-bb2.xmax()) ); + const double expected_diag = 0.45; // user parameter for scaling + const double scaling = expected_diag/diag; + + std::ofstream out("geodesic_polygon.polylines.txt"); + out << std::setprecision(17); + + K::Point_3 center_pt = PMP::construct_point(center, mesh); + std::cout << "center = " << center_pt << "\n"; + VGoS::Dual_geodesic_solver solver; + VGoS::init_geodesic_dual_solver(solver, mesh); + + for (const std::vector& polygon : polygons) + { + std::vector> polar_coords = + VGoS::convert_to_polar_coordinates(polygon, center_2); + if (polygon.front()==polygon.back()) polar_coords.pop_back(); + + std::vector directions; + std::vector lens; + lens.reserve(polar_coords.size()); + directions.reserve(polar_coords.size()); + + for (const std::pair& coord : polar_coords) + { + lens.push_back(scaling * coord.first); + directions.emplace_back(std::cos(coord.second), std::sin(coord.second)); + } + + // last point is duplicated + std::vector out_polygon_path = VGoS::trace_geodesic_polygon(center,directions,lens,mesh, solver); + std::vector poly; + poly.reserve(out_polygon_path.size()); + VGoS::convert_path_to_polyline(out_polygon_path, mesh, std::back_inserter(poly)); + + out << poly.size(); + for (auto p : poly) + out << " " << p; + out << std::endl; + } + + // second method + out.close(); + out.open("geodesic_polygons.polylines.txt"); + out << std::setprecision(17); + + std::vector> polygons_3 + = VGoS::trace_geodesic_polygons(center, polygons, scaling, mesh, solver); + + for (const auto& polygon : polygons_3) + { + out << polygon.size(); + for (auto p : polygon) + out << " " << PMP::construct_point(p, mesh); + out << std::endl; + } + + // third method + out.close(); + out.open("geodesic_label.polylines.txt"); + out << std::setprecision(17); + + polygons_3.clear(); + polygons_3 = VGoS::trace_geodesic_label(center, polygons, scaling, mesh, solver); + + for (const auto& polygon : polygons_3) + { + out << polygon.size(); + for (auto p : polygon) + out << " " << PMP::construct_point(p, mesh); + out << std::endl; + } + + // now refine the input mesh + std::vector cst_hedges; + auto vnm = mesh.add_property_map("vnm", K::Vector_3(0,0,0)).first; + auto fnm = mesh.add_property_map("fnm", K::Vector_3(0,0,0)).first; + using VNM = decltype(vnm); + + PMP::compute_normals(mesh, vnm, fnm); + VGoS::refine_mesh_along_paths(polygons_3, mesh, vnm, fnm, std::back_inserter(cst_hedges)); + + std::ofstream("mesh_refined.off") << std::setprecision(17) << mesh; + + std::ofstream cst_edges("refinement_edges.polylines.txt"); + cst_edges.precision(17); + for (Mesh::Halfedge_index h : cst_hedges) + cst_edges << "2 " << mesh.point(source(h,mesh)) << " " << mesh.point(target(h,mesh)) << "\n"; + + + auto ecm = mesh.add_property_map("ecm", false).first; + for (Mesh::Halfedge_index h : cst_hedges) + ecm[edge(h, mesh)]=true; + + + // face index for doing flood fill and mark inside-out + Mesh::Face_index out_face(2612); + std::vector in_out(num_faces(mesh), -1); + + bool inorout=false; + std::vector queue, next_queue; + queue.push_back(out_face); + + while(!queue.empty()) + { + Mesh::Face_index f = queue.back(); + queue.pop_back(); + if (in_out[f]==-1) + { + in_out[f]=inorout?1:0; + Mesh::Halfedge_index h=halfedge(f, mesh); + for (int i=0; i<3; ++i) + { + Mesh::Face_index nf = face(opposite(h, mesh), mesh); + if (nf!=boost::graph_traits::null_face() && in_out[nf]==-1) + { + if (ecm[edge(h,mesh)]) + next_queue.push_back(nf); + else + queue.push_back(nf); + } + h=next(h, mesh); + } + } + if (queue.empty()) + { + queue.swap(next_queue); + inorout=!inorout; + } + } + + + + struct Visitor + : public PMP::Corefinement::Default_visitor + { + VNM vnm; + Visitor(VNM vnm) : vnm(vnm) {} + + std::vector > hedge_map; + void after_edge_duplicated(Mesh::Halfedge_index h, Mesh::Halfedge_index new_hedge, const Mesh&) + { + hedge_map.emplace_back(h, new_hedge); + } + + void after_vertex_copy(Mesh::Vertex_index v, const Mesh&, Mesh::Vertex_index nv, const Mesh&) + { + put(vnm, nv, get(vnm, v)); + } + }; + Visitor visitor(vnm); + + PMP::internal::split_along_edges(mesh, ecm, mesh.points(), visitor); + + + double delta = -0.005; + for (const auto& ph : visitor.hedge_map) + { + Mesh::Halfedge_index h1 = ph.first; + Mesh::Halfedge_index h2 = ph.second; + if (is_border(h1, mesh)) h1=opposite(h1, mesh); + if (is_border(h2, mesh)) h2=opposite(h2, mesh); + Mesh::Halfedge_index h = in_out[face(h1, mesh)]==1 ? h1 : h2; + + Mesh::Vertex_index v = target(h, mesh); + K::Vector_3 n = get(vnm, v); + mesh.point(v) = mesh.point(v)+delta*n; + } + + // interior vertices + for (Mesh::Vertex_index v : vertices(mesh)) + { + bool skip=false; + Mesh::Halfedge_index h = halfedge(v, mesh); + for (Mesh::Halfedge_index h : CGAL::halfedges_around_target(v, mesh)) + { + if (is_border(h, mesh) || in_out[face(h, mesh)]==0) + { + skip=true; + break; + } + } + if (!skip) + { + K::Vector_3 n = get(vnm, v); + mesh.point(v) = mesh.point(v)+delta*n; + } + } + + std::vector b1(visitor.hedge_map.size()); + std::vector b2(visitor.hedge_map.size()); + for (std::size_t i=0; i +#include +#include + +#include + +namespace VGoS = CGAL::Vector_graphics_on_surfaces; +namespace PMP = CGAL::Polygon_mesh_processing; + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Surface_mesh Mesh; +typedef PMP::Face_location Face_location; +typedef PMP::Edge_location Edge_location; + +std::vector +get_supporting_curve(std::string filename, + const Mesh& mesh) +{ + std::ifstream in(filename); + if (!in) + { + std::cerr << "ERROR: cannot open " << filename << "\n"; + exit(1); + } + + int nbp; + in >> nbp; + std::cout << "Support polyline size: " << nbp << "\n"; + std::vector polyline(nbp); + + for (int i=0;i> polyline[i]; + + std::vector face_locations(nbp); + + using AABB_face_graph_primitive = CGAL::AABB_face_graph_triangle_primitive; + using AABB_face_graph_traits = CGAL::AABB_traits_3; + + CGAL::AABB_tree tree; + PMP::build_AABB_tree(mesh, tree); + + //Remarks on snapping_tolerance of locate: it is actually a value used as is in barycentric coordinates with no relation to the size of the triangle + // and snapping is only done when at least one coordinate is negative (clamping would be a better name, snapping is what I'm doing below). + //TODO: the following works only if we don't have a lfs smaller than snap_tol + double snap_tol = 0.00001; + for(int i=0; i primitives; + tree.all_intersected_primitives(query, std::back_inserter(primitives)); + + switch (primitives.size()) + { + case 0: + std::cerr <<"ERROR points too far!\n"; + exit(1); + break; + case 1: + face_locations[i] = PMP::locate_in_face(polyline[i],primitives[0], mesh); + break; + case 2: + { + // TODO: we assume we don't have any boundary to it's clear it's on an edge + Mesh::Halfedge_index h1 = halfedge(primitives[0], mesh), + h2 = halfedge(primitives[1], mesh); + bool found=false; + for (int j=0;j<3;++j) + { + for (int k=0;k<3; ++k) + { + if (h1==opposite(h2, mesh)) + { + found = true; + break; + } + h2=next(h2, mesh); + } + if (found) break; + h1=next(h1, mesh); + } + CGAL_assertion(h1 == opposite(h2, mesh)); + + const K::Point_3& src = mesh.point(source(h1, mesh)); + const K::Point_3& tgt= mesh.point(target(h1, mesh)); + double t = std::sqrt(CGAL::squared_distance(polyline[i], src)/CGAL::squared_distance(src,tgt)); + face_locations[i] = PMP::locate_on_halfedge(h1, t, mesh); + } + break; + default: + { + std::vector vrts; + Mesh::Halfedge_index h=halfedge(primitives[0], mesh); + vrts.push_back(source(h, mesh)); + vrts.push_back(target(h, mesh)); + vrts.push_back(target(next(h, mesh), mesh)); + std::sort(vrts.begin(), vrts.end()); + for(std::size_t k=1;k tmp; + h=halfedge(primitives[k], mesh); + for (int j=0;j<3;++j) + { + Mesh::Vertex_index v=target(h, mesh); + if (std::binary_search(vrts.begin(), vrts.end(), v)) + tmp.push_back(v); + h=next(h, mesh); + } + tmp.swap(vrts); + std::sort(vrts.begin(), vrts.end()); + } + CGAL_assertion(vrts.size()==1); + + int offset=0; + h=halfedge(primitives[0], mesh); + if (target(h, mesh)==vrts[0]) + offset=1; + else + if (target(next(h, mesh), mesh)==vrts[0]) + offset=2; + else + CGAL_assertion(source(h, mesh)==vrts[0]); + std::array bary = CGAL::make_array(0.,0.,0.); + bary[offset]=1.; + face_locations[i] = std::make_pair(primitives[0], bary); + } + } + // face_locations[i] = PMP::locate_with_AABB_tree(polyline[i], tree, mesh, CGAL::parameters::snapping_tolerance(snap_tol)); + for (int k=0;k<3;++k) + if (face_locations[i].second[k]<0) face_locations[i].second[k]=0; + } + + // DEBUG code to check that the polyline has correctly been converted to face locations (it's not a path as if vertex is hit we don't have the face continuity property) + for (int i=0;i 1) ? std::string(argv[1]) + : CGAL::data_file_path("meshes/elephant.off"); + + std::string support_filename = (argc > 2) ? std::string(argv[2]) + : CGAL::data_file_path("XXXXX"); + + std::string filename_poly = (argc > 3) ? std::string(argv[3]) + : CGAL::data_file_path("XXXXX"); + + std::size_t nb_copies = atoi(argv[4]); + + Mesh mesh; + if(!CGAL::IO::read_polygon_mesh(filename, mesh) || !CGAL::is_triangle_mesh(mesh)) + { + std::cerr << "Invalid input." << std::endl; + return 1; + } + + std::vector> polygons; + std::ifstream in(filename_poly); + if (!in) + { + std::cerr << "Error cannot open " << filename_poly << "\n"; + return 1; + } + + int nb_pt; + K::Point_3 pt; + CGAL::Bbox_2 bb2; + while (in >> nb_pt) + { + polygons.emplace_back(); + polygons.back().reserve(nb_pt-1); + for (int i=0; i> pt; + polygons.back().emplace_back(pt.x(), pt.y()); + bb2+=polygons.back().back().bbox(); + } + in >> pt; + if (!in) + { + std::cerr << "Error reading input polygons\n"; + return 1; + } + // check if last point is duplicated + if (polygons.back().back().x()!=pt.x() || polygons.back().back().y()!=pt.y()) + { + polygons.back().emplace_back(pt.x(), pt.y()); + bb2+=polygons.back().back().bbox(); + } + if (!in) break; + } + + std::cout << polygons.size() << " polygons read\n"; + + + // duplicate them a bit to get longer text while getting another one + CGAL::Bbox_2 gbox; + for (const auto& polygon : polygons) + gbox+=CGAL::bbox_2(polygon.begin(), polygon.end()); + K::Vector_2 delta(gbox.xmax()-gbox.xmin(), 0); + std::size_t nbpoly=polygons.size(); + polygons.reserve(nbpoly*(nb_copies+1)); + for (std::size_t c=0; c solver; + VGoS::init_geodesic_dual_solver(solver, mesh); + + // get supporting curve + std::vector supporting_curve = get_supporting_curve(support_filename, mesh); + if (supporting_curve.empty()) return 1; + + + std::ofstream support_out("support.polylines.txt"); + + std::vector support_poly; + support_poly.reserve(supporting_curve.size()); + VGoS::convert_path_to_polyline(supporting_curve, mesh, std::back_inserter(support_poly)); + + support_out << std::setprecision(17) << support_poly.size(); + for (auto p : support_poly) + support_out << " " << p; + support_out << "\n"; + support_out.close(); + std::cout <<"supporting_curve generated!\n"; + + // convert polygons to polar coordinates + const double expected_height = 0.025; // user parameter for scaling + const double scaling = expected_height/(gbox.ymax()-gbox.ymin()); + + + std::ofstream out("label_on_curve.polylines.txt"); + out << std::setprecision(17); + + std::vector> polygons_3; + polygons_3 = VGoS::trace_geodesic_label_along_curve(supporting_curve, polygons, scaling, 0., false, mesh, solver); + + for (const auto& polygon_path : polygons_3) + { + std::vector poly; + poly.reserve(polygon_path.size()); + VGoS::convert_path_to_polyline(polygon_path, mesh, std::back_inserter(poly)); + + out << poly.size(); + for (auto p : poly) + out << " " << p; + out << std::endl; + } + + return 0; +} diff --git a/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/trace_polygon_on_svg_curve_example.cpp b/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/trace_polygon_on_svg_curve_example.cpp new file mode 100644 index 00000000000..5acab114bf1 --- /dev/null +++ b/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/trace_polygon_on_svg_curve_example.cpp @@ -0,0 +1,231 @@ +#include +#include +#include + +#include + +#include + +namespace VGoS = CGAL::Vector_graphics_on_surfaces; +namespace PMP = CGAL::Polygon_mesh_processing; + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Surface_mesh Mesh; +typedef PMP::Face_location Face_location; +typedef PMP::Edge_location Edge_location; + +std::vector +get_supporting_curve(std::string svg_filename, + const Mesh& mesh, + Face_location center, const VGoS::Dual_geodesic_solver& solver) +{ + std::vector res; + + NSVGimage* g_image = nsvgParseFromFile(svg_filename.c_str(), "px", 96.0f); + if (g_image == NULL) { + printf("Could not open SVG image.\n"); + return res; + } + + // extract control points + std::vector< std::array > bezier_curves; + CGAL::Bbox_2 bb2; + + // in SVG's the y axis points downward, so we must take the opposite y coordinates + for (NSVGshape* shape = g_image->shapes; shape != NULL; shape = shape->next) + { + for (NSVGpath* path = shape->paths; path != NULL; path = path->next) + { + CGAL::Bbox_2 path_bbox(path->bounds[0], -path->bounds[1], + path->bounds[2], -path->bounds[3]); + bb2+=path_bbox; + + float* pts=path->pts; + int npts=path->npts; + + for (int i=0; i> directions; + std::vector> lengths; + directions.reserve(bezier_curves.size()); + lengths.reserve(bezier_curves.size()); + + CGAL_assertion_code(K::Point_2 prev = bezier_curves[0][0];) + for (const std::array& bezier : bezier_curves) + { + CGAL_assertion(bezier[0]==prev); + CGAL_assertion_code(prev=bezier[3];) + std::vector> polar_coords = + VGoS::convert_to_polar_coordinates(bezier, center_2); + + int start_id=directions.empty()?0:1; + + directions.emplace_back(); + lengths.emplace_back(); + + assert(polar_coords.size()==4); + + for (int i=start_id;i<4; ++i) + { + lengths.back()[i] = scaling * polar_coords[i].first; + directions.back()[i]=K::Vector_2(std::cos(polar_coords[i].second), std::sin(polar_coords[i].second)); + } + if (start_id==1) + { + lengths.back()[0] = lengths[lengths.size()-2][3]; + directions.back()[0] = directions[directions.size()-2][3]; + } + } + + // trace bezier curves + std::vector< std::vector > resi = + VGoS::trace_bezier_curves(center, directions, lengths, 6, mesh, solver); + + //TODO: here we assume that curves are parameterized in the same order and are consecutives + for (const std::vector& r : resi) + res.insert(res.end(), r.begin(), std::prev(r.end())); + res.push_back(resi.back().back()); + + std::reverse(res.begin(), res.end()); // TODO: should be an option! + + return res; +} + +int main(int argc, char** argv) +{ + std::string filename = (argc > 1) ? std::string(argv[1]) + : CGAL::data_file_path("meshes/elephant.off"); + + std::string svg_filename = (argc > 2) ? std::string(argv[2]) + : CGAL::data_file_path("polylines_2/Archimedean_spiral.svg"); + + std::string filename_poly = (argc > 3) ? std::string(argv[3]) + : CGAL::data_file_path("XXXXX"); + + Mesh mesh; + if(!CGAL::IO::read_polygon_mesh(filename, mesh) || !CGAL::is_triangle_mesh(mesh)) + { + std::cerr << "Invalid input." << std::endl; + return 1; + } + + std::vector> polygons; + std::ifstream in(filename_poly); + if (!in) + { + std::cerr << "Error cannot open " << filename_poly << "\n"; + return 1; + } + + int nb_pt; + K::Point_3 pt; + CGAL::Bbox_2 bb2; + while (in >> nb_pt) + { + polygons.emplace_back(); + polygons.back().reserve(nb_pt-1); + for (int i=0; i> pt; + polygons.back().emplace_back(pt.x(), pt.y()); + bb2+=polygons.back().back().bbox(); + } + in >> pt; + if (!in) + { + std::cerr << "Error reading input polygons\n"; + return 1; + } + // check if last point is duplicated + if (polygons.back().back().x()!=pt.x() || polygons.back().back().y()!=pt.y()) + { + polygons.back().emplace_back(pt.x(), pt.y()); + bb2+=polygons.back().back().bbox(); + } + if (!in) break; + } + + std::cout << polygons.size() << " polygons read\n"; + + // tracing center + std::size_t nb_faces = faces(mesh).size(); + Mesh::Face_index f = *std::next(faces(mesh).begin(), (2154)%nb_faces); + Face_location center(f, CGAL::make_array(0.3,0.3,0.4)); + + K::Point_3 center_pt = PMP::construct_point(center, mesh); + std::cout << "center = " << center_pt << "\n"; + VGoS::Dual_geodesic_solver solver; + VGoS::init_geodesic_dual_solver(solver, mesh); + + + // get supporting curve + std::vector supporting_curve = get_supporting_curve(svg_filename, mesh, center, solver); + if (supporting_curve.empty()) return 1; + + std::cout <<"supporting_curve generated!\n"; + std::ofstream support_out("svg_support.polylines.txt"); + + std::vector support_poly; + support_poly.reserve(supporting_curve.size()); + VGoS::convert_path_to_polyline(supporting_curve, mesh, std::back_inserter(support_poly)); + + support_out << std::setprecision(17) << support_poly.size(); + for (auto p : support_poly) + support_out << " " << p; + support_out << "\n"; + support_out.close(); + + // convert polygons to polar coordinates + typename K::Point_2 center_2((bb2.xmax()+bb2.xmin())/2., (bb2.ymax()+bb2.ymin())/2.); + double diag = std::sqrt( CGAL::square(bb2.xmin()-bb2.xmax()) + CGAL::square(bb2.xmin()-bb2.xmax()) ); + const double expected_diag = 2.1; // user parameter for scaling + const double scaling = expected_diag/diag; + + + std::ofstream out("label_on_curve.polylines.txt"); + out << std::setprecision(17); + + std::vector> polygons_3; + polygons_3 = VGoS::trace_geodesic_label_along_curve(supporting_curve, polygons, scaling, 0., true, mesh, solver); + + for (const auto& polygon_path : polygons_3) + { + std::vector poly; + poly.reserve(polygon_path.size()); + VGoS::convert_path_to_polyline(polygon_path, mesh, std::back_inserter(poly)); + + out << poly.size(); + for (auto p : poly) + out << " " << p; + out << std::endl; + } + + return 0; +} diff --git a/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/trace_regular_polygons_example.cpp b/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/trace_regular_polygons_example.cpp new file mode 100644 index 00000000000..b85ef4efd8e --- /dev/null +++ b/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/trace_regular_polygons_example.cpp @@ -0,0 +1,76 @@ +#include +#include +#include + +#include + +#if 0 +#include +#else + +#endif + +namespace VGoS = CGAL::Vector_graphics_on_surfaces; +namespace PMP = CGAL::Polygon_mesh_processing; + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Surface_mesh Mesh; +typedef PMP::Face_location Face_location; +typedef PMP::Edge_location Edge_location; + + +int main(int argc, char** argv) +{ + std::string filename = (argc > 1) ? std::string(argv[1]) + : CGAL::data_file_path("meshes/elephant.off"); + + Mesh mesh; + if(!CGAL::IO::read_polygon_mesh(filename, mesh) || !CGAL::is_triangle_mesh(mesh)) + { + std::cerr << "Invalid input." << std::endl; + return 1; + } + + std::size_t nb_faces = faces(mesh).size(); + int n_sides=6; + std::vector lenghts = { 0.03, 0.02, 0.01, 0.005 }; + + // take two random faces and pick the centroid + CGAL::Random rnd = CGAL::get_default_random(); + //CGAL::Random rnd(1707129825); + + std::cout << "seed " << rnd.get_seed() << std::endl; + Mesh::Face_index f = *std::next(faces(mesh).begin(), rnd.get_int(0, nb_faces)); + + Face_location center(f, CGAL::make_array(0.3,0.3,0.4)); + std::vector directions(n_sides); + double step=2*CGAL_PI/n_sides; + for(int i=0;i solver; + VGoS::init_geodesic_dual_solver(solver, mesh); + for (double len : lenghts) + { + std::vector lens(n_sides,len); + std::vector polygon_path = VGoS::trace_geodesic_polygon(center,directions,lens,mesh, solver); + + std::vector poly; + poly.reserve(polygon_path.size()); + VGoS::convert_path_to_polyline(polygon_path, mesh, std::back_inserter(poly)); + + out << poly.size(); + for (auto p : poly) + out << " " << p; + out << "\n"; + } + + return 0; +} diff --git a/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/trace_svg_example.cpp b/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/trace_svg_example.cpp new file mode 100644 index 00000000000..bf45a889ec1 --- /dev/null +++ b/Vector_graphics_on_surfaces/examples/Vector_graphics_on_surfaces/trace_svg_example.cpp @@ -0,0 +1,502 @@ +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include + +namespace VGoS = CGAL::Vector_graphics_on_surfaces; +namespace PMP = CGAL::Polygon_mesh_processing; + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Surface_mesh Mesh; +typedef PMP::Face_location Face_location; +typedef PMP::Edge_location Edge_location; + + +int main(int argc, char** argv) +{ + std::string mesh_filename = (argc > 1) ? std::string(argv[1]) + : CGAL::data_file_path("meshes/elephant.off"); + + std::string svg_filename = (argc > 2) ? std::string(argv[2]) + : CGAL::data_file_path("polylines_2/nano.svg"); + + const double expected_diag = (argc>3) ? atof(argv[3]) : 0.6; // user parameter for scaling + int face_index = (argc>4) ? atoi(argv[4]) : 2154; + double b0 = (argc>5) ? atof(argv[5]) : 0.3; + double b1 = (argc>6) ? atof(argv[6]) : 0.3; + double b2 = (argc>7) ? atof(argv[7]) : 0.4; + + Mesh mesh; + if(!CGAL::IO::read_polygon_mesh(mesh_filename, mesh) || !CGAL::is_triangle_mesh(mesh)) + { + std::cerr << "Invalid input mesh." << std::endl; + return 1; + } + + NSVGimage* g_image = nsvgParseFromFile(svg_filename.c_str(), "px", 96.0f); + if (g_image == NULL) { + printf("Could not open SVG image.\n"); + return 1; + } + + // extract control points + std::vector< std::array > bezier_control_points; + CGAL::Bbox_2 bb2; + + // in SVG's the y axis points downward, so we must take the opposite y coordinates + for (NSVGshape* shape = g_image->shapes; shape != NULL; shape = shape->next) + { + for (NSVGpath* path = shape->paths; path != NULL; path = path->next) + { + CGAL::Bbox_2 path_bbox(path->bounds[0], -path->bounds[1], + path->bounds[2], -path->bounds[3]); + bb2+=path_bbox; + + float* pts=path->pts; + int npts=path->npts; + + for (int i=0; i solver; + VGoS::init_geodesic_dual_solver(solver, mesh); + + std::size_t nb_faces = faces(mesh).size(); + Mesh::Face_index f = *std::next(faces(mesh).begin(), (face_index)%nb_faces); + Face_location center(f, CGAL::make_array(b0,b1,b2)); + + std::size_t num_subdiv = 3; + + // option 1: trace bezier curves with a single center + { + CGAL::Real_timer time; + time.start(); + std::vector> directions; + std::vector> lengths; + directions.reserve(bezier_control_points.size()); + lengths.reserve(bezier_control_points.size()); + + for (const std::array& bezier : bezier_control_points) + { + std::vector> polar_coords = + VGoS::convert_to_polar_coordinates(bezier, center_2); + + directions.emplace_back(); + lengths.emplace_back(); + + assert(polar_coords.size()==4); + + for (int i=0;i<4; ++i) + { + lengths.back()[i] = scaling * polar_coords[i].first; + directions.back()[i]=K::Vector_2(std::cos(polar_coords[i].second), std::sin(polar_coords[i].second)); + } + } + + std::vector< std::vector > res = + VGoS::trace_bezier_curves(center, directions, lengths, num_subdiv, mesh, solver); + + // write result + std::ofstream out("svg_option1.polylines.txt"); + out << std::setprecision(17); + for (const auto& b : res) + { + std::vector poly; + poly.reserve(b.size()); + VGoS::convert_path_to_polyline(b, mesh, std::back_inserter(poly)); + + + out << poly.size(); + for (const K::Point_3& p : poly) + out << " " << p; + out << "\n"; + } + time.stop(); + std::cout << "option 1 took: " << time.time() << "\n"; + } + + // option 2: regroup control polygon and trace bezier curves with a center per group + { + CGAL::Real_timer time; + time.start(); + + typedef boost::adjacency_list< boost::vecS, boost::vecS, + boost::undirectedS, + K::Point_2> Graph; + using Graph_vertex = Graph::vertex_descriptor; + + Graph graph; + + std::vector vrts; + std::map pt_map; + Graph_vertex null_vertex = boost::graph_traits::null_vertex(); + + // build bezier control segment graph + for (std::size_t eid=0; eid& bezier = bezier_control_points[eid]; + + auto insert_res_0 = pt_map.emplace(bezier[0], null_vertex); + if (insert_res_0.second) + { + vrts.push_back(add_vertex(graph)); + graph[vrts.back()]=bezier[0]; + insert_res_0.first->second = vrts.back(); + } + + auto insert_res_3 = pt_map.emplace(bezier[3], null_vertex); + if (insert_res_3.second) + { + vrts.push_back(add_vertex(graph)); + graph[vrts.back()]=bezier[3]; + insert_res_3.first->second = vrts.back(); + } + + Graph_vertex v1 = add_vertex(graph); + Graph_vertex v2 = add_vertex(graph); + graph[v1]=bezier[1]; + graph[v2]=bezier[2]; + + add_edge(insert_res_0.first->second, v1, graph); + add_edge(v1, v2, graph); + add_edge(v2, insert_res_3.first->second, graph); + } + + // visitor for split_graph_into_polylines that will also take care of + // the drawing and writing in the output file + struct Collect_visitor + { + const Graph& graph; + std::vector>& polylines_2; + std::vector current_polyline; + + Collect_visitor(const Graph& g, std::vector>& poly2) + : graph(g) + , polylines_2(poly2) + {} + + void start_new_polyline() + { + current_polyline.clear(); + } + + void add_node(Graph_vertex v) + { + current_polyline.push_back(v); + } + + void end_polyline() + { + // recover polyline of control points + polylines_2.emplace_back(); + polylines_2.back().reserve(current_polyline.size()); + for (Graph_vertex v : current_polyline) + polylines_2.back().push_back(graph[v]); + } + }; + + + std::vector> polylines_2; + Collect_visitor svisitor(graph, polylines_2); + CGAL::split_graph_into_polylines(graph, svisitor); + std::size_t nb_poly = polylines_2.size(); + + // regroup polylines if they have overlapping bboxes + std::vector bboxes_2(nb_poly); + using UF = CGAL::Union_find; + std::vector uf_handles(nb_poly); + UF uf; + for(std::size_t i=0;i partition_id(nb_poly,-1); + std::vector< std::vector > partition; + partition.reserve(uf.number_of_sets()); + + for (std::size_t i=0; i> polygons_3; + std::ofstream out("svg_option2.polylines.txt"); + out << std::setprecision(17); + bool all_closed=true; + for (const std::vector& pids : partition) + { + // center of the polyline + CGAL::Bbox_2 poly_bb; + for (std::size_t i : pids) + poly_bb+=bboxes_2[i]; + K::Point_2 poly_center((poly_bb.xmin()+poly_bb.xmax())/2., (poly_bb.ymin()+poly_bb.ymax())/2.); + + // path from the general center to the polyline center + K::Vector_2 dir(center_2, poly_center); + std::vector path_2_center = + VGoS::straightest_geodesic(center, dir, scaling * std::sqrt(dir.squared_length()), mesh); + Face_location poly_center_loc = path_2_center.back(); + + // update initial angle from the global center to the polyline center + using Impl=VGoS::internal::Locally_shortest_path_imp; + K::Vector_2 initial_dir(1,0);// TODO: this is arbitray and could be a user input or from PCA... + + + // TODO: avoid using the shortest path and directly use the straightest! + std::vector shortest_path; + VGoS::locally_shortest_path(center, poly_center_loc, mesh, shortest_path, solver); + + typename K::Vector_2 v = initial_dir; + // TODO: the code uses the straightest path but since the code expects Edge_location, + // it is not working directly. We need to extract the edge encoded by the face location + // for(std::size_t i=1;i directions; + std::vector lengths; + directions.reserve(polylines_2[pid].size()); + lengths.reserve(polylines_2[pid].size()); + + // polar coordinates from the polyline center and convertion to lengths and directions + std::vector> polar_coords = + VGoS::convert_to_polar_coordinates(polylines_2[pid], poly_center); + + for (const std::pair& polar_coord : polar_coords) + { + lengths.push_back(scaling * polar_coord.first); + directions.push_back(K::Vector_2(std::cos(polar_coord.second+theta), std::sin(polar_coord.second+theta))); + } + + bool is_closed=polylines_2[pid].front()==polylines_2[pid].back(); + if (is_closed) + { + lengths.pop_back(); + directions.pop_back(); + } + else + all_closed=false; + + + std::vector res = + VGoS::trace_polyline_of_bezier_curves(poly_center_loc, directions, lengths, + is_closed, // use [directions/lengths].front as last control point? + num_subdiv, mesh, solver); + // write result + std::vector poly3; + poly3.reserve(res.size()); + VGoS::convert_path_to_polyline(res, mesh, std::back_inserter(poly3)); + + out << poly3.size(); + for (const K::Point_3& p : poly3) + out << " " << p; + out << "\n"; + + polygons_3.push_back(std::move(res)); + } + } + out.close(); + + time.stop(); + std::cout << "option 2 took: " << time.time() << "\n"; + + + if (all_closed) + { + std::cout << "Now carve the input mesh\n"; + // copy/pasted from trace_polygon_example.cpp + // now refine the input mesh + std::vector cst_hedges; + auto vnm = mesh.add_property_map("vnm", K::Vector_3(0,0,0)).first; + auto fnm = mesh.add_property_map("fnm", K::Vector_3(0,0,0)).first; + using VNM = decltype(vnm); + + PMP::compute_normals(mesh, vnm, fnm); + VGoS::refine_mesh_along_paths(polygons_3, mesh, vnm, fnm, std::back_inserter(cst_hedges)); + + std::ofstream("mesh_refined.off") << std::setprecision(17) << mesh; + + std::ofstream cst_edges("refinement_edges.polylines.txt"); + cst_edges.precision(17); + for (Mesh::Halfedge_index h : cst_hedges) + cst_edges << "2 " << mesh.point(source(h,mesh)) << " " << mesh.point(target(h,mesh)) << "\n"; + + + auto ecm = mesh.add_property_map("ecm", false).first; + for (Mesh::Halfedge_index h : cst_hedges) + ecm[edge(h, mesh)]=true; + + + // face index for doing flood fill and mark inside-out + Mesh::Face_index out_face(2612); + std::vector in_out(num_faces(mesh), -1); + + bool inorout=false; + std::vector queue, next_queue; + queue.push_back(out_face); + + while(!queue.empty()) + { + Mesh::Face_index f = queue.back(); + queue.pop_back(); + if (in_out[f]==-1) + { + in_out[f]=inorout?1:0; + Mesh::Halfedge_index h=halfedge(f, mesh); + for (int i=0; i<3; ++i) + { + Mesh::Face_index nf = face(opposite(h, mesh), mesh); + if (nf!=boost::graph_traits::null_face() && in_out[nf]==-1) + { + if (ecm[edge(h,mesh)]) + next_queue.push_back(nf); + else + queue.push_back(nf); + } + h=next(h, mesh); + } + } + if (queue.empty()) + { + queue.swap(next_queue); + inorout=!inorout; + } + } + + struct Visitor + : public PMP::Corefinement::Default_visitor + { + VNM vnm; + Visitor(VNM vnm) : vnm(vnm) {} + + std::vector > hedge_map; + void after_edge_duplicated(Mesh::Halfedge_index h, Mesh::Halfedge_index new_hedge, const Mesh&) + { + hedge_map.emplace_back(h, new_hedge); + } + + void after_vertex_copy(Mesh::Vertex_index v, const Mesh&, Mesh::Vertex_index nv, const Mesh&) + { + put(vnm, nv, get(vnm, v)); + } + }; + Visitor visitor(vnm); + + PMP::internal::split_along_edges(mesh, ecm, mesh.points(), visitor); + + double delta = -0.0005; + for (const auto& ph : visitor.hedge_map) + { + Mesh::Halfedge_index h1 = ph.first; + Mesh::Halfedge_index h2 = ph.second; + if (is_border(h1, mesh)) h1=opposite(h1, mesh); + if (is_border(h2, mesh)) h2=opposite(h2, mesh); + Mesh::Halfedge_index h = in_out[face(h1, mesh)]==1 ? h1 : h2; + + Mesh::Vertex_index v = target(h, mesh); + K::Vector_3 n = get(vnm, v); + mesh.point(v) = mesh.point(v)+delta*n; + } + + // interior vertices + for (Mesh::Vertex_index v : vertices(mesh)) + { + bool skip=false; + Mesh::Halfedge_index h = halfedge(v, mesh); + for (Mesh::Halfedge_index h : CGAL::halfedges_around_target(v, mesh)) + { + if (is_border(h, mesh) || in_out[face(h, mesh)]==0) + { + skip=true; + break; + } + } + if (!skip) + { + K::Vector_3 n = get(vnm, v); + mesh.point(v) = mesh.point(v)+delta*n; + } + } + + std::vector b1(visitor.hedge_map.size()); + std::vector b2(visitor.hedge_map.size()); + for (std::size_t i=0; i + +#include +#include +#include +#include + +#include +//TODO: split to avoid redundant linking +#include +#include + +#include +#include +#include +#include +#include + +#ifdef CGAL_DEBUG_BSURF +#include +#endif + +namespace CGAL { +namespace Vector_graphics_on_surfaces { + +#ifndef CGAL_BSURF_USE_DIJKSTRA_SP +template +struct Dual_geodesic_solver; +#endif + +template +void locally_shortest_path(CGAL::Polygon_mesh_processing::Face_location src, + CGAL::Polygon_mesh_processing::Face_location tgt, + const TriangleMesh &tmesh, + EdgeLocationRange &edge_locations +#ifndef CGAL_BSURF_USE_DIJKSTRA_SP + , const Dual_geodesic_solver& solver = Dual_geodesic_solver() +#endif +); + +/*! + * \ingroup VGSMiscellaneous + * array containing the locations of the four control points of a Bézier segment on a triangle mesh + * `init_geodesic_dual_solver()`. + * \tparam FT floating point number type (float or double) + * \tparam TriangleMesh a model of `FaceListGraph` and `EdgeListGraph` + */ +template +using Bezier_segment = std::array, 4>; + +namespace internal { + +template +struct Locally_shortest_path_imp +{ + using face_descriptor = + typename boost::graph_traits::face_descriptor; + using vertex_descriptor = + typename boost::graph_traits::vertex_descriptor; + using halfedge_descriptor = + typename boost::graph_traits::halfedge_descriptor; + + using Point_2 = typename K::Point_2; + using Point_3 = typename K::Point_3; + using Vector_2 = typename K::Vector_2; + using Vector_3 = typename K::Vector_3; + using FT = typename K::FT; + +#ifdef CGAL_DEBUG_BSURF + static + void dump_path(const std::vector& path, + const std::vector& lerps, + const CGAL::Polygon_mesh_processing::Face_location& src, + const CGAL::Polygon_mesh_processing::Face_location& tgt, + const TriangleMesh& mesh) + { + namespace PMP = CGAL::Polygon_mesh_processing; + static int i = -1; + std::cout << "dump current path in path_"+std::to_string(i+1)+".polylines.txt\n"; + std::ofstream out("path_"+std::to_string(++i)+".polylines.txt"); + out << path.size()+2 << " " << PMP::construct_point(src, mesh); + for(std::size_t i=0; i(path[i], make_array(lerps[i], 1.-lerps[i])), mesh); + out << " " << PMP::construct_point(tgt, mesh) << "\n"; + } +#endif + + // TODO: recode using CGAL code? + static + Vector_2 + intersect_circles(const Vector_2 &c2, FT R2, + const Vector_2 &c1, FT R1) + { + auto R = (c2 - c1).squared_length(); + assert(R > 0); + auto invR = FT(1) / R; + Vector_2 result = c1+c2; + + result = result + (c2 - c1) * ((R1 - R2) * invR); + auto A = 2 * (R1 + R2) * invR; + auto B = (R1 - R2) * invR; + auto s = A - B * B - 1; + assert(s >= 0); + result = result + Vector_2(c2.y() - c1.y(), c1.x() - c2.x()) * sqrt(s); + return result / 2.; + } + + static std::array + init_flat_triangle( const halfedge_descriptor& h, + const VertexPointMap &vpm, const TriangleMesh &mesh) + { + std::array triangle_vertices = make_array( + source(h, mesh), target(h, mesh), target(next(h, mesh), mesh)); + + std::array tr2d; + tr2d[0] = Vector_2(0, 0); + tr2d[1] = + Vector_2(0, sqrt(squared_distance(get(vpm, triangle_vertices[0]), + get(vpm, triangle_vertices[1])))); + auto rx = squared_distance(get(vpm, triangle_vertices[0]), + get(vpm, triangle_vertices[2])); + auto ry = squared_distance(get(vpm, triangle_vertices[1]), + get(vpm, triangle_vertices[2])); + tr2d[2] = intersect_circles(tr2d[0], rx, tr2d[1], ry); + + return tr2d; + } + + static std::array + init_source_triangle(halfedge_descriptor hopp, + const VertexPointMap &vpm, + const TriangleMesh &mesh, + CGAL::Polygon_mesh_processing::Face_location src) + { + halfedge_descriptor h = opposite(hopp, mesh); + std::array triangle_vertices = make_array( + source(h, mesh), target(h, mesh), target(next(h, mesh), mesh)); + + std::array tr2d; + tr2d[0] = Vector_2(0, 0); + tr2d[1] = + Vector_2(0, sqrt(squared_distance(get(vpm, triangle_vertices[0]), + get(vpm, triangle_vertices[1])))); + FT rx = squared_distance(get(vpm, triangle_vertices[0]), + get(vpm, triangle_vertices[2])); + FT ry = squared_distance(get(vpm, triangle_vertices[1]), + get(vpm, triangle_vertices[2])); + tr2d[2] = intersect_circles(tr2d[0], rx, tr2d[1], ry); + + + halfedge_descriptor href = halfedge(src.first, mesh); + if (href!=h) + { + if (href==next(h, mesh)) + { + std::array tmp = CGAL::make_array(src.second[2], src.second[0], src.second[1]); + src.second = tmp; + } + else + { + CGAL_assertion(next(href, mesh)==h); + std::array tmp = CGAL::make_array(src.second[1], src.second[2], src.second[0]); + src.second = tmp; + } + } + + auto point_coords = tr2d[0] * src.second[0] + tr2d[1] * src.second[1] + + tr2d[2] * src.second[2]; +#ifdef CGAL_DEBUG_BSURF + std::cout << "4 " << tr2d[0] - point_coords << " 0 " + << tr2d[1] - point_coords << " 0 " + << tr2d[2] - point_coords << " 0 " + << tr2d[0] - point_coords << " 0\n"; +#endif + return make_array(tr2d[0] - point_coords, + tr2d[1] - point_coords); + } + + static std::array + init_target_triangle(halfedge_descriptor h, + const std::array& flat_tid, + const VertexPointMap &vpm, const TriangleMesh &mesh, + CGAL::Polygon_mesh_processing::Face_location tgt) + { + std::array triangle_vertices = make_array( + source(h, mesh), target(h, mesh), target(next(h, mesh), mesh)); + + std::array tr2d; + tr2d[0] = flat_tid[1]; + tr2d[1] = flat_tid[0]; + FT rx = squared_distance(get(vpm, triangle_vertices[0]), + get(vpm, triangle_vertices[2])); + FT ry = squared_distance(get(vpm, triangle_vertices[1]), + get(vpm, triangle_vertices[2])); + tr2d[2] = intersect_circles(tr2d[0], rx, tr2d[1], ry); + + + halfedge_descriptor href = halfedge(tgt.first, mesh); + if (href!=h) + { + if (href==next(h, mesh)) + { + std::array tmp = CGAL::make_array(tgt.second[2], tgt.second[0], tgt.second[1]); + tgt.second = tmp; + } + else + { + CGAL_assertion(next(href, mesh)==h); + std::array tmp = CGAL::make_array(tgt.second[1], tgt.second[2], tgt.second[0]); + tgt.second = tmp; + } + } + + auto point_coords = tr2d[0] * tgt.second[0] + tr2d[1] * tgt.second[1] + + tr2d[2] * tgt.second[2]; + +#ifdef CGAL_DEBUG_BSURF + std::cout << "4 " << tr2d[0] << " 0 " + << tr2d[1] << " 0 " + << tr2d[2] << " 0 " + << tr2d[0] << " 0\n"; +#endif + + return make_array(point_coords, + point_coords); + } + + //case k=0 assume that flat_tid has been flattened putting x-axis aligned with h_edge + static + std::array + unfold_face(const halfedge_descriptor& h_edge, + const VertexPointMap &vpm, const TriangleMesh &mesh, + const std::array& flat_tid,const int k=0) + { + halfedge_descriptor h_opp = opposite(h_edge, mesh); + + + vertex_descriptor v = target(next(h_opp,mesh),mesh); + vertex_descriptor a = source(h_edge, mesh); + vertex_descriptor b = target(h_edge,mesh); + FT r0 = squared_distance(get(vpm,v), get(vpm,a)); + FT r1 = squared_distance(get(vpm,v), get(vpm,b)); + + Vector_2 v2 = intersect_circles(flat_tid[(k+1)%3], r1, flat_tid[k], r0); + + halfedge_descriptor h_ref_opp = halfedge(face(h_opp,mesh),mesh); + + std::array res; + + if(h_ref_opp==h_opp) + { + res[0]=flat_tid[(k+1)%3]; + res[1]=flat_tid[k]; + res[2]=v2; + } else if(next(h_ref_opp,mesh)==h_opp) + { + res[0]=v2; + res[1]=flat_tid[(k+1)%3]; + res[2]=flat_tid[k]; + }else + { + assert(prev(h_ref_opp,mesh)==h_opp); + res[0]=flat_tid[k]; + res[1]=v2; + res[2]=flat_tid[(k+1)%3]; + } + + + + return res; + } + + static + std::array + unfold_face(halfedge_descriptor h_curr, halfedge_descriptor h_next, + const VertexPointMap &vpm, const TriangleMesh &mesh, + const std::array& flat_tid) + { + halfedge_descriptor h_next_opp = opposite(h_next, mesh); + CGAL_assertion(face( h_curr, mesh) == face(h_next_opp, mesh)); + + + vertex_descriptor v = target(next(h_curr,mesh),mesh); + vertex_descriptor a = target(h_curr,mesh); + vertex_descriptor b = source(h_curr, mesh); + FT r0 = squared_distance(get(vpm,v), get(vpm,a)); + FT r1 = squared_distance(get(vpm,v), get(vpm,b)); + + Vector_2 v2 = intersect_circles(flat_tid[1], r1, flat_tid[0], r0); + + + std::array res; + if(next(h_curr, mesh) == h_next_opp) + { + res[0]=flat_tid[0]; + //res[2]=flat_tid[1]; + res[1]=v2; +#ifdef CGAL_DEBUG_BSURF + std::cout << "4 " << res[0] << " 0 " + << res[1] << " 0 " + << flat_tid[1] << " 0 " + << res[0] << " 0\n"; +#endif + } + else + { + CGAL_assertion(prev(h_curr, mesh) == h_next_opp); + res[0]=v2; + res[1]=flat_tid[1]; + //res[2]=flat_tid[0]; +#ifdef CGAL_DEBUG_BSURF + std::cout << "4 " << res[0] << " 0 " + << res[1] << " 0 " + << flat_tid[0] << " 0 " + << res[0] << " 0\n"; +#endif + } + + return res; + } + static + std::vector< std::array> + unfold_strip(const std::vector& initial_path, + const CGAL::Polygon_mesh_processing::Face_location& src, + const CGAL::Polygon_mesh_processing::Face_location& tgt, + const VertexPointMap &vpm, const TriangleMesh &mesh) + { + std::size_t s=initial_path.size(); + std::vector> result(s+1); // Vector_2 should be Point_2 as they are funnel endpoints in 2D (after flattening) + result[0]=init_source_triangle(initial_path[0], vpm, mesh, src); +#ifdef CGAL_DEBUG_BSURF + std::cout << "unfolding faces\n"; +#endif + for(std::size_t i=1;i &path) { + // Among vertices around which the path curves, find the vertex + // with maximum angle. We are going to fix that vertex. Actually, max_index is + // the index of the first face containing that vertex. + std::size_t max_index = -1; + FT max_angle = 0.; + for (std::size_t i = 1; i < path.size()-1; ++i) { + Vector_2 pos = path[i].pos; + Vector_2 prev = path[i - 1].pos; + Vector_2 next = path[i + 1].pos; + Vector_2 v0 = pos - prev; + v0 = v0 / sqrt(v0.squared_length()); + Vector_2 v1 = next - pos; + v1 = v1 / sqrt(v1.squared_length()); + FT angle = 1 - scalar_product(v0, v1); + if (angle > max_angle) { + max_index = path[i].face; + max_angle = angle; + } + } + +#ifdef CGAL_DEBUG_BSURF + std::cout << "funnels ("<< max_index << ")"; + for (auto f : path) + std::cout << " " << f.pos << " |"; + std::cout << "\n"; +#endif + + return max_index; + } + + static + std::vector + funnel(const std::vector< std::array>& portals, std::size_t& max_index) + { + // Find straight path. + Vector_2 start(NULL_VECTOR); + int apex_index = 0; + int left_index = 0; + int right_index = 0; + Vector_2 apex = start; + Vector_2 left_bound = portals[0][0]; + Vector_2 right_bound = portals[0][1]; + + // Add start point. + std::vector points = std::vector{{apex_index, apex}}; + points.reserve(portals.size()); + // @Speed: is this slower than an inlined function? + auto area = [](const Vector_2 a, const Vector_2 b, const Vector_2 c) { // TODO replace with orientation predicate + return determinant(b - a, c - a); + }; + + for (std::size_t i = 0; i < portals.size(); ++i) + { + auto left = portals[i][0], right = portals[i][1]; + // Update right vertex. + if (area(apex, right_bound, right) <= 0) { + if (apex == right_bound || area(apex, left_bound, right) > 0) { + // Tighten the funnel. + right_bound = right; + right_index = i; + } else { + // Right over left, insert left to path and restart scan from + // portal left point. + if (left_bound != apex) { + points.push_back({left_index, left_bound}); + // Make current left the new apex. + apex = left_bound; + apex_index = left_index; + // Reset portal + left_bound = apex; + right_bound = apex; + left_index = apex_index; + right_index = apex_index; + // Restart scan + i = apex_index; + continue; + } + } + } + + // Update left vertex. + if (area(apex, left_bound, left) >= 0) { + if (apex == left_bound || area(apex, right_bound, left) < 0) { + // Tighten the funnel. + left_bound = left; + left_index = i; + } else { + if (right_bound != apex) { + points.push_back({right_index, right_bound}); + // Make current right the new apex. + apex = right_bound; + apex_index = right_index; + // Reset portal + left_bound = apex; + right_bound = apex; + left_index = apex_index; + right_index = apex_index; + // Restart scan + i = apex_index; + continue; + } + } + } + } + + // This happens when we got an apex on the last edge of the strip + if (points.back().pos != portals.back()[0]) { + points.push_back({(int)portals.size() - 1, portals.back()[0]}); + } + assert(points.back().pos == portals.back()[0]); + assert(points.back().pos == portals.back()[1]); + + std::vector lerps; + lerps.reserve(portals.size()); + for (std::size_t i = 0; i < points.size() - 1; i++) { + auto a = points[i].pos; + auto b = points[i + 1].pos; + for (auto k = points[i].face; k < points[i + 1].face; ++k) { + auto portal = portals[k]; + // assert(cross(b - a, portal.second - portal.first) > 0); + +#ifdef CGAL_DEBUG_BSURF + std::cout << "i=" << i << "\n"; + std::cout << "a=" << a << " b=" << b << " portal[0]=" << portal[0] << " portal[1]=" << portal[1] << "\n"; +#endif + FT s = intersect_segments(a, b, portal[0], portal[1]); + +#ifdef CGAL_DEBUG_BSURF + std::cout << "s=" << s << "\n"; +#endif + lerps.push_back(std::clamp(s, 0.0, 1.0)); + } + } + + auto index = 1; +#ifdef CGAL_DEBUG_BSURF + std::cout << "setting funnel_point indices\n"; +#endif + for (std::size_t i = 0; i < portals.size(); ++i) + { +#ifdef CGAL_DEBUG_BSURF + std::cout << " i=" << i << " index = " << index << "\n"; + std::cout << " portals[i][0]=" << portals[i][0] << " portals[i][1]=" << portals[i][1] << "\n"; + std::cout << " points[index].pos = " <>& portals, + std::vector& lerps, + std::vector& path, + CGAL::Polygon_mesh_processing::Face_location& src, + CGAL::Polygon_mesh_processing::Face_location& tgt, + const VertexPointMap &vpm, const TriangleMesh &mesh, std::size_t index) + { + namespace PMP = CGAL::Polygon_mesh_processing; + +#ifdef CGAL_DEBUG_BSURF + dump_path(path, lerps, src, tgt, mesh); +#endif + + typedef boost::graph_traits BGT; + vertex_descriptor vertex=BGT::null_vertex(); + + // TODO: use a while loop breaking when no apex vertices not already visited are available + for (std::size_t i = 0; i < portals.size() * 2 && index != std::size_t(-1); ++i) + { +#ifdef CGAL_DEBUG_BSURF + std::cout << "Improving path " << path.size() << " hedges\n"; + std::cout << "src = " << PMP::construct_point(src, mesh) << "\n"; + std::cout << "tgt = " << PMP::construct_point(tgt, mesh) << "\n"; +#endif + + vertex_descriptor new_vertex=BGT::null_vertex(); + halfedge_descriptor h_curr = path[index]; + + bool is_target = false; + if (lerps[index] == 0) { + new_vertex = target(h_curr,mesh); + is_target = true; + } else if (lerps[index] == 1) { + new_vertex = source(h_curr,mesh); + } + if (new_vertex == vertex) break; + vertex = new_vertex; + +#ifdef CGAL_DEBUG_BSURF + std::cout << " Current strip with Apex: " << get(vpm, new_vertex) << "\n"; + std::cout << " is_target? " << is_target << "\n"; + std::cout << " -- path --\n"; + for (auto h : path) + { + std::cout << " 4 " << get(vpm, source(h, mesh)) + << " " << get(vpm, target(h, mesh)) + << " " << get(vpm, target(next(h, mesh), mesh)) + << " " << get(vpm, source(h, mesh)) << std::endl; + std::cout << edge(h, mesh) << std::endl; + } + std::cout << " ----------\n"; +#endif + + // if I hit the source vertex v of h_curr, then h_next has v as source, thus we turn ccw around v in path + // Similarly, if I hit the target vertex v of h_curr, then h_next has v as target, thus we turn cw around v in path + std::size_t curr_index = index+1; + + std::vector new_hedges; + + if (is_target) + { + face_descriptor target_face; + + if (curr_index == path.size()) + { + target_face=tgt.first; + } + else + { + while (target(path[curr_index], mesh) == new_vertex) + { + if(curr_index==path.size()-1) + { + target_face=tgt.first; + curr_index=path.size(); + break; + } + ++curr_index; + } + if (curr_index != path.size()) + target_face = face(opposite(path[curr_index], mesh), mesh); + } + + halfedge_descriptor h_loop=opposite(prev(opposite(h_curr, mesh), mesh), mesh); + + while(target_face!=face(h_loop, mesh)) + { + new_hedges.push_back(h_loop); + h_loop=opposite(prev(h_loop,mesh), mesh); + } + new_hedges.push_back(h_loop); + } + else + { + face_descriptor target_face; + + if (curr_index == path.size()) + { + target_face=tgt.first; + } + else + { + while (source(path[curr_index], mesh) == new_vertex) + { + if(curr_index==path.size()-1) + { + target_face=tgt.first; + curr_index=path.size(); + break; + } + ++curr_index; + } + if (curr_index != path.size()) + target_face=face(opposite(path[curr_index], mesh), mesh); + } + halfedge_descriptor h_loop=opposite(next(opposite(h_curr, mesh), mesh), mesh); // skip the face before h_curr (as we won't remove it from path) + + while(target_face!=face(h_loop, mesh)) + { + new_hedges.push_back(h_loop); + h_loop=opposite(next(h_loop,mesh), mesh); + } + new_hedges.push_back(h_loop); + } + + // replace the halfedges incident to the apex vertex with the opposite part of the ring + std::vector new_path(path.begin(),path.begin()+index); + new_path.insert(new_path.end(), new_hedges.begin(), new_hedges.end()); + new_path.insert(new_path.end(), path.begin()+curr_index, path.end()); + path.swap(new_path); + + strip_path(mesh, src, tgt, path); + +#ifdef CGAL_DEBUG_BSURF + std::cout << " -- new path --\n"; + for (auto h : path) + { + std::cout << " 4 " << get(vpm, source(h, mesh)) + << " " << get(vpm, target(h, mesh)) + << " " << get(vpm, target(next(h, mesh), mesh)) + << " " << get(vpm, source(h, mesh)) << std::endl; + std::cout << face(h, mesh) << std::endl; + std::cout << edge(h, mesh) << std::endl; + } + std::cout << " ----------\n"; +#endif + + portals=unfold_strip(path,src,tgt,vpm,mesh); + lerps=funnel(portals,index); + CGAL_assertion(lerps.size()==path.size()); + +#ifdef CGAL_DEBUG_BSURF + dump_path(path, lerps, src, tgt, mesh); +#endif + } + +#ifdef CGAL_DEBUG_BSURF + std::cout << " Final strip\n"; + for (auto h : path) + { + std::cout << " 4 " << get(vpm, source(h, mesh)) + << " " << get(vpm, target(h, mesh)) + << " " << get(vpm, target(next(h, mesh), mesh)) + << " " << get(vpm, source(h, mesh)) << "\n"; + + } +#endif + + } + + + + #if 1 + //:::::::::::::::::::::Parallel Transport:::::::::::::::::::::::::::: + static + Vector_3 face_normal(const VertexPointMap &vpm, + const TriangleMesh &mesh, + const face_descriptor& f) + { halfedge_descriptor h=halfedge(f,mesh); + Point_3 p0=get(vpm,source(h,mesh)); + Point_3 p1=get(vpm,target(h,mesh)); + Point_3 p2=get(vpm,target(next(h,mesh),mesh)); + Vector_3 n = triangle_normal(p0,p1,p2); + n=n/sqrt(n.squared_length()); + + return n; + } + static + std::vector polar_basis(const VertexPointMap &vpm, + const TriangleMesh &mesh, + const face_descriptor& f, + const Vector_3& normal) + { + halfedge_descriptor h=halfedge(f,mesh); + Point_3 p0=get(vpm,source(h,mesh)); + Point_3 p1=get(vpm,target(h,mesh)); + Vector_3 e=p1-p0; + e=e/sqrt(e.squared_length()); + Vector_3 e_perp=cross_product(normal,e); + + return std::vector{e,e_perp,normal}; + + } + + + // //:::::::::::::::::::::Straightest Geodesic:::::::::::::::::::::::::::: + static + std::tuple> point_in_triangle(const VertexPointMap &vpm, + const TriangleMesh &mesh, + const face_descriptor& face, + const Point_3 point, + float tol=1e-5) { + // http://www.r-5.org/files/books/computers/algo-list/realtime-3d/Christer_Ericson-Real-Time_Collision_Detection-EN.pdf + // pag.48 + std::array b = make_array(0.,0.,0.); + halfedge_descriptor h=halfedge(face,mesh); + Point_3 v0 = get(vpm,source(h,mesh)); + Point_3 v1 = get(vpm,target(h,mesh)); + Point_3 v2 = get(vpm,target(next(h,mesh),mesh)); + + Vector_3 u = v1 - v0, v = v2 - v0, w = point - v0; + double d00 = u.squared_length(), d01 = u*v, d11 = v.squared_length(), d20 = w*u, + d21 = w*v, d = d00 * d11 - d01 * d01; + + if (d == 0) + return {false, make_array(0.,0.,0.)}; + + b[2] = (d00 * d21 - d01 * d20) / d; + assert(!isnan(b[2])); + b[1] = (d11 * d20 - d01 * d21) / d; + assert(!isnan(b[1])); + b[0] = 1 - b[1] - b[2]; + assert(!isnan(b[0])); + + for (auto i = 0; i < 3; ++i) { + if (b[i] < -tol || b[i] > 1.0 + tol) + return {false, make_array(0.,0.,0.)}; + } + + return {true, b}; +} + static + Eigen::Matrix3d rot_matrix (const FT& angle, const Vector_3& axis) + { + Eigen::Matrix3d result; + double c=std::cos(angle); + double s=std::sin(angle); + result < flat_from = init_flat_triangle(h,vpm,mesh); + std::array flat_to = unfold_face(h,vpm,mesh,flat_from); + Vector_2 bary = Vector_2{0.333, 0.333}; + Vector_2 c0 = 0.33*(flat_from[0]+flat_from[1]+flat_from[2]); + Vector_2 c1 = 0.33*(flat_to[0]+flat_to[1]+flat_to[2]); + Vector_2 e0 = flat_from[0] - c0; + Vector_2 e1 = flat_to[0] - c1; + + Vector_2 w = c1 - c0; + FT phi_ij = angle(e0, w); + if (e0.x()*w.y()-e0.y()*w.x() < 0) + phi_ij = 2 * CGAL_PI - phi_ij; + w *= -1; + FT phi_ji = angle(e1, w); + + if (e1.x()*w.y()-e1.y()*w.x() < 0) + phi_ji = 2 * CGAL_PI - phi_ji; + + Vector_3 n_from= face_normal(vpm, mesh, from); + std::vector e_from = polar_basis(vpm, mesh, from); + double teta = angle(e_from, v); + if (cross_product(e_from, v)*n_from < 0) + teta = 2 * CGAL_PI - teta; + + std::vector e_to = polar_basis(vpm, mesh, from ,to); + Vector_3 n_to = face_normal(vpm, mesh, to); + double rot = teta + phi_ji + CGAL_PI - phi_ij; + e_to =e_to*sqrt(v.squared_length()); + v = rotate_vector(e_to, n_to, rot); + + } + static + Eigen::Matrix3d transformation_matrix(const Vector_2 &x1, const Vector_2 &y1, + const Vector_2 &O1) + { + Eigen::Matrix3d T = Eigen::Matrix3d::Zero(); + T << x1.x(), y1.x(), O1.x(), x1.y(), y1.y(), O1.y(), 0, 0, 1; + + return T; + } + + //h_ref is the reference halfedge of the face we are in, h_edge is the halfedge along which we want to unfold + static + Vector_2 compute_new_dir(const halfedge_descriptor h_ref, + const halfedge_descriptor h_edge, + const std::array& prev_coords, + const std::array& curr_coords, + const VertexPointMap& vpm, + const TriangleMesh& mesh) + { + int k=0; + + if(h_edge==prev(h_ref,mesh)) + k=2; + else if(h_edge==next(h_ref,mesh)) + k=1; + else + assert(h_edge==h_ref); + + std::array flat_curr = init_flat_triangle(h_ref,vpm,mesh); + std::array flat_prev = unfold_face(h_edge,vpm,mesh,flat_curr,k); + + Vector_2 prev_flat_p=prev_coords[0]*flat_prev[0]+prev_coords[1]*flat_prev[1]+prev_coords[2]*flat_prev[2]; + Vector_2 curr_flat_p=curr_coords[0]*flat_curr[0]+curr_coords[1]*flat_curr[1]+curr_coords[2]*flat_curr[2]; + +#ifdef CGAL_DEBUG_BSURF + std::cout<<" k is "< flat_curr = init_flat_triangle(h_ref,vpm,mesh); + std::array flat_prev = unfold_face(h_edge,vpm,mesh,flat_curr,k); + + Vector_2 prev_origin = flat_prev[0]; + Vector_2 prev_y = flat_prev[1]-flat_prev[0]; + Vector_2 prev_x = Vector_2(prev_y.y(), - prev_y.x()); + + Eigen::Matrix3d T= transformation_matrix(prev_x, prev_y, prev_origin); + return switch_reference_frame_vector(T, prev_dir); + } + + static Vector_2 switch_reference_frame_vector(const Eigen::Matrix3d &T, const Vector_2 &p) { + Eigen::Vector3d V; + V << p.x(), p.y(), 0; + Eigen::Vector3d t = T * V; + return Vector_2(FT(t(0)), FT(t(1))); + } + + static + void parallel_transport_through_flattening(Vector_2& v, + const VertexPointMap &vpm, + const TriangleMesh &mesh, + const face_descriptor& from, + const face_descriptor& to) + { + + + } + + Vector_3 parallel_transport_along_path(const std::vector>& edge_locations, + const VertexPointMap &vpm, + const TriangleMesh& mesh, + const CGAL::Polygon_mesh_processing::Face_location& src, + const CGAL::Polygon_mesh_processing::Face_location& tgt, + const Vector_3& v) + { + if(edge_locations.size()==0) + return v; + Vector_3 w=v; + face_descriptor f_start=src.first; + for(std::size_t i=0;i + // polthier_condition_at_vert(const TriangleMesh& mesh, + // const VertexPointMap &vpm + // const vertex_descriptor& vid, + // const face_descriptor& tid, + // const int kv, + // const Vector_3 &dir) + // { + // FT total_angle=get_total_angle(vid,mesh,vpm); + // FT theta = 0.5 * total_angle; + // halfedge_descriptor h=halfedge(tid,mesh); + // while(target(h,mesh)!=vid) + // h=next(h,mesh); + + // Point_3 vert = get(vpm,vid); + // Point_3 vert_adj=get(vpm,source(h,mesh)); + // Vector_3 v = vert - vert_adj; + + // FT acc = angle(v, dir); + // FT prev_angle = acc; + // face_descriptor curr_tid = tid; + + // while (acc < theta) { + // h=prev(opposite(h.mesh),mesh); + // prev_angle = acc; + // Point_3 next_vert_adj=get(vpm,source(h,mesh)); + // acc += angle(vert_adj - vert,next_vert_adj - vert); + // vert_adj=next_vert_adj; + // } + // auto offset = theta - prev; + // Point_3 prev_vert_adj=get(vpm,target(next(h,mesh,mesh))); + + // Vector_3 result=Vector_3{prev_vert_adj.x()-vert.x(),prev_vert_adj.y()-vert.y(),prev_vert_adj.z()-vert.z()}; + // face_descriptor f=face(h,mesh); + // Vector_3 n=face_normal(vpm,mesh,f); + // result=rotate_vector(result,n,offset); + // return {result,f}; + // } + + + static + std::vector + get_3D_basis_at_point(const CGAL::Polygon_mesh_processing::Face_location& p, + const TriangleMesh& mesh, + const VertexPointMap &vpm) + { + halfedge_descriptor h=halfedge(p.first,mesh); + + Point_3 p0=get(vpm, source(h, mesh)); + Point_3 p1=get(vpm, target(h, mesh)); + Point_3 p2= get(vpm,target(next(h, mesh), mesh)); + + Vector_3 e=p1-p0; + e/=sqrt(e.squared_length()); + Vector_3 n= triangle_normal(p0,p1,p2); + Vector_3 e1= cross_product(n,e); + + return {e,e1,n}; + } + + + static + std::tuple + point_is_edge(const CGAL::Polygon_mesh_processing::Face_location& p, + const FT& tol=1e-5) + { + auto bary=p.second; + if (bary[0] > tol && bary[1] > tol && bary[2] <= tol) + return {true, 0}; + if (bary[1] > tol && bary[2] > tol && bary[0] <= tol) + return {true, 1}; + if (bary[2] > tol && bary[0] > tol && bary[1] <= tol) + return {true, 2}; + + return {false, -1}; + } + +//TODO get rid of tol or at least get it dependent on the input? +// + static + std::tuple + point_is_vert(const CGAL::Polygon_mesh_processing::Face_location& p,const FT& tol=1e-5) + { + auto bary=p.second; + if (bary[0] > tol && bary[1] <= tol && bary[2] <= tol) + return {true, 0}; + if (bary[1] > tol && bary[0] <= tol && bary[2] <= tol) + return {true, 1}; + if (bary[2] > tol && bary[0] <= tol && bary[1] <= tol) + return {true, 2}; + + return {false, -1}; + } + + static + //https://rootllama.wordpress.com/2014/06/20/ray-line-segment-intersection-test-in-2d/ + std::pair + intersect(const Vector_2 &direction, const Vector_2 &left, + const Vector_2 &right,const Vector_2& origin) + { + Vector_2 v1 = origin-left; + Vector_2 v2 = right - left; + Vector_2 v3(-direction.y(), direction.x()); + double t0 = (v2.x()*v1.y()-v2.y()*v1.x()) / (v2*v3); + double t1 = v1*v3/ ( v2*v3 ); + return std::make_pair(t0, t1); + }; + + static + std::tuple + segment_in_tri(const Vector_2& p, + const std::array& tri, + const Vector_2& dir,const int offset) + { + //rotated the triangle in order to test intersection at meaningful edges before + std::array rotated_tri=tri; + // TODO rotated_tri.begin() + offset ? + for(int k=0;k 0 && t1 >= -1e-4 && t1 <= 1 + 1e-4) + { +#ifdef CGAL_DEBUG_BSURF + int h=((i+1)%3 + offset)%3; + Vector_2 intersection_point=(1-t1)*tri[h]+t1*tri[(h+1)%3]; + std::cout<<" 2D intersection point: "<< intersection_point << std::endl; +#endif + return {((i+1)%3 + offset)%3, std::clamp(t1, 0., 1.)}; //return the offset w.r.t h_ref + } + } + CGAL_assertion(false); + return {-1,-1}; + } + + // return an angle in degrees + static + FT get_total_angle(const vertex_descriptor& vid, + const TriangleMesh& mesh, + const VertexPointMap &vpm) + { + halfedge_descriptor h_ref= halfedge(vid,mesh); + halfedge_descriptor h_start=h_ref; + halfedge_descriptor h_next = next(h_start,mesh); + + FT theta=approximate_angle(get(vpm,source(h_start,mesh))-get(vpm,target(h_start,mesh)), + get(vpm,target(h_next,mesh))-get(vpm,target(h_start,mesh))); + h_start=opposite(h_next,mesh); + h_next=next(h_start,mesh); + while(h_start!=h_ref) + { + theta+=approximate_angle(get(vpm,source(h_start,mesh))-get(vpm,target(h_start,mesh)), + get(vpm,target(h_next,mesh))-get(vpm,target(h_start,mesh))); + h_start=opposite(h_next,mesh); + h_next=next(h_start,mesh); + } + + return theta; + } + + static + std::tuple + polthier_condition_at_vert(const TriangleMesh& mesh, + const VertexPointMap &vpm, + const vertex_descriptor& vid, + const face_descriptor& tid, + const FT& init_angle) + { + auto get_vid_offset=[&mesh](const halfedge_descriptor& h_ref,const vertex_descriptor& vid) + { + if(source(h_ref,mesh)==vid) + return 0; + + if(target(h_ref,mesh)==vid) + return 1; + + if(target(next(h_ref,mesh),mesh)==vid) + return 2; + + std::cerr<<"Error! Halfedges are in different faces"< flat_tid = init_flat_triangle(halfedge(face(h,mesh),mesh),vpm,mesh); + int kv=get_vid_offset(halfedge(face(h,mesh),mesh),vid); + + + Vector_2 q = (1 - alpha) * flat_tid[(kv+1)%3] + alpha * flat_tid[(kv+2)%3]; + Vector_2 new_dir = q - flat_tid[kv]; +#ifdef CGAL_DEBUG_BSURF + Vector_2 flat_p=flat_tid[kv]; + std::cout<<"final h "<< edge(h,mesh)<> + straightest_geodesic(const CGAL::Polygon_mesh_processing::Face_location& start_loc, + const TriangleMesh& mesh, + const VertexPointMap &vpm, + Vector_2 dir,const FT& len) + { + namespace PMP = CGAL::Polygon_mesh_processing; + +#ifdef CGAL_DEBUG_BSURF + { + std::ofstream log("/tmp/log.txt"); + log << std::setprecision(17); + log << start_loc.first << " " << start_loc.second[0] << " " << start_loc.second[1] << " " << start_loc.second[2] << "\n"; + log << "dir = " << dir << "\n"; + log << "len = " << len << "\n"; + } +#endif + + if (len == 0) return {start_loc}; + + auto get_halfedge_offset=[&mesh](const halfedge_descriptor& h_ref,const halfedge_descriptor& h_curr) + { + if( h_ref == h_curr) return 0; + if( next(h_ref, mesh) == h_curr ) return 1; + if( prev(h_ref, mesh) == h_curr ) return 2; + + std::cout<<"Error! Halfedges are in different faces"<& bary_edge) + { + std::array bary_edge_in_face=make_array(0.,0.,0.); + if (h_face!=h_edge) + { + if (h_face==next(h_edge, mesh)) + { + bary_edge_in_face[0]=bary_edge[1]; + bary_edge_in_face[1]=0; + bary_edge_in_face[2]=bary_edge[0]; + } + else + { + assert(h_face==prev(h_edge,mesh)); + bary_edge_in_face[0]=0; + bary_edge_in_face[1]=bary_edge[0]; + bary_edge_in_face[2]=bary_edge[1]; + } + } + else + { + bary_edge_in_face[0]=bary_edge[0]; + bary_edge_in_face[1]=bary_edge[1]; + bary_edge_in_face[2]=0; + } + + return bary_edge_in_face; + }; + + std::vector> result; + FT accumulated=0.; + face_descriptor curr_tid=start_loc.first; + std::array curr_flat_tid=init_flat_triangle(halfedge(curr_tid,mesh),vpm,mesh); + Vector_2 flat_p= start_loc.second[0]*curr_flat_tid[0]+start_loc.second[1]*curr_flat_tid[1]+start_loc.second[2]*curr_flat_tid[2]; + PMP::Face_location curr_p=start_loc; + + PMP::descriptor_variant var_des = PMP::get_descriptor_from_location(start_loc,mesh); + switch(var_des.index()) + { + case 0: + { + // TODO: handle is src_vertex is a border vertex: does not necessarily mean early exist as it depends on the direction + vertex_descriptor src_vertex = std::get(var_des); + halfedge_descriptor href = halfedge(curr_tid, mesh); + int src_index = source(href, mesh)==src_vertex?0:(target(href,mesh)==src_vertex?1:2); + // first get the direction in the basis of the triangle to get its angle with an edge incident to the vertex + + Vector_2 v02 = curr_flat_tid[(src_index+2)%3] - curr_flat_tid[src_index]; + Vector_2 v01 = curr_flat_tid[(src_index+1)%3] - curr_flat_tid[src_index]; + double unnorm_cos_theta_1 = v02 * dir; + double unnorm_cos_theta_2 = v02 * v01; + + // 2 orientation tests to check if the angle is larger or smaller than Pi + // 2 orientation tests are needed as curr_flat_tid orientation could be inverted during flattening + Orientation dir_ori=orientation(v02,dir), v01_ori=orientation(v02,v01); + + // check if dir is the cone centered at curr_flat_tid[src_vertex] (TODO: should be a predicate!) + if (orientation(dir, v01)==LEFT_TURN || orientation(dir, v02)==RIGHT_TURN) + { + //TODO: should be made robust with snapping and predicates or something à la robust construction traits + + // compute all the angles around the vertex from where the path starts + std::vector angles; + double acc_angle=0; + halfedge_descriptor hloop=href; + while (target(hloop, mesh)!=src_vertex) hloop=next(hloop, mesh); + halfedge_descriptor hstart=hloop; + Vector_3 n(get(vpm,target(hloop, mesh)), get(vpm,source(hloop, mesh))); + n /= std::sqrt(n*n); + do{ + hloop=opposite(next(hloop, mesh), mesh); + Vector_3 nv(get(vpm,target(hloop, mesh)), get(vpm,source(hloop, mesh))); + nv /= std::sqrt(nv*nv); + angles.push_back( std::acos(n * nv) ); + acc_angle+=angles.back(); + n = nv; + }while(hloop!=hstart); + + CGAL_assertion( std::acos(unnorm_cos_theta_2 / std::sqrt((v01*v01) * (v02*v02))) == angles[0]); // will not be true because of rounding + + // normalise the angle to bring them in [0,1] + for (double& angle : angles) + angle /= acc_angle; + + // compute and normalise the target angle + double target_theta = std::acos(unnorm_cos_theta_1 / std::sqrt((dir*dir) * (v02*v02))); + if ( dir_ori!=v01_ori ) target_theta = 2 * CGAL_PI - target_theta; + target_theta /= acc_angle; + + double curr_angle = 0; + int ia = 0; + do{ + curr_angle+=angles[ia]; + if (curr_angle >= target_theta) break; + hloop=opposite(next(hloop, mesh), mesh); + ++ia; + }while(hloop!=hstart); + + CGAL_assertion(hloop!=hstart); + + // angle in target face wrt hloop + double delta = (target_theta - curr_angle + angles[ia]) * acc_angle; + + // using the law of the sinus we have: + CGAL_assertion(target(hloop, mesh) == src_vertex); + Vector_3 vloop(get(vpm, source(hloop, mesh)),get(vpm, target(hloop, mesh))); + Vector_3 vprev(get(vpm, source(hloop, mesh)),get(vpm, target(next(hloop, mesh), mesh))); + double dvloop = std::sqrt(vloop*vloop); + double dvprev = std::sqrt(vprev*vprev); + double theta = std::acos( vloop*vprev / dvloop / dvprev ); + double d_opp_delta = dvloop * sin(delta) / sin(CGAL_PI-delta-theta); + + double alpha = d_opp_delta/dvprev; + + if (is_border(hloop,mesh)) + { + // border case falling outside of the mesh TODO: check with Claudio + result.push_back(curr_p); + return result; + } + + halfedge_descriptor href = halfedge(face(hloop,mesh),mesh); + curr_flat_tid=init_flat_triangle(href,vpm,mesh); + int src_id = source(href, mesh)==src_vertex?0:(target(href,mesh)==src_vertex?1:2); + +#ifdef CGAL_DEBUG_BSURF + std::cout << "new_face= " << face(href, mesh) << "\n"; + std::array debug_pt; + debug_pt[0]=mesh.point(source(href,mesh)) - Point_3(0.,0.,0.); + debug_pt[1]=mesh.point(target(href,mesh)) - Point_3(0.,0.,0.); + debug_pt[2]=mesh.point(target(next(href,mesh),mesh)) - Point_3(0.,0.,0.); + std::cout << "direction inter pt: " << (alpha) * debug_pt[(src_id+1)%3] + (1-alpha) * debug_pt[(src_id+2)%3] << "\n"; +#endif + // point on the opposite edge of src_vertex in face(href, mesh) + Vector_2 ip = (alpha) * curr_flat_tid[(src_id+1)%3] + (1-alpha) * curr_flat_tid[(src_id+2)%3]; + dir = ip - curr_flat_tid[src_id]; + curr_tid=face(href, mesh); + flat_p = curr_flat_tid[src_id]; + curr_p.second=CGAL::make_array(0.,0.,0.); + curr_p.second[src_id]=1.; + curr_p.first=curr_tid; + } + } + break; + case 1: + { + halfedge_descriptor h = std::get(var_des); + halfedge_descriptor href = halfedge(face(h,mesh), mesh); + int opp_id = h==href?2:(next(href,mesh)==h?0:1); + Vector_2 h_2d = curr_flat_tid[(opp_id+2)%3] - curr_flat_tid[(opp_id+1)%3]; + + // check if the line starts in the current face (TODO should be a predicate) + if ( orientation(h_2d, dir) == RIGHT_TURN ) + { + if ( is_border(opposite(h, mesh), mesh) ) + { + result.push_back(curr_p); + return result; + } + + // take the same point but in the opposite face + halfedge_descriptor h_start=h; + h=opposite(h, mesh); + curr_tid=face(h, mesh); + curr_p.first=curr_tid; + href=halfedge(curr_tid, mesh); + CGAL_assertion(start_loc.second[opp_id]==0); + std::array bary2 = CGAL::make_array(start_loc.second[(opp_id+2)%3], start_loc.second[(opp_id+1)%3]); // swap done here + int opp_id = h==href?2:(next(href,mesh)==h?0:1); + curr_p.second[opp_id]=FT(0); + curr_p.second[(opp_id+1)%3]=bary2[0]; + curr_p.second[(opp_id+2)%3]=bary2[1]; + + + // dir must also be updated: + // unfold the new triangle into the basis of the input one + // compute the angle between the new triangle ref halfedge and dir + // compute the new dir thanks to that angle. + halfedge_descriptor start_href=halfedge(start_loc.first, mesh); + int k=h_start==prev(start_href,mesh)?2 :( h_start==next(start_href,mesh)?1:0); + std::array new_flat_tid_in_curr_basis=unfold_face(h_start,vpm, mesh,curr_flat_tid, k); + + Vector_2 ybase = new_flat_tid_in_curr_basis[1]-new_flat_tid_in_curr_basis[0]; + double cos_theta = ybase * dir / std::sqrt(ybase*ybase) / std::sqrt(dir*dir); + double theta = std::acos(cos_theta); + dir = Vector_2(-std::sin(theta), std::cos(theta)); + + curr_flat_tid=init_flat_triangle(halfedge(curr_tid,mesh),vpm,mesh); + flat_p= curr_p.second[0]*curr_flat_tid[0]+curr_p.second[1]*curr_flat_tid[1]+curr_p.second[2]*curr_flat_tid[2]; + } + } + break; + default: + break; + } + + PMP::Face_location prev_p; + Vector_2 curr_dir=dir; + halfedge_descriptor h_ref=halfedge(curr_tid,mesh); + halfedge_descriptor h_curr=h_ref; + + + + result.push_back(curr_p); +#ifdef CGAL_DEBUG_BSURF + std::cout << "p= " << PMP::construct_point(curr_p,mesh) << ")\n"; +#endif + + //TODO not sure why we need that + auto [is_vert, kv] = point_is_vert(curr_p); + auto [is_edge, ke] = point_is_edge(curr_p); + if (is_vert) + h_curr=get_halfedge(kv,h_ref); + else if (is_edge) + h_curr=get_halfedge(ke,h_ref); + +#ifdef CGAL_DEBUG_BSURF + std::cout << "Accumulated loop starts\n"; +#endif + while (accumulated < len) + { +#ifdef CGAL_DEBUG_BSURF + std::cout << "--->" << accumulated << " vs " << len << "\n"; +#endif + int curr_offset=get_halfedge_offset(h_ref,h_curr); + + auto [k, t1] = segment_in_tri(flat_p, curr_flat_tid, curr_dir,curr_offset); +#ifdef CGAL_DEBUG_BSURF + std::cout << " t1 = " << t1 << std::endl; +#endif + CGAL_assertion(k!=-1); + std::array new_bary=make_array(0.,0.,0.); + PMP::Face_location point_on_edge; + + new_bary[k] = 1 - t1; + new_bary[(k + 1) % 3] = t1; + point_on_edge.first=curr_tid; + point_on_edge.second=new_bary; + std::tie(is_vert, kv) = point_is_vert(point_on_edge); +#ifdef CGAL_DEBUG_BSURF + std::cout << " is_vert? " << is_vert << "\n"; +#endif + if (is_vert) + { + point_on_edge.second = make_array(0.,0.,0.); + point_on_edge.second[kv]=1; + vertex_descriptor vid = target(get_halfedge(kv, prev(h_ref,mesh)), mesh); +#ifdef CGAL_DEBUG_BSURF + std::cout<< "hit vertex "<< vid < len TODO: what about below for edges? + break; + } + if (border_reached) break; + + + // Point_3 vert = get(vpm,vid); + // Point_3 vert_adj=get(vpm,source(h,mesh)); + // Vector_3 v = vert - vert_adj; + + //TODO add a 2D version of approximate_angle in CGAL + // FT init_angle = approximate_angle(curr_flat_tid[(kv+2)%3]-curr_flat_tid[kv], curr_dir); + auto tmp =curr_flat_tid[kv]-curr_flat_tid[(kv+2)%3]; + FT init_angle = approximate_angle(Vector_3(tmp.x(), tmp.y(), 0), Vector_3(curr_dir.x(), curr_dir.y(), 0)); + std::tie(curr_dir, curr_tid, h_curr) = + polthier_condition_at_vert(mesh,vpm,vid,curr_tid, init_angle); + + h_ref=halfedge(curr_tid,mesh); + kv=get_vid_offset(h_ref,vid); + CGAL_assertion(kv!=-1); + std::array tmp_bary=make_array(0.,0.,0.); + tmp_bary[kv]=1; + + curr_flat_tid=init_flat_triangle(h_ref,vpm,mesh); + prev_p = curr_p; + + curr_p=PMP::Face_location(curr_tid,tmp_bary); + int k=get_vid_offset(h_ref,target(h_curr,mesh)); + flat_p=curr_flat_tid[k]; + } + else + { + h_curr=opposite(get_halfedge(k, h_ref),mesh); + + // break if hitting the boundary + if (is_border(h_curr, mesh)) + { + result.push_back(point_on_edge); + accumulated += sqrt(squared_distance(PMP::construct_point(point_on_edge,mesh), + PMP::construct_point(curr_p,mesh))); + prev_p=point_on_edge; //if accumulated > len we need to be in the common face + break; + } + + face_descriptor adj = face(h_curr,mesh); + std::array curr_alpha=make_array(t1,1-t1); //reversed because will switch face + new_bary=edge_barycentric_coordinate(h_curr,halfedge(adj,mesh),curr_alpha); + prev_p = curr_p; + curr_p.first=adj; + curr_p.second= new_bary; + accumulated += sqrt(squared_distance(PMP::construct_point(curr_p,mesh), + PMP::construct_point(prev_p,mesh))); + curr_tid = adj; + h_ref=halfedge(curr_tid,mesh); + //TODO curr_dir should be normalized everytime (to try to avoid numerical errors) + curr_dir = compute_new_dir(h_ref,h_curr,prev_p.second,curr_p.second,vpm,mesh); + curr_flat_tid=init_flat_triangle(h_ref,vpm,mesh); + flat_p= curr_p.second[0]*curr_flat_tid[0]+curr_p.second[1]*curr_flat_tid[1]+curr_p.second[2]*curr_flat_tid[2]; +#ifdef CGAL_DEBUG_BSURF + std::cout << " h_curr " << edge(h_curr, mesh)<<"\n"; + std::cout << " adj " << adj<<"\n"; + Vector_2 intersection_point=new_bary[0]*curr_flat_tid[0]+new_bary[1]*curr_flat_tid[1]+new_bary[2]*curr_flat_tid[2]; + std::cout << " New intersection point is "<< intersection_point< + // polthier_straightest_geodesic(const vector &triangles, const vector &positions, + // const vector &adjacencies, const vector> &v2t, + // const vector> &angles, const vector &total_angles, + // const mesh_point &p, const vec2f &dir, const float &len) + // { + // auto result = vector{}; + // auto accumulated = 0.f; + // auto curr_tid = p.face; + // auto curr_p = p; + // auto prev_p = mesh_point{}; + // auto curr_dir = dir; + + // result.push_back(p); + + // auto k_start = 0; + // auto [is_vert, kv] = point_is_vert(p); + // auto [is_edge, ke] = point_is_edge(p); + // if (is_vert) + // k_start = kv; + // else if (is_edge) + // k_start = ke; + + // auto count = 0; + // while (accumulated < len) + // { + // ++count; + // auto [k, t1] = straightest_path_in_tri(triangles, positions, curr_p, + // curr_dir, k_start); + + // auto new_bary = zero3f; + // auto point_on_edge = mesh_point{}; + // if (k != -1) { + // new_bary[k] = 1 - t1; + // new_bary[(k + 1) % 3] = t1; + // point_on_edge = mesh_point{curr_tid, vec2f{new_bary.y, new_bary.z}}; + // } else { + // std::tie(is_edge, ke) = point_is_edge(curr_p, 5e-3); + // std::tie(is_vert, kv) = point_is_vert(curr_p, 5e-3); + // auto bary = get_bary(curr_p.uv); + // if (is_edge) { + // k = ke; + // t1 = bary[(k + 1) % 3]; + // point_on_edge = curr_p; + // } else if (is_vert) { + // auto bary3 = zero3f; + // bary3[kv] = 1; + // point_on_edge = {curr_p.face, {bary3.y, bary3.z}}; + // } else { + // std::cout << "Error!This should not happen" << std::endl; + // return result; + // } + // } + + // std::tie(is_vert, kv) = point_is_vert(point_on_edge); + + // if (is_vert) { + // auto vid = triangles[curr_tid][kv]; + // if (angles[vid].size() == 0) + // return result; + + // accumulated += + // length(eval_position(triangles, positions, curr_p) - positions[vid]); + // auto dir3d = normalize(eval_position(triangles, positions, curr_p) - + // positions[vid]); + // prev_p = curr_p; + // std::tie(curr_dir, curr_p, k_start) = + // polthier_condition_at_vert(triangles, positions, adjacencies, + // total_angles, vid, curr_tid, dir3d); + // curr_tid = curr_p.face; + // if (curr_tid == -1) + // return result; + + // } else { + // auto adj = adjacencies[curr_tid][k]; + // if (adj == -1) + // return result; + // auto h = find(adjacencies[adj], curr_tid); + + // new_bary = zero3f; + // new_bary[h] = t1; + // new_bary[(h + 1) % 3] = 1 - t1; + + // prev_p = curr_p; + // curr_p = mesh_point{adj, vec2f{new_bary.y, new_bary.z}}; + // accumulated += length(eval_position(triangles, positions, curr_p) - + // eval_position(triangles, positions, prev_p)); + + // auto T = switch_reference_frame(triangles, positions, adj, curr_tid); + // curr_dir = switch_reference_frame_vector(T, curr_dir); + + // curr_tid = adj; + // k_start = h; + // } + + // result.push_back(curr_p); + // } + + // auto excess = accumulated - len; + // auto prev_pos = eval_position(triangles, positions, result.rbegin()[1]); + // auto last_pos = eval_position(triangles, positions, result.back()); + // auto alpha = excess / length(last_pos - prev_pos); + // auto pos = alpha * prev_pos + (1 - alpha) * last_pos; + + // auto [inside, bary] = + // point_in_triangle(triangles, positions, prev_p.face, pos); + // if (!inside) + // std::cout << "error!This point should be in the triangle" << std::endl; + + // result.pop_back(); + // result.push_back(mesh_point{prev_p.face, bary}); + + // return result; + // } +#endif + + static + void + strip_path(const TriangleMesh& tmesh, + CGAL::Polygon_mesh_processing::Face_location& src, + CGAL::Polygon_mesh_processing::Face_location& tgt, + std::vector& initial_path) + { + namespace PMP = CGAL::Polygon_mesh_processing; + + // retrieve the vertex id of a face location describing a vertex + auto is_vertex = [](const PMP::Face_location& fl) + { + if (fl.second[0]==0 && fl.second[1]==0) + return 2; + if (fl.second[1]==0 && fl.second[2]==0) + return 0; + if (fl.second[0]==0 && fl.second[2]==0) + return 1; + return -1; + }; + + // retrieve the source vertex id of a face location describing a point on a halfedge + auto is_edge = [](const PMP::Face_location& fl) + { + if (fl.second[0]==0 && fl.second[1]!=0 && fl.second[2]!=0) + return 2; + if (fl.second[1]==0 && fl.second[2]!=0 && fl.second[0]!=0) + return 0; + if (fl.second[2]==0 && fl.second[0]!=0 && fl.second[1]!=0) + return 1; + return -1; + }; + + // strip initial_path from extra halfedges after src and also update src + int vi = is_vertex(src); + if (vi!=-1) + { + halfedge_descriptor hsrc=prev(halfedge(src.first, tmesh), tmesh); + for (int i=0; i fl_src = {FT(0),FT(0),FT(0)}; + fl_src[vid]=FT(1); + src=std::make_pair(face(hsrc, tmesh), fl_src); + } + } + else + { + int ei = is_edge(src); + if (ei!=-1) + { + halfedge_descriptor hsrc=prev(halfedge(src.first, tmesh), tmesh); + for (int i=0; i fl_src = {FT(0),FT(0),FT(0)}; + fl_src[hid]=src.second[(ei+2)%3]; + fl_src[(hid+2)%3]=src.second[ei]; + src=std::make_pair(face(new_hsrc, tmesh), fl_src); + } + } + } + + // strip initial_path from extra halfedges before tgt and also update tgt + vi = is_vertex(tgt); + if (vi!=-1) + { + halfedge_descriptor htgt=prev(halfedge(tgt.first, tmesh), tmesh); + for (int i=0; i fl_tgt = {FT(0),FT(0),FT(0)}; + fl_tgt[vid]=FT(1); + tgt=std::make_pair(face(new_htgt, tmesh), fl_tgt); + } + } + else + { + int ei = is_edge(tgt); + if (ei!=-1) + { + halfedge_descriptor htgt=prev(halfedge(tgt.first, tmesh), tmesh); + for (int i=0; i fl_tgt = {FT(0),FT(0),FT(0)}; + fl_tgt[hid]=tgt.second[(ei+2)%3]; + fl_tgt[(hid+2)%3]=tgt.second[ei]; + tgt=std::make_pair(face(new_htgt, tmesh), fl_tgt); + } + } + } + } +}; // end of Locally_shortest_path_imp + +template +struct Bezier_tracing_impl +{ + using face_descriptor = + typename boost::graph_traits::face_descriptor; + using vertex_descriptor = + typename boost::graph_traits::vertex_descriptor; + using halfedge_descriptor = + typename boost::graph_traits::halfedge_descriptor; + + using Point_2 = typename K::Point_2; + using Point_3 = typename K::Point_3; + using Vector_2 = typename K::Vector_2; + using Vector_3 = typename K::Vector_3; + using FT = typename K::FT; + + using Face_location = Polygon_mesh_processing::Face_location; + using Edge_location = Polygon_mesh_processing::Edge_location; + + template + static + std::vector + get_positions(const EdgeLocationRange& edge_locations, + const TriangleMesh& mesh, + const Face_location& src, + const Face_location& tgt) + { + namespace PMP = CGAL::Polygon_mesh_processing; + + std::vector result; + result.reserve(edge_locations.size()+2); + result.push_back(PMP::construct_point(src,mesh)); + for(auto& e: edge_locations) + result.push_back(PMP::construct_point(e,mesh)); + + result.push_back(PMP::construct_point(tgt,mesh)); +//TODO: we must guarantee that result is sorted and unique (rounding issue?) + return result; + } + + template + static + std::vector + path_parameters(const EdgeLocationRange& edge_locations, + const TriangleMesh& mesh, + const Face_location& src, + const Face_location& tgt) + { + std::vector pos=get_positions(edge_locations,mesh,src,tgt); + FT L=0.; + std::vector result(pos.size()); + for(std::size_t i=0;i + static + Face_location + eval_point_on_geodesic(const EdgeLocationRange& edge_locations, + const TriangleMesh& mesh, + const Face_location& src, + const Face_location& tgt, + const std::vector& parameters,/// edge length parameterization of the path from src to tgt through edge_locations + const FT& t) + { + if (t==0) return src; + if (t==1) return tgt; + + if(src.first==tgt.first) + { + std::array bary; + bary[0]=(1-t)*src.second[0]+t*tgt.second[0]; + bary[1]=(1-t)*src.second[1]+t*tgt.second[1]; + bary[2]=(1-t)*src.second[2]+t*tgt.second[2]; + return {src.first,bary}; + } + + std::size_t i = 0; + for (; i < parameters.size() - 1; i++) + { + if (parameters[i + 1] >= t) break; + } + FT t_low = parameters[i]; + FT t_high = parameters[i + 1]; + CGAL_assertion(t_high!=t_low); + FT alpha = (t - t_low) / (t_high - t_low); + std::array bary_low; + std::array bary_high; + + // warning there is an offset of the index: parameters contains one extra element (src) at 0 + // while edge_locations does not + face_descriptor curr_tid = i==0?src.first:face(halfedge(edge_locations[i-1].first,mesh),mesh); + halfedge_descriptor h_face = halfedge(curr_tid, mesh); + auto edge_barycentric_coordinate = + [&mesh, h_face](halfedge_descriptor h_edge, + const std::array& bary_edge) + { + std::array bary_edge_in_face; + if (h_face!=h_edge) + { + if (h_face==next(h_edge, mesh)) + { + bary_edge_in_face[0]=bary_edge[1]; + bary_edge_in_face[1]=0; + bary_edge_in_face[2]=bary_edge[0]; + } + else + { + bary_edge_in_face[0]=0; + bary_edge_in_face[1]=bary_edge[0]; + bary_edge_in_face[2]=bary_edge[1]; + } + } + else + { + bary_edge_in_face[0]=bary_edge[0]; + bary_edge_in_face[1]=bary_edge[1]; + bary_edge_in_face[2]=0; + } + + return bary_edge_in_face; + }; + + if(i==0) + bary_low=src.second; + else + { + halfedge_descriptor h_low = halfedge(edge_locations[i-1].first, mesh); + bary_low = edge_barycentric_coordinate(h_low, edge_locations[i-1].second); + } + + if(i==parameters.size()-2) + bary_high=tgt.second; + else + { + halfedge_descriptor h_high = opposite(halfedge(edge_locations[i].first, mesh), mesh); + CGAL_assertion(face(h_high,mesh)==curr_tid); + std::array edge_bary_high=edge_locations[i].second; + std::swap(edge_bary_high[0],edge_bary_high[1]); + bary_high = edge_barycentric_coordinate(h_high, edge_bary_high); + } + + std::array bary; + bary[0]=(1-alpha)*bary_low[0]+alpha*bary_high[0]; + bary[1]=(1-alpha)*bary_low[1]+alpha*bary_high[1]; + bary[2]=(1-alpha)*bary_low[2]+alpha*bary_high[2]; + + return {curr_tid,bary}; + } + + static + Face_location + geodesic_lerp(const TriangleMesh &mesh, + const Face_location& src, + const Face_location& tgt,const FT& t +#ifndef CGAL_BSURF_USE_DIJKSTRA_SP + , const Dual_geodesic_solver& solver +#endif + ) + { + std::vector edge_locations; + locally_shortest_path(src,tgt,mesh, edge_locations, solver); + std::vector parameters=path_parameters(edge_locations,mesh,src,tgt); + Face_location point = eval_point_on_geodesic(edge_locations,mesh,src,tgt,parameters,t); + return point; + } + + + static + std::pair,Bezier_segment> + subdivide_bezier_polygon(const TriangleMesh& mesh, + const Bezier_segment& polygon, + const FT& t +#ifndef CGAL_BSURF_USE_DIJKSTRA_SP + , const Dual_geodesic_solver& solver +#endif + ) + { +#ifndef CGAL_BSURF_USE_DIJKSTRA_SP + Face_location Q0 = geodesic_lerp(mesh, polygon[0], polygon[1], t, solver); + Face_location Q1 = geodesic_lerp(mesh, polygon[1], polygon[2], t, solver); + Face_location Q2 = geodesic_lerp(mesh, polygon[2], polygon[3], t, solver); + Face_location R0 = geodesic_lerp(mesh, Q0, Q1, t, solver); + Face_location R1 = geodesic_lerp(mesh, Q1, Q2, t, solver); + Face_location S = geodesic_lerp(mesh, R0, R1, t, solver); +#else + Face_location Q0 = geodesic_lerp(mesh, polygon[0], polygon[1], t); + Face_location Q1 = geodesic_lerp(mesh, polygon[1], polygon[2], t); + Face_location Q2 = geodesic_lerp(mesh, polygon[2], polygon[3], t); + Face_location R0 = geodesic_lerp(mesh, Q0, Q1, t); + Face_location R1 = geodesic_lerp(mesh, Q1, Q2, t); + Face_location S = geodesic_lerp(mesh, R0, R1, t); +#endif + return {{polygon[0], Q0, R0, S}, {S, R1, Q2, polygon[3]}}; + } +}; + + +template +struct Geodesic_circle_impl +{ + using face_descriptor = + typename boost::graph_traits::face_descriptor; + using vertex_descriptor = + typename boost::graph_traits::vertex_descriptor; + using halfedge_descriptor = + typename boost::graph_traits::halfedge_descriptor; + + using Point_2 = typename K::Point_2; + using Point_3 = typename K::Point_3; + using Vector_2 = typename K::Vector_2; + using Vector_3 = typename K::Vector_3; + using FT = typename K::FT; + + using Face_location = Polygon_mesh_processing::Face_location; + using Edge_location = Polygon_mesh_processing::Edge_location; + + struct geodesic_solver { + struct graph_edge { + int node = -1; + double len=DBL_MAX; + }; + std::vector> graph; + }; + + static + void connect_nodes(geodesic_solver &solver, + const vertex_descriptor& a, + const vertex_descriptor& b, + const VertexIndexMap& vidmap, + const FT& len) + { + // TODO: avoid cast + unsigned int vida=get(vidmap,a); + unsigned int vidb=get(vidmap,b); + solver.graph[vida].push_back({static_cast(vidb), len}); + solver.graph[vidb].push_back({static_cast(vida), len}); + + } + + static + double opposite_nodes_arc_length(const VertexPointMap &vpm, + const vertex_descriptor& a, + const vertex_descriptor& c, + const vertex_descriptor& b, + const vertex_descriptor& d) + { + // Triangles (a, b, d) and (b, d, c) are connected by (b, d) edge + // Nodes a and c must be connected. + + Vector_3 ba = get(vpm,a) - get(vpm,b); + Vector_3 bc = get(vpm,c) - get(vpm,b); + Vector_3 bd = get(vpm,d) - get(vpm,b); + + Vector_3 ba_norm=ba/sqrt(ba.squared_length()); + Vector_3 bd_norm=bd/sqrt(bd.squared_length()); + Vector_3 bc_norm=bc/sqrt(bc.squared_length()); + + double cos_alpha = ba_norm * bd_norm; + double cos_beta = bc_norm * bd_norm; + double sin_alpha = sqrt((std::max)(0.0, 1 - cos_alpha * cos_alpha)); + double sin_beta = sqrt((std::max)(0.0, 1 - cos_beta * cos_beta)); + + + // cos(alpha + beta) + double cos_alpha_beta = cos_alpha * cos_beta - sin_alpha * sin_beta; + CGAL_assertion(cos_alpha_beta>-1); + if (cos_alpha_beta <= -1) return DBL_MAX; + + // law of cosines (generalized Pythagorean theorem) + double len = ba.squared_length() + bc.squared_length() - + sqrt(ba.squared_length()) * sqrt(bc.squared_length()) * 2 * cos_alpha_beta; + + CGAL_assertion(len>0); + if (len <= 0) + return DBL_MAX; + else + return sqrt(len); + } + + static + void connect_opposite_nodes(geodesic_solver& solver, + const VertexPointMap &vpm, + const TriangleMesh &mesh, + const VertexIndexMap& vidmap, + const vertex_descriptor& a, + const vertex_descriptor& b, + const halfedge_descriptor& h) + { + vertex_descriptor v0 =target(next(h,mesh),mesh); + vertex_descriptor v1= target(next(opposite(h,mesh),mesh),mesh); + auto len = opposite_nodes_arc_length(vpm, v0,v1,a,b); + //std::cout<0); + connect_nodes(solver,a,b,vidmap,len); + } + face_descriptor nei=face(opposite(h,mesh),mesh); + if(f + make_dual_geodesic_solver(const VertexPointMap &vpm, + const FaceIndexMap& tidmap, + const TriangleMesh &mesh) + { + auto compute_dual_weights=[&mesh,&vpm](const halfedge_descriptor& h) + { + using Impl = Locally_shortest_path_imp; + std::array flat_tid= Impl::init_flat_triangle(h,vpm,mesh); + std::array flat_nei= Impl::unfold_face(h,vpm,mesh,flat_tid); + + Vector_2 c0=0.33*(flat_tid[0]+flat_tid[1]+flat_tid[2]); + Vector_2 c1=0.33*(flat_nei[0]+flat_nei[1]+flat_nei[2]); + + return sqrt((c1 - c0).squared_length()); + }; + + Dual_geodesic_solver solver; + solver.graph.resize(faces(mesh).size()); + for (auto f : faces(mesh)) { + halfedge_descriptor h=halfedge(f,mesh); + int entry=get(tidmap,f); + for (auto i = 0; i < 3; ++i) { + solver.graph[entry][i].node = get(tidmap,face(opposite(h,mesh),mesh)); + solver.graph[entry][i].len = compute_dual_weights(h); + h=next(h,mesh); + } + } + return solver; + } + + // `update` is a function that is executed during expansion, every time a node + // is put into queue. `exit` is a function that tells whether to expand the + // current node or perform early exit. + template + static + void visit_geodesic_graph(std::vector &field, const geodesic_solver &solver, + const std::vector &sources, Update &&update, + Stop &&stop, Exit &&exit) + { + /* + This algorithm uses the heuristic Small Label Fisrt and Large Label Last + https://en.wikipedia.org/wiki/Shortest_Path_Faster_Algorithm + + Large Label Last (LLL): When extracting nodes from the queue, pick the + front one. If it weights more than the average weight of the queue, put + on the back and check the next node. Continue this way. + Sometimes average_weight is less than every value due to floating point + errors (doesn't happen with double precision). + + Small Label First (SLF): When adding a new node to queue, instead of + always pushing it to the end of the queue, if it weights less than the + front node of the queue, it is put on front. Otherwise the node is put at + the end of the queue. + */ + + auto in_queue = std::vector(solver.graph.size(), false); + + // Cumulative weights of elements in queue. Used to keep track of the + // average weight of the queue. + auto cumulative_weight = 0.0; + + // setup queue + auto queue = std::deque(); + for (auto source : sources) { + in_queue[source] = true; + cumulative_weight += field[source]; + queue.push_back(source); + } + + while (!queue.empty()) { + auto node = queue.front(); + auto average_weight = (cumulative_weight / queue.size()); + + // Large Label Last (see comment at the beginning) + for (std::size_t tries = 0; tries < queue.size() + 1; tries++) { + if (field[node] <= average_weight) + break; + queue.pop_front(); + queue.push_back(node); + node = queue.front(); + } + + // Remove node from queue. + queue.pop_front(); + in_queue[node] = false; + cumulative_weight -= field[node]; + + // Check early exit condition. + if (exit(node)) + break; + if (stop(node)) + continue; + + for (auto i = 0; i < (int)solver.graph[node].size(); i++) { + // Distance of neighbor through this node + double new_distance = field[node] + solver.graph[node][i].len; + + auto neighbor = solver.graph[node][i].node; + + auto old_distance = field[neighbor]; + if (new_distance >= old_distance) + continue; + + if (in_queue[neighbor]) { + // If neighbor already in queue, don't add it. + // Just update cumulative weight. + cumulative_weight += new_distance - old_distance; + } else { + // If neighbor not in queue, add node to queue using Small Label + // First (see comment at the beginning). + if (queue.empty() || (new_distance < field[queue.front()])) + queue.push_front(neighbor); + else + queue.push_back(neighbor); + + // Update queue information. + in_queue[neighbor] = true; + cumulative_weight += new_distance; + } + + // Update distance of neighbor. + field[neighbor] = new_distance; + update(node, neighbor, new_distance); + } + } + } + + template + static + void visit_dual_geodesic_graph(std::vector &field, + const Dual_geodesic_solver &solver, + const std::vector &sources, + Update &&update, + Stop &&stop, + Exit &&exit) + { + /* + This algorithm uses the heuristic Small Label Fisrt and Large Label Last + https://en.wikipedia.org/wiki/Shortest_Path_Faster_Algorithm + + Large Label Last (LLL): When extracting nodes from the queue, pick the + front one. If it weights more than the average weight of the queue, put + on the back and check the next node. Continue this way. + Sometimes average_weight is less than every value due to floating point + errors (doesn't happen with double precision). + + Small Label First (SLF): When adding a new node to queue, instead of + always pushing it to the end of the queue, if it weights less than the + front node of the queue, it is put on front. Otherwise the node is put at + the end of the queue. + */ + + auto in_queue = std::vector(solver.graph.size(), false); + + // Cumulative weights of elements in queue. Used to keep track of the + // average weight of the queue. + auto cumulative_weight = 0.0; + + // setup queue + auto queue = std::deque(); + for (auto source : sources) { + in_queue[source] = true; + cumulative_weight += field[source]; + queue.push_back(source); + } + + while (!queue.empty()) { + auto node = queue.front(); + auto average_weight = (cumulative_weight / queue.size()); + + // Large Label Last (see comment at the beginning) + for (std::size_t tries = 0; tries < queue.size() + 1; tries++) { + if (field[node] <= average_weight) + break; + queue.pop_front(); + queue.push_back(node); + node = queue.front(); + } + + // Remove node from queue. + queue.pop_front(); + in_queue[node] = false; + cumulative_weight -= field[node]; + + // Check early exit condition. + if (exit(node)) + break; + if (stop(node)) + continue; + + for (auto i = 0; i < (int)solver.graph[node].size(); i++) { + // Distance of neighbor through this node + auto new_distance = field[node] + solver.graph[node][i].len; + auto neighbor = solver.graph[node][i].node; + + auto old_distance = field[neighbor]; + if (new_distance >= old_distance) + continue; + + if (in_queue[neighbor]) { + // If neighbor already in queue, don't add it. + // Just update cumulative weight. + cumulative_weight += new_distance - old_distance; + } else { + // If neighbor not in queue, add node to queue using Small Label + // First (see comment at the beginning). + if (queue.empty() || (new_distance < field[queue.front()])) + queue.push_front(neighbor); + else + queue.push_back(neighbor); + + // Update queue information. + in_queue[neighbor] = true; + cumulative_weight += new_distance; + } + + // Update distance of neighbor. + field[neighbor] = new_distance; + update(node, neighbor, new_distance); + } + } + } + + static + std::vector + compute_geodesic_distances(const geodesic_solver &solver, + const VertexIndexMap& vidmap, + const std::vector> &sources_and_dist) + { + auto update = [](int /* node */, int /* neighbor */, double /* new_distance */) {}; + auto stop = [](int /* node */) { return false; }; + auto exit = [](int /* node */) { return false; }; + + auto distances = std::vector{}; + distances.assign(solver.graph.size(), DBL_MAX); + std::vectorsources_id((sources_and_dist.size())); + for (std::size_t i = 0; i < sources_and_dist.size(); ++i) { + sources_id[i] = get(vidmap,sources_and_dist[i].first); + distances[sources_id[i]] = sources_and_dist[i].second; + } + visit_geodesic_graph(distances, solver, sources_id, update, stop, exit); + + return distances; + } + + static + std::vector solve_with_targets(const geodesic_solver& solver, + const VertexIndexMap& vidmap, + const std::vector> &sources_and_dist, + const std::vector> &targets_and_dist) + { + auto update = [](int node, int neighbor, float new_distance) {}; + auto stop = [](int node) { return false; }; + double max_distance = DBL_MAX; + std::vector exit_verts(targets_and_dist.size()); + for (auto i = 0; i < targets_and_dist.size(); ++i) { + exit_verts[i]=get(vidmap,targets_and_dist[i].first); + } + auto exit = [&exit_verts](int node) { + auto it = find(exit_verts.begin(), exit_verts.end(), node); + if (it != exit_verts.end()) + exit_verts.erase(it); + + if (exit_verts.empty()) + return true; + + return false; + }; + + auto distances = std::vector(solver.graph.size(), DBL_MAX); + std::vectorsources_id((sources_and_dist.size())); + for (auto i = 0; i < sources_and_dist.size(); ++i) { + sources_id[i] = get(vidmap,sources_and_dist[i].first); + distances[sources_id[i]] = sources_and_dist[i].second; + } + + visit_geodesic_graph(distances, solver, sources_id, update, stop, exit); + return distances; + } + + //compute the length between two opposite vertices by flattening + //TODO: handle concave configurations + static + double length_by_flattening(const VertexPointMap &vpm, + const TriangleMesh &mesh, + const halfedge_descriptor& h) + { + using Impl = Locally_shortest_path_imp; + std::array flat_tid=Impl::init_flat_triangle(h,vpm,mesh); + std::array flat_nei=Impl::unfold_face(h,vpm,mesh,flat_tid); + return sqrt((flat_tid[2]-flat_nei[2]).squared_length()); + } + + static + Point_3 eval_position(const VertexPointMap &vpm, + const TriangleMesh &mesh, + const Face_location& p) + { + halfedge_descriptor h=halfedge(p.first,mesh); + return CGAL::barycenter(get(vpm, source(h,mesh)), p.second[0], get(vpm, target(h,mesh)), p.second[1], get(vpm, target(next(h,mesh),mesh)), p.second[2]); + } + + // compute the distance between a point p and some vertices around him + // TODO: consider to take more vertices (increase accuracy) + // TODO: handle concave configurations + static + std::vector> + nodes_around_point(const VertexPointMap &vpm, + const TriangleMesh &mesh, + const Face_location& p) + { + auto get_vid=[&mesh](const int k,const face_descriptor& tid) + { + halfedge_descriptor h=halfedge(tid,mesh); + switch(k) + { + case 0: + return source(h,mesh); + case 1: + return target(h,mesh); + default: + return target(next(h,mesh),mesh); + } + }; + + std::vector> nodes; + nodes.reserve(6); + auto [is_vert,offset]=Locally_shortest_path_imp::point_is_vert(p); + if (is_vert) { + vertex_descriptor vid = get_vid(offset,p.first); + nodes.push_back({vid, 0}); + } else { + face_descriptor tid = p.first; + Point_3 pos = eval_position(vpm, mesh, p); + halfedge_descriptor h=halfedge(tid,mesh); + for (auto i = 0; i < 3; ++i) { + vertex_descriptor p0 = source(h,mesh); + //connect to current vertex + double d = sqrt(squared_distance(get(vpm,p0),pos)); + nodes.push_back(std::make_pair(p0, d)); + + //connecting to opposite vertex w.r.t to current halfedge + vertex_descriptor opp = target(next(opposite(h,mesh),mesh),mesh); + double l = + length_by_flattening(vpm,mesh,h); + nodes.push_back(std::make_pair(opp, l)); + h=next(h,mesh); + } + } + + return nodes; + } + + //compute geodesic distance field from p + //TODO: can be easiliy extended to more than one source + static + std::vector + compute_geodesic_distances(const geodesic_solver& solver, + const VertexPointMap& vpm, + const VertexIndexMap& vim, + const TriangleMesh &mesh, + const Face_location& p) + { + std::vector> source_nodes=nodes_around_point(vpm,mesh,p); + + return compute_geodesic_distances(solver, vim, source_nodes); + } + + //compute the geodesic distance field from src, and stop the propagation + //once tgt is reached + static + std::vector + compute_pruned_geodesic_distances(const geodesic_solver &solver, + const VertexPointMap &vpm, + const VertexIndexMap& vidmap, + const TriangleMesh &mesh, + const Face_location& src, + const Face_location& tgt) + { + std::vector> + source_nodes = nodes_around_point(vpm,mesh,src); + + std::vector> + target_nodes = nodes_around_point(vpm,mesh,tgt); + + return solve_with_targets(solver, source_nodes, target_nodes); + } + + template + static + std::vector + strip_on_dual_graph(const Dual_geodesic_solver &solver, + const TriangleMesh &mesh, + const int src, + const int tgt) + { + if (src == tgt) + return {}; + + auto common_halfedge = [&mesh](face_descriptor f1, face_descriptor f2) { + halfedge_descriptor h = halfedge(f1, mesh); + for (int i = 0; i < 3; ++i) { + if (face(opposite(h, mesh), mesh) == f2) + return h; + h = next(h, mesh); + } + CGAL_assertion(!"faces do no share a common edge"); + return halfedge_descriptor(); + }; + // initialize once for all and sparsely cleanup at the end of every solve + std::vector parents(solver.graph.size(), -1); + std::vector field(solver.graph.size(), DBL_MAX); + std::vector id_to_face_map(faces(mesh).begin(), faces(mesh).end()); + + field[src]=0.0; + std::vector sources = {src}; + auto update = [&parents](int node, int neighbor, double) { + parents[neighbor] = node; + }; + auto stop = [](int) { return false; }; + auto exit = [tgt](int node) { return node==tgt; }; + + visit_dual_geodesic_graph(field,solver, sources, update, stop, exit); + + // extract_strip + std::vector strip; + int node = tgt; + CGAL_assertion(parents[tgt] != -1); + //update the result using id_to_face_map + strip.reserve((int)std::sqrt(parents.size())); + while (parents[node] != -1) { + strip.push_back(common_halfedge(id_to_face_map[node],id_to_face_map[parents[node]])); + node = parents[node]; + } + std::reverse(strip.begin(),strip.end()); + return strip; + } +}; + +#ifdef CGAL_BSURF_USE_DIJKSTRA_SP + class Dijkstra_end_exception : public std::exception + { + const char* what() const throw () + { + return "Dijkstra shortest path: reached the target vertex."; + } + }; + + template + class Stop_at_target_Dijkstra_visitor : boost::default_dijkstra_visitor + { + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + + vertex_descriptor destination_vd; + + public: + Stop_at_target_Dijkstra_visitor(vertex_descriptor destination_vd) + : destination_vd(destination_vd) + { } + + void initialize_vertex(const vertex_descriptor& /*s*/, const Graph& /*mesh*/) const { } + void examine_vertex(const vertex_descriptor& /*s*/, const Graph& /*mesh*/) const { } + void examine_edge(const edge_descriptor& /*e*/, const Graph& /*mesh*/) const { } + void edge_relaxed(const edge_descriptor& /*e*/, const Graph& /*mesh*/) const { } + void discover_vertex(const vertex_descriptor& /*s*/, const Graph& /*mesh*/) const { } + void edge_not_relaxed(const edge_descriptor& /*e*/, const Graph& /*mesh*/) const { } + void finish_vertex(const vertex_descriptor &vd, const Graph& /* mesh*/) const + { + if(vd == destination_vd) + throw Dijkstra_end_exception(); + } + }; +#endif + + +} // namespace internal + +#ifndef CGAL_BSURF_USE_DIJKSTRA_SP +/*! + * \ingroup VGSMiscellaneous + * Geodesic solver class used to store precomputed information of a given mesh for + * approximate geodesic compution. Those information are computed by the function + * `init_geodesic_dual_solver()`. + * \tparam FT floating point number type (float or double) + */ +template +struct Dual_geodesic_solver +{ + struct Edge { + int node = -1; + FT len = DBL_MAX; + }; + std::vector> graph = {}; +}; + +/*! + * \ingroup VGSMiscellaneous + * fills `solver` for a given mesh `tmesh`. It is the user responsability to + * call again this function if `tmesh` or the points of its vertices are modified. + * If `solver` was used in a previous call to this function, information will be overwritten. + * \tparam TriangleMesh a model of `FaceListGraph` and `EdgeListGraph` + * \tparam FT floating point number type (float or double) + * \param solver the container for the precomputed information + * \param tmesh triangle mesh to be considered for the precomputations + * \todo add named parameters + * \todo make sure solver.graph is cleared before filling it + */ +template +void init_geodesic_dual_solver(Dual_geodesic_solver& solver, const TriangleMesh& tmesh) +{ + //TODO replace with named parameter + using VPM = typename boost::property_map::const_type; + using K = typename Kernel_traits::value_type>::type; + VPM vpm = get(CGAL::vertex_point, tmesh); + typedef typename GetInitializedFaceIndexMap::const_type FIM; + typedef typename GetInitializedVertexIndexMap::const_type VIM; + const FIM fim = get_initialized_face_index_map(tmesh, parameters::default_values()); + + using Impl2 = typename internal::Geodesic_circle_impl; + solver=Impl2::make_dual_geodesic_solver(vpm, fim, tmesh); +} +#endif + +/*! + * \ingroup VGSFunctions + * computes an approximated geodesic shortest path between two locations on a + * `tmesh`. The points`src` and `tgt` must be on the same connected component. + * \tparam TriangleMesh a model of `FaceListGraph` and `EdgeListGraph` + * \tparam FT floating point number type (float or double) + * \tparam EdgeLocationRange a model of `BackInsertionSequence` whose value type `CGAL::Polygon_mesh_processing::Edge_location`. + * \param src source of the path + * \param tgt target of the path + * \param tmesh input triangle mesh to compute the path on + * \param edge_locations contains the path as a sequence of edge locations. + * Two consecutives edges `e_k` and `e_kp1` stored in `edge_locations` are such that + * `face(halfedge(e_k, tmesh), tmesh) == face(opposite(halfedge(e_kp1, tmesh), tmesh), tmesh))`. + * In parcular, it means that if the path goes through a vertex of `tmesh`, several + * edge locations will be reported to maintain this property. + * Additionally, if `src` is in the interior of a face `f`, then the first edge location `e_0` + * of `edge_locations` is such that `f == face(opposite(halfedge(e_0, tmesh), tmesh), tmesh))`. + * Similarly, if `tgt` is in the interior of a face `f`, then the last edge location `e_n` + * of `edge_locations` is such that `f == face(halfedge(e_n, tmesh), tmesh)`. + * \param solver container for the precomputed information. If not initialized, it will be initialized internally. + * \todo add named parameters + * \todo should we have halfedge location instead? + */ +template +void locally_shortest_path(CGAL::Polygon_mesh_processing::Face_location src, + CGAL::Polygon_mesh_processing::Face_location tgt, + const TriangleMesh &tmesh, + EdgeLocationRange &edge_locations +#ifndef CGAL_BSURF_USE_DIJKSTRA_SP + , const Dual_geodesic_solver& solver +#endif +) +{ + namespace PMP = CGAL::Polygon_mesh_processing; + + typedef boost::graph_traits BGT; + typedef typename BGT::halfedge_descriptor halfedge_descriptor; + typedef typename BGT::vertex_descriptor vertex_descriptor; + typedef typename BGT::face_descriptor face_descriptor; + + // start by checking if it is not a trivial path + if (src.first == tgt.first) return; + + auto variant_src = PMP::get_descriptor_from_location(src,tmesh); + auto variant_tgt = PMP::get_descriptor_from_location(tgt,tmesh); + + std::vector src_visible_face; + if (const face_descriptor* f_ptr = std::get_if(&variant_src)) + { + src_visible_face.push_back(*f_ptr); + } + else + { + if (const halfedge_descriptor* h_ptr = std::get_if(&variant_src)) + { + if (!is_border(*h_ptr, tmesh)) + src_visible_face.push_back(face(*h_ptr, tmesh)); + if (!is_border(opposite(*h_ptr, tmesh), tmesh)) + src_visible_face.push_back(face(opposite(*h_ptr, tmesh), tmesh)); + } + else + { + vertex_descriptor v = std::get(variant_src); + for(halfedge_descriptor h : halfedges_around_target(v, tmesh)) + { + if (!is_border(h, tmesh)) src_visible_face.push_back(face(h, tmesh)); + } + } + } + std::sort(src_visible_face.begin(), src_visible_face.end()); + if (const face_descriptor* f_ptr = std::get_if(&variant_tgt)) + { + if (std::find(src_visible_face.begin(), src_visible_face.end(), *f_ptr)!=src_visible_face.end()) + return; + } + else + { + auto is_visible = [&src_visible_face, &tmesh](halfedge_descriptor h) + { + return !is_border(h, tmesh) && + std::find(src_visible_face.begin(), + src_visible_face.end(), face(h, tmesh))!=src_visible_face.end(); + }; + if (const halfedge_descriptor* h_ptr = std::get_if(&variant_tgt)) + { + if (is_visible(*h_ptr) || is_visible(opposite(*h_ptr, tmesh))) return; + } + else + { + vertex_descriptor v = std::get(variant_tgt); + for(halfedge_descriptor h : halfedges_around_target(v, tmesh)) + { + if (is_visible(h)) return; + } + } + } + + + + //TODO replace with named parameter + using VPM = typename boost::property_map::const_type; + using K = typename Kernel_traits::value_type>::type; + using Impl = internal::Locally_shortest_path_imp; + VPM vpm = get(CGAL::vertex_point, tmesh); + + + +// TODO : if (edge(vsrc, vtgt, mesh) || tgt && src on the same edge ) return; + + +//TODO: handle cases of src and tgt not in the same connected component (assert?) + +#ifdef CGAL_BSURF_USE_DIJKSTRA_SP + typedef typename BGT::edge_descriptor edge_descriptor; + typename boost::property_map< + TriangleMesh, CGAL::dynamic_face_property_t>::const_type + predecessor_map = + get(CGAL::dynamic_face_property_t(), tmesh); + typename boost::property_map>::const_type + distance_map = get(CGAL::dynamic_face_property_t(), tmesh); + typename boost::property_map< + TriangleMesh, CGAL::dynamic_edge_property_t>::const_type weight_map = + get(CGAL::dynamic_edge_property_t(), tmesh); + + auto compute_dual_weights=[&tmesh,&vpm](const halfedge_descriptor& h) + { + auto flat_tid= Impl::init_flat_triangle(h,vpm,tmesh); + auto flat_nei= Impl::unfold_face(h,vpm,tmesh,flat_tid); + + Vector_2 c0=0.33*(flat_tid[0]+flat_tid[1]+flat_tid[2]); + Vector_2 c1=0.33*(flat_nei[0]+flat_nei[1]+flat_nei[2]); + + return sqrt((c1 - c0).squared_length()); + }; +//TODO: handle boundary edges + Dual dual(tmesh); + + // TODO: fill the weight map using something better than euclidean distance + // TODO: the edge map could be precomputed if we know that several queries will be done + // TODO: construct the dual graph once at the beginning and then use Dijkstra to + // to navigate it at every query + for (edge_descriptor ed : edges(tmesh)) + { + halfedge_descriptor h=halfedge(ed, tmesh), hopp=opposite(h, tmesh); + + // distance between centroids of unfolded triangles + put(weight_map, ed, compute_dual_weights(h)); + // Euclidean distance between centroids + // put(weight_map, ed, + // sqrt(squared_distance( + // centroid(get(vpm, source(h, tmesh)), get(vpm, target(h, tmesh)), get(vpm, target(next(h, tmesh), tmesh))), + // centroid(get(vpm, source(hopp, tmesh)), get(vpm, target(hopp, tmesh)), get(vpm, target(next(hopp, tmesh), tmesh))) + // ))); + } + + internal::Stop_at_target_Dijkstra_visitor> vis(tgt.first); + + try{ + boost::dijkstra_shortest_paths(dual, src.first, + boost::distance_map(distance_map) + .predecessor_map(predecessor_map) + .weight_map(weight_map) + .visitor(vis)); + } + catch(const internal::Dijkstra_end_exception&) + {} + + // TODO try stopping dijkstra as soon tgt is out of the queue. + + + std::vector initial_path; + + auto common_halfedge = [](face_descriptor f1, face_descriptor f2, + const TriangleMesh &tmesh) { + halfedge_descriptor h = halfedge(f1, tmesh); + for (int i = 0; i < 3; ++i) { + if (face(opposite(h, tmesh), tmesh) == f2) + return h; + h = next(h, tmesh); + } + CGAL_assertion(!"faces do no share a common edge"); + return halfedge_descriptor(); + }; + + face_descriptor current_face = tgt.first; + while (true) { + face_descriptor prev = get(predecessor_map, current_face); + halfedge_descriptor h = common_halfedge(current_face, prev, tmesh); + initial_path.push_back(h); + if (prev == src.first) + break; + current_face = prev; + } + std::reverse(initial_path.begin(), initial_path.end()); +#else + //TODO: VIM is not needed here + typedef typename GetInitializedVertexIndexMap::const_type VIM; + typedef typename GetInitializedFaceIndexMap::const_type FIM; + const FIM fim = get_initialized_face_index_map(tmesh, parameters::default_values()); + + using Impl2 = typename internal::Geodesic_circle_impl; + + std::vector initial_path + = (solver.graph.empty()) + ? Impl2::strip_on_dual_graph(Impl2::make_dual_geodesic_solver(vpm, fim, tmesh), tmesh, get(fim, src.first), get(fim,tgt.first)) + : Impl2::strip_on_dual_graph(solver, tmesh, get(fim, src.first), get(fim,tgt.first)); +#endif + + CGAL_assertion(face(opposite(initial_path.front(), tmesh), tmesh)==src.first); + CGAL_assertion(face(initial_path.back(), tmesh)==tgt.first); + + // remove extra halfedges from inital_path and update src/tgt to get a minimal sequence + // in case src and/or tgt are on a vertex or an edge +#ifdef CGAL_DEBUG_BSURF + std::size_t initial_path_size_before = initial_path.size(); +#endif + Impl::strip_path(tmesh, src, tgt, initial_path); + if (initial_path.empty()) return; + + CGAL_assertion(face(opposite(initial_path.front(), tmesh), tmesh)==src.first); + CGAL_assertion(face(initial_path.back(), tmesh)==tgt.first); +#ifdef CGAL_DEBUG_BSURF + if (initial_path_size_before != initial_path.size()) + { + std::cout << "initial_path cured: " << initial_path_size_before << " ---> " << initial_path.size() << "\n"; + for (halfedge_descriptor h : initial_path) + { + std::cout << " - " << edge(h,tmesh) << "\n"; + } + std::cout << " src=" << src.first << " ("<< src.second[0] << "," << src.second[1] << "," << src.second[2] << ")\n"; + std::cout << " tgt=" << tgt.first << " ("<< tgt.second[0] << "," << tgt.second[1] << "," << tgt.second[2] << ")\n"; + std::cout << " Updated src/tgt " << PMP::construct_point(src, tmesh) << " | " << PMP::construct_point(tgt, tmesh) << "\n"; + } +#endif + + // here portals contains 2D coordinates of endpoints of edges in initial_path + std::vector< std::array> portals=Impl::unfold_strip(initial_path,src,tgt,vpm,tmesh); + std::size_t max_index=0; + + // lerps are barycentric coordinates of the shortest path restricted to the unfolded strip + std::vector lerps=Impl::funnel(portals,max_index); + // TODO: if you comment this if you don't want to shorten the path (option?). + // but this part is really fast so maybe does not make sense. + Impl::straighten_path(portals,lerps,initial_path,src,tgt,vpm,tmesh,max_index); + CGAL_assertion(lerps.size()==initial_path.size()); + + edge_locations.reserve(initial_path.size()); + for(std::size_t i=0; i +std::vector> +recursive_de_Casteljau(const TriangleMesh& mesh, + const Bezier_segment& control_points, + const int num_subdiv +#ifndef CGAL_BSURF_USE_DIJKSTRA_SP + , const Dual_geodesic_solver& solver = Dual_geodesic_solver() +#endif + ) +{ + namespace PMP = CGAL::Polygon_mesh_processing; + + //TODO replace with named parameter + using VPM = typename boost::property_map::const_type; + using K = typename Kernel_traits::value_type>::type; + using Impl = internal::Bezier_tracing_impl; + + // init solver if empty + const Dual_geodesic_solver* solver_ptr=&solver; + Dual_geodesic_solver local_solver; + if (solver.graph.empty()) + { + solver_ptr = &local_solver; + init_geodesic_dual_solver(local_solver, mesh); + } + + std::vector> segments(1,control_points); + std::vector> result; + for (auto subdivision = 0; subdivision < num_subdiv; ++subdivision) + { + result.clear(); + result.reserve(segments.size() * 2); + for (std::size_t i = 0; i < segments.size(); ++i) + { + auto [split0, split1] = Impl::subdivide_bezier_polygon(mesh, segments[i], 0.5, *solver_ptr); + result.push_back(split0); + result.push_back(split1); + } + + CGAL_assertion( 2*segments.size()==result.size() ); + + std::swap(segments, result); + } + + std::vector> final_result; + final_result.reserve(result.size() * 3 + 1); + final_result.push_back(std::move(result.front()[0])); + for (Bezier_segment& bs : result) + { + for (int i=1;i<4;++i) + final_result.push_back(std::move(bs[i])); + } + +#if 0 + // nasty trick to build the vector from a pair of iterators + // using the fact that data in array and vector are contiguous + return {(PMP::Face_location*)segments.data(), + (PMP::Face_location*)segments.data() + segments.size() * 4}; +#endif + return final_result; +} + +/*! + * \ingroup VGSMiscellaneous + * computes the approximate geodesic distances of `center` to all the vertices of the mesh (only one connected component is expected) + * and put the distance in `distance_map` + */ +template +void approximate_geodesic_distance_field(const CGAL::Polygon_mesh_processing::Face_location& center, + VertexDistanceMap distance_map, + const TriangleMesh& tmesh) +{ + // TODO: the solver could be init once and used several times for different centers + // in particular, it can be tweaked to compute the Voronoi diagram of the initial centers + // or geodesic furthest point sampling. + // TODO: add a parameter for the link size to increase to precision of the approximation of the distance + // that is you increase the size of the neighborhood of each vertex and you connect in the graph each vertex to its neighbors + // with weight being the geodesic shortest path. + // the more neighbors you have, the better is the approximation + // TODO: graph construction can be done in parallel + // TODO: concave flattening should be handled to improve the approximation of the distance + // (shortest path is not a line in that case) + + //TODO replace with named parameter + using VPM = typename boost::property_map::const_type; + VPM vpm = get(CGAL::vertex_point, tmesh); + using K = typename Kernel_traits::value_type>::type; + + typedef typename GetInitializedVertexIndexMap::const_type VIM; + const VIM vim = get_initialized_vertex_index_map(tmesh, parameters::default_values()); + typedef typename GetInitializedFaceIndexMap::const_type FIM; + + using Impl = typename internal::Geodesic_circle_impl; + + typename Impl::geodesic_solver solver = Impl::make_geodesic_solver(vpm, vim,tmesh); + std::vector distances = Impl::compute_geodesic_distances(solver, vpm, vim, tmesh, center); + + for (typename Impl::vertex_descriptor v : vertices(tmesh)) + { + put(distance_map, v, distances[get(vim,v)]); + } +} + +/*! + * \ingroup VGSFunctions + * computes a path on a triangle mesh that is computed by starting a walk on `tmesh` + * given a direction and a maximum distance. The distance will not be achieved if a border edge + * is reached before. + * \tparam TriangleMesh a model of `FaceListGraph` and `EdgeListGraph` + * \tparam FT floating point number type (float or double) + * \param tmesh input triangle mesh to compute the path on + * \param src the source of the path + * \param len the distance to walk along the straightest + * \param dir the initial direction of the walk, given as a 2D vector in the face of src, the halfedge of the face being the y-axis. + * \return the straightest path (not containing `src`) + * \todo add named parameters + * \todo do we want to also have a way to return Bézier segments? The output is actually Bézier segments subdivided. + * \todo offer something better than a 2D vector for the direction + */ +template +std::vector> +straightest_geodesic(const CGAL::Polygon_mesh_processing::Face_location &src, + const typename K::Vector_2& dir, + const typename K::FT len, + const TriangleMesh &tmesh) +{ + //TODO replace with named parameter + using VPM = typename boost::property_map::const_type; + using Impl = internal::Locally_shortest_path_imp; + VPM vpm = get(CGAL::vertex_point, tmesh); + + + return Impl::straightest_geodesic(src, tmesh, vpm, dir, len); +} + + +/*! + * \ingroup VGSMiscellaneous + * converts the coordinates of a range of points into polar coordinates with respect to a given center + * \tparam K a model of `PolygonTraits_2` + * \tparam PointRange_2 a model of the concept `RandomAccessContainer` with `K::Point_2` as value type + * \param points the input points to convert. + * \param center the point of reference for the polar coordinates. + * If omitted, then centroid of `polygon` will be used. + * \return the polar coordinates of the points as a pair (distance, angle) + */ +template +std::vector> +convert_to_polar_coordinates(const PointRange_2& points, + std::optional center = std::nullopt) +{ + std::vector> result; + // compute naively the center + if (center==std::nullopt) + { + Bbox_2 bb = bbox_2(points.begin(), points.end()); + center = typename K::Point_2((bb.xmax()+bb.xmin())/2., (bb.ymax()+bb.ymin())/2); + } + + bool is_closed = points.front()==points.back(); + std::size_t nbp=points.size(); + if (is_closed) --nbp; + for (std::size_t i=0; iy())/* /d */, (pt.x()-center->x())/* /d */); + result.emplace_back(d, polar); + // std::cout << center->x()+d*std::cos(polar) << " " << center->x()+d*std::sin(polar) << " 0\n"; + } + if (is_closed) result.push_back(result.front()); + + return result; +} + + + +/*! + * \ingroup VGSFunctions + * computes the face location of each vertex of a 2D polygon on `tmesh`. + * The vertices of the polygon are given as a pair of direction and distance with respect to a point corresponding to `center` on `tmesh`. + * \tparam TriangleMesh a model of `FaceListGraph` and `EdgeListGraph` + * \tparam K a model of `Kernel` with `K::FT` being a floating point number type (float or double) + * \param center the location on `tmesh` used as reference. The y-axis used for coordinates is `halfedge(center.first, tmesh)`. + * \param directions contains the direction one need to move from `center` to reach each vertex of the polygon. + * \param lengths the distance one need to move from `center` along the direction at the same position in `directions` to reach each vertex of the polygon. + * \param tmesh input triangle mesh supporting the vertices of the output polygon + * \param solver container for the precomputed information. If not initialized, it will be initialized internally. + * \return a face location for each vertex of the polygon + * \todo add named parameters + * \todo polygon orientation is not handled in the function and should be done outside of the function for now + * \todo offer something better than a 2D vector for the direction + * \todo directly handle polyline? + * \todo why the first polygon vertex is duplicated by the function? (most probably for the example but it shouldn't be done here) + */ +template +std::vector> +trace_geodesic_polygon(const CGAL::Polygon_mesh_processing::Face_location ¢er, + const std::vector& directions, + const std::vector& lengths, + const TriangleMesh &tmesh, + const Dual_geodesic_solver& solver = {}) +{ + namespace PMP = CGAL::Polygon_mesh_processing; + + std::size_t n=directions.size(); + std::vector> result; + std::vector> vertices(n); +#ifdef CGAL_DEBUG_BSURF + std::cout << "trace_geodesic_polygon\n"; +#endif + for(std::size_t i=0;i(center,directions[i],lengths[i],tmesh).back(); +#ifdef CGAL_DEBUG_BSURF + std::cout << PMP::construct_point(vertices[i], tmesh) << "\n"; +#endif + } + std::vector> edge_locations; + + const Dual_geodesic_solver* solver_ptr=&solver; + Dual_geodesic_solver local_solver; + if (solver.graph.empty()) + { + solver_ptr = &local_solver; + init_geodesic_dual_solver(local_solver, tmesh); + } + + for(std::size_t i=0;i(vertices[i],vertices[(i+1)%n],tmesh, edge_locations, *solver_ptr); + result.push_back(vertices[i]); + for(auto& el : edge_locations) + { + result.push_back(PMP::to_face_location(el, tmesh)); + } + } + result.push_back(vertices[0]); + + return result; +} + +/*! + * \ingroup VGSFunctions + * computes for each vertex of each polygon in `polygons` a face location on `tmesh`, where `center` represents the center of the 2D bounding box of the polygons. + * This method computes the location of the center of the bounding box of each polygon on the mesh with respect to `center` and calls `trace_geodesic_polygon()` with that center with + * appropriate directions and distances to have a consistent orientation for the polygons. + * \tparam TriangleMesh a model of `FaceListGraph` and `EdgeListGraph` + * \tparam K a model of `Kernel` with `K::FT` being a floating point number type (float or double) + * \param center the location on `tmesh` corresponding to the center of the 2D bounding box of the polygons. + * \param polygons 2D polygons + * \param scaling a scaling factor to scale the polygons on `tmesh` (considering geodesic distances on `tmesh`) + * \param tmesh input triangle mesh supporting the vertices of the output polygon + * \param solver container for the precomputed information. If not initialized, it will be initialized internally. + * \return a face location for each vertex of each polygon + * \todo add named parameters + * \todo polygon orientation is not handled in the function and should be done outside of the function for now + * \todo for better rendering we can group polygons to have one center for the same group of polygon (useful for letters that are not simply connected) + * \todo what if boundary is reached + */ +template +std::vector>> +trace_geodesic_polygons(const CGAL::Polygon_mesh_processing::Face_location ¢er, + const std::vector>& polygons, + const typename K::FT scaling, + const TriangleMesh &tmesh, + const Dual_geodesic_solver& solver = {}) +{ + namespace PMP = CGAL::Polygon_mesh_processing; + + using VPM = typename boost::property_map::const_type; + using Impl = internal::Locally_shortest_path_imp; + VPM vpm = get(CGAL::vertex_point, tmesh); + + + using halfedge_descriptor = typename boost::graph_traits::halfedge_descriptor; + std::vector polygon_bboxes; + polygon_bboxes.reserve(polygons.size()); + Bbox_2 gbox; + for (const std::vector& polygon : polygons) + { + polygon_bboxes.push_back( bbox_2(polygon.begin(), polygon.end()) ); + gbox+=polygon_bboxes.back(); + } + + const Dual_geodesic_solver* solver_ptr=&solver; + Dual_geodesic_solver local_solver; + if (solver.graph.empty()) + { + solver_ptr = &local_solver; + init_geodesic_dual_solver(local_solver, tmesh); + } + + std::vector>> result(polygons.size()); + typename K::Vector_2 initial_dir(1,0);// TODO: input parameter or 2D PCA of the centers? + + for(std::size_t i=0;i > spath = + straightest_geodesic(center, dir, scaling * std::sqrt(dir.squared_length()),tmesh); + PMP::Face_location polygon_center = spath.back(); + + // TODO: avoid using the shortest path and directly use the straightest! + std::vector> shortest_path; + locally_shortest_path(center, polygon_center, tmesh, shortest_path, *solver_ptr); + + // update direction + typename K::Vector_2 v = initial_dir; + for(std::size_t i=0;i> polar_coords = + convert_to_polar_coordinates(polygons[i], + typename K::Point_2((polygon_bboxes[i].xmin()+polygon_bboxes[i].xmax())/2., + (polygon_bboxes[i].ymin()+polygon_bboxes[i].ymax())/2.)); + if (polygons[i].front()==polygons[i].back()) + polar_coords.pop_back(); + + double theta = atan2(v.y(),v.x()); + + std::vector directions; + std::vector lens; + lens.reserve(polar_coords.size()); + directions.reserve(polar_coords.size()); + + for (const std::pair& coord : polar_coords) + { + lens.push_back(scaling * coord.first); + directions.emplace_back(std::cos(coord.second+theta), std::sin(coord.second+theta)); + } + result[i] = trace_geodesic_polygon(polygon_center, directions, lens, tmesh, *solver_ptr); + } + + return result; +} + + + +/*! + * \ingroup VGSFunctions + * computes for each vertex of each polygon in `polygons` a face location on `tmesh`, where `center` represents the center of the 2D bounding box of the polygons. + * This method starts by considering the segment splitting in two halves along the y-axis the bounding box of the polygons. 2D centers for each polygon are + * computed on this segment as the intersection with the line splitting the bounding box of the polygon in two halves along the x-axis. + * The splitting segment is then drawn on `tmesh` and the face location of the 2D centers is found. + * `trace_geodesic_polygon()` is then called for each polygon and center, with appropriate directions and distances to have a consistent orientation for the polygons. + * \tparam TriangleMesh a model of `FaceListGraph` and `EdgeListGraph` + * \tparam K a model of `Kernel` with `K::FT` being a floating point number type (float or double) + * \param center the location on `tmesh` corresponding to the center of the 2D bounding box of the polygons. + * \param polygons 2D polygons + * \param scaling a scaling factor to scale the polygons on `tmesh` (considering geodesic distances on `tmesh`) + * \param tmesh input triangle mesh supporting the vertices of the output polygon + * \param solver container for the precomputed information. If not initialized, it will be initialized internally. + * \return a face location for each vertex of each polygon + * \todo add named parameters + * \todo polygon orientation is not handled in the function and should be done outside of the function for now + * \todo for better rendering we can group polygons to have one center for the same group of polygon (useful for letters that are not simply connected) + * \todo what if boundary is reached + */ +template +std::vector>> +trace_geodesic_label(const CGAL::Polygon_mesh_processing::Face_location ¢er, + const std::vector>& polygons, + const typename K::FT scaling, + const TriangleMesh &tmesh, + const Dual_geodesic_solver& solver = {}) +{ + namespace PMP = CGAL::Polygon_mesh_processing; + + using VPM = typename boost::property_map::const_type; + using Impl = internal::Locally_shortest_path_imp; + VPM vpm = get(CGAL::vertex_point, tmesh); + using Point_2 = typename K::Point_2; + using Vector_2 = typename K::Vector_2; + + + std::vector polygon_bboxes; + polygon_bboxes.reserve(polygons.size()); + Bbox_2 gbox; + for (const std::vector& polygon : polygons) + { + polygon_bboxes.push_back( bbox_2(polygon.begin(), polygon.end()) ); + gbox+=polygon_bboxes.back(); + } + + const Dual_geodesic_solver* solver_ptr=&solver; + Dual_geodesic_solver local_solver; + if (solver.graph.empty()) + { + solver_ptr = &local_solver; + init_geodesic_dual_solver(local_solver, tmesh); + } + + std::vector>> result(polygons.size()); + // Vector_2 initial_dir(1,0);// TODO: input parameter or 2D PCA of the centers? + + // 1D partition of the letters + Point_2 c2( (gbox.xmin()+gbox.xmax())/2., + (gbox.ymin()+gbox.ymax())/2. ); + Point_2 left_most(gbox.xmin(), c2.y()); + Point_2 right_most(gbox.xmax(), c2.y()); + double len = (gbox.xmax()-gbox.xmin())/2.; + + std::vector< PMP::Face_location > left_path = + straightest_geodesic(center, left_most-c2, scaling * len, tmesh); + std::vector< PMP::Face_location > right_path = + straightest_geodesic(center, right_most-c2, scaling * len, tmesh); + + CGAL_assertion(left_path.size() >=2); + CGAL_assertion(right_path.size() >=2); + + // TODO: precompute distances along supporting curve + stop when exceeding max distance needed + + for(std::size_t i=0;i polygon_center; + double xc = (polygon_bboxes[i].xmin()+polygon_bboxes[i].xmax())/2.; + + auto get_polygon_center = [&tmesh, &vpm](const std::vector>& path, + double targetd) + { + // use left + double acc=0.; + std::size_t k=0; + while(true) + { + double len = std::sqrt(squared_distance(PMP::construct_point(path[k], tmesh), + PMP::construct_point(path[k+1], tmesh))); + acc+=len; + if (acc == targetd) + { + // TODO: if you land here and path[k] is on an input vertex, it might be that loc_k==loc_k1 + + double theta=0; + PMP::Face_location loc_k=path[k], loc_k1=path[k+1]; + CGAL_assertion_code(bool OK=) + PMP::locate_in_common_face(loc_k, loc_k1, tmesh); + CGAL_assertion(OK); + std::array flat_triangle = + Impl::init_flat_triangle(halfedge(loc_k.first,tmesh),vpm,tmesh); + Vector_2 src = loc_k.second[0]*flat_triangle[0]+loc_k.second[1]*flat_triangle[1]+loc_k.second[2]*flat_triangle[2]; + Vector_2 tgt = loc_k1.second[0]*flat_triangle[0]+loc_k1.second[1]*flat_triangle[1]+loc_k1.second[2]*flat_triangle[2]; + + Vector_2 dir2 = tgt-src; + theta = atan2(dir2.y(), dir2.x()); + + return std::make_pair(path[k+1], theta); + } + + if (acc > targetd) + { + double excess = acc-targetd; + + PMP::Face_location loc_k=path[k], loc_k1=path[k+1]; + CGAL_assertion_code(bool OK=) + PMP::locate_in_common_face(loc_k, loc_k1, tmesh); + CGAL_assertion(OK); + + PMP::Face_location polygon_center; + polygon_center.first=loc_k.first; + double alpha = excess/len; + + for(int ii=0; ii<3;++ii) + polygon_center.second[ii] = loc_k.second[ii]*alpha+loc_k1.second[ii]*(1.-alpha); + + std::array flat_triangle = + Impl::init_flat_triangle(halfedge(polygon_center.first,tmesh),vpm,tmesh); + Vector_2 src = loc_k.second[0]*flat_triangle[0]+loc_k.second[1]*flat_triangle[1]+loc_k.second[2]*flat_triangle[2]; + Vector_2 tgt = loc_k1.second[0]*flat_triangle[0]+loc_k1.second[1]*flat_triangle[1]+loc_k1.second[2]*flat_triangle[2]; + + Vector_2 dir2 = tgt-src; + double theta = atan2(dir2.y(), dir2.x()); + + return std::make_pair(polygon_center, theta); + } + + if (++k==path.size()-1) + { + PMP::Face_location loc_k=path[k-1], loc_k1=path[k]; + CGAL_assertion_code(bool OK=) + PMP::locate_in_common_face(loc_k, loc_k1, tmesh); + CGAL_assertion(OK); + std::array flat_triangle = + Impl::init_flat_triangle(halfedge(loc_k.first,tmesh),vpm,tmesh); + Vector_2 src = loc_k.second[0]*flat_triangle[0]+loc_k.second[1]*flat_triangle[1]+loc_k.second[2]*flat_triangle[2]; + Vector_2 tgt = loc_k1.second[0]*flat_triangle[0]+loc_k1.second[1]*flat_triangle[1]+loc_k1.second[2]*flat_triangle[2]; + + Vector_2 dir2 = tgt-src; + double theta = atan2(dir2.y(), dir2.x()); + + return std::make_pair(path.back(), theta); + } + } + }; + double theta; + if (xc> shortest_path; + locally_shortest_path(center, polygon_center, tmesh, shortest_path, *solver_ptr); + + + std::vector> polar_coords = + convert_to_polar_coordinates(polygons[i], + typename K::Point_2((polygon_bboxes[i].xmin()+polygon_bboxes[i].xmax())/2., + (gbox.ymin()+gbox.ymax())/2.)); + + //already duplicated in trace_geodesic_polygon + if (polar_coords.front()==polar_coords.back()) polar_coords.pop_back(); + + std::vector directions; + std::vector lens; + lens.reserve(polar_coords.size()); + directions.reserve(polar_coords.size()); + + for (const std::pair& coord : polar_coords) + { + lens.push_back(scaling * coord.first); + directions.emplace_back(std::cos(coord.second+theta), std::sin(coord.second+theta)); + } + result[i] = trace_geodesic_polygon(polygon_center, directions, lens, tmesh, *solver_ptr); + } + + return result; +} + + +/*! + * \ingroup VGSMiscellaneous + * computes the length of a path on a triangle mesh. + * \tparam FT floating point number type (float or double) + * \tparam TriangleMesh a model of `FaceGraph` + * \param path a path described as a range of face locations, with the property that + for two consecutive face locations, there exists a face in `tmesh` containing the two corresponding points. + * \param tmesh the triangle mesh supporing the path + * \todo add named parameters + * \todo generic range + */ +template +FT path_length(const std::vector>& path, + const TriangleMesh &tmesh) +{ + namespace PMP = CGAL::Polygon_mesh_processing; + + std::size_t lpath = path.size(); + if(lpath<2) + return 0; + + using VPM = typename boost::property_map::const_type; + VPM vpm = get(CGAL::vertex_point, tmesh); + + FT len(0); + + for (std::size_t i=0; i +FT path_length(const std::vector>& path, + const CGAL::Polygon_mesh_processing::Face_location& src, + const CGAL::Polygon_mesh_processing::Face_location& tgt, + const TriangleMesh &tmesh) +{ + namespace PMP = CGAL::Polygon_mesh_processing; + + std::size_t lpath = path.size(); + if(lpath==0) + return sqrt(squared_distance(PMP::construct_point(src,tmesh), + PMP::construct_point(tgt,tmesh))); + + using VPM = typename boost::property_map::const_type; + VPM vpm = get(CGAL::vertex_point, tmesh); + + FT len=sqrt(squared_distance(PMP::construct_point(src,tmesh), + PMP::construct_point(path[0],tmesh))); + + for (std::size_t i=0; i +std::vector>> +trace_geodesic_label_along_curve(const std::vector>& supporting_curve, + const std::vector>& polygons, + const typename K::FT scaling, + const typename K::FT padding, + const bool is_centered, + const TriangleMesh &tmesh, + const Dual_geodesic_solver& solver = {}) +{ + namespace PMP = CGAL::Polygon_mesh_processing; + + using VPM = typename boost::property_map::const_type; + using Impl = internal::Locally_shortest_path_imp; + VPM vpm = get(CGAL::vertex_point, tmesh); + using Point_2 = typename K::Point_2; + using Vector_2 = typename K::Vector_2; + using Face_loc = PMP::Face_location; + + std::vector polygon_bboxes; + polygon_bboxes.reserve(polygons.size()); + Bbox_2 gbox; + for (const std::vector& polygon : polygons) + { + polygon_bboxes.push_back( bbox_2(polygon.begin(), polygon.end()) ); + gbox+=polygon_bboxes.back(); + } + + const Dual_geodesic_solver* solver_ptr=&solver; + Dual_geodesic_solver local_solver; + if (solver.graph.empty()) + { + solver_ptr = &local_solver; + init_geodesic_dual_solver(local_solver, tmesh); + } + + std::vector> result(polygons.size()); +// Vector_2 initial_dir(1,0);// TODO: input parameter or 2D PCA of the centers? + + // 1D partition of the letters + Point_2 c2( (gbox.xmin()+gbox.xmax())/2., + (gbox.ymin()+gbox.ymax())/2. ); + Point_2 left_most(gbox.xmin(), c2.y()); + Point_2 right_most(gbox.xmax(), c2.y()); + + std::vector support_len(supporting_curve.size(),0); + for (std::size_t i=0; ipadding + scaling * (gbox.xmax()-gbox.xmin())) || + (is_centered && (support_len.back()> scaling * (gbox.xmax()-gbox.xmin()))) ); + + typename K::FT pad = is_centered ? (support_len.back()-(scaling*(gbox.xmax()-gbox.xmin())))/2 + : padding; + + auto get_polygon_center = [&tmesh, &vpm, &supporting_curve, &support_len](double targetd) + { + // use left + std::size_t k=0; + while(true) // TODO get rid of the while and std::lower_bound + { + double acc = support_len[k]; + if (acc == targetd) + { + // note: if the supporting curve goes though supporting_curve[k], that vertex might be duplicated (if it's a path). + // But it is not an issue for the code as the 0 length segments do not contribute to the distance and + // points of loc_k and loc_k1 are expected to be always different. + + double theta=0; + //TODO: should pick k-1 if k==supporting_curve.size()-1 + PMP::Face_location loc_k=supporting_curve[k], loc_k1=supporting_curve[k+1]; + CGAL_assertion_code(bool OK=) + PMP::locate_in_common_face(loc_k, loc_k1, tmesh); + CGAL_assertion(OK); + std::array flat_triangle = + Impl::init_flat_triangle(halfedge(loc_k.first,tmesh),vpm,tmesh); + Vector_2 src = loc_k.second[0]*flat_triangle[0]+loc_k.second[1]*flat_triangle[1]+loc_k.second[2]*flat_triangle[2]; + Vector_2 tgt = loc_k1.second[0]*flat_triangle[0]+loc_k1.second[1]*flat_triangle[1]+loc_k1.second[2]*flat_triangle[2]; + + Vector_2 dir2 = tgt-src; + theta = atan2(dir2.y(), dir2.x()); + + return std::make_pair(supporting_curve[k], theta); + } + + if (acc > targetd) + { + double excess = acc-targetd; + + PMP::Face_location loc_k=supporting_curve[k-1], loc_k1=supporting_curve[k]; + CGAL_assertion_code(bool OK=) + PMP::locate_in_common_face(loc_k, loc_k1, tmesh); + CGAL_assertion(OK); + + PMP::Face_location polygon_center; + polygon_center.first=loc_k.first; + double alpha = excess/(support_len[k]-support_len[k-1]); + + for(int ii=0; ii<3;++ii) + polygon_center.second[ii] = loc_k.second[ii]*alpha+loc_k1.second[ii]*(1.-alpha); + + std::array flat_triangle = + Impl::init_flat_triangle(halfedge(polygon_center.first,tmesh),vpm,tmesh); + Vector_2 src = loc_k.second[0]*flat_triangle[0]+loc_k.second[1]*flat_triangle[1]+loc_k.second[2]*flat_triangle[2]; + Vector_2 tgt = loc_k1.second[0]*flat_triangle[0]+loc_k1.second[1]*flat_triangle[1]+loc_k1.second[2]*flat_triangle[2]; + + Vector_2 dir2 = tgt-src; + double theta = atan2(dir2.y(), dir2.x()); + + return std::make_pair(polygon_center, theta); + } + + if (++k==supporting_curve.size()) + { + //TODO: shall we throw an exception instead or return false? + + PMP::Face_location loc_k=supporting_curve[k-2], loc_k1=supporting_curve[k-1]; + CGAL_assertion_code(bool OK=) + PMP::locate_in_common_face(loc_k, loc_k1, tmesh); + CGAL_assertion(OK); + std::array flat_triangle = + Impl::init_flat_triangle(halfedge(loc_k.first,tmesh),vpm,tmesh); + Vector_2 src = loc_k.second[0]*flat_triangle[0]+loc_k.second[1]*flat_triangle[1]+loc_k.second[2]*flat_triangle[2]; + Vector_2 tgt = loc_k1.second[0]*flat_triangle[0]+loc_k1.second[1]*flat_triangle[1]+loc_k1.second[2]*flat_triangle[2]; + + Vector_2 dir2 = tgt-src; + double theta = atan2(dir2.y(), dir2.x()); + + return std::make_pair(supporting_curve.back(), theta); + } + } + }; + + for(std::size_t i=0;i polygon_center; + double xc = (polygon_bboxes[i].xmin()+polygon_bboxes[i].xmax())/2.; + + double theta; + std::tie(polygon_center, theta)=get_polygon_center(scaling * (xc-gbox.xmin())+pad); + + std::vector> polar_coords = + convert_to_polar_coordinates(polygons[i], + typename K::Point_2((polygon_bboxes[i].xmin()+polygon_bboxes[i].xmax())/2., + (gbox.ymin()+gbox.ymax())/2.)); + //already duplicated in trace_geodesic_polygon + if (polar_coords.front()==polar_coords.back()) polar_coords.pop_back(); + + std::vector directions; + std::vector lens; + lens.reserve(polar_coords.size()); + directions.reserve(polar_coords.size()); + + for (const std::pair& coord : polar_coords) + { + lens.push_back(scaling * coord.first); + directions.emplace_back(std::cos(coord.second+theta), std::sin(coord.second+theta)); + } + result[i] = trace_geodesic_polygon(polygon_center, directions, lens, tmesh, *solver_ptr); + } + + return result; +} + +/*! + * \ingroup VGSFunctions + * computes a path representing a Bézier curve defined by four control points. + * Control points are defined by the endpoints of straightest geodesic curves starting from `center` along given directions and distances. + * The iterative de Casteljau subdivision algorithm is applied to create more control points that are then connected with locally shortest paths. + * The output path is such that for two consecutive face locations, there must exist a face in `tmesh` containing the two corresponding points. + * \tparam TriangleMesh a model of `FaceListGraph` and `EdgeListGraph` + * \tparam K a model of `Kernel` with `K::FT` being a floating point number type (float or double) + * \param center the location on `tmesh` where straightest geodesic for the placement of control points starts. The y-axis used is `halfedge(center.first, tmesh)`. + * \param directions contains the direction of the straightest geodesic for each control point + * \param lengths contains the length of the straightest geodesic for each control point + * \param num_subdiv the number of iterations of the de Casteljau subdivision algorithm + * \param tmesh input triangle mesh supporting the vertices of the output polygon + * \param solver container for the precomputed information. If not initialized, it will be initialized internally. + * \return a face location for each vertex of each polygon + * \todo add named parameters + * \todo polygon orientation is not handled in the function and should be done outside of the function for now + * \todo for better rendering we can group polygons to have one center for the same group of polygon (useful for letters that are not simply connected) + * \todo what if boundary is reached + */ +template +std::vector< std::vector> > +trace_bezier_curves(const CGAL::Polygon_mesh_processing::Face_location ¢er, + const std::vector>& directions, + const std::vector>& lengths, + const int num_subdiv, + const TriangleMesh &tmesh +#ifndef CGAL_BSURF_USE_DIJKSTRA_SP + , const Dual_geodesic_solver& solver = {} +#endif +) +{ + namespace PMP = CGAL::Polygon_mesh_processing; + + using FT = typename K::FT; + + std::size_t n=directions.size(); + std::vector< std::vector> > result(n); + +#ifndef CGAL_BSURF_USE_DIJKSTRA_SP + const Dual_geodesic_solver* solver_ptr=&solver; + Dual_geodesic_solver local_solver; + if (solver.graph.empty()) + { + solver_ptr = &local_solver; + init_geodesic_dual_solver(local_solver, tmesh); + } +#endif + +#ifdef CGAL_DEBUG_BSURF + std::ofstream debug_cp("/tmp/control_points.xyz"); + std::ofstream debug_ep("/tmp/end_points.xyz"); + debug_cp << std::setprecision(17); + debug_ep << std::setprecision(17); +#endif + for (std::size_t i=0; i control_loc; + for (int k=0;k<4; ++k) + { + control_loc[k] = straightest_geodesic(center,directions[i][k],lengths[i][k],tmesh).back(); + } + + #ifdef CGAL_DEBUG_BSURF + debug_ep << PMP::construct_point(control_loc[0], tmesh) << "\n"; + debug_ep << PMP::construct_point(control_loc[3], tmesh) << "\n"; + debug_cp << PMP::construct_point(control_loc[1], tmesh) << "\n"; + debug_cp << PMP::construct_point(control_loc[2], tmesh) << "\n"; + #endif + + std::vector> bezier = + recursive_de_Casteljau(tmesh, control_loc, num_subdiv +#ifndef CGAL_BSURF_USE_DIJKSTRA_SP + , *solver_ptr +#endif + ); + + result[i].reserve(bezier.size()); + result[i].push_back(bezier[0]); + for(std::size_t b=0; b& loc = bezier[b]; + const PMP::Face_location& loc1 = bezier[b+1]; + + // connect the two face locations with shortest path is they are in different faces + if (loc.first!=loc1.first) + { + std::vector> edge_locations; + locally_shortest_path(loc, loc1, tmesh, edge_locations, solver); + result[i].reserve(result[i].size()+edge_locations.size()); + for (const PMP::Edge_location& e : edge_locations) + result[i].push_back(PMP::to_face_location(e, tmesh)); + } + result[i].push_back(loc1); + } + } + + return result; +} + + +/*! + * \ingroup VGSFunctions + * computes a path representing a Bézier polyline (a sequence of Bézier curves having an identical control points, + * that is the fourth control point of the nth curve is the first control point of the (n+1)th curve). + * Control points are defined by the endpoints of straightest geodesic curves starting from `center` along given directions and distances. + * The iterative de Casteljau subdivision algorithm is applied to create more control points that are then connected with locally shortest paths. + * The output path is such that for two consecutive face locations, there must exist a face in `tmesh` containing the two corresponding points. + * The first portion of the curve is defined by the 4 first values in `directions` and `lengths`. The second portion is defined by the 4'th value + * and the next 3, and so on until the end is reached. If `is_closed` is true, then the first value will be used with the last three to define the last portion. + * \tparam TriangleMesh a model of `FaceListGraph` and `EdgeListGraph` + * \tparam K a model of `Kernel` with `K::FT` being a floating point number type (float or double) + * \param center the location on `tmesh` where straightest geodesic for the placement of control points starts. The y-axis used is `halfedge(center.first, tmesh)`. + * \param directions contains the direction of the straightest geodesic for each control point + * \param lengths contains the length of the straightest geodesic for each control point + * \param num_subdiv the number of iterations of the de Casteljau subdivision algorithm + * \param is_closed if true [directions/lengths].front() will be used as additional last point, generating a closed path + * \param tmesh input triangle mesh supporting the vertices of the output polygon + * \param solver container for the precomputed information. If not initialized, it will be initialized internally. + * \return a face location for each vertex of each polygon + * \todo add named parameters + * \todo polygon orientation is not handled in the function and should be done outside of the function for now + * \todo for better rendering we can group polygons to have one center for the same group of polygon (useful for letters that are not simply connected) + * \todo what if boundary is reached + */ +//TODO: make sure it is consistent with the rest to not duplicate the last point if closed +template +std::vector> +trace_polyline_of_bezier_curves(const CGAL::Polygon_mesh_processing::Face_location ¢er, + const std::vector& directions, + const std::vector& lengths, + bool is_closed, // use [directions/lengths].front as last control point? + const int num_subdiv, + const TriangleMesh &tmesh +#ifndef CGAL_BSURF_USE_DIJKSTRA_SP + , const Dual_geodesic_solver& solver = {} +#endif +) +{ + namespace PMP = CGAL::Polygon_mesh_processing; + + using FT = typename K::FT; + + std::size_t n = (directions.size() - (is_closed?0:1))/3; + CGAL_assertion( n * 3 + (is_closed?0:1) == directions.size() ); + + std::vector> result; + + // n is the number of quadruple of control points + // After num_subdiv steps, we have 2^num_subdiv * n quadruples of control points + + + // even if closed we will duplicate the last point + // (this is a lower bound without taking into account shortest path between points) + result.reserve( (1<* solver_ptr=&solver; + Dual_geodesic_solver local_solver; + if (solver.graph.empty()) + { + solver_ptr = &local_solver; + init_geodesic_dual_solver(local_solver, tmesh); + } +#endif + +#ifdef CGAL_DEBUG_BSURF + std::ofstream debug_cp("/tmp/control_points.xyz"); + std::ofstream debug_ep("/tmp/end_points.xyz"); + debug_cp << std::setprecision(17); + debug_ep << std::setprecision(17); +#endif + + PMP::Face_location prev_loc = straightest_geodesic(center,directions[0],lengths[0],tmesh).back(), + first_loc = prev_loc; + + for (std::size_t i=0; i control_loc; + control_loc[0]=prev_loc; + for (int k=1;k<4; ++k) + { + if (k!=3 || !is_closed || 3*i+k!=directions.size()) + control_loc[k] = straightest_geodesic(center,directions[3*i+k],lengths[3*i+k],tmesh).back(); + else + control_loc[k] = first_loc; + } + prev_loc=control_loc[3]; + + #ifdef CGAL_DEBUG_BSURF + debug_ep << PMP::construct_point(control_loc[0], tmesh) << "\n"; + debug_ep << PMP::construct_point(control_loc[3], tmesh) << "\n"; + debug_cp << PMP::construct_point(control_loc[1], tmesh) << "\n"; + debug_cp << PMP::construct_point(control_loc[2], tmesh) << "\n"; + #endif + + std::vector> bezier = + recursive_de_Casteljau(tmesh, control_loc, num_subdiv +#ifndef CGAL_BSURF_USE_DIJKSTRA_SP + , *solver_ptr +#endif + ); + + if (i==0) + result.push_back(bezier[0]); + for(std::size_t b=0; b& loc = bezier[b]; + const PMP::Face_location& loc1 = bezier[b+1]; + + // connect the two face locations with shortest path is they are in different faces + if (loc.first!=loc1.first) + { + std::vector> edge_locations; + locally_shortest_path(loc, loc1, tmesh, edge_locations, solver); + for (const PMP::Edge_location& e : edge_locations) + result.push_back(PMP::to_face_location(e, tmesh)); + } + result.push_back(loc1); + } + } + + return result; +} + + + +/*! + * \ingroup VGSFunctions + * refines `tmesh` so that each path in `paths` corresponds to a set edges of `tmesh` after the call. + * Note that each path must be such that for two consecutive face locations, there exists a face in `tmesh` containing the two corresponding points. + * \tparam TriangleMesh a model of `MutableFaceGraph` + * \tparam K a model of `Kernel` with `K::FT` being a floating point number type (float or double) + * \tparam VNM a model of `ReadWritePropertyMap` with `boost::graph_traits::%vertex_descriptor` as key type and `Kernel::Vector_3` as value type. + * \tparam FNM a model of `ReadWritePropertyMap` with `boost::graph_traits::%face_descriptor` as key type and `Kernel::Vector_3` as value type. + * \tparam OutputIterator an output iterator accepting `boost::graph_traits::%halfedge_descriptor` to be put + * \param tmesh the triangle mesh to be refined + * \param paths a path described as a range of edge locations, with the property that + * for two consecutive edge locations, there exists a face in `tmesh` containing the two corresponding points. + * \param vnm property map associating a normal to each vertex of `tmesh` that is updated by this function + * \param fnm property map associating a normal to each face of `tmesh` that is updated by this function + * \param out output iterator where created halfedges are put + * \todo add named parameters + * \todo vnm and fnm are optional + * \todo intersection between path or self-intersections are not handle. Should it be? If not what do we do? + * \todo out should contain edges, and also existing edges already part of a path + * \todo shall we also have the edges in the order of the input rather than all at once + */ +template +// std::vector::vertex_descriptor> +OutputIterator +refine_mesh_along_paths(const std::vector>>& paths, + TriangleMesh& tmesh, + VNM vnm, + FNM fnm, + OutputIterator out) +{ + namespace PMP = CGAL::Polygon_mesh_processing; + + // TODO: nothing is done here to identify identical points + + using vertex_descriptor = typename boost::graph_traits::vertex_descriptor; + using edge_descriptor = typename boost::graph_traits::edge_descriptor; + using halfedge_descriptor = typename boost::graph_traits::halfedge_descriptor; + using face_descriptor = typename boost::graph_traits::face_descriptor; + using Face_loc = PMP::Face_location; + using dscrptr_vrnt = PMP::descriptor_variant; + using Point_3 = typename K::Point_3; + using EK = CGAL::Exact_predicates_exact_constructions_kernel; + + // 3 types of points: on edges, on faces, and existing vertices + typename boost::property_map>::type + fid_map = get(CGAL::dynamic_face_property_t(), tmesh, std::size_t(-1)); + typename boost::property_map>::type + eid_map = get(CGAL::dynamic_edge_property_t(), tmesh, std::size_t(-1)); + typename boost::property_map>::type + vid_map = get(CGAL::dynamic_vertex_property_t(), tmesh, std::size_t(-1)); + + std::vector< std::vector> > edges_per_face; + std::vector< std::vector > points_per_face; + std::vector< std::vector > points_per_edge; + std::vector< edge_descriptor > edges_to_refine; + std::vector< face_descriptor > faces_to_refine; + std::vector< std::pair > input_vertices; + + std::vector polyline_locations; + + auto register_point = [&](const dscrptr_vrnt& var, std::size_t id) + { + switch(var.index()) + { + case 1: + { + halfedge_descriptor hd=std::get(var); + edge_descriptor ed = edge(hd, tmesh); + if (get(eid_map, ed)==std::size_t(-1)) + { + put(eid_map, ed, edges_to_refine.size()); + points_per_edge.emplace_back(); + edges_to_refine.push_back(ed); + } + points_per_edge[get(eid_map,ed)].push_back(id); + return; + } + case 2: + { + face_descriptor fd=std::get(var); + if (get(fid_map, fd)==std::size_t(-1)) + { + put(fid_map, fd, faces_to_refine.size()); + points_per_face.emplace_back(); + edges_per_face.emplace_back(); + faces_to_refine.push_back(fd); + } + points_per_face[get(fid_map,fd)].push_back(id); + return; + } + default: + input_vertices.emplace_back(id, std::get(var)); + } + }; + + auto register_segment = [&](const dscrptr_vrnt& v1, const dscrptr_vrnt& v2, + const Face_loc& loc1, const Face_loc& loc2, + std::size_t i1, std::size_t i2) + { + if (v1.index()==0 && v2.index()==0) return; + if (v1.index()==1 && v2.index()==1 && v1==v2) return; + if (v1.index()==0 && v2.index()==1) + { + vertex_descriptor vd = std::get(v1); + halfedge_descriptor hd = std::get(v2); + if (vd==source(hd, tmesh) || vd==target(hd, tmesh)) + return; + } + if (v2.index()==0 && v1.index()==1) + { + vertex_descriptor vd = std::get(v2); + halfedge_descriptor hd = std::get(v1); + if (vd==source(hd, tmesh) || vd==target(hd, tmesh)) + return; + } + Face_loc copy1=loc1, copy2=loc2; + CGAL_assertion_code(bool OK=) + PMP::locate_in_common_face(copy1, copy2, tmesh); + CGAL_assertion(OK); + if (get(fid_map,copy1.first)==std::size_t(-1)) + { + put(fid_map, copy1.first, faces_to_refine.size()); + points_per_face.emplace_back(); + edges_per_face.emplace_back(); + faces_to_refine.push_back(copy1.first); + } + edges_per_face[get(fid_map, copy1.first)].emplace_back(i1,i2); + }; + + for (const std::vector& path : paths) + { + bool closed = PMP::are_locations_identical(path.front(), path.back(), tmesh); + std::size_t nbp = path.size(); + + if (closed) nbp-=1; + + std::size_t prev_id=polyline_locations.size(), first_id=prev_id; + dscrptr_vrnt prev_var = PMP::get_descriptor_from_location(path.front(), tmesh), + first_var = prev_var; + + polyline_locations.push_back( path.front() ); + register_point(prev_var, prev_id); + + for (std::size_t i=1; i face_normals(faces_to_refine.size()); + std::vector> face_input_vertices(faces_to_refine.size()); + for (std::size_t fid=0; fid exact_points(polyline_locations.size()); + std::vector polyline_vertices( + polyline_locations.size(), boost::graph_traits::null_vertex()); + + // collect split point on edges + std::vector hedges_to_refine(edges_to_refine.size()); + std::vector< std::vector > > points_on_hedges(edges_to_refine.size()); + for (std::size_t eid=0; eid > coordinates; + coordinates.reserve(points_per_edge[eid].size()); + for (std::size_t pid : points_per_edge[eid]) + { + halfedge_descriptor hd = halfedge(polyline_locations[pid].first, tmesh); + int src_id=0; + if (edge(hd, tmesh)!=e) + { + hd=next(hd, tmesh); ++src_id; + if (edge(hd, tmesh)!=e) + { + hd=next(hd, tmesh); + ++src_id; + } + } + CGAL_assertion(edge(hd, tmesh)==e); + coordinates.emplace_back(polyline_locations[pid].second[src_id], polyline_locations[pid].second[(src_id+1)%3]); + if (hd!=href) + { + CGAL_assertion( opposite(hd, tmesh)==href ); + std::swap( coordinates.back().first, coordinates.back().second ); + } + } + + // now sort coordinates + std::vector ids(coordinates.size()); + std::iota(ids.begin(), ids.end(), 0); + + std::sort(ids.begin(), ids.end(), [&coordinates](std::size_t i, std::size_t j) + { return coordinates[i].first < coordinates[j].first; }); + + points_on_hedges[eid].reserve(ids.size()); + for (std::size_t id : ids) + { + points_on_hedges[eid].emplace_back( + PMP::construct_point(polyline_locations[points_per_edge[eid][id]], tmesh), + points_per_edge[eid][id]); + + exact_points[points_per_edge[eid][id]] = + PMP::construct_point(polyline_locations[points_per_edge[eid][id]], tmesh, + parameters::vertex_point_map(make_cartesian_converter_property_map(tmesh.points()))); + + } + } + + CGAL::Cartesian_converter to_exact; + + // add new vertices per face + for (std::size_t fid=0; fid::null_vertex() ); + polyline_vertices[vid] = vd; + exact_points[vid] = PMP::construct_point(polyline_locations[vid], tmesh, + parameters::vertex_point_map(make_cartesian_converter_property_map(tmesh.points()))); + } + } + + // set existing vertex + for (const std::pair& id_and_vd : input_vertices) + { + polyline_vertices[id_and_vd.first]=id_and_vd.second; + put(vid_map, id_and_vd.second, id_and_vd.first); + exact_points[id_and_vd.first]=to_exact(tmesh.point(id_and_vd.second)); + } + + //Now split edges + for (std::size_t eid=0; eid& pt_and_id : points_on_hedges[eid]) + { + CGAL_assertion( polyline_vertices[pt_and_id.second]==boost::graph_traits::null_vertex() ); + href = ::CGAL::Euler::split_edge(href, tmesh); + polyline_vertices[pt_and_id.second]=target(href, tmesh); + // TODO: use VPM + tmesh.point(target(href,tmesh)) = std::move(pt_and_id.first); + put(vid_map, target(href,tmesh), pt_and_id.second); + put(vnm, target(href,tmesh), edge_normal); + } + } + + // triangulate faces + using CDT_traits = Projection_traits_3; + using Vb = Triangulation_vertex_base_with_info_2; + using Fb = Constrained_triangulation_face_base_2; + using TDS_2 = Triangulation_data_structure_2; + using CDT = Constrained_Delaunay_triangulation_2; + using CDT_Vertex_handle = typename CDT::Vertex_handle; + + std::size_t nb_verts = polyline_vertices.size(); + polyline_vertices.resize(polyline_vertices.size()+3); // for input triangle vertices + exact_points.resize(exact_points.size()+3); + + std::vector id_to_cdt_vhandles(nb_verts); + + for (std::size_t fid=0; fid tri_verts = + make_array(face_input_vertices[fid][0], face_input_vertices[fid][1], face_input_vertices[fid][2]); + + std::array vhandles; + vhandles[0]=cdt.insert_outside_affine_hull(to_exact(tmesh.point(tri_verts[0]))); + vhandles[1]=cdt.insert_outside_affine_hull(to_exact(tmesh.point(tri_verts[1]))); + vhandles[2] = cdt.tds().insert_dim_up(cdt.infinite_vertex(), false); + vhandles[2]->set_point(to_exact(tmesh.point(tri_verts[2]))); + + // assign ids to input triangle vertices + std::size_t offset=nb_verts; + for (int i=0;i<3;++i) + { + vhandles[i]->info()=tri_verts[i]; + if (get(vid_map, tri_verts[i])>=nb_verts) // and not simply -1 as previous triangulation could have set another dummy value + { + polyline_vertices[offset]=tri_verts[i]; + put(vid_map, tri_verts[i], offset); + exact_points[offset]=vhandles[i]->point(); + ++offset; + } + else + { + id_to_cdt_vhandles[get(vid_map, tri_verts[i])]=vhandles[i]; + } + } + + + face_descriptor fd = faces_to_refine[fid]; + std::array face_sides; + for (int i=0; i<3; ++i) + { + for (halfedge_descriptor h : halfedges_around_target(face_input_vertices[fid][i], tmesh)) + { + if (face(h, tmesh) == fd) + { + face_sides[i]=h; + break; + } + } + } + + // insert points on edges + CDT_Vertex_handle vinf = cdt.infinite_vertex(); + for (int k=0; k<3; ++k) + { + if (next(face_sides[k], tmesh) != face_sides[(k+1)%3]) + { + CDT_Vertex_handle src = vhandles[k], tgt = vhandles[(k+1)%3]; + halfedge_descriptor hcurr=next(face_sides[k], tmesh); + CGAL_assertion(src->info()==source(hcurr, tmesh)); + CGAL_assertion(tgt->info()==target(face_sides[(k+1)%3], tmesh)); + + CDT_Vertex_handle prev = src; + while(hcurr!=face_sides[(k+1)%3]) + { + typename CDT::Face_handle fh; + CGAL_assertion_code(bool ok =) + cdt.is_face(prev, tgt, vinf, fh); + CGAL_assertion(ok); + CGAL_assertion( get(vid_map, target(hcurr,tmesh)) != std::size_t(-1) ); + prev=cdt.insert_in_edge(exact_points[get(vid_map, target(hcurr,tmesh))], fh, fh->index(vinf)); + id_to_cdt_vhandles[get(vid_map, target(hcurr,tmesh))]=prev; + prev->info() = target(hcurr,tmesh); + cdt.restore_Delaunay(prev); // TODO maybe not each time but one global? + CGAL_assertion(cdt.is_valid()); + hcurr=next(hcurr, tmesh); + } + } + } + + // insert points in the face (TODO: we probably don't need spatial sorting) + for (std::size_t vid : points_per_face[fid]) + { + id_to_cdt_vhandles[vid] = cdt.insert(exact_points[vid]); + id_to_cdt_vhandles[vid]->info()=polyline_vertices[vid]; + } + + + + // insert constrained edges + for (const std::pair& pi : edges_per_face[fid]) + { + cdt.insert_constraint(id_to_cdt_vhandles[pi.first],id_to_cdt_vhandles[pi.second]); + } + + // register cdt edge -> halfedge + halfedge_descriptor hd = halfedge(fd, tmesh); + std::map, halfedge_descriptor> edge_to_hedge; + // triangle boundary first + for (halfedge_descriptor h : halfedges_around_face(hd, tmesh)) + { + edge_to_hedge[std::make_pair(get(vid_map,source(h,tmesh)), get(vid_map,target(h,tmesh)))]=h; + } + //grab edges that are not on the convex hull (these have already been created) + for (typename CDT::Finite_edges_iterator it=cdt.finite_edges_begin(); + it!=cdt.finite_edges_end(); ++it) + { + typename CDT::Vertex_handle cdt_v0=it->first->vertex( cdt.ccw(it->second) ); + typename CDT::Vertex_handle cdt_v1=it->first->vertex( cdt.cw(it->second) ); + + // consider edges not on the convex hull (not on the boundary of the face) + // and create the corresponding halfedges + if ( !cdt.is_infinite(it->first->vertex(it->second)) && + !cdt.is_infinite(cdt.mirror_vertex(it->first,it->second)) ) + { + edge_descriptor e=add_edge(tmesh); + halfedge_descriptor h=halfedge(e,tmesh), h_opp=opposite(h,tmesh); + std::size_t i0=get(vid_map,cdt_v0->info()), i1=get(vid_map,cdt_v1->info()); + vertex_descriptor v0=polyline_vertices[i0], v1=polyline_vertices[i1]; + + set_target(h,v0,tmesh); + set_target(h_opp,v1,tmesh); + set_halfedge(v0,h,tmesh); + set_halfedge(v1,h_opp,tmesh); + + edge_to_hedge[std::make_pair(i0,i1)]=h_opp; + edge_to_hedge[std::make_pair(i1,i0)]=h; + + if (it->first->is_constrained(it->second)) + *out++=h; + } + } + + //grab triangles. + face_descriptor current_face = fd; + typename K::Vector_3 face_normal = get(fnm, fd); + for (typename CDT::Finite_faces_iterator it=cdt.finite_faces_begin(), + it_end=cdt.finite_faces_end();;) + { + typename CDT::Vertex_handle cdt_v0=it->vertex(0); + typename CDT::Vertex_handle cdt_v1=it->vertex(1); + typename CDT::Vertex_handle cdt_v2=it->vertex(2); + + std::size_t i0=get(vid_map,cdt_v0->info()), i1=get(vid_map,cdt_v1->info()), i2=get(vid_map,cdt_v2->info()); + + CGAL_assertion(edge_to_hedge.count(std::make_pair(i0,i1))!= 0); + CGAL_assertion(edge_to_hedge.count(std::make_pair(i1,i2))!= 0); + CGAL_assertion(edge_to_hedge.count(std::make_pair(i2,i0))!= 0); + + halfedge_descriptor h01=edge_to_hedge[std::make_pair(i0,i1)]; + halfedge_descriptor h12=edge_to_hedge[std::make_pair(i1,i2)]; + halfedge_descriptor h20=edge_to_hedge[std::make_pair(i2,i0)]; + + set_next(h01,h12,tmesh); + set_next(h12,h20,tmesh); + set_next(h20,h01,tmesh); + + //update face halfedge + set_halfedge(current_face,h01,tmesh); + + //update face of halfedges + set_face(h01,current_face,tmesh); + set_face(h12,current_face,tmesh); + set_face(h20,current_face,tmesh); + + if ( ++it!=it_end ) + { + current_face=add_face(tmesh); + put(fnm, current_face, face_normal); + } + else + break; + } + } + + return out; +} + +/*! + * \ingroup VGSMiscellaneous + * converts a path on a triangle mesh to the corresponding polyline of points. + * If `path` contains identical consecutive vertices, only one point will be put in `poly_out` for this vertex. + * \tparam FT floating point number type (float or double) + * \tparam TriangleMesh a model of `FaceGraph` + * \tparam OutputIterator an output iterator accepting points from `tmesh` + * \param path a path described as a range of face locations, with the property that + for two consecutive face locations, there exists a face in `tmesh` containing the two corresponding points. + * \param tmesh the triangle mesh supporing the path + * \param poly_out output iterator where points of the polyline are put. + * \todo add named parameters + * \todo generic range + */ +template +OutputIterator +convert_path_to_polyline(const std::vector>& path, + const TriangleMesh& tmesh, + OutputIterator poly_out) +{ + namespace PMP = CGAL::Polygon_mesh_processing; + + using vertex_descriptor = typename boost::graph_traits::vertex_descriptor; + + vertex_descriptor vd = boost::graph_traits::null_vertex(); + for (const PMP::Face_location& loc : path) + { + std::optional ov = PMP::vertex_descriptor_from_location(loc, tmesh); + if (ov.has_value()) + { + if (vd==*ov) continue; + vd=*ov; + } + *poly_out++=PMP::construct_point(loc, tmesh); + } + return poly_out; +} + +/*! + * \ingroup VGSMiscellaneous + * converts a path on a triangle mesh to the corresponding polyline of points. + * If `path` contains identical consecutive vertices, only one point will be put in `poly_out` for this vertex. + * \tparam FT floating point number type (float or double) + * \tparam TriangleMesh a model of `FaceGraph` + * \tparam OutputIterator an output iterator accepting points from `tmesh` + * \param src source of the path + * \param tgt target of the path + * \param path a path described as a range of edge locations, with the property that + for two consecutive edge locations, there exists an edge in `tmesh` containing the two corresponding points. + * \param tmesh the triangle mesh supporing the path + * \param poly_out output iterator where points of the polyline are put. + * \todo add named parameters + * \todo generic range + */ +template +OutputIterator +convert_path_to_polyline(const CGAL::Polygon_mesh_processing::Face_location& src, + const std::vector>& path, + const CGAL::Polygon_mesh_processing::Face_location& tgt, + const TriangleMesh& tmesh, + OutputIterator poly_out) +{ + namespace PMP = CGAL::Polygon_mesh_processing; + + using vertex_descriptor = typename boost::graph_traits::vertex_descriptor; + + *poly_out++=PMP::construct_point(src, tmesh); + vertex_descriptor vd = boost::graph_traits::null_vertex(); + for (const PMP::Edge_location& loc : path) + { + std::optional ov = PMP::vertex_descriptor_from_location(loc, tmesh); + if (ov.has_value()) + { + if (vd==*ov) continue; + vd=*ov; + } + *poly_out++=PMP::construct_point(loc, tmesh); + } + *poly_out++=PMP::construct_point(tgt, tmesh); + return poly_out; +} + +// template +// std::vector +// tangent_path_direction(const std::vector>& path, +// const CGAL::Polygon_mesh_processing::Face_location& src, +// const CGAL::Polygon_mesh_processing::Face_location& tgt, +// const TriangleMesh &tmesh,const bool initial=true) +// { +// auto find = [](const std::array &vec, int x) { +// for (int i = 0; i < vec.size(); i++) +// if (vec[i] == x) +// return i; +// return -1; +// }; +// typename K::Vector_2 direction; +// using VPM = typename boost::property_map::const_type; +// VPM vpm = get(CGAL::vertex_point, tmesh); + + +// if(initial) +// { + +// halfedge_descriptor h_ref=halfedge(src.first,mesh); +// std::array flat_tid=init_flat_triangle(h_ref,vpm,tmesh); +// if(path.size()==0) +// { +// assert(src.first==tgt.first);//TODO:src and tgt may have different faces because we do not update them when cleaning the strip +// typename K::Vector_2 flat_src=src.second[0]*flat_tid[0]+src.second[1]*flat_tid[1]+src.second[2]*flat_tid[2]; +// typename K::Vector_2 flat_tgt=tgt.second[0]*flat_tid[0]+tgt.second[1]*flat_tid[1]+tgt.second[2]*flat_tid[2]; +// direction=normalize(flat_tgt-flat_src); +// }else{ +// halfedge_descriptor h_edge=halfedge(path[0].first,tmesh); +// int k=0; + +// if(h_edge==prev(h_ref,mesh)) +// k=2; +// else if(h_edge==next(h_ref,mesh)) +// k=1; +// else +// assert(h_edge==h_ref); +// } + + + + + + + + + + +// } + + + + + + +// } + +template +std::vector +trace_agap_polygon(const CGAL::Polygon_mesh_processing::Face_location ¢er, + const std::vector& directions, + const std::vector& lengths, + const TriangleMesh &tmesh, + const Dual_geodesic_solver& solver = {}) +{ + namespace PMP = CGAL::Polygon_mesh_processing; + + size_t n=directions.size(); + std::vector result; + std::vector> vertices(n); + for(std::size_t i=0;i(center,directions[i],lengths[i],tmesh).back(); + + std::vector> edge_locations; + + const Dual_geodesic_solver* solver_ptr=&solver; + Dual_geodesic_solver local_solver; + if (solver.graph.empty()) + { + solver_ptr = &local_solver; + init_geodesic_dual_solver(local_solver, tmesh); + } + + for(std::size_t i=0;i(vertices[i],vertices[(i+1)%n],tmesh, edge_locations, *solver_ptr); + result.push_back(PMP::construct_point(vertices[i],tmesh)); + for(auto& el : edge_locations) + { + result.push_back(PMP::construct_point(el, tmesh)); + } + } + result.push_back(PMP::construct_point(vertices[0],tmesh)); + + return result; + +} + +} // namespace Vector_graphics_on_surfaces +} // namespace CGAL + +#endif diff --git a/Vector_graphics_on_surfaces/package_info/Vector_graphics_on_surfaces/copyright b/Vector_graphics_on_surfaces/package_info/Vector_graphics_on_surfaces/copyright new file mode 100644 index 00000000000..2427333ef57 --- /dev/null +++ b/Vector_graphics_on_surfaces/package_info/Vector_graphics_on_surfaces/copyright @@ -0,0 +1 @@ +GeometryFactory diff --git a/Vector_graphics_on_surfaces/package_info/Vector_graphics_on_surfaces/dependencies b/Vector_graphics_on_surfaces/package_info/Vector_graphics_on_surfaces/dependencies new file mode 100644 index 00000000000..2bb950b1db4 --- /dev/null +++ b/Vector_graphics_on_surfaces/package_info/Vector_graphics_on_surfaces/dependencies @@ -0,0 +1,30 @@ +AABB_tree +Algebraic_foundations +Arithmetic_kernel +BGL +CGAL_Core +Cartesian_kernel +Circulator +Distance_2 +Distance_3 +Filtered_kernel +Hash_map +Homogeneous_kernel +Installation +Intersections_2 +Intersections_3 +Interval_support +Kernel_23 +Kernel_d +Modular_arithmetic +Number_types +Polygon_mesh_processing +Profiling_tools +Property_map +Random_numbers +STL_Extension +Spatial_searching +Spatial_sorting +Stream_support +TDS_2 +Triangulation_2 diff --git a/Vector_graphics_on_surfaces/package_info/Vector_graphics_on_surfaces/description.txt b/Vector_graphics_on_surfaces/package_info/Vector_graphics_on_surfaces/description.txt new file mode 100644 index 00000000000..d878bed185e --- /dev/null +++ b/Vector_graphics_on_surfaces/package_info/Vector_graphics_on_surfaces/description.txt @@ -0,0 +1,2 @@ +Package Vector Graphics on Triangulated Surface Meshes: +This package provides functions to draw basic geometric primitives such as geodesics, Bézier and B-Spline curves... \ No newline at end of file diff --git a/Vector_graphics_on_surfaces/package_info/Vector_graphics_on_surfaces/license.txt b/Vector_graphics_on_surfaces/package_info/Vector_graphics_on_surfaces/license.txt new file mode 100644 index 00000000000..8bb8efcb72b --- /dev/null +++ b/Vector_graphics_on_surfaces/package_info/Vector_graphics_on_surfaces/license.txt @@ -0,0 +1 @@ +GPL (v3 or later) diff --git a/Vector_graphics_on_surfaces/package_info/Vector_graphics_on_surfaces/maintainer b/Vector_graphics_on_surfaces/package_info/Vector_graphics_on_surfaces/maintainer new file mode 100644 index 00000000000..b02f658d427 --- /dev/null +++ b/Vector_graphics_on_surfaces/package_info/Vector_graphics_on_surfaces/maintainer @@ -0,0 +1 @@ +Claudio Mancinelli and Sébastien Loriot \ No newline at end of file diff --git a/Vector_graphics_on_surfaces/test/Vector_graphics_on_surfaces/CMakeLists.txt b/Vector_graphics_on_surfaces/test/Vector_graphics_on_surfaces/CMakeLists.txt new file mode 100644 index 00000000000..27a53ec8eca --- /dev/null +++ b/Vector_graphics_on_surfaces/test/Vector_graphics_on_surfaces/CMakeLists.txt @@ -0,0 +1,17 @@ +# Created by the script cgal_create_cmake_script +# This is the CMake script for compiling a CGAL application. + +cmake_minimum_required(VERSION 3.12...3.29) +project(Vector_graphics_on_surfaces_Tests) + +find_package(CGAL REQUIRED) +find_package(Eigen3 REQUIRED) + +include(CGAL_Eigen3_support) + +create_single_source_cgal_program( "test_straightest_geodesic.cpp" ) +target_link_libraries(test_straightest_geodesic PRIVATE CGAL::Eigen3_support) +create_single_source_cgal_program("test_Bsurf.cpp") +target_link_libraries(test_Bsurf PUBLIC CGAL::Eigen3_support) +create_single_source_cgal_program("test_Bsurf_grid.cpp") +target_link_libraries(test_Bsurf_grid PUBLIC CGAL::Eigen3_support) diff --git a/Vector_graphics_on_surfaces/test/Vector_graphics_on_surfaces/test_Bsurf.cpp b/Vector_graphics_on_surfaces/test/Vector_graphics_on_surfaces/test_Bsurf.cpp new file mode 100644 index 00000000000..e028de3bd23 --- /dev/null +++ b/Vector_graphics_on_surfaces/test/Vector_graphics_on_surfaces/test_Bsurf.cpp @@ -0,0 +1,287 @@ +#include +#include +#include +#include + +namespace VGoS = CGAL::Vector_graphics_on_surfaces; +namespace PMP = CGAL::Polygon_mesh_processing; + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Surface_mesh Mesh; +typedef PMP::Face_location Face_location; +typedef PMP::Edge_location Edge_location; + + +int main(int argc, char** argv) +{ + std::string filename = (argc > 1) ? std::string(argv[1]) + : CGAL::data_file_path("meshes/elephant.off"); + + Mesh mesh; + if(!CGAL::IO::read_polygon_mesh(filename, mesh) || !CGAL::is_triangle_mesh(mesh)) + { + std::cerr << "Invalid input." << std::endl; + return 1; + } + VGoS::Dual_geodesic_solver solver; + VGoS::init_geodesic_dual_solver(solver, mesh); + + + std::size_t nb_hedges = halfedges(mesh).size(); + + // take two random faces and pick the centroid + CGAL::Random rnd = CGAL::get_default_random(); + //CGAL::Random rnd(1706604646); + // CGAL::Random rnd(1695724381); + // CGAL::Random rnd(1695813638); + + std::cout << "seed " << rnd.get_seed() << std::endl; + Mesh::Halfedge_index h = *std::next(halfedges(mesh).begin(), rnd.get_int(0, nb_hedges)); + // Mesh::Halfedge_index h = *std::next(halfedges(mesh).begin(), 2*6178); // <---- interesting two different locally shortest paths + + std::cout << "h = " << h << "\n"; + + std::ofstream out("locally_shortest_path.polylines.txt"); + +// test src/tgt being 2 opposite vertices of an edge +#if 1 + for (int loop=0; loop<3; ++loop) + { + Mesh::Vertex_index v1 = target(next(h, mesh), mesh); + Mesh::Vertex_index v2 = target(next(opposite(h, mesh), mesh), mesh); + + bool first_run = true; + std::size_t expected_size=0; + for(Mesh::Halfedge_index h1 : CGAL::halfedges_around_target(v1, mesh)) + { + Mesh::Face_index f1 = face(h1, mesh); + Mesh::Halfedge_index hf1 = prev(halfedge(f1, mesh), mesh); + int i1=0; + while (target(hf1, mesh)!= v1) + { + hf1=next(hf1,mesh); + ++i1; + } + Face_location src(f1, CGAL::make_array(0.,0.,0.)); + src.second[i1]=1; + for(Mesh::Halfedge_index h2 : CGAL::halfedges_around_target(v2, mesh)) + { + Mesh::Face_index f2 = face(h2, mesh); + Mesh::Halfedge_index hf2 = prev(halfedge(f2, mesh), mesh); + int i2=0; + while (target(hf2, mesh)!= v2) + { + hf2=next(hf2,mesh); + ++i2; + } + Face_location tgt(f2, CGAL::make_array(0.,0.,0.)); + tgt.second[i2]=1; + + std::cout << "Running " << f1 << " " << v1 << " | " << f2 << " " << v2 << "\n"; + + std::vector edge_locations; + auto src_bk=src, tgt_bk=tgt; + VGoS::locally_shortest_path(src, tgt, mesh, edge_locations, solver); + assert(get(PMP::get_descriptor_from_location(src,mesh))==v1); + assert(get(PMP::get_descriptor_from_location(tgt,mesh))==v2); + + out << edge_locations.size()+2; + out << " " << PMP::construct_point(src, mesh); + for (auto el : edge_locations) + out << " " << PMP::construct_point(el, mesh); + out << " " << PMP::construct_point(tgt, mesh) << "\n"; + out << std::flush; + + if (first_run) + { + first_run=false; + expected_size=edge_locations.size(); + } + if(edge_locations.size() != expected_size) + { + std::cout << edge_locations.size() << " vs " << expected_size << "\n"; + } + CGAL_warning(edge_locations.size() == expected_size); + + src=src_bk; + tgt=tgt_bk; + edge_locations.clear(); + VGoS::locally_shortest_path(tgt, src, mesh, edge_locations, solver); + assert(get(PMP::get_descriptor_from_location(src,mesh))==v1); + assert(get(PMP::get_descriptor_from_location(tgt,mesh))==v2); + + out << edge_locations.size()+2; + out << " " << PMP::construct_point(tgt, mesh); + for (auto el : edge_locations) + out << " " << PMP::construct_point(el, mesh); + out << " " << PMP::construct_point(src, mesh) << "\n"; + out << std::flush; + CGAL_warning(edge_locations.size() == expected_size); + } + } + h = next(h, mesh); + } +#endif + +#if 1 + // test src is a vertex and tgt is on an edge + for (int i=0; i<3; ++i) + { + Mesh::Vertex_index v1 = target(next(h, mesh), mesh); + Mesh::Halfedge_index h2 = opposite(next(opposite(h, mesh), mesh), mesh); + bool first_run=true; + std::size_t expected_size=0; + + for(Mesh::Halfedge_index h1 : CGAL::halfedges_around_target(v1, mesh)) + { + Mesh::Face_index f1 = face(h1, mesh); + Mesh::Halfedge_index hf1 = prev(halfedge(f1, mesh), mesh); + int i1=0; + while (target(hf1, mesh)!= v1) + { + hf1=next(hf1,mesh); + ++i1; + } + Face_location src(f1, CGAL::make_array(0.,0.,0.)); + src.second[i1]=1; + + // define tgt + Mesh::Face_index f2 = face(h2, mesh); + Mesh::Halfedge_index hf2 = prev(halfedge(f2, mesh), mesh); + int k=0; + while(hf2!=h2) + { + hf2=next(hf2, mesh); + ++k; + } + + Face_location tgt(f2, CGAL::make_array(0.5,0.5,0.5)); + tgt.second[(k+1)%3]=0; + tgt.second[(k+2)%3]=0.25; + tgt.second[(k)%3]=0.75; + + assert(get(PMP::get_descriptor_from_location(src,mesh))==v1); + assert(edge(get(PMP::get_descriptor_from_location(tgt,mesh)), mesh)==edge(h2,mesh)); + + std::cout << "Running " << f1 << " " << v1 << " | " << f2 << " " << edge(h2,mesh) << "\n"; + // std::cout << " " << PMP::construct_point(src, mesh) << " | " << PMP::construct_point(tgt, mesh) << "\n"; + std::vector edge_locations; + auto src_bk=src, tgt_bk=tgt; + VGoS::locally_shortest_path(src, tgt, mesh, edge_locations, solver); + + out << edge_locations.size()+2; + out << " " << PMP::construct_point(src, mesh); + for (auto el : edge_locations) + out << " " << PMP::construct_point(el, mesh); + out << " " << PMP::construct_point(tgt, mesh) << "\n"; + out << std::flush; + + if (first_run) + { + first_run=false; + expected_size=edge_locations.size(); + } + CGAL_warning(edge_locations.size()==expected_size); + + assert(get(PMP::get_descriptor_from_location(src,mesh))==v1); + assert(edge(get(PMP::get_descriptor_from_location(tgt,mesh)), mesh)==edge(h2,mesh)); + + src=src_bk; + tgt=tgt_bk; + edge_locations.clear(); + VGoS::locally_shortest_path(tgt, src, mesh, edge_locations, solver); + + out << edge_locations.size()+2; + out << " " << PMP::construct_point(tgt, mesh); + for (auto el : edge_locations) + out << " " << PMP::construct_point(el, mesh); + out << " " << PMP::construct_point(src, mesh) << "\n"; + out << std::flush; + CGAL_warning(edge_locations.size()==expected_size); + + assert(get(PMP::get_descriptor_from_location(src,mesh))==v1); + assert(edge(get(PMP::get_descriptor_from_location(tgt,mesh)), mesh)==edge(h2,mesh)); + } + h = next(h, mesh); + } +#endif +#if 1 + // test src is on an edge and tgt is on an edge + for (int i=0; i<3; ++i) + { + Mesh::Halfedge_index h1 = opposite(next(h, mesh), mesh); + Mesh::Halfedge_index h2 = opposite(next(opposite(h, mesh), mesh), mesh); + + // define src + Mesh::Face_index f1 = face(h1, mesh); + Mesh::Halfedge_index hf1 = prev(halfedge(f1, mesh), mesh); + int k1=0; + while(hf1!=h1) + { + hf1=next(hf1, mesh); + ++k1; + } + + Face_location src(f1, CGAL::make_array(0.5,0.5,0.5)); + src.second[(k1+1)%3]=0; + src.second[(k1+2)%3]=0.25; + src.second[(k1)%3]=0.75; + // define tgt + Mesh::Face_index f2 = face(h2, mesh); + Mesh::Halfedge_index hf2 = prev(halfedge(f2, mesh), mesh); + int k2=0; + while(hf2!=h2) + { + hf2=next(hf2, mesh); + ++k2; + } + + Face_location tgt(f2, CGAL::make_array(0.5,0.5,0.5)); + tgt.second[(k2+1)%3]=0; + tgt.second[(k2+2)%3]=0.25; + tgt.second[(k2)%3]=0.75; + + assert(edge(get(PMP::get_descriptor_from_location(src,mesh)), mesh)==edge(h1,mesh)); + assert(edge(get(PMP::get_descriptor_from_location(tgt,mesh)), mesh)==edge(h2,mesh)); + + std::cout << "Running " << f1 << " " << edge(h1,mesh) << " | " << f2 << " " << edge(h2,mesh) << "\n"; + + std::vector edge_locations; + auto src_bk=src, tgt_bk=tgt; + VGoS::locally_shortest_path(src, tgt, mesh, edge_locations, solver); + + out << edge_locations.size()+2; + out << " " << PMP::construct_point(src, mesh); + for (auto el : edge_locations) + out << " " << PMP::construct_point(el, mesh); + out << " " << PMP::construct_point(tgt, mesh) << "\n"; + out << std::flush; + + std::size_t expected_size = edge_locations.size(); + + assert(edge(get(PMP::get_descriptor_from_location(src,mesh)), mesh)==edge(h1,mesh)); + assert(edge(get(PMP::get_descriptor_from_location(tgt,mesh)), mesh)==edge(h2,mesh)); + + src=src_bk; + tgt=tgt_bk; + edge_locations.clear(); + VGoS::locally_shortest_path(tgt, src, mesh, edge_locations, solver); + + out << edge_locations.size()+2; + out << " " << PMP::construct_point(tgt, mesh); + for (auto el : edge_locations) + out << " " << PMP::construct_point(el, mesh); + out << " " << PMP::construct_point(src, mesh) << "\n"; + out << std::flush; + + CGAL_warning(edge_locations.size()==expected_size); + + assert(edge(get(PMP::get_descriptor_from_location(src,mesh)), mesh)==edge(h1,mesh)); + assert(edge(get(PMP::get_descriptor_from_location(tgt,mesh)), mesh)==edge(h2,mesh)); + + h = next(h, mesh); + } + #endif + + return 0; +} diff --git a/Vector_graphics_on_surfaces/test/Vector_graphics_on_surfaces/test_Bsurf_grid.cpp b/Vector_graphics_on_surfaces/test/Vector_graphics_on_surfaces/test_Bsurf_grid.cpp new file mode 100644 index 00000000000..a7a2be5aca5 --- /dev/null +++ b/Vector_graphics_on_surfaces/test/Vector_graphics_on_surfaces/test_Bsurf_grid.cpp @@ -0,0 +1,164 @@ +#include +#include +#include +#include +#include + +#include + +namespace VGoS = CGAL::Vector_graphics_on_surfaces; +namespace PMP = CGAL::Polygon_mesh_processing; + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Surface_mesh Mesh; +typedef PMP::Face_location Face_location; +typedef PMP::Edge_location Edge_location; + + +//TODO: move to header and add doc +template +void simplify_sequence(std::vector>& edge_locations, const TriangleMesh& tm) +{ + std::cout << " simplify_sequence:\n"; + std::cout << " input.size() = " << edge_locations.size() << "\n"; + typedef boost::graph_traits BGT; + typedef typename BGT::vertex_descriptor vertex_descriptor; + + auto print_seq = [&]() + { + std::cout << " "; + typedef typename BGT::halfedge_descriptor halfedge_descriptor; + for (auto el : edge_locations) + { + auto var = PMP::get_descriptor_from_location(el, tm); + if (const vertex_descriptor* vd_ptr = std::get_if(&var)) + std::cout << " " << *vd_ptr; + else + { + if (const halfedge_descriptor* ed_ptr = std::get_if(&var)) + std::cout << " " << *ed_ptr; + else + std::cout << " ERROR"; + } + } + std::cout << "\n"; + }; + print_seq(); + + + std::vector> filtered_locations; + std::size_t end=edge_locations.size(); + filtered_locations.reserve(edge_locations.size()); + for(std::size_t curr=0; curr(&var)) + { + do{ + vertex_descriptor vd = *vd_ptr; + ++curr; + if (curr!=end) + { + var = PMP::get_descriptor_from_location(edge_locations[curr], tm); + vd_ptr = std::get_if(&var); + if (vd_ptr==nullptr || *vd_ptr!=vd) + break; + } + else + break; + } + while(true); + } + else + ++curr; + } + + if (filtered_locations.size()!=edge_locations.size()) + filtered_locations.swap(edge_locations); + + std::cout << " output.size() = " << edge_locations.size() << "\n"; + print_seq(); +} + +void test_vertex_pair(unsigned int vid1, unsigned int vid2, const Mesh& mesh) +{ + Mesh::Vertex_index v1(vid1), v2(vid2); + std::cout << "running " << v1 << " " << v2 << "\n"; + static int rid=-1; + std::string fname="grid_shortest_path"+std::to_string(++rid)+".polylines.txt"; + std::cout << " writing in " << fname << "\n"; + + std::ofstream out(fname); + out << std::setprecision(17); + + Mesh::Face_index f1=face(halfedge(v1, mesh), mesh), f2 = face(halfedge(v2, mesh), mesh); + + Face_location src(f1, CGAL::make_array(0.,0.,0.)); + Face_location tgt(f2, CGAL::make_array(0.,0.,0.)); + + Mesh::Halfedge_index hf1 = prev(halfedge(f1, mesh), mesh); + int i1=0; + while (target(hf1, mesh)!= v1) + { + hf1=next(hf1,mesh); + ++i1; + } + src.second[i1]=1; + Mesh::Halfedge_index hf2 = prev(halfedge(f2, mesh), mesh); + int i2=0; + while (target(hf2, mesh)!= v2) + { + hf2=next(hf2,mesh); + ++i2; + } + tgt.second[i2]=1; + + std::vector edge_locations; + auto src_bk=src, tgt_bk=tgt; + VGoS::locally_shortest_path(src, tgt, mesh, edge_locations); + assert(get(PMP::get_descriptor_from_location(src,mesh))==v1); + assert(get(PMP::get_descriptor_from_location(tgt,mesh))==v2); + + simplify_sequence(edge_locations, mesh); + + out << edge_locations.size()+2; + out << " " << PMP::construct_point(src, mesh); + for (auto el : edge_locations) + out << " " << PMP::construct_point(el, mesh); + out << " " << PMP::construct_point(tgt, mesh) << "\n"; + out << std::flush; + std::size_t expected_size = edge_locations.size(); + + /// other direction + + src=src_bk; + tgt=tgt_bk; + edge_locations.clear(); + VGoS::locally_shortest_path(tgt, src, mesh, edge_locations); + assert(get(PMP::get_descriptor_from_location(src,mesh))==v1); + assert(get(PMP::get_descriptor_from_location(tgt,mesh))==v2); + + out << edge_locations.size()+2; + out << " " << PMP::construct_point(tgt, mesh); + for (auto el : edge_locations) + out << " " << PMP::construct_point(el, mesh); + out << " " << PMP::construct_point(src, mesh) << "\n"; + out << std::flush; + CGAL_assertion(edge_locations.size() == expected_size); +} + +int main() +{ + Mesh tmp; + CGAL::make_grid(10, 10, tmp, [](int i, int j){return K::Point_3(i,j,0);}, true); + Mesh mesh; + CGAL::Polygon_mesh_processing::extrude_mesh(tmp, mesh, K::Vector_3(0,0,-1)); + std::ofstream("grid.off") << mesh; + + test_vertex_pair(84,36,mesh); + test_vertex_pair(40,36,mesh); + test_vertex_pair(78,28,mesh); + + return 0; +} diff --git a/Vector_graphics_on_surfaces/test/Vector_graphics_on_surfaces/test_straightest_geodesic.cpp b/Vector_graphics_on_surfaces/test/Vector_graphics_on_surfaces/test_straightest_geodesic.cpp new file mode 100644 index 00000000000..92297b74e69 --- /dev/null +++ b/Vector_graphics_on_surfaces/test/Vector_graphics_on_surfaces/test_straightest_geodesic.cpp @@ -0,0 +1,108 @@ +#include +#include +#include +#include +#include + +#include + +namespace PMP = CGAL::Polygon_mesh_processing; +namespace VGoS = CGAL::Vector_graphics_on_surfaces; + +using K = CGAL::Exact_predicates_inexact_constructions_kernel; +using Mesh = CGAL::Surface_mesh; +using Face_location = PMP::Face_location; +using Edge_location = PMP::Edge_location; + +int main() +{ + Mesh mesh; + Mesh::Halfedge_index h = CGAL::make_quad(K::Point_3(0,0,0), K::Point_3(1,0,0), + K::Point_3(1,1,0),K::Point_3(0,1,0), + mesh); + CGAL::Euler::split_face(h, next(next(h,mesh), mesh), mesh); + PMP::isotropic_remeshing(faces(mesh), 0.15, mesh); + + std::ofstream("square.off") << std::setprecision(17) << mesh; + + // test starting on edges + //int runid=0; + for (Mesh::Halfedge_index h : halfedges(mesh)) + { + if (is_border(h,mesh)) continue; + + Mesh::Face_index f = face(h, mesh); + Face_location src(f, CGAL::make_array(0.5, 0.5, 0)); + + // K::Point_3 src_pt = PMP::construct_point(src, mesh); + // std::cout << "src = " << src_pt << "\n"; + + double target_distance = 0.3; + + //std::ofstream out("straightest_geodesic_path_"+std::to_string(runid)+".polylines.txt"); + for (double n=0; n<8; ++n) + { + double theta = 2*CGAL_PI/8*n; + K::Vector_2 dir(std::cos(theta), std::sin(theta)); + std::vector path = VGoS::straightest_geodesic(src, dir, target_distance, mesh); + + //TODO: check the output is of correct length + +/* + std::vector poly; + poly.reserve(path.size()); + VGoS::convert_path_to_polyline(path, mesh, std::back_inserter(poly)); + + + out << path.size() << " "; + + for (auto p : poly) + out << " " << p; + out << "\n"; +*/ + } +// ++runid; + } + + // test starting on vertices + //int runid=0; + for (Mesh::Halfedge_index h : halfedges(mesh)) + { + if (is_border(h,mesh)) continue; + + Mesh::Face_index f = face(h, mesh); + Face_location src(f, CGAL::make_array(0., 1., 0.)); + + // K::Point_3 src_pt = PMP::construct_point(src, mesh); + // std::cout << "src = " << src_pt << "\n"; + + double target_distance = 0.3; + + //std::ofstream out("straightest_geodesic_path_"+std::to_string(runid)+".polylines.txt"); + for (double n=0; n<8; ++n) + { + double theta = 2*CGAL_PI/16*n; + K::Vector_2 dir(std::cos(theta), std::sin(theta)); + std::vector path = VGoS::straightest_geodesic(src, dir, target_distance, mesh); + + //TODO: check the output is of correct length + +/* + std::vector poly; + poly.reserve(path.size()); + VGoS::convert_path_to_polyline(path, mesh, std::back_inserter(poly)); + + + out << path.size() << " "; + + for (auto p : poly) + out << " " << p; + out << "\n"; +*/ + } +// ++runid; + } + + //TODO: we still need to check that the output is correct (which is currently not the case at vertex) + return 1; +}