Interpolating scattered 3 dimensional data
7 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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 comentarios
Respuestas (1)
Anton Semechko
el 4 de Jul. de 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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/189760/image.jpeg)
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 comentarios
Ver también
Categorías
Más información sobre Surface and Mesh Plots en Help Center y File Exchange.
Productos
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!