Tomography Geometry Matrices

Background

It is common to want to recover the distribution of light emission from the plasma in the poloidal plane, assuming toroidal symmetry, from cameras viewing approximately tangentially. The emission is reconstructed on a grid defined the poloidal \((R,Z)\) plane. The relationship between the brightness of emission from each grid element and the signal at each camera pixel is then a large system of linear equations, where the brightness of a pixel is given by a weighted sum of the brightness of each grid element along that pixel’s line of sight. This is most compactly expressed by a matrix multiplication \(\mathbf{Ax = b}\) where \(\mathbf{A}\) is the geometry matrix, \(\mathbf{x}\) is a vector of the brightnesses of each grid element (which we want to find), and \(\mathbf{b}\) is a vector of the pixel brightnesses. The calcam.gm module provides tools for generating the geometry matrix \(\mathbf{A}\).

For examples of using the features documented on this page, see the Examples page.

Model assumptions and limitations

In calculation of the geometry matrices, camera sight-lines are assumed to be infinitely thin pencil beams, i.e. finite etendue and depth-of-field effects are not included. If ray casting with binning = 1, each image pixel is characterised by a single pencil beam sight-line at the image centre, i.e. the finite size of the pixels is not accounted for. The values of matrix element \(i,j\) is given by the length, in metres, of the \(i^{th}\) sight line which passes through the \(j^{th}\) grid cell.

The Geometry Matrix class

class calcam.gm.GeometryMatrix(grid, raydata, pixel_order='C', calc_status_callback=calcam_status_printer)

Class to represent a scalar geometry matrix and associated metadata.

A new geometry matrix can be created by instantiating this class with the below parameters, or a saved geometry matrix can be loaded from disk using the fromfile() class method.

The matrix itself can be accessed in the data attribute, where it is stored as a sparse matrix using the scipy.sprase.csr_matrix class.

Parameters
  • grid (calcam.gm.PoloidalVolumeGrid) – Reconstruction grid to use

  • raydata (calcam.RayData) – Ray data for the camera to be inverted

  • pixel_order (str) – What pixel order to use when flattening the 2D image array in to the 1D data vector. Default ‘C’ goes left-to-right, then row-by-row (NumPy default), alternatively ‘F’ goes top-to-bottom, then column-by-column (MATLAB default).

  • calc_status_callback (callable) – Callable which takes a single argument, which will be called with status updates about the calculation. It will be called with either a string for textual status updates or a float from 0 to 1 specifying the progress of the calculation. By default, status updates are printed to stdout. If set to None, no status updates are issued.

data

The geometry matrix data itself.

Type

scipy.sparse.csr_matrix

format_image(image, coords=None)

Format a given 2D camera image in to a 1D data vector (i.e. \(b\) in \(Ax = b\)) appropriate for use with this geometry matrix. This will bin the image, remove any excluded pixels and reshape it to a 1D vector.

Parameters
  • image (numpy.ndarray) – Input image.

  • coords (str) – Either ‘Display’ or ‘Original’, specifies what orientation the input image is in. If not givwn, it will be auto-detected if possible.

Returns

1xN_pixels image data vector. Note that this is returned as a sparse matrix object despite its density being 100%; this is for consistency with the matrix itself.

Return type

scipy.sparse.csr_matrix

classmethod fromfile(filename)

Load a saved geometry matrix from disk.

Parameters

filename (str) – File name to load from. Can be a NumPy (.npz), MATLAB (.mat) or zipped ASCII (.zip) file.

Returns

Loaded geometry matrix.

Return type

calcam.GeometryMatrix

get_included_pixels()

Get a mask showing which image pixels are included in the geometry matrix.

Returns

Boolean array the same shape as the camera image after binning, where True indicates the corresponding pixel is included and False indicates the pixel is excluded.

Return type

numpy.ndarray

get_los_coverage()

Get the number of lines of sight viewing each grid element.

Returns

Vector with as many elements as there grid cells, with values being how many sight-lines interact with that grid element.

Return type

numpy.ndarray

grid

The inversion grid associated with the geometry matrix.

Type

calcam.gm.PoloidalVolumeGrid

save(filename)

Save the geometry matrix to a file.

Parameters

filename (str) – File name to save to, including file extension. The file extension determines the format to be saved: ‘.npz’ for compressed NumPy binary format, ‘.mat’ for MATLAB format or ‘.zip’ for Zipped collection of ASCII files.

Note: .npz is the recommended format; .mat is provided if compatibility with MATLAB is required but produces larger file sizes, and .zip is provided to make maximally compatible data files but is extremely slow to save and load.

set_binning(binning)

Set the level of image binning. Can be used to decrease the size of the matrix to reduce memory or computation requirements for inversions.

Parameters

binning (float) – Desired image binning. Must be larger than the existing binning value.

set_included_pixels(pixel_mask, coords=None)

Set which image pixels should be included, or not. Can be used to exclude image pixels which are known to have bad data or otherwise do not conform to the assumptions of the inversion.

Note: excluding pixels is a non-reversible process since their matrix rows will be removed. It is therefore recommended to keep a master copy of the matrix with all pixels included and then use this function on a transient copy of the matrix.

Parameters
  • pixel_mask (numpy.ndarray) – Boolean array the same shape as the un-binned camera image, where True or 1 indicates a pixel to be included and False or 0 represents a pixel to be excluded.

  • coords (str) – Either ‘Display’ or ‘Original’, specifies what orientation the input pixel mask is in. If not givwn, it will be auto-detected if possible.

Reconstruction grids

class calcam.gm.PoloidalVolumeGrid(vertices, cells, wall_contour=None, src=None)

Class for representing tomographic reconstruction grids with polygonal grid cells in the R, Z plane. Quantities are assumed to be uniform within the volume of each grid cell. Grid cells can be arbitrary polygons but all cells in the grid must all have the same number of sides.

Grids can be constructed by directly instantiating this class with the following parameters, or alternatively convenience functions such as calcam.gm.squaregrid() are provided for easily creating various types of grid.

Parameters
  • vertices (numpy.ndarray) – (N_verts x 2) array of floats containing the (R,Z) coordinates of the grid cell vertices.

  • cells (numpy.ndarray) – (N_cells x N_verts_per_cell) array of integers specifying which vertices (indexes in to the vertices array) define each grid cell. For each cell (array row), the vertices must be listed in order around the cell perimeter (in either direction).

  • wall_contour (numpy.ndaray) – Nx2 array containing the R,Z wall contour of the machine. If provided, this is used for plotting purposes.

  • src (str) – Human readable string describing where the grid came from, for data provenance purposes.

property extent

A 4-element tuple of floats containing the R,Z extent of the grid in metres (Rmin,Rmax,Zmin,Zmax)

get_cell_intersections(ray_start, ray_end, plot=False)

Get the intersections of a ray, i.e. a straight line in 3D space, with the grid cell boundaries.

Parameters
  • ray_start (sequence) – 3-element sequence containing the X,Y,Z coordinates of the ray’s start position.

  • ray_end (sequence) – 3-element sequence containing the X,Y,Z coordinates of the ray’s end position.

  • plot (bool) – If set to True, the function will also plot the ray in R,Z and red circles at each intersection location. Mainly intended for debugging and algorithm demonstration.

Returns

Intersection information:

  • NumPy array containing the lengths along the ray where intersections with grid cell boundaries lines were found (in order from start to end of the ray).

  • list of lists containing the indices of the grid cell(s) which were involved in each intersection. Each intersection will be associated with 1 or more grid cells.

Return type

tuple

interpolate(data, r_new, z_new, fill_value=nan)

Given a vector of data on the grid, return the data values at positions (r_new,z_new). Since quantities are assumed uniform within each grid cell, this is a “nearest neighbour” type interpolation where the value returned at new each point is the value of the grid cell that point lies within.

Note: this is not particularly optimised for speed, but at least provides the functionality for now.

Parameters
  • data (np.ndarray) – 1D array containing the data defined on the grid; must have the same number of elements as there are grid cells.

  • r_new (np.ndarray) – Arry of new R positions; must be the same shape as z_new

  • z_new (np.ndarray) – Array of new Z positions, must be the same shape as r_new

  • fill_value (float) – What value to return for (r_new,z_new) points which lie outside the grid. Default is np.nan.

Returns

Array the same shape as r_new and z_new containing the data values at positions (r_new,z_new).

Return type

np.ndarray

property n_cells

The number of cells in the grid.

property n_segments

The number of line segments (cell sides) in the grid.

property n_vertices

The number of vertices in the grid.

plot(data=None, clim=None, cmap=None, line_colour=0, 0, 0, cell_linewidth=None, cblabel='', axes=None)

Either plot a given data vector on the grid, or if no data vector is given, plot the grid itself.

Parameters
  • data (array) – Data to plot on the grid. Must be a 1D array with as many elements as there are grid cells. If not given, only the grid structure is plotted.

  • clim (sequence or None) – 2 element sequence giving the colourmap limits for the data [min,max]. If set to None, the data min and max are used.

  • cmap (str or None) – Name of the matplotlib colourmap to use to display the data. If not given, matplotlib’s default colourmap will be used.

  • line_colour (sequence or None) – 3-element sequence of values between 0 and 1 specifying the R,G,B colour with which to draw the wall contour and grid cell boundaries. If set to None, the wall and cell boundaries are not drawn.

  • cell_linewidth (float) – Line width to use to show the grid cell boundaries. If set to 0, the grid cell boundaries will not be drawn.

  • cblabel (str) – Label for the data colour bar, if plotting data.

  • axes (matplotlib.pyplot.Axes) – Matplotlib axes on which to plot. If not given, a new figure will be created.

Returns

MatPlotLib objects comprising the plot:

  • matplotlib.axes.Axes : The matplotlib axes containing the plot.

  • matplotlib.collections.PatchCollection : PatchCollection containing the patches used to show the data and grid cells. This object has useful methods for further adjusting the plot like set_cmap(), set_clim() etc.

  • list of matplotlib.lines.Line2D : List of matpltolib line objects making up the wall contour.

Return type

tuple

remove_cells(cell_inds)

Remove grid cells with the given indices from the grid.

Parameters

cell_inds (sequence) – Sequence of integers specifying which cell indices to remove.

Creating commonly used grid types

The following functions are provided for convenience to easily generate some commonly used types of grid.

calcam.gm.squaregrid(wall_contour, cell_size, rmin=None, rmax=None, zmin=None, zmax=None)

Create a reconstruction grid with square grid cells.

Parameters
  • wall_contour (str or numpy.ndarray) – Either the name of a Calcam CAD model from which to use the wall contour, or an N x 2 array of R,Z points defining the machine wall.

  • cell_size (float) – Side length of each grid cell in metres.

  • rmin, rmax, zmin, zmax (floats) – Optional limits of the grid extent in the R, Z plane. Any combination of these may or may not be given; if none are given the entire wall contour interior is gridded.

Returns

Generated grid.

Return type

calcam.PoloidalVolumeGrid

calcam.gm.trigrid(wall_contour, max_cell_scale, rmin=None, rmax=None, zmin=None, zmax=None, **kwargs)

Create a reconstruction grid with triangular grid cells conforming to the wall contour. Requires the MeshPy package (https://pypi.org/project/MeshPy/) to be installed, since it uses J. Shewchuk’s “triangle” via the MeshPy to generate the grid.

Parameters
  • wall_contour (str or numpy.ndarray) – Either the name of a Calcam CAD model from which to use the wall contour, or an N x 2 array of R,Z points defining the machine wall.

  • max_cell_scale (float) – Approximate maximum allowed length scale for a grid cell, in metres.

  • rmin, rmax, zmin, zmax (floats) – Optional limits of the grid extent in the R, Z plane. Any combination of these may or may not be given; if none are given the entire wall contour interior is gridded.

Any additional keyword arguments will be passed directly to the triangle mesher, meshpy.triangle.build(). This can be used to further control the meshing; see the MeshPy documentation for available arguments.

Returns

Generated grid.

Return type

calcam.PoloidalVolumeGrid