summaryrefslogtreecommitdiff
path: root/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h
blob: 7f96c94ca34d54b6658466ab1ffa848d09c54bab (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
/*    This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
 *    See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
 *    Author(s):       Vincent Rouvreau
 *
 *    Copyright (C) 2018 Inria
 *
 *    Modification(s):
 *      - 2019/08 Vincent Rouvreau: Fix issue #10 for CGAL and Eigen3
 *      - YYYY/MM Author: Description of the modification
 */

#ifndef ALPHA_COMPLEX_3D_H_
#define ALPHA_COMPLEX_3D_H_

#include <boost/version.hpp>
#include <boost/variant.hpp>

#include <gudhi/Debug_utils.h>
#include <gudhi/Alpha_complex_options.h>

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_3.h>
#include <CGAL/Periodic_3_regular_triangulation_traits_3.h>
#include <CGAL/Periodic_3_regular_triangulation_3.h>
#include <CGAL/Regular_triangulation_3.h>
#include <CGAL/Alpha_shape_3.h>
#include <CGAL/Alpha_shape_cell_base_3.h>
#include <CGAL/Alpha_shape_vertex_base_3.h>

#include <CGAL/Object.h>
#include <CGAL/tuple.h>
#include <CGAL/iterator.h>
#include <CGAL/version.h>  // for CGAL_VERSION_NR

#include <Eigen/src/Core/util/Macros.h>  // for EIGEN_VERSION_AT_LEAST

#include <boost/container/static_vector.hpp>

#include <iostream>
#include <vector>
#include <unordered_map>
#include <stdexcept>
#include <cstddef>  // for std::size_t
#include <memory>       // for std::unique_ptr
#include <type_traits>  // for std::conditional and std::enable_if
#include <limits>  // for numeric_limits<>

// Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10
#if CGAL_VERSION_NR < 1041101000
# error Alpha_complex_3d is only available for CGAL >= 4.11
#endif

#if !EIGEN_VERSION_AT_LEAST(3,1,0)
# error Alpha_complex_3d is only available for Eigen3 >= 3.1.0 installed with CGAL
#endif

namespace Gudhi {

namespace alpha_complex {

#ifdef GUDHI_CAN_USE_CXX11_THREAD_LOCAL
thread_local
#endif  // GUDHI_CAN_USE_CXX11_THREAD_LOCAL
    double RELATIVE_PRECISION_OF_TO_DOUBLE = 0.00001;

// Value_from_iterator returns the filtration value from an iterator on alpha shapes values
//
// FAST                         SAFE                         EXACT
// CGAL::to_double(*iterator)   CGAL::to_double(*iterator)   CGAL::to_double(iterator->exact())

template <complexity Complexity>
struct Value_from_iterator {
  template <typename Iterator>
  static double perform(Iterator it) {
    // Default behaviour
    return CGAL::to_double(*it);
  }
};

template <>
struct Value_from_iterator<complexity::EXACT> {
  template <typename Iterator>
  static double perform(Iterator it) {
    return CGAL::to_double(it->exact());
  }
};

/**
 * \class Alpha_complex_3d
 * \brief Alpha complex data structure for 3d specific case.
 *
 * \ingroup alpha_complex
 *
 * \details
 * The data structure is constructing a <a href="https://doc.cgal.org/latest/Alpha_shapes_3/index.html">CGAL 3D Alpha
 * Shapes</a> from a range of points (can be read from an OFF file, cf. Points_off_reader).
 * Duplicate points are inserted once in the Alpha_complex.
 *
 * \tparam Complexity shall be `Gudhi::alpha_complex::complexity` type. Default value is
 * `Gudhi::alpha_complex::complexity::SAFE`.
 *
 * \tparam Weighted Boolean used to set/unset the weighted version of Alpha_complex_3d. Default value is false.
 *
 * \tparam Periodic Boolean used to set/unset the periodic version of Alpha_complex_3d. Default value is false.
 *
 * For the weighted version, weights values are explained on CGAL
 * <a href="https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title0">Alpha shapes 3d</a> and
 * <a href="https://doc.cgal.org/latest/Triangulation_3/index.html#Triangulation3secclassRegulartriangulation">Regular
 * triangulation</a> documentation.
 *
 * For the periodic version, refer to the
 * <a href="https://doc.cgal.org/latest/Periodic_3_triangulation_3/index.html">CGAL’s 3D Periodic Triangulations User
 * Manual </a> for more details.
 * The periodicity is defined by an iso-oriented cuboid with diagonal opposite vertices (x_min, y_min, z_min) and
 * (x_max, y_max, z_max).
 *
 * Please refer to \ref alpha_complex for examples.
 *
 * \remark When Alpha_complex_3d is constructed with an infinite value of alpha (default value), the complex is a
 * 3d Delaunay complex.
 *
 */
template <complexity Complexity = complexity::SAFE, bool Weighted = false, bool Periodic = false>
class Alpha_complex_3d {
  // Epick = Exact_predicates_inexact_constructions_kernel
  // Epeck = Exact_predicates_exact_constructions_kernel
  // Exact_alpha_comparison_tag = exact version of CGAL Alpha_shape_3 and of its objects (Alpha_shape_vertex_base_3 and
  //                           Alpha_shape_cell_base_3). Not available if weighted or periodic.
  //                           Can be CGAL::Tag_false or CGAL::Tag_true. Default is False.
  //                           cf. https://doc.cgal.org/latest/Alpha_shapes_3/classCGAL_1_1Alpha__shape__3.html
  //
  // We could use Epick + CGAL::Tag_true for not weighted nor periodic, but during benchmark, we found a bug
  // https://github.com/CGAL/cgal/issues/3460
  // This is the reason we only use Epick + CGAL::Tag_false, or Epeck
  //
  // FAST                         SAFE                         EXACT
  // Epick + CGAL::Tag_false      Epeck                        Epeck
  using Predicates = typename std::conditional<(Complexity == complexity::FAST),
                                               CGAL::Exact_predicates_inexact_constructions_kernel,
                                               CGAL::Exact_predicates_exact_constructions_kernel>::type;

  // The other way to do a conditional type. Here there are 3 possibilities
  template <typename Predicates, bool Weighted_version, bool Periodic_version>
  struct Kernel_3 {};

  template <typename Predicates, bool Is_periodic>
  struct Kernel_3<Predicates, Is_periodic, false> {
    using Kernel = Predicates;
  };

  template <typename Predicates>
  struct Kernel_3<Predicates, false, true> {
    using Kernel = CGAL::Periodic_3_Delaunay_triangulation_traits_3<Predicates>;
  };
  template <typename Predicates>
  struct Kernel_3<Predicates, true, true> {
    using Kernel = CGAL::Periodic_3_regular_triangulation_traits_3<Predicates>;
  };

  using Kernel = typename Kernel_3<Predicates, Weighted, Periodic>::Kernel;

  using TdsVb = typename std::conditional<Periodic, CGAL::Periodic_3_triangulation_ds_vertex_base_3<>,
                                          CGAL::Triangulation_ds_vertex_base_3<>>::type;

  using Tvb = typename std::conditional<Weighted, CGAL::Regular_triangulation_vertex_base_3<Kernel, TdsVb>,
                                        CGAL::Triangulation_vertex_base_3<Kernel, TdsVb>>::type;

  using Vb = CGAL::Alpha_shape_vertex_base_3<Kernel, Tvb>;

  using TdsCb = typename std::conditional<Periodic, CGAL::Periodic_3_triangulation_ds_cell_base_3<>,
                                          CGAL::Triangulation_ds_cell_base_3<>>::type;

  using Tcb = typename std::conditional<Weighted, CGAL::Regular_triangulation_cell_base_3<Kernel, TdsCb>,
                                        CGAL::Triangulation_cell_base_3<Kernel, TdsCb>>::type;

  using Cb = CGAL::Alpha_shape_cell_base_3<Kernel, Tcb>;
  using Tds = CGAL::Triangulation_data_structure_3<Vb, Cb>;

  // The other way to do a conditional type. Here there 4 possibilities, cannot use std::conditional
  template <typename Kernel, typename Tds, bool Weighted_version, bool Periodic_version>
  struct Triangulation_3 {};

  template <typename Kernel, typename Tds>
  struct Triangulation_3<Kernel, Tds, false, false> {
    using Dt = CGAL::Delaunay_triangulation_3<Kernel, Tds>;
    using Weighted_point_3 = void;
  };
  template <typename Kernel, typename Tds>
  struct Triangulation_3<Kernel, Tds, true, false> {
    using Dt = CGAL::Regular_triangulation_3<Kernel, Tds>;
    using Weighted_point_3 = typename Dt::Weighted_point;
  };
  template <typename Kernel, typename Tds>
  struct Triangulation_3<Kernel, Tds, false, true> {
    using Dt = CGAL::Periodic_3_Delaunay_triangulation_3<Kernel, Tds>;
    using Weighted_point_3 = void;
  };
  template <typename Kernel, typename Tds>
  struct Triangulation_3<Kernel, Tds, true, true> {
    using Dt = CGAL::Periodic_3_regular_triangulation_3<Kernel, Tds>;
    using Weighted_point_3 = typename Dt::Weighted_point;
  };

  /** \brief Is either Delaunay_triangulation_3 (Weighted = false and Periodic = false),
   * Regular_triangulation_3 (Weighted = true and Periodic = false),
   * Periodic_3_Delaunay_triangulation_3 (Weighted = false and Periodic = true)
   * or Periodic_3_regular_triangulation_3 (Weighted = true and Periodic = true).
   *
   * This type is required by `Gudhi::alpha_complex::Alpha_complex_3d::Alpha_shape_3`.
   * */
  using Dt = typename Triangulation_3<Kernel, Tds, Weighted, Periodic>::Dt;

 public:
  /** \brief The <a href="https://doc.cgal.org/latest/Alpha_shapes_3/classCGAL_1_1Alpha__shape__3.html">CGAL 3D Alpha
   * Shapes</a> type.
   *
   *  The `Gudhi::alpha_complex::Alpha_complex_3d` is a wrapper on top of this class to ease the standard, weighted
   *  and/or periodic build of the Alpha complex 3d.*/
  using Alpha_shape_3 = CGAL::Alpha_shape_3<Dt>;

  /** \brief The alpha values type.
   * Must be compatible with double. */
  using FT = typename Alpha_shape_3::FT;

  /** \brief Gives public access to the Bare_point_3 (bare aka. unweighed) type.
   * Here is a Bare_point_3 constructor example:
\code{.cpp}
using Alpha_complex_3d = Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, false, false>;

// x0 = 1., y0 = -1.1, z0 = -1..
Alpha_complex_3d::Bare_point_3 p0(1., -1.1, -1.);
\endcode
   * */
  using Bare_point_3 = typename Kernel::Point_3;

  /** \brief Gives public access to the Weighted_point_3 type. A Weighted point can be constructed as follows:
\code{.cpp}
using Weighted_alpha_complex_3d = Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, true, false>;

// x0 = 1., y0 = -1.1, z0 = -1., weight = 4.
Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_point_3(1., -1.1, -1.), 4.);
\endcode
 *
 * Note: This type is defined to void if Alpha complex is not weighted.
 *
 * */
  using Weighted_point_3 = typename Triangulation_3<Kernel, Tds, Weighted, Periodic>::Weighted_point_3;

  /** \brief `Alpha_complex_3d::Point_3` type is either a `Alpha_complex_3d::Bare_point_3` (Weighted = false) or a
   * `Alpha_complex_3d::Weighted_point_3` (Weighted = true).
   */
  using Point_3 = typename Alpha_shape_3::Point;

 private:
  using Dispatch =
      CGAL::Dispatch_output_iterator<CGAL::cpp11::tuple<CGAL::Object, FT>,
                                     CGAL::cpp11::tuple<std::back_insert_iterator<std::vector<CGAL::Object>>,
                                                        std::back_insert_iterator<std::vector<FT>>>>;

  using Cell_handle = typename Alpha_shape_3::Cell_handle;
  using Facet = typename Alpha_shape_3::Facet;
  using Edge = typename Alpha_shape_3::Edge;
  using Alpha_vertex_handle = typename Alpha_shape_3::Vertex_handle;
  using Vertex_list = boost::container::static_vector<Alpha_vertex_handle, 4>;

 public:
  /** \brief Alpha_complex constructor from a list of points.
   *
   * @param[in] points Range of points to triangulate. Points must be in `Alpha_complex_3d::Point_3`.
   *
   * @pre Available if Alpha_complex_3d is not Periodic.
   *
   * The type InputPointRange must be a range for which std::begin and std::end return input iterators on a
   * `Alpha_complex_3d::Point_3`.
   */
  template <typename InputPointRange>
  Alpha_complex_3d(const InputPointRange& points) {
    static_assert(!Periodic, "This constructor is not available for periodic versions of Alpha_complex_3d");

    alpha_shape_3_ptr_ = std::unique_ptr<Alpha_shape_3>(
        new Alpha_shape_3(std::begin(points), std::end(points), 0, Alpha_shape_3::GENERAL));
  }

  /** \brief Alpha_complex constructor from a list of points and associated weights.
   *
   * @exception std::invalid_argument In debug mode, if points and weights do not have the same size.
   *
   * @param[in] points Range of points to triangulate. Points must be in `Alpha_complex_3d::Bare_point_3`.
   * @param[in] weights Range of weights on points. Weights shall be in double.
   *
   * @pre Available if Alpha_complex_3d is Weighted and not Periodic.
   *
   * The type InputPointRange must be a range for which std::begin and
   * std::end return input iterators on a `Alpha_complex_3d::Bare_point_3`.
   * The type WeightRange must be a range for which std::begin and
   * std::end return an input iterator on a double.
   */
  template <typename InputPointRange, typename WeightRange>
  Alpha_complex_3d(const InputPointRange& points, WeightRange weights) {
    static_assert(Weighted, "This constructor is not available for non-weighted versions of Alpha_complex_3d");
    static_assert(!Periodic, "This constructor is not available for periodic versions of Alpha_complex_3d");
    GUDHI_CHECK((weights.size() == points.size()),
                std::invalid_argument("Points number in range different from weights range number"));

    std::vector<Weighted_point_3> weighted_points_3;

    std::size_t index = 0;
    weighted_points_3.reserve(points.size());
    while ((index < weights.size()) && (index < points.size())) {
      weighted_points_3.push_back(Weighted_point_3(points[index], weights[index]));
      index++;
    }

    alpha_shape_3_ptr_ = std::unique_ptr<Alpha_shape_3>(
        new Alpha_shape_3(std::begin(weighted_points_3), std::end(weighted_points_3), 0, Alpha_shape_3::GENERAL));
  }

  /** \brief Alpha_complex constructor from a list of points and an iso-cuboid coordinates.
   *
   * @exception std::invalid_argument In debug mode, if the size of the cuboid in every directions is not the same.
   *
   * @param[in] points Range of points to triangulate. Points must be in `Alpha_complex_3d::Point_3`.
   * @param[in] x_min Iso-oriented cuboid x_min.
   * @param[in] y_min Iso-oriented cuboid y_min.
   * @param[in] z_min Iso-oriented cuboid z_min.
   * @param[in] x_max Iso-oriented cuboid x_max.
   * @param[in] y_max Iso-oriented cuboid y_max.
   * @param[in] z_max Iso-oriented cuboid z_max.
   *
   * @pre Available if Alpha_complex_3d is Periodic.
   *
   * The type InputPointRange must be a range for which std::begin and std::end return input iterators on a
   * `Alpha_complex_3d::Point_3`.
   *
   * @note In weighted version, please check weights are greater than zero, and lower than 1/64*cuboid length
   * squared.
   */
  template <typename InputPointRange>
  Alpha_complex_3d(const InputPointRange& points, FT x_min, FT y_min,
                   FT z_min, FT x_max, FT y_max, FT z_max) {
    static_assert(Periodic, "This constructor is not available for non-periodic versions of Alpha_complex_3d");
    // Checking if the cuboid is the same in x,y and z direction. If not, CGAL will not process it.
    GUDHI_CHECK(
        (x_max - x_min == y_max - y_min) && (x_max - x_min == z_max - z_min) && (z_max - z_min == y_max - y_min),
        std::invalid_argument("The size of the cuboid in every directions is not the same."));

    // Define the periodic cube
    Dt pdt(typename Kernel::Iso_cuboid_3(x_min, y_min, z_min, x_max, y_max, z_max));
    // Heuristic for inserting large point sets (if pts is reasonably large)
    pdt.insert(std::begin(points), std::end(points), true);
    // As pdt won't be modified anymore switch to 1-sheeted cover if possible
    if (!pdt.is_triangulation_in_1_sheet()) {
      throw std::invalid_argument("Unable to construct a triangulation within a single periodic domain.");
    }
    pdt.convert_to_1_sheeted_covering();

    // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode. This is the default mode
    // Maybe need to set it to GENERAL mode
    alpha_shape_3_ptr_ = std::unique_ptr<Alpha_shape_3>(new Alpha_shape_3(pdt, 0, Alpha_shape_3::GENERAL));
  }

  /** \brief Alpha_complex constructor from a list of points, associated weights and an iso-cuboid coordinates.
   *
   * @exception std::invalid_argument In debug mode, if points and weights do not have the same size.
   * @exception std::invalid_argument In debug mode, if the size of the cuboid in every directions is not the same.
   * @exception std::invalid_argument In debug mode, if a weight is negative, zero, or greater than 1/64*cuboid length
   * squared.
   *
   * @param[in] points Range of points to triangulate. Points must be in `Alpha_complex_3d::Bare_point_3`.
   * @param[in] weights Range of weights on points. Weights shall be in double.
   * @param[in] x_min Iso-oriented cuboid x_min.
   * @param[in] y_min Iso-oriented cuboid y_min.
   * @param[in] z_min Iso-oriented cuboid z_min.
   * @param[in] x_max Iso-oriented cuboid x_max.
   * @param[in] y_max Iso-oriented cuboid y_max.
   * @param[in] z_max Iso-oriented cuboid z_max.
   *
   * @pre Available if Alpha_complex_3d is Weighted and Periodic.
   *
   * The type InputPointRange must be a range for which std::begin and
   * std::end return input iterators on a `Alpha_complex_3d::Bare_point_3`.
   * The type WeightRange must be a range for which std::begin and
   * std::end return an input iterator on a double.
   * The type of x_min, y_min, z_min, x_max, y_max and z_max must be a double.
   */
  template <typename InputPointRange, typename WeightRange>
  Alpha_complex_3d(const InputPointRange& points, WeightRange weights, FT x_min, FT y_min,
                   FT z_min, FT x_max, FT y_max, FT z_max) {
    static_assert(Weighted, "This constructor is not available for non-weighted versions of Alpha_complex_3d");
    static_assert(Periodic, "This constructor is not available for non-periodic versions of Alpha_complex_3d");
    GUDHI_CHECK((weights.size() == points.size()),
                std::invalid_argument("Points number in range different from weights range number"));
    // Checking if the cuboid is the same in x,y and z direction. If not, CGAL will not process it.
    GUDHI_CHECK(
        (x_max - x_min == y_max - y_min) && (x_max - x_min == z_max - z_min) && (z_max - z_min == y_max - y_min),
        std::invalid_argument("The size of the cuboid in every directions is not the same."));

    std::vector<Weighted_point_3> weighted_points_3;

    std::size_t index = 0;
    weighted_points_3.reserve(points.size());

#ifdef GUDHI_DEBUG
    // Defined in GUDHI_DEBUG to avoid unused variable warning for GUDHI_CHECK
    FT maximal_possible_weight = 0.015625 * (x_max - x_min) * (x_max - x_min);
#endif

    while ((index < weights.size()) && (index < points.size())) {
      GUDHI_CHECK((weights[index] < maximal_possible_weight) && (weights[index] >= 0),
                  std::invalid_argument("Invalid weight at index " + std::to_string(index + 1) +
                                        ". Must be positive and less than maximal possible weight = 1/64*cuboid length "
                                        "squared, which is not an acceptable input."));
      weighted_points_3.push_back(Weighted_point_3(points[index], weights[index]));
      index++;
    }

    // Define the periodic cube
    Dt pdt(typename Kernel::Iso_cuboid_3(x_min, y_min, z_min, x_max, y_max, z_max));
    // Heuristic for inserting large point sets (if pts is reasonably large)
    pdt.insert(std::begin(weighted_points_3), std::end(weighted_points_3), true);
    // As pdt won't be modified anymore switch to 1-sheeted cover if possible
    if (!pdt.is_triangulation_in_1_sheet()) {
      throw std::invalid_argument("Unable to construct a triangulation within a single periodic domain.");
    }
    pdt.convert_to_1_sheeted_covering();

    // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode. This is the default mode
    // Maybe need to set it to GENERAL mode
    alpha_shape_3_ptr_ = std::unique_ptr<Alpha_shape_3>(new Alpha_shape_3(pdt, 0, Alpha_shape_3::GENERAL));
  }

  /** \brief Inserts all Delaunay triangulation into the simplicial complex.
   * It also computes the filtration values accordingly to the \ref createcomplexalgorithm
   *
   * \tparam SimplicialComplexForAlpha3d must meet `SimplicialComplexForAlpha3d` concept.
   *
   * @param[in] complex SimplicialComplexForAlpha3d to be created.
   * @param[in] max_alpha_square maximum for alpha square value. Default value is +\f$\infty\f$, and there is very
   * little point using anything else since it does not save time.
   *
   * @return true if creation succeeds, false otherwise.
   *
   * @pre The simplicial complex must be empty (no vertices).
   *
   */
  template <typename SimplicialComplexForAlpha3d,
            typename Filtration_value = typename SimplicialComplexForAlpha3d::Filtration_value>
  bool create_complex(SimplicialComplexForAlpha3d& complex,
                      Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity()) {
    if (complex.num_vertices() > 0) {
      std::cerr << "Alpha_complex_3d create_complex - complex is not empty\n";
      return false;  // ----- >>
    }

    using Complex_vertex_handle = typename SimplicialComplexForAlpha3d::Vertex_handle;
    using Simplex_tree_vector_vertex = std::vector<Complex_vertex_handle>;

#ifdef DEBUG_TRACES
    std::size_t count_vertices = 0;
    std::size_t count_edges = 0;
    std::size_t count_facets = 0;
    std::size_t count_cells = 0;
#endif  // DEBUG_TRACES
    std::vector<CGAL::Object> objects;
    std::vector<FT> alpha_values;

    Dispatch dispatcher = CGAL::dispatch_output<CGAL::Object, FT>(std::back_inserter(objects),
                                                                                std::back_inserter(alpha_values));

    alpha_shape_3_ptr_->filtration_with_alpha_values(dispatcher);
#ifdef DEBUG_TRACES
    std::cout << "filtration_with_alpha_values returns : " << objects.size() << " objects" << std::endl;
#endif  // DEBUG_TRACES

    using Alpha_value_iterator = typename std::vector<FT>::const_iterator;
    Alpha_value_iterator alpha_value_iterator = alpha_values.begin();
    for (auto object_iterator : objects) {
      Vertex_list vertex_list;

      // Retrieve Alpha shape vertex list from object
      if (const Cell_handle* cell = CGAL::object_cast<Cell_handle>(&object_iterator)) {
        for (auto i = 0; i < 4; i++) {
#ifdef DEBUG_TRACES
          std::cout << "from cell[" << i << "] - Point coordinates (" << (*cell)->vertex(i)->point() << ")"
                    << std::endl;
#endif  // DEBUG_TRACES
          vertex_list.push_back((*cell)->vertex(i));
        }
#ifdef DEBUG_TRACES
        count_cells++;
#endif  // DEBUG_TRACES
      } else if (const Facet* facet = CGAL::object_cast<Facet>(&object_iterator)) {
        for (auto i = 0; i < 4; i++) {
          if ((*facet).second != i) {
#ifdef DEBUG_TRACES
            std::cout << "from facet=[" << i << "] - Point coordinates (" << (*facet).first->vertex(i)->point() << ")"
                      << std::endl;
#endif  // DEBUG_TRACES
            vertex_list.push_back((*facet).first->vertex(i));
          }
        }
#ifdef DEBUG_TRACES
        count_facets++;
#endif  // DEBUG_TRACES
      } else if (const Edge* edge = CGAL::object_cast<Edge>(&object_iterator)) {
        for (auto i : {(*edge).second, (*edge).third}) {
#ifdef DEBUG_TRACES
          std::cout << "from edge[" << i << "] - Point coordinates (" << (*edge).first->vertex(i)->point() << ")"
                    << std::endl;
#endif  // DEBUG_TRACES
          vertex_list.push_back((*edge).first->vertex(i));
        }
#ifdef DEBUG_TRACES
        count_edges++;
#endif  // DEBUG_TRACES
      } else if (const Alpha_vertex_handle* vertex = CGAL::object_cast<Alpha_vertex_handle>(&object_iterator)) {
#ifdef DEBUG_TRACES
        count_vertices++;
        std::cout << "from vertex - Point coordinates (" << (*vertex)->point() << ")" << std::endl;
#endif  // DEBUG_TRACES
        vertex_list.push_back((*vertex));
      }
      // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex
      Simplex_tree_vector_vertex the_simplex;
      for (auto the_alpha_shape_vertex : vertex_list) {
        auto the_map_iterator = map_cgal_simplex_tree.find(the_alpha_shape_vertex);
        if (the_map_iterator == map_cgal_simplex_tree.end()) {
          // alpha shape not found
          Complex_vertex_handle vertex = map_cgal_simplex_tree.size();
#ifdef DEBUG_TRACES
          std::cout << "Point (" << the_alpha_shape_vertex->point() << ") not found - insert new vertex id " << vertex
                    << std::endl;
#endif  // DEBUG_TRACES
          the_simplex.push_back(vertex);
          map_cgal_simplex_tree.emplace(the_alpha_shape_vertex, vertex);
        } else {
          // alpha shape found
          Complex_vertex_handle vertex = the_map_iterator->second;
#ifdef DEBUG_TRACES
          std::cout << "Point (" << the_alpha_shape_vertex->point() << ") found as vertex id " << vertex << std::endl;
#endif  // DEBUG_TRACES
          the_simplex.push_back(vertex);
        }
      }
      // Construction of the simplex_tree
      Filtration_value filtr = Value_from_iterator<Complexity>::perform(alpha_value_iterator);

#ifdef DEBUG_TRACES
      std::cout << "filtration = " << filtr << std::endl;
#endif  // DEBUG_TRACES
      complex.insert_simplex(the_simplex, static_cast<Filtration_value>(filtr));
      GUDHI_CHECK(alpha_value_iterator != alpha_values.end(), "CGAL provided more simplices than values");
      ++alpha_value_iterator;
    }

#ifdef DEBUG_TRACES
    std::cout << "vertices \t" << count_vertices << std::endl;
    std::cout << "edges \t\t" << count_edges << std::endl;
    std::cout << "facets \t\t" << count_facets << std::endl;
    std::cout << "cells \t\t" << count_cells << std::endl;
#endif  // DEBUG_TRACES
    // --------------------------------------------------------------------------------------------
    // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension
    complex.make_filtration_non_decreasing();
    // Remove all simplices that have a filtration value greater than max_alpha_square
    complex.prune_above_filtration(max_alpha_square);
    // --------------------------------------------------------------------------------------------
    return true;
  }

  /** \brief get_point returns the point corresponding to the vertex given as parameter.
   *
   * @param[in] vertex Vertex handle of the point to retrieve.
   * @return The point found.
   * @exception std::out_of_range In case vertex is not found (cf. std::vector::at).
   */
  const Point_3& get_point(std::size_t vertex) {
    if (map_cgal_simplex_tree.size() != cgal_vertex_iterator_vector.size()) {
      cgal_vertex_iterator_vector.resize(map_cgal_simplex_tree.size());
      for (auto map_iterator : map_cgal_simplex_tree) {
        cgal_vertex_iterator_vector[map_iterator.second] = map_iterator.first;
      }

    }
    auto cgal_vertex_iterator = cgal_vertex_iterator_vector.at(vertex);
    return cgal_vertex_iterator->point();
  }

 private:
  // use of a unique_ptr on cgal Alpha_shape_3, as copy and default constructor is not available - no need to be freed
  std::unique_ptr<Alpha_shape_3> alpha_shape_3_ptr_;

  // Map type to switch from CGAL vertex iterator to simplex tree vertex handle.
  std::unordered_map<Alpha_vertex_handle, std::size_t> map_cgal_simplex_tree;
  // Vector type to switch from simplex tree vertex handle to CGAL vertex iterator.
  std::vector<Alpha_vertex_handle> cgal_vertex_iterator_vector;
};

}  // namespace alpha_complex

}  // namespace Gudhi

#endif  // ALPHA_COMPLEX_3D_H_