how to modify the variables for changing FOR loop to PARFOR loop??

I have a code described here which takes lot of time to run on my 32 cores machine. I tried to use all cores by changing for loop by parfor loop. I get an error about in final_lmatrix which ccannot be classied.

CC = bwconncomp(B_fill_slice);
L = labelmatrix(CC);
numObjects = CC.NumObjects;
prev_regions = 0; %initialize
new_regions = numObjects;
disp('init_regions')
disp(new_regions);
iter = 1;
while ((new_regions - prev_regions)>0)
    disp('sub');
    disp(new_regions - prev_regions);
    prev_regions = new_regions;
    disp('iter');
    disp(iter);
    new_num = 0;
    final_lmatrix = zeros(r,c,p, 'uint8'); % preallocation
    for n = 1:new_regions
        region = (L==n);
        l = water(region); %compute watershed and return labelmatrix
        if (n~=1)
            l = reassign_label(uint8(l),r,c,p,new_num);
        end
        final_lmatrix = l+final_lmatrix; %contains the updated label matrix to be used for the next iteration
        new_num = max(final_lmatrix(:));
    end  %end for loop
    new_regions = new_num;
    L = final_lmatrix;
    disp('new_regions');disp(new_regions);
    iter = iter+1;
end %end while loop

thanks

[Rik]

The entire matrix is required in each iteration. This setup is not possible to use in a parfor. You will have to edit it so that each iteration is independent of the others (e.g. by having a container array you index).

[Rizwan Khan]

could you please help in editing this as I am new in parallel computing??

[Rik]

Every iteration should independent. That means you can only use n as an index. If you don't have memory issues, that means you should use final_lmatrix(:,:,:,n) inside the loop.

Answer 1

[Jan]

What is the bottleneck of your code? Use the profiler to find it. If it is the water() command, parallelize this command:

lArray = zeros(r,c,p, new_regions, 'uint8')
parfor n = 1:new_regions
    lArray(:, :, :, n) = water(L == n);
end

Does this work? Then the re-assigning of the lables runs in a separate sequential loop, or you can create a more efficient method.

[Rizwan Khan]

Yes, water ( ) command is the bottleneck but issue I am getting with the function of re-assigning lables. We have to update the label numbers in the subsequest regions.

Purpose of reassign label function: I am applying my own watershed transform locally on the 3D image. Suppose region 1st has 5 labels numbering 1 to 5 and region 2nd has 3 labels which gives lables as 1 to 3 and therefore I ressign the lables to start as 6,7, and 8 to avoid overwriting and son on.

Answer 2

[Bruno Luong]

Could you try this inner loop (EDIT fix issue with 0 label):

final_lmatrix = 0;
parfor n = 1:new_regions
    region = (L==n);
    lbn = water(region); %compute watershed and return labelmatrix
    m = max(lbn);
    LUT = [0, (1-2*2^-n)+2^-n*((1:m)/m)];
    lbn = LUT(lbn); % I don't know what is your reassign_label function does
    final_lmatrix = final_lmatrix + lbn;
end  %end for loop
[~,~,final_lmatrix] = unique(final_lmatrix);
final_lmatrix = final_lmatrix - 1;

[Rizwan Khan]

reassign label function: I am applying my own watershed transform locally on the 3D image. Suppose region 1st has 5 labels numbering 1 to 5 and region 2nd has 3 labels which gives lables as 1 to 3 and therefore I ressign the lables to start as 6,7, and 8 to avoid overwriting.

[Bruno Luong]

In principle my method does the reassignment of labels without collision and doesn't need to be run sequential, so it can run in parallel.

Though the indexing but LUT might be not the most efficient as it is impemented. If this code works you can use the same idea but compute directly the reassignment without using LUT.

[Bruno Luong]

Perhaps a better code:

final_lmatrix = 0;
parfor n = 1:new_regions
    region = (L==n);
    lbn = water(region); %compute watershed and return labelmatrix
    m = max(lbn,[],'all');
    b = lbn > 0;
    lbn(b) = n + lbn(b)/m;
    final_lmatrix = final_lmatrix + lbn;
end  %end for loop
[~,~,final_lmatrix] = unique(final_lmatrix);
final_lmatrix = final_lmatrix - 1;

[Rizwan Khan]

sorry, its not working. Problem is in reassigning the labels.

% Please look into the reassining function I used but now using vectorized form for speed:

function o_lmat = reassign_label(o_lmat,h,w,p,n_prev)

for a = 1:h % for loops to reassign labels

for b = 1:w

for c = 1:p

if (o_lmat(a,b,c) ~= 0)

o_lmat(a,b,c) = o_lmat(a,b,c)+n_prev;

end

end

end

end

% vectorized form: I am using this now

function o_lmat = reassign_label(o_lmat, h, w,p,n_prev)

idx = o_lmat(:) ~=0;

o_lmat(idx) = o_lmat(idx) + n_prev;

end

[Bruno Luong]

So it looks like your label is higher for n small (pushed up by later iterations).

I did the opposite since I didn't know the detail of your function.

It's really not hard to change to match your labeling convention.