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