0

SART メソッドを使用して、断層撮影アルゴリズムに取り組んでいます。3D ボリュームを取得するために 2D プロジェクションを使用しています。

私の主な問題は、計算時間が非常に長いことです..

たとえば、サイズが 88*75 ピクセルのプロジェクションでは、プログラムが完了するまでに 5 日かかります。

matlab プロファイラーを使用すると、(重い) ラドン行列の乗算を含む行に非常に時間がかかっていることがわかります。

そのため、アルゴリズムを並列化したいのですが、よくわかりません。

これは、並列化したいアルゴリズムの一部です。

line which need a lot of time: Correction = mtimesx(facteur, Matrice_Radon_Transposee, 'SPEED');


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
4

2 に答える 2

0

このコードをベクトル化する方法が本当にわかりません..

すべてのコードをベクトル化する必要はありません。次の 2 行だけです。

Correction = mtimesx(facteur, Matrice_Radon_Transposee, 'SPEED'); 

Correction = Correction*Soustraction;

では、for ループを並列化するのは面白いでしょうか? :

for k = 1: ite

これを書くのは正しいですか：

result_Projection_old = zeros(1,ite);
result_Ik = zeros(1,ite);

parfor  k = 1 : ite                

    Somme = 0;      
    Projection_old = result_Projection_old(k);
    Ik = result_Ik(k);

    Soustraction = Projection2 - Projection_old;


    for z=1:LNproj
        for a=1:L_carre

            Somme = Somme + Matrice_Radon(z,a)*Matrice_Radon(z,a);

        end

        facteur = landa/Somme;
        Correction = mtimesx(facteur,Matrice_Radon_Transposee,'SPEED');    
        Correction = Correction*Soustraction;
        Ik1 = Ik + Correction;
        Ik=Ik1; 
        result_Ik(k) = Ik;                                                          

        Projection_old=Matrice_Radon*Ik;
        result_Projection_old(k) = Projection_old;                     

    end
end

このコードでは、次の反復で反復 k-1 から Ik の値を取得するために、result_Ik(k) = Ik を記述する代わりに、再構成が機能しません。Result_Ik(k+1) = Ik を記述する必要があります。 (前の繰り返し)

于 2014-07-08T07:38:51.963 に答える