From d0ffa69c9efb923f737c61a5a23e07377e2e6dc9 Mon Sep 17 00:00:00 2001 From: Wenzel Jakob Date: Thu, 23 Jun 2011 19:41:29 +0200 Subject: [PATCH] various cosmetic changes involving dielectric materials, still debugging roughglass.. --- include/mitsuba/core/util.h | 20 +++---- src/bsdfs/dielectric.cpp | 29 +++++----- src/bsdfs/roughglass.cpp | 100 ++++++++++++++++++++--------------- src/libcore/chisquare.cpp | 2 +- src/libcore/util.cpp | 33 ++++++------ src/tests/test_chisquare.cpp | 29 +++++++--- 6 files changed, 118 insertions(+), 95 deletions(-) diff --git a/include/mitsuba/core/util.h b/include/mitsuba/core/util.h index 9cf5cd93..3f181310 100644 --- a/include/mitsuba/core/util.h +++ b/include/mitsuba/core/util.h @@ -423,30 +423,30 @@ extern MTS_EXPORT_CORE Point2 squareToTriangle(const Point2 &sample); * Calculates the unpolarized fresnel reflection coefficient for a * dielectric material * - * \param cosTheta1 + * \param cosThetaI * Cosine of the angle between the normal and the incident ray - * \param cosTheta2 + * \param cosThetaT * Cosine of the angle between the normal and the transmitted ray - * \param etaExt - * Refraction coefficient outside of the material - * \param etaInt - * Refraction coefficient inside the material + * \param etaI + * Refraction coefficient at the incident direction + * \param etaT + * Refraction coefficient at the transmitted direction */ -extern MTS_EXPORT_CORE Float fresnelDielectric(Float cosTheta1, - Float cosTheta2, Float etaExt, Float etaInt); +extern MTS_EXPORT_CORE Float fresnelDielectric(Float cosThetaI, + Float cosThetaT, Float etaI, Float etaT); /** * Calculates the unpolarized fresnel reflection coefficient for a * dielectric material. Handles incidence from either sides. * - * \param cosTheta1 + * \param cosThetaI * Cosine of the angle between the normal and the incident ray * \param etaExt * Refraction coefficient outside of the material * \param etaInt * Refraction coefficient inside the material */ -extern MTS_EXPORT_CORE Float fresnel(Float cosTheta1, Float etaExt, +extern MTS_EXPORT_CORE Float fresnel(Float cosThetaI, Float etaExt, Float etaInt); /** diff --git a/src/bsdfs/dielectric.cpp b/src/bsdfs/dielectric.cpp index 43449f4b..9442584d 100644 --- a/src/bsdfs/dielectric.cpp +++ b/src/bsdfs/dielectric.cpp @@ -97,37 +97,36 @@ public: Float refract(Float intIOR, Float extIOR, const Vector &wi, Vector &wo, ETransportQuantity quantity) const { - Float cosTheta1 = Frame::cosTheta(wi); - bool entering = cosTheta1 > 0.0f; + Float cosThetaI = Frame::cosTheta(wi), + etaI = extIOR, etaT = intIOR; + bool entering = cosThetaI > 0.0f; - /* Swap the indices of refraction if the interaction starts - at the inside of the object */ if (!entering) - std::swap(intIOR, extIOR); + std::swap(etaT, etaI); - Float eta = extIOR/intIOR; + Float eta = etaI / etaT; /* Using Snell's law, calculate the squared sine of the angle between the normal and the transmitted ray */ - Float sinTheta2Sqr = eta*eta * Frame::sinTheta2(wi); + Float sinThetaTSqr = eta*eta * Frame::sinTheta2(wi); - if (sinTheta2Sqr > 1.0f) /* Total internal reflection! */ + if (sinThetaTSqr > 1.0f) /* Total internal reflection! */ return 0.0f; /* Compute the cosine, but guard against numerical imprecision */ - Float cosTheta2 = std::sqrt(std::max((Float) 0.0f, 1.0f - sinTheta2Sqr)); + Float cosThetaT = std::sqrt(1.0f - sinThetaTSqr); if (entering) - cosTheta2 = -cosTheta2; + cosThetaT = -cosThetaT; - /* With cos(N, transmittedRay) on tap, calculating the + /* With cos(N, transmittedRay) avilable, calculating the transmission direction is straightforward */ - wo = Vector(-eta*wi.x, -eta*wi.y, cosTheta2); + wo = Vector(-eta*wi.x, -eta*wi.y, cosThetaT); /* Finally compute transmission coefficient. When transporting - radiance, account for the change at boundaries with different - indices of refraction. */ + radiance, account for the solid angle change at boundaries + with different indices of refraction. */ if (quantity == ERadiance) - return (extIOR*extIOR)/(intIOR*intIOR); + return (etaI*etaI) / (etaT*etaT); else return 1.0f; } diff --git a/src/bsdfs/roughglass.cpp b/src/bsdfs/roughglass.cpp index 24f8a63c..f7418608 100644 --- a/src/bsdfs/roughglass.cpp +++ b/src/bsdfs/roughglass.cpp @@ -155,8 +155,8 @@ public: } /* Prevent numerical issues in other stages of the model */ - if (result < 1e-10) - result = 0; + //if (result < 1e-40) + // result = 0; return result; } @@ -169,7 +169,7 @@ public: * \param v An arbitrary direction */ Float smithG1(const Vector &v, const Vector &m) const { - const Float tanTheta = std::abs(Frame::tanTheta(v)); + const Float tanTheta = std::abs(Frame::tanTheta(v)); Float alpha = m_alpha; /* perpendicular incidence -- no shadowing/masking */ @@ -244,7 +244,9 @@ public: Log(EError, "Invalid distribution function!"); } - return Normal(sphericalDirection(thetaM, phiM)); + Normal result(sphericalDirection(thetaM, phiM)); + Assert(result.z >= 0); + return result; } @@ -256,43 +258,33 @@ public: return (value < 0) ? -1.0f : 1.0f; } - Float refract(const Vector &wi, Vector &wo, ETransportQuantity quantity) const { - Float cosTheta1 = Frame::cosTheta(wi); - Float intIOR = m_intIOR, extIOR = m_extIOR; - bool entering = cosTheta1 > 0.0f; + bool refract(const Vector &wi, Vector &wo) const { + Float cosThetaI = Frame::cosTheta(wi), + etaI = m_extIOR, etaT = m_intIOR; + bool entering = cosThetaI > 0.0f; - /* Swap the indices of refraction if the interaction starts - at the inside of the object */ if (!entering) - std::swap(intIOR, extIOR); + std::swap(etaT, etaI); - Float eta = extIOR/intIOR; + Float eta = etaI / etaT; /* Using Snell's law, calculate the squared sine of the angle between the normal and the transmitted ray */ - Float sinTheta2Sqr = eta*eta * Frame::sinTheta2(wi); + Float sinThetaTSqr = eta*eta * Frame::sinTheta2(wi); - if (sinTheta2Sqr > 1.0f) /* Total internal reflection! */ - return 0.0f; + if (sinThetaTSqr > 1.0f) /* Total internal reflection! */ + return false; - /* Use the sin^2+cos^2=1 identity - max() guards against - numerical imprecision*/ - Float cosTheta2 = std::sqrt(std::max((Float) 0.0f, 1.0f - sinTheta2Sqr)); + /* Compute the cosine, but guard against numerical imprecision */ + Float cosThetaT = std::sqrt(1.0f - sinThetaTSqr); if (entering) - cosTheta2 = -cosTheta2; + cosThetaT = -cosThetaT; - /* Having cos(N, transmittedRay), calculating the actual - direction becomes easy. */ - wo = Vector(-eta*wi.x, -eta*wi.y, cosTheta2); + /* With cos(N, transmittedRay) avilable, calculating the + transmission direction is straightforward */ + wo = Vector(-eta*wi.x, -eta*wi.y, cosThetaT); - /* Finally compute transmission coefficient. When transporting - radiance, account for the change at boundaries with different - indices of refraction. */ - if (quantity == ERadiance) - return (extIOR*extIOR)/(intIOR*intIOR) * - (1.0f - fresnel(Frame::cosTheta(wi), extIOR, intIOR)); - else - return 1.0f - fresnel(Frame::cosTheta(wi), extIOR, intIOR); + return true; } inline Spectrum fReflection(const BSDFQueryRecord &bRec) const { @@ -396,19 +388,19 @@ public: } inline Float pdfTransmission(const BSDFQueryRecord &bRec, Float alpha) const { - Float etaI = m_extIOR, etaO = m_intIOR; + Float etaI = m_extIOR, etaT = m_intIOR; if (Frame::cosTheta(bRec.wi) * Frame::cosTheta(bRec.wo) >= 0) return 0.0f; if (Frame::cosTheta(bRec.wi) < 0) - std::swap(etaI, etaO); + std::swap(etaI, etaT); - Vector Ht = -normalize(bRec.wi*etaI+bRec.wo*etaO); + Vector Ht = -normalize(bRec.wi*etaI + bRec.wo*etaT); /* Jacobian of the half-direction transform. */ - Float sqrtDenom = etaI * dot(bRec.wi, Ht) + etaO * dot(bRec.wo, Ht); - Float dwht_dwo = (etaO*etaO * absDot(bRec.wo, Ht)) / (sqrtDenom*sqrtDenom); + Float sqrtDenom = etaI * dot(bRec.wi, Ht) + etaT * dot(bRec.wo, Ht); + Float dwht_dwo = (etaT*etaT * absDot(bRec.wo, Ht)) / (sqrtDenom*sqrtDenom); return evalD(Ht, alpha) * std::abs(Frame::cosTheta(Ht)) * dwht_dwo; } @@ -422,11 +414,12 @@ public: /* Suggestion by Bruce Walter: sample using a slightly different value of alpha. This in practice limits the weights to values <= 4. See also \ref sample() */ - Float alpha = m_alpha * (1.2f - 0.2f * std::sqrt( - std::abs(Frame::cosTheta(bRec.wi)))); +// Float alpha = m_alpha * (1.2f - 0.2f * std::sqrt( +// std::abs(Frame::cosTheta(bRec.wi)))); + Float alpha = m_alpha; if (hasReflection && hasTransmission) { - /* PDF for importance sampling according to approximate F + /* PDF for importance sampling according to approximate Fresnel coefficients (approximate, because we don't know the microfacet normal at this point) */ Float fr = fresnel(Frame::cosTheta(bRec.wi), m_extIOR, m_intIOR); @@ -466,14 +459,31 @@ public: } inline Spectrum sampleTransmission(BSDFQueryRecord &bRec, Float alpha, const Point2 &sample) const { + Vector m = sampleD(sample, alpha); + cout << endl; +#if 1 + Float etaI = m_extIOR, etaT = m_intIOR; + + if (Frame::cosTheta(bRec.wi) < 0) + std::swap(etaI, etaT); + + Float eta = etaI / etaT, c = dot(bRec.wi, m); + + Vector wo2 = m * (eta*c - signum(bRec.wi.z) * std::sqrt(1 + eta * (c*c-1))) + - bRec.wi * eta; +#endif + cout << "wo2: " << wo2.toString() << endl; +#if 1 /* Sample M, the microsurface normal */ - Frame mFrame(sampleD(sample, alpha)); + Frame mFrame(m); /* Perfect specular reflection along the microsurface normal */ - if (refract(mFrame.toLocal(bRec.wi), bRec.wo, bRec.quantity) == 0) + if (!refract(mFrame.toLocal(bRec.wi), bRec.wo)) return Spectrum(0.0f); bRec.wo = mFrame.toWorld(bRec.wo); +#endif + cout << "bRec.wo: " << bRec.wo.toString() << endl; bRec.sampledComponent = 1; bRec.sampledType = EGlossyTransmission; @@ -502,11 +512,13 @@ public: value of alpha. This in practice limits the weights to values <= 4. The change is of course also accounted for in \ref pdf(), hence no error is introduced. */ - Float alpha = m_alpha * (1.2f - 0.2f * std::sqrt( - std::abs(Frame::cosTheta(bRec.wi)))); + //Float alpha = m_alpha * (1.2f - 0.2f * std::sqrt( + // std::abs(Frame::cosTheta(bRec.wi)))); + + Float alpha = m_alpha; if (hasReflection && hasTransmission) { - /* PDF for importance sampling according to approximate F + /* PDF for importance sampling according to approximate Fresnel coefficients (approximate, because we don't know the microfacet normal at this point) */ Float fr = fresnel(Frame::cosTheta(bRec.wi), m_extIOR, m_intIOR); @@ -515,7 +527,7 @@ public: sample.x /= fr; return sampleReflection(bRec, alpha, sample); } else { - sample.x = (sample.x - fr) / (1-fr); + sample.x = (sample.x - fr) / (1 - fr); return sampleTransmission(bRec, alpha, sample); } } else if (hasReflection) { diff --git a/src/libcore/chisquare.cpp b/src/libcore/chisquare.cpp index c5ba934e..29e14277 100644 --- a/src/libcore/chisquare.cpp +++ b/src/libcore/chisquare.cpp @@ -91,7 +91,7 @@ void ChiSquare::fill( Float min[2], max[2]; size_t idx = 0; - NDIntegrator integrator(1, 2, 100000, 0, 1e-6f); + NDIntegrator integrator(1, 2, 1000000, 1e-8f, 1e-8f); Float maxError = 0, integral = 0; for (int i=0; i 1.0f) + if (sinThetaT > 1.0f) return 1.0f; /* Total internal reflection! */ - /* Use the sin^2+cos^2=1 identity - max() guards against - numerical imprecision*/ - Float cosTheta2 = std::sqrt(std::max((Float) 0.0f, - 1.0f - sinTheta2*sinTheta2)); + Float cosThetaT = std::sqrt(1.0f - sinThetaT*sinThetaT); /* Finally compute the reflection coefficient */ - return fresnelDielectric(std::abs(cosTheta1), cosTheta2, - etaInt, etaExt); + return fresnelDielectric(std::abs(cosThetaI), + cosThetaT, etaI, etaT); } Float radicalInverse(int b, size_t i) { diff --git a/src/tests/test_chisquare.cpp b/src/tests/test_chisquare.cpp index 42afe007..509ab262 100644 --- a/src/tests/test_chisquare.cpp +++ b/src/tests/test_chisquare.cpp @@ -38,7 +38,8 @@ public: class BSDFAdapter { public: BSDFAdapter(const BSDF *bsdf, Random *random, const Vector &wi, int component = -1) - : m_bsdf(bsdf), m_random(random), m_wi(wi), m_component(component) { } + : m_bsdf(bsdf), m_random(random), m_wi(wi), m_component(component), + m_largestWeight(0) { } std::pair generateSample() { Point2 sample(m_random->nextFloat(), m_random->nextFloat()); @@ -69,15 +70,17 @@ public: for (int i=0; i= 0 && b >= 0); + Float min = std::min(a, b); + Float err = std::abs(a - b); + m_largestWeight = std::max(m_largestWeight, a); if (min < Epsilon && err > Epsilon) // absolute error threshold mismatch = true; else if (min > Epsilon && err/min > Epsilon) // relative error threshold mismatch = true; } - + if (mismatch) Log(EWarn, "Inconsistency: f=%s, pdf=%f, sampled f/pdf=%s", f.toString().c_str(), pdfVal, sampled.toString().c_str()); @@ -95,11 +98,14 @@ public: return 0.0f; return m_bsdf->pdf(bRec); } + + inline Float getLargestWeight() const { return m_largestWeight; } private: ref m_bsdf; ref m_random; Vector m_wi; int m_component; + Float m_largestWeight; }; void test01_BSDF() { @@ -116,8 +122,10 @@ public: continue; const BSDF *bsdf = static_cast(objects[i]); + Float largestWeight = 0; Log(EInfo, "Processing BSDF model %s", bsdf->toString().c_str()); +#if 0 Log(EInfo, "Checking the combined model for %i incident directions", wiSamples); /* Test for a number of different incident directions */ @@ -135,8 +143,8 @@ public: // Initialize the tables used by the chi-square test chiSqr->fill( - boost::bind(&BSDFAdapter::generateSample, adapter), - boost::bind(&BSDFAdapter::pdf, adapter, _1) + boost::bind(&BSDFAdapter::generateSample, &adapter), + boost::bind(&BSDFAdapter::pdf, &adapter, _1) ); // (the following assumes that the distribution has 1 parameter, e.g. exponent value) @@ -150,8 +158,10 @@ public: } else { succeed(); } + largestWeight = std::max(largestWeight, adapter.getLargestWeight()); ++testCount; } +#endif if (bsdf->getComponentCount() > 1) { for (int comp=0; compgetComponentCount(); ++comp) { @@ -167,13 +177,14 @@ public: wi = squareToHemispherePSA(Point2(random->nextFloat(), random->nextFloat())); BSDFAdapter adapter(bsdf, random, wi, comp); + ref chiSqr = new ChiSquare(thetaBins, 2*thetaBins, wiSamples); chiSqr->setLogLevel(EDebug); // Initialize the tables used by the chi-square test chiSqr->fill( - boost::bind(&BSDFAdapter::generateSample, adapter), - boost::bind(&BSDFAdapter::pdf, adapter, _1) + boost::bind(&BSDFAdapter::generateSample, &adapter), + boost::bind(&BSDFAdapter::pdf, &adapter, _1) ); // (the following assumes that the distribution has 1 parameter, e.g. exponent value) @@ -187,10 +198,12 @@ public: } else { succeed(); } + largestWeight = std::max(largestWeight, adapter.getLargestWeight()); ++testCount; } } } + Log(EInfo, "Done with this model. The largest encountered sampling weight was = %f", largestWeight); } Log(EInfo, "%i/%i BSDF checks succeeded", testCount-failureCount, testCount); }