diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index bbeae168..af832884 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -12,7 +12,7 @@ jobs: fail-fast: false matrix: numpy-version: ["numpy"] - python-version: ["3.9", "3.10", "3.11", "3.12"] + python-version: ["3.9", "3.10", "3.11", "3.12", "3.13"] steps: - name: Checkout code uses: actions/checkout@v4 @@ -21,7 +21,7 @@ jobs: with: python-version: ${{ matrix.python-version }} - name: Install Python dependencies - run: python -m pip install --prefer-binary meson-python meson ninja setuptools "cython>=0.28,<3.0" "matplotlib>=3,<4" ${{ matrix.numpy-version }} + run: python -m pip install --prefer-binary meson-python meson ninja setuptools "cython>=3.1" "matplotlib>=3,<4" ${{ matrix.numpy-version }} - name: Build and install Raysect run: dev/install_editable.sh - name: Run tests diff --git a/CHANGELOG.txt b/CHANGELOG.txt index 7dd3faf6..afb6b336 100644 --- a/CHANGELOG.txt +++ b/CHANGELOG.txt @@ -13,6 +13,9 @@ Build changes: - Meson-python performs the build out of the project folder, so the project source folders will no longer be polluted with build artefacts. - Cython build annotations are now always enabled, the annotation files can be found in the build folder under the build artefacts folder associated with each so file (*.so.p folder). +API changes: +* Corrected spelling of all classes, methods and functions where "targeted" was incorrectly spelt "targetted". + Release 0.8.1 (12 Feb 2023) --------------------------- diff --git a/demos/maths/plot_targetted_sampler.py b/demos/maths/plot_targeted_sampler.py similarity index 86% rename from demos/maths/plot_targetted_sampler.py rename to demos/maths/plot_targeted_sampler.py index 0978a647..d46f2a7b 100644 --- a/demos/maths/plot_targetted_sampler.py +++ b/demos/maths/plot_targeted_sampler.py @@ -2,7 +2,7 @@ import matplotlib.pyplot as plt from math import atan2, sqrt, degrees -from raysect.core.math import Point3D, TargettedHemisphereSampler, TargettedSphereSampler +from raysect.core.math import Point3D, TargetedHemisphereSampler, TargetedSphereSampler def display_samples(samples, title): @@ -34,8 +34,8 @@ def display_samples(samples, title): observation_point = Point3D(0, 0, 0) # generate samplers -hemisphere = TargettedHemisphereSampler(targets) -sphere = TargettedSphereSampler(targets) +hemisphere = TargetedHemisphereSampler(targets) +sphere = TargetedSphereSampler(targets) # sample for origin point and point at (0, 0, -10) h_samples = hemisphere(observation_point, samples=samples) diff --git a/demos/observers/cornell_box_cooke_triplet.py b/demos/observers/cornell_box_cooke_triplet.py index 2f7ee055..653632f2 100644 --- a/demos/observers/cornell_box_cooke_triplet.py +++ b/demos/observers/cornell_box_cooke_triplet.py @@ -4,7 +4,7 @@ from raysect.optical import World, Node, translate, rotate, Point3D from raysect.optical.material import Lambert, UniformSurfaceEmitter, AbsorbingSurface, Checkerboard from raysect.optical.library import * -from raysect.optical.observer import TargettedCCDArray +from raysect.optical.observer import TargetedCCDArray from raysect.optical.observer import RGBPipeline2D, BayerPipeline2D, PowerPipeline2D from raysect.optical.observer import RGBAdaptiveSampler2D from raysect.core.math import mm @@ -184,9 +184,9 @@ rgb = RGBPipeline2D(display_unsaturated_fraction=0.96, name="sRGB") sampler = RGBAdaptiveSampler2D(rgb, ratio=10, fraction=0.2, min_samples=1000, cutoff=0.01) -# CCD targetting all rays at last lens element for speed -ccd = TargettedCCDArray( - targetted_path_prob=1.0, targets=[l3], +# CCD targeting all rays at last lens element for speed +ccd = TargetedCCDArray( + targeted_path_prob=1.0, targets=[l3], width=mm(35), pixels=(512, 512), parent=image_plane, transform=translate(0, 0, 0)*rotate(0, 0, 180), pipelines=[rgb] diff --git a/demos/observers/cornell_box_real_pinhole.py b/demos/observers/cornell_box_real_pinhole.py index 64e43e38..d456cc13 100644 --- a/demos/observers/cornell_box_real_pinhole.py +++ b/demos/observers/cornell_box_real_pinhole.py @@ -5,7 +5,7 @@ from raysect.optical import World, Node, translate, rotate, Point3D from raysect.optical.material import Lambert, UniformSurfaceEmitter, NullMaterial from raysect.optical.library import * -from raysect.optical.observer import RGBPipeline2D, BayerPipeline2D, PowerPipeline2D, TargettedCCDArray +from raysect.optical.observer import RGBPipeline2D, BayerPipeline2D, PowerPipeline2D, TargetedCCDArray from raysect.optical.observer import RGBAdaptiveSampler2D @@ -144,7 +144,7 @@ camera = Node(parent=world, transform=translate(0, 0, -3.3)) pinhole = Sphere(0.0005, camera, transform=translate(0, 0, 0), material=NullMaterial()) -film = TargettedCCDArray(targetted_path_prob=1.0, targets=[pinhole], width=0.1, pixels=(512, 512), parent=camera, +film = TargetedCCDArray(targeted_path_prob=1.0, targets=[pinhole], width=0.1, pixels=(512, 512), parent=camera, transform=translate(0, 0, -0.1207), pipelines=pipelines) film.frame_sampler = sampler film.pixel_samples = 250 diff --git a/demos/observers/metal_with_lens.py b/demos/observers/metal_with_lens.py index e1bbaa8c..a5bcc6ce 100644 --- a/demos/observers/metal_with_lens.py +++ b/demos/observers/metal_with_lens.py @@ -7,7 +7,7 @@ from raysect.optical.library.metal import Gold, Silver, Copper, Titanium, Aluminium, Beryllium from raysect.optical.material import Lambert, UniformSurfaceEmitter, AbsorbingSurface, NullMaterial from raysect.optical.library import schott -from raysect.optical.observer import RGBPipeline2D, BayerPipeline2D, TargettedCCDArray, CCDArray, RGBAdaptiveSampler2D +from raysect.optical.observer import RGBPipeline2D, BayerPipeline2D, TargetedCCDArray, CCDArray, RGBAdaptiveSampler2D from raysect.optical.colour import ciexyz_x, ciexyz_y, ciexyz_z @@ -49,7 +49,7 @@ pipelines = [rgb, bayer] # ccd = CCDArray(parent=camera, pipelines=pipelines) -ccd = TargettedCCDArray(targets=[aperture], parent=camera, pipelines=pipelines) +ccd = TargetedCCDArray(targets=[aperture], parent=camera, pipelines=pipelines) ccd.frame_sampler = sampler ccd.pixels = (180*2, 120*2) # (360, 240) ccd.pixel_samples = 250 diff --git a/demos/observers/targetted_pixel.py b/demos/observers/targeted_pixel.py similarity index 76% rename from demos/observers/targetted_pixel.py rename to demos/observers/targeted_pixel.py index ff272219..8ca58ed1 100644 --- a/demos/observers/targetted_pixel.py +++ b/demos/observers/targeted_pixel.py @@ -3,7 +3,7 @@ from raysect.primitive import Box from raysect.optical import World, translate, Point3D, Node -from raysect.optical.observer import Pixel, TargettedPixel, PowerPipeline0D +from raysect.optical.observer import Pixel, TargetedPixel, PowerPipeline0D from raysect.optical.material import UnitySurfaceEmitter @@ -28,16 +28,16 @@ basic_pipeline = PowerPipeline0D(name="Basic Pixel Observer") basic_pixel = Pixel(parent=world, pixel_samples=SAMPLES, pipelines=[basic_pipeline]) -# setup targetted pixel -targetted_pipeline = PowerPipeline0D(name="Targeted Pixel Observer") -targetted_pixel = TargettedPixel(parent=world, targets=targets, pixel_samples=SAMPLES, pipelines=[targetted_pipeline]) -targetted_pixel.targetted_path_prob = 1 +# setup targeted pixel +targeted_pipeline = PowerPipeline0D(name="Targeted Pixel Observer") +targeted_pixel = TargetedPixel(parent=world, targets=targets, pixel_samples=SAMPLES, pipelines=[targeted_pipeline]) +targeted_pixel.targeted_path_prob = 1 # render ion() basic_pixel.observe() print() -targetted_pixel.observe() +targeted_pixel.observe() ioff() show() diff --git a/demos/optics/etendue_of_pinhole.py b/demos/optics/etendue_of_pinhole.py index 38a85b9b..e1daac54 100644 --- a/demos/optics/etendue_of_pinhole.py +++ b/demos/optics/etendue_of_pinhole.py @@ -3,7 +3,7 @@ import matplotlib.pyplot as plt from raysect.core import Point3D, Vector3D, rotate_basis, translate, Ray as CoreRay -from raysect.core.math.sampler import DiskSampler3D, RectangleSampler3D, TargettedHemisphereSampler +from raysect.core.math.sampler import DiskSampler3D, RectangleSampler3D, TargetedHemisphereSampler from raysect.optical import World from raysect.primitive import Box, Cylinder, Subtract from raysect.optical.material import AbsorbingSurface, NullMaterial @@ -35,8 +35,8 @@ def raytraced_etendue(distance, detector_radius=0.001, ray_count=100000, batches sphere = target.bounding_sphere() spheres = [(sphere.centre.transform(detector_transform), sphere.radius, 1.0)] - # instance targetted pixel sampler - targetted_sampler = TargettedHemisphereSampler(spheres) + # instance targeted pixel sampler + targeted_sampler = TargetedHemisphereSampler(spheres) point_sampler = DiskSampler3D(detector_radius) @@ -53,8 +53,8 @@ def raytraced_etendue(distance, detector_radius=0.001, ray_count=100000, batches passed = 0.0 for origin in origins: - # obtain targetted vector sample - direction, pdf = targetted_sampler(origin, pdf=True) + # obtain targeted vector sample + direction, pdf = targeted_sampler(origin, pdf=True) path_weight = R_2_PI * direction.z/pdf origin = origin.transform(detector_transform) diff --git a/dev/root-meson.build b/dev/root-meson.build index 526acc1a..98253b94 100644 --- a/dev/root-meson.build +++ b/dev/root-meson.build @@ -4,5 +4,6 @@ py = import('python').find_installation(pure: false) numpy = dependency('numpy', method: 'config-tool') fs = import('fs') -cython_args = ['--annotate'] +# disabling explicit noexcept until pycharm fixes their cython 3 support (causes a large number of build warnings) +cython_args = ['--annotate', '-X legacy_implicit_noexcept=True'] cython_dependencies = [numpy] diff --git a/meson.build b/meson.build index 072eed68..3a2b2cf0 100644 --- a/meson.build +++ b/meson.build @@ -7,7 +7,8 @@ py = import('python').find_installation(pure: false) numpy = dependency('numpy', method: 'config-tool') fs = import('fs') -cython_args = ['--annotate'] +# disabling explicit noexcept until pycharm fixes their cython 3 support (causes a large number of build warnings) +cython_args = ['--annotate', '-X legacy_implicit_noexcept=True'] cython_dependencies = [numpy] subdir('raysect') diff --git a/pyproject.toml b/pyproject.toml index a66e34b9..b16903e1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,5 +27,5 @@ Issues = "https://github.com/raysect/source/issues" Changelog = "https://github.com/raysect/source/blob/master/CHANGELOG.txt" [build-system] -requires = ["meson-python", "setuptools", "wheel", "numpy", "cython<3.0"] +requires = ["meson-python", "setuptools", "wheel", "numpy", "cython>=3.1"] build-backend = "mesonpy" diff --git a/raysect/core/acceleration/accelerator.pxd b/raysect/core/acceleration/accelerator.pxd index ca008e6b..495d3051 100644 --- a/raysect/core/acceleration/accelerator.pxd +++ b/raysect/core/acceleration/accelerator.pxd @@ -33,10 +33,9 @@ from raysect.core.ray cimport Ray from raysect.core.math cimport Point3D from raysect.core.intersection cimport Intersection + cdef class Accelerator: cpdef build(self, list primitives) - cpdef Intersection hit(self, Ray ray) - cpdef list contains(self, Point3D point) diff --git a/raysect/core/acceleration/accelerator.pyx b/raysect/core/acceleration/accelerator.pyx index 3ff5e3f4..e6279d91 100644 --- a/raysect/core/acceleration/accelerator.pyx +++ b/raysect/core/acceleration/accelerator.pyx @@ -32,13 +32,10 @@ cdef class Accelerator: cpdef build(self, list primitives): - pass cpdef Intersection hit(self, Ray ray): - raise NotImplementedError("Accelerator virtual method hit() has not been implemented.") cpdef list contains(self, Point3D point): - raise NotImplementedError("Accelerator virtual method contains() has not been implemented.") \ No newline at end of file diff --git a/raysect/core/acceleration/boundprimitive.pxd b/raysect/core/acceleration/boundprimitive.pxd index 8a1e35ed..cec01939 100644 --- a/raysect/core/acceleration/boundprimitive.pxd +++ b/raysect/core/acceleration/boundprimitive.pxd @@ -37,12 +37,11 @@ from raysect.core.intersection cimport Intersection cdef class BoundPrimitive: - cdef readonly Primitive primitive - cdef readonly BoundingBox3D box - cdef bint _primitive_tested + cdef: + readonly Primitive primitive + readonly BoundingBox3D box + bint _primitive_tested cdef Intersection hit(self, Ray ray) - cdef Intersection next_intersection(self) - cdef bint contains(self, Point3D point) \ No newline at end of file diff --git a/raysect/core/acceleration/kdtree.pxd b/raysect/core/acceleration/kdtree.pxd index 6616e2ef..16d0612d 100644 --- a/raysect/core/acceleration/kdtree.pxd +++ b/raysect/core/acceleration/kdtree.pxd @@ -33,12 +33,15 @@ from raysect.core.acceleration.accelerator cimport Accelerator as _Accelerator from raysect.core.math.spatial.kdtree3d cimport KDTree3DCore as _KDTreeCore from raysect.core.intersection cimport Intersection + cdef class _PrimitiveKDTree(_KDTreeCore): + cdef: list primitives Intersection hit_intersection cdef class KDTree(_Accelerator): + cdef _PrimitiveKDTree _kdtree diff --git a/raysect/core/acceleration/kdtree.pyx b/raysect/core/acceleration/kdtree.pyx index b2900780..2d7b28df 100644 --- a/raysect/core/acceleration/kdtree.pyx +++ b/raysect/core/acceleration/kdtree.pyx @@ -154,7 +154,6 @@ cdef class _PrimitiveKDTree(_KDTreeCore): # dereference the primitives and check if they contain the point enclosing_primitives = [] for item in range(count): - index = self._nodes[id].items[item] primitive = self.primitives[index] if primitive.contains(point): diff --git a/raysect/core/acceleration/unaccelerated.pxd b/raysect/core/acceleration/unaccelerated.pxd index 3d0149a5..d9ce0979 100644 --- a/raysect/core/acceleration/unaccelerated.pxd +++ b/raysect/core/acceleration/unaccelerated.pxd @@ -35,5 +35,6 @@ from raysect.core.boundingbox cimport BoundingBox3D cdef class Unaccelerated(Accelerator): - cdef list primitives - cdef BoundingBox3D world_box + cdef: + list primitives + BoundingBox3D world_box diff --git a/raysect/core/acceleration/unaccelerated.pyx b/raysect/core/acceleration/unaccelerated.pyx index 8695c549..c07e407f 100644 --- a/raysect/core/acceleration/unaccelerated.pyx +++ b/raysect/core/acceleration/unaccelerated.pyx @@ -37,6 +37,7 @@ from raysect.core.math cimport Point3D from raysect.core.intersection cimport Intersection from raysect.core.acceleration.boundprimitive cimport BoundPrimitive + cdef class Unaccelerated(Accelerator): def __init__(self): @@ -54,7 +55,6 @@ cdef class Unaccelerated(Accelerator): self.world_box = BoundingBox3D() for primitive in primitives: - accel_primitive = BoundPrimitive(primitive) self.primitives.append(accel_primitive) self.world_box.union(accel_primitive.box) @@ -70,23 +70,17 @@ cdef class Unaccelerated(Accelerator): # does the ray intersect the space containing the primitives if not self.world_box.hit(ray): - return None # find the closest primitive-ray intersection closest_intersection = None - # intial search distance is maximum possible ray extent + # initial search distance is maximum possible ray extent distance = ray.max_distance - for primitive in self.primitives: - intersection = primitive.hit(ray) - if intersection is not None: - if intersection.ray_distance < distance: - distance = intersection.ray_distance closest_intersection = intersection @@ -101,15 +95,12 @@ cdef class Unaccelerated(Accelerator): BoundPrimitive primitive if not self.world_box.contains(point): - return [] enclosing_primitives = [] for primitive in self.primitives: - if primitive.contains(point): - enclosing_primitives.append(primitive.primitive) return enclosing_primitives \ No newline at end of file diff --git a/raysect/core/boundingbox.pxd b/raysect/core/boundingbox.pxd index cd1e59b4..29167b96 100644 --- a/raysect/core/boundingbox.pxd +++ b/raysect/core/boundingbox.pxd @@ -40,37 +40,21 @@ cdef class BoundingBox3D: cdef Point3D upper cdef Point3D get_centre(self) - cpdef bint hit(self, Ray ray) - cpdef tuple full_intersection(self, Ray ray) - cdef bint intersect(self, Ray ray, double *front_intersection, double *back_intersection) - cdef void _slab(self, double origin, double direction, double lower, double upper, double *front_intersection, double *back_intersection) nogil - cpdef bint contains(self, Point3D point) - cpdef object union(self, BoundingBox3D box) - cpdef object extend(self, Point3D point, double padding=*) - cpdef double surface_area(self) - cpdef double volume(self) - cpdef list vertices(self) - cpdef double extent(self, int axis) except -1 - cpdef int largest_axis(self) - cpdef double largest_extent(self) - cpdef object pad(self, double padding) - cpdef object pad_axis(self, int axis, double padding) - cpdef BoundingSphere3D enclosing_sphere(self) @@ -95,23 +79,14 @@ cdef class BoundingBox2D: cdef Point2D upper cpdef bint contains(self, Point2D point) - cpdef object union(self, BoundingBox2D box) - cpdef object extend(self, Point2D point, double padding=*) - cpdef double surface_area(self) - cpdef list vertices(self) - cpdef double extent(self, int axis) except -1 - cpdef int largest_axis(self) - cpdef double largest_extent(self) - cpdef object pad(self, double padding) - cpdef object pad_axis(self, int axis, double padding) diff --git a/raysect/core/boundingbox.pyx b/raysect/core/boundingbox.pyx index a86f7fdb..5c54fc41 100644 --- a/raysect/core/boundingbox.pyx +++ b/raysect/core/boundingbox.pyx @@ -32,18 +32,19 @@ # TODO: add docstrings cimport cython +from libc.math cimport INFINITY from raysect.core.math cimport new_point3d, new_point2d -# cython doesn't have a built-in infinity constant, this compiles to +infinity -DEF INFINITY = 1e999 # axis defines -DEF X_AXIS = 0 -DEF Y_AXIS = 1 -DEF Z_AXIS = 2 +cdef enum: + X_AXIS = 0 + Y_AXIS = 1 + Z_AXIS = 2 + # defines the padding on the sphere which encloses the BoundingBox3D. -DEF SPHERE_PADDING = 1.000001 +cdef const double SPHERE_PADDING = 1.000001 @cython.freelist(256) @@ -258,6 +259,7 @@ cdef class BoundingBox3D: return False if (point.z < self.lower.z) or (point.z > self.upper.z): return False + return True cpdef object union(self, BoundingBox3D box): diff --git a/raysect/core/boundingsphere.pxd b/raysect/core/boundingsphere.pxd index 0026b721..821bea48 100644 --- a/raysect/core/boundingsphere.pxd +++ b/raysect/core/boundingsphere.pxd @@ -39,19 +39,11 @@ cdef class BoundingSphere3D: Point3D centre cpdef bint hit(self, Ray ray) - cpdef tuple full_intersection(self, Ray ray) - cdef bint intersect(self, Ray ray, double *front_intersection, double *back_intersection) - cpdef bint contains(self, Point3D point) - cpdef object union(self, BoundingSphere3D sphere) - cpdef object extend(self, Point3D point, double padding=*) - cpdef double surface_area(self) - cpdef double volume(self) - cpdef object pad(self, double padding) diff --git a/raysect/core/containers.pxd b/raysect/core/containers.pxd index add9dded..3de638c5 100644 --- a/raysect/core/containers.pxd +++ b/raysect/core/containers.pxd @@ -46,22 +46,16 @@ cdef class LinkedList: _Item last cpdef bint is_empty(self) - cpdef add(self, object value) - cpdef add_items(self, object iterable) - cpdef object get_index(self, int index) - cpdef insert(self, object value, int index) - cpdef object remove(self, int index) cdef class Stack(LinkedList): cpdef push(self, object value) - cpdef object pop(self) diff --git a/raysect/core/intersection.pyx b/raysect/core/intersection.pyx index 48ec1f7f..70cad491 100644 --- a/raysect/core/intersection.pyx +++ b/raysect/core/intersection.pyx @@ -31,6 +31,7 @@ cimport cython + @cython.freelist(256) cdef class Intersection: """ diff --git a/raysect/core/math/_mat4.pxd b/raysect/core/math/_mat4.pxd index 9993668c..40d34ff1 100644 --- a/raysect/core/math/_mat4.pxd +++ b/raysect/core/math/_mat4.pxd @@ -34,10 +34,7 @@ cdef class _Mat4: cdef double m[4][4] cdef double get_element(self, int row, int column) - cdef void set_element(self, int row, int column, double v) - cpdef bint is_identity(self, double tolerance=*) - cpdef bint is_close(self, _Mat4 other, double tolerance=*) diff --git a/raysect/core/math/_mat4.pyx b/raysect/core/math/_mat4.pyx index f8828744..4c5c775b 100644 --- a/raysect/core/math/_mat4.pyx +++ b/raysect/core/math/_mat4.pyx @@ -54,7 +54,6 @@ cdef class _Mat4: # special handling for _Mat4 if isinstance(v, _Mat4): - m = <_Mat4>v for i in range(0, 4): for j in range(0, 4): @@ -95,7 +94,6 @@ cdef class _Mat4: Expects a tuple (row, column) as the index. e.g. v = matrix[1, 2] - """ cdef int row, column diff --git a/raysect/core/math/_vec3.pxd b/raysect/core/math/_vec3.pxd index 70b913c7..df241982 100644 --- a/raysect/core/math/_vec3.pxd +++ b/raysect/core/math/_vec3.pxd @@ -32,16 +32,10 @@ cdef class _Vec3: cdef public double x, y, z - cpdef double dot(self, _Vec3 v) - cpdef double angle(self, _Vec3 v) - cdef double get_length(self) nogil - cdef object set_length(self, double v) - cdef double get_index(self, int index) nogil - cdef void set_index(self, int index, double value) nogil diff --git a/raysect/core/math/_vec3.pyx b/raysect/core/math/_vec3.pyx index 0a1e3d88..02c16bd4 100644 --- a/raysect/core/math/_vec3.pyx +++ b/raysect/core/math/_vec3.pyx @@ -58,8 +58,8 @@ cdef class _Vec3: return self.y elif i == 2: return self.z - else: - raise IndexError("Index out of range [0, 2].") + + raise IndexError("Index out of range [0, 2].") def __setitem__(self, int i, double value): """Sets the vector coordinates by index ([0,1,2] -> [x,y,z]).""" @@ -74,6 +74,7 @@ cdef class _Vec3: raise IndexError("Index out of range [0, 2].") def __iter__(self): + yield self.x yield self.y yield self.z @@ -189,8 +190,8 @@ cdef class _Vec3: return self.y elif index == 2: return self.z - else: - return NAN + + return NAN cdef void set_index(self, int index, double value) nogil: """ diff --git a/raysect/core/math/affinematrix.pxd b/raysect/core/math/affinematrix.pxd index a5de0d98..abbfe35d 100644 --- a/raysect/core/math/affinematrix.pxd +++ b/raysect/core/math/affinematrix.pxd @@ -35,7 +35,6 @@ from raysect.core.math._mat4 cimport _Mat4 cdef class AffineMatrix3D(_Mat4): cpdef AffineMatrix3D inverse(self) - cdef AffineMatrix3D mul(self, AffineMatrix3D m) diff --git a/raysect/core/math/affinematrix.pyx b/raysect/core/math/affinematrix.pyx index 4a67222c..f585d199 100644 --- a/raysect/core/math/affinematrix.pyx +++ b/raysect/core/math/affinematrix.pyx @@ -33,7 +33,6 @@ cimport cython from libc.math cimport fabs - cdef class AffineMatrix3D(_Mat4): """A 4x4 affine matrix. @@ -133,7 +132,7 @@ cdef class AffineMatrix3D(_Mat4): s += ", " return s + "])" - def __mul__(object x, object y): + def __mul__(self, object y): """Multiplication operator. >>> from raysect.core import translate, rotate_x @@ -146,26 +145,28 @@ cdef class AffineMatrix3D(_Mat4): cdef AffineMatrix3D mx, my - if isinstance(x, AffineMatrix3D) and isinstance(y, AffineMatrix3D): + if isinstance(y, AffineMatrix3D): - mx = x + mx = self my = y - return new_affinematrix3d(mx.m[0][0] * my.m[0][0] + mx.m[0][1] * my.m[1][0] + mx.m[0][2] * my.m[2][0] + mx.m[0][3] * my.m[3][0], - mx.m[0][0] * my.m[0][1] + mx.m[0][1] * my.m[1][1] + mx.m[0][2] * my.m[2][1] + mx.m[0][3] * my.m[3][1], - mx.m[0][0] * my.m[0][2] + mx.m[0][1] * my.m[1][2] + mx.m[0][2] * my.m[2][2] + mx.m[0][3] * my.m[3][2], - mx.m[0][0] * my.m[0][3] + mx.m[0][1] * my.m[1][3] + mx.m[0][2] * my.m[2][3] + mx.m[0][3] * my.m[3][3], - mx.m[1][0] * my.m[0][0] + mx.m[1][1] * my.m[1][0] + mx.m[1][2] * my.m[2][0] + mx.m[1][3] * my.m[3][0], - mx.m[1][0] * my.m[0][1] + mx.m[1][1] * my.m[1][1] + mx.m[1][2] * my.m[2][1] + mx.m[1][3] * my.m[3][1], - mx.m[1][0] * my.m[0][2] + mx.m[1][1] * my.m[1][2] + mx.m[1][2] * my.m[2][2] + mx.m[1][3] * my.m[3][2], - mx.m[1][0] * my.m[0][3] + mx.m[1][1] * my.m[1][3] + mx.m[1][2] * my.m[2][3] + mx.m[1][3] * my.m[3][3], - mx.m[2][0] * my.m[0][0] + mx.m[2][1] * my.m[1][0] + mx.m[2][2] * my.m[2][0] + mx.m[2][3] * my.m[3][0], - mx.m[2][0] * my.m[0][1] + mx.m[2][1] * my.m[1][1] + mx.m[2][2] * my.m[2][1] + mx.m[2][3] * my.m[3][1], - mx.m[2][0] * my.m[0][2] + mx.m[2][1] * my.m[1][2] + mx.m[2][2] * my.m[2][2] + mx.m[2][3] * my.m[3][2], - mx.m[2][0] * my.m[0][3] + mx.m[2][1] * my.m[1][3] + mx.m[2][2] * my.m[2][3] + mx.m[2][3] * my.m[3][3], - mx.m[3][0] * my.m[0][0] + mx.m[3][1] * my.m[1][0] + mx.m[3][2] * my.m[2][0] + mx.m[3][3] * my.m[3][0], - mx.m[3][0] * my.m[0][1] + mx.m[3][1] * my.m[1][1] + mx.m[3][2] * my.m[2][1] + mx.m[3][3] * my.m[3][1], - mx.m[3][0] * my.m[0][2] + mx.m[3][1] * my.m[1][2] + mx.m[3][2] * my.m[2][2] + mx.m[3][3] * my.m[3][2], - mx.m[3][0] * my.m[0][3] + mx.m[3][1] * my.m[1][3] + mx.m[3][2] * my.m[2][3] + mx.m[3][3] * my.m[3][3]) + return new_affinematrix3d( + mx.m[0][0] * my.m[0][0] + mx.m[0][1] * my.m[1][0] + mx.m[0][2] * my.m[2][0] + mx.m[0][3] * my.m[3][0], + mx.m[0][0] * my.m[0][1] + mx.m[0][1] * my.m[1][1] + mx.m[0][2] * my.m[2][1] + mx.m[0][3] * my.m[3][1], + mx.m[0][0] * my.m[0][2] + mx.m[0][1] * my.m[1][2] + mx.m[0][2] * my.m[2][2] + mx.m[0][3] * my.m[3][2], + mx.m[0][0] * my.m[0][3] + mx.m[0][1] * my.m[1][3] + mx.m[0][2] * my.m[2][3] + mx.m[0][3] * my.m[3][3], + mx.m[1][0] * my.m[0][0] + mx.m[1][1] * my.m[1][0] + mx.m[1][2] * my.m[2][0] + mx.m[1][3] * my.m[3][0], + mx.m[1][0] * my.m[0][1] + mx.m[1][1] * my.m[1][1] + mx.m[1][2] * my.m[2][1] + mx.m[1][3] * my.m[3][1], + mx.m[1][0] * my.m[0][2] + mx.m[1][1] * my.m[1][2] + mx.m[1][2] * my.m[2][2] + mx.m[1][3] * my.m[3][2], + mx.m[1][0] * my.m[0][3] + mx.m[1][1] * my.m[1][3] + mx.m[1][2] * my.m[2][3] + mx.m[1][3] * my.m[3][3], + mx.m[2][0] * my.m[0][0] + mx.m[2][1] * my.m[1][0] + mx.m[2][2] * my.m[2][0] + mx.m[2][3] * my.m[3][0], + mx.m[2][0] * my.m[0][1] + mx.m[2][1] * my.m[1][1] + mx.m[2][2] * my.m[2][1] + mx.m[2][3] * my.m[3][1], + mx.m[2][0] * my.m[0][2] + mx.m[2][1] * my.m[1][2] + mx.m[2][2] * my.m[2][2] + mx.m[2][3] * my.m[3][2], + mx.m[2][0] * my.m[0][3] + mx.m[2][1] * my.m[1][3] + mx.m[2][2] * my.m[2][3] + mx.m[2][3] * my.m[3][3], + mx.m[3][0] * my.m[0][0] + mx.m[3][1] * my.m[1][0] + mx.m[3][2] * my.m[2][0] + mx.m[3][3] * my.m[3][0], + mx.m[3][0] * my.m[0][1] + mx.m[3][1] * my.m[1][1] + mx.m[3][2] * my.m[2][1] + mx.m[3][3] * my.m[3][1], + mx.m[3][0] * my.m[0][2] + mx.m[3][1] * my.m[1][2] + mx.m[3][2] * my.m[2][2] + mx.m[3][3] * my.m[3][2], + mx.m[3][0] * my.m[0][3] + mx.m[3][1] * my.m[1][3] + mx.m[3][2] * my.m[2][3] + mx.m[3][3] * my.m[3][3] + ) return NotImplemented @@ -231,38 +232,42 @@ cdef class AffineMatrix3D(_Mat4): t[16] = self.m[1][1] * self.m[3][3] - self.m[1][3] * self.m[3][1] t[17] = self.m[1][2] * self.m[3][3] - self.m[1][3] * self.m[3][2] - return new_affinematrix3d((self.m[2][2] * t[16] - self.m[2][1] * t[17] - self.m[2][3] * t[15]) * idet, - (self.m[2][1] * t[11] - self.m[2][2] * t[10] + self.m[2][3] * t[ 9]) * idet, - (self.m[3][1] * t[ 5] - self.m[3][2] * t[ 4] + self.m[3][3] * t[ 3]) * idet, - -t[21] * idet, - (self.m[2][0] * t[17] - self.m[2][2] * t[14] + self.m[2][3] * t[13]) * idet, - (self.m[2][2] * t[ 8] - self.m[2][0] * t[11] - self.m[2][3] * t[ 7]) * idet, - (self.m[3][2] * t[ 2] - self.m[3][0] * t[ 5] - self.m[3][3] * t[ 1]) * idet, - t[20] * idet, - (self.m[2][1] * t[14] - self.m[2][0] * t[16] - self.m[2][3] * t[12]) * idet, - (self.m[2][0] * t[10] - self.m[2][1] * t[ 8] + self.m[2][3] * t[ 6]) * idet, - (self.m[3][0] * t[ 4] - self.m[3][1] * t[ 2] + self.m[3][3] * t[ 0]) * idet, - -t[19] * idet, - (self.m[2][0] * t[15] - self.m[2][1] * t[13] + self.m[2][2] * t[12]) * idet, - (self.m[2][1] * t[ 7] - self.m[2][0] * t[ 9] - self.m[2][2] * t[ 6]) * idet, - (self.m[3][1] * t[ 1] - self.m[3][0] * t[ 3] - self.m[3][2] * t[ 0]) * idet, - t[18] * idet) + return new_affinematrix3d( + (self.m[2][2] * t[16] - self.m[2][1] * t[17] - self.m[2][3] * t[15]) * idet, + (self.m[2][1] * t[11] - self.m[2][2] * t[10] + self.m[2][3] * t[ 9]) * idet, + (self.m[3][1] * t[ 5] - self.m[3][2] * t[ 4] + self.m[3][3] * t[ 3]) * idet, + -t[21] * idet, + (self.m[2][0] * t[17] - self.m[2][2] * t[14] + self.m[2][3] * t[13]) * idet, + (self.m[2][2] * t[ 8] - self.m[2][0] * t[11] - self.m[2][3] * t[ 7]) * idet, + (self.m[3][2] * t[ 2] - self.m[3][0] * t[ 5] - self.m[3][3] * t[ 1]) * idet, + t[20] * idet, + (self.m[2][1] * t[14] - self.m[2][0] * t[16] - self.m[2][3] * t[12]) * idet, + (self.m[2][0] * t[10] - self.m[2][1] * t[ 8] + self.m[2][3] * t[ 6]) * idet, + (self.m[3][0] * t[ 4] - self.m[3][1] * t[ 2] + self.m[3][3] * t[ 0]) * idet, + -t[19] * idet, + (self.m[2][0] * t[15] - self.m[2][1] * t[13] + self.m[2][2] * t[12]) * idet, + (self.m[2][1] * t[ 7] - self.m[2][0] * t[ 9] - self.m[2][2] * t[ 6]) * idet, + (self.m[3][1] * t[ 1] - self.m[3][0] * t[ 3] - self.m[3][2] * t[ 0]) * idet, + t[18] * idet + ) cdef AffineMatrix3D mul(self, AffineMatrix3D m): - return new_affinematrix3d(self.m[0][0] * m.m[0][0] + self.m[0][1] * m.m[1][0] + self.m[0][2] * m.m[2][0] + self.m[0][3] * m.m[3][0], - self.m[0][0] * m.m[0][1] + self.m[0][1] * m.m[1][1] + self.m[0][2] * m.m[2][1] + self.m[0][3] * m.m[3][1], - self.m[0][0] * m.m[0][2] + self.m[0][1] * m.m[1][2] + self.m[0][2] * m.m[2][2] + self.m[0][3] * m.m[3][2], - self.m[0][0] * m.m[0][3] + self.m[0][1] * m.m[1][3] + self.m[0][2] * m.m[2][3] + self.m[0][3] * m.m[3][3], - self.m[1][0] * m.m[0][0] + self.m[1][1] * m.m[1][0] + self.m[1][2] * m.m[2][0] + self.m[1][3] * m.m[3][0], - self.m[1][0] * m.m[0][1] + self.m[1][1] * m.m[1][1] + self.m[1][2] * m.m[2][1] + self.m[1][3] * m.m[3][1], - self.m[1][0] * m.m[0][2] + self.m[1][1] * m.m[1][2] + self.m[1][2] * m.m[2][2] + self.m[1][3] * m.m[3][2], - self.m[1][0] * m.m[0][3] + self.m[1][1] * m.m[1][3] + self.m[1][2] * m.m[2][3] + self.m[1][3] * m.m[3][3], - self.m[2][0] * m.m[0][0] + self.m[2][1] * m.m[1][0] + self.m[2][2] * m.m[2][0] + self.m[2][3] * m.m[3][0], - self.m[2][0] * m.m[0][1] + self.m[2][1] * m.m[1][1] + self.m[2][2] * m.m[2][1] + self.m[2][3] * m.m[3][1], - self.m[2][0] * m.m[0][2] + self.m[2][1] * m.m[1][2] + self.m[2][2] * m.m[2][2] + self.m[2][3] * m.m[3][2], - self.m[2][0] * m.m[0][3] + self.m[2][1] * m.m[1][3] + self.m[2][2] * m.m[2][3] + self.m[2][3] * m.m[3][3], - self.m[3][0] * m.m[0][0] + self.m[3][1] * m.m[1][0] + self.m[3][2] * m.m[2][0] + self.m[3][3] * m.m[3][0], - self.m[3][0] * m.m[0][1] + self.m[3][1] * m.m[1][1] + self.m[3][2] * m.m[2][1] + self.m[3][3] * m.m[3][1], - self.m[3][0] * m.m[0][2] + self.m[3][1] * m.m[1][2] + self.m[3][2] * m.m[2][2] + self.m[3][3] * m.m[3][2], - self.m[3][0] * m.m[0][3] + self.m[3][1] * m.m[1][3] + self.m[3][2] * m.m[2][3] + self.m[3][3] * m.m[3][3]) + return new_affinematrix3d( + self.m[0][0] * m.m[0][0] + self.m[0][1] * m.m[1][0] + self.m[0][2] * m.m[2][0] + self.m[0][3] * m.m[3][0], + self.m[0][0] * m.m[0][1] + self.m[0][1] * m.m[1][1] + self.m[0][2] * m.m[2][1] + self.m[0][3] * m.m[3][1], + self.m[0][0] * m.m[0][2] + self.m[0][1] * m.m[1][2] + self.m[0][2] * m.m[2][2] + self.m[0][3] * m.m[3][2], + self.m[0][0] * m.m[0][3] + self.m[0][1] * m.m[1][3] + self.m[0][2] * m.m[2][3] + self.m[0][3] * m.m[3][3], + self.m[1][0] * m.m[0][0] + self.m[1][1] * m.m[1][0] + self.m[1][2] * m.m[2][0] + self.m[1][3] * m.m[3][0], + self.m[1][0] * m.m[0][1] + self.m[1][1] * m.m[1][1] + self.m[1][2] * m.m[2][1] + self.m[1][3] * m.m[3][1], + self.m[1][0] * m.m[0][2] + self.m[1][1] * m.m[1][2] + self.m[1][2] * m.m[2][2] + self.m[1][3] * m.m[3][2], + self.m[1][0] * m.m[0][3] + self.m[1][1] * m.m[1][3] + self.m[1][2] * m.m[2][3] + self.m[1][3] * m.m[3][3], + self.m[2][0] * m.m[0][0] + self.m[2][1] * m.m[1][0] + self.m[2][2] * m.m[2][0] + self.m[2][3] * m.m[3][0], + self.m[2][0] * m.m[0][1] + self.m[2][1] * m.m[1][1] + self.m[2][2] * m.m[2][1] + self.m[2][3] * m.m[3][1], + self.m[2][0] * m.m[0][2] + self.m[2][1] * m.m[1][2] + self.m[2][2] * m.m[2][2] + self.m[2][3] * m.m[3][2], + self.m[2][0] * m.m[0][3] + self.m[2][1] * m.m[1][3] + self.m[2][2] * m.m[2][3] + self.m[2][3] * m.m[3][3], + self.m[3][0] * m.m[0][0] + self.m[3][1] * m.m[1][0] + self.m[3][2] * m.m[2][0] + self.m[3][3] * m.m[3][0], + self.m[3][0] * m.m[0][1] + self.m[3][1] * m.m[1][1] + self.m[3][2] * m.m[2][1] + self.m[3][3] * m.m[3][1], + self.m[3][0] * m.m[0][2] + self.m[3][1] * m.m[1][2] + self.m[3][2] * m.m[2][2] + self.m[3][3] * m.m[3][2], + self.m[3][0] * m.m[0][3] + self.m[3][1] * m.m[1][3] + self.m[3][2] * m.m[2][3] + self.m[3][3] * m.m[3][3] + ) diff --git a/raysect/core/math/cython/triangle.pxd b/raysect/core/math/cython/triangle.pxd index 61597ea7..1530798c 100644 --- a/raysect/core/math/cython/triangle.pxd +++ b/raysect/core/math/cython/triangle.pxd @@ -31,16 +31,16 @@ cdef bint inside_triangle(double v1x, double v1y, double v2x, double v2y, - double v3x, double v3y, double px, double py) nogil + double v3x, double v3y, double px, double py) nogil cdef void barycentric_coords(double v1x, double v1y, double v2x, double v2y, - double v3x, double v3y, double px, double py, - double *alpha, double *beta, double *gamma) nogil + double v3x, double v3y, double px, double py, + double *alpha, double *beta, double *gamma) nogil cdef bint barycentric_inside_triangle(double alpha, double beta, double gamma) nogil cdef double barycentric_interpolation(double alpha, double beta, double gamma, - double va, double vb, double vc) nogil + double va, double vb, double vc) nogil diff --git a/raysect/core/math/cython/triangle.pyx b/raysect/core/math/cython/triangle.pyx index a3e20a6f..8e3c5ecb 100644 --- a/raysect/core/math/cython/triangle.pyx +++ b/raysect/core/math/cython/triangle.pyx @@ -33,7 +33,7 @@ cimport cython cdef bint inside_triangle(double v1x, double v1y, double v2x, double v2y, - double v3x, double v3y, double px, double py) nogil: + double v3x, double v3y, double px, double py) nogil: """ Cython utility for testing if point is inside a triangle. @@ -52,8 +52,7 @@ cdef bint inside_triangle(double v1x, double v1y, double v2x, double v2y, :rtype: bool """ - cdef: - double ux, uy, vx, vy + cdef double ux, uy, vx, vy # calculate vectors ux = v2x - v1x @@ -158,7 +157,7 @@ cdef bint barycentric_inside_triangle(double alpha, double beta, double gamma) n cdef double barycentric_interpolation(double alpha, double beta, double gamma, - double va, double vb, double vc) nogil: + double va, double vb, double vc) nogil: """ Cython utility for interpolation of data at triangle vertices. diff --git a/raysect/core/math/cython/utility.pyx b/raysect/core/math/cython/utility.pyx index 9cc8147e..8f8bdcb9 100644 --- a/raysect/core/math/cython/utility.pyx +++ b/raysect/core/math/cython/utility.pyx @@ -114,8 +114,7 @@ cdef double interpolate(double[::1] x, double[::1] y, double p) nogil: :rtype: double """ - cdef: - int index, top_index + cdef int index, top_index index = find_index(x, p) @@ -202,7 +201,6 @@ cdef double integrate(double[::1] x, double[::1] y, double x0, double x1) nogil: else: integral_sum = 0.0 - if lower_index == 0: # add contribution from point below array diff --git a/raysect/core/math/function/float/function1d/base.pyx b/raysect/core/math/function/float/function1d/base.pyx index c078e6ee..7695fc99 100644 --- a/raysect/core/math/function/float/function1d/base.pyx +++ b/raysect/core/math/function/float/function1d/base.pyx @@ -65,110 +65,161 @@ cdef class Function1D(FloatFunction): def __repr__(self): return 'Function1D(x)' - def __add__(object a, object b): - if is_callable(a): - if is_callable(b): - # a() + b() - return AddFunction1D(a, b) - elif isinstance(b, numbers.Real): - # a() + B -> B + a() - return AddScalar1D( b, a) - elif isinstance(a, numbers.Real): - if is_callable(b): - # A + b() - return AddScalar1D( a, b) + def __add__(self, object b): + + if is_callable(b): + # a() + b() + return AddFunction1D(self, b) + + elif isinstance(b, numbers.Real): + # a() + B -> B + a() + return AddScalar1D( b, self) + return NotImplemented + + def __radd__(self, object a): + return self.__add__(a) + + def __sub__(self, object b): + + if is_callable(b): + # a() - b() + return SubtractFunction1D(self, b) + + elif isinstance(b, numbers.Real): + # a() - B -> -B + a() + return AddScalar1D(-( b), self) - def __sub__(object a, object b): - if is_callable(a): - if is_callable(b): - # a() - b() - return SubtractFunction1D(a, b) - elif isinstance(b, numbers.Real): - # a() - B -> -B + a() - return AddScalar1D(-( b), a) - elif isinstance(a, numbers.Real): - if is_callable(b): - # A - b() - return SubtractScalar1D( a, b) return NotImplemented + + def __rsub__(self, object a): - def __mul__(object a, object b): if is_callable(a): - if is_callable(b): - # a() * b() - return MultiplyFunction1D(a, b) - elif isinstance(b, numbers.Real): - # a() * B -> B * a() - return MultiplyScalar1D( b, a) + # a() - b() + return SubtractFunction1D(a, self) + elif isinstance(a, numbers.Real): - if is_callable(b): - # A * b() - return MultiplyScalar1D( a, b) + # A - b() + return SubtractScalar1D( a, self) + return NotImplemented + def __mul__(self, object b): + + if is_callable(b): + # a() * b() + return MultiplyFunction1D(self, b) + + elif isinstance(b, numbers.Real): + # a() * B -> B * a() + return MultiplyScalar1D( b, self) + + return NotImplemented + + def __rmul__(self, object a): + return self.__mul__(a) + @cython.cdivision(True) - def __truediv__(object a, object b): + def __truediv__(self, object b): + cdef double v + + if is_callable(b): + # a() / b() + return DivideFunction1D(self, b) + + elif isinstance(b, numbers.Real): + # a() / B -> 1/B * a() + v = b + if v == 0.0: + raise ZeroDivisionError("Scalar used as the denominator of the division is zero valued.") + return MultiplyScalar1D(1 / v, self) + + return NotImplemented + + def __rtruediv__(self, object a): + if is_callable(a): - if is_callable(b): - # a() / b() - return DivideFunction1D(a, b) - elif isinstance(b, numbers.Real): - # a() / B -> 1/B * a() - v = b - if v == 0.0: - raise ZeroDivisionError("Scalar used as the denominator of the division is zero valued.") - return MultiplyScalar1D(1/v, a) + # a() / b() + return DivideFunction1D(a, self) + elif isinstance(a, numbers.Real): - if is_callable(b): - # A * b() - return DivideScalar1D( a, b) + # A / b() + return DivideScalar1D( a, self) + return NotImplemented - def __mod__(object a, object b): + def __mod__(self, object b): + cdef double v + + if is_callable(b): + # a() % b() + return ModuloFunction1D(self, b) + + elif isinstance(b, numbers.Real): + # a() % B + v = b + if v == 0.0: + raise ZeroDivisionError("Scalar used as the divisor of the division is zero valued.") + return ModuloFunctionScalar1D(self, v) + + return NotImplemented + + def __rmod__(self, object a): + if is_callable(a): - if is_callable(b): - # a() % b() - return ModuloFunction1D(a, b) - elif isinstance(b, numbers.Real): - # a() % B - v = b - if v == 0.0: - raise ZeroDivisionError("Scalar used as the divisor of the division is zero valued.") - return ModuloFunctionScalar1D(a, v) + # a() % b() + return ModuloFunction1D(a, self) + elif isinstance(a, numbers.Real): - if is_callable(b): - # A % b() - return ModuloScalarFunction1D( a, b) + # A % b() + return ModuloScalarFunction1D( a, self) + return NotImplemented def __neg__(self): return MultiplyScalar1D(-1, self) - def __pow__(object a, object b, object c): + def __pow__(self, object b, object c): + if c is not None: # Optimised implementation of pow(a, b, c) not available: fall back # to general implementation - return (a ** b) % c + return PowFunction1D(self, b) % c + + if is_callable(b): + # a() ** b() + return PowFunction1D(self, b) + + elif isinstance(b, numbers.Real): + # a() ** b + return PowFunctionScalar1D(self, b) + + return NotImplemented + + def __rpow__(self, object a, object c): + + if c is not None: + # Optimised implementation of pow(a, b, c) not available: fall back + # to general implementation + return PowFunction1D(a, self) % c + if is_callable(a): - if is_callable(b): - # a() ** b() - return PowFunction1D(a, b) - elif isinstance(b, numbers.Real): - # a() ** b - return PowFunctionScalar1D(a, b) + # a() ** b() + return PowFunction1D(a, self) + elif isinstance(a, numbers.Real): - if is_callable(b): - # a ** b() - return PowScalarFunction1D( a, b) + # A ** b() + return PowScalarFunction1D( a, self) + return NotImplemented def __abs__(self): return AbsFunction1D(self) def __richcmp__(self, object other, int op): + if is_callable(other): if op == Py_EQ: return EqualsFunction1D(self, other) @@ -182,6 +233,7 @@ cdef class Function1D(FloatFunction): return LessEqualsFunction1D(self, other) if op == Py_GE: return GreaterEqualsFunction1D(self, other) + if isinstance(other, numbers.Real): if op == Py_EQ: return EqualsScalar1D( other, self) @@ -199,6 +251,7 @@ cdef class Function1D(FloatFunction): if op == Py_GE: # f() >= K -> K <= f return LessEqualsScalar1D( other, self) + return NotImplemented diff --git a/raysect/core/math/function/float/function1d/samplers.pyx b/raysect/core/math/function/float/function1d/samplers.pyx index 16c6649b..e17f5439 100644 --- a/raysect/core/math/function/float/function1d/samplers.pyx +++ b/raysect/core/math/function/float/function1d/samplers.pyx @@ -35,6 +35,7 @@ from .autowrap cimport autowrap_function1d cimport cython cimport numpy as np + @cython.boundscheck(False) @cython.wraparound(False) cpdef tuple sample1d(object function, double x_min, double x_max, int x_samples): @@ -58,7 +59,6 @@ cpdef tuple sample1d(object function, double x_min, double x_max, int x_samples) if x_samples < 1: raise ValueError("The argument x_samples must be >= 1") - # ensures that func is of type Function1D. I.e. if 'function' argument was Python function, it'll get autowrapped # into Function1D object func = autowrap_function1d(function) diff --git a/raysect/core/math/function/float/function2d/base.pyx b/raysect/core/math/function/float/function2d/base.pyx index 933fc8b4..a689415b 100644 --- a/raysect/core/math/function/float/function2d/base.pyx +++ b/raysect/core/math/function/float/function2d/base.pyx @@ -63,110 +63,162 @@ cdef class Function2D(FloatFunction): """ return self.evaluate(x, y) - def __add__(object a, object b): - if is_callable(a): - if is_callable(b): - # a() + b() - return AddFunction2D(a, b) - elif isinstance(b, numbers.Real): - # a() + B -> B + a() - return AddScalar2D( b, a) - elif isinstance(a, numbers.Real): - if is_callable(b): - # A + b() - return AddScalar2D( a, b) + def __add__(self, object b): + + if is_callable(b): + # a() + b() + return AddFunction2D(self, b) + + elif isinstance(b, numbers.Real): + # a() + B -> B + a() + return AddScalar2D( b, self) + return NotImplemented - def __sub__(object a, object b): - if is_callable(a): - if is_callable(b): - # a() - b() - return SubtractFunction2D(a, b) - elif isinstance(b, numbers.Real): - # a() - B -> -B + a() - return AddScalar2D(-( b), a) - elif isinstance(a, numbers.Real): - if is_callable(b): - # A - b() - return SubtractScalar2D( a, b) + def __radd__(self, object a): + return self.__add__(a) + + def __sub__(self, object b): + + if is_callable(b): + # a() - b() + return SubtractFunction2D(self, b) + + elif isinstance(b, numbers.Real): + # a() - B -> -B + a() + return AddScalar2D(-( b), self) + return NotImplemented - def __mul__(object a, object b): + def __rsub__(self, object a): + if is_callable(a): - if is_callable(b): - # a() * b() - return MultiplyFunction2D(a, b) - elif isinstance(b, numbers.Real): - # a() * B -> B * a() - return MultiplyScalar2D( b, a) + # a() - b() + return SubtractFunction2D(a, self) + elif isinstance(a, numbers.Real): - if is_callable(b): - # A * b() - return MultiplyScalar2D( a, b) + # A - b() + return SubtractScalar2D( a, self) + return NotImplemented + def __mul__(self, object b): + + if is_callable(b): + # a() * b() + return MultiplyFunction2D(self, b) + + elif isinstance(b, numbers.Real): + # a() * B -> B * a() + return MultiplyScalar2D( b, self) + + return NotImplemented + + def __rmul__(self, object a): + return self.__mul__(a) + @cython.cdivision(True) - def __truediv__(object a, object b): + def __truediv__(self, object b): + cdef double v + + if is_callable(b): + # a() / b() + return DivideFunction2D(self, b) + + elif isinstance(b, numbers.Real): + # a() / B -> 1/B * a() + v = b + if v == 0.0: + raise ZeroDivisionError("Scalar used as the denominator of the division is zero valued.") + return MultiplyScalar2D(1/v, self) + + return NotImplemented + + @cython.cdivision(True) + def __rtruediv__(self, object a): + if is_callable(a): - if is_callable(b): - # a() / b() - return DivideFunction2D(a, b) - elif isinstance(b, numbers.Real): - # a() / B -> 1/B * a() - v = b - if v == 0.0: - raise ZeroDivisionError("Scalar used as the denominator of the division is zero valued.") - return MultiplyScalar2D(1/v, a) + # a() / b() + return DivideFunction2D(a, self) + elif isinstance(a, numbers.Real): - if is_callable(b): - # A * b() - return DivideScalar2D( a, b) + # A / b() + return DivideScalar2D( a, self) + return NotImplemented - def __mod__(object a, object b): + def __mod__(self, object b): + cdef double v + + if is_callable(b): + # a() % b() + return ModuloFunction2D(self, b) + + elif isinstance(b, numbers.Real): + # a() % B + v = b + if v == 0.0: + raise ZeroDivisionError("Scalar used as the divisor of the division is zero valued.") + return ModuloFunctionScalar2D(self, v) + + return NotImplemented + + def __rmod__(self, object a): + if is_callable(a): - if is_callable(b): - # a() % b() - return ModuloFunction2D(a, b) - elif isinstance(b, numbers.Real): - # a() % B - v = b - if v == 0.0: - raise ZeroDivisionError("Scalar used as the divisor of the division is zero valued.") - return ModuloFunctionScalar2D(a, v) + # a() % b() + return ModuloFunction2D(a, self) + elif isinstance(a, numbers.Real): - if is_callable(b): - # A % b() - return ModuloScalarFunction2D( a, b) + # A % b() + return ModuloScalarFunction2D( a, self) + return NotImplemented def __neg__(self): return MultiplyScalar2D(-1, self) - def __pow__(object a, object b, object c): + def __pow__(self, object b, object c): + if c is not None: # Optimised implementation of pow(a, b, c) not available: fall back # to general implementation - return (a ** b) % c + return PowFunction2D(self, b) % c + + if is_callable(b): + # a() ** b() + return PowFunction2D(self, b) + + elif isinstance(b, numbers.Real): + # a() ** b + return PowFunctionScalar2D(self, b) + + return NotImplemented + + def __rpow__(self, object a, object c): + + if c is not None: + # Optimised implementation of pow(a, b, c) not available: fall back + # to general implementation + return PowFunction2D(a, self) % c + if is_callable(a): - if is_callable(b): - # a() ** b() - return PowFunction2D(a, b) - elif isinstance(b, numbers.Real): - # a() ** b - return PowFunctionScalar2D(a, b) + # a() ** b() + return PowFunction2D(a, self) + elif isinstance(a, numbers.Real): - if is_callable(b): - # a ** b() - return PowScalarFunction2D( a, b) + # A ** b() + return PowScalarFunction2D( a, self) + return NotImplemented def __abs__(self): return AbsFunction2D(self) def __richcmp__(self, object other, int op): + if is_callable(other): if op == Py_EQ: return EqualsFunction2D(self, other) @@ -180,6 +232,7 @@ cdef class Function2D(FloatFunction): return LessEqualsFunction2D(self, other) if op == Py_GE: return GreaterEqualsFunction2D(self, other) + if isinstance(other, numbers.Real): if op == Py_EQ: return EqualsScalar2D( other, self) @@ -197,6 +250,7 @@ cdef class Function2D(FloatFunction): if op == Py_GE: # f() >= K -> K <= f return LessEqualsScalar2D( other, self) + return NotImplemented diff --git a/raysect/core/math/function/float/function2d/interpolate/common.pyx b/raysect/core/math/function/float/function2d/interpolate/common.pyx index 4a4ae152..58fb27e0 100644 --- a/raysect/core/math/function/float/function2d/interpolate/common.pyx +++ b/raysect/core/math/function/float/function2d/interpolate/common.pyx @@ -38,15 +38,16 @@ from raysect.core.math.cython cimport barycentric_inside_triangle, barycentric_c cimport cython # bounding box is padded by a small amount to avoid numerical accuracy issues -DEF BOX_PADDING = 1e-6 +cdef const double BOX_PADDING = 1e-6 # convenience defines -DEF V1 = 0 -DEF V2 = 1 -DEF V3 = 2 +cdef enum: + V1 = 0 + V2 = 1 + V3 = 2 -DEF X = 0 -DEF Y = 1 + X = 0 + Y = 1 cdef class MeshKDTree2D(KDTree2DCore): diff --git a/raysect/core/math/function/float/function3d/base.pyx b/raysect/core/math/function/float/function3d/base.pyx index 80f80b95..f5cdf065 100644 --- a/raysect/core/math/function/float/function3d/base.pyx +++ b/raysect/core/math/function/float/function3d/base.pyx @@ -64,110 +64,162 @@ cdef class Function3D(FloatFunction): """ return self.evaluate(x, y, z) - def __add__(object a, object b): - if is_callable(a): - if is_callable(b): - # a() + b() - return AddFunction3D(a, b) - elif isinstance(b, numbers.Real): - # a() + B -> B + a() - return AddScalar3D( b, a) - elif isinstance(a, numbers.Real): - if is_callable(b): - # A + b() - return AddScalar3D( a, b) + def __add__(self, object b): + + if is_callable(b): + # a() + b() + return AddFunction3D(self, b) + + elif isinstance(b, numbers.Real): + # a() + B -> B + a() + return AddScalar3D( b, self) + return NotImplemented - def __sub__(object a, object b): - if is_callable(a): - if is_callable(b): - # a() - b() - return SubtractFunction3D(a, b) - elif isinstance(b, numbers.Real): - # a() - B -> -B + a() - return AddScalar3D(-( b), a) - elif isinstance(a, numbers.Real): - if is_callable(b): - # A - b() - return SubtractScalar3D( a, b) + def __radd__(self, object a): + return self.__add__(a) + + def __sub__(self, object b): + + if is_callable(b): + # a() - b() + return SubtractFunction3D(self, b) + + elif isinstance(b, numbers.Real): + # a() - B -> -B + a() + return AddScalar3D(-( b), self) + return NotImplemented + + def __rsub__(self, object a): - def __mul__(object a, object b): if is_callable(a): - if is_callable(b): - # a() * b() - return MultiplyFunction3D(a, b) - elif isinstance(b, numbers.Real): - # a() * B -> B * a() - return MultiplyScalar3D( b, a) + # a() - b() + return SubtractFunction3D(a, self) + elif isinstance(a, numbers.Real): - if is_callable(b): - # A * b() - return MultiplyScalar3D( a, b) + # A - b() + return SubtractScalar3D( a, self) + return NotImplemented + def __mul__(self, object b): + + if is_callable(b): + # a() * b() + return MultiplyFunction3D(self, b) + + elif isinstance(b, numbers.Real): + # a() * B -> B * a() + return MultiplyScalar3D( b, self) + + return NotImplemented + + def __rmul__(self, object a): + return self.__mul__(a) + @cython.cdivision(True) - def __truediv__(object a, object b): + def __truediv__(self, object b): + cdef double v + + if is_callable(b): + # a() / b() + return DivideFunction3D(self, b) + + elif isinstance(b, numbers.Real): + # a() / B -> 1/B * a() + v = b + if v == 0.0: + raise ZeroDivisionError("Scalar used as the denominator of the division is zero valued.") + return MultiplyScalar3D(1/v, self) + + return NotImplemented + + @cython.cdivision(True) + def __rtruediv__(self, object a): + if is_callable(a): - if is_callable(b): - # a() / b() - return DivideFunction3D(a, b) - elif isinstance(b, numbers.Real): - # a() / B -> 1/B * a() - v = b - if v == 0.0: - raise ZeroDivisionError("Scalar used as the denominator of the division is zero valued.") - return MultiplyScalar3D(1/v, a) + # a() / b() + return DivideFunction3D(a, self) + elif isinstance(a, numbers.Real): - if is_callable(b): - # A * b() - return DivideScalar3D( a, b) + # A / b() + return DivideScalar3D( a, self) + return NotImplemented - def __mod__(object a, object b): + def __mod__(self, object b): + cdef double v + + if is_callable(b): + # a() % b() + return ModuloFunction3D(self, b) + + elif isinstance(b, numbers.Real): + # a() % B + v = b + if v == 0.0: + raise ZeroDivisionError("Scalar used as the divisor of the division is zero valued.") + return ModuloFunctionScalar3D(self, v) + + return NotImplemented + + def __rmod__(self, object a): + if is_callable(a): - if is_callable(b): - # a() % b() - return ModuloFunction3D(a, b) - elif isinstance(b, numbers.Real): - # a() % B - v = b - if v == 0.0: - raise ZeroDivisionError("Scalar used as the divisor of the division is zero valued.") - return ModuloFunctionScalar3D(a, v) + # a() % b() + return ModuloFunction3D(a, self) + elif isinstance(a, numbers.Real): - if is_callable(b): - # A % b() - return ModuloScalarFunction3D( a, b) + # A % b() + return ModuloScalarFunction3D( a, self) + return NotImplemented def __neg__(self): return MultiplyScalar3D(-1, self) - def __pow__(object a, object b, object c): + def __pow__(self, object b, object c): + if c is not None: # Optimised implementation of pow(a, b, c) not available: fall back # to general implementation - return (a ** b) % c + return PowFunction3D(self, b) % c + + if is_callable(b): + # a() ** b() + return PowFunction3D(self, b) + + elif isinstance(b, numbers.Real): + # a() ** b + return PowFunctionScalar3D(self, b) + + return NotImplemented + + def __rpow__(self, object a, object c): + + if c is not None: + # Optimised implementation of pow(a, b, c) not available: fall back + # to general implementation + return PowFunction3D(a, self) % c + if is_callable(a): - if is_callable(b): - # a() ** b() - return PowFunction3D(a, b) - elif isinstance(b, numbers.Real): - # a() ** b - return PowFunctionScalar3D(a, b) + # a() ** b() + return PowFunction3D(a, self) + elif isinstance(a, numbers.Real): - if is_callable(b): - # a ** b() - return PowScalarFunction3D( a, b) + # A ** b() + return PowScalarFunction3D( a, self) + return NotImplemented def __abs__(self): return AbsFunction3D(self) def __richcmp__(self, object other, int op): + if is_callable(other): if op == Py_EQ: return EqualsFunction3D(self, other) @@ -181,6 +233,7 @@ cdef class Function3D(FloatFunction): return LessEqualsFunction3D(self, other) if op == Py_GE: return GreaterEqualsFunction3D(self, other) + if isinstance(other, numbers.Real): if op == Py_EQ: return EqualsScalar3D( other, self) @@ -198,6 +251,7 @@ cdef class Function3D(FloatFunction): if op == Py_GE: # f() >= K -> K <= f return LessEqualsScalar3D( other, self) + return NotImplemented diff --git a/raysect/core/math/function/float/function3d/interpolate/common.pyx b/raysect/core/math/function/float/function3d/interpolate/common.pyx index 9b65345a..4fa0cdf0 100644 --- a/raysect/core/math/function/float/function3d/interpolate/common.pyx +++ b/raysect/core/math/function/float/function3d/interpolate/common.pyx @@ -38,17 +38,18 @@ from raysect.core.math.cython cimport barycentric_inside_tetrahedra, barycentric cimport cython # bounding box is padded by a small amount to avoid numerical accuracy issues -DEF BOX_PADDING = 1e-6 +cdef const double BOX_PADDING = 1e-6 # convenience defines -DEF V1 = 0 -DEF V2 = 1 -DEF V3 = 2 -DEF V4 = 3 - -DEF X = 0 -DEF Y = 1 -DEF Z = 2 +cdef enum: + V1 = 0 + V2 = 1 + V3 = 2 + V4 = 3 + + X = 0 + Y = 1 + Z = 2 cdef class MeshKDTree3D(KDTree3DCore): diff --git a/raysect/core/math/function/vector3d/function1d/base.pyx b/raysect/core/math/function/vector3d/function1d/base.pyx index 0c1e371c..5b920431 100644 --- a/raysect/core/math/function/vector3d/function1d/base.pyx +++ b/raysect/core/math/function/vector3d/function1d/base.pyx @@ -65,31 +65,34 @@ cdef class Function1D(Vector3DFunction): """ return self.evaluate(x) - def __add__(object a, object b): - if is_callable(a) or isinstance(a, Vector3D): - if is_callable(b) or isinstance(b, Vector3D): - return AddFunction1D(a, b) + def __add__(self, object b): + if is_callable(b) or isinstance(b, Vector3D): + return AddFunction1D(self, b) return NotImplemented + + def __radd__(self, object a): + return self.__add__(a) - def __sub__(object a, object b): - if is_callable(a) or isinstance(a, Vector3D): - if is_callable(b) or isinstance(b, Vector3D): - return SubtractFunction1D(a, b) + def __sub__(self, object b): + if is_callable(b) or isinstance(b, Vector3D): + return SubtractFunction1D(self, b) return NotImplemented - - def __mul__(object a, object b): + + def __rsub__(self, object a): if is_callable(a) or isinstance(a, Vector3D): - if float_is_callable(b) or isinstance(b, numbers.Real): - return MultiplyFunction1D(a, b) - if is_callable(b) or isinstance(b, Vector3D): - if float_is_callable(a) or isinstance(a, numbers.Real): - return MultiplyFunction1D(b, a) + return SubtractFunction1D(a, self) + + def __mul__(self, object b): + if float_is_callable(b) or isinstance(b, numbers.Real): + return MultiplyFunction1D(self, b) return NotImplemented + + def __rmul__(self, object a): + return self.__mul__(a) - def __truediv__(object a, object b): - if is_callable(a) or isinstance(a, Vector3D): - if float_is_callable(b) or isinstance(b, numbers.Real): - return DivideFunction1D(a, b) + def __truediv__(self, object b): + if float_is_callable(b) or isinstance(b, numbers.Real): + return DivideFunction1D(self, b) return NotImplemented def __neg__(self): diff --git a/raysect/core/math/function/vector3d/function2d/base.pyx b/raysect/core/math/function/vector3d/function2d/base.pyx index 021dbf71..d7856199 100644 --- a/raysect/core/math/function/vector3d/function2d/base.pyx +++ b/raysect/core/math/function/vector3d/function2d/base.pyx @@ -66,31 +66,34 @@ cdef class Function2D(Vector3DFunction): """ return self.evaluate(x, y) - def __add__(object a, object b): - if is_callable(a) or isinstance(a, Vector3D): - if is_callable(b) or isinstance(b, Vector3D): - return AddFunction2D(a, b) + def __add__(self, object b): + if is_callable(b) or isinstance(b, Vector3D): + return AddFunction2D(self, b) return NotImplemented - def __sub__(object a, object b): - if is_callable(a) or isinstance(a, Vector3D): - if is_callable(b) or isinstance(b, Vector3D): - return SubtractFunction2D(a, b) - return NotImplemented + def __radd__(self, object a): + return self.__add__(a) - def __mul__(object a, object b): - if is_callable(a) or isinstance(a, Vector3D): - if float_is_callable(b) or isinstance(b, numbers.Real): - return MultiplyFunction2D(a, b) + def __sub__(self, object b): if is_callable(b) or isinstance(b, Vector3D): - if float_is_callable(a) or isinstance(a, numbers.Real): - return MultiplyFunction2D(b, a) + return SubtractFunction2D(self, b) return NotImplemented - - def __truediv__(object a, object b): + + def __rsub__(self, object a): if is_callable(a) or isinstance(a, Vector3D): - if float_is_callable(b) or isinstance(b, numbers.Real): - return DivideFunction2D(a, b) + return SubtractFunction2D(a, self) + + def __mul__(self, object b): + if float_is_callable(b) or isinstance(b, numbers.Real): + return MultiplyFunction2D(self, b) + return NotImplemented + + def __rmul__(self, object a): + return self.__mul__(a) + + def __truediv__(self, object b): + if float_is_callable(b) or isinstance(b, numbers.Real): + return DivideFunction2D(self, b) return NotImplemented def __neg__(self): diff --git a/raysect/core/math/function/vector3d/function3d/base.pyx b/raysect/core/math/function/vector3d/function3d/base.pyx index 453a49bd..83e0f604 100644 --- a/raysect/core/math/function/vector3d/function3d/base.pyx +++ b/raysect/core/math/function/vector3d/function3d/base.pyx @@ -66,31 +66,34 @@ cdef class Function3D(Vector3DFunction): """ return self.evaluate(x, y, z) - def __add__(object a, object b): - if is_callable(a) or isinstance(a, Vector3D): - if is_callable(b) or isinstance(b, Vector3D): - return AddFunction3D(a, b) + def __add__(self, object b): + if is_callable(b) or isinstance(b, Vector3D): + return AddFunction3D(self, b) return NotImplemented - def __sub__(object a, object b): - if is_callable(a) or isinstance(a, Vector3D): - if is_callable(b) or isinstance(b, Vector3D): - return SubtractFunction3D(a, b) - return NotImplemented + def __radd__(self, object a): + return self.__add__(a) - def __mul__(object a, object b): - if is_callable(a) or isinstance(a, Vector3D): - if float_is_callable(b) or isinstance(b, numbers.Real): - return MultiplyFunction3D(a, b) + def __sub__(self, object b): if is_callable(b) or isinstance(b, Vector3D): - if float_is_callable(a) or isinstance(a, numbers.Real): - return MultiplyFunction3D(b, a) + return SubtractFunction3D(self, b) return NotImplemented - def __truediv__(object a, object b): + def __rsub__(self, object a): if is_callable(a) or isinstance(a, Vector3D): - if float_is_callable(b) or isinstance(b, numbers.Real): - return DivideFunction3D(a, b) + return SubtractFunction3D(a, self) + + def __mul__(self, object b): + if float_is_callable(b) or isinstance(b, numbers.Real): + return MultiplyFunction3D(self, b) + return NotImplemented + + def __rmul__(self, object a): + return self.__mul__(a) + + def __truediv__(self, object b): + if float_is_callable(b) or isinstance(b, numbers.Real): + return DivideFunction3D(self, b) return NotImplemented def __neg__(self): diff --git a/raysect/core/math/normal.pxd b/raysect/core/math/normal.pxd index ab43f7af..98ceb5bc 100644 --- a/raysect/core/math/normal.pxd +++ b/raysect/core/math/normal.pxd @@ -33,30 +33,20 @@ from raysect.core.math._vec3 cimport _Vec3 from raysect.core.math.vector cimport Vector3D from raysect.core.math.affinematrix cimport AffineMatrix3D + cdef class Normal3D(_Vec3): cpdef Vector3D cross(self, _Vec3 v) - cpdef Normal3D normalise(self) - cpdef Normal3D transform(self, AffineMatrix3D m) - cpdef Normal3D transform_with_inverse(self, AffineMatrix3D m) - cdef Normal3D neg(self) - cdef Normal3D add(self, _Vec3 v) - cdef Normal3D sub(self, _Vec3 v) - cdef Normal3D mul(self, double m) - cdef Normal3D div(self, double m) - cpdef Normal3D copy(self) - cpdef Vector3D as_vector(self) - cpdef Vector3D orthogonal(self) diff --git a/raysect/core/math/normal.pyx b/raysect/core/math/normal.pyx index 78169bc0..eb680668 100644 --- a/raysect/core/math/normal.pyx +++ b/raysect/core/math/normal.pyx @@ -34,6 +34,7 @@ cimport cython from libc.math cimport sqrt, fabs from raysect.core.math.vector cimport new_vector3d + cdef class Normal3D(_Vec3): """ Represents a normal vector in 3D affine space. @@ -83,117 +84,103 @@ cdef class Normal3D(_Vec3): return self.x == n.x and self.y == n.y and self.z == n.z elif op == 3: # __ne__() return self.x != n.x or self.y != n.y or self.z != n.z - else: - return NotImplemented + + return NotImplemented def __neg__(self): """Returns a normal with the reverse orientation.""" - return new_normal3d(-self.x, - -self.y, - -self.z) + return new_normal3d(-self.x, -self.y, -self.z) - def __add__(object x, object y): + def __add__(self, object y): """Addition operator.""" - cdef _Vec3 vx, vy - - if isinstance(x, _Vec3) and isinstance(y, _Vec3): + cdef _Vec3 v - vx = <_Vec3>x - vy = <_Vec3>y + if isinstance(y, _Vec3): + v = <_Vec3>y + return new_normal3d(self.x + v.x, self.y + v.y, self.z + v.z) - return new_normal3d(vx.x + vy.x, - vx.y + vy.y, - vx.z + vy.z) - - else: - - return NotImplemented + return NotImplemented + + def __radd__(self, object x): + """Reverse addition operator.""" + + return self.__add__(x) - def __sub__(object x, object y): + def __sub__(self, object y): """Subtract operator.""" - cdef _Vec3 vx, vy + cdef _Vec3 v - if isinstance(x, _Vec3) and isinstance(y, _Vec3): + if isinstance(y, _Vec3): + v = <_Vec3>y + return new_normal3d(self.x - v.x, self.y - v.y, self.z - v.z) - vx = <_Vec3>x - vy = <_Vec3>y + return NotImplemented - return new_normal3d(vx.x - vy.x, - vx.y - vy.y, - vx.z - vy.z) + def __rsub__(self, object x): + """Reverse subtract operator.""" - else: + cdef _Vec3 v - return NotImplemented + if isinstance(x, _Vec3): + v = <_Vec3>x + return new_normal3d(v.x - self.x, v.y - self.y, v.z - self.z) + + return NotImplemented - def __mul__(object x, object y): + def __mul__(self, object y): """Multiply operator.""" cdef double s - cdef Normal3D v - cdef AffineMatrix3D m, minv - if isinstance(x, numbers.Real) and isinstance(y, Normal3D): - - s = x - v = y - - return new_normal3d(s * v.x, - s * v.y, - s * v.z) + if isinstance(y, numbers.Real): + s = y + return new_normal3d(s * self.x, s * self.y, s * self.z) - elif isinstance(x, Normal3D) and isinstance(y, numbers.Real): + return NotImplemented - s = y - v = x + def __rmul__(self, object x): + """Reverse multiply operator.""" - return new_normal3d(s * v.x, - s * v.y, - s * v.z) + cdef: + double s + AffineMatrix3D m, minv - elif isinstance(x, AffineMatrix3D) and isinstance(y, Normal3D): + if isinstance(x, numbers.Real): + s = x + return new_normal3d(s * self.x, s * self.y, s * self.z) + elif isinstance(x, AffineMatrix3D): m = x - v = y - minv = m.inverse() - return new_normal3d(minv.m[0][0] * v.x + minv.m[1][0] * v.y + minv.m[2][0] * v.z, - minv.m[0][1] * v.x + minv.m[1][1] * v.y + minv.m[2][1] * v.z, - minv.m[0][2] * v.x + minv.m[1][2] * v.y + minv.m[2][2] * v.z) + return new_normal3d( + minv.m[0][0] * self.x + minv.m[1][0] * self.y + minv.m[2][0] * self.z, + minv.m[0][1] * self.x + minv.m[1][1] * self.y + minv.m[2][1] * self.z, + minv.m[0][2] * self.x + minv.m[1][2] * self.y + minv.m[2][2] * self.z + ) - else: - - return NotImplemented + return NotImplemented @cython.cdivision(True) - def __truediv__(object x, object y): + def __truediv__(self, object y): """Division operator.""" cdef double d - cdef Normal3D v - if isinstance(x, Normal3D) and isinstance(y, numbers.Real): + if isinstance(y, numbers.Real): - d = y + d = y # prevent divide my zero if d == 0.0: - raise ZeroDivisionError("Cannot divide a vector by a zero scalar.") - v = x d = 1.0 / d + return new_normal3d(d * self.x, d * self.y, d * self.z) - return new_normal3d(d * v.x, - d * v.y, - d * v.z) - - else: - - return NotImplemented + return NotImplemented cpdef Vector3D cross(self, _Vec3 v): """ @@ -205,9 +192,11 @@ cdef class Normal3D(_Vec3): :rtype: Vector3D """ - return new_vector3d(self.y * v.z - v.y * self.z, - self.z * v.x - v.z * self.x, - self.x * v.y - v.x * self.y) + return new_vector3d( + self.y * v.z - v.y * self.z, + self.z * v.x - v.z * self.x, + self.x * v.y - v.x * self.y + ) @cython.cdivision(True) cpdef Normal3D normalise(self): @@ -224,15 +213,11 @@ cdef class Normal3D(_Vec3): # if current length is zero, problem is ill defined t = self.x * self.x + self.y * self.y + self.z * self.z if t == 0.0: - raise ZeroDivisionError("A zero length vector can not be normalised as the direction of a zero length vector is undefined.") # normalise and rescale vector t = 1.0 / sqrt(t) - - return new_normal3d(self.x * t, - self.y * t, - self.z * t) + return new_normal3d(self.x * t, self.y * t, self.z * t) cpdef Normal3D transform(self, AffineMatrix3D m): """ @@ -256,9 +241,11 @@ cdef class Normal3D(_Vec3): cdef AffineMatrix3D minv minv = m.inverse() - return new_normal3d(minv.m[0][0] * self.x + minv.m[1][0] * self.y + minv.m[2][0] * self.z, - minv.m[0][1] * self.x + minv.m[1][1] * self.y + minv.m[2][1] * self.z, - minv.m[0][2] * self.x + minv.m[1][2] * self.y + minv.m[2][2] * self.z) + return new_normal3d( + minv.m[0][0] * self.x + minv.m[1][0] * self.y + minv.m[2][0] * self.z, + minv.m[0][1] * self.x + minv.m[1][1] * self.y + minv.m[2][1] * self.z, + minv.m[0][2] * self.x + minv.m[1][2] * self.y + minv.m[2][2] * self.z + ) cpdef Normal3D transform_with_inverse(self, AffineMatrix3D m): """ @@ -277,9 +264,11 @@ cdef class Normal3D(_Vec3): :rtype: Normal3D """ - return new_normal3d(m.m[0][0] * self.x + m.m[1][0] * self.y + m.m[2][0] * self.z, - m.m[0][1] * self.x + m.m[1][1] * self.y + m.m[2][1] * self.z, - m.m[0][2] * self.x + m.m[1][2] * self.y + m.m[2][2] * self.z) + return new_normal3d( + m.m[0][0] * self.x + m.m[1][0] * self.y + m.m[2][0] * self.z, + m.m[0][1] * self.x + m.m[1][1] * self.y + m.m[2][1] * self.z, + m.m[0][2] * self.x + m.m[1][2] * self.y + m.m[2][2] * self.z + ) cdef Normal3D neg(self): """ @@ -289,9 +278,7 @@ cdef class Normal3D(_Vec3): to the equivalent python operator. """ - return new_normal3d(-self.x, - -self.y, - -self.z) + return new_normal3d(-self.x, -self.y, -self.z) cdef Normal3D add(self, _Vec3 v): """ @@ -301,9 +288,7 @@ cdef class Normal3D(_Vec3): to the equivalent python operator. """ - return new_normal3d(self.x + v.x, - self.y + v.y, - self.z + v.z) + return new_normal3d(self.x + v.x, self.y + v.y, self.z + v.z) cdef Normal3D sub(self, _Vec3 v): """ @@ -313,9 +298,7 @@ cdef class Normal3D(_Vec3): to the equivalent python operator. """ - return new_normal3d(self.x - v.x, - self.y - v.y, - self.z - v.z) + return new_normal3d(self.x - v.x, self.y - v.y, self.z - v.z) cdef Normal3D mul(self, double m): """ @@ -325,9 +308,7 @@ cdef class Normal3D(_Vec3): to the equivalent python operator. """ - return new_normal3d(self.x * m, - self.y * m, - self.z * m) + return new_normal3d(self.x * m, self.y * m, self.z * m) @cython.cdivision(True) cdef Normal3D div(self, double d): @@ -339,14 +320,10 @@ cdef class Normal3D(_Vec3): """ if d == 0.0: - raise ZeroDivisionError("Cannot divide a vector by a zero scalar.") d = 1.0 / d - - return new_normal3d(self.x * d, - self.y * d, - self.z * d) + return new_normal3d(self.x * d, self.y * d, self.z * d) cpdef Normal3D copy(self): """ @@ -355,9 +332,7 @@ cdef class Normal3D(_Vec3): :rtype: Normal3D """ - return new_normal3d(self.x, - self.y, - self.z) + return new_normal3d(self.x, self.y, self.z) cpdef Vector3D as_vector(self): """ @@ -366,9 +341,7 @@ cdef class Normal3D(_Vec3): :rtype: Vector3D """ - return new_vector3d(self.x, - self.y, - self.z) + return new_vector3d(self.x, self.y, self.z) cpdef Vector3D orthogonal(self): """ diff --git a/raysect/core/math/point.pxd b/raysect/core/math/point.pxd index ced255a4..3665eecc 100644 --- a/raysect/core/math/point.pxd +++ b/raysect/core/math/point.pxd @@ -38,19 +38,12 @@ cdef class Point3D: cdef public double x, y, z cpdef Vector3D vector_to(self, Point3D p) - cpdef double distance_to(self, Point3D p) - cpdef Point3D transform(self, AffineMatrix3D m) - cdef Point3D add(self, _Vec3 v) - cdef Point3D sub(self, _Vec3 v) - cpdef Point3D copy(self) - cdef double get_index(self, int index) nogil - cdef void set_index(self, int index, double value) nogil @@ -75,19 +68,12 @@ cdef class Point2D: cdef public double x, y cpdef Vector2D vector_to(self, Point2D p) - cpdef double distance_to(self, Point2D p) - # cpdef Point3D transform(self, AffineMatrix3D m) - cdef Point2D add(self, Vector2D v) - cdef Point2D sub(self, Vector2D v) - cpdef Point2D copy(self) - cdef double get_index(self, int index) nogil - cdef void set_index(self, int index, double value) nogil diff --git a/raysect/core/math/point.pyx b/raysect/core/math/point.pyx index 20294b54..c1089302 100644 --- a/raysect/core/math/point.pyx +++ b/raysect/core/math/point.pyx @@ -59,7 +59,6 @@ cdef class Point3D: >>> from raysect.core import Point3D >>> a = Point3D(0, 1, 2) - """ def __init__(self, double x=0.0, double y=0.0, double z=0.0): @@ -86,8 +85,8 @@ cdef class Point3D: return self.x == p.x and self.y == p.y and self.z == p.z elif op == 3: # __ne__() return self.x != p.x or self.y != p.y or self.z != p.z - else: - return NotImplemented + + return NotImplemented def __getitem__(self, int i): """Returns the point coordinates by index ([0,1,2] -> [x,y,z]). @@ -103,8 +102,8 @@ cdef class Point3D: return self.y elif i == 2: return self.z - else: - raise IndexError("Index out of range [0, 2].") + + raise IndexError("Index out of range [0, 2].") def __setitem__(self, int i, double value): """Sets the point coordinates by index ([0,1,2] -> [x,y,z]). @@ -132,87 +131,70 @@ cdef class Point3D: >>> x, y, z (0.0, 1.0, 2.0) """ + yield self.x yield self.y yield self.z - def __add__(object x, object y): + def __add__(self, object y): """Addition operator. >>> Point3D(1, 0, 0) + Vector3D(0, 1, 0) Point3D(1.0, 1.0, 0.0) """ - cdef Point3D p cdef _Vec3 v - if isinstance(x, Point3D) and isinstance(y, _Vec3): - - p = x - v = <_Vec3>y - - else: - - return NotImplemented + if isinstance(y, _Vec3): + v = <_Vec3> y + return new_point3d(self.x + v.x, self.y + v.y, self.z + v.z) - return new_point3d(p.x + v.x, - p.y + v.y, - p.z + v.z) + return NotImplemented - def __sub__(object x, object y): + def __sub__(self, object y): """Subtraction operator. >>> Point3D(1, 0, 0) - Vector3D(0, 1, 0) Point3D(1.0, -1.0, 0.0) """ - cdef Point3D p cdef _Vec3 v - if isinstance(x, Point3D) and isinstance(y, _Vec3): - - p = x + if isinstance(y, _Vec3): v = <_Vec3>y + return new_point3d(self.x - v.x, self.y - v.y, self.z - v.z) - return new_point3d(p.x - v.x, - p.y - v.y, - p.z - v.z) - - else: - - return NotImplemented + return NotImplemented @cython.cdivision(True) - def __mul__(object x, object y): + def __rmul__(self, object x): """Multiplication operator. :param AffineMatrix3D x: transformation matrix x - :param Point3D y: point to transform :return: Matrix multiplication of a 3D transformation matrix with the input point. :rtype: Point3D """ cdef AffineMatrix3D m - cdef Point3D v cdef double w - if isinstance(x, AffineMatrix3D) and isinstance(y, Point3D): + if isinstance(x, AffineMatrix3D): m = x - v = y # 4th element of homogeneous coordinate - w = m.m[3][0] * v.x + m.m[3][1] * v.y + m.m[3][2] * v.z + m.m[3][3] + w = m.m[3][0] * self.x + m.m[3][1] * self.y + m.m[3][2] * self.z + m.m[3][3] if w == 0.0: - raise ZeroDivisionError("Bad matrix transform, 4th element of homogeneous coordinate is zero.") # pre divide for speed (dividing is much slower than multiplying) w = 1.0 / w - return new_point3d((m.m[0][0] * v.x + m.m[0][1] * v.y + m.m[0][2] * v.z + m.m[0][3]) * w, - (m.m[1][0] * v.x + m.m[1][1] * v.y + m.m[1][2] * v.z + m.m[1][3]) * w, - (m.m[2][0] * v.x + m.m[2][1] * v.y + m.m[2][2] * v.z + m.m[2][3]) * w) + return new_point3d( + (m.m[0][0] * self.x + m.m[0][1] * self.y + m.m[0][2] * self.z + m.m[0][3]) * w, + (m.m[1][0] * self.x + m.m[1][1] * self.y + m.m[1][2] * self.z + m.m[1][3]) * w, + (m.m[2][0] * self.x + m.m[2][1] * self.y + m.m[2][2] * self.z + m.m[2][3]) * w + ) return NotImplemented @@ -244,9 +226,7 @@ cdef class Point3D: """ - return new_vector3d(p.x - self.x, - p.y - self.y, - p.z - self.z) + return new_vector3d(p.x - self.x, p.y - self.y, p.z - self.z) cpdef double distance_to(self, Point3D p): """ @@ -293,15 +273,15 @@ cdef class Point3D: # 4th element of homogeneous coordinate w = m.m[3][0] * self.x + m.m[3][1] * self.y + m.m[3][2] * self.z + m.m[3][3] if w == 0.0: - raise ZeroDivisionError("Bad matrix transform, 4th element of homogeneous coordinate is zero.") # pre divide for speed (dividing is much slower than multiplying) w = 1.0 / w - - return new_point3d((m.m[0][0] * self.x + m.m[0][1] * self.y + m.m[0][2] * self.z + m.m[0][3]) * w, - (m.m[1][0] * self.x + m.m[1][1] * self.y + m.m[1][2] * self.z + m.m[1][3]) * w, - (m.m[2][0] * self.x + m.m[2][1] * self.y + m.m[2][2] * self.z + m.m[2][3]) * w) + return new_point3d( + (m.m[0][0] * self.x + m.m[0][1] * self.y + m.m[0][2] * self.z + m.m[0][3]) * w, + (m.m[1][0] * self.x + m.m[1][1] * self.y + m.m[1][2] * self.z + m.m[1][3]) * w, + (m.m[2][0] * self.x + m.m[2][1] * self.y + m.m[2][2] * self.z + m.m[2][3]) * w + ) cdef Point3D add(self, _Vec3 v): """ @@ -311,9 +291,7 @@ cdef class Point3D: to the equivalent python operator. """ - return new_point3d(self.x + v.x, - self.y + v.y, - self.z + v.z) + return new_point3d(self.x + v.x, self.y + v.y, self.z + v.z) cdef Point3D sub(self, _Vec3 v): """ @@ -323,9 +301,7 @@ cdef class Point3D: to the equivalent python operator. """ - return new_point3d(self.x - v.x, - self.y - v.y, - self.z - v.z) + return new_point3d(self.x - v.x, self.y - v.y, self.z - v.z) cpdef Point3D copy(self): """ @@ -340,9 +316,7 @@ cdef class Point3D: Point3D(0.0, 1.0, 2.0) """ - return new_point3d(self.x, - self.y, - self.z) + return new_point3d(self.x, self.y, self.z) cdef double get_index(self, int index) nogil: """ @@ -404,7 +378,6 @@ cdef class Point2D: """ def __init__(self, double x=0.0, double y=0.0): - self.x = x self.y = y @@ -441,8 +414,8 @@ cdef class Point2D: return self.x elif i == 1: return self.y - else: - raise IndexError("Index out of range [0, 1].") + + raise IndexError("Index out of range [0, 1].") def __setitem__(self, int i, double value): """Sets the point coordinates by index ([0,1] -> [x,y]). @@ -469,67 +442,54 @@ cdef class Point2D: (1.0, 1.0) """ + yield self.x yield self.y - def __add__(object x, object y): + def __add__(self, object y): """Addition operator. >>> Point2D(1, 0) + Vector2D(0, 1) Point2D(1.0, 1.0) """ - cdef Point2D p cdef Vector2D v - if isinstance(x, Point2D) and isinstance(y, Vector2D): - - p = x + if isinstance(y, Vector2D): v = y + return new_point2d(self.x + v.x, self.y + v.y) - else: - - return NotImplemented - - return new_point2d(p.x + v.x, p.y + v.y) + return NotImplemented - def __sub__(object x, object y): + def __sub__(self, object y): """Subtraction operator. >>> Point2D(1, 0) - Vector2D(0, 1) Point2D(1.0, -1.0) """ - cdef Point2D p cdef Vector2D v - if isinstance(x, Point2D) and isinstance(y, Vector2D): - - p = x + if isinstance(y, Vector2D): v = y + return new_point2d(self.x - v.x, self.y - v.y) - return new_point2d(p.x - v.x, p.y - v.y) - - else: - - return NotImplemented + return NotImplemented @cython.cdivision(True) - def __mul__(object x, object y): + def __rmul__(self, object x): """Multiply operator.""" - raise NotImplemented + return NotImplemented # cdef AffineMatrix3D m - # cdef Point3D v # cdef double w # - # if isinstance(x, AffineMatrix3D) and isinstance(y, Point3D): + # if isinstance(x, AffineMatrix3D): # # m = x - # v = y # # # 4th element of homogeneous coordinate - # w = m.m[3][0] * v.x + m.m[3][1] * v.y + m.m[3][2] * v.z + m.m[3][3] + # w = m.m[3][0] * self.x + m.m[3][1] * self.y + m.m[3][2] * self.z + m.m[3][3] # if w == 0.0: # # raise ZeroDivisionError("Bad matrix transform, 4th element of homogeneous coordinate is zero.") @@ -537,9 +497,9 @@ cdef class Point2D: # # pre divide for speed (dividing is much slower than multiplying) # w = 1.0 / w # - # return new_point3d((m.m[0][0] * v.x + m.m[0][1] * v.y + m.m[0][2] * v.z + m.m[0][3]) * w, - # (m.m[1][0] * v.x + m.m[1][1] * v.y + m.m[1][2] * v.z + m.m[1][3]) * w, - # (m.m[2][0] * v.x + m.m[2][1] * v.y + m.m[2][2] * v.z + m.m[2][3]) * w) + # return new_point3d((m.m[0][0] * self.x + m.m[0][1] * self.y + m.m[0][2] * self.z + m.m[0][3]) * w, + # (m.m[1][0] * self.x + m.m[1][1] * self.y + m.m[1][2] * self.z + m.m[1][3]) * w, + # (m.m[2][0] * self.x + m.m[2][1] * self.y + m.m[2][2] * self.z + m.m[2][3]) * w) # # return NotImplemented diff --git a/raysect/core/math/quaternion.pxd b/raysect/core/math/quaternion.pxd index bb6d1339..2c9ec72e 100644 --- a/raysect/core/math/quaternion.pxd +++ b/raysect/core/math/quaternion.pxd @@ -38,41 +38,23 @@ cdef class Quaternion: cdef public double x, y, z, s cpdef Quaternion copy(self) - cpdef Quaternion conjugate(self) - cpdef Quaternion inverse(self) - cpdef Quaternion normalise(self) - cpdef bint is_unit(self, double tolerance=*) - cpdef Quaternion transform(self, AffineMatrix3D m) - cpdef AffineMatrix3D as_matrix(self) - cpdef Quaternion quaternion_to(self, Quaternion q) - cdef Quaternion neg(self) - cdef Quaternion add(self, Quaternion q) - cdef Quaternion sub(self, Quaternion q) - cdef Quaternion mul_quaternion(self, Quaternion q) - cdef Quaternion mul_scalar(self, double d) - cdef Quaternion div_quaternion(self, Quaternion q) - cdef Quaternion div_scalar(self, double d) - cdef Vector3D get_axis(self) - cdef double get_angle(self) - cdef double get_length(self) nogil - cdef object set_length(self, double v) @@ -96,4 +78,4 @@ cdef inline Quaternion new_quaternion(double x, double y, double z, double s): cdef Quaternion new_quaternion_from_matrix(AffineMatrix3D matrix) -cdef Quaternion new_quaternion_from_axis_angle(Vector3D axis, double angle) \ No newline at end of file +cdef Quaternion new_quaternion_from_axis_angle(Vector3D axis, double angle) diff --git a/raysect/core/math/quaternion.pyx b/raysect/core/math/quaternion.pyx index 26d85e06..dec92228 100644 --- a/raysect/core/math/quaternion.pyx +++ b/raysect/core/math/quaternion.pyx @@ -31,13 +31,14 @@ import numbers cimport cython -from libc.math cimport sqrt, sin, cos, asin, acos, atan2, fabs, M_PI, copysign +from libc.math cimport sqrt, sin, cos, acos, fabs from raysect.core.math.vector cimport new_vector3d from raysect.core.math.affinematrix cimport new_affinematrix3d, AffineMatrix3D -DEF RAD2DEG = 57.29577951308232000 # 180 / pi -DEF DEG2RAD = 0.017453292519943295 # pi / 180 + +cdef const double RAD2DEG = 57.29577951308232000 # 180 / pi +cdef const double DEG2RAD = 0.017453292519943295 # pi / 180 cdef class Quaternion: @@ -72,8 +73,8 @@ cdef class Quaternion: return self.z elif i == 3: return self.s - else: - raise IndexError("Index out of range [0, 3].") + + raise IndexError("Index out of range [0, 3].") def __setitem__(self, int i, double value): """Sets the quaternion coordinates by index ([0,1,2,3] -> [x,y,z,s]). @@ -105,6 +106,7 @@ cdef class Quaternion: >>> x, y, z, s (0.0, 1.0, 2.0, 3.0) """ + yield self.x yield self.y yield self.z @@ -127,7 +129,7 @@ cdef class Quaternion: return new_quaternion(-self.x, -self.y, -self.z, -self.s) - def __eq__(object x, object y): + def __eq__(self, object y): """ Equality operator. @@ -137,18 +139,15 @@ cdef class Quaternion: True """ - cdef Quaternion q1, q2 - - if isinstance(x, Quaternion) and isinstance(y, Quaternion): + cdef Quaternion q - q1 = x - q2 = y - return q1.x == q2.x and q1.y == q2.y and q1.z == q2.z and q1.s == q2.s + if isinstance(y, Quaternion): + q = y + return self.x == q.x and self.y == q.y and self.z == q.z and self.s == q.s - else: - raise TypeError('A quaternion can only be equality tested against another quaternion.') + raise TypeError('A quaternion can only be equality tested against another quaternion.') - def __add__(object x, object y): + def __add__(self, object y): """ Addition operator. @@ -158,18 +157,15 @@ cdef class Quaternion: Quaternion(0.0, 1.0, 0.0, 1.0) """ - cdef Quaternion q1, q2 + cdef Quaternion q - if isinstance(x, Quaternion) and isinstance(y, Quaternion): + if isinstance(y, Quaternion): + q = y + return new_quaternion(self.x + q.x, self.y + q.y, self.z + q.z, self.s + q.s) - q1 = x - q2 = y - return new_quaternion(q1.x + q2.x, q1.y + q2.y, q1.z + q2.z, q1.s + q2.s) + return NotImplemented - else: - return NotImplemented - - def __sub__(object x, object y): + def __sub__(self, object y): """Subtraction operator. .. code-block:: pycon @@ -178,18 +174,15 @@ cdef class Quaternion: Quaternion(0.0, -1.0, 0, 1.0) """ - cdef Quaternion q1, q2 + cdef Quaternion q - if isinstance(x, Quaternion) and isinstance(y, Quaternion): + if isinstance(y, Quaternion): + q = y + return new_quaternion(self.x - q.x, self.y - q.y, self.z - q.z, self.s - q.s) - q1 = x - q2 = y - return new_quaternion(q1.x - q2.x, q1.y - q2.y, q1.z - q2.z, q1.s - q2.s) + return NotImplemented - else: - return NotImplemented - - def __mul__(object x, object y): + def __mul__(self, object y): """Multiplication operator. .. code-block:: pycon @@ -201,31 +194,37 @@ cdef class Quaternion: """ cdef double s - cdef Quaternion q1, q2 + cdef Quaternion q + + if isinstance(y, numbers.Real): + s = y + return self.mul_scalar(s) - if isinstance(x, numbers.Real) and isinstance(y, Quaternion): + elif isinstance(y, Quaternion): + q = y + return self.mul_quaternion(q) - s = x - q1 = y - return q1.mul_scalar(s) + return NotImplemented - elif isinstance(x, Quaternion) and isinstance(y, numbers.Real): + def __rmul__(self, object x): + """Reverse multiplication operator. - q1 = x - s = y - return q1.mul_scalar(s) + .. code-block:: pycon - elif isinstance(x, Quaternion) and isinstance(y, Quaternion): + >>> 2 * Quaternion(0, 0, 1, 1) + Quaternion(0.0, 0.0, 2.0, 2.0) + """ - q1 = x - q2 = y - return q1.mul_quaternion(q2) + cdef double s - else: - return NotImplemented() + if isinstance(x, numbers.Real): + s = x + return self.mul_scalar(s) + + return NotImplemented @cython.cdivision(True) - def __truediv__(object x, object y): + def __truediv__(self, object y): """Division operator. .. code-block:: pycon @@ -237,22 +236,17 @@ cdef class Quaternion: """ cdef double d - cdef Quaternion q1, q2, q2_inv - - if isinstance(x, Quaternion) and isinstance(y, numbers.Real): + cdef Quaternion q + if isinstance(y, numbers.Real): d = y - q1 = x - return q1.div_scalar(d) + return self.div_scalar(d) - elif isinstance(x, Quaternion) and isinstance(y, Quaternion): + elif isinstance(y, Quaternion): + q = y + return self.div_quaternion(q) - q1 = x - q2 = y - return q1.div_quaternion(q2) - - else: - raise TypeError('Unsupported operand type. Expects a real number.') + raise TypeError('Unsupported operand type. Expects a real number.') @property def length(self): @@ -264,34 +258,40 @@ cdef class Quaternion: >>> Quaternion(1, 2, 3, 0).length 3.7416573867739413 """ + return self.get_length() @length.setter def length(self, value): self.set_length(value) - @property def axis(self): """ The axis around which this quaternion rotates. """ + return self.get_axis() @property def angle(self): - """The magnitude of rotation around this quaternion's rotation axis in degrees.""" + """ + The magnitude of rotation around this quaternion's rotation axis in degrees. + """ + return self.get_angle() cpdef Quaternion copy(self): - """Returns a copy of this quaternion.""" + """ + Returns a copy of this quaternion. + """ return new_quaternion(self.x, self.y, self.z, self.s) cpdef Quaternion conjugate(self): """ - Complex conjugate operator. - + Complex conjugate operator. + .. code-block:: pycon >>> Quaternion(1, 2, 3, 0).conjugate() @@ -322,7 +322,7 @@ cdef class Quaternion: The returned quaternion is normalised to have norm length 1.0 - a unit quaternion. .. code-block:: pycon - + >>> a = Quaternion(1, 2, 3, 0) >>> a.normalise() Quaternion(0.26726, 0.53452, 0.80178, 0.0) @@ -345,6 +345,7 @@ cdef class Quaternion: :param float tolerance: The numerical tolerance by which the quaternion norm can differ by 1.0. """ + return fabs(1.0 - self.get_length()) <= tolerance cpdef Quaternion transform(self, AffineMatrix3D m): @@ -409,10 +410,12 @@ cdef class Quaternion: m21 = 2*qy*qz + 2*qx*qs m22 = 1 - 2*qx2 - 2*qy2 - return new_affinematrix3d(m00, m01, m02, 0, - m10, m11, m12, 0, - m20, m21, m22, 0, - 0, 0, 0, 1) + return new_affinematrix3d( + m00, m01, m02, 0, + m10, m11, m12, 0, + m20, m21, m22, 0, + 0, 0, 0, 1 + ) cpdef Quaternion quaternion_to(self, Quaternion q): """ @@ -421,12 +424,12 @@ cdef class Quaternion: This method calculates the quaternion required to map this quaternion onto the supplied quaternion. Both quaternions will be normalised and a normalised quaternion will be returned. - + .. code-block:: pycon - + >>> from raysect.core.math import Quaternion >>> - >>> q1 = Quaternion.from_axis_angle(Vector3D(1,0,0), 10) + >>> q1 = Quaternion.from_axis_angle(Vector3D(1,0,0), 10) >>> q2 = Quaternion.from_axis_angle(Vector3D(1,0,0), 25) >>> d = q1.quaternion_to(q2) >>> d @@ -435,9 +438,9 @@ cdef class Quaternion: 15.000000000000027 >>> d.axis Vector3D(1.0, 0.0, 0.0) - + :param Quaternion q: The target quaternion. - :return: A new Quaternion object representing the specified rotation. + :return: A new Quaternion object representing the specified rotation. """ return q.normalise().mul_quaternion(self.normalise().conjugate()).normalise() @@ -565,7 +568,9 @@ cdef class Quaternion: @cython.cdivision(True) cdef double get_angle(self): - """The magnitude of rotation around this quaternion's rotation axis in degrees.""" + """ + The magnitude of rotation around this quaternion's rotation axis in degrees. + """ cdef Quaternion q = self.normalise() return 2 * acos(q.s) * RAD2DEG @@ -578,6 +583,7 @@ cdef class Quaternion: Use instead of Python attribute access in cython code. """ + return sqrt(self.x * self.x + self.y * self.y + self.z * self.z + self.s * self.s) @cython.cdivision(True) @@ -686,4 +692,4 @@ cdef Quaternion new_quaternion_from_axis_angle(Vector3D axis, double angle): qz = axis.z * sin(theta_2) qs = cos(theta_2) - return new_quaternion(qx, qy, qz, qs) \ No newline at end of file + return new_quaternion(qx, qy, qz, qs) diff --git a/raysect/core/math/random.pyx b/raysect/core/math/random.pyx index c47c024a..6029fbbf 100644 --- a/raysect/core/math/random.pyx +++ b/raysect/core/math/random.pyx @@ -96,8 +96,9 @@ from libc.math cimport cos, sin, asin, log, fabs, sqrt, M_PI as PI from libc.stdint cimport uint64_t, int64_t cimport cython -DEF NN = 312 -DEF MM = 156 +cdef enum: + NN = 312 + MM = 156 # The array for the state vector cdef uint64_t mt[NN] diff --git a/raysect/core/math/sampler/__init__.pxd b/raysect/core/math/sampler/__init__.pxd index 8931c646..98840aff 100644 --- a/raysect/core/math/sampler/__init__.pxd +++ b/raysect/core/math/sampler/__init__.pxd @@ -31,5 +31,5 @@ from raysect.core.math.sampler.solidangle cimport * from raysect.core.math.sampler.surface3d cimport * -from raysect.core.math.sampler.targetted cimport * +from raysect.core.math.sampler.targeted cimport * diff --git a/raysect/core/math/sampler/__init__.py b/raysect/core/math/sampler/__init__.py index 0181770d..39b18ae2 100644 --- a/raysect/core/math/sampler/__init__.py +++ b/raysect/core/math/sampler/__init__.py @@ -31,4 +31,4 @@ from .solidangle import * from .surface3d import * -from .targetted import * +from .targeted import * diff --git a/raysect/core/math/sampler/meson.build b/raysect/core/math/sampler/meson.build index c6b49aff..6fea0bca 100644 --- a/raysect/core/math/sampler/meson.build +++ b/raysect/core/math/sampler/meson.build @@ -5,8 +5,8 @@ target_path = 'raysect/core/math/sampler' # source files py_files = ['__init__.py'] -pyx_files = ['solidangle.pyx', 'surface3d.pyx', 'targetted.pyx'] -pxd_files = ['__init__.pxd', 'solidangle.pxd', 'surface3d.pxd', 'targetted.pxd'] +pyx_files = ['solidangle.pyx', 'surface3d.pyx', 'targeted.pyx'] +pxd_files = ['__init__.pxd', 'solidangle.pxd', 'surface3d.pxd', 'targeted.pxd'] data_files = [] # compile cython diff --git a/raysect/core/math/sampler/solidangle.pxd b/raysect/core/math/sampler/solidangle.pxd index ca9343e9..8a3131c3 100644 --- a/raysect/core/math/sampler/solidangle.pxd +++ b/raysect/core/math/sampler/solidangle.pxd @@ -29,25 +29,15 @@ # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. -from libc.math cimport M_PI, M_1_PI, sqrt, sin, cos -from raysect.core.math cimport Point2D, new_point2d, Point3D, new_point3d, Vector3D, new_vector3d -from raysect.core.math.random cimport uniform -from raysect.core.math.cython cimport barycentric_coords, barycentric_interpolation - -DEF R_2_PI = 0.15915494309189535 # 1 / (2 * pi) -DEF R_4_PI = 0.07957747154594767 # 1 / (4 * pi) +from raysect.core.math cimport Vector3D cdef class SolidAngleSampler: cpdef double pdf(self, Vector3D sample) - cdef Vector3D sample(self) - cdef tuple sample_with_pdf(self) - cdef list samples(self, int samples) - cdef list samples_with_pdfs(self, int samples) diff --git a/raysect/core/math/sampler/solidangle.pyx b/raysect/core/math/sampler/solidangle.pyx index 04068a42..c6a29719 100644 --- a/raysect/core/math/sampler/solidangle.pyx +++ b/raysect/core/math/sampler/solidangle.pyx @@ -29,14 +29,14 @@ # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. -from libc.math cimport M_PI, M_1_PI, sqrt, sin, cos, asin +from libc.math cimport M_PI, M_1_PI, sqrt, sin, cos from raysect.core.math cimport Vector3D, new_vector3d from raysect.core.math.random cimport uniform -# TODO: add tests - idea: solve the lighting equation with a uniform emitting surface with each sampler and check the mean radiance is unity -DEF R_2_PI = 0.15915494309189535 # 1 / (2 * pi) -DEF R_4_PI = 0.07957747154594767 # 1 / (4 * pi) +# TODO: add tests - idea: solve the lighting equation with a uniform emitting surface with each sampler and check the mean radiance is unity +cdef const double R_2_PI = 0.15915494309189535 # 1 / (2 * pi) +cdef const double R_4_PI = 0.07957747154594767 # 1 / (4 * pi) cdef class SolidAngleSampler: @@ -73,7 +73,7 @@ cdef class SolidAngleSampler: cpdef double pdf(self, Vector3D sample): """ Generates a pdf for a given sample value. - + Vectors *must* be normalised. :param Vector3D sample: The sample point at which to get the pdf. @@ -242,7 +242,7 @@ cdef class ConeUniformSampler(SolidAngleSampler): Generates a uniform weighted random vector from a cone. The cone is aligned along the z-axis. - + :param angle: Angle of the cone in degrees (default=45). .. code-block:: pycon diff --git a/raysect/core/math/sampler/surface3d.pxd b/raysect/core/math/sampler/surface3d.pxd index 72dc91bf..3b007989 100644 --- a/raysect/core/math/sampler/surface3d.pxd +++ b/raysect/core/math/sampler/surface3d.pxd @@ -35,11 +35,8 @@ from raysect.core.math cimport Point3D cdef class SurfaceSampler3D: cdef Point3D sample(self) - cdef tuple sample_with_pdf(self) - cdef list samples(self, int samples) - cdef list samples_with_pdfs(self, int samples) diff --git a/raysect/core/math/sampler/targetted.pxd b/raysect/core/math/sampler/targeted.pxd similarity index 93% rename from raysect/core/math/sampler/targetted.pxd rename to raysect/core/math/sampler/targeted.pxd index a938a571..b5633b4d 100644 --- a/raysect/core/math/sampler/targetted.pxd +++ b/raysect/core/math/sampler/targeted.pxd @@ -32,7 +32,8 @@ cimport numpy as np from raysect.core.math cimport Point3D, Vector3D -cdef class _TargettedSampler: + +cdef class _TargetedSampler: cdef: double _total_weight @@ -41,24 +42,17 @@ cdef class _TargettedSampler: double[::1] _cdf_mv cdef object _validate_targets(self) - cpdef double pdf(self, Point3D point, Vector3D sample) - cdef Vector3D sample(self, Point3D point) - cdef tuple sample_with_pdf(self, Point3D point) - cdef list samples(self, Point3D point, int samples) - cdef list samples_with_pdfs(self, Point3D point, int samples) - cdef object _calculate_cdf(self) - cdef tuple _pick_sphere(self) -cdef class TargettedHemisphereSampler(_TargettedSampler): +cdef class TargetedHemisphereSampler(_TargetedSampler): pass -cdef class TargettedSphereSampler(_TargettedSampler): +cdef class TargetedSphereSampler(_TargetedSampler): pass \ No newline at end of file diff --git a/raysect/core/math/sampler/targetted.pyx b/raysect/core/math/sampler/targeted.pyx similarity index 97% rename from raysect/core/math/sampler/targetted.pyx rename to raysect/core/math/sampler/targeted.pyx index 3aabb070..9dffa327 100644 --- a/raysect/core/math/sampler/targetted.pyx +++ b/raysect/core/math/sampler/targeted.pyx @@ -34,11 +34,11 @@ cimport numpy as np from raysect.core.math cimport Point3D, Vector3D, AffineMatrix3D from raysect.core.math.random cimport uniform, vector_sphere, vector_cone_uniform, vector_hemisphere_cosine from raysect.core.math.cython cimport find_index, rotate_basis -from libc.math cimport M_PI, M_1_PI, asin, acos, sqrt +from libc.math cimport M_PI, M_1_PI, asin, sqrt cimport cython -cdef class _TargettedSampler: +cdef class _TargetedSampler: def __init__(self, object targets): """ @@ -46,7 +46,7 @@ cdef class _TargettedSampler: Each target tuple consists of (Point3D sphere_centre, double sphere_radius, double weight). - :param list targets: A list of tuples describing spheres for targetted sampling. + :param list targets: A list of tuples describing spheres for targeted sampling. """ self._targets = tuple(targets) @@ -248,16 +248,16 @@ cdef class _TargettedSampler: return self._targets[index] -cdef class TargettedHemisphereSampler(_TargettedSampler): +cdef class TargetedHemisphereSampler(_TargetedSampler): """ - Generates vectors on a hemisphere targetting a set of target spheres. + Generates vectors on a hemisphere targeting a set of target spheres. This sampler takes a list of spheres and corresponding weighting factors. To generate a sample a sphere is randomly selected according the distribution of the sphere weights and launches a ray at the solid angle subtended by the target sphere. - If the targetted sphere intersects with the hemisphere's base plane, lies + If the targeted sphere intersects with the hemisphere's base plane, lies behind the plane, or the origin point from which samples are generated lies inside the target sphere, a sample is randomly selected from the full hemisphere using a cosine weighted distribution. @@ -389,9 +389,9 @@ cdef class TargettedHemisphereSampler(_TargettedSampler): return sample.transform(rotation) -cdef class TargettedSphereSampler(_TargettedSampler): +cdef class TargetedSphereSampler(_TargetedSampler): """ - Generates vectors targetting a set of target spheres. + Generates vectors targeting a set of target spheres. This sampler takes a list of spheres and corresponding weighting factors. To generate a sample a sphere is randomly selected according the diff --git a/raysect/core/math/spatial/kdtree2d.pxd b/raysect/core/math/spatial/kdtree2d.pxd index 1ebfb847..c2c7e414 100644 --- a/raysect/core/math/spatial/kdtree2d.pxd +++ b/raysect/core/math/spatial/kdtree2d.pxd @@ -30,10 +30,10 @@ # POSSIBILITY OF SUCH DAMAGE. from raysect.core.boundingbox cimport BoundingBox2D -from raysect.core.ray cimport Ray from raysect.core.math.point cimport Point2D from libc.stdint cimport int32_t + # c-structure that represent a kd-tree node cdef struct kdnode: @@ -69,56 +69,32 @@ cdef class KDTree2DCore: double _empty_bonus cdef int32_t _build(self, list items, BoundingBox2D bounds, int32_t depth=*) - cdef tuple _split(self, list items, BoundingBox2D bounds) - cdef void _get_edges(self, list items, int32_t axis, int32_t *num_edges, edge **edges_ptr) - cdef void _free_edges(self, edge **edges_ptr) - cdef BoundingBox2D _get_lower_bounds(self, BoundingBox2D bounds, double split, int32_t axis) - cdef BoundingBox2D _get_upper_bounds(self, BoundingBox2D bounds, double split, int32_t axis) - cdef int32_t _new_leaf(self, list ids) - cdef int32_t _new_branch(self, tuple split_solution, int32_t depth) - cdef int32_t _new_node(self) - cpdef bint is_contained(self, Point2D point) - cdef bint _is_contained(self, Point2D point) - cdef bint _is_contained_node(self, int32_t id, Point2D point) - cdef bint _is_contained_branch(self, int32_t id, Point2D point) - cdef bint _is_contained_leaf(self, int32_t id, Point2D point) - cpdef list items_containing(self, Point2D point) - cdef list _items_containing(self, Point2D point) - cdef list _items_containing_node(self, int32_t id, Point2D point) - cdef list _items_containing_branch(self, int32_t id, Point2D point) - cdef list _items_containing_leaf(self, int32_t id, Point2D point) - cdef void _reset(self) - cdef double _read_double(self, object file) - cdef int32_t _read_int32(self, object file) cdef class KDTree2D(KDTree2DCore): cdef bint _is_contained_leaf(self, int32_t id, Point2D point) - cpdef bint _is_contained_items(self, list items, Point2D point) - cdef list _items_containing_leaf(self, int32_t id, Point2D point) - cpdef list _items_containing_items(self, list items, Point2D point) \ No newline at end of file diff --git a/raysect/core/math/spatial/kdtree2d.pyx b/raysect/core/math/spatial/kdtree2d.pyx index 879f4a4d..526eee0d 100644 --- a/raysect/core/math/spatial/kdtree2d.pyx +++ b/raysect/core/math/spatial/kdtree2d.pyx @@ -40,16 +40,20 @@ from libc.stdint cimport int32_t from libc.math cimport log, ceil cimport cython -# this number of nodes will be pre-allocated when the kd-tree is initially created -DEF INITIAL_NODE_COUNT = 128 -# friendly name for first node -DEF ROOT_NODE = 0 +# constants +cdef enum: -# node types -DEF LEAF = -1 # leaf node -DEF X_AXIS = 0 # branch, x-axis split -DEF Y_AXIS = 1 # branch, y-axis split + # this number of nodes will be pre-allocated when the kd-tree is initially created + INITIAL_NODE_COUNT = 128 + + # friendly name for first node + ROOT_NODE = 0 + + # node types + LEAF = -1 # leaf node + X_AXIS = 0 # branch, x-axis split + Y_AXIS = 1 # branch, y-axis split cdef class Item2D: @@ -69,7 +73,6 @@ cdef class Item2D: """ def __init__(self, int32_t id, BoundingBox2D box): - self.id = id self.box = box diff --git a/raysect/core/math/spatial/kdtree3d.pxd b/raysect/core/math/spatial/kdtree3d.pxd index fd7d2745..cd80f059 100644 --- a/raysect/core/math/spatial/kdtree3d.pxd +++ b/raysect/core/math/spatial/kdtree3d.pxd @@ -69,66 +69,37 @@ cdef class KDTree3DCore: double _empty_bonus cdef int32_t _build(self, list items, BoundingBox3D bounds, int32_t depth=*) - cdef tuple _split(self, list items, BoundingBox3D bounds) - cdef void _get_edges(self, list items, int32_t axis, int32_t *num_edges, edge **edges_ptr) - cdef void _free_edges(self, edge **edges_ptr) - cdef BoundingBox3D _get_lower_bounds(self, BoundingBox3D bounds, double split, int32_t axis) - cdef BoundingBox3D _get_upper_bounds(self, BoundingBox3D bounds, double split, int32_t axis) - cdef int32_t _new_leaf(self, list ids) - cdef int32_t _new_branch(self, tuple split_solution, int32_t depth) - cdef int32_t _new_node(self) - cpdef bint is_contained(self, Point3D point) - cdef bint _is_contained(self, Point3D point) - cdef bint _is_contained_node(self, int32_t id, Point3D point) - cdef bint _is_contained_branch(self, int32_t id, Point3D point) - cdef bint _is_contained_leaf(self, int32_t id, Point3D point) - cpdef bint trace(self, Ray ray) - cdef bint _trace(self, Ray ray) - cdef bint _trace_node(self, int32_t id, Ray ray, double min_range, double max_range) - cdef bint _trace_branch(self, int32_t id, Ray ray, double min_range, double max_range) - cdef bint _trace_leaf(self, int32_t id, Ray ray, double max_range) - cpdef list items_containing(self, Point3D point) - cdef list _items_containing(self, Point3D point) - cdef list _items_containing_node(self, int32_t id, Point3D point) - cdef list _items_containing_branch(self, int32_t id, Point3D point) - cdef list _items_containing_leaf(self, int32_t id, Point3D point) - cdef void _reset(self) - cdef double _read_double(self, object file) - cdef int32_t _read_int32(self, object file) cdef class KDTree3D(KDTree3DCore): cdef bint _trace_leaf(self, int32_t id, Ray ray, double max_range) - cpdef bint _trace_items(self, list items, Ray ray, double max_range) - cdef list _items_containing_leaf(self, int32_t id, Point3D point) - cpdef list _items_containing_items(self, list items, Point3D point) \ No newline at end of file diff --git a/raysect/core/math/spatial/kdtree3d.pyx b/raysect/core/math/spatial/kdtree3d.pyx index f8a95850..c1a8a195 100644 --- a/raysect/core/math/spatial/kdtree3d.pyx +++ b/raysect/core/math/spatial/kdtree3d.pyx @@ -40,17 +40,21 @@ from libc.stdint cimport int32_t from libc.math cimport log, ceil cimport cython -# this number of nodes will be pre-allocated when the kd-tree is initially created -DEF INITIAL_NODE_COUNT = 128 -# friendly name for first node -DEF ROOT_NODE = 0 +# constants +cdef enum: -# node types -DEF LEAF = -1 # leaf node -DEF X_AXIS = 0 # branch, x-axis split -DEF Y_AXIS = 1 # branch, y-axis split -DEF Z_AXIS = 2 # branch, z-axis split + # this number of nodes will be pre-allocated when the kd-tree is initially created + INITIAL_NODE_COUNT = 128 + + # friendly name for first node + ROOT_NODE = 0 + + # node types + LEAF = -1 # leaf node + X_AXIS = 0 # branch, x-axis split + Y_AXIS = 1 # branch, y-axis split + Z_AXIS = 2 # branch, z-axis split cdef class Item3D: diff --git a/raysect/core/math/statsarray.pxd b/raysect/core/math/statsarray.pxd index 6ea7bf7e..3352331a 100644 --- a/raysect/core/math/statsarray.pxd +++ b/raysect/core/math/statsarray.pxd @@ -40,13 +40,9 @@ cdef class StatsBin: readonly int samples cpdef object clear(self) - cpdef StatsBin copy(self) - cpdef object add_sample(self, double sample) - cpdef object combine_samples(self, double mean, double variance, int sample_count) - cpdef double error(self) @@ -62,19 +58,12 @@ cdef class StatsArray1D: int[::1] samples_mv cpdef object clear(self) - cpdef StatsArray1D copy(self) - cpdef object add_sample(self, int x, double sample) - cpdef object combine_samples(self, int x, double mean, double variance, int sample_count) - cpdef double error(self, int x) - cpdef ndarray errors(self) - cdef void _new_buffers(self) - cdef object _bounds_check(self, int x) @@ -90,19 +79,12 @@ cdef class StatsArray2D: int[:,::1] samples_mv cpdef object clear(self) - cpdef StatsArray2D copy(self) - cpdef object add_sample(self, int x, int y, double sample) - cpdef object combine_samples(self, int x, int y, double mean, double variance, int sample_count) - cpdef double error(self, int x, int y) - cpdef ndarray errors(self) - cdef void _new_buffers(self) - cdef object _bounds_check(self, int x, int y) @@ -118,17 +100,10 @@ cdef class StatsArray3D: int[:,:,::1] samples_mv cpdef object clear(self) - cpdef StatsArray3D copy(self) - cpdef object add_sample(self, int x, int y, int z, double sample) - cpdef object combine_samples(self, int x, int y, int z, double mean, double variance, int sample_count) - cpdef double error(self, int x, int y, int z) - cpdef ndarray errors(self) - cdef void _new_buffers(self) - cdef object _bounds_check(self, int x, int y, int z) diff --git a/raysect/core/math/statsarray.pyx b/raysect/core/math/statsarray.pyx index 491467af..4b4a60b2 100644 --- a/raysect/core/math/statsarray.pyx +++ b/raysect/core/math/statsarray.pyx @@ -52,13 +52,19 @@ cdef class StatsBin: self.samples = 0 cpdef object clear(self): - """ Erase the current statistics stored in this StatsBin. """ + """ + Erase the current statistics stored in this StatsBin. + """ + self.mean = 0.0 self.variance = 0.0 self.samples = 0 cpdef StatsBin copy(self): - """ Instantiate a new StatsBin object with the same statistical results. """ + """ + Instantiate a new StatsBin object with the same statistical results. + """ + obj = StatsBin() obj.mean = self.mean obj.variance = self.variance @@ -71,6 +77,7 @@ cdef class StatsBin: :param float sample: The sample value to be added. """ + _add_sample(sample, &self.mean, &self.variance, &self.samples) cpdef object combine_samples(self, double mean, double variance, int sample_count): @@ -115,7 +122,10 @@ cdef class StatsBin: self.samples = nt cpdef double error(self): - """ Compute the standard error of this sample distribution. """ + """ + Compute the standard error of this sample distribution. + """ + return _std_error(self.variance, self.samples) @@ -154,17 +164,25 @@ cdef class StatsArray1D: @property def shape(self): - """ The numpy style array shape of the underlying StatsArray. """ + """ + The numpy style array shape of the underlying StatsArray. + """ return (self.length, ) cpdef object clear(self): - """ Erase the current statistics stored in this StatsArray. """ + """ + Erase the current statistics stored in this StatsArray. + """ + self._new_buffers() @cython.initializedcheck(False) cpdef StatsArray1D copy(self): - """ Instantiate a new StatsArray1D object with the same statistical results. """ + """ + Instantiate a new StatsArray1D object with the same statistical results. + """ + obj = StatsArray1D(self.length) obj.mean_mv[:] = self.mean_mv[:] obj.variance_mv[:] = self.variance_mv[:] @@ -335,16 +353,25 @@ cdef class StatsArray2D: @property def shape(self): - """ The numpy style array shape of the underlying StatsArray. """ + """ + The numpy style array shape of the underlying StatsArray. + """ + return self.nx, self.ny cpdef object clear(self): - """ Erase the current statistics stored in this StatsArray. """ + """ + Erase the current statistics stored in this StatsArray. + """ + self._new_buffers() @cython.initializedcheck(False) cpdef StatsArray2D copy(self): - """ Instantiate a new StatsArray2D object with the same statistical results. """ + """ + Instantiate a new StatsArray2D object with the same statistical results. + """ + obj = StatsArray2D(self.nx, self.ny) obj.mean_mv[:] = self.mean_mv[:] obj.variance_mv[:] = self.variance_mv[:] @@ -531,18 +558,27 @@ cdef class StatsArray3D: @property def shape(self): - """ The numpy style array shape of the underlying StatsArray. """ + """ + The numpy style array shape of the underlying StatsArray. + """ + return self.nx, self.ny, self.nz cpdef object clear(self): - """ Erase the current statistics stored in this StatsArray. """ + """ + Erase the current statistics stored in this StatsArray. + """ + self._new_buffers() @cython.boundscheck(False) @cython.wraparound(False) @cython.initializedcheck(False) cpdef StatsArray3D copy(self): - """ Instantiate a new StatsArray3D object with the same statistical results. """ + """ + Instantiate a new StatsArray3D object with the same statistical results. + """ + obj = StatsArray3D(self.nx, self.ny, self.nz) obj.mean_mv[:] = self.mean_mv[:] obj.variance_mv[:] = self.variance_mv[:] diff --git a/raysect/core/math/transform.pyx b/raysect/core/math/transform.pyx index bf431ea0..41c1cd0b 100644 --- a/raysect/core/math/transform.pyx +++ b/raysect/core/math/transform.pyx @@ -29,18 +29,14 @@ # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. -from libc.math cimport sin, cos, sqrt, asin, atan2, M_PI as pi +from libc.math cimport sin, cos, sqrt, asin, atan2 from raysect.core.math.affinematrix cimport new_affinematrix3d from raysect.core.math.point cimport new_point3d cimport cython -DEF RAD2DEG = 57.29577951308232000 # 180 / pi -DEF DEG2RAD = 0.017453292519943295 # pi / 180 - - -FORWARD = 'forward' -UP = 'up' +cdef const double RAD2DEG = 57.29577951308232000 # 180 / pi +cdef const double DEG2RAD = 0.017453292519943295 # pi / 180 cpdef AffineMatrix3D translate(double x, double y, double z): @@ -108,7 +104,7 @@ cpdef AffineMatrix3D rotate_x(double angle): cdef double r - r = pi * angle / 180.0 + r = DEG2RAD * angle return new_affinematrix3d(1, 0, 0, 0, 0, cos(r), -sin(r), 0, 0, sin(r), cos(r), 0, @@ -136,7 +132,7 @@ cpdef AffineMatrix3D rotate_y(double angle): cdef double r - r = pi * angle / 180.0 + r = DEG2RAD * angle return new_affinematrix3d(cos(r), 0, sin(r), 0, 0, 1, 0, 0, -sin(r), 0, cos(r), 0, @@ -164,7 +160,7 @@ cpdef AffineMatrix3D rotate_z(double angle): cdef double r - r = pi * angle / 180.0 + r = DEG2RAD * angle return new_affinematrix3d(cos(r), -sin(r), 0, 0, sin(r), cos(r), 0, 0, 0, 0, 1, 0, @@ -195,7 +191,7 @@ cpdef AffineMatrix3D rotate_vector(double angle, Vector3D v): cdef double r, s, c, ci vn = v.normalise() - r = pi * angle / 180.0 + r = DEG2RAD * angle s = sin(r) c = cos(r) ci = 1.0 - c @@ -292,16 +288,16 @@ cpdef AffineMatrix3D rotate_basis(Vector3D forward, Vector3D up): cpdef tuple to_cylindrical(Point3D point): """ - Convert the given 3D point in cartesian space to cylindrical coordinates. - + Convert the given 3D point in cartesian space to cylindrical coordinates. + :param Point3D point: The 3D point to be transformed into cylindrical coordinates. :rtype: tuple :return: A tuple of r, z, phi coordinates. - + .. code-block:: pycon - + >>> from raysect.core.math import to_cylindrical, Point3D - + >>> point = Point3D(1, 1, 1) >>> to_cylindrical(point) (1.4142135623730951, 1.0, 45.0) @@ -318,15 +314,15 @@ cpdef tuple to_cylindrical(Point3D point): cpdef Point3D from_cylindrical(double r, double z, double phi): """ Convert a 3D point in cylindrical coordinates to a point in cartesian coordinates. - + :param float r: The radial coordinate. :param float z: The z-axis height coordinate. :param float phi: The azimuthal coordinate in degrees. :rtype: Point3D :return: A Point3D in cartesian space. - + .. code-block:: pycon - + >>> from raysect.core.math import from_cylindrical >>> from_cylindrical(1, 0, 45) @@ -348,21 +344,21 @@ cpdef Point3D from_cylindrical(double r, double z, double phi): cpdef (double, double, double) extract_rotation(AffineMatrix3D m, bint z_up=False): """ Extracts the rotation component of the affine matrix. - + The yaw, pitch and roll can be extracted for two common coordinate conventions by specifying the z_axis orientation: - - forward: +ve z is forward, +ve y is up, +ve x is left + + forward: +ve z is forward, +ve y is up, +ve x is left up: +ve z is up, +ve y is left, +ve x is forward The Raysect default is z axis forward. This function can be switched - to z axis up by setting the z_up parameter to True. - + to z axis up by setting the z_up parameter to True. + The matrix must consist of only rotation and translation operations. - + :param AffineMatrix3D m: An affine matrix. - :param bint z_up: Is the z-axis pointed upwards (default=False). - :return: A tuple containing (yaw, pitch, roll). + :param bint z_up: Is the z-axis pointed upwards (default=False). + :return: A tuple containing (yaw, pitch, roll). """ cdef double yaw, pitch, roll @@ -385,11 +381,11 @@ cpdef (double, double, double) extract_rotation(AffineMatrix3D m, bint z_up=Fals cpdef (double, double, double) extract_translation(AffineMatrix3D m): """ Extracts the translation component of the affine matrix. - + The matrix must consist of only rotation and translation operations. - :param AffineMatrix3D m: An affine matrix. - :return: tuple containing the x, y and z components of the translation. + :param AffineMatrix3D m: An affine matrix. + :return: tuple containing the x, y and z components of the translation. """ return m.get_element(0, 3), m.get_element(1, 3), m.get_element(2, 3) diff --git a/raysect/core/math/units.pyx b/raysect/core/math/units.pyx index 07b8375d..8488725d 100644 --- a/raysect/core/math/units.pyx +++ b/raysect/core/math/units.pyx @@ -31,6 +31,7 @@ from libc.math cimport M_1_PI + cpdef double km(double v): """ Converts kilometers to meters. @@ -40,6 +41,7 @@ cpdef double km(double v): """ return v * 1e3 + cpdef double cm(double v): """ Converts centimeters to meters. @@ -49,6 +51,7 @@ cpdef double cm(double v): """ return v * 1e-2 + cpdef double mm(double v): """ Converts millimeters to meters. @@ -58,6 +61,7 @@ cpdef double mm(double v): """ return v * 1e-3 + cpdef double um(double v): """ Converts micrometers to meters. @@ -67,6 +71,7 @@ cpdef double um(double v): """ return v * 1e-6 + cpdef double nm(double v): """ Converts nanometers to meters. @@ -76,6 +81,7 @@ cpdef double nm(double v): """ return v * 1e-9 + cpdef double mile(double v): """ Converts miles to meters. @@ -85,6 +91,7 @@ cpdef double mile(double v): """ return v * 1609.34 + cpdef double yard(double v): """ Converts yards to meters. @@ -94,6 +101,7 @@ cpdef double yard(double v): """ return v * 0.9144 + cpdef double foot(double v): """ Converts feet to meters. @@ -103,6 +111,7 @@ cpdef double foot(double v): """ return v * 0.3048 + cpdef double inch(double v): """ Converts inches to meters. @@ -112,6 +121,7 @@ cpdef double inch(double v): """ return v * 0.0254 + cpdef double mil(double v): """ Converts mils (thousandths of an inch) to meters. @@ -121,6 +131,7 @@ cpdef double mil(double v): """ return v * 2.54e-5 + cpdef radian(double v): """ Converts radians to degrees. diff --git a/raysect/core/math/vector.pxd b/raysect/core/math/vector.pxd index 07a0f8a5..c0843002 100644 --- a/raysect/core/math/vector.pxd +++ b/raysect/core/math/vector.pxd @@ -32,30 +32,20 @@ from raysect.core.math._vec3 cimport _Vec3 from raysect.core.math.affinematrix cimport AffineMatrix3D + cdef class Vector3D(_Vec3): cpdef Vector3D cross(self, _Vec3 v) - cpdef Vector3D normalise(self) - cpdef Vector3D transform(self, AffineMatrix3D m) - cdef Vector3D neg(self) - cdef Vector3D add(self, _Vec3 v) - cdef Vector3D sub(self, _Vec3 v) - cdef Vector3D mul(self, double m) - cdef Vector3D div(self, double m) - cpdef Vector3D copy(self) - cpdef Vector3D orthogonal(self) - cpdef Vector3D lerp(self, Vector3D b, double t) - cpdef Vector3D slerp(self, Vector3D b, double t) @@ -80,33 +70,19 @@ cdef class Vector2D: cdef public double x, y cpdef double dot(self, Vector2D v) - cdef double get_length(self) nogil - cdef object set_length(self, double v) - cdef double get_index(self, int index) nogil - cdef void set_index(self, int index, double value) nogil - cpdef double cross(self, Vector2D v) - cpdef Vector2D normalise(self) - # cpdef Vector2D transform(self, AffineMatrix2D m): - cdef Vector2D neg(self) - cdef Vector2D add(self, Vector2D v) - cdef Vector2D sub(self, Vector2D v) - cdef Vector2D mul(self, double m) - cdef Vector3D div(self, double d) - cpdef Vector2D copy(self) - cpdef Vector2D orthogonal(self) diff --git a/raysect/core/math/vector.pyx b/raysect/core/math/vector.pyx index 054479b0..54e294a2 100644 --- a/raysect/core/math/vector.pyx +++ b/raysect/core/math/vector.pyx @@ -33,7 +33,8 @@ import numbers cimport cython from libc.math cimport sqrt, fabs, NAN, acos, cos, sin -DEF EPSILON = 1e-12 + +cdef const double EPSILON = 1e-12 cdef class Vector3D(_Vec3): @@ -68,12 +69,16 @@ cdef class Vector3D(_Vec3): self.z = z def __repr__(self): - """Returns a string representation of the Vector3D object.""" + """ + Returns a string representation of the Vector3D object. + """ return "Vector3D(" + str(self.x) + ", " + str(self.y) + ", " + str(self.z) + ")" def __richcmp__(self, object other, int op): - """Provides basic vector comparison operations.""" + """ + Provides basic vector comparison operations. + """ cdef Vector3D v @@ -85,11 +90,12 @@ cdef class Vector3D(_Vec3): return self.x == v.x and self.y == v.y and self.z == v.z elif op == 3: # __ne__() return self.x != v.x or self.y != v.y or self.z != v.z - else: - return NotImplemented + + return NotImplemented def __getitem__(self, int i): - """Returns the vector coordinates by index ([0,1,2] -> [x,y,z]). + """ + Returns the vector coordinates by index ([0,1,2] -> [x,y,z]). >>> a = Vector3D(1, 0, 0) >>> a[0] @@ -102,11 +108,12 @@ cdef class Vector3D(_Vec3): return self.y elif i == 2: return self.z - else: - raise IndexError("Index out of range [0, 2].") + + raise IndexError("Index out of range [0, 2].") def __setitem__(self, int i, double value): - """Sets the vector coordinates by index ([0,1,2] -> [x,y,z]). + """ + Sets the vector coordinates by index ([0,1,2] -> [x,y,z]). >>> a = Vector3D(1, 0, 0) >>> a[1] = 2 @@ -124,122 +131,140 @@ cdef class Vector3D(_Vec3): raise IndexError("Index out of range [0, 2].") def __iter__(self): - """Iterates over the vector coordinates (x, y, z) + """ + Iterates over the vector coordinates (x, y, z) >>> a = Vector3D(0, 1, 2) >>> x, y, z = a >>> x, y, z (0.0, 1.0, 2.0) """ + yield self.x yield self.y yield self.z def __neg__(self): - """Returns a vector with the reverse orientation (negation operator). + """ + Returns a vector with the reverse orientation (negation operator). >>> a = Vector3D(1, 0, 0) >>> -a Vector3D(-1.0, -0.0, -0.0) """ - return new_vector3d(-self.x, - -self.y, - -self.z) + return new_vector3d(-self.x, -self.y, -self.z) - def __add__(object x, object y): - """Addition operator. + def __add__(self, object y): + """ + Addition operator. >>> Vector3D(1, 0, 0) + Vector3D(0, 1, 0) Vector3D(1.0, 1.0, 0.0) """ - cdef _Vec3 vx, vy + cdef _Vec3 vy - if isinstance(x, _Vec3) and isinstance(y, _Vec3): + if isinstance(y, _Vec3): + vy = <_Vec3> y + return new_vector3d(self.x + vy.x, self.y + vy.y, self.z + vy.z) - vx = <_Vec3>x - vy = <_Vec3>y + return NotImplemented - return new_vector3d(vx.x + vy.x, - vx.y + vy.y, - vx.z + vy.z) + def __radd__(self, object x): + """ + Reverse addition operator. - else: + >>> Vector3D(1, 0, 0) + Vector3D(0, 1, 0) + Vector3D(1.0, 1.0, 0.0) + """ - return NotImplemented + return self.__add__(x) - def __sub__(object x, object y): - """Subtraction operator. + def __sub__(self, object y): + """ + Subtraction operator. >>> Vector3D(1, 0, 0) - Vector3D(0, 1, 0) Vector3D(1.0, -1.0, 0.0) """ - cdef _Vec3 vx, vy + cdef _Vec3 vy - if isinstance(x, _Vec3) and isinstance(y, _Vec3): + if isinstance(y, _Vec3): + vy = <_Vec3> y + return new_vector3d(self.x - vy.x, self.y - vy.y, self.z - vy.z) - vx = <_Vec3>x - vy = <_Vec3>y + return NotImplemented - return new_vector3d(vx.x - vy.x, - vx.y - vy.y, - vx.z - vy.z) + def __rsub__(self, object x): + """ + Reverse subtraction operator. - else: + >>> Vector3D(1, 0, 0) - Vector3D(0, 1, 0) + Vector3D(1.0, -1.0, 0.0) + """ - return NotImplemented + cdef _Vec3 vx - def __mul__(object x, object y): - """Multiplication operator. + if isinstance(x, _Vec3): + vx = <_Vec3> x + return new_vector3d(vx.x - self.x, vx.y - self.y, vx.z - self.z) + + return NotImplemented + + def __mul__(self, object y): + """ + Multiplication operator. 3D vectors can be multiplied with both scalars and transformation matrices. >>> from raysect.core import Vector3D, rotate_x - >>> 2 * Vector3D(1, 2, 3) + >>> Vector3D(1, 2, 3) * 2 Vector3D(2.0, 4.0, 6.0) - >>> rotate_x(90) * Vector3D(0, 0, 1) - Vector3D(0.0, -1.0, 0.0) """ cdef double s - cdef Vector3D v cdef AffineMatrix3D m - if isinstance(x, numbers.Real) and isinstance(y, Vector3D): - - s = x - v = y + if isinstance(y, numbers.Real): + s = y + return new_vector3d(s * self.x, s * self.y, s * self.z) - return new_vector3d(s * v.x, - s * v.y, - s * v.z) + return NotImplemented - elif isinstance(x, Vector3D) and isinstance(y, numbers.Real): - - s = y - v = x + def __rmul__(self, object x): + """ + Reverse multiplication operator. - return new_vector3d(s * v.x, - s * v.y, - s * v.z) + 3D vectors can be multiplied with both scalars and transformation matrices. - elif isinstance(x, AffineMatrix3D) and isinstance(y, Vector3D): + >>> from raysect.core import Vector3D, rotate_x + >>> 2 * Vector3D(1, 2, 3) + Vector3D(2.0, 4.0, 6.0) + >>> rotate_x(90) * Vector3D(0, 0, 1) + Vector3D(0.0, -1.0, 0.0) + """ - m = x - v = y + cdef double s + cdef AffineMatrix3D m - return new_vector3d(m.m[0][0] * v.x + m.m[0][1] * v.y + m.m[0][2] * v.z, - m.m[1][0] * v.x + m.m[1][1] * v.y + m.m[1][2] * v.z, - m.m[2][0] * v.x + m.m[2][1] * v.y + m.m[2][2] * v.z) + if isinstance(x, numbers.Real): + s = x + return new_vector3d(s * self.x, s * self.y, s * self.z) - else: + elif isinstance(x, AffineMatrix3D): + m = x + return new_vector3d( + m.m[0][0] * self.x + m.m[0][1] * self.y + m.m[0][2] * self.z, + m.m[1][0] * self.x + m.m[1][1] * self.y + m.m[1][2] * self.z, + m.m[2][0] * self.x + m.m[2][1] * self.y + m.m[2][2] * self.z + ) - return NotImplemented + return NotImplemented @cython.cdivision(True) - def __truediv__(object x, object y): + def __truediv__(self, object y): """Division operator. >>> Vector3D(1, 1, 1) / 2 @@ -247,27 +272,19 @@ cdef class Vector3D(_Vec3): """ cdef double d - cdef Vector3D v - if isinstance(x, Vector3D) and isinstance(y, numbers.Real): + if isinstance(y, numbers.Real): d = y # prevent divide my zero if d == 0.0: - raise ZeroDivisionError("Cannot divide a vector by a zero scalar.") - v = x d = 1.0 / d + return new_vector3d(d * self.x, d * self.y, d * self.z) - return new_vector3d(d * v.x, - d * v.y, - d * v.z) - - else: - - raise TypeError("Unsupported operand type. Expects a real number.") + raise TypeError("Unsupported operand type. Expects a real number.") cpdef Vector3D cross(self, _Vec3 v): """ @@ -286,9 +303,11 @@ cdef class Vector3D(_Vec3): Vector3D(0.0, 0.0, 1.0) """ - return new_vector3d(self.y * v.z - v.y * self.z, - self.z * v.x - v.z * self.x, - self.x * v.y - v.x * self.y) + return new_vector3d( + self.y * v.z - v.y * self.z, + self.z * v.x - v.z * self.x, + self.x * v.y - v.x * self.y + ) @cython.cdivision(True) cpdef Vector3D normalise(self): @@ -311,15 +330,11 @@ cdef class Vector3D(_Vec3): # if current length is zero, problem is ill defined t = self.x * self.x + self.y * self.y + self.z * self.z if t == 0.0: - raise ZeroDivisionError("A zero length vector can not be normalised as the direction of a zero length vector is undefined.") # normalise and rescale vector t = 1.0 / sqrt(t) - - return new_vector3d(self.x * t, - self.y * t, - self.z * t) + return new_vector3d(self.x * t, self.y * t, self.z * t) cpdef Vector3D transform(self, AffineMatrix3D m): """ @@ -347,9 +362,11 @@ cdef class Vector3D(_Vec3): Vector3D(0.0, -1.0, 6.123233995736766e-17) """ - return new_vector3d(m.m[0][0] * self.x + m.m[0][1] * self.y + m.m[0][2] * self.z, - m.m[1][0] * self.x + m.m[1][1] * self.y + m.m[1][2] * self.z, - m.m[2][0] * self.x + m.m[2][1] * self.y + m.m[2][2] * self.z) + return new_vector3d( + m.m[0][0] * self.x + m.m[0][1] * self.y + m.m[0][2] * self.z, + m.m[1][0] * self.x + m.m[1][1] * self.y + m.m[1][2] * self.z, + m.m[2][0] * self.x + m.m[2][1] * self.y + m.m[2][2] * self.z + ) cdef Vector3D neg(self): """ @@ -359,9 +376,7 @@ cdef class Vector3D(_Vec3): to the equivalent python operator. """ - return new_vector3d(-self.x, - -self.y, - -self.z) + return new_vector3d(-self.x, -self.y, -self.z) cdef Vector3D add(self, _Vec3 v): """ @@ -371,9 +386,7 @@ cdef class Vector3D(_Vec3): to the equivalent python operator. """ - return new_vector3d(self.x + v.x, - self.y + v.y, - self.z + v.z) + return new_vector3d(self.x + v.x, self.y + v.y, self.z + v.z) cdef Vector3D sub(self, _Vec3 v): """ @@ -383,9 +396,7 @@ cdef class Vector3D(_Vec3): to the equivalent python operator. """ - return new_vector3d(self.x - v.x, - self.y - v.y, - self.z - v.z) + return new_vector3d(self.x - v.x, self.y - v.y, self.z - v.z) cdef Vector3D mul(self, double m): """ @@ -395,9 +406,7 @@ cdef class Vector3D(_Vec3): to the equivalent python operator. """ - return new_vector3d(self.x * m, - self.y * m, - self.z * m) + return new_vector3d(self.x * m, self.y * m, self.z * m) @cython.cdivision(True) cdef Vector3D div(self, double d): @@ -409,14 +418,10 @@ cdef class Vector3D(_Vec3): """ if d == 0.0: - raise ZeroDivisionError("Cannot divide a vector by a zero scalar.") d = 1.0 / d - - return new_vector3d(self.x * d, - self.y * d, - self.z * d) + return new_vector3d(self.x * d, self.y * d, self.z * d) cpdef Vector3D copy(self): """ @@ -431,9 +436,7 @@ cdef class Vector3D(_Vec3): Vector3D(1.0, 1.0, 1.0) """ - return new_vector3d(self.x, - self.y, - self.z) + return new_vector3d(self.x, self.y, self.z) # todo: this is common code with normal, move into math.cython and call cpdef Vector3D orthogonal(self): @@ -493,7 +496,6 @@ cdef class Vector3D(_Vec3): raise ValueError("Vector lerp parameter t must be in range (0, 1).") t_minus = 1 - t - return new_vector3d(self.x * t_minus + b.x * t, self.y * t_minus + b.y * t, self.z * t_minus + b.z * t) cpdef Vector3D slerp(self, Vector3D b, double t): @@ -631,12 +633,16 @@ cdef class Vector2D: self.y = y def __repr__(self): - """Returns a string representation of the Vector2D object.""" + """ + Returns a string representation of the Vector2D object. + """ return "Vector2D(" + str(self.x) + ", " + str(self.y) + ")" def __richcmp__(self, object other, int op): - """Provides basic vector comparison operations.""" + """ + Provides basic vector comparison operations. + """ cdef Vector2D v @@ -648,11 +654,12 @@ cdef class Vector2D: return self.x == v.x and self.y == v.y elif op == 3: # __ne__() return self.x != v.x or self.y != v.y - else: - return NotImplemented + + return NotImplemented def __getitem__(self, int i): - """Returns the vector coordinates by index ([0,1] -> [x,y]). + """ + Returns the vector coordinates by index ([0,1] -> [x,y]). >>> a = Vector2D(1, 0) >>> a[0] @@ -663,11 +670,12 @@ cdef class Vector2D: return self.x elif i == 1: return self.y - else: - raise IndexError("Index out of range [0, 1].") + + raise IndexError("Index out of range [0, 1].") def __setitem__(self, int i, double value): - """Sets the vector coordinates by index ([0,1] -> [x,y]). + """ + Sets the vector coordinates by index ([0,1] -> [x,y]). >>> a = Vector2D(1, 0) >>> a[1] = 2 @@ -683,7 +691,8 @@ cdef class Vector2D: raise IndexError("Index out of range [0, 1].") def __iter__(self): - """Iterates over the vector coordinates (x, y) + """ + Iterates over the vector coordinates (x, y) >>> a = Vector2D(1, 0) >>> x, y = a @@ -691,11 +700,13 @@ cdef class Vector2D: (1.0, 0.0) """ + yield self.x yield self.y def __neg__(self): - """Returns a vector with the reverse orientation (negation operator). + """ + Returns a vector with the reverse orientation (negation operator). >>> a = Vector2D(1, 0) >>> -a @@ -705,113 +716,93 @@ cdef class Vector2D: return new_vector2d(-self.x, -self.y) - def __add__(object x, object y): - """Addition operator. + def __add__(self, object y): + """ + Addition operator. >>> Vector2D(1, 0) + Vector2D(0, 1) Vector2D(1.0, 1.0) """ - cdef Vector2D vx, vy - - if isinstance(x, Vector2D) and isinstance(y, Vector2D): + cdef Vector2D vy - vx = x + if isinstance(y, Vector2D): vy = y + return new_vector2d(self.x + vy.x, self.y + vy.y) - return new_vector2d(vx.x + vy.x, vx.y + vy.y) - - else: - - return NotImplemented + return NotImplemented - def __sub__(object x, object y): - """Subtraction operator. + def __sub__(self, object y): + """ + Subtraction operator. >>> Vector2D(1, 0) - Vector2D(0, 1) Vector2D(1.0, -1.0) """ - cdef Vector2D vx, vy - - if isinstance(x, Vector2D) and isinstance(y, Vector2D): + cdef Vector2D vy - vx = x + if isinstance(y, Vector2D): vy = y + return new_vector2d(self.x - vy.x, self.y - vy.y) - return new_vector2d(vx.x - vy.x, vx.y - vy.y) - - else: + return NotImplemented - return NotImplemented - - # TODO - add 2D affine transformations - def __mul__(object x, object y): - """Multiplication operator. + def __mul__(self, object y): + """ + Multiplication operator. - >>> 2 * Vector3D(1, 2) + >>> Vector3D(1, 2) * 2 Vector2D(2.0, 4.0) """ cdef double s - cdef Vector2D v - # cdef AffineMatrix2D m - - if isinstance(x, numbers.Real) and isinstance(y, Vector2D): - - s = x - v = y - return new_vector2d(s * v.x, s * v.y,) + if isinstance(y, numbers.Real): + s = y + return new_vector2d(s * self.x, s * self.y) - elif isinstance(x, Vector2D) and isinstance(y, numbers.Real): + raise TypeError("Unsupported operand type. Expects a real number.") - s = y - v = x + def __rmul__(self, object x): + """ + Reverse multiplication operator. - return new_vector2d(s * v.x, s * v.y) + >>> 2 * Vector3D(1, 2) + Vector2D(2.0, 4.0) + """ - # elif isinstance(x, AffineMatrix3D) and isinstance(y, Vector3D): - # - # m = x - # v = y - # - # return new_vector3d(m.m[0][0] * v.x + m.m[0][1] * v.y + m.m[0][2] * v.z, - # m.m[1][0] * v.x + m.m[1][1] * v.y + m.m[1][2] * v.z, - # m.m[2][0] * v.x + m.m[2][1] * v.y + m.m[2][2] * v.z) + cdef double s - else: + if isinstance(x, numbers.Real): + s = x + return new_vector2d(s * self.x, s * self.y) - return NotImplemented + raise TypeError("Unsupported operand type. Expects a real number.") @cython.cdivision(True) - def __truediv__(object x, object y): - """Division operator. + def __truediv__(self, object y): + """ + Division operator. >>> Vector2D(1, 1) / 2 Vector2D(0.5, 0.5) """ cdef double d - cdef Vector2D v - if isinstance(x, Vector2D) and isinstance(y, numbers.Real): + if isinstance(y, numbers.Real): d = y # prevent divide my zero if d == 0.0: - raise ZeroDivisionError("Cannot divide a vector by a zero scalar.") - v = x d = 1.0 / d + return new_vector2d(d * self.x, d * self.y) - return new_vector2d(d * v.x, d * v.y) - - else: - - raise TypeError("Unsupported operand type. Expects a real number.") + raise TypeError("Unsupported operand type. Expects a real number.") @property def length(self): @@ -827,6 +818,7 @@ cdef class Vector2D: 1.4142135623730951 """ + return self.get_length() @length.setter @@ -880,7 +872,6 @@ cdef class Vector2D: # normalise and rescale vector t = v / sqrt(t) - self.x = self.x * t self.y = self.y * t @@ -897,8 +888,8 @@ cdef class Vector2D: return self.x elif index == 1: return self.y - else: - return NAN + + return NAN cdef void set_index(self, int index, double value) nogil: """ @@ -959,7 +950,6 @@ cdef class Vector2D: # if current length is zero, problem is ill defined t = self.x * self.x + self.y * self.y if t == 0.0: - raise ZeroDivisionError("A zero length vector can not be normalised as the direction of a zero length vector is undefined.") # normalise and rescale vector @@ -1036,11 +1026,9 @@ cdef class Vector2D: """ if d == 0.0: - raise ZeroDivisionError("Cannot divide a vector by a zero scalar.") d = 1.0 / d - return new_vector2d(self.x * d, self.y * d) cpdef Vector2D copy(self): @@ -1073,11 +1061,7 @@ cdef class Vector2D: """ - cdef: - Vector2D n - - n = self.normalise() - + cdef Vector2D n = self.normalise() return new_vector2d(-n.y, n.x) diff --git a/raysect/core/ray.pxd b/raysect/core/ray.pxd index db89e3cf..4e02656b 100644 --- a/raysect/core/ray.pxd +++ b/raysect/core/ray.pxd @@ -39,7 +39,6 @@ cdef class Ray: cdef public double max_distance cpdef Point3D point_on(self, double t) - cpdef Ray copy(self, Point3D origin=*, Vector3D direction=*) diff --git a/raysect/core/ray.pyx b/raysect/core/ray.pyx index 8eca9339..e0549923 100644 --- a/raysect/core/ray.pyx +++ b/raysect/core/ray.pyx @@ -29,12 +29,10 @@ # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. +from libc.math cimport INFINITY from raysect.core.math cimport new_point3d cimport cython -# cython doesn't have a built-in infinity constant, this compiles to +infinity -DEF INFINITY = 1e999 - @cython.freelist(256) cdef class Ray: diff --git a/raysect/core/scenegraph/node.pxd b/raysect/core/scenegraph/node.pxd index 2a80246b..e0350249 100644 --- a/raysect/core/scenegraph/node.pxd +++ b/raysect/core/scenegraph/node.pxd @@ -32,10 +32,9 @@ from raysect.core.scenegraph._nodebase cimport _NodeBase from raysect.core.math.affinematrix cimport AffineMatrix3D + cdef class Node(_NodeBase): cpdef AffineMatrix3D to(self, _NodeBase node) - cpdef AffineMatrix3D to_local(self) - cpdef AffineMatrix3D to_root(self) diff --git a/raysect/core/scenegraph/observer.pxd b/raysect/core/scenegraph/observer.pxd index a8fccb8e..0fdde98e 100644 --- a/raysect/core/scenegraph/observer.pxd +++ b/raysect/core/scenegraph/observer.pxd @@ -31,6 +31,7 @@ from raysect.core.scenegraph.node cimport Node + cdef class Observer(Node): cpdef observe(self) \ No newline at end of file diff --git a/raysect/core/scenegraph/primitive.pxd b/raysect/core/scenegraph/primitive.pxd index ad9d03f4..52f46b9a 100644 --- a/raysect/core/scenegraph/primitive.pxd +++ b/raysect/core/scenegraph/primitive.pxd @@ -44,19 +44,11 @@ cdef class Primitive(Node): cdef Material _material cdef Material get_material(self) - cpdef Intersection hit(self, Ray ray) - cpdef Intersection next_intersection(self) - cpdef bint contains(self, Point3D p) except -1 - cpdef BoundingBox3D bounding_box(self) - cpdef BoundingSphere3D bounding_sphere(self) - cpdef object instance(self, object parent=*, AffineMatrix3D transform=*, Material material=*, str name=*) - cpdef object notify_geometry_change(self) - - cpdef object notify_material_change(self) \ No newline at end of file + cpdef object notify_material_change(self) diff --git a/raysect/core/scenegraph/utility.pyx b/raysect/core/scenegraph/utility.pyx index 6947beb0..0870aa3d 100644 --- a/raysect/core/scenegraph/utility.pyx +++ b/raysect/core/scenegraph/utility.pyx @@ -35,6 +35,7 @@ from raysect.core.scenegraph._nodebase cimport _NodeBase from raysect.core.scenegraph.node cimport Node from raysect.core.scenegraph.signal cimport ChangeSignal + cdef class BridgeNode(Node): """ Specialised scene-graph root node that propagates geometry notifications. diff --git a/raysect/core/scenegraph/world.pxd b/raysect/core/scenegraph/world.pxd index a9e697af..57068045 100644 --- a/raysect/core/scenegraph/world.pxd +++ b/raysect/core/scenegraph/world.pxd @@ -35,18 +35,17 @@ from raysect.core.acceleration.accelerator cimport Accelerator from raysect.core.math cimport Point3D, AffineMatrix3D from raysect.core.scenegraph._nodebase cimport _NodeBase + cdef class World(_NodeBase): - cdef bint _rebuild_accelerator - cdef Accelerator _accelerator - cdef list _primitives - cdef list _observers + cdef: + bint _rebuild_accelerator + Accelerator _accelerator + list _primitives + list _observers cpdef AffineMatrix3D to(self, _NodeBase node) - cpdef Intersection hit(self, Ray ray) - cpdef list contains(self, Point3D point) - cpdef build_accelerator(self, bint force=*) diff --git a/raysect/optical/loggingray.pyx b/raysect/optical/loggingray.pyx index edc79b01..71334121 100644 --- a/raysect/optical/loggingray.pyx +++ b/raysect/optical/loggingray.pyx @@ -29,6 +29,8 @@ # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. +from libc.math cimport INFINITY + from raysect.core cimport Intersection, Point3D, Vector3D from raysect.core.math.random cimport probability @@ -39,9 +41,6 @@ from raysect.optical.material.material cimport Material cimport cython -# cython doesn't have a built-in infinity constant, this compiles to +infinity -DEF INFINITY = 1e999 - cdef class LoggingRay(Ray): diff --git a/raysect/optical/material/conductor.pxd b/raysect/optical/material/conductor.pxd index 83416241..9eb5aecb 100644 --- a/raysect/optical/material/conductor.pxd +++ b/raysect/optical/material/conductor.pxd @@ -51,11 +51,7 @@ cdef class RoughConductor(ContinuousBSDF): double _roughness cdef double _d(self, Vector3D s_half) - cdef double _g(self, Vector3D s_incoming, Vector3D s_outgoing) - cdef double _g1(self, Vector3D v) - cdef Spectrum _f(self, Spectrum spectrum, Vector3D s_outgoing, Vector3D s_normal) - cdef double _fresnel_conductor(self, double ci, double n, double k) nogil diff --git a/raysect/optical/material/debug.pxd b/raysect/optical/material/debug.pxd index 369603c8..0ec21e36 100644 --- a/raysect/optical/material/debug.pxd +++ b/raysect/optical/material/debug.pxd @@ -42,5 +42,4 @@ cdef class Light(NullVolume): cdef class PerfectReflectingSurface(Material): - pass diff --git a/raysect/optical/material/dielectric.pxd b/raysect/optical/material/dielectric.pxd index d541c7d3..969e18b5 100644 --- a/raysect/optical/material/dielectric.pxd +++ b/raysect/optical/material/dielectric.pxd @@ -32,6 +32,7 @@ from raysect.optical cimport SpectralFunction, NumericallyIntegratedSF from raysect.optical.material cimport Material + cdef class Sellmeier(NumericallyIntegratedSF): cdef: diff --git a/raysect/optical/material/emitter/unity.pxd b/raysect/optical/material/emitter/unity.pxd index 1f91e5c0..33e20a9c 100644 --- a/raysect/optical/material/emitter/unity.pxd +++ b/raysect/optical/material/emitter/unity.pxd @@ -29,7 +29,6 @@ # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. -from raysect.optical cimport World, Primitive, Ray, Spectrum, Point3D, Vector3D, Normal3D, AffineMatrix3D, Intersection from raysect.optical.material.emitter.homogeneous cimport HomogeneousVolumeEmitter from raysect.optical.material.material cimport NullVolume diff --git a/raysect/optical/material/emitter/unity.pyx b/raysect/optical/material/emitter/unity.pyx index 145fa1c9..a3fc2bea 100644 --- a/raysect/optical/material/emitter/unity.pyx +++ b/raysect/optical/material/emitter/unity.pyx @@ -29,6 +29,8 @@ # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. +from raysect.optical cimport World, Primitive, Ray, Spectrum, Point3D, Vector3D, Normal3D, AffineMatrix3D, Intersection + cimport cython diff --git a/raysect/optical/material/lambert.pyx b/raysect/optical/material/lambert.pyx index 9e2c1a67..a3ecd08b 100644 --- a/raysect/optical/material/lambert.pyx +++ b/raysect/optical/material/lambert.pyx @@ -32,7 +32,6 @@ from raysect.core.math.sampler cimport HemisphereCosineSampler from raysect.optical cimport Point3D, Vector3D, AffineMatrix3D, Primitive, World, Ray, Spectrum, SpectralFunction, ConstantSF, Intersection from raysect.optical.material cimport ContinuousBSDF -from numpy cimport ndarray cdef HemisphereCosineSampler hemisphere_sampler = HemisphereCosineSampler() diff --git a/raysect/optical/material/modifiers/roughen.pyx b/raysect/optical/material/modifiers/roughen.pyx index 024b0de3..b95da76f 100644 --- a/raysect/optical/material/modifiers/roughen.pyx +++ b/raysect/optical/material/modifiers/roughen.pyx @@ -36,7 +36,8 @@ from raysect.optical.material cimport Material # sets the maximum number of attempts to find a valid perturbed normal # it is highly unlikely (REALLY!) this number will ever be reached, it is just there for my paranoia # in the worst case 50% of the random hemisphere will always generate a valid solution... so P(fail) < 0.5^50! -DEF SAMPLE_ATTEMPTS = 50 +cdef enum: + SAMPLE_ATTEMPTS = 50 cdef HemisphereCosineSampler hemisphere_sampler = HemisphereCosineSampler() diff --git a/raysect/optical/observer/base/observer.pxd b/raysect/optical/observer/base/observer.pxd index b354b5db..39602353 100644 --- a/raysect/optical/observer/base/observer.pxd +++ b/raysect/optical/observer/base/observer.pxd @@ -55,31 +55,18 @@ cdef class _ObserverBase(Observer): public bint quiet cpdef list _slice_spectrum(self) - cpdef list _generate_templates(self, list slices) - cpdef object _render_pixel(self, tuple task, int slice_id, Ray template) - cpdef object _update_state(self, tuple packed_result, int slice_id) - cpdef list _generate_tasks(self) - cpdef list _obtain_pixel_processors(self, tuple task, int slice_id) - cpdef object _initialise_pipelines(self, double min_wavelength, double max_wavelength, int spectral_bins, list slices, bint quiet) - cpdef object _update_pipelines(self, tuple task, list results, int slice_id) - cpdef object _finalise_pipelines(self) - cpdef object _initialise_statistics(self, list tasks) - cpdef object _update_statistics(self, uint64_t sample_ray_count) - cpdef object _finalise_statistics(self) - cpdef list _obtain_rays(self, tuple task, Ray template) - cpdef double _obtain_sensitivity(self, tuple task) @@ -91,7 +78,6 @@ cdef class Observer0D(_ObserverBase): int _samples_per_task cpdef list _generate_rays(self, Ray template, int ray_count) - cpdef double _pixel_sensitivity(self) @@ -104,7 +90,6 @@ cdef class Observer1D(_ObserverBase): int _pixel_samples cpdef list _generate_rays(self, int pixel, Ray template, int ray_count) - cpdef double _pixel_sensitivity(self, int pixel) @@ -117,7 +102,6 @@ cdef class Observer2D(_ObserverBase): int _pixel_samples cpdef list _generate_rays(self, int x, int y, Ray template, int ray_count) - cpdef double _pixel_sensitivity(self, int x, int y) diff --git a/raysect/optical/observer/base/pipeline.pxd b/raysect/optical/observer/base/pipeline.pxd index 4b0437f4..0b6be0c2 100644 --- a/raysect/optical/observer/base/pipeline.pxd +++ b/raysect/optical/observer/base/pipeline.pxd @@ -35,32 +35,23 @@ from raysect.optical.observer.base.processor cimport PixelProcessor cdef class Pipeline0D: cpdef object initialise(self, double min_wavelength, double max_wavelength, int spectral_bins, list spectral_slices, bint quiet) - cpdef PixelProcessor pixel_processor(self, int slice_id) - cpdef object update(self, int slice_id, tuple packed_result, int samples) - cpdef object finalise(self) cdef class Pipeline1D: cpdef object initialise(self, int pixels, int pixel_samples, double min_wavelength, double max_wavelength, int spectral_bins, list spectral_slices, bint quiet) - cpdef PixelProcessor pixel_processor(self, int pixel, int slice_id) - cpdef object update(self, int pixel, int slice_id, tuple packed_result) - cpdef object finalise(self) cdef class Pipeline2D: cpdef object initialise(self, tuple pixels, int pixel_samples, double min_wavelength, double max_wavelength, int spectral_bins, list spectral_slices, bint quiet) - cpdef PixelProcessor pixel_processor(self, int x, int y, int slice_id) - cpdef object update(self, int x, int y, int slice_id, tuple packed_result) - cpdef object finalise(self) diff --git a/raysect/optical/observer/base/processor.pxd b/raysect/optical/observer/base/processor.pxd index c99cd6b5..10dbc820 100644 --- a/raysect/optical/observer/base/processor.pxd +++ b/raysect/optical/observer/base/processor.pxd @@ -31,9 +31,9 @@ from raysect.optical cimport Spectrum + cdef class PixelProcessor: cpdef object add_sample(self, Spectrum spectrum, double sensitivity) - cpdef tuple pack_results(self) diff --git a/raysect/optical/observer/base/processor.pyx b/raysect/optical/observer/base/processor.pyx index 9d13ff0d..ecf98dd7 100644 --- a/raysect/optical/observer/base/processor.pyx +++ b/raysect/optical/observer/base/processor.pyx @@ -31,6 +31,7 @@ cimport cython + @cython.freelist(64) cdef class PixelProcessor: """ diff --git a/raysect/optical/observer/imaging/__init__.pxd b/raysect/optical/observer/imaging/__init__.pxd index 2a62327d..946850b4 100644 --- a/raysect/optical/observer/imaging/__init__.pxd +++ b/raysect/optical/observer/imaging/__init__.pxd @@ -30,7 +30,7 @@ # POSSIBILITY OF SUCH DAMAGE. from raysect.optical.observer.imaging.ccd cimport CCDArray -from raysect.optical.observer.imaging.targetted_ccd cimport TargettedCCDArray +from raysect.optical.observer.imaging.targeted_ccd cimport TargetedCCDArray from raysect.optical.observer.imaging.orthographic cimport OrthographicCamera from raysect.optical.observer.imaging.pinhole cimport PinholeCamera from raysect.optical.observer.imaging.vector cimport VectorCamera diff --git a/raysect/optical/observer/imaging/__init__.py b/raysect/optical/observer/imaging/__init__.py index b3dfd66e..0028eedc 100644 --- a/raysect/optical/observer/imaging/__init__.py +++ b/raysect/optical/observer/imaging/__init__.py @@ -30,7 +30,7 @@ # POSSIBILITY OF SUCH DAMAGE. from .ccd import CCDArray -from .targetted_ccd import TargettedCCDArray +from .targeted_ccd import TargetedCCDArray from .orthographic import OrthographicCamera from .pinhole import PinholeCamera from .vector import VectorCamera diff --git a/raysect/optical/observer/imaging/meson.build b/raysect/optical/observer/imaging/meson.build index f20ce2f6..85a7aaa6 100644 --- a/raysect/optical/observer/imaging/meson.build +++ b/raysect/optical/observer/imaging/meson.build @@ -5,8 +5,8 @@ target_path = 'raysect/optical/observer/imaging' # source files py_files = ['__init__.py'] -pyx_files = ['ccd.pyx', 'opencv.pyx', 'orthographic.pyx', 'pinhole.pyx', 'targetted_ccd.pyx', 'vector.pyx'] -pxd_files = ['__init__.pxd', 'ccd.pxd', 'opencv.pxd', 'orthographic.pxd', 'pinhole.pxd', 'targetted_ccd.pxd', 'vector.pxd'] +pyx_files = ['ccd.pyx', 'opencv.pyx', 'orthographic.pyx', 'pinhole.pyx', 'targeted_ccd.pyx', 'vector.pyx'] +pxd_files = ['__init__.pxd', 'ccd.pxd', 'opencv.pxd', 'orthographic.pxd', 'pinhole.pxd', 'targeted_ccd.pxd', 'vector.pxd'] data_files = [] # compile cython diff --git a/raysect/optical/observer/imaging/targetted_ccd.pxd b/raysect/optical/observer/imaging/targeted_ccd.pxd similarity index 91% rename from raysect/optical/observer/imaging/targetted_ccd.pxd rename to raysect/optical/observer/imaging/targeted_ccd.pxd index f8af58b1..4e7fb777 100644 --- a/raysect/optical/observer/imaging/targetted_ccd.pxd +++ b/raysect/optical/observer/imaging/targeted_ccd.pxd @@ -29,17 +29,17 @@ # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. -from raysect.core.math.sampler cimport RectangleSampler3D, HemisphereCosineSampler, TargettedHemisphereSampler +from raysect.core.math.sampler cimport RectangleSampler3D, HemisphereCosineSampler, TargetedHemisphereSampler from raysect.optical.observer.base cimport Observer2D -cdef class TargettedCCDArray(Observer2D): +cdef class TargetedCCDArray(Observer2D): cdef: - double _width, _pixel_area, _image_delta, _image_start_x, _image_start_y, _targetted_path_prob + double _width, _pixel_area, _image_delta, _image_start_x, _image_start_y, _targeted_path_prob tuple _targets RectangleSampler3D _point_sampler HemisphereCosineSampler _cosine_sampler - TargettedHemisphereSampler _targetted_sampler + TargetedHemisphereSampler _targeted_sampler cdef object _update_image_geometry(self) diff --git a/raysect/optical/observer/imaging/targetted_ccd.pyx b/raysect/optical/observer/imaging/targeted_ccd.pyx similarity index 87% rename from raysect/optical/observer/imaging/targetted_ccd.pyx rename to raysect/optical/observer/imaging/targeted_ccd.pyx index 9a0defa6..5daa7468 100644 --- a/raysect/optical/observer/imaging/targetted_ccd.pyx +++ b/raysect/optical/observer/imaging/targeted_ccd.pyx @@ -33,34 +33,34 @@ from raysect.core.math.random cimport probability from raysect.optical.observer.sampler2d import FullFrameSampler2D from raysect.optical.observer.pipeline import RGBPipeline2D -from raysect.core.math.sampler cimport RectangleSampler3D, HemisphereCosineSampler, TargettedHemisphereSampler +from raysect.core.math.sampler cimport RectangleSampler3D, HemisphereCosineSampler, TargetedHemisphereSampler from raysect.optical cimport Primitive, BoundingSphere3D, Ray, AffineMatrix3D, Point3D, Vector3D, translate from libc.math cimport M_PI from raysect.optical.observer.base cimport Observer2D cimport cython -DEF R_2_PI = 0.15915494309189535 # 1 / (2 * pi) +cdef const double R_2_PI = 0.15915494309189535 # 1 / (2 * pi) -cdef class TargettedCCDArray(Observer2D): +cdef class TargetedCCDArray(Observer2D): """ An ideal CCD-like imaging sensor that preferentially targets a given list of primitives. - The targetted CCD is a regular array of square pixels. Each pixel samples red, green + The targeted CCD is a regular array of square pixels. Each pixel samples red, green and blue channels (behaves like a Foveon imaging sensor). The CCD sensor width is specified with the width parameter. The CCD height is calculated from the width and the number of vertical and horizontal pixels. The default width and sensor ratio approximates a 35mm camera sensor. - The targetted CCD takes a list of target primitives. Each pixel will target the + The targeted CCD takes a list of target primitives. Each pixel will target the bounding spheres that encompass each target primitive. Therefore, for best performance, the target primitives should be split up such that their surfaces are closely wrapped by the bounding sphere. The sampling algorithm fires a proportion of rays at the targets, and a portion sampled from the full hemisphere. The proportion that is fired towards the targets is controlled - with the targetted_path_prob attribute. By default this attribute is set to 0.9, i.e. + with the targeted_path_prob attribute. By default this attribute is set to 0.9, i.e. 90% of the rays are fired towards the targets. .. Warning.. @@ -68,7 +68,7 @@ cdef class TargettedCCDArray(Observer2D): targets. The user must ensure there are no sources of radiance outside of the targeted directions, otherwise they will not be sampled and the result will be biased. - :param list targets: The list of primitives for targetted sampling. + :param list targets: The list of primitives for targeted sampling. :param tuple pixels: A tuple of pixel dimensions for the camera (default=(720, 480)). :param float width: The CCD sensor x-width in metres (default=35mm). :param list pipelines: The list of pipelines that will process the spectrum measured @@ -76,7 +76,7 @@ cdef class TargettedCCDArray(Observer2D): :param kwargs: **kwargs and properties from Observer2D and _ObserverBase. """ - def __init__(self, targets, pixels=(720, 480), width=0.035, targetted_path_prob=None, parent=None, transform=None, name=None, pipelines=None): + def __init__(self, targets, pixels=(720, 480), width=0.035, targeted_path_prob=None, parent=None, transform=None, name=None, pipelines=None): # initial values to prevent undefined behaviour when setting via self.width self._width = 0.035 @@ -88,12 +88,12 @@ cdef class TargettedCCDArray(Observer2D): parent=parent, transform=transform, name=name) self._cosine_sampler = HemisphereCosineSampler() - self._targetted_sampler = None + self._targeted_sampler = None # setting width triggers calculation of image geometry calculations self.width = width self.targets = targets - self.targetted_path_prob = targetted_path_prob or 0.9 + self.targeted_path_prob = targeted_path_prob or 0.9 @property def pixels(self): @@ -158,7 +158,7 @@ cdef class TargettedCCDArray(Observer2D): self._targets = value @property - def targetted_path_prob(self): + def targeted_path_prob(self): """ The probability that an individual sample will be fired at a target instead of a sample from the whole hemisphere. @@ -169,14 +169,14 @@ cdef class TargettedCCDArray(Observer2D): :rtype: float """ - return self._targetted_path_prob + return self._targeted_path_prob - @targetted_path_prob.setter - def targetted_path_prob(self, double value): + @targeted_path_prob.setter + def targeted_path_prob(self, double value): if value < 0 or value > 1: raise ValueError("Targeted path probability must lie in the range [0, 1].") - self._targetted_path_prob = value + self._targeted_path_prob = value cdef object _update_image_geometry(self): @@ -212,8 +212,8 @@ cdef class TargettedCCDArray(Observer2D): sphere = target.bounding_sphere() spheres.append((sphere.centre.transform(self.to_local()), sphere.radius, 1.0)) - # instance targetted pixel sampler - self._targetted_sampler = TargettedHemisphereSampler(spheres) + # instance targeted pixel sampler + self._targeted_sampler = TargetedHemisphereSampler(spheres) # generate pixel transform pixel_x = self._image_start_x - self._image_delta * (ix + 0.5) @@ -230,17 +230,17 @@ cdef class TargettedCCDArray(Observer2D): # transform to local space from pixel space origin = origin.transform(pixel_to_local) - if probability(self._targetted_path_prob): - # obtain targetted vector sample - direction = self._targetted_sampler.sample(origin) + if probability(self._targeted_path_prob): + # obtain targeted vector sample + direction = self._targeted_sampler.sample(origin) else: # obtain cosine weighted hemisphere sample direction = self._cosine_sampler.sample() # calculate combined pdf and ray weight - pdf = self._targetted_path_prob * self._targetted_sampler.pdf(origin, direction) + \ - (1-self._targetted_path_prob) * self._cosine_sampler.pdf(direction) + pdf = self._targeted_path_prob * self._targeted_sampler.pdf(origin, direction) + \ + (1-self._targeted_path_prob) * self._cosine_sampler.pdf(direction) if pdf <= 0: raise ValueError('Ray direction probability is zero. The target object extends beyond the pixel horizon.') diff --git a/raysect/optical/observer/nonimaging/__init__.pxd b/raysect/optical/observer/nonimaging/__init__.pxd index 4e2f420a..f550d902 100644 --- a/raysect/optical/observer/nonimaging/__init__.pxd +++ b/raysect/optical/observer/nonimaging/__init__.pxd @@ -32,7 +32,7 @@ from raysect.optical.observer.nonimaging.sightline import SightLine from raysect.optical.observer.nonimaging.fibreoptic import FibreOptic from raysect.optical.observer.nonimaging.pixel import Pixel -from raysect.optical.observer.nonimaging.targetted_pixel import TargettedPixel +from raysect.optical.observer.nonimaging.targeted_pixel import TargetedPixel from raysect.optical.observer.nonimaging.mesh_pixel import MeshPixel from raysect.optical.observer.nonimaging.mesh_camera import MeshCamera diff --git a/raysect/optical/observer/nonimaging/__init__.py b/raysect/optical/observer/nonimaging/__init__.py index d1a8ee26..2b04ab5f 100644 --- a/raysect/optical/observer/nonimaging/__init__.py +++ b/raysect/optical/observer/nonimaging/__init__.py @@ -32,7 +32,7 @@ from .sightline import SightLine from .fibreoptic import FibreOptic from .pixel import Pixel -from .targetted_pixel import TargettedPixel +from .targeted_pixel import TargetedPixel from .mesh_pixel import MeshPixel from .mesh_camera import MeshCamera diff --git a/raysect/optical/observer/nonimaging/fibreoptic.pyx b/raysect/optical/observer/nonimaging/fibreoptic.pyx index 87969674..b9810367 100644 --- a/raysect/optical/observer/nonimaging/fibreoptic.pyx +++ b/raysect/optical/observer/nonimaging/fibreoptic.pyx @@ -40,7 +40,7 @@ cimport cython # 1 / (2 * PI) -DEF RECIP_2_PI = 0.15915494309189535 +cdef const double RECIP_2_PI = 0.15915494309189535 # TODO - provide a function for angular fall off for collection, instead of acceptance cone. diff --git a/raysect/optical/observer/nonimaging/mesh_camera.pyx b/raysect/optical/observer/nonimaging/mesh_camera.pyx index f47f5b83..60b1af97 100644 --- a/raysect/optical/observer/nonimaging/mesh_camera.pyx +++ b/raysect/optical/observer/nonimaging/mesh_camera.pyx @@ -43,14 +43,19 @@ from raysect.optical.observer.sampler1d import FullFrameSampler1D from raysect.primitive.mesh.mesh cimport Mesh cimport cython -# convenience defines -DEF X = 0 -DEF Y = 1 -DEF Z = 2 - -DEF V1 = 0 -DEF V2 = 1 -DEF V3 = 2 + +# constants +cdef enum: + + # axis + X = 0 + Y = 1 + Z = 2 + + # vertex + V1 = 0 + V2 = 1 + V3 = 2 cdef class MeshCamera(Observer1D): diff --git a/raysect/optical/observer/nonimaging/mesh_pixel.pyx b/raysect/optical/observer/nonimaging/mesh_pixel.pyx index 6391e417..9ed81706 100644 --- a/raysect/optical/observer/nonimaging/mesh_pixel.pyx +++ b/raysect/optical/observer/nonimaging/mesh_pixel.pyx @@ -43,14 +43,19 @@ from raysect.optical.observer.pipeline.spectral import SpectralPowerPipeline0D from raysect.primitive.mesh.mesh cimport Mesh cimport cython -# convenience defines -DEF X = 0 -DEF Y = 1 -DEF Z = 2 -DEF V1 = 0 -DEF V2 = 1 -DEF V3 = 2 +# constants +cdef enum: + + # axis + X = 0 + Y = 1 + Z = 2 + + # vertex + V1 = 0 + V2 = 1 + V3 = 2 cdef class MeshPixel(Observer0D): @@ -167,8 +172,8 @@ cdef class MeshPixel(Observer0D): self._collection_area = 0.0 self._cdf = np.zeros(self._triangles_mv.shape[0]) self._cdf_mv = self._cdf - - # calculate cumulative and total area simultaneously + + # calculate cumulative and total area simultaneously for i in range(self._triangles_mv.shape[0]): # obtain vertex indices @@ -184,7 +189,7 @@ cdef class MeshPixel(Observer0D): ) self._collection_area += triangle_area self._cdf_mv[i] = self._collection_area - + # normalise cumulative area to make cdf self._cdf /= self._collection_area diff --git a/raysect/optical/observer/nonimaging/meson.build b/raysect/optical/observer/nonimaging/meson.build index 8df56834..6ba3a60d 100644 --- a/raysect/optical/observer/nonimaging/meson.build +++ b/raysect/optical/observer/nonimaging/meson.build @@ -5,8 +5,8 @@ target_path = 'raysect/optical/observer/nonimaging' # source files py_files = ['__init__.py'] -pyx_files = ['fibreoptic.pyx', 'mesh_camera.pyx', 'mesh_pixel.pyx', 'pixel.pyx', 'sightline.pyx', 'targetted_pixel.pyx'] -pxd_files = ['__init__.pxd', 'fibreoptic.pxd', 'mesh_camera.pxd', 'mesh_pixel.pxd', 'pixel.pxd', 'sightline.pxd', 'targetted_pixel.pxd'] +pyx_files = ['fibreoptic.pyx', 'mesh_camera.pyx', 'mesh_pixel.pyx', 'pixel.pyx', 'sightline.pyx', 'targeted_pixel.pyx'] +pxd_files = ['__init__.pxd', 'fibreoptic.pxd', 'mesh_camera.pxd', 'mesh_pixel.pxd', 'pixel.pxd', 'sightline.pxd', 'targeted_pixel.pxd'] data_files = [] # compile cython diff --git a/raysect/optical/observer/nonimaging/targetted_pixel.pxd b/raysect/optical/observer/nonimaging/targeted_pixel.pxd similarity index 92% rename from raysect/optical/observer/nonimaging/targetted_pixel.pxd rename to raysect/optical/observer/nonimaging/targeted_pixel.pxd index dd9b7682..378a6a12 100644 --- a/raysect/optical/observer/nonimaging/targetted_pixel.pxd +++ b/raysect/optical/observer/nonimaging/targeted_pixel.pxd @@ -29,15 +29,15 @@ # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. -from raysect.core.math.sampler cimport RectangleSampler3D, HemisphereCosineSampler, TargettedHemisphereSampler +from raysect.core.math.sampler cimport RectangleSampler3D, HemisphereCosineSampler, TargetedHemisphereSampler from raysect.optical.observer.base cimport Observer0D -cdef class TargettedPixel(Observer0D): +cdef class TargetedPixel(Observer0D): cdef: - double _x_width, _y_width, _solid_angle, _collection_area, _targetted_path_prob + double _x_width, _y_width, _solid_angle, _collection_area, _targeted_path_prob tuple _targets RectangleSampler3D _point_sampler HemisphereCosineSampler _cosine_sampler - TargettedHemisphereSampler _targetted_sampler + TargetedHemisphereSampler _targeted_sampler diff --git a/raysect/optical/observer/nonimaging/targetted_pixel.pyx b/raysect/optical/observer/nonimaging/targeted_pixel.pyx similarity index 85% rename from raysect/optical/observer/nonimaging/targetted_pixel.pyx rename to raysect/optical/observer/nonimaging/targeted_pixel.pyx index d6b5048b..4a8d2c07 100644 --- a/raysect/optical/observer/nonimaging/targetted_pixel.pyx +++ b/raysect/optical/observer/nonimaging/targeted_pixel.pyx @@ -32,28 +32,28 @@ from libc.math cimport M_PI as PI from raysect.core.math.random cimport probability -from raysect.core.math.sampler cimport RectangleSampler3D, HemisphereCosineSampler, TargettedHemisphereSampler +from raysect.core.math.sampler cimport RectangleSampler3D, HemisphereCosineSampler, TargetedHemisphereSampler from raysect.optical cimport Ray, Primitive, Point3D, Vector3D, BoundingSphere3D from raysect.optical.observer.base cimport Observer0D from raysect.optical.observer.pipeline.spectral import SpectralPowerPipeline0D cimport cython -DEF R_2_PI = 0.15915494309189535 # 1 / (2 * pi) +cdef const double R_2_PI = 0.15915494309189535 # 1 / (2 * pi) -cdef class TargettedPixel(Observer0D): +cdef class TargetedPixel(Observer0D): """ A pixel observer that preferentially targets rays towards a given list of primitives. - The targetted pixel takes a list of target primitives. The observer targets the + The targeted pixel takes a list of target primitives. The observer targets the bounding sphere that encompasses a target primitive. Therefore, for best performance, the target primitives should be split up such that their surfaces are closely wrapped by the bounding sphere. The sampling algorithm fires a proportion of rays at the targets, and a portion sampled from the full hemisphere. The proportion that is fired towards the targets is controlled - with the targetted_path_prob attribute. By default this attribute is set to 0.9, i.e. + with the targeted_path_prob attribute. By default this attribute is set to 0.9, i.e. 90% of the rays are fired towards the targets. .. Warning.. @@ -62,7 +62,7 @@ cdef class TargettedPixel(Observer0D): targeted directions, otherwise they will not be sampled and the result will be biased. :param list targets: The list of primitives for targeted sampling. - :param float targetted_path_prob: The probability of sampling a targeted primitive VS sampling over the whole hemisphere. + :param float targeted_path_prob: The probability of sampling a targeted primitive VS sampling over the whole hemisphere. :param list pipelines: The list of pipelines that will process the spectrum measured by this pixel (default=SpectralPipeline0D()). :param float x_width: The rectangular collection area's width along the @@ -73,21 +73,21 @@ cdef class TargettedPixel(Observer0D): .. code-block:: pycon - >>> from raysect.optical.observer import TargettedPixel, PowerPipeline0D + >>> from raysect.optical.observer import TargetedPixel, PowerPipeline0D >>> >>> # set-up scenegraph >>> world = World() >>> emitter = Sphere(radius=sphere_radius, parent=world) >>> emitter.material = UnityVolumeEmitter() >>> - >>> # setup targetted pixel observer - >>> targetted_pipeline = PowerPipeline0D(name="Targeted Pixel Observer") - >>> targetted_pixel = TargettedPixel(parent=world, targets=[emitter], - >>> pixel_samples=250, pipelines=[targetted_pipeline]) - >>> targetted_pixel.observe() + >>> # setup targeted pixel observer + >>> targeted_pipeline = PowerPipeline0D(name="Targeted Pixel Observer") + >>> targeted_pixel = TargetedPixel(parent=world, targets=[emitter], + >>> pixel_samples=250, pipelines=[targeted_pipeline]) + >>> targeted_pixel.observe() """ - def __init__(self, targets, targetted_path_prob=None, + def __init__(self, targets, targeted_path_prob=None, pipelines=None, x_width=None, y_width=None, parent=None, transform=None, name=None, render_engine=None, pixel_samples=None, samples_per_task=None, spectral_rays=None, spectral_bins=None, min_wavelength=None, max_wavelength=None, ray_extinction_prob=None, ray_extinction_min_depth=None, @@ -105,14 +105,14 @@ cdef class TargettedPixel(Observer0D): self._x_width = 0.01 self._y_width = 0.01 self._cosine_sampler = HemisphereCosineSampler() - self._targetted_sampler = None + self._targeted_sampler = None self._solid_angle = 2 * PI self.x_width = x_width or 0.01 self.y_width = y_width or 0.01 self.targets = targets - self.targetted_path_prob = targetted_path_prob or 0.9 + self.targeted_path_prob = targeted_path_prob or 0.9 @property def x_width(self): @@ -205,7 +205,7 @@ cdef class TargettedPixel(Observer0D): self._targets = value @property - def targetted_path_prob(self): + def targeted_path_prob(self): """ The probability that an individual sample will be fired at a target instead of a sample from the whole hemisphere. @@ -216,14 +216,14 @@ cdef class TargettedPixel(Observer0D): :rtype: float """ - return self._targetted_path_prob + return self._targeted_path_prob - @targetted_path_prob.setter - def targetted_path_prob(self, double value): + @targeted_path_prob.setter + def targeted_path_prob(self, double value): if value < 0 or value > 1: raise ValueError("Targeted path probability must lie in the range [0, 1].") - self._targetted_path_prob = value + self._targeted_path_prob = value @cython.boundscheck(False) @cython.wraparound(False) @@ -249,8 +249,8 @@ cdef class TargettedPixel(Observer0D): sphere = target.bounding_sphere() spheres.append((sphere.centre.transform(self.to_local()), sphere.radius, 1.0)) - # instance targetted pixel sampler - self._targetted_sampler = TargettedHemisphereSampler(spheres) + # instance targeted pixel sampler + self._targeted_sampler = TargetedHemisphereSampler(spheres) # sample pixel origins origins = self._point_sampler.samples(ray_count) @@ -258,17 +258,17 @@ cdef class TargettedPixel(Observer0D): rays = [] for origin in origins: - if probability(self._targetted_path_prob): - # obtain targetted vector sample - direction = self._targetted_sampler.sample(origin) + if probability(self._targeted_path_prob): + # obtain targeted vector sample + direction = self._targeted_sampler.sample(origin) else: # obtain cosine weighted hemisphere sample direction = self._cosine_sampler.sample() # calculate combined pdf and ray weight - pdf = self._targetted_path_prob * self._targetted_sampler.pdf(origin, direction) + \ - (1-self._targetted_path_prob) * self._cosine_sampler.pdf(direction) + pdf = self._targeted_path_prob * self._targeted_sampler.pdf(origin, direction) + \ + (1-self._targeted_path_prob) * self._cosine_sampler.pdf(direction) if pdf <= 0: raise ValueError('Ray direction probability is zero. The target object extends beyond the pixel horizon.') diff --git a/raysect/optical/observer/pipeline/bayer.pxd b/raysect/optical/observer/pipeline/bayer.pxd index c31c915f..41d50d31 100644 --- a/raysect/optical/observer/pipeline/bayer.pxd +++ b/raysect/optical/observer/pipeline/bayer.pxd @@ -32,6 +32,7 @@ from raysect.optical.spectralfunction cimport SpectralFunction from raysect.optical.observer.base cimport PixelProcessor, Pipeline2D from raysect.core.math cimport StatsArray2D + cdef class BayerPipeline2D(Pipeline2D): cdef: @@ -56,17 +57,10 @@ cdef class BayerPipeline2D(Pipeline2D): bint _quiet cpdef object _start_display(self) - cpdef object _update_display(self, int x, int y) - cpdef object _refresh_display(self) - cpdef object _render_display(self, StatsArray2D frame, str status=*) - cpdef np.ndarray _generate_display_image(self, StatsArray2D frame) - cpdef double _calculate_white_point(self, np.ndarray image) - cpdef object display(self) - cpdef object save(self, str filename) diff --git a/raysect/optical/observer/pipeline/mono/power.pxd b/raysect/optical/observer/pipeline/mono/power.pxd index 12129362..347fee1a 100644 --- a/raysect/optical/observer/pipeline/mono/power.pxd +++ b/raysect/optical/observer/pipeline/mono/power.pxd @@ -83,19 +83,12 @@ cdef class PowerPipeline2D(Pipeline2D): bint _quiet cpdef object _start_display(self) - cpdef object _update_display(self, int x, int y) - cpdef object _refresh_display(self) - cpdef object _render_display(self, StatsArray2D frame, str status=*) - cpdef np.ndarray _generate_display_image(self, StatsArray2D frame) - cpdef double _calculate_white_point(self, np.ndarray image) - cpdef object display(self) - cpdef object save(self, str filename) diff --git a/raysect/optical/observer/pipeline/rgb.pxd b/raysect/optical/observer/pipeline/rgb.pxd index 8f40c6e0..35326988 100644 --- a/raysect/optical/observer/pipeline/rgb.pxd +++ b/raysect/optical/observer/pipeline/rgb.pxd @@ -54,21 +54,13 @@ cdef class RGBPipeline2D(Pipeline2D): bint _quiet cpdef object _start_display(self) - cpdef object _update_display(self, int x, int y) - cpdef object _refresh_display(self) - cpdef object _render_display(self, StatsArray3D frame, str status=*) - cpdef np.ndarray _generate_display_image(self, StatsArray3D frame) - cpdef double _calculate_sensitivity(self, np.ndarray image) - cpdef np.ndarray _generate_srgb_image(self, double[:,:,::1] xyz_image_mv) - cpdef object display(self) - cpdef object save(self, str filename) diff --git a/raysect/optical/ray.pyx b/raysect/optical/ray.pyx index 499289d6..78f07c64 100644 --- a/raysect/optical/ray.pyx +++ b/raysect/optical/ray.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 M_PI as PI, asin, cos +from libc.math cimport M_PI as PI, asin, cos, INFINITY from raysect.core cimport Intersection from raysect.core.math.random cimport probability @@ -39,9 +39,6 @@ from raysect.optical.spectrum cimport new_spectrum from raysect.optical.scenegraph cimport Primitive cimport cython -# cython doesn't have a built-in infinity constant, this compiles to +infinity -DEF INFINITY = 1e999 - cdef class Ray(CoreRay): """ diff --git a/raysect/optical/scenegraph/world.pxd b/raysect/optical/scenegraph/world.pxd index be376b6b..0249c64d 100644 --- a/raysect/optical/scenegraph/world.pxd +++ b/raysect/optical/scenegraph/world.pxd @@ -42,27 +42,18 @@ cdef class ImportanceManager: list _spheres cdef object _process_primitives(self, list primitives) - cdef object _calculate_cdf(self) - cdef tuple _pick_sphere(self) - cpdef Vector3D sample(self, Point3D origin) - cpdef double pdf(self, Point3D origin, Vector3D direction) - cpdef bint has_primitives(self) cdef class World(CoreWorld): - cdef: - ImportanceManager _importance + cdef ImportanceManager _importance cpdef build_importance(self, bint force=*) - cpdef Vector3D important_direction_sample(self, Point3D origin) - cpdef double important_direction_pdf(self, Point3D origin, Vector3D direction) - cpdef bint has_important_primitives(self) diff --git a/raysect/optical/spectrum.pxd b/raysect/optical/spectrum.pxd index c59ce20f..7c465e56 100644 --- a/raysect/optical/spectrum.pxd +++ b/raysect/optical/spectrum.pxd @@ -32,6 +32,7 @@ from raysect.optical.spectralfunction cimport SpectralFunction from numpy cimport ndarray + cdef class Spectrum(SpectralFunction): cdef: diff --git a/raysect/optical/spectrum.pyx b/raysect/optical/spectrum.pyx index 3a9ecac6..bb8ad335 100644 --- a/raysect/optical/spectrum.pyx +++ b/raysect/optical/spectrum.pyx @@ -34,7 +34,7 @@ from raysect.core.math.cython cimport integrate, interpolate from numpy cimport PyArray_SimpleNew, PyArray_FILLWBYTE, NPY_FLOAT64, npy_intp, import_array # Plank's constant * speed of light in a vacuum -DEF CONSTANT_HC = 1.9864456832693028e-25 +cdef const double CONSTANT_HC = 1.9864456832693028e-25 # required by numpy c-api import_array() diff --git a/raysect/primitive/box.pyx b/raysect/primitive/box.pyx index 8aa92536..3f585b1b 100644 --- a/raysect/primitive/box.pyx +++ b/raysect/primitive/box.pyx @@ -30,28 +30,27 @@ # POSSIBILITY OF SUCH DAMAGE. from raysect.core cimport new_point3d, Normal3D, new_normal3d, AffineMatrix3D, Material, new_intersection, BoundingBox3D -from libc.math cimport fabs +from libc.math cimport fabs, INFINITY cimport cython -# cython doesn't have a built-in infinity constant, this compiles to +infinity -DEF INFINITY = 1e999 - # bounding box is padded by a small amount to avoid numerical accuracy issues -DEF BOX_PADDING = 1e-9 +cdef const double BOX_PADDING = 1e-9 # additional ray distance to avoid re-hitting the same surface point -DEF EPSILON = 1e-9 - -# slab face enumeration -DEF NO_FACE = -1 -DEF LOWER_FACE = 0 -DEF UPPER_FACE = 1 - -# axis enumeration -DEF NO_AXIS = -1 -DEF X_AXIS = 0 -DEF Y_AXIS = 1 -DEF Z_AXIS = 2 +cdef const double EPSILON = 1e-9 + +cdef enum: + + # slab face enumeration + NO_FACE = -1 + LOWER_FACE = 0 + UPPER_FACE = 1 + + # axis enumeration + NO_AXIS = -1 + X_AXIS = 0 + Y_AXIS = 1 + Z_AXIS = 2 cdef class Box(Primitive): diff --git a/raysect/primitive/cone.pyx b/raysect/primitive/cone.pyx index 17355c8d..cff9592d 100644 --- a/raysect/primitive/cone.pyx +++ b/raysect/primitive/cone.pyx @@ -31,22 +31,20 @@ from raysect.core cimport new_point3d, Point3D, new_normal3d, AffineMatrix3D, Material, new_intersection, BoundingBox3D from raysect.core.math.cython cimport solve_quadratic, swap_double, swap_int -from libc.math cimport sqrt +from libc.math cimport sqrt, INFINITY cimport cython -# cython doesn't have a built-in infinity constant, this compiles to +infinity -DEF INFINITY = 1e999 - # bounding box is padded by a small amount to avoid numerical accuracy issues -DEF BOX_PADDING = 1e-9 +cdef const double BOX_PADDING = 1e-9 # additional ray distance to avoid re-hitting the same surface point -DEF EPSILON = 1e-9 +cdef const double EPSILON = 1e-9 # object type enumeration -DEF NO_TYPE = -1 -DEF CONE = 0 -DEF BASE = 1 +cdef enum: + NO_TYPE = -1 + CONE = 0 + BASE = 1 cdef class Cone(Primitive): diff --git a/raysect/primitive/csg.pyx b/raysect/primitive/csg.pyx index 6b760085..e32cbd29 100644 --- a/raysect/primitive/csg.pyx +++ b/raysect/primitive/csg.pyx @@ -31,13 +31,12 @@ # TODO: add more advanced material handling +from libc.math cimport INFINITY + from raysect.core cimport _NodeBase, ChangeSignal, Material, new_ray, new_intersection, Point3D, AffineMatrix3D, BoundingBox3D # bounding box is padded by a small amount to avoid numerical accuracy issues -DEF BOX_PADDING = 1e-9 - -# cython doesn't have a built in infinity definition -DEF INFINITY = 1e999 +cdef const double BOX_PADDING = 1e-9 cdef class CSGPrimitive(Primitive): @@ -112,7 +111,7 @@ cdef class CSGPrimitive(Primitive): def primitive_b(self): """ Component primitive B of the compound CSG primitive. - + :rtype: Primitive """ return self._primitive_b.primitive diff --git a/raysect/primitive/cylinder.pyx b/raysect/primitive/cylinder.pyx index 5840cc9f..f7ddd8d7 100644 --- a/raysect/primitive/cylinder.pyx +++ b/raysect/primitive/cylinder.pyx @@ -31,27 +31,26 @@ from raysect.core cimport AffineMatrix3D, new_normal3d, new_point3d, new_vector3d, Material, new_intersection, BoundingBox3D from raysect.core.math.cython cimport solve_quadratic, swap_double -from libc.math cimport sqrt, fabs +from libc.math cimport sqrt, fabs, INFINITY cimport cython -# cython doesn't have a built-in infinity constant, this compiles to +infinity -DEF INFINITY = 1e999 - # bounding box is padded by a small amount to avoid numerical accuracy issues -DEF BOX_PADDING = 1e-9 +cdef const double BOX_PADDING = 1e-9 # additional ray distance to avoid re-hitting the same surface point -DEF EPSILON = 1e-9 +cdef const double EPSILON = 1e-9 + +cdef enum: -# object type enumeration -DEF NO_TYPE = -1 -DEF CYLINDER = 0 -DEF SLAB = 1 + # object type enumeration + NO_TYPE = -1 + CYLINDER = 0 + SLAB = 1 -# slab face enumeration -DEF NO_FACE = -1 -DEF LOWER_FACE = 0 -DEF UPPER_FACE = 1 + # slab face enumeration + NO_FACE = -1 + LOWER_FACE = 0 + UPPER_FACE = 1 cdef class Cylinder(Primitive): diff --git a/raysect/primitive/lens/spherical.pyx b/raysect/primitive/lens/spherical.pyx index 545bbf76..4445b2d3 100644 --- a/raysect/primitive/lens/spherical.pyx +++ b/raysect/primitive/lens/spherical.pyx @@ -40,7 +40,7 @@ from libc.math cimport sqrt Basic spherical lens primitives. """ -DEF PADDING = 0.000001 +cdef const double PADDING = 0.000001 cdef class BiConvex(EncapsulatedPrimitive): @@ -128,7 +128,7 @@ cdef class BiConvex(EncapsulatedPrimitive): cdef bint _is_short(self): """ - Do the facing spheres overlap sufficiently to build a lens using just their intersection? + Do the facing spheres overlap sufficiently to build a lens using just their intersection? """ cdef double available_thickness = min( @@ -326,7 +326,7 @@ cdef class PlanoConvex(EncapsulatedPrimitive): # attach to local root (performed in EncapsulatedPrimitive init) super().__init__(lens, parent, transform, material, name) - + cdef void _calc_geometry(self): cdef double radius, radius_sqr @@ -343,7 +343,7 @@ cdef class PlanoConvex(EncapsulatedPrimitive): cdef bint _is_short(self): """ - Does the front sphere have sufficient radius to build the lens with just an intersection? + Does the front sphere have sufficient radius to build the lens with just an intersection? """ cdef double available_thickness = 2 * (self.curvature - self.curve_thickness) @@ -548,7 +548,7 @@ cdef class Meniscus(EncapsulatedPrimitive): cdef bint _is_short(self): """ - Does the front sphere have sufficient radius to build the lens with just an intersection? + Does the front sphere have sufficient radius to build the lens with just an intersection? """ cdef double available_thickness = 2 * self.front_curvature - self.front_thickness diff --git a/raysect/primitive/mesh/mesh.pxd b/raysect/primitive/mesh/mesh.pxd index 8c209581..7deb725a 100644 --- a/raysect/primitive/mesh/mesh.pxd +++ b/raysect/primitive/mesh/mesh.pxd @@ -60,37 +60,21 @@ cdef class MeshData(KDTree3DCore): int32_t _i cpdef Point3D vertex(self, int index) - cpdef ndarray triangle(self, int index) - cpdef Normal3D vertex_normal(self, int index) - cpdef Normal3D face_normal(self, int index) - cdef object _filter_triangles(self) - cdef object _flip_normals(self) - cdef object _generate_face_normals(self) - cdef BoundingBox3D _generate_bounding_box(self, int32_t i) - cdef void _calc_rayspace_transform(self, Ray ray) - cdef bint _hit_triangle(self, int32_t i, Ray ray, float[4] hit_data) - cpdef Intersection calc_intersection(self, Ray ray) - cdef Normal3D _intersection_normal(self) - cpdef bint contains(self, Point3D p) - cpdef BoundingBox3D bounding_box(self, AffineMatrix3D to_world) - cdef uint8_t _read_uint8(self, object file) - cdef bint _read_bool(self, object file) - cdef double _read_float(self, object file) diff --git a/raysect/primitive/mesh/mesh.pyx b/raysect/primitive/mesh/mesh.pyx index fabb03b5..048e21c2 100644 --- a/raysect/primitive/mesh/mesh.pyx +++ b/raysect/primitive/mesh/mesh.pyx @@ -35,7 +35,7 @@ import struct from numpy import array, float32, int32, zeros from raysect.core cimport Primitive, AffineMatrix3D, Normal3D, new_normal3d, Point3D, new_point3d, Vector3D, new_vector3d, Material, Ray, new_ray, Intersection, new_intersection, BoundingBox3D, new_boundingbox3d from raysect.core.math.spatial cimport KDTree3DCore, Item3D -from libc.math cimport fabs +from libc.math cimport fabs, INFINITY from numpy cimport float32_t, int32_t, uint8_t from cpython.bytes cimport PyBytes_AsString cimport cython @@ -45,37 +45,41 @@ The ray-triangle intersection used for the Mesh primitive is an implementation o "Watertight Ray/Triangle Intersection", S.Woop, C.Benthin, I.Wald, Journal of Computer Graphics Techniques (2013), Vol.2, No. 1 """ -# cython doesn't have a built-in infinity constant, this compiles to +infinity -DEF INFINITY = 1e999 - # bounding box is padded by a small amount to avoid numerical accuracy issues -DEF BOX_PADDING = 1e-6 +cdef const double BOX_PADDING = 1e-6 # additional ray distance to avoid re-hitting the same surface point -DEF EPSILON = 1e-6 +cdef const double EPSILON = 1e-6 + +# constants +cdef enum: -# convenience defines -DEF X = 0 -DEF Y = 1 -DEF Z = 2 + # axis index + X = 0 + Y = 1 + Z = 2 -DEF U = 0 -DEF V = 1 -DEF W = 2 -DEF T = 3 + # intersection data + U = 0 + V = 1 + W = 2 + T = 3 -DEF V1 = 0 -DEF V2 = 1 -DEF V3 = 2 -DEF N1 = 3 -DEF N2 = 4 -DEF N3 = 5 + # vertex index + V1 = 0 + V2 = 1 + V3 = 2 -DEF NO_INTERSECTION = -1 + # normal index + N1 = 3 + N2 = 4 + N3 = 5 -# raysect mesh format constants -DEF RSM_VERSION_MAJOR = 1 -DEF RSM_VERSION_MINOR = 0 + NO_INTERSECTION = -1 + + # raysect mesh format constants + RSM_VERSION_MAJOR = 1 + RSM_VERSION_MINOR = 0 cdef class MeshIntersection(Intersection): @@ -279,9 +283,9 @@ cdef class MeshData(KDTree3DCore): cpdef Point3D vertex(self, int index): """ Returns the specified vertex. - + :param index: The vertex index. - :return: A Point3D object. + :return: A Point3D object. """ if index < 0 or index >= self.vertices_mv.shape[0]: @@ -296,13 +300,13 @@ cdef class MeshData(KDTree3DCore): cpdef ndarray triangle(self, int index): """ Returns the specified triangle. - - The returned data will either be a 3 or 6 element numpy array. The + + The returned data will either be a 3 or 6 element numpy array. The first three element are the triangle's vertex indices. If present, the last three elements are the triangle's vertex normal indices. - + :param index: The triangle index. - :return: A numpy array. + :return: A numpy array. """ if index < 0 or index >= self.vertices_mv.shape[0]: @@ -316,9 +320,9 @@ cdef class MeshData(KDTree3DCore): cpdef Normal3D vertex_normal(self, int index): """ Returns the specified vertex normal. - + :param index: The vertex normal's index. - :return: A Normal3D object. + :return: A Normal3D object. """ if self._vertex_normals is None: @@ -339,9 +343,9 @@ cdef class MeshData(KDTree3DCore): cpdef Normal3D face_normal(self, int index): """ Returns the specified face normal. - + :param index: The face normal's index. - :return: A Normal3D object. + :return: A Normal3D object. """ if index < 0 or index >= self.face_normals_mv.shape[0]: diff --git a/raysect/primitive/parabola.pyx b/raysect/primitive/parabola.pyx index ac1a2c0e..47cc257c 100644 --- a/raysect/primitive/parabola.pyx +++ b/raysect/primitive/parabola.pyx @@ -31,23 +31,21 @@ from raysect.core cimport new_point3d, Point3D, new_normal3d, AffineMatrix3D, Material, new_intersection, BoundingBox3D from raysect.core.math.cython cimport solve_quadratic, swap_double, swap_int -from libc.math cimport sqrt +from libc.math cimport sqrt, INFINITY cimport cython -# cython doesn't have a built-in infinity constant, this compiles to +infinity -DEF INFINITY = 1e999 - # bounding box is padded by a small amount to avoid numerical accuracy issues -DEF BOX_PADDING = 1e-9 +cdef const double BOX_PADDING = 1e-9 # TODO - Perhaps should be calculated based on primitive scale # additional ray distance to avoid re-hitting the same surface point -DEF EPSILON = 1e-9 +cdef const double EPSILON = 1e-9 # object type enumeration -DEF NO_TYPE = -1 -DEF PARABOLA = 0 -DEF BASE = 1 +cdef enum: + NO_TYPE = -1 + PARABOLA = 0 + BASE = 1 cdef class Parabola(Primitive): diff --git a/raysect/primitive/sphere.pyx b/raysect/primitive/sphere.pyx index d5bbd650..a1fa13f3 100644 --- a/raysect/primitive/sphere.pyx +++ b/raysect/primitive/sphere.pyx @@ -35,11 +35,11 @@ from raysect.core.math.cython cimport solve_quadratic, swap_double # bounding box and sphere are padded by small amounts to avoid numerical accuracy issues -DEF BOX_PADDING = 1e-9 -DEF SPHERE_PADDING = 1.000000001 +cdef const double BOX_PADDING = 1e-9 +cdef const double SPHERE_PADDING = 1.000000001 # additional ray distance to avoid re-hitting the same surface point -DEF EPSILON = 1e-9 +cdef const double EPSILON = 1e-9 cdef class Sphere(Primitive): diff --git a/raysect/primitive/utility.pxd b/raysect/primitive/utility.pxd index 947645ac..cce6cc05 100644 --- a/raysect/primitive/utility.pxd +++ b/raysect/primitive/utility.pxd @@ -38,9 +38,6 @@ cdef class EncapsulatedPrimitive(Primitive): Primitive _primitive cpdef Intersection hit(self, Ray ray) - cpdef Intersection next_intersection(self) - cpdef bint contains(self, Point3D p) except -1 - cpdef BoundingBox3D bounding_box(self)