<jupyter_start><jupyter_text>Quality Metrics and Reconstruction DemoDemonstrate the use of full reference metrics by comparing the reconstruction of a simulated phantom using SIRT, ART, and MLEM.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
from xdesign import *
NPIXEL = 128<jupyter_output><empty_output><jupyter_text>Generate a phantomUse one of XDesign's various pre-made and procedurally generated phantoms.<jupyter_code>np.random.seed(0)
soil_like_phantom = Softwood()<jupyter_output><empty_output><jupyter_text>Generate a figure showing the phantom and save the discretized ground truth map for later.<jupyter_code>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()<jupyter_output><empty_output><jupyter_text>Simulate data acquisitionGenerate a list of probe coordinates to simulate data acquisition for parallel beam around 180 degrees.<jupyter_code>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)<jupyter_output><empty_output><jupyter_text>Make a probe.<jupyter_code>probe = Probe(size=1/NPIXEL)<jupyter_output><empty_output><jupyter_text>Use the probe to measure the phantom.<jupyter_code>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()<jupyter_output><empty_output><jupyter_text>ReconstructReconstruct the phantom using 3 different techniques: ART, SIRT, and MLEM.<jupyter_code>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)
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()<jupyter_output><empty_output><jupyter_text>Quality MetricsCompute local quality for each reconstruction using MS-SSIM, a convolution based quality metric.<jupyter_code>quality = list()
for rec in [rec_art, rec_sirt, rec_mlem]:
    scales, mscore, qmap = msssim(discrete, rec)
    quality.append(mscore)<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>plt.figure()
plt.bar(["ART", "SIRT", "MLEM"], quality)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12701252236135957, 'fraction_numerical': 0.018187239117471676, 'mean_word_length': 7.0843373493975905, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 37, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '16038', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 1294, 'n_tokens_neox': 1199, 'n_words': 285, 'path': 'docs/source/demos/FullReferenceMetrics.ipynb', 'repo_name': 'carterbox/xdesign', 'size': '3354'}