1 /*
2 * Copyright (c) 2012-2014, Bruno Levy
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * * Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 * * Redistributions in binary form must reproduce the above copyright notice,
11 * this list of conditions and the following disclaimer in the documentation
12 * and/or other materials provided with the distribution.
13 * * Neither the name of the ALICE Project-Team nor the names of its
14 * contributors may be used to endorse or promote products derived from this
15 * software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 * POSSIBILITY OF SUCH DAMAGE.
28 *
29 * If you modify this software, you should include a notice giving the
30 * name of the person performing the modification, the date of modification,
31 * and the reason for such modification.
32 *
33 * Contact: Bruno Levy
34 *
35 * Bruno.Levy@inria.fr
36 * http://www.loria.fr/~levy
37 *
38 * ALICE Project
39 * LORIA, INRIA Lorraine,
40 * Campus Scientifique, BP 239
41 * 54506 VANDOEUVRE LES NANCY CEDEX
42 * FRANCE
43 *
44 */
45
46 #include <geogram_gfx/gui/simple_application.h>
47 #include <geogram_gfx/GLUP/GLUP_private.h>
48 #include <geogram/delaunay/periodic_delaunay_3d.h>
49
50
51 namespace {
52
53 using namespace GEO;
54
55 /**
56 * \brief An application that demonstrates both
57 * GLUP primitives and glup_viewer application
58 * framework.
59 */
60 class DemoDelaunay3dApplication : public SimpleApplication {
61 public:
62
63 /**
64 * \brief DemoDelaunay3dApplication constructor.
65 */
DemoDelaunay3dApplication()66 DemoDelaunay3dApplication() : SimpleApplication("Delaunay3d") {
67
68 periodic_ = false;
69
70 // Define the 3d region that we want to display
71 // (xmin, ymin, zmin, xmax, ymax, zmax)
72 set_region_of_interest(0.0, 0.0, 0.0, 1.0, 1.0, 1.0);
73
74 draw_points_ = false;
75 draw_cells_ = true;
76 point_size_ = 10.0f;
77 cells_shrink_ = 0.1f;
78 nb_points_ = 100;
79 draw_box_ = false;
80 draw_period_ = false;
81
82 start_animation();
83 }
84
geogram_initialize(int argc,char ** argv)85 void geogram_initialize(int argc, char** argv) override {
86 SimpleApplication::geogram_initialize(argc, argv);
87 delaunay_ = new PeriodicDelaunay3d(periodic_, 1.0);
88 // In non-periodic mode, we need to keep the infinite
89 // tetrahedra to be able to query the combinatorics of
90 // all cells (including the infinite ones).
91 if(!periodic_) {
92 delaunay_->set_keeps_infinite(true);
93 }
94 init_random_points(nb_points_);
95 }
96
97 /**
98 * \brief DemoDelaunay3dApplication destructor.
99 */
~DemoDelaunay3dApplication()100 ~DemoDelaunay3dApplication() override {
101 }
102
103 /**
104 * \brief Displays and handles the GUI for object properties.
105 * \details Overloads Application::draw_object_properties().
106 */
draw_object_properties()107 void draw_object_properties() override {
108 SimpleApplication::draw_object_properties();
109 ImGui::Checkbox("Animate", animate_ptr());
110
111 if(ImGui::Button("One iteration", ImVec2(-1.0, 0.0))) {
112 Lloyd_iteration();
113 }
114
115 if(
116 ImGui::Button(
117 "+",
118 ImVec2(-ImGui::GetContentRegionAvail().x/2.0f,0.0f)
119 ) && nb_points_ < 10000
120 ) {
121 ++nb_points_;
122 init_random_points(nb_points_);
123 }
124 ImGui::SameLine();
125 if(
126 ImGui::Button(
127 "-", ImVec2(-1.0f, 0.0f)
128 ) && nb_points_ > 4
129 ) {
130 --nb_points_;
131 init_random_points(nb_points_);
132 }
133
134 if(ImGui::SliderInt("N", &nb_points_, 4, 10000)) {
135 init_random_points(nb_points_);
136 }
137
138 if(ImGui::Button("Reset", ImVec2(-1.0, 0.0))) {
139 init_random_points(nb_points_);
140 }
141
142 if(ImGui::Checkbox("Periodic", &periodic_)) {
143 delaunay_ = new PeriodicDelaunay3d(periodic_, 1.0);
144 // In non-periodic mode, we need to keep the infinite
145 // tetrahedra to be able to query the combinatorics of
146 // all cells (including the infinite ones).
147 if(!periodic_) {
148 delaunay_->set_keeps_infinite(true);
149 }
150 delaunay_->set_vertices(points_.size()/3, points_.data());
151 delaunay_->compute();
152 }
153
154 ImGui::Spacing();
155 ImGui::Separator();
156 ImGui::Spacing();
157
158 ImGui::Checkbox("Box", &draw_box_);
159 if(ImGui::Checkbox("Period 3x3x3", &draw_period_)) {
160 if(draw_period_) {
161 set_region_of_interest(
162 -1.0, -1.0, -1.0, 2.0, 2.0, 2.0
163 );
164 } else {
165 set_region_of_interest(
166 0.0, 0.0, 0.0, 1.0, 1.0, 1.0
167 );
168 }
169 }
170 ImGui::Checkbox("Points", &draw_points_);
171 ImGui::SameLine();
172 ImGui::SliderFloat("##PtSz.", &point_size_, 1.0f, 50.0f, "%.1f");
173 ImGui::Checkbox("Cells", &draw_cells_);
174 ImGui::SameLine();
175 ImGui::SliderFloat("##Shrk.", &cells_shrink_, 0.0f, 1.0f, "%.2f");
176 }
177
178 /**
179 * \brief Gets a Voronoi cell.
180 * \details In non-periodic mode, the cell is clipped by the domain.
181 * \param[in] v the index of the vertex
182 * \param[out] C the cell
183 */
get_cell(index_t v,ConvexCell & C)184 void get_cell(index_t v, ConvexCell& C) {
185 delaunay_->copy_Laguerre_cell_from_Delaunay(v, C, W_);
186 if(!periodic_) {
187 C.clip_by_plane(vec4( 1.0, 0.0, 0.0, 0.0));
188 C.clip_by_plane(vec4(-1.0, 0.0, 0.0, 1.0));
189 C.clip_by_plane(vec4( 0.0, 1.0, 0.0, 0.0));
190 C.clip_by_plane(vec4( 0.0,-1.0, 0.0, 1.0));
191 C.clip_by_plane(vec4( 0.0, 0.0, 1.0, 0.0));
192 C.clip_by_plane(vec4( 0.0, 0.0,-1.0, 1.0));
193 }
194 C.compute_geometry();
195 }
196
197 /**
198 * \brief Draws a cell.
199 * \details Needs to be called between glupBegin(GLUP_TRIANGLES)
200 * and glupEnd().
201 */
draw_cell(ConvexCell & C,index_t instance=0)202 void draw_cell(ConvexCell& C, index_t instance = 0) {
203 double s = double(cells_shrink_);
204 double Tx = double(Periodic::translation[instance][0]);
205 double Ty = double(Periodic::translation[instance][1]);
206 double Tz = double(Periodic::translation[instance][2]);
207
208 vec3 g;
209 if(cells_shrink_ != 0.0f) {
210 g = C.barycenter();
211 }
212
213 // For each vertex of the Voronoi cell
214 // Start at 1 (vertex 0 is point at infinity)
215 for(index_t v=1; v<C.nb_v(); ++v) {
216 index_t t = C.vertex_triangle(v);
217
218 // Happens if a clipping plane did not
219 // clip anything.
220 if(t == VBW::END_OF_LIST) {
221 continue;
222 }
223
224 // Now iterate on the Voronoi vertices of the
225 // Voronoi facet. In dual form this means iterating
226 // on the triangles incident to vertex v
227
228 vec3 P[3];
229 index_t n=0;
230 do {
231
232 // Triangulate the Voronoi facet and send the
233 // triangles to OpenGL/GLUP.
234 if(n == 0) {
235 P[0] = C.triangle_point(VBW::ushort(t));
236 } else if(n == 1) {
237 P[1] = C.triangle_point(VBW::ushort(t));
238 } else {
239 P[2] = C.triangle_point(VBW::ushort(t));
240 if(s == 0.0) {
241 for(index_t i=0; i<3; ++i) {
242 glupPrivateVertex3d(
243 P[i].x + Tx, P[i].y + Ty, P[i].z + Tz
244 );
245 }
246 } else {
247 for(index_t i=0; i<3; ++i) {
248 glupPrivateVertex3d(
249 s*g.x + (1.0-s)*P[i].x + Tx,
250 s*g.y + (1.0-s)*P[i].y + Ty,
251 s*g.z + (1.0-s)*P[i].z + Tz
252 );
253 }
254 }
255 P[1] = P[2];
256 }
257
258 // Move to the next triangle incident to
259 // vertex v.
260 index_t lv = C.triangle_find_vertex(t,v);
261 t = C.triangle_adjacent(t, (lv + 1)%3);
262
263 ++n;
264 } while(t != C.vertex_triangle(v));
265 }
266 }
267
268 /**
269 * \brief Draws the scene according to currently set primitive and
270 * drawing modes.
271 */
draw_scene()272 void draw_scene() override {
273 // Avoid re-entry, for instance when a message is sent to
274 // the logger, this triggers a graphic redisplay.
275 static bool locked = false;
276 if(locked) {
277 return;
278 }
279 locked = true;
280 if(animate()) {
281 Lloyd_iteration();
282 }
283 glupSetColor3f(GLUP_FRONT_AND_BACK_COLOR, 0.5f, 0.5f, 0.5f);
284 glupEnable(GLUP_LIGHTING);
285
286 index_t max_instance = (draw_period_ ? 27 : 1);
287
288
289
290 if(draw_points_) {
291 double R = double(point_size_) / 250.0;
292
293 for(index_t i=0; i<max_instance; ++i) {
294 glupTranslated(
295 double(Periodic::translation[i][0]),
296 double(Periodic::translation[i][1]),
297 double(Periodic::translation[i][2])
298 );
299
300 glupBegin(GLUP_SPHERES);
301 for(index_t v=0; v<points_.size()/3; ++v) {
302 glupPrivateVertex4d(
303 points_[3*v],
304 points_[3*v+1],
305 points_[3*v+2],
306 R
307 );
308 }
309 glupEnd();
310
311 glupTranslated(
312 -double(Periodic::translation[i][0]),
313 -double(Periodic::translation[i][1]),
314 -double(Periodic::translation[i][2])
315 );
316 }
317
318 }
319
320 if(draw_cells_) {
321 ConvexCell C;
322 glupBegin(GLUP_TRIANGLES);
323 for(index_t v=0; v<points_.size()/3; ++v) {
324 get_cell(v, C);
325 if(draw_period_) {
326 for(index_t i=0; i<27; ++i) {
327 draw_cell(C,i);
328 }
329 } else {
330 draw_cell(C);
331 }
332 }
333 glupEnd();
334 }
335
336 if(draw_box_) {
337 static GLfloat cube[8][3] = {
338 {0.0f, 0.0f, 0.0f},
339 {0.0f, 0.0f, 1.0f},
340 {0.0f, 1.0f, 0.0f},
341 {0.0f, 1.0f, 1.0f},
342 {1.0f, 0.0f, 0.0f},
343 {1.0f, 0.0f, 1.0f},
344 {1.0f, 1.0f, 0.0f},
345 {1.0f, 1.0f, 1.0f}
346 };
347
348
349 glupEnable(GLUP_DRAW_MESH);
350 glupEnable(GLUP_ALPHA_DISCARD);
351 glupSetColor4f(
352 GLUP_FRONT_AND_BACK_COLOR, 1.0f, 1.0f, 1.0f, 0.0f
353 );
354 glupSetMeshWidth(10);
355 glupDisable(GLUP_LIGHTING);
356
357 for(index_t i=0; i<max_instance; ++i) {
358 glupTranslated(
359 double(Periodic::translation[i][0]),
360 double(Periodic::translation[i][1]),
361 double(Periodic::translation[i][2])
362 );
363
364 glupBegin(GLUP_QUADS);
365
366 glupVertex3fv(cube[0]);
367 glupVertex3fv(cube[1]);
368 glupVertex3fv(cube[3]);
369 glupVertex3fv(cube[2]);
370
371 glupVertex3fv(cube[4]);
372 glupVertex3fv(cube[5]);
373 glupVertex3fv(cube[7]);
374 glupVertex3fv(cube[6]);
375
376 glupVertex3fv(cube[0]);
377 glupVertex3fv(cube[1]);
378 glupVertex3fv(cube[5]);
379 glupVertex3fv(cube[4]);
380
381 glupVertex3fv(cube[1]);
382 glupVertex3fv(cube[3]);
383 glupVertex3fv(cube[7]);
384 glupVertex3fv(cube[5]);
385
386 glupVertex3fv(cube[3]);
387 glupVertex3fv(cube[2]);
388 glupVertex3fv(cube[6]);
389 glupVertex3fv(cube[7]);
390
391 glupVertex3fv(cube[2]);
392 glupVertex3fv(cube[0]);
393 glupVertex3fv(cube[4]);
394 glupVertex3fv(cube[6]);
395
396 glupEnd();
397
398 glupTranslated(
399 -double(Periodic::translation[i][0]),
400 -double(Periodic::translation[i][1]),
401 -double(Periodic::translation[i][2])
402 );
403 }
404
405
406
407 glupDisable(GLUP_DRAW_MESH);
408 glupDisable(GLUP_ALPHA_DISCARD);
409 }
410
411 locked = false;
412 }
413
414 /**
415 * \brief Initializes the Delaunay triangulation with
416 * random points.
417 * \param[in] nb_points_in number of points
418 */
init_random_points(int nb_points_in)419 void init_random_points(int nb_points_in) {
420 index_t nb_points = index_t(nb_points_in);
421 bool bkp = draw_cells_;
422 draw_cells_ = false;
423 points_.resize(3*nb_points);
424 for(index_t i=0; i<3*nb_points; ++i) {
425 points_[i] = Numeric::random_float64();
426 }
427 delaunay_->set_vertices(nb_points, points_.data());
428 delaunay_->compute();
429 draw_cells_ = bkp;
430 }
431
432 /**
433 * \brief Does one iteration of Lloyd relaxation.
434 * \details Moves each point to the centroid of its
435 * Voronoi cell. In non-periodic mode, the Voronoi
436 * cell is clipped by the domain.
437 */
Lloyd_iteration()438 void Lloyd_iteration() {
439 vector<double> new_points(points_.size());
440 ConvexCell C;
441 for(index_t v=0; v<points_.size()/3; ++v) {
442 get_cell(v, C);
443 vec3 g = C.barycenter();
444 new_points[3*v] = g.x;
445 new_points[3*v+1] = g.y;
446 new_points[3*v+2] = g.z;
447 }
448 // In periodic mode, points may escape out of
449 // the domain. Relocate them in [0,1]^3
450 for(index_t i=0; i<new_points.size(); ++i) {
451 if(new_points[i] < 0.0) {
452 new_points[i] += 1.0;
453 }
454 if(new_points[i] > 1.0) {
455 new_points[i] -= 1.0;
456 }
457 }
458 points_.swap(new_points);
459 delaunay_->set_vertices(points_.size()/3, points_.data());
460 delaunay_->compute();
461 }
462
463 private:
464 SmartPointer<PeriodicDelaunay3d> delaunay_;
465 bool periodic_;
466 bool draw_points_;
467 float point_size_;
468 bool draw_cells_;
469 bool draw_box_;
470 bool draw_period_;
471 float cells_shrink_;
472 vector<double> points_;
473 int nb_points_;
474 PeriodicDelaunay3d::IncidentTetrahedra W_;
475 };
476
477 }
478
main(int argc,char ** argv)479 int main(int argc, char** argv) {
480 DemoDelaunay3dApplication app;
481 app.start(argc, argv);
482 return 0;
483 }
484