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:

```python 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).