<jupyter_start><jupyter_text>Tomography reconstruction$\newcommand{\n}[1]{\left\|1 \right\|}$ $\renewcommand{\a}{\alpha}             $ $\renewcommand{\b}{\beta}              $ $\renewcommand{\c}{\gamma}             $ $\renewcommand{\d}{\delta}             $ $\newcommand{\D}{\Delta}               $ $\newcommand{\la}{\lambda}             $ $\renewcommand{\t}{\tau}               $ $\newcommand{\s}{\sigma}               $ $\newcommand{\e}{\varepsilon}          $ $\renewcommand{\th}{\theta}            $ $\newcommand{\x}{\bar x}               $ $\newcommand{\R}{\mathbb R}            $ $\newcommand{\N}{\mathbb N}            $ $\newcommand{\Z}{\mathbb Z}            $ $\newcommand{\E}{\mathcal E}           $ $\newcommand{\lr}[1]{\left\langle 1\right\rangle}$$\newcommand{\nf}[1]{\nabla f(1)}     $$\newcommand{\hx}{\hat x}               $$\newcommand{\hy}{\hat y}               $$\DeclareMathOperator{\prox}{prox}      $$\DeclareMathOperator{\argmin}{argmin}  $$\DeclareMathOperator{\dom}{dom}        $$\DeclareMathOperator{\id}{Id}          $$\DeclareMathOperator{\conv}{conv}      $We want to solve  $Ax = b$,where $A$ is a matrix, obtained from the projection tomography operator and $b$ is the observed sinogram.<jupyter_code>import scipy.sparse.linalg as spr_LA
import matplotlib as mpl
from skimage import data, transform, img_as_float, transform
from skimage.color import rgb2gray
from tomo_utils import generate_synthetic_data, build_projection_operator
from fixed_points import *
import matplotlib.pyplot as plt
import seaborn as sns

%load_ext autoreload
%autoreload 2
%matplotlib inline

sns.set()<jupyter_output>The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload<jupyter_text>Generate the data<jupyter_code># number of positions
n_pos = 256
# number or angles
n_angles = 128
#img = generate_synthetic_data(p)

x = rgb2gray(data.load('phantom.png'))
img = transform.resize(x, (n_pos, n_pos))

w = img.ravel()
A = build_projection_operator(n_pos, n_angles)
#A = build_projection_operator(n_pos, n_dir=n_angles, l_det=64)


n = n_pos * n_pos
m = n_pos * n_angles

x_true = w

# no noise
#b = A.dot(x_true)

# with noise
b = A.dot(x_true) + np.random.randn(m)

# starting point
x0 = np.zeros(n)

# define operator T:
norms = spr_LA.norm(A, axis=1)
T = lambda x: x - 1./m * A.T.dot((A.dot(x)-b)/norms**2)

J = lambda x: LA.norm(T(x)-x)


N = 1000

ans1 = krasn_mann(T, x0, 0, numb_iter=N)
ans2 = fixed_point_agraal(T, x0, numb_iter=N, phi=1.5, output=False)

x1 = ans1[1]
x2 = ans2[1]    
print("Fixed point residuals. KM and aGRAAL:", J(x1), J(x2))<jupyter_output>Fixed point residuals. KM and aGRAAL: 0.01206923787295626 4.7282722448181315e-05<jupyter_text>Show the results<jupyter_code>plt.plot(ans1[0], '--b', label="KM: $x^{k+1}=Tx^k$")

plt.plot(ans2[0], '#FFD700', label="aGRAAL")

plt.yscale('log')
plt.legend()


plt.xlabel(u'iterations, $k$')
plt.ylabel('residual')

plt.legend()
#plt.grid()
plt.savefig('figures/tomo-12-darkgrid.pdf',bbox_inches='tight')
plt.show()
plt.clf()<jupyter_output><empty_output><jupyter_text>What is the stepsize in aGRAAL?<jupyter_code>plt.plot(ans2[2], '.', color='#FFD700', label="aGRAAL")
plt.legend()


plt.xlabel(u'iterations, $k$')
plt.ylabel('stepsize $\lambda_k$')

plt.legend()
#plt.grid()
plt.savefig('figures/tomo-22-darkdrid.pdf',bbox_inches='tight')
plt.show()
plt.clf()<jupyter_output><empty_output><jupyter_text>Show the original image and reconstructed ones<jupyter_code>print("Original image and reconstructed ones")
img1 = x1.reshape(n_pos, n_pos)
img2 = x2.reshape(n_pos, n_pos)
fig, ax, = plt.subplots(nrows=1, ncols=3, figsize=(15, 5))
ax[0].imshow(img, cmap='gray')
ax[0].set_title("True")
ax[1].imshow(img1, cmap='gray')
ax[1].set_title("KM")
ax[2].imshow(img2, cmap='gray')
ax[2].set_title("aGRAAL")
plt.show(fig)<jupyter_output>Original image and reconstructed ones
{'fraction_non_alphanumeric': 0.17320598864223025, 'fraction_numerical': 0.028652555498193084, 'mean_word_length': 4.152925531914893, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 18, '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': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '619023', 'license': 'MIT', 'n_tokens_mistral': 1547, 'n_tokens_neox': 1443, 'n_words': 278, 'path': 'feasibility_problem/Tomography_reconstruction.ipynb', 'repo_name': 'ymalitsky/graal', 'size': '3874'}