Correctly wrap data for spherical interpolation.

Say I've got n scattered data points defined in spherical coordinates,

n=36;
az=random('unif',0,2*pi,n,1);
ele=random('unif',-pi/2,pi/2,n,1);
r=random('unif',3,7,n,1);

I want to interpolate to get a finer sampling of r. I've tried doing this using the TriScatteredInterp() function,

F = TriScatteredInterp(az,ele,r);
n_i=50;
az_i=linspace(0,2*pi,n_i); ele_i=linspace(-pi/2,pi/2,n_i);
[AZ_i ELE_i]=meshgrid(az_i,ele_i);
R_i=F(AZ_i,ELE_i);
[X Y Z]=sph2cart(AZ_i,ELE_i,R_i);
surf(X,Y,Z)

As the code is above, I get serious problems near az=0, and ele=pi/2. I believe this is because the data doesn't "wrap around" as it should in spherical coordinates. Thus the interpolation near az=0 is off and I get nan's at these locations. I tried to fix this by wrapping the data as follows.

n=10;
az=random('unif',0,2*pi,n,1); az=[az-2*pi; az; az+2*pi]; %wrap az
ele=random('unif',-pi/2,pi/2,n,1); ele=repmat(ele,[3 1]);%copy ele
r=random('unif',3,7,n,1); r=repmat(r,[3 1]);             %copy r
F = TriScatteredInterp(az,ele,r);
n_i=50;
az_i=linspace(0,2*pi,n_i); ele_i=linspace(-pi/2,pi/2,n_i);
[AZ_i ELE_i]=meshgrid(az_i,ele_i);
R_i=F(AZ_i,ELE_i);
[X Y Z]=sph2cart(AZ_i,ELE_i,R_i);
surf(X,Y,Z)

All I've done above is copied the data and changed the ele values on the copied version to go from -2*pi-0 and 2*pi04*pi. This improves the interpolation drastically however I still get nan's near the z-axis (i.e., where ele=pi/2). Basically, this puts a big hole in my data. Has anyone run into a similar problem. Is there a better "wrapping" technique I should use? Thanks in advance for any insight.

Justin

Answer 1

[Sven]

Justin,

You're in luck - I've had the exact same problem and worked reasonably hard to get an acceptable solution. See the code below... all I've changed is the variable names from to azimuth and elevation to theta and phi... I just had other variables with that naming convention so it all looked nicer.

The trick is grab the highest/lowest elevation point, and copy it across the "top" and "bottom" of the phi/theta space (ie, the poles). You'll see in the first figure which points are copied to where. The "best case scenario" is to have a point ON the poles (in which case no error at all is introduced when that point is copied to various THETA locations), but any point reasonably near to the pole should do just fine.

n=36;
sphSampTh = random('unif',0,2*pi,n,1);
sphSampPhi = random('unif',-pi/2,pi/2,n,1);
r=random('unif',3,7,n,1);

sphSampTh = mod(sphSampTh,2*pi); % Interpret theta from -pi:pi to 0:2pi
% "Wrap" the lateralmost thetas around to the opposite side
wrapSz = pi/4;
wrapThIdxs = find(abs(sphSampTh-pi)>(pi-wrapSz) & abs(sphSampTh-pi)<pi);
wrappedThs = sphSampTh(wrapThIdxs) - 2*pi*sign(sphSampTh(wrapThIdxs)-pi);
% Force the samples taken at the poles to be copied all along the thetaVec but still at the poles. 
maxPhiThVec = linspace(0,2*pi,20)';
[~,highPhiIdx] = max(sphSampPhi);
[~,lowPhiIdx] = min(sphSampPhi);
wrapPhiThIdxs = [maxPhiThVec repmat([pi/2 highPhiIdx],size(maxPhiThVec))];
wrapPhiThIdxs = cat(1, wrapPhiThIdxs, [maxPhiThVec repmat([-pi/2 lowPhiIdx],size(maxPhiThVec))]);
% Get a new set of ThPhi combinations that include these additional points, along with indices into
% the sampled coverage map.
gloThPhis = cat(1, [sphSampTh sphSampPhi], [wrappedThs sphSampPhi(wrapThIdxs)], wrapPhiThIdxs(:,1:2));
gloCoverageIdxs = cat(1, (1:n)', wrapThIdxs, wrapPhiThIdxs(:,3));

% NOW we can interp()
DT = DelaunayTri(gloThPhis(:,1),gloThPhis(:,2));
TSinterp = TriScatteredInterp(DT, r(gloCoverageIdxs));
sampThVec = 0:pi/60:2*pi; sampPhiVec = -pi/2:pi/60:pi/2;
[sampThMat, sampPhiMat] = meshgrid(sampThVec, sampPhiVec);
interpdR = TSinterp(sampThMat, sampPhiMat);

% Show the result(s)
figure, plot(sphSampTh,sphSampPhi,'b.',wrappedThs,sphSampPhi(wrapThIdxs),'ro',wrapPhiThIdxs(:,1), wrapPhiThIdxs(:,2),'go')
hold on, rectangle('Position',[0 -.5 2 1]*pi)
title('Scattered spherical sample locations')
xlabel 'Theta', ylabel 'Phi'

figure, imagesc(sampThVec,sampPhiVec,interpdR)

[X Y Z] = sph2cart(sampThMat,sampPhiMat,interpdR);
figure, surf(X,Y,Z)

[Justin Solomon]

Thanks a lot for your response! I took some of your ideas and implemented them in my code. The main one being the fact that you have to copy the same r value across all az values where ele is -pi/2 or pi/2.

Also, wrt wrapping, I found that in order to avoid artifacts along the z-axis, you should not only tile r in the az axis,

(az= [az-2*pi az az+2*pi]; r=[r r r])

but you should also "reflect" r above and below in the ele direction

(ele=[ele(1:end-1)-pi ele ele(2:end)+pi]; r=[flipud(r(2:end,:)) ;r ; flipud(r(1:end-1,:))];

Thus we have the correctly copied the data so the interpolation on the "edges" (i.e., where az= 0 or 2*pi and where ele=-pi/2 or pi.2) wraps around. Anyways, thanks again for your help! Justin

[Sven]

No problems, glad it helped.

And yes, "reflection" is needed in your case. In my original data, I actually had a point sitting directly at the poles, so my copying along the top/bottom was enough to cover the required area without the need for subsequent reflection.

Glad it's all sorted.