1 #include <CGAL/Three/Polyhedron_demo_plugin_interface.h>
2 #include <CGAL/Three/Scene_interface.h>
3 #include <CGAL/Three/Scene_item_rendering_helper.h>
4 #include <CGAL/Three/Triangle_container.h>
5 #include <CGAL/Three/Edge_container.h>
6 #include <CGAL/Three/Three.h>
7 #include <QApplication>
8 #include <QObject>
9 #include <QAction>
10 #include <QMainWindow>
11 #include <QInputDialog>
12 #include <QMessageBox>
13 #include <QMap>
14 #include "Messages_interface.h"
15 #include "Kernel_type.h"
16 #include "Scene_surface_mesh_item.h"
17 #include "Color_ramp.h"
18 #include "triangulate_primitive.h"
19 #include <CGAL/Polygon_mesh_processing/bbox.h>
20 #include <CGAL/Polygon_mesh_processing/distance.h>
21 #include <CGAL/Polygon_mesh_processing/compute_normal.h>
22 #include <CGAL/boost/iterator/counting_iterator.hpp>
23 #include <CGAL/Search_traits_3.h>
24 #include <CGAL/Spatial_sort_traits_adapter_3.h>
25 #include <CGAL/property_map.h>
26 #include <boost/container/flat_map.hpp>
27 using namespace CGAL::Three;
28 namespace PMP = CGAL::Polygon_mesh_processing;
29 typedef Viewer_interface Vi;
30 typedef Triangle_container Tc;
31 typedef Edge_container Ec;
32 
33 typedef Scene_surface_mesh_item Scene_face_graph_item;
34 typedef Scene_face_graph_item::Face_graph Face_graph;
35 
36 #if defined(CGAL_LINKED_WITH_TBB)
37 #include <tbb/parallel_for.h>
38 template <class AABB_tree, class Point_3>
39 struct Distance_computation{
40   const AABB_tree& tree;
41   const std::vector<Point_3>& sample_points;
42   const Point_3 initial_hint;
43   std::vector<double>& output;
44 
Distance_computationDistance_computation45   Distance_computation(const AABB_tree& tree,
46                        const Point_3 p,
47                        const std::vector<Point_3>& sample_points,
48                        std::vector<double>& out )
49     : tree(tree)
50     , sample_points(sample_points)
51     , initial_hint(p)
52     , output(out)
53   {
54   }
55   void
operator ()Distance_computation56   operator()(const tbb::blocked_range<std::size_t>& range) const
57   {
58     Point_3 hint = initial_hint;
59     for( std::size_t i = range.begin(); i != range.end(); ++i)
60     {
61       hint = tree.closest_point(sample_points[i], hint);
62       Kernel::FT dist = squared_distance(hint,sample_points[i]);
63       double d = CGAL::sqrt(dist);
64       output[i] = d;
65     }
66   }
67 };
68 #endif
69 
70 class Scene_distance_polyhedron_item: public Scene_item_rendering_helper
71 {
72   Q_OBJECT
73 public:
Scene_distance_polyhedron_item(Face_graph * poly,Face_graph * polyB,QString other_name,int sampling_pts)74   Scene_distance_polyhedron_item(Face_graph* poly, Face_graph* polyB, QString other_name, int sampling_pts)
75     : poly(poly),
76       poly_B(polyB),
77       other_poly(other_name)
78   {
79     setTriangleContainer(0, new Tc(Vi::PROGRAM_WITH_LIGHT,
80                                    false));
81     setEdgeContainer(0, new Ec(Vi::PROGRAM_NO_SELECTION, false));
82     nb_pts_per_face = sampling_pts;
83     this->setRenderingMode(FlatPlusEdges);
84     thermal_ramp.build_thermal();
85   }
supportsRenderingMode(RenderingMode m) const86   bool supportsRenderingMode(RenderingMode m) const {
87     return (m == Flat || m == FlatPlusEdges);
88   }
clone() const89   Scene_item* clone() const {return nullptr;}
toolTip() const90   QString toolTip() const {return QString("Item %1 with color indicating distance with %2").arg(this->name()).arg(other_poly);}
draw(Viewer_interface * viewer) const91   void draw(Viewer_interface *viewer) const
92   {
93     if(!isInit(viewer))
94       initGL(viewer);
95     if ( getBuffersFilled() &&
96          ! getBuffersInit(viewer))
97     {
98       initializeBuffers(viewer);
99       setBuffersInit(viewer, true);
100     }
101     getTriangleContainer(0)->draw(viewer, false);
102   }
103 
drawEdges(Viewer_interface * viewer) const104   void drawEdges(Viewer_interface* viewer) const
105   {
106     getEdgeContainer(0)->setColor(QColor(Qt::black));
107     getEdgeContainer(0)->draw(viewer, true);
108   }
109 
compute_bbox() const110   void compute_bbox() const {
111     setBbox(PMP::bbox(*poly));
112   }
113 
114 private:
115   Face_graph* poly;
116   Face_graph* poly_B;
117   QString other_poly;
118   mutable std::vector<float> m_vertices;
119   mutable std::vector<float> edge_vertices;
120   mutable std::vector<float> normals;
121   mutable std::vector<float> colors;
122   Color_ramp thermal_ramp;
123   int nb_pts_per_face;
124 
125   mutable std::size_t nb_pos;
126   mutable std::size_t nb_edge_pos;
127   //fills 'out' and returns the hausdorff distance for calibration of the color_ramp.
128 
compute_distances(const Face_graph & m,const std::vector<Kernel::Point_3> & sample_points,std::vector<double> & out) const129   double compute_distances(const Face_graph& m, const std::vector<Kernel::Point_3>& sample_points,
130                            std::vector<double>& out)const
131   {
132     typedef CGAL::AABB_face_graph_triangle_primitive<Face_graph> Primitive;
133     typedef CGAL::AABB_traits<Kernel, Primitive> Traits;
134     typedef CGAL::AABB_tree< Traits > Tree;
135 
136     Tree tree( faces(m).first, faces(m).second, m);
137     tree.build();
138     boost::graph_traits<Face_graph>::vertex_descriptor vd = *(vertices(m).first);
139     Traits::Point_3 hint = get(CGAL::vertex_point,m, vd);
140 
141 #if !defined(CGAL_LINKED_WITH_TBB)
142     double hdist = 0;
143     for(std::size_t i = 0; i<sample_points.size(); ++i)
144     {
145       hint = tree.closest_point(sample_points[i], hint);
146       Kernel::FT dist = squared_distance(hint,sample_points[i]);
147       double d = CGAL::sqrt(dist);
148       out[i]= d;
149       if (d>hdist) hdist=d;
150     }
151       return hdist;
152 #else
153     double distance=0;
154     Distance_computation<Tree, Kernel::Point_3> f(tree, hint, sample_points, out);
155     tbb::parallel_for(tbb::blocked_range<std::size_t>(0, sample_points.size()), f);
156     for(std::size_t i = 0; i< out.size(); ++i){
157       if(out[i] > distance)
158         distance = out[i];
159     }
160     return distance;
161 #endif
162   }
163 
computeElements() const164   void computeElements()const
165   {
166     QApplication::setOverrideCursor(Qt::WaitCursor);
167     m_vertices.resize(0);
168     edge_vertices.resize(0);
169     normals.resize(0);
170     colors.resize(0);
171 
172     typedef Kernel::Vector_3            Vector;
173     typedef boost::graph_traits<Face_graph>::face_descriptor   face_descriptor;
174     typedef boost::graph_traits<Face_graph>::vertex_descriptor vertex_descriptor;
175 
176     typedef boost::property_map<Face_graph,CGAL::vertex_point_t>::type VPmap;
177     VPmap vpmap = get(CGAL::vertex_point,*poly);
178 
179     //facets
180     {
181       boost::container::flat_map<face_descriptor, Vector> face_normals_map;
182       boost::associative_property_map< boost::container::flat_map<face_descriptor, Vector> >
183           nf_pmap(face_normals_map);
184       boost::container::flat_map<vertex_descriptor, Vector> vertex_normals_map;
185       boost::associative_property_map< boost::container::flat_map<vertex_descriptor, Vector> >
186           nv_pmap(vertex_normals_map);
187 
188       PMP::compute_normals(*poly, nv_pmap, nf_pmap);
189       std::vector<Kernel::Point_3> total_points(0);
190 
191       for(boost::graph_traits<Face_graph>::face_descriptor f : faces(*poly)) {
192         Vector nf = get(nf_pmap, f);
193         typedef FacetTriangulator<Face_graph, Kernel, boost::graph_traits<Face_graph>::vertex_descriptor> FT;
194 
195         //compute distance with other polyhedron
196         //sample facet
197         std::vector<Kernel::Point_3> sampled_points;
198         std::size_t nb_points =  (std::max)((int)std::ceil(nb_pts_per_face * PMP::face_area(f,*poly,PMP::parameters::geom_traits(Kernel()))),
199                                             1);
200         Kernel::Point_3 &p = get(vpmap,target(halfedge(f,*poly),*poly));
201         Kernel::Point_3 &q = get(vpmap,target(next(halfedge(f,*poly),*poly),*poly));
202         Kernel::Point_3 &r = get(vpmap,target(next(next(halfedge(f,*poly),*poly),*poly),*poly));
203         CGAL::Random_points_in_triangle_3<Kernel::Point_3> g(p, q, r);
204         std::copy_n(g, nb_points, std::back_inserter(sampled_points));
205         sampled_points.push_back(p);
206         sampled_points.push_back(q);
207         sampled_points.push_back(r);
208 
209         //triangle facets with sample points for color display
210         FT triangulation(f,sampled_points,nf,poly);
211 
212         if(triangulation.cdt->dimension() != 2 )
213         {
214           qDebug()<<"Error : cdt not right (dimension != 2). Facet not displayed";
215           continue;
216         }
217 
218         //iterates on the internal faces to add the vertices to the positions
219         //and the normals to the appropriate vectors
220 
221         for(FT::CDT::Finite_faces_iterator
222             ffit = triangulation.cdt->finite_faces_begin(),
223             end = triangulation.cdt->finite_faces_end();
224             ffit != end; ++ffit)
225         {
226           if(ffit->info().is_external)
227             continue;
228 
229           for (int i = 0; i<3; ++i)
230           {
231             total_points.push_back(ffit->vertex(i)->point());
232             m_vertices.push_back(ffit->vertex(i)->point().x());
233             m_vertices.push_back(ffit->vertex(i)->point().y());
234             m_vertices.push_back(ffit->vertex(i)->point().z());
235 
236             normals.push_back(nf.x());
237             normals.push_back(nf.y());
238             normals.push_back(nf.z());
239           }
240         }
241       }
242       //compute the distances
243       typedef CGAL::Spatial_sort_traits_adapter_3<Kernel,
244                 CGAL::Pointer_property_map<Kernel::Point_3>::type > Search_traits_3;
245 
246       std::vector<double> distances(total_points.size());
247       std::vector<std::size_t> indices;
248       indices.reserve(total_points.size());
249       std::copy(boost::counting_iterator<std::size_t>(0),
250                 boost::counting_iterator<std::size_t>(total_points.size()),
251                 std::back_inserter(indices));
252       spatial_sort(indices.begin(),
253                    indices.end(),
254                   Search_traits_3(CGAL::make_property_map(total_points)));
255       std::vector<Kernel::Point_3> sorted_points(total_points.size());
256       for(std::size_t i = 0; i < sorted_points.size(); ++i)
257       {
258         sorted_points[i] = total_points[indices[i]];
259       }
260 
261       double hausdorff = compute_distances(*poly_B,
262                                            sorted_points,
263                                            distances);
264       if(hausdorff == 0)
265         hausdorff++;
266       //compute the colors
267       colors.resize(sorted_points.size()*3);
268       for(std::size_t i=0; i<sorted_points.size(); ++i)
269       {
270         std::size_t k = indices[i];
271         double d = distances[i]/hausdorff;
272         colors[3*k]=thermal_ramp.r(d);
273         colors[3*k+1]=thermal_ramp.g(d);
274         colors[3*k+2]=thermal_ramp.b(d);
275       }
276     }
277 
278     //edges
279     {
280       //Lines
281       typedef Kernel::Point_3                Point;
282       typedef boost::graph_traits<Face_graph>::edge_descriptor        edge_descriptor;
283 
284       for(edge_descriptor he : edges(*poly)){
285         const Point& a = get(vpmap,target(he,*poly));
286         const Point& b = get(vpmap,source(he,*poly));
287         {
288 
289           edge_vertices.push_back(a.x());
290           edge_vertices.push_back(a.y());
291           edge_vertices.push_back(a.z());
292 
293           edge_vertices.push_back(b.x());
294           edge_vertices.push_back(b.y());
295           edge_vertices.push_back(b.z());
296         }
297       }
298     }
299 
300     Tc* tc = getTriangleContainer(0);
301     Ec* ec = getEdgeContainer(0);
302 
303     tc->allocate(
304           Tc::Flat_vertices,
305           m_vertices.data(),
306           static_cast<GLsizei>(m_vertices.size()*sizeof(float)));
307     tc->allocate(Tc::Flat_normals,
308                  normals.data(),
309                  static_cast<GLsizei>(normals.size()*sizeof(float)));
310     tc->allocate(Tc::FColors,
311                  colors.data(),
312                  static_cast<GLsizei>(colors.size()*sizeof(float)));
313 
314     ec->allocate(Ec::Vertices,
315                  edge_vertices.data(),
316                  static_cast<GLsizei>(edge_vertices.size()*sizeof(float)));
317     nb_pos = m_vertices.size();
318     nb_edge_pos = edge_vertices.size();
319     compute_bbox();
320     setBuffersFilled(true);
321     QApplication::restoreOverrideCursor();
322   }
323 
initializeBuffers(Viewer_interface * viewer) const324   void initializeBuffers(Viewer_interface *viewer)const
325   {
326     Tc* tc = getTriangleContainer(0);
327     Ec* ec = getEdgeContainer(0);
328     tc->initializeBuffers(viewer);
329     ec->initializeBuffers(viewer);
330     tc->setFlatDataSize(nb_pos);
331     ec->setFlatDataSize(nb_edge_pos);
332     m_vertices.clear();
333     edge_vertices.clear();
334     normals.clear();
335     colors.clear();
336     m_vertices.shrink_to_fit();
337     edge_vertices.shrink_to_fit();
338     normals.shrink_to_fit();
339     colors.shrink_to_fit();
340   }
341 };
342 class DistancePlugin :
343     public QObject,
344     public Polyhedron_demo_plugin_interface
345 {
346   Q_OBJECT
347   Q_INTERFACES(CGAL::Three::Polyhedron_demo_plugin_interface)
348   Q_PLUGIN_METADATA(IID "com.geometryfactory.PolyhedronDemo.PluginInterface/1.0" FILE "distance_plugin.json")
349 
350   typedef Kernel::Point_3 Point_3;
351 public:
352   //decides if the plugin's actions will be displayed or not.
applicable(QAction *) const353   bool applicable(QAction*) const
354   {
355     return scene->selectionIndices().size() == 2 &&
356         qobject_cast<Scene_face_graph_item*>(scene->item(scene->selectionIndices().first())) &&
357         qobject_cast<Scene_face_graph_item*>(scene->item(scene->selectionIndices().last()));
358   }
359   //the list of the actions of the plugin.
actions() const360   QList<QAction*> actions() const
361   {
362     return _actions;
363   }
364   //this acts like a constructor for the plugin. It gets the references to the mainwindow and the scene, and connects the action.
init(QMainWindow * mw,Scene_interface * sc,Messages_interface * mi)365   void init(QMainWindow* mw, Scene_interface* sc, Messages_interface* mi)
366   {
367     //gets the reference to the message interface, to display text in the console widget
368     this->messageInterface = mi;
369     //get the references
370     this->scene = sc;
371     this->mw = mw;
372     //creates the action
373     QAction *actionComputeDistance= new QAction(QString("Compute Distance Between Polyhedra"), mw);
374     //specifies the subMenu
375     actionComputeDistance->setProperty("subMenuName", "Polygon Mesh Processing");
376     //links the action
377     if(actionComputeDistance) {
378       connect(actionComputeDistance, SIGNAL(triggered()),
379               this, SLOT(createDistanceItems()));
380       _actions << actionComputeDistance;
381     }
382   }
383 public Q_SLOTS:
createDistanceItems()384   void createDistanceItems()
385   {
386     bool ok = false;
387     nb_pts_per_face = QInputDialog::getInt(mw, tr("Sampling"),
388                                                tr("Number of points per face:"),40, 1,2147483647,1, &ok);
389     if (!ok)
390       return;
391 
392     //check the initial conditions
393     Scene_face_graph_item* itemA = qobject_cast<Scene_face_graph_item*>(scene->item(scene->selectionIndices().first()));
394     Scene_face_graph_item* itemB = qobject_cast<Scene_face_graph_item*>(scene->item(scene->selectionIndices().last()));
395     if(! CGAL::is_triangle_mesh(*itemA->polyhedron()) ||
396        !CGAL::is_triangle_mesh(*itemB->polyhedron()) ){
397       CGAL::Three::Three::error(QString("Distance not computed. (Both polyhedra must be triangulated)"));
398       return;
399     }
400     QApplication::setOverrideCursor(Qt::WaitCursor);
401     Scene_distance_polyhedron_item* new_itemA = new Scene_distance_polyhedron_item(itemA->polyhedron(),itemB->polyhedron(), itemB->name(), nb_pts_per_face);
402     Scene_distance_polyhedron_item* new_itemB = new Scene_distance_polyhedron_item(itemB->polyhedron(),itemA->polyhedron(), itemA->name(), nb_pts_per_face);
403     itemA->setVisible(false);
404     itemB->setVisible(false);
405     new_itemA->setName(QString("%1 to %2").arg(itemA->name()).arg(itemB->name()));
406     new_itemB->setName(QString("%1 to %2").arg(itemB->name()).arg(itemA->name()));
407     scene->addItem(new_itemA);
408     scene->addItem(new_itemB);
409     QApplication::restoreOverrideCursor();
410   }
411 private:
412   int nb_pts_per_face;
413   QList<QAction*> _actions;
414   Messages_interface* messageInterface;
415   //The reference to the scene
416   Scene_interface* scene;
417   //The reference to the main window
418   QMainWindow* mw;
419 };
420 #include "Distance_plugin.moc"
421