<p>I'm building a software renderer but I think there might be some problem with my ray-triangle intersection accuracy.</p>

<p>I implemented the algorithm referring to <code>pbrt-v3</code> without some components which I think have little influence on the intersection test.</p>

<p>Below is the code perform intersection test:</p>

<pre><code>bool Triangle::intersectP(const Ray &amp;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 -&gt; p2
    float e1 = p2_.x()*p0_.y() - p0_.x()*p2_.y();   // p2 -&gt; p0
    float e2 = p0_.x()*p1_.y() - p1_.x()*p0_.y();   // p0 -&gt; 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 &lt; 0 || e1 &lt; 0 ||  e2 &lt; 0) &amp;&amp; (e0 &gt; 0 || e2 &gt; 0 || e1 &gt; 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 &lt; 0 &amp;&amp; (tScaled &gt;= 0 || tScaled &lt; r._tmax*det) ) return false;
    else if (det &gt; 0 &amp;&amp; (tScaled &lt;= 0 || tScaled &gt; r._tmax*det) ) return false;

    return true;
}
</code></pre>

<p>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.</p>

<p>Most of the time, this function works fine. However, artifacts appears when I do occlusion test to the light sources, where <code>r._tmax</code> has a specific value (the distance between intersection point and light sample point).</p>

<p>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.</p>

<pre><code>Vec3f PointEmitter::evalDirect(std::shared_ptr&lt;Scene&gt; scene, const IntersectInfo &amp;isect_info, const Point2f &amp;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-&gt;f(isect_info._wo, wi)*std::abs(isect_info._shading._n.dot(wi));

    return (scene-&gt;intersectP(IntersectInfo(_p).spawnRayTo(isect_info))) *100;
}
</code></pre>

<p><a href="https://i.sstatic.net/dFLQd.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/dFLQd.jpg" alt="enter image description here"></a></p>

<p>And another open box-like polygon, with same light:
<a href="https://i.sstatic.net/8EziQ.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/8EziQ.jpg" alt="enter image description here"></a>
There's a circle-like region obviously stands out on the surface.</p>

<p>I consider this comes from the float point error when comparing <code>r_tmax</code> and the distance. But the code already uses double-precision. Can anyone help?</p>

CC BY-SA 4.0
2018-11-10T21:22:48.347Z
Have you added a bias to your rays? Adding a bias means that you offset the ray position by the normal of the surface multiplied by a very small number, so that the ray is just off the surface. This means that the triangle that the starting point of a ray is on, is actually a bit behind the ray meaning that $t<0$. If you do not have the bias, the triangle would be at $t=0$, however we are working with limited precision and then sometimes it is below zero and sometimes above zero. It would not matter how much precision you have, you would always get this. Every render engine implements this.
bram0101
CC BY-SA 4.0
2018-11-10T23:19:26.577Z
@bram0101 Hi. Thanks for the advice. However I'm not sure if this solves the problem, though I didn't add the bias. This is because, in my code, the occlusion test ray starts from the point light and I don't think phenomenon is caused by the $t<0$ issue. But I'll definitely add the bias to make the code robust.
jinglei
CC BY-SA 4.0
2018-11-12T15:15:09.903Z
+1 to Rahul.  Definitely looks like "acne". With floating point, it's always worth remembering:  "Floating point numbers are like piles of sand; every time you move them around, you lose a little sand and pick up a little dirt. — Brian Kernighan and P.J. Plauger"
Simon F
CC BY-SA 4.0
2018-11-12T16:29:25.163Z
@Rahul Thank you for the reference link. I've added the bias and the phenomena seems gone away. However, one more question on how I choose the bias value. I used $1e-6$ at first and some of the "acnes" still exists, only after I used $1e-5$ did them disappear. But I'm afraid the bias value is too high?
jinglei
CC BY-SA 4.0
2018-11-19T06:42:00.807Z
It's been a long time, but you possibly need to make the bias relative to the magnitude to the input floats.
Simon F
<raytracing><c++><geometry><ray>
ray-triangle intersection precision - c++
jinglei
# ray-triangle intersection precision - c++

(Asked by: username_0 on 2018-11-10)

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;
}

```

[![enter image description here](https://i.sstatic.net/dFLQd.jpg)](https://i.sstatic.net/dFLQd.jpg)

And another open box-like polygon, with same light:
[![enter image description here](https://i.sstatic.net/8EziQ.jpg)](https://i.sstatic.net/8EziQ.jpg)
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?

**Tags:** raytracing, c++, geometry, ray

> **Comments:**
> **username_1**: Have you added a bias to your rays? Adding a bias means that you offset the ray position by the normal of the surface multiplied by a very small number, so that the ray is just off the surface. This means that the triangle that the starting point of a ray is on, is actually a bit behind the ray meaning that $t<0$. If you do not have the bias, the triangle would be at $t=0$, however we are working with limited precision and then sometimes it is below zero and sometimes above zero. It would not matter how much precision you have, you would always get this. Every render engine implements this.  
> **username_0**: @username_1 Hi. Thanks for the advice. However I'm not sure if this solves the problem, though I didn't add the bias. This is because, in my code, the occlusion test ray starts from the point light and I don't think phenomenon is caused by the $t<0$ issue. But I'll definitely add the bias to make the code robust.  
> **username_2**: +1 to Rahul.  Definitely looks like "acne". With floating point, it's always worth remembering:  "Floating point numbers are like piles of sand; every time you move them around, you lose a little sand and pick up a little dirt. — Brian Kernighan and P.J. Plauger"  
> **username_0**: @Rahul Thank you for the reference link. I've added the bias and the phenomena seems gone away. However, one more question on how I choose the bias value. I used $1e-6$ at first and some of the "acnes" still exists, only after I used $1e-5$ did them disappear. But I'm afraid the bias value is too high?  
> **username_2**: It's been a long time, but you possibly need to make the bias relative to the magnitude to the input floats.  

---

## 0 answers

