diff --git a/demos/primitives/raysect_primitives.py b/demos/primitives/raysect_primitives.py new file mode 100644 index 00000000..9e6825fd --- /dev/null +++ b/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() diff --git a/demos/primitives/simple_torus.py b/demos/primitives/simple_torus.py new file mode 100644 index 00000000..91c0d66e --- /dev/null +++ b/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() diff --git a/docs/source/api_reference/primitives/geometric_primitives.rst b/docs/source/api_reference/primitives/geometric_primitives.rst index d0b833cb..a21d2b87 100644 --- a/docs/source/api_reference/primitives/geometric_primitives.rst +++ b/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: diff --git a/docs/source/how_it_works.rst b/docs/source/how_it_works.rst index d64f5916..ef6994ff 100644 --- a/docs/source/how_it_works.rst +++ b/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. diff --git a/docs/source/images/raysect_primitives.png b/docs/source/images/raysect_primitives.png index 0439ab93..cd2301f3 100644 Binary files a/docs/source/images/raysect_primitives.png and b/docs/source/images/raysect_primitives.png differ diff --git a/docs/source/primitives.rst b/docs/source/primitives.rst index ac997661..1a94b7ae 100644 --- a/docs/source/primitives.rst +++ b/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 diff --git a/raysect/core/math/cython/utility.pxd b/raysect/core/math/cython/utility.pxd index 55251357..71cc0303 100644 --- a/raysect/core/math/cython/utility.pxd +++ b/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) diff --git a/raysect/core/math/cython/utility.pyx b/raysect/core/math/cython/utility.pyx index 8f8bdcb9..693ea8c7 100644 --- a/raysect/core/math/cython/utility.pyx +++ b/raysect/core/math/cython/utility.pyx @@ -29,7 +29,7 @@ # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. -from libc.math cimport sqrt +from libc.math cimport sqrt, fabs, cbrt, acos, cos, M_PI cimport cython #TODO: Write unit tests! @@ -380,16 +380,16 @@ cdef bint solve_quadratic(double a, double b, double c, double *t0, double *t1) The a, b and c arguments are the three constants of the quadratic equation: f = a.x^2 + b.x^2 + c - + If the quadratic equation has 1 or 2 real roots, this function will return True. If there are no real roots this method will return False. - + The values of the real roots, are returned by setting the values of the memory locations pointed to by t0 and t1. In the case of a single root, both t0 and t1 will have the same value. If there are not roots, the values - of t0 and t1 will be undefined. + of t0 and t1 will be undefined. - :param double a: Quadratic constant. + :param double a: Quadratic constant. :param double b: Quadratic constant. :param double c: Quadratic constant. :param double t0: 1st root of the quadratic. @@ -419,6 +419,334 @@ cdef bint solve_quadratic(double a, double b, double c, double *t0, double *t1) return True +@cython.cdivision(True) +cdef int solve_cubic(double a, double b, double c, double d, double *t0, double *t1, double *t2) nogil: + """ + Calculates the roots of a cubic equation. + + The a, b, c and d arguments are the four constants of the cubic equation: + + f = a.x^3 + b.x^2 + c.x + d + + The cubic equation has 1, 2 or 3 real roots, and this function returns either 1 or 3, but in a special case, + 2 of the 3 roots found by this function will be equal to each other up to machine epsilon. + + + The values of the roots, are returned by setting the values of the memory locations pointed to by t0, t1, and t2. + + In the case of three real roots, the roots themselves back in t0, t1 and t2. + In the case of two real roots, a pair in t0, t1, and t2 will have the same value, and the other is a different one. + In the case of one real root, t0 will be the real root, and t1 +- i * t2 will a pair of complex-conjugated roots. + + The practical algorithm is followed by https://quarticequations.com + + :param double a: Cubic constant. + :param double b: Cubic constant. + :param double c: Cubic constant. + :param double d: Cubic constant. + :param double t0: 1st real root. + :param double t1: either 2nd real root or real part of coplex-conjugated roots. + :param double t2: either 3rd real root or imaginary part of coplex-conjugated roots. + :return: Number of real roots. + :rtype: int + """ + cdef: + double q, r, sq_b, cb_q, D, A, z0, phi, u + + # normal form: x^3 + bx^2 + cx + d = 0 + b /= a + c /= a + d /= a + + # convert depressed cubic: y^3 + 3qy - 2r = 0 + sq_b = b * b + q = (3.0 * c - sq_b) / 9.0 + r = (c * b - 3.0 * d) / 6.0 - b * sq_b / 27.0 + + # calculate discriminant + cb_q = q * q * q + D = cb_q + r * r + + # one real root and a pair of complex-conjugate roots + if D > 0: + A = cbrt(fabs(r) + sqrt(D)) + if r < 0: + z0 = q / A - A + else: + z0 = A - q / A + + t0[0] = z0 - b / 3.0 + t1[0] = -0.5 * z0 - b / 3.0 + t2[0] = 0.5 * sqrt(3.0) * (A + q / A) + + return 1 + + # Trigonometric solution for three real roots + else: + if is_zero(q): + phi = 0.0 + + # otherwise q < 0 because of D = q^3 + r^2 < 0 + else: + phi = acos(r / sqrt(-cb_q)) / 3.0 + + u = 2.0 * sqrt(-q) + + t0[0] = u * cos(phi) - b / 3.0 + t1[0] = -u * cos(phi + M_PI / 3.0) - b / 3.0 + t2[0] = -u * cos(phi - M_PI / 3.0) - b / 3.0 + + return 3 + + +cdef int solve_biquadratic(double a, double c, double e, double *t0, double *t1, double *t2, double *t3) nogil: + """ + Calculate the real roots of a bi quadratic equation. + + The a, c, and e arguments are the 3 constants of the biquadratic equation: + + f = a.x^4 + c.x^2 + e + + The biquadratic equation has 0, 2, or 4 real roots, and this function will return the number of real roots. + + The values of the real roots, are returned by setting the values of the + memory locations pointed to by t0, t1, t2, and t3. If there are two or four real roots, + the values of t2 and t3 will be undefined. If there is no real root, + all values will be undefined. + + :param double a: Biquadratic constant. + :param double c: Biquadratic constant. + :param double e: Biquadratic constant. + :param double t0: 1st root of the biquadratic. + :param double t1: 2nd root of the biquadratic. + :param double t2: 3rd root of the biquadratic. + :param double t3: 4th root of the biquadratic. + :return: Number of real roots. + :rtype: int + """ + cdef double s0, s1, sx0, sx1 + + # solve quadratic for x^2 + if not solve_quadratic(a, c, e, &s0, &s1): + return 0 + + # ensure s0 < s1 + if s0 > s1: + swap_double(&s0, &s1) + + # 0 <= s0 <= s1, 4 real roots + if s0 >= 0: + sx0 = sqrt(s0) + sx1 = sqrt(s1) + t0[0] = -sx1 + t1[0] = -sx0 + t2[0] = sx0 + t3[0] = sx1 + return 4 + + # s0 < 0 <= s1, 2 real roots + elif s1 >= 0: + sx1 = sqrt(s1) + t0[0] = -sx1 + t1[0] = sx1 + return 2 + + # s0 < s1 <= 0, no real root + else: + return 0 + + +cdef int _solve_depressed_quartic(double p, double q, double r, double *t0, double *t1, double *t2, double *t3) nogil: + """ + Solve depressed quartic: x^4 + p.x^2 + q.x + r + """ + cdef: + int num + double s0, sigma, A, B, sq_A, sq_B + + if q > 0: + sigma = 1.0 + else: + sigma = -1.0 + + # q = 0 => x^4 + p.x^2 + r = 0 + if is_zero(q): + return solve_biquadratic(1.0, p, r, t0, t1, t2, t3) + + # solve resolvent cubic: t^3 - 2pt^2 + (p^2-4r)t + q^2 = 0 + # using Van der Waerden's method + num = solve_cubic(1.0, -2.0 * p, p * p - 4.0 * r, q * q, t0, t1, t2) + + if num > 1: + # sort roots to t0 < t1 < t2 + sort_three_doubles(t0, t1, t2) + + # t0 <= 0 => t1*t2 >= 0 because vieta's therem: -t0*t1*t2 = q^2 + if t0[0] <= 0: + s0 = sqrt(-t0[0]) + A = -t1[0] - t2[0] - 2.0 * sigma * sqrt(t1[0] * t2[0]) + B = -t1[0] - t2[0] + 2.0 * sigma * sqrt(t1[0] * t2[0]) + + # four real roots + if A >= 0 and B >= 0: + sq_A = sqrt(A) + sq_B = sqrt(B) + t0[0] = 0.5 * (s0 + sq_A) + t1[0] = 0.5 * (s0 - sq_A) + t2[0] = 0.5 * (-s0 + sq_B) + t3[0] = 0.5 * (-s0 - sq_B) + return 4 + + # two real roots + elif A < 0 and B >= 0: + sq_B = sqrt(B) + t0[0] = 0.5 * (-s0 + sq_B) + t1[0] = 0.5 * (-s0 - sq_B) + return 2 + + # two real roots + elif A >= 0 and B < 0: + sq_A = sqrt(A) + t0[0] = 0.5 * (s0 + sq_A) + t1[0] = 0.5 * (s0 - sq_A) + return 2 + + # no real root + else: + return 0 + + # if resolvent cubic solutions have only one real root t0 + else: + if t0[0] <= 0: + s0 = sqrt(-t0[0]) + A = -2.0 * t1[0] - 2.0 * sigma * sqrt(t1[0] * t1[0] + t2[0] * t2[0]) + B = -2.0 * t1[0] + 2.0 * sigma * sqrt(t1[0] * t1[0] + t2[0] * t2[0]) + + # four real roots + if A >= 0 and B >= 0: + sq_A = sqrt(A) + sq_B = sqrt(B) + t0[0] = 0.5 * (s0 + sq_A) + t1[0] = 0.5 * (s0 - sq_A) + t2[0] = 0.5 * (-s0 + sq_B) + t3[0] = 0.5 * (-s0 - sq_B) + return 4 + + # two real roots + elif A < 0 and B >= 0: + sq_B = sqrt(B) + t0[0] = 0.5 * (-s0 + sq_B) + t1[0] = 0.5 * (-s0 - sq_B) + return 2 + + # two real roots + elif A >= 0 and B < 0: + sq_A = sqrt(A) + t0[0] = 0.5 * (s0 + sq_A) + t1[0] = 0.5 * (s0 - sq_A) + return 2 + + # no real root + else: + return 0 + + # no real root if -t0 < 0 + return 0 + + +@cython.cdivision(True) +cdef void one_newton_step(double b, double c, double d, double e, double *x) nogil: + """ + One step Newton's method for monic quartic polinomial: x^4 + b.x^3 + c.x^2 * d.x + e = 0 + + :param double b: Qurtic constant. + :param double c: Qurtic constant. + :param double d: Qurtic constant. + :param double e: Qurtic constant. + :param double x: one root of the quartic. + """ + cdef double fx, dfx + dfx = ((4.0 * x[0] + 3 * b) * x[0] + 2.0 * d) * x[0] + e + if not is_zero(dfx): + fx = (((x[0] + b) * x[0] + c) * x[0] + d) * x[0] + e + x[0] = x[0] - fx / dfx + + +@cython.cdivision(True) +cdef int solve_quartic(double a, double b, double c, double d, double e, + double *t0, double *t1, double *t2, double *t3) nogil: + """ + Calculates the real roots of a quartic equation with Van der Waerden method. + + The a, b, c, d and e arguments are the five constants of the quartic equation: + + f = a.x^4 + b.x^3 + c.x^2 + d.x + e + + The quartic equation has 0, 1, 2, 3 or 4 real roots, and this function will return the number of real roots. + + The values of the real roots, are returned by setting the values of the + memory locations pointed to by t0, t1, t2, and t3. If there are one or two real roots, + the values of t2 and t3 will be undefined. If there is no real root, + all values will be undefined. + + The practical algorithm of Van der Waerden method is followed by https://quarticequations.com + + :param double a: Qurtic constant. + :param double b: Qurtic constant. + :param double c: Qurtic constant. + :param double d: Qurtic constant. + :param double e: Qurtic constant. + :param double t0: 1st root of the quartic. + :param double t1: 2nd root of the quartic. + :param double t2: 3rd root of the quartic. + :param double t3: 4th root of the quartic. + :return: Number of real roots. + :rtype: int + """ + cdef: + double p, q, r, sq_b + int cubic_num, num = 0 + + # normal form: x^4 + bx^3 + cx^2 + dx + e = 0 + b /= a + c /= a + d /= a + e /= a + + # substitute x = y - b / 4 to eliminate quadric term: y^4 + py^2 + qy + r = 0 + sq_b = b * b + p = c - 3 * sq_b / 8.0 + q = sq_b * b / 8.0 - 0.5 * b * c + d + r = -3.0 * sq_b * sq_b / 256.0 + sq_b * c / 16.0 - b * d / 4.0 + e + + if is_zero(r): + # no absolute term: y(y^3 + py + q) = 0 + t0[0] = 0 + cubic_num = solve_cubic(1.0, 0.0, p, q, t1, t2, t3) + num = 1 + cubic_num + + else: + # solve depressed quartic + num = _solve_depressed_quartic(p, q, r, t0, t1, t2, t3) + + # resubstitute + t0[0] -= b / 4.0 + t1[0] -= b / 4.0 + t2[0] -= b / 4.0 + t3[0] -= b / 4.0 + + # One newton step for each real root + if num > 0: + one_newton_step(b, c, d, e, t0) + one_newton_step(b, c, d, e, t1) + if num > 2: + one_newton_step(b, c, d, e, t2) + if num > 3: + one_newton_step(b, c, d, e, t3) + + return num + + @cython.boundscheck(False) @cython.wraparound(False) cdef bint winding2d(double[:,::1] vertices) nogil: @@ -538,3 +866,24 @@ def _test_winding2d(p): def _point_inside_polygon(vertices, ptx, pty): """Expose cython function for testing.""" return point_inside_polygon(vertices, ptx, pty) + +def _solve_cubic(a, b, c, d): + """Expose cython function for testing.""" + t0 = 0.0 + t1 = 0.0 + t2 = 0.0 + num = 0.0 + num = solve_cubic(a, b, c, d, &t0, &t1, &t2) + + return (t0, t1, t2, num) + +def _solve_quartic(a, b, c, d, e): + """Expose cython function for testing.""" + t0 = 0.0 + t1 = 0.0 + t2 = 0.0 + t3 = 0.0 + num = 0.0 + num = solve_quartic(a, b, c, d, e, &t0, &t1, &t2, &t3) + + return (t0, t1, t2, t3, num) diff --git a/raysect/primitive/__init__.pxd b/raysect/primitive/__init__.pxd index 2b236180..cdcf2fae 100644 --- a/raysect/primitive/__init__.pxd +++ b/raysect/primitive/__init__.pxd @@ -29,6 +29,7 @@ from raysect.primitive.box cimport Box from raysect.primitive.sphere cimport Sphere +from raysect.primitive.torus cimport Torus from raysect.primitive.cylinder cimport Cylinder from raysect.primitive.csg cimport Union, Intersect, Subtract from raysect.primitive.mesh cimport Mesh diff --git a/raysect/primitive/__init__.py b/raysect/primitive/__init__.py index fc571c09..fa3c6656 100644 --- a/raysect/primitive/__init__.py +++ b/raysect/primitive/__init__.py @@ -28,6 +28,7 @@ # POSSIBILITY OF SUCH DAMAGE. from .box import Box +from .torus import Torus from .sphere import Sphere from .cylinder import Cylinder from .csg import Union, Intersect, Subtract diff --git a/raysect/primitive/meson.build b/raysect/primitive/meson.build index 6dab2462..1f14a347 100644 --- a/raysect/primitive/meson.build +++ b/raysect/primitive/meson.build @@ -5,8 +5,8 @@ target_path = 'raysect/primitive' # source files py_files = ['__init__.py'] -pyx_files = ['box.pyx', 'cone.pyx', 'csg.pyx', 'cylinder.pyx', 'parabola.pyx', 'sphere.pyx', 'utility.pyx'] -pxd_files = ['__init__.pxd', 'box.pxd', 'cone.pxd', 'csg.pxd', 'cylinder.pxd', 'parabola.pxd', 'sphere.pxd', 'utility.pxd'] +pyx_files = ['box.pyx', 'cone.pyx', 'csg.pyx', 'cylinder.pyx', 'parabola.pyx', 'sphere.pyx', 'torus.pyx', 'utility.pyx'] +pxd_files = ['__init__.pxd', 'box.pxd', 'cone.pxd', 'csg.pxd', 'cylinder.pxd', 'parabola.pxd', 'sphere.pxd', 'torus.pxd', 'utility.pxd'] data_files = [] # compile cython diff --git a/raysect/primitive/torus.pxd b/raysect/primitive/torus.pxd new file mode 100644 index 00000000..9643552b --- /dev/null +++ b/raysect/primitive/torus.pxd @@ -0,0 +1,46 @@ +# cython: language_level=3 + +# Copyright (c) 2014-2021, Dr Alex Meakins, Raysect Project +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in the +# documentation and/or other materials provided with the distribution. +# +# 3. Neither the name of the Raysect Project nor the names of its +# contributors may be used to endorse or promote products derived from +# this software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. + +from raysect.core cimport Primitive, Point3D, Vector3D, Ray, Intersection + +cdef class Torus(Primitive): + + cdef: + double _major_radius, _minor_radius + bint _further_intersection + int _next_t_index + int _num_t + double[4] _cached_t + Point3D _cached_origin + Vector3D _cached_direction + Ray _cached_ray + + cdef Intersection _generate_intersection(self, Ray ray, Point3D origin, Vector3D direction, double ray_distance) diff --git a/raysect/primitive/torus.pyx b/raysect/primitive/torus.pyx new file mode 100644 index 00000000..0b9818aa --- /dev/null +++ b/raysect/primitive/torus.pyx @@ -0,0 +1,371 @@ +# cython: language_level=3 + +# Copyright (c) 2014-2021, Dr Alex Meakins, Raysect Project +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in the +# documentation and/or other materials provided with the distribution. +# +# 3. Neither the name of the Raysect Project nor the names of its +# contributors may be used to endorse or promote products derived from +# this software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. + + +from raysect.core cimport Material, new_intersection, BoundingBox3D, BoundingSphere3D, new_point3d, new_normal3d, Normal3D, AffineMatrix3D +from raysect.core.math.cython cimport solve_quartic, swap_double, sort_three_doubles, sort_four_doubles +from libc.math cimport hypot + + +# bounding box and sphere are padded by small amounts to avoid numerical accuracy issues +DEF BOX_PADDING = 1e-9 +DEF SPHERE_PADDING = 1.000000001 + +# additional ray distance to avoid re-hitting the same surface point +DEF EPSILON = 1e-9 + + +cdef class Torus(Primitive): + """ + A torus primitive. + + The torus is defined by major and minor radius. + The major radius is the distance from the center of the tube to the center of the torus. + The minor radius is the radius of the tube. + The center of the torus corresponds to the origin of the local coordinate system. + The axis of revolution coincides with the z-axis, and The center of the torus tube lies + on the x-y plane. + + :param float major_radius: Major radius of the torus in meters (default = 1.0). + :param float minor_radius: Minor radius of the torus in meters (default = 0.5). + :param Node parent: Scene-graph parent node or None (default = None). + :param AffineMatrix3D transform: An AffineMatrix3D defining the local coordinate system relative to the scene-graph parent (default = identity matrix). + :param Material material: A Material object defining the torus's material (default = None). + :param str name: A string specifying a user-friendly name for the torus (default = ""). + + :ivar float major_radius: The major radius of the torus in meters. + :ivar float minor_radius: The minor radius of the torus in meters. + + .. code-block:: pycon + + >>> from raysect.core import translate + >>> from raysect.primitive import Torus + >>> from raysect.optical import World + >>> from raysect.optical.material import UniformSurfaceEmitter + >>> from raysect.optical.library.spectra.colours import orange + >>> + >>> world = World() + >>> + >>> torus = Torus(1.0, 0.5, parent=world, transform=translate(3, 0, 0), + material=UniformSurfaceEmitter(orange), name="orange torus") + """ + + def __init__(self, double major_radius=1.0, double minor_radius=0.5, object parent=None, + AffineMatrix3D transform=None, Material material=None, str name=None): + + super().__init__(parent, transform, material, name) + + if major_radius < minor_radius or minor_radius < 0.0: + raise ValueError("Torus minor radius cannot be less than zero and greater than major radius.") + + self._major_radius = major_radius + self._minor_radius = minor_radius + + # initialise next intersection caching and control attributes + self._further_intersection = False + self._next_t_index = 0 + self._num_t = 0 + self._cached_origin = None + self._cached_direction = None + self._cached_ray = None + + @property + def major_radius(self): + """ + The major radius of this torus. + + :rtype: float + """ + return self._major_radius + + @major_radius.setter + def major_radius(self, double major_radius): + + # don't do anything if the value is unchanged + if major_radius == self._major_radius: + return + + if major_radius < 0.0: + raise ValueError("Torus major radius cannot be less than zero.") + if major_radius < self._minor_radius: + raise ValueError("Torus major radius cannot be less than minor radius.") + self._major_radius = major_radius + + # the next intersection cache has been invalidated by the major radius change + self._further_intersection = False + + # any geometry caching in the root node is now invalid, inform root + self.notify_geometry_change() + + @property + def minor_radius(self): + """ + The minor radius of this torus. + + :rtype: float + """ + return self._minor_radius + + @minor_radius.setter + def minor_radius(self, double minor_radius): + + # don't do anything if the value is unchanged + if minor_radius == self._minor_radius: + return + + if minor_radius < 0.0: + raise ValueError("Torus minor radius cannot be less than zero.") + if minor_radius > self._major_radius: + raise ValueError("Torus minor radius cannot be greater than major radius.") + self._minor_radius = minor_radius + + # the next intersection cache has been invalidated by the minor radius change + self._further_intersection = False + + # any geometry caching in the root node is now invalid, inform root + self.notify_geometry_change() + + cpdef Intersection hit(self, Ray ray): + + cdef: + Point3D origin + Vector3D direction + double sq_origin_xy, sq_direction_xy, sq_origin, sq_direction + double origin_direction_xy, origin_dot_direction, sq_r, sq_R, R2_r2 + double a, b, c, d, e, t_closest + double[4] t + int num, i + + # reset further intersection state + self._further_intersection = False + + # convert ray parameters to local space + origin = ray.origin.transform(self.to_local()) + direction = ray.direction.transform(self.to_local()) + + # calculate temporary values + sq_origin_xy = origin.x * origin.x + origin.y * origin.y + sq_direction_xy = direction.x * direction.x + direction.y * direction.y + + sq_origin = sq_origin_xy + origin.z * origin.z + sq_direction = sq_direction_xy + direction.z * direction.z + + origin_direction_xy = origin.x * direction.x + origin.y * direction.y + origin_dot_direction = origin_direction_xy + origin.z * direction.z + + sq_r = self._minor_radius * self._minor_radius + sq_R = self._major_radius * self._major_radius + R2_r2 = sq_R - sq_r + + # coefficients of quartic equation + a = sq_direction * sq_direction + b = 4.0 * sq_direction * origin_dot_direction + c = 2.0 * (2.0 * origin_dot_direction * origin_dot_direction + sq_direction * (sq_origin + R2_r2)) - 4.0 * sq_R * sq_direction_xy + d = 4.0 * origin_dot_direction * (sq_origin + R2_r2) - 8.0 * sq_R * origin_direction_xy + e = (sq_origin + R2_r2) * (sq_origin + R2_r2) - 4.0 * sq_R * sq_origin_xy + + # calculate intersection distances by solving the quartic equation + # ray misses if there are no real roots of the quartic + num = solve_quartic(a, b, c, d, e, &t[0], &t[1], &t[2], &t[3]) + + if num == 0: + return None + + elif num == 1: + # test the intersection points inside the ray search range [0, max_distance] + if t[0] > ray.max_distance or t[0] < 0.0: + return None + else: + t_closest = t[0] + + else: + # sorting solutions in each number of them + if num == 2: + # ensure t0 < t1 + if t[0] > t[1]: + swap_double(&t[0], &t[1]) + + # substitute the last value into undefined variables + t[2] = t[1] + t[3] = t[1] + + elif num == 3: + # ensure t0 < t1 < t2 + sort_three_doubles(&t[0], &t[1], &t[2]) + + # substitute the last value into undefined variables + t[3] = t[2] + + elif num == 4: + # ensure t0 < t1 < t2 < t3 + sort_four_doubles(&t[0], &t[1], &t[2], &t[3]) + else: + return None + + # test the intersection points inside the ray search range [0, max_distance] + if t[0] > ray.max_distance or t[3] < 0.0: + return None + + # cache the all intersection points + self._num_t = num + self._cached_t = t + + for i in range(num - 1): + if t[i] >= 0.0: + t_closest = t[i] + if t[i + 1] <= ray.max_distance: + self._further_intersection = True + self._next_t_index = i + 1 + self._cached_ray = ray + self._cached_origin = origin + self._cached_direction = direction + + return self._generate_intersection(ray, origin, direction, t_closest) + + if t[num - 1] <= ray.max_distance: + t_closest = t[num - 1] + else: + return None + + return self._generate_intersection(ray, origin, direction, t_closest) + + cpdef Intersection next_intersection(self): + + cdef: + double next_t + Intersection new_intersection + + if not self._further_intersection: + return None + + next_t = self._cached_t[self._next_t_index] + + # generate the next intersection. + new_intersection = self._generate_intersection( + self._cached_ray, + self._cached_origin, + self._cached_direction, + next_t, + ) + + # update the next intersection index + # if there are no more intersections, disable further intersection state + # and reset the next_t_index to 0 + if self._next_t_index < self._num_t - 1: + self._next_t_index += 1 + else: + self._further_intersection = False + self._next_t_index = 0 + + return new_intersection + + cdef Intersection _generate_intersection(self, Ray ray, Point3D origin, Vector3D direction, double ray_distance): + + cdef Point3D hit_point, inside_point, outside_point + cdef Normal3D normal + cdef double alpha, delta_x, delta_y, delta_z + cdef bint exiting + + # point of surface intersection in local space + hit_point = new_point3d(origin.x + ray_distance * direction.x, + origin.y + ray_distance * direction.y, + origin.z + ray_distance * direction.z) + + # normal is normalised vector from torus tube centre to hit_point + alpha = self._major_radius / hypot(hit_point.x, hit_point.y) + normal = new_normal3d((1.0 - alpha) * hit_point.x, (1.0 - alpha) * hit_point.y, hit_point.z) + normal = normal.normalise() + + # calculate points inside and outside of surface for daughter rays to + # spawn from - these points are displaced from the surface to avoid + # re-hitting the same surface + delta_x = EPSILON * normal.x + delta_y = EPSILON * normal.y + delta_z = EPSILON * normal.z + + inside_point = new_point3d(hit_point.x - delta_x, hit_point.y - delta_y, hit_point.z - delta_z) + outside_point = new_point3d(hit_point.x + delta_x, hit_point.y + delta_y, hit_point.z + delta_z) + + # is ray exiting surface + exiting = direction.dot(normal) >= 0.0 + + return new_intersection(ray, ray_distance, self, hit_point, inside_point, outside_point, + normal, exiting, self.to_local(), self.to_root()) + + cpdef bint contains(self, Point3D point) except -1: + + cdef Point3D local_point + cdef double discriminant, distance_xy, distance_sqr, sq_R, R2_r2 + + # convert world space point to local space + local_point = point.transform(self.to_local()) + + # calculate the interior discriminant + distance_xy = local_point.x * local_point.x + local_point.y * local_point.y + distance_sqr = distance_xy + local_point.z * local_point.z + sq_R = self._major_radius * self._major_radius + R2_r2 = sq_R - self._minor_radius * self._minor_radius + discriminant = distance_sqr * distance_sqr + 2.0 * distance_sqr * R2_r2 + R2_r2 * R2_r2 - 4.0 * sq_R * distance_xy + + # point is outside torus if discriminant is greater than 0 + return discriminant <= 0.0 + + cpdef BoundingBox3D bounding_box(self): + + cdef: + double extent + list points + Point3D point + BoundingBox3D box + + box = BoundingBox3D() + + # calculate local bounds + extent = self._major_radius + self._minor_radius + BOX_PADDING + box.lower = new_point3d(-extent, -extent, -self._minor_radius - BOX_PADDING) + box.upper = new_point3d(extent, extent, self._minor_radius + BOX_PADDING) + + # obtain local space vertices + points = box.vertices() + + # convert points to world space + box = BoundingBox3D() + for point in points: + box.extend(point.transform(self.to_root())) + + return box + + cpdef BoundingSphere3D bounding_sphere(self): + cdef Point3D centre = new_point3d(0, 0, 0).transform(self.to_root()) + return BoundingSphere3D(centre, (self._major_radius + self._minor_radius) * SPHERE_PADDING) + + cpdef object instance(self, object parent=None, AffineMatrix3D transform=None, Material material=None, str name=None): + return Torus(self._major_radius, self._minor_radius, parent, transform, material, name) \ No newline at end of file