find index on 3d matrix

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

Answer 1

[Sven]

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?

[gianluca]

thanks Sven, I try to do your advices