fully switched over to new spline code
parent
9708dbb3d6
commit
eb4a823ead
|
@ -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 <tt>extrapolate=false</tt>tt>
|
||||
* 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 <tt>extrapolate=false</tt>tt>
|
||||
* 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 <tt>size.x*size.y</tt> cells containing regularly
|
||||
* spaced evaluations of the function to be interpolated on the domain <tt>[min, max]</tt>.
|
||||
* 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 <tt>extrapolate=false</tt>tt> 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 <tt>i</tt>-th array must have
|
||||
* size <tt>size[i]</tt> and contain position values in \a increasing order.
|
||||
* \param data
|
||||
* A 2D floating point array of <tt>size.x*size.y</tt> cells containing irregularly
|
||||
* spaced evaluations of the function to be interpolated on the domain <tt>[min, max]</tt>.
|
||||
* 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 <tt>extrapolate=false</tt>tt> 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 <tt>size.x*size.y*size.z</tt> cells containing regularly
|
||||
* spaced evaluations of the function to be interpolated on the domain <tt>[min, max]</tt>.
|
||||
* 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 <tt>extrapolate=false</tt>tt> 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 <tt>i</tt>-th array must have
|
||||
* size <tt>size[i]</tt> and contain position values in \a increasing order.
|
||||
* \param data
|
||||
* A 2D floating point array of <tt>size.x*size.y</tt> cells containing irregularly
|
||||
* spaced evaluations of the function to be interpolated on the domain <tt>[min, max]</tt>.
|
||||
* 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 <tt>extrapolate=false</tt>tt> 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); }
|
||||
|
||||
|
|
|
@ -20,6 +20,7 @@
|
|||
#define __ROUGH_TRANSMITTANCE_H
|
||||
|
||||
#include <mitsuba/core/fstream.h>
|
||||
#include <mitsuba/core/spline.h>
|
||||
#include <mitsuba/core/fresolver.h>
|
||||
#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<m_alphaSamples; ++i) {
|
||||
for (size_t j=0; j<m_thetaSamples; ++j)
|
||||
newTrans[i*m_thetaSamples + j] = interpCubic3D(
|
||||
for (size_t j=0; j<m_thetaSamples; ++j) {
|
||||
newTrans[i*m_thetaSamples + j] = evalCubicInterp3D(
|
||||
Point3(j*dTheta, i*dAlpha, warpedEta),
|
||||
trans, Point3(0.0f), Point3(1.0f),
|
||||
Size3(m_thetaSamples, m_alphaSamples, m_etaSamples));
|
||||
trans, Size3(m_thetaSamples, m_alphaSamples, m_etaSamples),
|
||||
Point3(0.0f), Point3(1.0f));
|
||||
}
|
||||
|
||||
newDiffTrans[i] = interpCubic2D(
|
||||
newDiffTrans[i] = evalCubicInterp2D(
|
||||
Point2(i*dAlpha, warpedEta),
|
||||
diffTrans, Point2(0.0f), Point2(1.0f),
|
||||
Size2(m_alphaSamples, m_etaSamples));
|
||||
diffTrans, Size2(m_alphaSamples, m_etaSamples),
|
||||
Point2(0.0f), Point2(1.0f));
|
||||
}
|
||||
|
||||
delete[] m_trans;
|
||||
|
@ -364,13 +367,13 @@ public:
|
|||
Float dTheta = 1.0f / (m_thetaSamples - 1);
|
||||
|
||||
for (size_t i=0; i<m_thetaSamples; ++i)
|
||||
newTrans[i] = interpCubic2D(
|
||||
newTrans[i] = evalCubicInterp2D(
|
||||
Point2(i*dTheta, warpedAlpha),
|
||||
m_trans, Point2(0.0f), Point2(1.0f),
|
||||
Size2(m_thetaSamples, m_alphaSamples));
|
||||
m_trans, Size2(m_thetaSamples, m_alphaSamples),
|
||||
Point2(0.0f), Point2(1.0f));
|
||||
|
||||
newDiffTrans[0] = interpCubic1D(warpedAlpha,
|
||||
m_diffTrans, 0.0f, 1.0f, m_alphaSamples);
|
||||
newDiffTrans[0] = evalCubicInterp1D(warpedAlpha,
|
||||
m_diffTrans, m_alphaSamples, 0.0f, 1.0f);
|
||||
|
||||
delete[] m_trans;
|
||||
delete[] m_diffTrans;
|
||||
|
|
|
@ -469,377 +469,6 @@ bool solveLinearSystem2x2(const Float a[2][2], const Float b[2], Float x[2]) {
|
|||
return true;
|
||||
}
|
||||
|
||||
Float interpCubic1D(Float x, const Float *data, Float min, Float max, size_t size, bool extrapolate) {
|
||||
/* Give up when given an out-of-range or NaN argument */
|
||||
if (!(x >= 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;
|
||||
|
||||
|
|
|
@ -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<nPts; ++i)
|
||||
out[i] = 2 * in[i] * interpCubic1D(in[i], data, 0, 1, resolution);
|
||||
out[i] = 2 * in[i] * evalCubicInterp1D(in[i], data, resolution, 0, 1);
|
||||
}
|
||||
|
||||
class TestRoughTransmittance : public TestCase {
|
||||
|
|
Loading…
Reference in New Issue