1 // Copyright (c) 2005 INRIA Sophia-Antipolis (France).
2 // All rights reserved.
3 //
4 // This file is part of CGAL (www.cgal.org).
5 //
6 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Stream_lines_2/include/CGAL/Triangular_field_2.h $
7 // $Id: Triangular_field_2.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot
8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
9 //
10 //
11 // Author(s) : Abdelkrim Mebarki <Abdelkrim.Mebarki@sophia.inria.fr>
12
13 #ifndef CGAL_TRIANGULAR_FIELD_2_H_
14 #define CGAL_TRIANGULAR_FIELD_2_H_
15
16 #include <CGAL/license/Stream_lines_2.h>
17
18
19 #include <CGAL/Delaunay_triangulation_2.h>
20 #include <CGAL/Polygon_2_algorithms.h>
21 #include <CGAL/squared_distance_2.h>
22
23 #include <cfloat>
24 #include <cstdlib>
25
26 #include <iostream>
27 #include <fstream>
28 #include <list>
29 #include <cmath>
30 #include <string>
31
32 #include <CGAL/Triangulation_face_base_with_info_2.h>
33
34 namespace CGAL {
35
36 template <class StreamLinesTraits_2>
37 class Triangular_field_2
38 {
39 public:
40 typedef Triangular_field_2<StreamLinesTraits_2> Vector_field_2;
41 typedef StreamLinesTraits_2 Geom_traits;
42 typedef typename StreamLinesTraits_2::FT FT;
43 typedef typename StreamLinesTraits_2::Point_2 Point_2;
44 typedef typename StreamLinesTraits_2::Vector_2 Vector_2;
45 protected:
46 typedef CGAL::Triangulation_vertex_base_2<StreamLinesTraits_2> Vb;
47 typedef CGAL::Triangulation_face_base_with_info_2<Vector_2, StreamLinesTraits_2> Fb;
48 typedef CGAL::Triangulation_data_structure_2<Vb,Fb> TDS;
49 typedef CGAL::Delaunay_triangulation_2<StreamLinesTraits_2,TDS> D_Ttr;
50 typedef typename D_Ttr::Vertex_handle Vertex_handle;
51 typedef typename TDS::Face_handle Face_handle;
52 public:
53 typedef typename D_Ttr::Vertex_iterator Vertex_iterator;
54 D_Ttr m_D_Ttr;
55 protected:
56 Vector_2 get_vector_field(const Point_2 & p) const;
57
58 FT get_density_field(const Point_2 & p) const;
59
60 template <class PointInputIterator, class VectorInputIterator>
fill(PointInputIterator pBegin,PointInputIterator pEnd,VectorInputIterator vBegin)61 void fill(PointInputIterator pBegin, PointInputIterator pEnd, VectorInputIterator vBegin)
62 {
63 std::cout << "reading file...\n";
64 while(pBegin != pEnd)
65 {
66 Point_2 p;
67 Vector_2 v;
68 p = (*pBegin);
69 v = (*vBegin);
70 m_D_Ttr.insert(p);
71 field_map[p] = v;
72 if (m_D_Ttr.number_of_vertices() == 1)
73 {
74 maxx = minx = p.x();
75 maxy = miny = p.y();
76 }
77 if(p.x()<minx)
78 {
79 minx = p.x();
80 }
81 if(p.y()<miny)
82 {
83 miny = p.y();
84 }
85 if(p.x()>maxx)
86 {
87 maxx = p.x();
88 }
89 if(p.y()>maxy)
90 {
91 maxy = p.y();
92 }
93 ++pBegin;
94 ++vBegin;
95 }
96 std::cout << "number of samples " << m_D_Ttr.number_of_vertices() << "\n";
97 }
98 public:
99 template <class PointInputIterator, class VectorInputIterator>
Triangular_field_2(PointInputIterator pBegin,PointInputIterator pEnd,VectorInputIterator vBegin)100 Triangular_field_2(PointInputIterator pBegin, PointInputIterator
101 pEnd, VectorInputIterator vBegin)
102 {
103 fill(pBegin, pEnd, vBegin);
104 }
105
106 inline typename Geom_traits::Iso_rectangle_2 bbox() const;
107
get_field(const Point_2 & p)108 std::pair<Vector_2,FT> get_field(const Point_2 & p) const
109 {
110 CGAL_assertion(is_in_domain(p));
111 Vector_2 v = get_vector_field(p);
112 FT density = get_density_field(p);
113 return std::make_pair(v, density);
114 }
115
116 bool is_in_domain(const Point_2 & p) const;
117
118 FT get_integration_step(const Point_2 &) const;
119
120 FT get_integration_step() const;
121 protected:
122 FT minx;
123 FT miny;
124 FT maxx;
125 FT maxy;
126 protected:
127 mutable std::map<Point_2, Vector_2> field_map;
128
distance(const Point_2 & p,const Point_2 & q)129 FT distance(const Point_2 & p, const Point_2 & q)
130 {
131 return sqrt( CGAL::squared_distance(p, q) );
132 }
133 };
134
135 template <class StreamLinesTraits_2>
136 inline
137 typename Triangular_field_2<StreamLinesTraits_2>::Geom_traits::Iso_rectangle_2
bbox()138 Triangular_field_2<StreamLinesTraits_2>::bbox() const
139 {
140 return typename Geom_traits::Iso_rectangle_2(minx, miny, maxx,
141 maxy);
142 }
143
144 template <class StreamLinesTraits_2>
145 bool
is_in_domain(const Point_2 & p)146 Triangular_field_2<StreamLinesTraits_2>::is_in_domain(const Point_2 &
147 p) const
148 {
149 Face_handle f = m_D_Ttr.locate(p);
150 return !m_D_Ttr.is_infinite(f);
151 }
152
153 template <class StreamLinesTraits_2>
154 typename Triangular_field_2<StreamLinesTraits_2>::Vector_2
get_vector_field(const Point_2 & p)155 Triangular_field_2<StreamLinesTraits_2>::get_vector_field(const
156 Point_2 & p)
157 const
158 {
159 Face_handle m_Face_handle = m_D_Ttr.locate(p);
160 CGAL_assertion(is_in_domain(p));
161 Vertex_handle v0 = m_Face_handle->vertex(0);
162 Vertex_handle v1 = m_Face_handle->vertex(1);
163 Vertex_handle v2 = m_Face_handle->vertex(2);
164 const Point_2 & p0 = v0->point();
165 const Point_2 & p1 = v1->point();
166 const Point_2 & p2 = v2->point();
167 FT s0,s1,s2,s;
168 std::vector<Point_2> vec;
169 vec.push_back(p0); vec.push_back(p1); vec.push_back(p2);
170 s = polygon_area_2(vec.begin(), vec.end(), m_D_Ttr.geom_traits());
171 vec.clear();
172 vec.push_back(p); vec.push_back(p1); vec.push_back(p2);
173 s0 = polygon_area_2(vec.begin(), vec.end(), m_D_Ttr.geom_traits());
174 vec.clear();
175 vec.push_back(p0); vec.push_back(p); vec.push_back(p2);
176 s1 = polygon_area_2(vec.begin(), vec.end(), m_D_Ttr.geom_traits());
177 vec.clear();
178 vec.push_back(p0); vec.push_back(p1); vec.push_back(p);
179 s2 = polygon_area_2(vec.begin(), vec.end(), m_D_Ttr.geom_traits());
180 vec.clear();
181 s0 = s0 / s; s1 = s1 / s; s2 = s2 / s;
182 Vector_2 v_0 = field_map[p0];
183 Vector_2 v_1 = field_map[p1];
184 Vector_2 v_2 = field_map[p2];
185 FT x = ((v_0.x()*s0)+(v_1.x()*s1)+(v_2.x()*s2));
186 FT y = ((v_0.y()*s0)+(v_1.y()*s1)+(v_2.y()*s2));
187 FT normal = sqrt((x)*(x) + (y)*(y));
188 if (normal != 0)
189 {
190 x = x / normal;
191 y = y / normal;
192 }
193 Vector_2 v = Vector_2(x, y);
194 return v;
195 }
196
197 template <class StreamLinesTraits_2>
198 typename Triangular_field_2<StreamLinesTraits_2>::FT
get_density_field(const Point_2 & p)199 Triangular_field_2<StreamLinesTraits_2>::get_density_field(const
200 Point_2 &
201 p) const
202 {
203 return p.x();
204 }
205
206 template<class StreamLinesTraits_2>
207 typename Triangular_field_2<StreamLinesTraits_2>::FT
get_integration_step(const Point_2 &)208 Triangular_field_2<StreamLinesTraits_2>::get_integration_step(const
209 Point_2
210 &) const
211 {
212 return 1.0;
213 }
214
215 template<class StreamLinesTraits_2>
216 typename Triangular_field_2<StreamLinesTraits_2>::FT
get_integration_step()217 Triangular_field_2<StreamLinesTraits_2>::get_integration_step() const
218 {
219 return 1.0;
220 }
221
222 } //namespace CGAL
223
224 #endif
225