mingzailao
/
mitsuba
mirror of https://github.com/Quartenia/mitsuba.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

cleaned up brent solver

metadata
Wenzel Jakob 2011-04-13 02:21:21 +02:00
parent 9baca24e5e
commit 9aa5c13daa
6 changed files with 343 additions and 167 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

189
include/mitsuba/core/brent.h
View File

@ -20,6 +20,7 @@
#define __BRENT_H
#include <mitsuba/mitsuba.h>
#include <boost/function.hpp>
MTS_NAMESPACE_BEGIN
@ -27,24 +28,22 @@ MTS_NAMESPACE_BEGIN
* \brief Brent's method nonlinear zero finder
*
* The implementation is transcribed from the Apache Commons
* Java implementation. The function \a Func is required to be
* Java implementation. The supplied function is required to be
* continuous, but not necessarily smooth.
*
* \ingroup libcore
*/
template <typename T, typename Func> class BrentSolver {
class BrentSolver {
public:
typedef T value_type;
/// Return value of \ref BrentSolver::solve()
struct Result {
bool success;
int iterations;
value_type x;
value_type y;
Float x;
Float y;
/// Create a new result instance
inline Result(bool success, int iterations, value_type x, value_type y)
inline Result(bool success, int iterations, Float x, Float y)
: success(success), iterations(iterations), x(x), y(y) { }
/// Return a string representation of the result
@ -75,10 +74,10 @@ public:
* Absolute accuracy requirement of the position -- the
* iterations will stop when |minX-maxX|/minX < relAccuracyPos.
*/
inline BrentSolver(int maxIterations = 100,
value_type absAccuracy = 1e-6f,
value_type absAccuracyPos = 1e-6f,
value_type relAccuracyPos = 1e-6f)
inline BrentSolver(size_t maxIterations = 100,
Float absAccuracy = 1e-6f,
Float absAccuracyPos = 1e-6f,
Float relAccuracyPos = 1e-6f)
: m_maxIterations(maxIterations),
m_absAccuracy(absAccuracy),
m_absAccuracyPos(absAccuracyPos),
@ -94,27 +93,8 @@ public:
* \param max the upper bound for the interval.
* \return the value where the function is zero
*/
Result solve(Func &f, value_type min, value_type max) const {
// return the first endpoint if it is good enough
value_type yMin = f(min);
if (std::abs(yMin) <= m_absAccuracy)
return Result(true, 0, min, yMin);
// return the second endpoint if it is good enough
value_type yMax = f(max);
if (std::abs(yMax) <= m_absAccuracy)
return Result(true, 0, max, yMax);
value_type sign = yMin * yMax;
if (sign > 0) {
SLog(EWarn, "BrentSolver: Function values at the endpoints do not have different signs -- "
"endpoints: [%f, %f], values: [%f, %f]", min, max, yMin, yMax);
return Result(false, 0, 0, 0);
} else {
// solve using only the first endpoint as initial guess
return solve(f, min, yMin, max, yMax, min, yMin);
}
}
Result solve(const boost::function<Float (Float)> &func,
Float min, Float max) const;
/**
* \brief Find a zero in the given interval with an initial guess
@ -130,41 +110,8 @@ public:
* if no initial point is known)
* \return the value where the function is zero
*/
Result solve(Func &f, value_type min, value_type max, value_type initial) const {
if (initial < min || initial > max) {
SLog(EWarn, "BrentSolver: Invalid interval: lower=%f, initial=%f, upper=%f",
min, max, initial);
return Result(false, 0, 0, 0);
}
// return the initial guess if it is good enough
value_type yInitial = f(initial);
if (std::abs(yInitial) <= m_absAccuracy)
return Result(true, 0, initial, yInitial);
// return the first endpoint if it is good enough
value_type yMin = f(min);
if (std::abs(yMin) <= m_absAccuracy)
return Result(true, 0, min, yMin);
// reduce interval if min and initial bracket the root
if (yInitial * yMin < 0)
return solve(f, min, yMin, initial, yInitial, min, yMin);
// return the second endpoint if it is good enough
value_type yMax = f(max);
if (std::abs(yMax) <= m_absAccuracy)
return Result(true, 0, max, yMax);
// reduce interval if initial and max bracket the root
if (yInitial * yMax < 0)
return solve(f, initial, yInitial, max, yMax, initial, yInitial);
SLog(EWarn, "BrentSolver: Function values at the endpoints do not have different signs -- "
"endpoints: [%f, %f], values: [%f, %f]", min, max, yMin, yMax);
return Result(false, 0, 0, 0);
}
Result solve(const boost::function<Float (Float)> &func,
Float min, Float max, Float initial) const;
/**
* Find a zero starting search according to the three provided points.
@ -179,105 +126,15 @@ public:
* \param y2 function value at the bracket point.
* \return the value where the function is zero
*/
Result solve(Func &f,
value_type x0, value_type y0,
value_type x1, value_type y1,
value_type x2, value_type y2) const {
value_type delta = x1 - x0;
value_type oldDelta = delta;
int i = 0;
while (i < m_maxIterations) {
if (std::abs(y2) < std::abs(y1)) {
// use the bracket point if is better than last approximation
x0 = x1;
x1 = x2;
x2 = x0;
y0 = y1;
y1 = y2;
y2 = y0;
}
if (std::abs(y1) <= m_absAccuracy) {
// Avoid division by very small values. Assume
// the iteration has converged (the problem may
// still be ill conditioned)
return Result(true, i, x1, y1);
}
value_type dx = x2 - x1;
value_type tolerance =
std::max(m_relAccuracyPos * std::abs(x1), m_absAccuracyPos);
if (std::abs(dx) <= tolerance)
return Result(true, i, x1, y1);
if ((std::abs(oldDelta) < tolerance) ||
(std::abs(y0) <= std::abs(y1))) {
// Force bisection.
delta = (value_type) 0.5f * dx;
oldDelta = delta;
} else {
value_type r3 = y1 / y0;
value_type p;
value_type p1;
// the equality test (x0 == x2) is intentional,
// it is part of the original Brent's method,
// it should NOT be replaced by proximity test
if (x0 == x2) {
// Linear interpolation.
p = dx * r3;
p1 = 1 - r3;
} else {
// Inverse quadratic interpolation.
value_type r1 = y0 / y2;
value_type r2 = y1 / y2;
p = r3 * (dx * r1 * (r1 - r2) - (x1 - x0) * (r2 - 1));
p1 = (r1 - 1) * (r2 - 1) * (r3 - 1);
}
if (p > 0) {
p1 = -p1;
} else {
p = -p;
}
if (2 * p >= (value_type) 1.5f * dx * p1 - std::abs(tolerance * p1) ||
p >= std::abs((value_type) 0.5f * oldDelta * p1)) {
// Inverse quadratic interpolation gives a value
// in the wrong direction, or progress is slow.
// Fall back to bisection.
delta = (value_type) 0.5f * dx;
oldDelta = delta;
} else {
oldDelta = delta;
delta = p / p1;
}
}
// Save old X1, Y1
x0 = x1;
y0 = y1;
// Compute new X1, Y1
if (std::abs(delta) > tolerance) {
x1 = x1 + delta;
} else if (dx > 0) {
x1 = x1 + (value_type) 0.5f * tolerance;
} else if (dx <= 0) {
x1 = x1 - (value_type) 0.5f * tolerance;
}
y1 = f(x1);
if ((y1 > 0) == (y2 > 0)) {
x2 = x0;
y2 = y0;
delta = x1 - x0;
oldDelta = delta;
}
i++;
}
SLog(EWarn, "BrentSolver: Maximum number of iterations (%i) exceeded!",
m_maxIterations);
return Result(false, i, x1, y1);
}
private:
int m_maxIterations;
value_type m_absAccuracy;
value_type m_absAccuracyPos;
value_type m_relAccuracyPos;
Result solve(const boost::function<Float (Float)> &func,
Float x0, Float y0,
Float x1, Float y1,
Float x2, Float y2) const;
protected:
size_t m_maxIterations;
Float m_absAccuracy;
Float m_absAccuracyPos;
Float m_relAccuracyPos;
};
MTS_NAMESPACE_END

2
src/libcore/SConscript
View File

@ -25,7 +25,7 @@ if coreEnv.has_key('JPEGLIB'):
coreEnv.Prepend(CPPDEFINES = [['MTS_BUILD_MODULE', 'MTS_MODULE_CORE']])
libcore_objects = [
'class.cpp', 'object.cpp', 'statistics.cpp', 'thread.cpp',
'class.cpp', 'object.cpp', 'statistics.cpp', 'thread.cpp', 'brent.cpp',
'logger.cpp', 'appender.cpp', 'formatter.cpp', 'lock.cpp',
'random.cpp', 'timer.cpp', 'util.cpp', 'properties.cpp',
'transform.cpp', 'spectrum.cpp', 'aabb.cpp', 'stream.cpp',

186
src/libcore/brent.cpp Normal file
View File

@ -0,0 +1,186 @@
/*
This file is part of Mitsuba, a physically based rendering system.
Copyright (c) 2007-2010 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/core/brent.h>
/**
* \file brent.cpp
* \brief Brent's method nonlinear zero finder
*
* The implementation is transcribed from the Apache Commons
* Java implementation.
*/
MTS_NAMESPACE_BEGIN
BrentSolver::Result BrentSolver::solve(const boost::function<Float (Float)> &f,
Float min, Float max) const {
// return the first endpoint if it is good enough
Float yMin = f(min);
if (std::abs(yMin) <= m_absAccuracy)
return Result(true, 0, min, yMin);
// return the second endpoint if it is good enough
Float yMax = f(max);
if (std::abs(yMax) <= m_absAccuracy)
return Result(true, 0, max, yMax);
Float sign = yMin * yMax;
if (sign > 0) {
SLog(EWarn, "BrentSolver: Function values at the endpoints do not have different signs -- "
"endpoints: [%f, %f], values: [%f, %f]", min, max, yMin, yMax);
return Result(false, 0, 0, 0);
} else {
// solve using only the first endpoint as initial guess
return solve(f, min, yMin, max, yMax, min, yMin);
}
}
BrentSolver::Result BrentSolver::solve(const boost::function<Float (Float)> &f,
Float min, Float max, Float initial) const {
if (initial < min || initial > max) {
SLog(EWarn, "BrentSolver: Invalid interval: lower=%f, initial=%f, upper=%f",
min, max, initial);
return Result(false, 0, 0, 0);
}
// return the initial guess if it is good enough
Float yInitial = f(initial);
if (std::abs(yInitial) <= m_absAccuracy)
return Result(true, 0, initial, yInitial);
// return the first endpoint if it is good enough
Float yMin = f(min);
if (std::abs(yMin) <= m_absAccuracy)
return Result(true, 0, min, yMin);
// reduce interval if min and initial bracket the root
if (yInitial * yMin < 0)
return solve(f, min, yMin, initial, yInitial, min, yMin);
// return the second endpoint if it is good enough
Float yMax = f(max);
if (std::abs(yMax) <= m_absAccuracy)
return Result(true, 0, max, yMax);
// reduce interval if initial and max bracket the root
if (yInitial * yMax < 0)
return solve(f, initial, yInitial, max, yMax, initial, yInitial);
SLog(EWarn, "BrentSolver: Function values at the endpoints do not have different signs -- "
"endpoints: [%f, %f], values: [%f, %f]", min, max, yMin, yMax);
return Result(false, 0, 0, 0);
}
BrentSolver::Result BrentSolver::solve(const boost::function<Float (Float)> &f,
Float x0, Float y0,
Float x1, Float y1,
Float x2, Float y2) const {
Float delta = x1 - x0;
Float oldDelta = delta;
size_t i = 0;
while (i < m_maxIterations) {
if (std::abs(y2) < std::abs(y1)) {
// use the bracket point if is better than last approximation
x0 = x1;
x1 = x2;
x2 = x0;
y0 = y1;
y1 = y2;
y2 = y0;
}
if (std::abs(y1) <= m_absAccuracy) {
// Avoid division by very small values. Assume
// the iteration has converged (the problem may
// still be ill conditioned)
return Result(true, i, x1, y1);
}
Float dx = x2 - x1;
Float tolerance =
std::max(m_relAccuracyPos * std::abs(x1), m_absAccuracyPos);
if (std::abs(dx) <= tolerance)
return Result(true, i, x1, y1);
if ((std::abs(oldDelta) < tolerance) ||
(std::abs(y0) <= std::abs(y1))) {
// Force bisection.
delta = (Float) 0.5f * dx;
oldDelta = delta;
} else {
Float r3 = y1 / y0;
Float p;
Float p1;
// the equality test (x0 == x2) is intentional,
// it is part of the original Brent's method,
// it should NOT be replaced by proximity test
if (x0 == x2) {
// Linear interpolation.
p = dx * r3;
p1 = 1 - r3;
} else {
// Inverse quadratic interpolation.
Float r1 = y0 / y2;
Float r2 = y1 / y2;
p = r3 * (dx * r1 * (r1 - r2) - (x1 - x0) * (r2 - 1));
p1 = (r1 - 1) * (r2 - 1) * (r3 - 1);
}
if (p > 0) {
p1 = -p1;
} else {
p = -p;
}
if (2 * p >= (Float) 1.5f * dx * p1 - std::abs(tolerance * p1) ||
p >= std::abs((Float) 0.5f * oldDelta * p1)) {
// Inverse quadratic interpolation gives a value
// in the wrong direction, or progress is slow.
// Fall back to bisection.
delta = (Float) 0.5f * dx;
oldDelta = delta;
} else {
oldDelta = delta;
delta = p / p1;
}
}
// Save old X1, Y1
x0 = x1;
y0 = y1;
// Compute new X1, Y1
if (std::abs(delta) > tolerance) {
x1 = x1 + delta;
} else if (dx > 0) {
x1 = x1 + (Float) 0.5f * tolerance;
} else if (dx <= 0) {
x1 = x1 - (Float) 0.5f * tolerance;
}
y1 = f(x1);
if ((y1 > 0) == (y2 > 0)) {
x2 = x0;
y2 = y0;
delta = x1 - x0;
oldDelta = delta;
}
i++;
}
SLog(EWarn, "BrentSolver: Maximum number of iterations (" SIZE_T_FMT ") exceeded!",
m_maxIterations);
return Result(false, i, x1, y1);
}
MTS_NAMESPACE_END

1
src/phase/SConscript
View File

@ -3,5 +3,6 @@ Import('env', 'plugins')
plugins += env.SharedLibrary('#plugins/isotropic', ['isotropic.cpp'])
plugins += env.SharedLibrary('#plugins/hg', ['hg.cpp'])
plugins += env.SharedLibrary('#plugins/kkay', ['kkay.cpp'])
plugins += env.SharedLibrary('#plugins/microflake', ['microflake.cpp'])
Export('plugins')

72
src/phase/microflake.cpp Normal file
View File

@ -0,0 +1,72 @@
/*
This file is part of Mitsuba, a physically based rendering system.
Copyright (c) 2007-2010 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/render/phase.h>
#include <mitsuba/render/sampler.h>
#include "microflake.h"
MTS_NAMESPACE_BEGIN
class MicroflakePhaseFunction : public PhaseFunction {
public:
MicroflakePhaseFunction(const Properties &props)
: PhaseFunction(props) {
}
MicroflakePhaseFunction(Stream *stream, InstanceManager *manager)
: PhaseFunction(stream, manager) {
}
virtual ~MicroflakePhaseFunction() { }
void serialize(Stream *stream, InstanceManager *manager) const {
PhaseFunction::serialize(stream, manager);
}
void configure() {
}
Spectrum sample(PhaseFunctionQueryRecord &pRec,
Sampler *sampler) const {
Point2 sample(sampler->next2D());
pRec.wo = squareToSphere(sample);
return Spectrum(1.0f);
}
Spectrum sample(PhaseFunctionQueryRecord &pRec,
Float &pdf, Sampler *sampler) const {
pRec.wo = squareToSphere(sampler->next2D());
pdf = 1/(4 * (Float) M_PI);
return Spectrum(pdf);
}
Spectrum f(const PhaseFunctionQueryRecord &pRec) const {
return Spectrum(1/(4 * (Float) M_PI));
}
std::string toString() const {
return "MicroflakePhaseFunction[]";
}
MTS_DECLARE_CLASS()
};
MTS_IMPLEMENT_CLASS_S(MicroflakePhaseFunction, false, PhaseFunction)
MTS_EXPORT_PLUGIN(MicroflakePhaseFunction, "Microflake phase function");
MTS_NAMESPACE_END

60
src/phase/microflake.h
View File

@ -19,6 +19,11 @@
#if !defined(__MICROFLAKE_SIGMA_T_H)
#define __MICROFLAKE_SIGMA_T_H
#include <mitsuba/core/statistics.h>
#include <mitsuba/core/brent.h>
#include <boost/math/special_functions/erf.hpp>
#include <boost/bind.hpp>
MTS_NAMESPACE_BEGIN
#define SIGMA_T_RANGE_MIN 4e-08
@ -166,6 +171,61 @@ double sigmaT_fiberGaussian(double stddev, double sinTheta) {
return result;
}
static StatsCounter brentSolves("Micro-flake model",
"Brent solver calls");
static StatsCounter avgBrentFunEvals("Micro-flake model",
"Average Brent solver function evaluations", EAverage);
struct GaussianFiberDistribution {
public:
inline GaussianFiberDistribution(Float stddev) : m_stddev(stddev) {
m_normalization = 1/(std::sqrt(2*M_PI) * m_stddev *
boost::math::erf<Float>(1/(SQRT_TWO * m_stddev)));
m_c1 = 1.0f/boost::math::erf<Float>(1/(SQRT_TWO * m_stddev));
}
/// Evaluate \sigma_t as a function of \cos\theta
Float sigmaT(Float cosTheta) const {
return sigmaT_fiberGaussian(m_stddev,
std::sqrt(std::max((Float) 0, 1-cosTheta*cosTheta)));
}
/// Evaluate the longitudinal density as a function of \cos\theta
Float pdfCosTheta(Float cosTheta) const {
return std::exp(-cosTheta*cosTheta/(2*m_stddev*m_stddev)) * m_normalization;
}
/// Evaluate the longitudinal CDF as a function of \cos\theta
inline Float cdf(Float theta) const {
return 0.5f * (1.0f -
boost::math::erf<Float>(std::cos(theta)/(SQRT_TWO * m_stddev))*m_c1);
}
/**
* Apply the inversion method to sample \theta given
* a uniformly distributed \xi on [0, 1]
*/
Float sample(Float xi) const {
BrentSolver brentSolver(100, 1e-6f);
BrentSolver::Result result = brentSolver.solve(
boost::bind(&GaussianFiberDistribution::cdfFunctor, this, xi, _1), 0, M_PI);
SAssert(result.success);
avgBrentFunEvals.incrementBase();
++brentSolves;
return result.x;
}
private:
Float cdfFunctor(Float theta, Float xi) const {
++avgBrentFunEvals;
return cdf(theta)-xi;
}
Float m_stddev;
Float m_normalization;
Float m_c1;
};
MTS_NAMESPACE_END
#endif /* __MICROFLAKE_SIGMA_T_H */