MVE - Multi-View Environment mve-devel
Loading...
Searching...
No Matches
basis_function.h
Go to the documentation of this file.
1/*
2 * Copyright (C) 2015, Simon Fuhrmann
3 * TU Darmstadt - Graphics, Capture and Massively Parallel Computing
4 * All rights reserved.
5 *
6 * This software may be modified and distributed under the terms
7 * of the BSD 3-Clause license. See the LICENSE.txt file for details.
8 */
9
10#ifndef FSSR_BASIS_FUNCTION_HEADER
11#define FSSR_BASIS_FUNCTION_HEADER
12
13#include "math/vector.h"
14#include "math/matrix.h"
15#include "fssr/defines.h"
16#include "fssr/sample.h"
17
19
20/* ---------------------- Gaussian functions ---------------------- */
21
23template <typename T>
24T
25gaussian (T const& sigma, math::Vector<T, 3> const& x);
26
28template <typename T>
29T
30gaussian_normalized (T const& sigma, math::Vector<T, 3> const& x);
31
32/* ---------------------- FSSR basis function --------------------- */
33
42template <typename T>
43T
44fssr_basis (T const& scale, math::Vector<T, 3> const& pos,
45 math::Vector<T, 3>* deriv = nullptr);
46
47/* --------------------- FSSR weight function --------------------- */
48
55template <typename T>
56T
57fssr_weight (T const& scale, math::Vector<T, 3> const& pos,
58 math::Vector<T, 3>* deriv = nullptr);
59
60/* -------------------------- Helper functions --------------------------- */
61
67void
68evaluate (math::Vec3f const& pos, Sample const& sample,
69 double* value, double* weight,
70 math::Vector<double, 3>* value_deriv,
71 math::Vector<double, 3>* weight_deriv);
72
75transform_position (math::Vec3f const& pos, Sample const& sample);
76
78void
80
82void
84
86
87/*
88 * ========================= Implementation ============================
89 */
90
92
93template <typename T>
94inline T
95gaussian (T const& sigma, math::Vector<T, 3> const& x)
96{
97 return std::exp(-x.dot(x) / (T(2) * MATH_POW2(sigma)));
98}
99
100template <typename T>
101inline T
102gaussian_normalized (T const& sigma, math::Vector<T, 3> const& x)
103{
104 return gaussian(sigma, x) / (sigma * MATH_SQRT_2PI);
105}
106
107/* -------------------------------------------------------------------- */
108
109template <typename T>
110inline T
111fssr_basis (T const& scale, math::Vector<T, 3> const& pos,
112 math::Vector<T, 3>* deriv)
113{
114 double const gaussian_value = gaussian(scale, pos);
115 double const value_norm = T(2) * MATH_PI * MATH_POW4(scale);
116
117 if (deriv != nullptr)
118 {
119 double const deriv_norm = value_norm * T(2) * MATH_POW2(scale);
120 (*deriv)[0] = T(2) * (MATH_POW2(scale) - MATH_POW2(pos[0]))
121 * gaussian_value / deriv_norm;
122 (*deriv)[1] = (-T(2) * pos[0] * pos[1]) * gaussian_value / deriv_norm;
123 (*deriv)[2] = (-T(2) * pos[0] * pos[2]) * gaussian_value / deriv_norm;
124 }
125
126 return pos[0] * gaussian_value / value_norm;
127}
128
129/* -------------------------------------------------------------------- */
130
131#if FSSR_NEW_WEIGHT_FUNCTION
132
133template <typename T>
134T
135fssr_weight (T const& scale, math::Vector<T, 3> const& pos,
136 math::Vector<T, 3>* deriv)
137{
138 T const square_radius = pos.square_norm() / MATH_POW2(scale);
139 if (square_radius >= T(9))
140 {
141 if (deriv != nullptr)
142 deriv->fill(T(0));
143 return T(0);
144 }
145
146 if (deriv != nullptr)
147 {
148 T const deriv_factor = -T(4) / T(3)
149 + T(48) / T(54) * std::sqrt(square_radius)
150 - T(4) / T(27) * square_radius;
151
152 (*deriv)[0] = deriv_factor * pos[0] / scale;
153 (*deriv)[1] = deriv_factor * pos[1] / scale;
154 (*deriv)[2] = deriv_factor * pos[2] / scale;
155 }
156
157 /*
158 * w(r > 0) = 1.0 - 2.0/3.0 * r**2 + 8.0/27.0 * r**3 - 1.0/27.0 * r**4
159 * w(r < 0) = 1.0 - 2.0/3.0 * r**2 - 8.0/27.0 * r**3 - 1.0/27.0 * r**4
160 */
161 return T(1) - T(2) / T(3) * square_radius
162 + T(8) / T(27) * std::pow(square_radius, T(1.5))
163 - T(1) / T(27) * MATH_POW2(square_radius);
164}
165
166#else
167
168template <typename T>
169T
170fssr_weight (T const& scale, math::Vector<T, 3> const& pos,
171 math::Vector<T, 3>* deriv)
172{
173 T const x = pos[0] / scale;
174 T const y = pos[1] / scale;
175 T const z = pos[2] / scale;
176 T const square_radius = MATH_POW2(y) + MATH_POW2(z);
177
178 /*
179 * wx(x > 0) = 1.0 - 1.0/3.0 * r**2 + 2.0/27.0 * r**3
180 * wx(x < 0) = 1.0 + 2.0/3.0 * r + 1.0/9.0 * r**2
181 */
182 T weight_x = T(0);
183 if (x > T(-3) && x < T(0))
184 {
185 weight_x = T(1)
186 + T(2) / T(3) * x
187 + T(1) / T(9) * MATH_POW2(x);
188 }
189 else if (x >= T(0) && x < T(3))
190 {
191 weight_x = T(1)
192 - T(1) / T(3) * MATH_POW2(x)
193 + T(2) / T(27) * MATH_POW3(x);
194 }
195
196 /*
197 * wyz(r) = 1.0 - 1.0/3.0 r**2 / s**2 + 2.0/27.0 r**3 / s**3
198 */
199 T weight_yz = T(0);
200 if (square_radius < T(9))
201 {
202 weight_yz = T(1)
203 - T(1) / T(3) * square_radius
204 + T(2) / T(27) * std::pow(square_radius, T(1.5));
205 }
206
207 if (deriv != nullptr)
208 {
209 /*
210 * wx'(x < 0) = 2.0/9.0 * x / s^2 + 2.0/3.0 / s
211 * wx'(x > 0) = 6.0/27.0 * x^2 / s^3 - 2.0/3.0 * x / s^2
212 */
213 T deriv_x = T(0);
214 if (x > T(-3) && x <= T(0))
215 {
216 deriv_x = (T(2) / T(3) + T(2) / T(9) * x) / scale;
217 }
218 else if (x > T(0) && x < T(3))
219 {
220 deriv_x = (-T(2) / T(3) * x + T(6) / T(27) * MATH_POW2(x)) / scale;
221 }
222
223 /*
224 * w(y,z) = 2/27 (y^2 + z^2)^(3/2) / s^3 - 1/3 (y^2 + z^2) / s^2 + 1
225 * d/dy w'(y,z) = y/s * (12 / (54 s^2) * (y^2 + z^2)^(1/2) - 2 / (3s))
226 * d/dz w'(y,z) = z/s * (12 / (54 s^2) * (y^2 + z^2)^(1/2) - 2 / (3s))
227 */
228 T deriv_y = T(0);
229 T deriv_z = T(0);
230 if (square_radius < T(9))
231 {
232 T const factor = (T(12) / (T(54) * scale)
233 * std::sqrt(MATH_POW2(pos[1]) + MATH_POW2(pos[2]))
234 - T(2) / T(3)) / MATH_POW2(scale);
235 deriv_y = factor * pos[1];
236 deriv_z = factor * pos[2];
237 }
238
239 (*deriv)[0] = deriv_x * weight_yz;
240 (*deriv)[1] = deriv_y * weight_x;
241 (*deriv)[2] = deriv_z * weight_x;
242 }
243
244 return weight_x * weight_yz;
245}
246
247#endif
248
249/* -------------------------------------------------------------------- */
250
251inline math::Vec3f
252transform_position (math::Vec3f const& pos, Sample const& sample)
253{
254 math::Matrix3f rot;
255 rotation_from_normal(sample.normal, &rot);
256 return rot * (pos - sample.pos);
257}
258
259/* -------------------------------------------------------------------- */
260
261#if 0 /* Following is dead but potentially still useful code. */
267template <typename T>
268inline T
269linear_ramp (Sample const& sample, math::Vector<T, 3> const& pos)
270{
271 return (pos - sample.pos).dot(sample.normal);
272}
273
278template <typename T>
279T
280weighting_function_mpu (T const& x)
281{
282 if (x <= T(-3) || x >= T(3))
283 return T(0);
284
285 T xf = (x + T(3)) / T(2);
286 if (xf <= T(1))
287 return MATH_POW2(xf) / T(2);
288 if (xf <= T(2))
289 return (T(-2) * MATH_POW2(xf) + T(6) * xf - T(3)) / T(2);
290 return MATH_POW2(T(3) - xf) / T(2);
291}
292
298template <typename T>
299inline T
300weighting_function_mpu (T const& sample_scale, math::Vector<T, 3> const& pos)
301{
302 return weighting_function_mpu(pos.norm() / sample_scale);
303}
304#endif
305
307
308#endif // FSSR_BASIS_FUNCTION_HEADER
Matrix class for arbitrary dimensions and types.
Definition matrix.h:54
Vector class for arbitrary dimensions and types.
Definition vector.h:87
T square_norm(void) const
Computes the squared norm of the vector (much cheaper).
Definition vector.h:441
T norm(void) const
Computes the norm (length) of the vector.
Definition vector.h:434
T dot(Vector< T, N > const &other) const
Dot (or scalar) product between self and another vector.
Definition vector.h:542
Vector< T, N > & fill(T const &value)
Fills all vector elements with the given value.
Definition vector.h:381
#define FSSR_NAMESPACE_END
Definition defines.h:14
#define FSSR_NAMESPACE_BEGIN
Definition defines.h:13
#define MATH_POW3(x)
Definition defines.h:69
#define MATH_POW4(x)
Definition defines.h:70
#define MATH_SQRT_2PI
Definition defines.h:54
#define MATH_PI
Definition defines.h:42
#define MATH_POW2(x)
Definition defines.h:68
T gaussian(T const &sigma, math::Vector< T, 3 > const &x)
The Gaussian function in 3D.
T fssr_weight(T const &scale, math::Vector< T, 3 > const &pos, math::Vector< T, 3 > *deriv=nullptr)
Evaluates the FSSR weight function and directional derivatives (optional).
T gaussian_normalized(T const &sigma, math::Vector< T, 3 > const &x)
The normalized Gaussian function in 3D.
math::Vec3f transform_position(math::Vec3f const &pos, Sample const &sample)
Transforms 'pos' according to the samples position and normal.
void rotation_from_normal(math::Vec3f const &normal, math::Matrix3f *rot)
Generates a rotation matrix that transforms in the FSSR LCS.
T fssr_basis(T const &scale, math::Vector< T, 3 > const &pos, math::Vector< T, 3 > *deriv=nullptr)
Evaluates the FSSR basis function and directional derivatives (optional).
void evaluate(math::Vec3f const &pos, Sample const &sample, double *value, double *weight, math::Vector< double, 3 > *value_deriv, math::Vector< double, 3 > *weight_deriv)
Rotates the given point in the LCS of the sample, evaluates the basis and weight functions and their ...
Representation of a point sample.
Definition sample.h:22
math::Vec3f pos
Definition sample.h:23
math::Vec3f normal
Definition sample.h:24