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', trim_rows=True, trim_columns=True, 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).
trim_rows (bool) – Whether to automatically remove all-zero matrix rows, i.e. rows corresponding to pixels which don’t see any of the cells in the inversion grid.
trim_columns (bool) – Whether to automatically remove all-zero matrix columns, i.e. columns corresponding to grid cells not seen by any pixel. This will also remove these cells from the grid object.
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.

Note

.npz is the recommended file 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.

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.

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.

unformat_image(im_vector, coords='Native', fill_value=nan)

Formats a given a 1D data vector of image pixel values (i.e. \(b\) in \(Ax = b\)) back in to a 2D image array.

Parameters

im_vector (numpy.ndarray or sparse matrix) – 1D vector of image data, must have length equal to the number of rows in the geometry matrix.
coords (str) – What image orientation to return the image: ‘Original’, ‘Display’ or ‘Native’ (whichever was used when creating the geometry matrix).
fill_value (float) – Value to return in any image pixels which are not included in the geometry matrix.

Returns

2D array containing the image.

Return type

numpy.ndarray

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 (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.
rmax – 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.
zmin – 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.
zmax – 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.gm.PoloidalVolumeGrid

calcam.gm.trigrid(wall_contour, max_cell_scale, rmin=None, rmax=None, zmin=None, zmax=None, triopts=[])

Create a reconstruction grid with triangular grid cells conforming to the wall contour, using Jonathan Richard Shewchuk’s “triangle” triangulation library.

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 (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.
rmax – 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.
zmin – 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.
zmax – 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.
trioppts (list of str) – List of string options to pass to triangle.triangulate(), see triangle trignale documentation page for possible options.

Returns

Generated grid.

Return type

calcam.gm.PoloidalVolumeGrid

calcam.gm.solps_grid(solps_file, wall_contour, rmin=None, rmax=None, zmin=None, zmax=None, cell_attrs=False)

Load a grid from a SOLPS-ITER DivGeo-Carre file.

Parameters

solps_file (str) – Filename of DivGeo-Carre .geo file containing grid coordinates to load.
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.
rmin (float) – 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 SOLPS-ITER grid is loaded.
rmax (float) – 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 SOLPS-ITER grid is loaded.
zmin (float) – 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 SOLPS-ITER grid is loaded.
zmax (float) – 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 SOLPS-ITER grid is loaded.
cell_attrs (bool) – Whether to return the cell centre coordinates and cell indices as provided by SOLPS-ITER.

Returns

The reconstruction grid.

if cell_attrs is set to true, also returns an Nx8 NumPy array where N is the number of grid cells. The 8 columns of this array contain for each cell: Poloidal index, Radial index, R centre, Z centre, and indices in to the resulting grid object’s vertex list of the 4 vertices of the cell.

Return type

calcam.gm.PoloidalVolumeGrid