Aide pour vectoriser une boucle
Bonjour !
Je debute en Python, et vectoriser des boucles for ca n'a jamais ete mon fort (mouhahaha).
Bref, je voudrais rendre mon code plus rapide en vectorisant la boucle suivante :
Code:
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
|
def Projection(t,theta,P):
x1 = P[0]
y1 = P[1]
A = P[2]
B = P[3]
alpha = np.pi*P[4]/180
rho = P[5]
s = np.sqrt(x1**2+y1**2)
gamma = np.arctan2(y1,x1)
thp = theta - alpha
tp = t - s*np.cos(gamma-theta)
a2 = A**2*(np.cos(thp))**2+B**2*(np.sin(thp))**2
if np.abs(tp)<=np.sqrt(a2):
projection = (2*rho*A*B/a2)*np.sqrt(a2-tp**2)
else:
projection = 0
return projection
N = 256 #size of the matrix
Proj = np.zeros((10,N,N))
T = np.linspace(-1,1,N)
for e in range(10): #for all the ellipse
for thet in range(N): #for all the angles
theta = np.pi*thet/N
for tt in range(N): #for all the samples
t = T[tt]
Proj[e,thet,tt] = Projection(t,theta,phantom[e,:])
Sino = Proj.sum(axis=0) |
Cette boucle n'est pas top et elle prend du temps pour de grand N.
J'ai essaye de vectoriser de cette maniere :
Code:
1
2
3
4
5
6
7
8
9
10
|
theta = np.linspace(0,math.pi,N)
t = np.linspace(-1,1,N)
T,TH = pyl.meshgrid(t,theta)
T = T.flatten(1)
TH = TH.flatten(1)
projMat = np.zeros((10,N,N))
for e in range(10):
projMat[e,:,:] = Projection(T,TH,phantom[e,:]) |
Mais le probleme c'est l'evaluation de ma condition "if" dans la definition de Projection (enfin je pense), j'ai donc modifie en remplacant tout ce qui etait dans le "if,else" par :
Code:
1
2
|
projection[(t%2) <= np.sqrt(thp)] |
Quand je lance mon script j'obtiens l'erreur suivante :
Traceback (most recent call last):
File "hw2v3.py", line 46, in <module>
projMat[e,:,:] = Projection(T,TH,phantom[e,:])
ValueError: shape mismatch: objects cannot be broadcast to a single shape
Je me doute bien que ma matrice projMat n'a pas la bonne dimension pour faire ce que je veux faire, mais je bloque un peu...
Merci pour votre aide !
