broadcast 2D 3D, bsxfun, , :
[m,n,r] = size(I);
out = I(bsxfun(@plus,linInd,(m*n*(0:r-1))'))
%// ---------------- 2D Case ---------------------
im = randi(9,10,10);
imSize = size(im);
rowPnts = [3,6,8,4];
colPnts = [6,3,8,5];
linInd = sub2ind(imSize,rowPnts,colPnts);
%// ---------------- 3D Case ---------------------
I = randi(9,10,10,4);
%// BSXFUN solution
[m,n,r] = size(I);
out = I(bsxfun(@plus,linInd,(m*n*(0:r-1))')); %//'
%// Tedious work of splitting
Ir = I(:,:,1);
Ig = I(:,:,2);
Ib = I(:,:,3);
Ia = I(:,:,4);
>> Ir(linInd)
ans =
8 9 1 6
>> Ig(linInd)
ans =
1 5 9 8
>> Ib(linInd)
ans =
8 5 3 8
>> Ia(linInd)
ans =
8 8 3 3
>> out
out =
8 9 1 6
1 5 9 8
8 5 3 8
8 8 3 3