10.6. Tomography

Tomographic inversion problems often arise in the study of plasma diagnostics, especially when making line integrated measurements. Tomographic inversion techniques allow us to try recovering the underlying plasma properties from a reduced set of measurements. In Cherab we implement a basic set of tomography algorithms, for a wider range of options please consult the dedicated tomography libraries such as ToFu.

In general, recovering a plasma emission profile with tomography is an ill-posed problem. It is customary to describe the system in terms of a sensitivity matrix \(\mathbf{W}\). The elements \(W_{k,l}\) describe the coupling between the \(N_s\) plasma emission voxels \(x_l\) and measured power \(\Phi_k\) at \(N_d\) detectors. The whole detector set is typically represented as the matrix equation

\[\mathbf{\Phi} = \mathbf{W} \mathbf{x}.\]

The power for the k th detector can be expressed as

\[\Phi_k = \sum_{l=1}^{N_s} W_{k,l} \, x_l,\]

where \(k\) and \(l\) are the indices for the detectors and source voxels respectively.

In this module we implement a number of basic inversion algorithms for recovering the emissivity vector \(\mathbf{x}\) from a set of measurements \(\mathbf{\Phi}\) and sensitivity matrix \(\mathbf{W}\).

10.6.1. Inversion Methods

cherab.tools.inversions.sart.invert_sart(geometry_matrix, measurement_vector, initial_guess, max_iterations, relaxation, conv_tol)

Performs a SART inversion on the specified measurement vector.

This function implements the Simultaneous Algebraic Reconstruction Technique (SART), as published in A. Andersen, and A. Kak, Ultrasonic imaging 6, 81 (1984). The SART method is an iterative inversion scheme where the source cells are updated with the formula

\[\begin{split}x_l^{(i+1)} = f_{sart}(x_l^{(i)}) = x_l^{(i)} + \\frac{\omega}{W_{\oplus,l}} \sum_{k=1}^{N_d} \\frac{W_{k,l}}{W_{k,\oplus}} (\Phi_k - \hat{\Phi}_k),\end{split}\]

where

\[W_{k,\oplus} = \sum_{l=1}^{N_s} W_{k,l}, \quad W_{\oplus, l} = \sum_{k=1}^{N_d} W_{k,l}.\]

Here \(x_l^{(i)}\) is the previous estimate for the emission at voxel \(l\) in iteration \(i\). The SART method effectively updates each cell by the weighted average error between the forward modelled \(\hat{\Phi}_k\) and observed \(\Phi_k\) measurements. The observed errors are weighted by both their proportion of the total ray length (\(W_{k,\oplus}\)) and the sum of the effective ray paths crossing that cell (\(W_{\oplus, l}\)).

Parameters:

  • geometry_matrix (np.ndarray) – The sensitivity matrix describing the coupling between the detectors and the voxels. Must be an array with shape \((N_d, N_s)\).

  • measurement_vector (np.ndarray) – The measured power/radiance vector with shape \((N_d)\).

  • initial_guess – An optional initial guess, can be an array of shape \((N_s)\) or a constant value that will be used to seed the algorithm.

  • max_iterations (int) – The maximum number of iterations to run the SART algorithm before returning a result, defaults to max_iterations=250.

  • relaxation (float) – The relaxation hyperparameter, defaults to relaxation=1. Consult the reference papers for more information on this hyperparameter.

  • conv_tol (float) – The convergence limit at which the algorithm will be terminated, unless the maximum number of iterations has been reached. The convergence is calculated as the normalised squared difference between the measurement and solution vectors.

Returns:

A tuple with the inverted solution vector \(\mathbf{x}\) as an ndarray with shape \((N_s)\), and the convergence achieved as a float.

>>> from cherab.tools.inversions import invert_sart
>>> inverted_solution, conv = invert_sart(weight_matrix, observations, max_iterations=100)
cherab.tools.inversions.sart.invert_constrained_sart(geometry_matrix, laplacian_matrix, measurement_vector, initial_guess, max_iterations, relaxation, beta_laplace, conv_tol)

Performs a constrained SART inversion on the specified measurement vector.

The core of the constrained SART algorithm is identical to the basic SART algorithm implemented in invert_sart(). The only difference is that now the iterative update formula includes a regularisation operator.

\[x_l^{(i+1)} = f_{sart}(x_l^{(i)}) - \hat{\mathcal{L}}_{iso}(x_l^{(i)}).\]

In this particular function we have implemented a isotropic Laplacian smoothness operator,

\[\begin{split}\hat{\mathcal{L}}_{iso}(x_l^{(i)}) = \\beta_L (Cx_l^{(i)} - \sum_{c=1}^C x_c^{(i)}).\end{split}\]

Here, \(c\) is the index for the sum over the neighbouring voxels. The regularisation hyperparameter \(\\beta_L\) determines the amount of local smoothness imposed on the solution vector. When \(\\beta_L = 0\), the solution is fully determined by the measurements, and as \(\\beta_L \\rightarrow 1\), the solution is dominated by the smoothness operator.

Parameters:

  • geometry_matrix (np.ndarray) – The sensitivity matrix describing the coupling between the detectors and the voxels. Must be an array with shape \((N_d, N_s)\).

  • laplacian_matrix (np.ndarray) – The laplacian regularisation matrix of shape \((N_s, N_s)\).

  • measurement_vector (np.ndarray) – The measured power/radiance vector with shape \((N_d)\).

  • initial_guess – An optional initial guess, can be an array of shape \((N_s)\) or a constant value that will be used to seed the algorithm.

  • max_iterations (int) – The maximum number of iterations to run the SART algorithm before returning a result, defaults to max_iterations=250.

  • relaxation (float) – The relaxation hyperparameter, defaults to relaxation=1. Consult the reference papers for more information on this hyperparameter.

  • beta_laplace (float) – The regularisation hyperparameter in the range [0, 1]. Defaults to beta_laplace=0.01.

  • conv_tol (float) – The convergence limit at which the algorithm will be terminated, unless the maximum number of iterations has been reached. The convergence is calculated as the normalised squared difference between the measurement and solution vectors.

Returns:

A tuple with the inverted solution vector \(\mathbf{x}\) as an ndarray with shape \((N_s)\), and the convergence achieved as a float.

>>> from cherab.tools.inversions import invert_constrained_sart
>>> inverted_solution, conv = invert_constrained_sart(weight_matrix, laplacian, observations)
class cherab.tools.inversions.opencl.sart_opencl.SartOpencl(geometry_matrix, laplacian_matrix=None, device=None, block_size=256, copy_column_major=True, block_size_row_maj=64, use_atomic=True, steps_per_thread=64, verbose=False)

A GPU-accelerated version of SART inversion. The geometry matrix and Laplacian matrix are provided on initialisation because they must be copied to GPU memory, which takes time. Inversions may be performed multiple times for different measurement vectors without copying the matrices each time. If required, the Laplacian matrix can be updated by calling update_laplacian_matrix(new_laplacian_matrix) method. Note that computations are performed with single precision.

Parameters:

  • geometry_matrix (np.ndarray) – The sensitivity matrix describing the coupling between the detectors and the voxels. Must be an array with shape \((N_d, N_s)\).

  • laplacian_matrix (np.ndarray) – The laplacian regularisation matrix of shape \((N_s, N_s)\). Default value: laplacian_matrix=None.

  • device (pyopencl.Device) – OpenCL device which will be used for computations. Default value: device=None (autoselect).

  • block_size (int) – Number of GPU threads per block. Must be the power of 2. For the best performance try from 256 to 1024 for Nvidia (use 1024 on high-end GPUs), from 64 to 256 for AMD and from 16 to 64 for Intel GPUs. Default value: block_size=256.

  • copy_column_major (bool) – If True, the two copies of geometry matrix will be stored in GPU memory. One in row-major order and the other one in column-major order. This provides much better performance of the inversions but requires twice as much GPU memory. Default value: copy_column_major=True.

  • block_size_row_maj (int) – If copy_column_major is set to False, this parameter defines the number of GPU threads per block in mat_vec_mult_row_maj() kernel used to calculate y_hat. Must be lower than block_size. Default value: block_size_row_maj=64 (optimal value for Nvidia GPUs).

  • use_atomic (bool) – If True, increases the number of thread blocks that can run in parallel with the help of atomic operations (custom atomic add on floats). Set this to False, if the atomic operations are running slow on your device (Nvidia GPUs before Kepler, some AMD APUs, some Intel GPUs). Default value: use_atomic=True.

  • steps_per_thread (int) – If use_atomic is set to True, this parameters defines the maximum number of loop steps performed by the parallel threads in a single thread block. Default value: steps_per_thread=64 (optimal for Nvidia GPUs).

  • verbose (bool) – Verbose output, defaults to verbose=False.

>>> with SartOpencl(geometry_matrix, block_size=1024) as invert_sart:
>>>     solution, residual = invert_sart(measurement_vector)
>>> ### or ###
>>> inv_sart = SartOpencl(geometry_matrix, block_size=1024)
>>> solution, residual = inv_sart(measurement_vector)
>>> inv_sart.clean()
__call__(measurement_vector, initial_guess=None, max_iterations=250, relaxation=1.0, beta_laplace=0.01, conv_tol=0.0001, time_limit=None)

Performs the inversion for a given measurement vector.

Parameters:

  • measurement_vector (np.ndarray) – The measured power/radiance vector with shape \((N_d)\).

  • initial_guess – An optional initial guess, can be an array of shape \((N_s)\) or a constant value that will be used to seed the algorithm.

  • max_iterations (int) – The maximum number of iterations to run the SART algorithm before returning a result, defaults to max_iterations=250.

  • relaxation (float) – The relaxation hyperparameter, defaults to relaxation=1. Consult the reference papers for more information on this hyperparameter.

  • beta_laplace (float) – The regularisation hyperparameter in the range [0, 1]. Defaults to beta_laplace=0.01.

  • conv_tol (float) – The convergence limit at which the algorithm will be terminated, unless the maximum number of iterations has been reached. The convergence is calculated as the normalised squared difference between the measurement and solution vectors. Note that reaching convergence lower than 1.e-6 is hardly possible on GPUs due to single precision calculations and relaxed math.

  • time_limit (float) – If set, the iterations will stop after this time limit (in seconds) is reached. Default value: time_limit=None.

Returns:

A tuple with the inverted solution vector \(\mathbf{x}\) as an ndarray with shape \((N_s)\), and the list of convergence values achieved after each iteration step.

clean()

Releases GPU buffers

update_laplacian_matrix(laplacian_matrix)

Updates the Laplacian matrix in GPU memory

Parameters:

laplacian_matrix (np.ndarray) – The laplacian regularisation matrix of shape \((N_s, N_s)\).

cherab.tools.inversions.nnls.invert_regularised_nnls(w_matrix, b_vector, alpha=0.01, tikhonov_matrix=None, **kwargs)

Solves \(\mathbf{b} = \mathbf{W} \mathbf{x}\) for the vector \(\mathbf{x}\), using Tikhonov regulariastion.

This is a thin wrapper around scipy.optimize.nnls, which modifies the arguments to include the supplied Tikhonov regularisation matrix.

The values of w_matrix, b_vector and alpha * tikhonov_matrix are notmalised by max(b_vector) before passing them to scipy.optimize.nnls().

Parameters:

  • w_matrix (np.ndarray) – The sensitivity matrix describing the coupling between the detectors and the voxels. Must be an array with shape \((N_d, N_s)\).

  • b_vector (np.ndarray) – The measured power/radiance vector with shape \((N_d)\).

  • alpha (float) – The regularisation hyperparameter \(\alpha\) which determines the regularisation strength of the tikhonov matrix.

  • tikhonov_matrix (np.ndarray) – The tikhonov regularisation matrix operator, an array with shape \((N_s, N_s)\). If None, the identity matrix is used.

  • **kwargs – Keyword arguments passed to scipy.optimize.nnls.

Returns:

(x, norm), the solution vector and the residual norm.

>>> from cherab.tools.inversions import invert_regularised_nnls
>>> x, norm = invert_regularised_nnls(w_matrix, b_vector, tikhonov_matrix=tikhonov_matrix)
cherab.tools.inversions.svd.invert_svd(w_matrix, b_vector)

Performs a Singular Value Decomposition (SVD) operation inversion.

Parameters:

  • w_matrix (np.ndarray) – The sensitivity matrix describing the coupling between the detectors and the voxels. Must be an array with shape \((N_d, N_s)\).

  • b_vector (np.ndarray) – The measured power/radiance vector with shape \((N_d)\).

Returns:

The solution vector x as an ndarray.

10.6.2. Voxels

class cherab.tools.inversions.voxels.Voxel

A Voxel base class.

Each Voxel is a Node in the scenegraph. Each Voxel type that inherits from this class defines its own geometry.

Variables:

volume (float) – The geometric volume of this voxel.

class cherab.tools.inversions.voxels.AxisymmetricVoxel

An axis-symmetric Voxel.

This Voxel is symmetric about the vertical z-axis. The cross section of the voxel can be arbitrarily defined by a polygon in the r-z plane. The type of geometric primitive used to define the geometric extent of this Voxel can be selected by the user and either of type Mesh or CSG. The two representations should approximately the same geometry but have different performance goals. The CSG representation uses lower memory and is a better choice when large numbers of Voxels will be present in a single scene. The Mesh representation is split into smaller components and better for cases where multiple importance sampling is important, such as weight matrices including reflection effects.

Parameters:

  • vertices – An Nx2 array specifying the voxel’s polygon outline in the r-z plane.

  • parent (Node) – The scenegraph to which this Voxel is attached.

  • material (Material) – The emission material of this Voxel, defaults to a UnityVolumeEmitter() for weight matrix calculations.

  • primitive_type (str) – Specifies the primitive type, can be either ‘mesh’ or ‘csg’. Defaults to the CSG representation.

Variables:

  • volume (float) – The geometric volume of this voxel.

  • cross_sectional_area (float) – The cross sectional area of the voxel in the r-z plane.

  • cross_section_centroid (Point2D) – The centroid of the voxel in the r-z plane.

emissivity_from_function(emission_function, grid_samples)

Calculate the average emissivity in the voxel.

Parameters:

  • emission_function (callable) – a function defining the emissivity in (r, ϕ, z) space

  • grid_samples (int) – the number of samples of the emissivitiy to use to calculate the average

Return float emissivity:

the average emissivity in the voxel cross section

Note that while the emissivity function is a 3D function, for Axisymmetric voxels the return value should be independent of toroidal angle ϕ.

class cherab.tools.inversions.voxels.VoxelCollection

The base class for collections of voxels.

Used for managing a collection of voxels when calculating a weight matrix for example.

Variables:

  • count (float) – The number of voxels in this collection.

  • total_volume (float) – The total volume of all voxels.

emissivities_from_function(emission_function, grid_samples=10)

Returns an array of sampled emissivities at each voxel location.

Note that the results will be nonsense if you mix an emission function and VoxelCollection with incompatible symmetries.

Parameters:

  • emission_function (Function3D) – Emission function to sample over.

  • grid_samples (int) – Number of emission samples to average over.

Return type:

np.ndarray

parent_all_voxels()

Add all voxels in this collection to the scenegraph.

set_active(item)

Set the ith voxel as an active emitter.

Parameters:

item – If item is an int, the ith voxel will be configured as an active emitter, all the others will be turned off. If item is the string ‘all’, all voxels will be active emitters.

unparent_all_voxels()

Remove all voxels in this collection from the scenegraph.

class cherab.tools.inversions.voxels.ToroidalVoxelGrid(voxel_coordinates, name='', parent=None, transform=None, active='all', primitive_type='csg')

A collection of axis-symmetric toroidal voxels.

This object manages a collection of voxels, where each voxel in the collection is an AxisymmetricVoxel object.

Parameters:

  • voxel_coordinates – An array/list of voxels, where each voxel element is defined by a list of 2D points.

  • name (str) – The name of this voxel collection.

  • parent (Node) – The parent scenegraph to which these voxels belong.

  • transform (AffineMatrix3D) – The coordinate transformation of this local coordinate system relative to the scenegraph parent, defaults to the identity transform.

  • active – Selects which voxels are active emitters in the initialised state. If active is an int, the ith voxel will be configured as an active emitter, all the others will be turned off. If active is the string ‘all’, all voxels will be active emitters.

  • primitive_type (str) – The geometry type to use for the AxisymmetricVoxel instances, can be [‘mesh’, ‘csg’]. See their documentation for more information. Defaults to primitive_type=’csg’.

plot(title=None, voxel_values=None, ax=None, vmin=None, vmax=None, cmap=None)

Plots a voxel grid.

If no voxel data values are provided, the plot is an outline of the grid in the r-z plane. If voxel values are provided, this method plots the voxel grid coloured by the voxel intensities.

Parameters:

  • title (str) – The title of the plot.

  • voxel_values (np.ndarray) – A 1D numpy array with length equal to the number of voxels in the collection.

  • ax – The matplotlib Axes object on which the plot will be made. If None, this function generates a new plot.

  • vmin (float) – The minimum value for the colour map.

  • vmax (float) – The maximum value for the colour map.

  • cmap – The matplotlib colour map to use for colouring the voxel intensities.

set_active(item)

Set the ith voxel as an active emitter.

Parameters:

item – If item is an int, the ith voxel will be configured as an active emitter, all the others will be turned off. If item is the string ‘all’, all voxels will be active emitters.

10.6.3. Ray Transfer Objects

Ray transfer objects accelerate the calculation of geometry matrices (or Ray Transfer Matrices as they were called in S. Kajita, et al. Contrib. Plasma Phys., 2016, 1-9) in the case of regular spatial grids. As in the case of Voxels, the spectral array is used to store the data for individual light sources (in this case the grid cells or their unions), however no voxels are created at all. Instead, a custom integration along the ray is implemented. Ray transfer objects allow to calculate geometry matrices for a single value of wavelength. Use RayTransferBox class for Cartesian grids and RayTransferCylinder class for cylindrical grids (3D or axisymmetrical).

Performance tips:

  • The best performance is achieved when Ray Transfer Objects are used with special pipelines and optimised materials (currently only rough metals are optimised, see the demos).

  • When the number of individual light sources and respective bins in the spectral array is higher than ~50-70 thousands, the lack of CPU cache memory becomes a serious factor affecting performance. Therefore, it is not recommended to use hyper-threading when calculating geometry matrices for a large number of light sources. It is also recommended to divide the calculation into several parts and to calculate partial geometry matrices for not more than ~50-70 thousands of light sources in a single run. Partial geometry matrices can easily be combined into one when all computations are complete.

class cherab.tools.raytransfer.raytransfer.RayTransferObject(primitive)

Basic class for ray transfer objects.

Variables:

  • voxel_map (np.ndarray) – An array containing the indices of the light sources.

  • ~.mask (np.ndarray) – A boolean mask array showing active (True) and inactive (False) gird cells.

  • parent (Node) – Scene-graph parent node.

  • transform (AffineMatrix3D) – An AffineMatrix3D defining the local co-ordinate system relative to the scene-graph parent.

  • step (float) – Integration step of volume integrator.

  • bins (int) – Number of light sources (the size of spectral array must be equal to this value).

invert_voxel_map()

Returns a list of arrays of cell indices belonging to each light source. This list is an inversion of voxel_map array.

class cherab.tools.raytransfer.raytransfer.RayTransferBox(xmax, ymax, zmax, nx, ny, nz, step=None, voxel_map=None, mask=None, parent=None, transform=None)

Ray transfer object for rectangular emitter defined on a regular 3D \((X, Y, Z)\) grid. The grid starts at (0, 0, 0). Use transform parameter to move it.

Parameters:

  • xmax (float) – Upper bound of grid in X direction (in meters).

  • ymax (float) – Upper bound of grid in Y direction (in meters).

  • zmax (float) – Upper bound of grid in Z direction (in meters).

  • nx (int) – Number of grid points in X direction.

  • ny (int) – Number of grid points in Y direction.

  • nz (int) – Number of grid points in Z direction.

  • step (float) – The step of integration along the ray (in meters), defaults to step = 0.1 * min(xmax / nx, ymax / ny, zmax / nz).

  • voxel_map (np.ndarray) – An array with shape (nx, ny, nz) containing the indices of the light sources. This array maps the cells in \((X, Y, Z)\) space to the respective voxels (light sources). The cells with identical indices in voxel_map array form a single voxel (light source). If voxel_map[ix, iy, iz] == -1, the cell with index (ix, iy, iz) will not be mapped to any light source. This parameters allows to apply a custom geometry (pixelated though) to the light sources. Default value: voxel_map=None.

  • mask (np.ndarray) – A boolean mask array with shape (nx, ny, nz). Allows to include (mask[ix, iy, iz] == True) or exclude (mask[ix, iy, iz] == False) the cells from the calculation. The ray tranfer matrix will be calculated only for those cells for which mask is True. This parameter is ignored if voxel_map is provided, defaults to mask=None (all cells are included).

  • parent (Node) – Scene-graph parent node or None (default = None).

  • transform (AffineMatrix3D) – An AffineMatrix3D defining the local co-ordinate system relative to the scene-graph parent (default = identity matrix).

>>> from raysect.optical import World, translate
>>> from cherab.tools.raytransfer import RayTransferBox
>>> world = World()
>>> rtb = RayTransferBox(xmax=1., ymax=1., zmax=1., nx=100, ny=100, nz=100)
>>> rtb.parent = world
>>> rtb.transform = translate(-0.5, -0.5, -0.5)
>>> ### cutting out a sphere of radius 0.5 ###
>>> x = np.linspace(-0.495, 0.495, 100)
>>> xsqr = x * x
>>> ### mask is a bollean array of shape (100, 100, 100) ###
>>> mask = xsqr[:, None, None] + xsqr[None, :, None] + xsqr[None, None, :] < 0.25
>>> rtb.mask = mask  # all cells outside this sphere are excluded
...
>>> camera.spectral_bins = rtb.bins
>>> # ray transfer matrix will be calculated for 600.5 nm
>>> camera.min_wavelength = 600.
>>> camera.max_wavelength = 601.
class cherab.tools.raytransfer.raytransfer.RayTransferCylinder(radius_outer, height, n_radius, n_height, radius_inner=0, n_polar=1, period=360.0, step=None, voxel_map=None, mask=None, parent=None, transform=None)

Ray transfer object for cylindrical emitter defined on a regular 3D \((R, \phi, Z)\) grid. This emitter is periodic in \(\phi\) direction. The base of the cylinder is located at Z = 0 plane. Use transform parameter to move it.

Parameters:

  • radius_outer (float) – Radius of the outer cylinder and the upper bound of grid in R direction (in meters).

  • height (float) – Height of the cylinder and the length of grid in Z direction (in meters).

  • n_radius (int) – Number of grid points in R direction.

  • n_height (int) – Number of grid points in Z direction.

  • radius_inner (float) – Radius of the inner cylinder and the lower bound of grid in R direction (in meters), defaults to radius_inner=0.

  • n_polar (int) – Number of grid points in \(\phi\) direction, defaults to n_polar=1 (axisymmetric case).

  • period (float) – A period in \(\phi\) direction (in degree), defaults to period=360.

  • step (float) – The step of integration along the ray (in meters), defaults to step = 0.1 * min((radius_outer - radius_inner)/n_radius, height/n_height).

  • voxel_map (np.ndarray) – An array with shape (n_radius, n_polar, n_height) containing the indices of the light sources. This array maps the cells in \((R, \phi, Z)\) space to the respective voxels (light sources). The cells with identical indices in voxel_map array form a single voxel (light source). If voxel_map[ir, iphi, iz] == -1, the cell with index (ir, iphi, iz) will not be mapped to any light source. This parameters allows to apply a custom geometry (pixelated though) to the light sources. Default value: voxel_map=None. Convert 2D (axisymmetric) voxel_map to 3D with voxel_map = voxel_map[:, None, :].

  • mask (np.ndarray) – A boolean mask array with shape (n_radius, n_polar, n_height). Allows to include (mask[ir, iphi, iz] == True) or exclude (mask[ir, iphi, iz] == False) the cells from the calculation. The ray tranfer matrix will be calculated only for those cells for which mask is True. This parameter is ignored if voxel_map is provided, defaults to mask=None (all cells are included). Convert 2D (axisymmetric) mask to 3D with mask = mask[:, None, :].

  • parent (Node) – Scene-graph parent node or None (default = None).

  • transform (AffineMatrix3D) – An AffineMatrix3D defining the local co-ordinate system relative to the scene-graph parent (default = identity matrix).

>>> from raysect.optical import World, translate
>>> from cherab.tools.raytransfer import RayTransferCylinder
>>> world = World()
>>> rtc = RayTransferCylinder(radius_outer=8., height=10., n_radius=400, n_height=1000,
                              radius_inner=4.)
>>> rtc.parent = world
>>> rtc.transform = translate(0, 0, -5.)
...
>>> camera.spectral_bins = rtc.bins
>>> # ray transfer matrix will be calculated for 600.5 nm
>>> camera.min_wavelength = 600.
>>> camera.max_wavelength = 601.

Emitters and integrators

The following emitters and integrators are used in ray transfer objects. Note that these emitters support other integrators as well, however high performance with other integrators is not guaranteed.

class cherab.tools.raytransfer.emitters.RayTransferEmitter

Basic class for ray transfer emitters defined on a regular 3D grid. Ray transfer emitters are used to calculate ray transfer matrices (geometry matrices) for a single value of wavelength.

Parameters:

  • grid_shape (tuple) – The shape of regular grid (the number of grid cells along each direction).

  • grid_steps (tuple) – The sizes of grid cells along each direction.

  • voxel_map (np.ndarray) – An array with shape grid_shape containing the indices of the light sources. This array maps the cells of regular grid to the respective voxels (light sources). The cells with identical indices in voxel_map array form a single voxel (light source). If voxel_map[i1, i2, i3] == -1, the cell with indices (i1, i2, i3) will not be mapped to any light source. This parameters allows to apply a custom geometry (pixelated though) to the light sources. Default value: voxel_map=None.

  • mask (np.ndarray) – A boolean mask array with shape grid_shape. Allows to include (mask[i1, i2, i3] == True) or exclude (mask[i1, i2, i3] == False) the cells from the calculation. The ray tranfer matrix will be calculated only for those cells for which mask is True. This parameter is ignored if voxel_map is provided, defaults to mask=None (all cells are included).

  • integrator (raysect.optical.material.VolumeIntegrator) – Volume integrator, defaults to integrator=NumericalVolumeIntegrator()

Variables:

  • grid_shape (tuple) – The shape of regular 3D grid.

  • grid_steps (tuple) – The sizes of grid cells along each direction.

  • voxel_map (np.ndarray) – An array containing the indices of the light sources.

  • ~.mask (np.ndarray) – A boolean mask array showing active (True) and inactive (False) gird cells.

  • bins (int) – Number of light sources (the size of spectral array must be equal to this value).

class cherab.tools.raytransfer.emitters.RayTransferIntegrator

Basic class for ray transfer integrators that calculate distances traveled by the ray through the voxels defined on a regular grid.

Parameters:

  • step (float) – Integration step (in meters), defaults to step=0.001.

  • min_samples (int) – The minimum number of samples to use over integration range, defaults to min_samples=2.

Variables:

  • step (float) – Integration step.

  • min_samples (int) – The minimum number of samples to use over integration range.

class cherab.tools.raytransfer.emitters.CartesianRayTransferEmitter

A unit emitter defined on a regular 3D \((X, Y, Z)\) grid, which can be used to calculate ray transfer matrices (geometry matrices). Note that for performance reason there are no boundary checks in emission_function(), or in CartesianRayTranferIntegrator, so this emitter must be placed inside a bounding box.

Parameters:

  • grid_shape (tuple) – The shape of regular \((X, Y, Z)\) grid. The number of points in X, Y and Z directions.

  • grid_steps (tuple) – The sizes of grid cells in X, Y and Z directions (in meters).

  • voxel_map (np.ndarray) – An array with shape grid_shape containing the indices of the light sources. This array maps the cells in \((X, Y, Z)\) space to the respective voxels (light sources). The cells with identical indices in voxel_map array form a single voxel (light source). If voxel_map[ix, iy, iz] == -1, the cell with indices (ix, iy, iz) will not be mapped to any light source. This parameters allows to apply a custom geometry (pixelated though) to the light sources. Default value: voxel_map=None.

  • mask (np.ndarray) – A boolean mask array with shape grid_shape. Allows to include (mask[ix, iy, iz] == True) or exclude (mask[ix, iy, iz] == False) the cells from the calculation. The ray tranfer matrix will be calculated only for those cells for which mask is True. This parameter is ignored if voxel_map is provided, defaults to mask=None (all cells are included).

  • integrator (raysect.optical.material.VolumeIntegrator) – Volume integrator, defaults to integrator=CartesianRayTransferIntegrator(step=0.1 * min(grid_steps))

Variables:

  • dx (float) – The size of grid cell in X direction (equals to grid_shape[0]).

  • dy (float) – The size of grid cell in Y direction (equals to grid_shape[1]).

  • dz (float) –

    The size of grid cell in Z direction (equals to grid_shape[2]).

    >>> from raysect.optical import World, translate, Point3D
>>> from raysect.primitive import Box
>>> from cherab.tools.raytransfer import CartesianRayTransferEmitter
>>> world = World()
>>> grid_shape = (10, 10, 10)
>>> grid_steps = (0.5, 0.5, 0.5)
>>> material = CartesianRayTransferEmitter(grid_shape, grid_steps)
>>> eps = 1.e-6  # ray must never leave the grid when passing through the volume
>>> upper = Point3D(grid_shape[0] * grid_steps[0] - eps,
                    grid_shape[1] * grid_steps[1] - eps,
                    grid_shape[2] * grid_steps[2] - eps)
>>> bounding_box = Box(lower=Point3D(0, 0, 0), upper=upper, material=material,
                       parent=world)
>>> bounding_box.transform = translate(-2.5, -2.5, -2.5)
...
>>> camera.spectral_bins = material.bins
>>> # ray transfer matrix will be calculated for 600.5 nm
>>> camera.min_wavelength = 600.
>>> camera.max_wavelength = 601.

class cherab.tools.raytransfer.emitters.CartesianRayTransferIntegrator

Calculates the distances traveled by the ray through the voxels defined on a regular grid in Cartesian coordinate system: \((X, Y, Z)\). This integrator is used with the CartesianRayTransferEmitter material to calculate ray transfer matrices (geometry matrices). The value for each voxel is stored in respective bin of the spectral array. The distances traveled by the ray through the voxel is calculated approximately and the accuracy depends on the integration step.

class cherab.tools.raytransfer.emitters.CylindricalRayTransferEmitter

A unit emitter defined on a regular 3D \((R, \phi, Z)\) grid, which can be used to calculate ray transfer matrices (geometry matrices) for a single value of wavelength. This emitter is periodic in \(\phi\) direction. Note that for performance reason there are no boundary checks in emission_function(), or in CylindricalRayTranferIntegrator, so this emitter must be placed between a couple of coaxial cylinders that act like a bounding box.

Parameters:

  • grid_shape (tuple) – The shape of regular \((R, \phi, Z)\) 3D grid. If grid_shape[1] = 1, the emitter is axisymmetric.

  • grid_steps (tuple) – The sizes of grid cells in R, \(\phi\) and Z directions. The size in \(\phi\) must be provided in degrees (sizes in R and Z are provided in meters). The period in \(\phi\) direction is defined as grid_shape[1] * grid_steps[1]. Note that the period must be a multiple of 360.

  • voxel_map (np.ndarray) – An array with shape grid_shape containing the indices of the light sources. This array maps the cells in \((R, \phi, Z)\) space to the respective voxels (light sources). The cells with identical indices in voxel_map array form a single voxel (light source). If voxel_map[ir, iphi, iz] == -1, the cell with indices (ir, iphi, iz) will not be mapped to any light source. This parameters allows to apply a custom geometry (pixelated though) to the light sources. Default value: voxel_map=None.

  • mask (np.ndarray) – A boolean mask array with shape grid_shape. Allows to include (mask[ir, iphi, iz] == True) or exclude (mask[ir, iphi, iz] == False) the cells from the calculation. The ray tranfer matrix will be calculated only for those cells for which mask is True. This parameter is ignored if voxel_map is provided, defaults to mask=None (all cells are included).

  • integrator (raysect.optical.material.VolumeIntegrator) – Volume integrator, defaults to integrator=CylindricalRayTransferIntegrator(step=0.1*min(grid_shape[0], grid_shape[-1])).

  • rmin (float) – Lower bound of grid in R direction (in meters), defaults to rmin=0.

Variables:

  • period (float) – The period in \(\phi\) direction (equals to grid_shape[1] * grid_steps[1]).

  • rmin (float) – Lower bound of grid in R direction.

  • dr (float) – The size of grid cell in R direction (equals to grid_shape[0]).

  • dphi (float) – The size of grid cell in \(\phi\) direction (equals to grid_shape[1]).

  • dz (float) – The size of grid cell in Z direction (equals to grid_shape[2]).

>>> from raysect.optical import World, translate
>>> from raysect.primitive import Cylinder, Subtract
>>> from cherab.tools.raytransfer import CylindricalRayTransferEmitter
>>> world = World()
>>> grid_shape = (10, 1, 10)  # axisymmetric case
>>> grid_steps = (0.5, 360, 0.5)
>>> rmin = 2.5
>>> material = CylindricalRayTransferEmitter(grid_shape, grid_steps, rmin=rmin)
>>> eps = 1.e-6  # ray must never leave the grid when passing through the volume
>>> radius_outer = grid_shape[0] * grid_steps[0] - eps
>>> height = grid_shape[2] * grid_steps[2] - eps
>>> radius_inner = rmin + eps
>>> bounding_box = Subtract(Cylinder(radius_outer, height), Cylinder(radius_inner, height),
                            material=material, parent=world)  # bounding primitive
>>> bounding_box.transform = translate(0, 0, -2.5)
...
>>> camera.spectral_bins = material.bins
>>> # ray transfer matrix will be calculated for 600.5 nm
>>> camera.min_wavelength = 600.
>>> camera.max_wavelength = 601.
class cherab.tools.raytransfer.emitters.CylindricalRayTransferIntegrator

Calculates the distances traveled by the ray through the voxels defined on a regular grid in cylindrical coordinate system: \((R, \phi, Z)\). This integrator is used with the CylindricalRayTransferEmitter material class to calculate ray transfer matrices (geometry matrices). The value for each voxel is stored in respective bin of the spectral array. It is assumed that the emitter is periodic in \(\phi\) direction with a period equal to material.period. The distances traveled by the ray through the voxel is calculated approximately and the accuracy depends on the integration step.

Pipelines

Very simple but fast pipelines for ray transfer matrix (geometry matrix) calculation. When calculating the ray transfer matrix, the spectral array is used to store the radiance from individual unit light sources and not the actual spectrum. In this case the spectral array may contain ~ 10,000 spectral bins but the wavelengths for all of them are equal. Spectral pipelines from Raysect still can be used, but they are slower compared to ray transfer pipelines. Note that the standard error is not calculated in these pipelines, only the mean value. Dispersive rendering and adaptive sampling features are removed to improve the performance. Use spectral pipelines from Raysect if you need these features.

class cherab.tools.raytransfer.pipelines.RayTransferPipeline0D(name='RayTransferPipeline0D', kind='power')

Simple 0D pipeline for ray transfer matrix (geometry matrix) calculation.

Parameters:

  • name (str) – The name of the pipeline. Default is ‘RayTransferPipeline0D’.

  • kind (str) – The kind of the pipeline. Can be ‘power’ (default) or ‘radiance’. In the case of ‘power’, the resulting matrix is multiplied by the sensitivity of the detector, and the units of the matrix are [m^3 sr], which gives the units of power [W] for the product of the ray transfer matrix and the emission profile. In case of ‘radiance’, the sensitivity is not taken into account and the matrix is calculated in [m], which gives the units of radiance [W m^-2 sr^-1] for the product of the ray transfer matrix and the emission profile. Note that if the sensitivity of the detector is 1 (e.g. PinholeCamera, VectorCamera), the ‘power’ and ‘radiance’ give the same results.

Variables:

matrix (np.ndarray) – Ray transfer matrix, a 1D array of size \(N_{bin}\).

>>> from cherab.tools.raytransfer import RayTransferPipeline0D
>>> pipeline = RayTransferPipeline0D(kind='radiance')
class cherab.tools.raytransfer.pipelines.RayTransferPipeline1D(name='RayTransferPipeline1D', kind='power')

Simple 1D pipeline for ray transfer matrix (geometry matrix) calculation.

Parameters:

  • name (str) – The name of the pipeline. Default is ‘RayTransferPipeline0D’.

  • kind (str) – The kind of the pipeline. Can be ‘power’ (default) or ‘radiance’. In the case of ‘power’, the resulting matrix is multiplied by the sensitivity of the detector, and the units of the matrix are [m^3 sr], which gives the units of power [W] for the product of the ray transfer matrix and the emission profile. In case of ‘radiance’, the sensitivity is not taken into account and the matrix is calculated in [m], which gives the units of radiance [W m^-2 sr^-1] for the product of the ray transfer matrix and the emission profile. Note that if the sensitivity of the detector is 1 (e.g. PinholeCamera, VectorCamera), the ‘power’ and ‘radiance’ give the same results.

Variables:

matrix (np.ndarray) – Ray transfer matrix, a 2D array of shape \((N_{pixel}, N_{bin})\).

>>> from cherab.tools.raytransfer import RayTransferPipeline1D
>>> pipeline = RayTransferPipeline1D(kind='radiance')
class cherab.tools.raytransfer.pipelines.RayTransferPipeline2D(name='RayTransferPipeline2D', kind='power')

Simple 2D pipeline for ray transfer matrix (geometry matrix) calculation.

Parameters:

  • name (str) – The name of the pipeline. Default is ‘RayTransferPipeline0D’.

  • kind (str) – The kind of the pipeline. Can be ‘power’ (default) or ‘radiance’. In the case of ‘power’, the resulting matrix is multiplied by the sensitivity of the detector, and the units of the matrix are [m^3 sr], which gives the units of power [W] for the product of the ray transfer matrix and the emission profile. In case of ‘radiance’, the sensitivity is not taken into account and the matrix is calculated in [m], which gives the units of radiance [W m^-2 sr^-1] for the product of the ray transfer matrix and the emission profile. Note that if the sensitivity of the detector is 1 (e.g. PinholeCamera, VectorCamera), the ‘power’ and ‘radiance’ give the same results.

Variables:

matrix (np.ndarray) – Ray transfer matrix, a 3D array of shape \((N_x, N_y, N_{bin})\).

>>> from cherab.tools.raytransfer import RayTransferPipeline2D
>>> pipeline = RayTransferPipeline2D(kind='radiance')