Mitsuba Renderer  0.5.0
shvector.h
Go to the documentation of this file.
1 /*
2  This file is part of Mitsuba, a physically based rendering system.
3
4  Copyright (c) 2007-2014 by Wenzel Jakob and others.
5
6  Mitsuba is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License Version 3
9
10  Mitsuba is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU General Public License for more details.
14
15  You should have received a copy of the GNU General Public License
16  along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19 #pragma once
20 #if !defined(__MITSUBA_CORE_SHVECTOR_H_)
21 #define __MITSUBA_CORE_SHVECTOR_H_
22
23 #include <mitsuba/mitsuba.h>
25 #include <Eigen/Core>
26
28
29 /* Precompute normalization coefficients for the first 10 bands */
30 #define SH_NORMTBL_SIZE 10
31
32 struct SHVector;
33
34 /**
35  * \brief Stores the diagonal blocks of a spherical harmonic
36  * rotation matrix
37  *
38  * \ingroup libcore
39  * \ingroup libpython
40  */
42  typedef Eigen::Matrix<Float, Eigen::Dynamic, Eigen::Dynamic> Matrix;
43
44  std::vector<Matrix> blocks;
45
46  /// Construct a new rotation storage for the given number of bands
47  inline SHRotation(int bands) : blocks(bands) {
48  for (int i=0; i<bands; ++i) {
49  int dim = 2*i+1;
50  blocks[i] = Matrix(dim, dim);
51  }
52  }
53
54  /**
55  * \brief Transform a coefficient vector and store the result into
56  * the given target vector.
57  *
58  * The source and target must have the same number of bands.
59  */
60  void operator()(const SHVector &source, SHVector &target) const;
61 };
62
63 /**
64  * \brief Stores a truncated real spherical harmonics representation of
65  * an L2-integrable function.
66  *
67  * Also provides some other useful functionality, such as evaluation,
68  * projection and rotation.
69  *
70  * The Mathematica equivalent of the basis functions implemented here is:
71  *
72  * \code
73  * SphericalHarmonicQ[l_, m_, \[Theta]_, \[Phi]_] :=
74  * Piecewise[{
75  * {SphericalHarmonicY[l, m, \[Theta], \[Phi]], m == 0},
76  * {Sqrt[2]*Re[SphericalHarmonicY[l, m, \[Theta], \[Phi]]], m > 0},
77  * {Sqrt[2]*Im[SphericalHarmonicY[l, -m, \[Theta], \[Phi]]], m < 0}
78  * }]
79  * \endcode
80  *
81  * \ingroup libcore
82  * \ingroup libpython
83  */
85 public:
86  /// Construct an invalid SH vector
87  inline SHVector()
88  : m_bands(0) {
89  }
90
91  /// Construct a new SH vector (initialized to zero)
92  inline SHVector(int bands)
93  : m_bands(bands), m_coeffs(bands*bands) {
94  clear();
95  }
96
97  /// Unserialize a SH vector to a binary data stream
98  SHVector(Stream *stream);
99
100  /// Copy constructor
101  inline SHVector(const SHVector &v) : m_bands(v.m_bands),
102  m_coeffs(v.m_coeffs) {
103  }
104
105  /// Return the number of stored SH coefficient bands
106  inline int getBands() const {
107  return m_bands;
108  }
109
110  /// Serialize a SH vector to a binary data stream
111  void serialize(Stream *stream) const;
112
113  /// Get the energy per band
114  inline Float energy(int band) const {
115  Float result = 0;
116  for (int m=-band; m<=band; ++m)
117  result += std::abs(operator()(band,m));
118  return result;
119  }
120
121  /// Assignment
122  inline SHVector &operator=(const SHVector &v) {
123  m_bands = v.m_bands;
124  m_coeffs = v.m_coeffs;
125  return *this;
126  }
127
128  /// Set all coefficients to zero
129  inline void clear() {
130  m_coeffs.setZero();
131  }
132
134  inline SHVector& operator+=(const SHVector &v) {
135  ptrdiff_t extendBy = v.m_coeffs.size() - m_coeffs.size();
136  if (extendBy > 0) {
137  m_coeffs.conservativeResize(v.m_coeffs.rows());
138  m_coeffs.tail(extendBy).setZero();
139  m_bands = v.m_bands;
140  }
142  return *this;
143  }
144
146  inline SHVector operator+(const SHVector &v) const {
147  SHVector vec(std::max(m_bands, v.m_bands));
148  if (m_bands > v.m_bands) {
149  vec.m_coeffs = m_coeffs;
151  } else {
152  vec.m_coeffs = v.m_coeffs;
154  }
155  return vec;
156  }
157
158  /// Component-wise subtraction
159  inline SHVector& operator-=(const SHVector &v) {
160  ptrdiff_t extendBy = v.m_coeffs.size() - m_coeffs.size();
161  if (extendBy > 0) {
162  m_coeffs.conservativeResize(v.m_coeffs.rows());
163  m_coeffs.tail(extendBy).setZero();
164  m_bands = v.m_bands;
165  }
167  return *this;
168  }
169
170  /// Component-wise subtraction
171  inline SHVector operator-(const SHVector &v) const {
172  SHVector vec(std::max(m_bands, v.m_bands));
173  if (m_bands > v.m_bands) {
174  vec.m_coeffs = m_coeffs;
176  } else {
177  vec.m_coeffs = -v.m_coeffs;
179  }
180  return vec;
181  }
182
183  /// Add a scalar multiple of another vector
184  inline SHVector& madd(Float f, const SHVector &v) {
185  ptrdiff_t extendBy = v.m_coeffs.size() - m_coeffs.size();
186  if (extendBy > 0) {
187  m_coeffs.conservativeResize(v.m_coeffs.rows());
188  m_coeffs.tail(extendBy).setZero();
189  m_bands = v.m_bands;
190  }
192
193  return *this;
194  }
195
196  /// Scalar multiplication
197  inline SHVector &operator*=(Float f) {
198  m_coeffs *= f;
199  return *this;
200  }
201
202  /// Scalar multiplication
203  inline SHVector operator*(Float f) const {
204  SHVector vec(m_bands);
205  vec.m_coeffs = m_coeffs * f;
206  return vec;
207  }
208
209  /// Scalar division
210  inline SHVector &operator/=(Float f) {
211  m_coeffs *= (Float) 1 / f;
212  return *this;
213  }
214
215  /// Scalar division
216  inline SHVector operator/(Float f) const {
217  SHVector vec(m_bands);
218  vec.m_coeffs = m_coeffs * (1/f);
219  return vec;
220  }
221
222  /// Negation operator
223  inline SHVector operator-() const {
224  SHVector vec(m_bands);
225  vec.m_coeffs = -m_coeffs;
226  return vec;
227  }
228
229  /// Access coefficient m (in {-l, ..., l}) on band l
230  inline Float &operator()(int l, int m) {
231  return m_coeffs[l*(l+1) + m];
232  }
233
234  /// Access coefficient m (in {-l, ..., l}) on band l
235  inline const Float &operator()(int l, int m) const {
236  return m_coeffs[l*(l+1) + m];
237  }
238
239  /// Evaluate for a direction given in spherical coordinates
240  Float eval(Float theta, Float phi) const;
241
242  /// Evaluate for a direction given in Cartesian coordinates
243  Float eval(const Vector &v) const;
244
245  /**
246  * \brief Evaluate for a direction given in spherical coordinates.
247  *
248  * This function is much faster but only works for azimuthally
249  * invariant functions
250  */
251  Float evalAzimuthallyInvariant(Float theta, Float phi) const;
252
253  /**
254  * \brief Evaluate for a direction given in cartesian coordinates.
255  *
256  * This function is much faster but only works for azimuthally
257  * invariant functions
258  */
259  Float evalAzimuthallyInvariant(const Vector &v) const;
260
261  /// Check if this function is azumuthally invariant
262  bool isAzimuthallyInvariant() const;
263
264  /// Equality comparison operator
265  inline bool operator==(const SHVector &v) const {
266  return m_bands == v.m_bands && m_coeffs == v.m_coeffs;
267  }
268
269  /// Equality comparison operator
270  inline bool operator!=(const SHVector &v) const {
271  return !operator==(v);
272  }
273
274  /// Dot product
275  inline friend Float dot(const SHVector &v1, const SHVector &v2);
276
277  /// Normalize so that the represented function becomes a valid distribution
278  void normalize();
279
280  /// Compute the second spherical moment (analytic)
281  Matrix3x3 mu2() const;
282
283  /// Brute-force search for the minimum value over the sphere
284  Float findMinimum(int res) const;
285
286  /// Add a constant value
288
289  /**
290  * \brief Convolve the SH representation with the supplied kernel.
291  *
292  * Based on the Funk-Hecke theorem -- the kernel must be rotationally
293  * symmetric around the Z-axis.
294  */
295  void convolve(const SHVector &kernel);
296
297  /// Project the given function onto a SH basis (using a 2D composite Simpson's rule)
298  template<typename Functor> void project(const Functor &f, int res = 32) {
299  SAssert(res % 2 == 0);
300  /* Nested composite Simpson's rule */
301  Float hExt = M_PI / res,
302  hInt = (2*M_PI)/(res*2);
303
304  for (int l=0; l<m_bands; ++l)
305  for (int m=-l; m<=l; ++m)
306  operator()(l,m) = 0;
307
308  Float *sinPhi = (Float *) alloca(sizeof(Float)*m_bands),
309  *cosPhi = (Float *) alloca(sizeof(Float)*m_bands);
310
311  for (int i=0; i<=res; ++i) {
312  Float theta = hExt*i, cosTheta = std::cos(theta);
313  Float weightExt = (i & 1) ? 4.0f : 2.0f;
314  if (i == 0 || i == res)
315  weightExt = 1.0f;
316
317  for (int j=0; j<=res*2; ++j) {
318  Float phi = hInt*j;
319  Float weightInt = (j & 1) ? 4.0f : 2.0f;
320  if (j == 0 || j == 2*res)
321  weightInt = 1.0f;
322
323  for (int m=0; m<m_bands; ++m) {
324  sinPhi[m] = std::sin((m+1)*phi);
325  cosPhi[m] = std::cos((m+1)*phi);
326  }
327
328  Float value = f(sphericalDirection(theta, phi))*std::sin(theta)
329  * weightExt*weightInt;
330
331  for (int l=0; l<m_bands; ++l) {
332  for (int m=1; m<=l; ++m) {
333  Float L = legendreP(l, m, cosTheta) * normalization(l, m);
334  operator()(l, -m) += value * SQRT_TWO * sinPhi[m-1] * L;
335  operator()(l, m) += value * SQRT_TWO * cosPhi[m-1] * L;
336  }
337
338  operator()(l, 0) += value * legendreP(l, 0, cosTheta) * normalization(l, 0);
339  }
340  }
341  }
342
343  for (int l=0; l<m_bands; ++l)
344  for (int m=-l; m<=l; ++m)
345  operator()(l,m) *= hExt*hInt/9;
346  }
347
348  /// Compute the relative L2 error
349  template<typename Functor> Float l2Error(const Functor &f, int res = 32) const {
350  SAssert(res % 2 == 0);
351  /* Nested composite Simpson's rule */
352  Float hExt = M_PI / res,
353  hInt = (2*M_PI)/(res*2);
354  Float error = 0.0f, denom=0.0f;
355
356  for (int i=0; i<=res; ++i) {
357  Float theta = hExt*i;
358  Float weightExt = (i & 1) ? 4.0f : 2.0f;
359  if (i == 0 || i == res)
360  weightExt = 1.0f;
361
362  for (int j=0; j<=res*2; ++j) {
363  Float phi = hInt*j;
364  Float weightInt = (j & 1) ? 4.0f : 2.0f;
365  if (j == 0 || j == 2*res)
366  weightInt = 1.0f;
367
368  Float value1 = f(sphericalDirection(theta, phi));
369  Float value2 = eval(theta, phi);
370  Float diff = value1-value2;
371  Float weight = std::sin(theta)*weightInt*weightExt;
372
373  error += diff*diff*weight;
374  denom += value1*value1*weight;
375  }
376  }
377
378  return error/denom;
379  }
380
381  /// Turn into a string representation
382  std::string toString() const;
383
384  /// Return a normalization coefficient
385  inline static Float normalization(int l, int m) {
386  if (l < SH_NORMTBL_SIZE)
387  return m_normalization[l*(l+1)/2 + m];
388  else
389  return computeNormalization(l, m);
390  }
391
392  /**
393  * \brief Recursively computes rotation matrices for each band of SH coefficients.
394  *
395  * Based on 'Rotation Matrices for Real Spherical Harmonics. Direct Determination by Recursion'
396  * by Ivanic and Ruedenberg. The implemented tables follow the notation in
397  * 'Spherical Harmonic Lighting: The Gritty Details' by Robin Green.
398  */
399  static void rotation(const Transform &t, SHRotation &rot);
400
401  /// Precomputes normalization coefficients for the first few bands
402  static void staticInitialization();
403
404  /// Free the memory taken up by staticInitialization()
405  static void staticShutdown();
406 protected:
407  /// Helper function for rotation() -- computes a diagonal block based on the previous level
408  static void rotationBlock(const SHRotation::Matrix &M1, const SHRotation::Matrix &Mp, SHRotation::Matrix &Mn);
409
410  /// Compute a normalization coefficient
411  static Float computeNormalization(int l, int m);
412 private:
413  int m_bands;
414  Eigen::Matrix<Float, Eigen::Dynamic, 1> m_coeffs;
415  static Float *m_normalization;
416 };
417
418 inline Float dot(const SHVector &v1, const SHVector &v2) {
419  const size_t size = std::min(v1.m_coeffs.size(), v2.m_coeffs.size());
421 }
422
423 /**
424  * \brief Implementation of 'Importance Sampling Spherical Harmonics'
425  * by W. Jarsz, N. Carr and H. W. Jensen (EUROGRAPHICS 2009)
426  *
427  * \ingroup libcore
428  * \ingroup libpython
429  */
431 public:
432  /**
433  * \brief Precompute a spherical harmonics sampler object
434  *
435  * \param bands Number of SH coefficient bands to support
436  * \param depth Number of recursive sample warping steps.
437  */
438  SHSampler(int bands, int depth);
439
440  /**
441  * \brief Warp a uniform sample in [0,1]^2 to one that is
442  * approximately proportional to the specified function.
443  *
444  * The resulting sample will have spherical coordinates
445  * [0,pi]x[0,2pi] and its actual PDF (which might be
446  * slightly different from the function evaluated at the
447  * sample, even if $f$ is a distribution) will be returned.
448  */
449  Float warp(const SHVector &f, Point2 &sample) const;
450
451  /// Return information on the size of the precomputed tables
452  std::string toString() const;
453
455 protected:
456  /// Virtual destructor
457  virtual ~SHSampler();
458
459  /* Index into the assoc. legendre polynomial table */
460  inline int I(int l, int m) const { return l*(l+1)/2 + m; }
461
462  /* Index into the phi table */
463  inline int P(int m) const { return m + m_bands; }
464
465  inline Float lookupIntegral(int depth, int zBlock, int phiBlock, int l, int m) const {
466  return -m_phiMap[depth][phiBlock][P(m)] * m_legendreMap[depth][zBlock][I(l, std::abs(m))];
467  }
468
469  /// Recursively compute assoc. legendre & phi integrals
470  Float *legendreIntegrals(Float a, Float b);
471  Float *phiIntegrals(Float a, Float b);
472
473  /// Integrate a SH expansion over the specified mip-map region
474  Float integrate(int depth, int zBlock, int phiBlock, const SHVector &f) const;
475 protected:
476  int m_bands;
477  int m_depth;
482 };
483
485
486 #endif /* __MITSUBA_CORE_SHVECTOR_H_ */
Implementation of &#39;Importance Sampling Spherical Harmonics&#39; by W. Jarsz, N. Carr and H...
Definition: shvector.h:430
SHVector & operator*=(Float f)
Scalar multiplication.
Definition: shvector.h:197
int m_depth
Definition: shvector.h:477
SHRotation(int bands)
Construct a new rotation storage for the given number of bands.
Definition: shvector.h:47
int getBands() const
Return the number of stored SH coefficient bands.
Definition: shvector.h:106
SHVector & madd(Float f, const SHVector &v)
Add a scalar multiple of another vector.
Definition: shvector.h:184
Float *** m_legendreMap
Definition: shvector.h:479
int P(int m) const
Definition: shvector.h:463
Float lookupIntegral(int depth, int zBlock, int phiBlock, int l, int m) const
Definition: shvector.h:465
Definition: fwd.h:48
#define SQRT_TWO
Definition: constants.h:94
Stores the diagonal blocks of a spherical harmonic rotation matrix.
Definition: shvector.h:41
Float *** m_phiMap
Definition: shvector.h:478
SHVector(const SHVector &v)
Copy constructor.
Definition: shvector.h:101
int m_bands
Definition: shvector.h:476
Eigen::Matrix< Float, Eigen::Dynamic, Eigen::Dynamic > Matrix
Definition: shvector.h:42
Float * m_normalization
Definition: shvector.h:481
#define MTS_EXPORT_CORE
Definition: getopt.h:29
int m_dataSize
Definition: shvector.h:480
SHVector & operator+=(const SHVector &v)
Definition: shvector.h:134
Normal normalize(const Normal &n)
Definition: normal.h:73
#define M_PI
Definition: constants.h:90
float legendreP(int l, float x)
Evaluate the l-th Legendre polynomial using recurrence (single precision)
#define MTS_NAMESPACE_BEGIN
Definition: platform.h:137
SHVector operator-() const
Negation operator.
Definition: shvector.h:223
SHVector()
Construct an invalid SH vector.
Definition: shvector.h:87
bool operator==(const SHVector &v) const
Equality comparison operator.
Definition: shvector.h:265
SHVector & operator-=(const SHVector &v)
Component-wise subtraction.
Definition: shvector.h:159
SHVector & operator/=(Float f)
Scalar division.
Definition: shvector.h:210
SHVector(int bands)
Construct a new SH vector (initialized to zero)
Definition: shvector.h:92
Float energy(int band) const
Get the energy per band.
Definition: shvector.h:114
SHVector & operator=(const SHVector &v)
Assignment.
Definition: shvector.h:122
Float & operator()(int l, int m)
Access coefficient m (in {-l, ..., l}) on band l.
Definition: shvector.h:230
#define SAssert(cond)
Static&#39;&#39; assertion (to be used outside of classes that derive from Object)
Definition: logger.h:79
Abstract seekable stream class.
Definition: stream.h:58
#define MTS_DECLARE_CLASS()
This macro must be used in the initial definition in classes that derive from Object.
Definition: class.h:158
bool operator!=(const SHVector &v) const
Equality comparison operator.
Definition: shvector.h:270
Definition: fwd.h:99
void clear()
Set all coefficients to zero.
Definition: shvector.h:129
Encapsulates a 4x4 linear transformation and its inverse.
Definition: transform.h:33
Float l2Error(const Functor &f, int res=32) const
Compute the relative L2 error.
Definition: shvector.h:349
MTS_EXPORT_CORE Vector sphericalDirection(Float theta, Float phi)
Convert spherical coordinates to a direction.
SHVector operator+(const SHVector &v) const
Definition: shvector.h:146
SHVector operator*(Float f) const
Scalar multiplication.
Definition: shvector.h:203
Definition: fwd.h:96
SHVector operator-(const SHVector &v) const
Component-wise subtraction.
Definition: shvector.h:171
static Float normalization(int l, int m)
Return a normalization coefficient.
Definition: shvector.h:385
Float dot(const SHVector &v1, const SHVector &v2)
Definition: shvector.h:418
const Float & operator()(int l, int m) const
Access coefficient m (in {-l, ..., l}) on band l.
Definition: shvector.h:235
Parent of all Mitsuba classes.
Definition: object.h:38
virtual std::string toString() const
Return a human-readable string representation of the object&#39;s contents.
void project(const Functor &f, int res=32)
Project the given function onto a SH basis (using a 2D composite Simpson&#39;s rule)
Definition: shvector.h:298
#define MTS_NAMESPACE_END
Definition: platform.h:138
SHVector operator/(Float f) const
Scalar division.
Definition: shvector.h:216
Basic 3x3 matrix data type.
Definition: matrix.h:560
Stores a truncated real spherical harmonics representation of an L2-integrable function.
Definition: shvector.h:84
std::vector< Matrix > blocks
Definition: shvector.h:44
#define SH_NORMTBL_SIZE
Definition: shvector.h:30