Logical Indexing for a slice of a 3D-Matrix

I have two matrices, on which I do some work:

% MyMat1: n x m x 1
% MyMat2: n x m x 3

This is what I do until now:

ValsToBeProcessed = (MyMat1 > 5);

MyMat2_1 = MyMat2(:,:,1);
MyMat2_2 = MyMat2(:,:,2);
MyMat2_3 = MyMat2(:,:,3);

MyMat2_2(ValsToBeProcessed) = MyMat1(ValsToBeProcessed);

MyMat2 = cat(3, MyMat2_1, MyMat2_2, MyMat2_3);

Is there any simpler way to do this? Something like:

ValsToBeProcessed = (MyMat1 > 5);

MyMat2(ValsToBeProcessed, 2) = MyMat1(ValsToBeProcessed);

Answer 1

[Lemar DeSalis]

Thanks for all your answers!

In the meantime I found a different solution:

ValsToBeProcessed3D = false(size(MyMat2));
ValsToBeProcessed3D(:,:,2) = ValsToBeProcessed;

MyMat2(ValsToBeProcessed3D) = MyMat1(ValsToBeProcessed);

This is more "readable" to me, so I'll use that version (although your methods may be faster).

[Sean de Wolski]

That's actually a winner in timings :) Congrats!

+1

[Teja Muppirala]

Nice work. Easy to read, and efficient

Answer 2

[Sean de Wolski]

idx = find(MyMat1>5);  %indices in 2d matrix
MyMat2(idx+numel(MyMat1)) = MyMat1(idx);
%indices in 3d matrix and extracted from 2d matrix

NOTE:

%I can't guarantee that this will be any faster than your original way compacted to:
slice2 = MyMat2(:,:,2)
idx = MyMat1>5;
slice2(idx) = MyMat1(idx);
MyMat2(:,:,2) = slice2;

Answer 3

[Andrei Bobrov]

bsxfun(@times,MyMat2,(MyMat1>5));

EDIT

M5 = MyMat1>5;
MyMat2(find(M5)+numel(MyMat2)/3)=MyMat1(M5);

[Sean de Wolski]

They only want the second slice, not all three. And you'd have to convert the logical expression to the class of MyMat for the multiplication op to work.

cast(MyMat1>5,class(MyMat2))

[Andrei Bobrov]

Hi Sean de! Excuse me, my carelessness, correct...

Answer 4

[Teja Muppirala]

Probably a bit faster:

IDX = cast(MyMat1>5,class(MyMat1));
MyMat2(:,:,2) = MyMat1.*IDX+MyMat2(:,:,2).*(1-IDX);

Answer 5

[Sean de Wolski]

Obligatory timings :)

MyMat2store = rand(1000,1000,3);
MyMat1 = rand(1000)*10;
t1 = 0;
t2 = 0;
t3 = 0;
for ii = 1:30;
  MyMat2 = MyMat2store;
  tic
  idx = find(MyMat1>5);
  MyMat2(idx+numel(MyMat1)) = MyMat1(idx);
  t1 = t1+toc;
  MyMat2 = MyMat2store;
  tic
  slice2 = MyMat2(:,:,2);
  idx = MyMat1>5;
  slice2(idx) = MyMat1(idx);
  MyMat2(:,:,2) = slice2;
  t2 = t2+toc;
  MyMat2 = MyMat2store;
  tic
  IDX = cast(MyMat1>5,class(MyMat1));
  MyMat2(:,:,2) = MyMat1.*IDX+MyMat2(:,:,2).*(1-IDX);
  t3 = t3+toc;
end
disp([t1 t2 t3]);
%with MyMat1 = rand(1000)*10
%   2.9776       2.8133       2.6096
%   2.9078        2.772       2.5669
%with MyMat1 = rand(1000)*7
%   2.4349       2.5948       2.6454
%   2.4591       2.6174        2.724
%with MyMat1 = rand(1000)*20
%   3.5691       3.0149       2.6538
%   3.5123       2.9869       2.6282
%with MyMat1 = rand(1000)*5.5
%   1.7049       2.1195       2.6007
%   1.7177       2.1276        2.604

So Teja's method is relatively constant regardless of the percentage of MyMat1 >5. The two extraction methods that I proposed are faster when the percentage of MyMat1>5 is small and slower when that percentage is large. These findings make sense considering the operations involved.

Answer 6

[Teja Muppirala]

After just trying it out of curiousity, I find that a simple FOR loop actually works better than all of these vectorized solutions (though again it depends a little on the composition of MyMat1) The JIT has really become incredibly good over the last couple of years.

N = numel(MyMat1);
for n = 1:N
    if MyMat1(n) > 5
        MyMat2(N+n) = MyMat1(n);
    end
end

[Jan Siegmund]

Loops work wonderful on double. Try using a self written numeric classes and loops are devastatingly slow.