How to accelerate a function ?

Hello !

I'm currently working on a segmentation problem, where I have to recreate the skeleton of a blood vessel network. I currently have two skeleton, one without any false branches, but with a lot of missing connections, and an other with all the actual connections and then some, which are false. My idea to have a correct AND complete skeleton was to link the loose end on the first skeleton using the second one.

In order to do so, I use this function :

 function skel = Relie(skelsec,squeltouffu,graph,lcon)
% skelsec is the correct but incomplete skeleton, 
% squeltouffu is the overdone skeleton with false branches, 
% graph is the connectivity graph of the overdone squeleton, where each voxel (3D pixel) is connected to his 
% neighbours and represented by it's linear indice
% lcon is the list of all the loose end in the correct-yet-incomplete skeleton
disp('Calcul du nombre de fins a relier')
skel = skelsec ;
skel(:) = 0 ;
contouf = bwconncomp(squeltouffu) ;
i = 1 ;
u = lcon.' ;
while (i<=length(lcon(1,:)))&&(length(lcon)>1)
    tic
    % premiere fin de vaisseau
    z1 = lcon(1,i) ;
    x1 = lcon(2,i) ;
    y1 = lcon(3,i) ;
    ind1 = sub2ind(size(skelsec),z1,x1,y1) ;
    % On la retire de la liste des fins de vaisseau restantes
    lcon2 = setdiff(u,[z1,x1,y1],'rows') ;
    lcon2 = lcon2.' ;
    % On compare les fins
    for j = 1:length(lcon2(1,:))
        % deuxieme fin de vaisseau
        z2 = lcon2(1,j) ;
        x2 = lcon2(2,j) ;
        y2 = lcon2(3,j) ;
        ind2 = sub2ind(size(skelsec),z2,x2,y2) ;
        %%on verifie si les fins sont dans un objet
        k = 1 ;
        verif = 0 ;
        while (k <=contouf.NumObjects)&&(verif == 0)
            vec = contouf.PixelIdxList{k} ;
            if (~isequal(size(find(ismember(vec,ind1))),[1,0]))&&(~isequal(size(find(ismember(vec,ind2))),[1,0]))
                path = shortestpath(graph,ind1,ind2) ;
                skel(path) = 1 ;
                v = setdiff(u,[z1,x1,y1],'rows') ;
                v = setdiff(u,[z2,x2,y2],'rows') ;
                lcon = v.' ;
                verif = 1 ;
            end
        end
    end
    i = i+1 ;
    toc
end
end
  end

It basically test all the loose ends in the skeleton and see if it can reconnect them. However, I have in my current skeleton 929 loose ends, and each test takes 30sec, which means I won't have any answer before hours ! It wouldn't normally bother me, but I will soon have to test this with way bigger skeletons, and I can't really afford to wait for 3 days for this function to finish.

Would anyone have some tips to increase the speed here ? I work on a very powerful machine with 32 GPU processor, but I don't think I can use them in parallel, as I use the result of the previous iteration in the new one.

Thanks for the help ! :)

Answer 1

[Guillaume]

The best way to speed up your code is to get rid of the loops. You have three nested loops. In your second loop, you keep converting almost the whole array into linear indices. You're basically repeating the same operation time and time again. You could simply calculate all the linear indices once before the loops, and just loop over them.

In essence, what you are doing is computing the cartesian product of the set of linear indices with itself, and see if the elements of that set belong to the pixels of each connected component. That can be done much faster:

%lcon: is a 3xN array

%convert all columns of lcon into linear indices at once:
indices = sub2ind(size(skelsec), lcon(1, :), lcon(2, :), lcon(3, :));  %convert all columns into linear indices at once
%compute cartesian product of indices with itself:
[ind1, ind2] = ndgrid(indices);
cartprodind = [ind1(:), ind2(:)];  %each row of cartprodind is a pair of indices
%remove pairs with identical indices:
cartprodind(cartprodind(:, 1) == cartprodind(:, 2), :) = [];
%I don't really understand what you're doing in your most inner loop
%the code you've written contain bugs (k doesn't change, you're overwriting v)
%so I assume it's not correct
%In any case, to check if a pair belongs to a connected component:
for rowpair = 1 : size(cartprodind, 1)
   for component = 1 : contouf.NumObjects
      if all(ismember(cartprodind(rowpair, :), contouf.PixelIdxList{component}))
         %do something
         break;
      end
   end
end

Answer 2

[Delphine Pommier]

Thanks for all of your answers ! Indeed, there was several errors that I missed, and I corrected them as soon as you pointed them ^^'

I can now do those calculations in a few minutes only, so thanks a lot ! I really appreciate you taking the time to answer.