Parallélisation d'un algorithme de Tomographie en vue d'une utilisation de GPUs
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:
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.
