Add new primitive: TORUS 🍩

demos/primitives/raysect_primitives.py

@@ -0,0 +1,126 @@
+# External imports
+from time import strftime
+
+from matplotlib import pyplot as plt
+
+# Raysect imports
+from raysect.optical import Point3D, World, d65_white, rotate, translate
+from raysect.optical.library import schott
+from raysect.optical.material import Checkerboard, Lambert
+from raysect.optical.observer import PinholeCamera, RGBAdaptiveSampler2D, RGBPipeline2D
+from raysect.primitive import Box, Cone, Cylinder, Parabola, Sphere, Torus
+
+# 1. Create Primitives
+# --------------------
+
+# Box defining the ground plane
+ground = Box(
+    lower=Point3D(-50, -0.01, -50), upper=Point3D(50, 0.0, 50), material=Lambert()
+)
+
+# checker board wall that acts as emitter
+emitter = Box(
+    lower=Point3D(-100, -100, 10),
+    upper=Point3D(100, 100, 10.1),
+    material=Checkerboard(4, d65_white, d65_white, 0.1, 2.0),
+)
+
+# Primitive showcasing all geometric features
+# Note that the primitives must be displaced slightly above the ground plane to prevent numerically issues that could
+# cause a light leak at the intersection between the objects and the ground.
+cylinder = Cylinder(
+    radius=1.5,
+    height=3.0,
+    transform=translate(1.5 * 3 + 1.0, 0.0001, 0) * rotate(0, 90, 0),
+    material=schott("N-BK7"),
+)
+cone = Cone(
+    radius=1.5,
+    height=3.0,
+    transform=translate(1.5 + 0.2, 0.0001, 0) * rotate(0, 90, 0),
+    material=schott("N-BK7"),
+)
+sphere = Sphere(
+    radius=1.5,
+    transform=translate(-1.5 - 0.2, 1.5 + 0.0001, 0),
+    material=schott("N-BK7"),
+)
+box = Box(
+    lower=Point3D(-1.5, 0.0, -1.5),
+    upper=Point3D(1.5, 3.0, 1.5),
+    transform=translate(-1.5 * 3 - 1.0, 0.0001, 0),
+    material=schott("N-BK7"),
+)
+parabola = Parabola(
+    radius=2.0,
+    height=1.0,
+    transform=translate(2.5, 1.0 + 0.0001, -5.0) * rotate(0, -90, 0),
+    material=schott("N-BK7"),
+)
+torus = Torus(
+    major_radius=1.0,
+    minor_radius=0.5,
+    transform=translate(-2.5, 0.5 + 0.0001, -5.0) * rotate(0, 90, 0),
+    material=schott("N-BK7"),
+)
+
+
+# 2. Add Observer
+# ---------------
+
+# Process the ray-traced spectra with the RGB pipeline.
+rgb = RGBPipeline2D(display_unsaturated_fraction=0.96)
+sampler = RGBAdaptiveSampler2D(
+    rgb, ratio=10, fraction=0.2, min_samples=2000, cutoff=0.01
+)
+
+# camera
+camera = PinholeCamera(
+    (512, 512), pipelines=[rgb], transform=translate(-7, 12, -15) * rotate(-25, -40, 0)
+)
+
+# camera - pixel sampling settings
+camera.fov = 45
+camera.pixel_samples = 250
+
+# camera - ray sampling settings
+camera.spectral_rays = 15
+camera.spectral_bins = 15
+camera.ray_max_depth = 100
+camera.ray_extinction_prob = 0.1
+camera.min_wavelength = 375.0
+camera.max_wavelength = 740.0
+
+
+# 3. Build Scenegraph
+# -------------------
+
+world = World()
+
+ground.parent = world
+emitter.parent = world
+camera.parent = world
+cylinder.parent = world
+cone.parent = world
+sphere.parent = world
+box.parent = world
+parabola.parent = world
+torus.parent = world
+
+# 4. Observe()
+# ------------
+name = "raysect_primitives"
+timestamp = strftime("%Y-%m-%d_%H-%M-%S")
+render_pass = 1
+plt.ion()
+while not camera.render_complete:
+    print(f"Rendering pass {render_pass}...")
+    camera.observe()
+    rgb.save(f"{name}_{timestamp}_pass_{render_pass}.png")
+    render_pass += 1
+    print()
+
+# display final result
+plt.ioff()
+rgb.display()
+plt.show()

demos/primitives/simple_torus.py

@@ -0,0 +1,64 @@
+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.observer import PinholeCamera, RGBAdaptiveSampler2D, RGBPipeline2D
+from raysect.primitive import Box, Cylinder, Torus
+
+world = World()
+
+# Torus
+torus = Torus(
+    1.0,
+    0.5,
+    world,
+    transform=translate(0, 0.0, 0.6),
+    material=Copper(),
+)
+
+# 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()

docs/source/api_reference/primitives/geometric_primitives.rst

@@ -21,3 +21,7 @@ Geometric Primitives
 .. autoclass:: raysect.primitive.Parabola
    :members: radius, height
    :show-inheritance:
+
+.. autoclass:: raysect.primitive.Torus
+   :members: major_radius, minor_radius
+   :show-inheritance:

docs/source/how_it_works.rst

@@ -58,7 +58,7 @@ Primitives
 Scenes in Raysect consist of primitives, which are the basic objects making up the scene. These are objects that rays
 can interact with, e.g light sources, lenses, mirrors, diffuse surfaces. Types of primitives:
 
-* Mathematically defined surfaces and solids (i.e. sphere, box, cylinder, cone).
+* Mathematically defined surfaces and solids (i.e. sphere, box, cylinder, cone, parabola, torus).
 * Constructive Solid Geometry Operators (union, intersect, subtract).
 * Tri-poly meshes optimised for millions of polygons, support instancing. Importers for STL and OBJ.
 

docs/source/primitives.rst

@@ -3,7 +3,7 @@
 Primitives
 **********
 
-The raysect primitives: sphere; box; cylinder; and cone.
+The raysect primitives: sphere; box; cylinder; cone; parabola; and torus.
 
 .. image:: images/raysect_primitives.png
    :align: center
@@ -30,6 +30,15 @@ Cone
 
 .. autoclass:: raysect.primitive.Cone
 
+Parabola
+~~~~~~~~
+
+.. autoclass:: raysect.primitive.Parabola
+
+Torus
+~~~~~
+
+.. autoclass:: raysect.primitive.Torus
 
 ==============
 CSG Operations

raysect/core/math/cython/utility.pxd

@@ -31,6 +31,8 @@
 
 cimport cython
 
+DEF EQN_EPS = 1.0e-9
+
 cdef int find_index(double[::1] x, double v) nogil
 
 cdef double interpolate(double[::1] x, double[::1] y, double p) nogil
@@ -64,12 +66,42 @@ cdef inline void swap_int(int *a, int *b) nogil:
     a[0] = b[0]
     b[0] = temp
 
+cdef inline void sort_three_doubles(double *a, double *b, double *c) nogil:
+    if a[0] > b[0]:
+        swap_double(a, b)
+    if b[0] > c[0]:
+        swap_double(b, c)
+        if a[0] > c[0]:
+            swap_double(a, c)
+
+cdef inline void sort_four_doubles(double *a, double *b, double *c, double *d) nogil:
+    if a[0] > b[0]:
+        swap_double(a, b)
+    if b[0] > c[0]:
+        swap_double(b, c)
+    if c[0] > d[0]:
+        swap_double(c, d)
+    if a[0] > b[0]:
+        swap_double(a, b)
+    if b[0] > c[0]:
+        swap_double(b, c)
+    if a[0] > b[0]:
+        swap_double(a, b)
+
+cdef inline bint is_zero(double v) nogil:
+    return v < EQN_EPS and v > -EQN_EPS
+
 @cython.cdivision(True)
 cdef inline double lerp(double x0, double x1, double y0, double y1, double x) nogil:
     return ((y1 - y0) / (x1 - x0)) * (x - x0) + y0
 
 cdef bint solve_quadratic(double a, double b, double c, double *t0, double *t1) nogil
 
+cdef int solve_cubic(double a, double b, double c, double d, double *t0, double *t1, double *t2) nogil
+
+cdef int solve_quartic(double a, double b, double c, double d, double e,
+                       double *t0, double *t1, double *t2, double *t3) nogil
+
 cdef bint winding2d(double[:,::1] vertices) nogil
 
 cdef bint point_inside_polygon(double[:,::1] vertices, double ptx, double pty)