Interpolating scattered 3 dimensional data

I have measured electric field data in three dimensions of the following form:

 pos = [x y z]
 ef = [e_x e_y e_z]

The matrices are 1000x3 in size, and the positions are located in a half sphere (cartesian coordinates). I want to be able to interpolate the electric field at some point in space, so that I receive all the three values of the electric field components, not just the norm. interp3 won't work as the points are not in a grid. scatteredInterpolant needs the norm as the input, so it does not fit my needs. Any suggestions on what function to use?

[Aarti Dwivedi]

Just to be clear, you have data in cartesian coordinates, right?

[Markus Toivonen]

Yes, I do. I will add that information to the first post. Thanks.

[KSSV]

Scatteredinterpolant should work..why you think it will not work?

[Markus Toivonen]

Initially I thought that it would not work, because if it counts the components separately, it would provide the wrong answer. But this does not seem to be the case, and you can individually calculate them:

 e_x_new = scatterInterpolant(X,Y,Z,e(:,1),x,y,z)

and same for the two different components.

[KSSV]

Hope your problem solved..

Answer 1

[Anton Semechko]

Hi, Markus, here is a demo of how to perform linear interpolation of vector fields on a unit half-sphere:

half_sphere_interpolation_demo

function half_sphere_interpolation_demo
% Demo of how to perform linear interpolation of scalar and vector fields 
% defined on a unit half-sphere.

% PART 1: Simulate data sampled on a half-sphere 
% =========================================================================

% Generate a set of N randomly distributed points on the northern hemisphere
N=1E3;
X=RandSplHalfSphere(N);

% Triangulate points
x=StereoProj(X);
Tri=delaunay(x);
TR=triangulation(Tri,X);

% Generate vector field sampled at X
t=X(:,3);
t(t>1)=1;
t=acos(t);

f=@(t) sin([4*t 8*t 12*t]);
F=f(t);

figure('color','w')

TTL={'Fx' 'Fy' 'Fz'};
CLim=zeros(3,2);
for i=1:3
    ha=subtightplot(2,3,i,[0.1 0.05],[0.05 0.05]);
    h=trimesh(TR);
    set(h,'EdgeColor','k','EdgeAlpha',0.25,'FaceColor','interp','FaceVertexCData',F(:,i))
    axis equal
    colorbar(ha,'Location','SouthOutside')    
    set(get(ha,'Title'),'String',TTL{i},'FontSize',20)
    CLim(i,:)=get(ha,'CLim');
end
drawnow

% PART 2: Linear interpolation
% =========================================================================

% Generate new set of M points on a half-sphere
M=1E4;
Y=RandSplHalfSphere(M);

% Find faces of TR and containing Y  
[T,uvw]=SphericalTriangleIntersection(TR,Y);
id_val=~isnan(T); % points in Y that intersect with TR

% IMPORTANT: In this demo, TR does not completely cover the entire 
%            northern hemisphere, and thus it will not be possible to 
%            interpolate F at some of the points in Y (close 
%            to the equator) that do not intersect with TR. Indices 
%            corresponding to these points are isnan(T).

% Only retain points in Y that intersect with TR
Y=Y(id_val,:);
T=T(id_val,:);
uvw=uvw(id_val,:);

% Interpolate F at Y
T=Tri(T,:); % triangles containing Y
F_int=bsxfun(@times,uvw(:,1),F(T(:,1),:)) + ... 
      bsxfun(@times,uvw(:,2),F(T(:,2),:)) + ...
      bsxfun(@times,uvw(:,3),F(T(:,3),:));

% Visualize interpolation errors for each component of F
t=Y(:,3);
t(t>1)=1;
t=acos(t);
F_ref=f(t);
dF=F_int-F_ref;

TR2=triangulation(delaunay(StereoProj(Y)),Y);
TTL={'Fx error' 'Fy error' 'Fz error'};
for i=1:3
    ha=subtightplot(2,3,3+i,[0.1 0.05],[0.05 0.05]);
    h=trimesh(TR2);    
    set(h,'EdgeColor','k','EdgeAlpha',0.25,'FaceColor','interp','FaceVertexCData',dF(:,i))
    axis equal
    colorbar(ha,'Location','SouthOutside')
    set(get(ha,'Title'),'String',TTL{i},'FontSize',20)
    set(ha,'CLim',CLim(i,:))
    colormap(ha,'jet')
end

function [T,uvw]=SphericalTriangleIntersection(TR,P)
% Find barycentric coordinates of the points of intersection between a
% hemispherical mesh, TR, and set of rays, P. Note that rays emanating from
% the origin = points on a unit sphere.

% Use stereographic projection (with north pole as the origin) to 
% identify triangles intersected by the rays
x=StereoProj(TR.Points);
p=StereoProj(P);
tr=triangulation(TR.ConnectivityList,x);
[T,uvw]=pointLocation(tr,p);

function x=StereoProj(X)
% Stereographic projection from a unit sphere onto a plane tangent to the
% north pole. 

x=bsxfun(@rdivide,X(:,1:2),1+X(:,3));

function X=RandSplHalfSphere(N)
% Use rejection sampling to generate N random points on the northern 
% hemisphere 

X=[];
while size(X,1)<N
    Xi=randn(N,3);
    Xi(Xi(:,3)<0,:)=[];
    X=cat(1,X,Xi);
end
X=X(1:N,:);
X=bsxfun(@rdivide,X,sqrt(sum(X.^2,2)));  

% -----------------------------------------------------------------

'subtightplot' function used in this demo can be downloaded from here

% -----------------------------------------------------------------