summaryrefslogtreecommitdiff
path: root/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h
blob: 21707e5e5eb57c07f8fc68d8044341bae195de2b (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
/*    This file is part of the Gudhi Library. The Gudhi library
 *    (Geometric Understanding in Higher Dimensions) is a generic C++
 *    library for computational topology.
 *
 *    Author(s):       Vincent Rouvreau
 *
 *    Copyright (C) 2018 Inria
 *
 *    This program is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 3 of the License, or
 *    (at your option) any later version.
 *
 *    This program is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#ifndef ALPHA_COMPLEX_3D_H_
#define ALPHA_COMPLEX_3D_H_


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

#include <CGAL/Object.h>
#include <CGAL/tuple.h>
#include <CGAL/iterator.h>


namespace Gudhi {

namespace alpha_complex {

template<typename AlphaComplex3dOptions>
class Alpha_complex_3d {
private:
  using Alpha_shape_3 = typename AlphaComplex3dOptions::Alpha_shape_3;
  // filtration with alpha values needed type definition
  using Alpha_value_type = typename Alpha_shape_3::FT;
  using Object = CGAL::Object;
  using Dispatch = CGAL::Dispatch_output_iterator<CGAL::cpp11::tuple<Object, Alpha_value_type>,
      CGAL::cpp11::tuple<std::back_insert_iterator<std::vector<Object> >,
          std::back_insert_iterator<std::vector<Alpha_value_type> > > >;
  using Cell_handle = typename Alpha_shape_3::Cell_handle;
  using Facet = typename Alpha_shape_3::Facet;
  using Edge_3 = typename Alpha_shape_3::Edge;
  using Vertex_handle = typename Alpha_shape_3::Vertex_handle;

public:
  /** \brief Alpha_complex constructor from a list of points.
  *
  * Duplicate points are inserted once in the Alpha_complex. This is the reason why the vertices may be not contiguous.
  *
  * @param[in] points Range of points to triangulate. Points must be in AlphaComplex3dOptions::Point_3
  *
  * The type InputPointRange must be a range for which std::begin and
  * std::end return input iterators on a AlphaComplex3dOptions::Point_3.
  */
  template<typename InputPointRange >
  Alpha_complex_3d(const InputPointRange& points) {
    static_assert(!AlphaComplex3dOptions::weighted, "This constructor is not available for weighted versions of Alpha_complex_3d");
    static_assert(!AlphaComplex3dOptions::periodic, "This constructor is not available for periodic versions of Alpha_complex_3d");
    std::cout << points[0] << std::endl;
    Alpha_shape_3 alpha_shape_3(std::begin(points), std::end(points), 0, Alpha_shape_3::GENERAL);

    Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects),
                                                                    std::back_inserter(the_alpha_values));

    alpha_shape_3.filtration_with_alpha_values(disp);
#ifdef DEBUG_TRACES
    std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl;
#endif  // DEBUG_TRACES

  }

  /** \brief Alpha_complex constructor from a list of points.
*
* Duplicate points are inserted once in the Alpha_complex. This is the reason why the vertices may be not contiguous.
*
* @param[in] points Range of points to triangulate. Points must be in Kernel::Point_d
*
* The type InputPointRange must be a range for which std::begin and
* std::end return input iterators on a Kernel::Point_d.
*/
  template<typename InputPointRange , typename WeightRange>
  Alpha_complex_3d(const InputPointRange& points, WeightRange weights) {
    static_assert(AlphaComplex3dOptions::weighted,
                  "This constructor is not available for non-weighted versions of Alpha_complex_3d");
    static_assert(!AlphaComplex3dOptions::periodic,
                  "This constructor is not available for periodic versions of Alpha_complex_3d");
    GUDHI_CHECK((weights.size() == points.size()),
                std::invalid_argument("Alpha_complex_3d constructor with weights requires points number to be equal with points number"));

    using Weighted_point_3 = typename AlphaComplex3dOptions::Weighted_point_3;
    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 alpha_shape_3(std::begin(weighted_points_3), std::end(weighted_points_3), 0, Alpha_shape_3::GENERAL);

    Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects),
                                                                    std::back_inserter(the_alpha_values));

    alpha_shape_3.filtration_with_alpha_values(disp);
#ifdef DEBUG_TRACES
    std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl;
#endif  // DEBUG_TRACES
  }

  /** \brief Alpha_complex constructor from a list of points.
*
* Duplicate points are inserted once in the Alpha_complex. This is the reason why the vertices may be not contiguous.
*
* @param[in] points Range of points to triangulate. Points must be in Kernel::Point_d
*
* The type InputPointRange must be a range for which std::begin and
* std::end return input iterators on a Kernel::Point_d.
*/
  template<typename InputPointRange>
  Alpha_complex_3d(const InputPointRange& points,
                   Alpha_value_type x_min, Alpha_value_type y_min, Alpha_value_type z_min,
                   Alpha_value_type x_max, Alpha_value_type y_max, Alpha_value_type z_max) {
    static_assert(!AlphaComplex3dOptions::weighted,
                  "This constructor is not available for weighted versions of Alpha_complex_3d");
    static_assert(AlphaComplex3dOptions::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."));

    using Periodic_delaunay_triangulation_3 = typename AlphaComplex3dOptions::Periodic_delaunay_triangulation_3;
    using Iso_cuboid_3 = typename AlphaComplex3dOptions::Iso_cuboid_3;
    // Define the periodic cube
    Periodic_delaunay_triangulation_3 pdt(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
    GUDHI_CHECK(pdt.is_triangulation_in_1_sheet(),
                std::invalid_argument("Uable to construct a triangulation within a single periodic domain."));

    // 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 as(pdt, 0, Alpha_shape_3::GENERAL);

    // filtration with alpha values from alpha shape
    std::vector<Object> the_objects;
    std::vector<Alpha_value_type> the_alpha_values;

    Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects),
                                                                    std::back_inserter(the_alpha_values));

    as.filtration_with_alpha_values(disp);
#ifdef DEBUG_TRACES
    std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl;
#endif  // DEBUG_TRACES

  }

private:
  std::vector<Object> the_objects;
  std::vector<Alpha_value_type> the_alpha_values;

};

}  // namespace alpha_complex

}  // namespace Gudhi

#endif  // ALPHA_COMPLEX_3D_H_