/*
This file is part of Mitsuba, a physically based rendering system.
Copyright (c) 2007-2014 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 <http://www.gnu.org/licenses/>.
*/
#if !defined(__MICROFACET_H)
#define __MICROFACET_H
#include <mitsuba/mitsuba.h>
#include <mitsuba/core/frame.h>
#include <mitsuba/core/properties.h>
#include <boost/algorithm/string.hpp>
MTS_NAMESPACE_BEGIN
/**
* \brief Implementation of the Beckman and GGX / Trowbridge-Reitz microfacet
* distributions and various useful sampling routines
*
* Based on the papers
*
* "Microfacet Models for Refraction through Rough Surfaces"
* by Bruce Walter, Stephen R. Marschner, Hongsong Li, and Kenneth E. Torrance
*
* and
*
* "Importance Sampling Microfacet-Based BSDFs using the Distribution of Visible Normals"
* by Eric Heitz and Eugene D'Eon
*
* The visible normal sampling code was provided by Eric Heitz and Eugene D'Eon.
*/
class MicrofacetDistribution {
public:
/// Supported distribution types
enum EType {
/// Beckmann distribution derived from Gaussian random surfaces
EBeckmann = 0,
/// GGX: Long-tailed distribution for very rough surfaces (aka. Trowbridge-Reitz distr.)
EGGX = 1,
/// Phong distribution (with the anisotropic extension by Ashikhmin and Shirley)
EPhong = 2
};
/**
* Create an isotropic microfacet distribution of the specified type
*
* \param type
* The desired type of microfacet distribution
* \param alpha
* The surface roughness
*/
inline MicrofacetDistribution(EType type, Float alpha, bool sampleVisible = true)
: m_type(type), m_alphaU(alpha), m_alphaV(alpha), m_sampleVisible(sampleVisible),
m_exponentU(0.0f), m_exponentV(0.0f) {
if (m_type == EPhong)
computePhongExponent();
}
/**
* Create an anisotropic microfacet distribution of the specified type
*
* \param type
* The desired type of microfacet distribution
* \param alphaU
* The surface roughness in the tangent direction
* \param alphaV
* The surface roughness in the bitangent direction
*/
inline MicrofacetDistribution(EType type, Float alphaU, Float alphaV, bool sampleVisible = true)
: m_type(type), m_alphaU(alphaU), m_alphaV(alphaV), m_sampleVisible(sampleVisible),
m_exponentU(0.0f), m_exponentV(0.0f) {
if (m_type == EPhong)
computePhongExponent();
}
/**
* \brief Create a microfacet distribution from a Property data
* structure
*/
MicrofacetDistribution(const Properties &props, EType type = EBeckmann,
Float alphaU = 0.1f, Float alphaV = 0.1f, bool sampleVisible = true)
: m_type(type), m_alphaU(alphaU), m_alphaV(alphaV), m_exponentU(0.0f),
m_exponentV(0.0f) {
if (props.hasProperty("distribution")) {
std::string distr = boost::to_lower_copy(props.getString("distribution"));
if (distr == "beckmann")
m_type = EBeckmann;
else if (distr == "ggx")
m_type = EGGX;
else if (distr == "phong" || distr == "as")
m_type = EPhong;
else
SLog(EError, "Specified an invalid distribution \"%s\", must be "
"\"beckmann\", \"ggx\", or \"phong\"/\"as\"!", distr.c_str());
}
if (props.hasProperty("alpha")) {
m_alphaU = m_alphaV = props.getFloat("alpha");
if (props.hasProperty("alphaU") || props.hasProperty("alphaV"))
SLog(EError, "Microfacet model: please specify either 'alpha' or 'alphaU'/'alphaV'.");
} else if (props.hasProperty("alphaU") || props.hasProperty("alphaV")) {
if (!props.hasProperty("alphaU") || !props.hasProperty("alphaV"))
SLog(EError, "Microfacet model: both 'alphaU' and 'alphaV' must be specified.");
if (props.hasProperty("alpha"))
SLog(EError, "Microfacet model: please specify either 'alpha' or 'alphaU'/'alphaV'.");
m_alphaU = props.getFloat("alphaU");
m_alphaV = props.getFloat("alphaV");
}
m_sampleVisible = props.getBoolean("sampleVisible", sampleVisible);
/* Visible normal sampling is not supported for the Phong / Ashikhmin-Shirley distribution */
if (m_type == EPhong) {
m_sampleVisible = false;
computePhongExponent();
}
}
/// Return the distribution type
inline EType getType() const { return m_type; }
/// Return the roughness (isotropic case)
inline Float getAlpha() const { return m_alphaU; }
/// Return the roughness along the tangent direction
inline Float getAlphaU() const { return m_alphaU; }
/// Return the roughness along the bitangent direction
inline Float getAlphaV() const { return m_alphaV; }
/// Return the Phong exponent (isotropic case)
inline Float getExponent() const { return m_exponentU; }
/// Return the Phong exponent along the tangent direction
inline Float getExponentU() const { return m_exponentU; }
/// Return the Phong exponent along the bitangent direction
inline Float getExponentV() const { return m_exponentV; }
/// Return whether or not only visible normals are sampled?
inline bool getSampleVisible() const { return m_sampleVisible; }
/// Is this an anisotropic microfacet distribution?
inline bool isAnisotropic() const { return m_alphaU != m_alphaV; }
/// Is this an anisotropic microfacet distribution?
inline bool isIsotropic() const { return m_alphaU == m_alphaV; }
/// Scale the roughness values by some constant
inline void scaleAlpha(Float value) {
m_alphaU *= value;
m_alphaV *= value;
if (m_type == EPhong)
computePhongExponent();
}
/**
* \brief Evaluate the microfacet distribution function
*
* \param m
* The microfacet normal
*/
inline Float eval(const Vector &m) const {
if (Frame::cosTheta(m) <= 0)
return 0.0f;
Float cosTheta2 = Frame::cosTheta2(m);
Float beckmannExponent = ((m.x*m.x) / (m_alphaU * m_alphaU)
+ (m.y*m.y) / (m_alphaV * m_alphaV)) / cosTheta2;
Float result;
switch (m_type) {
case EBeckmann: {
/* Beckmann distribution function for Gaussian random surfaces - [Walter 2005] evaluation */
result = math::fastexp(-beckmannExponent) /
(M_PI * m_alphaU * m_alphaV * cosTheta2 * cosTheta2);
}
break;
case EGGX: {
/* GGX / Trowbridge-Reitz distribution function for rough surfaces */
Float root = ((Float) 1 + beckmannExponent) * cosTheta2;
result = (Float) 1 / (M_PI * m_alphaU * m_alphaV * root * root);
}
break;
case EPhong: {
/* Isotropic case: Phong distribution. Anisotropic case: Ashikhmin-Shirley distribution */
Float exponent = interpolatePhongExponent(m);
result = std::sqrt((m_exponentU + 2) * (m_exponentV + 2))
* INV_TWOPI * std::pow(Frame::cosTheta(m), exponent);
}
break;
default:
SLog(EError, "Invalid distribution type!");
return -1;
}
/* Prevent potential numerical issues in other stages of the model */
if (result * Frame::cosTheta(m) < 1e-20f)
result = 0;
return result;
}
/**
* \brief Wrapper function which calls \ref sampleAll() or \ref sampleVisible()
* depending on the parameters of this class
*/
inline Normal sample(const Vector &wi, const Point2 &sample, Float &pdf) const {
Normal m;
if (m_sampleVisible) {
m = sampleVisible(wi, sample);
pdf = pdfVisible(wi, m);
} else {
m = sampleAll(sample, pdf);
}
return m;
}
/**
* \brief Wrapper function which calls \ref sampleAll() or \ref sampleVisible()
* depending on the parameters of this class
*/
inline Normal sample(const Vector &wi, const Point2 &sample) const {
Normal m;
if (m_sampleVisible) {
m = sampleVisible(wi, sample);
} else {
Float pdf;
m = sampleAll(sample, pdf);
}
return m;
}
/**
* \brief Wrapper function which calls \ref pdfAll() or \ref pdfVisible()
* depending on the parameters of this class
*/
inline Float pdf(const Vector &wi, const Vector &m) const {
if (m_sampleVisible)
return pdfVisible(wi, m);
else
return pdfAll(m);
}
/**
* \brief Draw a sample from the microfacet normal distribution
* (including *all* normals) and return the associated
* probability density
*
* \param sample
* A uniformly distributed 2D sample
* \param pdf
* The probability density wrt. solid angles
*/
inline Normal sampleAll(const Point2 &sample, Float &pdf) const {
/* The azimuthal component is always selected
uniformly regardless of the distribution */
Float cosThetaM = 0.0f;
Float sinPhiM, cosPhiM;
Float alphaSqr;
switch (m_type) {
case EBeckmann: {
/* Beckmann distribution function for Gaussian random surfaces */
if (isIsotropic()) {
/* Sample phi component (isotropic case) */
math::sincos((2.0f * M_PI) * sample.y, &sinPhiM, &cosPhiM);
alphaSqr = m_alphaU * m_alphaU;
} else {
/* Sample phi component (anisotropic case) */
Float phiM = std::atan(m_alphaV / m_alphaU *
std::tan(M_PI + 2*M_PI*sample.y)) + M_PI * std::floor(2*sample.y + 0.5f);
math::sincos(phiM, &sinPhiM, &cosPhiM);
Float cosSc = cosPhiM / m_alphaU, sinSc = sinPhiM / m_alphaV;
alphaSqr = 1.0f / (cosSc*cosSc + sinSc*sinSc);
}
/* Sample theta component */
Float tanThetaMSqr = alphaSqr * -math::fastlog(1.0f - sample.x);
cosThetaM = 1.0f / std::sqrt(1.0f + tanThetaMSqr);
/* Compute probability density of the sampled position */
pdf = (1.0f - sample.x) / (M_PI*m_alphaU*m_alphaV*cosThetaM*cosThetaM*cosThetaM);
}
break;
case EGGX: {
/* GGX / Trowbridge-Reitz distribution function for rough surfaces */
if (isIsotropic()) {
/* Sample phi component (isotropic case) */
math::sincos((2.0f * M_PI) * sample.y, &sinPhiM, &cosPhiM);
/* Sample theta component */
alphaSqr = m_alphaU*m_alphaU;
} else {
/* Sample phi component (anisotropic case) */
Float phiM = std::atan(m_alphaV / m_alphaU *
std::tan(M_PI + 2*M_PI*sample.y)) + M_PI * std::floor(2*sample.y + 0.5f);
math::sincos(phiM, &sinPhiM, &cosPhiM);
Float cosSc = cosPhiM / m_alphaU, sinSc = sinPhiM / m_alphaV;
alphaSqr = 1.0f / (cosSc*cosSc + sinSc*sinSc);
}
/* Sample theta component */
Float tanThetaMSqr = alphaSqr * sample.x / (1.0f - sample.x);
cosThetaM = 1.0f / std::sqrt(1.0f + tanThetaMSqr);
/* Compute probability density of the sampled position */
Float temp = 1+tanThetaMSqr/alphaSqr;
pdf = INV_PI / (m_alphaU*m_alphaV*cosThetaM*cosThetaM*cosThetaM*temp*temp);
}
break;
case EPhong: {
Float phiM;
Float exponent;
if (isIsotropic()) {
phiM = (2.0f * M_PI) * sample.y;
exponent = m_exponentU;
} else {
/* Sampling method based on code from PBRT */
if (sample.y < 0.25f) {
sampleFirstQuadrant(4 * sample.y, phiM, exponent);
} else if (sample.y < 0.5f) {
sampleFirstQuadrant(4 * (0.5f - sample.y), phiM, exponent);
phiM = M_PI - phiM;
} else if (sample.y < 0.75f) {
sampleFirstQuadrant(4 * (sample.y - 0.5f), phiM, exponent);
phiM += M_PI;
} else {
sampleFirstQuadrant(4 * (1 - sample.y), phiM, exponent);
phiM = 2 * M_PI - phiM;
}
}
math::sincos(phiM, &sinPhiM, &cosPhiM);
cosThetaM = std::pow(sample.x, 1.0f / (exponent + 2.0f));
pdf = std::sqrt((m_exponentU + 2.0f) * (m_exponentV + 2.0f))
* INV_TWOPI * std::pow(cosThetaM, exponent + 1.0f);
}
break;
default:
SLog(EError, "Invalid distribution type!");
pdf = -1;
return Vector(-1);
}
/* Prevent potential numerical issues in other stages of the model */
if (pdf < 1e-20f)
pdf = 0;
Float sinThetaM = std::sqrt(
std::max((Float) 0, 1 - cosThetaM*cosThetaM));
return Vector(
sinThetaM * cosPhiM,
sinThetaM * sinPhiM,
cosThetaM
);
}
/**
* \brief Returns the density function associated with
* the \ref sampleAll() function.
*
* \param m
* The microfacet normal
*/
inline Float pdfAll(const Vector &m) const {
/* PDF is just D(m) * cos(theta_M) */
return eval(m) * Frame::cosTheta(m);
}
/**
* \brief Draw a sample from the distribution of visible normals
* and return the associated probability density
*
* \param _wi
* A reference direction that defines the set of visible normals
* \param sample
* A uniformly distributed 2D sample
* \param pdf
* The probability density wrt. solid angles
*/
inline Normal sampleVisible(const Vector &_wi, const Point2 &sample) const {
/* Step 1: stretch wi */
Vector wi = normalize(Vector(
m_alphaU * _wi.x,
m_alphaV * _wi.y,
_wi.z
));
/* Get polar coordinates */
Float theta = 0, phi = 0;
if (wi.z < (Float) 0.99999) {
theta = std::acos(wi.z);
phi = std::atan2(wi.y, wi.x);
}
Float sinPhi, cosPhi;
math::sincos(phi, &sinPhi, &cosPhi);
/* Step 2: simulate P22_{wi}(slope.x, slope.y, 1, 1) */
Vector2 slope = sampleVisible11(theta, sample);
/* Step 3: rotate */
slope = Vector2(
cosPhi * slope.x - sinPhi * slope.y,
sinPhi * slope.x + cosPhi * slope.y);
/* Step 4: unstretch */
slope.x *= m_alphaU;
slope.y *= m_alphaV;
/* Step 5: compute normal */
Float normalization = (Float) 1 / std::sqrt(slope.x*slope.x
+ slope.y*slope.y + (Float) 1.0);
return Normal(
-slope.x * normalization,
-slope.y * normalization,
normalization
);
}
/// Implements the probability density of the function \ref sampleVisible()
Float pdfVisible(const Vector &wi, const Vector &m) const {
return smithG1(wi, m) * absDot(wi, m) * eval(m) / std::abs(Frame::cosTheta(wi));
}
/**
* \brief Smith's shadowing-masking function G1 for each
* of the supported microfacet distributions
*
* \param v
* An arbitrary direction
* \param m
* The microfacet normal
*/
Float smithG1(const Vector &v, const Vector &m) const {
/* Ensure consistent orientation (can't see the back
of the microfacet from the front and vice versa) */
if (dot(v, m) * Frame::cosTheta(v) <= 0)
return 0.0f;
/* Perpendicular incidence -- no shadowing/masking */
Float tanTheta = std::abs(Frame::tanTheta(v));
if (tanTheta == 0.0f)
return 1.0f;
Float alpha = projectRoughness(v);
switch (m_type) {
case EPhong:
case EBeckmann: {
Float a = 1.0f / (alpha * tanTheta);
if (a >= 1.6f)
return 1.0f;
/* Use a fast and accurate (<0.35% rel. error) rational
approximation to the shadowing-masking function */
Float aSqr = a*a;
return (3.535f * a + 2.181f * aSqr)
/ (1.0f + 2.276f * a + 2.577f * aSqr);
}
break;
case EGGX: {
Float root = alpha * tanTheta;
return 2.0f / (1.0f + std::sqrt(1.0f + root*root));
}
break;
default:
SLog(EError, "Invalid distribution type!");
return -1.0f;
}
}
/**
* \brief Separable shadow-masking function based on Smith's
* one-dimensional masking model
*/
inline Float G(const Vector &wi, const Vector &wo, const Vector &m) const {
return smithG1(wi, m) * smithG1(wo, m);
}
/// Return a string representation of the name of a distribution
inline static std::string distributionName(EType type) {
switch (type) {
case EBeckmann: return "beckmann"; break;
case EGGX: return "ggx"; break;
case EPhong: return "phong"; break;
default: return "invalid"; break;
}
}
/// Return a string representation of the contents of this instance
std::string toString() const {
return formatString("MicrofacetDistribution[type=\"%s\", alphaU=%f, alphaV=%f]",
distributionName(m_type).c_str(), m_alphaU, m_alphaV);
}
protected:
/// Compute the effective roughness projected on direction \c v
inline Float projectRoughness(const Vector &v) const {
Float invSinTheta2 = 1 / Frame::sinTheta2(v);
if (isIsotropic() || invSinTheta2 <= 0)
return m_alphaU;
Float cosPhi2 = v.x * v.x * invSinTheta2;
Float sinPhi2 = v.y * v.y * invSinTheta2;
return std::sqrt(cosPhi2 * m_alphaU * m_alphaU + sinPhi2 * m_alphaV * m_alphaV);
}
/// Compute the interpolated roughness for the Phong model
inline Float interpolatePhongExponent(const Vector &v) const {
Float invSinTheta2 = 1 / Frame::sinTheta2(v);
if (isIsotropic() || invSinTheta2 <= 0)
return m_exponentU;
Float cosPhi2 = v.x * v.x * invSinTheta2;
Float sinPhi2 = v.y * v.y * invSinTheta2;
return m_exponentU * cosPhi2 + m_exponentV * sinPhi2;
}
/**
* \brief Visible normal sampling code for the alpha=1 case
*
* Source: supplemental material of "Importance Sampling
* Microfacet-Based BSDFs using the Distribution of Visible Normals"
*/
Vector2 sampleVisible11(Float thetaI, Point2 sample) const {
const Float SQRT_PI_INV = 1 / std::sqrt(M_PI);
Vector2 slope;
switch (m_type) {
case EBeckmann: {
/* Special case (normal incidence) */
if (thetaI < 1e-4f) {
Float sinPhi, cosPhi;
Float r = std::sqrt(-math::fastlog(1.0f-sample.x));
math::sincos(2 * M_PI * sample.y, &sinPhi, &cosPhi);
return Vector2(r * cosPhi, r * sinPhi);
}
/* The original inversion routine from the paper contained
discontinuities, which causes issues for QMC integration
and techniques like Kelemen-style MLT. The following code
performs a numerical inversion with better behavior */
Float tanThetaI = std::tan(thetaI);
Float cotThetaI = 1 / tanThetaI;
/* Search interval -- everything is parameterized
in the erf() domain */
Float a = -1, c = math::erf(cotThetaI);
Float sample_x = std::max(sample.x, (Float) 1e-6f);
/* Start with a good initial guess */
//Float b = (1-sample_x) * a + sample_x * c;
/* We can do better (inverse of an approximation computed in Mathematica) */
Float fit = 1 + thetaI*(-0.876f + thetaI * (0.4265f - 0.0594f*thetaI));
Float b = c - (1+c) * std::pow(1-sample_x, fit);
/* Normalization factor for the CDF */
Float normalization = 1 / (1 + c + SQRT_PI_INV*
tanThetaI*std::exp(-cotThetaI*cotThetaI));
int it = 0;
while (++it < 10) {
/* Bisection criterion -- the oddly-looking
boolean expression are intentional to check
for NaNs at little additional cost */
if (!(b >= a && b <= c))
b = 0.5f * (a + c);
/* Evaluate the CDF and its derivative
(i.e. the density function) */
Float invErf = math::erfinv(b);
Float value = normalization*(1 + b + SQRT_PI_INV*
tanThetaI*std::exp(-invErf*invErf)) - sample_x;
Float derivative = normalization * (1
- invErf*tanThetaI);
if (std::abs(value) < 1e-5f)
break;
/* Update bisection intervals */
if (value > 0)
c = b;
else
a = b;
b -= value / derivative;
}
/* Now convert back into a slope value */
slope.x = math::erfinv(b);
/* Simulate Y component */
slope.y = math::erfinv(2.0f*std::max(sample.y, (Float) 1e-6f) - 1.0f);
};
break;
case EGGX: {
/* Special case (normal incidence) */
if (thetaI < 1e-4f) {
Float sinPhi, cosPhi;
Float r = math::safe_sqrt(sample.x / (1 - sample.x));
math::sincos(2 * M_PI * sample.y, &sinPhi, &cosPhi);
return Vector2(r * cosPhi, r * sinPhi);
}
/* Precomputations */
Float tanThetaI = std::tan(thetaI);
Float a = 1 / tanThetaI;
Float G1 = 2.0f / (1.0f + math::safe_sqrt(1.0f + 1.0f / (a*a)));
/* Simulate X component */
Float A = 2.0f * sample.x / G1 - 1.0f;
if (std::abs(A) == 1)
A -= math::signum(A)*Epsilon;
Float tmp = 1.0f / (A*A - 1.0f);
Float B = tanThetaI;
Float D = math::safe_sqrt(B*B*tmp*tmp - (A*A - B*B) * tmp);
Float slope_x_1 = B * tmp - D;
Float slope_x_2 = B * tmp + D;
slope.x = (A < 0.0f || slope_x_2 > 1.0f / tanThetaI) ? slope_x_1 : slope_x_2;
/* Simulate Y component */
Float S;
if (sample.y > 0.5f) {
S = 1.0f;
sample.y = 2.0f * (sample.y - 0.5f);
} else {
S = -1.0f;
sample.y = 2.0f * (0.5f - sample.y);
}
/* Improved fit */
Float z =
(sample.y * (sample.y * (sample.y * (-(Float) 0.365728915865723) + (Float) 0.790235037209296) -
(Float) 0.424965825137544) + (Float) 0.000152998850436920) /
(sample.y * (sample.y * (sample.y * (sample.y * (Float) 0.169507819808272 - (Float) 0.397203533833404) -
(Float) 0.232500544458471) + (Float) 1) - (Float) 0.539825872510702);
slope.y = S * z * std::sqrt(1.0f + slope.x*slope.x);
};
break;
default:
SLog(EError, "Invalid distribution type!");
return Vector2(-1);
};
return slope;
}
/// Helper routine: convert from Beckmann-style roughness values to Phong exponents (Walter et al.)
void computePhongExponent() {
m_exponentU = std::max(2.0f / (m_alphaU * m_alphaU) - 2.0f, (Float) 0.0f);
m_exponentV = std::max(2.0f / (m_alphaV * m_alphaV) - 2.0f, (Float) 0.0f);
}
/// Helper routine: sample the azimuthal part of the first quadrant of the A&S distribution
void sampleFirstQuadrant(Float u1, Float &phi, Float &exponent) const {
Float cosPhi, sinPhi;
phi = std::atan(
std::sqrt((m_exponentU + 2.0f) / (m_exponentV + 2.0f)) *
std::tan(M_PI * u1 * 0.5f));
math::sincos(phi, &sinPhi, &cosPhi);
/* Return the interpolated roughness */
exponent = m_exponentU * cosPhi * cosPhi + m_exponentV * sinPhi * sinPhi;
}
protected:
EType m_type;
Float m_alphaU, m_alphaV;
bool m_sampleVisible;
Float m_exponentU, m_exponentV;
};
MTS_NAMESPACE_END
#endif /* __MICROFACET_H */