From eb4a823ead0a588bd99a2d11ae1da314fbbfa4f2 Mon Sep 17 00:00:00 2001 From: Wenzel Jakob Date: Tue, 27 Nov 2012 11:47:25 -0500 Subject: [PATCH] fully switched over to new spline code --- include/mitsuba/core/util.h | 174 ----------------- src/bsdfs/rtrans.h | 51 ++--- src/libcore/util.cpp | 371 ------------------------------------ src/tests/test_rtrans.cpp | 2 +- 4 files changed, 28 insertions(+), 570 deletions(-) diff --git a/include/mitsuba/core/util.h b/include/mitsuba/core/util.h index c31e6df7..2a02fd2e 100644 --- a/include/mitsuba/core/util.h +++ b/include/mitsuba/core/util.h @@ -340,180 +340,6 @@ extern MTS_EXPORT_CORE bool solveQuadratic(Float a, Float b, extern MTS_EXPORT_CORE bool solveQuadraticDouble(double a, double b, double c, double &x0, double &x1); -/** - * \brief Evaluate a cubic spline interpolant of a regularly sampled 1D function - * - * This implementation relies on Catmull-Rom splines, i.e. it uses finite - * differences to approximate the derivatives at the endpoints of each spline - * segment. - * - * \param x - * Evaluation point - * \param data - * Floating point array containing \c size regularly spaced evaluations - * in the range [\c min,\c max] of the function to be approximated. - * \param min - * Position of the first knot - * \param max - * Position of the last knot - * \param size - * Denotes the size of the \c data array - * \param extrapolate - * Extrapolate data values when \c x is out of range? (default: \c false) - * \return - * The interpolated value or zero when extrapolate=falsett> - * and \c x lies outside of [\c min, \c max] - */ -extern MTS_EXPORT_CORE Float interpCubic1D(Float x, const Float *data, - Float min, Float max, size_t size, bool extrapolate = false); - -/** - * \brief Evaluate a cubic spline interpolant of an \a irregularly sampled 1D function - * - * This implementation relies on Catmull-Rom splines, i.e. it uses finite - * differences to approximate the derivatives at the endpoints of each spline - * segment. - * - * \param x - * Evaluation point - * \param nodes - * Floating point array containing \c size irregularly spaced values - * denoting positions the where the function to be interpolated was evaluated. - * They must be provided in \a increasing order. - * \param data - * Floating point array containing interpolant values matched to - * the entries of \c nodes. - * \param size - * Denotes the size of the \c data array - * \param extrapolate - * Extrapolate data values when \c x is out of range? (default: \c false) - * \return - * The interpolated value or zero when extrapolate=falsett> - * and \c x lies outside of \a [\c min, \c max] - */ -extern MTS_EXPORT Float interpCubic1DIrregular(Float x, const Float *nodes, - const Float *data, size_t size, bool extrapolate = false); - -/** - * \brief Evaluate a cubic spline interpolant of a regularly sampled 2D function - * - * This implementation relies on a tensor product of Catmull-Rom splines, i.e. it uses - * finite differences to approximate the derivatives at the endpoints of each spline - * patch. - * - * \param p - * Evaluation point - * \param data - * A 2D floating point array of size.x*size.y cells containing regularly - * spaced evaluations of the function to be interpolated on the domain [min, max]. - * Consecutive entries of this array correspond to increments in the 'x' coordinate. - * \param min - * Position of the first knot on each dimension - * \param max - * Position of the last knot on each dimension - * \param size - * Denotes the size of the \c data array (along each dimension) - * \param extrapolate - * Extrapolate data values when \c p is out of range? (default: \c false) - * \return - * The interpolated value or zero when extrapolate=falsett> and - * \c p lies outside of the knot range - */ -extern MTS_EXPORT_CORE Float interpCubic2D(const Point2 &p, const Float *data, - const Point2 &min, const Point2 &max, const Size2 &size, bool extrapolate = false); - -/** - * \brief Evaluate a cubic spline interpolant of an \a irregularly sampled 2D function - * - * This implementation relies on a tensor product of Catmull-Rom splines, i.e. it uses - * finite differences to approximate the derivatives at the endpoints of each spline - * region. - * - * When the underlying function is sampled on a regular grid, \ref interpCubic2D() - * should be preferred, since data lookups will be considerably faster. - * - * \param p - * Evaluation point - * \param nodes - * Pointer to a list for each dimension denoting the positions where the function - * to be interpolated was evaluated. The i-th array must have - * size size[i] and contain position values in \a increasing order. - * \param data - * A 2D floating point array of size.x*size.y cells containing irregularly - * spaced evaluations of the function to be interpolated on the domain [min, max]. - * Consecutive entries of this array correspond to increments in the 'x' coordinate. - * \param size - * Denotes the size of the \c data array (along each dimension) - * \param extrapolate - * Extrapolate data values when \c p is out of range? (default: \c false) - * \return - * The interpolated value or zero when extrapolate=falsett> and - * \c p lies outside of the knot range - */ -extern MTS_EXPORT_CORE Float interpCubic2DIrregular(const Point2 &p, const Float **nodes, - const Float *data, const Size2 &size, bool extrapolate = false); - -/** - * \brief Evaluate a cubic spline interpolant of a regularly sampled 3D function - * - * This implementation relies on a tensor product of Catmull-Rom splines, i.e. it uses - * finite differences to approximate the derivatives at the endpoints of each spline - * region. - * - * \param p - * Evaluation point of the interpolant - * \param data - * A 3D floating point array of size.x*size.y*size.z cells containing regularly - * spaced evaluations of the function to be interpolated on the domain [min, max]. - * Consecutive entries of this array correspond to increments in the 'x' coordinate, - * then 'y', and finally 'z' increments. - * \param min - * Position of the first knot on each dimension - * \param max - * Position of the last knot on each dimension - * \param size - * Denotes the size of the \c data array (along each dimension) - * \param extrapolate - * Extrapolate data values when \c p is out of range? (default: \c false) - * \return - * The interpolated value or zero when extrapolate=falsett> and - * \c p lies outside of the knot range - */ -extern MTS_EXPORT_CORE Float interpCubic3D(const Point3 &p, const Float *data, - const Point3 &min, const Point3 &max, const Size3 &size, bool extrapolate = false); - -/** - * \brief Evaluate a cubic spline interpolant of an \a irregularly sampled 3D function - * - * This implementation relies on a tensor product of Catmull-Rom splines, i.e. it uses - * finite differences to approximate the derivatives at the endpoints of each spline - * region. - * - * When the underlying function is sampled on a regular grid, \ref interpCubic3D() - * should be preferred, since data lookups will be considerably faster. - * - * \param p - * Evaluation point - * \param nodes - * Pointer to a list for each dimension denoting the positions where the function - * to be interpolated was evaluated. The i-th array must have - * size size[i] and contain position values in \a increasing order. - * \param data - * A 2D floating point array of size.x*size.y cells containing irregularly - * spaced evaluations of the function to be interpolated on the domain [min, max]. - * Consecutive entries of this array correspond to increments in the 'x' coordinate, - * then 'y', and finally 'z' increments. - * \param size - * Denotes the size of the \c data array (along each dimension) - * \param extrapolate - * Extrapolate data values when \c p is out of range? (default: \c false) - * \return - * The interpolated value or zero when extrapolate=falsett> and - * \c p lies outside of the knot range - */ -extern MTS_EXPORT_CORE Float interpCubic3DIrregular(const Point3 &p, const Float **nodes, - const Float *data, const Size3 &size, bool extrapolate = false); - //// Convert radians to degrees inline Float radToDeg(Float value) { return value * (180.0f / M_PI); } diff --git a/src/bsdfs/rtrans.h b/src/bsdfs/rtrans.h index 14fbb4cb..1c4c248f 100644 --- a/src/bsdfs/rtrans.h +++ b/src/bsdfs/rtrans.h @@ -20,6 +20,7 @@ #define __ROUGH_TRANSMITTANCE_H #include +#include #include #include "microfacet.h" @@ -186,17 +187,17 @@ public: if (m_alphaFixed && m_etaFixed) { SAssert(cosTheta >= 0); - result = interpCubic1D(warpedCosTheta, - m_trans, 0.0f, 1.0f, m_thetaSamples); + result = evalCubicInterp1D(warpedCosTheta, + m_trans, m_thetaSamples, 0.0f, 1.0f); } else if (m_etaFixed) { SAssert(cosTheta >= 0); Float warpedAlpha = std::pow((alpha - m_alphaMin) / (m_alphaMax-m_alphaMin), (Float) 0.25f); - result = interpCubic2D(Point2(warpedCosTheta, warpedAlpha), - m_trans, Point2(0.0f), Point2(1.0f), - Size2(m_thetaSamples, m_alphaSamples)); + result = evalCubicInterp2D(Point2(warpedCosTheta, warpedAlpha), + m_trans, Size2(m_thetaSamples, m_alphaSamples), + Point2(0.0f), Point2(1.0f)); } else { if (cosTheta < 0) { cosTheta = -cosTheta; @@ -220,9 +221,9 @@ public: Float warpedEta = std::pow((eta - m_etaMin) / (m_etaMax-m_etaMin), (Float) 0.25f); - result = interpCubic3D(Point3(warpedCosTheta, warpedAlpha, warpedEta), - data, Point3(0.0f), Point3(1.0f), - Size3(m_thetaSamples, m_alphaSamples, m_etaSamples)); + result = evalCubicInterp3D(Point3(warpedCosTheta, warpedAlpha, warpedEta), + data, Size3(m_thetaSamples, m_alphaSamples, m_etaSamples), + Point3(0.0f), Point3(1.0f)); } return std::min((Float) 1.0f, std::max((Float) 0.0f, result)); @@ -250,8 +251,8 @@ public: Float warpedAlpha = std::pow((alpha - m_alphaMin) / (m_alphaMax-m_alphaMin), (Float) 0.25f); - result = interpCubic1D(warpedAlpha, m_diffTrans, - 0.0f, 1.0f, m_alphaSamples); + result = evalCubicInterp1D(warpedAlpha, + m_diffTrans, m_alphaSamples, 0.0f, 1.0f); } else { Float *data = m_diffTrans; if (eta < 1) { @@ -270,8 +271,9 @@ public: Float warpedEta = std::pow((eta - m_etaMin) / (m_etaMax-m_etaMin), (Float) 0.25f); - result = interpCubic2D(Point2(warpedAlpha, warpedEta), data, - Point2(0.0f), Point2(1.0f), Size2(m_alphaSamples, m_etaSamples)); + result = evalCubicInterp2D(Point2(warpedAlpha, warpedEta), data, + Size2(m_alphaSamples, m_etaSamples), Point2(0.0f), Point2(1.0f)); + } return std::min((Float) 1.0f, std::max((Float) 0.0f, result)); @@ -317,16 +319,17 @@ public: dTheta = 1.0f / (m_thetaSamples - 1); for (size_t i=0; i= min && x <= max) && !extrapolate) - return 0.0f; - - /* Transform 'x' so that knots lie at integer positions */ - Float t = ((x - min) * (size - 1)) / (max - min); - - /* Find the index of the left knot in the queried subinterval, be - robust to cases where 't' lies exactly on the right endpoint */ - size_t k = std::max((size_t) 0, std::min((size_t) t, size - 2)); - - Float f0 = data[k], - f1 = data[k+1], - d0, d1; - - /* Approximate the derivatives */ - if (k > 0) - d0 = 0.5f * (data[k+1] - data[k-1]); - else - d0 = data[k+1] - data[k]; - - if (k + 2 < size) - d1 = 0.5f * (data[k+2] - data[k]); - else - d1 = data[k+1] - data[k]; - - /* Compute the relative position within the interval */ - t = t - (Float) k; - - Float t2 = t*t, t3 = t2*t; - - return - ( 2*t3 - 3*t2 + 1) * f0 + - (-2*t3 + 3*t2) * f1 + - ( t3 - 2*t2 + t) * d0 + - ( t3 - t2) * d1; -} - -Float interpCubic1DIrregular(Float x, const Float *nodes, const Float *data, size_t size, bool extrapolate) { - /* Give up when given an out-of-range or NaN argument */ - if (!(x >= nodes[0] && x <= nodes[size-1]) && !extrapolate) - return 0.0f; - - size_t k = (size_t) std::max((ptrdiff_t) 0, std::min((ptrdiff_t) size - 2, - std::lower_bound(nodes, nodes + size, x) - nodes - 1)); - - Float f0 = data[k], - f1 = data[k+1], - width = nodes[k+1] - nodes[k], - invWidth = 1.0f / width, - d0, d1; - - /* Approximate the derivatives */ - if (k > 0) - d0 = (f1 - data[k-1]) / (nodes[k+1] - nodes[k-1]); - else - d0 = (f1 - f0) * invWidth; - - if (k + 2 < size) - d1 = (data[k+2] - f0) / (nodes[k+2] - nodes[k]); - else - d1 = (f1 - f0) * invWidth; - - Float t = (x - nodes[k]) * invWidth; - Float t2 = t*t, t3 = t2*t; - - return - ( 2*t3 - 3*t2 + 1) * f0 + - (-2*t3 + 3*t2) * f1 + - (( t3 - 2*t2 + t) * d0 + - ( t3 - t2) * d1) * width; -} - - -Float interpCubic2D(const Point2 &p, const Float *data, - const Point2 &min, const Point2 &max, const Size2 &size, bool extrapolate) { - Float knotWeights[2][4]; - Size2 knot; - - /* Compute interpolation weights separately for each dimension */ - for (int dim=0; dim<2; ++dim) { - Float *weights = knotWeights[dim]; - /* Give up when given an out-of-range or NaN argument */ - if (!(p[dim] >= min[dim] && p[dim] <= max[dim]) && !extrapolate) - return 0.0f; - - /* Transform 'p' so that knots lie at integer positions */ - Float t = ((p[dim] - min[dim]) * (size[dim] - 1)) - / (max[dim]-min[dim]); - - /* Find the index of the left knot in the queried subinterval, be - robust to cases where 't' lies exactly on the right endpoint */ - knot[dim] = std::min((size_t) t, size[dim] - 2); - - /* Compute the relative position within the interval */ - t = t - (Float) knot[dim]; - - /* Compute node weights */ - Float t2 = t*t, t3 = t2*t; - weights[0] = 0.0f; - weights[1] = 2*t3 - 3*t2 + 1; - weights[2] = -2*t3 + 3*t2; - weights[3] = 0.0f; - - /* Derivative weights */ - Float d0 = t3 - 2*t2 + t, - d1 = t3 - t2; - - /* Turn derivative weights into node weights using - an appropriate chosen finite differences stencil */ - if (knot[dim] > 0) { - weights[2] += 0.5f * d0; - weights[0] -= 0.5f * d0; - } else { - weights[2] += d0; - weights[1] -= d0; - } - - if (knot[dim] + 2 < size[dim]) { - weights[3] += 0.5f * d1; - weights[1] -= 0.5f * d1; - } else { - weights[2] += d1; - weights[1] -= d1; - } - } - - Float result = 0.0f; - for (int y=-1; y<=2; ++y) { - Float wy = knotWeights[1][y+1]; - for (int x=-1; x<=2; ++x) { - Float wxy = knotWeights[0][x+1] * wy; - - if (wxy == 0) - continue; - - size_t pos = (knot[1] + y) * size[0] + knot[0] + x; - - result += data[pos] * wxy; - } - } - return result; -} - -Float interpCubic2DIrregular(const Point2 &p, const Float **nodes_, - const Float *data, const Size2 &size, bool extrapolate) { - Float knotWeights[2][4]; - Size2 knot; - - /* Compute interpolation weights separately for each dimension */ - for (int dim=0; dim<2; ++dim) { - const Float *nodes = nodes_[dim]; - Float *weights = knotWeights[dim]; - - /* Give up when given an out-of-range or NaN argument */ - if (!(p[dim] >= nodes[0] && p[dim] <= nodes[size[dim]-1]) && !extrapolate) - return 0.0f; - - /* Find the index of the left knot in the queried subinterval, be - robust to cases where 't' lies exactly on the right endpoint */ - size_t k = (size_t) std::max((ptrdiff_t) 0, std::min((ptrdiff_t) size[dim] - 2, - std::lower_bound(nodes, nodes + size[dim], p[dim]) - nodes - 1)); - knot[dim] = k; - - Float width = nodes[k+1] - nodes[k], invWidth = 1 / width; - - /* Compute the relative position within the interval */ - Float t = (p[dim] - nodes[k]) * invWidth, - t2 = t*t, t3 = t2*t; - - /* Compute node weights */ - weights[0] = 0.0f; - weights[1] = 2*t3 - 3*t2 + 1; - weights[2] = -2*t3 + 3*t2; - weights[3] = 0.0f; - - /* Derivative weights */ - Float d0 = (t3 - 2*t2 + t) * width, - d1 = (t3 - t2) * width; - - /* Turn derivative weights into node weights using - an appropriate chosen finite differences stencil */ - if (k > 0) { - Float factor = 1 / (nodes[k+1]-nodes[k-1]); - weights[2] += d0 * factor; - weights[0] -= d0 * factor; - } else { - weights[2] += d0 * invWidth; - weights[1] -= d0 * invWidth; - } - - if (k + 2 < size[dim]) { - Float factor = 1 / (nodes[k+2]-nodes[k]); - weights[3] += d1 * factor; - weights[1] -= d1 * factor; - } else { - weights[2] += d1 * invWidth; - weights[1] -= d1 * invWidth; - } - } - - Float result = 0.0f; - for (int y=-1; y<=2; ++y) { - Float wy = knotWeights[1][y+1]; - for (int x=-1; x<=2; ++x) { - Float wxy = knotWeights[0][x+1] * wy; - - if (wxy == 0) - continue; - - size_t pos = (knot[1] + y) * size[0] + knot[0] + x; - - result += data[pos] * wxy; - } - } - return result; -} - -Float interpCubic3D(const Point3 &p, const Float *data, - const Point3 &min, const Point3 &max, const Size3 &size, bool extrapolate) { - Float knotWeights[3][4]; - Size3 knot; - - /* Compute interpolation weights separately for each dimension */ - for (int dim=0; dim<3; ++dim) { - Float *weights = knotWeights[dim]; - /* Give up when given an out-of-range or NaN argument */ - if (!(p[dim] >= min[dim] && p[dim] <= max[dim]) && !extrapolate) - return 0.0f; - - /* Transform 'p' so that knots lie at integer positions */ - Float t = ((p[dim] - min[dim]) * (size[dim] - 1)) - / (max[dim]-min[dim]); - - /* Find the index of the left knot in the queried subinterval, be - robust to cases where 't' lies exactly on the right endpoint */ - knot[dim] = std::min((size_t) t, size[dim] - 2); - - /* Compute the relative position within the interval */ - t = t - (Float) knot[dim]; - - /* Compute node weights */ - Float t2 = t*t, t3 = t2*t; - weights[0] = 0.0f; - weights[1] = 2*t3 - 3*t2 + 1; - weights[2] = -2*t3 + 3*t2; - weights[3] = 0.0f; - - /* Derivative weights */ - Float d0 = t3 - 2*t2 + t, - d1 = t3 - t2; - - /* Turn derivative weights into node weights using - an appropriate chosen finite differences stencil */ - if (knot[dim] > 0) { - weights[2] += 0.5f * d0; - weights[0] -= 0.5f * d0; - } else { - weights[2] += d0; - weights[1] -= d0; - } - - if (knot[dim] + 2 < size[dim]) { - weights[3] += 0.5f * d1; - weights[1] -= 0.5f * d1; - } else { - weights[2] += d1; - weights[1] -= d1; - } - } - - Float result = 0.0f; - for (int z=-1; z<=2; ++z) { - Float wz = knotWeights[2][z+1]; - for (int y=-1; y<=2; ++y) { - Float wyz = knotWeights[1][y+1] * wz; - for (int x=-1; x<=2; ++x) { - Float wxyz = knotWeights[0][x+1] * wyz; - - if (wxyz == 0) - continue; - - size_t pos = ((knot[2] + z) * size[1] + (knot[1] + y)) - * size[0] + knot[0] + x; - - result += data[pos] * wxyz; - } - } - } - return result; -} - -Float interpCubic3DIrregular(const Point3 &p, const Float **nodes_, - const Float *data, const Size3 &size, bool extrapolate) { - Float knotWeights[3][4]; - Size3 knot; - - /* Compute interpolation weights separately for each dimension */ - for (int dim=0; dim<3; ++dim) { - const Float *nodes = nodes_[dim]; - Float *weights = knotWeights[dim]; - - /* Give up when given an out-of-range or NaN argument */ - if (!(p[dim] >= nodes[0] && p[dim] <= nodes[size[dim]-1]) && !extrapolate) - return 0.0f; - - /* Find the index of the left knot in the queried subinterval, be - robust to cases where 't' lies exactly on the right endpoint */ - size_t k = (size_t) std::max((ptrdiff_t) 0, std::min((ptrdiff_t) size[dim] - 2, - std::lower_bound(nodes, nodes + size[dim], p[dim]) - nodes - 1)); - knot[dim] = k; - - Float width = nodes[k+1] - nodes[k], invWidth = 1 / width; - - /* Compute the relative position within the interval */ - Float t = (p[dim] - nodes[k]) * invWidth, - t2 = t*t, t3 = t2*t; - - /* Compute node weights */ - weights[0] = 0.0f; - weights[1] = 2*t3 - 3*t2 + 1; - weights[2] = -2*t3 + 3*t2; - weights[3] = 0.0f; - - /* Derivative weights */ - Float d0 = (t3 - 2*t2 + t) * width, - d1 = (t3 - t2) * width; - - /* Turn derivative weights into node weights using - an appropriate chosen finite differences stencil */ - if (k > 0) { - Float factor = 1 / (nodes[k+1]-nodes[k-1]); - weights[2] += d0 * factor; - weights[0] -= d0 * factor; - } else { - weights[2] += d0 * invWidth; - weights[1] -= d0 * invWidth; - } - - if (k + 2 < size[dim]) { - Float factor = 1 / (nodes[k+2]-nodes[k]); - weights[3] += d1 * factor; - weights[1] -= d1 * factor; - } else { - weights[2] += d1 * invWidth; - weights[1] -= d1 * invWidth; - } - } - - Float result = 0.0f; - for (int z=-1; z<=2; ++z) { - Float wz = knotWeights[2][z+1]; - for (int y=-1; y<=2; ++y) { - Float wyz = knotWeights[1][y+1] * wz; - for (int x=-1; x<=2; ++x) { - Float wxyz = knotWeights[0][x+1] * wyz; - - if (wxyz == 0) - continue; - - size_t pos = ((knot[2] + z) * size[1] + (knot[1] + y)) - * size[0] + knot[0] + x; - - result += data[pos] * wxyz; - } - } - } - return result; -} - void stratifiedSample1D(Random *random, Float *dest, int count, bool jitter) { Float invCount = 1.0f / count; diff --git a/src/tests/test_rtrans.cpp b/src/tests/test_rtrans.cpp index 49c98637..cd9332a4 100644 --- a/src/tests/test_rtrans.cpp +++ b/src/tests/test_rtrans.cpp @@ -37,7 +37,7 @@ void transmittanceIntegrand(const BSDF *bsdf, const Vector &wi, size_t nPts, con void diffTransmittanceIntegrand(Float *data, size_t resolution, size_t nPts, const Float *in, Float *out) { for (size_t i=0; i