#!/usr/bin/env python
import time
import numpy as np
import tomograpy
import lo
# object
obj = tomograpy.centered_cubic_map(10, 64)
obj[:] = tomograpy.phantom.shepp_logan(obj.shape)
# data 
radius = 200
a = tomograpy.fov(obj, radius)
data = tomograpy.centered_stack(a, 128, n_images=60, radius=radius, max_lon=np.pi)
# model
kwargs = {"pb":"pb", "obj_rmin":1.5, "data_rmin":1.5}
P, D, obj_mask, data_mask = tomograpy.models.thomson(data, obj, u=.5, **kwargs)
# projection
t = time.time()
data[:] = (P * obj.ravel()).reshape(data.shape)
print("projection time : " + str(time.time() - t))
# data
# backprojection
t = time.time()
x0 = P.T * data.ravel()
bpj = x0.reshape(obj.shape)
print("backprojection time : " + str(time.time() - t))
# inversion using scipy.sparse.linalg
t = time.time()
sol = lo.acg(P, data.ravel(), D, 1e-3 * np.ones(3),  maxiter=100, tol=1e-8)
sol = sol.reshape(obj.shape)
print("inversion time : " + str(time.time() - t))

[((95, 131), 'tomograpy.centered_cubic_map', 'tomograpy.centered_cubic_map', (['(10)', '(64)'], {}), '(10, 64)\n', (123, 131), False, 'import tomograpy\n'), ((141, 181), 'tomograpy.phantom.shepp_logan', 'tomograpy.phantom.shepp_logan', (['obj.shape'], {}), '(obj.shape)\n', (170, 181), False, 'import tomograpy\n'), ((207, 233), 'tomograpy.fov', 'tomograpy.fov', (['obj', 'radius'], {}), '(obj, radius)\n', (220, 233), False, 'import tomograpy\n'), ((241, 316), 'tomograpy.centered_stack', 'tomograpy.centered_stack', (['a', '(128)'], {'n_images': '(60)', 'radius': 'radius', 'max_lon': 'np.pi'}), '(a, 128, n_images=60, radius=radius, max_lon=np.pi)\n', (265, 316), False, 'import tomograpy\n'), ((407, 459), 'tomograpy.models.thomson', 'tomograpy.models.thomson', (['data', 'obj'], {'u': '(0.5)'}), '(data, obj, u=0.5, **kwargs)\n', (431, 459), False, 'import tomograpy\n'), ((476, 487), 'time.time', 'time.time', ([], {}), '()\n', (485, 487), False, 'import time\n'), ((615, 626), 'time.time', 'time.time', ([], {}), '()\n', (624, 626), False, 'import time\n'), ((776, 787), 'time.time', 'time.time', ([], {}), '()\n', (785, 787), False, 'import time\n'), ((828, 838), 'numpy.ones', 'np.ones', (['(3)'], {}), '(3)\n', (835, 838), True, 'import numpy as np\n'), ((569, 580), 'time.time', 'time.time', ([], {}), '()\n', (578, 580), False, 'import time\n'), ((716, 727), 'time.time', 'time.time', ([], {}), '()\n', (725, 727), False, 'import time\n'), ((925, 936), 'time.time', 'time.time', ([], {}), '()\n', (934, 936), False, 'import time\n')]