<jupyter_start><jupyter_text>Error Analysis<jupyter_code>import numpy as np
import matplotlib.pylab as plt
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits import mplot3d

from skimage.transform import iradon, iradon_sart

from scipy.ndimage import sobel
from scipy.ndimage import binary_dilation

%matplotlib inline
# import of the phantom library which contain the phantom class
import exact_sinogram as esin
def normalize(x):
    return (x-x.min())/(x.max() - x.min())
def gradient(A):
    return np.sqrt(sobel(A, axis=0)**2 + sobel(A, axis=-1)**2)
def ndilate(A,n):
    for i in range(n):
        A = binary_dilation(A)
    return A
def error(X,Targ,p = np.inf):
    Mask = ndilate(gradient(P)>.5,2)
    Error = X-Targ
    Error[Mask==1] = 0
    Error[P==0] = 0
    return np.linalg.norm(Error,p)<jupyter_output><empty_output><jupyter_text>Creation of the Phantom<jupyter_code>circle = True
phantom_type = 'modified_shepp_logan'
#  phantom_type : 'ellipses' (or 'shepp_logan') , 'modified_shepp_logan', 'squares', 'rectangles'
n_points = 300;   # number of pixels

Phm = esin.Phantom(phantom_type = phantom_type, circle = circle)    

P = Phm.get_phantom(N = n_points)<jupyter_output><empty_output><jupyter_text>Calculation of the Sinograms<jupyter_code># Creation of the vector theta of the angles on which the Radon is calculated
theta_vec_deg = np.linspace(0, 359, 360)
theta_vec_rad = np.deg2rad(theta_vec_deg)
analytical_sinogram = Phm.get_sinogram(N = n_points, theta_vec = theta_vec_rad)

# Calculation of the Sinogram with the Aproximated formula
sinogram = esin.my_radon(P, theta=theta_vec_deg)<jupyter_output><empty_output><jupyter_text>Backprojection with the iradon function of the library Skimage<jupyter_code>P_an = normalize(iradon(analytical_sinogram, theta = theta_vec_deg, circle=circle))

P_inv = normalize(iradon(sinogram, theta= theta_vec_deg , circle=circle))<jupyter_output><empty_output><jupyter_text>Analysis of the results<jupyter_code>error(P_an, P,2)
error(P_inv, P,2)
Mask = ndilate(gradient(P)>.1,4)
##############################################################################
# Comparison of the recalculated Phantom images
deltaP = np.abs(P-P_an)
deltaP[Mask] = 0
deltaP[P==0] = 0
plt.imshow(deltaP, cmap = cm.gray)
plt.title('|P - P_an|')
#errors:
print('max(|P-P_an|) = ', np.max(deltaP))
##############################################################################
##############################################################################
# Comparison of the recalculated Phantom images
deltaP = np.abs(P-P_inv)
deltaP[Mask] = 0
deltaP[P==0] = 0
plt.imshow(deltaP, cmap = cm.gray)
plt.title('|P - P_inv|')
#errors:
print('max(|P-P_inv|) = ', np.max(deltaP))
##############################################################################
##############################################################################
# 3D Plot of the error of the recalculated Phantom with the approximated radon/iradon
X, Y = np.meshgrid(np.arange(300), np.arange(300))
Z = np.abs(P-P_inv)
Z[Mask] = 0
Z[P==0] = 0
fig = plt.figure();
ax = plt.axes(projection='3d')
ax.contour3D(X, Y, Z, 50, cmap='binary')
ax.set_xlabel('x');
ax.set_ylabel('y');
ax.set_zlabel('z');
plt.title('|P-P_inv|')
plt.show()
##############################################################################
##############################################################################
# 3D Plot of the error of the recalculated Phantom with the analytical radon/iradon
X, Y = np.meshgrid(np.arange(300), np.arange(300))
Z = np.abs(P_an-P)
Z[Mask] = 0
Z[P==0] = 0
plt.figure();
ax = plt.axes(projection='3d')
ax.contour3D(X, Y, Z, 50, cmap='binary')
ax.set_xlabel('x');
ax.set_ylabel('y');
ax.set_zlabel('z');
plt.title('|P-P_an|')
plt.show()
##############################################################################<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.2818134981968058, 'fraction_numerical': 0.015455950540958269, 'mean_word_length': 6.956967213114754, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 23, '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': '268307', 'license': 'MIT', 'n_tokens_mistral': 1463, 'n_tokens_neox': 1345, 'n_words': 289, 'path': 'Error Analysis.ipynb', 'repo_name': 'francescat93/Exact_sinogram', 'size': '3882'}