Skip to content
GitLab
Commit 50aef922 authored by Philip Trettner's avatar Philip Trettner
Browse files

Merge branch 'agrabowy' into 'develop' 

Added uniform over tube, cylinder, hemisphere, capsule, and cone. Added project onto triangle, disk, circle, hemisphere2, and cone. Added contains for cone

See merge request !49
parents db2bfb53 2b7c0814
Hide whitespace changes
Inline Side-by-side
src/typed-geometry/functions/angle.hh
View file @ 50aef922
......@@ -26,8 +26,8 @@ TG_NODISCARD constexpr angle_t<fractional_result<ScalarT>> angle_between(vec<D,
template <int D, class ScalarT>
TG_NODISCARD constexpr angle_t<fractional_result<ScalarT>> angle_between(dir<D, ScalarT> const& a, dir<D, ScalarT> const& b)
{
constexpr auto lower = decltype(dot(normal(a), normal(b)))(-1);
constexpr auto upper = decltype(dot(normal(a), normal(b)))(1);
constexpr auto lower = decltype(dot(a, b))(-1);
constexpr auto upper = decltype(dot(a, b))(1);
return acos(clamp(dot(a, b), lower, upper));
}
......
src/typed-geometry/functions/contains.hh
View file @ 50aef922
......@@ -38,7 +38,7 @@ TG_NODISCARD constexpr auto contains(A const& a, pos<D, ScalarT> const& p, Scala
template <int D, class ScalarT>
TG_NODISCARD constexpr bool contains(pos<D, ScalarT> const& b, pos<D, ScalarT> const& o, ScalarT eps = ScalarT(0))
{
if (eps > 0)
if (eps > ScalarT(0))
return distance_sqr(b, o) < eps * eps;
return b == o;
}
......@@ -119,7 +119,7 @@ TG_NODISCARD constexpr bool contains(triangle<2, ScalarT> const& t, pos<2, Scala
auto A1 = cross(pv2, pv0);
auto A2 = cross(pv0, pv1);
if (eps > 0)
if (eps > ScalarT(0))
return ((A0 >= -std::copysign(eps, A0)) == (A1 >= -std::copysign(eps, A0))) && //
((A1 >= -std::copysign(eps, A0)) == (A2 >= -std::copysign(eps, A0)));
......@@ -138,17 +138,17 @@ TG_NODISCARD constexpr bool contains(triangle<3, ScalarT> const& t, pos<3, Scala
// checking whether point lies on right side of any edge
auto e = t.pos1 - t.pos0;
auto C = cross(e, p - t.pos0);
if (dot(n, C) < 0)
if (dot(n, C) < ScalarT(0))
return false;
e = t.pos2 - t.pos1;
C = cross(e, p - t.pos1);
if (dot(n, C) < 0)
if (dot(n, C) < ScalarT(0))
return false;
e = t.pos0 - t.pos2;
C = cross(e, p - t.pos2);
if (dot(n, C) < 0)
if (dot(n, C) < ScalarT(0))
return false;
// point always on left side
......@@ -168,7 +168,7 @@ TG_NODISCARD constexpr bool contains(cylinder<3, ScalarT> const& c, pos<3, Scala
auto hsqd = length_sqr(ad);
auto rsqd = pow2(c.radius);
if (d0 < 0 || d0 > hsqd) // behind a cap
if (d0 < ScalarT(0) || d0 > hsqd) // behind a cap
return false;
// check whether distance from p to axis is less or equal to radius
......@@ -191,4 +191,19 @@ TG_NODISCARD constexpr bool contains(disk<2, ScalarT> const& d, pos<2, ScalarT>
return distance_sqr(d.center, p) <= r * r;
}
template <class ScalarT>
TG_NODISCARD constexpr bool contains(cone<3, ScalarT> const& c, pos<3, ScalarT> const& p, ScalarT eps = ScalarT(0))
{
auto center = c.base.center - eps * c.base.normal;
if (dot(p - center, c.base.normal) < ScalarT(0))
return false; // Not inside if on the other side of the base
auto apex = c.base.center + (c.height + eps) * c.base.normal;
auto pRing = c.base.center + (c.base.radius + eps) * any_normal(c.base.normal);
// Inside iff the point is closer to the axis (in terms of angle wrt. the apex) than some point on the outer boundary
return dot(-c.base.normal, normalize(p - apex)) > dot(-c.base.normal, normalize(pRing - apex));
}
} // namespace tg
src/typed-geometry/functions/project.hh
View file @ 50aef922
......@@ -89,6 +89,51 @@ TG_NODISCARD constexpr pos<D, ScalarT> project(pos<D, ScalarT> const& p, halfspa
return p - pl.normal * tg::max(ScalarT(0), dot(p, pl.normal) - pl.dis);
}
template <class ScalarT>
TG_NODISCARD constexpr pos<3, ScalarT> project(pos<3, ScalarT> const& p, triangle<3, ScalarT> const& t)
{
auto pPlane = project(p, hyperplane<3, ScalarT>(normal(t), t.pos0));
if (contains(t, pPlane))
return pPlane;
auto p0 = project(pPlane, segment<3, ScalarT>(t.pos0, t.pos1));
auto p1 = project(pPlane, segment<3, ScalarT>(t.pos0, t.pos2));
auto p2 = project(pPlane, segment<3, ScalarT>(t.pos1, t.pos2));
auto d0 = distance_sqr(p0, pPlane);
auto d1 = distance_sqr(p1, pPlane);
auto d2 = distance_sqr(p2, pPlane);
if (d0 <= d1 && d0 <= d2)
return p0;
else if (d1 <= d2)
return p1;
else
return p2;
}
template <class ScalarT>
TG_NODISCARD constexpr pos<2, ScalarT> project(pos<2, ScalarT> const& p, triangle<2, ScalarT> const& t)
{
if (contains(t, p))
return p;
auto p0 = project(p, segment<2, ScalarT>(t.pos0, t.pos1));
auto p1 = project(p, segment<2, ScalarT>(t.pos0, t.pos2));
auto p2 = project(p, segment<2, ScalarT>(t.pos1, t.pos2));
auto d0 = distance_sqr(p0, p);
auto d1 = distance_sqr(p1, p);
auto d2 = distance_sqr(p2, p);
if (d0 <= d1 && d0 <= d2)
return p0;
else if (d1 <= d2)
return p1;
else
return p2;
}
template <int D, class ScalarT>
TG_NODISCARD constexpr pos<D, ScalarT> project(pos<D, ScalarT> const& p, sphere<D, ScalarT> const& sp)
{
......@@ -98,18 +143,32 @@ TG_NODISCARD constexpr pos<D, ScalarT> project(pos<D, ScalarT> const& p, sphere<
return sp.center + dir_to_p * sp.radius;
}
template <int D, class ScalarT>
TG_NODISCARD constexpr pos<D, ScalarT> project(pos<D, ScalarT> const& p, hemisphere<D, ScalarT> const& h)
template <class ScalarT>
TG_NODISCARD constexpr pos<3, ScalarT> project(pos<3, ScalarT> const& p, hemisphere<3, ScalarT> const& h) // boundary, including caps
{
auto dir_to_p = tg::normalize_safe(p - h.center);
if (is_zero_vector(dir_to_p))
return h.center + h.normal * h.radius;
if (dot(dir_to_p, h.normal) < 0)
return project(p, disk<D, ScalarT>(h.center, h.radius, h.normal));
if (dot(dir_to_p, h.normal) >= ScalarT(0))
return h.center + dir_to_p * h.radius;
return h.center + dir_to_p * h.radius;
return project(p, disk<3, ScalarT>(h.center, h.radius, h.normal));
}
template <class ScalarT>
TG_NODISCARD constexpr pos<2, ScalarT> project(pos<2, ScalarT> const& p, hemisphere<2, ScalarT> const& h) // boundary, including caps
{
auto dir_to_p = tg::normalize_safe(p - h.center);
if (is_zero_vector(dir_to_p))
return h.center + h.normal * h.radius;
if (dot(dir_to_p, h.normal) >= ScalarT(0))
return h.center + dir_to_p * h.radius;
auto v = perpendicular(h.normal) * h.radius;
return project(p, segment<2, ScalarT>(h.center - v, h.center + v));
}
template <int D, class ScalarT>
......@@ -122,7 +181,7 @@ TG_NODISCARD constexpr pos<D, ScalarT> project(pos<D, ScalarT> const& p, ball<D,
}
template <class ScalarT>
TG_NODISCARD constexpr pos<3, ScalarT> project(pos<3, ScalarT> const& p, tube<3, ScalarT> const& t)
TG_NODISCARD constexpr pos<3, ScalarT> project(pos<3, ScalarT> const& p, tube<3, ScalarT> const& t) // boundary
{
auto lp = project(p, t.axis);
auto dir = normalize_safe(p - lp);
......@@ -147,6 +206,18 @@ TG_NODISCARD constexpr pos<3, ScalarT> project(pos<3, ScalarT> const& p, disk<3,
return d.center + dir * d.radius;
}
template <class ScalarT>
TG_NODISCARD constexpr pos<2, ScalarT> project(pos<2, ScalarT> const& p, disk<2, ScalarT> const& d)
{
if (distance_sqr(p, d.center) <= d.radius * d.radius)
return p;
auto dir = normalize_safe(p - d.center);
if (is_zero_vector(dir))
dir = tg::dir<2, ScalarT>::pos_x;
return d.center + dir * d.radius;
}
template <class ScalarT>
TG_NODISCARD constexpr pos<3, ScalarT> project(pos<3, ScalarT> const& p, circle<3, ScalarT> const& c)
......@@ -159,9 +230,18 @@ TG_NODISCARD constexpr pos<3, ScalarT> project(pos<3, ScalarT> const& p, circle<
return c.center + dir * c.radius;
}
template <class ScalarT>
TG_NODISCARD constexpr pos<2, ScalarT> project(pos<2, ScalarT> const& p, circle<2, ScalarT> const& c)
{
auto dir = normalize_safe(p - c.center);
if (is_zero_vector(dir))
dir = tg::dir<2, ScalarT>::pos_x;
return c.center + dir * c.radius;
}
template <class ScalarT>
TG_NODISCARD constexpr pos<3, ScalarT> project(pos<3, ScalarT> const& p, cylinder<3, ScalarT> const& c)
TG_NODISCARD constexpr pos<3, ScalarT> project(pos<3, ScalarT> const& p, cylinder<3, ScalarT> const& c) // boundary, including caps
{
auto dir = direction(c);
......@@ -182,19 +262,30 @@ TG_NODISCARD constexpr pos<3, ScalarT> project(pos<3, ScalarT> const& p, cylinde
}
template <class ScalarT>
TG_NODISCARD constexpr pos<3, ScalarT> project(pos<3, ScalarT> const& p, capsule<3, ScalarT> const& c)
TG_NODISCARD constexpr pos<3, ScalarT> project(pos<3, ScalarT> const& p, capsule<3, ScalarT> const& c) // boundary, including caps
{
auto t = coordinates(c.axis, p);
if (t < 0)
if (t < ScalarT(0))
return project(p, sphere<3, ScalarT>(c.axis.pos0, c.radius));
if (t > 1)
if (t > ScalarT(1))
return project(p, sphere<3, ScalarT>(c.axis.pos1, c.radius));
return project(p, tube<3, ScalarT>(c.axis, c.radius));
}
template <class ScalarT>
TG_NODISCARD constexpr pos<3, ScalarT> project(pos<3, ScalarT> const& p, cone<3, ScalarT> const& c)
{
auto closestOnBase = project(p, c.base);
auto apex = c.base.center + c.height * c.base.normal;
if (dot(p - closestOnBase, closestOnBase - apex) >= ScalarT(0))
return closestOnBase;
return project(p, inf_cone<3, ScalarT>(apex, -c.base.normal, ScalarT(2) * angle_between(normalize(closestOnBase - apex), -c.base.normal)));
}
template <class ScalarT>
TG_NODISCARD constexpr pos<3, ScalarT> project(pos<3, ScalarT> const& p, inf_cone<3, ScalarT> const& icone)
{
......@@ -211,13 +302,13 @@ TG_NODISCARD constexpr pos<3, ScalarT> project(pos<3, ScalarT> const& p, inf_con
if (tg::are_collinear(p_apex_dir, static_cast<vec<3, ScalarT>>(icone.opening_dir)))
{
// p is "above" the apex
if (dot(p_apex_dir, icone.opening_dir) < 0)
if (dot(p_apex_dir, icone.opening_dir) < ScalarT(0))
return icone.apex;
// any point on the cone in normal direction from p is the closest point
auto h = tg::length(p_apex);
auto l = tg::cos(icone.opening_angle / 2) * h;
auto r = tan(icone.opening_angle / 2);
auto l = tg::cos(icone.opening_angle / ScalarT(2)) * h;
auto r = tan(icone.opening_angle / ScalarT(2));
dir_t ortho_dir = tg::any_normal(icone.opening_dir);
auto pt_on_cone = icone.apex + icone.opening_dir + ortho_dir * r;
dir_t on_surface_dir = normalize(pt_on_cone - icone.apex);
......@@ -231,21 +322,21 @@ TG_NODISCARD constexpr pos<3, ScalarT> project(pos<3, ScalarT> const& p, inf_con
dir_t y_axis = -icone.opening_dir;
dir_t plane_normal = normalize(cross(p - c, vec<3, ScalarT>(y_axis)));
dir_t x_axis = normalize(cross(y_axis, plane_normal));
if (dot(p - c, x_axis) < 0)
if (dot(p - c, x_axis) < ScalarT(0))
x_axis = -x_axis;
// construct the 2D surface normal of the cone in the plane
auto r = tan(icone.opening_angle / 2);
vec2_t r_ = {r, 0};
auto r = tan(icone.opening_angle / ScalarT(2));
vec2_t r_ = {r, ScalarT(0)};
vec2_t p_ = {dot(p - c, x_axis), dot(p - c, y_axis)};
vec2_t peak_ = {0, 1};
vec2_t peak_ = {ScalarT(0), ScalarT(1)};
dir2_t r_vec = normalize(r_ - peak_);
dir2_t n_ = tg::perpendicular(r_vec);
if (n_.y < 0)
if (n_.y < ScalarT(0))
n_ = -n_;
// reconstruct 3D closest point
if (dot(r_vec, p_ - peak_) > 0)
if (dot(r_vec, p_ - peak_) > ScalarT(0))
{
auto d = dot(p_ - peak_, n_);
auto proj_p2 = p_ - d * n_;
......
src/typed-geometry/functions/uniform.hh
View file @ 50aef922
......@@ -220,6 +220,50 @@ TG_NODISCARD constexpr pos<3, ScalarT> uniform(Rng& rng, disk<3, ScalarT> const&
return d.center + direction.x * x + direction.y * y;
}
template <class ScalarT, class Rng>
TG_NODISCARD constexpr pos<3, ScalarT> uniform(Rng& rng, tube<3, ScalarT> const& t) // boundary
{
auto c = circle<3, ScalarT>(pos<3, ScalarT>::zero, t.radius, normalize(t.axis.pos1 - t.axis.pos0));
return uniform(rng, t.axis) + vec<3, ScalarT>(uniform(rng, c));
}
template <class ScalarT, class Rng>
TG_NODISCARD constexpr pos<3, ScalarT> uniform(Rng& rng, cylinder<3, ScalarT> const& c) // boundary, including caps
{
auto x = c.axis.pos1 - c.axis.pos0;
auto h = length(x);
auto sideArea = ScalarT(2) * c.radius * h; // * Pi, but that does not matter here
auto capArea = c.radius * c.radius; // * Pi
auto totalArea = ScalarT(2) * capArea + sideArea;
auto part = detail::uniform01<ScalarT>(rng) * totalArea;
if (part < sideArea) // Uniform sampling on cylinder side
return uniform(rng, tube<3, ScalarT>(c.axis, c.radius));
// Otherwise sampling on one of the caps
auto capDisk = disk<3, ScalarT>(part < sideArea + capArea ? c.axis.pos0 : c.axis.pos1, c.radius, normalize(x));
return uniform(rng, capDisk);
}
template <class ScalarT, class Rng>
TG_NODISCARD constexpr pos<3, ScalarT> uniform(Rng& rng, capsule<3, ScalarT> const& c) // boundary, including caps (does no_caps make sense here?)
{
auto x = c.axis.pos1 - c.axis.pos0;
auto h = length(x);
auto sideArea = ScalarT(2) * c.radius * h; // * Pi, but that does not matter here
auto capArea = ScalarT(2) * c.radius * c.radius; // * Pi
auto totalArea = ScalarT(2) * capArea + sideArea;
auto part = detail::uniform01<ScalarT>(rng) * totalArea;
if (part < sideArea) // Uniform sampling on capsule side
return uniform(rng, tube<3, ScalarT>(c.axis, c.radius));
// Otherwise sampling on one of the caps
auto capHemi = hemisphere<3, ScalarT>();
capHemi.radius = c.radius;
capHemi.center = part < sideArea + capArea ? c.axis.pos0 : c.axis.pos1;
capHemi.normal = part < sideArea + capArea ? -normalize(x) : normalize(x);
return uniform(rng, capHemi);
}
template <int D, class ScalarT, class Rng>
TG_NODISCARD constexpr pos<D, ScalarT> uniform(Rng& rng, sphere<D, ScalarT> const& s)
{
......@@ -246,6 +290,35 @@ TG_NODISCARD constexpr pos<D, ScalarT> uniform(Rng& rng, ball<D, ScalarT> const&
}
}
template <class ScalarT, class Rng>
TG_NODISCARD constexpr pos<3, ScalarT> uniform(Rng& rng, cone<3, ScalarT> const& c) // boundary, no_caps (not on base)
{
auto ub = tg::aabb<2, ScalarT>::minus_one_to_one;
while (true)
{
auto p = uniform_vec(rng, ub);
auto l = length_sqr(p);
if (l <= ScalarT(1))
{
p *= c.base.radius;
auto x = any_normal(c.base.normal);
auto y = cross(c.base.normal, x);
return c.base.center + p.x * x + p.y * y + (ScalarT(1) - sqrt(l)) * c.base.normal * c.height;
}
}
}
template <int D, class ScalarT, class Rng>
TG_NODISCARD constexpr pos<D, ScalarT> uniform(Rng& rng, hemisphere<D, ScalarT> const& h) // boundary, no_caps (not on base)
{
auto p = uniform(rng, sphere<D, ScalarT>(h.center, h.radius));
auto v = p - h.center;
if (dot(v, h.normal) >= ScalarT(0))
return p;
else
return h.center - v;
}
template <int D, class ScalarT, class Rng, class = enable_if<is_floating_point<ScalarT>>>
[[deprecated("potentially misleading operation. use uniform_vec(rng, tg::aabb3(..)) or uniform_vec(rng, tg::segment3(..)) depending on your intended "
"semantics")]] TG_NODISCARD constexpr vec<D, ScalarT>
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment