球坐标中的插值曲线

Question

我有点,属于3D空间中的一些非参数曲线.我需要在球坐标中工作并在给定azimuth和elevation值中插入给定的曲线.

你能用这种插值算法帮助我吗？

PS我不确定,从我在1D空间中绘制的一维曲线的问题想法是明确的.所以,以防下面的代码创建这样的,但参数化.

x=linspace(-2*pi,2*pi,10^3);
y=sin(x);
z=sinh(x);

更新:

感谢@jodag关于上采样的非常好的帖子.slerp是完美的!然而,我的问题是不同的,或者我不明白如何使用上采样来解决它.假设我将我的点定义为数组(只是一些虚拟数据用于演示):

x_given = rand(9, 1) - 0.5; 
y_given = rand(9, 1) - 0.5; 
z_given = rand(9, 1) - 0.5; 
x_given(end+1) = x_given(1); % should work for closed curves
y_given(end+1) = y_given(1); % should work for closed curves
z_given(end+1) = x_given(1); % should work for closed curves

我想找到x_i, y_i, z_i给定的所有点()elevation_q和azimuth_q:

elevation_q = pi./4;
azimuth_q = pi./7;
[x_i, y_i ,z_i] = interpolateCurve(x_given, y_given, y_given, elevation_q, azimuth_q);

这就是我真正需要的.我明白,我们不需要前后从笛卡尔坐标到球面坐标.但是,即使我将取代x_given, y_given, y_given于elevation_given, azimuth_given, r_given,我不知道怎么写interpolateCurve()功能.

Answer 1

这是高维空间中曲线的线性插值的示例.请注意,我引入了t与点的索引相对应的参数化变量.

function [xq,yq,zq] = interp_sph(x,y,z,upsample)
    [az,el,r] = cart2sph(x,y,z);
    t = 1:numel(x);
    tq = 1:(1/upsample):numel(x);
    azq = interp1(t,az,tq);
    elq = interp1(t,el,tq);
    rq = interp1(t,r,tq);
    [xq,yq,zq] = sph2cart(azq,elq,rq);
end

从问题中不清楚这是否是您实际需要的.在球坐标中存在直接插值的问题.例如,如果两个点彼此接近但一个方位角= -179度而另一个方位角= 179而不是采用最短的2度路径,则插值将产生跨越[-179,179]的样本.您可能感兴趣的是slerp,它提供了更好的插值结果.

我写了以下slerp示例来演示.请注意,slerp要求1/sin(theta)定义theta函数中两点之间的角度.如果theta == pi那时我们遇到问题并且slerp将无法工作,那么请注意您的数据是否接近原点.

function [xq,yq,zq] = interp_slerp(x,y,z,upsample)
    r = sqrt(x.^2+y.^2+z.^2);
    vq = zeros(3,(upsample-1)*(numel(x)-1)+1);
    for idx = 1:numel(x)-1
        t = 0:1;
        tq = linspace(0,1,upsample);
        % compute interpolated direction
        dq = slerp([x(idx),y(idx),z(idx)],[x(idx+1),y(idx+1),z(idx+1)],tq);
        % linearly interpolate radius
        rq = interp1(t,r(idx:(idx+1)),tq);
        % get interpolated path
        idxq = (idx-1)*upsample + 1;
        vq(:,idxq:(idxq+(upsample-1))) = bsxfun(@times,rq,dq);
    end
    xq = vq(1,:);
    yq = vq(2,:);
    zq = vq(3,:);
end

function pq = slerp(p0,p1,t)
    p0 = p0(:) / norm(p0);
    p1 = p1(:) / norm(p1);
    theta = acos( dot(p0,p1) );
    if(0==theta)
        pq = a;
    elseif(pi==theta)
        error('Angle between points cannot be exactly pi.');
    else
        pq = bsxfun(@times,(sin((1.0-t)*theta)/sin(theta)),p0) + bsxfun(@times,(sin(t*theta)/sin(theta)),p1);
    end
end

要演示此问题,请interp_sph考虑以下示例.

% create test path
az = pi+sin(linspace(pi/2,3*pi/2,10));
el = cos(linspace(-pi,pi,10));
r = 1 + 1/5*sin(6*az).*sin(5*el);
[x,y,z] = sph2cart(az,el,r);

% test both spherical interpolation and slerp
upsample = 20;
[xq_sph,yq_sph,zq_sph] = interp_sph(x,y,z,upsample);
[xq_slerp,yq_slerp,zq_slerp] = interp_slerp(x,y,z,upsample);

% plot results
plot3(x,y,z,'LineWidth',2); hold on;
plot3(xq_sph,yq_sph,zq_sph,'LineWidth',2);
plot3(xq_slerp,yq_slerp,zq_slerp,'LineWidth',2);
axis('vis3d');
axis([-1,1,-1,1,-1,1]);
grid('on');
legend('Cartisian', 'Spherical', 'SLERP');

这将产生下图.请注意,当方位角从正变为负时,球形插值需要很长时间,但是slerp会采用预期的路径.