Discussion :

[Alex3434]

Bonjour à tous,

Je viens vers vous car travaillant actuellement sur la conception d'un imageur 3D, j'ai réalisé un algorithme permettant par le biais de la méthode itérative SART d'effectuer la reconstruction tomographique d'un objet en 3D.

Le problème est que le temps de calcul est extrêmement long (~ 5 jours pour des projection en 88*75 pixels), ceci venant du fait que la matrice de Radon (matrice permettant de passer du domaine image au domaine projection et vice versa) est de grande taille.

En utilisant le profiler, j'ai pu cibler les lignes gourmandes en temps de calcul, et c'est bien la multiplication avec la matrice de Radon qui est mise en cause.

Je me suis donc renseigné sur la parallélisation en utilisant par exemple des boucles parfor au lieu de boucles for, le problème est que j'ai vraiment du mal avec la logique de la parallélisation.

Voici la partie de mon algorithme que j'aimerai paralléliser :

De la ligne 34 (inclus) à la ligne 75 (inclus)

Code MATLAB :

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
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
%% Méthode SART
 
% Dans cette partie, nous allons utiliser une méthode de correction
% itérative de type additive.
 
 
fprintf('Méthode SART : Start\n');
 
Temporaire=zeros(L,L);
Erreur=zeros(H,ite);
pos2=1;
 
 for ligne=1:H                                                              % Balayage de l'ensemble des sinogrammes
 
    Ik1=zeros(L*L,1);                                                       % Initialisation des matrices images aux itérations k+1 et k
    Ik=zeros(L*L,1);                                                        %
    Ik2=zeros(L,L);
    Projection(:,:)=double(IMG_Sino{ligne});                                % Stockage du sinogramme = Projection_mesurée
    pos=1;
 
    for incr0=1:Nproj
        for incr=1:L
 
            Projection2(pos,1)=Projection(incr,incr0);                      %Vecteur colonne proj
 
            pos=pos+1;
        end
    end
 
 
    Projection_old=zeros(L*Nproj,1);                                        % Initialisation de la matrice projection_itération_précédente                                      
    Somme = 0;
    Err = 100;
    for  k=1:ite                                                            % Boucle définie par le nombre d'itérations choisi par l'utilisateur
 
        Soustraction = Projection2 - Projection_old;
 
        Erreur(ligne,k)=mean(Soustraction);
        Err=Erreur(ligne,k);
 
        for z=1:L*Nproj
            for a=1:L*L
 
                Somme = Somme + Matrice_Radon(z,a)*Matrice_Radon(z,a);
 
            end
 
            facteur = landa/Somme;
 
            Correction = mtimesx(facteur,Matrice_Radon_Transposee,'SPEED');                       % Formule de correction additive de l'image
 
            Correction = Correction*Soustraction;
            Ik1 = Ik + Correction;
            Ik=Ik1;                                                             % Stockage de l'image corrigée dans l'image_itération_précédente
 
             Projection_old=Matrice_Radon*Ik;                          % Calcul de la projection à partir de l'image_itération_précédente
 
        end
 
 
        pos3=1;
 
        for incr2=1:L
            for incr3=1:L
 
                Ik2(incr2,incr3)=Ik1(pos3,1);                                    % Transformation vecteur en matrice
                pos3=pos3+1;
 
            end
        end
        Coupe_SIRT(ligne,:,:)=Ik2;                                              % Stockage des coupes générées
 
    end
    fprintf('Progression : %2d/%d\n',ligne,H);                                  % Affichage de l'avancement du traitement
 end

Merci par avance, et j'espère que mon code ne vous rebutera pas trop.

Cordialement,

Alexandre.

[ToTo13]

Bonjour,
tout d'abord, il faut éviter les calculs dans les tests des boucles (L*Nproj, L*L, etc.) : ce sont des valeurs fixes pré-calculable.
Est ce que c'est du MatLab ? Si oui, c'est le langage le plus lent qui soit. Re-développe donc cela en C/C++ et tu auras quelque chose qui ira au moins 10 fois plus vite, voire plus si tu parviens à utiliser des optimisations SSE2/auto-vectorisation.
Si je suis pas familier avec les opérations que tu fais (corrige moi si je me trompe) :
- les soustractions/additions de matrices se font en temps linéaires, donc inutile de paralléliser, surtout si tu optimises.
- les multiplications de matrices sont parallélisables, il suffit pour cela de découper (virtuellement) la matrice résultat en bandes et de faire faire le calcul/remplissage de chaque bande à un thread différent. Cela étant possible/facile car on ne fait que des lectures.

[Alex3434]

Bonjour et merci pour votre réponse,

Oui il s'agit bien d'un code Matlab.

Je suis bien conscient que Matlab n'est pas optimal d'un point de vue de vitesse étant donné que c'est du script et non du code compilé, cependant cela fait partie d'une première phase afin de valider des méthodes, dans l'avenir un développement en C est envisageable.

En ce qui concerne les soustractions/additions : En utilisant le profiler, j'ai pu cibler les lignes gourmandes en temps de calcul, et il est ressortit que c'était bien les opérations concernant la matrice de Radon qui prenaient le plus de temps, mais il est vrai qu'en ce qui concerne les soustractions/additions le temps de calcul est négligeable, ce n'est pas sur ce point qu'il y a besoin de paralléliser.

Pour la multiplication c'est une autre paire de manche, la matrice de Radon étant très grande, le nombre de calcul demandé l'est aussi. J'ai déjà optimisé un peu le code en utilisant un mexfile appelé mtimesx dédié à la multiplication de matrices et censé accélérer l'algo, il en résulte que j'ai pu (en plus d'autres petites optimisations) diviser le temps d'exécution par 2 (Avant il me fallait plus d'une semaine pour faire le meme calcul).

Je vais essayer d'intégrer la découpe virtuelle, mais j'ai du mal a visualiser quelles bandes des matrices je dois découper pour les multiplier entres elles.

En fait quand je parle de parallélisation, j'avais pensé paralléliser la grande boucle : for ligne=1:H en utilisant la boucle parfor au lieu d'une for, car cette boucle peut voir sa valeur H (représentant le nombre de ligne de pixels sur l'image) être très grande.

Ainsi chaque coeur effectuerait la reconstruction d'une coupe.

Ou alors si c'est trop complexe, paralléliser simplement la boucle : for k=1:ite afin que chaque coeur calcule une itération. Le problème dans ce cas vient du fait qu'on a besoin des valeurs de l'itération précédente pour effectuer le calcul, donc je pense que cela va poser des problèmes.

Edit : Voici une version que j'ai réalisé ou je parallélise le calcul le plus gourmand en temps, n'ayant que deux coeurs, cela rallonge mon temps d'exécution par rapport à la version du code non parallélisé.

Code MATLAB :

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

1
2
3
4
5
6
 Correction = mtimesx(facteur,Matrice_Radon_Transposee,'SPEED');    % Formule de correction additive de l'image
            parfor ll = 1 : L_carre
            Correction2 (ll)= Correction(ll,:)*Soustraction;
            end
            Ik1 = Ik + Correction2;
            Ik=Ik1;

Pensez vous que cette parallélisation puisse vraiment apporter un gain ?

[ToTo13]

Pour la parallélisation de la multiplication, c'est assez simple.
- disons que tu as nbCPU (ou thread utilisables)
- tu veux remplir une matrice résultat de N lignes et M colonnes.
- tu considères ta matrice résultat comme N / nbCPU (= X) bandes (0... X-1, X... 2X-1, ... N).
- chaque thread/CPU s'occupe de remplir une bande différente

Je rajoute à nouveau que ce genre de calcul est plus certainement optimisable par de l'auto-vectorisation en C/C++, ce qui ferait un gain monstrueux en performances. A froid il suffirait juste de transposer la deuxième matrice et ainsi de faire des multiplications ligne ligne et non ligne colonne.

[Alex3434]

Bon je vais programmer ca en C/C++ histoire de ne pas perdre plus de temps.

Merci beaucoup pour votre aide et pour m'avoir accordé un peu de votre temps.

Cordialement,

Alexandre

[ToTo13]

Tu peux peut être commencer par ne faire que la partie multiplication / Radon.

Surtout pense auto-vectorisation en plus de parallélisme, c'est vraiment important.

[Alex3434]

Oui c'est exactement ce que j'ai pensé faire, histoire de ne pas avoir fait mon code matlab pour rien, et puis comme seul le calcul avec la matrice de Radon est long, autant juste calculer ça en C.

[Jean Dumoncel]

Bonjour,

même si matlab n'est pas réputé pour être une bête de course, quelques optimisations apportent parfois de bonnes surprises. Et ton code comporte quelques éléments à améliorer, en voici deux qui profitent du fait que matlab, c'est du calcul matriciel :
Tu peux changer ceci :

Code MATLAB :

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

1
2
3
4
5
6
7
8
for incr0=1:Nproj
    for incr=1:L
 
        Projection2(pos,1)=Projection(incr,incr0);                      %Vecteur colonne proj
 
        pos=pos+1;
    end
end

Par :

Code MATLAB :

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

Projection2=Projection(:);

De même, tu dois pouvoir supprimer ces boucles :

Code MATLAB :

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

1
2
3
4
5
6
7
8
for incr2=1:L
    for incr3=1:L
 
        Ik2(incr2,incr3)=Ik1(pos3,1);                                    % Transformation vecteur en matrice
        pos3=pos3+1;
 
    end
end

en utilisant reshape à la place.


Et à la place de :

Code MATLAB :

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

1
2
3
4
5
for a=1:L*L
 
    Somme = Somme + Matrice_Radon(z,a)*Matrice_Radon(z,a);
 
end

Tu peux écrire :

Code MATLAB :

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

Somme = Somme + sum(Matrice_Radon(z,:).*Matrice_Radon(z,:));

[Alex3434]

Merci pour les petites astuces, pour de petites images ca a bien augmenté la vitesse mais pour de plus grande c'est anecdotique..

N'ayant pas d'autres solutions pour accélérer les multiplications avec la matrice de Radon, je continue à développer en C.

[Alex3434]

Bonjour,


J'ai trouvé une petite astuce qui m'a permis de passer de 5 jours d'exécution à 1 minute ! En fait la matrice de Radon comportant beaucoups de 0 (et donc des calculs inutiles), en définissant la matrice comme une matrice sparse, on ne garde que les cases comportant une valeur différente de 0 sans perdre l'info sur son emplacement dans la matrice.

Ainsi le nombre de multiplications se voit réduit de façon importante.

En espérant que cette astuce puisse servir à quelqu'un !