Discussion :

[lammouch]

Salut,
voici un code que j'ai fait, en fait il permet de remplir une matrice 3D, element par element , il s'execute en 3 min pour une matrice de 9*9*9mais si j'ai besoin de travailler avec une matrice de 200*200*200 par exemple, ca va prendre enormement de temps!!
je crois qu'il faut trouver un moyen de remplir la matrice d'un facon plus rapide en evitant les boucles...mais j'arrives pas...

voici la partie du code que je veux modifier:

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
 
V=zeros(fix(d),fix(d),fix(d));
F=pos_voxel1(d,u); 
fprintf('%.0f\n',d);
dprime=power(d,3);
MIN=min(F(:,2)); 
nb=1000; 
tt=1;
 
%%%%%%%%%%% calcul de l'intensité de chaque voxel %%%%%%%%%%%%%%%%%%%%%
for i=1:1:dprime
   X=d+(MIN-F(i,2))/u;
   Y=d+(MIN+F(i,1))/u;
   Z=d+(MIN+F(i,3))/u;
   for j=1:1:nprime
      S=sqrt((F(i,1)-C(j,1))^2+(F(i,2)-C(j,2))^2+(F(i,3)-C(j,3))^2);
      Sprime=S./a;
      if (i==nb*tt) 
         fprintf('%.0f %% de blob2voxel ont été exécutés en %.2f secondes, soit %.2f minutes\n',i*100/dprime,etime(clock,t1),etime(clock,t1)/60) 
         tt=tt+1;
      end
      if real(S)<=c 
         h=8.63*sqrt(1-power(Sprime,2));
         t=power((1-Sprime.^2),M);
         b=abs((real(besseli(m,h))./real((besseli(m,alpha)))).*real(t));
         V(X,Y,Z)= V(X,Y,Z)+ Q(j,1)*b;
      else V(X,Y,Z)=V(X,Y,Z);
      end
   end
end

est ce que vous pouvez me trouver une solution pour ameliorer la vitesse d'execution de mon programme??

merci=)
Lama

[tonylataupe]

pourrais-tu donner les valeurs de u,d que tu utilise ainsi que la fonction pos_voxel1 que l'on puisse exécuter ton code. Ca sera plus facile pour l'optimiser.
Sinon, je pense que pour gagner du temps, il faut calculer directement toute ta matrice (ou chaque matrice 2D), au lieu de remplir chaque élément. Je ne sais pas si c'est possible, je n'ai pas trop compris les calculs effectués.

[FR119492]

Tiens !
Une matrice à 3 dimensions. C'est nouveau, ça !
Jean-Marc Blanc

[mr_samurai]

Salut,
Voici un début de piste d'optimisation :

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
 
Tmp = d + (MIN-F)/u;
for i=1:dprime
    S = sqrt(sum((repmat(F(i,:),nprime,1)-C).^2));
    Sprime=S./a;
    if (i==nb*tt)
        fprintf('%.0f %% de blob2voTmp(i,2)el ont été eTmp(i,2)écutés en %.2f secondes, soit %.2f minutes\n',i*100/dprime,etime(clock,t1),etime(clock,t1)/60)
        tt=tt+1;
    end
 
    [J,I] = find(real(S)<=c);
 
    if % ....
        % ....
    end
end

bonne chance , +++

[lammouch]

je vous donne le code en entier, peut etre ca sera plus facile..

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
33
34
function blob2voxel(nom,a,n,g,m,alpha,u)
 
dz=g*(sqrt(3))/2; 
M=m/2;
nprime=power(n,3);
c=a+u;
MT=ones(n,n,n);
Q=mat2vect(MT)
C=pos_blob(g,n,dz);
 
d=21;
V=zeros(fix(d),fix(d),fix(d));
F=pos_voxel1(d,u); 
fprintf('%.0f\n',d);
dprime=power(d,3);
MIN=min(F(:,2)); 
 
for i=1:1:dprime
    X=d+(MIN-F(i,2))/u;
    Y=d+(MIN+F(i,1))/u;
    Z=d+(MIN+F(i,3))/u;
    for j=1:1:nprime
        S=sqrt((F(i,1)-C(j,1))^2+(F(i,2)-C(j,2))^2+(F(i,3)-C(j,3))^2);
        Sprime=S./a;
        if real(S)<=c 
            h=8.63*sqrt(1-power(Sprime,2));
            t=power((1-Sprime.^2),M);
            b=abs((real(besseli(m,h))./real((besseli(m,alpha)))).*real(t));
            V(X,Y,Z)= V(X,Y,Z)+ Q(j,1)*b;
        else
            V(X,Y,Z)=V(X,Y,Z);
        end
    end
end

le code de la foction posvoxel1 est:

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
function A=pos_voxel1(d,u)
F=zeros(power(d,3),3);
 
k=1;
 
for c=1:1:d 
    for j=1:1:d
        for i=1:1:d
            l=-1*((j-0.5)*u-(u*d)./2);
            s=1*((i-0.5)*u-(u*d)./2); 
            w=1*((c-0.5)*u-(u*d)./2);
            F(k,:)=[s l w];
            k=k+1;
        end
    end
end
 
A=F;

et le code de la fonction pos_blob(g,n,dz) est:

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
function B=pos_blob(g,n,dz)
C=zeros(power(n,3),3);
 
k=1;
 
for i =1:1:n
    if mod(i,2)~=0
        for f=-fix(n/2):1:fix(n/2)
            for c=-fix(n/2):1:fix(n/2)
            x=-g*f ;
            y=-g*c ;
            z=(n-1)*(dz/2)-(i-1)*(dz);
            C(k,:)=[x y z];
            k=k+1;
            end
        end
    else 
        for f=-fix(n/2):1:fix(n/2)
            for c=-fix(n/2):1:fix(n/2)
                x=-1*(g*f+g/2) ;
                y= -1*(g*c+g/2) ;
                z=(n-1)*(dz/2)-(i-1)*(dz);
                C(k,:)=[x y z];
                k=k+1;
            end
        end
    end
end
C;

je ne sais pas si l'utilisation des matrices sparses peut m'apporter une aide...qu'estce que vous en pensez?
merci :)

[lammouch]

salut a nouveau, et désolée pour n'avoir pas ajouté les balises au code...

en fait je n'ai pas tres bien compris a quoi sert la comande suivante:
[J,I] = find(real(S)<=c)
ca remplace la condition if??
et est ce que vous savez une methode qui me permet de remplir une matrice d'un seul coup?

Merci,
Lama

[mr_samurai]

find permet d'avoir les indices qui vérifient une condition donnée. pratique pour remplacer les test 'IF'

Exemple :

Code :

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

1
2
3
4
5
6
7
A = [1 2 3 4 5; 6 7 8 9 10];
 
% Indices de A tel que A(i,j)>8
[I,J] = find(A>8)
 
% Les valeurs : 
Aval = A(sub2ind(size(A),I,J))

L'execution donne :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
I =
     2
     2
 
J =
     4
     5
 
Aval =
     9
    10

Plus d'informations : help find , help sub2ind

++