Interpolating scattered 3 dimensional data

10 次查看（过去 30 天）
Markus Toivonen
Markus Toivonen 2018-7-4
回答： Anton Semechko 2018-7-4
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 个评论
显示 3更早的评论
Markus Toivonen
Markus Toivonen 2018-7-4
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Anton Semechko
Anton Semechko 2018-7-4
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
% -----------------------------------------------------------------

类别

在 Help Center 和 File Exchange 中查找有关 Surface and Mesh Plots 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by