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
|
#include "triangle2d.hpp"
#include "distance.hpp"
#include "robust_orientation.hpp"
#include "segment2d.hpp"
#include "base/math.hpp"
#include "std/chrono.hpp"
#include "std/random.hpp"
using namespace m2::robust;
namespace m2
{
bool IsPointInsideTriangle(m2::PointD const & pt, m2::PointD const & p1,
m2::PointD const & p2, m2::PointD const & p3)
{
double const s1 = OrientedS(p1, p2, pt);
double const s2 = OrientedS(p2, p3, pt);
double const s3 = OrientedS(p3, p1, pt);
// In the case of degenerate triangles we need to check that pt lies
// on (p1, p2), (p2, p3) or (p3, p1).
if (s1 == 0.0 && s2 == 0.0 && s3 == 0.0)
{
return IsPointOnSegment(pt, p1, p2) || IsPointOnSegment(pt, p2, p3) ||
IsPointOnSegment(pt, p3, p1);
}
return ((s1 >= 0.0 && s2 >= 0.0 && s3 >= 0.0) ||
(s1 <= 0.0 && s2 <= 0.0 && s3 <= 0.0));
}
bool IsPointStrictlyInsideTriangle(m2::PointD const & pt, m2::PointD const & p1,
m2::PointD const & p2, m2::PointD const & p3)
{
double const s1 = OrientedS(p1, p2, pt);
double const s2 = OrientedS(p2, p3, pt);
double const s3 = OrientedS(p3, p1, pt);
return ((s1 > 0.0 && s2 > 0.0 && s3 > 0.0) ||
(s1 < 0.0 && s2 < 0.0 && s3 < 0.0));
}
bool IsPointInsideTriangles(m2::PointD const & pt, vector<m2::TriangleD> const & v)
{
for (auto const & triangle : v)
{
if (IsPointInsideTriangle(pt, triangle.p1(), triangle.p2(), triangle.p3()))
return true;
}
return false;
}
m2::PointD GetRandomPointInsideTriangle(m2::TriangleD const & t)
{
size_t kDistribMax = 1000;
default_random_engine engine(system_clock::now().time_since_epoch().count());
uniform_int_distribution<> distrib(0, kDistribMax);
double const r1 = sqrt(static_cast<double>(distrib(engine)) / kDistribMax);
double const r2 = static_cast<double>(distrib(engine)) / kDistribMax;
return t.m_points[0] * (1.0 - r1) + t.m_points[1] * r1 * (1.0 - r2) + t.m_points[2] * r2 * r1;
}
m2::PointD GetRandomPointInsideTriangles(vector<m2::TriangleD> const & v)
{
if (v.empty())
return m2::PointD();
default_random_engine engine(system_clock::now().time_since_epoch().count());
uniform_int_distribution<> distrib(0, v.size() - 1);
return GetRandomPointInsideTriangle(v[distrib(engine)]);
}
m2::PointD ProjectPointToTriangles(m2::PointD const & pt, vector<m2::TriangleD> const & v)
{
if (v.empty())
return pt;
auto distToLine = m2::DistanceToLineSquare<m2::PointD>();
int minT = -1;
int minI = -1;
double minDist = numeric_limits<double>::max();
for (int t = 0; t < static_cast<int>(v.size()); t++)
{
for (int i = 0; i < 3; i++)
{
distToLine.SetBounds(v[t].m_points[i], v[t].m_points[(i + 1) % 3]);
double const dist = distToLine(pt);
if (dist < minDist)
{
minDist = dist;
minT = t;
minI = i;
}
}
}
auto projectToLine = m2::ProjectionToSection<m2::PointD>();
projectToLine.SetBounds(v[minT].m_points[minI], v[minT].m_points[(minI + 1) % 3]);
return projectToLine(pt);
}
} // namespace m2;
|
