MATLAB: Find index on 3d matrix

findgriddatamesh

Hi, I've a grid
[X,Y,Z] = meshgrid(XI,YI,ZI);
where
XI = 1:1:4;
YI = 1:1:4;
ZI = -500:100:5000;
and a 3d surface defined by points xi, yi, zi
[A,B] = meshgrid(XI,YI);
S = griddata(xi,yi,zi,A,B);
I've to performe simple index calculations (linear extrapolation of temperature with depth) on my grid [X,Y,Z] but with different coefficients depending on the node i,j,k position: above or below the surface S. I think I've to find k index of my grid that are above and below the surface.
for i = 1:4
for j = 1:4
for k = 1:(length(Z)-1)
ind(i,j,k) = find(Z(i,j,k) >= S)
end
end
end
The problem is that the dimension of the grid and surface matrix are different! Any idea how to solve this problem?
Gianluca

Best Answer

Hi Gianluca,
The MATLAB docs say "TriScatteredInterp is the recommended alternative to griddata as it is generally more efficient", so I will use the TriScatteredInterp example.
What you've got is some surface S, defined by points that aren't necessarily on your [XI,YI] grid. If you can interpolate your surface S at the [XI,YI] grid locations, then you can get to the next part of your question. So let's do that interpolation:
% Create a data set:
x = rand(100,1)*4 - 2;
y = rand(100,1)*4 - 2;
S = x.*exp(-x.^2-y.^2) * 1000;
% Construct the interpolant:
F = TriScatteredInterp(x,y,S);
% Evaluate the interpolant at the locations [XI,YI].
XI = -2:0.25:2;
YI = -2:0.1:2;
[XImat,YImat] = meshgrid(XI,YI);
ZImat = F(XImat,YImat);
Now, you've got a variable ZImat that lines up perfectly with your [XI,YI] grid and represents your original surface S. Next, you've got a set of ZI locations:
ZI = -500:100:500;
Let's make a logical matrix showing which [XI,YI,ZI] grid locations are above your surface S (which is now represented by ZImat):
BW = false(length(YI),length(XI),length(ZI));
for i = 1:length(YI)
    for j = 1:length(XI)
        BW(i,j,:) = ZImat(i,j)>ZI;
    end
end
And, just because it's always interesting to see a different way of doing the same thing, here's a very efficient, vectorised "single-line" way of getting the same BW matrix:
BW = bsxfun(@gt, ZImat, reshape(ZI,1,1,[]));
Now, to "find the first k index of the grid that is above the surface", you can either go with a loop:
KImat = zeros(size(ZImat));
for i = 1:length(YI)
    for j = 1:length(XI)
        firstIndex = find(BW(i,j,:),1);
        if ~isempty(firstIndex)
            KImat(i,j) = firstIndex;
        end
    end
end
Or, with the help of a file-exchange entry find_ndim
KImat = find_ndim(BW,3,'first');
Does this help you out?
Related Question