Extract a specified interval of numbers from a multidimensional array

Hello,

I have a multidimensional array, let's call it error, of the size (37,4,4,3) which is basically filled with numbers ranging from 0-1 describing the amount of error for a certain application. What I want to do is extract intervals for each vector (:,i,j,k) of the total array, when the numbers are within a certain tolerance of 0.03. For example if we take error(:,1,1,1), I want to be able to get a strict interval with the values from 0-0.03.

So far I have been using the find function in a for-loop to determine the 'first' position which is a value below 0.03 and the 'last' position which is a value below 0.03. The problem with this is that sometimes there is a value exceeding 0.03 within those positions, which I don't want.

for k = 1:3
    for j = 1:4
        for i = 1:4
            pos_first(:,i,j,k) = find(error(:,i,j,k) <= 0.3,1,'first');
            pos_last(:,i,j,k) = find(error(:,i,j,k) <= 0.3 & error(:,i,j,k) > 0,1,'last');
        end
    end
end

I also tried to write error(error < 0.03 & error > 0), but that bunched all of the numbers into one single vector, which I don't want. I need to have it in a format where I can distinguish between the dimensions (:,4,4,3).

I therefore wonder if anyone know a smart way to accomplish this.

/ Jonatan

Answer 1

[Guillaume]

Note: calling your variable error may not be wise, as it'll prevent you from using the error function.

error <= 0.3 is a logical array of 0 and 1. What you're basically asking is how to detect the first continuous sequence of 1s. This is easily achieved by taking the diff of that sequence. Any transition from 0 to 1 will result in 1 in the diff and transition from 1 to 0 will result in -1, while unchanged values result in 0. So you're simply left with detecting the first 1 and -1 pair:

for k = 1:size(error, 4)  %don't hardcode endpoints!
    for j = 1:size(error, 3)
        for i = 1:size(error, 2)
            diffseq = diff([0; error(:, i,j,k) <= 0.3; 0]);  %borders with 0 to make sure we have transition from 0 to 1 and 1 to 0 even if the sequence starts or ends with 1.

              pos_first(i,j,k) = find(diffseq == 1, 1);
              pos_last(i,j,k) = find(diffseq == -1, 1) - 1;
          end
      end
  end

Note: you could achieve the same without a loop:

sz = size(error);
sz(1) = 1;
diffseq = diff([zeros(sz); error <= 0.3; zeros(sz)]);  %padding accross all dimensions but rows
[i, j, k, l] = ind2sub(size(error) + [1, 0, 0, 0], find(diffseq == 1));
pos_first = accumarray(i, [j, k, l], [], @min);
[i, j, k, l] = ind2sub(size(error) + [1, 0, 0, 0], find(diffseq == -1));
pos_last = accummaray(i, [j, k, l], [], @min) - 1;

No idea if it's faster. It's certainly not clearer.

[Jonatan Boberg]

Hello Guillaume,

Thank you very much for your reply. I tried the diff function as you mentioned and it worked well. However, in some cases the diff function gives me the following output:

ans =

       1
       0
      -1
       0
       0
       0
       1
       0
       0
       0
       0
       0
       0
       0
       0
       0
       0
       0
       0
       0
       0
       0
       0
       0
       0
       0
       0
       0
       0
       0
       0
       0
       0
      -1
       0
       0
       0
       0

And for that case I would like to get the interval between the second instance of 1 to -1. In short, I always want to get the interval with the most values in it. Do you know of any way to always get the largest interval of values which lies between 1 and -1? Much appreciated!

EDIT: I mean that I want to get the positions of the longest interval.

/ Jonatan

[Guillaume]

I can't think of a (fast) way without a loop. It's trivial with the loop, the difference between find(diffseq == 1) and find(diffseq == -1) is the length of the sequences. Just pick the longest one and its corresponding indices:

for k = 1:size(error, 4)  %don't hardcode endpoints!
   for j = 1:size(error, 3)
      for i = 1:size(error, 2)
          diffseq = diff([0; error(:, i,j,k) <= 0.3; 0]);  %borders with 0 to make sure we have transition from 0 to 1 and 1 to 0 even if the sequence starts or ends with 1.
          runstarts = find(diffseq == 1);
          runends = find(diffseq == 1) - 1;
          [~, longest] = max(runends - runstarts);
          pos_first(i,j,k) = runstarts(longest);
          pos_last(i,j,k) = runends(longest);
      end
   end
end