-
Notifications
You must be signed in to change notification settings - Fork 16
/
radon.pyx
120 lines (85 loc) · 3.05 KB
/
radon.pyx
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
# -*- coding: utf-8 -*-
import numpy as np
cimport numpy as np
from scipy.sparse import csr_matrix
cimport cradon
try: range=xrange
except: pass
def radon2d(data, theta):
if np.min(theta) < 0.0 or np.max(theta) >= np.pi:
raise ValueError('theta should be within [0 pi)')
nx, ny = data.shape
if nx != ny:
raise RuntimeError('data should be a square array')
nw = nx
sinogram = np.zeros((nw,theta.size))
Tx = np.zeros((nw,2))
Rx = np.zeros((nw,2))
cdef size_t nTx = Tx.shape[0]
for nt in range(theta.size):
if theta[nt] == 0.0:
for n in range(nw):
Tx[n,0] = n
Tx[n,1] = 0
Rx[n,0] = n
Rx[n,1] = ny-1
elif np.abs(theta[nt] - np.pi/2) < 1.e-10:
for n in range(nw):
Tx[n,0] = 0
Tx[n,1] = n
Rx[n,0] = nx-1
Rx[n,1] = n
elif theta[nt] < np.pi/2:
xs = nw/2 * (1. - np.cos(theta[nt]) )
ys = nw/2 * (1. - np.sin(theta[nt]) )
for n in range(nw):
xk = xs + n * np.cos(theta[nt])
yk = ys + n * np.sin(theta[nt])
x0 = 0.0
y0 = yk + (xk-x0) / np.tan(theta[nt])
if y0 > ny-1:
y0 = ny-1
x0 = xk - (y0-yk) * np.tan(theta[nt])
x1 = nx-1
y1 = yk - (x1-xk) / np.tan(theta[nt])
if y1 < 0.0:
y1 = 0.0
x1 = xk + (yk-y1) * np.tan(theta[nt])
Tx[n,0] = x0
Tx[n,1] = y0
Rx[n,0] = x1
Rx[n,1] = y1
else:
xs = nw/2 * (1. - np.cos(theta[nt]) )
ys = nw/2 * (1. - np.sin(theta[nt]) )
for n in range(nw):
xk = xs + n * np.cos(theta[nt])
yk = ys + n * np.sin(theta[nt])
x0 = nx-1
y0 = yk + (x0-xk) * np.tan(theta[nt]-np.pi/2.0)
if y0 > ny-1:
y0 = ny-1
x0 = xk + (y0-yk) / np.tan(theta[nt]-np.pi/2.0)
x1 = 0.0
y1 = yk - (xk-x1) * np.tan(theta[nt]-np.pi/2.0)
if y1 < 0.0:
y1 = 0.0
x1 = xk - (yk-y1) / np.tan(theta[nt]-np.pi/2.0)
Tx[n,0] = x0
Tx[n,1] = y0
Rx[n,0] = x1
Rx[n,1] = y1
if np.max(Rx) > nx-1:
print('Rx '+str(n))
print(Rx)
if np.max(Tx) > nx-1:
print('Tx '+str(n))
print(Tx)
Ldata = ([0.0], [0.0], [0.0])
cradon.radon2d(<double*> np.PyArray_DATA(Tx), <double*> np.PyArray_DATA(Rx), nTx, nx, ny, Ldata)
M = nTx
N = nx * ny
L = csr_matrix(Ldata, shape=(M,N))
p = L.dot(data.flatten())
sinogram[:,nt] = p
return sinogram