# Interpolating scattered 3 dimensional data

11 vues (au cours des 30 derniers jours)
Markus Toivonen le 4 Juil 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 commentairesAfficher 3 commentaires plus anciensMasquer 3 commentaires plus anciens
Markus Toivonen le 4 Juil 2018
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 le 4 Juil 2018

Connectez-vous pour commenter.

### Réponses (1)

Anton Semechko le 4 Juil 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)));
% -----------------------------------------------------------------
% -----------------------------------------------------------------
##### 0 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

### Catégories

En savoir plus sur Surface and Mesh Plots dans Help Center et File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by