Add new primitive: TORUS 🍩

[munechika-koyo]

Hello,

I would like to create a PR about a new Torus primitive. Additionally, I implemented required functionalities to solve quartic, cubic equations, etc.

If you agree to this PR, then I will add or modify the documentations written about primitives' sections. The example rendering of a torus primitive is below, the script of which is demos/primitives/simple_torus.py.

I would appreciate it if you would review my codes and comment this PR.

[munechika-koyo]

@vsnever @mattngc @CnlPepper This feature is not prioritized over other issues, but how do you feel about adding this?

[CnlPepper]

I've not had a chance to look at this, but, assuming you haven't, could you make test demo that checks it works correctly with the CSG operations (tests next intersection)?

[munechika-koyo]

@CnlPepper I tried to calculate csg oprations by:

• Torus + Sphere(right) - Sphere(left) - Box(front) (Copper matrial)
• Torus * Cylinder (Glass material)

Here is the script I used:

from matplotlib import pyplot as plt

from raysect.optical import ConstantSF, Point3D, World, d65_white, rotate, translate
from raysect.optical.library.metal import Copper
from raysect.optical.material import Lambert, UniformSurfaceEmitter
from raysect.optical.library import schott
from raysect.optical.observer import PinholeCamera, RGBAdaptiveSampler2D, RGBPipeline2D
from raysect.primitive import (
    Box,
    Cylinder,
    Torus,
    Sphere,
    Union,
    Intersect,
    Subtract,
)


# glass matrial
glass = schott("N-BK7")

world = World()

# Toruses
torus1 = Torus(1.0, 0.5)
torus2 = Torus(1.0, 0.5)

# Spheres
sphere1 = Sphere(0.6, transform=translate(1.0, 0.0, 0.0))
sphere2 = Sphere(0.6, transform=translate(-1.0, 0.0, 0.0))

# Box
sqrt2 = 2 ** 0.5
box = Box(
    Point3D(-1.6, -1.6, -1.6),
    Point3D(1.6, 1.6, 1.6),
    transform=translate(0.0, 2.21, 0.0),
)

# cylinder
cylinder = Cylinder(0.6, 2.0, transform=translate(0.0, 1.0, 0.0))

# Torus1 - Sphere1 + Sphere2 - Box
Subtract(
    Union(Subtract(torus1, sphere1), sphere2),
    box,
    world,
    transform=translate(0.0, 0.0, 0.6),
    material=Copper(),
)

# Torus2 * Cylinder
Intersect(
    torus2,
    cylinder,
    world,
    transform=translate(0.0, 0.3, 0.6),
    material=glass,
)


# floor
Box(
    Point3D(-100, -100, -10),
    Point3D(100, 100, 0),
    parent=world,
    material=Lambert(ConstantSF(1.0)),
)

# emitter
Cylinder(
    3.0,
    100.0,
    parent=world,
    transform=translate(0, 0, 8) * rotate(90, 0, 0) * translate(0, 0, -50),
    material=UniformSurfaceEmitter(d65_white, 1.0),
)

# camera
rgb = RGBPipeline2D(display_unsaturated_fraction=0.995)
sampler = RGBAdaptiveSampler2D(rgb, min_samples=500, fraction=0.1, cutoff=0.01)
camera = PinholeCamera(
    (512, 512),
    parent=world,
    transform=rotate(0, 45, 0) * translate(0, 0, 5) * rotate(0, -180, 0),
    pipelines=[rgb],
    frame_sampler=sampler,
)
camera.spectral_bins = 21
camera.spectral_rays = 1
camera.pixel_samples = 250
camera.ray_max_depth = 10000
camera.ray_extinction_min_depth = 3
camera.ray_extinction_prob = 0.01


# start ray tracing
plt.ion()
for p in range(0, 1000):
    print(f"Rendering pass {p}...")
    camera.observe()
    print()

plt.ioff()
rgb.display()
plt.show()

[vsnever]

I think the implementation of the next_intersection() method is inconsistent with its intended behaviour. The description of this method says: https://github.com/raysect/source/blob/20f57259136ded2db692a967e486ebec88066b8a/raysect/core/scenegraph/primitive.pyx#L105-L119 Unlike any other primitive in Raysect, there are up to four possible intersections of a ray with a torus, but the method always stops after the second intersection.

    cpdef Intersection next_intersection(self):

        if not self._further_intersection:
            return None

        # this is the 2nd intersection
        self._further_intersection = False
        return self._generate_intersection(self._cached_ray, self._cached_origin, self._cached_direction, self._next_t)

The method is used only in the CSGPrimitive._identify_intersection() method.

As a result, extra intersections are generated in some cases. The output of this script:

from raysect.optical import Point3D, Vector3D, World, translate, InterpolatedSF
from raysect.optical.material import Dielectric
from raysect.primitive import Union, Torus, Cylinder
from raysect.optical.loggingray import LoggingRay


glass = Dielectric(index=InterpolatedSF([10, 10000], [1., 1.]),
                   transmission=InterpolatedSF([10, 10000], [1., 1.]),
                   transmission_only=True)

world = World()

torus = Torus(1.0, 0.5)

cylinder = Cylinder(1.0, 1.0, transform=translate(0, 0, -0.5))

union = Union(torus, cylinder, material=glass, parent=world)

ray = LoggingRay(origin=Point3D(2., 0, 0), direction=Vector3D(-1, 0, 0))
ray.trace(world)
for intersection in ray.log:
    print(intersection.hit_point, intersection.exiting)
    print()

is

Point3D(1.5, 0.0, 0.0) False

Point3D(1.4999997776219791, 0.0, 0.0) False

Point3D(-1.0, 0.0, 0.0) True

Point3D(-1.5, 0.0, 0.0) True

The intersection at Point3D(-1.0, 0.0, 0.0) is an extra one because the ray is not exciting the primitive at (-1, 0, 0).