XDesign

XDesign is an open-source Python package for generating configurable x-ray imaging phantoms, simulating data acquisition, and benchmarking x-ray tomographic image reconstruction.

Goals

• Assist faster development of new generation tomographic reconstruction methods

• Allow quantitative comparison of different reconstruction methods

• Create a framework for designing x-ray imaging experiments

Current Scope

• Customizable 2D phantoms constructed from circles and convex polygons

• Quantitative reconstruction quality and probe coverage metrics

• Attenuation interactions with X-ray probes of uniform flux

• Use of analytic (exact) solutions for algorithms and computation

Contribute

• Issue Tracker: https://github.com/tomography/xdesign/issues

• Documentation: https://github.com/tomography/xdesign/tree/master/docs

• Source Code: https://github.com/tomography/xdesign/tree/master/xdesign

• Tests: https://github.com/tomography/xdesign/tree/master/tests

License

The project is licensed under the BSD-3 license.

Install

Since version 0.5, XDesign is available on the conda-forge channel. Install it in the usual way:

$ conda install xdesign -c conda-forge

API Documentation

XDesign aims to provide tools for designing xray-imaging experiments.

These tools include a computational geometry library for simulating x-ray phantoms and data acquisition, quality metrics for quantitatively rating the image reconstructions and scanning procedures, and other helpful functions for specifying motion and coding aperatures for computational imaging.

xdesign.acquisition

Define objects and methods for simulated data acquisition.

This not only includes physical things like Probe, detectors, turntables, and lenses, but non-physical things such as scanning patterns.

Module author: Doga Gursoy <dgursoy@aps.anl.gov>

Module author: Daniel J Ching <carterbox@users.noreply.github.com>

Classes:

Probe([p1, p2, size, intensity, energy])

A square cross-section x-ray beam for probing Phantoms.

Functions:

raster_scan2D(sa, st[, meta])

Return the coordinates of a 2D raster scan.

class xdesign.acquisition.Probe(p1=None, p2=None, size=0.0, intensity=1.0, energy=15.0)

Bases: xdesign.geometry.line.Line

A square cross-section x-ray beam for probing Phantoms.

p1, p2

Deprecated since version 0.4: Measure now uses theta, h, v coordinates instead.

Type

xdesign.geometry.Point

size

The size of probe in centimeters.

Type

float, cm (default: 0.0 cm)

intensity

The intensity of the beam in candela.

Type

float, cd (default: 1.0 cd)

energy

The energy of the probe in eV.

Type

float, eV (default: 15 eV)

.. todo::

Implement additional attributes for Probe such as wavelength, etc.

property cross_section

Return the cross-sectional area of a square beam.

distance(other)

Return the closest distance between entities.

half_space()

Return the half space polytope respresentation of the probe.

intersect(polygon)

Return the intersection with polygon.

measure(phantom, theta, h, perc=None)

Measure the phantom from the given position.

Parameters

theta, h – The coordinates of the Probe.

xdesign.acquisition.raster_scan2D(sa, st, meta=False)

Return the coordinates of a 2D raster scan.

Parameters

• sa (int) – The number of projeciton angles in [0, 2PI).

• st (int) – The number of Probe steps at each projection angle. [-0.5, 0.5)

• nmeta (int >= 0) – The number of meta steps. Meta steps are the offset from the starting Probe position after each rotation.

Returns

theta, h, v (np.array (M,)) – Probe positions for scan

xdesign.codes

Generate codes for space- and time-coded apertures.

Module author: Daniel Ching

Code Generating Functions:

mura_1d(L)	Return the longest MURA whose length is less than or equal to L.
mura_2d(M[, N])	Return the largest 2D MURA whose lengths are less than M and N.
raskar(npool)	Return the coded mask from Raskar et al.

xdesign.codes.mura_1d(L)

Return the longest MURA whose length is less than or equal to L.

From Wikipedia: A Modified uniformly redundant array (MURA) can be generated in any length L that is prime and of the form:

L = 4m + 1, m = 1, 2, 3, ...,

the first six such values being L = 5, 13, 17, 29, 37. The binary sequence of a linear MURA is given by A[0:L] where:

A[i] = {
    0 if i = 0,
    1 if i is a quadratic residue modulo L, i != 0,
    0 otherwise,
}

xdesign.codes.mura_2d(M, N=None)

Return the largest 2D MURA whose lengths are less than M and N.

From Wikipedia: A rectangular MURA, A[0:M, 0:N], is defined as follows:

A[i, j] = {
    0 if i = 0,
    1 if j = 0, i != 0,
    1 if C[i] * C[j] = 1,
    0 othewise,
}

C[i] = {
    1 if i is a quadratic residue modulo p,
    -1 otherwise,
}

where p is the length of the matching side M, N.

xdesign.codes.raskar(npool)

Return the coded mask from Raskar et al.

xdesign.constants

Constants in cgs units.

xdesign.constants.AVOGADRO_NUMBER

Avagadro constant [1/mol]

Type

float

xdesign.constants.BOLTZMANN_CONSTANT

Boltzmann constant [erg/k]

Type

float

xdesign.constants.CLASSICAL_ELECTRON_RADIUS

Classical electron radius [cm]

Type

float

xdesign.constants.ELECTRONIC_CHARGE

Electronic charge [esu]

Type

float

xdesign.constants.ELECTRON_VOLT

Electron volt (keV) [erg]

Type

float

xdesign.constants.ELECTRON_MASS

Electron mass [g]

Type

float

xdesign.constants.FINE_STRUCTURE_CONSTANT

Fine structure constant

Type

float

xdesign.constants.PLANCK_CONSTANT

Reduced planck’s constant [keV*s]

Type

float

xdesign.constants.PROTON_MASS

Proton mass [g]

Type

float

xdesign.constants.SPEED_OF_LIGHT

Speed of light in vacuum [cm/s]

Type

float

xdesign.constants.THOMPSON_CROSS_SECTION

Thomson cross section [cm^2]

Type

float

xdesign.constants.PI

Ratio of a circle’s circumference to its diameter

Type

float

xdesign.constants.wavelength(energy)

Return wavelength [cm] of light given energy [keV].

xdesign.formats

Functions:

xdesign.geometry

Define geometric objects for compuational geometry operations.

Module author: Doga Gursoy <dgursoy@aps.anl.gov>

Module author: Daniel J Ching <carterbox@users.noreply.github.com>

Zero dimensional entities

Define zero dimensional geometric entities.

Classes:

Point(x)

Define a point in ND cartesian space.

class xdesign.geometry.point.Point(x)

Bases: xdesign.geometry.entity.Entity

Define a point in ND cartesian space.

Parameters

x (ndarray, list) – ND coordinates of the point.

Raises

TypeError – If x is not a list or ndarray.

collision(other)

Return True if this Point collides with another entity.

contains(other)

Return wether the other is within the bounds of the Point.

Points can only contain other Points.

distance(other)

Return the closest distance this and the other.

property norm

Euclidian (L2) norm of the vector to the point.

rotate(theta, point=None, axis=None)

Rotates theta radians around an axis.

scale(vector)

Scale the ambient space in each dimension according to vector.

Scaling is centered on the origin.

translate(vector)

Translate the point along the given vector.

property x

Dimension 0 of the point.

property y

Dimension 1 of the point.

property z

Dimension 2 of the point.

xdesign.geometry.point.calc_standard(A)

Return the standard equation for the row-wise points in A.

The coefficents (c_{0}*x + … = c_{1}) describe the hyper-plane defined by the row-wise N-dimensional points A.

Parameters

A (np.array (…, N, N)) – Each row is an N-dimensional point on the plane.

Returns

• c0 (np.array (…, N)) – The first N coefficients for the hyper-plane

• c1 (np.array (…, 1)) – The last coefficient for the hyper-plane

xdesign.geometry.point.dim(self)

Return the dimensionality of the ambient space.

xdesign.geometry.point.distance(self, other)

Return the closest distance this and the other.

xdesign.geometry.point.norm(self)

Euclidian (L2) norm of the vector.

xdesign.geometry.point.rotated(self, theta, center=None, axis=None)

Rotates theta radians around an axis.

One dimensional entities

Define one dimensional geometric entities.

Classes:

Line(p1, p2)	Line in 2D cartesian space.
Segment(p1, p2)	Defines a finite line segment from two unique points.

class xdesign.geometry.line.Line(p1, p2)

Bases: xdesign.geometry.line.LinearEntity

Line in 2D cartesian space.

The constructor takes two unique Point.

p1

Type

Point

p2

Type

Point

distance(other)

Returns the closest distance between entities.

intercept(n)

Calculates the intercept for the nth dimension.

property standard

Returns coeffients for the first N-1 standard equation coefficients. The Nth is returned separately.

property xintercept

Return the x-intercept.

property yintercept

Return the y-intercept.

class xdesign.geometry.line.Segment(p1, p2)

Bases: xdesign.geometry.line.Line

Defines a finite line segment from two unique points.

distance(other)

Return the distance to the other.

property midpoint

Return the midpoint of the line segment.

Two dimensional entities

Define two dimensional geometric entities.

Classes:

Curve(center)	The base class for closed manifolds defined by a single equation.
Circle(center, radius[, sign])	Circle in 2D cartesian space.
Polygon(vertices[, sign])	A convex polygon in 2D cartesian space.
RegularPolygon(center, radius, order[, ...])	A regular polygon in 2D cartesian space.
Triangle(p1, p2, p3)	Triangle in 2D cartesian space.
Rectangle(center, side_lengths)	Rectangle in 2D cartesian space.
Square(center[, side_length, radius])	Square in 2D cartesian space.
Mesh([obj, faces])	A collection of Entities

class xdesign.geometry.area.Circle(center, radius, sign=1)

Bases: xdesign.geometry.area.Curve

Circle in 2D cartesian space.

center

The center point of the circle.

Type

Point

radius

The radius of the circle.

Type

scalar

sign

The sign of the area

Type

int (-1 or 1)

property area

Return the area.

property bounding_box

Return the axis-aligned bounding box as two numpy vectors.

property circumference

Returns the circumference.

contains(other)

Return whether other is a proper subset.

Return one boolean for all geometric entities. Return an array of boolean for array input.

property diameter

Returns the diameter.

property list

Return list representation for saving to files.

property patch

Returns a matplotlib patch.

class xdesign.geometry.area.Curve(center)

Bases: xdesign.geometry.entity.Entity

The base class for closed manifolds defined by a single equation. e.g. Circle, Sphere, or Torus.

center

Type

Point

rotate(theta, point=None, axis=None)

Rotates the Curve by theta radians around an axis which passes through a point radians.

translate(vector)

Translates the Curve along a vector.

class xdesign.geometry.area.Mesh(obj=None, faces=[])

Bases: xdesign.geometry.entity.Entity

A collection of Entities

faces

A list of the Entities

Type

list

area

The total area of the Entities

Type

float

population

The number entities in the Mesh

Type

int

radius

The radius of a bounding circle

Type

float

append(t)

Add a triangle to the mesh.

property bounding_box

Return the axis-aligned bounding box as two numpy vectors.

property center

contains(other)

Return whether this Mesh contains other.

FOR ALL x, THERE EXISTS a face of the Mesh that contains x AND (ALL cut outs that contain x or THERE DOES NOT EXIST a cut out).

import_triangle(obj)

Loads mesh data from a Python Triangle dict.

property patch

pop(i=- 1)

Pop i-th triangle from the mesh.

rotate(theta, point=None, axis=None)

Rotate entity around an axis which passes through a point by theta radians.

scale(vector)

Scale entity.

translate(vector)

Translate entity.

class xdesign.geometry.area.Polygon(vertices, sign=1)

Bases: xdesign.geometry.entity.Entity

A convex polygon in 2D cartesian space.

It is defined by a number of distinct vertices of class Point. Superclasses include Square, Triangle, etc.

vertices

Type

List of Points

sign

The sign of the area

Type

int (-1 or 1)

:raises ValueError : If the number of vertices is less than three.:

area = <MagicMock name='mock()' id='140346947267856'>

property bounding_box

Return the axis-aligned bounding box as two numpy vectors.

property bounds

Returns a 4-tuple (xmin, ymin, xmax, ymax) representing the bounding rectangle for the Polygon.

center = <MagicMock name='mock()' id='140346947267856'>

contains(other)

Return whether this Polygon contains the other.

property edges

Return a list of lines connecting the points of the Polygon.

half_space = <MagicMock name='mock()' id='140346947267856'>

property list

Return list representation.

property numpy

Return Numpy representation.

property numverts

property patch

Returns a matplotlib patch.

perimeter = <MagicMock name='mock()' id='140346947267856'>

radius = <MagicMock name='mock()' id='140346947267856'>

rotate(theta, point=None, axis=None)

Rotates the Polygon around an axis which passes through a point by theta radians.

translate(vector)

Translates the polygon by a vector.

class xdesign.geometry.area.Rectangle(center, side_lengths)

Bases: xdesign.geometry.area.Polygon

Rectangle in 2D cartesian space.

Defined by a point and a vector to enforce perpendicular sides.

Parameters

side_lengths (array) – The lengths of the sides

area = <MagicMock name='mock()' id='140346947267856'>

class xdesign.geometry.area.RegularPolygon(center, radius, order, angle=0, sign=1)

Bases: xdesign.geometry.area.Polygon

A regular polygon in 2D cartesian space.

It is defined by the polynomial center, order, and radius.

By default (i.e. when the angle parameter is zero), the regular polygon is oriented so that one of the vertices is at coordinates \((x + r, x)\) where \(x\) is the x-coordinate of center and \(r\) = radius. The angle parameter is only meaningful modulo \(2\pi /\) order since rotation by \(2\pi /\) order gives a result equivalent to no rotation.

Parameters

• center (Point) – The center of the polygon

• radius (float) – Distance from polygon center to vertices

• order (int) – Order of the polygon (e.g. order 6 is a hexagon).

• angle (float) – Optional rotation angle in radians.

• sign (int (-1 or 1)) – Optional sign of the area (see Polygon)

class xdesign.geometry.area.Square(center, side_length=None, radius=None)

Bases: xdesign.geometry.area.Rectangle

Square in 2D cartesian space.

Defined by a point and a length to enforce perpendicular sides.

class xdesign.geometry.area.Triangle(p1, p2, p3)

Bases: xdesign.geometry.area.Polygon

Triangle in 2D cartesian space.

It is defined by three distinct points.

area = <MagicMock name='mock()' id='140346947267856'>

center = <MagicMock name='mock()' id='140346947267856'>

Intersect

Define algorithms to support intersection calculation.

Functions:

clip_SH(clipEdges, polygon)	Clip a polygon using the Sutherland-Hodgeman algorithm.
halfspacecirc(d, r)	Return the area of intersection between a circle and half-plane.

xdesign.geometry.intersect.clip_SH(clipEdges, polygon)

Clip a polygon using the Sutherland-Hodgeman algorithm.

Parameters

• clipEdges [[A, b], …] – half-spaces defined by coefficients

• polygon

xdesign.geometry.intersect.halfspacecirc(d, r)

Return the area of intersection between a circle and half-plane.

Returns the smaller fraction of a circle split by a line d units away from the center of the circle.

Parameters

• d (scalar) – The distance from the line to the center of the circle

• r (scalar) – The radius of the circle

Returns

f (scalar) – The proportion of the circle in the smaller half-space

References

Glassner, A. S. (Ed.). (2013). Graphics gems. Elsevier.

xdesign.material

Defines objects which auto-generate a parameterized Phantom.

Module author: Daniel J Ching <carterbox@users.noreply.github.com>

Module author: Doga Gursoy <dgursoy@aps.anl.gov>

Classes:

SimpleMaterial([mass_attenuation])	Simple material with constant mass_attenuation parameter only.
XraylibMaterial(compound, density)	Materials which use xraylib data to automatically calculate material properties based on beam energy in KeV.

class xdesign.material.SimpleMaterial(mass_attenuation=1.0)

Bases: xdesign.material.Material

Simple material with constant mass_attenuation parameter only.

density

The mass density of the material

Type

float [g/cm^3] (default: 1.0)

linear_attenuation(energy)

linear x-ray attenuation [1/cm] for the energy [KeV].

mass_attenuation(energy)

mass x-ray attenuation [1/cm] for the energy [KeV].

class xdesign.material.XraylibMaterial(compound, density)

Bases: xdesign.material.Material

Materials which use xraylib data to automatically calculate material properties based on beam energy in KeV.

compound

Molecular formula of the material.

Type

string

density

The mass density of the material

Type

float [g/cm^3] (default: 1.0)

beta()

x.__getitem__(y) <==> x[y]

delta()

x.__getitem__(y) <==> x[y]

xdesign.metrics

Objects and methods for computing the quality of reconstructions.

Module author: Daniel J Ching <carterbox@users.noreply.github.com>

Coverage metrics

Objects and methods for computing coverage based quality metrics.

These methods are based on the scanning trajectory only.

Module author: Daniel J Ching <carterbox@users.noreply.github.com>

Functions:

coverage_approx(gmin, gsize, ngrid, ...[, ...])

Approximate procedure coverage with a Riemann sum.

xdesign.metrics.coverage.coverage_approx(gmin, gsize, ngrid, probe_size, theta, h, v, weights=None, anisotropy=1, num_rays=16)

Approximate procedure coverage with a Riemann sum.

The intersection between the beam and each pixel is approximated by using a Reimann sum of n rectangles: width beam.size / n and length dist where dist is the length of segment of the line alpha which passes through the pixel parallel to the beam.

If anisotropy is True, then coverage_map.shape is (M, N, 2, 2), where the two extra dimensions contain coverage anisotopy information as a second order tensor.

Parameters

• procedure (Probe generator) – A generator which defines a scanning procedure by returning a sequence of Probe objects.

• region (np.array [cm]) – A rectangle in which to map the coverage. Specify the bounds as [[min_x, max_x], [min_y, max_y]]. i.e. column vectors pointing to the min and max corner.

• pixel_size (float [cm]) – The edge length of the pixels in the coverage map in centimeters.

• n (int) – The number of lines per beam

• anisotropy (bool) – Whether the coverage map includes anisotropy information

Returns

coverage_map (numpy.ndarray) – A discretized map of the Probe coverage.

See also

plot.plot_coverage_anisotropy()

Full-reference metrics

Defines full-referene image quality metricsself.

These methods require a ground truth in order to make a quality assessment.

Module author: Daniel J Ching <carterbox@users.noreply.github.com>

Functions:

pcc(A, B[, masks])	Return the Pearson product-moment correlation coefficients (PCC).
ssim(img1, img2[, sigma, L, K, scale, ...])	Return the Structural SIMilarity index (SSIM) of two images.
msssim(img0, img1[, nlevels, sigma, L, K, ...])	Multi-Scale Structural SIMilarity index (MS-SSIM).

class xdesign.metrics.fullref.ImageQuality(original, reconstruction)

Bases: object

Store information about image quality.

img0

Type

array

img1

Stacks of reference and deformed images.

Type

array

metrics

A dictionary with image quality information organized by scale. metric[scale] = (mean_quality, quality_map)

Type

dict

method

The metric used to calculate the quality

Type

string

quality(method='MSSSIM', L=1.0, **kwargs)

Compute the full-reference image quality of each image pair.

Available methods include SSIM [8], MSSSIM [9], VIFp [7]

Parameters

• method (string, optional, (default: MSSSIM)) – The quality metric desired for this comparison. Options include: SSIM, MSSSIM, VIFp

• L (scalar) – The dynamic range of the data. This value is 1 for float representations and 2^bitdepth for integer representations.

xdesign.metrics.fullref.msssim(img0, img1, nlevels=5, sigma=1.2, L=1.0, K=(0.01, 0.03), alpha=4, beta_gamma=None)

Multi-Scale Structural SIMilarity index (MS-SSIM).

Parameters

• img0 (array)

• img1 (array) – Two images for comparison.

• nlevels (int) – The max number of levels to analyze

• sigma (float) – Sets the standard deviation of the gaussian filter. This setting determines the minimum scale at which quality is assessed.

• L (scalar) – The dynamic range of the data. This value is 1 for float representations and 2^bitdepth for integer representations.

• K (2-tuple) – A list of two constants which help prevent division by zero.

• alpha (float) – The exponent which weights the contribution of the luminance term.

• beta_gamma (list) – The exponent which weights the contribution of the contrast and structure terms at each level.

Returns

metrics (dict) – A dictionary with image quality information organized by scale. metric[scale] = (mean_quality, quality_map) The valid range for SSIM is [-1, 1].

References

Multi-scale Structural Similarity Index (MS-SSIM) Z. Wang, E. P. Simoncelli and A. C. Bovik, “Multi-scale structural similarity for image quality assessment,” Invited Paper, IEEE Asilomar Conference on Signals, Systems and Computers, Nov. 2003

xdesign.metrics.fullref.pcc(A, B, masks=None)

Return the Pearson product-moment correlation coefficients (PCC).

Parameters

• A, B (ndarray) – The two images to be compared

• masks (list of ndarrays, optional) – If supplied, the data under each mask is computed separately.

Returns

covariances (array, list of arrays)

xdesign.metrics.fullref.ssim(img1, img2, sigma=1.2, L=1, K=(0.01, 0.03), scale=None, alpha=4, beta_gamma=4)

Return the Structural SIMilarity index (SSIM) of two images.

A modified version of the Structural SIMilarity index (SSIM) based on an implementation by Helder C. R. de Oliveira, based on the implementation by Antoine Vacavant, ISIT lab, antoine.vacavant@iut.u-clermont1.fr http://isit.u-clermont1.fr/~anvacava

xdesign.metrics.fullref.img1

Type

array

xdesign.metrics.fullref.img2

Two images for comparison.

Type

array

xdesign.metrics.fullref.sigma

Sets the standard deviation of the gaussian filter. This setting determines the minimum scale at which quality is assessed.

Type

float

xdesign.metrics.fullref.L

The dynamic range of the data. The difference between the minimum and maximum of the data: 2^bitdepth for integer representations.

Type

scalar

xdesign.metrics.fullref.K

A list of two constants which help prevent division by zero.

Type

2-tuple

xdesign.metrics.fullref.alpha

The exponent which weights the contribution of the luminance term.

Type

float

xdesign.metrics.fullref.beta_gamma

The exponent which weights the contribution of the contrast and structure terms at each level.

Type

list

Returns

metrics (dict) – A dictionary with image quality information organized by scale. metric[scale] = (mean_quality, quality_map) The valid range for SSIM is [-1, 1].

References

Z. Wang, A. C. Bovik, H. R. Sheikh and E. P. Simoncelli. Image quality assessment: From error visibility to structural similarity. IEEE Transactions on Image Processing, 13(4):600–612, 2004.

Z. Wang and A. C. Bovik. Mean squared error: Love it or leave it? - A new look at signal fidelity measures. IEEE Signal Processing Magazine, 26(1):98–117, 2009.

Silvestre-Blanes, J., & Pérez-Lloréns, R. (2011, September). SSIM and their dynamic range for image quality assessment. In ELMAR, 2011 Proceedings (pp. 93-96). IEEE.

Standards-based metrics

Defines standards based image quality metrics.

These methods require the reconstructed image to be of a specifically shaped standard object such as a siemens star or a zone plate.

Module author: Daniel J Ching <carterbox@users.noreply.github.com>

Functions:

compute_mtf(phantom, image)	Approximate the modulation tranfer function using the HyperbolicCocentric phantom.
compute_mtf_ffst(phantom, image[, Ntheta])	Calculate the MTF using the method described in [2].
compute_mtf_lwkj(phantom, image)	Calculate the MTF using the modulated Siemens Star method in [6].
compute_nps_ffst(phantom, A[, B, plot_type])	Calculate the noise power spectrum from a unit circle image using the method from [2].
compute_neq_d(phantom, A, B)	Calculate the NEQ according to recommendations by [1].

xdesign.metrics.standards.compute_mtf(phantom, image)

Approximate the modulation tranfer function using the HyperbolicCocentric phantom. Calculate the MTF from the modulation depth at each edge on the line from (0.5,0.5) to (0.5,1). MTF = (hi-lo)/(hi+lo)

Deprecated since version 0.3: This method rapidly becomes inaccurate at small wavelenths because the measurement gets out of phase with the waves due to rounding error. Use another one of the MTF functions instead. This method will be removed in xdesign 0.6.

See also

compute_mtf_ffst() compute_mtf_lwkj()

Parameters

• phantom (HyperbolicConcentric) – Predefined phantom of cocentric rings whose widths decay parabolically.

• image (ndarray) – The reconstruction of the above phantom.

Returns

• wavelength (list) – wavelenth in the scale of the original phantom

• MTF (list) – MTF values

xdesign.metrics.standards.compute_mtf_ffst(phantom, image, Ntheta=4)

Calculate the MTF using the method described in [2].

See also

compute_mtf_lwkj()

Parameters

• phantom (UnitCircle) – Predefined phantom with single circle whose radius is less than 0.5.

• image (ndarray) – The reconstruction of the phantom.

• Ntheta (scalar) – The number of directions at which to calculate the MTF.

Returns

• wavenumber (ndarray) – wavelenth in the scale of the original phantom

• MTF (ndarray) – MTF values

• bin_centers (ndarray) – the center of the bins if Ntheta >= 1

xdesign.metrics.standards.compute_mtf_lwkj(phantom, image)

Calculate the MTF using the modulated Siemens Star method in [6].

See also

compute_mtf_ffst()

Parameters

• phantom (SiemensStar)

• image (ndarray) – The reconstruciton of the SiemensStar

Returns

• frequency (array) – The spatial frequency in cycles per unit length

• M (array) – The MTF values for each frequency

xdesign.metrics.standards.compute_neq_d(phantom, A, B)

Calculate the NEQ according to recommendations by [1].

Parameters

• phantom (UnitCircle) – The unit circle class with radius less than 0.5

• A (ndarray) – The reconstruction of the above phantom.

• B (ndarray) – The reconstruction of the above phantom with different noise. This second reconstruction enables allows use of trend subtraction instead of zero mean normalization.

Returns

• mu_b – The spatial frequencies

• NEQ – the Noise Equivalent Quanta

xdesign.metrics.standards.compute_nps_ffst(phantom, A, B=None, plot_type='frequency')

Calculate the noise power spectrum from a unit circle image using the method from [2].

Parameters

• phantom (UnitCircle) – The unit circle phantom.

• A (ndarray) – The reconstruction of the above phantom.

• B (ndarray) – The reconstruction of the above phantom with different noise. This second reconstruction enables allows use of trend subtraction instead of zero mean normalization.

• plot_type (string) – ‘histogram’ returns a plot binned by radial coordinate wavenumber ‘frequency’ returns a wavenumber vs wavenumber plot

Returns

• bins – Bins for the radially binned NPS

• counts – NPS values for the radially binned NPS

• X, Y – Frequencies for the 2D frequency plot NPS

• NPS (2Darray) – the NPS for the 2D frequency plot

xdesign.phantom

Objects and methods for computing the quality of reconstructions.

Module author: Daniel J Ching <carterbox@users.noreply.github.com>

Phantoms

Defines an object for simulating X-ray phantoms.

Module author: Daniel J Ching <carterbox@users.noreply.github.com>

Phantom([geometry, children, material])	An object for the purpose of evaluating X-ray imaging methods.
save_phantom(phantom, filename)	Save phantom to file as a python repr.
load_phantom(filename)	Load phantom from file containing a python repr.
pickle_phantom(phantom, filename)	Save phantom to file as a python pickle.
unpickle_phantom(filename)	Load phantom from file as a python pickle.

class xdesign.phantom.phantom.Phantom(geometry=None, children=[], material=None)

Bases: object

An object for the purpose of evaluating X-ray imaging methods.

Phantoms may be hierarchical structures with children that are contained within and/or a parent which contains them. They have two parts: a geometry and properties. The geometry defines the spatial extent over which the properties are valid. Properties are parameters which a Probe uses to measure the Phantom.

All Phantoms must fit within the geometry of their ancestors. Phantoms whose geometry is None act as containers.

geometry

The spatial boundary of the Phantom; may be None.

Type

Entity

children

A list of Phantoms contained in this Phantom.

parent

The Phantom containing this Phantom.

material

The mass_attenuation of the phantom.

population

The number of decendents of this phantom.

append(child)

Add a child to the Phantom.

Only add the child if it is contained within the geometry of its ancestors.

property center

Return the centroid of the Phantom.

property density

Return the geometric density of the Phantom.

property geometry

Return the geometry of the Phantom.

property is_leaf

Return whether the Phantom is a leaf node.

pop(i=- 1)

Pop the i-th child from the Phantom.

property radius

Return the radius of the smallest boundary sphere.

rotate(theta, point=Point([]), axis=None)

Rotate around an axis that passes through the given point.

sprinkle(counts, radius, gap=0, region=None, material=None, max_density=1, shape=<class 'xdesign.geometry.area.Circle'>)

Sprinkle a number of Circle shaped Phantoms around the Phantom. Uses various termination criteria to determine when to stop trying to add circles.

Parameters

• counts (int) – The number of circles to be added.

• radius (scalar or list) – The radius of the circles to be added.

• gap (float, optional) – The minimum distance between circle boundaries. A negative value allows overlapping edges.

• region (Entity, optional) – The new circles are confined to this shape. None if the circles are allowed anywhere.

• max_density (scalar, optional) – Stops adding circles when the geometric density of the phantom reaches this ratio.

• material (scalar, optional) – A mass attenuation parameter passed to the circles.

Returns

counts (scalar) – The number of circles successfully added.

translate(vector)

Translate the Phantom.

property volume

Return the volume of the Phantom

xdesign.phantom.phantom.load_phantom(filename)

Load phantom from file containing a python repr.

xdesign.phantom.phantom.pickle_phantom(phantom, filename)

Save phantom to file as a python pickle.

xdesign.phantom.phantom.save_phantom(phantom, filename)

Save phantom to file as a python repr.

xdesign.phantom.phantom.unpickle_phantom(filename)

Load phantom from file as a python pickle.

Standard phantoms

Defines an object for simulating X-ray phantoms.

Module author: Daniel J Ching <carterbox@users.noreply.github.com>

XDesignDefault()	Generates a Phantom for internal testing of XDesign.
HyperbolicConcentric([min_width, exponent])	Generates a series of cocentric alternating black and white circles whose radii are changing at a parabolic rate.
DynamicRange([steps, jitter, geometry])	Generates a phantom of randomly placed circles for determining dynamic range.
DogaCircles([n_sizes, size_ratio, n_shuffles])	Rows of increasingly smaller circles.
SlantedSquares([count, angle, gap])	Generates a collection of slanted squares.
UnitCircle([radius, material])	Generates a phantom with a single circle in its center.
SiemensStar([n_sectors, center, radius])	Generates a Siemens star.

class xdesign.phantom.standards.DogaCircles(n_sizes=5, size_ratio=0.5, n_shuffles=5)

Bases: xdesign.phantom.phantom.Phantom

Rows of increasingly smaller circles. Initally arranged in an ordered Latin square, the inital arrangement can be randomly shuffled.

radii

radii of circles

Type

ndarray

x

x position of circles

Type

ndarray

y

y position of circles

Type

ndarray

class xdesign.phantom.standards.DynamicRange(steps=10, jitter=True, geometry=Rectangle(<MagicMock name='mock()' id='140346947267856'>, <MagicMock name='mock().tolist()' id='140346817705296'>))

Bases: xdesign.phantom.phantom.Phantom

Generates a phantom of randomly placed circles for determining dynamic range.

Parameters

• steps (scalar, optional) – The orders of magnitude (base 2) that the colors of the circles cover.

• jitter (bool, optional) – True : circles are placed in a jittered grid False : circles are randomly placed

• shape (string, optional)

class xdesign.phantom.standards.HyperbolicConcentric(min_width=0.1, exponent=0.5)

Bases: xdesign.phantom.phantom.Phantom

Generates a series of cocentric alternating black and white circles whose radii are changing at a parabolic rate. These line spacings cover a range of scales and can be used to estimate the Modulation Transfer Function. The radii change according to this function: r(n) = r0*(n+1)^k.

radii

The list of radii of the circles

Type

list

widths

The list of the widths of the bands

Type

list

class xdesign.phantom.standards.SiemensStar(n_sectors=4, center=Point([]), radius=0.5)

Bases: xdesign.phantom.phantom.Phantom

Generates a Siemens star.

ratio

The spatial frequency times the proportional radius. e.g to get the frequency, f, divide this ratio by some fraction of the maximum radius: f = ratio/radius_fraction

Type

scalar

class xdesign.phantom.standards.SlantedSquares(count=10, angle=0.08726646259972222, gap=0)

Bases: xdesign.phantom.phantom.Phantom

Generates a collection of slanted squares. Squares are arranged in concentric circles such that the space between squares is at least gap. The size of the squares is adaptive such that they all remain within the unit circle.

angle

the angle of slant in radians

Type

scalar

count

the total number of squares

Type

scalar

gap

the minimum space between squares

Type

scalar

side_length

the size of the squares

Type

scalar

squares_per_level

the number of squares at each level

Type

list

radius_per_level

the radius at each level

Type

list

n_levels

the number of levels

Type

scalar

class xdesign.phantom.standards.UnitCircle(radius=0.5, material=SimpleMaterial(mass_attenuation=1.0))

Bases: xdesign.phantom.phantom.Phantom

Generates a phantom with a single circle in its center.

class xdesign.phantom.standards.XDesignDefault

Bases: xdesign.phantom.phantom.Phantom

Generates a Phantom for internal testing of XDesign.

The default phantom is: (1) nested, it contains phantoms within phantoms; (2) geometrically simple, the sinogram can be verified visually; and (3) representative, it contains the three main geometric elements: circle, polygon, and mesh.

Custom phantoms

Defines an object for simulating X-ray phantoms.

Module author: Daniel J Ching <carterbox@users.noreply.github.com>

Soil([porosity])	Generates a phantom with structure similar to soil.
WetCircles()
Foam([size_range, gap, porosity])	Generates a phantom with structure similar to foam.
Softwood()	Generate a Phantom with structure similar to wood.

class xdesign.phantom.custom.Foam(size_range=[0.05, 0.01], gap=0, porosity=1)

Bases: xdesign.phantom.standards.UnitCircle

Generates a phantom with structure similar to foam.

class xdesign.phantom.custom.Softwood

Bases: xdesign.phantom.phantom.Phantom

Generate a Phantom with structure similar to wood.

Parameters

• ringsize (float [cm]) – The thickness of the annual rings in cm.

• latewood_fraction (float) – The volume ratio of latewood cells to earlywood cells

• ray_fraction (float) – The ratio of rows of ray cells to rows of tracheids

• ray_height (float [cm]) – The height of the ray cells

• cell_width, cell_height (float [cm]) – The shape of the earlywood cells

• cell_thickness (float [cm]) – The thickness of the earlywood cell walls

• frame (arraylike [cm]) – A bounding box for the cells

class xdesign.phantom.custom.Soil(porosity=0.412)

Bases: xdesign.phantom.standards.UnitCircle

Generates a phantom with structure similar to soil.

References

Schlüter, S., Sheppard, A., Brown, K., & Wildenschild, D. (2014). Image processing of multiphase images obtained via X‐ray microtomography: a review. Water Resources Research, 50(4), 3615-3639.

class xdesign.phantom.custom.WetCircles

Bases: xdesign.phantom.standards.UnitCircle

xdesign.plot

Contains functions for visualizing Phantom and ImageQuality metrics.

DEFAULT_COLOR_MAPmatplotlib.colors.Colormap

The color map used to choose property colors.

DEFAULT_COLORmatplotlib.colors

The face color of geometry.

POLY_COLORmatplotlib.colors

The face color of polygons.

DEFAULT_EDGE_COLORmatplotlib.colors

The color of geometry edges.

POLY_EDGE_COLORmatplotlib.colors

The color of polygon edges.

LABEL_COLORmatplotlib.colors

The color of number labels on phantom plots.

POLY_LINEWIDTHfloat

The edge width for polygons. See matplotlib.patches.Patch.set_linewidth().

CURVE_LINEWIDTHfloat

The edge width for curves. See matplotlib.patches.Patch.set_linewidth().

PLOT_STYLES :

A list of 126 unique line styles.

Module author: Daniel J Ching <carterbox@users.noreply.github.com>

Classes:

Phantom Plotting Functions:

sidebyside(p[, size, labels, prop, figsize, dpi])	Displays the geometry and the discrete property function of the given Phantom side by side.
discrete_phantom(phantom, size[, ratio, ...])	Return a discrete map of the property in the phantom.
plot_phantom(phantom[, axis, labels, ...])	Plot a Phantom to the given axis.
plot_mesh(mesh[, axis, alpha, c])	Plot a Mesh to the given axis.
plot_polygon(polygon[, axis, alpha, c])	Plot a Polygon to the given axis.
plot_curve(curve[, axis, alpha, c])	Plot a Curve to the given axis.

Metrics Plotting Functions:

plot_coverage_anisotropy(coverage_map, **kwargs)	Plot the coverage anisotropy using 2D glyphs.
plot_metrics(imqual)	Plot full reference metrics of ImageQuality data.
plot_mtf(faxis, MTF[, labels])	Plot the MTF.
plot_nps(X, Y, NPS)	Plot the 2D frequency plot for the NPS.
plot_neq(freq, NEQ)	Plot the NEQ.
get_pie_glyphs(xy, values[, color, trace_normal])	Returns a list of pie glyphs at coordinates xy representing values

xdesign.plot.combine_grid(Amin, A, Bmin, B)

Add grid B to grid A by aligning min corners and clipping B

Parameters

• Amin, Bmin (int tuple) – The coordinates of the minimum corner of A and B

• A, B (numpy.ndarray) – The two arrays to add to each other

Returns

AB (numpy.ndarray) – The combined grid

Raises

ValueError – If A and B are do not have the same number of dimensions

xdesign.plot.discrete_geometry(geometry, psize, ratio=9)

Draw the geometry onto a patch the size of its bounding box.

Parameters

• geometry (geometry.Entity) – A geometric object with dim, bounding_box, and contains methods

• psize (float [cm]) – The real size of the pixels in the discrete image

• ratio (int (default: 9)) – The supersampling ratio for antialiasing. 1 means no antialiasing

Returns

• corner (1darray [cm]) – The min corner of the patch

• patch (ndarray) – The discretized geometry in it’s bounding box

Raises

ValueError – If ratio is less than 1 or psize is less than or equal to 0.

xdesign.plot.discrete_phantom(phantom, size, ratio=9, uniform=True, prop='linear_attenuation')

Return a discrete map of the property in the phantom.

The values of overlapping phantom.Phantom are additive.

Parameters

• phantom (phantom.Phantom)

• size (scalar) – The side length in pixels of the resulting 1 by 1 cm image.

• ratio (scalar, optional (default: 9)) – The antialiasing works by supersampling. This parameter controls how many pixels in the larger representation are averaged for the final representation. e.g. if ratio = 9, then the final pixel values are the average of 81 pixels.

• uniform (boolean, optional (default: True)) – When set to False, changes the way pixels are averaged from a uniform weights to gaussian weigths.

• prop (str, optional (default: linear_attenuation)) – The name of the property to discretize

Returns

image (numpy.ndarray) – The discrete representation of the Phantom that is size x size. 0 if phantom has no geometry or material property.

Raises

ValueError – If size is less than or equal to 0

xdesign.plot.get_pie_glyphs(xy, values, color='coverage', trace_normal=1, **kwargs)

Returns a list of pie glyphs at coordinates xy representing values

The areas of the pie sectors are proportional to the elements of the vector that the glyph represents. The default color of the Glyph is determined by the sum of the values divided by the trace_normal and the plot.DEFAULT_COLOR_MAP.

Parameters

• xy ((M, 2) float) – Locations of glyph centers

• values ((M, N) float) – Bin sizes for each glyph

• color (string) –

  coverage

  The color is determined by the sum of the values.

  standard deviation

  The color is standard deviation of the bins

  Kullback-Leibler

  The color is the Kullback-Leibler devation

  random

  The color is randomly assigned from the DEFAULT_COLOR_MAP.

• trace_normal (float) – A scalar used to normalize the trace for coloring the glyph.

• kwargs (dict) – Arguments passed to the patches constructor

See also

plot_coverage_anisotropy()

xdesign.plot.multiroll(x, shift, axis=None)

Roll an array along each axis.

Parameters

• x (array_like) – Array to be rolled.

• shift (sequence of int) – Number of indices by which to shift each axis.

• axis (sequence of int, optional) – The axes to be rolled. If not given, all axes is assumed, and len(shift) must equal the number of dimensions of x.

Returns

y (numpy array, with the same type and size as x) – The rolled array.

Notes

The length of x along each axis must be positive. The function does not handle arrays that have axes with length 0.

See also

numpy.roll()

Example

Here’s a two-dimensional array:

>>> x = np.arange(20).reshape(4,5)
>>> x
array([[ 0,  1,  2,  3,  4],
       [ 5,  6,  7,  8,  9],
       [10, 11, 12, 13, 14],
       [15, 16, 17, 18, 19]])

Roll the first axis one step and the second axis three steps:

>>> multiroll(x, [1, 3])
array([[17, 18, 19, 15, 16],
       [ 2,  3,  4,  0,  1],
       [ 7,  8,  9,  5,  6],
       [12, 13, 14, 10, 11]])

That’s equivalent to:

>>> np.roll(np.roll(x, 1, axis=0), 3, axis=1)
array([[17, 18, 19, 15, 16],
       [ 2,  3,  4,  0,  1],
       [ 7,  8,  9,  5,  6],
       [12, 13, 14, 10, 11]])

Not all the axes must be rolled. The following uses the axis argument to roll just the second axis:

>>> multiroll(x, [2], axis=[1])
array([[ 3,  4,  0,  1,  2],
       [ 8,  9,  5,  6,  7],
       [13, 14, 10, 11, 12],
       [18, 19, 15, 16, 17]])

which is equivalent to:

>>> np.roll(x, 2, axis=1)
array([[ 3,  4,  0,  1,  2],
       [ 8,  9,  5,  6,  7],
       [13, 14, 10, 11, 12],
       [18, 19, 15, 16, 17]])

References

Warren Weckesser

xdesign.plot.plot_angle_intensity(angle, intensity, background_color=(0.3, 0.3, 0.3))

Plot the phase angle and intensity in the same plot using 2D glyphs.

The glyphs are 120 degree pie glyphs whose orientation and hue are determined by the angle and whose brightness are determined by the intensity.

xdesign.plot.plot_coverage_anisotropy(coverage_map, **kwargs)

Plot the coverage anisotropy using 2D glyphs.

Parameters

kwargs – Keyword arguments for the Glyphs.

See also

metrics.coverage_approx(), Glyph

xdesign.plot.plot_curve(curve, axis=None, alpha=None, c=None)

Plot a Curve to the given axis.

Parameters

• curve (Curve) – A Curve to plot on the given axis.

• axis (matplotlib.axis.Axis, optional) – The axis where the Curve should be plotted. None creates a new axis.

• alpha (float, optional) – The plot opaqueness. 0 is transparent. 1 is opaque.

• c (matplotlib.colors, optional) – The color of the plotted curve.

xdesign.plot.plot_geometry(geometry, axis=None, alpha=None, c=None, z=0.0, t=0.0001)

Plot a Entity on the given axis.

Parameters

• geometry (Entity) – A geometry to plot on the given axis.

• axis (matplotlib.axis.Axis, optional) – The axis where the geometry should be plotted. None creates a new axis.

• alpha (float, optional) – The plot opaqueness. 0 is transparent. 1 is opaque.

• c (matplotlib.colors, optional) – The color of the plotted geometry.

xdesign.plot.plot_mesh(mesh, axis=None, alpha=None, c=None)

Plot a Mesh to the given axis.

Parameters

• mesh (Mesh) – A Mesh to plot on the given axis.

• axis (matplotlib.axis.Axis, optional) – The axis where the Mesh should be plotted. None creates a new axis.

• alpha (float, optional) – The plot opaqueness. 0 is transparent. 1 is opaque.

• c (matplotlib.colors, optional) – The color of the plotted Mesh.

xdesign.plot.plot_metrics(imqual)

Plot full reference metrics of ImageQuality data.

Parameters

imqual (ImageQuality) – The data to plot.

References

Colors taken from this gist

xdesign.plot.plot_mtf(faxis, MTF, labels=None)

Plot the MTF. Return the figure reference.

xdesign.plot.plot_neq(freq, NEQ)

Plot the NEQ. Return the figure reference.

xdesign.plot.plot_nps(X, Y, NPS)

Plot the 2D frequency plot for the NPS. Return the figure reference.

xdesign.plot.plot_phantom(phantom, axis=None, labels=None, c_props=[], c_map=None, i=- 1, z=0.0, t=0.0001)

Plot a Phantom to the given axis.

Parameters

• phantom (Phantom) – A phantom to be plotted.

• axis (matplotlib.axis.Axis) – The axis where the phantom should be plotted. None creates a new axis.

• labels (bool, optional) – True : Each Phantom given a unique number.

• c_props (list of str, optional) – List of Phantom properties to use for colormapping the geometries. [] colors the geometries by type.

• c_map (function, optional) – A function which takes the list of prop(s) for a Phantom as input and returns a matplolib color specifier. [5]

xdesign.plot.plot_polygon(polygon, axis=None, alpha=None, c=None)

Plot a Polygon to the given axis.

Parameters

• polygon (Polygon) – A Polygon to plot on the given axis.

• axis (matplotlib.axis.Axis, optional) – The axis where the Polygon should be plotted. None creates a new axis.

• alpha (float, optional) – The plot opaqueness. 0 is transparent. 1 is opaque.

• c (matplotlib.colors, optional) – The color of the plotted Polygon.

xdesign.plot.sidebyside(p, size=100, labels=None, prop='mass_attenuation', figsize=(6, 3), dpi=100, **kwargs)

Displays the geometry and the discrete property function of the given Phantom side by side.

xdesign.recon

Defines methods for reconstructing data from the acquisition module.

The algorithm module contains methods for reconstructing tomographic data including gridrec, SIRT, ART, and MLEM. These methods can be used as benchmarks for custom reconstruction methods or as an easy way to access reconstruction algorithms for developing other methods such as noise correction.

Note

Using tomopy is recommended instead of these functions for heavy computation.

Module author: Doga Gursoy <dgursoy@aps.anl.gov>

Functions:

art(gmin, gsize, data, theta, h, init[, ...])	Reconstruct data using ART algorithm.
sirt(gmin, gsize, data, theta, h, init[, ...])	Reconstruct data using SIRT algorithm.
mlem(gmin, gsize, data, theta, h, init[, niter])	Reconstruct data using MLEM algorithm.
update_progress(progress)	Draw a process bar in the terminal.

xdesign.recon.art(gmin, gsize, data, theta, h, init, niter=10, weights=None, save_interval=None)

Reconstruct data using ART algorithm. [4].

xdesign.recon.mlem(gmin, gsize, data, theta, h, init, niter=10)

Reconstruct data using MLEM algorithm.

xdesign.recon.sirt(gmin, gsize, data, theta, h, init, niter=10, weights=None, save_interval=None)

Reconstruct data using SIRT algorithm. [3].

xdesign.recon.update_progress(progress)

Draw a process bar in the terminal.

Parameters

process (float) – The percentage completed e.g. 0.10 for 10%

Examples

This section contains Jupyter Notebooks examples for various XDesign functions.

To run these examples in a notebook install Jupyter and run the notebooks from their source

This page was generated from /home/docs/checkouts/readthedocs.org/user_builds/xdesign/checkouts/0.5.x/docs/source/demos/Shepp.ipynb. Interactive online version:

Simple How-to Explaining Phantoms

Demonstrate simple basic custom phantom and sinogram generation with XDesign.

import matplotlib.pyplot as plt
import numpy as np
from xdesign import *

###################|###################|###################|###################|

Phantom creation

Create various Phantoms each with unique geometry. Make non-convex polygons by meshing together convex polygons.

# Make a circle with a triangle cut out
m = Mesh()
m.append(Circle(Point([0.0, 0.0]), radius=0.5))
m.append(-Triangle(Point([-0.3, -0.2]),
                   Point([0.0, -0.3]),
                   Point([0.3, -0.2])))


head = Phantom(geometry=m)

# Make two eyes separately
eyeL = Phantom(geometry=Circle(Point([-0.2, 0.0]), radius=0.1))
eyeR = Phantom(geometry=Circle(Point([0.2, 0.0]), radius=0.1))

Define materials to use in the phantom. Assigning multiple phantoms the same material saves memory.

material = SimpleMaterial(mass_attenuation=1.0)

head.material = material
eyeL.material = material
eyeR.material = material

Collect the phantoms together by making the eyes and mouth children of the head Phantom.

head.append(eyeL)
head.append(eyeR)

print(repr(head))

Phantom(geometry=Mesh(faces=[Circle(center=Point([0.0, 0.0]), radius=0.5, sign=1), Triangle(Point([-0.3, -0.2]), Point([0.0, -0.3]), Point([0.3, -0.2]))]), children=[Phantom(geometry=Circle(center=Point([-0.2, 0.0]), radius=0.1, sign=1), children=[], material=SimpleMaterial(mass_attenuation=1.0)), Phantom(geometry=Circle(center=Point([0.2, 0.0]), radius=0.1, sign=1), children=[], material=SimpleMaterial(mass_attenuation=1.0))], material=SimpleMaterial(mass_attenuation=1.0))

/home/user/python/venvs/py373/lib/python3.7/site-packages/xdesign/geometry/area.py:724: UserWarning: Didn't check that Mesh contains Circle.
  warnings.warn("Didn't check that Mesh contains Circle.")

Viewing phantom geometry and properties

Plot the Phantom geometry and properties with a colorbar.

fig = plt.figure(figsize=(7, 3), dpi=100)

# plot geometry
axis = fig.add_subplot(121, aspect='equal')
plt.grid()
plot_phantom(head, axis=axis, labels=False)
plt.xlim([-.5, .5])
plt.ylim([-.5, .5])

# plot property
plt.subplot(1, 2, 2)
im = plt.imshow(discrete_phantom(head, 100, prop='mass_attenuation'),
                interpolation='none', cmap=plt.cm.inferno, origin='lower')

# plot colorbar
fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.16, 0.05, 0.7])
fig.colorbar(im, cax=cbar_ax)

# save the figure
plt.savefig('Shepp_sidebyside.png', dpi=600,
        orientation='landscape', papertype=None, format=None,
        transparent=True, bbox_inches='tight', pad_inches=0.0,
        frameon=False)
plt.show()

Simulate data acquisition

Simulate data acquisition for parallel beam around 180 degrees.

NPIXEL = 100
theta, h = np.meshgrid(np.linspace(0, np.pi, NPIXEL, endpoint=False),
                       np.linspace(0, 1, NPIXEL, endpoint=False) - 0.5 + 1/NPIXEL/2)
theta = theta.flatten()
h = h.flatten()
v = h*0

probe = Probe(size=1/NPIXEL)

sino = probe.measure(head, theta, h)
sino = -np.log(sino)

/home/user/python/venvs/py373/lib/python3.7/site-packages/xdesign/geometry/algorithms.py:54: RuntimeWarning: halfspacecirc was out of bounds, -1.6769698407692601e-09
  RuntimeWarning)

Plot the sinogram.

plt.figure(figsize=(8, 8), dpi=100)
plt.imshow(np.reshape(sino, (NPIXEL, NPIXEL)), cmap='inferno',
           interpolation='nearest')
plt.savefig('Shepp_sinogram.png', dpi=600,
        orientation='landscape', papertype=None, format=None,
        transparent=True, bbox_inches='tight', pad_inches=0.0,
        frameon=False)
plt.show()

This page was generated from /home/docs/checkouts/readthedocs.org/user_builds/xdesign/checkouts/0.5.x/docs/source/demos/StandardPatterns.ipynb. Interactive online version:

Standard Test Patterns

Generates sidebyside plots of all the standard test patterns in xdesign.

from xdesign import *
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

p = SlantedSquares(count=16, angle=5/360*2*np.pi, gap=0.01)
sidebyside(p)
plt.savefig('SlantedSquares_sidebyside.png', dpi='figure',
        orientation='landscape', papertype=None, format=None,
        transparent=True, bbox_inches='tight', pad_inches=0.0,
        frameon=False)
plt.show()

h = HyperbolicConcentric()
sidebyside(h)
plt.savefig('HyperbolicConcentric_sidebyside.png', dpi='figure',
        orientation='landscape', papertype=None, format=None,
        transparent=True, bbox_inches='tight', pad_inches=0.0,
        frameon=False)
plt.show()

u = UnitCircle(radius=0.4)
sidebyside(u)
plt.savefig('UnitCircle_sidebyside.png', dpi='figure',
        orientation='landscape', papertype=None, format=None,
        transparent=True, bbox_inches='tight', pad_inches=0.0,
        frameon=False)
plt.show()

d = DynamicRange(steps=16, jitter=True)
sidebyside(d)
plt.savefig('DynamicRange_sidebyside.png', dpi='figure',
        orientation='landscape', papertype=None, format=None,
        transparent=True, bbox_inches='tight', pad_inches=0.0,
        frameon=False)
plt.show()

l = DogaCircles(n_sizes=8, size_ratio=0.5, n_shuffles=0)
l.rotate(np.pi/2, Point([0.0, 0.0]))
sidebyside(l)
plt.savefig('DogaCircles_sidebyside.png', dpi='figure',
        orientation='landscape', papertype=None, format=None,
        transparent=True, bbox_inches='tight', pad_inches=0.0,
        frameon=False)
plt.show()

s = SiemensStar(32)
sidebyside(s)
plt.savefig('SiemensStar_sidebyside.png', dpi='figure',
        orientation='landscape', papertype=None, format=None,
        transparent=True, bbox_inches='tight', pad_inches=0.0,
        frameon=False)
plt.show()

fig = plt.figure(figsize=(8, 6), dpi=200)
gs1 = gridspec.GridSpec(3, 4)
gs1.update(wspace=0.4, hspace=0.4) # set the spacing between axes.
phantoms = [l, d, u, h, p, s]
letters = ['a','b','c','d','e','f','g']
for i in range(0, len(phantoms)):

    axis = plt.subplot(gs1[2*i], aspect=1)
    plot_phantom(phantoms[i], axis=axis)
    plt.grid()
#     axis.invert_yaxis()
    axis.set_xticks(np.linspace(-0.5, 0.5, 6, True))
    axis.set_yticks(np.linspace(-0.5, 0.5, 6, True))
    plt.xlim([-.5, .5])
    plt.ylim([-.5, .5])
    plt.title('('+ letters[i] +')')

    axis = plt.subplot(gs1[2*i+1], aspect=1)
    plt.imshow(discrete_phantom(phantoms[i], 200), cmap='inferno', origin='lower')
    axis.set_xticks(np.linspace(0, 200, 6, True))
    axis.set_yticks(np.linspace(0, 200, 6, True))

plt.savefig('standard_patterns.png', dpi='figure',
        orientation='landscape', papertype=None, format=None,
        transparent=True, bbox_inches='tight', pad_inches=0.0,
        frameon=False)
plt.show()

/home/user/python/venvs/py373/lib/python3.7/site-packages/matplotlib/cbook/__init__.py:424: MatplotlibDeprecationWarning:
Passing one of 'on', 'true', 'off', 'false' as a boolean is deprecated; use an actual boolean (True/False) instead.
  warn_deprecated("2.2", "Passing one of 'on', 'true', 'off', 'false' as a "

This page was generated from /home/docs/checkouts/readthedocs.org/user_builds/xdesign/checkouts/0.5.x/docs/source/demos/Parameterized.ipynb. Interactive online version:

Parameterized Phantom Generation

Demonstrates how a parameterized function can generate 4 different phantoms from the same parameterized class.

from xdesign import *
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.gridspec as gridspec

SIZE = 600

np.random.seed(0) # random seed for repeatability
p1 = Foam(size_range=[0.05, 0.01], gap=0, porosity=1)
d1 = discrete_phantom(p1, SIZE)
plt.imshow(d1, cmap='viridis')
plt.show()

/home/user/python/venvs/py373/lib/python3.7/site-packages/xdesign/phantom.py:368: RuntimeWarning: Reached termination criteria of 500 attempts before adding all of the circles.
  kTERM_CRIT), RuntimeWarning)

np.random.seed(0) # random seed for repeatability
p2 = Foam(size_range=[0.07, 0.01], gap=0, porosity=0.75)
d2 = discrete_phantom(p2, SIZE)
plt.imshow(d2, cmap='viridis')
plt.show()

np.random.seed(0) # random seed for repeatability
p3 = Foam(size_range=[0.1, 0.01], gap=0, porosity=0.5)
d3 = discrete_phantom(p3, SIZE)
plt.imshow(d3, cmap='viridis')
plt.show()

np.random.seed(0) # random seed for repeatability
p4 = Foam(size_range=[0.1, 0.01], gap=0.015, porosity=0.5)
d4 = discrete_phantom(p4, SIZE)
plt.imshow(d4, cmap='viridis')
plt.show()

Create a composite figure of all four discrete phantoms.

fig = plt.figure(dpi=100)
gs1 = gridspec.GridSpec(2, 2)
gs1.update(wspace=0.1, hspace=0.1) # set the spacing between axes.

plt.subplot(gs1[0])
plt.title('(a)')
plt.axis('off')
plt.imshow(d1, interpolation='none', cmap=plt.cm.gray)
plt.subplot(gs1[1])
plt.title('(b)')
plt.axis('off')
plt.imshow(d2, interpolation='none', cmap=plt.cm.gray)
plt.subplot(gs1[2])
plt.title('(c)')
plt.axis('off')
plt.imshow(d3, interpolation='none', cmap=plt.cm.gray)
plt.subplot(gs1[3])
plt.title('(d)')
plt.axis('off')
plt.imshow(d4, interpolation='none', cmap=plt.cm.gray)
fig.set_size_inches(6, 6)
plt.savefig('Foam_parameterized.png', dpi='figure',
        orientation='landscape', papertype=None, format=None,
        transparent=True, bbox_inches='tight', pad_inches=0.0,
        frameon=False)
plt.show()

This page was generated from /home/docs/checkouts/readthedocs.org/user_builds/xdesign/checkouts/0.5.x/docs/source/demos/NoReferenceMetrics.ipynb. Interactive online version:

No Reference Metrics

Demonstrate the use of the no-reference metrics: noise power spectrum (NPS), modulation transfer function (MTF), and noise equivalent quanta (NEQ).

import matplotlib.pylab as plt
import numpy as np
import tomopy

from xdesign import *

Simulate data aqcuisition

Generate a UnitCircle standards phantom. For the modulation transfer function (MTF), the radius must be less than 0.5, otherwise the circle touches the edges of the field of view.

NPIXEL = 100

p = UnitCircle(radius=0.35, material=SimpleMaterial(7.23))
sidebyside(p, NPIXEL)
plt.show()

Noise power spectrum (NPS) and Noise Equivalent Quanta (NEQ) are meaningless without noise, so add some poisson noise to the simulated data using np.random.poisson. Collecting two sinograms allows us to isolate the noise by subtracting out the circle.

# Define the scan positions using theta and horizontal coordinates
theta, h = np.meshgrid(np.linspace(0, np.pi, NPIXEL, endpoint=False),
                       np.linspace(0, 1, NPIXEL, endpoint=False) - 0.5 + 1/NPIXEL/2)
# Reshape the returned arrays into vectors
theta = theta.flatten()
h = h.flatten()

num_photons = 1e4
# Define a probe that sends 1e4 photons per exposure
probe = Probe(size=1/NPIXEL, intensity=num_photons)
# Measure the phantom
data = probe.measure(p, theta, h)

/home/beams0/DCHING/Documents/xdesign/src/xdesign/geometry/algorithms.py:54: RuntimeWarning: halfspacecirc was out of bounds, -1.1249850495609337e-09
  RuntimeWarning)
/home/beams0/DCHING/Documents/xdesign/src/xdesign/geometry/algorithms.py:54: RuntimeWarning: halfspacecirc was out of bounds, -1.597134646758036e-10
  RuntimeWarning)
/home/beams0/DCHING/Documents/xdesign/src/xdesign/geometry/algorithms.py:54: RuntimeWarning: halfspacecirc was out of bounds, -7.460398965264403e-11
  RuntimeWarning)
/home/beams0/DCHING/Documents/xdesign/src/xdesign/geometry/algorithms.py:54: RuntimeWarning: halfspacecirc was out of bounds, -8.34076696598629e-11
  RuntimeWarning)
/home/beams0/DCHING/Documents/xdesign/src/xdesign/geometry/algorithms.py:54: RuntimeWarning: halfspacecirc was out of bounds, -1.9854995425561128e-10
  RuntimeWarning)

# Add poisson noise to the result.
noisy_data_0 = np.random.poisson(data)
noisy_data_1 = np.random.poisson(data)
# Normalize the data by the incident intensity
normalized_data_0 = noisy_data_0 / num_photons
normalized_data_1 = noisy_data_1 / num_photons

# Linearize the data by taking the negative log
sinoA = -np.log(normalized_data_0).reshape(NPIXEL, NPIXEL).T
sinoB = -np.log(normalized_data_1).reshape(NPIXEL, NPIXEL).T

plt.imshow(sinoA, cmap='viridis', origin='lower')
plt.colorbar()
plt.show()

# Reconstruct the data using TomoPy
recA = tomopy.recon(np.expand_dims(sinoA, 1), theta,
                    algorithm='gridrec', center=(sinoA.shape[1]-1)/2.) * NPIXEL
recB = tomopy.recon(np.expand_dims(sinoB, 1), theta,
                    algorithm='gridrec', center=(sinoB.shape[1]-1)/2.) * NPIXEL

plt.imshow(recA[0], cmap='inferno', interpolation="none")
plt.colorbar()
plt.savefig('UnitCircle_noise0.png', dpi=600,
        orientation='landscape', papertype=None, format=None,
        transparent=True, bbox_inches='tight', pad_inches=0.0,
        frameon=False)
plt.show()

plt.imshow(recB[0], cmap='inferno', interpolation="none")
plt.colorbar()
plt.show()

Calculate MTF

Friedman’s method

Use Friedman et al’s method for computing the MTF. You can separate the MTF into multiple directions or average them all together using the Ntheta argument.

mtf_freq, mtf_value, labels = compute_mtf_ffst(p, recA[0], Ntheta=4)

The MTF is really a symmetric function around zero frequency, so usually people just show the positive portion. Sometimes, there is a peak at the higher spatial frequencies instead of the MTF approaching zero. This is probably because of aliasing noise content with frequencies higher than the Nyquist frequency.

plot_mtf(mtf_freq, mtf_value, labels)
plt.gca().set_xlim([0,50]) # hide negative portion of MTF
plt.show()

Siemens star method

You can also use a Siemens Star to calculate the MTF using a fitted sinusoidal method instead of the slanted edges that the above method uses.

s = SiemensStar(n_sectors=32, center=Point([0, 0]), radius=0.5)
d = sidebyside(s, 100)
plt.show()

Here we are using the discreet verison of the phantom (without noise), so we are only limited by the resolution of the image.

mtf_freq, mtf_value = compute_mtf_lwkj(s, d)

plot_mtf(mtf_freq, mtf_value, labels=None)
plt.gca().set_xlim([0,50]) # hide portion of MTF beyond Nyquist frequency
plt.show()

Calculate NPS

Calculate the radial or 2D frequency plot of the NPS.

X, Y, NPS = compute_nps_ffst(p, recA[0], plot_type='frequency',B=recB[0])

plot_nps(X, Y, NPS)
plt.show()

bins, counts = compute_nps_ffst(p, recA[0], plot_type='histogram',B=recB[0])

plt.figure()
plt.bar(bins, counts)
plt.xlabel('spatial frequency [cycles/length]')
plt.title('Noise Power Spectrum')
plt.show()

Calculate NEQ

freq, NEQ = compute_neq_d(p, recA[0], recB[0])

plt.figure()
plt.plot(freq.flatten(), NEQ.flatten())
plt.xlabel('spatial frequency [cycles/length]')
plt.title('Noise Equivalent Quanta')
plt.show()

This page was generated from /home/docs/checkouts/readthedocs.org/user_builds/xdesign/checkouts/0.5.x/docs/source/demos/FullReferenceMetrics.ipynb. Interactive online version:

Quality Metrics and Reconstruction Demo

Demonstrate the use of full reference metrics by comparing the reconstruction of a simulated phantom using SIRT, ART, and MLEM.

import numpy as np
import matplotlib.pyplot as plt
from xdesign import *

NPIXEL = 128

Generate a phantom

Use one of XDesign’s various pre-made and procedurally generated phantoms.

np.random.seed(0)
soil_like_phantom = Softwood()

Generate a figure showing the phantom and save the discretized ground truth map for later.

discrete = sidebyside(soil_like_phantom, NPIXEL)
if False:
    plt.savefig('Soil_sidebyside.png', dpi='figure',
        orientation='landscape', papertype=None, format=None,
        transparent=True, bbox_inches='tight', pad_inches=0.0,
        frameon=False)
plt.show()

Simulate data acquisition

Generate a list of probe coordinates to simulate data acquisition for parallel beam around 180 degrees.

angles = np.linspace(0, np.pi, NPIXEL, endpoint=False),
positions = np.linspace(0, 1, NPIXEL, endpoint=False) - 0.5 + 1/NPIXEL/2
theta, h = np.meshgrid(angles, positions)

Make a probe.

probe = Probe(size=1/NPIXEL)

Use the probe to measure the phantom.

sino = probe.measure(soil_like_phantom, theta, h)

# Transform data from attenuated intensity to attenuation coefficient
sino = -np.log(sino)

plt.imshow(sino.reshape(NPIXEL, NPIXEL), cmap='viridis', origin='lower')
plt.show()

Reconstruct

Reconstruct the phantom using 3 different techniques: ART, SIRT, and MLEM.

niter = 32  # number of iterations
gmin = [-0.5, -0.5]
gsize = [1, 1]
data = sino

init = np.full((NPIXEL, NPIXEL), 1e-12)
rec_art = art(gmin, gsize, data, theta, h, init, niter)


init = np.full((NPIXEL, NPIXEL), 1e-12)
rec_sirt = sirt(gmin, gsize, data, theta, h, init, niter)


init = np.full((NPIXEL, NPIXEL), 1e-12)
rec_mlem = mlem(gmin, gsize, data, theta, h, init, niter)

[##########] 100.00%
[##########] 100.00%
[##########] 100.00%

plt.figure(figsize=(12,4), dpi=100)
plt.subplot(141)
plt.imshow(rec_art, cmap='gray', origin='lower', vmin=0, vmax=1)
plt.title('ART')
plt.subplot(142)
plt.imshow(rec_sirt, cmap='gray', interpolation='none', origin='lower', vmin=0, vmax=1)
plt.title('SIRT')
plt.subplot(143)
plt.imshow(rec_mlem, cmap='gray', interpolation='none', origin='lower', vmin=0, vmax=1)
plt.title('MLEM')
plt.subplot(144)
plt.imshow(discrete, cmap='gray', interpolation='none', origin='lower', vmin=0, vmax=1)
plt.title('truth')
plt.show()

Quality Metrics

Compute local quality for each reconstruction using MS-SSIM, a convolution based quality metric.

quality = list()
for rec in [rec_art, rec_sirt, rec_mlem]:
    scales, mscore, qmap = msssim(discrete, rec)
    quality.append(mscore)

Plot the average quality at for each reconstruction. Then display the local quality map for each reconstruction to see why certain reconstructions are ranked higher than others.

plt.figure()
plt.bar(["ART", "SIRT", "MLEM"], quality)
plt.show()

References

1

James T. Dobbins. Effects of undersampling on the proper interpretation of modulation transfer function, noise power spectra, and noise equivalent quanta of digital imaging systems. Medical Physics, 22(2):171–181, 1995. URL: http://scitation.aip.org/content/aapm/journal/medphys/22/2/10.1118/1.597600, doi:http://dx.doi.org/10.1118/1.597600.

2

Saul N. Friedman, George S. K. Fung, Jeffrey H. Siewerdsen, and Benjamin M. W. Tsui. A simple approach to measure computed tomography (ct) modulation transfer function (mtf) and noise-power spectrum (nps) using the american college of radiology (acr) accreditation phantom. Medical Physics, 40(5):, 2013. URL: http://scitation.aip.org/content/aapm/journal/medphys/40/5/10.1118/1.4800795, doi:http://dx.doi.org/10.1118/1.4800795.

3

Peter Gilbert. Iterative methods for the three-dimensional reconstruction of an object from projections. Journal of Theoretical Biology, 1972. URL: https://doi.org/10.1016/0022-5193(72)90180-4, doi:10.1016/0022-5193(72)90180-4.

4

Richard Gordon, Robert Bender, and Gabor T. Herman. Algebraic Reconstruction Techniques (ART) for three-dimensional electron microscopy and X-ray photography. Journal of Theoretical Biology, 29(3):471–481, dec 1970. URL: https://doi.org/10.1016/0022-5193(70)90109-8, doi:10.1016/0022-5193(70)90109-8.

5

J. D. Hunter. Matplotlib: a 2d graphics environment. Computing In Science & Engineering, 9(3):90–95, 2007.

6

Christian Loebich, Dietmar Wueller, Bruno Klingen, and Anke Jaeger. Digital camera resolution measurement using sinusoidal siemens stars. In Electronic Imaging 2007, 65020N–65020N. International Society for Optics and Photonics, 2007.

7

H. R. Sheikh and A. C. Bovik. Image information and visual quality. IEEE Transactions on Image Processing, 15(2):430–444, Feb 2006. doi:10.1109/TIP.2005.859378.

8

Zhou Wang and Alan C Bovik. A universal image quality index. IEEE signal processing letters, 9(3):81–84, 2002.

9

Zhou Wang, Eero P Simoncelli, and Alan C Bovik. Multiscale structural similarity for image quality assessment. In Signals, Systems and Computers, 2004. Conference Record of the Thirty-Seventh Asilomar Conference on, volume 2, 1398–1402. Ieee, 2003.

Credits

We kindly request that you cite the following article(s) [A1] if you use XDesign.

A1

Daniel J. Ching and Doǧa Gürsoy. XDesign: an open-source software package for designing X-ray imaging phantoms and experiments. Journal of Synchrotron Radiation, 24(2):537–544, Mar 2017. URL: https://doi.org/10.1107/S1600577517001928, doi:10.1107/S1600577517001928.