Add new primitive: TORUS 🍩

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 = 21
+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()

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)