XDesign¶

xdesign.acquisition¶

Functions:

class xdesign.acquisition.Beam(p1, p2, size=0)[source]¶

Bases: xdesign.geometry.Line

Beam (thick line) in 2-D cartesian space.

It is defined by two distinct points.

p1¶

Point

p2¶

Point

size¶

scalar, optional

Size of the beam.

half_space¶

Returns the half space polytope respresentation of the infinite beam.

class xdesign.acquisition.Probe(p1, p2, size=0)[source]¶

Bases: xdesign.acquisition.Beam

measure(phantom, noise=False)[source]¶

Return the probe measurement given phantom. When noise is > 0, poisson noise is added to the returned measurement.

record()[source]¶

translate(dx)[source]¶

Translates beam along its normal direction.

xdesign.acquisition.sinogram(sx, sy, phantom, noise=False)[source]¶

Generates sinogram given a phantom.

Parameters:	sx (int) – Number of rotation angles. sy (int) – Number of detection pixels (or sample translations). phantom (Phantom)
Returns:	ndarray – Sinogram.

xdesign.acquisition.angleogram(sx, sy, phantom, noise=False)[source]¶

Generates angleogram given a phantom.

Parameters:	sx (int) – Number of rotation angles. sy (int) – Number of detection pixels (or sample translations). phantom (Phantom)
Returns:	ndarray – Angleogram.

xdesign.algorithms¶

Functions:

xdesign.algorithms.reconstruct(probe, shape, data, rec)[source]¶

Reconstruct single x-ray beam data.

Parameters:	probe (Probe) shape (list) – shape of the reconstruction grid. data (scalar) – Probe measurement dtaa. rec (ndarray) – Initialization matrix.

xdesign.algorithms.art(probe, data, init, niter=10)[source]¶

Reconstruct data.

xdesign.algorithms.sirt(probe, data, init, niter=10)[source]¶

Reconstruct data.

xdesign.algorithms.mlem(probe, data, init, niter=10)[source]¶

Reconstruct data.

xdesign.algorithms.stream(probe, data, init)[source]¶

Reconstruct data.

xdesign.algorithms.update_progress(progress)[source]¶

Draws a process bar in the terminal.

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

xdesign.feature¶

Functions:

class xdesign.feature.Feature(geometry, mass_atten=1)[source]¶

Bases: object

A geometric region(s) and associated materials properti(es). Properties and geometry can be manipulated from this level, but rendering must access the geometry directly.

geometry¶

Entity

Defines a region where the properties are valid.

mass_atten¶

scalar

The mass attenuation coefficient of the Feature.

add_property(name, function)[source]¶

Adds a property by name to the Feature. NOTE: Properties added here are not cached because they are probably never called with the same parameters ever or the property is static so it doesnt need caching.

area¶

Returns the total surface area of the feature

center¶

Returns the centroid of the feature.

radius¶

Returns the radius of the smallest boundary circle

rotate(theta, p, axis=None)[source]¶

Rotate feature geometry around a line. Rotating property functions is not supported.

translate(x, y)[source]¶

Translate feature geometry. Translating property functions is not supported.

volume¶

Returns the volume of the feature

xdesign.geometry¶

Functions:

class xdesign.geometry.Entity[source]¶

Bases: object

Base class for all geometric entities. All geometric entities should have these methods.

collision(other)[source]¶

Returns True if this entity collides with another entity.

contains(point)[source]¶

Returns True if the point(s) is within the bounds of the entity. point is either a Point or an N points listed as an NxD array.

dim¶

The dimensionality of the entity

distance(other)[source]¶

Return the closest distance between entities.

midpoint(other)[source]¶

Return the midpoint between entities.

rotate(theta, point=None, axis=None)[source]¶

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

scale(vector)[source]¶

Scale entity in each dimension according to vector. Scaling is centered on the origin.

translate(vector)[source]¶

Translate entity along vector.

class xdesign.geometry.Point(x)[source]¶

Bases: xdesign.geometry.Entity

A point in ND cartesian space.

_x¶

1darray

ND coordinates of the point.

x¶

scalar

dimension 0 of the point.

y¶

scalar

dimension 1 of the point.

z¶

scalar

dimension 2 of the point.

norm¶

scalar

The L2/euclidian norm of the point.

collision(other)[source]¶

Returns True if this entity collides with another entity.

contains(other)[source]¶

Returns True if the other is within the bounds of the entity.

dim¶

distance(other)[source]¶

Return the closest distance between entities.

norm

Reference: http://stackoverflow.com/a/23576322

rotate(theta, point=None, axis=None)[source]¶

Rotate entity around an axis by theta radians.

scale(vector)[source]¶

Scale entity in each dimension according to vector. Scaling is centered on the origin.

translate(vector)[source]¶

Translate point along vector.

x

y

z

class xdesign.geometry.Circle(center, radius)[source]¶

Bases: xdesign.geometry.Curve

Circle in 2-D cartesian space.

center¶

Point

Defines the center point of the circle.

radius¶

scalar

Radius of the circle.

area¶

Return area.

circumference¶

Return circumference.

contains(points)[source]¶

diameter¶

Return diameter.

list¶

Return list representation for saving to files.

patch¶

scale(val)[source]¶

Scale.

class xdesign.geometry.Line(p1, p2)[source]¶

Bases: xdesign.geometry.LinearEntity

Line in 2-D cartesian space.

It is defined by two distinct points.

p1¶

Point

p2¶

Point

distance(other)[source]¶

Return the closest distance between entities. REF: http://geomalgorithms.com/a02-_lines.html

intercept(n)[source]¶

Calculates the intercept for the nth dimension.

standard¶

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

xintercept¶

Return the x-intercept.

yintercept¶

Return the y-intercept.

class xdesign.geometry.Polygon(vertices)[source]¶

Bases: xdesign.geometry.Entity

A convex polygon in 2-D cartesian space.

It is defined by a number of distinct points.

vertices¶

sequence of Points

area¶

Return the area of the entity.

bounds¶

Return a tuple (xmin, ymin, xmax, ymax) representing the bounding rectangle for the geometric figure.

contains(points)[source]¶

half_space¶

list¶

Return list representation.

numpy¶

Return Numpy representation.

numverts¶

patch¶

perimeter¶

Return the perimeter of the entity.

rotate(theta, point=None, axis=None)[source]¶

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

translate(vector)[source]¶

Translate polygon.

class xdesign.geometry.Triangle(p1, p2, p3)[source]¶

Bases: xdesign.geometry.Polygon

Triangle in 2-D cartesian space.

It is defined by three distinct points.

area¶

center¶

class xdesign.geometry.Rectangle(p1, p2, p3, p4)[source]¶

Bases: xdesign.geometry.Polygon

Rectangle in 2-D cartesian space.

It is defined by four distinct points.

area¶

center¶

class xdesign.geometry.Square(center, side_length)[source]¶

Bases: xdesign.geometry.Rectangle

Square in 2-D cartesian space.

It is defined by a center and a side length.

class xdesign.geometry.Mesh(obj=None)[source]¶

Bases: xdesign.geometry.Entity

A mesh object. It is a collection of polygons

append(t)[source]¶

Add a triangle to the mesh.

center¶

contains(points)[source]¶

half_space¶

import_triangle(obj)[source]¶

Loads mesh data from a Python Triangle dict.

patch¶

pop(i=-1)[source]¶

Pop i-th triangle from the mesh.

rotate(theta, point=None, axis=None)[source]¶

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

scale(vector)[source]¶

Scale entity.

translate(vector)[source]¶

Translate entity.

xdesign.material¶

Functions:

class xdesign.material.HyperbolicConcentric(min_width=0.1, exponent=0.5)[source]¶

Bases: xdesign.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¶

list

The list of radii of the circles

widths¶

list

The list of the widths of the bands

class xdesign.material.DynamicRange(steps=10, jitter=True, shape=u'square')[source]¶

Bases: xdesign.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.material.DogaCircles(n_sizes=5, size_ratio=0.5, n_shuffles=5)[source]¶

Bases: xdesign.phantom.Phantom

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

radii¶

ndarray

radii of circles

x¶

ndarray

x position of circles

y¶

ndarray

y position of circles

class xdesign.material.SlantedSquares(count=10, angle=0.08726646259972222, gap=0)[source]¶

Bases: xdesign.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¶

scalar

the angle of slant in radians

count¶

scalar

the total number of squares

gap¶

scalar

the minimum space between squares

side_length¶

scalar

the size of the squares

squares_per_level¶

list

the number of squares at each level

radius_per_level¶

list

the radius at each level

n_levels¶

scalar

the number of levels

class xdesign.material.UnitCircle(radius=0.5, mass_atten=1)[source]¶

Bases: xdesign.phantom.Phantom

Generates a phantom with a single circle in its center.

class xdesign.material.WetCircles[source]¶

Bases: xdesign.phantom.Phantom

class xdesign.material.SiemensStar(n_sectors=4, center=<xdesign.geometry.Point object>, radius=0.5)[source]¶

Bases: xdesign.phantom.Phantom

Generates a Siemens star.

ratio¶

scalar

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

class xdesign.material.Soil(porosity=0.412)[source]¶

Bases: xdesign.phantom.Phantom

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.material.Foam(size_range=[0.05, 0.01], gap=0, porosity=1)[source]¶

Bases: xdesign.phantom.Phantom

Generates a phantom with structure similar to foam.

xdesign.metrics¶

Functions:

class xdesign.metrics.ImageQuality(original, reconstruction, method=u'')[source]¶

Bases: object

Stores information about image quality.

orig¶

numpy.ndarray

recon¶

numpy.ndarray

qualities¶

list of scalars

maps¶

list of numpy.ndarray

scales¶

list of scalars

add_quality(quality, scale, maps=None)[source]¶

Parameters:	quality (scalar, list) – The average quality for the image map (array, list of arrays, optional) – the local quality rating across the image scale (scalar, list) – the size scale at which the quality was calculated

sort()[source]¶

Sorts the qualities by scale

xdesign.metrics.probability_mask(phantom, size, ratio=8, uniform=True)[source]¶

Returns the probability mask for each phase in the phantom.

Parameters:

size (scalar) – The side length in pixels of the resulting square image. ratio (scalar, optional) – The discretization works by drawing the shapes in a larger space then averaging and downsampling. This parameter controls how many pixels in the larger representation are averaged for the final representation. e.g. if ratio = 8, then the final pixel values are the average of 64 pixels. uniform (boolean, optional) – When set to False, changes the way pixels are averaged from a uniform weights to gaussian weights.

Returns:

image (list of numpy.ndarray) – A list of float masks for each phase in the phantom.

xdesign.metrics.compute_quality(reference, reconstructions, method=u'MSSSIM', L=1)[source]¶

Computes image quality metrics for each of the reconstructions.

Parameters:

reference (array) – the discrete reference image. In a future release, we will determine the best way to compare a continuous domain to a discrete reconstruction. reconstructions (list of arrays) – A list of discrete reconstructions method (string, optional) – The quality metric desired for this comparison. Options include: SSIM, MSSSIM L (scalar) – The dynamic range of the data. This value is 1 for float representations and 2^bitdepth for integer representations.

Returns:

metrics (list of ImageQuality)

xdesign.metrics.compute_PCC(A, B, masks=None)[source]¶

Computes the Pearson product-moment correlation coefficients (PCC) for the two images.

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.compute_likeness(A, B, masks)[source]¶

Predicts the likelihood that each pixel in B belongs to a phase based on the histogram of A.

Parameters:	A (ndarray) B (ndarray) masks (list of ndarrays)
Returns:	likelihoods (list of ndarrays)

xdesign.metrics.compute_background_ttest(image, masks)[source]¶

Determines whether the background has significantly different luminance than the other phases.

Parameters:	image (ndarray) masks (list of ndarrays) – Masks for the background and any other phases. Does not autogenerate the non-background mask because maybe you want to compare only two phases.
Returns:	tstat (scalar) pvalue (scalar)

xdesign.metrics.compute_mtf(phantom, image, Ntheta=4)[source]¶

Uses method described in Friedman et al to calculate the MTF.

Parameters:	phantom (UnitCircle) – Predefined phantom with single circle whose radius is less than 0.5. image (ndarray) – The reconstruction of the above 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

References

S. N. Friedman, G. S. K. Fung, J. H. Siewerdsen, and B. 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,” Med. Phys. 40, 051907-1 - 051907-9 (2013). http://dx.doi.org/10.1118/1.4800795

xdesign.metrics.compute_nps(phantom, A, B=None, plot_type=u'frequency')[source]¶

Calculates the noise power spectrum from a unit circle image. The peak at low spatial frequency is probably due to aliasing. Invesigation into supressing this peak is necessary.

References

S. N. Friedman, G. S. K. Fung, J. H. Siewerdsen, and B. 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,” Med. Phys. 40, 051907-1 - 051907-9 (2013). http://dx.doi.org/10.1118/1.4800795

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.metrics.compute_neq(phantom, A, B)[source]¶

Calculates the NEQ according to recommendations by JT Dobbins.

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.compute_mtf_siemens(phantom, image)[source]¶

Calculates the MTF using the modulated Siemens Star method in Loebich et al. (2007).

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.phantom¶

Functions:

class xdesign.phantom.Phantom(shape=u'circle')[source]¶

Bases: object

Phantom generation class.

shape¶

string

Shape of the phantom. Available options: circle, square.

population¶

scalar

Number of generated features in the phantom.

area¶

scalar

Area of the features in the phantom.

feature¶

list

List of features.

append(feature)[source]¶

Add a feature to the top of the phantom.

density¶

Returns the area density of the phantom. (Does not acount for mass_atten.)

insert(i, feature)[source]¶

Insert a feature at a given depth.

list¶

load(filename)[source]¶

Load phantom from file.

numpy()[source]¶

Returns the Numpy representation.

pop(i=-1)[source]¶

Pop i-th feature from the phantom.

reverse()[source]¶

Reverse the features of the phantom

rotate(theta, origin=<xdesign.geometry.Point object>, axis=None)[source]¶

Rotate phantom around a point.

save(filename)[source]¶

Save phantom to file.

sort(param=u'mass_atten', reverse=False)[source]¶

Sorts the features by mass_atten or size

sprinkle(counts, radius, gap=0, region=None, mass_atten=1, max_density=1)[source]¶

Sprinkle a number of circles. 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. mass_atten (scalar, optional) – A mass attenuation parameter passed to the circles.

Returns:

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

translate(dx, dy)[source]¶

Translate phantom.

xdesign.plot¶

Functions:

xdesign.plot.plot_phantom(phantom, axis=None, labels=None, c_props=[], c_map=None)[source]¶

Plots a phantom to the given axis.

Parameters:	labels (bool, optional) – Each feature is numbered according to its index in the phantom. c_props (list of str, optional) – List of feature properties to use for colormapping the geometries. c_map (function, optional) – A function which takes the a list of prop(s) for a Feature as input and returns a matplolib color specifier.

xdesign.plot.plot_feature(feature, axis=None, alpha=None, c=None)[source]¶

Plots a feature on the given axis.

Parameters:	alpha (float) – The plot opaqueness. 0 is transparent. 1 is opaque. c (matplotlib color specifier) – See http://matplotlib.org/api/colors_api.html

xdesign.plot.plot_mesh(mesh, axis=None, alpha=None, c=None)[source]¶

Plots a mesh to the given axis.

Parameters:	alpha (float) – The plot opaqueness. 0 is transparent. 1 is opaque. c (matplotlib color specifier) – See http://matplotlib.org/api/colors_api.html

xdesign.plot.plot_polygon(polygon, axis=None, alpha=None, c=None)[source]¶

Plots a polygon to the given axis.

Parameters:	alpha (float) – The plot opaqueness. 0 is transparent. 1 is opaque. c (matplotlib color specifier) – See http://matplotlib.org/api/colors_api.html

xdesign.plot.plot_curve(curve, axis=None, alpha=None, c=None)[source]¶

Plots a curve to the given axis.

Parameters:	alpha (float) – The plot opaqueness. 0 is transparent. 1 is opaque. c (matplotlib color specifier) – See http://matplotlib.org/api/colors_api.html

xdesign.plot.discrete_phantom(phantom, size, ratio=8, uniform=True, prop=u'mass_atten')[source]¶

Returns discrete representation of the property function, prop, in the phantom. The values of overlapping features are additive.

Parameters:

size (scalar) – The side length in pixels of the resulting square image. ratio (scalar, optional) – 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 = 8, then the final pixel values are the average of 64 pixels. uniform (boolean, optional) – When set to False, changes the way pixels are averaged from a uniform weights to gaussian weigths. prop (str, optional) – The name of the property function to discretize

Returns:

image (numpy.ndarray) – The discrete representation of the phantom that is size x size.

xdesign.plot.sidebyside(p, size=100, labels=None, prop=u'mass_atten')[source]¶

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

xdesign.plot.plot_metrics(imqual)[source]¶

Plots full reference metrics of ImageQuality data.

Parameters:	imqual (ImageQuality) – The data to plot.

References

Colors taken from this gist <https://gist.github.com/thriveth/8560036>

xdesign.plot.plot_mtf(faxis, MTF, labels=None)[source]¶

Plots the MTF. Returns the figure reference.

xdesign.plot.plot_nps(X, Y, NPS)[source]¶

Plots the 2D frequency plot for the NPS. Returns the figure reference.

xdesign.plot.plot_neq(freq, NEQ)[source]¶

Plots the NEQ. Returns the figure reference.

Quality Metrics and Reconstruction Demo¶

Demonstrates the use of full reference metrics.

import numpy as np
import matplotlib.pyplot as plt
from skimage.exposure import adjust_gamma, rescale_intensity
from xdesign import *


def rescale(reconstruction, hi):
    I = rescale_intensity(reconstruction, out_range=(0., 1.))
    return adjust_gamma(I, 1, hi)

Generate a predetermined phantom that resembles soil.

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

Generate a figure showing the phantom and save the discrete conversion for later.

discrete = sidebyside(soil_like_phantom, 100)
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 for parallel beam around 180 degrees.

sx, sy = 100, 100
step = 1. / sy
prb = Probe(Point([step / 2., -10]), Point([step / 2., 10]), step)
theta = np.pi / sx
sino = np.zeros(sx * sy)
a = 0
for m in range(sx):
    for n in range(sy):
        update_progress((m*sy + n)/(sx*sy))
        sino[a] = prb.measure(soil_like_phantom)
        a += 1
        prb.translate(step)
    prb.translate(-1)
    prb.rotate(theta, Point([0.5, 0.5]))
update_progress(1)

plt.imshow(np.reshape(sino, (sx, sy)), cmap='gray', interpolation='nearest')
# plt.hist(sino)
plt.show(block=True)

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

hi = 1  # highest expected value in reconstruction (for rescaling)
niter = 10  # number of iterations

init = 1e-12 * np.ones((sx, sy))
rec_art = art(prb, sino, init, niter)
rec_art = rescale(np.rot90(rec_art)[::-1], hi)

init = 1e-12 * np.ones((sx, sy))
rec_sirt = sirt(prb, sino, init, niter)
rec_sirt = rescale(np.rot90(rec_sirt)[::-1], hi)

init = 1e-12 * np.ones((sx, sy))
rec_mlem = mlem(prb, sino, init, niter)
rec_mlem = rescale(np.rot90(rec_mlem)[::-1], hi)

plt.figure(figsize=(12,4))
plt.subplot(131)
plt.imshow(rec_art, cmap='gray', interpolation='none')
plt.title('ART')
plt.subplot(132)
plt.imshow(rec_sirt, cmap='gray', interpolation='none')
plt.title('SIRT')
plt.subplot(133)
plt.imshow(rec_mlem, cmap='gray', interpolation='none')
plt.title('MLEM')
plt.show()

Compute quality metrics using the MSSSIM method.

metrics = compute_quality(discrete, [rec_art, rec_sirt, rec_mlem], method="MSSSIM")

Plotting the results shows the average quality at each level for each reconstruction in a line plot. Then it produces the local quality maps for each reconstruction so we might see why certain reconstructions are ranked higher than others.

In this case, it’s clear that ART is ranking higer than MLEM and SIRT at smaller scales because the small dark particles are visible; whereas for SIRT and MLEM they are unresolved. We can also see that the large yellow circles have a more accurately rendered luminance for SIRT and MLEM which is what causes these methods to be ranked higher at larger scales.

plot_metrics(metrics)
plt.show()

No Reference Metrics¶

Demonstrates the use of the no-reference metrics: NPS, MTF, and NEQ.

from xdesign import *
import tomopy
import numpy as np
import matplotlib.pylab as plt

Generate a UnitCircle test phantom. For the MTF, the radius must be less than 0.5, otherwise the circle touches the edges of the field of view.

p = UnitCircle(mass_atten=4, radius=0.35)
sidebyside(p, 100)
plt.show()

Generate two sinograms and reconstruct. Noise power spectrum and Noise Equivalent Quanta are meaningless withouth noise so add some poisson noise to the reconstruction process with the noise argument. Collecting two sinograms allows us to isolate the noise by subtracting out the circle.

np.random.seed(0)
sinoA = sinogram(100, 100, p, noise=0.1)
sinoB = sinogram(100, 100, p, noise=0.1)
theta = np.arange(0, np.pi, np.pi / 100.)

recA = tomopy.recon(np.expand_dims(sinoA, 1), theta, algorithm='gridrec', center=(sinoA.shape[1]-1)/2.)
recB = tomopy.recon(np.expand_dims(sinoB, 1), theta, algorithm='gridrec', center=(sinoB.shape[1]-1)/2.)

Take a look at the two noisy reconstructions.

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¶

This metric is meaningful without noise. You can separate the MTF into multiple directions or average them all together using the Ntheta argument.

mtf_freq, mtf_value, labels = compute_mtf(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()

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.5, 0.5]), 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_siemens(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¶

You can also calculate the radial or 2D frequency plot of the NPS.

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

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

bins, counts = compute_nps(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(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()

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()

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=600)
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()

Simple Phantom Construction Demo¶

Demonstrates simple basic custom phantom and sinogram generation.

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

Create various Features with geometries and assign attenuation values to each of the Features.

head = Feature(Circle(Point([0.5, 0.5]), radius=0.5))
head.mass_atten = 1
eyeL = Feature(Circle(Point([0.3, 0.5]), radius=0.1))
eyeL.mass_atten = 1
eyeR = Feature(Circle(Point([0.7, 0.5]), radius=0.1))
eyeR.mass_atten = 1
mouth = Feature(Triangle(Point([0.2, 0.7]), Point([0.5, 0.8]), Point([0.8, 0.7])))
mouth.mass_atten = -1

Add ‘Features’ to ‘Phantom’.

Shepp = Phantom()
Shepp.append(head)
Shepp.append(eyeL)
Shepp.append(eyeR)
Shepp.append(mouth)

Plot the Phantom geometry and properties with a colorbar.

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

# plot geometry
axis = fig.add_subplot(121, aspect='equal')
plt.grid('on')
plt.gca().invert_yaxis()
plot_phantom(Shepp, axis=axis, labels=False)

# plot property
plt.subplot(1, 2, 2)
im = plt.imshow(discrete_phantom(Shepp, 100, prop='mass_atten'), interpolation='none', cmap=plt.cm.inferno)

# 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 for parallel beam around 180 degrees.

sx, sy = 100, 100
step = 1. / sy
prb = Probe(Point([step / 2., -10]), Point([step / 2., 10]), step)
theta = np.pi / sx
sino = np.zeros(sx * sy)

a = 0
for m in range(sx):
    for n in range(sy):
        update_progress((m*sy + n)/(sx*sy))
        sino[a] = prb.measure(Shepp)
        a += 1
        prb.translate(step)
    prb.translate(-1)
    prb.rotate(theta, Point([0.5, 0.5]))
update_progress(1)

Plot the sinogram.

plt.figure(figsize=(8, 8))
plt.imshow(np.reshape(sino, (sx, sy)), 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()

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, mass_atten=1)
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, shape='square')
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.5, 0.5]))
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=600)
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)
    plt.grid('on')
    plt.gca().invert_yaxis()
    plot_phantom(phantoms[i], axis=axis)
    plt.title('('+ letters[i] +')')
    plt.subplot(gs1[2*i+1], aspect=1)
    plt.imshow(discrete_phantom(phantoms[i], 200), cmap='inferno')

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()

Wet Circles¶

import numpy as np
from scipy.spatial import Delaunay
import matplotlib.pyplot as plt
from xdesign import *
from skimage.exposure import adjust_gamma, rescale_intensity


def rescale(reconstruction, hi):
    I = rescale_intensity(reconstruction, out_range=(0., 1.))
    return adjust_gamma(I, 1, hi)

wet = WetCircles()
sidebyside(wet, size=200)
plt.savefig('Wet_sidebyside.png', dpi='figure',
        orientation='landscape', papertype=None, format=None,
        transparent=True, bbox_inches='tight', pad_inches=0.0,
        frameon=False)
plt.show(block=True)

sx, sy = 100, 100
step = 1. / sy
prb = Probe(Point([step / 2., -10]), Point([step / 2., 10]), step)
theta = np.pi / sx
sino = np.zeros(sx * sy)
a = 0
for m in range(sx):
    for n in range(sy):
        update_progress((m*sy + n)/(sx*sy))
        sino[a] = prb.measure(wet)
        a += 1
        prb.translate(step)
    prb.translate(-1)
    prb.rotate(theta, Point([0.5, 0.5]))


plt.figure(figsize=(8, 8))
plt.imshow(np.reshape(sino, (sx, sy)), cmap='gray', interpolation='nearest')
plt.show(block=True)

hi = 1
niter = 20
# Reconstruct object.
init = 1e-12 * np.ones((sx, sy))
rec_art = art(prb, sino, init, niter)
rec_art = rescale(np.rot90(rec_art)[::-1], hi)
plt.figure(figsize=(8, 8))
plt.imshow(rec_art, cmap='gray', interpolation='none')
plt.title('ART')

init = 1e-12 * np.ones((sx, sy))
rec_sirt = sirt(prb, sino, init, niter)
rec_sirt = rescale(np.rot90(rec_sirt)[::-1], hi)
plt.figure(figsize=(8, 8))
plt.imshow(rec_sirt, cmap='gray', interpolation='none')
plt.title('SIRT')

init = 1e-12 * np.ones((sx, sy))
rec_mlem = mlem(prb, sino, init, niter)
rec_mlem = rescale(np.rot90(rec_mlem)[::-1], hi)
plt.figure(figsize=(8, 8))
plt.imshow(rec_mlem, cmap='gray', interpolation='none')
plt.title('MLEM')
plt.show()