0

我正在使用 SART 方法研究断层扫描算法。我正在使用 2D 投影来获得 3D 体积。

我的主要问题是计算时间真的很长..

例如，对于大小为 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

我在这里看到很多嵌套的 for 循环。您应该在 MATLAB 中尽可能使用矢量化代码而不是 for 循环，以加快您的算法。也不能使用嵌套 parfor。

只要您的循环彼此独立（当然您需要安装 Parallel Computing Toolbox），就可以进行并行化。在这种情况下，使用它非常简单：

  result = zeros(1,n); % preallocate result storing array
  parfor i=1:n % parallel for loop
    %your code here, e.g.:
    result(i) = ...;
  end

可以手动设置一个 MATALB 工作池，但也可以只使用 parfor 并让 MATLAB 自动处理该池（至少在 r2014a 中）。

请注意，并行化并不总能加快速度。

于 2014-07-07T16:48:28.967 回答
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

使用此代码，重建不起作用，而不是编写 result_Ik(k) = Ik 我需要编写 Result_Ik(k+1) = Ik 以便在下一次迭代中从迭代 k-1 中获得 Ik 的值（上一次迭代）

于 2014-07-08T07:38:51.963 回答