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

with the fixed ray epsilon issue, sphere intersections can be switched back to single precision

metadata
Wenzel Jakob 2010-11-03 22:14:26 +01:00
parent cce9775c0d
commit 41efe48e01
3 changed files with 14 additions and 80 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -96,12 +96,6 @@ extern MTS_EXPORT_CORE Float lanczosSinc(Float t, Float tau = 2);
*/
extern MTS_EXPORT_CORE bool solveQuadratic(Float a, Float b, Float c, Float &x0, Float &x1);
/**
* \brief Similar to solveQuadratic(), but always uses double precision independent
* of the chosen compile-time precision.
*/
extern MTS_EXPORT_CORE bool solveQuadraticDouble(double a, double b, double c, double &x0, double &x1);
/**
* \brief Solve a 2x2 linear equation system using basic linear algebra
*/

40
src/libcore/util.cpp
View File

@ -387,46 +387,6 @@ bool solveQuadratic(Float a, Float b, Float c, Float &x0, Float &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];

48
src/shapes/sphere.cpp
View File

@ -87,23 +87,13 @@ public:
}
bool rayIntersect(const Ray &ray, Float mint, Float maxt, Float &t, void *tmp) const {
/* 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;
Vector o = ray.o - m_center;
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*o.x + ray.d.y*o.y + ray.d.z*o.z);
Float C = o.x*o.x + o.y*o.y + o.z*o.z - m_radius*m_radius;
/* Transform into the local coordinate system and normalize */
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))
Float nearT, farT;
if (!solveQuadratic(A, B, C, nearT, farT))
return false;
if (nearT > maxt || farT < mint)
@ -111,32 +101,22 @@ public:
if (nearT < mint) {
if (farT > maxt)
return false;
t = (Float) farT;
t = farT;
} else {
t = (Float) nearT;
t = nearT;
}
return true;
}
bool rayIntersect(const Ray &ray, Float mint, Float maxt) const {
/* Do the following in double precision. This helps to avoid
self-intersections when approximating planes using giant spheres */
Vector o = ray.o - m_center;
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*o.x + ray.d.y*o.y + ray.d.z*o.z);
Float C = o.x*o.x + o.y*o.y + o.z*o.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;
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))
Float nearT, farT;
if (!solveQuadratic(A, B, C, nearT, farT))
return false;
if (nearT > maxt || farT < mint)