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

fixed some accuracy issues

metadata
Wenzel Jakob 2010-10-24 14:14:12 +02:00
parent c747d91263
commit 0013e696f7
6 changed files with 54 additions and 37 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
include/mitsuba/core/constants.h
View File

@ -22,17 +22,17 @@
/* Choice of precision */
#ifdef DOUBLE_PRECISION
#define Float double
#define Epsilon 1e-8
#define Epsilon 1e-6
#else
#ifndef SINGLE_PRECISION
#define SINGLE_PRECISION
#endif
#define Float float
#define Epsilon 1e-5f
#define Epsilon 1e-4f
#endif
/// relative eps. for shadow rays
#define ShadowEpsilon 1e-3
#define ShadowEpsilon 1e-3f
#ifdef M_PI
#undef M_PI

5
include/mitsuba/render/gkdtree.h
View File

@ -40,9 +40,6 @@
#define MTS_KD_BLOCKSIZE_KD (512*1024/sizeof(KDNode))
#define MTS_KD_BLOCKSIZE_IDX (512*1024/sizeof(uint32_t))
/// 32 byte temporary storage for intersection computations
#define MTS_KD_INTERSECTION_TEMP 32
/// Use a simple hashed 8-entry mailbox per thread
//#define MTS_KD_MAILBOX_ENABLED 1
//#define MTS_KD_MAILBOX_SIZE 8
@ -3117,7 +3114,7 @@ template<typename Derived> template <bool shadowRay> FINLINE bool
template <typename Derived> void GenericKDTree<Derived>::findCosts(
Float &traversalCost, Float &intersectionCost) {
ref<Random> random = new Random();
uint8_t temp[MTS_KD_INTERSECTION_TEMP];
uint8_t temp[128];
BSphere bsphere = m_aabb.getBSphere();
int nRays = 10000000, warmup = nRays/4;
Vector *A = new Vector[nRays-warmup];

7
include/mitsuba/render/kdtree.h
View File

@ -29,6 +29,13 @@
#endif
#endif
#if defined(SINGLE_PRECISION)
/// 32 byte temporary storage for intersection computations
#define MTS_KD_INTERSECTION_TEMP 32
#else
#define MTS_KD_INTERSECTION_TEMP 64
#endif
MTS_NAMESPACE_BEGIN
/**

2
src/librender/kdtree.cpp
View File

@ -112,7 +112,7 @@ bool KDTree::rayIntersect(const Ray &ray, Intersection &its) const {
if (m_aabb.rayIntersect(ray, mint, maxt)) {
/* Use an adaptive ray epsilon */
Float rayMinT = ray.mint;
if (rayMinT == Epsilon)
if (rayMinT == Epsilon)
rayMinT *= std::max(std::max(std::abs(ray.o.x),
std::abs(ray.o.y)), std::abs(ray.o.z));

23
src/librender/preview.cpp
View File

@ -65,7 +65,7 @@ void PreviewWorker::processIncoherent(const WorkUnit *workUnit, WorkResult *work
int pos = 0;
Intersection its;
Spectrum value, bsdfVal;
Vector toIts;
Vector toVPL;
Ray primary, secondary;
int numRays = 0;
@ -87,33 +87,32 @@ void PreviewWorker::processIncoherent(const WorkUnit *workUnit, WorkResult *work
else
value = Spectrum(0.0f);
toIts = its.p - m_vpl.its.p;
secondary = Ray(m_vpl.its.p, toIts, 0.001, 0.999);
toVPL = m_vpl.its.p - its.p;
secondary = Ray(its.p, toVPL, ShadowEpsilon, 1-ShadowEpsilon);
++numRays;
if (m_kdtree->rayIntersect(secondary)) {
block->setPixel(pos++, value);
continue;
}
Float length = toVPL.length();
toVPL/=length;
Float length = toIts.length();
toIts /= length;
BSDFQueryRecord rr(its, -its.toLocal(toIts));
BSDFQueryRecord rr(its, its.toLocal(toVPL));
rr.wi = normalize(rr.wi);
bsdfVal = its.shape->getBSDF()->fCos(rr);
length = std::max(length, m_minDist);
if (m_vpl.type == ESurfaceVPL) {
BSDFQueryRecord bRec(m_vpl.its, m_vpl.its.toLocal(toIts));
BSDFQueryRecord bRec(m_vpl.its, -m_vpl.its.toLocal(toVPL));
bRec.quantity = EImportance;
value += m_vpl.P * bsdfVal * m_vpl.its.shape->getBSDF()->fCos(bRec) / (length*length);
} else {
EmissionRecord eRec(m_vpl.luminaire,
ShapeSamplingRecord(m_vpl.its.p, m_vpl.its.shFrame.n), toIts);
ShapeSamplingRecord(m_vpl.its.p, m_vpl.its.shFrame.n), -toVPL);
eRec.type = EmissionRecord::EPreview;
value += m_vpl.P * bsdfVal * m_vpl.luminaire->f(eRec)
* ((m_vpl.luminaire->getType() & Luminaire::EOnSurface ?
dot(m_vpl.its.shFrame.n, toIts) : (Float) 1)
dot(m_vpl.its.shFrame.n, -toVPL) : (Float) 1)
/ (length*length));
}
block->setPixel(pos++, value);
@ -198,7 +197,7 @@ void PreviewWorker::processCoherent(const WorkUnit *workUnit, WorkResult *workRe
primRay4.o[0].ps = _mm_set1_ps(m_cameraO.x);
primRay4.o[1].ps = _mm_set1_ps(m_cameraO.y);
primRay4.o[2].ps = _mm_set1_ps(m_cameraO.z);
secItv4.mint.ps = _mm_set1_ps(0.001);
secItv4.mint.ps = _mm_set1_ps(ShadowEpsilon);
/* Work on 2x2 sub-blocks */
for (int y=sy; y<ey; y += 2, pos += width) {
@ -272,7 +271,7 @@ void PreviewWorker::processCoherent(const WorkUnit *workUnit, WorkResult *workRe
_mm_mul_ps(nSecD[0].ps, lumDir[0]),
_mm_mul_ps(nSecD[1].ps, lumDir[1])),
_mm_mul_ps(nSecD[2].ps, lumDir[2])));
secItv4.maxt.ps = _mm_set1_ps(0.999);
secItv4.maxt.ps = _mm_set1_ps(1-ShadowEpsilon);
/* Shading (scalar) --- this is way too much work and should be
rewritten to be smarter in special cases */

48
src/shapes/sphere.cpp
View File

@ -86,16 +86,23 @@ public:
}
bool rayIntersect(const Ray &ray, Float mint, Float maxt, Float &t, void *tmp) const {
const Vector ro = ray.o - m_center;
/* Do the following in double precision. This helps to avoid
self-intersections when approximating planes using giant spheres */
const double
rox = (double) (ray.o.x - m_center.x),
roy = (double) (ray.o.y - m_center.y),
roz = (double) (ray.o.z - m_center.z),
rdx = (double) ray.d.x,
rdy = (double) ray.d.y,
rdz = (double) ray.d.z;
/* Transform into the local coordinate system and normalize */
Float A = ray.d.x*ray.d.x + ray.d.y*ray.d.y + ray.d.z*ray.d.z;
Float B = 2 * (ray.d.x*ro.x + ray.d.y*ro.y + ray.d.z*ro.z);
Float C = ro.x*ro.x + ro.y*ro.y +
ro.z*ro.z - m_radius*m_radius;
double A = rdx*rdx + rdy*rdy + rdz*rdz;
double B = 2 * (rdx*rox + rdy*roy + rdz*roz);
double C = rox*rox + roy*roy + roz*roz - m_radius*m_radius;
Float nearT, farT;
if (!solveQuadratic(A, B, C, nearT, farT))
double nearT, farT;
if (!solveQuadraticDouble(A, B, C, nearT, farT))
return false;
if (nearT > maxt || farT < mint)
@ -103,25 +110,32 @@ public:
if (nearT < mint) {
if (farT > maxt)
return false;
t = farT;
t = (Float) farT;
} else {
t = nearT;
t = (Float) nearT;
}
return true;
}
bool rayIntersect(const Ray &ray, Float mint, Float maxt) const {
const Vector ro = ray.o - m_center;
/* Do the following in double precision. This helps to avoid
self-intersections when approximating planes using giant spheres */
/* Transform into the local coordinate system and normalize */
Float A = ray.d.x*ray.d.x + ray.d.y*ray.d.y + ray.d.z*ray.d.z;
Float B = 2 * (ray.d.x*ro.x + ray.d.y*ro.y + ray.d.z*ro.z);
Float C = ro.x*ro.x + ro.y*ro.y +
ro.z*ro.z - m_radius*m_radius;
const double
rox = (double) (ray.o.x - m_center.x),
roy = (double) (ray.o.y - m_center.y),
roz = (double) (ray.o.z - m_center.z),
rdx = (double) ray.d.x,
rdy = (double) ray.d.y,
rdz = (double) ray.d.z;
Float nearT, farT;
if (!solveQuadratic(A, B, C, nearT, farT))
double A = rdx*rdx + rdy*rdy + rdz*rdz;
double B = 2 * (rdx*rox + rdy*roy + rdz*roz);
double C = rox*rox + roy*roy + roz*roz - m_radius*m_radius;
double nearT, farT;
if (!solveQuadraticDouble(A, B, C, nearT, farT))
return false;
if (nearT > maxt || farT < mint)