Skip to content

Commit e1f4050

Browse files
davidbrownedavidbrowne
committed
more examples of basic geometric algorithms
1 parent 6749e0d commit e1f4050

File tree

4 files changed

+257
-7
lines changed

4 files changed

+257
-7
lines changed

‎README.md‎

Lines changed: 9 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -20,19 +20,21 @@ using namespace dsga;
2020

2121
// get a 2D vector that is perpendicular (rotated 90 degrees counter-clockwise)
2222
// to a 2D vector in the plane
23-
constexpr auto get_perpendicular1(const vec2 &some_vec) noexcept
23+
template <floating_point_scalar T>
24+
constexpr auto get_perpendicular1(const basic_vector<T, 2> &some_vec) noexcept
2425
{
25-
auto cos90 = 0.0f;
26-
auto sin90 = 1.0f;
26+
auto cos90 = 0.0f;
27+
auto sin90 = 1.0f;
2728

28-
// rotation matrix -- components in column major order
29-
return mat2(cos90, sin90, -sin90, cos90) * some_vec;
29+
// rotation matrix -- components in column major order
30+
return basic_matrix<T, 2, 2>(cos90, sin90, -sin90, cos90) * some_vec;
3031
}
3132

3233
// same as above, different implementation
33-
constexpr auto get_perpendicular2(const vec2 &some_vec) noexcept
34+
template <floating_point_scalar T>
35+
constexpr auto get_perpendicular2(const basic_vector<T, 2> &some_vec) noexcept
3436
{
35-
return vec2(-1, 1) * some_vec.yx;
37+
return basic_vector<T, 2>(-1, 1) * some_vec.yx;
3638
}
3739
```
3840

‎VS2022/dsga.vcxproj‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -169,6 +169,7 @@
169169
<ClInclude Include="..\dev_3rd\doctest.h" />
170170
<ClInclude Include="..\dev_3rd\nanobench.h" />
171171
<ClInclude Include="..\examples\angle.hxx" />
172+
<ClInclude Include="..\examples\basic.hxx" />
172173
<ClInclude Include="..\examples\bezier.hxx" />
173174
<ClInclude Include="..\examples\format_output.hxx" />
174175
<ClInclude Include="..\examples\ostream_output.hxx" />

‎VS2022/dsga.vcxproj.filters‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -42,6 +42,9 @@
4242
<ClInclude Include="..\examples\stl.hxx">
4343
<Filter>Header Files</Filter>
4444
</ClInclude>
45+
<ClInclude Include="..\examples\basic.hxx">
46+
<Filter>Header Files</Filter>
47+
</ClInclude>
4548
</ItemGroup>
4649
<ItemGroup>
4750
<None Include="..\README.md" />

‎examples/basic.hxx‎

Lines changed: 244 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,244 @@
1+
2+
// Copyright David Browne 2020-2023.
3+
// Distributed under the Boost Software License, Version 1.0.
4+
// (See accompanying file LICENSE_1_0.txt or copy at
5+
// https://www.boost.org/LICENSE_1_0.txt)
6+
7+
#include "dsga.hxx"
8+
9+
// get a 2D vector that is perpendicular (rotated 90 degrees counter-clockwise)
10+
// to a 2D vector in the plane
11+
template <dsga::floating_point_scalar T>
12+
constexpr auto get_perpendicular1(const dsga::basic_vector<T, 2> &some_vec) noexcept
13+
{
14+
auto cos90 = 0.0f;
15+
auto sin90 = 1.0f;
16+
17+
// rotation matrix -- components in column major order
18+
return dsga::basic_matrix<T, 2, 2>(cos90, sin90, -sin90, cos90) * some_vec;
19+
}
20+
21+
// same as above, different implementation
22+
template <dsga::floating_point_scalar T>
23+
constexpr auto get_perpendicular2(const dsga::basic_vector<T, 2> &some_vec) noexcept
24+
{
25+
return dsga::basic_vector<T, 2>(-1, 1) * some_vec.yx;
26+
}
27+
28+
// if p1 == p2 == p3, then there is a singularity -- we will have 0/0 problem, when real answer should be p1 or p2 or p3.
29+
// the return value c is the center point of a circle inscribed in a triangle represented by the vertices p1, p2, and p3.
30+
// a line segment from c to any of the vertices bisects the angles at the vertices.
31+
template <bool W1, dsga::floating_point_scalar T, typename D1, bool W2, typename D2, bool W3, typename D3>
32+
constexpr dsga::basic_vector<T, 3u> triangle_incenter(const dsga::vector_base<W1, T, 3u, D1> &p1,
33+
const dsga::vector_base<W2, T, 3u, D2> &p2,
34+
const dsga::vector_base<W3, T, 3u, D3> &p3)
35+
{
36+
auto mag1 = dsga::distance(p2, p3);
37+
auto mag2 = dsga::distance(p3, p1);
38+
auto mag3 = dsga::distance(p1, p2);
39+
40+
return (p1 * mag1 + p2 * mag2 + p3 * mag3) / (mag1 + mag2 + mag3);
41+
}
42+
43+
// the return value c is the center point of the biggest sphere inscribed in a tetrahedron represented by the vertices p1,
44+
// p2, p3, and the implicit origin. c is equidistant from the four planes of the triangle faces of the tetrahedron.
45+
template <bool W1, dsga::floating_point_scalar T, typename D1, bool W2, typename D2, bool W3, typename D3>
46+
constexpr dsga::basic_vector<T, 3u> tetrahedron_incenter(const dsga::vector_base<W1, T, 3u, D1> &p1,
47+
const dsga::vector_base<W2, T, 3u, D2> &p2,
48+
const dsga::vector_base<W3, T, 3u, D3> &p3)
49+
{
50+
auto mag1 = dsga::length(dsga::cross_matrix(p2) * p3);
51+
auto mag2 = dsga::length(dsga::cross_matrix(p3) * p1);
52+
auto mag3 = dsga::length(dsga::cross_matrix(p1) * p2);
53+
auto mag4 = dsga::length(dsga::cross_matrix(p2 - p1) * (p3 - p1));
54+
55+
return (p1 * mag1 + p2 * mag2 + p3 * mag3) / (mag1 + mag2 + mag3 + mag4);
56+
}
57+
58+
// find center of circle that goes through the three points
59+
template <bool W1, dsga::floating_point_scalar T, typename D1, bool W2, typename D2, bool W3, typename D3>
60+
constexpr auto three_point_circle_center(const dsga::vector_base<W1, T, 3u, D1> &p1,
61+
const dsga::vector_base<W2, T, 3u, D2> &p2,
62+
const dsga::vector_base<W3, T, 3u, D3> &p3) noexcept
63+
{
64+
auto v = p2 - p1;
65+
auto u = dsga::basic_vector<T, 3u>(p2);
66+
auto w = p3 - p2;
67+
68+
//auto u = p1;
69+
//auto v = p2;
70+
//auto w = p3;
71+
72+
auto cross_term = dsga::cross_matrix(v) * w;
73+
74+
return u + T(0.5) * (dsga::dot(v, v) * dsga::outerProduct(w, w) - dsga::dot(w, w) * dsga::outerProduct(v, v)) * (v + w) / dsga::dot(cross_term, cross_term);
75+
}
76+
77+
// find radius of circle that goes through the three points
78+
template <bool W1, dsga::floating_point_scalar T, typename D1, bool W2, typename D2, bool W3, typename D3>
79+
constexpr auto three_point_circle_radius(const dsga::vector_base<W1, T, 3u, D1> &p1,
80+
const dsga::vector_base<W2, T, 3u, D2> &p2,
81+
const dsga::vector_base<W3, T, 3u, D3> &p3) noexcept
82+
{
83+
auto v = p2 - p1;
84+
[[maybe_unused]] auto u = dsga::basic_vector<T, 3u>(p2);
85+
auto w = p3 - p2;
86+
87+
//auto u = p1;
88+
//auto v = p2;
89+
//auto w = p3;
90+
91+
auto cross_term = dsga::cross_matrix(v) * w;
92+
93+
return T(0.5) * dsga::length(v) * dsga::length(w) * dsga::length(v + w) / dsga::length(cross_term);
94+
}
95+
96+
// gives closest projection point from point to a line made from line segment p1 <=> p2
97+
template <bool W1, dsga::floating_point_scalar T, typename D1, bool W2, typename D2, bool W3, typename D3>
98+
constexpr auto project_to_line1(const dsga::vector_base<W1, T, 3u, D1> &point,
99+
const dsga::vector_base<W2, T, 3u, D2> &p1,
100+
const dsga::vector_base<W3, T, 3u, D3> &p2) noexcept
101+
{
102+
auto hyp = point - p1;
103+
auto v1 = p2 - p1;
104+
auto t = dsga::dot(hyp, v1) / dsga::dot(v1, v1);
105+
106+
return p1 + (t * v1);
107+
}
108+
109+
// same as above, different implementation
110+
template <bool W1, dsga::floating_point_scalar T, typename D1, bool W2, typename D2, bool W3, typename D3>
111+
constexpr auto project_to_line2(const dsga::vector_base<W1, T, 3u, D1> &point,
112+
const dsga::vector_base<W2, T, 3u, D2> &p1,
113+
const dsga::vector_base<W3, T, 3u, D3> &p2) noexcept
114+
{
115+
auto hyp = point - p1;
116+
auto v1 = p2 - p1;
117+
return p1 + dsga::outerProduct(v1, v1) * hyp / dsga::dot(v1, v1);
118+
}
119+
120+
// same as above, different implementation, paying more attention to attenuating roundoff
121+
template <bool W1, dsga::floating_point_scalar T, typename D1, bool W2, typename D2, bool W3, typename D3>
122+
constexpr auto project_to_line3(const dsga::vector_base<W1, T, 3u, D1> &point,
123+
const dsga::vector_base<W2, T, 3u, D2> &p1,
124+
const dsga::vector_base<W3, T, 3u, D3> &p2) noexcept
125+
{
126+
auto hyp1 = point - p1;
127+
auto v1 = p2 - p1;
128+
auto hyp2 = point - p2;
129+
auto v2 = p1 - p2;
130+
131+
auto hyp = hyp1;
132+
auto v = v1;
133+
auto u = dsga::basic_vector(p1);
134+
135+
if (dsga::dot(hyp1, hyp1) > dsga::dot(hyp2, hyp2))
136+
{
137+
hyp = hyp2;
138+
v = v2;
139+
u = dsga::basic_vector(p2);
140+
}
141+
142+
return u + dsga::outerProduct(v, v) * hyp / dsga::dot(v, v);
143+
}
144+
145+
// gives minimum distance from point to a line made from line segment p1 <=> p2
146+
template <bool W1, dsga::floating_point_scalar T, typename D1, bool W2, typename D2, bool W3, typename D3>
147+
constexpr T distance_to_line(const dsga::vector_base<W1, T, 3u, D1> &point,
148+
const dsga::vector_base<W2, T, 3u, D2> &p1,
149+
const dsga::vector_base<W3, T, 3u, D3> &p2) noexcept
150+
{
151+
auto hyp = point - p1;
152+
auto v1 = p2 - p1;
153+
auto t = dsga::dot(hyp, v1) / dsga::dot(v1, v1);
154+
155+
return dsga::length(hyp - (t * v1));
156+
}
157+
158+
// project a point in 3D space to the closest point on a plane, where plane defined by 3 CCW points
159+
template <bool W1, dsga::floating_point_scalar T, typename D1, bool W2, typename D2, bool W3, typename D3, bool W4, typename D4>
160+
constexpr auto project_to_plane1(const dsga::vector_base<W1, T, 3u, D1> &point,
161+
const dsga::vector_base<W2, T, 3u, D2> &p1,
162+
const dsga::vector_base<W3, T, 3u, D3> &p2,
163+
const dsga::vector_base<W4, T, 3u, D4> &p3) noexcept
164+
{
165+
auto p = [](const auto &u, auto &v, const auto &w) { return dsga::cross_matrix(v - u) * (w - u); };
166+
auto p_val = p(p1, p2, p3);
167+
auto p_cross = dsga::cross_matrix(p_val);
168+
169+
return p1 - p_cross * p_cross * (point - p1) / dsga::dot(p_val, p_val);
170+
}
171+
172+
// same as above, different implementation
173+
template <bool W1, dsga::floating_point_scalar T, typename D1, bool W2, typename D2, bool W3, typename D3, bool W4, typename D4>
174+
constexpr auto project_to_plane2(const dsga::vector_base<W1, T, 3u, D1> &point,
175+
const dsga::vector_base<W2, T, 3u, D2> &p1,
176+
const dsga::vector_base<W3, T, 3u, D3> &p2,
177+
const dsga::vector_base<W4, T, 3u, D4> &p3) noexcept
178+
{
179+
auto triangle_norm = [](const auto &u, auto &v, const auto &w) { return dsga::cross_matrix(v - u) * (w - u); };
180+
auto N = triangle_norm(p1, p2, p3);
181+
auto d = dsga::dot(N, p1);
182+
183+
return point - ((dsga::dot(N, point) - d) / dsga::dot(N, N)) * N;
184+
}
185+
186+
#if ATTENUATE_ROUNDOFF
187+
template <bool W1, dsga::floating_point_scalar T, typename D1, bool W2, typename D2, bool W3, typename D3, bool W4, typename D4>
188+
constexpr auto project_to_plane(const dsga::vector_base<W1, T, 3u, D1> &point,
189+
const dsga::vector_base<W2, T, 3u, D2> &p1,
190+
const dsga::vector_base<W3, T, 3u, D3> &p2,
191+
const dsga::vector_base<W4, T, 3u, D4> &p3) noexcept
192+
{
193+
auto p = [](const auto &u, auto &v, const auto &w) { return dsga::cross_matrix(v - u) * (w - u); };
194+
auto delta_u = p3 - p2;
195+
auto delta_v = p1 - p3;
196+
auto delta_w = p2 - p1;
197+
198+
auto md2u = dsga::dot(delta_u, delta_u);
199+
auto md2v = dsga::dot(delta_v, delta_v);
200+
auto md2w = dsga::dot(delta_w, delta_w);
201+
dsga::basic_vector<T, 3u> p_val{0};
202+
203+
// for p_val, choose u that maximizes (v - w).length()
204+
if (md2u > md2v)
205+
{
206+
if (md2u > md2w)
207+
p_val = p(p1, p2, p3);
208+
else
209+
p_val = p(p3, p1, p2);
210+
}
211+
else
212+
{
213+
if (md2v > md2w)
214+
p_val = p(p2, p3, p1);
215+
else
216+
p_val = p(p3, p1, p2);
217+
}
218+
219+
auto p_cross = dsga::cross_matrix(p_val);
220+
221+
dsga::basic_vector<T, 3u> anchor{0};
222+
auto m2u = dsga::dot(p1, p1);
223+
auto m2v = dsga::dot(p2, p2);
224+
auto m2w = dsga::dot(p3, p3);
225+
226+
// for anchor point, choose u that minimizes u.length()
227+
if (m2u < m2v)
228+
{
229+
if (m2u < m2w)
230+
anchor = p1;
231+
else
232+
anchor = p3;
233+
}
234+
else
235+
{
236+
if (m2u < m2v)
237+
anchor = p1;
238+
else
239+
anchor = p2;
240+
}
241+
242+
return anchor - p_cross * p_cross * (point - anchor) / dsga::dot(p_val, p_val);
243+
}
244+
#endif

0 commit comments

Comments
 (0)