core icon indicating copy to clipboard operation
core copied to clipboard

Add RayTransferEmitter that uses a function to map a point in space to a light source

Open vsnever opened this issue 4 years ago • 0 comments

I would like to add a ray transfer emitter, which would use a custom function to map a point in 3D space to a unit light source in the array. This ray transfer emitter will be compatible with Discrete2DMesh, Discrete3DMesh, SOLPSFunction2D, etc.

Here is an implementation:

cdef class FunctionRayTransferEmitter(InhomogeneousVolumeEmitter):
    """
    Ray transfer emitter defined with a function, which returns an index of the light source
    for a given point in space.

    :param Function3D raytransfer_function: A function that maps a point in 3D space to
        an array of volumetric light sources with unit emissivity.
        Must return the index of the light source in the array starting with 1.
        Zero means that the point is outside the volume for all light sources
        in the array.
    :param raysect.optical.material.VolumeIntegrator integrator: Volume integrator,
        defaults to `integrator=NumericalIntegrator()`.
    """

    cdef:
        readonly Function3D raytransfer_function

    def __init__(self, raytransfer_function, integrator=None):

        super().__init__(integrator)
        self.raytransfer_function = autowrap_function3d(raytransfer_function)

    @cython.boundscheck(False)
    @cython.wraparound(False)
    @cython.initializedcheck(False)
    cpdef Spectrum emission_function(self, Point3D point, Vector3D direction, Spectrum spectrum,
                                     World world, Ray ray, Primitive primitive,
                                     AffineMatrix3D world_to_local, AffineMatrix3D local_to_world):

        cdef int isource

        isource = <int> libc_round(self.raytransfer_function.evaluate(point.x, point.y, point.z))
        if isource > 0:
            spectrum.samples_mv[isource] += 1.

        return spectrum

Usage case:

import numpy as np
from Py3DViewer import Tetmesh
from raysect.optical import World, translate, Point3D, Vector3D, rotate_basis
from raysect.core.math.function.float.function3d.interpolate import Discrete3DMesh
from raysect.optical.observer import PinholeCamera, FullFrameSampler2D
from raysect.primitive import Box
from cherab.tools.raytransfer import FunctionRayTransferEmitter, RayTransferPipeline2D, RoughTitanium

mesh = Tetmesh("path_to_raysect_demos/resources/stanford_bunny.mesh")
try:
    tetra_index = Discrete3DMesh(mesh.vertices, mesh.polys, np.arange(1, mesh.num_polys + 1, dtype=float), False, 0)
    bins = mesh.num_polys
except AttributeError:  # for backward compatibility with older versions of Py3DViewer
    tetra_index = Discrete3DMesh(mesh.vertices, mesh.tets, np.arange(1, mesh.num_tets + 1, dtype=float), False, 0)
    bins = mesh.num_tets

world = World()
lower = Point3D(mesh.bbox[0, 0], mesh.bbox[0, 1], mesh.bbox[0, 2])
upper = Point3D(mesh.bbox[1, 0], mesh.bbox[1, 1], mesh.bbox[1, 2])
emitter = Box(lower, upper, material=FunctionRayTransferEmitter(tetra_index), parent=world)
floor = Box(Point3D(-100, -0.1, -100), Point3D(100, -0.01, 100), world, material=RoughTitanium(0.1))

camera_pos = Point3D(-0.015, 0.204, 0.25)
camera_translate = translate(camera_pos.x, camera_pos.y, camera_pos.z)
camera_rotate = rotate_basis(camera_pos.vector_to(Point3D(0, 0.06, 0)), Vector3D(0, 1, 0))
pipeline = RayTransferPipeline2D()
camera = PinholeCamera((32, 32), parent=world, transform=camera_translate * camera_rotate,
                                   pipelines=[pipeline], frame_sampler=FullFrameSampler2D())
camera.fov = 50
camera.pixel_samples = 1000
camera.spectral_bins = bins
# ray transfer matrix will be calculated for 600.5 nm
camera.min_wavelength = 600.
camera.max_wavelength = 601.

# integration resolution
emitter.material.integrator.step = 0.001

camera.observe()

However, for the observers with the large number of pixels, the combination of FunctionRayTransferEmitter and Discrete3(2)DMesh cannot replace the combination of Cartesian(Cylindrical)RayTransferEmitter and voxel_map due to significantly lower computational speed.

I will create a PR when #307 is merged to avoid merging conflicts.

vsnever avatar Oct 21 '21 15:10 vsnever