one line takes ages, but is it possible to optimise any more?

Need a little help, I'm running this code a lot and small differences lead to big time savings.

I have a large matrix (512x512), if you imagine this matrix is a square and then draw a circle fully inside this square, the pixels in this circle need to be identified quickly. So far I've done this with a pre-calculated logical 512x512 matrix that has a 1 if the pixel is in the circle, and a 0 otherwise, then I do:

for i = 1:512^2 if circ(i) == 1 ... end end

Still, this "circ(i) == 1" test is using over 50% of my time(there are for loops outside this, so we're searching this matrix many millions of times). I tried other things, like:

cir = 1:512^2; cir = cir(circ); for i = cir end

^^ In this case I'm pre-calculating even the exact indices I need, but this takes even longer. I've tried using geometry and defining a radial distance to the centre of the matrix but this takes a little longer as well.

I guess it must be getting slowed down by the RAM reading speed, but I don't know how to do anything about this.

Anyway thanks for any suggestions.

edit- curious if people think it's a memory limit... because if I convert it from logical to double, there is a small time advantage. I have no idea how that works.

Answer 1

[Andrew Newell]

If circ is your pre-defined logical matrix and c is the matrix to test, you could use

incircle = circ & cir;

to get a logical matrix with ones for the entries you want; and

I = find(circ & cir);

if you need the indices.

[Michael]

Thanks for your reply. I'm not sure what this does, I can't test it because the two matrices are different sizes (I think cir is a little bit truncated because the first and last few elements of circ are zero);

But I already have the indices in cir, it just takes really long to say "for i = cir"; for a reason I don't understand, "for i = 1:512^2; if circ(i) == 1" is much faster.

Mike

Answer 2

[Hin Kwan Wong]

Your code is running slow because MATLAB is not good for running serial for loops. It's good for vectorised calculations.

While there are probably more efficient way in terms of image processing to your problem, ny proposed method is simple to understand:

If the radius of the circle is always touching the matrix boundary, you have the equation of the circle and it is just a matter of quantisation, you can divise formula for every rows of the matrix:

You use the circle equation X^2 + Y^2 = r, so X=plus or minus sqrt(r-y^2) For each Y value (row numbers of the matrix offset from centre) you and find X Which is offset from the centre of the matrix. From this you fill up that row of the matrix with one.

matrix=zeros(512);
N=length(matrix) %=512
for row = 1:N
Xoffset = (sqrt(abs(r-(Y-N/2)^2))
Xmin = round(N/2-Xoffset);
Xmax = round(N/2+Xoffset);
matrix(Xmin:Xmax,row)=ones(Xmax-Xmin,1);
end

Of course the above is simplified example, you need to apply checks for Xmin never goes to 0, Xmax doesn't get out of range, the appropiate quantisation using round(), ceil() and floor(), etc.

Answer 3

[Hin Kwan Wong]

Sorry I may have misinterpreted your question

But then it seems easy to me

You use use Data(circ) to extract a linear vector containing all data points on a matrx fitting inside that circle

for example Intersect = Data(circ) Then you can do traditional tricks with logical indexing:

Intersect(some criteria) = something else eg. Intersect(Intersect>10)=0; or find(Intersect>10)

[Michael]

Thanks for the reply, I've thought about vectorizing it, but it's become problematic, I'm either not capable enough as a programmer to do it or it's not possible because there are further for/if loops inside this one, and some of these are criteria which depend on i, so I struggle to figure how how to put this as a logical argument without having any variable "i" to work with. I'm giving it another go, but it's a good 50 lines of text!

[Hin Kwan Wong]

The for loop with i will get into the way, if you vectorize you need to vectorize everything including the operations in that for loop.

Data(circ) contains the intersected data. But find(circ) gives all the index (i) values of them.

Search for Techniques for Improving Performance in MATLAB help

There are a lot of good guidance in there into writing faster code, and vectorization techniques.

You can try to play around with functions like arrayfun