<p><em>Many</em> years ago I worked on a ray tracer that handled parametric surfaces, so this is unlikely to be state of the art, but, IIRC, I used a combination of interval arithmetic with (binary?) subdivision and Newton-Rhapson.</p>

<p>The interval arithmetic + subdivision constructed (conservative) bounding boxes which could be used for intersection rejection. I think I <em>may</em> also have used interval arithmetic on the 1st derivative to help determine when it was safe to launch into Newton-Rhapson. Typically the Newton steps converged very rapidly. </p>

<p>You can also pre-dice your parametric surfaces and put them into an acceleration structure (Voxel grid/BVH) to speed up the process.</p>

<p>This work was based on/inspired by Daniel Toth's 1985 paper <a href="https://www.cs.drexel.edu/~david/Classes/CS586/Papers/p171-toth.pdf" rel="noreferrer">"On Ray Tracing Parametric Surfaces"</a></p>

CC BY-SA 3.0
2017-03-14T17:30:32.940Z
Thanks for the answer but I have to admit that I don't understand everything. Could you put some concrete example or code so that I can't see more clearly please ?
Telokis
CC BY-SA 3.0
2017-03-15T12:10:07.650Z
Ok. I haven't got time today, but will try to add some more info soon(ish).
Simon F
CC BY-SA 3.0
2017-03-16T17:00:50.857Z
What you do suggest is that I keep using my method of constructing triangles but that I should use a spatial subdivision algorithm to "fasten" the calculation?
Telokis
CC BY-SA 3.0
2017-03-17T10:19:44.837Z
It seems to me you have two basic options though, in some senses, they can both converge to an 'equivalent' result: (a) Interval arithmetic, to weed out regions followed by Newton-Rhapson to solve for the ray intersection or (b) dice (i.e. tessellate)  your object into "sufficiently small" triangles that you don't see discontinuity artefacts (as in your first image). In a sense these are "similar" since in (a) each Newton-R iteration approximates the surface locally as a plane (like a small triangle). I suspect b will be much easier to implement but tricky to estimate ideal subdivision level.
Simon F
CC BY-SA 3.0
2017-03-14T09:44:21.260Z
2017-03-14T09:44:21.260Z
Simon F
<p>In general when doing parametric in a ray tracer you need a solution for </p>

<p>$\begin{cases}  P = v * t + C \\
 P = f(u,v)
\end{cases}$</p>

<p>for the lowest $t$ where $f(u,v)$ is your parametric function $C$ is the camera position and $v$ is the ray direction.</p>

<p>There are a few general ways, for example if you can tell whether 2 points are on the same side or not you can take 2 points on the ray on opposite sides and binary search until you have the depth. You can find the first pair by marching along the ray.</p>

CC BY-SA 3.0
2017-03-14T17:36:20.940Z
What do you mean by "the same side"? I know the general idea of "how to do it" but I don't know how it happens in practice. Povray gives such an accurate and pretty result, I don't understand how it can possibly do it.
Telokis
CC BY-SA 3.0
2017-03-14T09:47:48.797Z
2017-03-14T09:47:48.797Z
ratchet freak
<p>After taking a look at the <a href="https://www.wikiwand.com/en/M%C3%B6bius_strip#/Geometry_and_topology" rel="nofollow noreferrer">Mobius strip</a>, I noticed its equation is really simple and tried to add it into my Raytracer.</p>

<p>I tried a "naive" way by simply generating <code>N</code> triangles attached to each other to obtain the desired shape. While this approach works, the result it not really pretty:</p>

<p><a href="https://i.sstatic.net/1w5Kx.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/1w5Kx.png" alt="IMG"></a></p>

<p>(By the way I probably have an issue with my normals but I don't know where it comes from.)</p>

<p>I tried it with PovRay and the result was astonishing. Perfectly smooth strip made in a far far FAR smaller time than mine. I'm pretty sure Povray is <strong>well</strong> optimized but I also think it won't generate triangles like I did.</p>

<p><a href="https://i.sstatic.net/DN5p4.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/DN5p4.png" alt="povray"></a></p>

<p>In case that might help, here is the actual code used (C++) :</p>

<pre><code>float step = .1f;
float halW = 0.5f;

_facets.clear();

auto lambda = [this] (float v, float t) {
  Vec_t p;

  float cdv = Tools::Cos(2 * v);
  float sdv = Tools::Sin(2 * v);
  float ctv = Tools::Cos(v);
  float stv = Tools::Sin(v);

  float c = 2 + t * ctv;

  p.x = c * cdv;
  p.z = c * sdv;
  p.y = t * stv;

  return p;
};

for (float v = 0.f; v &lt; Globals::PI; v += step)
{
  if (v &gt; Globals::PI)
    v = Globals::PI;

  for (float t = -halW; t &lt; halW; t += step)
  {
    if (t &gt; halW)
      t = halW;

    Vec3 p1 = lambda(v, t);
    Vec3 p2 = lambda(v + step, t);
    Vec3 p3 = lambda(v, t + step);
    Vec3 p4 = lambda(v + step, t + step);
    _facets.emplace_back(p1, p2, p3);
    _facets.emplace_back(p3, p2, p4);
  }
}
</code></pre>

<p><strong>TL;DR</strong></p>

<p>How can I handle parametric surfaces like <a href="https://www.wikiwand.com/en/M%C3%B6bius_strip#/Geometry_and_topology" rel="nofollow noreferrer">this one</a> in raytracing?</p>

<p><strong>Edit</strong>
After letting the above algorithm run for about 20 hours, I got a way prettier result (with 3 torsions instead of 1)</p>

<p><a href="https://cdn.discordapp.com/attachments/148687431728627713/291972504245895178/unknown.png" rel="nofollow noreferrer"><img src="https://cdn.discordapp.com/attachments/148687431728627713/291972504245895178/unknown.png" alt="Better result"></a></p>

CC BY-SA 3.0
2017-03-19T01:14:16.007Z
For the specific problem outlined in red in the first image, it looks like you might need to check which side of the surface you are approaching before deciding which way the normal should point, otherwise some regions will be shaded incorrectly. This shows up in particular for a Möbius strip, since there must be a point at which the normals switch direction due to the twist. There must be some adjacent triangles that have almost opposite normal directions (like the light and dark adjacent triangles in the image).
trichoplax is on Codidact now
CC BY-SA 3.0
2017-03-19T02:33:45.410Z
Yes, I thought that as well but I am setting a `isInside` flag in the triangle intersection method. If it is true I then negate the normal vector. The flag is set if `det < 0` using an intersection method I don't recall the name for now
Telokis
CC BY-SA 3.0
2017-03-19T03:04:25.550Z
I am using the [Möller–Trumbore intersection algorithm](https://www.wikiwand.com/en/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm)
Telokis
<raytracing><rendering><c++>
How to handle a parametric equation in raytracing?
Telokis
# How to handle a parametric equation in raytracing?

(Asked by: username_0 on 2017-03-14)

After taking a look at the [Mobius strip](https://www.wikiwand.com/en/M%C3%B6bius_strip#/Geometry_and_topology), I noticed its equation is really simple and tried to add it into my Raytracer.

I tried a "naive" way by simply generating `N` triangles attached to each other to obtain the desired shape. While this approach works, the result it not really pretty:

[![IMG](https://i.sstatic.net/1w5Kx.png)](https://i.sstatic.net/1w5Kx.png)

(By the way I probably have an issue with my normals but I don't know where it comes from.)

I tried it with PovRay and the result was astonishing. Perfectly smooth strip made in a far far FAR smaller time than mine. I'm pretty sure Povray is **well** optimized but I also think it won't generate triangles like I did.

[![povray](https://i.sstatic.net/DN5p4.png)](https://i.sstatic.net/DN5p4.png)

In case that might help, here is the actual code used (C++) :

```
float step = .1f;
float halW = 0.5f;

_facets.clear();

auto lambda = [this] (float v, float t) {
  Vec_t p;

  float cdv = Tools::Cos(2 * v);
  float sdv = Tools::Sin(2 * v);
  float ctv = Tools::Cos(v);
  float stv = Tools::Sin(v);

  float c = 2 + t * ctv;

  p.x = c * cdv;
  p.z = c * sdv;
  p.y = t * stv;

  return p;
};

for (float v = 0.f; v < Globals::PI; v += step)
{
  if (v > Globals::PI)
    v = Globals::PI;

  for (float t = -halW; t < halW; t += step)
  {
    if (t > halW)
      t = halW;

    Vec3 p1 = lambda(v, t);
    Vec3 p2 = lambda(v + step, t);
    Vec3 p3 = lambda(v, t + step);
    Vec3 p4 = lambda(v + step, t + step);
    _facets.emplace_back(p1, p2, p3);
    _facets.emplace_back(p3, p2, p4);
  }
}

```

**TL;DR**

How can I handle parametric surfaces like [this one](https://www.wikiwand.com/en/M%C3%B6bius_strip#/Geometry_and_topology) in raytracing?

**Edit**
After letting the above algorithm run for about 20 hours, I got a way prettier result (with 3 torsions instead of 1)

[![Better result](https://cdn.discordapp.com/attachments/148687431728627713/291972504245895178/unknown.png)](https://cdn.discordapp.com/attachments/148687431728627713/291972504245895178/unknown.png)

**Tags:** raytracing, rendering, c++

> **Comments:**
> **username_1**: For the specific problem outlined in red in the first image, it looks like you might need to check which side of the surface you are approaching before deciding which way the normal should point, otherwise some regions will be shaded incorrectly. This shows up in particular for a Möbius strip, since there must be a point at which the normals switch direction due to the twist. There must be some adjacent triangles that have almost opposite normal directions (like the light and dark adjacent triangles in the image).  
> **username_0**: Yes, I thought that as well but I am setting a `isInside` flag in the triangle intersection method. If it is true I then negate the normal vector. The flag is set if `det < 0` using an intersection method I don't recall the name for now  
> **username_0**: I am using the [Möller–Trumbore intersection algorithm](https://www.wikiwand.com/en/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm)  

---

## 2 answers

### Answer by username_2 on 2017-03-14. Score: 6

*Many* years ago I worked on a ray tracer that handled parametric surfaces, so this is unlikely to be state of the art, but, IIRC, I used a combination of interval arithmetic with (binary?) subdivision and Newton-Rhapson.

The interval arithmetic + subdivision constructed (conservative) bounding boxes which could be used for intersection rejection. I think I *may* also have used interval arithmetic on the 1st derivative to help determine when it was safe to launch into Newton-Rhapson. Typically the Newton steps converged very rapidly.

You can also pre-dice your parametric surfaces and put them into an acceleration structure (Voxel grid/BVH) to speed up the process.

This work was based on/inspired by Daniel Toth's 1985 paper ["On Ray Tracing Parametric Surfaces"](https://www.cs.drexel.edu/~david/Classes/CS586/Papers/p171-toth.pdf)

> **Comments:**
> **username_0**: Thanks for the answer but I have to admit that I don't understand everything. Could you put some concrete example or code so that I can't see more clearly please ?  
> **username_2**: Ok. I haven't got time today, but will try to add some more info soon(ish).  
> **username_0**: What you do suggest is that I keep using my method of constructing triangles but that I should use a spatial subdivision algorithm to "fasten" the calculation?  
> **username_2**: It seems to me you have two basic options though, in some senses, they can both converge to an 'equivalent' result: (a) Interval arithmetic, to weed out regions followed by Newton-Rhapson to solve for the ray intersection or (b) dice (i.e. tessellate)  your object into "sufficiently small" triangles that you don't see discontinuity artefacts (as in your first image). In a sense these are "similar" since in (a) each Newton-R iteration approximates the surface locally as a plane (like a small triangle). I suspect b will be much easier to implement but tricky to estimate ideal subdivision level.  

---

### Answer by username_3 on 2017-03-14. Score: 2

In general when doing parametric in a ray tracer you need a solution for

$\begin{cases} P = v \* t + C \\
P = f(u,v)
\end{cases}$

for the lowest $t$ where $f(u,v)$ is your parametric function $C$ is the camera position and $v$ is the ray direction.

There are a few general ways, for example if you can tell whether 2 points are on the same side or not you can take 2 points on the ray on opposite sides and binary search until you have the depth. You can find the first pair by marching along the ray.

> **Comments:**
> **username_0**: What do you mean by "the same side"? I know the general idea of "how to do it" but I don't know how it happens in practice. Povray gives such an accurate and pretty result, I don't understand how it can possibly do it.  

---

