Apply 2D threshold to 3D array

Question

With a geographical grid the size 20x30, I have two (temperature) variables:

The data A with the size 20x30x100
and a threshold the size 20x30

I'd like to apply the threshold to the data, i.e. to cut out the values in A that are above threshold, with every grid point having its own threshold. Since that will give a different number of values for each grid point, I thought to pad the rest with zeros, so that the resulting variable, let's call it B, will also be of the size 20x30x100.

I was thinking to do something like this, but there's something wrong with the loop:

B = sort(A,3); %// sort third dimension in ascending order
threshold_3d = repmat(threshold,1,1,100); %// make threshold into same size as B

for i=1:20
    for j=1:30
        if B(i,j,:) > threshold_3d(i,j,:); %// if B is above threshold
          B(i,j,:); %// keep values
        else
          B(i,j,:) = 0; %// otherwise set to zero
        end
    end
end

What is the correct way to do the loop?
What other options to do this are there?

Thanks for any help!

Answer 1

You can use bsxfun for a much more efficient solution which will internally take care of the replication being done with repmat, like so -

B = bsxfun(@times,B,bsxfun(@gt,B,threshold))

A more efficient solution might be to use logical indexing to set the False elements from the mask created by bsxfun(gt, i.e. True ones with the use of bsxfun(@le in B to zeros thereby avoiding the use of bsxfun(@times, which for huge multidimensional arrays could be a bit expensive, like so -

B(bsxfun(@le,B,threshold)) = 0

Note on efficiency : Being a relational operation, going with the vectorized operation with bsxfun would provide both memory and runtime efficiency. The memory efficiency part has been discussed here - BSXFUN on memory efficiency with relational operations and the performance numbers have been researched here - Comparing BSXFUN and REPMAT.

Sample run -

 >> B
 B(:,:,1) =
      8     3     9
      2     8     3
 B(:,:,2) =
      4     1     8
      4     5     6
 B(:,:,3) =
      4     8     5
      5     6     5
 >> threshold
 threshold =
      1     3     9
      1     9     1
 >> B(bsxfun(@le,B,threshold)) = 0
 B(:,:,1) =
      8     0     0
      2     0     3
 B(:,:,2) =
      4     0     0
      4     0     6
 B(:,:,3) =
      4     8     0
      5     0     5