Logical Indexing instead of loops in 3D-Matrices

Dear All Community members,

I have two 3D matrices, AB1/AB2=1250x441x451 and two arrays, C=1x5 and D=1x1250.

Instead of running loops I want to use a vectorization code to do the following operation;

If D(i)<=C (for each value in C)

ABC(i)=AB1(i)

else

ABC(i)=AB2(i).

As an end-result I want five different 3D matrices ("ABC1", "ABC2"..."ABC5"), one for each value in C.

In order to use vectorization I could repeat the 3D matrices ("AB1" and "AB2") and the 1D array ("D") 5 times using repmat and repeat the 1D array ("C") 1250 times using repelem. But I am not sure how to proceed any further or even if this is the best approach..

How do I apply the logical condition (D<C) to the 3D matrices?

All help is greatly appreciated.

Kind regards

[Rik]

Can you post working code that does what you need? That way it is easier to see what you mean. You can use the code below to create the variables you describe.

AB1=rand(1250,441,451);
AB2=rand(1250,441,451);
C=rand(1,5);
D=rand(1,1250);
ABC=cell(size(C));%don't use dynamically numbered variable names

[Askeladden2]

Thank you Rik for the swift reply.

I dont want to use multiple loops, so I have not created a code. But I can give a simplified example instead;

AB1=4x3x3

AB2=4x3x3

C=1x4

D=1x3

AB1=[1 2 3;4 5 6; 7 8 9; 10 11 12]; AB2=[13 14 15; 16 17 18; 19 20 21; 22 23 24]; C=[1:1:4]; D=[0.5:1:2.5];

Third dimension of AB1 and AB2= [0.1:0.1:0.3]

Result:

ABC1=AB1

ABC2=[13 14 15; 4 5 6; 7 8 9;10 11 12] and third dimension 0.1 0.2 0.3

ABC3=[13 14 15;16 17 18; 7 8 9; 10 11 12] and third dimension 0.1 0.2 0.3.

Thanks in advance

[Bruno Luong]

Your example is incoherent with the question, you said in the question that length(D) matches size(AB1,1) and size(AB2,2).

But then in your example C is match, not D.

You also seem to reverse the comparison.

If you run this it gives the results you expect, but it's not exactly as with your question. Let you sort it out:

AB1=[1 2 3;4 5 6; 7 8 9; 10 11 12];
AB2=[13 14 15; 16 17 18; 19 20 21; 22 23 24]; 
D=[1:1:4];
C=[0.5:1:2.5];

k=double((D(:)<=C(:).'));
AB12=[AB1(:),AB2(:)];
n = size(AB12,1);
k = reshape(k,size(AB1,1),1,1,[]);
i = reshape(1:n,size(AB1));
ABC12=AB12(i+n*k);
ABC12_cell=num2cell(ABC12,[1 2 3]);
[ABC1,ABC2,ABC3]=deal(ABC12_cell{:})

[J. Alex Lee]

what does it mean that the third dimension is (0.1:0.1:0.3)? For your purposes does it matter at all?

Answer 1

[Bruno Luong]

This is code according to the specification stated in the question, NOT the example that is incoherent.

% Dummy test data
AB1=rand(250,441,451);
AB2=rand(250,441,451);
D=rand(1,size(AB1,1));
C=rand(1,5);

k=double((D(:)>C(:).'));
AB12=[AB1(:),AB2(:)];
n = size(AB12,1);
k = reshape(k,size(AB1,1),1,1,[]);
i = reshape(1:n,size(AB1));
ABC12=AB12(i+n*k);
ABC12_cell=num2cell(ABC12,[1 2 3]);
[ABC1,ABC2,ABC3,ABC4,ABC5]=deal(ABC12_cell{:});

[Askeladden2]

Yes, this is the solution I was looking for. I was too quick creating the example.

Thank you very much Bruno!

[Rik]

You shouldn't use numbered variables like that last line. That way you force yourself to never change the number of elements in C, because otherwise you'll have to write more code to account for that. If you use ABC12_cell instead, you can trivially loop over all elements.

[Bruno Luong]

Take Rik's advise seriously.

I would suggest using 4D array ABC12 rather than the cell or - worse of all - the separate variable names.

[J. Alex Lee]

+1 for neat vectorized code and >2D thinking (bow down)

[Askeladden2]

Thank you very much for the tip Rik!

Answer 2

[J. Alex Lee]

What is the purpose of vectorization? Is it for imagined speed gains, or just an "in principle" question?

Bruno's answer is fascinating in principle, but is undecipherable for the likes of me, but also signifcantly slower than loops (if I'm doing the right thing)

% Dummy test data
AB1=rand(250,441,451);
AB2=rand(250,441,451);
C=rand(1,5);
D=rand(1,size(AB1,1));

%% Bruno's answer
tic

k=double((D(:)>C(:).'));
AB12=[AB1(:),AB2(:)];
n = size(AB12,1);
k = reshape(k,size(AB1,1),1,1,[]);
i = reshape(1:n,size(AB1));
ABC12=AB12(i+n*k);
ABC12_cell=num2cell(ABC12,[1 2 3]);
[ABC1,ABC2,ABC3,ABC4,ABC5]=deal(ABC12_cell{:});

%% loops to achieve same thing

tic

ABCD = cell(size(C));
for k = 1:length(C)
    ABCD{k} = AB1;
    mask = D > C(k);
    ABCD{k}(mask,:,:) = AB2(mask,:,:);
end

toc

%% check answers between methods
%  as Rik points out, this is why naming variables serially is silly.
isequal(ABCD{1},ABC1)
isequal(ABCD{2},ABC2)
isequal(ABCD{3},ABC3)
isequal(ABCD{4},ABC4)
isequal(ABCD{5},ABC5)

The result is

Elapsed time is 11.659369 seconds.
Elapsed time is 3.062867 seconds.

[Bruno Luong]

+1, simple and efficient for-loop.

[Askeladden2]

For this specific example it was just speed gain imagination!

I also voted +1 for this simple and neat code.