diff --git a/include/mitsuba/core/frame.h b/include/mitsuba/core/frame.h index 3410ebed..ac01156f 100644 --- a/include/mitsuba/core/frame.h +++ b/include/mitsuba/core/frame.h @@ -135,7 +135,7 @@ struct Frame { Float sinTheta = Frame::sinTheta(v); if (sinTheta == 0.0f) return 1.0f; - return clamp(v.y / sinTheta, (Float) -1.0f, (Float) 1.0f); + return math::clamp(v.y / sinTheta, (Float) -1.0f, (Float) 1.0f); } /** \brief Assuming that the given direction is in the local coordinate @@ -144,21 +144,21 @@ struct Frame { Float sinTheta = Frame::sinTheta(v); if (sinTheta == 0.0f) return 1.0f; - return clamp(v.x / sinTheta, (Float) -1.0f, (Float) 1.0f); + return math::clamp(v.x / sinTheta, (Float) -1.0f, (Float) 1.0f); } /** \brief Assuming that the given direction is in the local coordinate * system, return the squared sine of the phi parameter in spherical * coordinates */ inline static Float sinPhi2(const Vector &v) { - return clamp(v.y * v.y / sinTheta2(v), (Float) 0.0f, (Float) 1.0f); + return math::clamp(v.y * v.y / sinTheta2(v), (Float) 0.0f, (Float) 1.0f); } /** \brief Assuming that the given direction is in the local coordinate * system, return the squared cosine of the phi parameter in spherical * coordinates */ inline static Float cosPhi2(const Vector &v) { - return clamp(v.x * v.x / sinTheta2(v), (Float) 0.0f, (Float) 1.0f); + return math::clamp(v.x * v.x / sinTheta2(v), (Float) 0.0f, (Float) 1.0f); } /// Equality test diff --git a/include/mitsuba/core/math.h b/include/mitsuba/core/math.h index 95da452e..0cec8dec 100644 --- a/include/mitsuba/core/math.h +++ b/include/mitsuba/core/math.h @@ -23,7 +23,143 @@ MTS_NAMESPACE_BEGIN +/** + * Contains elementary 1D math functions that were either not provided by the standard, + * or which are not consistently provided on all platforms/compilers + */ namespace math { + +/// Cross-platform implementation of the error function +extern MTS_EXPORT_CORE Float erf(Float x); + +/// Cross-platform implementation of the inverse error function +extern MTS_EXPORT_CORE Float erfinv(Float x); + +/// sqrt(a^2 + b^2) without range issues (like 'hypot' on compilers that support C99, single precision) +extern MTS_EXPORT_CORE float hypot2(float a, float b); + +/// sqrt(a^2 + b^2) without range issues (like 'hypot' on compilers that support C99, double precision) +extern MTS_EXPORT_CORE double hypot2(double a, double b); + +/// Base-2 logarithm (single precision) +extern MTS_EXPORT_CORE float log2(float value); + +/// Base-2 logarithm (double precision) +extern MTS_EXPORT_CORE double log2(double value); + +/// Generic clamping function +template inline Scalar clamp(Scalar value, Scalar min, Scalar max) { + return std::min(max, std::max(min, value)); +} + +/// Linearly interpolate between two values +template inline Scalar lerp(Scalar t, Scalar v1, Scalar v2) { + return ((Scalar) 1 - t) * v1 + t * v2; +} + +/// S-shaped smoothly varying interpolation between two values +template inline Scalar smoothStep(Scalar min, Scalar max, Scalar value) { + Scalar v = clamp((value - min) / (max - min), (Scalar) 0, (Scalar) 1); + return v * v * (-2 * v + 3); +} + +/// Always-positive modulo function (assumes b > 0) +inline int32_t modulo(int32_t a, int32_t b) { + int32_t r = a % b; + return (r < 0) ? r+b : r; +} + +/// Always-positive modulo function (assumes b > 0) +inline int64_t modulo(int64_t a, int64_t b) { + int64_t r = a % b; + return (r < 0) ? r+b : r; +} + +#if defined(MTS_AMBIGUOUS_SIZE_T) +inline ssize_t modulo(ssize_t a, ssize_t b) { + if (sizeof(ssize_t) == 8) + return modulo((int64_t) a, (int64_t) b); + else + return modulo((int32_t) a, (int32_t) b); +} +#endif + +/// Always-positive modulo function, single precision version (assumes b > 0) +inline float modulo(float a, float b) { + float r = std::fmod(a, b); + return (r < 0.0f) ? r+b : r; +} + +/// Always-positive modulo function, double precision version (assumes b > 0) +inline double modulo(double a, double b) { + double r = std::fmod(a, b); + return (r < 0.0) ? r+b : r; +} + +/// Integer floor function (single precision) +template inline int floorToInt(Scalar value) { return (int) std::floor(value); } + +/// Integer ceil function (single precision) +template inline int ceilToInt(Scalar value) { return (int) std::ceil(value); } + +/// Integer round function (single precision) +inline int roundToInt(float value) { return (int) ::roundf(value); } + +/// Integer round function (double precision) +inline int roundToInt(double value) { return (int) ::round(value); } + +/// Base-2 logarithm (32-bit integer version) +extern MTS_EXPORT_CORE int log2i(uint32_t value); + +/// Base-2 logarithm (64-bit integer version) +extern MTS_EXPORT_CORE int log2i(uint64_t value); + +#if defined(MTS_AMBIGUOUS_SIZE_T) +inline int log2i(size_t value) { + if (sizeof(size_t) == 8) + return log2i((uint64_t) value); + else + return log2i((uint32_t) value); +} +#endif + +/// Check if an integer is a power of two (unsigned 32 bit version) +inline bool isPowerOfTwo(uint32_t i) { return (i & (i-1)) == 0; } + +/// Check if an integer is a power of two (signed 32 bit version) +inline bool isPowerOfTwo(int32_t i) { return i > 0 && (i & (i-1)) == 0; } + +/// Check if an integer is a power of two (64 bit version) +inline bool isPowerOfTwo(uint64_t i) { return (i & (i-1)) == 0; } + +/// Check if an integer is a power of two (signed 64 bit version) +inline bool isPowerOfTwo(int64_t i) { return i > 0 && (i & (i-1)) == 0; } + +#if defined(MTS_AMBIGUOUS_SIZE_T) +inline bool isPowerOfTwo(size_t value) { + if (sizeof(size_t) == 8) /// will be optimized away + return isPowerOfTwo((uint64_t) value); + else + return isPowerOfTwo((uint32_t) value); +} +#endif + +/// Round an integer to the next power of two +extern MTS_EXPORT_CORE uint32_t roundToPowerOfTwo(uint32_t i); + +/// Round an integer to the next power of two (64 bit version) +extern MTS_EXPORT_CORE uint64_t roundToPowerOfTwo(uint64_t i); + +#if defined(MTS_AMBIGUOUS_SIZE_T) +/// Round an integer to the next power of two +inline size_t roundToPowerOfTwo(size_t value) { + if (sizeof(size_t) == 8) /// will be optimized away + return (size_t) roundToPowerOfTwo((uint64_t) value); + else + return (size_t) roundToPowerOfTwo((uint32_t) value); +} +#endif + #if defined(__LINUX__) && defined(__x86_64__) /* The Linux/x86_64 single precision implementations of 'exp' diff --git a/include/mitsuba/core/matrix.inl b/include/mitsuba/core/matrix.inl index 4c71d427..388fcc4e 100644 --- a/include/mitsuba/core/matrix.inl +++ b/include/mitsuba/core/matrix.inl @@ -341,7 +341,7 @@ template void // Compute implicit shift T g = d[l]; T p = (d[l + 1] - g) / (2.0f * e[l]); - T r = hypot2(1, p); + T r = math::hypot2((T) 1, p); if (p < 0) r = -r; @@ -368,7 +368,7 @@ template void s2 = s; g = c * e[i]; h = c * p; - r = hypot2(p, e[i]); + r = math::hypot2(p, e[i]); e[i + 1] = s * r; s = e[i] / r; c = p / r; diff --git a/include/mitsuba/core/quat.h b/include/mitsuba/core/quat.h index bd8577fd..6ae59141 100644 --- a/include/mitsuba/core/quat.h +++ b/include/mitsuba/core/quat.h @@ -367,7 +367,7 @@ template inline TQuaternion slerp(const TQuaternion &q1, // Revert to plain linear interpolation return normalize(q1 * (1.0f - t) + q2 * t); } else { - Float theta = math::safe_acos(clamp(cosTheta, (Float) -1.0f, (Float) 1.0f)); + Float theta = math::safe_acos(math::clamp(cosTheta, (Float) -1.0f, (Float) 1.0f)); Float thetap = theta * t; TQuaternion qperp = normalize(q2 - q1 * cosTheta); return q1 * std::cos(thetap) + qperp * std::sin(thetap); diff --git a/include/mitsuba/core/rfilter.h b/include/mitsuba/core/rfilter.h index 1c3dab41..c3632dfb 100644 --- a/include/mitsuba/core/rfilter.h +++ b/include/mitsuba/core/rfilter.h @@ -133,7 +133,7 @@ template struct Resampler { filterRadius *= scale; } - m_taps = ceilToInt(filterRadius * 2); + m_taps = math::ceilToInt(filterRadius * 2); if (sourceRes == targetRes && (m_taps % 2) != 1) --m_taps; m_halfTaps = m_taps / 2; @@ -149,7 +149,7 @@ template struct Resampler { Float center = (i + (Float) 0.5f) / targetRes * sourceRes; /* Determine the index of the first original sample that might contribute */ - m_start[i] = floorToInt(center - filterRadius + (Float) 0.5f); + m_start[i] = math::floorToInt(center - filterRadius + (Float) 0.5f); /* Determine the size of center region, on which to run fast non condition-aware code */ if (m_start[i] < 0) @@ -438,13 +438,13 @@ private: if (EXPECT_NOT_TAKEN(pos < 0 || pos >= m_sourceRes)) { switch (m_bc) { case ReconstructionFilter::EClamp: - pos = clamp(pos, 0, m_sourceRes - 1); + pos = math::clamp(pos, 0, m_sourceRes - 1); break; case ReconstructionFilter::ERepeat: - pos = modulo(pos, m_sourceRes); + pos = math::modulo(pos, m_sourceRes); break; case ReconstructionFilter::EMirror: - pos = modulo(pos, 2*m_sourceRes); + pos = math::modulo(pos, 2*m_sourceRes); if (pos >= m_sourceRes) pos = 2*m_sourceRes - pos - 1; break; diff --git a/include/mitsuba/core/ssemath.h b/include/mitsuba/core/ssemath.h index 32b58d99..3db56ddf 100644 --- a/include/mitsuba/core/ssemath.h +++ b/include/mitsuba/core/ssemath.h @@ -81,7 +81,7 @@ namespace math { * holds column 1 of the original matrix, and so on. * \author Intel Intrinsics Guide for AVX2 */ - FINLINE void transpose(__m128& row0, __m128& row1, + FINLINE void transpose_ps(__m128& row0, __m128& row1, __m128& row2, __m128& row3) { __m128 tmp3, tmp2, tmp1, tmp0; tmp0 = _mm_unpacklo_ps(row0, row1); @@ -96,12 +96,12 @@ namespace math { } /// Component-wise clamp: max(min(x, maxVal), minVal) - inline __m128 clamp(__m128 x, __m128 minVal, __m128 maxVal) { + inline __m128 clamp_ps(__m128 x, __m128 minVal, __m128 maxVal) { return _mm_max_ps(_mm_min_ps(x, maxVal), minVal); } /// Sum of all elements in the vector - inline float hsum(__m128 vec) { + inline float hsum_ps(__m128 vec) { __m128 tmp = _mm_shuffle_ps(vec, vec, _MM_SHUFFLE(1,0,3,2)); __m128 sum_tmp = _mm_add_ps(vec, tmp); tmp = _mm_shuffle_ps(sum_tmp, sum_tmp, _MM_SHUFFLE(2,3,0,1)); @@ -110,7 +110,7 @@ namespace math { } /// Maximum across all the elements of a vector - inline float hmax(__m128 vec) { + inline float hmax_ps(__m128 vec) { __m128 tmp = _mm_shuffle_ps(vec, vec, _MM_SHUFFLE(1,0,3,2)); __m128 tmp_max = _mm_max_ps(vec, tmp); tmp = _mm_shuffle_ps(tmp_max, tmp_max, _MM_SHUFFLE(2,3,0,1)); @@ -119,7 +119,7 @@ namespace math { } /// Minimum across all the elements of a vector - inline float hmin(__m128 vec) { + inline float hmin_ps(__m128 vec) { __m128 tmp = _mm_shuffle_ps(vec, vec, _MM_SHUFFLE(1,0,3,2)); __m128 tmp_min = _mm_min_ps(vec, tmp); tmp = _mm_shuffle_ps(tmp_min, tmp_min, _MM_SHUFFLE(2,3,0,1)); diff --git a/include/mitsuba/core/ssevector.h b/include/mitsuba/core/ssevector.h index d895b810..168ab881 100644 --- a/include/mitsuba/core/ssevector.h +++ b/include/mitsuba/core/ssevector.h @@ -477,7 +477,7 @@ FINLINE void transpose(SSEVector4f& row0, SSEVector4f& row1, row3 = _mm_movehl_ps(tmp3, tmp1); } -} // namespace math +} // namespace sse MTS_NAMESPACE_END diff --git a/include/mitsuba/core/util.h b/include/mitsuba/core/util.h index 621d5673..f3e8ae14 100644 --- a/include/mitsuba/core/util.h +++ b/include/mitsuba/core/util.h @@ -271,121 +271,6 @@ template void permute_inplace( //! @{ \name Numerical utility functions // ----------------------------------------------------------------------- -/// sqrt(a^2 + b^2) without underflow (like 'hypot' on compilers that support C99) -extern MTS_EXPORT_CORE Float hypot2(Float a, Float b); - -/// Base-2 logarithm -extern MTS_EXPORT_CORE Float log2(Float value); - -/// Always-positive modulo function (assumes b > 0) -inline int32_t modulo(int32_t a, int32_t b) { - int32_t r = a % b; - return (r < 0) ? r+b : r; -} - -/// Always-positive modulo function (assumes b > 0) -inline int64_t modulo(int64_t a, int64_t b) { - int64_t r = a % b; - return (r < 0) ? r+b : r; -} - -#if defined(MTS_AMBIGUOUS_SIZE_T) -inline ssize_t modulo(ssize_t a, ssize_t b) { - if (sizeof(ssize_t) == 8) - return modulo((int64_t) a, (int64_t) b); - else - return modulo((int32_t) a, (int32_t) b); -} -#endif - -/// Always-positive modulo function, float version (assumes b > 0) -inline Float modulo(Float a, Float b) { - Float r = std::fmod(a, b); - return (r < 0) ? r+b : r; -} - -/// Compute the signum (a.k.a. "sign") function -inline Float signum(Float value) { - if (value < 0) - return -1; - else if (value > 0) - return 1; - else - return 0; -} - -/// Compute the signum (a.k.a. "sign") function, and return an integer value -inline int signumToInt(Float value) { - if (value < 0) - return -1; - else if (value > 0) - return 1; - else - return 0; -} - -/// Integer floor function -inline int floorToInt(Float value) { - return (int) std::floor(value); -} - -/// Integer ceil function -inline int ceilToInt(Float value) { - return (int) std::ceil(value); -} - -/// Base-2 logarithm (32-bit integer version) -extern MTS_EXPORT_CORE int log2i(uint32_t value); - -/// Base-2 logarithm (64-bit integer version) -extern MTS_EXPORT_CORE int log2i(uint64_t value); - -#if defined(MTS_AMBIGUOUS_SIZE_T) -inline int log2i(size_t value) { - if (sizeof(size_t) == 8) - return log2i((uint64_t) value); - else - return log2i((uint32_t) value); -} -#endif - -/// Check if an integer is a power of two (unsigned 32 bit version) -inline bool isPowerOfTwo(uint32_t i) { return (i & (i-1)) == 0; } - -/// Check if an integer is a power of two (signed 32 bit version) -inline bool isPowerOfTwo(int32_t i) { return i > 0 && (i & (i-1)) == 0; } - -/// Check if an integer is a power of two (64 bit version) -inline bool isPowerOfTwo(uint64_t i) { return (i & (i-1)) == 0; } - -/// Check if an integer is a power of two (signed 64 bit version) -inline bool isPowerOfTwo(int64_t i) { return i > 0 && (i & (i-1)) == 0; } - -#if defined(MTS_AMBIGUOUS_SIZE_T) -inline bool isPowerOfTwo(size_t value) { - if (sizeof(size_t) == 8) /// will be optimized away - return isPowerOfTwo((uint64_t) value); - else - return isPowerOfTwo((uint32_t) value); -} -#endif - -/// Round an integer to the next power of two -extern MTS_EXPORT_CORE uint32_t roundToPowerOfTwo(uint32_t i); - -/// Round an integer to the next power of two (64 bit version) -extern MTS_EXPORT_CORE uint64_t roundToPowerOfTwo(uint64_t i); - -#if defined(MTS_AMBIGUOUS_SIZE_T) -/// Round an integer to the next power of two -inline size_t roundToPowerOfTwo(size_t value) { - if (sizeof(size_t) == 8) /// will be optimized away - return (size_t) roundToPowerOfTwo((uint64_t) value); - else - return (size_t) roundToPowerOfTwo((uint32_t) value); -} -#endif - /** * \brief Solve a quadratic equation of the form a*x^2 + b*x + c = 0. * \return \c true if a solution could be found @@ -407,22 +292,6 @@ inline Float radToDeg(Float value) { return value * (180.0f / M_PI); } /// Convert degrees to radians inline Float degToRad(Float value) { return value * (M_PI / 180.0f); } -/// Generic clamping function -template inline Scalar clamp(Scalar value, Scalar min, Scalar max) { - return std::min(max, std::max(min, value)); -} - -/// Linearly interpolate between two values -inline Float lerp(Float t, Float v1, Float v2) { - return ((Float) 1 - t) * v1 + t * v2; -} - -/// S-shaped smoothly varying interpolation between two values -inline Float smoothStep(Float min, Float max, Float value) { - Float v = clamp((value - min) / (max - min), (Float) 0, (Float) 1); - return v * v * (-2 * v + 3); -} - /** * \brief Numerically well-behaved routine for computing the angle * between two unit direction vectors diff --git a/include/mitsuba/core/warp.h b/include/mitsuba/core/warp.h index ad6b2275..da0e4a12 100644 --- a/include/mitsuba/core/warp.h +++ b/include/mitsuba/core/warp.h @@ -31,29 +31,28 @@ MTS_NAMESPACE_BEGIN * The main application of this class is to generate uniformly * distributed or weighted point sets in certain common target domains. */ -class MTS_EXPORT_CORE Warp { -public: +namespace warp { // ============================================================= //! @{ \name Warping techniques related to spheres and subsets // ============================================================= /// Uniformly sample a vector on the unit sphere with respect to solid angles - static Vector squareToUniformSphere(const Point2 &sample); + extern MTS_EXPORT_CORE Vector squareToUniformSphere(const Point2 &sample); /// Density of \ref squareToUniformSphere() with respect to solid angles - static inline Float squareToUniformSpherePdf() { return INV_FOURPI; } + extern MTS_EXPORT_CORE inline Float squareToUniformSpherePdf() { return INV_FOURPI; } /// Uniformly sample a vector on the unit hemisphere with respect to solid angles - static Vector squareToUniformHemisphere(const Point2 &sample); + extern MTS_EXPORT_CORE Vector squareToUniformHemisphere(const Point2 &sample); /// Density of \ref squareToUniformHemisphere() with respect to solid angles - static inline Float squareToUniformHemispherePdf() { return INV_TWOPI; } + extern MTS_EXPORT_CORE inline Float squareToUniformHemispherePdf() { return INV_TWOPI; } /// Sample a cosine-weighted vector on the unit hemisphere with respect to solid angles - static Vector squareToCosineHemisphere(const Point2 &sample); + extern MTS_EXPORT_CORE Vector squareToCosineHemisphere(const Point2 &sample); /// Density of \ref squareToCosineHemisphere() with respect to solid angles - static inline Float squareToCosineHemispherePdf(const Vector &d) + extern MTS_EXPORT_CORE inline Float squareToCosineHemispherePdf(const Vector &d) { return INV_PI * Frame::cosTheta(d); } /** @@ -63,7 +62,7 @@ public: * \param cosCutoff Cosine of the cutoff angle * \param sample A uniformly distributed sample on \f$[0,1]^2\f$ */ - static Vector squareToUniformCone(Float cosCutoff, const Point2 &sample); + extern MTS_EXPORT_CORE Vector squareToUniformCone(Float cosCutoff, const Point2 &sample); /** * \brief Uniformly sample a vector that lies within a given @@ -72,7 +71,7 @@ public: * \param cosCutoff Cosine of the cutoff angle * \param sample A uniformly distributed sample on \f$[0,1]^2\f$ */ - static inline Float squareToUniformConePdf(Float cosCutoff) { + extern MTS_EXPORT_CORE inline Float squareToUniformConePdf(Float cosCutoff) { return INV_TWOPI / (1-cosCutoff); } @@ -84,41 +83,41 @@ public: // ============================================================= /// Uniformly sample a vector on a 2D disk - static Point2 squareToUniformDisk(const Point2 &sample); + extern MTS_EXPORT_CORE Point2 squareToUniformDisk(const Point2 &sample); /// Density of \ref squareToUniformDisk per unit area - static inline Float squareToUniformDiskPdf() { return INV_PI; } + extern MTS_EXPORT_CORE inline Float squareToUniformDiskPdf() { return INV_PI; } /// Low-distortion concentric square to disk mapping by Peter Shirley (PDF: 1/PI) - static Point2 squareToUniformDiskConcentric(const Point2 &sample); + extern MTS_EXPORT_CORE Point2 squareToUniformDiskConcentric(const Point2 &sample); /// Inverse of the mapping \ref squareToUniformDiskConcentric - static Point2 uniformDiskToSquareConcentric(const Point2 &p); + extern MTS_EXPORT_CORE Point2 uniformDiskToSquareConcentric(const Point2 &p); /// Density of \ref squareToUniformDisk per unit area - static inline Float squareToUniformDiskConcentricPdf() { return INV_PI; } + extern MTS_EXPORT_CORE inline Float squareToUniformDiskConcentricPdf() { return INV_PI; } /// Convert an uniformly distributed square sample into barycentric coordinates - static Point2 squareToUniformTriangle(const Point2 &sample); + extern MTS_EXPORT_CORE Point2 squareToUniformTriangle(const Point2 &sample); /** * \brief Sample a point on a 2D standard normal distribution * * Internally uses the Box-Muller transformation */ - static Point2 squareToStdNormal(const Point2 &sample); + extern MTS_EXPORT_CORE Point2 squareToStdNormal(const Point2 &sample); /// Density of \ref squareToStdNormal per unit area - static Float squareToStdNormalPdf(const Point2 &pos); + extern MTS_EXPORT_CORE Float squareToStdNormalPdf(const Point2 &pos); /// Warp a uniformly distributed square sample to a 2D tent distribution - static Point2 squareToTent(const Point2 &sample); + extern MTS_EXPORT_CORE Point2 squareToTent(const Point2 &sample); /** * \brief Warp a uniformly distributed sample on [0, 1] to a nonuniform * tent distribution with nodes {a, b, c} */ - static Float intervalToNonuniformTent(Float a, Float b, Float c, Float sample); + extern MTS_EXPORT_CORE Float intervalToNonuniformTent(Float a, Float b, Float c, Float sample); //! @} // ============================================================= diff --git a/include/mitsuba/render/gkdtree.h b/include/mitsuba/render/gkdtree.h index 52744844..d482533f 100644 --- a/include/mitsuba/render/gkdtree.h +++ b/include/mitsuba/render/gkdtree.h @@ -978,7 +978,7 @@ protected: /* Establish an ad-hoc depth cutoff value (Formula from PBRT) */ if (m_maxDepth == 0) - m_maxDepth = (int) (8 + 1.3f * log2i(primCount)); + m_maxDepth = (int) (8 + 1.3f * math::log2i(primCount)); m_maxDepth = std::min(m_maxDepth, (SizeType) MTS_KD_MAXDEPTH); KDLog(m_logLevel, "Creating a preliminary index list (%s)", diff --git a/include/mitsuba/render/mipmap.h b/include/mitsuba/render/mipmap.h index d5e9a959..5a878e35 100644 --- a/include/mitsuba/render/mipmap.h +++ b/include/mitsuba/render/mipmap.h @@ -508,15 +508,15 @@ public: switch (m_bcu) { case ReconstructionFilter::ERepeat: // Assume that the input repeats in a periodic fashion - x = modulo(x, size.x); + x = math::modulo(x, size.x); break; case ReconstructionFilter::EClamp: // Clamp to the outermost sample position - x = clamp(x, 0, size.x - 1); + x = math::clamp(x, 0, size.x - 1); break; case ReconstructionFilter::EMirror: // Assume that the input is mirrored along the boundary - x = modulo(x, 2*size.x); + x = math::modulo(x, 2*size.x); if (x >= size.x) x = 2*size.x - x - 1; break; @@ -536,15 +536,15 @@ public: switch (m_bcv) { case ReconstructionFilter::ERepeat: // Assume that the input repeats in a periodic fashion - y = modulo(y, size.y); + y = math::modulo(y, size.y); break; case ReconstructionFilter::EClamp: // Clamp to the outermost sample position - y = clamp(y, 0, size.y - 1); + y = math::clamp(y, 0, size.y - 1); break; case ReconstructionFilter::EMirror: // Assume that the input is mirrored along the boundary - y = modulo(y, 2*size.y); + y = math::modulo(y, 2*size.y); if (y >= size.y) y = 2*size.y - y - 1; break; @@ -565,7 +565,7 @@ public: /// Evaluate the texture at the given resolution using a box filter inline Value evalBox(int level, const Point2 &uv) const { const Vector2i &size = m_pyramid[level].getSize(); - return evalTexel(level, floorToInt(uv.x*size.x), floorToInt(uv.y*size.y)); + return evalTexel(level, math::floorToInt(uv.x*size.x), math::floorToInt(uv.y*size.y)); } /** @@ -585,7 +585,7 @@ public: const Vector2i &size = m_pyramid[level].getSize(); Float u = uv.x * size.x - 0.5f, v = uv.y * size.y - 0.5f; - int xPos = floorToInt(u), yPos = floorToInt(v); + int xPos = math::floorToInt(u), yPos = math::floorToInt(v); Float dx1 = u - xPos, dx2 = 1.0f - dx1, dy1 = v - yPos, dy2 = 1.0f - dy1; @@ -612,7 +612,7 @@ public: const Vector2i &size = m_pyramid[level].getSize(); Float u = uv.x * size.x - 0.5f, v = uv.y * size.y - 0.5f; - int xPos = floorToInt(u), yPos = floorToInt(v); + int xPos = math::floorToInt(u), yPos = math::floorToInt(v); Float dx = u - xPos, dy = v - yPos; const Value p00 = evalTexel(level, xPos, yPos); @@ -645,7 +645,7 @@ public: F = A*C - B*B*0.25f; /* Compute the major and minor radii */ - Float root = hypot2(A-C, B), + Float root = math::hypot2(A-C, B), Aprime = 0.5f * (A + C - root), Cprime = 0.5f * (A + C + root), majorRadius = Aprime != 0 ? std::sqrt(F / Aprime) : 0, @@ -654,8 +654,8 @@ public: if (m_filterType == ETrilinear || !(minorRadius > 0) || !(majorRadius > 0) || F < 0) { /* Determine a suitable mip map level, while preferring blurring over aliasing */ - Float level = log2(std::max(majorRadius, Epsilon)); - int ilevel = floorToInt(level); + Float level = math::log2(std::max(majorRadius, Epsilon)); + int ilevel = math::floorToInt(level); if (ilevel < 0) { /* Bilinear interpolation (lookup is smaller than 1 pixel) */ @@ -701,7 +701,7 @@ public: /* Determine a suitable MIP map level, such that the filter covers a reasonable amount of pixels */ - Float level = std::max((Float) 0.0f, log2(minorRadius)); + Float level = std::max((Float) 0.0f, math::log2(minorRadius)); int ilevel = (int) level; Float a = level - ilevel; @@ -786,8 +786,8 @@ protected: Float invDet = 1.0f / (-B*B + 4.0f*A*C), deltaU = 2.0f * std::sqrt(C * invDet), deltaV = 2.0f * std::sqrt(A * invDet); - int u0 = ceilToInt(u - deltaU), u1 = floorToInt(u + deltaU); - int v0 = ceilToInt(v - deltaV), v1 = floorToInt(v + deltaV); + int u0 = math::ceilToInt(u - deltaU), u1 = math::floorToInt(u + deltaU); + int v0 = math::ceilToInt(v - deltaV), v1 = math::floorToInt(v + deltaV); /* Scale the coefficients by the size of the Gaussian lookup table */ Float As = A * MTS_MIPMAP_LUT_SIZE, diff --git a/include/mitsuba/render/sahkdtree3.h b/include/mitsuba/render/sahkdtree3.h index 8053d121..a0070ab5 100644 --- a/include/mitsuba/render/sahkdtree3.h +++ b/include/mitsuba/render/sahkdtree3.h @@ -518,8 +518,8 @@ public: for (int i=0; inextFloat(), random->nextFloat()), sample2(random->nextFloat(), random->nextFloat()); - Point p1 = bsphere.center + Warp::squareToUniformSphere(sample1) * bsphere.radius; - Point p2 = bsphere.center + Warp::squareToUniformSphere(sample2) * bsphere.radius; + Point p1 = bsphere.center + warp::squareToUniformSphere(sample1) * bsphere.radius; + Point p2 = bsphere.center + warp::squareToUniformSphere(sample2) * bsphere.radius; Ray ray(p1, normalize(p2-p1), 0.0f); Float mint, maxt, t; if (m_aabb.rayIntersect(ray, mint, maxt)) { diff --git a/src/bsdfs/difftrans.cpp b/src/bsdfs/difftrans.cpp index f976d449..de8d37a2 100644 --- a/src/bsdfs/difftrans.cpp +++ b/src/bsdfs/difftrans.cpp @@ -93,7 +93,7 @@ public: Spectrum sample(BSDFSamplingRecord &bRec, const Point2 &sample) const { if (!(bRec.typeMask & EDiffuseTransmission)) return Spectrum(0.0f); - bRec.wo = Warp::squareToCosineHemisphere(sample); + bRec.wo = warp::squareToCosineHemisphere(sample); if (Frame::cosTheta(bRec.wi) > 0) bRec.wo.z *= -1; bRec.eta = 1.0f; @@ -105,7 +105,7 @@ public: Spectrum sample(BSDFSamplingRecord &bRec, Float &pdf, const Point2 &sample) const { if (!(bRec.typeMask & m_combinedType)) return Spectrum(0.0f); - bRec.wo = Warp::squareToCosineHemisphere(sample); + bRec.wo = warp::squareToCosineHemisphere(sample); if (Frame::cosTheta(bRec.wi) > 0) bRec.wo.z *= -1; bRec.eta = 1.0f; diff --git a/src/bsdfs/diffuse.cpp b/src/bsdfs/diffuse.cpp index 7cfe9a65..a12ffa46 100644 --- a/src/bsdfs/diffuse.cpp +++ b/src/bsdfs/diffuse.cpp @@ -123,14 +123,14 @@ public: || Frame::cosTheta(bRec.wo) <= 0) return 0.0f; - return Warp::squareToCosineHemispherePdf(bRec.wo); + return warp::squareToCosineHemispherePdf(bRec.wo); } Spectrum sample(BSDFSamplingRecord &bRec, const Point2 &sample) const { if (!(bRec.typeMask & EDiffuseReflection) || Frame::cosTheta(bRec.wi) <= 0) return Spectrum(0.0f); - bRec.wo = Warp::squareToCosineHemisphere(sample); + bRec.wo = warp::squareToCosineHemisphere(sample); bRec.eta = 1.0f; bRec.sampledComponent = 0; bRec.sampledType = EDiffuseReflection; @@ -141,11 +141,11 @@ public: if (!(bRec.typeMask & EDiffuseReflection) || Frame::cosTheta(bRec.wi) <= 0) return Spectrum(0.0f); - bRec.wo = Warp::squareToCosineHemisphere(sample); + bRec.wo = warp::squareToCosineHemisphere(sample); bRec.eta = 1.0f; bRec.sampledComponent = 0; bRec.sampledType = EDiffuseReflection; - pdf = Warp::squareToCosineHemispherePdf(bRec.wo); + pdf = warp::squareToCosineHemispherePdf(bRec.wo); return m_reflectance->eval(bRec.its); } diff --git a/src/bsdfs/irawan.cpp b/src/bsdfs/irawan.cpp index 53d6cf86..2f03ec71 100644 --- a/src/bsdfs/irawan.cpp +++ b/src/bsdfs/irawan.cpp @@ -155,8 +155,8 @@ public: Spectrum result(0.0f); m_initialization = true; for (size_t i=0; inextFloat(), random->nextFloat())); - bRec.wo = Warp::squareToCosineHemisphere(Point2(random->nextFloat(), random->nextFloat())); + bRec.wi = warp::squareToCosineHemisphere(Point2(random->nextFloat(), random->nextFloat())); + bRec.wo = warp::squareToCosineHemisphere(Point2(random->nextFloat(), random->nextFloat())); its.uv = Point2(random->nextFloat(), random->nextFloat()); result += eval(bRec, ESolidAngle) / Frame::cosTheta(bRec.wo); @@ -177,8 +177,8 @@ public: (1 - its.uv.y) * m_repeatV); Point2 xy(uv.x * m_pattern.tileWidth, uv.y * m_pattern.tileHeight); Point2i lookup( - modulo((int) xy.x, m_pattern.tileWidth), - modulo((int) xy.y, m_pattern.tileHeight)); + math::modulo((int) xy.x, m_pattern.tileWidth), + math::modulo((int) xy.y, m_pattern.tileHeight)); int yarnID = m_pattern.pattern[lookup.x + lookup.y * m_pattern.tileWidth] - 1; const Yarn &yarn = m_pattern.yarns.at(yarnID); @@ -203,8 +203,8 @@ public: Point2 xy(uv.x * m_pattern.tileWidth, uv.y * m_pattern.tileHeight); Point2i lookup( - modulo((int) xy.x, m_pattern.tileWidth), - modulo((int) xy.y, m_pattern.tileHeight)); + math::modulo((int) xy.x, m_pattern.tileWidth), + math::modulo((int) xy.y, m_pattern.tileHeight)); int yarnID = m_pattern.pattern[lookup.x + lookup.y * m_pattern.tileWidth] - 1; @@ -330,7 +330,7 @@ public: measure != ESolidAngle) return 0.0f; - return Warp::squareToCosineHemispherePdf(bRec.wo); + return warp::squareToCosineHemispherePdf(bRec.wo); } Spectrum sample(BSDFSamplingRecord &bRec, const Point2 &sample) const { @@ -344,7 +344,7 @@ public: return Spectrum(0.0f); /* Lacking a better sampling method, generate cosine-weighted directions */ - bRec.wo = Warp::squareToCosineHemisphere(sample); + bRec.wo = warp::squareToCosineHemisphere(sample); bRec.eta = 1.0f; bRec.sampledComponent = 0; bRec.sampledType = EGlossyReflection; @@ -362,11 +362,11 @@ public: return Spectrum(0.0f); /* Lacking a better sampling method, generate cosine-weighted directions */ - bRec.wo = Warp::squareToCosineHemisphere(sample); + bRec.wo = warp::squareToCosineHemisphere(sample); bRec.eta = 1.0f; bRec.sampledComponent = 0; bRec.sampledType = EGlossyReflection; - pdf = Warp::squareToCosineHemispherePdf(bRec.wo); + pdf = warp::squareToCosineHemispherePdf(bRec.wo); return eval(bRec, ESolidAngle) / pdf; } @@ -437,7 +437,7 @@ public: if (ss == 0.0f) As = A; else - As = A * (1.0f - smoothStep(0, 1, (std::abs(u_of_v) + As = A * (1.0f - math::smoothStep((Float) 0, (Float) 1, (std::abs(u_of_v) - (1.0f - ss) * umax) / (ss * umax))); // fs is scattering function. diff --git a/src/bsdfs/phong.cpp b/src/bsdfs/phong.cpp index 911b5702..8131a24e 100644 --- a/src/bsdfs/phong.cpp +++ b/src/bsdfs/phong.cpp @@ -159,7 +159,7 @@ public: Float diffuseProb = 0.0f, specProb = 0.0f; if (hasDiffuse) - diffuseProb = Warp::squareToCosineHemispherePdf(bRec.wo); + diffuseProb = warp::squareToCosineHemispherePdf(bRec.wo); if (hasSpecular) { Float alpha = dot(bRec.wo, reflect(bRec.wi)), @@ -225,7 +225,7 @@ public: if (Frame::cosTheta(bRec.wo) <= 0) return Spectrum(0.0f); } else { - bRec.wo = Warp::squareToCosineHemisphere(sample); + bRec.wo = warp::squareToCosineHemisphere(sample); bRec.sampledComponent = 0; bRec.sampledType = EDiffuseReflection; } diff --git a/src/bsdfs/plastic.cpp b/src/bsdfs/plastic.cpp index aa99bd0f..b048435e 100644 --- a/src/bsdfs/plastic.cpp +++ b/src/bsdfs/plastic.cpp @@ -270,7 +270,7 @@ public: else diff /= 1 - m_fdrInt; - return diff * (Warp::squareToCosineHemispherePdf(bRec.wo) + return diff * (warp::squareToCosineHemispherePdf(bRec.wo) * m_invEta2 * (1-Fi) * (1-Fo)); } @@ -300,7 +300,7 @@ public: if (std::abs(dot(reflect(bRec.wi), bRec.wo)-1) < DeltaEpsilon) return probSpecular; } else if (hasDiffuse && measure == ESolidAngle) { - return Warp::squareToCosineHemispherePdf(bRec.wo) * (1-probSpecular); + return warp::squareToCosineHemispherePdf(bRec.wo) * (1-probSpecular); } return 0.0f; @@ -334,7 +334,7 @@ public: } else { bRec.sampledComponent = 1; bRec.sampledType = EDiffuseReflection; - bRec.wo = Warp::squareToCosineHemisphere(Point2( + bRec.wo = warp::squareToCosineHemisphere(Point2( (sample.x - probSpecular) / (1 - probSpecular), sample.y )); @@ -356,7 +356,7 @@ public: } else { bRec.sampledComponent = 1; bRec.sampledType = EDiffuseReflection; - bRec.wo = Warp::squareToCosineHemisphere(sample); + bRec.wo = warp::squareToCosineHemisphere(sample); Float Fo = fresnelDielectricExt(Frame::cosTheta(bRec.wo), m_eta); Spectrum diff = m_diffuseReflectance->eval(bRec.its); @@ -398,7 +398,7 @@ public: } else { bRec.sampledComponent = 1; bRec.sampledType = EDiffuseReflection; - bRec.wo = Warp::squareToCosineHemisphere(Point2( + bRec.wo = warp::squareToCosineHemisphere(Point2( (sample.x - probSpecular) / (1 - probSpecular), sample.y )); @@ -411,7 +411,7 @@ public: diff /= 1 - m_fdrInt; pdf = (1-probSpecular) * - Warp::squareToCosineHemispherePdf(bRec.wo); + warp::squareToCosineHemispherePdf(bRec.wo); return diff * (m_invEta2 * (1-Fi) * (1-Fo) / (1-probSpecular)); } @@ -424,7 +424,7 @@ public: } else { bRec.sampledComponent = 1; bRec.sampledType = EDiffuseReflection; - bRec.wo = Warp::squareToCosineHemisphere(sample); + bRec.wo = warp::squareToCosineHemisphere(sample); Float Fo = fresnelDielectricExt(Frame::cosTheta(bRec.wo), m_eta); Spectrum diff = m_diffuseReflectance->eval(bRec.its); @@ -433,7 +433,7 @@ public: else diff /= 1 - m_fdrInt; - pdf = Warp::squareToCosineHemispherePdf(bRec.wo); + pdf = warp::squareToCosineHemispherePdf(bRec.wo); return diff * (m_invEta2 * (1-Fi) * (1-Fo)); } diff --git a/src/bsdfs/roughcoating.cpp b/src/bsdfs/roughcoating.cpp index 5f30e7bb..ea315d41 100644 --- a/src/bsdfs/roughcoating.cpp +++ b/src/bsdfs/roughcoating.cpp @@ -262,7 +262,7 @@ public: if (hasSpecular && Frame::cosTheta(bRec.wo) * Frame::cosTheta(bRec.wi) > 0) { /* Calculate the reflection half-vector */ const Vector H = normalize(bRec.wo+bRec.wi) - * signum(Frame::cosTheta(bRec.wo)); + * math::signum(Frame::cosTheta(bRec.wo)); /* Evaluate the microsurface normal distribution */ const Float D = m_distribution.eval(H, alphaT); @@ -317,7 +317,7 @@ public: /* Calculate the reflection half-vector */ const Vector H = normalize(bRec.wo+bRec.wi) - * signum(Frame::cosTheta(bRec.wo)); + * math::signum(Frame::cosTheta(bRec.wo)); /* Evaluate the roughness texture */ Float alpha = m_alpha->eval(bRec.its).average(); diff --git a/src/bsdfs/roughdiffuse.cpp b/src/bsdfs/roughdiffuse.cpp index a107a1f8..91d4f40c 100644 --- a/src/bsdfs/roughdiffuse.cpp +++ b/src/bsdfs/roughdiffuse.cpp @@ -223,30 +223,30 @@ public: || Frame::cosTheta(bRec.wo) <= 0) return 0.0f; - return Warp::squareToCosineHemispherePdf(bRec.wo); + return warp::squareToCosineHemispherePdf(bRec.wo); } Spectrum sample(BSDFSamplingRecord &bRec, const Point2 &sample) const { if (!(bRec.typeMask & EGlossyReflection) || Frame::cosTheta(bRec.wi) <= 0) return Spectrum(0.0f); - bRec.wo = Warp::squareToCosineHemisphere(sample); + bRec.wo = warp::squareToCosineHemisphere(sample); bRec.eta = 1.0f; bRec.sampledComponent = 0; bRec.sampledType = EGlossyReflection; return eval(bRec, ESolidAngle) / - Warp::squareToCosineHemispherePdf(bRec.wo); + warp::squareToCosineHemispherePdf(bRec.wo); } Spectrum sample(BSDFSamplingRecord &bRec, Float &pdf, const Point2 &sample) const { if (!(bRec.typeMask & EGlossyReflection) || Frame::cosTheta(bRec.wi) <= 0) return Spectrum(0.0f); - bRec.wo = Warp::squareToCosineHemisphere(sample); + bRec.wo = warp::squareToCosineHemisphere(sample); bRec.eta = 1.0f; bRec.sampledComponent = 0; bRec.sampledType = EGlossyReflection; - pdf = Warp::squareToCosineHemispherePdf(bRec.wo); + pdf = warp::squareToCosineHemispherePdf(bRec.wo); return eval(bRec, ESolidAngle) / pdf; } diff --git a/src/bsdfs/roughplastic.cpp b/src/bsdfs/roughplastic.cpp index 27d8128c..b08b97d3 100644 --- a/src/bsdfs/roughplastic.cpp +++ b/src/bsdfs/roughplastic.cpp @@ -420,7 +420,7 @@ public: } if (hasDiffuse) - result += probDiffuse * Warp::squareToCosineHemispherePdf(bRec.wo); + result += probDiffuse * warp::squareToCosineHemispherePdf(bRec.wo); return result; } @@ -472,7 +472,7 @@ public: } else { bRec.sampledComponent = 1; bRec.sampledType = EDiffuseReflection; - bRec.wo = Warp::squareToCosineHemisphere(sample); + bRec.wo = warp::squareToCosineHemisphere(sample); } bRec.eta = 1.0f; diff --git a/src/bsdfs/ward.cpp b/src/bsdfs/ward.cpp index 9e2cde9c..e8cd1a16 100644 --- a/src/bsdfs/ward.cpp +++ b/src/bsdfs/ward.cpp @@ -250,7 +250,7 @@ public: } if (hasDiffuse) - diffuseProb = Warp::squareToCosineHemispherePdf(bRec.wo); + diffuseProb = warp::squareToCosineHemispherePdf(bRec.wo); if (hasDiffuse && hasSpecular) return m_specularSamplingWeight * specProb + @@ -311,7 +311,7 @@ public: if (Frame::cosTheta(bRec.wo) <= 0.0f) return Spectrum(0.0f); } else { - bRec.wo = Warp::squareToCosineHemisphere(sample); + bRec.wo = warp::squareToCosineHemisphere(sample); bRec.sampledComponent = 0; bRec.sampledType = EDiffuseReflection; } diff --git a/src/emitters/area.cpp b/src/emitters/area.cpp index 4f5b782a..397b9fd3 100644 --- a/src/emitters/area.cpp +++ b/src/emitters/area.cpp @@ -115,9 +115,9 @@ public: Spectrum sampleDirection(DirectionSamplingRecord &dRec, PositionSamplingRecord &pRec, const Point2 &sample, const Point2 *extra) const { - Vector local = Warp::squareToCosineHemisphere(sample); + Vector local = warp::squareToCosineHemisphere(sample); dRec.d = Frame(pRec.n).toWorld(local); - dRec.pdf = Warp::squareToCosineHemispherePdf(local); + dRec.pdf = warp::squareToCosineHemispherePdf(local); dRec.measure = ESolidAngle; return Spectrum(1.0f); } @@ -148,7 +148,7 @@ public: Float time) const { PositionSamplingRecord pRec(time); m_shape->samplePosition(pRec, spatialSample); - Vector local = Warp::squareToCosineHemisphere(directionalSample); + Vector local = warp::squareToCosineHemisphere(directionalSample); ray.setTime(time); ray.setOrigin(pRec.p); ray.setDirection(Frame(pRec.n).toWorld(local)); diff --git a/src/emitters/constant.cpp b/src/emitters/constant.cpp index 7a7c4c32..eefd32ac 100644 --- a/src/emitters/constant.cpp +++ b/src/emitters/constant.cpp @@ -108,7 +108,7 @@ public: Spectrum samplePosition(PositionSamplingRecord &pRec, const Point2 &sample, const Point2 *extra) const { - Vector d = Warp::squareToUniformSphere(sample); + Vector d = warp::squareToUniformSphere(sample); pRec.p = m_sceneBSphere.center + d * m_sceneBSphere.radius; pRec.n = -d; @@ -129,9 +129,9 @@ public: Spectrum sampleDirection(DirectionSamplingRecord &dRec, PositionSamplingRecord &pRec, const Point2 &sample, const Point2 *extra) const { - Vector local = Warp::squareToCosineHemisphere(sample); + Vector local = warp::squareToCosineHemisphere(sample); dRec.d = Frame(pRec.n).toWorld(local); - dRec.pdf = Warp::squareToCosineHemispherePdf(local); + dRec.pdf = warp::squareToCosineHemispherePdf(local); dRec.measure = ESolidAngle; return Spectrum(1.0f); } @@ -160,8 +160,8 @@ public: const Point2 &spatialSample, const Point2 &directionalSample, Float time) const { - Vector v0 = Warp::squareToUniformSphere(spatialSample); - Vector v1 = Warp::squareToCosineHemisphere(directionalSample); + Vector v0 = warp::squareToUniformSphere(spatialSample); + Vector v1 = warp::squareToCosineHemisphere(directionalSample); ray.setOrigin(m_geoBSphere.center + v0 * m_geoBSphere.radius); ray.setDirection(Frame(-v0).toWorld(v1)); @@ -176,12 +176,12 @@ public: Float pdf; if (!dRec.refN.isZero()) { - d = Warp::squareToCosineHemisphere(sample); - pdf = Warp::squareToCosineHemispherePdf(d); + d = warp::squareToCosineHemisphere(sample); + pdf = warp::squareToCosineHemispherePdf(d); d = Frame(dRec.refN).toWorld(d); } else { - d = Warp::squareToUniformSphere(sample); - pdf = Warp::squareToUniformSpherePdf(); + d = warp::squareToUniformSphere(sample); + pdf = warp::squareToUniformSpherePdf(); } /* Intersect against the bounding sphere. This is not really @@ -220,7 +220,7 @@ public: if (!dRec.refN.isZero()) pdfSA = INV_PI * std::max((Float) 0.0f, dot(dRec.d, dRec.refN)); else - pdfSA = Warp::squareToUniformSpherePdf(); + pdfSA = warp::squareToUniformSpherePdf(); if (dRec.measure == ESolidAngle) return pdfSA; diff --git a/src/emitters/directional.cpp b/src/emitters/directional.cpp index b00db66b..9091c574 100644 --- a/src/emitters/directional.cpp +++ b/src/emitters/directional.cpp @@ -103,7 +103,7 @@ public: Spectrum samplePosition(PositionSamplingRecord &pRec, const Point2 &sample, const Point2 *extra) const { const Transform &trafo = m_worldTransform->eval(pRec.time); - Point2 p = Warp::squareToUniformDiskConcentric(sample); + Point2 p = warp::squareToUniformDiskConcentric(sample); Vector perpOffset = trafo(Vector(p.x, p.y, 0) * m_bsphere.radius); Vector d = trafo(Vector(0, 0, 1)); @@ -147,7 +147,7 @@ public: const Point2 &directionalSample, Float time) const { const Transform &trafo = m_worldTransform->eval(time); - Point2 p = Warp::squareToUniformDiskConcentric(spatialSample); + Point2 p = warp::squareToUniformDiskConcentric(spatialSample); Vector perpOffset = trafo(Vector(p.x, p.y, 0) * m_bsphere.radius); Vector d = trafo(Vector(0, 0, 1)); ray.setOrigin(m_bsphere.center - d*m_bsphere.radius + perpOffset); diff --git a/src/emitters/envmap.cpp b/src/emitters/envmap.cpp index f76cd905..e4d81f53 100644 --- a/src/emitters/envmap.cpp +++ b/src/emitters/envmap.cpp @@ -411,7 +411,7 @@ public: Spectrum samplePosition(PositionSamplingRecord &pRec, const Point2 &sample, const Point2 *extra) const { - Vector d = Warp::squareToUniformSphere(sample); + Vector d = warp::squareToUniformSphere(sample); pRec.p = m_sceneBSphere.center + d * m_sceneBSphere.radius; pRec.n = -d; @@ -503,7 +503,7 @@ public: Vector d; Spectrum value; Float pdf; internalSampleDirection(directionalSample, d, value, pdf); d = -trafo(d); - Point2 offset = Warp::squareToUniformDiskConcentric(spatialSample); + Point2 offset = warp::squareToUniformDiskConcentric(spatialSample); Vector perpOffset = Frame(d).toWorld(Vector(offset.x, offset.y, 0)); ray.setOrigin(m_geoBSphere.center + (perpOffset - d) * m_geoBSphere.radius); @@ -572,10 +572,10 @@ public: /* Using the remaining bits of precision to shift the sample by an offset drawn from a tent function. This effectively creates a sampling strategy for a linearly interpolated environment map */ - Point2 pos = Point2((Float) col, (Float) row) + Warp::squareToTent(sample); + Point2 pos = Point2((Float) col, (Float) row) + warp::squareToTent(sample); /* Bilinearly interpolate colors from the adjacent four neighbors */ - int xPos = floorToInt(pos.x), yPos = floorToInt(pos.y); + int xPos = math::floorToInt(pos.x), yPos = math::floorToInt(pos.y); Float dx1 = pos.x - xPos, dx2 = 1.0f - dx1, dy1 = pos.y - yPos, dy2 = 1.0f - dy1; @@ -587,8 +587,8 @@ public: /* Compute the final color and probability density of the sample */ value = (value1 + value2) * m_scale; - pdf = (value1.getLuminance() * m_rowWeights[clamp(yPos, 0, m_size.y-1)] + - value2.getLuminance() * m_rowWeights[clamp(yPos+1, 0, m_size.y-1)]) * m_normalization; + pdf = (value1.getLuminance() * m_rowWeights[math::clamp(yPos, 0, m_size.y-1)] + + value2.getLuminance() * m_rowWeights[math::clamp(yPos+1, 0, m_size.y-1)]) * m_normalization; /* Turn into a proper direction on the sphere */ Float sinPhi, cosPhi, sinTheta, cosTheta; @@ -616,7 +616,7 @@ public: Float u = uv.x * m_size.x - 0.5f, v = uv.y * m_size.y - 0.5f; /* Bilinearly interpolate colors from the adjacent four neighbors */ - int xPos = floorToInt(u), yPos = floorToInt(v); + int xPos = math::floorToInt(u), yPos = math::floorToInt(v); Float dx1 = u - xPos, dx2 = 1.0f - dx1, dy1 = v - yPos, dy2 = 1.0f - dy1; @@ -627,8 +627,8 @@ public: stats::filteredLookups.incrementBase(); Float sinTheta = math::safe_sqrt(1-d.y*d.y); - return (value1.getLuminance() * m_rowWeights[clamp(yPos, 0, m_size.y-1)] + - value2.getLuminance() * m_rowWeights[clamp(yPos+1, 0, m_size.y-1)]) + return (value1.getLuminance() * m_rowWeights[math::clamp(yPos, 0, m_size.y-1)] + + value2.getLuminance() * m_rowWeights[math::clamp(yPos+1, 0, m_size.y-1)]) * m_normalization / std::max(std::abs(sinTheta), Epsilon); } diff --git a/src/emitters/point.cpp b/src/emitters/point.cpp index e73cc3d7..00ffe1f8 100644 --- a/src/emitters/point.cpp +++ b/src/emitters/point.cpp @@ -101,7 +101,7 @@ public: PositionSamplingRecord &pRec, const Point2 &sample, const Point2 *extra) const { - dRec.d = Warp::squareToUniformSphere(sample); + dRec.d = warp::squareToUniformSphere(sample); dRec.pdf = INV_FOURPI; dRec.measure = ESolidAngle; return Spectrum(1.0f); @@ -124,7 +124,7 @@ public: const Transform &trafo = m_worldTransform->eval(time); ray.setTime(time); ray.setOrigin(trafo(Point(0.0f))); - ray.setDirection(Warp::squareToUniformSphere(directionalSample)); + ray.setDirection(warp::squareToUniformSphere(directionalSample)); return m_intensity * (4 * M_PI); } diff --git a/src/emitters/sky.cpp b/src/emitters/sky.cpp index 28b26583..d91d043e 100644 --- a/src/emitters/sky.cpp +++ b/src/emitters/sky.cpp @@ -442,7 +442,7 @@ protected: result.clampNegative(); if (m_extend) - result *= smoothStep(0, 1, 2 - 2*coords.elevation*INV_PI); + result *= math::smoothStep((Float) 0, (Float) 1, 2 - 2*coords.elevation*INV_PI); return result * m_scale; } diff --git a/src/emitters/spot.cpp b/src/emitters/spot.cpp index 8e1808b4..1356ccb4 100644 --- a/src/emitters/spot.cpp +++ b/src/emitters/spot.cpp @@ -147,16 +147,16 @@ public: const Point2 &sample, const Point2 *extra) const { const Transform &trafo = m_worldTransform->eval(pRec.time); - Vector d = Warp::squareToUniformCone(m_cosCutoffAngle, sample); + Vector d = warp::squareToUniformCone(m_cosCutoffAngle, sample); dRec.d = trafo(d); - dRec.pdf = Warp::squareToUniformConePdf(m_cosCutoffAngle); + dRec.pdf = warp::squareToUniformConePdf(m_cosCutoffAngle); dRec.measure = ESolidAngle; return evalDirection(dRec, pRec)/dRec.pdf; } Float pdfDirection(const DirectionSamplingRecord &dRec, const PositionSamplingRecord &pRec) const { - return (dRec.measure == ESolidAngle) ? Warp::squareToUniformConePdf(m_cosCutoffAngle) : 0.0f; + return (dRec.measure == ESolidAngle) ? warp::squareToUniformConePdf(m_cosCutoffAngle) : 0.0f; } Spectrum evalDirection(const DirectionSamplingRecord &dRec, @@ -172,12 +172,12 @@ public: Float time) const { const Transform &trafo = m_worldTransform->eval(time); - Vector local = Warp::squareToUniformCone( + Vector local = warp::squareToUniformCone( m_cosCutoffAngle, directionalSample); ray.setTime(time); ray.setOrigin(trafo.transformAffine(Point(0.0f))); ray.setDirection(trafo(local)); - Float dirPdf = Warp::squareToUniformConePdf(m_cosCutoffAngle); + Float dirPdf = warp::squareToUniformConePdf(m_cosCutoffAngle); return m_intensity * falloffCurve(local) / dirPdf; } diff --git a/src/emitters/sun.cpp b/src/emitters/sun.cpp index 615dcb1c..21f94d20 100644 --- a/src/emitters/sun.cpp +++ b/src/emitters/sun.cpp @@ -203,7 +203,7 @@ public: for (size_t i=0; irayIntersect(shadowRay)) diff --git a/src/libbidir/mut_manifold.cpp b/src/libbidir/mut_manifold.cpp index f7c82007..50dbb69a 100644 --- a/src/libbidir/mut_manifold.cpp +++ b/src/libbidir/mut_manifold.cpp @@ -440,7 +440,7 @@ bool ManifoldPerturbation::sampleMutation( cosTheta = dot(wo_old, n_old); - Float dTheta = Warp::squareToStdNormal(m_sampler->next2D()).x + Float dTheta = warp::squareToStdNormal(m_sampler->next2D()).x * 0.5f * M_PI / m_probFactor; math::sincos(dTheta, &sinPhi, &cosPhi); diff --git a/src/libbidir/vertex.cpp b/src/libbidir/vertex.cpp index 88a038ee..f62c1481 100644 --- a/src/libbidir/vertex.cpp +++ b/src/libbidir/vertex.cpp @@ -383,7 +383,7 @@ int PathVertex::sampleSensor(const Scene *scene, Sampler *sampler, } bool PathVertex::perturbPosition(const Scene *scene, Sampler *sampler, Float stddev) { - Point2 step = Warp::squareToStdNormal(sampler->next2D()) * stddev; + Point2 step = warp::squareToStdNormal(sampler->next2D()) * stddev; EVertexType type = (EVertexType) this->type; Ray ray; @@ -462,7 +462,7 @@ Float PathVertex::perturbPositionPdf(const PathVertex *target, Float stddev) con Vector rel = itsOld.geoFrame.toLocal(itsOld.p - itsNew.p); Point2 rel2 = Point2(rel.x, rel.y) / stddev; - return Warp::squareToStdNormalPdf(rel2) * absDot(itsOld.geoFrame.n, itsNew.geoFrame.n) / (stddev*stddev); + return warp::squareToStdNormalPdf(rel2) * absDot(itsOld.geoFrame.n, itsNew.geoFrame.n) / (stddev*stddev); } break; @@ -473,7 +473,7 @@ Float PathVertex::perturbPositionPdf(const PathVertex *target, Float stddev) con Vector rel = Frame(prOld.n).toLocal(prOld.p - prNew.p); Point2 rel2 = Point2(rel.x, rel.y) / stddev; - return Warp::squareToStdNormalPdf(rel2) * absDot(prOld.n, prNew.n) / (stddev*stddev); + return warp::squareToStdNormalPdf(rel2) * absDot(prOld.n, prNew.n) / (stddev*stddev); } break; diff --git a/src/libcore/SConscript b/src/libcore/SConscript index 7559d368..b6a741fa 100644 --- a/src/libcore/SConscript +++ b/src/libcore/SConscript @@ -39,7 +39,7 @@ libcore_objects = [ 'class.cpp', 'object.cpp', 'statistics.cpp', 'thread.cpp', 'brent.cpp', 'logger.cpp', 'appender.cpp', 'formatter.cpp', 'lock.cpp', 'qmc.cpp', 'random.cpp', 'timer.cpp', 'util.cpp', 'properties.cpp', 'half.cpp', - 'transform.cpp', 'spectrum.cpp', 'aabb.cpp', 'stream.cpp', + 'transform.cpp', 'spectrum.cpp', 'aabb.cpp', 'stream.cpp', 'math.cpp', 'fstream.cpp', 'plugin.cpp', 'triangle.cpp', 'bitmap.cpp', 'fmtconv.cpp', 'serialization.cpp', 'sstream.cpp', 'cstream.cpp', 'mstream.cpp', 'sched.cpp', 'sched_remote.cpp', 'sshstream.cpp', diff --git a/src/libcore/bitmap.cpp b/src/libcore/bitmap.cpp index 9e5c48cb..8f9f1ef4 100644 --- a/src/libcore/bitmap.cpp +++ b/src/libcore/bitmap.cpp @@ -813,9 +813,9 @@ void Bitmap::convolve(const Bitmap *_kernel) { /* Copy and zero-pad the convolution kernel in a wraparound fashion */ if (ch < channelCountKernel) { for (size_t y=0; ygetFloat16Data()[(x+y*kernelSize)*channelCountKernel+ch]; } } @@ -831,9 +831,9 @@ void Bitmap::convolve(const Bitmap *_kernel) { /* Copy and zero-pad the convolution kernel in a wraparound fashion */ if (ch < channelCountKernel) { for (size_t y=0; ygetFloat32Data()[(x+y*kernelSize)*channelCountKernel+ch]; } } @@ -849,9 +849,9 @@ void Bitmap::convolve(const Bitmap *_kernel) { /* Copy and zero-pad the convolution kernel in a wraparound fashion */ if (ch < channelCountKernel) { for (size_t y=0; ygetFloat64Data()[(x+y*kernelSize)*channelCountKernel+ch]; } } diff --git a/src/libcore/chisquare.cpp b/src/libcore/chisquare.cpp index 98481051..052b15bb 100644 --- a/src/libcore/chisquare.cpp +++ b/src/libcore/chisquare.cpp @@ -105,9 +105,9 @@ void ChiSquare::fill( Point2 sphCoords = toSphericalCoordinates(boost::get<0>(sample)); int thetaBin = std::min(std::max(0, - floorToInt(sphCoords.x * factor.x)), m_thetaBins-1); + math::floorToInt(sphCoords.x * factor.x)), m_thetaBins-1); int phiBin = std::min(std::max(0, - floorToInt(sphCoords.y * factor.y)), m_phiBins-1); + math::floorToInt(sphCoords.y * factor.y)), m_phiBins-1); m_table[thetaBin * m_phiBins + phiBin] += boost::get<1>(sample); if (boost::get<1>(sample) > 0 && boost::get<2>(sample) == EDiscrete) discreteDirections.insert(boost::get<0>(sample)); @@ -123,9 +123,9 @@ void ChiSquare::fill( Float pdf = pdfFn(direction, EDiscrete); int thetaBin = std::min(std::max(0, - floorToInt(sphCoords.x * factor.x)), m_thetaBins-1); + math::floorToInt(sphCoords.x * factor.x)), m_thetaBins-1); int phiBin = std::min(std::max(0, - floorToInt(sphCoords.y * factor.y)), m_phiBins-1); + math::floorToInt(sphCoords.y * factor.y)), m_phiBins-1); m_refTable[thetaBin * m_phiBins + phiBin] += pdf * m_sampleCount; } diff --git a/src/libcore/math.cpp b/src/libcore/math.cpp new file mode 100644 index 00000000..a1709e86 --- /dev/null +++ b/src/libcore/math.cpp @@ -0,0 +1,146 @@ +/* + This file is part of Mitsuba, a physically based rendering system. + + Copyright (c) 2007-2012 by Wenzel Jakob and others. + + Mitsuba is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License Version 3 + as published by the Free Software Foundation. + + Mitsuba is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +*/ + +#include + +MTS_NAMESPACE_BEGIN + +namespace math { + +Float erfinv(Float x) { + // Based on "Approximating the erfinv function" by Mark Giles + Float w = -math::fastlog(((Float) 1 - x)*((Float) 1 + x)); + Float p; + if (w < (Float) 5) { + w = w - (Float) 2.5; + p = (Float) 2.81022636e-08; + p = (Float) 3.43273939e-07 + p*w; + p = (Float) -3.5233877e-06 + p*w; + p = (Float) -4.39150654e-06 + p*w; + p = (Float) 0.00021858087 + p*w; + p = (Float) -0.00125372503 + p*w; + p = (Float) -0.00417768164 + p*w; + p = (Float) 0.246640727 + p*w; + p = (Float) 1.50140941 + p*w; + } else { + w = std::sqrt(w) - (Float) 3; + p = (Float) -0.000200214257; + p = (Float) 0.000100950558 + p*w; + p = (Float) 0.00134934322 + p*w; + p = (Float) -0.00367342844 + p*w; + p = (Float) 0.00573950773 + p*w; + p = (Float) -0.0076224613 + p*w; + p = (Float) 0.00943887047 + p*w; + p = (Float) 1.00167406 + p*w; + p = (Float) 2.83297682 + p*w; + } + return p*x; +} + +Float erf(Float x) { + Float a1 = (Float) 0.254829592; + Float a2 = (Float) -0.284496736; + Float a3 = (Float) 1.421413741; + Float a4 = (Float) -1.453152027; + Float a5 = (Float) 1.061405429; + Float p = (Float) 0.3275911; + + // Save the sign of x + Float sign = math::signum(x); + x = std::abs(x); + + // A&S formula 7.1.26 + Float t = (Float) 1.0 / ((Float) 1.0 + p*x); + Float y = (Float) 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*math::fastexp(-x*x); + + return sign*y; +} + +float hypot2(float a, float b) { + float r; + if (std::abs(a) > std::abs(b)) { + r = b / a; + r = std::abs(a) * std::sqrt(1.0f + r*r); + } else if (b != 0.0f) { + r = a / b; + r = std::abs(b) * std::sqrt(1.0f + r*r); + } else { + r = 0.0f; + } + return r; +} + + +double hypot2(double a, double b) { + double r; + if (std::abs(a) > std::abs(b)) { + r = b / a; + r = std::abs(a) * std::sqrt(1.0 + r*r); + } else if (b != 0.0) { + r = a / b; + r = std::abs(b) * std::sqrt(1.0 + r*r); + } else { + r = 0.0; + } + return r; +} + +float log2(float value) { + const float invLn2 = 1.0f / std::log(2.0f); + return fastlog(value) * invLn2; +} + +double log2(double value) { + const double invLn2 = 1.0 / std::log(2.0); + return fastlog(value) * invLn2; +} + +int log2i(uint32_t value) { + int r = 0; + while ((value >> r) != 0) + r++; + return r-1; +} + +int log2i(uint64_t value) { + int r = 0; + while ((value >> r) != 0) + r++; + return r-1; +} + +/* Fast rounding & power-of-two test algorithms from PBRT */ +uint32_t roundToPowerOfTwo(uint32_t i) { + i--; + i |= i >> 1; i |= i >> 2; + i |= i >> 4; i |= i >> 8; + i |= i >> 16; + return i+1; +} + +uint64_t roundToPowerOfTwo(uint64_t i) { + i--; + i |= i >> 1; i |= i >> 2; + i |= i >> 4; i |= i >> 8; + i |= i >> 16; i |= i >> 32; + return i+1; +} + +}; + +MTS_NAMESPACE_END diff --git a/src/libcore/spectrum.cpp b/src/libcore/spectrum.cpp index f4466f32..9acdfb3b 100644 --- a/src/libcore/spectrum.cpp +++ b/src/libcore/spectrum.cpp @@ -196,7 +196,7 @@ Float Spectrum::eval(Float lambda) const { "is configured for RGB-based rendering"); return 0.0f; #else - int index = floorToInt((lambda - SPECTRUM_MIN_WAVELENGTH) * + int index = math::floorToInt((lambda - SPECTRUM_MIN_WAVELENGTH) * ((Float) SPECTRUM_SAMPLES / (Float) SPECTRUM_RANGE)); if (index < 0 || index >= SPECTRUM_SAMPLES) @@ -676,8 +676,8 @@ Float InterpolatedSpectrum::average(Float lambdaMin, Float lambdaMax) const { if (cb <= ca) continue; - Float interpA = lerp((ca - a) * invAB, fa, fb); - Float interpB = lerp((cb - a) * invAB, fa, fb); + Float interpA = math::lerp((ca - a) * invAB, fa, fb); + Float interpB = math::lerp((cb - a) * invAB, fa, fb); result += 0.5f * (interpA + interpB) * (cb-ca); } @@ -703,7 +703,7 @@ Float InterpolatedSpectrum::eval(Float lambda) const { b = m_wavelengths[idx1], fa = m_values[idx1-1], fb = m_values[idx1]; - return lerp((lambda - a) / (b-a), fb, fa); + return math::lerp((lambda - a) / (b-a), fb, fa); } else if (idx2 == idx1+1) { /* Hit a value exactly */ return m_values[idx1]; diff --git a/src/libcore/triangle.cpp b/src/libcore/triangle.cpp index a1965bd2..3aa40c13 100644 --- a/src/libcore/triangle.cpp +++ b/src/libcore/triangle.cpp @@ -27,7 +27,7 @@ Point Triangle::sample(const Point *positions, const Normal *normals, const Point &p1 = positions[idx[1]]; const Point &p2 = positions[idx[2]]; - Point2 bary = Warp::squareToUniformTriangle(sample); + Point2 bary = warp::squareToUniformTriangle(sample); Vector sideA = p1 - p0, sideB = p2 - p0; Point p = p0 + (sideA * bary.x) + (sideB * bary.y); diff --git a/src/libcore/util.cpp b/src/libcore/util.cpp index 20179d57..401abce5 100644 --- a/src/libcore/util.cpp +++ b/src/libcore/util.cpp @@ -399,11 +399,6 @@ std::string getFQDN() { return fqdn; } -Float log2(Float value) { - const Float invLn2 = (Float) 1.0f / math::fastlog((Float) 2.0f); - return math::fastlog(value) * invLn2; -} - std::string formatString(const char *fmt, ...) { char tmp[512]; va_list iterator; @@ -444,37 +439,6 @@ std::string formatString(const char *fmt, ...) { return std::string(tmp); } -int log2i(uint32_t value) { - int r = 0; - while ((value >> r) != 0) - r++; - return r-1; -} - -int log2i(uint64_t value) { - int r = 0; - while ((value >> r) != 0) - r++; - return r-1; -} - -/* Fast rounding & power-of-two test algorithms from PBRT */ -uint32_t roundToPowerOfTwo(uint32_t i) { - i--; - i |= i >> 1; i |= i >> 2; - i |= i >> 4; i |= i >> 8; - i |= i >> 16; - return i+1; -} - -uint64_t roundToPowerOfTwo(uint64_t i) { - i--; - i |= i >> 1; i |= i >> 2; - i |= i >> 4; i |= i >> 8; - i |= i >> 16; i |= i >> 32; - return i+1; -} - // ----------------------------------------------------------------------- // Numerical utility functions // ----------------------------------------------------------------------- @@ -910,18 +874,4 @@ std::string memString(size_t size, bool precise) { return os.str(); } -Float hypot2(Float a, Float b) { - Float r; - if (std::abs(a) > std::abs(b)) { - r = b / a; - r = std::abs(a) * std::sqrt(1 + r*r); - } else if (b != 0) { - r = a / b; - r = std::abs(b) * std::sqrt(1 + r*r); - } else { - r = 0; - } - return r; -} - MTS_NAMESPACE_END diff --git a/src/libcore/vmf.cpp b/src/libcore/vmf.cpp index 59f2dcf7..a2a6aaed 100644 --- a/src/libcore/vmf.cpp +++ b/src/libcore/vmf.cpp @@ -38,7 +38,7 @@ Float VonMisesFisherDistr::eval(Float cosTheta) const { Vector VonMisesFisherDistr::sample(const Point2 &sample) const { if (m_kappa == 0) - return Warp::squareToUniformSphere(sample); + return warp::squareToUniformSphere(sample); #if 0 Float cosTheta = math::fastlog(math::fastexp(-m_kappa) + 2 * diff --git a/src/libcore/warp.cpp b/src/libcore/warp.cpp index ef54171b..8cb07fb2 100644 --- a/src/libcore/warp.cpp +++ b/src/libcore/warp.cpp @@ -20,7 +20,9 @@ MTS_NAMESPACE_BEGIN -Vector Warp::squareToUniformSphere(const Point2 &sample) { +namespace warp { + +Vector squareToUniformSphere(const Point2 &sample) { Float z = 1.0f - 2.0f * sample.y; Float r = math::safe_sqrt(1.0f - z*z); Float sinPhi, cosPhi; @@ -28,7 +30,7 @@ Vector Warp::squareToUniformSphere(const Point2 &sample) { return Vector(r * cosPhi, r * sinPhi, z); } -Vector Warp::squareToUniformHemisphere(const Point2 &sample) { +Vector squareToUniformHemisphere(const Point2 &sample) { Float z = sample.x; Float tmp = math::safe_sqrt(1.0f - z*z); @@ -38,8 +40,8 @@ Vector Warp::squareToUniformHemisphere(const Point2 &sample) { return Vector(cosPhi * tmp, sinPhi * tmp, z); } -Vector Warp::squareToCosineHemisphere(const Point2 &sample) { - Point2 p = Warp::squareToUniformDiskConcentric(sample); +Vector squareToCosineHemisphere(const Point2 &sample) { + Point2 p = squareToUniformDiskConcentric(sample); Float z = math::safe_sqrt(1.0f - p.x*p.x - p.y*p.y); /* Guard against numerical imprecisions */ @@ -49,7 +51,7 @@ Vector Warp::squareToCosineHemisphere(const Point2 &sample) { return Vector(p.x, p.y, z); } -Vector Warp::squareToUniformCone(Float cosCutoff, const Point2 &sample) { +Vector squareToUniformCone(Float cosCutoff, const Point2 &sample) { Float cosTheta = (1-sample.x) + sample.x * cosCutoff; Float sinTheta = math::safe_sqrt(1.0f - cosTheta * cosTheta); @@ -60,7 +62,7 @@ Vector Warp::squareToUniformCone(Float cosCutoff, const Point2 &sample) { sinPhi * sinTheta, cosTheta); } -Point2 Warp::squareToUniformDisk(const Point2 &sample) { +Point2 squareToUniformDisk(const Point2 &sample) { Float r = std::sqrt(sample.x); Float sinPhi, cosPhi; math::sincos(2.0f * M_PI * sample.y, &sinPhi, &cosPhi); @@ -71,12 +73,12 @@ Point2 Warp::squareToUniformDisk(const Point2 &sample) { ); } -Point2 Warp::squareToUniformTriangle(const Point2 &sample) { +Point2 squareToUniformTriangle(const Point2 &sample) { Float a = math::safe_sqrt(1.0f - sample.x); return Point2(1 - a, a * sample.y); } -Point2 Warp::squareToUniformDiskConcentric(const Point2 &sample) { +Point2 squareToUniformDiskConcentric(const Point2 &sample) { Float r1 = 2.0f*sample.x - 1.0f; Float r2 = 2.0f*sample.y - 1.0f; @@ -99,7 +101,7 @@ Point2 Warp::squareToUniformDiskConcentric(const Point2 &sample) { return Point2(r * cosPhi, r * sinPhi); } -Point2 Warp::uniformDiskToSquareConcentric(const Point2 &p) { +Point2 uniformDiskToSquareConcentric(const Point2 &p) { Float r = std::sqrt(p.x * p.x + p.y * p.y), phi = std::atan2(p.y, p.x), a, b; @@ -126,7 +128,7 @@ Point2 Warp::uniformDiskToSquareConcentric(const Point2 &p) { return Point2(0.5f * (a+1), 0.5f * (b+1)); } -Point2 Warp::squareToStdNormal(const Point2 &sample) { +Point2 squareToStdNormal(const Point2 &sample) { Float r = std::sqrt(-2 * math::fastlog(1-sample.x)), phi = 2 * M_PI * sample.y; Point2 result; @@ -134,7 +136,7 @@ Point2 Warp::squareToStdNormal(const Point2 &sample) { return result * r; } -Float Warp::squareToStdNormalPdf(const Point2 &pos) { +Float squareToStdNormalPdf(const Point2 &pos) { return INV_TWOPI * math::fastexp(-(pos.x*pos.x + pos.y*pos.y)/2.0f); } @@ -152,14 +154,14 @@ static Float intervalToTent(Float sample) { return sign * (1 - std::sqrt(sample)); } -Point2 Warp::squareToTent(const Point2 &sample) { +Point2 squareToTent(const Point2 &sample) { return Point2( intervalToTent(sample.x), intervalToTent(sample.y) ); } -Float Warp::intervalToNonuniformTent(Float a, Float b, Float c, Float sample) { +Float intervalToNonuniformTent(Float a, Float b, Float c, Float sample) { Float factor; if (sample * (c-a) < b-a) { @@ -173,4 +175,6 @@ Float Warp::intervalToNonuniformTent(Float a, Float b, Float c, Float sample) { return b + factor * (1-math::safe_sqrt(sample)); } +}; + MTS_NAMESPACE_END diff --git a/src/libhw/gltexture.cpp b/src/libhw/gltexture.cpp index 26d4783e..d357eb84 100644 --- a/src/libhw/gltexture.cpp +++ b/src/libhw/gltexture.cpp @@ -121,7 +121,7 @@ void GLTexture::init() { glFramebufferTexture2DEXT(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, m_glType, m_id, 0); } else if (m_type == ETextureCubeMap) { - Assert(m_size.x == m_size.y && isPowerOfTwo(m_size.x)); + Assert(m_size.x == m_size.y && math::isPowerOfTwo(m_size.x)); Assert(m_fbType == EColorBuffer); Assert(m_samples == 1); @@ -178,7 +178,7 @@ void GLTexture::init() { glFramebufferTexture2DEXT(GL_FRAMEBUFFER_EXT, GL_DEPTH_ATTACHMENT_EXT, m_glType, m_id, 0); } else if (m_type == ETextureCubeMap) { - Assert(m_size.x == m_size.y && isPowerOfTwo(m_size.x)); + Assert(m_size.x == m_size.y && math::isPowerOfTwo(m_size.x)); for (int i=0; i<6; i++) glTexImage2D(GL_TEXTURE_CUBE_MAP_POSITIVE_X + i, 0, m_internalFormat, m_size.x, m_size.y, 0, m_format, m_dataFormat, NULL); @@ -260,8 +260,8 @@ void GLTexture::refresh() { glPixelStorei(GL_UNPACK_ALIGNMENT, 1); if (m_type == ETexture1D) { - Assert((isPowerOfTwo(m_size.x) && m_size.y == 1) - || (isPowerOfTwo(m_size.y) && m_size.x == 1)); + Assert((math::isPowerOfTwo(m_size.x) && m_size.y == 1) + || (math::isPowerOfTwo(m_size.y) && m_size.x == 1)); if (isMipMapped()) { /* Let GLU generate mipmaps for us */ @@ -282,7 +282,7 @@ void GLTexture::refresh() { } else if (m_type == ETextureCubeMap) { Assert(bitmap != NULL); Assert(bitmap->getWidth() == bitmap->getHeight()); - Assert(isPowerOfTwo(bitmap->getWidth())); + Assert(math::isPowerOfTwo(bitmap->getWidth())); if (isMipMapped()) glTexParameteri(m_glType, GL_GENERATE_MIPMAP, GL_TRUE); diff --git a/src/libhw/vpl.cpp b/src/libhw/vpl.cpp index cf62fe76..e5e5bfee 100644 --- a/src/libhw/vpl.cpp +++ b/src/libhw/vpl.cpp @@ -319,7 +319,7 @@ void VPLShaderManager::setVPL(const VPL &vpl) { if (vpl.type == ESurfaceVPL || (vpl.type == EPointEmitterVPL && vpl.emitter->getType() & Emitter::EOnSurface)) { ray = Ray(vpl.its.p, vpl.its.shFrame.toWorld( - Warp::squareToCosineHemisphere(sample)), 0); + warp::squareToCosineHemisphere(sample)), 0); #if defined(MTS_VPL_USE_PARABOLOID_MAPS) m_shadowMapType = ShadowMapGenerator::EParaboloid; @@ -329,7 +329,7 @@ void VPLShaderManager::setVPL(const VPL &vpl) { m_shadowMapType = ShadowMapGenerator::EHemicube; #endif } else if (vpl.type == EPointEmitterVPL) { - ray = Ray(vpl.its.p, Warp::squareToUniformSphere(sample), 0); + ray = Ray(vpl.its.p, warp::squareToUniformSphere(sample), 0); #if defined(MTS_VPL_USE_SINGLE_PASS) m_shadowMapType = ShadowMapGenerator::ECubeSinglePass; #else diff --git a/src/libpython/core.cpp b/src/libpython/core.cpp index 74e2473f..648c925f 100644 --- a/src/libpython/core.cpp +++ b/src/libpython/core.cpp @@ -2352,7 +2352,6 @@ void export_core() { bp::def("getTotalSystemMemory", &getTotalSystemMemory); bp::def("getFQDN", &getFQDN); bp::def("rdtsc", &rdtsc); - bp::def("hypot2", &hypot2); /* Functions from qmc.h */ bp::def("radicalInverse2Single", radicalInverse2Single); diff --git a/src/librender/integrator.cpp b/src/librender/integrator.cpp index 815987ca..8cb9d9c4 100644 --- a/src/librender/integrator.cpp +++ b/src/librender/integrator.cpp @@ -75,7 +75,7 @@ Spectrum SamplingIntegrator::E(const Scene *scene, const Intersection &its, /* Sample the indirect illumination component */ if (handleIndirect) { query.newQuery(RadianceQueryRecord::ERadianceNoEmission, medium); - Vector d = frame.toWorld(Warp::squareToCosineHemisphere(query.nextSample2D())); + Vector d = frame.toWorld(warp::squareToCosineHemisphere(query.nextSample2D())); ++query.depth; query.medium = medium; E += Li(RayDifferential(its.p, d, its.time), query) * M_PI; diff --git a/src/librender/noise.cpp b/src/librender/noise.cpp index a012b303..b593e34a 100644 --- a/src/librender/noise.cpp +++ b/src/librender/noise.cpp @@ -63,9 +63,9 @@ inline static Float noiseWeight(Float t) { Float Noise::perlinNoise(const Point &p) { // Compute noise cell coordinates and offsets - int ix = floorToInt(p.x), - iy = floorToInt(p.y), - iz = floorToInt(p.z); + int ix = math::floorToInt(p.x), + iy = math::floorToInt(p.y), + iz = math::floorToInt(p.z); Float dx = p.x - ix, dy = p.y - iy, @@ -90,21 +90,21 @@ Float Noise::perlinNoise(const Point &p) { wy = noiseWeight(dy), wz = noiseWeight(dz); - Float x00 = lerp(wx, w000, w100), - x10 = lerp(wx, w010, w110), - x01 = lerp(wx, w001, w101), - x11 = lerp(wx, w011, w111), - y0 = lerp(wy, x00, x10), - y1 = lerp(wy, x01, x11); + Float x00 = math::lerp(wx, w000, w100), + x10 = math::lerp(wx, w010, w110), + x01 = math::lerp(wx, w001, w101), + x11 = math::lerp(wx, w011, w111), + y0 = math::lerp(wy, x00, x10), + y1 = math::lerp(wy, x01, x11); - return lerp(wz, y0, y1); + return math::lerp(wz, y0, y1); } Float Noise::fbm(const Point &p, const Vector &dpdx, const Vector &dpdy, Float omega, int maxOctaves) { // Compute number of octaves for antialiased FBm Float s2 = std::max(dpdx.lengthSquared(), dpdy.lengthSquared()); - Float foctaves = std::min((Float) maxOctaves, 1.f - .5f * log2(s2)); + Float foctaves = std::min((Float) maxOctaves, 1.f - .5f * math::log2(s2)); int octaves = (int) foctaves; // Compute sum of octaves of noise for FBm @@ -115,7 +115,7 @@ Float Noise::fbm(const Point &p, const Vector &dpdx, o *= omega; } Float partialOctave = foctaves - octaves; - sum += o * smoothStep(.3f, .7f, partialOctave) + sum += o * math::smoothStep((Float) .3f, (Float) .7f, partialOctave) * perlinNoise(lambda * p); return sum; } @@ -124,7 +124,7 @@ Float Noise::turbulence(const Point &p, const Vector &dpdx, const Vector &dpdy, Float omega, int maxOctaves) { // Compute number of octaves for antialiased FBm Float s2 = std::max(dpdx.lengthSquared(), dpdy.lengthSquared()); - Float foctaves = std::min((Float) maxOctaves, 1.f - .5f * log2(s2)); + Float foctaves = std::min((Float) maxOctaves, 1.f - .5f * math::log2(s2)); int octaves = (int) foctaves; // Compute sum of octaves of noise for turbulence @@ -135,7 +135,7 @@ Float Noise::turbulence(const Point &p, const Vector &dpdx, o *= omega; } Float partialOctave = foctaves - octaves; - sum += o * smoothStep(.3f, .7f, partialOctave) + sum += o * math::smoothStep((Float) .3f, (Float) .7f, partialOctave) * std::abs(perlinNoise(lambda * p)); return sum; } diff --git a/src/librender/vpl.cpp b/src/librender/vpl.cpp index 474248fc..dac763a6 100644 --- a/src/librender/vpl.cpp +++ b/src/librender/vpl.cpp @@ -141,7 +141,7 @@ size_t generateVPLs(const Scene *scene, Random *random, appendVPL(scene, random, lumVPL, false, vpls); dRec.d = -diRec.d; - Point2 offset = Warp::squareToUniformDiskConcentric(sampler->next2D()); + Point2 offset = warp::squareToUniformDiskConcentric(sampler->next2D()); Vector perpOffset = Frame(diRec.d).toWorld(Vector(offset.x, offset.y, 0)); BSphere geoBSphere = scene->getKDTree()->getAABB().getBSphere(); pRec.p = geoBSphere.center + (perpOffset - dRec.d) * geoBSphere.radius; diff --git a/src/mtsgui/mainwindow.cpp b/src/mtsgui/mainwindow.cpp index e11f2f18..70c21329 100644 --- a/src/mtsgui/mainwindow.cpp +++ b/src/mtsgui/mainwindow.cpp @@ -1080,17 +1080,17 @@ void MainWindow::on_glView_crop(int type, int x, int y, int width, int height) { int maxNewSize = std::max(width, height); Float magnification = maxOldSize / (Float) maxNewSize; - width = floorToInt(width*magnification); - height = floorToInt(height*magnification); + width = math::floorToInt(width*magnification); + height = math::floorToInt(height*magnification); - filmProps.setInteger("cropOffsetX", floorToInt(magnification + filmProps.setInteger("cropOffsetX", math::floorToInt(magnification * (x + oldCropOffset.x)), false); - filmProps.setInteger("cropOffsetY", floorToInt(magnification + filmProps.setInteger("cropOffsetY", math::floorToInt(magnification * (y + oldCropOffset.y)), false); filmProps.setInteger("cropWidth", width); filmProps.setInteger("cropHeight", height); - filmProps.setInteger("width", ceilToInt(oldSize.x*magnification), false); - filmProps.setInteger("height", ceilToInt(oldSize.y*magnification), false); + filmProps.setInteger("width", math::ceilToInt(oldSize.x*magnification), false); + filmProps.setInteger("height", math::ceilToInt(oldSize.y*magnification), false); } else { width = context->originalSize.x; height = context->originalSize.y; diff --git a/src/mtsgui/rendersettingsdlg.cpp b/src/mtsgui/rendersettingsdlg.cpp index c99f3d8a..d7cdf63a 100644 --- a/src/mtsgui/rendersettingsdlg.cpp +++ b/src/mtsgui/rendersettingsdlg.cpp @@ -474,11 +474,11 @@ void RenderSettingsDialog::apply(SceneContext *ctx) { Vector2i oldCropSize = oldFilm->getCropSize(); Point2i oldCropOffset = oldFilm->getCropOffset(); - Vector2i size(std::ceil((oldSize.x * cropSize.x) / (Float) oldCropSize.x), - std::ceil((oldSize.y * cropSize.y) / (Float) oldCropSize.y)); + Vector2i size(math::roundToInt((oldSize.x * cropSize.x / (Float) oldCropSize.x)), + math::roundToInt((oldSize.y * cropSize.y / (Float) oldCropSize.y))); - Point2i cropOffset(std::floor((oldCropOffset.x * cropSize.x) / (Float) oldCropSize.x), - std::floor((oldCropOffset.y * cropSize.y) / (Float) oldCropSize.y)); + Point2i cropOffset(math::roundToInt((oldCropOffset.x * cropSize.x / (Float) oldCropSize.x)), + math::roundToInt((oldCropOffset.y * cropSize.y / (Float) oldCropSize.y))); filmProps.setInteger("width", size.x, false); filmProps.setInteger("height", size.y, false); diff --git a/src/mtsgui/simdtonemap.cpp b/src/mtsgui/simdtonemap.cpp index 21bf4c25..89979682 100644 --- a/src/mtsgui/simdtonemap.cpp +++ b/src/mtsgui/simdtonemap.cpp @@ -93,6 +93,8 @@ inline mitsuba::math::SSEVector4f am_log(const mitsuba::math::SSEVector4f& x) { typedef mitsuba::math::SSEVector4f v4f; typedef mitsuba::math::SSEVector4i v4i; using mitsuba::math::castAsFloat; + using mitsuba::math::clamp_ps; + using mitsuba::math::fastpow_ps; // Constants const v4f min_normal(castAsFloat(v4i::constant<0x00800000>())); @@ -280,7 +282,7 @@ void tonemap(const PixelRGBA32FGroup* const begin, buffer[aOffset+i] = it->p[3]; transpose(buffer[rOffset+i], buffer[gOffset+i], buffer[bOffset+i], buffer[aOffset+i]); - buffer[aOffset+i] = clamp(buffer[aOffset+i], + buffer[aOffset+i] = clamp_ps(buffer[aOffset+i], V4f::zero(), const_1f); } @@ -299,12 +301,12 @@ void tonemap(const PixelRGBA32FGroup* const begin, // Clamp and nonlinear display transform of the color components if (displayMethod == ESRGB) { for (ptrdiff_t i = 0; i != numColorIter; ++i) { - const V4f s = clamp(buffer[i], V4f::zero(), const_1f); + const V4f s = clamp_ps(buffer[i], V4f::zero(), const_1f); buffer[i] = sRGB_Remez43(s); } } else if (displayMethod == EGamma) { for (ptrdiff_t i = 0; i != numColorIter; ++i) { - const V4f s = clamp(buffer[i], V4f::zero(), const_1f); + const V4f s = clamp_ps(buffer[i], V4f::zero(), const_1f); buffer[i] = fastpow_ps(s, invGamma); } } @@ -534,9 +536,13 @@ bool luminance(const mitsuba::Bitmap* source, const float multiplier, result.maxLuminance = max(rTail.maxLuminance, result.maxLuminance); } - outMaxLuminance = hmax(result.maxLuminance); - const float sumLogLuminance = hsum(result.sumLogLuminance); - outAvgLogLuminance = mitsuba::math::fastexp(sumLogLuminance / pixelCount); + using mitsuba::math::hmax_ps; + using mitsuba::math::hsum_ps; + using mitsuba::math::fastexp; + + outMaxLuminance = hmax_ps(result.maxLuminance); + const float sumLogLuminance = hsum_ps(result.sumLogLuminance); + outAvgLogLuminance = fastexp(sumLogLuminance / pixelCount); return true; }; diff --git a/src/phase/isotropic.cpp b/src/phase/isotropic.cpp index e753598e..63d6f667 100644 --- a/src/phase/isotropic.cpp +++ b/src/phase/isotropic.cpp @@ -62,19 +62,19 @@ public: Float sample(PhaseFunctionSamplingRecord &pRec, Sampler *sampler) const { Point2 sample(sampler->next2D()); - pRec.wo = Warp::squareToUniformSphere(sample); + pRec.wo = warp::squareToUniformSphere(sample); return 1.0f; } Float sample(PhaseFunctionSamplingRecord &pRec, Float &pdf, Sampler *sampler) const { - pRec.wo = Warp::squareToUniformSphere(sampler->next2D()); - pdf = Warp::squareToUniformSpherePdf(); + pRec.wo = warp::squareToUniformSphere(sampler->next2D()); + pdf = warp::squareToUniformSpherePdf(); return 1.0f; } Float eval(const PhaseFunctionSamplingRecord &pRec) const { - return Warp::squareToUniformSpherePdf(); + return warp::squareToUniformSpherePdf(); } Float getMeanCosine() const { diff --git a/src/phase/kkay.cpp b/src/phase/kkay.cpp index 42c01169..d34ba60e 100644 --- a/src/phase/kkay.cpp +++ b/src/phase/kkay.cpp @@ -84,19 +84,19 @@ public: Float sample(PhaseFunctionSamplingRecord &pRec, Sampler *sampler) const { - pRec.wo = Warp::squareToUniformSphere(sampler->next2D()); + pRec.wo = warp::squareToUniformSphere(sampler->next2D()); return eval(pRec) * (4 * M_PI); } Float sample(PhaseFunctionSamplingRecord &pRec, Float &pdf, Sampler *sampler) const { - pRec.wo = Warp::squareToUniformSphere(sampler->next2D()); - pdf = Warp::squareToUniformSpherePdf(); + pRec.wo = warp::squareToUniformSphere(sampler->next2D()); + pdf = warp::squareToUniformSpherePdf(); return eval(pRec) * (4 * M_PI); } Float pdf(const PhaseFunctionSamplingRecord &pRec) const { - return Warp::squareToUniformSpherePdf(); + return warp::squareToUniformSpherePdf(); } Float eval(const PhaseFunctionSamplingRecord &pRec) const { diff --git a/src/phase/microflake.cpp b/src/phase/microflake.cpp index 52318ae3..9e0518d9 100644 --- a/src/phase/microflake.cpp +++ b/src/phase/microflake.cpp @@ -129,7 +129,7 @@ public: if (pRec.mRec.orientation.isZero()) { /* What to do when the local orientation is undefined */ #if 0 - pRec.wo = Warp::squareToUniformSphere(sampler->next2D()); + pRec.wo = warp::squareToUniformSphere(sampler->next2D()); return 1.0f; #else return 0.0f; diff --git a/src/samplers/halton.cpp b/src/samplers/halton.cpp index c44d05df..cf5da2f8 100644 --- a/src/samplers/halton.cpp +++ b/src/samplers/halton.cpp @@ -226,21 +226,15 @@ public: y = x_ - (int64_t) (d * y_); } - /// Positive modulo for 64-bit integers - uint64_t modulo(int64_t a, int64_t b) { - int64_t result = a - (a/b) * b; - return (uint64_t) ((result < 0) ? result+b : result); - } - /** - * \brief Compute the multiplicative inverse of a modulo n, + * \brief Compute the multiplicative inverse of a math::modulo n, * where a and n a re relative prime. The result is in the * range 1, ..., n - 1. */ uint64_t multiplicativeInverse(int64_t a, int64_t n) { int64_t x, y; extendedGCD(a, n, x, y); - return modulo(x, n); + return math::modulo(x, n); } void setFilmResolution(const Vector2i &res, bool blocked) { diff --git a/src/samplers/hammersley.cpp b/src/samplers/hammersley.cpp index cd8b9490..a0454a6e 100644 --- a/src/samplers/hammersley.cpp +++ b/src/samplers/hammersley.cpp @@ -184,8 +184,8 @@ public: dimensions. This is required to support blocked rendering. */ for (int i=0; i<2; ++i) m_resolution[i] = std::min((uint32_t) MAX_RESOLUTION, - roundToPowerOfTwo((uint32_t) res[i])); - m_logHeight = log2i((uint32_t) m_resolution.y); + math::roundToPowerOfTwo((uint32_t) res[i])); + m_logHeight = math::log2i((uint32_t) m_resolution.y); m_samplesPerBatch = m_sampleCount; m_factor = (Float) 1.0f / (m_sampleCount * diff --git a/src/samplers/ldsampler.cpp b/src/samplers/ldsampler.cpp index 3a155d4f..3ab50fc3 100644 --- a/src/samplers/ldsampler.cpp +++ b/src/samplers/ldsampler.cpp @@ -80,8 +80,8 @@ public: /* Dimension, up to which which low discrepancy samples are guaranteed to be available. */ m_maxDimension = props.getInteger("dimension", 4); - if (!isPowerOfTwo(m_sampleCount)) { - m_sampleCount = roundToPowerOfTwo(m_sampleCount); + if (!math::isPowerOfTwo(m_sampleCount)) { + m_sampleCount = math::roundToPowerOfTwo(m_sampleCount); Log(EWarn, "Sample count should be a power of two -- rounding to " SIZE_T_FMT, m_sampleCount); } diff --git a/src/samplers/sobol.cpp b/src/samplers/sobol.cpp index 0d4d2cdd..ac3908b3 100644 --- a/src/samplers/sobol.cpp +++ b/src/samplers/sobol.cpp @@ -149,11 +149,11 @@ public: m_resolution = 1; m_logResolution = 0; } else { - uint32_t resolution = roundToPowerOfTwo( + uint32_t resolution = math::roundToPowerOfTwo( (uint32_t) std::max(res.x, res.y)); m_resolution = (Float) resolution; - m_logResolution = log2i(resolution); + m_logResolution = math::log2i(resolution); } } diff --git a/src/sensors/fluencemeter.cpp b/src/sensors/fluencemeter.cpp index 6bd5a084..21bb4c0b 100644 --- a/src/sensors/fluencemeter.cpp +++ b/src/sensors/fluencemeter.cpp @@ -92,7 +92,7 @@ public: const Transform &trafo = m_worldTransform->eval(ray.time); ray.setOrigin(trafo(Point(0.0f))); - ray.setDirection(trafo(Warp::squareToUniformSphere(pixelSample))); + ray.setDirection(trafo(warp::squareToUniformSphere(pixelSample))); return Spectrum(1.0f); } @@ -118,7 +118,7 @@ public: PositionSamplingRecord &pRec, const Point2 &sample, const Point2 *extra) const { - dRec.d = Warp::squareToUniformSphere(sample); + dRec.d = warp::squareToUniformSphere(sample); dRec.pdf = INV_FOURPI; dRec.measure = ESolidAngle; return Spectrum(1.0f); diff --git a/src/sensors/irradiancemeter.cpp b/src/sensors/irradiancemeter.cpp index 9ef414e5..c136291e 100644 --- a/src/sensors/irradiancemeter.cpp +++ b/src/sensors/irradiancemeter.cpp @@ -115,7 +115,7 @@ public: ray.setOrigin(pRec.p); ray.setDirection(Frame(pRec.n).toWorld( - Warp::squareToCosineHemisphere(otherSample))); + warp::squareToCosineHemisphere(otherSample))); return Spectrum(M_PI); } @@ -149,9 +149,9 @@ public: Spectrum sampleDirection(DirectionSamplingRecord &dRec, PositionSamplingRecord &pRec, const Point2 &sample, const Point2 *extra) const { - Vector local = Warp::squareToCosineHemisphere(sample); + Vector local = warp::squareToCosineHemisphere(sample); dRec.d = Frame(pRec.n).toWorld(local); - dRec.pdf = Warp::squareToCosineHemispherePdf(local); + dRec.pdf = warp::squareToCosineHemispherePdf(local); dRec.measure = ESolidAngle; return Spectrum(1.0f); } diff --git a/src/sensors/spherical.cpp b/src/sensors/spherical.cpp index e39c2143..3abc473d 100644 --- a/src/sensors/spherical.cpp +++ b/src/sensors/spherical.cpp @@ -155,7 +155,7 @@ public: Vector d = normalize(m_worldTransform->eval(pRec.time).inverse()(dRec.d)); samplePosition = Point2( - modulo(std::atan2(d.x, -d.z) * INV_TWOPI, (Float) 1) * m_resolution.x, + math::modulo(std::atan2(d.x, -d.z) * INV_TWOPI, (Float) 1) * m_resolution.x, math::safe_acos(d.y) * INV_PI * m_resolution.y ); @@ -173,7 +173,7 @@ public: d *= invDist; dRec.uv = Point2( - modulo(std::atan2(d.x, -d.z) * INV_TWOPI, (Float) 1) * m_resolution.x, + math::modulo(std::atan2(d.x, -d.z) * INV_TWOPI, (Float) 1) * m_resolution.x, math::safe_acos(d.y) * INV_PI * m_resolution.y ); diff --git a/src/sensors/telecentric.cpp b/src/sensors/telecentric.cpp index 02bf566d..d1a55224 100644 --- a/src/sensors/telecentric.cpp +++ b/src/sensors/telecentric.cpp @@ -194,7 +194,7 @@ public: Spectrum sampleRay(Ray &ray, const Point2 &pixelSample, const Point2 &otherSample, Float timeSample) const { - Point2 diskSample = Warp::squareToUniformDiskConcentric(otherSample) + Point2 diskSample = warp::squareToUniformDiskConcentric(otherSample) * (m_apertureRadius / m_scale.x); ray.time = sampleTime(timeSample); @@ -223,7 +223,7 @@ public: Spectrum sampleRayDifferential(RayDifferential &ray, const Point2 &pixelSample, const Point2 &otherSample, Float timeSample) const { - Point2 diskSample = Warp::squareToUniformDiskConcentric(otherSample) + Point2 diskSample = warp::squareToUniformDiskConcentric(otherSample) * (m_apertureRadius / m_scale.x); ray.time = sampleTime(timeSample); @@ -285,7 +285,7 @@ public: double rand4 = (tmp2 & 0x3FFFFFF) * (1.0 / 0x3FFFFFF); #endif - Point2 aperturePos = Warp::squareToUniformDiskConcentric(Point2(rand1, rand2)) + Point2 aperturePos = warp::squareToUniformDiskConcentric(Point2(rand1, rand2)) * (m_apertureRadius / m_scale.x); Point2 samplePos(rand3, rand4); @@ -382,7 +382,7 @@ public: radius += apertureRadius; /* Sample the ray origin */ - Point2 disk = Warp::squareToUniformDiskConcentric(sample); + Point2 disk = warp::squareToUniformDiskConcentric(sample); Point diskP(disk.x*radius+localP.x, disk.y*radius+localP.y, 0.0f); /* Compute the intersection with the focal plane */ diff --git a/src/sensors/thinlens.cpp b/src/sensors/thinlens.cpp index 79cb1eb5..6045f9f0 100644 --- a/src/sensors/thinlens.cpp +++ b/src/sensors/thinlens.cpp @@ -292,7 +292,7 @@ public: Spectrum sampleRay(Ray &ray, const Point2 &pixelSample, const Point2 &otherSample, Float timeSample) const { - Point2 tmp = Warp::squareToUniformDiskConcentric(otherSample) + Point2 tmp = warp::squareToUniformDiskConcentric(otherSample) * m_apertureRadius; ray.time = sampleTime(timeSample); @@ -323,7 +323,7 @@ public: Spectrum sampleRayDifferential(RayDifferential &ray, const Point2 &pixelSample, const Point2 &otherSample, Float timeSample) const { - Point2 tmp = Warp::squareToUniformDiskConcentric(otherSample) + Point2 tmp = warp::squareToUniformDiskConcentric(otherSample) * m_apertureRadius; ray.time = sampleTime(timeSample); @@ -364,7 +364,7 @@ public: const Point2 &sample, const Point2 *extra) const { const Transform &trafo = m_worldTransform->eval(pRec.time); - Point2 aperturePos = Warp::squareToUniformDiskConcentric(sample) + Point2 aperturePos = warp::squareToUniformDiskConcentric(sample) * m_apertureRadius; pRec.p = trafo.transformAffine( @@ -449,7 +449,7 @@ public: } /* Sample a position on the aperture (in local coordinates) */ - Point2 tmp = Warp::squareToUniformDiskConcentric(sample) + Point2 tmp = warp::squareToUniformDiskConcentric(sample) * m_apertureRadius; Point apertureP(tmp.x, tmp.y, 0); @@ -495,7 +495,7 @@ public: const Point2 &aaSample) const { Float right = std::tan(m_xfov * M_PI/360) * m_nearClip, left = -right; Float top = right / m_aspect, bottom = -top; - Point2 apertureP = Warp::squareToUniformDiskConcentric(apertureSample) + Point2 apertureP = warp::squareToUniformDiskConcentric(apertureSample) * m_apertureRadius; Vector2 jitterScale( diff --git a/src/shapes/disk.cpp b/src/shapes/disk.cpp index 922f8e83..b2a45cd6 100644 --- a/src/shapes/disk.cpp +++ b/src/shapes/disk.cpp @@ -249,7 +249,7 @@ public: void samplePosition(PositionSamplingRecord &pRec, const Point2 &sample) const { const Transform &trafo = m_objectToWorld->eval(pRec.time); - Point2 p = Warp::squareToUniformDiskConcentric(sample); + Point2 p = warp::squareToUniformDiskConcentric(sample); pRec.p = trafo(Point3(p.x, p.y, 0)); pRec.n = trafo(normalize(Normal(0,0,1))); diff --git a/src/shapes/heightfield.cpp b/src/shapes/heightfield.cpp index fb6c61f8..72bb86e4 100644 --- a/src/shapes/heightfield.cpp +++ b/src/shapes/heightfield.cpp @@ -199,6 +199,15 @@ public: return (size_t) m_levelSize[0].x * (size_t) m_levelSize[0].y; } + inline static int signumToInt(Float value) { + if (value < 0) + return -1; + else if (value > 0) + return 1; + else + return 0; + } + bool rayIntersect(const Ray &_ray, Float mint, Float maxt, Float &t, void *tmp) const { StackEntry stack[MTS_QTREE_MAXDEPTH]; @@ -233,14 +242,14 @@ public: maxY = (int) std::max(enterPt.y, exitPt.y); /* Determine quadtree level */ - int level = clamp(1 + log2i( + int level = math::clamp(1 + math::log2i( std::max((uint32_t) (minX ^ maxX), (uint32_t) (minY ^ maxY))), 0, m_levelCount-1); /* Compute X and Y coordinates at that level */ const Vector2i &blockSize = m_blockSize[level]; - int x = clamp(minX / blockSize.x, 0, m_levelSize[level].x-1), - y = clamp(minY / blockSize.y, 0, m_levelSize[level].y-1); + int x = math::clamp(minX / blockSize.x, 0, m_levelSize[level].x-1), + y = math::clamp(minY / blockSize.y, 0, m_levelSize[level].y-1); stack[stackIdx].level = level; stack[stackIdx].x = x; @@ -504,8 +513,8 @@ public: m_dataSize = m_bitmap->getSize(); if (m_dataSize.x < 2) m_dataSize.x = 2; if (m_dataSize.y < 2) m_dataSize.y = 2; - if (!isPowerOfTwo(m_dataSize.x - 1)) m_dataSize.x = (int) roundToPowerOfTwo((uint32_t) m_dataSize.x - 1) + 1; - if (!isPowerOfTwo(m_dataSize.y - 1)) m_dataSize.y = (int) roundToPowerOfTwo((uint32_t) m_dataSize.y - 1) + 1; + if (!math::isPowerOfTwo(m_dataSize.x - 1)) m_dataSize.x = (int) math::roundToPowerOfTwo((uint32_t) m_dataSize.x - 1) + 1; + if (!math::isPowerOfTwo(m_dataSize.y - 1)) m_dataSize.y = (int) math::roundToPowerOfTwo((uint32_t) m_dataSize.y - 1) + 1; if (m_bitmap->getSize() != m_dataSize) { m_bitmap = m_bitmap->convert(Bitmap::ELuminance, Bitmap::EFloat); @@ -533,7 +542,7 @@ public: Log(EInfo, "Building acceleration data structure for %ix%i height field ..", m_dataSize.x, m_dataSize.y); ref timer = new Timer(); - m_levelCount = (int) std::max(log2i((uint32_t) m_dataSize.x-1), log2i((uint32_t) m_dataSize.y-1)) + 1; + m_levelCount = (int) std::max(math::log2i((uint32_t) m_dataSize.x-1), math::log2i((uint32_t) m_dataSize.y-1)) + 1; m_levelSize = new Vector2i[m_levelCount]; m_numChildren = new Vector2i[m_levelCount]; diff --git a/src/shapes/sphere.cpp b/src/shapes/sphere.cpp index d35ad11c..b255da6d 100644 --- a/src/shapes/sphere.cpp +++ b/src/shapes/sphere.cpp @@ -256,7 +256,7 @@ public: } void samplePosition(PositionSamplingRecord &pRec, const Point2 &sample) const { - Vector v = Warp::squareToUniformSphere(sample); + Vector v = warp::squareToUniformSphere(sample); pRec.p = Point(v * m_radius) + m_center; pRec.n = Normal(v); @@ -300,8 +300,8 @@ public: Float cosAlpha = math::safe_sqrt(1.0f - sinAlpha * sinAlpha); dRec.d = Frame(refToCenter * invRefDist).toWorld( - Warp::squareToUniformCone(cosAlpha, sample)); - dRec.pdf = Warp::squareToUniformConePdf(cosAlpha); + warp::squareToUniformCone(cosAlpha, sample)); + dRec.pdf = warp::squareToUniformConePdf(cosAlpha); /* Distance to the projection of the sphere center onto the ray (dRec.ref, dRec.d) */ @@ -336,7 +336,7 @@ public: } else { /* The reference point lies inside the sphere => use uniform sphere sampling. */ - Vector d = Warp::squareToUniformSphere(sample); + Vector d = warp::squareToUniformSphere(sample); dRec.p = m_center + d * m_radius; dRec.n = Normal(d); @@ -366,7 +366,7 @@ public: if (sinAlpha < 1-Epsilon) { /* The reference point lies outside the sphere */ Float cosAlpha = math::safe_sqrt(1 - sinAlpha*sinAlpha); - Float pdfSA = Warp::squareToUniformConePdf(cosAlpha); + Float pdfSA = warp::squareToUniformConePdf(cosAlpha); if (dRec.measure == ESolidAngle) return pdfSA; diff --git a/src/subsurface/bluenoise.cpp b/src/subsurface/bluenoise.cpp index 20ff403c..5efa405c 100644 --- a/src/subsurface/bluenoise.cpp +++ b/src/subsurface/bluenoise.cpp @@ -100,7 +100,7 @@ void blueNoisePointSet(const Scene *scene, const std::vector &shapes, sa = areaDistr.getSum(); /* Start with a fairly dense set of white noise points */ - int nsamples = ceilToInt(15 * sa / (M_PI * radius * radius)); + int nsamples = math::ceilToInt(15 * sa / (M_PI * radius * radius)); SLog(EInfo, "Creating a blue noise point set (radius=%f, " "surface area = %f)", radius, sa); @@ -145,7 +145,7 @@ void blueNoisePointSet(const Scene *scene, const std::vector &shapes, Vector3i cellCount; for (int i=0; i<3; ++i) - cellCount[i] = std::max(1, ceilToInt(extents[i] * invCellWidth)); + cellCount[i] = std::max(1, math::ceilToInt(extents[i] * invCellWidth)); SLog(EInfo, " phase 2: computing cell indices .."); #if defined(MTS_OPENMP) diff --git a/src/tests/test_chisquare.cpp b/src/tests/test_chisquare.cpp index 381853ca..ae2c580e 100644 --- a/src/tests/test_chisquare.cpp +++ b/src/tests/test_chisquare.cpp @@ -419,9 +419,9 @@ public: Vector wi; if (bsdf->getType() & BSDF::EBackSide) - wi = Warp::squareToUniformSphere(sampler->next2D()); + wi = warp::squareToUniformSphere(sampler->next2D()); else - wi = Warp::squareToCosineHemisphere(sampler->next2D()); + wi = warp::squareToCosineHemisphere(sampler->next2D()); BSDFAdapter adapter(bsdf, sampler, wi, -1); ref chiSqr = new ChiSquare(thetaBins, 2*thetaBins, wiSamples); @@ -466,9 +466,9 @@ public: Vector wi; if (bsdf->getType(comp) & BSDF::EBackSide) - wi = Warp::squareToUniformSphere(sampler->next2D()); + wi = warp::squareToUniformSphere(sampler->next2D()); else - wi = Warp::squareToCosineHemisphere(sampler->next2D()); + wi = warp::squareToCosineHemisphere(sampler->next2D()); BSDFAdapter adapter(bsdf, sampler, wi, comp); @@ -533,11 +533,11 @@ public: MediumSamplingRecord mRec; /* Sampler fiber/particle orientation */ - mRec.orientation = Warp::squareToUniformSphere(sampler->next2D()); + mRec.orientation = warp::squareToUniformSphere(sampler->next2D()); /* Test for a number of different incident directions */ for (int j=0; jnext2D()); + Vector wi = warp::squareToUniformSphere(sampler->next2D()); PhaseFunctionAdapter adapter(mRec, phase, sampler, wi); ref chiSqr = new ChiSquare(thetaBins, 2*thetaBins, wiSamples); diff --git a/src/tests/test_kd.cpp b/src/tests/test_kd.cpp index eccfb8f8..3a75a545 100644 --- a/src/tests/test_kd.cpp +++ b/src/tests/test_kd.cpp @@ -112,8 +112,8 @@ public: for (size_t i=0; inextFloat(), random->nextFloat()), sample2(random->nextFloat(), random->nextFloat()); - Point p1 = bsphere.center + Warp::squareToUniformSphere(sample1) * bsphere.radius; - Point p2 = bsphere.center + Warp::squareToUniformSphere(sample2) * bsphere.radius; + Point p1 = bsphere.center + warp::squareToUniformSphere(sample1) * bsphere.radius; + Point p2 = bsphere.center + warp::squareToUniformSphere(sample2) * bsphere.radius; Ray r(p1, normalize(p2-p1), 0.0f); Intersection its; diff --git a/src/tests/test_sh.cpp b/src/tests/test_sh.cpp index 17f471ca..883c0b67 100644 --- a/src/tests/test_sh.cpp +++ b/src/tests/test_sh.cpp @@ -41,7 +41,7 @@ public: for (int m=-l; m<=l; ++m) vec1(l, m) = random->nextFloat(); - Vector axis(Warp::squareToUniformSphere(Point2(random->nextFloat(), random->nextFloat()))); + Vector axis(warp::squareToUniformSphere(Point2(random->nextFloat(), random->nextFloat()))); Transform trafo = Transform::rotate(axis, random->nextFloat()*360); SHRotation rot(vec1.getBands()); @@ -51,7 +51,7 @@ public: rot(vec1, vec2); for (int i=0; i<100; ++i) { - Vector dir1(Warp::squareToUniformSphere(Point2(random->nextFloat(), random->nextFloat()))), dir2; + Vector dir1(warp::squareToUniformSphere(Point2(random->nextFloat(), random->nextFloat()))), dir2; trafo(dir1, dir2); Float value1 = vec1.eval(dir2); diff --git a/src/textures/checkerboard.cpp b/src/textures/checkerboard.cpp index dfe93516..776729a9 100644 --- a/src/textures/checkerboard.cpp +++ b/src/textures/checkerboard.cpp @@ -64,8 +64,8 @@ public: } inline Spectrum eval(const Point2 &uv) const { - int x = 2*modulo((int) (uv.x * 2), 2) - 1, - y = 2*modulo((int) (uv.y * 2), 2) - 1; + int x = 2*math::modulo((int) (uv.x * 2), 2) - 1, + y = 2*math::modulo((int) (uv.y * 2), 2) - 1; if (x*y == 1) return m_color0; diff --git a/src/textures/gridtexture.cpp b/src/textures/gridtexture.cpp index 557d0b91..1b700d34 100644 --- a/src/textures/gridtexture.cpp +++ b/src/textures/gridtexture.cpp @@ -73,8 +73,8 @@ public: } inline Spectrum eval(const Point2 &uv) const { - Float x = uv.x - floorToInt(uv.x); - Float y = uv.y - floorToInt(uv.y); + Float x = uv.x - math::floorToInt(uv.x); + Float y = uv.y - math::floorToInt(uv.y); if (x > .5) x -= 1; diff --git a/src/textures/wireframe.cpp b/src/textures/wireframe.cpp index c6867b89..c4a756b4 100644 --- a/src/textures/wireframe.cpp +++ b/src/textures/wireframe.cpp @@ -118,7 +118,7 @@ public: minDist = std::min(minDist, (cur + d1 * dot(d1, d2) - its.p).lengthSquared()); } - Float a = smoothStep(m_lineWidth*(1.f-m_stepWidth), m_lineWidth, std::sqrt(minDist)); + Float a = math::smoothStep(m_lineWidth*(1.f-m_stepWidth), m_lineWidth, std::sqrt(minDist)); return m_edgeColor*(1-a) + m_interiorColor*a; } diff --git a/src/utils/kdbench.cpp b/src/utils/kdbench.cpp index 73f6653e..ad5b3c39 100644 --- a/src/utils/kdbench.cpp +++ b/src/utils/kdbench.cpp @@ -223,8 +223,8 @@ public: for (size_t i=0; inextFloat(), random->nextFloat()), sample2(random->nextFloat(), random->nextFloat()); - Point p1 = bsphere.center + Warp::squareToUniformSphere(sample1) * bsphere.radius; - Point p2 = bsphere.center + Warp::squareToUniformSphere(sample2) * bsphere.radius; + Point p1 = bsphere.center + warp::squareToUniformSphere(sample1) * bsphere.radius; + Point p2 = bsphere.center + warp::squareToUniformSphere(sample2) * bsphere.radius; Ray r(p1, normalize(p2-p1), 0.0f); Intersection its; diff --git a/src/volume/gridvolume.cpp b/src/volume/gridvolume.cpp index 23d2ccf6..8aefdded 100644 --- a/src/volume/gridvolume.cpp +++ b/src/volume/gridvolume.cpp @@ -336,9 +336,9 @@ public: Float lookupFloat(const Point &_p) const { const Point p = m_worldToGrid.transformAffine(_p); - const int x1 = floorToInt(p.x), - y1 = floorToInt(p.y), - z1 = floorToInt(p.z), + const int x1 = math::floorToInt(p.x), + y1 = math::floorToInt(p.y), + z1 = math::floorToInt(p.z), x2 = x1+1, y2 = y1+1, z2 = z1+1; if (x1 < 0 || y1 < 0 || z1 < 0 || x2 >= m_res.x || @@ -389,9 +389,9 @@ public: Spectrum lookupSpectrum(const Point &_p) const { const Point p = m_worldToGrid.transformAffine(_p); - const int x1 = floorToInt(p.x), - y1 = floorToInt(p.y), - z1 = floorToInt(p.z), + const int x1 = math::floorToInt(p.x), + y1 = math::floorToInt(p.y), + z1 = math::floorToInt(p.z), x2 = x1+1, y2 = y1+1, z2 = z1+1; if (x1 < 0 || y1 < 0 || z1 < 0 || x2 >= m_res.x || @@ -466,9 +466,9 @@ public: Vector lookupVector(const Point &_p) const { const Point p = m_worldToGrid.transformAffine(_p); - const int x1 = floorToInt(p.x), - y1 = floorToInt(p.y), - z1 = floorToInt(p.z), + const int x1 = math::floorToInt(p.x), + y1 = math::floorToInt(p.y), + z1 = math::floorToInt(p.z), x2 = x1+1, y2 = y1+1, z2 = z1+1; if (x1 < 0 || y1 < 0 || z1 < 0 || x2 >= m_res.x || diff --git a/src/volume/hgridvolume.cpp b/src/volume/hgridvolume.cpp index b449c0c0..5e8d8f1c 100644 --- a/src/volume/hgridvolume.cpp +++ b/src/volume/hgridvolume.cpp @@ -143,9 +143,9 @@ public: Float lookupFloat(const Point &_p) const { const Point p = m_worldToGrid.transformAffine(_p); - const int x = floorToInt(p.x), - y = floorToInt(p.y), - z = floorToInt(p.z); + const int x = math::floorToInt(p.x), + y = math::floorToInt(p.y), + z = math::floorToInt(p.z); if (x < 0 || x >= m_res.x || y < 0 || y >= m_res.y || z < 0 || z >= m_res.z) @@ -160,9 +160,9 @@ public: Spectrum lookupSpectrum(const Point &_p) const { const Point p = m_worldToGrid.transformAffine(_p); - const int x = floorToInt(p.x), - y = floorToInt(p.y), - z = floorToInt(p.z); + const int x = math::floorToInt(p.x), + y = math::floorToInt(p.y), + z = math::floorToInt(p.z); if (x < 0 || x >= m_res.x || y < 0 || y >= m_res.y || z < 0 || z >= m_res.z) @@ -177,9 +177,9 @@ public: Vector lookupVector(const Point &_p) const { const Point p = m_worldToGrid.transformAffine(_p); - const int x = floorToInt(p.x), - y = floorToInt(p.y), - z = floorToInt(p.z); + const int x = math::floorToInt(p.x), + y = math::floorToInt(p.y), + z = math::floorToInt(p.z); if (x < 0 || x >= m_res.x || y < 0 || y >= m_res.y || z < 0 || z >= m_res.z) diff --git a/src/volume/volcache.cpp b/src/volume/volcache.cpp index 31d896d8..7cf9b7fe 100644 --- a/src/volume/volcache.cpp +++ b/src/volume/volcache.cpp @@ -86,7 +86,7 @@ public: /// Size of an individual block (must be a power of 2) m_blockSize = props.getInteger("blockSize", 8); - if (!isPowerOfTwo(m_blockSize)) + if (!math::isPowerOfTwo(m_blockSize)) Log(EError, "Block size must be a power of two!"); /* Width of an individual voxel. Will use the step size of the @@ -150,7 +150,7 @@ public: * Transform::translate(-Vector(m_aabb.min)) * m_worldToVolume; m_voxelMask = m_blockSize-1; m_blockMask = ~(m_blockSize-1); - m_blockShift = log2i((uint32_t) m_blockSize); + m_blockShift = math::log2i((uint32_t) m_blockSize); Log(EInfo, "Volume cache configuration"); Log(EInfo, " Block size in voxels = %i", m_blockSize);