Parcours d'un cube de donnees

Version imprimable

Bonjour tout le monde,

Je vous expose mon probleme :

Je possede deux cubes de donnees, par exemple de taille 1000x500x500 et je fait parcourir ces gros cubes par un petit cube 10x8x8 par exemple et pour chaque petit cube, je fait ensuite un traitement via les donnees recueillies entre le donnes_1 et le donnes_2.

Mon code ressemple a plus ou moins a cela :
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 % taille des inputs et de l'output, donc des cubes donnees_1 et donnes_2 nz=1000; nx=500; ny=500; % taille du petit cube wz=10; wx=8; wy=8; for i=1:nz for j=1:nx for k=1:ny output(i+floor(wz/2),j+floor(wx/2),k+floor(wy/2))=... treatment(donnees_1(i:i+wz-1,j:j+wx-1,k:k+wy-1),... donnees_2(i:i+wz-1,j:j+wx-1,k:k+wy-1))); end end end
En faisant cela, le temps est..... long.... savez-vous s'il existe une astuce pour faire un cube glissant sans boucle for ?

Merci beaucoup de votre aide.

14/03/2013, 15h10
FLB

Salut,
tout dépend de ce que tu fais dans ton traitement?
Si tu ne peux pas vectoriser les calculs, tu peux regarder du côté des mex functions, qui permettent de coder en c, et donc d'accélérer les calculs en boucles.
Une autre solution est d'utiliser la parallele toolbox, et de transformer les "for" en "parfor".
14/03/2013, 20h28
Jerome Briot

Il faudrait dans un premier temps préallouer la mémoire pour la variable output :

Code:

output = zeros(nz,nx,ny);

Voir la FAQ : Qu'est-ce que la préallocation de mémoire ?

Bonjour a tous et merci pour ces reponses,

J'alloue deja mon output avant la boucle, le code que j'avais donne etait un exemple de mon code, voila le vrai :
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 G=zeros(nz,nx,ny); for i=1:nz-WL+1 for j=1:nx-WX+1 for k=1:ny-WY+1 T1=ones(WL,WX,WY);T2=ones(WL,WX,WY);T3=ones(WL,WX,WY); waitbar((i*j*k)/((nz-WL+1)*(nx-WX+1)*(ny-WY+1)),h,['Processing ... ' num2str(((i*j*k)/((nz-WL+1)*(nx-WX+1)*(ny-WY+1)))*100) '%']); num=sum(sum(sum( ((X(i:i+WL-1,j:j+WX-1,k:k+WY-1)-(T1*mean(mean(mean(X(i:i+WL-1,j:j+WX-1,k:k+WY-1)))))) .* ... (Y(i:i+WL-1,j:j+WX-1,k:k+WY-1)-(T2*mean(mean(mean(Y(i:i+WL-1,j:j+WX-1,k:k+WY-1)))))))))); den=sum(sum(sum((X(i:i+WL-1,j:j+WX-1,k:k+WY-1)-(T3*mean(mean(mean(X(i:i+WL-1,j:j+WX-1,k:k+WY-1)))))))).^2); G(i+floor(WL/2),j+floor(WX/2),k+floor(WY/2))=num./den; end end end
En gros a chaque boucle j'extrait les donnees de mes deux gros cubes input (X et Y) le petit cube et je fait une regression lineaire grace aux points contenues dans ces deux petits cubes et je sauve dans G (output) la pente de la droite ainsi obtenue.
J'ai essaye avec un parfor sur la premiere boucle (Du coup en mettant G(i,... au lieu de G(i+floor....))) mais le temps de calcul est toujours aussi long j'ai l'impression.

Par rapport aux mex, j'en ai fait un en Fortran qui fonctionne (il marche tres vite pour de petits cubes), cependant pour des grands cube, la fonction est lente, mais pas pour realiser le calcul (j'ai mis un mexPrintf dans les boucles pour voir les etapes du calcul) mais l'etape qui est tres longue est le lancement de la fonction.
C'est a dire que des que je vois mon premier mexPrintf dans ma premiere boucle le calcul est tres rapide, mais ce premier mex est tres long a arriver, et ne realisant que ces boucles dans la fonction, j'en deduis que ce qui prend le plus de temps est le passage des arguments de matlab a Fortran. Bizare non? oO?

Merci beaucoup

Il faut absolument éviter d'utiliser waitbar.
Fais une recherche sur ce forum pour savoir pourquoi !

Sinon, il y a quelques optimisation possible dans le code :

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
G = zeros(nz,nx,ny);
 
T1 = ones(WL,WX,WY);
T2 = ones(WL,WX,WY);
T3 = ones(WL,WX,WY);
 
wl = floor(WL/2);
wx = floor(WX/2);
wy = floor(WY/2);
 
for i = 1:nz-WL+1
    for j = 1:nx-WX+1
        for k = 1:ny-WY+1            
 
            num = sum(sum(sum(((X(i:i+WL-1,j:j+WX-1,k:k+WY-1)-(T1*mean(mean(mean(X(i:i+WL-1,j:j+WX-1,k:k+WY-1)))))) .* ...
            (Y(i:i+WL-1,j:j+WX-1,k:k+WY-1)-(T2*mean(mean(mean(Y(i:i+WL-1,j:j+WX-1,k:k+WY-1))))))))));
 
            den = sum(sum(sum((X(i:i+WL-1,j:j+WX-1,k:k+WY-1)-(T3*mean(mean(mean(X(i:i+WL-1,j:j+WX-1,k:k+WY-1)))))))).^2);
 
            G(i+wl,j+wx,k+wy)=num./den;
 
        end
    end
end

Je regarde le problème de plus près...

Les matrices T1, T2 et T3 ne servent à rien :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
G = zeros(nz,nx,ny);
 
wl = floor(WL/2);
wx = floor(WX/2);
wy = floor(WY/2);
 
for i = 1:nz-WL+1
    for j = 1:nx-WX+1
        for k = 1:ny-WY+1            
 
            num = sum(sum(sum(((X(i:i+WL-1,j:j+WX-1,k:k+WY-1)-mean(mean(mean(X(i:i+WL-1,j:j+WX-1,k:k+WY-1))))) .* ...
            (Y(i:i+WL-1,j:j+WX-1,k:k+WY-1)-mean(mean(mean(Y(i:i+WL-1,j:j+WX-1,k:k+WY-1)))))))));
 
            den = sum(sum(sum((X(i:i+WL-1,j:j+WX-1,k:k+WY-1)-mean(mean(mean(X(i:i+WL-1,j:j+WX-1,k:k+WY-1))))))).^2);
 
            G(i+wl,j+wx,k+wy)=num./den;
 
        end
    end
end

Voici une version améliorée :
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 for i = 1:nz-WL+1 for j = 1:nx-WX+1 for k = 1:ny-WY+1 x = X(i:i+WL-1,j:j+WX-1,k:k+WY-1); y = Y(i:i+WL-1,j:j+WX-1,k:k+WY-1); num = sum( ( x(:) - mean( x(:) ) ) .* ( y(:) - mean( y(:) ) ) ) ; den = sum( sum( sum( x - mean( x(:) ) ) ) .^2); G(i+wl,j+wx,k+wy) = num./den; end end end
Comme les opérations élémentaires ne se font pas dans le
même ordre que dans la solution précédente, il se peut que les deux codes ne donne pas exactement le même résultat

18/03/2013, 17h46
Glouglou69

Bonjour,

Merci de votre reponse, en effet le code est bien accelere avec ca, mais il prend toujours 2 heures au lieux de 7...
Alors que par exemple, si je prend une fenetre verticale, donc WZ=15, par exemple et WX=WY=1, ainsi je supprime les deux autres boucles, le calcul prend 10 minutes grace a la vectorisation.

Du coup, j'ai deux questions :
- Est-ce possible dans ce cas d'eviter les boucles ?
- Pourquoi lors de l'appel de mon mex en Fortran matlab met tant de temps a transmettre les inputs ?

Merci beaucoup en tout cas.