w = 40;
h = 50;
n = 1000;
ch = rand( [h,w,n] );
mask = rand( [h,w] ) < .95;

disp('Shai')
tic;
r = reshape( ch, [], size(ch,3) );
Mean_res = mean( r(mask(:),: ), 2 );
toc

disp('Divakar 1')
tic;
Mean_bsxfun = sum(reshape(bsxfun(@times,ch,mask),[],size(ch,3)),1)./sum(mask(:));
toc
disp('Divakar 2')
tic;
Mean_matmult  = mask(:).'*reshape(ch,[],size(ch,3))./sum(mask(:));
toc