#!/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))

[((7, 6, 7, 42), 'tomograpy.centered_cubic_map', 'tomograpy.centered_cubic_map', ({(7, 35, 7, 37): '10', (7, 39, 7, 41): '64'}, {}), '(10, 64)', False, 'import tomograpy\n'), ((8, 9, 8, 49), 'tomograpy.phantom.shepp_logan', 'tomograpy.phantom.shepp_logan', ({(8, 39, 8, 48): 'obj.shape'}, {}), '(obj.shape)', False, 'import tomograpy\n'), ((11, 4, 11, 30), 'tomograpy.fov', 'tomograpy.fov', ({(11, 18, 11, 21): 'obj', (11, 23, 11, 29): 'radius'}, {}), '(obj, radius)', False, 'import tomograpy\n'), ((12, 7, 12, 82), 'tomograpy.centered_stack', 'tomograpy.centered_stack', (), '', False, 'import tomograpy\n'), ((15, 28, 15, 79), 'tomograpy.models.thomson', 'tomograpy.models.thomson', (), '', False, 'import tomograpy\n'), ((17, 4, 17, 15), 'time.time', 'time.time', ({}, {}), '()', False, 'import time\n'), ((22, 4, 22, 15), 'time.time', 'time.time', ({}, {}), '()', False, 'import time\n'), ((27, 4, 27, 15), 'time.time', 'time.time', ({}, {}), '()', False, 'import time\n'), ((28, 40, 28, 50), 'numpy.ones', 'np.ones', ({(28, 48, 28, 49): '3'}, {}), '(3)', True, 'import numpy as np\n'), ((19, 33, 19, 44), 'time.time', 'time.time', ({}, {}), '()', False, 'import time\n'), ((25, 37, 25, 48), 'time.time', 'time.time', ({}, {}), '()', False, 'import time\n'), ((30, 32, 30, 43), 'time.time', 'time.time', ({}, {}), '()', False, 'import time\n')]