提交 d5ee6219 编写于 作者: M Matt Pharr

Update from book source. No functional changes.

上级 db69a235
......@@ -215,6 +215,7 @@ PixelSensor *PixelSensor::Create(const ParameterDictionary &parameters,
if (sensorName != "cie1931" && whiteBalanceTemp == 0)
whiteBalanceTemp = 6500;
const Float K_m = 683;
Float imagingRatio = Pi * exposureTime * ISO * K_m / (C * fNumber * fNumber);
if (sensorName == "cie1931") {
......
......@@ -177,6 +177,23 @@ inline RGB Lerp(Float t, const RGB &s1, const RGB &s2) {
class XYZ {
public:
// XYZ Public Methods
XYZ() = default;
PBRT_CPU_GPU
XYZ(Float X, Float Y, Float Z) : X(X), Y(Y), Z(Z) {}
PBRT_CPU_GPU
Float Average() const { return (X + Y + Z) / 3; }
PBRT_CPU_GPU
Point2f xy() const { return Point2f(X / (X + Y + Z), Y / (X + Y + Z)); }
PBRT_CPU_GPU
static XYZ FromxyY(Point2f xy, Float Y = 1) {
if (xy.y == 0)
return XYZ(0, 0, 0);
return XYZ(xy.x * Y / xy.y, Y, (1 - xy.x - xy.y) * Y / xy.y);
}
PBRT_CPU_GPU
XYZ &operator+=(const XYZ &s) {
X += s.X;
......@@ -286,23 +303,6 @@ class XYZ {
std::string ToString() const;
XYZ() = default;
PBRT_CPU_GPU
XYZ(Float X, Float Y, Float Z) : X(X), Y(Y), Z(Z) {}
PBRT_CPU_GPU
Float Average() const { return (X + Y + Z) / 3; }
PBRT_CPU_GPU
Point2f xy() const { return Point2f(X / (X + Y + Z), Y / (X + Y + Z)); }
PBRT_CPU_GPU
static XYZ FromxyY(Point2f xy, Float Y = 1) {
if (xy.y == 0)
return XYZ(0, 0, 0);
return XYZ(xy.x * Y / xy.y, Y, (1 - xy.x - xy.y) * Y / xy.y);
}
// XYZ Public Members
Float X = 0, Y = 0, Z = 0;
};
......
......@@ -105,10 +105,6 @@ template <int N>
class SquareMatrix;
// Math Inline Functions
PBRT_CPU_GPU inline Float Lerp(Float x, Float a, Float b) {
return (1 - x) * a + x * b;
}
// http://www.plunk.org/~hatch/rightway.php
PBRT_CPU_GPU
inline Float SinXOverX(Float x) {
......@@ -134,6 +130,10 @@ inline Float WindowedSinc(Float x, Float radius, Float tau) {
return Sinc(x) * Sinc(x / tau);
}
PBRT_CPU_GPU inline Float Lerp(Float x, Float a, Float b) {
return (1 - x) * a + x * b;
}
#ifdef PBRT_IS_MSVC
#pragma warning(push)
#pragma warning(disable : 4018) // signed/unsigned mismatch
......
......@@ -25,15 +25,6 @@
namespace pbrt {
// Sampling Function Definitions
Point2f RejectionSampleDisk(RNG &rng) {
Point2f p;
do {
p.x = 1 - 2 * rng.Uniform<Float>();
p.y = 1 - 2 * rng.Uniform<Float>();
} while (p.x * p.x + p.y * p.y > 1);
return p;
}
pstd::array<Float, 3> SampleSphericalTriangle(const pstd::array<Point3f, 3> &v,
const Point3f &p, const Point2f &u,
Float *pdf) {
......@@ -390,6 +381,15 @@ Vector3f SampleHenyeyGreenstein(Vector3f wo, Float g, Point2f u, Float *pdf) {
return wi;
}
Point2f RejectionSampleDisk(RNG &rng) {
Point2f p;
do {
p.x = 1 - 2 * rng.Uniform<Float>();
p.y = 1 - 2 * rng.Uniform<Float>();
} while (p.x * p.x + p.y * p.y > 1);
return p;
}
Float SampleCatmullRom(pstd::span<const Float> nodes, pstd::span<const Float> f,
pstd::span<const Float> F, Float u, Float *fval, Float *pdf) {
CHECK_EQ(nodes.size(), f.size());
......@@ -560,11 +560,11 @@ void PiecewiseConstant2D::TestCompareDistributions(const PiecewiseConstant2D &da
Float eps) {
PiecewiseConstant1D::TestCompareDistributions(da.pMarginal, db.pMarginal, eps);
ASSERT_EQ(da.pConditionalY.size(), db.pConditionalY.size());
ASSERT_EQ(da.pConditionalV.size(), db.pConditionalV.size());
ASSERT_EQ(da.domain, db.domain);
for (size_t i = 0; i < da.pConditionalY.size(); ++i)
PiecewiseConstant1D::TestCompareDistributions(da.pConditionalY[i],
db.pConditionalY[i], eps);
for (size_t i = 0; i < da.pConditionalV.size(); ++i)
PiecewiseConstant1D::TestCompareDistributions(da.pConditionalV[i],
db.pConditionalV[i], eps);
}
// AliasTable Method Definitions
......
......@@ -49,6 +49,12 @@ Point2f InvertSphericalRectangleSample(const Point3f &pRef, const Point3f &v00,
PBRT_CPU_GPU
Vector3f SampleHenyeyGreenstein(Vector3f wo, Float g, Point2f u, Float *pdf = nullptr);
PBRT_CPU_GPU inline int SampleDiscrete(pstd::span<const Float> weights, Float u,
Float *pmf = nullptr, Float *uRemapped = nullptr);
PBRT_CPU_GPU inline Float SampleLinear(Float u, Float a, Float b);
PBRT_CPU_GPU inline Float InvertLinearSample(Float x, Float a, Float b);
PBRT_CPU_GPU
Float SampleCatmullRom(pstd::span<const Float> nodes, pstd::span<const Float> f,
pstd::span<const Float> cdf, Float sample, Float *fval = nullptr,
......@@ -61,8 +67,115 @@ Float SampleCatmullRom2D(pstd::span<const Float> nodes1, pstd::span<const Float>
Float *pdf = nullptr);
// Sampling Inline Functions
PBRT_CPU_GPU
inline Float XYZMatchingPDF(Float lambda) {
if (lambda < 360 || lambda > 830)
return 0;
return 0.0039398042f / Sqr(std::cosh(0.0072f * (lambda - 538)));
}
PBRT_CPU_GPU
inline Float SampleXYZMatching(Float u) {
return 538 - 138.888889f * std::atanh(0.85691062f - 1.82750197f * u);
}
PBRT_CPU_GPU inline pstd::array<Float, 3> SampleUniformTriangle(const Point2f &u) {
Float b0, b1;
if (u[0] < u[1]) {
b0 = u[0] / 2;
b1 = u[1] - b0;
} else {
b1 = u[1] / 2;
b0 = u[0] - b1;
}
return {b0, b1, 1 - b0 - b1};
}
PBRT_CPU_GPU
inline Point2f InvertUniformTriangleSample(const pstd::array<Float, 3> &b) {
if (b[0] > b[1]) {
// b0 = u[0] - u[1] / 2, b1 = u[1] / 2
return {b[0] + b[1], 2 * b[1]};
} else {
// b1 = u[1] - u[0] / 2, b0 = u[0] / 2
return {2 * b[0], b[1] + b[0]};
}
}
PBRT_CPU_GPU
inline Float TentPDF(Float x, Float radius) {
if (std::abs(x) >= radius)
return 0;
return 1 / radius - std::abs(x) / Sqr(radius);
}
PBRT_CPU_GPU
inline Float SampleTent(Float u, Float radius) {
if (SampleDiscrete({0.5f, 0.5f}, u, nullptr, &u) == 0)
return -radius + radius * SampleLinear(u, 0, 1);
else
return radius * SampleLinear(u, 1, 0);
}
PBRT_CPU_GPU
inline Float InvertTentSample(Float x, Float radius) {
if (x <= 0)
return (1 - InvertLinearSample(-x / radius, 1, 0)) / 2;
else
return 0.5f + InvertLinearSample(x / radius, 1, 0) / 2;
}
PBRT_CPU_GPU inline Vector3f SampleTrowbridgeReitz(Float alpha_x, Float alpha_y,
Point2f u) {
Float cosTheta, phi;
if (alpha_x == alpha_y) {
// Sample $\cos \theta$ for isotropic Trowbridge--Reitz distribution
Float tanTheta2 = alpha_x * alpha_x * u[0] / (1 - u[0]);
cosTheta = 1 / std::sqrt(1 + tanTheta2);
phi = 2 * Pi * u[1];
} else {
// Sample $\cos \theta$ for anisotropic Trowbridge--Reitz distribution
phi = std::atan(alpha_y / alpha_x * std::tan(2 * Pi * u[1] + .5f * Pi));
if (u[1] > .5f)
phi += Pi;
Float sinPhi = std::sin(phi), cosPhi = std::cos(phi);
Float alpha2 = 1 / (Sqr(cosPhi / alpha_x) + Sqr(sinPhi / alpha_y));
Float tanTheta2 = alpha2 * u[0] / (1 - u[0]);
cosTheta = 1 / std::sqrt(1 + tanTheta2);
}
Float sinTheta = SafeSqrt(1 - Sqr(cosTheta));
return SphericalDirection(sinTheta, cosTheta, phi);
}
// Via Eric Heitz's jcgt sample code...
PBRT_CPU_GPU inline Vector3f SampleTrowbridgeReitzVisibleArea(Vector3f w, Float alpha_x,
Float alpha_y, Point2f u) {
// Transform _w_ to hemispherical configuration for visible area sampling
Vector3f wh = Normalize(Vector3f(alpha_x * w.x, alpha_y * w.y, w.z));
// Find orthonormal basis for visible area microfacet sampling
Vector3f T1 =
(wh.z < 0.99999f) ? Normalize(Cross(Vector3f(0, 0, 1), wh)) : Vector3f(1, 0, 0);
Vector3f T2 = Cross(wh, T1);
// Sample parameterization of projected microfacet area
Float r = std::sqrt(u[0]);
Float phi = 2 * Pi * u[1];
Float t1 = r * std::cos(phi), t2 = r * std::sin(phi);
Float s = 0.5f * (1 + wh.z);
t2 = (1 - s) * std::sqrt(1 - t1 * t1) + s * t2;
// Reproject to hemisphere and transform normal to ellipsoid configuration
Vector3f nh =
t1 * T1 + t2 * T2 + std::sqrt(std::max<Float>(0, 1 - t1 * t1 - t2 * t2)) * wh;
CHECK_RARE(1e-5f, nh.z == 0);
return Normalize(
Vector3f(alpha_x * nh.x, alpha_y * nh.y, std::max<Float>(1e-6f, nh.z)));
}
PBRT_CPU_GPU inline int SampleDiscrete(pstd::span<const Float> weights, Float u,
Float *pmf = nullptr, Float *uRemapped = nullptr) {
Float *pmf, Float *uRemapped) {
// Handle empty _weights_ for discrete sampling
if (weights.empty()) {
if (pmf != nullptr)
......@@ -144,8 +257,9 @@ inline Float InvertNormalSample(Float x, Float mu = 0, Float sigma = 1) {
PBRT_CPU_GPU inline Point2f SampleTwoNormal(const Point2f &u, Float mu = 0,
Float sigma = 1) {
return {mu + sigma * std::sqrt(-2 * std::log(1 - u[0])) * std::cos(2 * Pi * u[1]),
mu + sigma * std::sqrt(-2 * std::log(1 - u[0])) * std::sin(2 * Pi * u[1])};
Float r2 = -2 * std::log(1 - u[0]);
return {mu + sigma * std::sqrt(r2 * std::cos(2 * Pi * u[1])),
mu + sigma * std::sqrt(r2 * std::sin(2 * Pi * u[1]))};
}
PBRT_CPU_GPU inline Float LogisticPDF(Float x, Float s) {
......@@ -350,8 +464,7 @@ PBRT_CPU_GPU inline Point2f SampleBilinear(Point2f u, pstd::span<const Float> w)
DCHECK_EQ(4, w.size());
Point2f p;
// Sample $v$ for bilinear marginal distribution
Float v0 = w[0] + w[1], v1 = w[2] + w[3];
p[1] = SampleLinear(u[1], v0, v1);
p[1] = SampleLinear(u[1], w[0] + w[1], w[2] + w[3]);
// Sample $u$ for bilinear conditional distribution
p[0] = SampleLinear(u[0], Lerp(p[1], w[0], w[2]), Lerp(p[1], w[1], w[3]));
......@@ -370,114 +483,7 @@ PBRT_CPU_GPU inline Float BalanceHeuristic(int nf, Float fPdf, int ng, Float gPd
PBRT_CPU_GPU inline Float PowerHeuristic(int nf, Float fPdf, int ng, Float gPdf) {
Float f = nf * fPdf, g = ng * gPdf;
return (f * f) / (f * f + g * g);
}
PBRT_CPU_GPU
inline Float XYZMatchingPDF(Float lambda) {
if (lambda < 360 || lambda > 830)
return 0;
return 0.0039398042f / Sqr(std::cosh(0.0072f * (lambda - 538)));
}
PBRT_CPU_GPU
inline Float SampleXYZMatching(Float u) {
return 538 - 138.888889f * std::atanh(0.85691062f - 1.82750197f * u);
}
PBRT_CPU_GPU inline pstd::array<Float, 3> SampleUniformTriangle(const Point2f &u) {
Float b0, b1;
if (u[0] < u[1]) {
b0 = u[0] / 2;
b1 = u[1] - b0;
} else {
b1 = u[1] / 2;
b0 = u[0] - b1;
}
return {b0, b1, 1 - b0 - b1};
}
PBRT_CPU_GPU
inline Point2f InvertUniformTriangleSample(const pstd::array<Float, 3> &b) {
if (b[0] > b[1]) {
// b0 = u[0] - u[1] / 2, b1 = u[1] / 2
return {b[0] + b[1], 2 * b[1]};
} else {
// b1 = u[1] - u[0] / 2, b0 = u[0] / 2
return {2 * b[0], b[1] + b[0]};
}
}
PBRT_CPU_GPU
inline Float TentPDF(Float x, Float radius) {
if (std::abs(x) >= radius)
return 0;
return 1 / radius - std::abs(x) / Sqr(radius);
}
PBRT_CPU_GPU
inline Float SampleTent(Float u, Float radius) {
if (SampleDiscrete({0.5f, 0.5f}, u, nullptr, &u) == 0)
return -radius + radius * SampleLinear(u, 0, 1);
else
return radius * SampleLinear(u, 1, 0);
}
PBRT_CPU_GPU
inline Float InvertTentSample(Float x, Float radius) {
if (x <= 0)
return (1 - InvertLinearSample(-x / radius, 1, 0)) / 2;
else
return 0.5f + InvertLinearSample(x / radius, 1, 0) / 2;
}
PBRT_CPU_GPU inline Vector3f SampleTrowbridgeReitz(Float alpha_x, Float alpha_y,
Point2f u) {
Float cosTheta, phi;
if (alpha_x == alpha_y) {
// Sample $\cos \theta$ for isotropic Trowbridge--Reitz distribution
Float tanTheta2 = alpha_x * alpha_x * u[0] / (1 - u[0]);
cosTheta = 1 / std::sqrt(1 + tanTheta2);
phi = 2 * Pi * u[1];
} else {
// Sample $\cos \theta$ for anisotropic Trowbridge--Reitz distribution
phi = std::atan(alpha_y / alpha_x * std::tan(2 * Pi * u[1] + .5f * Pi));
if (u[1] > .5f)
phi += Pi;
Float sinPhi = std::sin(phi), cosPhi = std::cos(phi);
Float alpha2 = 1 / (Sqr(cosPhi / alpha_x) + Sqr(sinPhi / alpha_y));
Float tanTheta2 = alpha2 * u[0] / (1 - u[0]);
cosTheta = 1 / std::sqrt(1 + tanTheta2);
}
Float sinTheta = SafeSqrt(1 - Sqr(cosTheta));
return SphericalDirection(sinTheta, cosTheta, phi);
}
// Via Eric Heitz's jcgt sample code...
PBRT_CPU_GPU inline Vector3f SampleTrowbridgeReitzVisibleArea(Vector3f w, Float alpha_x,
Float alpha_y, Point2f u) {
// Transform _w_ to hemispherical configuration for visible area sampling
Vector3f wh = Normalize(Vector3f(alpha_x * w.x, alpha_y * w.y, w.z));
// Find orthonormal basis for visible area microfacet sampling
Vector3f T1 =
(wh.z < 0.99999f) ? Normalize(Cross(Vector3f(0, 0, 1), wh)) : Vector3f(1, 0, 0);
Vector3f T2 = Cross(wh, T1);
// Sample parameterization of projected microfacet area
Float r = std::sqrt(u[0]);
Float phi = 2 * Pi * u[1];
Float t1 = r * std::cos(phi), t2 = r * std::sin(phi);
Float s = 0.5f * (1 + wh.z);
t2 = (1 - s) * std::sqrt(1 - t1 * t1) + s * t2;
// Reproject to hemisphere and transform normal to ellipsoid configuration
Vector3f nh =
t1 * T1 + t2 * T2 + std::sqrt(std::max<Float>(0, 1 - t1 * t1 - t2 * t2)) * wh;
CHECK_RARE(1e-5f, nh.z == 0);
return Normalize(
Vector3f(alpha_x * nh.x, alpha_y * nh.y, std::max<Float>(1e-6f, nh.z)));
return Sqr(f) / (Sqr(f) + Sqr(g));
}
// Sample from e^(-c x), x from 0 to xMax
......@@ -740,7 +746,7 @@ class PiecewiseConstant2D {
public:
// PiecewiseConstant2D Public Methods
PiecewiseConstant2D() = default;
PiecewiseConstant2D(Allocator alloc) : pConditionalY(alloc), pMarginal(alloc) {}
PiecewiseConstant2D(Allocator alloc) : pConditionalV(alloc), pMarginal(alloc) {}
PiecewiseConstant2D(pstd::span<const Float> data, int nx, int ny,
Allocator alloc = {})
: PiecewiseConstant2D(data, nx, ny, Bounds2f(Point2f(0, 0), Point2f(1, 1)),
......@@ -754,8 +760,8 @@ class PiecewiseConstant2D {
PBRT_CPU_GPU
size_t BytesUsed() const {
return pConditionalY.size() *
(pConditionalY[0].BytesUsed() + sizeof(pConditionalY[0])) +
return pConditionalV.size() *
(pConditionalV[0].BytesUsed() + sizeof(pConditionalV[0])) +
pMarginal.BytesUsed();
}
......@@ -764,33 +770,32 @@ class PiecewiseConstant2D {
PBRT_CPU_GPU
Point2i Resolution() const {
return {int(pConditionalY[0].size()), int(pMarginal.size())};
return {int(pConditionalV[0].size()), int(pMarginal.size())};
}
std::string ToString() const {
return StringPrintf("[ PiecewiseConstant2D domain: %s pConditionalY: %s "
return StringPrintf("[ PiecewiseConstant2D domain: %s pConditionalV: %s "
"pMarginal: %s ]",
domain, pConditionalY, pMarginal);
domain, pConditionalV, pMarginal);
}
static void TestCompareDistributions(const PiecewiseConstant2D &da,
const PiecewiseConstant2D &db, Float eps = 1e-5);
PiecewiseConstant2D(pstd::span<const Float> func, int nx, int ny, Bounds2f domain,
PiecewiseConstant2D(pstd::span<const Float> func, int nu, int nv, Bounds2f domain,
Allocator alloc = {})
: domain(domain), pConditionalY(alloc), pMarginal(alloc) {
CHECK_EQ(func.size(), (size_t)nx * (size_t)ny);
pConditionalY.reserve(ny);
for (int y = 0; y < ny; ++y)
// Compute conditional sampling distribution for $\tilde{y}$
// TODO: emplace_back is key so the alloc sticks. WHY?
pConditionalY.emplace_back(func.subspan(y * nx, nx), domain.pMin[0],
: domain(domain), pConditionalV(alloc), pMarginal(alloc) {
CHECK_EQ(func.size(), (size_t)nu * (size_t)nv);
pConditionalV.reserve(nv);
for (int v = 0; v < nv; ++v)
// Compute conditional sampling distribution for $\tilde{v}$
pConditionalV.emplace_back(func.subspan(v * nu, nu), domain.pMin[0],
domain.pMax[0], alloc);
// Compute marginal sampling distribution $p[\tilde{y}]$
// Compute marginal sampling distribution $p[\tilde{v}]$
std::vector<Float> marginalFunc;
marginalFunc.reserve(ny);
for (int y = 0; y < ny; ++y)
marginalFunc.push_back(pConditionalY[y].funcInt);
marginalFunc.reserve(nv);
for (int v = 0; v < nv; ++v)
marginalFunc.push_back(pConditionalV[v].Integral());
pMarginal =
PiecewiseConstant1D(marginalFunc, domain.pMin[1], domain.pMax[1], alloc);
}
......@@ -801,9 +806,9 @@ class PiecewiseConstant2D {
PBRT_CPU_GPU
Point2f Sample(const Point2f &u, Float *pdf = nullptr) const {
Float pdfs[2];
int y;
Float d1 = pMarginal.Sample(u[1], &pdfs[1], &y);
Float d0 = pConditionalY[y].Sample(u[0], &pdfs[0]);
int v;
Float d1 = pMarginal.Sample(u[1], &pdfs[1], &v);
Float d0 = pConditionalV[v].Sample(u[0], &pdfs[0]);
if (pdf != nullptr)
*pdf = pdfs[0] * pdfs[1];
return Point2f(d0, d1);
......@@ -812,10 +817,10 @@ class PiecewiseConstant2D {
PBRT_CPU_GPU
Float PDF(const Point2f &pr) const {
Point2f p = Point2f(domain.Offset(pr));
int ix =
Clamp(int(p[0] * pConditionalY[0].size()), 0, pConditionalY[0].size() - 1);
int iy = Clamp(int(p[1] * pMarginal.size()), 0, pMarginal.size() - 1);
return pConditionalY[iy].func[ix] / pMarginal.funcInt;
int iu =
Clamp(int(p[0] * pConditionalV[0].size()), 0, pConditionalV[0].size() - 1);
int iv = Clamp(int(p[1] * pMarginal.size()), 0, pMarginal.size() - 1);
return pConditionalV[iv].func[iu] / pMarginal.Integral();
}
PBRT_CPU_GPU
......@@ -826,8 +831,8 @@ class PiecewiseConstant2D {
Float p1o = (p[1] - domain.pMin[1]) / (domain.pMax[1] - domain.pMin[1]);
if (p1o < 0 || p1o > 1)
return {};
int offset = Clamp(p1o * pConditionalY.size(), 0, pConditionalY.size() - 1);
pstd::optional<Float> cInv = pConditionalY[offset].Invert(p[0]);
int offset = Clamp(p1o * pConditionalV.size(), 0, pConditionalV.size() - 1);
pstd::optional<Float> cInv = pConditionalV[offset].Invert(p[0]);
if (!cInv)
return {};
return Point2f(*cInv, *mInv);
......@@ -836,7 +841,7 @@ class PiecewiseConstant2D {
private:
// PiecewiseConstant2D Private Members
Bounds2f domain;
pstd::vector<PiecewiseConstant1D> pConditionalY;
pstd::vector<PiecewiseConstant1D> pConditionalV;
PiecewiseConstant1D pMarginal;
};
......
......@@ -44,6 +44,7 @@ Float SpectrumToPhotometric(SpectrumHandle s) {
for (Float lambda = Lambda_min; lambda <= Lambda_max; ++lambda)
y += Spectra::Y()(lambda) * s(lambda);
const Float K_m = 683;
return y * K_m;
}
......@@ -67,11 +68,10 @@ Float PiecewiseLinearSpectrum::operator()(Float lambda) const {
return 0;
// Find offset to largest _lambdas_ below _lambda_ and interpolate
int offset =
FindInterval(lambdas.size(), [&](int index) { return lambdas[index] <= lambda; });
DCHECK(lambda >= lambdas[offset] && lambda <= lambdas[offset + 1]);
Float t = (lambda - lambdas[offset]) / (lambdas[offset + 1] - lambdas[offset]);
return Lerp(t, values[offset], values[offset + 1]);
int o = FindInterval(lambdas.size(), [&](int i) { return lambdas[i] <= lambda; });
DCHECK(lambda >= lambdas[o] && lambda <= lambdas[o + 1]);
Float t = (lambda - lambdas[o]) / (lambdas[o + 1] - lambdas[o]);
return Lerp(t, values[o], values[o + 1]);
}
Float PiecewiseLinearSpectrum::MaxValue() const {
......@@ -173,17 +173,6 @@ std::string DenselySampledSpectrum::ToString() const {
return s;
}
std::string SampledWavelengths::ToString() const {
std::string r = "[ SampledWavelengths lambda: [";
for (size_t i = 0; i < lambda.size(); ++i)
r += StringPrintf(" %f%c", lambda[i], i != lambda.size() - 1 ? ',' : ' ');
r += "] pdf: [";
for (size_t i = 0; i < lambda.size(); ++i)
r += StringPrintf(" %f%c", pdf[i], i != pdf.size() - 1 ? ',' : ' ');
r += "] ]";
return r;
}
std::string SampledSpectrum::ToString() const {
std::string str = "[ ";
for (int i = 0; i < NSpectrumSamples; ++i) {
......@@ -195,12 +184,25 @@ std::string SampledSpectrum::ToString() const {
return str;
}
std::string SampledWavelengths::ToString() const {
std::string r = "[ SampledWavelengths lambda: [";
for (size_t i = 0; i < lambda.size(); ++i)
r += StringPrintf(" %f%c", lambda[i], i != lambda.size() - 1 ? ',' : ' ');
r += "] pdf: [";
for (size_t i = 0; i < lambda.size(); ++i)
r += StringPrintf(" %f%c", pdf[i], i != pdf.size() - 1 ? ',' : ' ');
r += "] ]";
return r;
}
XYZ SampledSpectrum::ToXYZ(const SampledWavelengths &lambda) const {
// Sample the $X$, $Y$, and $Z$ matching curves at _lambda_
SampledSpectrum X = Spectra::X().Sample(lambda);
SampledSpectrum Y = Spectra::Y().Sample(lambda);
SampledSpectrum Z = Spectra::Z().Sample(lambda);
SampledSpectrum pdf = lambda.PDF();
// Evaluate estimator to compute $(x,y,z)$ coefficients
SampledSpectrum pdf = lambda.PDF();
return XYZ(SafeDiv(X * *this, pdf).Average(), SafeDiv(Y * *this, pdf).Average(),
SafeDiv(Z * *this, pdf).Average()) /
CIE_Y_integral;
......
......@@ -33,7 +33,6 @@ constexpr Float Lambda_min = 360, Lambda_max = 830;
static constexpr int NSpectrumSamples = 4;
static constexpr Float CIE_Y_integral = 106.856895;
static constexpr Float K_m = 683;
// SpectrumHandle Definition
class BlackbodySpectrum;
......@@ -49,7 +48,7 @@ class SpectrumHandle : public TaggedPointer<ConstantSpectrum, DenselySampledSpec
RGBUnboundedSpectrum, RGBIlluminantSpectrum,
BlackbodySpectrum> {
public:
// SpectrumHandle Public Methods
// Spectrum Interface
using TaggedPointer::TaggedPointer;
std::string ToString() const;
......@@ -70,7 +69,7 @@ PBRT_CPU_GPU inline Float Blackbody(Float lambda, Float T) {
const Float c = 299792458;
const Float h = 6.62606957e-34;
const Float kb = 1.3806488e-23;
// Return emitted radiance for blackbody at wavelength _lambda[i]_
// Return emitted radiance for blackbody at wavelength _lambda_
Float l = lambda * 1e-9f;
Float Le = (2 * h * c * c) / (Pow<5>(l) * (FastExp((h * c) / (l * kb * T)) - 1));
CHECK(!IsNaN(Le));
......@@ -78,8 +77,9 @@ PBRT_CPU_GPU inline Float Blackbody(Float lambda, Float T) {
}
namespace Spectra {
DenselySampledSpectrum D(Float temperature, Allocator alloc);
DenselySampledSpectrum D(Float T, Allocator alloc);
} // namespace Spectra
Float SpectrumToPhotometric(SpectrumHandle s);
Float SpectrumToPhotometric(SpectrumHandle s);
......@@ -260,7 +260,6 @@ class SampledSpectrum {
private:
friend class SOA<SampledSpectrum>;
// SampledSpectrum Private Members
pstd::array<Float, NSpectrumSamples> values;
};
......@@ -344,13 +343,16 @@ class SampledWavelengths {
private:
// SampledWavelengths Private Members
friend class SOA<SampledWavelengths>;
pstd::array<Float, NSpectrumSamples> lambda;
pstd::array<Float, NSpectrumSamples> pdf;
pstd::array<Float, NSpectrumSamples> lambda, pdf;
};
// Spectrum Definitions
class ConstantSpectrum {
public:
PBRT_CPU_GPU
ConstantSpectrum(Float c) : c(c) {}
PBRT_CPU_GPU
Float operator()(Float lambda) const { return c; }
// ConstantSpectrum Public Methods
PBRT_CPU_GPU
SampledSpectrum Sample(const SampledWavelengths &) const;
......@@ -360,13 +362,7 @@ class ConstantSpectrum {
std::string ToString() const;
ConstantSpectrum(Float c) : c(c) {}
PBRT_CPU_GPU
Float operator()(Float lambda) const { return c; }
private:
// ConstantSpectrum Private Members
Float c;
};
......@@ -481,8 +477,8 @@ class BlackbodySpectrum {
PBRT_CPU_GPU
BlackbodySpectrum(Float T) : T(T) {
// Compute blackbody normalization constant for given temperature
Float lambdaMax = Float(2.8977721e-3 / T * 1e9);
normalizationFactor = 1 / Blackbody(lambdaMax, T);
Float lambdaMax = Float(2.8977721e-3 / T);
normalizationFactor = 1 / Blackbody(lambdaMax * 1e9f, T);
}
PBRT_CPU_GPU
......@@ -600,11 +596,10 @@ class RGBIlluminantSpectrum {
};
// SampledSpectrum Inline Functions
PBRT_CPU_GPU
inline SampledSpectrum SafeDiv(const SampledSpectrum &s1, const SampledSpectrum &s2) {
PBRT_CPU_GPU inline SampledSpectrum SafeDiv(SampledSpectrum a, SampledSpectrum b) {
SampledSpectrum r;
for (int i = 0; i < NSpectrumSamples; ++i)
r[i] = (s2[i] != 0) ? s1[i] / s2[i] : 0.;
r[i] = (b[i] != 0) ? a[i] / b[i] : 0.;
return r;
}
......
Markdown is supported
0% .
You are about to add 0 people to the discussion. Proceed with caution.
先完成此消息的编辑!
想要评论请 注册