MATLAB: Correctly wrap data for spherical interpolation.

spherical interpolationwrapping

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

Best Answer

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)
Related Question