Discussion :

[Hew]

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 !

[josmiley]

tu peux replacer

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
 if np.abs(tp)<=np.sqrt(a2):
        projection = (2*rho*A*B/a2)*np.sqrt(a2-tp**2)
    else:
        projection = 0

par

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

projection = (2*rho*A*B/a2)*np.sqrt(a2-tp**2)*(np.abs(tp)<=np.sqrt(a2))

[Hew]

Merci ! J'ai fini par y arriver

[jlg_47]

Merci ! J'ai fini par y arriver

Bonjour,

Est ce que tu pourrais, stp, mettre ton bout de code?
Ce serait sympa, cela m'intéresse de voir comment tu as fait.

Je débute avec Numpy et moi aussi vectoriser les boucle FOR n'est pas mon fort!!

Bonne journée