Can this function work any faster?

Hello All, Well I'm making a function that takes about 0.022 seconds to complete, I need to try to make it even less time then that, just getting peoples opinions. The end result needs to run at 120 Hz, I have tried most matlab built in functions, but there too slow for what I need, I know because I've tried. This is designed to be used for an accelerator beam image, it is just a "spot" on the screen and the image coming in has about 606x1182 pixels in it. The inputs are an image and a threshold(which is usually == 0). I pretty much have striped down to bare essentials for this code, however hope there was something I missed that will make it faster. I was just hoping for peoples comments or options about any way to make this code faster(or if there is a better way):

function [clipedImg,tfr,tfc] = ImageClipper(useImg,threshold)          
useImg = double(useImg); %????????????????????????

mx = max(useImg);
my = max(useImg,[],2);

sorted=sort(my(:));
dist=diff([sorted; sorted(end)-1]);
idx=find(dist);
num=[idx(1); diff(idx)];
n=max(num);
ty=sorted(idx(num==n));

v = mx > ty;
na = length(v);
b = [1 -1];
[k,zf] = filter(b, 1, v);
k(na+1:na+2-1) = zf;
P = find(k==1);
L = find(k==-1)-P;
LP = [max(L) P(L==max(L))];
tfc = [false(1,LP(1,2)-threshold) true(1,LP(1,1)+2*threshold) false(1,length(mx)-sum(LP)-threshold)];

v = my > ty;
na = length(v);
b = [1 -1];
[k,zf] = filter(b, 1, v);
k(na+1:na+2-1) = zf;
P = find(k==1);
L = find(k==-1)-P;
LP = [max(L) P(L==max(L))];
tfr = [false(1,LP(1,2)-threshold) true(1,LP(1,1)+2*threshold) false(1,length(my)-sum(LP)-threshold)];

clipedImg = useImg(tfr,tfc);

Thanks!

[Matt Kindig]

Can you post a sample image as well for us to test with?

[dpb]

Looks like lots of sort and find operations--they're certainly time-consuming.

BUT, have you profiled it to confirm where the bottleneck really are? There's no point in just random attempts at optimization by code changes until you know where the time is being spent.

After that, in general the real gains in speed are virtually always "algorithm, algorithm, algorithm". Suggest you post what the actual processing objective is in definitive descriptive prose -- that has better chance of promoting a more useful direction of discussion on how to implement that more efficiently than having somebody just look at undocumented code.

[Chris E.]

For Matt's qustion,

I wish I could! I'm really not sure how and it may loose its information if I do. But the next best thing I can do is give you a "simulated" image:

function ret = customgauss(gsize, sigmax, sigmay, theta, offset, factor, center)
ret     = zeros(gsize);
rbegin  = -round(gsize(1) / 2);
cbegin  = -round(gsize(2) / 2);
for r=1:gsize(1)
    for c=1:gsize(2)
        ret(r,c) = rotgauss(rbegin+r,cbegin+c, theta, sigmax, sigmay, offset, factor, center);
    end
end
end

function val = rotgauss(x, y, theta, sigmax, sigmay, offset, factor, center)
xc      = center(1);
yc      = center(2);
theta   = (theta/180)*pi;
xm      = (x-xc)*cos(theta) - (y-yc)*sin(theta);
ym      = (x-xc)*sin(theta) + (y-yc)*cos(theta);
u       = (xm/sigmax)^2 + (ym/sigmay)^2;
val     = offset + factor*exp(-u/2);
end

To use this run the following:

useImg = customgauss([1182 1182], 20, 20, 3, 2, 100, [30 50]);

That will give you the image you need (or is close enough) and just use 0 for the threshold.

[Chris E.]

As for dpb's comment, Well I used the profiler and it appears to be all focused on the max function. More particularly:

mx = max(useImg);
my = max(useImg,[],2);

I did dry to use sum instead of max and came up with something quicker and really ugly, it seems to go faster but is not so great(I think) and I'm sure it can be cleaned up a little:

function [clipedImg,tfr,tfc] = ImageClipper(useImg,border,threshold)   
[m n] = size(useImg);
mx = sum(useImg)*1/m;
my = sum(useImg,2)*1/n;

sorted=sort(my);
idx=find(diff([sorted; sorted(end)-1]));
num=[idx(1); diff(idx)];
ty=sorted(idx(num==max(num)));
ty = ty(end);

[k,zf] = filter([1 -1], 1, (mx>ty));
k(n+1:n+2-1) = zf;
P = find(k==1);
[L posP]= max(find(k==-1)-P);
P = P(posP);
tfc = [zeros(1,P-threshold) ones(1,L+2*threshold) zeros(1,n-(L+P)-threshold)];

[k,zf] = filter([1 -1], 1, (my>ty));
k(m+1:m+2-1) = zf;
P = find(k==1);
[L posP]= max(find(k==-1)-P);
P = P(posP);
tfr = [zeros(1,P-threshold) ones(1,L+2*threshold) zeros(1,m-(L+P)-threshold)];

useImg = useImg+(tfr'*tfc);

mx = sum(useImg)*1/m;
my = sum(useImg,2)*1/n;

[k,zf] = filter([1 -1], 1, (mx>ty));
k(n+1:n+2-1) = zf;
P = find(k==1);
[L posP]= max(find(k==-1)-P);
P = P(posP);
tfc = [false(1,P-border) true(1,L+2*border) false(1,length(mx)-(L+P)-border)];

[k,zf] = filter([1 -1], 1, (my>ty));
k(m+1:m+2-1) = zf;
P = find(k==1);
[L posP]= max(find(k==-1)-P);
P = P(posP);
tfr = [false(1,P-border) true(1,L+2*border) false(1,length(mx)-(L+P)-border)];

clipedImg = useImg(tfr,tfc);

[Chris E.]

I also made the orignal post with max a tad faster, but still not 120 Hz fast, here is the updated code:

function [clipedImg,tfr,tfc] = ImageClipper(useImg,threshold)        
[m n] = size(useImg);

mx = max(useImg);
my = max(useImg,[],2);

sorted=sort(my);
idx=find(diff([sorted; sorted(end)-1]));
num=[idx(1); diff(idx)];
ty=sorted(idx(num==max(num)));
% ty = ty(end);

[k,zf] = filter([1 -1], 1, (mx>ty));
k(n+1:n+2-1) = zf;
P = find(k==1);
[L posP]= max(find(k==-1)-P);
P = P(posP);
tfc = [false(1,P-threshold) true(1,L+2*threshold) false(1,n-(L+P)-threshold)];

[k,zf] = filter([1 -1], 1, (my>ty));
k(m+1:m+2-1) = zf;
P = find(k==1);
[L posP]= max(find(k==-1)-P);
P = P(posP);
tfr = [false(1,P-threshold) true(1,L+2*threshold) false(1,m-(L+P)-threshold)];

clipedImg = useImg(tfr,tfc);

[dpb]

I still think it would be more productive to describe the desired algorithm instead...

[Chris E.]

Well all I need is to choose the "beam spot" in the image really quickly, then crop it with a border of x pixels, that is basically it. I need to do this for a stream of images coming in at 120 Hz, that is the hard part. The needed information is the beam centroid and the image of the beam cropped its self. I hope that helps!

[Walter Roberson]

http://www.mathworks.com/matlabcentral/answers/7924-where-can-i-upload-images-and-files-for-use-on-matlab-answers

http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup

[dpb]

I'm no image processor...would seem if you can't figure out alternate algorithm to minimize the amount of processing further the way to get speed would be

a) faster processor,

b) mex what you have

(not necessarily in that order)

What's the specific operational definition of " the image of the beam cropped"?

But, isn't the centroid always going to be

[mx,imx]=max(image(:));
[i,j]=sub2ind(size(useImg),imx);

?

That's still a time-consuming MAX, but only one, not multiple and no additional SORTs, FINDs, etc., ...

Then, the more specifics you can provide re: the bounds on what is the range in the image over which you have to operate the better.

Basically, given that image, (again as no imager; never done any "in anger" in 40+ yr career) I'd use the above starting point then look for first crossing below some threshold to the l/r and u/d to find the range boundaries and use that instead of the whole thing. How much more you need will depend on how noisy the input is, basically.

Does the beam size/shape change markedly or only the centroid position? And can it be anywhere in the frame from one to the next or is it going to be somewhere near where was last time once find it the first image? IOW can you not bound the search significantly and update location incrementally rather than start from scratch every cycle?

[Chris E.]

Hello dpb,

Well as for the centroid, it changes depending on where "more" of the beam will be, it is sometimes lopsided or overpowered in one particular area. As for the "jitter" of the beam, it will move, but usually in a smaller area. It usually has some kind of movement from frame to frame, but usually is small. I think there is a way to cut the image to an area of interest if needed and I'm looking it to that. The bounding box just needs to enclose the 95% of the beam, a Gaussian distribution in 2d usually is a good mimic of how the beam looks and acts. Thank you for the thoughts I will defiantly think of implementing them to see how it goes. I did try mex and found it made things faster, as for the faster processor, I'm kind of limited to my machinery I need to use, however it will be the only thing that runs on it. Please if you have any more questions I would love to answer them. Any speed increase is going to be helpful.

[Walter Roberson]

http://www.mathworks.com/matlabcentral/answers/7924-where-can-i-upload-images-and-files-for-use-on-matlab-answers

Answer 1

[Sean de Wolski]

Another husky!

You could avoid the max operation here as well:

[L posP]= max(find(k==-1)-P);

Instead:

find(k==-1,1,'last')-P

Now it only needs to find the first occurence of a -1 and not all of them allowing it to stop prematurely.

Avoid operations like the following that change the size or class of an array:

P = P(posP);

Also, to address DPB's point: do you have MATLAB Coder? If so, you could try generating code from this and then compiling it into a MEX file that might run faster.

Answer 2

[dpb]

My simple-minded solution -- undoubtedly will need some robustness added, but try adding sophistication to something as trivial as the following for a first cut...

[mx,imx]=max(image(:));               % find the max point as center
[i,j]=sub2ind(size(useImg),imx);        % turn into row, column
lvl=0.05*mx;                            % a cutoff threshold
% left/right and upper/lower bounds enclosing that level around max 
L=find(useImg(i,1:j)>lvl,1,'first');R=find(useImg(i,j:end)<lvl,1,'first')+j;
U=find(useImg(1:i,j)>lvl,1,'first');B=find(useImg(i:end,j)<lvl,1,'first')+i;
imClip=useImg(U:B,L:R);                 % the bounded image
image(imClip)

As my previous suggestion, I would think that even if take more sophistication the first time, after that if the beam can't be just anywhere you could limit the area over which you have to do the max to a fairly small region and thereby make it quite a lot smaller area for the search and thus be able to add more robustness over a smaller area to gain compute speed.

Then again, if the signal is pretty robust, maybe something as trivial as just doing the bounding on a strip P pixels wide centered on the max location instead of a single row/column as here would do it. Again, I'm no ImageAnalyst by any stretch...

[Chris E.]

Thank you for helping, this could be a possible solution and I do think it may help. I have used the simple version, i think your right in terms of robustness, however if I take a few of the max points and then try it out (so the max and then 2nd place and 3 rd place) it should be a decent and fast way of doing this. I will look into it for sure! Thanks

[dpb]

If the size of the beam image is relatively constant I'd think you could take the overall max, then look at an area around it of fixed size and pick a group of points within some relatively small tolerance of that value. From those you ought to be able to get a band in each direction from which you can do the bounding instead of the single line as I illustrated. Overall, however, this could be only a very small fraction of the overall image size which has to help speed.

Then, one thing you might note -- if you do use the diff() approach, it works well if you first subtract a threshold between the noise levels of the min and max and use sign() on that result. Then diff(sign(x-th)) leaves you w/ +/-2 at the zero crossings.