MATLAB Answers

0

Interpolating scattered 3 dimensional data

Asked by Markus Toivonen on 4 Jul 2018
Latest activity Answered by Anton Semechko on 4 Jul 2018
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?

  5 Comments

KSSV
on 4 Jul 2018
Scatteredinterpolant should work..why you think it will not work?
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
on 4 Jul 2018
Hope your problem solved..

Sign in to comment.

1 Answer

Answer by Anton Semechko on 4 Jul 2018

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
% -----------------------------------------------------------------

  0 Comments

Sign in to comment.