Thread Subject: N-dimensional find

Hi there,

I've noticed myself running the following style of loop quite a bit when trying to find the 'first' pixel in an image along a dimension. This is somewhat like measuring from the top row of an image towards the bottom row of the image, and recording the first collision with an "on" pixel.

% BW is some binary image

firstIdxs = zeros(size(BW,2);
for i=1:size(BW,2)
  idx = find(BW(:,i),1,'first');
  if ~isempty(idx)
    firstIdxs(i) = idx;
  end
end

I see this as analogous to using, say, sum(BW, 2) to sum along the 2nd dimension, whereas I'm trying to do something like find(BW, dim).

Is there a faster solution? Or, perhaps more relevant: is there a solution that doesn't clutter up my .m files as much as this implementation? If not, I'll just have to function-ize this one.

Cheers,
Sven.

"Sven" <sven.holcombe@gmail.deleteme.com> writes:

> Hi there,
>
> I've noticed myself running the following style of loop quite a bit
> when trying to find the 'first' pixel in an image along a
> dimension. This is somewhat like measuring from the top row of an
> image towards the bottom row of the image, and recording the first
> collision with an "on" pixel.
>
> % BW is some binary image
>
> firstIdxs = zeros(size(BW,2);
> for i=1:size(BW,2)
> idx = find(BW(:,i),1,'first');
> if ~isempty(idx)
> firstIdxs(i) = idx;
> end
> end
>
> I see this as analogous to using, say, sum(BW, 2) to sum along the 2nd
> dimension, whereas I'm trying to do something like find(BW, dim).
>
> Is there a faster solution? Or, perhaps more relevant: is there a
> solution that doesn't clutter up my .m files as much as this
> implementation? If not, I'll just have to function-ize this one.


Not sure about faster, but you could abuse the behavior of max:

[dummy firstIdxs] = max(logical(BW));

That is, 1 is the max value, and the location of that max value is the
one you're looking for. "help max" says that the first maximum is the
one that counts.

If this is faster, it's only because the loop overhead is bad. The find
solution technically only has to look at the values up to the first 1 in
each column, where the max solution has to look at every single value in
the matrix.

-Peter

"Sven" <sven.holcombe@gmail.deleteme.com> wrote in message <gnhibi$s1m$1@fred.mathworks.com>...
> Hi there,
>
> I've noticed myself running the following style of loop quite a bit when trying to find the 'first' pixel in an image along a dimension. This is somewhat like measuring from the top row of an image towards the bottom row of the image, and recording the first collision with an "on" pixel.
>
> % BW is some binary image
>
> firstIdxs = zeros(size(BW,2);
> for i=1:size(BW,2)
> idx = find(BW(:,i),1,'first');
> if ~isempty(idx)
> firstIdxs(i) = idx;
> end
> end
>
> I see this as analogous to using, say, sum(BW, 2) to sum along the 2nd dimension, whereas I'm trying to do something like find(BW, dim).
>
> Is there a faster solution? Or, perhaps more relevant: is there a solution that doesn't clutter up my .m files as much as this implementation? If not, I'll just have to function-ize this one.
>
> Cheers,
> Sven.

Something like this?

bw = rand(10) > .8 ;

firstidx = zeros(1,size(bw,2)) ;
[R,C] = find(cumsum(cumsum(bw))==1) ;
firstidx(C) = R % we need some care for all-zero columns

although I think the double cumsum makes it slower than a for-loop ;-)

hth
Jos

If it is speed you are after (and it sounds like maybe it isn't), it might be worth it to make a sub function with not one, but two loops (blasphemy!). Try this out, the speed gain will depend on the size of BW and it's sparsity. Run it a couple of times.



%--------------------------------------------------------------------------%
function [] = find_idx_dim()
% Demonstrate three different techniques for finding the index of the first
% non-zero element of each column in a binary matrix.
BW = round(rand(1000)); % A binary matrix.
% BW = rand(1000)>.9; % A sparse binary matrix. Try this out too.
t = zeros(1,3); % For timings.

tic
[FI1,FI1] = max(BW); % Technique 1.
t(1) = toc;

tic
FI2 = zeros(1,size(BW,2)); % Technique 2.
[R,C] = find(cumsum(cumsum(BW))==1);
FI2(C) = R;
t(2) = toc;

tic
FI3 = getidx(BW); % Technique 3, Call subfunction.
t(3) = toc;

% Show the relative times and if all are equal.
fprintf('\n%s%3.1f\t%3.1f\t%3.1f','Relative times: ',t/min(t))
if all(FI1==FI2 & FI1==FI3') % check for equality.
    fprintf('%s\n\n',' and all are equal')
else
    fprintf('%s\n\n',' and all are not equal')
end

function idx = getidx(BW)
% Subfunction that finds the first 1 index in each column.
[r,c] = size(BW);
idx = zeros(r,1);
for ii = 1:c
    for kk = 1:r
        if BW(kk,ii)==1
            idx(ii) = kk;
            break
        end
    end
end

%--------------------------------------------------------------------------%

Thanks all,

The short-term solution I've chosen which has cleaned up my code quite a bit is a combination of suggestions:
BW = rand(10,20)>.8 % Make the image
[A, B] = max(BW) % Hijack the max function
B(BW(B==1)==0) = 0 % Account for all-zero columns


I think that Matt's double for-loop probably is the fastest, and I am thinking of adding something like it to the file exchange (or Matt, you're more than welcome to add it if you'd like!). I would probably modify it a little to allow you to input the direction ala the SUM(BW,DIR), or MAX(X,[],DIM) functions.

I can think of how to make it generic enough for a 2D image called by
find_2d(BW,1)
OR
find_2d(BW,2)

Perhaps when I get a moment I'll try and work out how to make it work in ND. Unless of course, someone already knows ;)

Cheers,
Sven.


"Matt Fig" <spamanon@yahoo.com> wrote in message <gnhrhd$llc$1@fred.mathworks.com>...
> If it is speed you are after (and it sounds like maybe it isn't), it might be worth it to make a sub function with not one, but two loops (blasphemy!). Try this out, the speed gain will depend on the size of BW and it's sparsity. Run it a couple of times.
>
>
>
> %--------------------------------------------------------------------------%
> function [] = find_idx_dim()
> % Demonstrate three different techniques for finding the index of the first
> % non-zero element of each column in a binary matrix.
> BW = round(rand(1000)); % A binary matrix.
> % BW = rand(1000)>.9; % A sparse binary matrix. Try this out too.
> t = zeros(1,3); % For timings.
>
> tic
> [FI1,FI1] = max(BW); % Technique 1.
> t(1) = toc;
>
> tic
> FI2 = zeros(1,size(BW,2)); % Technique 2.
> [R,C] = find(cumsum(cumsum(BW))==1);
> FI2(C) = R;
> t(2) = toc;
>
> tic
> FI3 = getidx(BW); % Technique 3, Call subfunction.
> t(3) = toc;
>
> % Show the relative times and if all are equal.
> fprintf('\n%s%3.1f\t%3.1f\t%3.1f','Relative times: ',t/min(t))
> if all(FI1==FI2 & FI1==FI3') % check for equality.
> fprintf('%s\n\n',' and all are equal')
> else
> fprintf('%s\n\n',' and all are not equal')
> end
>
> function idx = getidx(BW)
> % Subfunction that finds the first 1 index in each column.
> [r,c] = size(BW);
> idx = zeros(r,1);
> for ii = 1:c
> for kk = 1:r
> if BW(kk,ii)==1
> idx(ii) = kk;
> break
> end
> end
> end
>
> %--------------------------------------------------------------------------%

"Sven" <sven.holcombe@gmail.deleteme.com> wrote in message <gni2rp$3ja$1@fred.mathworks.com>...
> Thanks all,
>
> The short-term solution I've chosen which has cleaned up my code quite a bit is a combination of suggestions:
> BW = rand(10,20)>.8 % Make the image
> [A, B] = max(BW) % Hijack the max function
> B(BW(B==1)==0) = 0 % Account for all-zero columns
>
...

If anyone is trying to use this neat trick with *sparse* logical matrices in R2009a, be aware that a MATLAB bug causes a crash:

http://www.mathworks.com/matlabcentral/newsreader/view_thread/246673

So, if your matrix is sparse and you want your code to be R2009a-compliant, you'll have to do something like:

[A, B] = max(spones(BW)) % converts BW to double to prevent bug

Took me a while to work out what was causing new crashes when moving to R2009a!

AJ

"Andrew Jackson" <andrewj@eml.cc> wrote in message
news:gu171n$rg$1@fred.mathworks.com...
> "Sven" <sven.holcombe@gmail.deleteme.com> wrote in message
> <gni2rp$3ja$1@fred.mathworks.com>...
>> Thanks all,
>>
>> The short-term solution I've chosen which has cleaned up my code quite a
>> bit is a combination of suggestions:
>> BW = rand(10,20)>.8 % Make the image
>> [A, B] = max(BW) % Hijack the max function
>> B(BW(B==1)==0) = 0 % Account for all-zero columns
>>
> ...
>
> If anyone is trying to use this neat trick with *sparse* logical matrices
> in R2009a, be aware that a MATLAB bug causes a crash:
>
> http://www.mathworks.com/matlabcentral/newsreader/view_thread/246673

FYI, that's bug report 529803 if you want to track it.

http://www.mathworks.com/support/bugreports/

*snip*

--
Steve Lord
slord@mathworks.com