280 lines
9.5 KiB
C++
280 lines
9.5 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/core/statistics.h>
|
|
#include <mitsuba/render/particleproc.h>
|
|
#include <mitsuba/render/medium.h>
|
|
#include <mitsuba/render/phase.h>
|
|
#include <mitsuba/render/range.h>
|
|
|
|
MTS_NAMESPACE_BEGIN
|
|
|
|
ParticleProcess::ParticleProcess(EMode mode, size_t workCount, size_t granularity,
|
|
const std::string &progressText, const void *progressReporterPayload)
|
|
: m_mode(mode), m_workCount(workCount), m_numGenerated(0),
|
|
m_granularity(granularity), m_receivedResultCount(0) {
|
|
|
|
/* Choose a suitable work unit granularity if none was specified */
|
|
if (m_granularity == 0)
|
|
m_granularity = std::max((size_t) 1, workCount /
|
|
(16 * Scheduler::getInstance()->getWorkerCount()));
|
|
|
|
/* Create a visual progress reporter */
|
|
m_progress = new ProgressReporter(progressText, workCount,
|
|
progressReporterPayload);
|
|
m_resultMutex = new Mutex();
|
|
}
|
|
|
|
ParticleProcess::~ParticleProcess() {
|
|
delete m_progress;
|
|
}
|
|
|
|
ParallelProcess::EStatus ParticleProcess::generateWork(WorkUnit *unit, int worker) {
|
|
RangeWorkUnit *range = static_cast<RangeWorkUnit *>(unit);
|
|
size_t workUnitSize;
|
|
|
|
if (m_mode == ETrace) {
|
|
if (m_numGenerated == m_workCount)
|
|
return EFailure; // There is no more work
|
|
|
|
workUnitSize = std::min(m_granularity, m_workCount - m_numGenerated);
|
|
} else {
|
|
if (m_receivedResultCount >= m_workCount)
|
|
return EFailure; // There is no more work
|
|
|
|
workUnitSize = m_granularity;
|
|
}
|
|
|
|
range->setRange(m_numGenerated, m_numGenerated + workUnitSize - 1);
|
|
m_numGenerated += workUnitSize;
|
|
|
|
return ESuccess;
|
|
}
|
|
|
|
void ParticleProcess::increaseResultCount(size_t resultCount) {
|
|
LockGuard lock(m_resultMutex);
|
|
m_receivedResultCount += resultCount;
|
|
m_progress->update(m_receivedResultCount);
|
|
}
|
|
|
|
ParticleTracer::ParticleTracer(int maxDepth, int rrDepth, bool emissionEvents)
|
|
: m_maxDepth(maxDepth), m_rrDepth(rrDepth), m_emissionEvents(emissionEvents) { }
|
|
|
|
ParticleTracer::ParticleTracer(Stream *stream, InstanceManager *manager)
|
|
: WorkProcessor(stream, manager) {
|
|
|
|
m_maxDepth = stream->readInt();
|
|
m_rrDepth = stream->readInt();
|
|
m_emissionEvents = stream->readBool();
|
|
}
|
|
|
|
void ParticleTracer::serialize(Stream *stream, InstanceManager *manager) const {
|
|
stream->writeInt(m_maxDepth);
|
|
stream->writeInt(m_rrDepth);
|
|
stream->writeBool(m_emissionEvents);
|
|
}
|
|
|
|
ref<WorkUnit> ParticleTracer::createWorkUnit() const {
|
|
return new RangeWorkUnit();
|
|
}
|
|
|
|
void ParticleTracer::prepare() {
|
|
Scene *scene = static_cast<Scene *>(getResource("scene"));
|
|
m_scene = new Scene(scene);
|
|
m_sampler = static_cast<Sampler *>(getResource("sampler"));
|
|
Sensor *newSensor = static_cast<Sensor *>(getResource("sensor"));
|
|
m_scene->removeSensor(scene->getSensor());
|
|
m_scene->addSensor(newSensor);
|
|
m_scene->setSensor(newSensor);
|
|
m_scene->initializeBidirectional();
|
|
}
|
|
|
|
void ParticleTracer::process(const WorkUnit *workUnit, WorkResult *workResult,
|
|
const bool &stop) {
|
|
const RangeWorkUnit *range = static_cast<const RangeWorkUnit *>(workUnit);
|
|
MediumSamplingRecord mRec;
|
|
Intersection its;
|
|
ref<Sensor> sensor = m_scene->getSensor();
|
|
bool needsTimeSample = sensor->needsTimeSample();
|
|
PositionSamplingRecord pRec(sensor->getShutterOpen()
|
|
+ 0.5f * sensor->getShutterOpenTime());
|
|
Ray ray;
|
|
|
|
m_sampler->generate(Point2i(0));
|
|
|
|
for (size_t index = range->getRangeStart(); index <= range->getRangeEnd() && !stop; ++index) {
|
|
m_sampler->setSampleIndex(index);
|
|
|
|
/* Sample an emission */
|
|
if (needsTimeSample)
|
|
pRec.time = sensor->sampleTime(m_sampler->next1D());
|
|
|
|
const Emitter *emitter = NULL;
|
|
const Medium *medium;
|
|
|
|
Spectrum power;
|
|
Ray ray;
|
|
|
|
if (m_emissionEvents) {
|
|
/* Sample the position and direction component separately to
|
|
generate emission events */
|
|
power = m_scene->sampleEmitterPosition(pRec, m_sampler->next2D());
|
|
emitter = static_cast<const Emitter *>(pRec.object);
|
|
medium = emitter->getMedium();
|
|
|
|
/* Forward the sampling event to the attached handler */
|
|
handleEmission(pRec, medium, power);
|
|
|
|
DirectionSamplingRecord dRec;
|
|
power *= emitter->sampleDirection(dRec, pRec,
|
|
emitter->needsDirectionSample() ? m_sampler->next2D() : Point2(0.5f));
|
|
ray.setTime(pRec.time);
|
|
ray.setOrigin(pRec.p);
|
|
ray.setDirection(dRec.d);
|
|
} else {
|
|
/* Sample both components together, which is potentially
|
|
faster / uses a better sampling strategy */
|
|
|
|
power = m_scene->sampleEmitterRay(ray, emitter,
|
|
m_sampler->next2D(), m_sampler->next2D(), pRec.time);
|
|
medium = emitter->getMedium();
|
|
handleNewParticle();
|
|
}
|
|
|
|
int depth = 1;
|
|
bool delta = false;
|
|
|
|
Spectrum throughput(1.0f); // unitless path throughput (used for russian roulette)
|
|
while (!throughput.isZero() && (depth <= m_maxDepth || m_maxDepth < 0)) {
|
|
m_scene->rayIntersectAll(ray, its);
|
|
|
|
/* ==================================================================== */
|
|
/* Radiative Transfer Equation sampling */
|
|
/* ==================================================================== */
|
|
if (medium && medium->sampleDistance(Ray(ray, 0, its.t), mRec, m_sampler)) {
|
|
/* Sample the integral
|
|
\int_x^y tau(x, x') [ \sigma_s \int_{S^2} \rho(\omega,\omega') L(x,\omega') d\omega' ] dx'
|
|
*/
|
|
|
|
throughput *= mRec.sigmaS * mRec.transmittance / mRec.pdfSuccess;
|
|
|
|
/* Forward the medium scattering event to the attached handler */
|
|
handleMediumInteraction(depth,
|
|
delta, mRec, medium, -ray.d, throughput*power);
|
|
|
|
PhaseFunctionSamplingRecord pRec(mRec, -ray.d, EImportance);
|
|
|
|
throughput *= medium->getPhaseFunction()->sample(pRec, m_sampler);
|
|
delta = false;
|
|
|
|
ray = Ray(mRec.p, pRec.wo, ray.time);
|
|
ray.mint = 0;
|
|
} else if (its.t == std::numeric_limits<Float>::infinity()) {
|
|
/* There is no surface in this direction */
|
|
break;
|
|
} else {
|
|
/* Sample
|
|
tau(x, y) (Surface integral). This happens with probability mRec.pdfFailure
|
|
Account for this and multiply by the proper per-color-channel transmittance.
|
|
*/
|
|
if (medium)
|
|
throughput *= mRec.transmittance / mRec.pdfFailure;
|
|
|
|
const BSDF *bsdf = its.getBSDF();
|
|
|
|
/* Forward the surface scattering event to the attached handler */
|
|
handleSurfaceInteraction(depth, delta, its, medium, throughput*power);
|
|
|
|
BSDFSamplingRecord bRec(its, m_sampler, EImportance);
|
|
Spectrum bsdfWeight = bsdf->sample(bRec, m_sampler->next2D());
|
|
if (bsdfWeight.isZero())
|
|
break;
|
|
|
|
/* Prevent light leaks due to the use of shading normals -- [Veach, p. 158] */
|
|
Vector wi = -ray.d, wo = its.toWorld(bRec.wo);
|
|
Float wiDotGeoN = dot(its.geoFrame.n, wi),
|
|
woDotGeoN = dot(its.geoFrame.n, wo);
|
|
if (wiDotGeoN * Frame::cosTheta(bRec.wi) <= 0 ||
|
|
woDotGeoN * Frame::cosTheta(bRec.wo) <= 0)
|
|
break;
|
|
|
|
/* Keep track of the weight, medium and relative
|
|
refractive index along the path */
|
|
throughput *= bsdfWeight;
|
|
if (its.isMediumTransition())
|
|
medium = its.getTargetMedium(woDotGeoN);
|
|
|
|
#if 0
|
|
/* This is somewhat unfortunate: for accuracy, we'd really want the
|
|
correction factor below to match the path tracing interpretation
|
|
of a scene with shading normals. However, this factor can become
|
|
extremely large, which adds unacceptable variance to output
|
|
renderings.
|
|
|
|
So for now, it is disabled. The adjoint particle tracer and the
|
|
photon mapping variants still use this factor for the last
|
|
bounce -- just not for the intermediate ones, which introduces
|
|
a small (though in practice not noticeable) amount of error. This
|
|
is also what the implementation of SPPM by Toshiya Hachisuka does.
|
|
|
|
Ultimately, we'll need better adjoint BSDF sampling strategies
|
|
that incorporate these extra terms */
|
|
|
|
/* Adjoint BSDF for shading normals -- [Veach, p. 155] */
|
|
throughput *= std::abs(
|
|
(Frame::cosTheta(bRec.wi) * woDotGeoN)/
|
|
(Frame::cosTheta(bRec.wo) * wiDotGeoN));
|
|
#endif
|
|
|
|
ray.setOrigin(its.p);
|
|
ray.setDirection(wo);
|
|
ray.mint = Epsilon;
|
|
delta = (bRec.sampledType & (BSDF::EDelta & ~BSDF::ENull));
|
|
}
|
|
|
|
if (depth++ >= m_rrDepth) {
|
|
/* Russian roulette: try to keep path weights equal to one,
|
|
Stop with at least some probability to avoid
|
|
getting stuck (e.g. due to total internal reflection) */
|
|
|
|
Float q = std::min(throughput.max(), (Float) 0.95f);
|
|
if (m_sampler->next1D() >= q)
|
|
break;
|
|
throughput /= q;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
void ParticleTracer::handleEmission(const PositionSamplingRecord &pRec,
|
|
const Medium *medium, const Spectrum &weight) { }
|
|
|
|
void ParticleTracer::handleNewParticle() { }
|
|
|
|
void ParticleTracer::handleSurfaceInteraction(int depth, bool delta,
|
|
const Intersection &its, const Medium *medium, const Spectrum &weight) { }
|
|
|
|
void ParticleTracer::handleMediumInteraction(int depth, bool delta,
|
|
const MediumSamplingRecord &mRec, const Medium *medium,
|
|
const Vector &wi, const Spectrum &weight) { }
|
|
|
|
MTS_IMPLEMENT_CLASS(RangeWorkUnit, false, WorkUnit)
|
|
MTS_IMPLEMENT_CLASS(ParticleProcess, true, ParallelProcess)
|
|
MTS_IMPLEMENT_CLASS(ParticleTracer, true, WorkProcessor)
|
|
MTS_NAMESPACE_END
|