#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
#include <CGAL/Barycentric_coordinates_3/Mean_value_coordinates_3.h>
using Surface_mesh = CGAL::Surface_mesh<Point_3>;
namespace PMP = CGAL::Polygon_mesh_processing;
int main(){
Surface_mesh concave;
const Point_3 p0(0, 3, 0);
const Point_3 p1(1, 1, 0);
const Point_3 p2(3, 0, 0);
const Point_3 p3(0, 0, 0);
const Point_3 p4(0, 0, 3);
const Point_3 p5(0, 3, 3);
const Point_3 p6(1, 1, 3);
const Point_3 p7(3, 0, 3);
PMP::triangulate_faces(faces(concave), concave);
std::vector<FT> coords;
std::vector<Point_3> queries{
Point_3(FT(1)/FT(2), FT(1)/FT(2), FT(1)), Point_3(FT(1)/FT(3), FT(1)/FT(3), FT(2)),
Point_3(FT(4)/FT(3), FT(1)/FT(3), FT(1)), Point_3(FT(4)/FT(3), FT(1)/FT(3), FT(2)),
Point_3(FT(1)/FT(3), FT(4)/FT(3), FT(1)), Point_3(FT(1)/FT(3), FT(4)/FT(3), FT(2))};
std::cout << std::endl << "Mean value coordinates : " << std::endl << std::endl;
for (const auto& query : queries){
coords.clear();
concave, query, std::back_inserter(coords), CP3::FAST);
for (std::size_t i = 0; i < coords.size() -1; ++i) {
std::cout << coords[i] << ", ";
}
std::cout << coords[coords.size() -1] << std::endl;
}
std::cout << std::endl;
return EXIT_SUCCESS;
}