808 lines
20 KiB
C++
808 lines
20 KiB
C++
/*
|
|
This file is part of Mitsuba, a physically based rendering system.
|
|
|
|
Copyright (c) 2007-2012 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/>.
|
|
*/
|
|
|
|
#include <mitsuba/mitsuba.h>
|
|
#include <mitsuba/core/util.h>
|
|
#include <mitsuba/core/random.h>
|
|
#include <mitsuba/core/quad.h>
|
|
#include <mitsuba/core/sse.h>
|
|
#include <boost/bind.hpp>
|
|
#include <stdarg.h>
|
|
#include <iomanip>
|
|
#include <errno.h>
|
|
|
|
#if defined(__OSX__)
|
|
#include <sys/sysctl.h>
|
|
#elif defined(WIN32)
|
|
#include <direct.h>
|
|
#else
|
|
#include <malloc.h>
|
|
#endif
|
|
|
|
#if defined(WIN32)
|
|
# include <windows.h>
|
|
# include <winsock2.h>
|
|
# include <ws2tcpip.h>
|
|
#else
|
|
# include <sys/types.h>
|
|
# include <sys/socket.h>
|
|
# include <netdb.h>
|
|
# include <fenv.h>
|
|
#endif
|
|
|
|
// SSE is not enabled in general when using double precision, however it is
|
|
// required in OS X for FP exception handling
|
|
#if defined(__OSX__) && !defined(MTS_SSE)
|
|
#include <xmmintrin.h>
|
|
#undef enable_fpexcept_sse
|
|
#undef query_fpexcept_sse
|
|
#undef disable_fpexcept_sse
|
|
|
|
namespace {
|
|
inline void enable_fpexcept_sse() {
|
|
_MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() &
|
|
~(_MM_MASK_INVALID | _MM_MASK_DIV_ZERO));
|
|
}
|
|
inline unsigned int query_fpexcept_sse() {
|
|
return (~_MM_GET_EXCEPTION_MASK() &
|
|
(_MM_MASK_INVALID | _MM_MASK_DIV_ZERO));
|
|
}
|
|
inline void disable_fpexcept_sse() {
|
|
_MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() |
|
|
_MM_MASK_INVALID | _MM_MASK_DIV_ZERO);
|
|
}
|
|
} // namespace
|
|
|
|
#endif
|
|
|
|
MTS_NAMESPACE_BEGIN
|
|
|
|
// -----------------------------------------------------------------------
|
|
// General utility functions
|
|
// -----------------------------------------------------------------------
|
|
|
|
std::vector<std::string> tokenize(const std::string &string, const std::string &delim) {
|
|
std::string::size_type lastPos = string.find_first_not_of(delim, 0);
|
|
std::string::size_type pos = string.find_first_of(delim, lastPos);
|
|
std::vector<std::string> tokens;
|
|
|
|
while (std::string::npos != pos || std::string::npos != lastPos) {
|
|
tokens.push_back(string.substr(lastPos, pos - lastPos));
|
|
lastPos = string.find_first_not_of(delim, pos);
|
|
pos = string.find_first_of(delim, lastPos);
|
|
}
|
|
|
|
return tokens;
|
|
}
|
|
|
|
std::string trim(const std::string& str) {
|
|
std::string::size_type
|
|
start = str.find_first_not_of(" \t\r\n"),
|
|
end = str.find_last_not_of(" \t\r\n");
|
|
|
|
return str.substr(start == std::string::npos ? 0 : start,
|
|
end == std::string::npos ? str.length() - 1 : end - start + 1);
|
|
}
|
|
|
|
std::string indent(const std::string &string, int amount) {
|
|
/* This could probably be done faster (is not
|
|
really speed-critical though) */
|
|
std::istringstream iss(string);
|
|
std::ostringstream oss;
|
|
std::string str;
|
|
bool firstLine = true;
|
|
while (!iss.eof()) {
|
|
std::getline(iss, str);
|
|
if (!firstLine) {
|
|
for (int i=0; i<amount; ++i)
|
|
oss << " ";
|
|
}
|
|
oss << str;
|
|
if (!iss.eof())
|
|
oss << endl;
|
|
firstLine = false;
|
|
}
|
|
return oss.str();
|
|
}
|
|
|
|
std::string memString(size_t size) {
|
|
Float value = (Float) size;
|
|
const char *prefixes[] = {
|
|
"B", "KiB", "MiB", "GiB", "TiB", "PiB"
|
|
};
|
|
int prefix = 0;
|
|
while (prefix < 5 && value > 1024.0f) {
|
|
value /= 1024.0f; ++prefix;
|
|
}
|
|
return formatString(prefix == 0 ?
|
|
"%.0f %s" : "%.2f %s", value, prefixes[prefix]);
|
|
}
|
|
|
|
void * __restrict allocAligned(size_t size) {
|
|
#if defined(WIN32)
|
|
return _aligned_malloc(size, L1_CACHE_LINE_SIZE);
|
|
#elif defined(__OSX__)
|
|
/* OSX malloc already returns 16-byte aligned data suitable
|
|
for AltiVec and SSE computations */
|
|
return malloc(size);
|
|
#else
|
|
return memalign(L1_CACHE_LINE_SIZE, size);
|
|
#endif
|
|
}
|
|
|
|
void freeAligned(void *ptr) {
|
|
#if defined(WIN32)
|
|
_aligned_free(ptr);
|
|
#else
|
|
free(ptr);
|
|
#endif
|
|
}
|
|
|
|
int getCoreCount() {
|
|
#if defined(WIN32)
|
|
SYSTEM_INFO sys_info;
|
|
GetSystemInfo(&sys_info);
|
|
return sys_info.dwNumberOfProcessors;
|
|
#elif defined(__OSX__)
|
|
int nprocs;
|
|
size_t nprocsSize = sizeof(int);
|
|
if (sysctlbyname("hw.activecpu", &nprocs, &nprocsSize, NULL, 0))
|
|
SLog(EError, "Could not detect the number of processors!");
|
|
return (int) nprocs;
|
|
#else
|
|
return sysconf(_SC_NPROCESSORS_CONF);
|
|
#endif
|
|
}
|
|
|
|
#if defined(WIN32)
|
|
std::string lastErrorText() {
|
|
DWORD errCode = GetLastError();
|
|
char *errorText = NULL;
|
|
if (!FormatMessage(FORMAT_MESSAGE_ALLOCATE_BUFFER
|
|
| FORMAT_MESSAGE_FROM_SYSTEM
|
|
| FORMAT_MESSAGE_IGNORE_INSERTS,
|
|
NULL,
|
|
errCode,
|
|
MAKELANGID(LANG_NEUTRAL, SUBLANG_DEFAULT),
|
|
(LPTSTR) &errorText,
|
|
0,
|
|
NULL)) {
|
|
return "Internal error while looking up an error code";
|
|
}
|
|
std::string result(errorText);
|
|
LocalFree(errorText);
|
|
return result;
|
|
}
|
|
#endif
|
|
|
|
bool enableFPExceptions() {
|
|
bool exceptionsWereEnabled = false;
|
|
#if defined(WIN32)
|
|
_clearfp();
|
|
uint32_t cw = _controlfp(0, 0);
|
|
exceptionsWereEnabled = ~cw & (_EM_INVALID | _EM_ZERODIVIDE | _EM_OVERFLOW);
|
|
cw &= ~(_EM_INVALID | _EM_ZERODIVIDE | _EM_OVERFLOW);
|
|
_controlfp(cw, _MCW_EM);
|
|
#elif defined(__OSX__)
|
|
exceptionsWereEnabled = query_fpexcept_sse() != 0;
|
|
#else
|
|
exceptionsWereEnabled =
|
|
fegetexcept() & (FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW);
|
|
feenableexcept(FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW);
|
|
#endif
|
|
enable_fpexcept_sse();
|
|
return exceptionsWereEnabled;
|
|
}
|
|
|
|
bool disableFPExceptions() {
|
|
bool exceptionsWereEnabled = false;
|
|
#if defined(WIN32)
|
|
_clearfp();
|
|
uint32_t cw = _controlfp(0, 0);
|
|
exceptionsWereEnabled = ~cw & (_EM_INVALID | _EM_ZERODIVIDE | _EM_OVERFLOW);
|
|
cw |= _EM_INVALID | _EM_ZERODIVIDE | _EM_OVERFLOW;
|
|
_controlfp(cw, _MCW_EM);
|
|
#elif defined(__OSX__)
|
|
exceptionsWereEnabled = query_fpexcept_sse() != 0;
|
|
#else
|
|
exceptionsWereEnabled =
|
|
fegetexcept() & (FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW);
|
|
fedisableexcept(FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW);
|
|
#endif
|
|
disable_fpexcept_sse();
|
|
return exceptionsWereEnabled;
|
|
}
|
|
|
|
void restoreFPExceptions(bool oldState) {
|
|
bool currentState;
|
|
#if defined(WIN32)
|
|
uint32_t cw = _controlfp(0, 0);
|
|
currentState = ~cw & (_EM_INVALID | _EM_ZERODIVIDE | _EM_OVERFLOW);
|
|
#elif defined(__OSX__)
|
|
currentState = query_fpexcept_sse() != 0;
|
|
#else
|
|
currentState = fegetexcept() & (FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW);
|
|
#endif
|
|
if (oldState != currentState) {
|
|
if (oldState)
|
|
enableFPExceptions();
|
|
else
|
|
disableFPExceptions();
|
|
}
|
|
}
|
|
|
|
std::string getHostName() {
|
|
char hostName[128];
|
|
if (gethostname(hostName, sizeof(hostName)) != 0)
|
|
#if defined(WIN32)
|
|
SLog(EError, "Could not retrieve the computer's host name: %s!",
|
|
lastErrorText().c_str());
|
|
#else
|
|
SLog(EError, "Could not retrieve the computer's host name : %s!",
|
|
strerror(errno));
|
|
#endif
|
|
return hostName;
|
|
}
|
|
|
|
std::string getFQDN() {
|
|
struct addrinfo *addrInfo = NULL, hints;
|
|
memset(&hints, 0, sizeof(addrinfo));
|
|
// Only look for IPv4 addresses -- perhaps revisit this later
|
|
hints.ai_family = AF_INET; // AF_UNSPEC
|
|
hints.ai_socktype = SOCK_STREAM;
|
|
hints.ai_protocol = IPPROTO_TCP;
|
|
|
|
int retVal = getaddrinfo(getHostName().c_str(), NULL, &hints, &addrInfo);
|
|
if (addrInfo == NULL || retVal != 0) {
|
|
SLog(EWarn, "Could not retrieve the computer's fully "
|
|
"qualified domain name: could not resolve host address \"%s\"!",
|
|
getHostName().c_str());
|
|
return getHostName();
|
|
}
|
|
|
|
char fqdn[NI_MAXHOST];
|
|
retVal = getnameinfo(addrInfo->ai_addr, sizeof(struct sockaddr_in),
|
|
fqdn, NI_MAXHOST, NULL, 0, 0);
|
|
if (retVal != 0) {
|
|
freeaddrinfo(addrInfo);
|
|
#if defined(WIN32)
|
|
SLog(EWarn, "Could not retrieve the computer's fully "
|
|
"qualified domain name: error %i!", WSAGetLastError());
|
|
#else
|
|
SLog(EWarn, "Could not retrieve the computer's fully "
|
|
"qualified domain name: error %i!", gai_strerror(retVal));
|
|
#endif
|
|
return getHostName();
|
|
}
|
|
|
|
freeaddrinfo(addrInfo);
|
|
|
|
return fqdn;
|
|
}
|
|
|
|
Float log2(Float value) {
|
|
const Float invLn2 = (Float) 1.0f / math::fastlog((Float) 2.0f);
|
|
return math::fastlog(value) * invLn2;
|
|
}
|
|
|
|
std::string formatString(const char *fmt, ...) {
|
|
char tmp[512];
|
|
va_list iterator;
|
|
|
|
#if defined(WIN32)
|
|
va_start(iterator, fmt);
|
|
size_t size = _vscprintf(fmt, iterator) + 1;
|
|
|
|
if (size >= sizeof(tmp)) {
|
|
char *dest = new char[size];
|
|
vsnprintf_s(dest, size, size-1, fmt, iterator);
|
|
va_end(iterator);
|
|
std::string result(dest);
|
|
delete[] dest;
|
|
return result;
|
|
}
|
|
|
|
vsnprintf_s(tmp, size, size-1, fmt, iterator);
|
|
va_end(iterator);
|
|
#else
|
|
va_start(iterator, fmt);
|
|
size_t size = vsnprintf(tmp, sizeof(tmp), fmt, iterator);
|
|
va_end(iterator);
|
|
|
|
if (size >= sizeof(tmp)) {
|
|
/* Overflow! -- dynamically allocate memory */
|
|
char *dest = new char[size+1];
|
|
va_start(iterator, fmt);
|
|
vsnprintf(dest, size+1, fmt, iterator);
|
|
va_end(iterator);
|
|
|
|
std::string result(dest);
|
|
delete[] dest;
|
|
return result;
|
|
}
|
|
#endif
|
|
|
|
return std::string(tmp);
|
|
}
|
|
|
|
int log2i(uint32_t value) {
|
|
int r = 0;
|
|
while ((value >> r) != 0)
|
|
r++;
|
|
return r-1;
|
|
}
|
|
|
|
int log2i(uint64_t value) {
|
|
int r = 0;
|
|
while ((value >> r) != 0)
|
|
r++;
|
|
return r-1;
|
|
}
|
|
|
|
/* Fast rounding & power-of-two test algorithms from PBRT */
|
|
uint32_t roundToPowerOfTwo(uint32_t i) {
|
|
i--;
|
|
i |= i >> 1; i |= i >> 2;
|
|
i |= i >> 4; i |= i >> 8;
|
|
i |= i >> 16;
|
|
return i+1;
|
|
}
|
|
|
|
uint64_t roundToPowerOfTwo(uint64_t i) {
|
|
i--;
|
|
i |= i >> 1; i |= i >> 2;
|
|
i |= i >> 4; i |= i >> 8;
|
|
i |= i >> 16; i |= i >> 32;
|
|
return i+1;
|
|
}
|
|
|
|
// -----------------------------------------------------------------------
|
|
// Numerical utility functions
|
|
// -----------------------------------------------------------------------
|
|
|
|
bool solveQuadratic(Float a, Float b, Float c, Float &x0, Float &x1) {
|
|
/* Linear case */
|
|
if (a == 0) {
|
|
if (b != 0) {
|
|
x0 = x1 = -c / b;
|
|
return true;
|
|
}
|
|
return false;
|
|
}
|
|
|
|
Float discrim = b*b - 4.0f*a*c;
|
|
|
|
/* Leave if there is no solution */
|
|
if (discrim < 0)
|
|
return false;
|
|
|
|
Float temp, sqrtDiscrim = std::sqrt(discrim);
|
|
|
|
/* Numerically stable version of (-b (+/-) sqrtDiscrim) / (2 * a)
|
|
*
|
|
* Based on the observation that one solution is always
|
|
* accurate while the other is not. Finds the solution of
|
|
* greater magnitude which does not suffer from loss of
|
|
* precision and then uses the identity x1 * x2 = c / a
|
|
*/
|
|
if (b < 0)
|
|
temp = -0.5f * (b - sqrtDiscrim);
|
|
else
|
|
temp = -0.5f * (b + sqrtDiscrim);
|
|
|
|
x0 = temp / a;
|
|
x1 = c / temp;
|
|
|
|
/* Return the results so that x0 < x1 */
|
|
if (x0 > x1)
|
|
std::swap(x0, x1);
|
|
|
|
return true;
|
|
}
|
|
|
|
bool solveQuadraticDouble(double a, double b, double c, double &x0, double &x1) {
|
|
/* Linear case */
|
|
if (a == 0) {
|
|
if (b != 0) {
|
|
x0 = x1 = -c / b;
|
|
return true;
|
|
}
|
|
return false;
|
|
}
|
|
|
|
double discrim = b*b - 4.0f*a*c;
|
|
|
|
/* Leave if there is no solution */
|
|
if (discrim < 0)
|
|
return false;
|
|
|
|
double temp, sqrtDiscrim = std::sqrt(discrim);
|
|
|
|
/* Numerically stable version of (-b (+/-) sqrtDiscrim) / (2 * a)
|
|
*
|
|
* Based on the observation that one solution is always
|
|
* accurate while the other is not. Finds the solution of
|
|
* greater magnitude which does not suffer from loss of
|
|
* precision and then uses the identity x1 * x2 = c / a
|
|
*/
|
|
if (b < 0)
|
|
temp = -0.5f * (b - sqrtDiscrim);
|
|
else
|
|
temp = -0.5f * (b + sqrtDiscrim);
|
|
|
|
x0 = temp / a;
|
|
x1 = c / temp;
|
|
|
|
/* Return the results so that x0 < x1 */
|
|
if (x0 > x1)
|
|
std::swap(x0, x1);
|
|
|
|
return true;
|
|
}
|
|
|
|
bool solveLinearSystem2x2(const Float a[2][2], const Float b[2], Float x[2]) {
|
|
Float det = a[0][0] * a[1][1] - a[0][1] * a[1][0];
|
|
|
|
if (det == 0)
|
|
return false;
|
|
|
|
Float inverse = (Float) 1.0f / det;
|
|
|
|
x[0] = (a[1][1] * b[0] - a[0][1] * b[1]) * inverse;
|
|
x[1] = (a[0][0] * b[1] - a[1][0] * b[0]) * inverse;
|
|
|
|
return true;
|
|
}
|
|
|
|
void stratifiedSample1D(Random *random, Float *dest, int count, bool jitter) {
|
|
Float invCount = 1.0f / count;
|
|
|
|
for (int i=0; i<count; i++) {
|
|
Float offset = jitter ? random->nextFloat() : 0.5f;
|
|
*dest++ = (i + offset) * invCount;
|
|
}
|
|
}
|
|
|
|
void stratifiedSample2D(Random *random, Point2 *dest, int countX, int countY, bool jitter) {
|
|
Float invCountX = 1.0f / countX;
|
|
Float invCountY = 1.0f / countY;
|
|
|
|
for (int x=0; x<countX; x++) {
|
|
for (int y=0; y<countY; y++) {
|
|
Float offsetX = jitter ? random->nextFloat() : 0.5f;
|
|
Float offsetY = jitter ? random->nextFloat() : 0.5f;
|
|
*dest++ = Point2(
|
|
(x + offsetX) * invCountX,
|
|
(y + offsetY) * invCountY
|
|
);
|
|
}
|
|
}
|
|
}
|
|
|
|
void latinHypercube(Random *random, Float *dest, size_t nSamples, size_t nDim) {
|
|
Float delta = 1 / (Float) nSamples;
|
|
for (size_t i = 0; i < nSamples; ++i)
|
|
for (size_t j = 0; j < nDim; ++j)
|
|
dest[nDim * i + j] = (i + random->nextFloat()) * delta;
|
|
for (size_t i = 0; i < nDim; ++i) {
|
|
for (size_t j = 0; j < nSamples; ++j) {
|
|
size_t other = random->nextSize(nSamples);
|
|
std::swap(dest[nDim * j + i], dest[nDim * other + i]);
|
|
}
|
|
}
|
|
}
|
|
|
|
Vector sphericalDirection(Float theta, Float phi) {
|
|
Float sinTheta, cosTheta, sinPhi, cosPhi;
|
|
|
|
math::sincos(theta, &sinTheta, &cosTheta);
|
|
math::sincos(phi, &sinPhi, &cosPhi);
|
|
|
|
return Vector(
|
|
sinTheta * cosPhi,
|
|
sinTheta * sinPhi,
|
|
cosTheta
|
|
);
|
|
}
|
|
|
|
void coordinateSystem(const Vector &a, Vector &b, Vector &c) {
|
|
if (std::abs(a.x) > std::abs(a.y)) {
|
|
Float invLen = 1.0f / std::sqrt(a.x * a.x + a.z * a.z);
|
|
c = Vector(a.z * invLen, 0.0f, -a.x * invLen);
|
|
} else {
|
|
Float invLen = 1.0f / std::sqrt(a.y * a.y + a.z * a.z);
|
|
c = Vector(0.0f, a.z * invLen, -a.y * invLen);
|
|
}
|
|
b = cross(c, a);
|
|
}
|
|
|
|
Point2 toSphericalCoordinates(const Vector &v) {
|
|
Point2 result(
|
|
std::acos(v.z),
|
|
std::atan2(v.y, v.x)
|
|
);
|
|
if (result.y < 0)
|
|
result.y += 2*M_PI;
|
|
return result;
|
|
}
|
|
|
|
Float fresnelDielectric(Float cosThetaI, Float cosThetaT, Float eta) {
|
|
if (EXPECT_NOT_TAKEN(eta == 1))
|
|
return 0.0f;
|
|
|
|
Float Rs = (cosThetaI - eta * cosThetaT)
|
|
/ (cosThetaI + eta * cosThetaT);
|
|
Float Rp = (eta * cosThetaI - cosThetaT)
|
|
/ (eta * cosThetaI + cosThetaT);
|
|
|
|
/* No polarization -- return the unpolarized reflectance */
|
|
return 0.5f * (Rs * Rs + Rp * Rp);
|
|
}
|
|
|
|
Float fresnelDielectricExt(Float cosThetaI_, Float &cosThetaT_, Float eta) {
|
|
if (EXPECT_NOT_TAKEN(eta == 1)) {
|
|
cosThetaT_ = -cosThetaI_;
|
|
return 0.0f;
|
|
}
|
|
|
|
/* Using Snell's law, calculate the squared sine of the
|
|
angle between the normal and the transmitted ray */
|
|
Float scale = (cosThetaI_ > 0) ? 1/eta : eta,
|
|
cosThetaTSqr = 1 - (1-cosThetaI_*cosThetaI_) * (scale*scale);
|
|
|
|
/* Check for total internal reflection */
|
|
if (cosThetaTSqr <= 0.0f) {
|
|
cosThetaT_ = 0.0f;
|
|
return 1.0f;
|
|
}
|
|
|
|
/* Find the absolute cosines of the incident/transmitted rays */
|
|
Float cosThetaI = std::abs(cosThetaI_);
|
|
Float cosThetaT = std::sqrt(cosThetaTSqr);
|
|
|
|
Float Rs = (cosThetaI - eta * cosThetaT)
|
|
/ (cosThetaI + eta * cosThetaT);
|
|
Float Rp = (eta * cosThetaI - cosThetaT)
|
|
/ (eta * cosThetaI + cosThetaT);
|
|
|
|
cosThetaT_ = (cosThetaI_ > 0) ? -cosThetaT : cosThetaT;
|
|
|
|
/* No polarization -- return the unpolarized reflectance */
|
|
return 0.5f * (Rs * Rs + Rp * Rp);
|
|
}
|
|
|
|
Float fresnelConductorApprox(Float cosThetaI, Float eta, Float k) {
|
|
Float cosThetaI2 = cosThetaI*cosThetaI;
|
|
|
|
Float tmp = (eta*eta + k*k) * cosThetaI2;
|
|
|
|
Float Rp2 = (tmp - (eta * (2 * cosThetaI)) + 1)
|
|
/ (tmp + (eta * (2 * cosThetaI)) + 1);
|
|
|
|
Float tmpF = eta*eta + k*k;
|
|
|
|
Float Rs2 = (tmpF - (eta * (2 * cosThetaI)) + cosThetaI2) /
|
|
(tmpF + (eta * (2 * cosThetaI)) + cosThetaI2);
|
|
|
|
return 0.5f * (Rp2 + Rs2);
|
|
}
|
|
|
|
Spectrum fresnelConductorApprox(Float cosThetaI, const Spectrum &eta, const Spectrum &k) {
|
|
Float cosThetaI2 = cosThetaI*cosThetaI;
|
|
|
|
Spectrum tmp = (eta*eta + k*k) * cosThetaI2;
|
|
|
|
Spectrum Rp2 = (tmp - (eta * (2 * cosThetaI)) + Spectrum(1.0f))
|
|
/ (tmp + (eta * (2 * cosThetaI)) + Spectrum(1.0f));
|
|
|
|
Spectrum tmpF = eta*eta + k*k;
|
|
|
|
Spectrum Rs2 = (tmpF - (eta * (2 * cosThetaI)) + Spectrum(cosThetaI2)) /
|
|
(tmpF + (eta * (2 * cosThetaI)) + Spectrum(cosThetaI2));
|
|
|
|
return 0.5f * (Rp2 + Rs2);
|
|
}
|
|
|
|
Float fresnelConductorExact(Float cosThetaI, Float eta, Float k) {
|
|
/* Modified from "Optics" by K.D. Moeller, University Science Books, 1988 */
|
|
|
|
Float cosThetaI2 = cosThetaI*cosThetaI,
|
|
sinThetaI2 = 1-cosThetaI2,
|
|
sinThetaI4 = sinThetaI2*sinThetaI2;
|
|
|
|
Float temp1 = eta*eta - k*k - sinThetaI2,
|
|
a2pb2 = math::safe_sqrt(temp1*temp1 + 4*k*k*eta*eta),
|
|
a = math::safe_sqrt(0.5f * (a2pb2 + temp1));
|
|
|
|
Float term1 = a2pb2 + cosThetaI2,
|
|
term2 = 2*a*cosThetaI;
|
|
|
|
Float Rs2 = (term1 - term2) / (term1 + term2);
|
|
|
|
Float term3 = a2pb2*cosThetaI2 + sinThetaI4,
|
|
term4 = term2*sinThetaI2;
|
|
|
|
Float Rp2 = Rs2 * (term3 - term4) / (term3 + term4);
|
|
|
|
return 0.5f * (Rp2 + Rs2);
|
|
}
|
|
|
|
Spectrum fresnelConductorExact(Float cosThetaI, const Spectrum &eta, const Spectrum &k) {
|
|
/* Modified from "Optics" by K.D. Moeller, University Science Books, 1988 */
|
|
|
|
Float cosThetaI2 = cosThetaI*cosThetaI,
|
|
sinThetaI2 = 1-cosThetaI2,
|
|
sinThetaI4 = sinThetaI2*sinThetaI2;
|
|
|
|
Spectrum temp1 = eta*eta - k*k - Spectrum(sinThetaI2),
|
|
a2pb2 = (temp1*temp1 + k*k*eta*eta*4).safe_sqrt(),
|
|
a = ((a2pb2 + temp1) * 0.5f).safe_sqrt();
|
|
|
|
Spectrum term1 = a2pb2 + Spectrum(cosThetaI2),
|
|
term2 = a*(2*cosThetaI);
|
|
|
|
Spectrum Rs2 = (term1 - term2) / (term1 + term2);
|
|
|
|
Spectrum term3 = a2pb2*cosThetaI2 + Spectrum(sinThetaI4),
|
|
term4 = term2*sinThetaI2;
|
|
|
|
Spectrum Rp2 = Rs2 * (term3 - term4) / (term3 + term4);
|
|
|
|
return 0.5f * (Rp2 + Rs2);
|
|
}
|
|
|
|
Vector reflect(const Vector &wi, const Normal &n) {
|
|
return 2 * dot(wi, n) * Vector(n) - wi;
|
|
}
|
|
|
|
Vector refract(const Vector &wi, const Normal &n, Float eta, Float cosThetaT) {
|
|
if (cosThetaT < 0)
|
|
eta = 1 / eta;
|
|
|
|
return n * (dot(wi, n) * eta + cosThetaT) - wi * eta;
|
|
}
|
|
|
|
Vector refract(const Vector &wi, const Normal &n, Float eta) {
|
|
if (EXPECT_NOT_TAKEN(eta == 1))
|
|
return -wi;
|
|
|
|
Float cosThetaI = dot(wi, n);
|
|
if (cosThetaI > 0)
|
|
eta = 1 / eta;
|
|
|
|
/* Using Snell's law, calculate the squared sine of the
|
|
angle between the normal and the transmitted ray */
|
|
Float cosThetaTSqr = 1 - (1-cosThetaI*cosThetaI) * (eta*eta);
|
|
|
|
/* Check for total internal reflection */
|
|
if (cosThetaTSqr <= 0.0f)
|
|
return Vector(0.0f);
|
|
|
|
return n * (cosThetaI * eta - math::signum(cosThetaI)
|
|
* std::sqrt(cosThetaTSqr)) - wi * eta;
|
|
}
|
|
|
|
Vector refract(const Vector &wi, const Normal &n, Float eta, Float &cosThetaT, Float &F) {
|
|
Float cosThetaI = dot(wi, n);
|
|
F = fresnelDielectricExt(cosThetaI, cosThetaT, eta);
|
|
|
|
if (F == 1.0f) /* Total internal reflection */
|
|
return Vector(0.0f);
|
|
|
|
if (cosThetaT < 0)
|
|
eta = 1 / eta;
|
|
|
|
return n * (eta * cosThetaI + cosThetaT) - wi * eta;
|
|
}
|
|
|
|
namespace {
|
|
/// Integrand used by fresnelDiffuseReflectance
|
|
inline Float fresnelDiffuseIntegrand(Float eta, Float xi) {
|
|
return fresnelDielectricExt(std::sqrt(xi), eta);
|
|
}
|
|
};
|
|
|
|
Float fresnelDiffuseReflectance(Float eta, bool fast) {
|
|
if (fast) {
|
|
/* Fast mode: the following code approximates the
|
|
* diffuse Frensel reflectance for the eta<1 and
|
|
* eta>1 cases. An evalution of the accuracy led
|
|
* to the following scheme, which cherry-picks
|
|
* fits from two papers where they are best.
|
|
*/
|
|
if (eta < 1) {
|
|
/* Fit by Egan and Hilgeman (1973). Works
|
|
reasonably well for "normal" IOR values (<2).
|
|
|
|
Max rel. error in 1.0 - 1.5 : 0.1%
|
|
Max rel. error in 1.5 - 2 : 0.6%
|
|
Max rel. error in 2.0 - 5 : 9.5%
|
|
*/
|
|
return -1.4399f * (eta * eta)
|
|
+ 0.7099f * eta
|
|
+ 0.6681f
|
|
+ 0.0636f / eta;
|
|
} else {
|
|
/* Fit by d'Eon and Irving (2011)
|
|
*
|
|
* Maintains a good accuracy even for
|
|
* unrealistic IOR values.
|
|
*
|
|
* Max rel. error in 1.0 - 2.0 : 0.1%
|
|
* Max rel. error in 2.0 - 10.0 : 0.2%
|
|
*/
|
|
Float invEta = 1.0f / eta,
|
|
invEta2 = invEta*invEta,
|
|
invEta3 = invEta2*invEta,
|
|
invEta4 = invEta3*invEta,
|
|
invEta5 = invEta4*invEta;
|
|
|
|
return 0.919317f - 3.4793f * invEta
|
|
+ 6.75335f * invEta2
|
|
- 7.80989f * invEta3
|
|
+ 4.98554f * invEta4
|
|
- 1.36881f * invEta5;
|
|
}
|
|
} else {
|
|
GaussLobattoIntegrator quad(1024, 0, 1e-5f);
|
|
return quad.integrate(
|
|
boost::bind(&fresnelDiffuseIntegrand, eta, _1), 0, 1);
|
|
}
|
|
|
|
return 0.0f;
|
|
}
|
|
|
|
std::string timeString(Float time, bool precise) {
|
|
std::ostringstream os;
|
|
|
|
if (std::isnan(time) || std::isinf(time))
|
|
return "inf";
|
|
|
|
char suffix = 's';
|
|
if (time > 60) {
|
|
time /= 60; suffix = 'm';
|
|
if (time > 60) {
|
|
time /= 60; suffix = 'h';
|
|
if (time > 12) {
|
|
time /= 12; suffix = 'd';
|
|
}
|
|
}
|
|
}
|
|
|
|
os << std::setprecision(precise ? 4 : 1)
|
|
<< std::fixed << time << suffix;
|
|
|
|
return os.str();
|
|
}
|
|
|
|
Float hypot2(Float a, Float b) {
|
|
Float r;
|
|
if (std::abs(a) > std::abs(b)) {
|
|
r = b / a;
|
|
r = std::abs(a) * std::sqrt(1 + r*r);
|
|
} else if (b != 0) {
|
|
r = a / b;
|
|
r = std::abs(b) * std::sqrt(1 + r*r);
|
|
} else {
|
|
r = 0;
|
|
}
|
|
return r;
|
|
}
|
|
|
|
MTS_NAMESPACE_END
|