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?

[vsnever]

I think, this is an excellent implementation, @munechika-koyo. I have only minor comments.

Some plasma diagnostics, such as ITER CXRS, use toric lenses and mirrors, so having the Torus primitive in Raysect would be very useful for Cherab as well.

If @CnlPepper approves this PR, could you please update both primitives.rst and api_reference/primitives/geometric_primitives.rst?

[vsnever]

As above, the function returns only 0, 2 or 4. The special cases of 1 and 3 real roots are not treated by this function separately from other cases due to machine epsilon.

[vsnever]

I think you don't need to add BOX_PADDING twice here. It's enough to add padding when you expand the bounding box.

[munechika-koyo]

I agree with you, and will remove it here.

[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]

Just checked on the hollow toroids

torus1 = Subtract(Torus(1.0, 0.5), Torus(1.0, 0.47))
torus2 = Subtract(Torus(1.0, 0.5), Torus(1.0, 0.47))

that everything works fine when rays hit the surface from the inside.

[vsnever]

I think the implementation of the next_intersection() method is inconsistent with its intended behaviour. The description of this method says:

source/raysect/core/scenegraph/primitive.pyx

Lines 105 to 119 in 20f5725

	cpdef Intersection next_intersection(self):
	"""
	Virtual method - to be implemented by derived classes.
	Returns the next intersection of the ray with the primitive along the
	ray path.
	This method may only be called following a call to hit(). If the ray
	has further intersections with the primitive, these may be obtained by
	repeatedly calling the next_intersection() method. Each call to
	next_intersection() will return the next ray-primitive intersection
	along the ray's path. If no further intersections are found or
	intersections lie outside the ray parameters then next_intersection()
	will return None.

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

[CnlPepper]

Oh this looks fun! I'd need to check this over with a fine tooth comb.... has the comment from @vsnever regarding the next_intersection() implementation been addressed?

[munechika-koyo]

@CnlPepper Thank you for commenting on this stale PR.
Since the implementation of the 4-times intersection is very specific to the torus primitive, I wonder whether we should modify other related codes.
Also, I haven't yet found a workaround for the issues mentioned above.
Could you please help make some changes to resolve them?

[CnlPepper]

The next_intersection() implementation just needs to return each interection found in order. For most of our current primitives (except mesh) there are only two intersections. So hit() returns the first, with the second cached. A subsequent call to next_intersection() returns the second, cached one. For the torus you'll need to cache up to 3 addition intersections on the object and return them one by one, in ray distance order for each call of next intersection.

[munechika-koyo]

Thank you for giving me informative advice!
I will try to implement the up to 3 hit-point caching functionality, and confirm if the above issue will be resolved.

[munechika-koyo]

The result of executing the above script that @vsnever proposed is as follows:

Point3D(1.5, 0.0, 0.0) False

Point3D(1.4999997776219791, 0.0, 0.0) False

Point3D(-1.5, 0.0, 0.0) True

It seems to be an appropriate behavior.
I am wondering a bit why the second one appeared, and the hit points were not symmetrically represented.

[munechika-koyo]

Note that some lines are still incompatible with Cython 3.1, and we should address them at #447 if this PR is merged beforehand.

[munechika-koyo]

I would like to modify the document about the primitive page by rendering the top image, adding a torus primitive.
Does anyone know where the script that renders its image is located? Or should I write it from scratch?

[CnlPepper]

I would like to modify the document about the primitive page by rendering the top image, adding a torus primitive. Does anyone know where the script that renders its image is located? Or should I write it from scratch?

I think that was a one-off script I created in the early Raysect development. I have a quick look but couldn't find it. If you recreate it, feel free to add it to demos/primitives if you like.

I'm pretty sure the glass material I used was N-Bk7. You'll need to hover the primitives slightly off the ground plane (1e-6 m or so) to prevent surface fighting with the ground plane.

[munechika-koyo]

I created a primitive image below. I will replace the previous image and add the script file, rendering it.

[munechika-koyo]

When executing the generate_meson_files.py script, almost all meson.build files were modified, which seems to reorder the elements of the source files lists.
The pathlib.Path.iterdir() here:

source/dev/generate_meson_files.py

Line 197 in 5897707

for child in path.iterdir():

couldn't return a path object in a fixed order. (It perhaps depends on the OS platform?)
To consistently yield the same result without any source change, we likely need a workaround, such as using sorting or globbing.

[CnlPepper]

Interesting, I'd not encountered that thus far. I'll make all the lists sorted to avoid the issue in the future. Thanks for spotting that.

[munechika-koyo]

I believe this PR is ready to merge!