
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.
Citation¶
Ching, Daniel J., and Doga Gürsoy. “XDesign: an open-source software package for designing x-ray imaging phantoms and experiments.” Journal of synchrotron radiation 24, no. 2 (2017): 537-544.
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.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:
|
A square cross-section x-ray beam for probing Phantoms. |
Functions:
|
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)[source]¶
Bases:
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
- .. todo::
Implement additional attributes for Probe such as wavelength, etc.
- property cross_section¶
Return the cross-sectional area of a square beam.
- xdesign.acquisition.raster_scan2D(sa, st, meta=False)[source]¶
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:
|
Return the longest MURA whose length is less than or equal to L. |
|
Return the largest 2D MURA whose lengths are less than M and N. |
|
Return the coded mask from Raskar et al. |
- xdesign.codes.mura_1d(L)[source]¶
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 byA[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)[source]¶
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.constants
¶
Constants in cgs units.
xdesign.formats
¶
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:
|
Define a point in ND cartesian space. |
- class xdesign.geometry.point.Point(x)[source]¶
Bases:
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.
- contains(other)[source]¶
Return wether the other is within the bounds of the Point.
Points can only contain other Points.
- property norm¶
Euclidian (L2) norm of the vector to the point.
- scale(vector)[source]¶
Scale the ambient space in each dimension according to vector.
Scaling is centered on the origin.
- 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)[source]¶
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-planec1 (
np.array
(…, 1)) – The last coefficient for the hyper-plane
One dimensional entities¶
Define one dimensional geometric entities.
Classes:
|
Line in 2D cartesian space. |
|
Defines a finite line segment from two unique points. |
- class xdesign.geometry.line.Line(p1, p2)[source]¶
Bases:
LinearEntity
Line in 2D cartesian space.
The constructor takes two unique
Point
.- 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.
Two dimensional entities¶
Define two dimensional geometric entities.
Classes:
|
The base class for closed manifolds defined by a single equation. |
|
Circle in 2D cartesian space. |
|
A convex polygon in 2D cartesian space. |
|
A regular polygon in 2D cartesian space. |
|
Triangle in 2D cartesian space. |
|
Rectangle in 2D cartesian space. |
|
Square in 2D cartesian space. |
|
A collection of Entities |
- class xdesign.geometry.area.Circle(center, radius, sign=1)[source]¶
Bases:
Curve
Circle in 2D cartesian space.
- radius¶
The radius of the circle.
- Type
scalar
- 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)[source]¶
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)[source]¶
Bases:
Entity
The base class for closed manifolds defined by a single equation. e.g.
Circle
,Sphere
, orTorus
.
- class xdesign.geometry.area.Mesh(obj=None, faces=[])[source]¶
Bases:
Entity
A collection of Entities
- property bounding_box¶
Return the axis-aligned bounding box as two numpy vectors.
- property center¶
- contains(other)[source]¶
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).
- property patch¶
- class xdesign.geometry.area.Polygon(vertices, sign=1)[source]¶
Bases:
Entity
A convex polygon in 2D cartesian space.
It is defined by a number of distinct vertices of class
Point
. Superclasses includeSquare
,Triangle
, etc.- vertices¶
- Type
List of Points
:raises ValueError : If the number of vertices is less than three.:
- area = <MagicMock name='mock()' id='139671082956112'>¶
- 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='139671082956112'>¶
- property edges¶
Return a list of lines connecting the points of the Polygon.
- half_space = <MagicMock name='mock()' id='139671082956112'>¶
- property list¶
Return list representation.
- property numpy¶
Return Numpy representation.
- property numverts¶
- property patch¶
Returns a matplotlib patch.
- perimeter = <MagicMock name='mock()' id='139671082956112'>¶
- radius = <MagicMock name='mock()' id='139671082956112'>¶
- class xdesign.geometry.area.Rectangle(center, side_lengths)[source]¶
Bases:
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='139671082956112'>¶
- class xdesign.geometry.area.RegularPolygon(center, radius, order, angle=0, sign=1)[source]¶
Bases:
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 ofcenter
and \(r\) =radius
. Theangle
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 polygonradius (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
)
Intersect¶
Define algorithms to support intersection calculation.
Functions:
|
Clip a polygon using the Sutherland-Hodgeman algorithm. |
|
Return the area of intersection between a circle and half-plane. |
- xdesign.geometry.intersect.clip_SH(clipEdges, polygon)[source]¶
Clip a polygon using the Sutherland-Hodgeman algorithm.
- Parameters
clipEdges [[A, b], …] – half-spaces defined by coefficients
polygon
- xdesign.geometry.intersect.halfspacecirc(d, r)[source]¶
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:
|
Simple material with constant mass_attenuation parameter only. |
|
Materials which use xraylib data to automatically calculate material properties based on beam energy in KeV. |
- class xdesign.material.SimpleMaterial(mass_attenuation=1.0)[source]¶
Bases:
Material
Simple material with constant mass_attenuation parameter only.
- class xdesign.material.XraylibMaterial(compound, density)[source]¶
Bases:
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
- beta()¶
x.__getitem__(y) <==> x[y]
- delta()¶
x.__getitem__(y) <==> x[y]
xdesign.metrics
¶
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:
|
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)[source]¶
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:
|
Return the Pearson product-moment correlation coefficients (PCC). |
|
Return the Structural SIMilarity index (SSIM) of two images. |
|
Multi-Scale Structural SIMilarity index (MS-SSIM). |
- class xdesign.metrics.fullref.ImageQuality(original, reconstruction)[source]¶
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
- method¶
The metric used to calculate the quality
- Type
string
- quality(method='MSSSIM', L=1.0, **kwargs)[source]¶
Compute the full-reference image quality of each image pair.
Available methods include SSIM [4], MSSSIM [5], VIFp [3]
- 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)[source]¶
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)[source]¶
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)[source]¶
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
- 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
- xdesign.metrics.fullref.beta_gamma¶
The exponent which weights the contribution of the contrast and structure terms at each level.
- Type
- 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¶
xdesign.phantom
¶
Phantoms¶
Standard phantoms¶
Custom phantoms¶
xdesign.plot
¶
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:
|
Reconstruct data using ART algorithm. |
|
Reconstruct data using SIRT algorithm. |
|
Reconstruct data using MLEM algorithm. |
|
Draw a process bar in the terminal. |
- xdesign.recon.art(gmin, gsize, data, theta, h, init, niter=10, weights=None, save_interval=None)[source]¶
Reconstruct data using ART algorithm. [2].
- xdesign.recon.mlem(gmin, gsize, data, theta, h, init, niter=10)[source]¶
Reconstruct data using MLEM algorithm.
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/latest/docs/source/demos/Shepp.ipynb.
Interactive online version:
Simple How-to Explaining Phantoms¶
Demonstrate simple basic custom phantom and sinogram generation with XDesign.
[1]:
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.
[2]:
# 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.
[3]:
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
.
[4]:
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.
[11]:
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.
[6]:
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
[7]:
probe = Probe(size=1/NPIXEL)
[8]:
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.
[9]:
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/latest/docs/source/demos/StandardPatterns.ipynb.
Interactive online version:
Standard Test Patterns¶
Generates sidebyside plots of all the standard test patterns in xdesign.
[1]:
from xdesign import *
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
[2]:
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()

[3]:
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()

[4]:
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()

[5]:
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()

[6]:
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()

[7]:
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()

[8]:
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/latest/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.
[1]:
from xdesign import *
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.gridspec as gridspec
SIZE = 600
[2]:
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)

[3]:
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()

[4]:
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()

[5]:
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.
[6]:
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/latest/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).
[1]:
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.
[2]:
NPIXEL = 100
[3]:
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.
[4]:
# 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()
[5]:
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/intersect.py:57: RuntimeWarning: halfspacecirc was out of bounds, -1.1249850495609337e-09
"halfspacecirc was out of bounds, {}".format(f), RuntimeWarning
/home/beams0/DCHING/Documents/xdesign/src/xdesign/geometry/intersect.py:57: RuntimeWarning: halfspacecirc was out of bounds, -1.597134646758036e-10
"halfspacecirc was out of bounds, {}".format(f), RuntimeWarning
/home/beams0/DCHING/Documents/xdesign/src/xdesign/geometry/intersect.py:57: RuntimeWarning: halfspacecirc was out of bounds, -7.460398965264403e-11
"halfspacecirc was out of bounds, {}".format(f), RuntimeWarning
/home/beams0/DCHING/Documents/xdesign/src/xdesign/geometry/intersect.py:57: RuntimeWarning: halfspacecirc was out of bounds, -8.34076696598629e-11
"halfspacecirc was out of bounds, {}".format(f), RuntimeWarning
/home/beams0/DCHING/Documents/xdesign/src/xdesign/geometry/intersect.py:57: RuntimeWarning: halfspacecirc was out of bounds, -1.9854995425561128e-10
"halfspacecirc was out of bounds, {}".format(f), RuntimeWarning
[6]:
# 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
[7]:
# 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
[8]:
plt.imshow(sinoA, cmap='viridis', origin='lower')
plt.colorbar()
plt.show()

[9]:
# 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
[10]:
plt.imshow(recA[0], cmap='inferno', interpolation="none")
plt.colorbar()
plt.show()

[11]:
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.
[12]:
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.
[13]:
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.
[14]:
s = SiemensStar(n_sectors=16, center=Point([0, 0]), radius=0.5)
d = sidebyside(s, 256)
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.
[15]:
mtf_freq, mtf_value = compute_mtf_lwkj(d, s.n_sectors)
[16]:
plot_mtf(mtf_freq, mtf_value, labels=None)
plt.show()

Calculate NPS¶
Calculate the radial or 2D frequency plot of the NPS.
[17]:
X, Y, NPS = compute_nps_ffst(p, recA[0], plot_type='frequency',B=recB[0])
[18]:
plot_nps(X, Y, NPS)
plt.show()

[19]:
bins, counts = compute_nps_ffst(p, recA[0], plot_type='histogram',B=recB[0])
[20]:
plt.figure()
plt.bar(bins, counts)
plt.xlabel('spatial frequency [cycles/length]')
plt.title('Noise Power Spectrum')
plt.show()

Calculate NEQ¶
[21]:
freq, NEQ = compute_neq_d(p, recA[0], recB[0])
[22]:
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/latest/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.
[1]:
import numpy as np
import matplotlib.pyplot as plt
from xdesign import *
[2]:
NPIXEL = 128
Generate a phantom¶
Use one of XDesign’s various pre-made and procedurally generated phantoms.
[3]:
np.random.seed(0)
soil_like_phantom = Softwood()
Generate a figure showing the phantom and save the discretized ground truth map for later.
[4]:
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.
[5]:
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.
[6]:
probe = Probe(size=1/NPIXEL)
Use the probe to measure the phantom.
[7]:
sino = probe.measure(soil_like_phantom, theta, h)
[8]:
# Transform data from attenuated intensity to attenuation coefficient
sino = -np.log(sino)
[9]:
plt.imshow(sino.reshape(NPIXEL, NPIXEL), cmap='viridis', origin='lower')
plt.show()

Reconstruct¶
Reconstruct the phantom using 3 different techniques: ART, SIRT, and MLEM.
[10]:
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%
[11]:
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.
[12]:
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.
[13]:
plt.figure()
plt.bar(["ART", "SIRT", "MLEM"], quality)
plt.show()

[ ]:
[ ]:
References¶
- 1
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.
- 2
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.
- 3
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.
- 4
Zhou Wang and Alan C Bovik. A universal image quality index. IEEE signal processing letters, 9(3):81–84, 2002.
- 5
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.