1 
2 //
3 // This source file is part of appleseed.
4 // Visit https://appleseedhq.net/ for additional information and resources.
5 //
6 // This software is released under the MIT license.
7 //
8 // Copyright (c) 2019 Esteban Tovagliari, The appleseedhq Organization
9 //
10 // Permission is hereby granted, free of charge, to any person obtaining a copy
11 // of this software and associated documentation files (the "Software"), to deal
12 // in the Software without restriction, including without limitation the rights
13 // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
14 // copies of the Software, and to permit persons to whom the Software is
15 // furnished to do so, subject to the following conditions:
16 //
17 // The above copyright notice and this permission notice shall be included in
18 // all copies or substantial portions of the Software.
19 //
20 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
21 // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
23 // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
25 // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
26 // THE SOFTWARE.
27 //
28 
29 #pragma once
30 
31 // appleseed.foundation headers.
32 #include "foundation/math/scalar.h"
33 #include "foundation/math/vector.h"
34 
35 // Standard headers.
36 #include <cassert>
37 #include <cmath>
38 
39 namespace foundation
40 {
41 
42 //
43 // References:
44 //
45 //   An Area-Preserving Parametrization for Spherical Rectangles.
46 //   https://www.arnoldrenderer.com/research/egsr2013_spherical_rectangle.pdf
47 //
48 
49 template <typename T>
50 class SphericalRectangleSampler
51 {
52   public:
SphericalRectangleSampler(const Vector<T,3> & origin,const Vector<T,3> & x_axis,const Vector<T,3> & y_axis,const Vector<T,3> & normal,const Vector<T,3> & center)53     SphericalRectangleSampler(
54         const Vector<T, 3>& origin,
55         const Vector<T, 3>& x_axis,
56         const Vector<T, 3>& y_axis,
57         const Vector<T, 3>& normal,
58         const Vector<T, 3>& center)
59       : m_x(x_axis)
60       , m_y(y_axis)
61       , m_z(normal)
62       , m_center(center)
63     {
64         const Vector<T, 3> d = origin - center;
65 
66         const T width = norm(x_axis);
67         const T height = norm(y_axis);
68 
69         const T x0 = dot(d, x_axis);
70         const T x1 = x0 + width;
71 
72         m_y0 = dot(d, y_axis);
73         m_y1 = m_y0 + height;
74         m_z0 = dot(d, normal);
75 
76         // z flip
77         if (m_z0 > T(0.0))
78         {
79             m_z0 *= T(-1.0);
80             m_z  *= T(-1.0);
81         }
82 
83         m_n0 = Vector<T, 3>(T(0.0), m_z0, -m_y0);
84         m_n2 = Vector<T, 3>(T(0.0), -m_z0, m_y1);
85         const Vector<T, 3> n1 = Vector<T, 3>(-m_z0, T(0.0), x1);
86         const Vector<T, 3> n3 = Vector<T, 3>(m_z0, T(0.0), -x0);
87 
88         m_y0y0 = square(m_y0);
89         m_y1y1 = square(m_y1);
90         m_z0z0 = square(m_z0);
91 
92         m_n0.z /= std::sqrt(m_z0z0 + m_y0y0);
93         m_n2.z /= std::sqrt(m_z0z0 + m_y1y1);
94         n1.z /= std::sqrt(m_z0z0 + square(x1));
95         n3.z /= std::sqrt(m_z0z0 + square(x0));
96 
97         const T g0 = std::acos(-m_n0.z * n1.z);
98         const T g1 = std::acos(-n1.z * m_n2.z);
99         m_g2 = std::acos(-m_n2.z * n3.z);
100         m_g3 = std::acos(-n3.z * m_n0.z);
101 
102         m_area = g0 + g1 + m_g2 + m_g3 - TwoPi<T>();
103     }
104 
get_area()105     T get_area() const
106     {
107         return m_area;
108     }
109 
sample(const Vector<T,2> & s)110     Vector<T, 3> sample(const Vector<T, 2>& s) const
111     {
112         const T au = s[0] * m_area - m_g2 - m_g3 + TwoPi<T>();
113 
114         const T b0 = m_n0.z;
115         const T b1 = m_n2.z;
116         const T fu = (std::cos(au) * b0 - b1) / std::sin(au);
117 
118         const T cu =
119             clamp(
120                 std::copysign(T(1.0), fu) / std::sqrt(square(fu) + square(b0)),
121                 T(-1.0),
122                 T(1.0));
123 
124         const T xu = (-cu * m_z0) / safe_sqrt(T(1.0) - square(cu));
125         const T d = std::sqrt(square(xu) + m_z0z0);
126         const T d2 = square(d);
127 
128         const T h0 = m_y0 / std::sqrt(d2 + m_y0y0);
129         const T h1 = m_y1 / std::sqrt(d2 + m_y1y1);
130 
131         const T hv = lerp(h0, h1, s[1]);
132         const T hv2 = square(hv);
133         const T yv = (hv * d) / safe_sqrt(T(1.0) - hv2);
134 
135         return m_center + xu * m_x + yv * m_y + m_z0 * m_z;
136     }
137 
138   private:
139     const Vector<T, 3>  m_x;
140     const Vector<T, 3>  m_y;
141     Vector<T, 3>        m_z;
142     const Vector<T, 3>  m_center;
143 
144     T                   m_y0;
145     T                   m_y1;
146     T                   m_z0;
147     T                   m_z0z0;
148     T                   m_y0y0;
149     T                   m_y1y1;
150 
151     T                   m_area;
152 
153     Vector<T, 3>        m_n0;
154     Vector<T, 3>        m_n2;
155     T                   m_g2;
156     T                   m_g3;
157 };
158 
159 }  // namespace foundation
160