How to optimize this code to run faster?

Hello,

I have a code that is attached that I am trying to optimize in order for it to run faster.

It contains a lot of for loops (which are very slow, I know). In the beginning this was the only way for me to code this as a beginner.

This is a small part of my original code, but it lacks in speed when I analyze it in MATLAB. It is very very slow, especially when I work with a large 3D set of images.

Can anyone help me with this, on how to optimize it. I would really appreciate it.

Thanks in advance!

Answer 1

[Adam Danz]

Your variable names induce anxiety in me :D

I'll tidy up the first for-loop so you can see how that's done.

I assume ndims(M) == 3 (ie, 'M' is 3D array)

If that's the case, 'Mmat' is a matrix and 'min1' is a column vector and

for ij1=1:numImages
    Mmin1(:,ij1)=Mmat(:,ij1)-min1(ij1,:);
end

can be replaced by

Mmin1 = Mmat - min1';

This works because the number of columns in Mmat equals the number of columns in min1' (transposed).

You could test that these two methods produced the same results by performing both of them and storing the results under different variable names. Then use isequal()

isequal(Mmin1, Mmin1_b)

Most of your for-loops follow similar logic so why don't you give it a shot and follow-up with specific areas where you got stuck, stating what you tried, etc.

[Mario]

Haha, I know, I do not know hat I was thinking while naming those variables. :D

Thanks for your help.

My Mmat is 222684x404 and min1 is 404x1. So that is why I put it in a for loop.

So your suggestion gives me an error "Error using - Matrix dimensions must agree. " as they are not the same size.

Any suggestions on how to solve this?

EDIT: Size of variable M is 462x482x404

[Mario]

Adam, I found a solution to my question above.

I am using currently:

 Mmin1=bsxfun(@minus, Mmat, min1');

It gives me the same result, but much faster. I will go through my code and if I get stuck at something, I will let you know.

[Walter Roberson]

bsxfun is a good solution to that.

Starting from R2016b, MATLAB would automatically detect the situation on a plain subtraction, and would automatically operate like bsxfun (though sometimes a little slower than bsxfun.)

[Mario]

Walter, do you have some faster alternative to regionprops, or reshape? How can I avoid for loop (if possible for a faster solution), on these examples:

s=struct('Centroid',[],'WeightedCentroid',[]);
for ij6=1:numImages
s = regionprops(logical(bw(:,:,ij6)), M2(:,:,ij6), 'Centroid'); %%bw and M2 are 3d sets of images
strCentroids(ij6).Centroids = [s.Centroid];
end

or for reshape:

for ij3=1:numImages
M2(:,:,ij3) = reshape(M1(:,ij3),Xsize, Ysize); %%M1 variable in this case is M x no. of images in 3D set - 22684x404 to be exact
end

Any ideas? I am using r2013b version.

[Walter Roberson]

reshape() is one of the very fastest MATLAB operations. reshape() does not copy data at all: it keeps the exact same internal pointer to data but creates a new header for it with the new size, and increments the "in use" count for the data block.

That said, you do not need to use a loop for that code: provided you are using all of M1, then the code can be just

M2 = reshape(m1, XSize, YSize, []);

[Mario]

Thanks Walter for the clarification. And you are right, I didn't needed for loop for reshape().

[Walter Roberson]

I do not think you are going to be able to do much better than looping with regionprops there, due to the need to identify isolated regions and find the centroid of each.

You are passing in the bw array to regionprops to identify regions, but because you are using logical values, regionprops detects regions out of the bw values: if you had passed in a double array then all of the 1 values together would have been treated as a single image, where with logical it means that the regions have to be located first.

The one thing you can do is preallocate:

 s(numImages) = struct('Centroid',[],'WeightedCentroid',[]);

This prevents strCentroids having to grow inside the loop like you are doing now.

[Mario]

Thanks again Walter for another hint.

I think I found my "black hole" where computational time is drastically slow.

It is this section of code, and I admit that it is not well written:

for ij8=1:numImages
ji = regionprops(logical(bw(:,:,ij8)), M2(:,:,ij8), {'Centroid','PixelValues','BoundingBox'});

jj(ij8).Centroids = [ji.Centroid];
jj(ij8).Centroids=reshape(jj(ij8).Centroids, 2, []);
jj(ij8).Centroids=transpose(jj(ij8).Centroids);

jj(ij8).PixelValues = {ji().PixelValues}.';
jj(ij8).PixelValues =cell2table(jj(ij8).PixelValues);
jj(ij8).PixelValues =table2struct(jj(ij8).PixelValues);

jj(ij8).BoundingBox = [ji.BoundingBox];
jj(ij8).BoundingBox=reshape(jj(ij8).BoundingBox, 4, []);
jj(ij8).BoundingBox=transpose(jj(ij8).BoundingBox);

for iji8=1:numImages
sP = regionprops(logical(bw2(:,:,iji8)), 'PixelList');
jj(iji8).PixelList={sP().PixelList}.';
jj(iji8).PixelList =cell2table(jj(iji8).PixelList);
jj(iji8).PixelList =table2struct(jj(iji8).PixelList);
end
end

I did not use tic toc to measure the time, but it is very slow.

[Walter Roberson]

jj(ij8).PixelValues = {ji().PixelValues}.';
jj(ij8).PixelValues =cell2table(jj(ij8).PixelValues);
jj(ij8).PixelValues =table2struct(jj(ij8).PixelValues);

could be replaced with

jj(ij8).PixelValues = struct('Var1', {ji.PixelValues}.');

and

jj(ij8).BoundingBox = [ji.BoundingBox];
jj(ij8).BoundingBox=reshape(jj(ij8).BoundingBox, 4, []);
jj(ij8).BoundingBox=transpose(jj(ij8).BoundingBox);

can be replaced with

jj(ij8).BoundingBox = vertcat(ji.BoundingBox);

It is not clear to me why you have a nested for loop on iji8 that continually overwrites jj(:).PixelList .

[Mario]

To be honest, I do not know either. :D

I was a beginner in MATLAB when I wrote this code which I now revisited in hopes to speed it up.

Can you help me sort the last part of this puzzle:

numPix=zeros(numImages,1);
for ii=1:numImages
numPix(ii) = numel(jj(ii).PixelValues);
end

numPixix=table(numPix);
numPixx=table2cell(numPixix);
for qq=1:ii
numPixxx{qq, 1}=[1:(numPixx{qq, 1})]';
end

%Extraction of standard deviation information
OOO1= {jj.PixelValues}.';
OOO2= {jj.PixelList}.';
OOOO= [OOO1,OOO2];

   for kk = 1:numel(numPixxx)
       for iiii=1:numel(numPixxx{kk,1}) 
       jj(kk).StandardDeviation(iiii).Var1 = std(double([jj(kk).PixelValues(iiii).Var1]));
       end
   end

sStd=cell2table(OOOO);
sStd=table2struct(sStd);
field = {'BoundingBox','PixelValues'};
SEW = rmfield(jj,field);

for uu=1:length(numPixx)
sStd(uu).OOOO = [SEW(uu).StandardDeviation]';
end

A1=(sStd);
A3=struct2table(A1);
N1 = size(numPix);

for ii2=1:N1
A22(ii2).Var1 = struct2table(A3.OOOO{ii2});
end
A3a=struct2table(A22);

for ii3=1:N1
A33(ii3).Var1 = table2array(A3a.Var1{ii3});
end 

for ii4 = 1:N1
MaxStD(ii4, :) = [max(A33(ii4).Var1)];
end

MaxStd=array2table(MaxStD);

for ii5 = 1:N1
[jj(ii5).MaxStd] = (MaxStd.MaxStD(ii5));
end

field1 = {'BoundingBox','PixelValues'};
SEW = rmfield(jj,field1);

for ii1=1:N1
SEW(ii1).Centroids = array2table(SEW(ii1).Centroids);
SEW(ii1).Centroids = table2struct(SEW(ii1).Centroids);
SEW(ii1).PixelList = array2table(SEW(ii1).PixelList);
SEW(ii1).PixelList = table2struct(SEW(ii1).PixelList);
end

for ik = 1:numel(numPixxx)
    for i7=1:numel(numPixxx{ik, 1})
if SEW(ik).MaxStd==SEW(ik).StandardDeviation(i7).Var1

       SEW(ik).Centr = SEW(ik).Centroids(i7);
       SEW(ik).Pixel = SEW(ik).PixelList(i7);
  end
      end
  end

for iki = 1:numel(numPixxx)
    for ii7=1:numel(numPixxx{iki, 1})
SEW(iki).Koord_Semena =[SEW(iki).Centr.Var1,SEW(iki).Centr.Var2];
SEW(iki).ListPix =[SEW(iki).Pixel.Var1.Var1];
    end
end

for ikI = 1:numel(numPixxx)
Koord_Semena(ikI,:) = SEW(ikI).Koord_Semena;

end

I am not sure what can be sorted/speed up in this section. This is the last part of this code that I am trying to modify/speed up.

[Adam Danz]

Wow, Walter helped you out a lot! I'll chip in two more examples.

numPix=zeros(numImages,1);
for ii=1:numImages
   numPix(ii) = numel(jj(ii).PixelValues);
end

can be replaced by

numPix = cellfun(@numel, {jj.PixelValues}');

As for

numPixix=table(numPix);
numPixx=table2cell(numPixix);
for qq=1:ii
    numPixxx{qq, 1}=[1:(numPixx{qq, 1})]';
end

assuming numPix is a vector (double, integers), it's replacement would be

numPixxx = cellfun(@(x)1:x, num2cell(numPix), 'UniformOutput', false);

Given the tools you learned from Walter and from these two examples, why don't you try the rest on your own. Google and this forum are great primary resources when getting stuck. If you've tried several methods and can't understand why they aren't working, follow-up here and include the methods you've tried so we can help set things straight.

[Mario]

Hi Adam, thank you for your valuable help and comments also.

I will try to complete the rest of the code myself and if I have any questions, I will post here.

Thanks!