I'm building a software renderer but I think there might be some problem with my ray-triangle intersection accuracy.
I implemented the algorithm referring to pbrt-v3
without some components which I think have little influence on the intersection test.
Below is the code perform intersection test:
bool Triangle::intersectP(const Ray &r, Point3f p0, Point3f p1, Point3f p2) const {
// translate
Point3f p0_ = p0 - r._o;
Point3f p1_ = p1 - r._o;
Point3f p2_ = p2 - r._o;
// permute, to make sure z-dimension is not zero
// permutation is just like to circular move right/left until z is largest dimemsion.
int kz = r._d.maxDim();
int kx = kz+1; if (kx == 3) kx = 0;
int ky = kx+1; if (ky == 3) ky = 0;
// swap coord
Vec3f d_ = r._d.permute(kx, ky, kz);
p0_ = p0_.permute(kx, ky, kz);
p1_ = p1_.permute(kx, ky, kz);
p2_ = p2_.permute(kx, ky, kz);
// shear
float Sx = -d_.x()/d_.z();
float Sy = -d_.y()/d_.z();
float Sz = 1.0f/d_.z();
// perform transform (p.z doesn't matter now)
p0_[0] += Sx*p0_[2];
p0_[1] += Sy*p0_[2];
p1_[0] += Sx*p1_[2];
p1_[1] += Sy*p1_[2];
p2_[0] += Sx*p2_[2];
p2_[1] += Sy*p2_[2];
// signed edge function to judge whether origin's projection in projected triangle
float e0 = p1_.x()*p2_.y() - p2_.x()*p1_.y(); // p1 -> p2
float e1 = p2_.x()*p0_.y() - p0_.x()*p2_.y(); // p2 -> p0
float e2 = p0_.x()*p1_.y() - p1_.x()*p0_.y(); // p0 -> p1
// perform high-precision at edges
if (e0 == 0.0f || e1 == 0.0f || e2 == 0.0f) {
double p0xp1y = (double)p0_.x() * (double)p1_.y();
double p1xp0y = (double)p1_.x() * (double)p0_.y();
e2 = (float)(p0xp1y - p1xp0y);
double p1xp2y = (double)p1_.x() * (double)p2_.y();
double p2xp1y = (double)p2_.x() * (double)p1_.y();
e0 = (float)(p1xp2y - p2xp1y);
double p2xp0y = (double)p2_.x() * (double)p0_.y();
double p0xp2y = (double)p0_.x() * p2_.y();
e1 = (float)(p2xp0y - p0xp2y);
}
// edge detection
if ( (e0 < 0 || e1 < 0 || e2 < 0) && (e0 > 0 || e2 > 0 || e1 > 0) )
return false;
float det = e0 + e1 + e2;
if (det == 0) // on edge
return false;
// compute z value (i.e. t)
p0_[2] *= Sz;
p1_[2] *= Sz;
p2_[2] *= Sz;
float tScaled = e0*p0_.z() + e1*p1_.z() + e2*p2_.z();
if (det < 0 && (tScaled >= 0 || tScaled < r._tmax*det) ) return false;
else if (det > 0 && (tScaled <= 0 || tScaled > r._tmax*det) ) return false;
return true;
}
Unlike pbrt, my implement use right-hand coordinate and row-major matrices, but I think this do not influence the correctness of this piece of code.
Most of the time, this function works fine. However, artifacts appears when I do occlusion test to the light sources, where r._tmax
has a specific value (the distance between intersection point and light sample point).
I tested the algorithm with a plane and a point light above it, however just use the position of the light, radiance is not computed, and this is why I think the problem is caused by intersection implementation. If unoccluded, a value is returned as radiance, otherwise zero.
Vec3f PointEmitter::evalDirect(std::shared_ptr<Scene> scene, const IntersectInfo &isect_info, const Point2f &u) const {
Vec3f wi;
float light_pdf;
VisibilityTester vt;
Vec3f Li_sampled = sample_Li(isect_info, wi, light_pdf, vt, u);
Vec3f f = isect_info._bsdf->f(isect_info._wo, wi)*std::abs(isect_info._shading._n.dot(wi));
return (scene->intersectP(IntersectInfo(_p).spawnRayTo(isect_info))) *100;
}
And another open box-like polygon, with same light:
There's a circle-like region obviously stands out on the surface.
I consider this comes from the float point error when comparing r_tmax
and the distance. But the code already uses double-precision. Can anyone help?