Commit 703e9987 authored by Julius NehringWirxel's avatar Julius NehringWirxel

Merge remote-tracking branch 'origin/develop' into feature/eigenvaluedecomposition

parents 05db9ddb 9e2b6e80
This diff is collapsed.
This diff is collapsed.
#pragma once
#include <typed-geometry/types/vec.hh>
namespace tg
{
namespace noise
{
// Helper functions for simplex.hh and perlin.hh
template <class ScalarT>
ScalarT mod289(ScalarT const x)
{
return x - floor(x * (ScalarT(1) / ScalarT(289))) * ScalarT(289);
}
template <int D, class ScalarT>
pos<D, ScalarT> mod289(pos<D, ScalarT> const& x)
{
return pos<D, ScalarT>(x - floor(x * (ScalarT(1) / ScalarT(289))) * ScalarT(289));
}
template <class ScalarT>
ScalarT permute(ScalarT const x)
{
return mod289(((x * ScalarT(34)) + ScalarT(1)) * x);
}
template <int D, class ScalarT>
pos<D, ScalarT> permute(const pos<D, ScalarT>& x)
{
auto ret = comp<D, ScalarT>(x) * ((comp<D, ScalarT>(x) * ScalarT(34)) + ScalarT(1));
return mod289(pos<D, ScalarT>(ret));
}
template <class ScalarT>
ScalarT taylorInvSqrt(ScalarT const r)
{
return ScalarT(1.79284291400159) - ScalarT(0.85373472095314) * r;
}
template <int D, class ScalarT>
pos<D, ScalarT> taylorInvSqrt(pos<D, ScalarT> const& r)
{
return ScalarT(1.79284291400159) - ScalarT(0.85373472095314) * r;
}
template <int D, class ScalarT>
pos<D, ScalarT> fade(const vec<D, ScalarT>& t)
{
auto ret = comp<D, ScalarT>(t);
return pos<D, ScalarT>(ret * ret * ret * (ret * (ret * ScalarT(6) - ScalarT(15)) + ScalarT(10)));
}
// See https://www.khronos.org/registry/OpenGL-Refpages/gl4/html/step.xhtml
template <int D, class ScalarT>
vec<D, ScalarT> step(const vec<D, ScalarT>& edge, const vec<D, ScalarT>& x)
{
auto ret = vec<D, ScalarT>::zero;
for (auto i = 0; i < D; i++)
{
if (x[i] >= edge[i])
ret[i] = ScalarT(1);
// else: 0
}
return ret;
}
// See https://www.khronos.org/registry/OpenGL-Refpages/gl4/html/lessThan.xhtml
template <class ScalarT>
vec<4, ScalarT> lessThan(vec<4, ScalarT> const& x, vec<4, ScalarT> const& y)
{
vec<4, ScalarT> ret;
for (auto i = 0; i < 4; i++)
{
ret[i] = x[i] < y[i] ? ScalarT(1) : ScalarT(0);
}
return ret;
}
template <class ScalarT>
vec<4, ScalarT> grad4(ScalarT j, vec<4, ScalarT> ip)
{
auto ones = vec<4, ScalarT>::one;
ones.w = -1;
vec<4, ScalarT> p, s;
auto res = floor(fract(vec<3, ScalarT>(comp<3, ScalarT>(j) * comp<3, ScalarT>(ip.x, ip.y, ip.z))) * ScalarT(7)) * ip.z - ScalarT(1);
p.x = res.x;
p.y = res.y;
p.z = res.z;
p.w = ScalarT(1.5) - dot(abs(vec<3, ScalarT>(p.x, p.y, p.z)), vec<3, ScalarT>(ones.x, ones.y, ones.z));
s = vec<4, ScalarT>(lessThan(p, vec<4, ScalarT>::zero));
auto xyz = comp<3, ScalarT>(p.x, p.y, p.z) + (comp<3, ScalarT>(s.x, s.y, s.z) * ScalarT(2) - comp<3, ScalarT>(1)) * comp<3, ScalarT>(s.w);
p.x = xyz[0];
p.y = xyz[1];
p.z = xyz[2];
return p;
}
template <int D, class ScalarT>
vec<D, ScalarT> clamp(vec<D, ScalarT> const& v, const ScalarT minval, const ScalarT maxval)
{
vec<D, ScalarT> ret;
for (auto i = 0; i < D; i++)
{
ret[i] = max(min(v[i], maxval), minval);
}
return ret;
}
// Helper functions (and data) for simplex noise https://github.com/SRombauts/SimplexNoise/blob/master/src/SimplexNoise.cpp
/**
* Permutation table. This is just a random jumble of all numbers 0-255.
*
* This produce a repeatable pattern of 256, but Ken Perlin stated
* that it is not a problem for graphic texture as the noise features disappear
* at a distance far enough to be able to see a repeatable pattern of 256.
*
* This needs to be exactly the same for all instances on all platforms,
* so it's easiest to just keep it as static explicit data.
* This also removes the need for any initialisation of this class.
*
* Note that making this an uint32_t[] instead of a uint8_t[] might make the
* code run faster on platforms with a high penalty for unaligned single
* byte addressing. Intel x86 is generally single-byte-friendly, but
* some other CPUs are faster with 4-aligned reads.
* However, a char[] is smaller, which avoids cache trashing, and that
* is probably the most important aspect on most architectures.
* This array is accessed a *lot* by the noise functions.
* A vector-valued noise over 3D accesses it 96 times, and a
* float-valued 4D noise 64 times. We want this to fit in the cache!
*/
static const tg::u8 perm[256] = {
151, 160, 137, 91, 90, 15, 131, 13, 201, 95, 96, 53, 194, 233, 7, 225, 140, 36, 103, 30, 69, 142, 8, 99, 37, 240, 21, 10, 23,
190, 6, 148, 247, 120, 234, 75, 0, 26, 197, 62, 94, 252, 219, 203, 117, 35, 11, 32, 57, 177, 33, 88, 237, 149, 56, 87, 174, 20,
125, 136, 171, 168, 68, 175, 74, 165, 71, 134, 139, 48, 27, 166, 77, 146, 158, 231, 83, 111, 229, 122, 60, 211, 133, 230, 220, 105, 92,
41, 55, 46, 245, 40, 244, 102, 143, 54, 65, 25, 63, 161, 1, 216, 80, 73, 209, 76, 132, 187, 208, 89, 18, 169, 200, 196, 135, 130,
116, 188, 159, 86, 164, 100, 109, 198, 173, 186, 3, 64, 52, 217, 226, 250, 124, 123, 5, 202, 38, 147, 118, 126, 255, 82, 85, 212, 207,
206, 59, 227, 47, 16, 58, 17, 182, 189, 28, 42, 223, 183, 170, 213, 119, 248, 152, 2, 44, 154, 163, 70, 221, 153, 101, 155, 167, 43,
172, 9, 129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224, 232, 178, 185, 112, 104, 218, 246, 97, 228, 251, 34, 242, 193, 238, 210, 144,
12, 191, 179, 162, 241, 81, 51, 145, 235, 249, 14, 239, 107, 49, 192, 214, 31, 181, 199, 106, 157, 184, 84, 204, 176, 115, 121, 50, 45,
127, 4, 150, 254, 138, 236, 205, 93, 222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 215, 61, 156, 180};
/**
* Helper function to hash an integer using the above permutation table
*
* This inline function costs around 1ns, and is called N+1 times for a noise of N dimension.
*
* Using a real hash function would be better to improve the "repeatability of 256" of the above permutation table,
* but fast integer Hash functions uses more time and have bad random properties.
*
* @param[in] i Integer value to hash
*
* @return 8-bits hashed value
*/
static inline tg::u8 hash(tg::i32 i) { return perm[static_cast<tg::u8>(i)]; }
/**
* Helper function to compute gradients-dot-residual vectors (1D)
*
* @note that these generate gradients of more than unit length. To make
* a close match with the value range of classic Perlin noise, the final
* noise values need to be rescaled to fit nicely within [-1,1].
* (The simplex noise functions as such also have different scaling.)
* Note also that these noise functions are the most practical and useful
* signed version of Perlin noise.
*
* @param[in] hash hash value
* @param[in] x distance to the corner
*
* @return gradient value
*/
template <class ScalarT>
static ScalarT grad(tg::i32 hash, ScalarT x)
{
const tg::i32 h = hash & 0x0F; // Convert low 4 bits of hash code
ScalarT grad = ScalarT(1) + (h & 7); // Gradient value 1.0, 2.0, ..., 8.0
if ((h & 8) != 0)
grad = -grad; // Set a random sign for the gradient
// float grad = gradients1D[h]; // NOTE : Test of Gradient look-up table instead of the above
return (grad * x); // Multiply the gradient with the distance
}
/**
* Helper functions to compute gradients-dot-residual vectors (2D)
*
* @param[in] hash hash value
* @param[in] x x coord of the distance to the corner
* @param[in] y y coord of the distance to the corner
*
* @return gradient value
*/
template <class ScalarT>
static ScalarT grad(tg::i32 hash, ScalarT x, ScalarT y)
{
const tg::i32 h = hash & 0x3F; // Convert low 3 bits of hash code
const ScalarT u = h < 4 ? x : y; // into 8 simple gradient directions,
const ScalarT v = h < 4 ? y : x;
return ((h & 1) ? -u : u) + ((h & 2) ? -ScalarT(2) * v : ScalarT(2) * v); // and compute the dot product with (x,y).
}
/**
* Helper functions to compute gradients-dot-residual vectors (3D)
*
* @param[in] hash hash value
* @param[in] x x coord of the distance to the corner
* @param[in] y y coord of the distance to the corner
* @param[in] z z coord of the distance to the corner
*
* @return gradient value
*/
template <class ScalarT>
static ScalarT grad(tg::i32 hash, ScalarT x, ScalarT y, ScalarT z)
{
tg::i32 h = hash & 15; // Convert low 4 bits of hash code into 12 simple
ScalarT u = h < 8 ? x : y; // gradient directions, and compute dot product.
ScalarT v = h < 4 ? y : h == 12 || h == 14 ? x : z; // Fix repeats at h = 12 to 15
return ((h & 1) ? -u : u) + ((h & 2) ? -v : v);
}
/**
* Computes the largest integer value not greater than the float one
*
* This method is faster than using (int32_t)std::floor(fp).
*
* I measured it to be approximately twice as fast:
* float: ~18.4ns instead of ~39.6ns on an AMD APU),
* double: ~20.6ns instead of ~36.6ns on an AMD APU),
* Reference: http://www.codeproject.com/Tips/700780/Fast-floor-ceiling-functions
*
* @param[in] fp float input value
*
* @return largest integer value not greater than fp
*/
template <class ScalarT>
static inline tg::i32 fastfloor(ScalarT fp)
{
auto i = static_cast<tg::i32>(fp);
return (fp < i) ? (i - 1) : (i);
}
} // namespace noise
} // namespace tg
......@@ -39,11 +39,13 @@ template <int D, class ScalarT, class TraitsT>
return this->min + (this->max - this->min) * size(c);
}
/// Note that the box goes from -1 to 1 instead of the usual 0 to 1
template <int D, class ScalarT, class TraitsT>
[[nodiscard]] constexpr pos<D, ScalarT> box<D, ScalarT, D, TraitsT>::operator[](comp<D, ScalarT> const& c) const
{
return this->center + this->half_extents * vec(c);
}
/// Note that the box goes from -1 to 1 instead of the usual 0 to 1
template <class ScalarT, class TraitsT>
[[nodiscard]] constexpr pos<3, ScalarT> box<2, ScalarT, 3, TraitsT>::operator[](comp<2, ScalarT> const& c) const
{
......
......@@ -20,9 +20,16 @@ struct sum_of_pos
constexpr sum_of_pos() = default;
/* implicit */ constexpr sum_of_pos(pos_t const& p) : accum(p) {}
constexpr pos<D, div_t> operator/(div_t const& rhs) const { return accum / rhs; }
constexpr pos<D, div_t> operator*(div_t const& rhs) const { return accum * rhs; }
constexpr auto operator/(div_t const& rhs) const
{
using T = std::common_type_t<ScalarT, div_t>;
return pos<D, T>(accum) / T(rhs);
}
constexpr auto operator*(div_t const& rhs) const
{
using T = std::common_type_t<ScalarT, div_t>;
return pos<D, T>(accum) * T(rhs);
}
constexpr sum_of_pos operator+(sum_of_pos const& rhs) const { return accum + vec(rhs); }
constexpr sum_of_pos operator+(pos_t const& rhs) const { return accum + vec(rhs); }
......
#pragma once
#include <typed-geometry/detail/noise/perlin.hh>
#include <typed-geometry/detail/noise/simplex.hh>
......@@ -4,6 +4,9 @@
#include <typed-geometry/types/objects/aabb.hh>
#include <typed-geometry/types/objects/box.hh>
#include <typed-geometry/types/objects/capsule.hh>
#include <typed-geometry/types/objects/hemisphere.hh>
#include <typed-geometry/types/objects/pyramid.hh>
#include <typed-geometry/types/objects/quad.hh>
#include <typed-geometry/types/objects/segment.hh>
#include <typed-geometry/types/objects/sphere.hh>
......@@ -29,6 +32,34 @@ template <int D, class ScalarT, class TraitsT>
return {s.center - s.radius, s.center + s.radius};
}
template <class ScalarT, class TraitsT>
[[nodiscard]] constexpr aabb<3, ScalarT> aabb_of(sphere<2, ScalarT, 3, TraitsT> const& s)
{
// See http://www.iquilezles.org/www/articles/diskbbox/diskbbox.htm
const auto e = abs(s.radius) * sqrt(ScalarT(1) - comp(s.normal) * comp(s.normal));
return {s.center - e, s.center + e};
}
template <class ScalarT, class TraitsT>
[[nodiscard]] constexpr aabb<1, ScalarT> aabb_of(hemisphere<1, ScalarT, TraitsT> const& h)
{
return aabb_of(h.center, h.center + h.normal * h.radius);
}
template <class ScalarT, class TraitsT>
[[nodiscard]] constexpr aabb<2, ScalarT> aabb_of(hemisphere<2, ScalarT, TraitsT> const& h)
{
const auto baseVec = h.radius * perpendicular(h.normal);
const auto sphereCorner = h.center + h.radius * sign(vec(h.normal));
return aabb_of(h.center - baseVec, h.center + baseVec, sphereCorner);
}
template <class ScalarT, class TraitsT>
[[nodiscard]] constexpr aabb<3, ScalarT> aabb_of(hemisphere<3, ScalarT, TraitsT> const& h)
{
const auto disk = sphere<2, ScalarT, 3>(h.center, h.radius, h.normal);
const auto sphereCorner = h.center + h.radius * sign(vec(h.normal));
return aabb_of(disk, sphereCorner);
}
template <int D, class ScalarT>
[[nodiscard]] constexpr aabb<D, ScalarT> aabb_of(segment<D, ScalarT> const& s)
{
......@@ -42,9 +73,9 @@ template <int D, class ScalarT>
}
template <int D, class ScalarT>
[[nodiscard]] constexpr aabb<D, ScalarT> aabb_of(quad<D, ScalarT> const& t)
[[nodiscard]] constexpr aabb<D, ScalarT> aabb_of(quad<D, ScalarT> const& q)
{
return aabb_of(t.pos00, t.pos10, t.pos11, t.pos01);
return aabb_of(q.pos00, q.pos10, q.pos11, q.pos01);
}
template <int ObjectD, class ScalarT, int DomainD, class TraitsT>
......@@ -57,6 +88,26 @@ template <int ObjectD, class ScalarT, int DomainD, class TraitsT>
return {b.center - diag, b.center + diag};
}
template <int D, class ScalarT, class TraitsT>
[[nodiscard]] constexpr aabb<D, ScalarT> aabb_of(capsule<D, ScalarT, TraitsT> const& c)
{
return aabb_of(sphere<D, ScalarT>(c.axis.pos0, c.radius), sphere<D, ScalarT>(c.axis.pos1, c.radius));
}
template <int D, class ScalarT, class TraitsT>
[[nodiscard]] constexpr aabb<D, ScalarT> aabb_of(cylinder<D, ScalarT, TraitsT> const& c)
{
const auto n = normalize(c.axis.pos1 - c.axis.pos0);
return aabb_of(sphere<2, ScalarT, 3>(c.axis.pos0, c.radius, n), sphere<2, ScalarT, 3>(c.axis.pos1, c.radius, n));
}
template <class BaseT, class TraitsT, class ScalarT = typename BaseT::scalar_t>
[[nodiscard]] constexpr aabb<3, ScalarT> aabb_of(pyramid<BaseT, TraitsT> const& p)
{
const auto apex = centroid(p.base) + normal(p.base) * p.height;
return aabb_of(p.base, apex);
}
template <class PrimA, class PrimB, class... PrimsT>
[[nodiscard]] constexpr auto aabb_of(PrimA const& pa, PrimB const& pb, PrimsT const&... prims) -> decltype(aabb_of(pa))
{
......
#pragma once
#include <typed-geometry/detail/scalar_traits.hh>
#include <typed-geometry/detail/special_values.hh>
#include <typed-geometry/types/objects/aabb.hh>
#include <typed-geometry/types/objects/capsule.hh>
#include <typed-geometry/types/objects/cylinder.hh>
#include <typed-geometry/types/objects/quad.hh>
#include <typed-geometry/types/objects/segment.hh>
#include <typed-geometry/types/objects/sphere.hh>
#include <typed-geometry/types/objects/triangle.hh>
......@@ -12,8 +12,6 @@
#include <typed-geometry/detail/operators/ops_pos.hh>
#include <typed-geometry/detail/operators/ops_vec.hh>
#include <typed-geometry/functions/basic/mix.hh>
// returns the arithmetic mean of all points contained in an object
// has variadic versions
......@@ -26,38 +24,44 @@ template <int D, class ScalarT>
}
template <int D, class ScalarT>
[[nodiscard]] constexpr pos<D, ScalarT> centroid(pos<D, ScalarT> const& a, pos<D, ScalarT> const& b)
[[nodiscard]] constexpr pos<D, fractional_result<ScalarT>> centroid(pos<D, ScalarT> const& a, pos<D, ScalarT> const& b)
{
return mix(a, b, ScalarT(0.5));
return (a + b) * fractional_result<ScalarT>(0.5);
}
template <int D, class ScalarT>
[[nodiscard]] constexpr pos<D, ScalarT> centroid(box<D, ScalarT> const& b)
template <int ObjectD, class ScalarT, int DomainD>
[[nodiscard]] constexpr pos<DomainD, ScalarT> centroid(box<ObjectD, ScalarT, DomainD> const& b)
{
return b.center;
}
template <int D, class ScalarT>
[[nodiscard]] constexpr pos<D, fractional_result<ScalarT>> centroid(aabb<D, ScalarT> const& p)
[[nodiscard]] constexpr pos<D, fractional_result<ScalarT>> centroid(aabb<D, ScalarT> const& b)
{
return mix(p.min, p.max, ScalarT(0.5));
return (b.min + b.max) * fractional_result<ScalarT>(0.5);
}
template <int D, class ScalarT>
[[nodiscard]] constexpr pos<D, fractional_result<ScalarT>> centroid(triangle<D, ScalarT> const& p)
[[nodiscard]] constexpr pos<D, fractional_result<ScalarT>> centroid(segment<D, ScalarT> const& s)
{
auto z = tg::zero<pos<D, ScalarT>>();
return z + ((p.pos0 - z) + (p.pos1 - z) + (p.pos2 - z)) / ScalarT(3);
return (s.pos0 + s.pos1) * fractional_result<ScalarT>(0.5);
}
template <int D, class ScalarT>
[[nodiscard]] constexpr pos<D, fractional_result<ScalarT>> centroid(segment<D, ScalarT> const& p)
[[nodiscard]] constexpr pos<D, fractional_result<ScalarT>> centroid(triangle<D, ScalarT> const& t)
{
return mix(p.pos0, p.pos1, fractional_result<ScalarT>(0.5));
return (t.pos0 + t.pos1 + t.pos2) / ScalarT(3);
}
template <int D, class ScalarT>
[[nodiscard]] constexpr pos<D, ScalarT> centroid(sphere<D, ScalarT> const& p)
[[nodiscard]] constexpr pos<D, fractional_result<ScalarT>> centroid(quad<D, ScalarT> const& q)
{
// TODO: The average of the four corners is not the same as the average of all points in the quad
return (q.pos00 + q.pos01 + q.pos10 + q.pos11) * fractional_result<ScalarT>(0.25);
}
template <int ObjectD, class ScalarT, int DomainD>
[[nodiscard]] constexpr pos<DomainD, ScalarT> centroid(sphere<ObjectD, ScalarT, DomainD> const& p)
{
return p.center;
}
......
......@@ -279,6 +279,13 @@ template <class ScalarT>
}
}
template <class ScalarT>
[[nodiscard]] constexpr bool contains(inf_cylinder<3, ScalarT> const& c, pos<3, ScalarT> const& p, dont_deduce<ScalarT> eps = ScalarT(0))
{
auto d = distance(c.axis, p);
return d <= c.radius + eps;
}
template <class ScalarT>
[[nodiscard]] constexpr bool contains(inf_cone<3, ScalarT> const& c, pos<3, ScalarT> const& p, dont_deduce<ScalarT> eps = ScalarT(0))
{
......@@ -286,7 +293,7 @@ template <class ScalarT>
auto apexToP = normalize_safe(p - apex);
if (apexToP == vec<3, ScalarT>::zero)
return true;
return angle_between(dir<3, ScalarT>(apexToP), c.opening_dir) <= c.opening_angle;
return angle_between(dir<3, ScalarT>(apexToP), c.opening_dir) <= ScalarT(0.5) * c.opening_angle; // opening_angle is between the cone surfaces
}
template <class ScalarT>
[[nodiscard]] constexpr bool contains(inf_cone_boundary<3, ScalarT> const& c, pos<3, ScalarT> const& p, dont_deduce<ScalarT> eps = ScalarT(0))
......@@ -297,8 +304,64 @@ template <class ScalarT>
auto apexInnerToP = normalize_safe(p - apexInner);
if (apexOuterToP == vec<3, ScalarT>::zero || apexInnerToP == vec<3, ScalarT>::zero)
return true;
return angle_between(dir<3, ScalarT>(apexOuterToP), c.opening_dir) <= c.opening_angle
&& angle_between(dir<3, ScalarT>(apexInnerToP), c.opening_dir) >= c.opening_angle;
return angle_between(dir<3, ScalarT>(apexOuterToP), c.opening_dir) <= ScalarT(0.5) * c.opening_angle
&& angle_between(dir<3, ScalarT>(apexInnerToP), c.opening_dir) >= ScalarT(0.5) * c.opening_angle;
}
template <class ScalarT>
[[nodiscard]] constexpr bool contains(pyramid<triangle<3, ScalarT>> const& py, pos<3, ScalarT> const& p, dont_deduce<ScalarT> eps = ScalarT(0))
{
const auto c = centroid(py.base);
auto n = normal(py.base);
if (dot(p - c + eps * n, n) < ScalarT(0))
return false; // Not inside if on the other side of the base
// Check if inside for each pyramid side
const auto apex = c + n * py.height;
n = normalize(cross(apex - py.base.pos0, py.base.pos1 - py.base.pos0));
if (dot(p - apex + eps * n, n) < ScalarT(0))
return false;
n = normalize(cross(apex - py.base.pos1, py.base.pos2 - py.base.pos1));
if (dot(p - apex + eps * n, n) < ScalarT(0))
return false;
n = normalize(cross(apex - py.base.pos2, py.base.pos0 - py.base.pos2));
if (dot(p - apex + eps * n, n) < ScalarT(0))
return false;
return true;
}
template <class ScalarT>
[[nodiscard]] constexpr bool contains(pyramid<box<2, ScalarT, 3>> const& py, pos<3, ScalarT> const& p, dont_deduce<ScalarT> eps = ScalarT(0))
{
const auto c = centroid(py.base);
auto n = normal(py.base);
if (dot(p - c + eps * n, n) < ScalarT(0))
return false; // Not inside if on the other side of the base
// Check if inside for each pyramid side
const auto apex = c + n * py.height;
const auto p0 = py.base[comp<2, ScalarT>(-1, -1)];
const auto p1 = py.base[comp<2, ScalarT>(1, -1)];
const auto p2 = py.base[comp<2, ScalarT>(1, 1)];
const auto p3 = py.base[comp<2, ScalarT>(-1, 1)];
n = normalize(cross(apex - p0, p1 - p0));
if (dot(p - apex + eps * n, n) < ScalarT(0))
return false;
n = normalize(cross(apex - p1, p2 - p1));
if (dot(p - apex + eps * n, n) < ScalarT(0))
return false;
n = normalize(cross(apex - p2, p3 - p2));
if (dot(p - apex + eps * n, n) < ScalarT(0))
return false;
n = normalize(cross(apex - p3, p0 - p3));
if (dot(p - apex + eps * n, n) < ScalarT(0))
return false;
return true;
}
} // namespace tg
......@@ -4,13 +4,17 @@
#include <typed-geometry/types/objects/halfspace.hh>
#include <typed-geometry/types/objects/line.hh>
#include <typed-geometry/types/objects/plane.hh>
#include <typed-geometry/types/objects/quad.hh>
#include <typed-geometry/types/objects/ray.hh>
#include <typed-geometry/types/objects/segment.hh>
#include <typed-geometry/types/objects/sphere.hh>
#include <typed-geometry/types/objects/triangle.hh>
#include <typed-geometry/functions/vector/normalize.hh>
#include <typed-geometry/functions/vector/perpendicular.hh>
#include "typed-geometry/functions/tests/vec_tests.hh"
// Computes the normal at the surface of an object
// Some objects have a fixed normal everywhere, some only at defined positions
// Evaluating the normal at other positions might be undefined or wrong
......@@ -29,6 +33,12 @@ template <int D, class ScalarT>
return h.normal;
}
template <class ScalarT>
[[nodiscard]] constexpr dir<3, ScalarT> normal(sphere<2, ScalarT, 3> const& d)
{
return d.normal;
}
template <class ScalarT>
[[nodiscard]] constexpr dir<2, ScalarT> normal(line<2, ScalarT> const& l)
{
......@@ -53,6 +63,15 @@ template <class ScalarT>
return normalize(cross(t.pos1 - t.pos0, t.pos2 - t.pos0));
}
template <class ScalarT>
[[nodiscard]] constexpr dir<3, fractional_result<ScalarT>> normal(quad<3, ScalarT> const& q)
{
// Assumes the quad is planar, as it is a requirement for pyramid<quad>
const auto res = normalize(cross(q.pos01 - q.pos00, q.pos10 - q.pos00));
TG_ASSERT(tg::are_orthogonal(q.pos11 - q.pos00, res)); // Checks that the four points are indeed coplanar
return res;
}
template <class ScalarT, class TraitsT>
[[nodiscard]] constexpr dir<3, fractional_result<ScalarT>> normal(box<2, ScalarT, 3, TraitsT> const& b)