MATLAB: Correct Interpolation on a pre-calculated grid that is a vector function of a vector

interpolation

I have calculated vector data, such that the data can be written and represented as , where each component is a function of all three coordinates; i.e. . I need to interpolate to be able to represent the inputs of a different vector anywhere, and have it be a vector output. This different vector is . Written out fully, you can see that each component of the vector M is an array which depends on all three x, y, and z coordinates. As a result, I can't simply use Interp3 three times. This means that I need the G's as a grid, so that I can query at points (Gx_q, Gy_q, Gz_q) and find their (x,y,z) inputs that give G.
Interp3 repeatedly gives errors of "Input is not a Meshgrid."
My vector data arrays (, , and ) are sized 31x31x31. I'm so lost; I've tried so many tricks that have failed. Any tricks? Will scatteredInterpolant() work?
The situation is analogous to this user's question (https://www.mathworks.com/matlabcentral/answers/408712-interpolating-scattered-3-dimensional-data), but distinctly different. The answers provided in () also are not quite related – something there is wrong, as well.
Previously, I've poorly or incorrectly explained my problem, and so those questions (and, as a result, answers) were phrased incorrectly. This is fixed.
To avoid confusion, I'll present the case where the M vector depends only on one coordinate, each, such that , which lets us write a scalar G as , and just do the process three times. The code would look like:
load('data.mat')
x = 0:1000:30000;
y = 0:1000:30000;
z = 0:1000:30000;
Gx_q = 721.0085; %These are example query points
Gy_q = 15.1313;
Gz_q = 45.9302;
Gx = x - squeeze(Mx(1,1,:))'; %These should be 1D arrays; for simplicity, in this case I would only need one column of each M_i array.
Gy = y - squeeze(My(1,1,:))';
Gz = z - squeeze(Mz(1,1,:))';
Gx_out = interp1(Gx,x,Gx_q,'linear','extrap');
Gy_out = interp1(Gy,y,Gy_q,'linear','extrap');
Gz_out = interp1(Gz,z,Gz_q,'linear','extrap');

Best Answer

so that I can query at points (Gx_q, Gy_q, Gz_q) and find their (x,y,z) inputs that give G.
But there may exist no such point or many of them. Suppose for example that Gx=1,Gy=1,Gz=1 everywhere in x,y,z. Then if Gx_q=5, Gx_q=6,Gz_q=7 there will exist no (x,y,z) agreeing with this query point for any reasonable interpolator. If, on the other hand Gx_q=1, Gx_q=1,Gz_q=1 then all (x,y,z) agree with it.
The best I think you can do is solve for one coordinate (x,y,z) at which the interpolation of G is closest to (Gx_q, Gy_q, Gz_q) according to some error metric. That can be approached as an inverse problem as in the following example:
[x,y,z]=ndgrid(0:1000:30000);
Gx=x-Mx;
Gy=y-My;
Gz=z-Mz;
%Find best grid point (x,y,z) - later refine that with fminsearch
Err = abs(Gx(:)-Gx_q)+abs(Gy(:)-Gy_q)+abs(Gz(:)-Gz_q);
[~,i]=min(Err);
xyzInitial = [x(i), y(i), z(i)];
%Refine:
Fx=griddedInterpolant(x,y,z,Gx);
Fy=griddedInterpolant(x,y,z,Gy);
Fz=griddedInterpolant(x,y,z,Gz);
xyz=fminsearch(@(xyz)lookupError(xyz,Fx,Fy,Fz),  xyzInitial);
x_out=xyz(1);
y_out=xyz(2);
z_out=xyz(3);
    function  Err=lookupError(xyz,Fx,Fy,Fz)
         gx=Fx(xyz);
         gy=Fy(xyz);
         gz=Fz(xyz);
        Err= abs(gx-Gx_q)+abs(gy-Gy_q)+abs(gz-Gz_q);
    end
Related Question