Discussion :

[Glouglou69]

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 :

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
 
% 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.

[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".

[Jerome Briot]

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

Code :

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

output = zeros(nz,nx,ny);

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

[Glouglou69]

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 :

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
 
 
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

[Jerome Briot]

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 :

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
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...

[Jerome Briot]

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

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
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

[Jerome Briot]

Voici une version améliorée :

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
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

[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.