Nabla/include/nbl/builtin/hlsl/shapes/spherical_triangle.hlsl at 62c79b40e00edcda8a3599c149198c2fb14fe88c · Devsh-Graphics-Programming/Nabla · GitHub

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
106
107
108
109
// Copyright (C) 2018-2023 - DevSH Graphics Programming Sp. z O.O.
// This file is part of the "Nabla Engine".
// For conditions of distribution and use, see copyright notice in nabla.h

#ifndef _NBL_BUILTIN_HLSL_SHAPES_SPHERICAL_TRIANGLE_INCLUDED_
#define _NBL_BUILTIN_HLSL_SHAPES_SPHERICAL_TRIANGLE_INCLUDED_

#include <nbl/builtin/hlsl/tgmath.hlsl>
#include <nbl/builtin/hlsl/cpp_compat.hlsl>
#include <nbl/builtin/hlsl/limits.hlsl>
#include <nbl/builtin/hlsl/math/functions.hlsl>
#include <nbl/builtin/hlsl/math/angle_adding.hlsl>
#include <nbl/builtin/hlsl/numbers.hlsl>

namespace nbl
{
namespace hlsl
{
namespace shapes
{

template<typename T>
struct SphericalTriangle
{
    using scalar_type = T;
    using vector3_type = vector<T, 3>;

    static SphericalTriangle<T> create(const vector3_type vertex0, const vector3_type vertex1, const vector3_type vertex2, const vector3_type origin)
    {
        SphericalTriangle<T> retval;
        retval.vertex0 = nbl::hlsl::normalize(vertex0 - origin);
        retval.vertex1 = nbl::hlsl::normalize(vertex1 - origin);
        retval.vertex2 = nbl::hlsl::normalize(vertex2 - origin);
        retval.cos_sides = vector3_type(hlsl::dot(retval.vertex1, retval.vertex2), hlsl::dot(retval.vertex2, retval.vertex0), hlsl::dot(retval.vertex0, retval.vertex1));
        const vector3_type csc_sides2 = hlsl::promote<vector3_type>(1.0) - retval.cos_sides * retval.cos_sides;
        retval.csc_sides.x = hlsl::rsqrt<scalar_type>(csc_sides2.x);
        retval.csc_sides.y = hlsl::rsqrt<scalar_type>(csc_sides2.y);
        retval.csc_sides.z = hlsl::rsqrt<scalar_type>(csc_sides2.z);
        return retval;
    }

    bool pyramidAngles()
    {
        return hlsl::any<vector<bool, 3> >(csc_sides >= (vector3_type)(numeric_limits<scalar_type>::max));
    }

    scalar_type solidAngleOfTriangle(NBL_REF_ARG(vector3_type) cos_vertices, NBL_REF_ARG(vector3_type) sin_vertices, NBL_REF_ARG(scalar_type) cos_a, NBL_REF_ARG(scalar_type) cos_c, NBL_REF_ARG(scalar_type) csc_b, NBL_REF_ARG(scalar_type) csc_c)
    {
        if (pyramidAngles())
            return 0.f;

        // these variables might eventually get optimized out
        cos_a = cos_sides[0];
        cos_c = cos_sides[2];
        csc_b = csc_sides[1];
        csc_c = csc_sides[2];

        // Both vertices and angles at the vertices are denoted by the same upper case letters A, B, and C. The angles A, B, C of the triangle are equal to the angles between the planes that intersect the surface of the sphere or, equivalently, the angles between the tangent vectors of the great circle arcs where they meet at the vertices. Angles are in radians. The angles of proper spherical triangles are (by convention) less than PI
        cos_vertices = hlsl::clamp((cos_sides - cos_sides.yzx * cos_sides.zxy) * csc_sides.yzx * csc_sides.zxy, hlsl::promote<vector3_type>(-1.0), hlsl::promote<vector3_type>(1.0)); // using Spherical Law of Cosines (TODO: do we need to clamp anymore? since the pyramid angles method introduction?)
        sin_vertices = hlsl::sqrt(hlsl::promote<vector3_type>(1.0) - cos_vertices * cos_vertices);

        math::sincos_accumulator<scalar_type> angle_adder = math::sincos_accumulator<scalar_type>::create(cos_vertices[0], sin_vertices[0]);
        angle_adder.addAngle(cos_vertices[1], sin_vertices[1]);
        angle_adder.addAngle(cos_vertices[2], sin_vertices[2]);
        return angle_adder.getSumofArccos() - numbers::pi<scalar_type>;
    }

    scalar_type solidAngleOfTriangle()
    {
        vector3_type dummy0,dummy1;
        scalar_type dummy2,dummy3,dummy4,dummy5;
        return solidAngleOfTriangle(dummy0,dummy1,dummy2,dummy3,dummy4,dummy5);
    }

    scalar_type projectedSolidAngleOfTriangle(const vector3_type receiverNormal, NBL_REF_ARG(vector3_type) cos_sides, NBL_REF_ARG(vector3_type) csc_sides, NBL_REF_ARG(vector3_type) cos_vertices)
    {
        if (pyramidAngles())
            return 0.f;

        vector3_type awayFromEdgePlane0 = hlsl::cross<vector3_type>(vertex1, vertex2) * csc_sides[0];
        vector3_type awayFromEdgePlane1 = hlsl::cross<vector3_type>(vertex2, vertex0) * csc_sides[1];
        vector3_type awayFromEdgePlane2 = hlsl::cross<vector3_type>(vertex0, vertex1) * csc_sides[2];

        // useless here but could be useful somewhere else
        cos_vertices[0] = hlsl::dot<vector3_type>(awayFromEdgePlane1, awayFromEdgePlane2);
        cos_vertices[1] = hlsl::dot<vector3_type>(awayFromEdgePlane2, awayFromEdgePlane0);
        cos_vertices[2] = hlsl::dot<vector3_type>(awayFromEdgePlane0, awayFromEdgePlane1);
        // TODO: above dot products are in the wrong order, either work out which is which, or try all 6 permutations till it works
        cos_vertices = hlsl::clamp<vector3_type>((cos_sides - cos_sides.yzx * cos_sides.zxy) * csc_sides.yzx * csc_sides.zxy, hlsl::promote<vector3_type>(-1.0), hlsl::promote<vector3_type>(1.0));

        matrix<scalar_type, 3, 3> awayFromEdgePlane = matrix<scalar_type, 3, 3>(awayFromEdgePlane0, awayFromEdgePlane1, awayFromEdgePlane2);
        const vector3_type externalProducts = hlsl::abs(hlsl::mul(/* transposed already */awayFromEdgePlane, receiverNormal));

        const vector3_type pyramidAngles = acos<scalar_type>(cos_sides);
        return hlsl::dot<vector3_type>(pyramidAngles, externalProducts) / (2.f * numbers::pi<scalar_type>);
    }

    vector3_type vertex0;
    vector3_type vertex1;
    vector3_type vertex2;
    vector3_type cos_sides;
    vector3_type csc_sides;
};

}
}
}

#endif